PLUTO
Macros | Functions
hlld.c File Reference

HLLD Riemann solver for the classical MHD equations. More...

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

Go to the source code of this file.

Macros

#define VERIFY_CONSISTENCY_CONDITION   NO
 

Functions

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

Detailed Description

HLLD Riemann solver for the classical MHD equations.

Solve the Riemann problem for the classical MHD equations with an ideal or isothermal equation of state. For the ideal case, we implement the four-state (of five-wave) HLLD solver of Miyoshi & Kusano (2005). For the isothermal case, we use the three-state approximation described by Mignone (2007).

The macro VERIFY_CONSISTENCY_CONDITION can be turned to YES to verify the correctness of the implementation.

On input, this function 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:

  • "A muLti-state HLL approximate Riemann Solver for ideal MHD", Miyoshi, T., Kusano, K., JCP 2005
  • "A simple and accurate Riemann solver for isothermal MHD", Mignone, JCP (2007) 225, 1427.
Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t) C. Zanni (zanni.nosp@m.@oat.nosp@m.o.ina.nosp@m.f.it)
Date
April 14, 2014

Definition in file hlld.c.

Macro Definition Documentation

#define VERIFY_CONSISTENCY_CONDITION   NO

Definition at line 37 of file hlld.c.

Function Documentation

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

Solve Riemann problem for the adiabatic MHD equations using the four-state HLLD Riemann solver of Miyoshi & Kusano (2005).

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.

Solve Riemann problem for the isothermal MHD equations using the three-state HLLD Riemann solver of Mignone (2007).

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 41 of file hlld.c.

54 {
55  int nv, i;
56  int revert_to_hllc;
57  double scrh, Uhll[NFLX];
58 
59  double usL[NFLX], ussl[NFLX];
60  double usR[NFLX], ussr[NFLX];
61 
62  double vsL, wsL, scrhL, S1L, sqrL, duL;
63  double vsR, wsR, scrhR, S1R, sqrR, duR;
64  double Bx, Bx1, SM, sBx, pts;
65  double vss, wss;
66  double *vL, *vR, *uL, *uR, *SL, *SR;
67  static double *ptL, *ptR;
68  static double **fL, **fR;
69  static double **VL, **VR, **UL, **UR;
70  static double *a2L, *a2R;
71  #if BACKGROUND_FIELD == YES
72  double B0n, B0t, B0b;
73  #endif
74  double **bgf;
75 
76  if (fL == NULL){
77  fL = ARRAY_2D(NMAX_POINT, NFLX, double);
78  fR = ARRAY_2D(NMAX_POINT, NFLX, double);
79 
80  ptR = ARRAY_1D(NMAX_POINT, double);
81  ptL = ARRAY_1D(NMAX_POINT, double);
82 
83  a2R = ARRAY_1D(NMAX_POINT, double);
84  a2L = ARRAY_1D(NMAX_POINT, double);
85  #ifdef GLM_MHD
86  VL = ARRAY_2D(NMAX_POINT, NVAR, double);
87  VR = ARRAY_2D(NMAX_POINT, NVAR, double);
88  UL = ARRAY_2D(NMAX_POINT, NVAR, double);
89  UR = ARRAY_2D(NMAX_POINT, NVAR, double);
90  #endif
91  }
92 
93  #if BACKGROUND_FIELD == YES
94  bgf = GetBackgroundField (beg, end, FACE_CENTER, grid);
95  #endif
96 
97  #if DIVB_CONTROL == EIGHT_WAVES
98  print ("! hlld Riemann solver does not work with Powell's 8-wave\n");
99  QUIT_PLUTO(1);
100  #endif
101 
102  #ifdef GLM_MHD
103 /*
104 double **uL0, **uR0;
105 static double *BxL, *BxR, *psiL, *psiR;
106 if (BxL == NULL){
107  BxL = ARRAY_1D(NMAX_POINT, double);
108  BxR = ARRAY_1D(NMAX_POINT, double);
109  psiL = ARRAY_1D(NMAX_POINT, double);
110  psiR = ARRAY_1D(NMAX_POINT, double);
111 }
112 uL0 = state->uL;
113 uR0 = state->uR;
114 for (i = beg; i <= end; i++){
115  BxL[i] = state->vL[i][BXn];
116  BxR[i] = state->vR[i][BXn];
117  psiL[i] = state->vL[i][PSI_GLM];
118  psiR[i] = state->vL[i][PSI_GLM];
119 }
120 state->uL = UL;
121 state->uR = UR;
122 
123  GLM_Solve (state, beg, end, grid);
124 */
125 
126  GLM_Solve (state, VL, VR, beg, end, grid);
127  PrimToCons (VL, UL, beg, end);
128  PrimToCons (VR, UR, beg, end);
129 
130  #else
131  VL = state->vL; UL = state->uL;
132  VR = state->vR; UR = state->uR;
133  #endif
134 
135 /* ----------------------------------------------------
136  Compute sound speed & fluxes at zone interfaces
137  ---------------------------------------------------- */
138 
139  SoundSpeed2 (VL, a2L, NULL, beg, end, FACE_CENTER, grid);
140  SoundSpeed2 (VR, a2R, NULL, beg, end, FACE_CENTER, grid);
141 
142  Flux (UL, VL, a2L, bgf, fL, ptL, beg, end);
143  Flux (UR, VR, a2R, bgf, fR, ptR, beg, end);
144 
145 /* ----------------------------------------
146  get max and min signal velocities
147  ---------------------------------------- */
148 
149  SL = state->SL; SR = state->SR;
150  HLL_Speed (VL, VR, a2L, a2R, bgf, SL, SR, beg, end);
151 
152  for (i = beg; i <= end; i++) {
153 
154  #if BACKGROUND_FIELD == YES
155  EXPAND (B0n = bgf[i][BXn]; ,
156  B0t = bgf[i][BXt]; ,
157  B0b = bgf[i][BXb];)
158  #endif
159 
160  /* ----------------------------------------
161  get max propagation speed for dt comp.
162  ---------------------------------------- */
163 
164  scrh = MAX(fabs(SL[i]), fabs(SR[i]));
165  cmax[i] = scrh;
166 
167  vL = VL[i]; uL = UL[i];
168  vR = VR[i]; uR = UR[i];
169 
170 /* ----------------------------------------------------------
171  COMPUTE FLUXES and STATES
172  ---------------------------------------------------------- */
173 
174  if (SL[i] >= 0.0){ /* ---- Region L ---- */
175 
176  for (nv = NFLX; nv--; ) state->flux[i][nv] = fL[i][nv];
177  state->press[i] = ptL[i];
178 
179  }else if (SR[i] <= 0.0) { /* ---- Region R ---- */
180 
181  for (nv = NFLX; nv--; ) state->flux[i][nv] = fR[i][nv];
182  state->press[i] = ptR[i];
183 
184  } else {
185 
186 #if SHOCK_FLATTENING == MULTID
187  if ((state->flag[i] & FLAG_HLL) || (state->flag[i+1] & FLAG_HLL)){
188  scrh = 1.0/(SR[i] - SL[i]);
189  for (nv = NFLX; nv--; ){
190  state->flux[i][nv] = SR[i]*SL[i]*(uR[nv] - uL[nv])
191  + SR[i]*fL[i][nv] - SL[i]*fR[i][nv];
192  state->flux[i][nv] *= scrh;
193  }
194  state->press[i] = (SR[i]*ptL[i] - SL[i]*ptR[i])*scrh;
195  continue;
196  }
197 #endif
198 
199  /* ---------------------------
200  Compute U*
201  --------------------------- */
202 
203  scrh = 1.0/(SR[i] - SL[i]);
204  Bx1 = Bx = (SR[i]*vR[BXn] - SL[i]*vL[BXn])*scrh;
205  #if BACKGROUND_FIELD == YES
206  Bx += B0n; /* Bx will be now the (normal) total field */
207  #endif
208  sBx = (Bx > 0.0 ? 1.0 : -1.0);
209 
210  duL = SL[i] - vL[VXn];
211  duR = SR[i] - vR[VXn];
212 
213  scrh = 1.0/(duR*uR[RHO] - duL*uL[RHO]);
214  SM = (duR*uR[MXn] - duL*uL[MXn] - ptR[i] + ptL[i])*scrh;
215 
216  pts = duR*uR[RHO]*ptL[i] - duL*uL[RHO]*ptR[i] +
217  vL[RHO]*vR[RHO]*duR*duL*(vR[VXn]- vL[VXn]);
218  pts *= scrh;
219 
220  usL[RHO] = uL[RHO]*duL/(SL[i] - SM);
221  usR[RHO] = uR[RHO]*duR/(SR[i] - SM);
222 
223  sqrL = sqrt(usL[RHO]);
224  sqrR = sqrt(usR[RHO]);
225 
226  S1L = SM - fabs(Bx)/sqrL;
227  S1R = SM + fabs(Bx)/sqrR;
228 
229  /* -------------------------------------------------------------
230  When S1L -> SL or S1R -> SR a degeneracy occurs.
231  Although Miyoshi & Kusano say that no jump exists, we don't
232  think this is actually true.
233  Indeed, vy*, vz*, By*, Bz* cannot be solved independently.
234  In this case we revert to the HLLC solver of Li (2005),
235  except for the term v.B in the region, which we compute in
236  our own way.
237  Note, that by comparing the expressions of Li (2005) and
238  Miyoshi & Kusano (2005), the only change involves a
239  re-definition of By* and Bz* in terms of By(HLL), Bz(HLL).
240  ------------------------------------------------------------- */
241 
242  revert_to_hllc = 0;
243 
244  if ( (S1L - SL[i]) < 1.e-4*(SM - SL[i]) ) revert_to_hllc = 1;
245  if ( (S1R - SR[i]) > -1.e-4*(SR[i] - SM) ) revert_to_hllc = 1;
246 
247  if (revert_to_hllc){
248 
249  scrh = 1.0/(SR[i] - SL[i]);
250  for (nv = NFLX; nv--; ){
251  Uhll[nv] = SR[i]*uR[nv] - SL[i]*uL[nv] + fL[i][nv] - fR[i][nv];
252  Uhll[nv] *= scrh;
253  }
254 
255  EXPAND(usL[BXn] = usR[BXn] = Uhll[BXn]; ,
256  usL[BXt] = usR[BXt] = Uhll[BXt]; ,
257  usL[BXb] = usR[BXb] = Uhll[BXb];)
258 
259  S1L = S1R = SM; /* region ** should never be computed since */
260  /* fluxes are given in terms of UL* and UR* */
261 
262  }else{
263 
264  scrhL = (uL[RHO]*duL*duL - Bx*Bx)/(uL[RHO]*duL*(SL[i] - SM) - Bx*Bx);
265  scrhR = (uR[RHO]*duR*duR - Bx*Bx)/(uR[RHO]*duR*(SR[i] - SM) - Bx*Bx);
266 
267  EXPAND(usL[BXn] = Bx1; ,
268  usL[BXt] = uL[BXt]*scrhL; ,
269  usL[BXb] = uL[BXb]*scrhL;)
270 
271  #if BACKGROUND_FIELD == YES
272  EXPAND( ; ,
273  usL[BXt] += B0t*(scrhL - 1.0); , /* Eq. [40] of */
274  usL[BXb] += B0b*(scrhL - 1.0);) /* Miyoshi etal (2010) */
275  #endif
276 
277  EXPAND(usR[BXn] = Bx1; ,
278  usR[BXt] = uR[BXt]*scrhR; ,
279  usR[BXb] = uR[BXb]*scrhR;)
280 
281  #if BACKGROUND_FIELD == YES
282  EXPAND( ; ,
283  usR[BXt] += B0t*(scrhR - 1.0); , /* Eq. [40] of */
284  usR[BXb] += B0b*(scrhR - 1.0);) /* Miyoshi etal (2010) */
285  #endif
286 
287  }
288 
289  scrhL = Bx/(uL[RHO]*duL);
290  scrhR = Bx/(uR[RHO]*duR);
291 
292  EXPAND( ; ,
293  vsL = vL[VXt] - scrhL*(usL[BXt] - uL[BXt]);
294  vsR = vR[VXt] - scrhR*(usR[BXt] - uR[BXt]); ,
295 
296  wsL = vL[VXb] - scrhL*(usL[BXb] - uL[BXb]);
297  wsR = vR[VXb] - scrhR*(usR[BXb] - uR[BXb]); )
298 
299  EXPAND(usL[MXn] = usL[RHO]*SM;
300  usR[MXn] = usR[RHO]*SM; ,
301 
302  usL[MXt] = usL[RHO]*vsL;
303  usR[MXt] = usR[RHO]*vsR; ,
304 
305  usL[MXb] = usL[RHO]*wsL;
306  usR[MXb] = usR[RHO]*wsR;)
307 
308  scrhL = EXPAND(vL[VXn]*Bx1, + vL[VXt]*uL[BXt], + vL[VXb]*uL[BXb]);
309  scrhL -= EXPAND( SM*Bx1, + vsL*usL[BXt], + wsL*usL[BXb]);
310 
311  usL[ENG] = duL*uL[ENG] - ptL[i]*vL[VXn] + pts*SM + Bx*scrhL;
312  usL[ENG] /= SL[i] - SM;
313 
314  scrhR = EXPAND(vR[VXn]*Bx1, + vR[VXt]*uR[BXt], + vR[VXb]*uR[BXb]);
315  scrhR -= EXPAND( SM*Bx1, + vsR*usR[BXt], + wsR*usR[BXb]);
316 
317  usR[ENG] = duR*uR[ENG] - ptR[i]*vR[VXn] + pts*SM + Bx*scrhR;
318  usR[ENG] /= SR[i] - SM;
319 
320  #ifdef GLM_MHD
321  usL[PSI_GLM] = usR[PSI_GLM] = vL[PSI_GLM];
322  #endif
323 
324  /* ------------------------------
325  compute HLLD flux
326  ------------------------------ */
327 
328  if (S1L >= 0.0){ /* ---- Region L* ---- */
329 
330  for (nv = NFLX; nv--; ){
331  state->flux[i][nv] = fL[i][nv] + SL[i]*(usL[nv] - uL[nv]);
332  }
333  state->press[i] = ptL[i];
334 
335  }else if (S1R <= 0.0) { /* ---- Region R* ---- */
336 
337  for (nv = NFLX; nv--; ){
338  state->flux[i][nv] = fR[i][nv] + SR[i]*(usR[nv] - uR[nv]);
339  }
340  state->press[i] = ptR[i];
341 
342  } else { /* -- This state exists only if B_x != 0 -- */
343 
344  /* ---------------------------
345  Compute U**
346  --------------------------- */
347 
348  ussl[RHO] = usL[RHO];
349  ussr[RHO] = usR[RHO];
350 
351  EXPAND( ,
352 
353  vss = sqrL*vsL + sqrR*vsR + (usR[BXt] - usL[BXt])*sBx;
354  vss /= sqrL + sqrR; ,
355 
356  wss = sqrL*wsL + sqrR*wsR + (usR[BXb] - usL[BXb])*sBx;
357  wss /= sqrL + sqrR;)
358 
359  EXPAND(ussl[MXn] = ussl[RHO]*SM;
360  ussr[MXn] = ussr[RHO]*SM; ,
361 
362  ussl[MXt] = ussl[RHO]*vss;
363  ussr[MXt] = ussr[RHO]*vss; ,
364 
365  ussl[MXb] = ussl[RHO]*wss;
366  ussr[MXb] = ussr[RHO]*wss;)
367 
368  EXPAND(ussl[BXn] = ussr[BXn] = Bx1; ,
369 
370  ussl[BXt] = sqrL*usR[BXt] + sqrR*usL[BXt] + sqrL*sqrR*(vsR - vsL)*sBx;
371  ussl[BXt] /= sqrL + sqrR;
372  ussr[BXt] = ussl[BXt]; ,
373 
374  ussl[BXb] = sqrL*usR[BXb] + sqrR*usL[BXb] + sqrL*sqrR*(wsR - wsL)*sBx;
375  ussl[BXb] /= sqrL + sqrR;
376  ussr[BXb] = ussl[BXb];)
377 
378  scrhL = EXPAND(SM*Bx1, + vsL*usL [BXt], + wsL*usL [BXb]);
379  scrhL -= EXPAND(SM*Bx1, + vss*ussl[BXt], + wss*ussl[BXb]);
380 
381  scrhR = EXPAND(SM*Bx1, + vsR*usR [BXt], + wsR*usR [BXb]);
382  scrhR -= EXPAND(SM*Bx1, + vss*ussr[BXt], + wss*ussr[BXb]);
383 
384  ussl[ENG] = usL[ENG] - sqrL*scrhL*sBx;
385  ussr[ENG] = usR[ENG] + sqrR*scrhR*sBx;
386 
387  #ifdef GLM_MHD
388  ussl[PSI_GLM] = ussr[PSI_GLM] = vL[PSI_GLM];
389  #endif
390 
391  /* --------------------------------------
392  verify consistency condition
393  -------------------------------------- */
394 
395 /*
396  for (nv = 0; nv < NFLX; nv++){
397  scrh = (SR[i] - S1R)*usR[nv] + (S1R - SM)*ussr[nv] +
398  (SM - S1L)*ussl[nv] + (S1L - SL[i])*usL[nv] -
399  SR[i]*UR[nv] + SL[i]*UL[nv] + fr[i][nv] - fl[i][nv];
400 
401  if (fabs(scrh) > 1.e-2){
402  printf (" ! Consistency condition violated, pt %d, nv %d, %12.6e \n",
403  i,nv,scrh);
404  printf ("%f %f %f %f %f %f\n",UL[BXn], usL[BXn], ussl[BXn],
405  ussr[BXn], usR[BXn], UR[BXn]);
406  printf ("%f %f %f %f %f\n",SL[i], S1L, SM, S1R, SR[i]);
407  printf ("%f %f %d\n",fr[i][nv],fl[i][nv],BXn);
408  exit(1);
409  }
410  }
411 */
412 
413  if (SM >= 0.0){ /* ---- Region L** ---- */
414 
415  for (nv = NFLX; nv--; ){
416  state->flux[i][nv] = fL[i][nv] + S1L*(ussl[nv] - usL[nv])
417  + SL[i]*(usL[nv] - uL[nv]);
418  }
419  state->press[i] = ptL[i];
420  }else{ /* ---- Region R** ---- */
421 
422  for (nv = NFLX; nv--; ){
423  state->flux[i][nv] = fR[i][nv] + S1R*(ussr[nv] - usR[nv])
424  + SR[i]*(usR[nv] - uR[nv]);
425  }
426  state->press[i] = ptR[i];
427  }
428  }
429  }
430  }
431 }
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
GetBackgroundField
double ** GetBackgroundField(int beg, int end, int where, Grid *grid)
Definition: bckgrnd_field.c:5
VXb
int VXb
Definition: globals.h:73
SoundSpeed2
void SoundSpeed2(double **v, double *cs2, double *h, int beg, int end, int pos, Grid *grid)
Definition: eos.c:16
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
VXt
int VXt
Definition: globals.h:73
STATE_1D::flag
unsigned char * flag
Definition: structs.h:168
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
MXn
int MXn
Definition: globals.h:74
VXn
int VXn
Definition: globals.h:73
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
BACKGROUND_FIELD
void BACKGROUND_FIELD(real x1, real x2, real x3, real *B0)
Definition: init.c:95
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
QUIT_PLUTO
#define QUIT_PLUTO(e_code)
Definition: macros.h:125
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: