8 #define EX(k,j,i) (vz[k][j][i]*By[k][j][i] - vy[k][j][i]*Bz[k][j][i])
9 #define EY(k,j,i) (vx[k][j][i]*Bz[k][j][i] - vz[k][j][i]*Bx[k][j][i])
10 #define EZ(k,j,i) (vy[k][j][i]*Bx[k][j][i] - vx[k][j][i]*By[k][j][i])
43 for (k = Z1->
kbeg; k <= Z1->kend; k++){
44 for (j = Z1->
jbeg; j <= Z1->jend; j++){
45 for (i = Z1->
ibeg; i <= Z1->iend; i++){
81 #define dEx_dyp(k,j,i) (emf->exj[k][j][i] - EX(k,j,i))
82 #define dEx_dzp(k,j,i) (emf->exk[k][j][i] - EX(k,j,i))
84 #define dEy_dxp(k,j,i) (emf->eyi[k][j][i] - EY(k,j,i))
85 #define dEy_dzp(k,j,i) (emf->eyk[k][j][i] - EY(k,j,i))
87 #define dEz_dxp(k,j,i) (emf->ezi[k][j][i] - EZ(k,j,i))
88 #define dEz_dyp(k,j,i) (emf->ezj[k][j][i] - EZ(k,j,i))
90 #define dEx_dym(k,j,i) (EX(k,j,i) - emf->exj[k][j-1][i])
91 #define dEx_dzm(k,j,i) (EX(k,j,i) - emf->exk[k-1][j][i])
93 #define dEy_dxm(k,j,i) (EY(k,j,i) - emf->eyi[k][j][i-1])
94 #define dEy_dzm(k,j,i) (EY(k,j,i) - emf->eyk[k-1][j][i])
96 #define dEz_dxm(k,j,i) (EZ(k,j,i) - emf->ezi[k][j][i-1])
97 #define dEz_dym(k,j,i) (EZ(k,j,i) - emf->ezj[k][j-1][i])
101 signed char sx, sy, sz;
102 double ***vx, ***vy, ***vz;
103 double ***
Bx, ***By, ***Bz;
113 for (k = emf->
kbeg; k <= emf->kend + KOFFSET; k++){
114 for (j = emf->
jbeg; j <= emf->jend + 1 ; j++){
115 for (i = emf->
ibeg; i <= emf->iend + 1 ; i++){
118 sy = emf->
svy[k][j][i]; ,
119 sz = emf->
svz[k][j][i];)
122 ju = sy > 0 ? j:j+1; ,
123 ku = sz > 0 ? k:k+1;)
233 #define dPP(a,x,y,i,j,k) (a[k][j][i] + 0.5*(d##a##_##d##x[k][j][i] + \
234 d##a##_##d##y[k][j][i]))
236 #define dPM(a,x,y,i,j,k) (a[k][j][i] + 0.5*(d##a##_##d##x[k][j][i] - \
237 d##a##_##d##y[k][j][i]))
239 #define dMM(a,x,y,i,j,k) (a[k][j][i] - 0.5*(d##a##_##d##x[k][j][i] + \
240 d##a##_##d##y[k][j][i]))
242 #define dMP(a,x,y,i,j,k) (a[k][j][i] - 0.5*(d##a##_##d##x[k][j][i] - \
243 d##a##_##d##y[k][j][i]))
247 #define dP(a,x,i,j,k) (a[k][j][i] + 0.5*(d##a##_##d##x[k][j][i]))
248 #define dM(a,x,i,j,k) (a[k][j][i] - 0.5*(d##a##_##d##x[k][j][i]))
253 double a_xp, a_yp, a_zp;
254 double a_xm, a_ym, a_zm;
256 double eSE, eSW, eNW, eNE;
257 double bS, bN, bW, bE;
260 double *xr, *yr, *zr;
263 double ***bx, ***by, ***bz;
264 double ***vx, ***vy, ***vz;
265 double ***dvx_dx, ***dvx_dy, ***dvx_dz;
266 double ***dvy_dx, ***dvy_dy, ***dvy_dz;
267 double ***dvz_dx, ***dvz_dy, ***dvz_dz;
269 double ***dbx_dy, ***dbx_dz;
270 double ***dby_dx, ***dby_dz;
271 double ***dbz_dx, ***dbz_dy;
288 for (k = emf->
kbeg; k <= emf->kend; k++){ kp = k + 1;
289 for (j = emf->
jbeg; j <= emf->jend; j++){ jp = j + 1;
290 for (i = emf->
ibeg; i <= emf->iend; i++){ ip = i + 1;
296 a_xp =
MAX(emf->
SxR[k][j][i], emf->
SxR[k][jp][i]);
297 a_xm =
MAX(emf->
SxL[k][j][i], emf->
SxL[k][jp][i]);
298 a_yp =
MAX(emf->
SyR[k][j][i], emf->
SyR[k][j][ip]);
299 a_ym =
MAX(emf->
SyL[k][j][i], emf->
SyL[k][j][ip]);
301 bS =
dP(bx, y, i , j , k); bW =
dP(by, x, i , j , k);
302 bN =
dM(bx, y, i , jp, k); bE =
dM(by, x, ip, j , k);
304 #if BACKGROUND_FIELD == YES
306 bS += B0[0]; bW += B0[1];
307 bN += B0[0]; bE += B0[1];
310 eSW =
dPP(vy,x,y,i ,j ,k)*bS -
dPP(vx,x,y,i ,j ,k)*bW;
311 eSE =
dMP(vy,x,y,ip,j ,k)*bS -
dMP(vx,x,y,ip,j ,k)*bE;
312 eNE =
dMM(vy,x,y,ip,jp,k)*bN -
dMM(vx,x,y,ip,jp,k)*bE;
313 eNW =
dPM(vy,x,y,i ,jp,k)*bN -
dPM(vx,x,y,i ,jp,k)*bW;
315 emf->
ez[
k][
j][
i] = a_xp*a_yp*eSW + a_xm*a_yp*eSE
316 + a_xm*a_ym*eNE + a_xp*a_ym*eNW;
317 emf->
ez[
k][
j][
i] /= (a_xp + a_xm)*(a_yp + a_ym);
318 emf->
ez[
k][
j][
i] -= a_yp*a_ym*(bN - bS)/(a_yp + a_ym);
319 emf->
ez[
k][
j][
i] += a_xp*a_xm*(bE - bW)/(a_xp + a_xm);
326 a_xp =
MAX(emf->
SyR[k][j][i], emf->
SyR[kp][j][i]);
327 a_xm =
MAX(emf->
SyL[k][j][i], emf->
SyL[kp][j][i]);
328 a_yp =
MAX(emf->
SzR[k][j][i], emf->
SzR[k][jp][i]);
329 a_ym =
MAX(emf->
SzL[k][j][i], emf->
SzL[k][jp][i]);
331 bS =
dP(by, z, i, j, k ); bW =
dP(bz, y, i, j , k);
332 bN =
dM(by, z, i, j, kp); bE =
dM(bz, y, i, jp, k);
334 #if BACKGROUND_FIELD == YES
336 bS += B0[1]; bW += B0[2];
337 bN += B0[1]; bE += B0[2];
340 eSW =
dPP(vz,y,z,i,j ,k )*bS -
dPP(vy,y,z,i ,j ,k )*bW;
341 eSE =
dMP(vz,y,z,i,jp,k )*bS -
dMP(vy,y,z,i ,jp,k )*bE;
342 eNE =
dMM(vz,y,z,i,jp,kp)*bN -
dMM(vy,y,z,i ,jp,kp)*bE;
343 eNW =
dPM(vz,y,z,i,j ,kp)*bN -
dPM(vy,y,z,i ,j ,kp)*bW;
345 emf->
ex[
k][
j][
i] = a_xp*a_yp*eSW + a_xm*a_yp*eSE
346 + a_xm*a_ym*eNE + a_xp*a_ym*eNW;
347 emf->
ex[
k][
j][
i] /= (a_xp + a_xm)*(a_yp + a_ym);
348 emf->
ex[
k][
j][
i] -= a_yp*a_ym*(bN - bS)/(a_yp + a_ym);
349 emf->
ex[
k][
j][
i] += a_xp*a_xm*(bE - bW)/(a_xp + a_xm);
355 a_xp =
MAX(emf->
SzR[k][j][i], emf->
SzR[k][j][ip]);
356 a_xm =
MAX(emf->
SzL[k][j][i], emf->
SzL[k][j][ip]);
357 a_yp =
MAX(emf->
SxR[k][j][i], emf->
SxR[kp][j][i]);
358 a_ym =
MAX(emf->
SxL[k][j][i], emf->
SxL[kp][j][i]);
360 bS =
dP(bz, x, i , j, k); bW =
dP(bx, z, i, j, k);
361 bN =
dM(bz, x, ip, j, k); bE =
dM(bx, z, i, j, kp);
363 #if BACKGROUND_FIELD == YES
365 bS += B0[2]; bW += B0[0];
366 bN += B0[2]; bE += B0[0];
369 eSW =
dPP(vx,z,x,i, j,k )*bS -
dPP(vz,z,x,i ,j ,k )*bW;
370 eSE =
dMP(vx,z,x,i, j,kp)*bS -
dMP(vz,z,x,i ,j ,kp)*bE;
371 eNE =
dMM(vx,z,x,ip,j,kp)*bN -
dMM(vz,z,x,ip,j ,kp)*bE;
372 eNW =
dPM(vx,z,x,ip,j,k )*bN -
dPM(vz,z,x,ip,j ,k )*bW;
374 emf->
ey[
k][
j][
i] = a_xp*a_yp*eSW + a_xm*a_yp*eSE
375 + a_xm*a_ym*eNE + a_xp*a_ym*eNW;
376 emf->
ey[
k][
j][
i] /= (a_xp + a_xm)*(a_yp + a_ym);
377 emf->
ey[
k][
j][
i] -= a_yp*a_ym*(bN - bS)/(a_yp + a_ym);
378 emf->
ey[
k][
j][
i] += a_xp*a_xm*(bE - bW)/(a_xp + a_xm);
413 double bx2, by2, bz2;
414 double vx4, vy4, vz4;
415 double ***vx, ***vy, ***vz;
416 double ***bx, ***by, ***bz;
422 for (k = emf->
kbeg; k <= emf->kend; k++){ kp = k + 1;
423 for (j = emf->
jbeg; j <= emf->jend; j++){ jp = j + 1;
424 for (i = emf->
ibeg; i <= emf->iend; i++){ ip = i + 1;
430 vx4 = 0.25*(vx[
k][
j][
i] + vx[
k][
j][ip] + vx[
k][jp][
i] + vx[
k][jp][ip]);
431 vy4 = 0.25*(vy[
k][
j][
i] + vy[
k][
j][ip] + vy[
k][jp][
i] + vy[
k][jp][ip]);
433 bx2 = 0.5*(bx[
k][
j][
i] + bx[
k][jp][
i]);
434 by2 = 0.5*(by[
k][
j][
i] + by[
k][
j][ip]);
436 emf->
ez[
k][
j][
i] = vy4*bx2 - vx4*by2;
443 vy4 = 0.25*(vy[
k][
j][
i] + vy[
k][jp][
i] + vy[kp][
j][
i] + vy[kp][jp][
i]);
444 vz4 = 0.25*(vz[
k][
j][
i] + vz[
k][jp][
i] + vz[kp][
j][
i] + vz[kp][jp][
i]);
446 by2 = 0.5*(by[
k][
j][
i] + by[kp][
j][
i]);
447 bz2 = 0.5*(bz[
k][
j][
i] + bz[
k][jp][
i]);
449 emf->
ex[
k][
j][
i] = vz4*by2 - vy4*bz2;
455 vz4 = 0.25*(vz[
k][
j][
i] + vz[kp][
j][
i] + vz[
k][
j][ip] + vz[kp][
j][ip]);
456 vx4 = 0.25*(vx[
k][
j][
i] + vx[kp][
j][
i] + vx[
k][
j][ip] + vx[kp][
j][ip]);
458 bz2 = 0.5*(bz[
k][
j][
i] + bz[
k][
j][ip]);
459 bx2 = 0.5*(bx[
k][
j][
i] + bx[kp][
j][
i]);
461 emf->
ey[
k][
j][
i] = vx4*bz2 - vz4*bx2;
double **** Vs
The main four-index data array used for face-centered staggered magnetic fields.
#define dP(a, x, i, j, k)
#define dPP(a, x, y, i, j, k)
void BackgroundField(double x1, double x2, double x3, double *B0)
int g_intStage
Gives the current integration stage of the time stepping method (predictor = 0, 1st corrector = 1...
double **** Vc
The main four-index data array used for cell-centered primitive variables.
void CT_EMF_ArithmeticAverage(const EMF *Z1, const double w)
Compute arithmetic average of EMF at cell edges.
void CT_EMF_HLL_Solver(const Data *d, const EMF *emf, Grid *grid)
Solve 2D Riemann problem for induction equation.
double *** ezi
Ez available at x-faces (i+1/2);.
double *** exk
Ex available at z-faces (k+1/2);.
void CT_EMF_CMUSCL_Average(const Data *d, const EMF *emf, Grid *grid)
–
double *** eyi
Ey available at x-faces (i+1/2);.
void CT_EMF_IntegrateToCorner(const Data *d, const EMF *emf, Grid *grid)
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;)
double *** eyk
Ey available at z-faces (k+1/2);.
double *** ezj
Ez available at y-faces (j+1/2);.
#define dM(a, x, i, j, k)
#define dPM(a, x, y, i, j, k)
double *** exj
Ex available at y-faces (j+1/2);.
#define dMP(a, x, y, i, j, k)
#define dMM(a, x, y, i, j, k)