PLUTO
Functions
init.c File Reference

Reconnection test (Harris sheet) in 2D. More...

#include "pluto.h"
Include dependency graph for init.c:

Go to the source code of this file.

Functions

void Init (double *v, double x, double y, double x3)
 
void Analysis (const Data *d, Grid *grid)
 
void BackgroundField (double x1, double x2, double x3, double *B0)
 
void UserDefBoundary (const Data *d, RBox *box, int side, Grid *grid)
 

Detailed Description

Reconnection test (Harris sheet) in 2D.

In this setup - see also Section 5.3 of Mignone et al., ApJS (2012) 198:7 - we reproduce a 2D Harris current sheet with magnetic field profile given by

\[ B_x(y) = B_0 \tanh(y/l) \]

where l is the half thickness of the layer. The density profile is given by

\[ \rho(y) = \rho_0 \cosh^{-2}(y/l) + \rho_{\infty} \]

We use $ \rho_0=1 $ and $ \rho_{\infty} = 0.2 $, following the guidelines of Birn et al., 2001, while l is user supplied.
In order to achieve equilibrium with the magnetic pressure, the thermal pressure is chosen to be $ p = c_s^2 \rho $, where $ c_s^2 = \frac{B_0^2}{2\rho_0} $. The initial equilibrium is pertubed by an additional magnetic field defined as

\[ \begin{array}{lcl} B_x(x,y) &=& \DS -\Psi_0\frac{\pi}{L_y}\cos\left(\frac{2\pi x}{L_x}\right) \sin\left(\frac{\pi y}{L_y}\right), \\ \noalign{\medskip} B_y(x,y) &=& \DS +\Psi_0 \frac{2\pi}{L_x}\sin\left(\frac{2\pi x}{L_x}\right) \cos\left(\frac{\pi y}{L_y}\right). \end{array} \]

The Lundquist number $ S $ of a plasma is defined as

\[ S = \frac{v_A L}{\eta} \]

where $ v_A $ is the Alfvén velocity, $ v_A = \DS \frac{B}{\sqrt{\rho}}$, $ L $ is a typical lenght scale, and $ \eta $ the plasma resistivity. The reconnection rate $\mathcal{E} = \DS \frac{v_{in}}{v_{out}}$, with $ v_{in} $ and $ v_{out}$ the plasma inflow and outflow velocities, follows the Sweet-Parker scaling law $\mathcal{E} \sim \frac{\delta}{L} \sim \frac{1}{\sqrt{S}}$. In this example several values of the resitivity $ \eta $, that correspond to different values of the Lundquist number $ S $, are provided. The reconnection rate, calculated as the ratio $ \frac{\delta}{L} $ (see Mignone et al., 2012) verifies the Sweet-Parker scaling in the range $ \eta = 10^{-2} - 10^{-4} $ (see the first figure below).

The input parameters (read from pluto.ini) for this test problem are:

  • g_inputParam[ETA]: sets the value of resistivity $ \eta $;
  • g_inputParam[WIDTH]: sets the layer width l;
  • g_inputParam[PSI0]: sets the amplitude of perturbation $\Psi_0$.
Note
  • Configuration #02 employs a small width (l -> 0, current-sheet) large resistivity (test passes only with the new implementation of the resistive-CT module in PLUTO 4.1. Crash with PLUTO 4.0).
  • Configuratation #09 employs adaptive mesh refinement as in the original PLUTO-Chombo paper (ApJS 2012).
fig1.png
Computed Sweet-Parker scaling for different values of eta with a resolution of 512x256.
vx_rho_plot__137.png
Density map and magnetic field lines for eta = 2.e-3 at t = 137.
Authors
E. Striani (edoar.nosp@m.do.s.nosp@m.trian.nosp@m.i@ia.nosp@m.ps.in.nosp@m.af.i.nosp@m.t)
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
July 29, 2014

Reference

  • "The PLUTO Code for Adaptive Mesh Computations in Astrophysical FLuid Dynamics" Mignone et al., ApJS (2012) 198,7
  • "Geospace Environmental Modeling (GEM) magnetic reconnection challenge" Birn et al., JGR (2001) 106, 3715

Definition in file init.c.

Function Documentation

void Analysis ( const Data *  d,
Grid *  grid 
)

Perform runtime data analysis.

Parameters
[in]dthe PLUTO Data structure
[in]gridpointer to array of Grid structures

Definition at line 113 of file init.c.

115 {
116 
117 }
void BackgroundField ( double  x1,
double  x2,
double  x3,
double *  B0 
)

Definition at line 120 of file init.c.

122 {
123  B0[0] = 0.0;
124  B0[1] = 0.0;
125  B0[2] = 0.0;
126 }
void Init ( double *  v,
double  x,
double  y,
double  x3 
)

The Init() function can be used to assign initial conditions as as a function of spatial position.

Parameters
[out]va pointer to a vector of primitive variables
[in]x1coordinate point in the 1st dimension
[in]x2coordinate point in the 2nd dimension
[in]x3coordinate point in the 3rdt dimension

The meaning of x1, x2 and x3 depends on the geometry:

\[ \begin{array}{cccl} x_1 & x_2 & x_3 & \mathrm{Geometry} \\ \noalign{\medskip} \hline x & y & z & \mathrm{Cartesian} \\ \noalign{\medskip} R & z & - & \mathrm{cylindrical} \\ \noalign{\medskip} R & \phi & z & \mathrm{polar} \\ \noalign{\medskip} r & \theta & \phi & \mathrm{spherical} \end{array} \]

Variable names are accessed by means of an index v[nv], where nv = RHO is density, nv = PRS is pressure, nv = (VX1, VX2, VX3) are the three components of velocity, and so forth.

Definition at line 82 of file init.c.

85 {
86 #if 1
87  double cs2 = 0.5, b0 = 1.0, l, Psi0;
88  double Lx, Ly, kx, ky;
89 
90  l = g_inputParam[WIDTH];
91  v[RHO] = 0.2 + 1.0/(cosh(y/l)*(cosh(y/l)));
92  v[PRS] = cs2*v[RHO]; /* v[PRS] = 0.5 in the original PLUTO (2012) paper. */
93  v[VX1] = 0.0;
94  v[VX2] = 0.0;
95 
96  #if PHYSICS == MHD || PHYSICS == RMHD
97  v[BX1] = b0*tanh(y/l);
98  v[BX2] = 0.0;
99 
100  Lx = g_domEnd[IDIR] - g_domBeg[IDIR]; kx = CONST_PI/Lx;
101  Ly = g_domEnd[JDIR] - g_domBeg[JDIR]; ky = CONST_PI/Ly;
102 
103 
104  Psi0 = g_inputParam[PSI0];
105  v[BX1] += -Psi0*ky*sin(ky*y)*cos(2.0*kx*x);
106  v[BX2] += Psi0*2.0*kx*sin(2.0*kx*x)*cos(ky*y);
107  #endif
108 #endif
109 
110 
111 }
VX2
#define VX2
Definition: mod_defs.h:29
RHO
#define RHO
Definition: mod_defs.h:19
PSI0
#define PSI0
Definition: definitions_01.h:28
VX1
#define VX1
Definition: mod_defs.h:28
IDIR
#define IDIR
Definition: pluto.h:193
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
g_domBeg
double g_domBeg[3]
Lower limits of the computational domain.
Definition: globals.h:125
GRID::x
double * x
Definition: structs.h:80
WIDTH
#define WIDTH
Definition: definitions_01.h:27
BX1
#define BX1
Definition: mod_defs.h:25
CONST_PI
#define CONST_PI
.
Definition: pluto.h:269
BX2
#define BX2
Definition: mod_defs.h:26
g_domEnd
double g_domEnd[3]
Upper limits of the computational domain.
Definition: globals.h:126
JDIR
#define JDIR
Definition: pluto.h:194
void UserDefBoundary ( const Data *  d,
RBox *  box,
int  side,
Grid *  grid 
)

Assign user-defined boundary conditions.

Parameters
[in,out]dpointer to the PLUTO data structure containing cell-centered primitive quantities (d->Vc) and staggered magnetic fields (d->Vs, when used) to be filled.
[in]boxpointer to a RBox structure containing the lower and upper indices of the ghost zone-centers/nodes or edges at which data values should be assigned.
[in]sidespecifies the boundary side where ghost zones need to be filled. It can assume the following pre-definite values: X1_BEG, X1_END, X2_BEG, X2_END, X3_BEG, X3_END. The special value side == 0 is used to control a region inside the computational domain.
[in]gridpointer to an array of Grid structures.

Assign user-defined boundary conditions in the lower boundary ghost zones. The profile is top-hat:

\[ V_{ij} = \left\{\begin{array}{ll} V_{\rm jet} & \quad\mathrm{for}\quad r_i < 1 \\ \noalign{\medskip} \mathrm{Reflect}(V) & \quad\mathrm{otherwise} \end{array}\right. \]

where $ V_{\rm jet} = (\rho,v,p)_{\rm jet} = (1,M,1/\Gamma)$ and M is the flow Mach number (the unit velocity is the jet sound speed, so $ v = M$).

Assign user-defined boundary conditions:

  • left side (x beg): constant shocked values
  • bottom side (y beg): constant shocked values for x < 1/6 and reflective boundary otherwise.
  • top side (y end): time-dependent boundary: for $ x < x_s(t) = 10 t/\sin\alpha + 1/6 + 1.0/\tan\alpha $ we use fixed (post-shock) values. Unperturbed values otherwise.

Assign user-defined boundary conditions at inner and outer radial boundaries. Reflective conditions are applied except for the azimuthal velocity which is fixed.

Assign user-defined boundary conditions.

Parameters
[in/out]d pointer to the PLUTO data structure containing cell-centered primitive quantities (d->Vc) and staggered magnetic fields (d->Vs, when used) to be filled.
[in]boxpointer to a RBox structure containing the lower and upper indices of the ghost zone-centers/nodes or edges at which data values should be assigned.
[in]sidespecifies on which side boundary conditions need to be assigned. side can assume the following pre-definite values: X1_BEG, X1_END, X2_BEG, X2_END, X3_BEG, X3_END. The special value side == 0 is used to control a region inside the computational domain.
[in]gridpointer to an array of Grid structures.

Set the injection boundary condition at the lower z-boundary (X2-beg must be set to userdef in pluto.ini). For $ R \le 1 $ we set constant input values (given by the GetJetValues() function while for $ R > 1 $ the solution has equatorial symmetry with respect to the z=0 plane. To avoid numerical problems with a "top-hat" discontinuous jump, we smoothly merge the inlet and reflected value using a profile function Profile().

Definition at line 130 of file init.c.

132 {
133  int i, j, k, nv;
134  double *x1, *x2, *x3;
135 
136  x1 = grid[IDIR].x;
137  x2 = grid[JDIR].x;
138  x3 = grid[KDIR].x;
139 
140  if (side == 0) { /* -- check solution inside domain -- */
141  DOM_LOOP(k,j,i){}
142  }
143 
144 }
DOM_LOOP
DOM_LOOP(k, j, i)
Definition: analysis.c:22
KDIR
#define KDIR
Definition: pluto.h:195
IDIR
#define IDIR
Definition: pluto.h:193
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
i
int i
Definition: analysis.c:2
JDIR
#define JDIR
Definition: pluto.h:194

Here is the call graph for this function: