37 double b2_old, b2_new, bx_ave, by_ave, bz_ave;
39 double ***bx, ***by, ***bz;
40 double *dx, *dy, *A1, *A2, *dV1, *dV2;
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++){
66 #if GEOMETRY == CARTESIAN
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]); )
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]);
81 #elif GEOMETRY == POLAR
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]);)
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]; ,
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]; ,
98 bz_ave = 0.5*(bz[
k][
j][
i] + bz[k - 1][
j][
i]);)
107 + UU[k][j][i][
BX2]*UU[k][j][i][
BX2],
108 + UU[k][j][i][
BX3]*UU[k][j][i][
BX3]);
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; )
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);
147 double ***
Bx, ***By, ***Bz;
148 double ***bx, ***by, ***bz;
155 Ap = grid[
IDIR].
A; Am = Ap - 1;
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]);
168 }
else if (side ==
X1_END){
169 Ap = grid[
IDIR].
A; Am = Ap - 1;
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]);
189 Ap = grid[
JDIR].
A; Am = Ap - 1;
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]);
202 }
else if (side ==
X2_END){
203 Ap = grid[
JDIR].
A; Am = Ap - 1;
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]);
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]);
235 }
else if (side ==
X3_END){
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]);
272 #define IEXT_LOOP(i) for (i = IOFFSET; i < NX1_TOT - IOFFSET; i++)
273 #define JEXT_LOOP(j) for (j = JOFFSET; j < NX2_TOT - JOFFSET; j++)
274 #define KEXT_LOOP(k) for (k = KOFFSET; k < NX3_TOT - KOFFSET; k++)
277 double *A1, *A2, *A3;
278 double ***
Bx, ***By, ***Bz;
279 double ***bx, ***by, ***bz;
295 #if GEOMETRY != SPHERICAL
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]); )
301 By[k][j][i] = (A2[j]*by[k][j][i] + A2[j-1]*by[k][j - 1][i])/
303 Bz[
k][
j][
i] = 0.5*(bz[
k][
j][
i] + bz[k - 1][
j][
i]);)
306 }
else if (side ==
X1_END){
308 #if GEOMETRY != SPHERICAL
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]); )
314 By[k][j][i] = (A2[j]*by[k][j][i] + A2[j-1]*by[k][j - 1][i])/
316 Bz[
k][
j][
i] = 0.5*(bz[
k][
j][
i] + bz[k - 1][
j][
i]);)
327 #if GEOMETRY == CARTESIAN
328 D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
330 Bz[
k][
j][
i] = 0.5*(bz[
k][
j][
i] + bz[k - 1][
j][
i]); )
332 D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
334 Bz[
k][
j][
i] = 0.5*(bz[
k][
j][
i] + bz[k - 1][
j][
i]); )
336 D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
339 Bz[
k][
j][
i] = 0.5*(bz[
k][
j][
i] + bz[k - 1][
j][
i]);)
342 }
else if (side ==
X2_END){
344 #if GEOMETRY == CARTESIAN
345 D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
347 Bz[
k][
j][
i] = 0.5*(bz[
k][
j][
i] + bz[k - 1][
j][
i]); )
349 D_EXPAND( Bx[k][j][i] = 0.5*(bx[k][j][i] + bx[k][j][i - 1]); ,
351 Bz[
k][
j][
i] = 0.5*(bz[
k][
j][
i] + bz[k - 1][
j][
i]); )
353 D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
356 Bz[
k][
j][
i] = 0.5*(bz[
k][
j][
i] + bz[k - 1][
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]); ,
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]); ,
377 D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
379 By[
k][
j][
i] = 0.5*(by[
k][
j][
i] + by[
k][j - 1][
i]); ,
383 D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
385 By[
k][
j][
i] = (A2[
j]*by[
k][
j][
i] + A2[j-1]*by[
k][j - 1][
i])/
390 }
else if (side ==
X3_END){
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]); ,
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]); ,
402 D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
404 By[
k][
j][
i] = 0.5*(by[
k][
j][
i] + by[
k][j - 1][
i]); ,
408 D_EXPAND(Bx[k][j][i] = (A1[i]*bx[k][j][i] + A1[i-1]*bx[k][j][i - 1])/
410 By[
k][
j][
i] = (A2[
j]*by[
k][
j][
i] + A2[j-1]*by[
k][j - 1][
i])/
#define X3_BEG
Boundary region at X3 beg.
#define X1_BEG
Boundary region at X1 beg.
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
long int NX2_TOT
Total number of zones in the X2 direction (boundaries included) for the local processor.
double **** Vc
The main four-index data array used for cell-centered primitive variables.
#define X3_END_LOOP(k, j, i)
#define X1_END
Boundary region at X1 end.
#define X2_BEG_LOOP(k, j, i)
#define X2_END_LOOP(k, j, i)
#define X2_END
Boundary region at X2 end.
long int NX3_TOT
Total number of zones in the X3 direction (boundaries included) for the local processor.
long int IEND
Upper grid index of the computational domain in the the X1 direction for the local processor...
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 th...
#define X1_END_LOOP(k, j, i)
void CT_AverageTransverseMagField(const Data *d, int side, Grid *grid)
Compute the transverse component of the volume-average magnetic field from the staggered components i...
#define X3_END
Boundary region at X3 end.
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;)
#define X1_BEG_LOOP(k, j, i)
long int KBEG
Lower grid index of the computational domain in the the X3 direction for the local processor...
#define X2_BEG
Boundary region at X2 beg.
long int KEND
Upper grid index of the computational domain in the the X3 direction for the local processor...
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
long int JEND
Upper grid index of the computational domain in the the X2 direction for the local processor...
#define X3_BEG_LOOP(k, j, i)
long int NX1_TOT
Total number of zones in the X1 direction (boundaries included) for the local processor.
long int IBEG
Lower grid index of the computational domain in the the X1 direction for the local processor...
double * A
Right interface area, A[i] = .
void CT_AverageMagneticField(double ****bf, double ****UU, Grid *grid)