Solve the Riemann problem for the relativistic MHD equations using the dominant-wave Riemann solver.
17 double *uL, *uR, *flux, cL, cR, lambda;
18 static double **fL, *pL, *a2L, *hL, *cmin_L, *cmax_L;
19 static double **fR, *pR, *a2R, *hR, *cmin_R, *cmax_R;
20 static double **VL, **VR, **UL, **UR;
21 double num, den, dU, dF, vn;
50 GLM_Solve (state, VL, VR, beg, end, grid);
54 VL = state->
vL; UL = state->
uL;
55 VR = state->
vR; UR = state->
uR;
65 Flux (UL, VL, hL, fL, pL, beg, end);
66 Flux (UR, VR, hR, fR, pR, beg, end);
75 for (i = beg; i <= end; i++) {
79 flux = state->
flux[
i];
83 cL =
MAX(fabs(cmax_L[i]), fabs(cmin_L[i]));
84 cR =
MAX(fabs(cmax_R[i]), fabs(cmin_R[i]));
85 cmax[
i] =
MAX(cL, cR);
86 state->
SL[
i] = -cmax[
i];
87 state->
SR[
i] = cmax[
i];
91 for (nv = 0; nv <
NFLX; nv++) {
93 dF = fR[
i][nv] - fL[
i][nv] + (nv ==
MXn ? (pR[
i] - pL[
i]):0.0);
98 lambda = fabs(num)/(fabs(den) + 1.e-12);
100 if (lambda > cmax[i]) lambda = cmax[
i];
102 vn = 0.5*( fabs(state->
vR[i][
VXn]) + fabs(state->
vL[i][
VXn]) );
103 vn =
MAX( fabs(state->
vR[i][
VXn]), fabs(state->
vL[i][VXn]) );
104 if (lambda < vn) lambda = vn;
113 for (nv =
NFLX; nv--; ) {
114 flux[nv] = 0.5*(fL[
i][nv] + fR[
i][nv] - lambda*(uR[nv] - uL[nv]));
116 state->
press[
i] = 0.5*(pL[
i] + pR[
i]);
123 #if DIVB_CONTROL == EIGHT_WAVES
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 SoundSpeed2(double **v, double *cs2, double *h, int beg, int end, int pos, Grid *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)
double ** vL
Primitive variables to the left of the interface, .
double * press
Upwind pressure term computed with the Riemann solver.
void MaxSignalSpeed(double **v, double *cs2, double *cmin, double *cmax, int beg, int end)
#define ARRAY_2D(nx, ny, type)
void POWELL_DIVB_SOURCE(const State_1D *, int, int, Grid *)
double ** uL
same as vL, in conservative vars