PLUTO
Macros | Functions
ct_field_average.c File Reference

Performs various magnetic field averaging operations. More...

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

Go to the source code of this file.

Macros

#define IEXT_LOOP(i)   for (i = IOFFSET; i < NX1_TOT - IOFFSET; i++)
 
#define JEXT_LOOP(j)   for (j = JOFFSET; j < NX2_TOT - JOFFSET; j++)
 
#define KEXT_LOOP(k)   for (k = KOFFSET; k < NX3_TOT - KOFFSET; k++)
 

Functions

void CT_AverageMagneticField (double ****bf, double ****UU, Grid *grid)
 
void CT_AverageNormalMagField (const Data *d, int side, Grid *grid)
 Compute the normal component of the volume-average magnetic field from the staggered components in the ghost zones. More...
 
void CT_AverageTransverseMagField (const Data *d, int side, Grid *grid)
 Compute the transverse component of the volume-average magnetic field from the staggered components in the ghost zones. More...
 

Detailed Description

Performs various magnetic field averaging operations.

Author
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Aug 16, 2012

Definition in file ct_field_average.c.

Macro Definition Documentation

#define IEXT_LOOP (   i)    for (i = IOFFSET; i < NX1_TOT - IOFFSET; i++)
#define JEXT_LOOP (   j)    for (j = JOFFSET; j < NX2_TOT - JOFFSET; j++)
#define KEXT_LOOP (   k)    for (k = KOFFSET; k < NX3_TOT - KOFFSET; k++)

Function Documentation

void CT_AverageMagneticField ( double ****  bf,
double ****  UU,
Grid *  grid 
)

Average staggered magnetic field components to form a zone-centered field, e.g., $ <B_i> = (b_{i+1/2} + b_{i-1/2})/2$

The averaging procedure is performed inside the computational domain augmented by an additional row/plane of ghost zones in the boundary regions. Inclusion of boundary regions is useful only during the predictor step of the CT-CTU algorithm, whereas is useless for RK time stepping.

When the CT_EN_CORRECTION flag is enabled, we also redefine the zone total energy using the newly formed cell-centered field, i.e.,

\[ E_i \to E_i - \frac{\mathbf{B}_i^2}{2} + \frac{<\mathbf{B}_i^2>}{2} \]

Parameters
[in]bfarray of staggered fields
[out]UUarray of conservative variables
[in]gridpointer to Grid structure

Definition at line 13 of file ct_field_average.c.

35 {
36  int i, j, k;
37  double b2_old, b2_new, bx_ave, by_ave, bz_ave;
38  double rp, rm;
39  double ***bx, ***by, ***bz;
40  double *dx, *dy, *A1, *A2, *dV1, *dV2;
41  double *r;
42 
43  D_EXPAND(bx = bf[BX1s]; ,
44  by = bf[BX2s]; ,
45  bz = bf[BX3s]; )
46 
47  dx = grid[IDIR].dx;
48  dy = grid[JDIR].dx;
49  A1 = grid[IDIR].A;
50  A2 = grid[JDIR].A;
51  dV1 = grid[IDIR].dV;
52  dV2 = grid[JDIR].dV;
53  r = grid[IDIR].x;
54 
55 /* ---------------------------------------------------------
56  Loop over all zones in the domain.
57  We include also one set of boundary zones which
58  is useful only during the predictor step of the CT-CTU
59  algorithm, whereas is useless for RK time stepping.
60  --------------------------------------------------------- */
61 
62  for (k = KBEG - KOFFSET; k <= KEND + KOFFSET; k++){
63  for (j = JBEG - JOFFSET; j <= JEND + JOFFSET; j++){
64  for (i = IBEG - IOFFSET; i <= IEND + IOFFSET; i++){
65 
66  #if GEOMETRY == CARTESIAN
67 
68  D_EXPAND( bx_ave = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
69  by_ave = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
70  bz_ave = 0.5*(bz[k][j][i] + bz[k - 1][j][i]); )
71 
72  #elif GEOMETRY == CYLINDRICAL
73 /*
74  bx_ave = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/(A1[i] + A1[i-1]);
75  by_ave = 0.5*( by[k][j][i] + by[k][j - 1][i]);
76 */
77 
78  bx_ave = 0.5*(bx[k][j][i] + bx[k][j][i - 1]);
79  by_ave = 0.5*(by[k][j][i] + by[k][j - 1][i]);
80 
81  #elif GEOMETRY == POLAR
82 
83  rp = grid[IDIR].xr[i];
84  rm = grid[IDIR].xl[i];
85 
86  D_EXPAND(bx_ave = (rp*bx[k][j][i] + rm*bx[k][j][i - 1])/(rp + rm); ,
87  by_ave = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
88  bz_ave = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);)
89 
90  #elif GEOMETRY == SPHERICAL
91 
92  D_EXPAND(bx_ave = 0.5*(A1[i]*bx[k][j][i] + A1[i - 1]*bx[k][j][i - 1]);
93  bx_ave *= dx[i]/dV1[i]; ,
94 
95  by_ave = 0.5*(A2[j]*by[k][j][i] + A2[j - 1]*by[k][j - 1][i]);
96  by_ave *= dy[j]/dV2[j]; ,
97 
98  bz_ave = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);)
99  #endif
100 
101  /* ---------------------------------------------
102  apply energy correction if necessary
103  --------------------------------------------- */
104 
105  #if CT_EN_CORRECTION == YES && HAVE_ENERGY
106  b2_old = D_EXPAND( UU[k][j][i][BX1]*UU[k][j][i][BX1],
107  + UU[k][j][i][BX2]*UU[k][j][i][BX2],
108  + UU[k][j][i][BX3]*UU[k][j][i][BX3]);
109  #endif
110 
111  D_EXPAND( UU[k][j][i][BX1] = bx_ave; ,
112  UU[k][j][i][BX2] = by_ave; ,
113  UU[k][j][i][BX3] = bz_ave; )
114 
115  #if CT_EN_CORRECTION == YES && HAVE_ENERGY
116  b2_new = D_EXPAND(bx_ave*bx_ave, + by_ave*by_ave, + bz_ave*bz_ave);
117  UU[k][j][i][ENG] += 0.5*(b2_new - b2_old);
118  #endif
119  }}}
120 
121 }
GRID::xr
double * xr
Definition: structs.h:81
YES
#define YES
Definition: pluto.h:25
SPHERICAL
#define SPHERICAL
Definition: pluto.h:36
GRID::dV
double * dV
Cell volume.
Definition: structs.h:86
BX3s
#define BX3s
Definition: ct.h:29
GRID::dx
double * dx
Definition: structs.h:83
CT_EN_CORRECTION
#define CT_EN_CORRECTION
Definition: definitions_01.h:49
GRID::xl
double * xl
Definition: structs.h:82
IDIR
#define IDIR
Definition: pluto.h:193
j
int j
Definition: analysis.c:2
IEND
long int IEND
Upper grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:37
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
CYLINDRICAL
#define CYLINDRICAL
Definition: pluto.h:34
BX3
#define BX3
Definition: mod_defs.h:27
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
i
int i
Definition: analysis.c:2
BX1
#define BX1
Definition: mod_defs.h:25
KBEG
long int KBEG
Lower grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:43
BX2s
#define BX2s
Definition: ct.h:28
BX2
#define BX2
Definition: mod_defs.h:26
KEND
long int KEND
Upper grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:45
JDIR
#define JDIR
Definition: pluto.h:194
JBEG
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:39
JEND
long int JEND
Upper grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:41
IBEG
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87
HAVE_ENERGY
#define HAVE_ENERGY
Definition: pluto.h:395

Here is the call graph for this function:

Here is the caller graph for this function:

void CT_AverageNormalMagField ( const Data *  d,
int  side,
Grid *  grid 
)

Compute the normal component of the volume-average magnetic field from the staggered components in the ghost zones.

For a given "side" of the boundary, average the staggered magnetic field components normal to that boundary into a cell-centered field. For instance, at X3_BEG we only average Bz. Transverse components should be averaged using a slightly different stencil. This function is automatically called from Boundary() for consistency between staggered and zone-centered fields.

Parameters
[in,out]dpointer to PLUTO Data structure
[in]sidethe side
[in]gridpointer to PLUTO Grid structure
Todo:
replace the loops with more compact macro, such as X1_BEG_LOOP()...

Definition at line 123 of file ct_field_average.c.

144 {
145  int i, j, k;
146  double *Ap, *Am;
147  double ***Bx, ***By, ***Bz;
148  double ***bx, ***by, ***bz;
149 
150 /* ------------------------------------------------------
151  X1 boundaries
152  ------------------------------------------------------ */
153 
154  if (side == X1_BEG){
155  Ap = grid[IDIR].A; Am = Ap - 1;
156  Bx = d->Vc[BX1]; bx = d->Vs[BX1s];
157  KTOT_LOOP(k) JTOT_LOOP(j) for (i = 0; i < IBEG; i++) {
158  #if GEOMETRY == CARTESIAN
159  Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]);
160  #elif GEOMETRY == CYLINDRICAL
161  Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]);
162  #elif GEOMETRY == POLAR
163  Bx[k][j][i] = (Ap[i]*bx[k][j][i] + Am[i]*bx[k][j][i - 1])/(Ap[i] + Am[i]);
164  #elif GEOMETRY == SPHERICAL
165  Bx[k][j][i] = (Ap[i]*bx[k][j][i] + Am[i]*bx[k][j][i - 1])/(Ap[i] + Am[i]);
166  #endif
167  }
168  }else if (side == X1_END){
169  Ap = grid[IDIR].A; Am = Ap - 1;
170  Bx = d->Vc[BX1]; bx = d->Vs[BX1s];
171  KTOT_LOOP(k) JTOT_LOOP(j) for (i = IEND+1; i < NX1_TOT; i++) {
172  #if GEOMETRY == CARTESIAN
173  Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]);
174  #elif GEOMETRY == CYLINDRICAL
175  Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]);
176  #elif GEOMETRY == POLAR
177  Bx[k][j][i] = (Ap[i]*bx[k][j][i] + Am[i]*bx[k][j][i - 1])/(Ap[i] + Am[i]);
178  #elif GEOMETRY == SPHERICAL
179  Bx[k][j][i] = (Ap[i]*bx[k][j][i] + Am[i]*bx[k][j][i - 1])/(Ap[i] + Am[i]);
180  #endif
181  }
182  }
183 
184 /* ------------------------------------------------------
185  X2 boundaries
186  ------------------------------------------------------ */
187 
188  if (side == X2_BEG){
189  Ap = grid[JDIR].A; Am = Ap - 1;
190  By = d->Vc[BX2]; by = d->Vs[BX2s];
191  KTOT_LOOP(k) for (j = 0; j < JBEG; j++) ITOT_LOOP(i){
192  #if GEOMETRY == CARTESIAN
193  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]);
194  #elif GEOMETRY == CYLINDRICAL
195  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]);
196  #elif GEOMETRY == POLAR
197  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]);
198  #elif GEOMETRY == SPHERICAL
199  By[k][j][i] = (Ap[j]*by[k][j][i] + Am[j]*by[k][j - 1][i])/(Ap[j] + Am[j]);
200  #endif
201  }
202  }else if (side == X2_END){
203  Ap = grid[JDIR].A; Am = Ap - 1;
204  By = d->Vc[BX2]; by = d->Vs[BX2s];
205  KTOT_LOOP(k) for (j = JEND+1; j < NX2_TOT; j++) ITOT_LOOP(i){
206  #if GEOMETRY == CARTESIAN
207  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]);
208  #elif GEOMETRY == CYLINDRICAL
209  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]);
210  #elif GEOMETRY == POLAR
211  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]);
212  #elif GEOMETRY == SPHERICAL
213  By[k][j][i] = (Ap[j]*by[k][j][i] + Am[j]*by[k][j - 1][i])/(Ap[j] + Am[j]);
214  #endif
215  }
216  }
217 
218 /* ------------------------------------------------------
219  X3 boundaries
220  ------------------------------------------------------ */
221 
222  if (side == X3_BEG){
223  Bz = d->Vc[BX3]; bz = d->Vs[BX3s];
224  for (k = 0; k < KBEG; k++) JTOT_LOOP(j) ITOT_LOOP(i){
225  #if GEOMETRY == CARTESIAN
226  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);
227  #elif GEOMETRY == CYLINDRICAL
228  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);
229  #elif GEOMETRY == POLAR
230  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);
231  #elif GEOMETRY == SPHERICAL
232  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);
233  #endif
234  }
235  }else if (side == X3_END){
236  Bz = d->Vc[BX3]; bz = d->Vs[BX3s];
237  for (k = KEND+1; k < NX3_TOT; k++) JTOT_LOOP(j) ITOT_LOOP(i){
238  #if GEOMETRY == CARTESIAN
239  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);
240  #elif GEOMETRY == CYLINDRICAL
241  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);
242  #elif GEOMETRY == POLAR
243  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);
244  #elif GEOMETRY == SPHERICAL
245  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);
246  #endif
247  }
248  }
249 }
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
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
KTOT_LOOP
#define KTOT_LOOP(k)
Definition: macros.h:40
BX3s
#define BX3s
Definition: ct.h:29
X1_END
#define X1_END
Boundary region at X1 end.
Definition: pluto.h:147
IDIR
#define IDIR
Definition: pluto.h:193
X2_END
#define X2_END
Boundary region at X2 end.
Definition: pluto.h:149
Bx
static double Bx
Definition: hlld.c:62
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
IEND
long int IEND
Upper grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:37
k
int k
Definition: analysis.c:2
BX1s
#define BX1s
Definition: ct.h:27
BX3
#define BX3
Definition: mod_defs.h:27
X3_END
#define X3_END
Boundary region at X3 end.
Definition: pluto.h:151
ITOT_LOOP
#define ITOT_LOOP(i)
Definition: macros.h:38
JTOT_LOOP
#define JTOT_LOOP(j)
Definition: macros.h:39
i
int i
Definition: analysis.c:2
BX1
#define BX1
Definition: mod_defs.h:25
KBEG
long int KBEG
Lower grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:43
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
KEND
long int KEND
Upper grid index of the computational domain in the the X3 direction for the local processor...
Definition: globals.h:45
JDIR
#define JDIR
Definition: pluto.h:194
JBEG
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:39
JEND
long int JEND
Upper grid index of the computational domain in the the X2 direction for the local processor...
Definition: globals.h:41
NX1_TOT
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
Definition: globals.h:55
IBEG
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
Definition: globals.h:35
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87

Here is the caller graph for this function:

void CT_AverageTransverseMagField ( const Data *  d,
int  side,
Grid *  grid 
)

Compute the transverse component of the volume-average magnetic field from the staggered components in the ghost zones.

For a given "side" of the boundary, average the staggered magnetic field components transverse to that boundary into a cell-centered field. For instance, at X2_BEG we average Bx and Bz. This call is necessary only if one requires consistency between staggered and cell-centered magnetic field components.

Parameters
[in,out]dpointer to PLUTO Data structure
[in]sidethe side
[in]gridpointer to PLUTO Grid structure
Todo:
replace the loops using more compact notations, such as X1_BEG_LOOP()...

Definition at line 252 of file ct_field_average.c.

275 {
276  int i, j, k;
277  double *A1, *A2, *A3;
278  double ***Bx, ***By, ***Bz;
279  double ***bx, ***by, ***bz;
280 
281  A1 = grid[IDIR].A;
282  A2 = grid[JDIR].A;
283  A3 = grid[KDIR].A;
284 
285  D_EXPAND(Bx = d->Vc[BX1]; bx = d->Vs[BX1s]; ,
286  By = d->Vc[BX2]; by = d->Vs[BX2s]; ,
287  Bz = d->Vc[BX3]; bz = d->Vs[BX3s]; )
288 
289 /* ------------------------------------------------------
290  X1 boundaries
291  ------------------------------------------------------ */
292 
293  if (side == X1_BEG){
294  X1_BEG_LOOP(k,j,i){
295  #if GEOMETRY != SPHERICAL
296  D_EXPAND( ,
297  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
298  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]); )
299  #else
300  D_EXPAND( ,
301  By[k][j][i] = (A2[j]*by[k][j][i] + A2[j-1]*by[k][j - 1][i])/
302  (A2[j] + A2[j-1]); ,
303  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);)
304  #endif
305  }
306  }else if (side == X1_END){
307  X1_END_LOOP(k,j,i){
308  #if GEOMETRY != SPHERICAL
309  D_EXPAND( ,
310  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
311  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]); )
312  #else
313  D_EXPAND( ,
314  By[k][j][i] = (A2[j]*by[k][j][i] + A2[j-1]*by[k][j - 1][i])/
315  (A2[j] + A2[j-1]); ,
316  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);)
317  #endif
318  }
319  }
320 
321 /* ------------------------------------------------------
322  X2 boundaries
323  ------------------------------------------------------ */
324 
325  if (side == X2_BEG){
326  X2_BEG_LOOP(k,j,i){
327  #if GEOMETRY == CARTESIAN
328  D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
329  ,
330  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]); )
331  #elif GEOMETRY == CYLINDRICAL
332  D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
333  ,
334  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]); )
335  #elif GEOMETRY == POLAR || GEOMETRY == SPHERICAL
336  D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
337  (A1[i] + A1[i-1]); ,
338  ,
339  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);)
340  #endif
341  }
342  }else if (side == X2_END){
343  X2_END_LOOP(k,j,i){
344  #if GEOMETRY == CARTESIAN
345  D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
346  ,
347  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]); )
348  #elif GEOMETRY == CYLINDRICAL
349  D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
350  ,
351  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]); )
352  #elif GEOMETRY == POLAR || GEOMETRY == SPHERICAL
353  D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
354  (A1[i] + A1[i-1]); ,
355  ,
356  Bz[k][j][i] = 0.5*(bz[k][j][i] + bz[k - 1][j][i]);)
357  #endif
358  }
359  }
360 
361 /* ------------------------------------------------------
362  X3 boundaries
363  ------------------------------------------------------ */
364 
365  if (side == X3_BEG){
366  X3_BEG_LOOP(k,j,i){
367  #if GEOMETRY == CARTESIAN
368  D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
369  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
370  )
371 
372  #elif GEOMETRY == CYLINDRICAL
373  D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
374  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
375  )
376  #elif GEOMETRY == POLAR
377  D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
378  (A1[i] + A1[i-1]); ,
379  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
380  )
381  #elif GEOMETRY == SPHERICAL
382 
383  D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
384  (A1[i] + A1[i-1]); ,
385  By[k][j][i] = (A2[j]*by[k][j][i] + A2[j-1]*by[k][j - 1][i])/
386  (A2[j] + A2[j-1]); ,
387  )
388  #endif
389  }
390  }else if (side == X3_END){
391  X3_END_LOOP(k,j,i){
392  #if GEOMETRY == CARTESIAN
393  D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
394  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
395  )
396 
397  #elif GEOMETRY == CYLINDRICAL
398  D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
399  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
400  )
401  #elif GEOMETRY == POLAR
402  D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
403  (A1[i] + A1[i-1]); ,
404  By[k][j][i] = 0.5*(by[k][j][i] + by[k][j - 1][i]); ,
405  )
406  #elif GEOMETRY == SPHERICAL
407 
408  D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
409  (A1[i] + A1[i-1]); ,
410  By[k][j][i] = (A2[j]*by[k][j][i] + A2[j-1]*by[k][j - 1][i])/
411  (A2[j] + A2[j-1]); ,
412  )
413  #endif
414  }
415  }
416 }
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
DATA::Vc
double **** Vc
The main four-index data array used for cell-centered primitive variables.
Definition: structs.h:31
SPHERICAL
#define SPHERICAL
Definition: pluto.h:36
X3_END_LOOP
#define X3_END_LOOP(k, j, i)
Definition: macros.h:52
BX3s
#define BX3s
Definition: ct.h:29
X1_END
#define X1_END
Boundary region at X1 end.
Definition: pluto.h:147
X2_BEG_LOOP
#define X2_BEG_LOOP(k, j, i)
Definition: macros.h:47
KDIR
#define KDIR
Definition: pluto.h:195
POLAR
#define POLAR
Definition: pluto.h:35
X2_END_LOOP
#define X2_END_LOOP(k, j, i)
Definition: macros.h:51
IDIR
#define IDIR
Definition: pluto.h:193
X2_END
#define X2_END
Boundary region at X2 end.
Definition: pluto.h:149
Bx
static double Bx
Definition: hlld.c:62
j
int j
Definition: analysis.c:2
if
if(divB==NULL)
Definition: analysis.c:10
k
int k
Definition: analysis.c:2
BX1s
#define BX1s
Definition: ct.h:27
X1_END_LOOP
#define X1_END_LOOP(k, j, i)
Definition: macros.h:50
GEOMETRY
#define GEOMETRY
Definition: definitions_01.h:4
CYLINDRICAL
#define CYLINDRICAL
Definition: pluto.h:34
BX3
#define BX3
Definition: mod_defs.h:27
X3_END
#define X3_END
Boundary region at X3 end.
Definition: pluto.h:151
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
X1_BEG_LOOP
#define X1_BEG_LOOP(k, j, i)
Definition: macros.h:46
i
int i
Definition: analysis.c:2
BX1
#define BX1
Definition: mod_defs.h:25
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
X3_BEG_LOOP
#define X3_BEG_LOOP(k, j, i)
Definition: macros.h:48
GRID::A
double * A
Right interface area, A[i] = .
Definition: structs.h:87

Here is the call graph for this function:

Here is the caller graph for this function: