25 int beg,
int end,
Grid *grid)
65 static double *wp, *wm, *Bxp, *Bxm;
75 #if PARABOLIC_FLUX != EXPLICIT
88 for (i = beg-1; i <= end+1; i++){
94 for (i = beg; i <= end+1; i++){
95 dB =
MC(wm[i], wm[i-1]);
96 dB +=
MC(wp[i], wp[i-1]);
109 for (i = beg; i <= end; i++){
110 for (nv =
NVAR; nv--; ){
111 VL[
i][nv] = state->
vL[
i][nv];
112 VR[
i][nv] = state->
vR[
i][nv];
118 VR[
i][
BXn] = Bxm[i+1];
132 #if COMPUTE_DIVB == YES
150 double dx, dy, dz, dtdx;
169 int beg,
int end,
Grid *grid)
183 double *Ar, *Ath, **bgf;
184 double *rhs, *v, *Bm, *pm;
185 static double *
divB, *dpsi, *Bp, *pp;
206 for (i = beg - 1; i <= end; i++) {
219 #if GEOMETRY == CARTESIAN
221 for (i = beg; i <= end; i++) {
222 divB[
i] = (Bp[
i] - Bm[
i])/GG->
dx[i];
223 dpsi[i] = (pp[i] - pm[i])/GG->
dx[
i];
226 #elif GEOMETRY == CYLINDRICAL
230 for (i = beg; i <= end; i++) {
231 divB[
i] = (Bp[
i]*Ar[
i] - Bm[
i]*Ar[i - 1])/GG->
dV[i];
232 dpsi[i] = (pp[i] - pm[i])/GG->
dx[
i];
234 }
else if (g_dir ==
JDIR){
235 for (i = beg; i <= end; i++) {
236 divB[
i] = (Bp[
i] - Bm[
i])/GG->
dx[i];
237 dpsi[i] = (pp[i] - pm[i])/GG->
dx[
i];
241 #elif GEOMETRY == POLAR
245 for (i = beg; i <= end; i++) {
246 divB[
i] = (Bp[
i]*Ar[
i] - Bm[
i]*Ar[i - 1])/GG->
dV[i];
247 dpsi[i] = (pp[i] - pm[i])/GG->
dx[
i];
249 }
else if (g_dir ==
JDIR){
251 for (i = beg; i <= end; i++) {
252 divB[
i] = (Bp[
i] - Bm[
i])/(r*GG->
dx[i]);
253 dpsi[
i] = (pp[
i] - pm[
i])/(r*GG->
dx[i]);
255 }
else if (g_dir ==
KDIR){
256 for (i = beg; i <= end; i++) {
257 divB[
i] = (Bp[
i] - Bm[
i])/GG->
dx[i];
258 dpsi[i] = (pp[i] - pm[i])/GG->
dx[
i];
262 #elif GEOMETRY == SPHERICAL
266 for (i = beg; i <= end; i++) {
267 divB[
i] = (Bp[
i]*Ar[
i] - Bm[
i]*Ar[i - 1])/GG->
dV[i];
268 dpsi[i] = (pp[i] - pm[i])/GG->
dx[
i];
270 }
else if (g_dir ==
JDIR){
273 for (i = beg; i <= end; i++) {
274 divB[
i] = (Bp[
i]*Ath[
i] - Bm[
i]*Ath[i - 1]) /
276 dpsi[
i] = (pp[
i] - pm[
i])/(r*GG->
dx[i]);
279 }
else if (g_dir ==
KDIR){
282 for (i = beg; i <= end; i++) {
283 divB[
i] = (Bp[
i] - Bm[
i])/(r*s*GG->
dx[i]);
284 dpsi[
i] = (pp[
i] - pm[
i])/(r*s*GG->
dx[i]);
290 #if BACKGROUND_FIELD == YES
298 for (i = beg; i <= end; i++) {
303 EXPAND(bx = v[
BX1]; ,
308 EXPAND(bx += bgf[i][
BX1]; ,
313 EXPAND (rhs[
MX1] -= dt*divB[i]*bx; ,
314 rhs[
MX2] -= dt*divB[
i]*by; ,
315 rhs[
MX3] -= dt*divB[
i]*bz;)
318 rhs[ENG] -= dt*v[
BXn]*dpsi[i];
336 double *x1, *x2, *x3;
337 static double **v, **lambda, *cs2;
398 #elif PHYSICS == RMHD
403 MPI_Allreduce (&
glm_ch, &gmaxc, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
408 #if COMPUTE_DIVB == YES
420 static int old_nstep = -1;
double ** v
Cell-centered primitive varables at the base time level, v[i] = .
void GLM_Solve(const State_1D *state, double **VL, double **VR, int beg, int end, Grid *grid)
double ** vh
Primitive state at n+1/2 (only for one step method)
double ** flux
upwind flux computed with the Riemann solver
void Eigenvalues(double **v, double *csound2, double **lambda, int beg, int end)
void GLM_ExtendedSource(const State_1D *state, double dt, int beg, int end, Grid *grid)
double ** rhs
Conservative right hand side.
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...
double **** Vc
The main four-index data array used for cell-centered primitive variables.
double ** vR
Primitive variables to the right of the interface, .
#define ARRAY_3D(nx, ny, nz, type)
double ** GetBackgroundField(int beg, int end, int where, Grid *grid)
void SoundSpeed2(double **v, double *cs2, double *h, int beg, int end, int pos, Grid *grid)
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
long int g_stepNumber
Gives the current integration step number.
int g_dir
Specifies the current sweep or direction of integration.
int g_j
x2 grid index when sweeping along the x1 or x3 direction.
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 SetIndexes(Index *indx, Grid *grid)
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
#define ARRAY_1D(nx, type)
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
double *** GLM_GetDivB(void)
double * bn
Face magentic field, bn = bx(i+1/2)
void BACKGROUND_FIELD(real x1, real x2, real x3, real *B0)
double ** vL
Primitive variables to the left of the interface, .
long int KBEG
Lower grid index of the computational domain in the the X3 direction for the local processor...
#define ARRAY_2D(nx, ny, type)
long int KEND
Upper grid index of the computational domain in the the X3 direction for the local processor...
double glm_ch
The propagation speed of divergence error.
long int JBEG
Lower grid index of the computational domain in the the X2 direction for the local processor...
#define GLM_ALPHA
Sets the damping rate of monopoles.
void GLM_Source(const Data_Arr Q, double dt, Grid *grid)
void GLM_Init(const Data *d, const Time_Step *Dts, Grid *grid)
long int JEND
Upper grid index of the computational domain in the the X2 direction for the local processor...
void GLM_ComputeDivB(const State_1D *state, Grid *grid)
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] = .