PLUTO
hlld.c
Go to the documentation of this file.
1 /* ///////////////////////////////////////////////////////////////////// */
2 /*!
3  \file
4  \brief HLLD Riemann solver for the classical MHD equations.
5 
6  Solve the Riemann problem for the classical MHD equations with
7  an ideal or isothermal equation of state.
8  For the ideal case, we implement the four-state (of five-wave)
9  HLLD solver of Miyoshi & Kusano (2005).
10  For the isothermal case, we use the three-state approximation
11  described by Mignone (2007).
12 
13  The macro VERIFY_CONSISTENCY_CONDITION can be turned to YES to verify
14  the correctness of the implementation.
15 
16  On input, this function takes left and right primitive state vectors
17  \c state->vL and \c state->vR at zone edge i+1/2;
18  On output, return flux and pressure vectors at the same interface
19  \c i+1/2 (note that the \c i refers to \c i+1/2).
20 
21  Also during this step, compute maximum wave propagation speed (cmax)
22  for explicit time step computation.
23 
24  \b Reference:
25  - "A muLti-state HLL approximate Riemann Solver for ideal MHD",
26  Miyoshi, T., Kusano, K., JCP 2005
27  - "A simple and accurate Riemann solver for isothermal MHD",
28  Mignone, JCP (2007) 225, 1427.
29 
30  \authors A. Mignone (mignone@ph.unito.it)
31  C. Zanni (zanni@oato.inaf.it)
32  \date April 14, 2014
33 */
34 /* ///////////////////////////////////////////////////////////////////// */
35 #include"pluto.h"
36 
37 #define VERIFY_CONSISTENCY_CONDITION NO
38 
39 #if EOS == IDEAL || EOS == PVTE_LAW
40 /* ********************************************************************* */
41 void HLLD_Solver (const State_1D *state, int beg, int end,
42  double *cmax, Grid *grid)
43 /*!
44  * Solve Riemann problem for the adiabatic MHD equations using the
45  * four-state HLLD Riemann solver of Miyoshi & Kusano (2005).
46  *
47  * \param[in,out] state pointer to State_1D structure
48  * \param[in] beg initial grid index
49  * \param[out] end final grid index
50  * \param[out] cmax 1D array of maximum characteristic speeds
51  * \param[in] grid pointer to array of Grid structures.
52  *
53  *********************************************************************** */
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 }
432 #endif
433 
434 #if EOS == ISOTHERMAL
435 /* ********************************************************************* */
436 void HLLD_Solver (const State_1D *state, int beg, int end,
437  double *cmax, Grid *grid)
438 /*!
439  * Solve Riemann problem for the isothermal MHD equations using the
440  * three-state HLLD Riemann solver of Mignone (2007).
441  *
442  * \param[in,out] state pointer to State_1D structure
443  * \param[in] beg initial grid index
444  * \param[out] end final grid index
445  * \param[out] cmax 1D array of maximum characteristic speeds
446  * \param[in] grid pointer to array of Grid structures.
447  *
448  *********************************************************************** */
449 {
450  int nv, i;
451  int revert_to_hll;
452  double scrh;
453  double usL[NFLX], *SL;
454  double usR[NFLX], *SR;
455  double usc[NFLX];
456 
457  double scrhL, S1L, duL;
458  double scrhR, S1R, duR;
459  double Bx, Bx1, SM, sBx, rho, sqrho;
460 
461  double *vL, *vR, *uL, *uR;
462  static double *ptL, *ptR, *a2L, *a2R;
463  static double **fL, **fR;
464  static double **VL, **VR, **UL, **UR;
465  #if BACKGROUND_FIELD == YES
466  double B0n, B0t, B0b;
467  #endif
468  double **bgf;
469 
470  if (fL == NULL){
471  fL = ARRAY_2D(NMAX_POINT, NFLX, double);
472  fR = ARRAY_2D(NMAX_POINT, NFLX, double);
473 
474  ptR = ARRAY_1D(NMAX_POINT, double);
475  ptL = ARRAY_1D(NMAX_POINT, double);
476 
477  a2R = ARRAY_1D(NMAX_POINT, double);
478  a2L = ARRAY_1D(NMAX_POINT, double);
479 
480  #ifdef GLM_MHD
481  VL = ARRAY_2D(NMAX_POINT, NVAR, double);
482  VR = ARRAY_2D(NMAX_POINT, NVAR, double);
483  UL = ARRAY_2D(NMAX_POINT, NVAR, double);
484  UR = ARRAY_2D(NMAX_POINT, NVAR, double);
485  #endif
486  }
487 
488  #if BACKGROUND_FIELD == YES
489  bgf = GetBackgroundField (beg, end, FACE_CENTER, grid);
490  #endif
491 
492  #if DIVB_CONTROL == EIGHT_WAVES
493  print ("! hlld Riemann solver does not work with Powell\n");
494  QUIT_PLUTO(1);
495  #endif
496 
497  #ifdef GLM_MHD
498  GLM_Solve (state, VL, VR, beg, end, grid);
499  PrimToCons (VL, UL, beg, end);
500  PrimToCons (VR, UR, beg, end);
501  #else
502  VL = state->vL; UL = state->uL;
503  VR = state->vR; UR = state->uR;
504  #endif
505 
506 /* ----------------------------------------------------
507  compute sound speed & fluxes at zone interfaces
508  ---------------------------------------------------- */
509 
510  SoundSpeed2 (VL, a2L, NULL, beg, end, FACE_CENTER, grid);
511  SoundSpeed2 (VR, a2R, NULL, beg, end, FACE_CENTER, grid);
512 
513  Flux (UL, VL, a2L, bgf, fL, ptL, beg, end);
514  Flux (UR, VR, a2R, bgf, fR, ptR, beg, end);
515 
516 /* ----------------------------------------
517  get max and min signal velocities
518  ---------------------------------------- */
519 
520  SL = state->SL; SR = state->SR;
521  HLL_Speed (VL, VR, a2L, a2R, bgf, SL, SR, beg, end);
522 
523  for (i = beg; i <= end; i++) {
524 
525  #if BACKGROUND_FIELD == YES
526  EXPAND (B0n = bgf[i][BXn]; ,
527  B0t = bgf[i][BXt]; ,
528  B0b = bgf[i][BXb];)
529  #endif
530 
531  /* ----------------------------------------
532  get max propagation speed for dt comp.
533  ---------------------------------------- */
534 
535  scrh = MAX(fabs(SL[i]), fabs(SR[i]));
536  cmax[i] = scrh;
537 
538  vL = VL[i]; uL = UL[i];
539  vR = VR[i]; uR = UR[i];
540 
541 /* ----------------------------------------------------------
542  COMPUTE FLUXES and STATES
543  ---------------------------------------------------------- */
544 
545  if (SL[i] >= 0.0){ /* ---- Region L ---- */
546 
547  for (nv = NFLX; nv--; ) state->flux[i][nv] = fL[i][nv];
548  state->press[i] = ptL[i];
549 
550  }else if (SR[i] <= 0.0) { /* ---- Region R ---- */
551 
552  for (nv = NFLX; nv--; ) state->flux[i][nv] = fR[i][nv];
553  state->press[i] = ptR[i];
554 
555  } else {
556 
557  scrh = 1.0/(SR[i] - SL[i]);
558  duL = SL[i] - vL[VXn];
559  duR = SR[i] - vR[VXn];
560 
561  Bx1 = Bx = (SR[i]*vR[BXn] - SL[i]*vL[BXn])*scrh;
562  #if BACKGROUND_FIELD == YES
563  Bx += B0n; /* total field */
564  #endif
565 
566  rho = (uR[RHO]*duR - uL[RHO]*duL)*scrh;
567  state->flux[i][RHO] = (SL[i]*uR[RHO]*duR - SR[i]*uL[RHO]*duL)*scrh;
568 
569  /* ---------------------------
570  compute S*
571  --------------------------- */
572 
573  sqrho = sqrt(rho);
574 
575  SM = state->flux[i][RHO]/rho;
576  S1L = SM - fabs(Bx)/sqrho;
577  S1R = SM + fabs(Bx)/sqrho;
578 
579  /* ---------------------------------------------
580  Prevent degeneracies when S1L -> SL or
581  S1R -> SR. Revert to HLL if necessary.
582  --------------------------------------------- */
583 
584  revert_to_hll = 0;
585 
586  if ( (S1L - SL[i]) < 1.e-4*(SR[i] - SL[i]) ) revert_to_hll = 1;
587  if ( (S1R - SR[i]) > -1.e-4*(SR[i] - SL[i]) ) revert_to_hll = 1;
588 
589  if (revert_to_hll){
590  scrh = 1.0/(SR[i] - SL[i]);
591  for (nv = NFLX; nv--; ){
592  state->flux[i][nv] = SL[i]*SR[i]*(uR[nv] - uL[nv]) +
593  SR[i]*fL[i][nv] - SL[i]*fR[i][nv];
594  state->flux[i][nv] *= scrh;
595  }
596  state->press[i] = (SR[i]*ptL[i] - SL[i]*ptR[i])*scrh;
597  continue;
598  }
599 
600  state->flux[i][MXn] = (SR[i]*fL[i][MXn] - SL[i]*fR[i][MXn]
601  + SR[i]*SL[i]*(uR[MXn] - uL[MXn]))*scrh;
602 
603  state->press[i] = (SR[i]*ptL[i] - SL[i]*ptR[i])*scrh;
604  #ifdef GLM_MHD
605  state->flux[i][BXn] = fL[i][BXn];
606  state->flux[i][PSI_GLM] = fL[i][PSI_GLM];
607  #else
608  state->flux[i][BXn] = SR[i]*SL[i]*(uR[BXn] - uL[BXn])*scrh;
609  #endif
610 
611  /* ---------------------------
612  Compute U*
613  --------------------------- */
614 
615  scrhL = 1.0/((SL[i] - S1L)*(SL[i] - S1R));
616  scrhR = 1.0/((SR[i] - S1L)*(SR[i] - S1R));
617 
618  #if BACKGROUND_FIELD == YES
619  EXPAND( ; ,
620  usL[MXt] = rho*vL[VXt] - Bx*(uL[BXt]+B0t)*(SM - vL[VXn])*scrhL;
621  usR[MXt] = rho*vR[VXt] - Bx*(uR[BXt]+B0t)*(SM - vR[VXn])*scrhR; ,
622 
623  usL[MXb] = rho*vL[VXb] - Bx*(uL[BXb]+B0b)*(SM - vL[VXn])*scrhL;
624  usR[MXb] = rho*vR[VXb] - Bx*(uR[BXb]+B0b)*(SM - vR[VXn])*scrhR;)
625  #else
626  EXPAND( ; ,
627  usL[MXt] = rho*vL[VXt] - Bx*uL[BXt]*(SM - vL[VXn])*scrhL;
628  usR[MXt] = rho*vR[VXt] - Bx*uR[BXt]*(SM - vR[VXn])*scrhR; ,
629 
630  usL[MXb] = rho*vL[VXb] - Bx*uL[BXb]*(SM - vL[VXn])*scrhL;
631  usR[MXb] = rho*vR[VXb] - Bx*uR[BXb]*(SM - vR[VXn])*scrhR;)
632  #endif
633 
634  #if BACKGROUND_FIELD == YES
635  EXPAND( ; ,
636  usL[BXt] = (uL[BXt]+B0t)/rho*(uL[RHO]*duL*duL - Bx*Bx)*scrhL-B0t;
637  usR[BXt] = (uR[BXt]+B0t)/rho*(uR[RHO]*duR*duR - Bx*Bx)*scrhR-B0t; ,
638 
639  usL[BXb] = (uL[BXb]+B0b)/rho*(uL[RHO]*duL*duL - Bx*Bx)*scrhL-B0b;
640  usR[BXb] = (uR[BXb]+B0b)/rho*(uR[RHO]*duR*duR - Bx*Bx)*scrhR-B0b;)
641  #else
642  EXPAND( ; ,
643  usL[BXt] = uL[BXt]/rho*(uL[RHO]*duL*duL - Bx*Bx)*scrhL;
644  usR[BXt] = uR[BXt]/rho*(uR[RHO]*duR*duR - Bx*Bx)*scrhR; ,
645 
646  usL[BXb] = uL[BXb]/rho*(uL[RHO]*duL*duL - Bx*Bx)*scrhL;
647  usR[BXb] = uR[BXb]/rho*(uR[RHO]*duR*duR - Bx*Bx)*scrhR;)
648  #endif
649 
650  if (S1L >= 0.0){ /* ---- Region L* ---- */
651 
652  EXPAND( ; ,
653  state->flux[i][MXt] = fL[i][MXt] + SL[i]*(usL[MXt] - uL[MXt]); ,
654  state->flux[i][MXb] = fL[i][MXb] + SL[i]*(usL[MXb] - uL[MXb]);
655  )
656  EXPAND( ; ,
657  state->flux[i][BXt] = fL[i][BXt] + SL[i]*(usL[BXt] - uL[BXt]); ,
658  state->flux[i][BXb] = fL[i][BXb] + SL[i]*(usL[BXb] - uL[BXb]);
659  )
660 
661  }else if (S1R <= 0.0) { /* ---- Region R* ---- */
662 
663  EXPAND( ; ,
664  state->flux[i][MXt] = fR[i][MXt] + SR[i]*(usR[MXt] - uR[MXt]); ,
665  state->flux[i][MXb] = fR[i][MXb] + SR[i]*(usR[MXb] - uR[MXb]);
666  )
667  EXPAND( ; ,
668  state->flux[i][BXt] = fR[i][BXt] + SR[i]*(usR[BXt] - uR[BXt]); ,
669  state->flux[i][BXb] = fR[i][BXb] + SR[i]*(usR[BXb] - uR[BXb]);
670  )
671 
672  } else {
673 
674  /* ---------------------------
675  Compute U** = Uc
676  --------------------------- */
677 
678  sBx = (Bx > 0.0 ? 1.0 : -1.0);
679 
680  EXPAND( ; ,
681  usc[MXt] = 0.5*(usR[MXt] + usL[MXt]
682  + (usR[BXt] - usL[BXt])*sBx*sqrho); ,
683  usc[MXb] = 0.5*( usR[MXb] + usL[MXb]
684  + (usR[BXb] - usL[BXb])*sBx*sqrho);)
685 
686  EXPAND( ; ,
687  usc[BXt] = 0.5*( usR[BXt] + usL[BXt]
688  + (usR[MXt] - usL[MXt])*sBx/sqrho); ,
689  usc[BXb] = 0.5*( usR[BXb] + usL[BXb]
690  + (usR[MXb] - usL[MXb])*sBx/sqrho);)
691 
692  EXPAND( ; ,
693  state->flux[i][MXt] = usc[MXt]*SM - Bx*usc[BXt]; ,
694  state->flux[i][MXb] = usc[MXb]*SM - Bx*usc[BXb]; )
695  #if BACKGROUND_FIELD == YES
696  EXPAND( ; ,
697  state->flux[i][MXt] -= Bx1*B0t; ,
698  state->flux[i][MXb] -= Bx1*B0b; )
699  #endif
700 
701  EXPAND( ; ,
702  state->flux[i][BXt] = usc[BXt]*SM - Bx*usc[MXt]/rho; ,
703  state->flux[i][BXb] = usc[BXb]*SM - Bx*usc[MXb]/rho;)
704  #if BACKGROUND_FIELD == YES
705  EXPAND( ; ,
706  state->flux[i][BXt] += B0t*SM; ,
707  state->flux[i][BXb] += B0b*SM;)
708  #endif
709 
710  /* --------------------------------------
711  verify consistency condition
712  -------------------------------------- */
713 
714  #if VERIFY_CONSISTENCY_CONDITION == YES
715  for (nv = NFLX; nv--; ){
716  if (nv == RHO || nv == MXn || nv == BXn) continue;
717  scrh = (S1L - SL[i])*usL[nv] + (S1R - S1L)*usc[nv] +
718  (SR[i] - S1R)*usR[nv] -
719  SR[i]*uR[nv] + SL[i]*uL[nv] + fR[i][nv] - fL[i][nv];
720 
721  if (fabs(scrh) > 1.e-6){
722  printf (" ! Consistency condition violated, pt %d, nv %d, %12.6e \n",
723  i,nv,scrh);
724  printf (" scrhL = %12.6e scrhR = %12.6e\n",scrhL, scrhR);
725  printf (" SL = %12.6e, S1L = %12.6e, S1R = %12.6e, SR = %12.6e\n",
726  SL[i],S1L,S1R, SR[i]);
727  Show(state->vL,i);
728  Show(state->vR,i);
729 
730  exit(1);
731  }
732  }
733  #endif
734 
735  }
736  }
737  }
738 }
739 #endif /* end #if on EOS */
740 
741 #undef VERIFY_CONSISTENCY_CONDITION
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
HLLD_Solver
void HLLD_Solver(const State_1D *state, int beg, int end, double *cmax, Grid *grid)
Definition: hlld.c:41
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
GRID
Definition: structs.h:78
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
pluto.h
PLUTO main header file.
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
VERIFY_CONSISTENCY_CONDITION
#define VERIFY_CONSISTENCY_CONDITION
Definition: hlld.c:37
Show
void Show(double **, int)
Definition: tools.c:132
NMAX_POINT
long int NMAX_POINT
Maximum number of points among the three directions, boundaries excluded.
Definition: globals.h:62
STATE_1D
Definition: structs.h:133
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