37 #define COUNT_FAILURES NO
69 double,
double *,
double *);
74 double *cmax,
Grid *grid)
89 static double **fluxL, **fluxR;
90 static double *pL, *pR, *a2L, *a2R, *hL, *hR;
91 static double **Uhll, **Fhll, **Vhll;
92 double *vL, *vR, *fL, *fR, *uL, *uR, *SL, *SR;
93 static double **VL, **VR, **UL, **UR;
95 double p0, f0, p, f, dp, dS_1;
105 #if COUNT_FAILURES == YES
106 static double totfail=0.0,totzones=1.0;
107 static int last_step = -1;
111 fp = fopen(
"hlld_fails.dat",
"w");
112 fprintf (fp,
"# Step Failures (x 100) \n");
113 fprintf (fp,
"# --------------------------\n");
118 MPI_Allreduce (&totfail, &scrh, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
120 MPI_Allreduce (&totzones, &scrh, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
124 fp = fopen(
"hlld_fails.dat",
"a+");
125 fprintf (fp,
" %d %f\n",
g_stepNumber, totfail/totzones*100.);
128 totfail = totzones = 0.0;
164 GLM_Solve (state, VL, VR, beg, end, grid);
168 VL = state->
vL; UL = state->
uL;
169 VR = state->
vR; UR = state->
uR;
179 Flux (UL, VL, hL, fluxL, pL, beg, end);
180 Flux (UR, VR, hR, fluxR, pR, beg, end);
184 SL = state->
SL; SR = state->
SR;
185 HLL_Speed (VL, VR, a2L, a2R, hL, hR, SL, SR, beg, end);
191 for (i = beg; i <= end; i++) {
194 if (!(grid[
IDIR].x[i] < 0.5 && grid[
IDIR].x[i+1] > 0.5))
continue;
197 #if COUNT_FAILURES == YES
201 scrh =
MAX(fabs(SL[i]), fabs(SR[i]));
204 vL = VL[
i]; vR = VR[
i]; fR = fluxR[
i];
205 uL = UL[
i]; uR = UR[
i]; fL = fluxL[
i];
209 for (nv =
NFLX; nv--; ) state->
flux[i][nv] = fL[nv];
210 state->
press[i] = pL[i];
212 }
else if (SR[i] <= 0.0){
214 for (nv =
NFLX; nv--; ) state->
flux[i][nv] = fR[nv];
215 state->
press[i] = pR[i];
221 dS_1 = 1.0/(SR[
i] - SL[
i]);
222 for (nv =
NFLX; nv--; ){
224 Uhll[
i][nv] = SR[
i]*uR[nv] - SL[
i]*uL[nv] + fL[nv] - fR[nv];
227 Fhll[
i][nv] = SR[
i]*fL[nv] - SL[
i]*fR[nv] + SL[
i]*SR[
i]*(uR[nv] - uL[nv]);
230 Uhll[
i][
MXn] += (pL[
i] - pR[
i])*dS_1;
233 double vxR = vR[
VXn];
234 double vxL = vL[
VXn];
235 Uhll[
i][nv] = (SR[
i] - vxR)*uR[nv] - (SL[i] - vxL)*uL[nv];
238 Fhll[
i][nv] = SR[
i]*uL[nv]*vL[
VXn] - SL[
i]*uR[nv]*vR[
VXn]
239 + SL[
i]*SR[
i]*(uR[nv] - uL[nv]);
246 #if SHOCK_FLATTENING == MULTID
248 for (nv =
NFLX; nv--; ) state->
flux[i][nv] = Fhll[i][nv];
249 state->
press[i] = (SR[i]*pL[i] - SL[i]*pR[i])*dS_1;
259 for (nv = 0; nv <
NFLX; nv++){
260 PaL.
R[nv] = SL[
i]*uL[nv] - fL[nv];
261 PaR.
R[nv] = SR[
i]*uR[nv] - fR[nv];
265 #if RMHD_REDUCED_ENERGY == YES
266 PaL.
R[ENG] += PaL.
R[
RHO];
267 PaR.
R[ENG] += PaR.
R[
RHO];
272 scrh =
MAX(pL[i], pR[i]);
273 if (
Bx*
Bx/scrh < 0.01) {
277 b = PaR.
R[ENG] - PaL.
R[ENG] + SR[
i]*PaL.
R[
MXn] - SL[
i]*PaR.
R[
MXn];
278 c = PaL.
R[
MXn]*PaR.
R[ENG] - PaR.
R[
MXn]*PaL.
R[ENG];
279 scrh = b*b - 4.0*a*
c;
280 scrh =
MAX(scrh,0.0);
281 p0 = 0.5*(- b + sqrt(scrh))*dS_1;
294 if (f0 != f0 || PaL.
fail) switch_to_hll = 1;
299 if (fabs(f0) > 1.e-12 && !switch_to_hll){
300 p = 1.025*p0; f = f0;
304 printf (
"k = %d, p = %12.6e, f = %12.6e\n",k,p,f);
308 if (f != f || PaL.
fail || (k > 7) || (fabs(f) > fabs(f0) && k > 4)) {
313 dp = (p - p0)/(f - f0)*f;
317 if (p < 0.0) p = 1.e-6;
318 if (fabs(dp) < 1.e-5*p || fabs(f) < 1.e-6)
break;
328 if (PaL.
fail) switch_to_hll = 1;
331 #if COUNT_FAILURES == YES
335 for (nv = NFLX; nv--; ) state->
flux[i][nv] = Fhll[i][nv];
336 state->
press[i] = (SR[i]*pL[i] - SL[i]*pR[i])*dS_1;
348 if (PaL.
Sa >= -1.e-6){
351 #if RMHD_REDUCED_ENERGY == YES
352 PaL.
u[ENG] -= PaL.
u[
RHO];
354 for (nv = NFLX; nv--; ) {
355 state->
flux[
i][nv] = fL[nv] + SL[
i]*(PaL.
u[nv] - uL[nv]);
359 }
else if (PaR.
Sa <= 1.e-6){
362 #if RMHD_REDUCED_ENERGY == YES
363 PaR.
u[ENG] -= PaR.
u[
RHO];
365 for (nv = NFLX; nv--; ) {
366 state->
flux[
i][nv] = fR[nv] + SR[
i]*(PaR.
u[nv] - uR[nv]);
374 #if RMHD_REDUCED_ENERGY == YES
375 PaL.
u[ENG] -= PaL.
u[
RHO];
378 for (nv = NFLX; nv--; ) {
379 state->
flux[
i][nv] = fL[nv] + SL[
i]*(PaL.
u[nv] - uL[nv])
380 + PaL.
Sa*(Uc[nv] - PaL.
u[nv]);
385 #if RMHD_REDUCED_ENERGY == YES
386 PaR.
u[ENG] -= PaR.
u[
RHO];
389 for (nv = NFLX; nv--; ) {
390 state->
flux[
i][nv] = fR[nv] + SR[
i]*(PaR.
u[nv] - uR[nv])
391 + PaR.
Sa*(Uc[nv] - PaR.
u[nv]);
403 #if DIVB_CONTROL == EIGHT_WAVES
409 #undef COUNT_FAILURES
420 double dK, Bxc, Byc, Bzc;
425 sBx =
Bx > 0.0 ? 1.0:-1.0;
432 dK = PaR->
Kx - PaL->
Kx + 1.e-12;
434 EXPAND(Bxc =
Bx*dK; ,
436 Byc = PaR->
By*(PaR->
Kx - PaR->
vx)
437 - PaL->
By*(PaL->
Kx - PaL->
vx)
438 +
Bx*(PaR->
vy - PaL->
vy); ,
440 Bzc = PaR->
Bz*(PaR->
Kx - PaR->
vx)
441 - PaL->
Bz*(PaL->
Kx - PaL->
vx)
442 +
Bx*(PaR->
vz - PaL->
vz);)
444 KLBc = EXPAND(PaL->
Kx*Bxc, + PaL->
Ky*Byc, + PaL->
Kz*Bzc);
445 KRBc = EXPAND(PaR->
Kx*Bxc, + PaR->
Ky*Byc, + PaR->
Kz*Bzc);
447 vxcL = PaL->
Kx - dK*
Bx*(1.0 - PaL->
K2)/(PaL->
sw*dK - KLBc);
448 vxcR = PaR->
Kx - dK*Bx*(1.0 - PaR->
K2)/(PaR->
sw*dK - KRBc);
452 Sc = 0.5*(vxcL + vxcR);
462 success = (vxcL - PaL->
Kx) > -1.e-6;
463 success *= (PaR->
Kx - vxcR) > -1.e-6;
465 success *= (PaL->
S - PaL->
vx) < 0.0;
466 success *= (PaR->
S - PaR->
vx) > 0.0;
468 success *= (PaR->
w - p) > 0.0;
469 success *= (PaL->
w - p) > 0.0;
470 success *= (PaL->
Sa - PaL->
S) > -1.e-6;
471 success *= (PaR->
S - PaR->
Sa) > -1.e-6;
473 PaL->
fail = !success;
499 double A, C, G,
X,
s;
500 double vx, vy, vz,
scrh, S, *R;
505 A = R[
MXn] + p*(1.0 - S*S) - S*R[ENG];
507 C = EXPAND(0.0, + R[BXt]*R[
MXt], + R[BXb]*R[
MXb]);
508 X =
Bx*(A*S*
Bx + C) - (A + G)*(S*p + R[ENG]);
512 EXPAND(vx = (
Bx*(A*
Bx + C*S) - (R[
MXn] + p)*(G + A) ); ,
514 vy = ( - (A + G -
Bx*
Bx*(1.0 - S*S))*R[
MXt]
515 + R[
BXt]*(C +
Bx*(S*R[
MXn] - R[ENG])) ); ,
517 vz = ( - (A + G -
Bx*
Bx*(1.0 - S*S))*R[
MXb]
518 + R[
BXb]*(C +
Bx*(S*R[
MXn] - R[ENG])) );)
520 scrh = EXPAND(vx*R[
MXn], + vy*R[MXt], + vz*R[MXb]);
521 scrh = X*R[ENG] -
scrh;
522 Pv->
w = p + scrh/(X*S - vx);
524 if (Pv->
w < 0.0)
return(0);
526 EXPAND(Pv->
vx = vx/X; ,
556 EXPAND(Pv->
Bx =
Bx; ,
557 Pv->
By = -(R[BXt]*(S*p + R[ENG]) -
Bx*R[MXt])/A; ,
558 Pv->
Bz = -(R[BXb]*(S*p + R[ENG]) -
Bx*R[MXb])/A;)
560 s =
Bx > 0.0 ? 1.0:-1.0;
561 if (side < 0) s *= -1.0;
563 Pv->
sw = s*sqrt(Pv->
w);
565 scrh = 1.0/(S*p + R[ENG] +
Bx*Pv->
sw);
566 EXPAND(Pv->
Kx = scrh*(R[MXn] + p + R[
BXn]*Pv->
sw); ,
567 Pv->
Ky = scrh*(R[MXt] + R[BXt]*Pv->
sw); ,
568 Pv->
Kz = scrh*(R[MXb] + R[BXb]*Pv->
sw);)
580 double vB, *ua, *R, S;
587 scrh = 1.0/(S - Pa->
vx);
590 EXPAND(ua[
BXn] =
Bx; ,
595 ua[ENG] = (R[ENG] + p*Pa->
vx - vB*
Bx)*scrh;
597 EXPAND(ua[
MXn] = (ua[ENG] + p)*Pa->
vx - vB*ua[BXn]; ,
598 ua[
MXt] = (ua[ENG] + p)*Pa->
vy - vB*ua[BXt]; ,
599 ua[
MXb] = (ua[ENG] + p)*Pa->
vz - vB*ua[BXb];)
611 double vxcL, vycL, vzcL, KLBc;
612 double vxcR, vycR, vzcR, KRBc;
613 double vxc, vyc, vzc, vBc;
614 double Bxc, Byc, Bzc, Sa, vxa;
619 dK = (PaR->
Kx - PaL->
Kx) + 1.e-12;
621 EXPAND(Bxc =
Bx*dK; ,
623 Byc = PaR->
By*(PaR->
Kx - PaR->
vx)
624 - PaL->
By*(PaL->
Kx - PaL->
vx)
625 +
Bx*(PaR->
vy - PaL->
vy); ,
627 Bzc = PaR->
Bz*(PaR->
Kx - PaR->
vx)
628 - PaL->
Bz*(PaL->
Kx - PaL->
vx)
629 +
Bx*(PaR->
vz - PaL->
vz);)
635 EXPAND(Uc[
BXn] = Bxc; ,
639 KLBc = EXPAND(PaL->
Kx*Bxc, + PaL->
Ky*Byc, + PaL->
Kz*Bzc);
640 KRBc = EXPAND(PaR->
Kx*Bxc, + PaR->
Ky*Byc, + PaR->
Kz*Bzc);
642 scrhL = (1.0 - PaL->
K2)/(PaL->
sw - KLBc);
643 scrhR = (1.0 - PaR->
K2)/(PaR->
sw - KRBc);
645 EXPAND(vxcL = PaL->
Kx - Uc[
BXn]*scrhL;
646 vxcR = PaR->
Kx - Uc[
BXn]*scrhR; ,
648 vycL = PaL->
Ky - Uc[
BXt]*scrhL;
649 vycR = PaR->
Ky - Uc[
BXt]*scrhR; ,
651 vzcL = PaL->
Kz - Uc[
BXb]*scrhL;
652 vzcR = PaR->
Kz - Uc[
BXb]*scrhR;)
654 EXPAND(vxc = 0.5*(vxcL + vxcR); ,
655 vyc = 0.5*(vycL + vycR); ,
656 vzc = 0.5*(vzcL + vzcR);)
670 vBc = EXPAND(vxc*Uc[
BXn], + vyc*Uc[
BXt], + vzc*Uc[
BXb]);
672 Uc[
RHO] = ua[
RHO]*(Sa - vxa)/(Sa - vxc);
673 Uc[ENG] = (Sa*ua[ENG] - ua[
MXn] + p*vxc - vBc*
Bx)/(Sa - vxc);
675 EXPAND(Uc[
MXn] = (Uc[ENG] + p)*vxc - vBc*
Bx; ,
676 Uc[
MXt] = (Uc[ENG] + p)*vyc - vBc*Uc[BXt]; ,
677 Uc[
MXb] = (Uc[ENG] + p)*vzc - vBc*Uc[BXb];)
687 double vel2, Bmag2, vB, lor2;
692 vB = EXPAND(v[VX1]*v[BX1], + v[VX2]*v[BX2], + v[VX3]*v[BX3]);
693 pt = v[PRS] + 0.5*(Bmag2*(1.0 -
vel2) + vB*vB);
701 printf (
"RHO_LEFT %18.9e\n",vL[
RHO]);
702 printf (
"VX_LEFT %18.9e\n",vL[
VXn]);
703 printf (
"VY_LEFT %18.9e\n",vL[
VXt]);
705 printf (
"VZ_LEFT %18.9e \n" ,vL[
VXb]);
707 printf (
"VZ_LEFT %18.9e \n" ,0.0);
709 printf (
"BY_LEFT %18.9e\n",vL[
BXt]);
711 printf (
"BZ_LEFT %18.9e \n" ,vL[
BXb]);
713 printf (
"BZ_LEFT %18.9e \n" ,0.0);
715 printf (
"PR_LEFT %18.9e \n" ,vL[PRS]);
717 printf (
"RHO_RIGHT %18.9e\n",vR[
RHO]);
718 printf (
"VX_RIGHT %18.9e\n",vR[
VXn]);
719 printf (
"VY_RIGHT %18.9e\n",vR[
VXt]);
721 printf (
"VZ_RIGHT %18.9e \n", vR[
VXb]);
723 printf (
"VZ_RIGHT %18.9e \n", 0.0);
725 printf (
"BY_RIGHT %18.9e\n", vR[
BXt]);
727 printf (
"BZ_RIGHT %18.9e\n", vR[
BXb]);
729 printf (
"BZ_RIGHT %18.9e\n",0.0);
731 printf (
"PR_RIGHT %18.9e \n",vR[PRS]);
732 printf (
"BX_CONST %18.9e \n",vR[
BXn]);
737 double pguess,
double *vL,
double *vR)
747 printf (
"pguess = %f, F = %f\n",pguess,
HLLD_Fstar(PaL, PaR, pguess));
748 printf (
"S = %f %f %f %f %f\n", PaL->
S, PaL->
Sa,
Sc, PaR->
Sa, PaR->
S);
750 fp = fopen(
"new.dat",
"w");
751 for (p = 0.01*pguess; p <= 10.0*pguess; p += 1.e-4*pguess){
753 fprintf (fp,
"%f %f %f %f\n", p, f,
Sc - PaL->
Sa, PaR->
Sa -
Sc);
struct RIEMANN_STATE Riemann_State
void GLM_Solve(const State_1D *state, double **VL, double **VR, int beg, int end, Grid *grid)
static void HLLD_GetCState(Riemann_State *, Riemann_State *, double, double *)
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)
void HLLD_Solver(const State_1D *state, int beg, int end, double *cmax, Grid *grid)
double * SR
Rightmost velocity in the Riemann fan at i+1/2.
int g_maxRiemannIter
Maximum number of iterations for iterative Riemann Solver.
double ** vR
Primitive variables to the right of the interface, .
static double HLLD_TotalPressure(double *)
void SoundSpeed2(double **v, double *cs2, double *h, int beg, int end, int pos, Grid *grid)
static void HLLD_PrintStates(double *, double *)
void HLL_DIVB_SOURCE(const State_1D *, double **, int, int, Grid *)
static int HLLD_GetRiemannState(Riemann_State *, double, int)
long int g_stepNumber
Gives the current integration step number.
int ConsToPrim(double **ucons, double **uprim, int ibeg, int iend, unsigned char *flag)
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)
static void HLLD_GetAState(Riemann_State *, double p)
static void HLLD_PrintWhatsWrong(Riemann_State *, Riemann_State *, double, double *, double *)
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)
static double HLLD_Fstar(Riemann_State *, Riemann_State *, double)
double ** uL
same as vL, in conservative vars