PLUTO
Functions
init.c File Reference

Magnetic field diffusion in 2D and 3D. More...

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

Go to the source code of this file.

Functions

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

Detailed Description

Magnetic field diffusion in 2D and 3D.

Sets the initial conditions for a magnetic field diffusion problem in 2 or 3 dimensions. This is a useful test to check the ability of the code to solve standard diffusion problems. The magnetic field has initially a Gaussian profile, and an anisotropic resistivity is possible. This problem has an analytical solution given, in 2D, by

\[ B_x(y,t) = \exp(-y^2/4\eta_zt)/\sqrt{t} \quad\quad\quad B_y(x,t) = \exp(-x^2/4\eta_zt)/\sqrt{t} \quad\quad\quad B_z(x,y,t) = \exp(-x^2/4\eta_yt)\exp(-y^2/4\eta_xt)/t \]

and in 3D by

\[ \begin{array}{lcl} B_x(y,z,t) &=& \exp(-y^2/4\eta_zt)\exp(-z^2/4\eta_yt)/t \\ \noalign{\medskip} B_y(x,z,t) &=& \exp(-x^2/4\eta_zt)\exp(-z^2/4\eta_xt)/t \\ \noalign{\medskip} B_z(x,y,t) &=& \exp(-x^2/4\eta_yt)\exp(-y^2/4\eta_xt)/t \end{array} \]

The initial condition is simply set using the previous profiles with t=1. In order to solve only the parabolic term in the induction equation ( $ \nabla\times\vec{J}=\nabla\times(\eta\nabla\times\vec{B})$) we give to the fluid a very large intertia using a high density value. Moreover, to avoid any fluid motion, the velocity is reset to zero at each time step by using the INTERNAL_BOUNDARY.

The runtime parameters that are read from pluto.ini are

  • g_inputParam[ETAX]: sets the resistivity along the $x$ direction;
  • g_inputParam[ETAY]: sets the resistivity along the $y$ direction;
  • g_inputParam[ETAZ]: sets the resistivity along the $z$ direction;
3D-cart-0.jpg
Initial and final profiles of the numerical (points) and analytical (lines) solutions for a component of the magnetic field.

The configurations use both EXPLICIT and STS time integrators and different geometries are explored:

Conf.GEOMETRY DIMT.STEPPING divBRESISTIVITY
#01 CARTESIAN 3 RK2 8W EXPLICIT
#02 CARTESIAN 3 HANCOCK GLMSTS
#03 CARTESIAN 3 RK3 CT EXPLICIT
#04 CARTESIAN 3 HANCOCK GLMEXPLICIT
#05 POLAR 3 RK2 8W EXPLICIT
#06 POLAR 3 RK2 8W STS
#07 SPHERICAL 3 RK2 8W EXPLICIT
#08 SPHERICAL 3 RK2 8W STS
#09 SPHERICAL 3 RK3 CT EXPLICIT
#10 CARTESIAN 2 HANCOCK CT EXPLICIT
#11 CARTESIAN 3 HANCOCK CT EXPLICIT
#12 CARTESIAN 2 HANCOCK CT STS
#13 SPHERICAL 3 RK2 CT STS
Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
T. Matsakos (titos.nosp@m.@odd.nosp@m.job.u.nosp@m.chic.nosp@m.ago.e.nosp@m.du)
Date
Sept 17, 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 84 of file init.c.

88 {}
void BoundValues ( double *  v,
double  x1,
double  x2,
double  x3,
double  t 
)

Definition at line 158 of file init.c.

164 {
165  double Bx, By, Bz;
166  double x, y, z;
167  double eta_x, eta_y, eta_z, st;
168 
169  eta_x = g_inputParam[ETAX];
170  eta_y = g_inputParam[ETAY];
171  eta_z = g_inputParam[ETAZ];
172 
173 /* -- find Cartesian coordinates from (x1,x2,x3) -- */
174 
175  #if GEOMETRY == CARTESIAN
176  x = x1; y = x2; z = x3;
177  #elif GEOMETRY == POLAR
178  x = x1*cos(x2) - 5.0;
179  y = x1*sin(x2) - 5.0;
180  z = x3 - 5.0;
181  #elif GEOMETRY == SPHERICAL
182  x = x1*cos(x3)*sin(x2) - 5.0;
183  y = x1*sin(x3)*sin(x2) - 5.0;
184  z = x1*cos(x2) - 5.0;
185  #else
186  #error geometry not valid
187  #endif
188 
189  v[RHO] = 1.0e9;
190  v[VX1] = 0.0;
191  v[VX2] = 0.0;
192  v[VX3] = 0.0;
193  #if EOS != ISOTHERMAL
194  v[PRS] = 1.0;
195  #endif
196  #ifdef GLM_MHD
197  v[PSI_GLM] = 0.0;
198  #endif
199 
200 /* -- find Cartesian components of magnetic field -- */
201 /*
202  Bx = exp(-0.25*(y*y)/g_inputParam[ETAZ]/t)*exp(-0.25*(z*z)/g_inputParam[ETAY]/t)/t;
203  By = exp(-0.25*(x*x)/g_inputParam[ETAZ]/t)*exp(-0.25*(z*z)/g_inputParam[ETAX]/t)/t;
204  Bz = exp(-0.25*(x*x)/g_inputParam[ETAY]/t)*exp(-0.25*(y*y)/g_inputParam[ETAX]/t)/t;
205 */
206 
207  #if DIMENSIONS == 2
208  st = sqrt(t);
209 
210  Bx = exp(-0.25*y*y/(eta_z*t))/st;
211  By = exp(-0.25*x*x/(eta_z*t))/st;
212  Bz = exp(-0.25*x*x/(eta_y*t))*exp(-0.25*y*y/(eta_x*t))/t;
213  #elif DIMENSIONS == 3
214  Bx = exp(-0.25*y*y/(eta_z*t))*exp(-0.25*z*z/(eta_y*t))/t;
215  By = exp(-0.25*x*x/(eta_z*t))*exp(-0.25*z*z/(eta_x*t))/t;
216  Bz = exp(-0.25*x*x/(eta_y*t))*exp(-0.25*y*y/(eta_x*t))/t;
217  #endif
218 
219 /* -----------------------------------------------------------
220  Transform back to original coordinate system
221  ----------------------------------------------------------- */
222 
223  #if GEOMETRY == CARTESIAN
224  v[BX1] = Bx;
225  v[BX2] = By;
226  v[BX3] = Bz;
227  #elif GEOMETRY == POLAR
228  v[BX1] = cos(x2)*Bx + sin(x2)*By;
229  v[BX2] = -sin(x2)*Bx + cos(x2)*By;
230  v[BX3] = Bz;
231  #elif GEOMETRY == SPHERICAL
232  v[BX1] = cos(x3)*sin(x2)*Bx + sin(x3)*sin(x2)*By + cos(x2)*Bz;
233  v[BX2] = cos(x3)*cos(x2)*Bx + sin(x3)*cos(x2)*By - sin(x2)*Bz;
234  v[BX3] = -sin(x3)*Bx + cos(x3)*By;
235  #else
236  print1 ("! BoundValues: GEOMETRY not defined\n");
237  QUIT_PLUTO(1);
238  #endif
239 
240 }
VX2
#define VX2
Definition: mod_defs.h:29
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
RHO
#define RHO
Definition: mod_defs.h:19
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
ETAY
#define ETAY
Definition: definitions_01.h:27
VX1
#define VX1
Definition: mod_defs.h:28
Bx
static double Bx
Definition: hlld.c:62
g_inputParam
double g_inputParam[32]
Array containing the user-defined parameters.
Definition: globals.h:131
GRID::x
double * x
Definition: structs.h:80
BX3
#define BX3
Definition: mod_defs.h:27
ETAX
#define ETAX
Definition: definitions_01.h:26
BX1
#define BX1
Definition: mod_defs.h:25
VX3
#define VX3
Definition: mod_defs.h:30
BX2
#define BX2
Definition: mod_defs.h:26
ETAZ
#define ETAZ
Definition: definitions_01.h:28
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125

Here is the call graph for this function:

Here is the caller graph for this function:

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

79 {
80  BoundValues(us, x1, x2, x3, 1.0);
81 }
BoundValues
void BoundValues(double *v, double x1, double x2, double x3, double t)
Definition: init.c:158

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

95 {
96  int i, j, k, nv;
97  double *x, *y, *z;
98  double *xp, *yp, *zp;
99  double t, dt, vb[256];
100 
101  if (side == 0){
102  #if GEOMETRY == CARTESIAN
103  TOT_LOOP(k,j,i){
104  EXPAND(d->Vc[VX1][k][j][i] = 0.0; ,
105  d->Vc[VX2][k][j][i] = 0.0; ,
106  d->Vc[VX3][k][j][i] = 0.0;)
107  #if ENTROPY_SWITCH == YES
108  d->flag[k][j][i] |= FLAG_ENTROPY;
109  #endif
110  }
111  #endif
112  }
113 
114  t = g_time + 1.0;
115 
116  x = grid[IDIR].x; xp = grid[IDIR].xr;
117  y = grid[JDIR].x; yp = grid[JDIR].xr;
118  z = grid[KDIR].x; zp = grid[KDIR].xr;
119 
120  if (side == X1_BEG || side == X2_BEG || side == X3_BEG ||
121  side == X1_END || side == X2_END || side == X3_END){ /* -- All Boundaries -- */
122  if (box->vpos == CENTER) {
123  BOX_LOOP(box,k,j,i){
124  BoundValues (vb, x[i], y[j], z[k], t);
125  for (nv = NVAR; nv--; ) d->Vc[nv][k][j][i] = vb[nv];
126  }
127  }else if (box->vpos == X1FACE){
128  #ifdef STAGGERED_MHD
129  if (side == X1_BEG || side == X1_END) return;
130  BOX_LOOP(box,k,j,i){
131  BoundValues (vb, xp[i], y[j], z[k], t);
132  d->Vs[BX1s][k][j][i] = vb[BX1];
133  }
134  #endif
135  }else if (box->vpos == X2FACE){
136  if (side == X2_BEG || side == X2_END) return;
137  #ifdef STAGGERED_MHD
138  BOX_LOOP(box,k,j,i){
139  BoundValues (vb, x[i], yp[j], z[k], t);
140  d->Vs[BX2s][k][j][i] = vb[BX2];
141  }
142  #endif
143 
144  }else if (box->vpos == X3FACE){
145  if (side == X3_BEG || side == X3_END) return;
146  #ifdef STAGGERED_MHD
147  BOX_LOOP(box,k,j,i) {
148  BoundValues (vb, x[i], y[j], zp[k], t);
149  d->Vs[BX3s][k][j][i] = vb[BX3];
150  }
151  #endif
152 
153  }
154  }
155 }
X3_BEG
#define X3_BEG
Boundary region at X3 beg.
Definition: pluto.h:150
X1_BEG
#define X1_BEG
Boundary region at X1 beg.
Definition: pluto.h:146
DATA::Vs
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
Definition: structs.h:43
VX2
#define VX2
Definition: mod_defs.h:29
CENTER
#define CENTER
Definition: pluto.h:200
RBOX::vpos
int vpos
Location of the variable inside the cell.
Definition: structs.h:359
BoundValues
void BoundValues(double *v, double x1, double x2, double x3, double t)
Definition: init.c:158
GRID::xr
double * xr
Definition: structs.h:81
X3FACE
#define X3FACE
Definition: pluto.h:203
BOX_LOOP
#define BOX_LOOP(B, k, j, i)
Definition: macros.h:70
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
BX3s
#define BX3s
Definition: ct.h:29
X1_END
#define X1_END
Boundary region at X1 end.
Definition: pluto.h:147
VX1
#define VX1
Definition: mod_defs.h:28
KDIR
#define KDIR
Definition: pluto.h:195
DATA::flag
unsigned char *** flag
Pointer to a 3D array setting useful integration flags that are retrieved during integration.
Definition: structs.h:55
X1FACE
#define X1FACE
Definition: pluto.h:201
TOT_LOOP
#define TOT_LOOP(k, j, i)
Definition: macros.h:44
FLAG_ENTROPY
#define FLAG_ENTROPY
Update pressure using entropy equation.
Definition: pluto.h:182
IDIR
#define IDIR
Definition: pluto.h:193
X2_END
#define X2_END
Boundary region at X2 end.
Definition: pluto.h:149
j
int j
Definition: analysis.c:2
k
int k
Definition: analysis.c:2
BX1s
#define BX1s
Definition: ct.h:27
GRID::x
double * x
Definition: structs.h:80
BX3
#define BX3
Definition: mod_defs.h:27
X3_END
#define X3_END
Boundary region at X3 end.
Definition: pluto.h:151
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
X2_BEG
#define X2_BEG
Boundary region at X2 beg.
Definition: pluto.h:148
BX2
#define BX2
Definition: mod_defs.h:26
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609
X2FACE
#define X2FACE
Definition: pluto.h:202

Here is the call graph for this function: