PLUTO
Macros | Functions
hllc.c File Reference

Implement the HLLC Riemann solver for relativistic MHD. More...

#include "pluto.h"
Include dependency graph for hllc.c:

Go to the source code of this file.

Macros

#define BX_MIN   1.e-6
 

Functions

void HLLC_Solver (const State_1D *state, int beg, int end, double *cmax, Grid *grid)
 

Detailed Description

Implement the HLLC Riemann solver for relativistic MHD.

Solve the Riemann problem for the relativistic MHD (RMHD) equations using the HLLC solver of Mignone & Bodo (2006).

On input, it takes left and right primitive state vectors state->vL and state->vR at zone edge i+1/2; On output, return flux and pressure vectors at the same interface i+1/2 (note that the i refers to i+1/2).

Also during this step, compute maximum wave propagation speed (cmax) for explicit time step computation.

Reference:

  • "An HLLC Riemann solver for relativistic flows - II. Magnetohydrodynamics", Mignone and Bodo, MNRAS (2006) 368, 1040.
Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
June 10, 2015

Definition in file hllc.c.

Macro Definition Documentation

#define BX_MIN   1.e-6

Definition at line 26 of file hllc.c.

Function Documentation

void HLLC_Solver ( const State_1D *  state,
int  beg,
int  end,
double *  cmax,
Grid *  grid 
)

Solve the RMHD Riemann problem using the HLLC Riemann solver.

Parameters
[in,out]statepointer to State_1D structure
[in]beginitial grid index
[out]endfinal grid index
[out]cmax1D array of maximum characteristic speeds
[in]gridpointer to array of Grid structures.

Definition at line 29 of file hllc.c.

41 {
42  int nv, i;
43  double scrh;
44  double usL[NFLX], usR[NFLX];
45  double Bx, Bys, Bzs;
46  double BtFBt, Bt2, FBt2;
47  double a, b, c;
48  double ps, vxs, vys, vzs, gammas_2, vBs;
49  double vxl, vxr, alpha_l, alpha_r;
50 
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;
56 
57  if (fL == NULL){
58  fL = ARRAY_2D(NMAX_POINT, NFLX, double);
59  fR = ARRAY_2D(NMAX_POINT, NFLX, double);
60 
61  Uhll = ARRAY_2D(NMAX_POINT, NVAR, double);
62  Fhll = ARRAY_2D(NMAX_POINT, NVAR, double);
63 
64  pL = ARRAY_1D(NMAX_POINT, double);
65  pR = ARRAY_1D(NMAX_POINT, double);
66  #ifdef GLM_MHD
67  VL = ARRAY_2D(NMAX_POINT, NVAR, double);
68  VR = ARRAY_2D(NMAX_POINT, NVAR, double);
69  UL = ARRAY_2D(NMAX_POINT, NVAR, double);
70  UR = ARRAY_2D(NMAX_POINT, NVAR, double);
71  #endif
72 
73  a2L = ARRAY_1D(NMAX_POINT, double);
74  a2R = ARRAY_1D(NMAX_POINT, double);
75  hL = ARRAY_1D(NMAX_POINT, double);
76  hR = ARRAY_1D(NMAX_POINT, double);
77  }
78 
79  #ifdef GLM_MHD
80  GLM_Solve (state, VL, VR, beg, end, grid);
81  PrimToCons (VL, UL, beg, end);
82  PrimToCons (VR, UR, beg, end);
83  #else
84  VL = state->vL; UL = state->uL;
85  VR = state->vR; UR = state->uR;
86  #endif
87 
88 /* ----------------------------------------------------
89  compute sound speed & fluxes at zone interfaces
90  ---------------------------------------------------- */
91 
92  SoundSpeed2 (VL, a2L, hL, beg, end, FACE_CENTER, grid);
93  SoundSpeed2 (VR, a2R, hR, beg, end, FACE_CENTER, grid);
94 
95  Flux (UL, VL, hL, fL, pL, beg, end);
96  Flux (UR, VR, hR, fR, pR, beg, end);
97 
98 /* ----------------------------------------
99  get max and min signal velocities
100  ---------------------------------------- */
101 
102  SL = state->SL; SR = state->SR;
103  HLL_Speed (VL, VR, a2L, a2R, hL, hR, SL, SR, beg, end);
104 
105 /* ----------------------------------------------
106  compute HLL state and flux
107  ---------------------------------------------- */
108 
109  for (i = beg; i <= end; i++) {
110 
111  scrh = MAX(fabs(SL[i]), fabs(SR[i]));
112  cmax[i] = scrh;
113 
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];
119  Uhll[i][nv] *= scrh;
120 
121  Fhll[i][nv] = SL[i]*SR[i]*(uR[nv] - uL[nv])
122  + SR[i]*fL[i][nv] - SL[i]*fR[i][nv];
123  Fhll[i][nv] *= scrh;
124  }
125  Uhll[i][MXn] += (pL[i] - pR[i])*scrh;
126  Fhll[i][MXn] += (SR[i]*pL[i] - SL[i]*pR[i])*scrh;
127  }
128 
129 /* ----------------------------------------------
130  Solve Riemann problem
131  ---------------------------------------------- */
132 
133  for (i = beg; i <= end; i++) {
134 
135  if (SL[i] >= 0.0){
136  for (nv = NFLX; nv--; ){
137  state->flux[i][nv] = fL[i][nv];
138  }
139  state->press[i] = pL[i];
140  }else if (SR[i] <= 0.0){
141  for (nv = NFLX; nv--; ){
142  state->flux[i][nv] = fR[i][nv];
143  }
144  state->press[i] = pR[i];
145  }else{
146 
147  uL = UL[i]; uR = UR[i];
148 
149 #if SHOCK_FLATTENING == MULTID
150  if ((state->flag[i] & FLAG_HLL) || (state->flag[i+1] & FLAG_HLL)){
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];
155  state->flux[i][nv] *= scrh;
156  }
157  state->press[i] = (SR[i]*pL[i] - SL[i]*pR[i])*scrh;
158  continue;
159  }
160 #endif
161 
162  vxl = VL[i][VXn];
163  vxr = VR[i][VXn];
164 
165  EXPAND(Bx = Uhll[i][BXn]; ,
166  Bys = Uhll[i][BXt]; ,
167  Bzs = Uhll[i][BXb];)
168 
169  #if RMHD_REDUCED_ENERGY == YES
170  Uhll[i][ENG] += Uhll[i][RHO];
171  Fhll[i][ENG] += Fhll[i][RHO];
172  #endif
173 
174  if (fabs(Bx) < BX_MIN) {
175 
176  BtFBt = Bt2 = FBt2 = 0.0;
177  a = Fhll[i][ENG];
178  b = - (Fhll[i][MXn] + Uhll[i][ENG]);
179  c = Uhll[i][MXn];
180 
181  }else{
182 
183  BtFBt = EXPAND(0.0, + Uhll[i][BXt]*Fhll[i][BXt],
184  + Uhll[i][BXb]*Fhll[i][BXb]);
185 
186  Bt2 = EXPAND(0.0, + Uhll[i][BXt]*Uhll[i][BXt],
187  + Uhll[i][BXb]*Uhll[i][BXb]);
188 
189  FBt2 = EXPAND(0.0, + Fhll[i][BXt]*Fhll[i][BXt],
190  + Fhll[i][BXb]*Fhll[i][BXb]);
191 
192  a = Fhll[i][ENG] - BtFBt;
193  b = Bt2 + FBt2 - (Fhll[i][MXn] + Uhll[i][ENG]);
194  c = Uhll[i][MXn] - BtFBt;
195  }
196 
197  /* -------------------------
198  solve quadratic, get v*
199  ------------------------- */
200 /*
201  if (fabs(a) > 1.e-12){
202  vxs = 0.5*(- b - sqrt(b*b - 4.0*a*c))/a;
203  }else{
204  vxs = -c/b;
205  }
206 */
207 
208 /*
209  scrh = -0.5*(b + DSIGN(b)*sqrt(b*b - 4.0*a*c));
210  vxs = c/scrh;
211 */
212 
213  scrh = 1.0 + sqrt(1.0 - 4.0*a*c/(b*b));
214  vxs = - 2.0*c/(b*scrh);
215 
216  if (fabs(Bx) < BX_MIN) {
217 
218  /* -------------------------------
219  the value of vy and vz
220  is irrelevant in this case
221  ------------------------------- */
222 
223  ps = Fhll[i][MXn] - Fhll[i][ENG]*vxs;
224 
225  alpha_l = (SL[i] - vxl)/(SL[i] - vxs);
226  alpha_r = (SR[i] - vxr)/(SR[i] - vxs);
227 
228  usL[RHO] = uL[RHO]*alpha_l;
229  usR[RHO] = uR[RHO]*alpha_r;
230 
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);
233 
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; )
240 
241  #if RMHD_REDUCED_ENERGY == YES
242  usL[MXn] += usL[RHO]*vxs;
243  usR[MXn] += usR[RHO]*vxs;
244  #endif
245 
246  EXPAND( usL[BXn] = Bx;
247  usR[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; )
252 
253  }else{
254 
255  EXPAND( ,
256  vys = (Bys*vxs - Fhll[i][BXt])/Bx; ,
257  vzs = (Bzs*vxs - Fhll[i][BXb])/Bx; )
258 
259  ps = (BtFBt - Fhll[i][ENG])*vxs
260  + Bx*Bx + Fhll[i][MXn] - FBt2;
261 
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);
265 
266  alpha_l = (SL[i] - vxl)/(SL[i] - vxs);
267  alpha_r = (SR[i] - vxr)/(SR[i] - vxs);
268 
269  usL[RHO] = uL[RHO]*alpha_l;
270  usR[RHO] = uR[RHO]*alpha_r;
271 
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);
274 
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); )
281 
282  #if RMHD_REDUCED_ENERGY == YES
283  usL[MXn] += usL[RHO]*vxs;
284  usR[MXn] += usR[RHO]*vxs;
285  #endif
286 
287  EXPAND( usL[BXn] = usR[BXn] = Bx; ,
288  usL[BXt] = usR[BXt] = Bys; ,
289  usL[BXb] = usR[BXb] = Bzs; )
290  }
291 
292  #ifdef GLM_MHD
293  usL[PSI_GLM] = usR[PSI_GLM] = uL[PSI_GLM];
294  #endif
295 
296 /* -------------------------------------------
297  Check consistency condition
298  ------------------------------------------- */
299 /*
300  for (nv = 0; nv < NVAR; nv++) {
301  scrh = (vxs - SL[i])*usL[nv] + (SR[i] - vxs)*usR[nv];
302  scrh -= SR[i]*uR[i][nv] - SL[i]*uL[i][nv] +
303  fL[i][nv] - fR[i][nv];
304  if (fabs(scrh/(SR[i]-SL[i])) > 1.e-5){
305  printf (" Consistency condition violated\n");
306  printf ("%d %d %12.6e\n",i,nv, scrh);
307  }
308  }
309 */
310 /* ---- Compute HLLC flux ---- */
311 
312  if (vxs > 0.0) {
313  for (nv = NFLX; nv--; ) {
314  state->flux[i][nv] = fL[i][nv] + SL[i]*(usL[nv] - uL[nv]);
315  }
316  state->press[i] = pL[i];
317  }else {
318  for (nv = NFLX; nv--; ) {
319  state->flux[i][nv] = fR[i][nv] + SR[i]*(usR[nv] - uR[nv]);
320  }
321  state->press[i] = pR[i];
322  }
323  }
324 
325  } /* -- end loop on points -- */
326 
327 /* --------------------------------------------------------
328  initialize souRce term
329  -------------------------------------------------------- */
330 
331  #if DIVB_CONTROL == EIGHT_WAVES
332 /*
333  POWELL_DIVB_SOURCE (state, beg, end, grid);
334 */
335  /* ----------------------------------------------------
336  to avoid conversion problems in HLL_DIVB_SOURCE,
337  we use the HLL average provided by SR = -SL = 1
338  ---------------------------------------------------- */
339 /*
340  for (i = beg; i <= end; i++) {
341  uL = state->uL[i]; uR = state->uR[i];
342  for (nv = 0; nv < NFLX; nv++) {
343  Uhll[i][nv] = 0.5*(uR[nv] + uL[nv] + fL[i][nv] - fR[i][nv]);
344  }
345  Uhll[i][MXn] += (pL[i] - pR[i])*0.5;
346  for (nv = NFLX; nv < NVAR; nv++) Uhll[i][nv] = 0.0;
347  }
348 */
349  HLL_DIVB_SOURCE (state, Uhll, beg+1, end, grid);
350  #endif
351 
352 }
RMHD_REDUCED_ENERGY
#define RMHD_REDUCED_ENERGY
By turning RMHD_REDUCED_ENERGY to YES, we let PLUTO evolve the total energy minus the mass density co...
Definition: mod_defs.h:193
a
static double a
Definition: init.c:135
MAX
#define MAX(a, b)
Definition: macros.h:101
GLM_Solve
void GLM_Solve(const State_1D *state, double **VL, double **VR, int beg, int end, Grid *grid)
Definition: glm.c:24
STATE_1D::flux
double ** flux
upwind flux computed with the Riemann solver
Definition: structs.h:149
Flux
void Flux(double **u, double **w, double *a2, double **fx, double *p, int beg, int end)
Definition: fluxes.c:23
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
YES
#define YES
Definition: pluto.h:25
STATE_1D::SR
double * SR
Rightmost velocity in the Riemann fan at i+1/2.
Definition: structs.h:167
PSI_GLM
#define PSI_GLM
Definition: mod_defs.h:34
STATE_1D::vR
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
SoundSpeed2
void SoundSpeed2(double **v, double *cs2, double *h, int beg, int end, int pos, Grid *grid)
Definition: eos.c:16
HLL_DIVB_SOURCE
void HLL_DIVB_SOURCE(const State_1D *, double **, int, int, Grid *)
Definition: source.c:53
BXn
int BXn
Definition: globals.h:75
MXt
int MXt
Definition: globals.h:74
NFLX
#define NFLX
Definition: mod_defs.h:32
FACE_CENTER
#define FACE_CENTER
Definition: pluto.h:206
Bx
static double Bx
Definition: hlld.c:62
STATE_1D::SL
double * SL
Leftmost velocity in the Riemann fan at i+1/2.
Definition: structs.h:166
STATE_1D::flag
unsigned char * flag
Definition: structs.h:168
MXn
int MXn
Definition: globals.h:74
menu.c
tuple c
Definition: menu.py:375
VXn
int VXn
Definition: globals.h:73
BX_MIN
#define BX_MIN
Definition: hllc.c:26
STATE_1D::uR
double ** uR
same as vR, in conservative vars
Definition: structs.h:145
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
PrimToCons
void PrimToCons(double **uprim, double **ucons, int ibeg, int iend)
Definition: mappers.c:26
i
int i
Definition: analysis.c:2
HLL_Speed
void HLL_Speed(double **vL, double **vR, double *a2L, double *a2R, double *SL, double *SR, int beg, int end)
Definition: hll_speed.c:24
STATE_1D::vL
double ** vL
Primitive variables to the left of the interface, .
Definition: structs.h:136
STATE_1D::press
double * press
Upwind pressure term computed with the Riemann solver.
Definition: structs.h:164
FLAG_HLL
#define FLAG_HLL
Use HLL Riemann solver.
Definition: pluto.h:181
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
BXt
int BXt
Definition: globals.h:75
NVAR
#define NVAR
Definition: pluto.h:609
BXb
int BXb
Definition: globals.h:75
GLM_MHD
#define GLM_MHD
Definition: glm.h:25
MXb
int MXb
Definition: globals.h:74
STATE_1D::uL
double ** uL
same as vL, in conservative vars
Definition: structs.h:144

Here is the call graph for this function: