30 double *cmax,
Grid *grid)
46 double BtFBt, Bt2, FBt2;
48 double ps, vxs, vys, vzs, gammas_2, vBs;
49 double vxl, vxr, alpha_l, alpha_r;
51 double *uL, *uR, *SL, *SR;
52 static double **fL, **fR;
53 static double *pR, *pL, *a2L, *a2R, *hL, *hR;
54 static double **Uhll, **Fhll;
55 static double **VL, **VR, **UL, **UR;
80 GLM_Solve (state, VL, VR, beg, end, grid);
84 VL = state->
vL; UL = state->
uL;
85 VR = state->
vR; UR = state->
uR;
95 Flux (UL, VL, hL, fL, pL, beg, end);
96 Flux (UR, VR, hR, fR, pR, beg, end);
102 SL = state->
SL; SR = state->
SR;
103 HLL_Speed (VL, VR, a2L, a2R, hL, hR, SL, SR, beg, end);
109 for (i = beg; i <= end; i++) {
111 scrh =
MAX(fabs(SL[i]), fabs(SR[i]));
114 uL = UL[
i]; uR = UR[
i];
115 scrh = 1.0/(SR[
i] - SL[
i]);
116 for (nv =
NFLX; nv--; ){
117 Uhll[
i][nv] = SR[
i]*uR[nv] - SL[
i]*uL[nv]
118 + fL[
i][nv] - fR[
i][nv];
121 Fhll[
i][nv] = SL[
i]*SR[
i]*(uR[nv] - uL[nv])
122 + SR[i]*fL[i][nv] - SL[i]*fR[i][nv];
125 Uhll[
i][
MXn] += (pL[
i] - pR[
i])*scrh;
126 Fhll[
i][
MXn] += (SR[
i]*pL[
i] - SL[
i]*pR[
i])*scrh;
133 for (i = beg; i <= end; i++) {
136 for (nv =
NFLX; nv--; ){
137 state->
flux[
i][nv] = fL[
i][nv];
140 }
else if (SR[i] <= 0.0){
141 for (nv =
NFLX; nv--; ){
142 state->
flux[
i][nv] = fR[
i][nv];
147 uL = UL[
i]; uR = UR[
i];
149 #if SHOCK_FLATTENING == MULTID
151 scrh = 1.0/(SR[
i] - SL[
i]);
152 for (nv =
NFLX; nv--; ){
153 state->
flux[
i][nv] = SL[
i]*SR[
i]*(uR[nv] - uL[nv])
154 + SR[i]*fL[i][nv] - SL[i]*fR[i][nv];
157 state->
press[
i] = (SR[
i]*pL[
i] - SL[
i]*pR[
i])*scrh;
165 EXPAND(Bx = Uhll[i][
BXn]; ,
166 Bys = Uhll[
i][
BXt]; ,
170 Uhll[i][ENG] += Uhll[i][
RHO];
171 Fhll[
i][ENG] += Fhll[
i][
RHO];
176 BtFBt = Bt2 = FBt2 = 0.0;
178 b = - (Fhll[
i][
MXn] + Uhll[
i][ENG]);
183 BtFBt = EXPAND(0.0, + Uhll[i][
BXt]*Fhll[i][
BXt],
184 + Uhll[i][
BXb]*Fhll[i][
BXb]);
186 Bt2 = EXPAND(0.0, + Uhll[i][BXt]*Uhll[i][BXt],
187 + Uhll[i][BXb]*Uhll[i][BXb]);
189 FBt2 = EXPAND(0.0, + Fhll[i][BXt]*Fhll[i][BXt],
190 + Fhll[i][BXb]*Fhll[i][BXb]);
192 a = Fhll[
i][ENG] - BtFBt;
193 b = Bt2 + FBt2 - (Fhll[
i][
MXn] + Uhll[
i][ENG]);
194 c = Uhll[
i][
MXn] - BtFBt;
213 scrh = 1.0 + sqrt(1.0 - 4.0*a*c/(b*b));
214 vxs = - 2.0*c/(b*
scrh);
223 ps = Fhll[
i][
MXn] - Fhll[
i][ENG]*vxs;
225 alpha_l = (SL[
i] - vxl)/(SL[i] - vxs);
226 alpha_r = (SR[
i] - vxr)/(SR[i] - vxs);
228 usL[
RHO] = uL[
RHO]*alpha_l;
229 usR[
RHO] = uR[
RHO]*alpha_r;
231 usL[ENG] = (SL[
i]*uL[ENG] - fL[
i][ENG] + ps*vxs)/(SL[i] - vxs);
232 usR[ENG] = (SR[
i]*uR[ENG] - fR[
i][ENG] + ps*vxs)/(SR[i] - vxs);
234 EXPAND( usL[
MXn] = (usL[ENG] + ps)*vxs;
235 usR[
MXn] = (usR[ENG] + ps)*vxs; ,
236 usL[
MXt] = uL[
MXt]*alpha_l;
237 usR[
MXt] = uR[
MXt]*alpha_r; ,
238 usL[
MXb] = uL[
MXb]*alpha_l;
239 usR[
MXb] = uR[
MXb]*alpha_r; )
246 EXPAND( usL[
BXn] = Bx;
248 usL[
BXt] = uL[
BXt]*alpha_l;
249 usR[
BXt] = uR[
BXt]*alpha_r; ,
250 usL[
BXb] = uL[
BXb]*alpha_l;
251 usR[
BXb] = uR[
BXb]*alpha_r; )
256 vys = (Bys*vxs - Fhll[i][
BXt])/Bx; ,
257 vzs = (Bzs*vxs - Fhll[
i][
BXb])/Bx; )
259 ps = (BtFBt - Fhll[i][ENG])*vxs
260 + Bx*Bx + Fhll[
i][
MXn] - FBt2;
262 gammas_2 = EXPAND(vxs*vxs, + vys*vys, + vzs*vzs);
263 gammas_2 = 1.0 - gammas_2;
264 vBs = EXPAND(vxs*Bx, + vys*Bys, + vzs*Bzs);
266 alpha_l = (SL[
i] - vxl)/(SL[i] - vxs);
267 alpha_r = (SR[
i] - vxr)/(SR[i] - vxs);
269 usL[
RHO] = uL[
RHO]*alpha_l;
270 usR[
RHO] = uR[
RHO]*alpha_r;
272 usL[ENG] = (SL[
i]*uL[ENG] - fL[
i][ENG] + ps*vxs - vBs*
Bx)/(SL[i] - vxs);
273 usR[ENG] = (SR[
i]*uR[ENG] - fR[
i][ENG] + ps*vxs - vBs*
Bx)/(SR[i] - vxs);
275 EXPAND( usL[
MXn] = (usL[ENG] + ps)*vxs - vBs*Bx;
276 usR[
MXn] = (usR[ENG] + ps)*vxs - vBs*Bx; ,
277 usL[
MXt] = (SL[
i]*uL[
MXt] - fL[
i][
MXt] - Bx*(Bys*gammas_2 + vBs*vys))/(SL[
i] - vxs);
278 usR[
MXt] = (SR[
i]*uR[
MXt] - fR[
i][
MXt] - Bx*(Bys*gammas_2 + vBs*vys))/(SR[
i] - vxs); ,
279 usL[
MXb] = (SL[
i]*uL[
MXb] - fL[
i][
MXb] - Bx*(Bzs*gammas_2 + vBs*vzs))/(SL[
i] - vxs);
280 usR[
MXb] = (SR[
i]*uR[
MXb] - fR[
i][
MXb] - Bx*(Bzs*gammas_2 + vBs*vzs))/(SR[
i] - vxs); )
287 EXPAND( usL[
BXn] = usR[
BXn] = Bx; ,
288 usL[
BXt] = usR[
BXt] = Bys; ,
289 usL[
BXb] = usR[
BXb] = Bzs; )
313 for (nv =
NFLX; nv--; ) {
314 state->
flux[
i][nv] = fL[
i][nv] + SL[
i]*(usL[nv] - uL[nv]);
318 for (nv =
NFLX; nv--; ) {
319 state->
flux[
i][nv] = fR[
i][nv] + SR[
i]*(usR[nv] - uR[nv]);
331 #if DIVB_CONTROL == EIGHT_WAVES
#define RMHD_REDUCED_ENERGY
By turning RMHD_REDUCED_ENERGY to YES, we let PLUTO evolve the total energy minus the mass density co...
void GLM_Solve(const State_1D *state, double **VL, double **VR, int beg, int end, Grid *grid)
double ** flux
upwind flux computed with the Riemann solver
void Flux(double **u, double **w, double *a2, double **fx, double *p, int beg, int end)
double * SR
Rightmost velocity in the Riemann fan at i+1/2.
double ** vR
Primitive variables to the right of the interface, .
void HLLC_Solver(const State_1D *state, int beg, int end, real *cmax, Grid *grid)
void SoundSpeed2(double **v, double *cs2, double *h, int beg, int end, int pos, Grid *grid)
void HLL_DIVB_SOURCE(const State_1D *, double **, int, int, Grid *)
double * SL
Leftmost velocity in the Riemann fan at i+1/2.
double ** uR
same as vR, in conservative vars
#define ARRAY_1D(nx, type)
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
void PrimToCons(double **uprim, double **ucons, int ibeg, int iend)
void HLL_Speed(double **vL, double **vR, double *a2L, double *a2R, double *SL, double *SR, int beg, int end)
double ** vL
Primitive variables to the left of the interface, .
double * press
Upwind pressure term computed with the Riemann solver.
#define FLAG_HLL
Use HLL Riemann solver.
#define ARRAY_2D(nx, ny, type)
double ** uL
same as vL, in conservative vars