78 #define USE_PR_GRADIENT YES
80 #ifdef FINITE_DIFFERENCE
81 #define USE_PR_GRADIENT NO
83 #define USE_PR_GRADIENT YES
89 int beg,
int end,
double dt,
Grid *grid)
105 double dtdx, dtdV,
scrh, rhog;
110 double **rhs = state->
rhs;
111 double **flux = state->
flux, *p = state->
press;
112 double **vh = state->
vh, **vp = state->
vp, **vm = state->
vm;
116 double w, wp, vphi, phi_c;
118 static double **fA, *phi_p;
123 double **visc_src = state->
visc_src;
124 double **tc_flux = state->
tc_flux;
125 double **res_flux = state->
res_flux;
129 #if GEOMETRY != CARTESIAN
143 Dust_Solver(state, beg - 1, end, Dts->
cmax, grid);
155 #if USE_PR_GRADIENT == NO
156 for (i = beg - 1; i <= end; i++) flux[i][
MXn] += p[i];
163 #if (BODY_FORCE & POTENTIAL)
165 for (i = beg-1; i <= end; i++) {
169 for (j = beg-1; j <= end; j++) {
173 for (k = beg-1; k <= end; k++) {
183 #if (defined FARGO && !defined SHEARINGBOX) ||\
184 (ROTATING_FRAME == YES) || (BODY_FORCE & POTENTIAL)
185 TotalFlux(state, phi_p, beg-1, end, grid);
188 #if GEOMETRY == CARTESIAN
192 for (i = beg; i <= end; i++) {
197 NVAR_LOOP(nv) rhs[
i][nv] = -dtdx*(flux[
i][nv] - flux[i-1][nv]);
198 #if USE_PR_GRADIENT == YES
199 rhs[
i][
MX1] -= dtdx*(p[
i] - p[i-1]);
204 #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
209 #if VISCOSITY == EXPLICIT
210 rhs_entr += (visc_flux[
i][ENG] - visc_flux[i-1][ENG]);
211 rhs_entr -= EXPAND( vh[i][
VX1]*(visc_flux[i][
MX1] - visc_flux[i-1][
MX1]) ,
212 + vh[i][
VX2]*(visc_flux[i][
MX2] - visc_flux[i-1][
MX2]) ,
213 + vh[i][
VX3]*(visc_flux[i][
MX3] - visc_flux[i-1][
MX3]));
216 #if THERMAL_CONDUCTION == EXPLICIT
217 rhs_entr += (tc_flux[
i][ENG] - tc_flux[i-1][ENG]);
219 rhs[
i][ENTR] += rhs_entr*dtdx*rhog;
225 for (j = beg; j <= end; j++) {
230 NVAR_LOOP(nv) rhs[
j][nv] = -dtdx*(flux[
j][nv] - flux[j-1][nv]);
231 #if USE_PR_GRADIENT == YES
232 rhs[
j][
MX2] -= dtdx*(p[
j] - p[j-1]);
237 #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
242 #if VISCOSITY == EXPLICIT
243 rhs_entr += (visc_flux[
j][ENG] - visc_flux[j-1][ENG]);
244 rhs_entr -= EXPAND( vh[j][
VX1]*(visc_flux[j][
MX1] - visc_flux[j-1][
MX1]) ,
245 + vh[j][
VX2]*(visc_flux[j][
MX2] - visc_flux[j-1][
MX2]) ,
246 + vh[j][
VX3]*(visc_flux[j][
MX3] - visc_flux[j-1][
MX3]));
249 #if THERMAL_CONDUCTION == EXPLICIT
250 rhs_entr += (tc_flux[
j][ENG] - tc_flux[j-1][ENG]);
252 rhs[
j][ENTR] += rhs_entr*dtdx*rhog;
259 for (k = beg; k <= end; k++) {
264 NVAR_LOOP(nv) rhs[
k][nv] = -dtdx*(flux[
k][nv] - flux[k-1][nv]);
265 #if USE_PR_GRADIENT == YES
266 rhs[
k][
MX3] -= dtdx*(p[
k] - p[k-1]);
271 #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
275 #if VISCOSITY == EXPLICIT
276 rhs_entr += (visc_flux[
k][ENG] - visc_flux[k-1][ENG]);
277 rhs_entr -= EXPAND( vh[k][
VX1]*(visc_flux[k][
MX1] - visc_flux[k-1][
MX1]) ,
278 + vh[k][
VX2]*(visc_flux[k][
MX2] - visc_flux[k-1][
MX2]) ,
279 + vh[k][
VX3]*(visc_flux[k][
MX3] - visc_flux[k-1][
MX3]));
282 #if THERMAL_CONDUCTION == EXPLICIT
283 rhs_entr += (tc_flux[
k][ENG] - tc_flux[k-1][ENG]);
285 rhs[
k][ENTR] += rhs_entr*dtdx*rhog;
290 #elif GEOMETRY == CYLINDRICAL
293 double R, z, phi, R_1;
296 #error "DUST not implemented in CYLINDRICAL coordinates"
303 for (i = beg; i <= end; i++){
312 EXPAND(rhs[i][
iMR] = - dtdV*(fA[i][
iMR] - fA[i-1][
iMR]); ,
313 rhs[
i][
iMZ] = - dtdV*(fA[
i][
iMZ] - fA[i-1][
iMZ]); ,
316 rhs[i][
iMR] -= dtdx*(p[i] - p[i-1]);
319 EXPAND(rhs[i][
iBR] = - dtdV*(fA[i][
iBR] - fA[i-1][
iBR]); ,
320 rhs[
i][
iBZ] = - dtdV*(fA[
i][
iBZ] - fA[i-1][
iBZ]); ,
323 rhs[i][
iBR] = - dtdx*(flux[i][
iBR] - flux[i-1][
iBR]);
327 IF_ENERGY(rhs[i][ENG] = -dtdV*(fA[i][ENG] - fA[i-1][ENG]);)
328 NSCL_LOOP(nv) rhs[
i][nv] = -dtdV*(fA[
i][nv] - fA[i-1][nv]);
332 #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
333 NVAR_LOOP(nv) vc[nv] = 0.5*(vp[
i][nv] + vm[
i][nv]);
338 #if VISCOSITY == EXPLICIT
339 rhs_entr += (fvA[
i][ENG] - fvA[i-1][ENG]);
340 rhs_entr -= EXPAND( vc[
VX1]*(fvA[i][
MX1] - fvA[i-1][
MX1]) ,
341 + vc[
VX2]*(fvA[i][
MX2] - fvA[i-1][
MX2]) ,
342 + vc[
VX3]*(fvA[i][
MX3] - fvA[i-1][
MX3])*fabs(R_1));
344 rhs_entr -= EXPAND( vc[
VX1]*visc_src[i][
MX1],
345 + vc[
VX2]*visc_src[i][
MX2],
346 + vc[
VX3]*visc_src[i][
MX3]);
349 #if THERMAL_CONDUCTION == EXPLICIT
350 rhs_entr += (A[
i]*tc_flux[
i][ENG] - A[i-1]*tc_flux[i-1][ENG]);
352 rhs[
i][ENTR] += rhs_entr*dtdV*rhog;
358 for (j = beg; j <= end; j++){
363 NVAR_LOOP(nv) rhs[
j][nv] = -dtdx*(flux[
j][nv] - flux[j-1][nv]);
364 #if USE_PR_GRADIENT == YES
365 rhs[
j][
iMZ] += - dtdx*(p[
j] - p[j-1]);
370 #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
375 #if VISCOSITY == EXPLICIT
376 rhs_entr += (visc_flux[
j][ENG] - visc_flux[j-1][ENG]);
377 rhs_entr -= EXPAND( vc[
VX1]*(visc_flux[j][
MX1] - visc_flux[j-1][
MX1]) ,
378 + vc[
VX2]*(visc_flux[j][
MX2] - visc_flux[j-1][
MX2]) ,
379 + vc[
VX3]*(visc_flux[j][
MX3] - visc_flux[j-1][
MX3]));
381 #if THERMAL_CONDUCTION == EXPLICIT
382 rhs_entr += (tc_flux[
j][ENG] - tc_flux[j-1][ENG]);
384 rhs[
j][ENTR] += rhs_entr*dtdx*rhog;
390 #elif GEOMETRY == POLAR
399 for (i = beg; i <= end; i++) {
408 EXPAND(rhs[i][
iMR] = - dtdV*(fA[i][
iMR] - fA[i-1][
iMR])
409 - dtdx*(p[i] - p[i-1]); ,
413 EXPAND(rhs[i][
iBR] = -dtdV*(fA[i][
iBR] - fA[i-1][
iBR]); ,
417 rhs[i][
iBR] = -dtdx*(flux[i][
iBR] - flux[i-1][
iBR]);
421 IF_ENERGY(rhs[i][ENG] = -dtdV*(fA[i][ENG] - fA[i-1][ENG]);)
423 NSCL_LOOP(nv) rhs[
i][nv] = -dtdV*(fA[
i][nv] - fA[i-1][nv]);
426 rhs[
i][MX2_D] *= R_1;)
431 for (nv =
NVAR; nv--; ) vc[nv] = 0.5*(vp[i][nv] + vm[i][nv]);
436 #if VISCOSITY == EXPLICIT
437 rhs_entr += (fvA[
i][ENG] - fvA[i-1][ENG]);
438 rhs_entr -= EXPAND( vc[
VX1]*(fvA[i][
iMR] - fvA[i-1][
iMR]) ,
440 + vc[
VX3]*(fvA[i][
iMZ] - fvA[i-1][
iMZ]));
442 rhs_entr -= EXPAND( vc[
VX1]*visc_src[i][
MX1],
443 + vc[
VX2]*visc_src[i][
MX2],
444 + vc[
VX3]*visc_src[i][
MX3]);
447 #if THERMAL_CONDUCTION == EXPLICIT
448 rhs_entr += (A[
i]*tc_flux[
i][ENG] - A[i-1]*tc_flux[i-1][ENG]);
450 rhs[
i][ENTR] += rhs_entr*dtdV*rhog;
458 for (j = beg; j <= end; j++){
463 NVAR_LOOP(nv) rhs[
j][nv] = -dtdx*(flux[
j][nv] - flux[j-1][nv]);
464 rhs[
j][
iMPHI] -= dtdx*(p[
j] - p[j-1]);
468 #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
473 #if VISCOSITY == EXPLICIT
474 rhs_entr += (visc_flux[
j][ENG] - visc_flux[j-1][ENG]);
475 rhs_entr -= EXPAND( vh[j][
VX1]*(visc_flux[j][
MX1] - visc_flux[j-1][
MX1]) ,
476 + vh[j][
VX2]*(visc_flux[j][
MX2] - visc_flux[j-1][
MX2]) ,
477 + vh[j][
VX3]*(visc_flux[j][
MX3] - visc_flux[j-1][
MX3]));
479 rhs_entr -= EXPAND( vh[j][
VX1]*visc_src[j][
MX1],
480 + vh[j][
VX2]*visc_src[j][
MX2],
481 + vh[j][
VX3]*visc_src[j][
MX3]);
483 #if THERMAL_CONDUCTION == EXPLICIT
484 rhs_entr += (tc_flux[
j][ENG] - tc_flux[j-1][ENG]);
486 rhs[
j][ENTR] += rhs_entr*dtdx*rhog;
493 for (k = beg; k <= end; k++){
498 VAR_LOOP(nv) rhs[
k][nv] = -dtdx*(flux[
k][nv] - flux[k-1][nv]);
499 rhs[
k][
MX3] -= dtdx*(p[
k] - p[k-1]);
503 #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
508 #if VISCOSITY == EXPLICIT
509 rhs_entr += (visc_flux[
k][ENG] - visc_flux[k-1][ENG]);
510 rhs_entr -= EXPAND( vh[j][
VX1]*(visc_flux[k][
MX1] - visc_flux[k-1][
MX1]) ,
511 + vh[j][
VX2]*(visc_flux[k][
MX2] - visc_flux[k-1][
MX2]) ,
512 + vh[j][
VX3]*(visc_flux[k][
MX3] - visc_flux[k-1][
MX3]));
515 #if THERMAL_CONDUCTION == EXPLICIT
516 rhs_entr += (tc_flux[
k][ENG] - tc_flux[k-1][ENG]);
518 rhs[
k][ENTR] += rhs_entr*dtdx*rhog;
524 #elif GEOMETRY == SPHERICAL
526 double r_1, th,
s, s_1;
529 double vc[
NVAR], dVdx;
532 for (i = beg; i <= end; i++) {
559 rhs[i][
iMR] = - dtdV*(fA[i][
iMR] - fA[i-1][
iMR])
560 - dtdx*(p[i] - p[i-1]); ,
566 rhs[i][
iBR] = -dtdV*(fA[i][
iBR] - fA[i-1][
iBR]); ,
571 rhs[i][
iBR] = -dtdx*(flux[i][
iBR] - flux[i-1][
iBR]);
575 IF_ENERGY(rhs[i][ENG] = -dtdV*(fA[i][ENG] - fA[i-1][ENG]);)
577 NSCL_LOOP(nv) rhs[
i][nv] = -dtdV*(fA[
i][nv] - fA[i-1][nv]);
579 rhs[
i][MX3_D] *= r_1;)
584 NVAR_LOOP(nv) vc[nv] = 0.5*(vp[
i][nv] + vm[
i][nv]);
590 #if VISCOSITY == EXPLICIT
591 rhs_entr += (fvA[
i][ENG] - fvA[i-1][ENG]);
592 rhs_entr -= EXPAND( vc[
VX1]*(fvA[i][
iMR] - fvA[i-1][
iMR]) ,
596 rhs_entr -= EXPAND( vc[
VX1]*visc_src[i][
MX1],
597 + vc[
VX2]*visc_src[i][
MX2],
598 + vc[
VX3]*visc_src[i][
MX3]);
600 #if THERMAL_CONDUCTION == EXPLICIT
601 rhs_entr += (A[
i]*tc_flux[
i][ENG] - A[i-1]*tc_flux[i-1][ENG]);
603 rhs[
i][ENTR] += rhs_entr*dtdV*rhog;
610 r_1 = 0.5*(x1p[
i]*x1p[
i] - x1p[i-1]*x1p[i-1])/dV1[i];
612 for (j = beg; j <= end; j++){
623 rhs[j][
iMR] = - dtdV*(fA[j][
iMR] - fA[j-1][
iMR]); ,
625 - dtdx*(p[j] - p[j-1]); ,
630 rhs[j][
iBR] = -dtdV*(fA[j][
iBR] - fA[j-1][
iBR]); ,
639 IF_ENERGY(rhs[j][ENG] = -dtdV*(fA[j][ENG] - fA[j-1][ENG]);)
641 NSCL_LOOP(nv) rhs[
j][nv] = -dtdV*(fA[
j][nv] - fA[j-1][nv]);
644 rhs[
j][MX3_D] *= fabs(s_1);)
653 #if VISCOSITY == EXPLICIT
654 rhs_entr += (fvA[
j][ENG] - fvA[j-1][ENG]);
655 rhs_entr -= EXPAND( vc[
VX1]*(fvA[j][
iMR] - fvA[j-1][
iMR]) ,
659 rhs_entr -= EXPAND( vc[
VX1]*visc_src[j][
MX1],
660 + vc[
VX2]*visc_src[j][
MX2],
661 + vc[
VX3]*visc_src[j][
MX3]);
663 #if THERMAL_CONDUCTION == EXPLICIT
664 rhs_entr += (A[
j]*tc_flux[
j][ENG] - A[j-1]*tc_flux[j-1][ENG]);
666 rhs[
j][ENTR] += rhs_entr*dtdV*rhog;
673 r_1 = 0.5*(x1p[
i]*x1p[
i] - x1p[i-1]*x1p[i-1])/dV1[i];
674 scrh = dt*r_1*dx2[
j]/dV2[
j];
675 for (k = beg; k <= end; k++) {
680 VAR_LOOP(nv) rhs[
k][nv] = -dtdx*(flux[
k][nv] - flux[k-1][nv]);
681 rhs[
k][
iMPHI] -= dtdx*(p[
k] - p[k-1]);
685 #if (ENTROPY_SWITCH) && (PARABOLIC_FLUX & EXPLICIT)
690 #if VISCOSITY == EXPLICIT
691 rhs_entr += (visc_flux[
k][ENG] - visc_flux[k-1][ENG]);
692 rhs_entr -= EXPAND( vc[
VX1]*(visc_flux[k][
MX1] - visc_flux[k-1][
MX1]) ,
693 + vc[
VX2]*(visc_flux[k][
MX2] - visc_flux[k-1][
MX2]) ,
694 + vc[
VX3]*(visc_flux[k][
MX3] - visc_flux[k-1][
MX3]));
696 #if THERMAL_CONDUCTION == EXPLICIT
697 rhs_entr += (tc_flux[
k][ENG] - tc_flux[k-1][ENG]);
699 rhs[
k][ENTR] += rhs_entr*dtdx*rhog;
716 #if INTERNAL_BOUNDARY == YES
729 for (i = beg-1; i <= end; i++) {
733 #if GEOMETRY == POLAR || GEOMETRY == SPHERICAL
735 cl /= fabs(grid[
IDIR].xgc[
g_i]);
737 #if GEOMETRY == SPHERICAL
748 int beg,
int end,
Grid *grid)
793 double *x1, *x2, *x3;
794 double *x1p, *x2p, *x3p;
810 for (i = beg; i <= end; i++){
818 #if (defined FARGO && !defined SHEARINGBOX) || (ROTATING_FRAME == YES)
820 #if GEOMETRY == SPHERICAL
821 IF_FARGO(wp = 0.5*(wA[j][i] + wA[j][i+1]);)
822 R = x1p[i]*sin(x2[j]);
824 IF_FARGO(wp = 0.5*(wA[k][i] + wA[k][i+1]);)
829 flux[i][iMPHI] += wp*flux[i][
RHO];
834 #if (BODY_FORCE & POTENTIAL) && HAVE_ENERGY
835 flux[
i][ENG] += flux[
i][
RHO]*phi_p[
i];
839 for (j = beg; j <= end; j++){
845 #if GEOMETRY == SPHERICAL
846 #if (defined FARGO) || (ROTATING_FRAME == YES)
848 R = x1[
i]*sin(x2p[j]);
849 IF_FARGO (wp += 0.5*(wA[j][i] + wA[j+1][i]);)
852 flux[j][iMPHI] += wp*flux[j][
RHO];
856 #if (BODY_FORCE & POTENTIAL) && HAVE_ENERGY
857 flux[
j][ENG] += flux[
j][
RHO]*phi_p[
j];
862 for (k = beg; k <= end; k++) {
869 #if (GEOMETRY != SPHERICAL) && (defined FARGO && !defined SHEARINGBOX)
870 wp = 0.5*(wA[
k][
i] + wA[k+1][
i]);
872 flux[k][iMPHI] += wp*flux[k][
RHO];
875 #if (BODY_FORCE & POTENTIAL) && HAVE_ENERGY
876 flux[
k][ENG] += flux[
k][
RHO]*phi_p[
k];
#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 RightHandSideSource(const State_1D *state, Time_Step *Dts, int beg, int end, double dt, double *phi_p, Grid *grid)
double ** rhs
Conservative right hand side.
static void TotalFlux(const State_1D *, double *, int, int, Grid *)
double ** res_flux
Resistive flux (current)
void GetAreaFlux(const State_1D *state, double **fA, double **fvA, int beg, int end, Grid *grid)
double * cmax
Maximum signal velocity for hyperbolic eqns.
double ** visc_flux
Viscosity flux.
double ** vp
prim vars at i+1/2 edge, vp[i] = vL(i+1/2)
double inv_dta
Inverse of advection (hyperbolic) time step, .
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.
double ** tc_flux
Thermal conduction flux.
double BodyForcePotential(double, double, double)
void RightHandSide(const State_1D *state, Time_Step *Dts, int beg, int end, double dt, Grid *grid)
int g_k
x3 grid index when sweeping along the x1 or x2 direction.
double * r_1
Geometrical factor 1/r.
double ** vm
prim vars at i-1/2 edge, vm[i] = vR(i-1/2)
#define ARRAY_1D(nx, type)
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
double ** FARGO_GetVelocity(void)
void InternalBoundaryReset(const State_1D *, Time_Step *, int, int, Grid *)
double * press
Upwind pressure term computed with the Riemann solver.
#define ARRAY_2D(nx, ny, type)
double * inv_dxi
inverse of the distance between the center of two cells, inv_dxi = .
void AdvectFlux(const State_1D *state, int beg, int end, Grid *grid)
double * A
Right interface area, A[i] = .