37 #if DIVB_CONTROL == EIGHT_WAVES
47 double btx, bty, btz, bx, by, bz, vx, vy, vz;
52 static double *
divB, *vp;
65 for (i = is - 1; i <= ie; i++) {
76 #if GEOMETRY == CARTESIAN
78 for (i = is; i <= ie; i++) {
79 divB[
i] = (vp[
i] - vm[
i])/GG->
dx[i];
82 #elif GEOMETRY == CYLINDRICAL
86 for (i = is; i <= ie; i++) {
87 divB[
i] = (vp[
i]*Ar[
i] - vm[
i]*Ar[i - 1])/GG->
dV[i];
89 }
else if (g_dir ==
JDIR){
90 for (i = is; i <= ie; i++) {
91 divB[
i] = (vp[
i] - vm[
i])/GG->
dx[i];
95 #elif GEOMETRY == POLAR
99 for (i = is; i <= ie; i++) {
100 divB[
i] = (vp[
i]*Ar[
i] - vm[
i]*Ar[i - 1])/GG->
dV[i];
102 }
else if (g_dir ==
JDIR){
104 for (i = is; i <= ie; i++) {
105 divB[
i] = (vp[
i] - vm[
i])/(r*GG->
dx[i]);
107 }
else if (g_dir ==
KDIR){
108 for (i = is; i <= ie; i++) {
109 divB[
i] = (vp[
i] - vm[
i])/GG->
dx[i];
113 #elif GEOMETRY == SPHERICAL
117 for (i = is; i <= ie; i++) {
118 divB[
i] = (vp[
i]*Ar[
i] - vm[
i]*Ar[i - 1])/GG->
dV[i];
120 }
else if (g_dir ==
JDIR){
123 for (i = is; i <= ie; i++) {
124 divB[
i] = (vp[
i]*Ath[
i] - vm[
i]*Ath[i - 1]) /
127 }
else if (g_dir ==
KDIR){
130 for (i = is; i <= ie; i++) {
131 divB[
i] = (vp[
i] - vm[
i])/(r*s*GG->
dx[i]);
137 #if BACKGROUND_FIELD == YES
145 for (i = is; i <= ie; i++) {
150 EXPAND (vx = v[
VX1]; ,
154 EXPAND (bx = btx = v[
BX1]; ,
165 EXPAND (src[
MX1] = -divB[i]*btx; ,
166 src[
MX2] = -divB[
i]*bty; ,
167 src[
MX3] = -divB[
i]*btz;)
170 src[ENG] = -divB[i]*(EXPAND(vx*bx, +vy*by, +vz*bz));
172 EXPAND (src[
BX1] = -divB[i]*vx; ,
173 src[
BX2] = -divB[
i]*vy; ,
174 src[
BX3] = -divB[
i]*vz;)
180 int beg,
int end,
Grid *grid)
188 double vc[
NVAR], *A, *src, *vm;
190 static double *
divB, *vp;
206 for (i = beg - 1; i <= end; i++) {
215 #if GEOMETRY == CARTESIAN
217 for (i = beg; i <= end; i++) {
218 divB[
i] = (vp[
i] - vm[
i])/GG->
dx[i];
221 #elif GEOMETRY == CYLINDRICAL
224 for (i = beg; i <= end; i++) {
225 divB[
i] = (vp[
i]*A[
i] - vm[
i]*A[i - 1])/GG->
dV[i];
227 }
else if (g_dir ==
JDIR){
228 for (i = beg; i <= end; i++) {
229 divB[
i] = (vp[
i] - vm[
i])/GG->
dx[i];
233 #elif GEOMETRY == POLAR
236 for (i = beg; i <= end; i++) {
237 divB[
i] = (vp[
i]*A[
i] - vm[
i]*A[i - 1])/GG->
dV[i];
239 }
else if (g_dir ==
JDIR){
241 for (i = beg; i <= end; i++) {
242 divB[
i] = (vp[
i] - vm[
i])/(r*GG->
dx[i]);
244 }
else if (g_dir ==
KDIR){
245 for (i = beg; i <= end; i++) {
246 divB[
i] = (vp[
i] - vm[
i])/GG->
dx[i];
250 #elif GEOMETRY == SPHERICAL
253 for (i = beg; i <= end; i++) {
254 divB[
i] = (vp[
i]*A[
i] - vm[
i]*A[i - 1])/GG->
dV[i];
256 }
else if (g_dir ==
JDIR){
258 for (i = beg; i <= end; i++) {
259 divB[
i] = (vp[
i]*A[
i] - vm[
i]*A[i - 1])/(r*GG->
dV[i]);
261 }
else if (g_dir ==
KDIR){
264 for (i = beg; i <= end; i++) {
265 divB[
i] = (vp[
i] - vm[
i])/(r*s*GG->
dx[i]);
275 for (i = beg; i <= end; i++) {
278 for (nv =
NFLX; nv--; ) vc[nv] = state->
vh[i][nv];
282 EXPAND(src[
MX1] = -vc[BX1]*divB[i]; ,
286 EXPAND(src[BX1] = -vc[
VX1]*divB[i]; ,
291 src[ENG] = -vB*divB[i];
double ** vh
Primitive state at n+1/2 (only for one step method)
void HLL_DivBSource(const State_1D *state, double **Uhll, int beg, int end, Grid *grid)
double ** flux
upwind flux computed with the Riemann solver
double ** vR
Primitive variables to the right of the interface, .
double ** GetBackgroundField(int beg, int end, int where, Grid *grid)
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.
void Roe_DivBSource(const State_1D *state, int is, int ie, Grid *grid)
#define ARRAY_1D(nx, type)
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
void BACKGROUND_FIELD(real x1, real x2, real x3, real *B0)
double ** vL
Primitive variables to the left of the interface, .
double * A
Right interface area, A[i] = .