86 double dtdx,
scrh, th, ct,
s;
93 double **rhs = state->
rhs;
94 double **flux = state->
flux;
95 double **vh = state->
vh;
96 double **vp = state->
vp;
97 double **vm = state->
vm;
101 double **Bg0, **
wA, w, wp, vphi, phi_c;
103 double vc[
NVAR], *vg;
109 #if (PHYSICS == MHD) && (BACKGROUND_FIELD == YES)
127 #if GEOMETRY == SPHERICAL
132 for (i = beg; i <= end; i++) {
139 #if GEOMETRY == CARTESIAN
146 #if (defined FARGO && !defined SHEARINGBOX)
150 #elif GEOMETRY == CYLINDRICAL
154 NFLX_LOOP(nv) vc[nv] = 0.5*(vp[i][nv] + vm[i][nv]);
162 #elif GEOMETRY == POLAR
166 NFLX_LOOP(nv) vc[nv] = 0.5*(vp[i][nv] + vm[i][nv]);
168 #if (defined FARGO) || (ROTATING_FRAME == YES)
174 rhs[
i][
MX1] += dt*( vc[
RHO]*vphi*vphi
177 IF_DUST(vphi_d = SELECT(0.0, vc[VX3_D], 0.0);
178 #if (defined FARGO) || (ROTATING_FRAME == YES)
181 rhs[
i][MX1_D] += dt*vc[RHO_D]*vphi_d*vphi_d*r_1;)
187 NFLX_LOOP(nv) vc[nv] = 0.5*(vp[i][nv] + vm[i][nv]);
188 vphi = SELECT(0.0, 0.0, vc[iVPHI]);
189 #if (defined FARGO) || (ROTATING_FRAME == YES)
195 Sm = vc[
RHO]*(EXPAND( 0.0, + vc[
VX2]*vc[
VX2], + vphi*vphi));
198 rhs[
i][
MX1] += dt*Sm*r_1;
200 IF_DUST(vphi_d = SELECT(0.0, 0.0, vc[VX3_D]);
201 #if (defined FARGO) || (ROTATING_FRAME == YES)
204 Sm_d = vc[RHO_D]*(EXPAND( 0.0, + vc[VX2_D]*vc[VX2_D], + vphi_d*vphi_d));
205 rhs[
i][MX1_D] += dt*Sm_d*r_1;)
214 IF_DUST (rhs[i][iMPHI_D] -= w*rhs[i][RHO_D];)
222 #if (BODY_FORCE & VECTOR)
226 IF_ENERGY(rhs[i][ENG] += dt*0.5*(flux[i][RHO] + flux[i-1][RHO])*g[
IDIR];)
230 rhs[
i][
MX1] -= dtdx*vg[
RHO]*(phi_p[
i] - phi_p[i-1]);
231 IF_DUST (rhs[i][MX1_D] -= dtdx*vg[RHO_D]*(phi_p[i] - phi_p[i-1]);)
233 rhs[
i][ENG] -= phi_c*rhs[
i][
RHO];)
242 #if GEOMETRY == POLAR
244 #elif GEOMETRY == SPHERICAL
245 scrh = r_1 = 0.5*(x1p[
i]*x1p[
i] - x1p[i-1]*x1p[i-1])/dV1[i];
247 for (j = beg; j <= end; j++) {
248 dtdx = dt/dx2[
j]*
scrh;
254 #if GEOMETRY != SPHERICAL
267 #elif GEOMETRY == SPHERICAL
275 vphi = SELECT(0.0, 0.0, vc[iVPHI]);
276 #if (defined FARGO) || (ROTATING_FRAME == YES)
282 Sm = vc[
RHO]*(EXPAND( 0.0, - vc[
iVTH]*vc[
iVR], + ct*vphi*vphi));
285 rhs[
j][
MX2] += dt*Sm*r_1;
287 IF_DUST(vphi_d = SELECT(0.0, 0.0, vc[VX3_D]);
288 #if (defined FARGO) || (ROTATING_FRAME == YES)
291 Sm_d = vc[RHO_D]*(EXPAND( 0.0, - vc[VX2_D]*vc[VX1_D], + ct*vphi_d*vphi_d));
292 rhs[
j][MX2_D] += dt*Sm_d*r_1;)
299 rhs[j][
MX3] -= w*rhs[j][RHO];
300 IF_DUST (rhs[j][MX3_D] -= w*rhs[j][RHO_D];)
301 IF_ENERGY(rhs[j][ENG] -= w*(rhs[j][
MX3] + 0.5*w*rhs[j][RHO]);)
314 IF_ENERGY(rhs[j][ENG] += dt*0.5*(flux[j][RHO] + flux[j-1][RHO])*g[
JDIR];)
318 rhs[
j][
MX2] -= dtdx*vg[
RHO]*(phi_p[
j] - phi_p[j-1]);
319 IF_DUST (rhs[j][MX2_D] -= dtdx*vg[RHO_D]*(phi_p[j] - phi_p[j-1]);)
321 rhs[
j][ENG] -= phi_c*rhs[
j][
RHO];)
328 #if GEOMETRY == SPHERICAL
329 r_1 = 0.5*(x1p[
i]*x1p[
i] - x1p[i-1]*x1p[i-1])/dV1[i];
330 scrh = dt*r_1*dx2[
j]/dV2[
j];
333 for (k = beg; k <= end; k++) {
337 #if (GEOMETRY == CARTESIAN) || (GEOMETRY == POLAR)
344 #if (defined FARGO && !defined SHEARINGBOX)
347 IF_DUST (rhs[k][MX2_D] -= w*rhs[k][RHO_D];)
348 IF_ENERGY(rhs[k][ENG] -= w*(rhs[k][
MX2] + 0.5*w*rhs[k][RHO]);)
360 IF_ENERGY(rhs[k][ENG] += dt*0.5*(flux[k][RHO] + flux[k-1][RHO])*g[
KDIR];)
364 rhs[
k][
MX3] -= dtdx*vg[
RHO]*(phi_p[
k] - phi_p[k-1]);
365 IF_DUST (rhs[k][MX3_D] -= dtdx*vg[RHO_D]*(phi_p[k] - phi_p[k-1]);)
367 rhs[
k][ENG] -= phi_c*rhs[
k][
RHO];)
378 for (i = beg; i <= end; i++) {
390 Dust_DragForce(state, beg, end, dt, grid);
397 #if (PHYSICS == MHD) && (DIVB_CONTROL == EIGHT_WAVES)
398 for (i = beg; i <= end; i++) {
399 EXPAND(rhs[i][
MX1] += dt*state->
src[i][
MX1]; ,
403 EXPAND(rhs[i][
BX1] += dt*state->
src[i][
BX1]; ,
407 rhs[
i][ENG] += dt*state->
src[
i][ENG];
416 #if (defined GLM_MHD) && (GLM_EXTENDED == YES)
#define IF_ROTATING_FRAME(a)
double ** vh
Primitive state at n+1/2 (only for one step method)
double ** flux
upwind flux computed with the Riemann solver
double ** visc_src
Viscosity source term.
void GLM_ExtendedSource(const State_1D *state, double dt, int beg, int end, Grid *grid)
double ** rhs
Conservative right hand side.
double ** GetBackgroundField(int beg, int end, int where, Grid *grid)
double ** vp
prim vars at i+1/2 edge, vp[i] = vL(i+1/2)
int g_i
x1 grid index when sweeping along the x2 or x3 direction.
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.
#define SB_OMEGA
Disk local orbital frequency .
double BodyForcePotential(double, double, double)
#define TotBB(v, b0, a, b)
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
#define SB_Q
The shear parameter, .
double * r_1
Geometrical factor 1/r.
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
double ** FARGO_GetVelocity(void)
void BodyForceVector(double *, double *, double, double, double)
double * ct
Geometrical factor cot(theta).
double * press
Upwind pressure term computed with the Riemann solver.
double * A
Right interface area, A[i] = .