PLUTO
Functions
init.c File Reference

Diffusion of a linear force-free magnetic field in cylindrical coordinates. More...

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

Go to the source code of this file.

Functions

void Init (double *us, double x1, double x2, double x3)
 
void UserDefBoundary (const Data *d, RBox *box, int side, Grid *grid)
 
void Analysis (const Data *d, Grid *grid)
 

Detailed Description

Diffusion of a linear force-free magnetic field in cylindrical coordinates.

Solve the diffusion equation for a linear force-free magnetic field in cylindrical coordinates. A force-free field satisfies

\[ \vec{J} \times \vec{B} = 0 \qquad\Longrightarrow\qquad \vec{J} = \nabla\times\vec{B} = \mu \vec{B} \]

which means that $\vec{J}$ and $\vec{B}$ must be parallel. If $ \mu $ is constant then the solution is given by the Bessel functions:

\[ B_z(r) = B_0 J_1(r) \,,\qquad B_\phi(r) = B_0 J_0(r) \]

For constant resistivity and zero velocity the induction equation simplifies as follows:

\[ \frac{d\vec{B}}{dt} = - \nabla\times (\eta \vec{J}) = - \nabla\times (\eta \mu \vec{B}) = -\eta \mu^2 \vec{B} \]

which admits the exact analytical solution $\vec{B}(r,t) = \vec{B}(r,0) \exp(-\eta\mu^2 t)$ meaning that the field remains force free also at subsequent times.

Note that if pressure and density are initially constant and the velocity is also initially zero everywhere, the previous solution is also an exact solution of the isothermal MHD equations but not of the adiabatic MHD equations because of the Ohmic dissipation term (magnetic energy transforms into heat). However, using a large density makes pressure effects negligible.

  • Configurations #01-03 solve the problem on a cylindrical domain away from the origin.
  • Configurations #04-06 extends the integration up to the origin and tests the quality of the discretization at r = 0
Author
T. Matsakos
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
April 13, 2014

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 149 of file init.c.

154 {
155 }
void Init ( double *  us,
double  x1,
double  x2,
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 52 of file init.c.

56 {
57  us[RHO] = g_inputParam[DENSITY];
58  us[VX1] = 0.0;
59  us[VX2] = 0.0;
60  us[VX3] = 0.0;
61  #if HAVE_ENERGY
62  us[PRS] = 1.0;
63  #endif
64  us[TRC] = 0.0;
65 
66  us[BX1] = 0.0;
67  us[BX2] = BesselJ0(x1);
68  us[BX3] = BesselJ1(x1);
69  #ifdef GLM_MHD
70  us[PSI_GLM] = 0.0;
71  #endif
72 }
VX2
#define VX2
Definition: mod_defs.h:29
DENSITY
#define DENSITY
Definition: definitions_01.h:26
BesselJ0
double BesselJ0(double x)
Definition: math_misc.c:13
RHO
#define RHO
Definition: mod_defs.h:19
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
TRC
#define TRC
Definition: pluto.h:581
VX1
#define VX1
Definition: mod_defs.h:28
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
BesselJ1
double BesselJ1(double x)
Definition: math_misc.c:45
BX3
#define BX3
Definition: mod_defs.h:27
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
BX2
#define BX2
Definition: mod_defs.h:26

Here is the call graph for this function:

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 74 of file init.c.

78 {
79  int i, j, k, nv;
80  double *r, e_t;
81  double t, dt, vin[256];
82  static float *j0, *j1; /* store bessel function to avoid overhead */
83 
84  t = g_time;
85  r = grid[IDIR].x;
86  e_t = exp(-t*1.0);
87 
88 /* -------------------------------------------------
89  Pre-compute Bessel functions to avoid overhead
90  ------------------------------------------------- */
91 
92  if (j0 == NULL){
93  j0 = ARRAY_1D(NX1_TOT, float);
94  j1 = ARRAY_1D(NX1_TOT, float);
95  ITOT_LOOP(i) {
96  j0[i] = BesselJ0(r[i]);
97  j1[i] = BesselJ1(r[i]);
98  }
99  }
100 
101 /* ----------------------------------------------------
102  Set default boundary values
103  ---------------------------------------------------- */
104 
105  vin[RHO] = g_inputParam[DENSITY];
106  vin[VX1] = vin[VX2] = vin[VX3] = 0.0;
107  vin[BX1] = vin[BX2] = vin[BX3] = 0.0;
108  #if HAVE_ENERGY
109  vin[PRS] = 1.0;
110  #endif
111  #ifdef GLM_MHD
112  vin[PSI_GLM] = 0.0;
113  #endif
114 
115 /* ---------------------------------------------------
116  Set boundary conditions
117  --------------------------------------------------- */
118 
119  if (side == X1_BEG){ /* -- X1_BEG boundary -- */
120  if (box->vpos == CENTER){
121  BOX_LOOP(box,k,j,i){
122  VAR_LOOP(nv) d->Vc[nv][k][j][i] = vin[nv];
123  d->Vc[BX2][k][j][i] = j0[i]*e_t;
124  d->Vc[BX3][k][j][i] = j1[i]*e_t;
125  }
126  }else if (box->vpos == X2FACE){
127  #ifdef STAGGERED_MHD
128  BOX_LOOP(box,k,j,i) d->Vs[BX2s][k][j][i] = j0[i]*e_t;
129  #endif
130  }
131  }
132 
133  if (side == X1_END){ /* -- X1_END boundary -- */
134  if (box->vpos == CENTER){
135  BOX_LOOP(box,k,j,i){
136  VAR_LOOP(nv) d->Vc[nv][k][j][i] = vin[nv];
137  d->Vc[BX2][k][j][i] = j0[i]*e_t;
138  d->Vc[BX3][k][j][i] = j1[i]*e_t;
139  }
140  }else if (box->vpos == X2FACE){
141  #ifdef STAGGERED_MHD
142  BOX_LOOP(box,k,j,i) d->Vs[BX2s][k][j][i] = j0[i]*e_t;
143  #endif
144  }
145  }
146 }
X1_BEG
#define X1_BEG
Boundary region at X1 beg.
Definition: pluto.h:146
VX2
#define VX2
Definition: mod_defs.h:29
DENSITY
#define DENSITY
Definition: definitions_01.h:26
CENTER
#define CENTER
Definition: pluto.h:200
BesselJ0
double BesselJ0(double x)
Definition: math_misc.c:13
RBOX::vpos
int vpos
Location of the variable inside the cell.
Definition: structs.h:359
RHO
#define RHO
Definition: mod_defs.h:19
BOX_LOOP
#define BOX_LOOP(B, k, j, i)
Definition: macros.h:70
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
X1_END
#define X1_END
Boundary region at X1 end.
Definition: pluto.h:147
VX1
#define VX1
Definition: mod_defs.h:28
VAR_LOOP
#define VAR_LOOP(n)
Definition: macros.h:226
IDIR
#define IDIR
Definition: pluto.h:193
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
j
int j
Definition: analysis.c:2
if
if(divB==NULL)
Definition: analysis.c:10
k
int k
Definition: analysis.c:2
GRID::x
double * x
Definition: structs.h:80
BesselJ1
double BesselJ1(double x)
Definition: math_misc.c:45
BX3
#define BX3
Definition: mod_defs.h:27
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
ITOT_LOOP
#define ITOT_LOOP(i)
Definition: macros.h:38
i
int i
Definition: analysis.c:2
g_time
double g_time
The current integration time.
Definition: globals.h:117
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
BX2s
#define BX2s
Definition: ct.h:28
BX2
#define BX2
Definition: mod_defs.h:26
X2FACE
#define X2FACE
Definition: pluto.h:202
NX1_TOT
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55

Here is the call graph for this function: