PLUTO
res_functions.c
Go to the documentation of this file.
1 /* ///////////////////////////////////////////////////////////////////// */
2 /*!
3  \file
4  \brief Compute the curl of magnetic field.
5 
6  Compute the electric current (defined as J = curl(B))
7  for the induction and the total energy equations.
8 
9  For constrained transport MHD, J has the same staggered location of
10  the electric field and the three components (Jx, Jy, Jz) are placed
11  at different locations inside the cell:
12 
13  - Jx at (i, j+1/2, k+1/2)
14  - Jy at (i+1/2, j, k+1/2)
15  - Jz at (i+1/2, j+1/2, k)
16 
17  The same rule apply to the components of resistivity \c eta which
18  are computed and stored inside this function.
19 
20  For cell-centered MHD, the three components of J are computed
21  during each sweep direction at cell interfaces, that is,
22 
23  - IDIR: (Jx, Jy, Jz) at (i+1/2, j, k)
24  - JDIR: (Jx, Jy, Jz) at (i, j+1/2, k)
25  - KDIR: (Jx, Jy, Jz) at (i, j, k+1/2)
26 
27  Although only the two transverse components of J are actually needed
28  during the update step, we compute also the normal component
29  since the resistivity coefficients eta may depend on the total current.
30 
31  For a compact implementation, we note that the curl of a vector in
32  the three system of coordinates normally adopted may be written as
33  \f[
34  \begin{array}{lcl}
35  \left(\nabla\times\vec{B}\right)_{x_1} &=& \DS
36  \frac{1}{d_{23}}\pd{}{x_2}(a_{23}B_{x_3})
37  - \frac{1}{d_{32}}\pd{B_{x_2}}{x_3}
38  \\ \noalign{\bigskip}
39  \left(\nabla\times\vec{B}\right)_{x_2} &=& \DS
40  \frac{1}{d_{31}}\pd{B_{x_1}}{x_3}
41  - \frac{1}{d_{13}}\pd{}{x_1}(a_{13}B_{x_3})
42  \\ \noalign{\bigskip}
43  \left(\nabla\times\vec{B}\right)_{x_3} &=& \DS
44  \frac{1}{d_{12}}\pd{}{x_1}(a_{12}B_{x_2})
45  - \frac{1}{d_{21}}\pd{B_{x_1}}{x_2}
46  \end{array}
47  \f]
48  where the coefficients \f$d_{nm}=1\f$ and \f$a_{nm}=1\f$ except:
49  \f[
50  \begin{array}{lll}
51  d_{nm} = 1; &
52  \quad a_{nm} = 1; &
53  \qquad ({\rm Cartesian})
54  \\ \noalign{\medskip}
55  d_{13} = r\,,\qquad &
56  d_{32} = d_{31} = \infty;
57  \quad a_{13} = r; &
58  \qquad ({\rm Cylindrical})
59  \\ \noalign{\medskip}
60  d_{23} = d_{12} = d_{21} = r; &
61  \quad a_{12} = r; &
62  \qquad ({\rm Polar})
63  \\ \noalign{\medskip}
64  d_{23} = d_{32} = d_{31} = r\sin\theta \,,\qquad
65  d_{13} = d_{12} = d_{21} = r; &
66  \quad a_{23} = \sin\theta\,,\qquad a_{13} = a_{12} = r; &
67  \qquad ({\rm Spherical})
68  \end{array}
69  \f]
70  In the actual implementation we use \f$d_{nm} \to 1/(d_{nm}\Delta x_n)\f$.
71 
72  \note In cylindrical (or polar) geometry, there's a potential
73  singular behavior with the definition of d13 (d12) at r=0.
74  An equivalent formulation can be obtained by using absolute
75  values in the geometrical coefficient and thus replacing
76  d13 = r[i+1/2]*dr with (r[i+1] |r[i+1]| - r[i] |r[i]|)/2 and
77  a13 = r[i] with abs(r[i]).
78  This is obviously the desired derivative away from the axis
79  (where r = |r|) and it gives the correct behavior at the
80  axis, where Bphi --> -Bphi.
81  Bphi = 0 --> Bphi \sim alpha*r --> d(r Bphi)/(r dr) = 2*alpha
82 
83 
84  \attention Do NOT use this function to compute terms like (curl V)^2 !!
85 
86  \authors A. Mignone (mignone@ph.unito.it)\n
87  T. Matsakos\n
88 
89  \date March 10, 2014
90 */
91 /* ///////////////////////////////////////////////////////////////////// */
92 #include "pluto.h"
93 
94 static double ***eta[3];
95 
96 /* ********************************************************************* */
97 void GetCurrent (const Data *d, int dir, Grid *grid)
98 /*!
99  * Compute the curl of magnetic field for constrained transport MHD or
100  * cell-centered MHD.
101  *
102  * \param [in,out] d pointer to the PLUTO data structure
103  * \param [in] dir sweep direction (useless for constrained
104  * transport MHD)
105  * \param [in] grid pointer to an array of Grid structures
106  *
107  *********************************************************************** */
108 {
109  int i, j, k;
110  static int first_call = 1;
111  double *dx1, *dx2, *dx3;
112  double *r, *rp, *th, *thp;
113  double ***Bx1, ***Bx2, ***Bx3;
114  double ***Jx1, ***Jx2, ***Jx3;
115  double dx1_Bx2 = 0.0, dx2_Bx1 = 0.0;
116  double dx2_Bx3 = 0.0, dx3_Bx2 = 0.0;
117  double dx1_Bx3 = 0.0, dx3_Bx1 = 0.0;
118  double d12, d21, d13, d31, d23, d32;
119  double h3;
120  static double ***a23Bx3, ***a13Bx3, ***a12Bx2;
121  RBox box;
122 
123 /* ------------------------------------------------------------------
124  1. Set pointer shortcuts.
125  The primary magnetic field will be the staggered one
126  (if CONSTRAINTED_TRANSPORT is used) or the cell-centered field
127  otherwise.
128 
129  Note: in 2+1/2 dimensions, we also need Jx1 and Jx2 which
130  contribute to the total energy equation but not
131  to the induction equation.
132  For this reason, we set Bx3 to be equal to the \b
133  cell centered value.
134  ---------------------------------------------------------------- */
135 
136  #ifdef STAGGERED_MHD
137  D_EXPAND(Bx1 = d->Vs[BX1s]; ,
138  Bx2 = d->Vs[BX2s]; ,
139  Bx3 = d->Vs[BX3s];)
140  #if DIMENSIONS == 2 && COMPONENTS == 3
141  Bx3 = d->Vc[BX3];
142  #endif
143  #else
144  EXPAND(Bx1 = d->Vc[BX1]; ,
145  Bx2 = d->Vc[BX2]; ,
146  Bx3 = d->Vc[BX3];)
147  #endif
148 
149  Jx1 = d->J[IDIR];
150  Jx2 = d->J[JDIR];
151  Jx3 = d->J[KDIR];
152 
153  dx1 = grid[IDIR].dx;
154  dx2 = grid[JDIR].dx;
155  dx3 = grid[KDIR].dx;
156 
157 /* ----------------------------------------------------------------
158  2. Allocate static memory areas for the three arrays
159  a12Bx2, a23Bx3, a13Bx3 if necessary. Otherwise, they will
160  just be shortcuts to the default magnetic field arrays.
161  ---------------------------------------------------------------- */
162 
163  if (first_call){
164  a12Bx2 = Bx2;
165  a23Bx3 = Bx3;
166  a13Bx3 = Bx3;
167  #if GEOMETRY == CYLINDRICAL && COMPONENTS == 3
168  a13Bx3 = ARRAY_3D(NX3_MAX, NX2_MAX, NX1_MAX, double);
169  #elif GEOMETRY == POLAR && COMPONENTS >= 2
170  a12Bx2 = ARRAY_3D(NX3_MAX, NX2_MAX, NX1_MAX, double);
171  #elif GEOMETRY == SPHERICAL
172  EXPAND( ; ,
173  a12Bx2 = ARRAY_3D(NX3_MAX, NX2_MAX, NX1_MAX, double); ,
174  a13Bx3 = ARRAY_3D(NX3_MAX, NX2_MAX, NX1_MAX, double);
175  a23Bx3 = ARRAY_3D(NX3_MAX, NX2_MAX, NX1_MAX, double);)
176  #endif
177  first_call = 0;
178  }
179 
180 /* --------------------------------------------------------------
181  3. Construct goemetrical coefficients and arrays
182  -------------------------------------------------------------- */
183 
184  #if GEOMETRY == CYLINDRICAL
185  r = grid[IDIR].x; rp = grid[IDIR].xr;
186  #if COMPONENTS == 3
187  TOT_LOOP(k,j,i) a13Bx3[k][j][i] = Bx3[k][j][i]*fabs(r[i]);
188  #endif
189  #elif GEOMETRY == POLAR
190  r = grid[IDIR].x; rp = grid[IDIR].xr;
191  #if COMPONENTS >= 2
192  TOT_LOOP(k,j,i) a12Bx2[k][j][i] = Bx2[k][j][i]*r[i];
193  #endif
194  #elif GEOMETRY == SPHERICAL
195  r = grid[IDIR].x; rp = grid[IDIR].xr;
196  th = grid[JDIR].x; thp = grid[JDIR].xr;
197  TOT_LOOP(k,j,i) {
198  EXPAND( ; ,
199  a12Bx2[k][j][i] = Bx2[k][j][i]*r[i]; ,
200  a13Bx3[k][j][i] = Bx3[k][j][i]*r[i];
201  a23Bx3[k][j][i] = Bx3[k][j][i]*sin(th[j]); )
202  }
203  #endif
204 
205 /* -------------------------------------------------------------
206  4a. For staggered MHD, we compute the three components of
207  currents at cell edges.
208  In a staggered formulation the different components of J
209  have different spatial location (same as electric field)
210  and one needs
211 
212  Jx at (i , j+1/2, k+1/2)
213  Jy at (i+1/2, j , k+1/2)
214  Jz at (i+1/2, j+1/2, k )
215  -------------------------------------------------------------- */
216 
217 #ifdef STAGGERED_MHD
218 
219  for (k = 0; k < NX3_TOT-KOFFSET; k++){
220  for (j = 0; j < NX2_TOT-JOFFSET; j++){
221  for (i = 0; i < NX1_TOT-IOFFSET; i++){
222 
223  D_EXPAND(d12 = d13 = 1.0/dx1[i]; ,
224  d21 = d23 = 1.0/dx2[j]; ,
225  d31 = d32 = 1.0/dx3[k];)
226 
227  #if GEOMETRY == CYLINDRICAL
228  d32 = d31 = 0.0;
229  d13 = 2.0/(fabs(r[i+1])*r[i+1] - fabs(r[i])*r[i]); /* = 1.0/(rp*dr) */
230  #elif GEOMETRY == POLAR
231  d23 /= r[i];
232  d21 /= rp[i];
233  d12 = 2.0/(fabs(r[i+1])*r[i+1] - fabs(r[i])*r[i]); /* = 1.0/(rp*dr) */
234  #elif GEOMETRY == SPHERICAL
235  D_EXPAND(d12 /= rp[i]; d13 /= rp[i]; ,
236  d21 /= rp[i]; d23 /= r[i]*sin(thp[j]); ,
237  d32 /= r[i]*sin(thp[j]); d31 /= rp[i]*sin(th[j]);)
238  #endif
239 
240  #if COMPONENTS == 3
241  Jx1[k][j][i] = FDIFF_X2(a23Bx3,k,j,i)*d23 - FDIFF_X3( Bx2,k,j,i)*d32;
242  Jx2[k][j][i] = FDIFF_X3( Bx1,k,j,i)*d31 - FDIFF_X1(a13Bx3,k,j,i)*d13;
243  #endif
244  Jx3[k][j][i] = FDIFF_X1(a12Bx2,k,j,i)*d12 - FDIFF_X2(Bx1,k,j,i)*d21;
245  }}}
246  ComputeStaggeredEta(d, grid);
247 
248 #else
249 
250 /* ------------------------------------------------------------- */
251 /*!
252  4b. For cell-centered MHD, we compute the three components of
253  currents at a cell interface.
254  -------------------------------------------------------------- */
255 
256  if (dir == IDIR){
257 
258  /* ----------------------------------------------------
259  IDIR: Compute {Jx1, Jx2, Jx3} at i+1/2,j,k faces.
260  ---------------------------------------------------- */
261 
262  box.ib = 0; box.ie = NX1_TOT-1-IOFFSET;
263  box.jb = JOFFSET; box.je = NX2_TOT-1-JOFFSET;
264  box.kb = KOFFSET; box.ke = NX3_TOT-1-KOFFSET;
265 
266  BOX_LOOP(&box,k,j,i){
267 
268  d12 = d13 = 1.0/dx1[i];
269  d21 = d23 = 1.0/dx2[j];
270  d31 = d32 = 1.0/dx3[k];
271 
272  #if GEOMETRY == CYLINDRICAL
273  d32 = d31 = 0.0;
274  d13 = 2.0/(fabs(r[i+1])*r[i+1] - fabs(r[i])*r[i]);
275  #elif GEOMETRY == POLAR
276  d23 = d21 = 1.0/(rp[i]*dx2[j]);
277  d12 = 2.0/(fabs(r[i+1])*r[i+1] - fabs(r[i])*r[i]);
278  #elif GEOMETRY == SPHERICAL
279  h3 = rp[i]*sin(th[j]);
280 
281  D_EXPAND(d12 = 1.0/(rp[i]*dx1[i]); d13 = 1.0/(rp[i]*dx1[i]); ,
282  d21 = 1.0/(rp[i]*dx2[j]); d23 = 1.0/(h3*dx2[j]); ,
283  d32 = 1.0/(h3*dx3[k]); d31 = 1.0/(h3*dx3[k]);)
284  #endif
285 
286  #if COMPONENTS == 3
287  dx2_Bx3 = 0.5*(CDIFF_X2(a23Bx3,k,j,i) + CDIFF_X2(a23Bx3,k,j,i+1))*d23;
288  dx3_Bx2 = 0.5*(CDIFF_X3(Bx2,k,j,i) + CDIFF_X3(Bx2,k,j,i+1) )*d32;
289  Jx1[k][j][i] = (dx2_Bx3 - dx3_Bx2);
290 
291  dx3_Bx1 = 0.5*(CDIFF_X3(Bx1,k,j,i) + CDIFF_X3(Bx1,k,j,i+1))*d31;
292  dx1_Bx3 = FDIFF_X1(a13Bx3,k,j,i)*d13;
293  Jx2[k][j][i] = (dx3_Bx1 - dx1_Bx3);
294  #endif
295 
296  dx1_Bx2 = FDIFF_X1(a12Bx2,k,j,i)*d12;
297  dx2_Bx1 = 0.5*(CDIFF_X2(Bx1,k,j,i) + CDIFF_X2(Bx1,k,j,i+1))*d21;
298  Jx3[k][j][i] = (dx1_Bx2 - dx2_Bx1);
299  }
300 
301  }else if (dir == JDIR){
302 
303  /* ----------------------------------------------------
304  JDIR: Compute {Jx1, Jx2, Jx3} at i,j+1/2,k faces.
305  ---------------------------------------------------- */
306 
307  box.ib = IOFFSET; box.ie = NX1_TOT-1-IOFFSET;
308  box.jb = 0; box.je = NX2_TOT-1-JOFFSET;
309  box.kb = KOFFSET; box.ke = NX3_TOT-1-KOFFSET;
310 
311  BOX_LOOP(&box,k,j,i){
312  d12 = d13 = 1.0/dx1[i];
313  d21 = d23 = 1.0/dx2[j];
314  d31 = d32 = 1.0/dx3[k];
315 
316  #if GEOMETRY == CYLINDRICAL
317  d32 = d31 = 0.0; /* axisymmetry */
318  d13 = 1.0/(r[i]*dx1[i]);
319  #elif GEOMETRY == POLAR
320  d23 = d21 = 1.0/(r[i]*dx2[j]);
321  d12 = 1.0/(r[i]*dx1[i]);
322  #elif GEOMETRY == SPHERICAL
323  h3 = r[i]*sin(thp[j]);
324 
325 /* !!! check singularity !!! */
326 
327  D_EXPAND(d12 = 1.0/(r[i]*dx1[i]); d13 = 1.0/(r[i]*dx1[i]); ,
328  d21 = 1.0/(r[i]*dx2[j]); d23 = 1.0/(h3*dx2[j]); ,
329  d32 = 1.0/(h3*dx3[k]); d31 = 1.0/(h3*dx3[k]);)
330  #endif
331 
332  #if COMPONENTS == 3
333  dx2_Bx3 = FDIFF_X2(a23Bx3,k,j,i)*d23;
334  dx3_Bx2 = 0.5*(CDIFF_X3(Bx2,k,j,i) + CDIFF_X3(Bx2,k,j+1,i))*d32;
335  Jx1[k][j][i] = (dx2_Bx3 - dx3_Bx2);
336 
337  dx3_Bx1 = 0.5*(CDIFF_X3(Bx1,k,j,i) + CDIFF_X3(Bx1,k,j+1,i))*d31;
338  dx1_Bx3 = 0.5*(CDIFF_X1(a13Bx3,k,j,i) + CDIFF_X1(a13Bx3,k,j+1,i))*d13;
339  Jx2[k][j][i] = (dx3_Bx1 - dx1_Bx3);
340  #endif
341  dx1_Bx2 = 0.5*(CDIFF_X1(a12Bx2,k,j,i) + CDIFF_X1(a12Bx2,k,j+1,i))*d12;
342  dx2_Bx1 = FDIFF_X2(Bx1,k,j,i)*d21;
343  Jx3[k][j][i] = (dx1_Bx2 - dx2_Bx1);
344  }
345 
346  }else if (dir == KDIR) {
347 
348  /* ----------------------------------------------------
349  KDIR: Compute {Jx1, Jx2, Jx3} at i,j,k+1/2 faces.
350  ---------------------------------------------------- */
351 
352  box.ib = IOFFSET; box.ie = NX1_TOT-1-IOFFSET;
353  box.jb = JOFFSET; box.je = NX2_TOT-1-JOFFSET;
354  box.kb = 0; box.ke = NX3_TOT-1-KOFFSET;
355 
356  BOX_LOOP(&box,k,j,i){
357  d12 = d13 = 1.0/dx1[i];
358  d21 = d23 = 1.0/dx2[j];
359  d31 = d32 = 1.0/dx3[k];
360 
361  #if GEOMETRY == POLAR
362  d23 = d21 = 1.0/(rp[i]*dx2[j]);
363  d12 = 1.0/(rp[i]*dx1[i]);
364  #elif GEOMETRY == SPHERICAL
365  h3 = r[i]*sin(th[j]);
366 
367  D_EXPAND(d12 = 1.0/(r[i]*dx1[i]); d13 = 1.0/(r[i]*dx1[i]); ,
368  d21 = 1.0/(r[i]*dx2[j]); d23 = 1.0/(h3*dx2[j]); ,
369  d32 = 1.0/(h3*dx3[k]); d31 = 1.0/(h3*dx3[k]);)
370  #endif
371 
372  #if COMPONENTS == 3
373  dx2_Bx3 = 0.5*(CDIFF_X2(a23Bx3,k,j,i) + CDIFF_X2(a23Bx3,k+1,j,i))*d23;
374  dx3_Bx2 = FDIFF_X3(Bx2,k,j,i)*d32;
375  Jx1[k][j][i] = (dx2_Bx3 - dx3_Bx2);
376 
377  dx3_Bx1 = FDIFF_X3(Bx1,k,j,i)*d31;
378  dx1_Bx3 = 0.5*(CDIFF_X1(a13Bx3,k,j,i) + CDIFF_X1(a13Bx3,k+1,j,i))*d13;
379  Jx2[k][j][i] = (dx3_Bx1 - dx1_Bx3);
380  #endif
381  dx1_Bx2 = 0.5*(CDIFF_X1(a12Bx2,k,j,i) + CDIFF_X1(a12Bx2,k+1,j,i))*d12;
382  dx2_Bx1 = 0.5*(CDIFF_X2(Bx1,k,j,i) + CDIFF_X2(Bx1,k+1,j,i))*d21;
383  Jx3[k][j][i] = (dx1_Bx2 - dx2_Bx1);
384  }
385  }
386 
387 #endif /* STAGGERED_MHD */
388 }
389 
390 #ifdef STAGGERED_MHD
391 /* ********************************************************************* */
392 void ComputeStaggeredEta(const Data *d, Grid *grid)
393 /*!
394  * Compute eta coefficients using the same staggering of the current:
395  *
396  * - eta_x1 at (i, j+1/2, k+1/2)
397  * - eta_x2 at (i+1/2, j, k+1/2)
398  * - eta_x3 at (i+1/2, j+1/2, k)
399  *
400  * Since each eta_x may depend on the total current, we also need
401  * to interpolate different current components at the same place.
402  *
403  *********************************************************************** */
404 {
405  int i, j, k, nv;
406  double Js[3], etas[3], vs[NVAR];
407  double ***Jx1, ***Jx2, ***Jx3;
408  double *x1, *x2, *x3;
409  double *x1r, *x2r, *x3r;
410 
411  if (eta[0] == NULL){
412  eta[IDIR] = ARRAY_3D(NX3_TOT,NX2_TOT,NX1_TOT,double);
413  eta[JDIR] = ARRAY_3D(NX3_TOT,NX2_TOT,NX1_TOT,double);
414  eta[KDIR] = ARRAY_3D(NX3_TOT,NX2_TOT,NX1_TOT,double);
415  }
416 
417  Jx1 = d->J[IDIR]; x1 = grid[IDIR].x; x1r = grid[IDIR].xr;
418  Jx2 = d->J[JDIR]; x2 = grid[JDIR].x; x2r = grid[JDIR].xr;
419  Jx3 = d->J[KDIR]; x3 = grid[KDIR].x; x3r = grid[KDIR].xr;
420 
421  for (k = 0; k < NX3_TOT-KOFFSET; k++){
422  for (j = 0; j < NX2_TOT-JOFFSET; j++){
423  for (i = 0; i < NX1_TOT-IOFFSET; i++){
424 
425  #if COMPONENTS == 2
426 
427  /* ---- Compute J and eta_x3 at i+1/2, j+1/2, k ---- */
428 
429  Js[IDIR] = 0.0;
430  Js[JDIR] = 0.0;
431  Js[KDIR] = Jx3[k][j][i];
432  VAR_LOOP(nv) vs[nv] = AVERAGE_XY(d->Vc[nv],k,j,i);
433  Resistive_eta (vs, x1r[i], x2r[j], x3[k], Js, etas);
434  eta[KDIR][k][j][i] = etas[KDIR];
435 
436  #elif COMPONENTS == 3
437 
438  /* ---- Compute J and eta_x1 at i, j+1/2, k+1/2 ---- */
439 
440  if (i > 0) {
441  Js[IDIR] = Jx1[k][j][i];
442  Js[JDIR] = AVERAGE_XY(Jx2,k,j,i-1);
443  Js[KDIR] = AVERAGE_XZ(Jx3,k,j,i-1);
444  VAR_LOOP(nv) vs[nv] = AVERAGE_YZ(d->Vc[nv],k,j,i);
445  Resistive_eta (vs, x1[i], x2r[j], x3r[k], Js, etas);
446  eta[IDIR][k][j][i] = etas[IDIR];
447  }
448 
449  /* ---- Compute J and eta_x2 at i+1/2, j, k+1/2 ---- */
450 
451  if (j > 0) {
452  Js[IDIR] = AVERAGE_XY(Jx1,k,j-1,i);
453  Js[JDIR] = Jx2[k][j][i];
454  Js[KDIR] = AVERAGE_YZ(Jx3,k,j-1,i);
455  VAR_LOOP(nv) vs[nv] = AVERAGE_XZ(d->Vc[nv],k,j,i);
456  Resistive_eta (vs, x1r[i], x2[j], x3r[k], Js, etas);
457  eta[JDIR][k][j][i] = etas[JDIR];
458  }
459 
460  /* ---- Compute J and eta_x3 at i+1/2, j+1/2, k ---- */
461 
462  if (k > 0 || KOFFSET == 0){
463  Js[IDIR] = AVERAGE_XZ(Jx1,k-1,j,i);
464  Js[JDIR] = AVERAGE_YZ(Jx2,k-1,j,i);
465  Js[KDIR] = Jx3[k][j][i];
466  VAR_LOOP(nv) vs[nv] = AVERAGE_XY(d->Vc[nv],k,j,i);
467  Resistive_eta (vs, x1r[i], x2r[j], x3[k], Js, etas);
468  eta[KDIR][k][j][i] = etas[KDIR];
469  }
470 
471  #endif
472 
473  }}}
474 
475 }
476 
477 /* ********************************************************************* */
478 Data_Arr GetStaggeredEta()
479 /*
480  *********************************************************************** */
481 {
482  return eta;
483 }
484 #endif
DATA::Vs
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
Definition: structs.h:43
FDIFF_X2
#define FDIFF_X2(b, k, j, i)
Definition: macros.h:237
RBOX::jb
int jb
Lower corner index in the x2 direction.
Definition: structs.h:349
Data_Arr
double **** Data_Arr
Definition: pluto.h:492
DATA::J
double **** J
Electric current defined as curl(B).
Definition: structs.h:54
GetCurrent
void GetCurrent(const Data *d, int dir, Grid *grid)
Definition: res_functions.c:97
FDIFF_X1
#define FDIFF_X1(b, k, j, i)
Definition: macros.h:236
AVERAGE_YZ
#define AVERAGE_YZ(q, k, j, i)
Definition: macros.h:264
CDIFF_X2
#define CDIFF_X2(b, k, j, i)
Definition: macros.h:241
GRID::xr
double * xr
Definition: structs.h:81
eta
static double *** eta[3]
Definition: res_functions.c:94
BOX_LOOP
#define BOX_LOOP(B, k, j, i)
Definition: macros.h:70
NX2_TOT
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
Definition: globals.h:57
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
RBOX::kb
int kb
Lower corner index in the x3 direction.
Definition: structs.h:351
BX3s
#define BX3s
Definition: ct.h:29
ARRAY_3D
#define ARRAY_3D(nx, ny, nz, type)
Definition: prototypes.h:172
GRID::dx
double * dx
Definition: structs.h:83
KDIR
#define KDIR
Definition: pluto.h:195
RBOX::ib
int ib
Lower corner index in the x1 direction.
Definition: structs.h:347
NX3_MAX
#define NX3_MAX
Definition: pluto.h:743
STAGGERED_MHD
#define STAGGERED_MHD
Definition: ct.h:15
FDIFF_X3
#define FDIFF_X3(b, k, j, i)
Definition: macros.h:238
VAR_LOOP
#define VAR_LOOP(n)
Definition: macros.h:226
NX2_MAX
#define NX2_MAX
Definition: pluto.h:742
TOT_LOOP
#define TOT_LOOP(k, j, i)
Definition: macros.h:44
IDIR
#define IDIR
Definition: pluto.h:193
GRID
Definition: structs.h:78
NX3_TOT
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
Definition: globals.h:59
j
int j
Definition: analysis.c:2
AVERAGE_XY
#define AVERAGE_XY(q, k, j, i)
Definition: macros.h:260
k
int k
Definition: analysis.c:2
BX1s
#define BX1s
Definition: ct.h:27
GRID::x
double * x
Definition: structs.h:80
GEOMETRY
#define GEOMETRY
Definition: definitions_01.h:4
Resistive_eta
void Resistive_eta(double *v, double x1, double x2, double x3, double *J, double *eta)
Definition: res_eta.c:17
CYLINDRICAL
#define CYLINDRICAL
Definition: pluto.h:34
BX3
#define BX3
Definition: mod_defs.h:27
COMPONENTS
#define COMPONENTS
Definition: definitions_01.h:3
CDIFF_X1
#define CDIFF_X1(b, k, j, i)
Definition: macros.h:240
pluto.h
PLUTO main header file.
D_EXPAND
D_EXPAND(tot/[n]=(double) grid[IDIR].np_int_glob;, tot/[n]=(double) grid[JDIR].np_int_glob;, tot/[n]=(double) grid[KDIR].np_int_glob;)
Definition: analysis.c:27
DATA
Definition: structs.h:30
i
int i
Definition: analysis.c:2
RBOX::ie
int ie
Upper corner index in the x1 direction.
Definition: structs.h:348
RBOX::ke
int ke
Upper corner index in the x3 direction.
Definition: structs.h:352
RBOX::je
int je
Upper corner index in the x2 direction.
Definition: structs.h:350
BX1
#define BX1
Definition: mod_defs.h:25
AVERAGE_XZ
#define AVERAGE_XZ(q, k, j, i)
Definition: macros.h:262
BX2s
#define BX2s
Definition: ct.h:28
BX2
#define BX2
Definition: mod_defs.h:26
CDIFF_X3
#define CDIFF_X3(b, k, j, i)
Definition: macros.h:242
JDIR
#define JDIR
Definition: pluto.h:194
NVAR
#define NVAR
Definition: pluto.h:609
RBOX
Definition: structs.h:346
NX1_MAX
#define NX1_MAX
Definition: pluto.h:741
NX1_TOT
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55