47 double rhs_x, rhs_y, rhs_z;
48 double *dx, *dy, *dz, *A1, *A2, *dV1, *dV2;
49 double *r, *rp, *th, r_2;
50 double ***ex, ***ey, ***ez;
55 #if CHECK_DIVB_CONDITION == YES
61 #if UPDATE_VECTOR_POTENTIAL == YES
76 if (
g_intStage == 2) SB_CorrectEMF (emf, Bs, grid);
78 SB_CorrectEMF (emf, d->
Vs, grid);
99 for (k = emf->
kbeg + KOFFSET; k <= emf->kend; k++){
100 for (j = emf->
jbeg + 1 ; j <= emf->jend; j++){
101 for (i = emf->
ibeg ; i <= emf->iend; i++){
103 #if GEOMETRY == CARTESIAN
106 - dt/dy[j]*(ez[k][j][i] - ez[k][j - 1][i]) ,
107 + dt/dz[k]*(ey[k][j][i] - ey[k - 1][j][i]) );
109 #elif GEOMETRY == CYLINDRICAL
111 rhs_x = - dt/dy[
j]*(ez[
k][
j][
i] - ez[
k][j - 1][
i]);
113 #elif GEOMETRY == POLAR
116 - dt/(A1[i]*dy[j])*(ez[k][j][i] - ez[k][j - 1][i]) ,
117 + dt/dz[k] *(ey[k][j][i] - ey[k - 1][j][i]));
119 #elif GEOMETRY == SPHERICAL
123 - dt/(rp[i]*dV2[j])*(A2[j]*ez[k][j][i] - A2[j - 1]*ez[k][j - 1][i]) ,
124 + dt*dy[j]/(rp[i]*dV2[j]*dz[k])*(ey[k][j][i] - ey[k - 1][j][i]));
136 for (k = emf->
kbeg + KOFFSET; k <= emf->kend; k++){
137 for (j = emf->
jbeg ; j <= emf->jend; j++){
138 for (i = emf->
ibeg + 1 ; i <= emf->iend; i++){
140 #if GEOMETRY == CARTESIAN
142 rhs_y =
D_EXPAND( dt/dx[i]*(ez[k][j][i] - ez[k][j][i - 1]) ,
144 - dt/dz[k]*(ex[k][j][i] - ex[k - 1][j][i]) );
146 #elif GEOMETRY == CYLINDRICAL
148 rhs_y = dt/dV1[
i]*(A1[
i]*ez[
k][
j][
i] - A1[i - 1]*ez[
k][
j][i - 1]);
150 #elif GEOMETRY == POLAR
152 rhs_y =
D_EXPAND( dt/dx[i]*(ez[k][j][i] - ez[k][j][i - 1]) ,
154 - dt/dz[k]*(ex[k][j][i] - ex[k - 1][j][i]));
156 #elif GEOMETRY == SPHERICAL
159 + dt/(r[i]*dx[i])*(rp[i]*ez[k][j][i] - rp[i - 1]*ez[k][j][i - 1]) ,
160 , - dt/(r[i]*A2[j]*dz[k])*(ex[k][j][i] - ex[k - 1][j][i]));
173 for (k = emf->
kbeg ; k <= emf->kend; k++){
174 for (j = emf->
jbeg + 1; j <= emf->jend; j++){
175 for (i = emf->
ibeg + 1; i <= emf->iend; i++){
177 #if GEOMETRY == CARTESIAN
179 rhs_z = - dt/dx[
i]*(ey[
k][
j][
i] - ey[
k][
j][i - 1])
180 + dt/dy[j]*(ex[k][j][i] - ex[k][j - 1][i]);
182 #elif GEOMETRY == POLAR
184 rhs_z = - dt/dV1[
i]*(A1[
i]*ey[
k][
j][
i] - A1[i - 1]*ey[
k][
j][i - 1])
185 + dt/(r[i]*dy[j])*(ex[
k][
j][
i] - ex[
k][j - 1][
i]);
187 #elif GEOMETRY == SPHERICAL
189 rhs_z = - dt/(r[
i]*dx[
i])*(rp[i]*ey[k][j][i] - rp[i - 1]*ey[k][j][i - 1])
190 + dt/(r[
i]*dy[
j])*(ex[k][j][i] - ex[k][j - 1][i]);
210 double ***bx, ***by, ***bz;
211 double *dx, *dy, *dz;
212 double *Ar, *
s, *dmu, *r, *rp, *rm;
259 #if GEOMETRY == CARTESIAN
261 divB =
D_EXPAND( (bx[k][j][i] - bx[k][j][i-1])*dy[j]*dz[k] ,
262 + (by[k][j][i] - by[k][j-1][i])*dx[i]*dz[k] ,
263 + (bz[k][j][i] - bz[k-1][j][i])*dx[i]*dy[j] );
265 #elif GEOMETRY == CYLINDRICAL
267 divB = dy[
j]*(bx[
k][
j][
i]*Ar[
i] - bx[
k][
j][i-1]*Ar[i-1])
268 + fabs(r[i])*dx[
i]*(by[
k][
j][
i] - by[
k][j-1][
i]);
274 #elif GEOMETRY == POLAR
277 dz[k]*dy[j]*(bx[k][j][i]*Ar[i] - bx[k][j][i-1]*Ar[i-1]) ,
278 + dx[i]*dz[k]*(by[k][j][i] - by[k][j-1][i]) ,
279 + r[i]*dx[i]*dy[j]*(bz[k][j][i] - bz[k-1][j][i])
282 #elif GEOMETRY == SPHERICAL
285 dmu[j] *dz[k]*(bx[k][j][i]*Ar[i] - bx[k][j][i-1]*Ar[i-1]) ,
286 + r[i]*dx[i]*dz[k]*(by[k][j][i]*s[j] - by[k][j-1][i]*s[j-1]) ,
287 + r[i]*dx[i]*dy[j]*(bz[k][j][i] - bz[k-1][j][i])
292 dbmax =
MAX(dbmax,fabs(divB));
293 if (fabs(divB) > 1.e-6) {
294 print (
"! CHECK_DIVB: div(B) = %12.6e != 0, rank = %d, ijk = %d %d %d \n",
296 print (
"! rank =%d, g_intStage = %d, Beg, End = [%d, %d] [%d, %d] [%d, %d]\n",
299 D_EXPAND(
print (
"b1: %12.6e %12.6e\n",bx[k][j][i],bx[k][j][i-1]); ,
300 print (
"b2: %12.6e %12.6e\n",by[k][j][i],by[k][j-1][i]); ,
301 print (
"b3: %12.6e %12.6e\n",bz[k][j][i],bz[k-1][j][i]); )
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.
int g_intStage
Gives the current integration stage of the time stepping method (predictor = 0, 1st corrector = 1...
EMF * CT_GetEMF(const Data *d, Grid *grid)
void CT_Update(const Data *d, Data_Arr Bs, double dt, Grid *grid)
void CT_CheckDivB(double ***bf[], Grid *grid)
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 print(const char *fmt,...)
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;)
long int KBEG
Lower grid index of the computational domain in the the X3 direction for the local processor...
long int KEND
Upper grid index of the computational domain in the the X3 direction for the local processor...
void VectorPotentialUpdate(const Data *d, const void *vp, const State_1D *state, const Grid *grid)
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
#define QUIT_PLUTO(e_code)
long int JEND
Upper grid index of the computational domain in the the X2 direction for the local processor...
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] = .