PLUTO
Macros | Functions
roe.c File Reference

Implementation of the Roe Riemann solver for the MHD equations. More...

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

Go to the source code of this file.

Macros

#define sqrt_1_2   (0.70710678118654752440)
 
#define HLL_HYBRIDIZATION   NO
 
#define CHECK_ROE_MATRIX   NO
 

Functions

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

Detailed Description

Implementation of the Roe Riemann solver for the MHD equations.

Solve the Riemann problem for the adiabatic and isothermal MHD equations using the linearized Riemann solver of Roe. The implementation follows the approach outlined by Cargo & Gallice (1997).

The isothermal version is recovered by taking the limit of $ \bar{a}^2 $ for $ \gamma \to 1$ which gives (page 451) $\bar{a}^2 \to {\rm g\_isoSoundSpeed2} + X$ where X is defined as in the adiabatic case. This follows by imposing zero jump across the entropy wave (first of Eq. 4.20, page 452) giving $ \Delta p = (\bar{a}^2 - X)\Delta\rho$. Then all the terms like $ [X\Delta\rho + \Delta p] \to \bar{a}^2\Delta\rho $.

When the macro HLL_HYBRIDIZATION is set to YES, we revert to the HLL solver whenever an unphysical state appear in the solution.

The macro CHECK_ROE_MATRIX can be used to verify that the characteristic decomposition reproduces the Roe matrix. This can be done also when BACKGROUND_FIELD is set to YES.

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:

  • "Roe Matrices for Ideal MHD and Systematic Construction of Roe Matrices for Systems of Conservation Laws", P. Cargo, G. Gallice, JCP (1997), 136, 446.
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 11, 2014

Definition in file roe.c.

Macro Definition Documentation

#define CHECK_ROE_MATRIX   NO

Definition at line 50 of file roe.c.

#define HLL_HYBRIDIZATION   NO

Definition at line 49 of file roe.c.

#define sqrt_1_2   (0.70710678118654752440)

Definition at line 48 of file roe.c.

Function Documentation

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

Solve Riemann problem for the adiabatic MHD equations using the Roe Riemann solver of Cargo & Gallice (1997).

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 53 of file roe.c.

66 {
67  int nv, i, j, k;
68  int ifail;
69  double rho, u, v, w, vel2, bx, by, bz, pr;
70  double a2, a, ca2, cf2, cs2;
71  double cs, ca, cf, b2;
72  double alpha_f, alpha_s, beta_y, beta_z, beta_v, scrh, sBx;
73  double dV[NFLX], dU[NFLX], *vL, *vR, *uL, *uR;
74  double *SL, *SR;
75  double Rc[NFLX][NFLX], eta[NFLX], lambda[NFLX];
76  double alambda[NFLX], Uv[NFLX];
77 
78  double sqrt_rho;
79  double delta, delta_inv;
80 
81  double g1, sl, sr, H, Hgas, HL, HR, Bx, By, Bz, X;
82  double vdm, BdB, beta_dv, beta_dB;
83  double bt2, Btmag, sqr_rho_L, sqr_rho_R;
84 
85  static double *pR, *pL, *a2L, *a2R;
86  static double **fL, **fR;
87  static double **VL, **VR, **UL, **UR;
88  double **bgf;
89  #if BACKGROUND_FIELD == YES
90  double B0x, B0y, B0z, B1x, B1y, B1z;
91  #endif
92  double Us[NFLX];
93  delta = 1.e-6;
94 
95 /* -----------------------------------------------------------
96  1. Allocate static memory areas
97  ----------------------------------------------------------- */
98 
99  if (fL == NULL){
100 
101  fL = ARRAY_2D(NMAX_POINT, NFLX, double);
102  fR = ARRAY_2D(NMAX_POINT, NFLX, double);
103 
104  pR = ARRAY_1D(NMAX_POINT, double);
105  pL = ARRAY_1D(NMAX_POINT, double);
106 
107  a2R = ARRAY_1D(NMAX_POINT, double);
108  a2L = ARRAY_1D(NMAX_POINT, double);
109 
110  #ifdef GLM_MHD
111  VL = ARRAY_2D(NMAX_POINT, NVAR, double);
112  VR = ARRAY_2D(NMAX_POINT, NVAR, double);
113  UL = ARRAY_2D(NMAX_POINT, NVAR, double);
114  UR = ARRAY_2D(NMAX_POINT, NVAR, double);
115  #endif
116  }
117 
118 /* -----------------------------------------------------------
119  2. Background field
120  ----------------------------------------------------------- */
121 
122  #if BACKGROUND_FIELD == YES
123  bgf = GetBackgroundField (beg, end, FACE_CENTER, grid);
124  #if (DIVB_CONTROL == EIGHT_WAVES)
125  print1 ("! Roe_Solver: background field and 8wave not tested\n");
126  QUIT_PLUTO(1);
127  #endif
128  #endif
129 
130 /* -----------------------------------------------------------
131  3. GLM pre-Rieman solver
132  ----------------------------------------------------------- */
133 
134  #ifdef GLM_MHD
135  GLM_Solve (state, VL, VR, beg, end, grid);
136  PrimToCons (VL, UL, beg, end);
137  PrimToCons (VR, UR, beg, end);
138  #else
139  VL = state->vL; UL = state->uL;
140  VR = state->vR; UR = state->uR;
141  #endif
142 
143 /* ----------------------------------------------------
144  4. Compute sound speed & fluxes at zone interfaces
145  ---------------------------------------------------- */
146 
147  SoundSpeed2 (VL, a2L, NULL, beg, end, FACE_CENTER, grid);
148  SoundSpeed2 (VR, a2R, NULL, beg, end, FACE_CENTER, grid);
149 
150  Flux (UL, VL, a2L, bgf, fL, pL, beg, end);
151  Flux (UR, VR, a2R, bgf, fR, pR, beg, end);
152 
153  SL = state->SL; SR = state->SR;
154 
155  #if EOS == IDEAL
156  g1 = g_gamma - 1.0;
157  #endif
158 
159 /* -------------------------------------------------
160  5. Some eigenvectors components will always be
161  zero so set Rc = 0 initially
162  -------------------------------------------------- */
163 
164  for (k = NFLX; k--; ) {
165  for (j = NFLX; j--; ) {
166  Rc[k][j] = 0.0;
167  }}
168 
169 /* ---------------------------------------------------------------------
170  6. Begin main loop
171  --------------------------------------------------------------------- */
172 
173  for (i = beg; i <= end; i++) {
174 
175  vL = VL[i]; uL = UL[i];
176  vR = VR[i]; uR = UR[i];
177 
178  /* -----------------------------------------------------------
179  6a. switch to HLL in proximity of strong shock
180  ----------------------------------------------------------- */
181 
182  #if SHOCK_FLATTENING == MULTID
183  if ((state->flag[i] & FLAG_HLL) || (state->flag[i+1] & FLAG_HLL)){
184  HLL_Speed (VL, VR, a2L, a2R, NULL, SL, SR, i, i);
185 
186  scrh = MAX(fabs(SL[i]), fabs(SR[i]));
187  cmax[i] = scrh;
188 
189  if (SL[i] > 0.0) {
190  for (nv = NFLX; nv--; ) state->flux[i][nv] = fL[i][nv];
191  state->press[i] = pL[i];
192  } else if (SR[i] < 0.0) {
193  for (nv = NFLX; nv--; ) state->flux[i][nv] = fR[i][nv];
194  state->press[i] = pR[i];
195  }else{
196  scrh = 1.0/(SR[i] - SL[i]);
197  for (nv = NFLX; nv--; ){
198  state->flux[i][nv] = SR[i]*SL[i]*(uR[nv] - uL[nv])
199  + SR[i]*fL[i][nv] - SL[i]*fR[i][nv];
200  state->flux[i][nv] *= scrh;
201  }
202  state->press[i] = (SR[i]*pL[i] - SL[i]*pR[i])*scrh;
203  }
204  continue;
205  }
206  #endif
207 
208  /* ----------------------------------------------------------------
209  6b. Compute jumps in primitive and conservative variables.
210 
211  Note on the velocity jump:
212  Formally the jump in velocity should be written using
213  conservative variables:
214 
215  Delta(v) = Delta(m/rho) = (Delta(m) - v*Delta(rho))/rho
216 
217  with rho and v defined by the second and first of [4.6]
218  Still, since we reconstruct v and not m, the following MAPLE
219  script shows that the previous expression simplifies to
220  Delta(v) = vR - vL.
221 
222  restart;
223  sR := sqrt(rho[R])/(sqrt(rho[R]) + sqrt(rho[L])):
224  sL := sqrt(rho[L])/(sqrt(rho[R]) + sqrt(rho[L])):
225  rhoa := sL*rho[R] + sR*rho[L]:
226  va := sL*v[L] + sR*v[R]:
227  dV := (rho[R]*v[R] - rho[L]*v[L] - va*(rho[R]-rho[L]))/rhoa:
228  simplify(dV);
229  ---------------------------------------------------------------- */
230 
231  for (nv = 0; nv < NFLX; nv++) {
232  dV[nv] = vR[nv] - vL[nv];
233  dU[nv] = uR[nv] - uL[nv];
234  }
235 
236  /* ---------------------------------
237  6c. Compute Roe averages
238  --------------------------------- */
239 
240  sqr_rho_L = sqrt(vL[RHO]);
241  sqr_rho_R = sqrt(vR[RHO]);
242 
243  sl = sqr_rho_L/(sqr_rho_L + sqr_rho_R);
244  sr = sqr_rho_R/(sqr_rho_L + sqr_rho_R);
245 
246 /* sl = sr = 0.5; */
247 
248  rho = sr*vL[RHO] + sl*vR[RHO];
249 
250  sqrt_rho = sqrt(rho);
251 
252  EXPAND (u = sl*vL[VXn] + sr*vR[VXn]; ,
253  v = sl*vL[VXt] + sr*vR[VXt]; ,
254  w = sl*vL[VXb] + sr*vR[VXb];)
255 
256  EXPAND (Bx = sr*vL[BXn] + sl*vR[BXn]; ,
257  By = sr*vL[BXt] + sl*vR[BXt]; ,
258  Bz = sr*vL[BXb] + sl*vR[BXb];)
259 
260  #if BACKGROUND_FIELD == YES
261  /* -- Define field B0 and total B. B1 is the deviation -- */
262  EXPAND (B0x = bgf[i][BXn]; B1x = sr*vL[BXn] + sl*vR[BXn]; Bx = B0x + B1x; ,
263  B0y = bgf[i][BXt]; B1y = sr*vL[BXt] + sl*vR[BXt]; By = B0y + B1y; ,
264  B0z = bgf[i][BXb]; B1z = sr*vL[BXb] + sl*vR[BXb]; Bz = B0z + B1z;)
265  #else
266  EXPAND (Bx = sr*vL[BXn] + sl*vR[BXn]; ,
267  By = sr*vL[BXt] + sl*vR[BXt]; ,
268  Bz = sr*vL[BXb] + sl*vR[BXb];)
269  #endif
270 
271  sBx = (Bx >= 0.0 ? 1.0 : -1.0);
272 
273  EXPAND(bx = Bx/sqrt_rho; ,
274  by = By/sqrt_rho; ,
275  bz = Bz/sqrt_rho; )
276 
277  bt2 = EXPAND(0.0 , + by*by, + bz*bz);
278  b2 = bx*bx + bt2;
279  Btmag = sqrt(bt2*rho);
280 
281  X = EXPAND(dV[BXn]*dV[BXn], + dV[BXt]*dV[BXt], + dV[BXb]*dV[BXb]);
282  X /= (sqr_rho_L + sqr_rho_R)*(sqr_rho_L + sqr_rho_R)*2.0;
283 
284  vdm = EXPAND(u*dU[MXn], + v*dU[MXt], + w*dU[MXb]);
285  #if BACKGROUND_FIELD == YES /* BdB = B1.dB1 (deviation only) */
286  BdB = EXPAND(B1x*dU[BXn], + B1y*dU[BXt], + B1z*dU[BXb]);
287  #else
288  BdB = EXPAND(Bx*dU[BXn], + By*dU[BXt], + Bz*dU[BXb]);
289  #endif
290 
291  /* ---------------------------------------
292  6d. Compute enthalpy and sound speed.
293  --------------------------------------- */
294 
295  #if EOS == ISOTHERMAL
296  a2 = 0.5*(a2L[i] + a2R[i]) + X; /* in most cases a2L = a2R
297  for isothermal MHD */
298  #elif EOS == BAROTROPIC
299  print ("! Roe_Solver: not implemented for barotropic EOS\n");
300  QUIT_PLUTO(1);
301  #elif EOS == IDEAL
302  vel2 = EXPAND(u*u, + v*v, + w*w);
303  dV[PRS] = g1*((0.5*vel2 - X)*dV[RHO] - vdm + dU[ENG] - BdB);
304 
305  HL = (uL[ENG] + pL[i])/vL[RHO];
306  HR = (uR[ENG] + pR[i])/vR[RHO];
307  H = sl*HL + sr*HR; /* total enthalpy */
308 
309  #if BACKGROUND_FIELD == YES
310  scrh = EXPAND(B1x*Bx, + B1y*By, + B1z*Bz);
311  Hgas = H - scrh/rho; /* gas enthalpy */
312  #else
313  Hgas = H - b2; /* gas enthalpy */
314  #endif
315 
316  a2 = (2.0 - g_gamma)*X + g1*(Hgas - 0.5*vel2);
317  if (a2 < 0.0) {
318  printf ("! Roe: a2 = %12.6e < 0.0 !! \n",a2);
319  Show(VL,i);
320  Show(VR,i);
321  QUIT_PLUTO(1);
322  }
323  #endif /* EOS == IDEAL */
324 
325  /* ------------------------------------------------------------
326  6e. Compute fast and slow magnetosonic speeds.
327 
328  The following expression appearing in the definitions
329  of the fast magnetosonic speed
330 
331  (a^2 - b^2)^2 + 4*a^2*bt^2 = (a^2 + b^2)^2 - 4*a^2*bx^2
332 
333  is always positive and avoids round-off errors.
334 
335  Note that we always use the total field to compute the
336  characteristic speeds.
337  ------------------------------------------------------------ */
338 
339  scrh = a2 - b2;
340  ca2 = bx*bx;
341  scrh = scrh*scrh + 4.0*bt2*a2;
342  scrh = sqrt(scrh);
343 
344  cf2 = 0.5*(a2 + b2 + scrh);
345  cs2 = a2*ca2/cf2; /* -- same as 0.5*(a2 + b2 - scrh) -- */
346 
347  cf = sqrt(cf2);
348  cs = sqrt(cs2);
349  ca = sqrt(ca2);
350  a = sqrt(a2);
351 
352  if (cf == cs) {
353  alpha_f = 1.0;
354  alpha_s = 0.0;
355  }else if (a <= cs) {
356  alpha_f = 0.0;
357  alpha_s = 1.0;
358  }else if (cf <= a){
359  alpha_f = 1.0;
360  alpha_s = 0.0;
361  }else{
362  scrh = 1.0/(cf2 - cs2);
363  alpha_f = (a2 - cs2)*scrh;
364  alpha_s = (cf2 - a2)*scrh;
365  alpha_f = MAX(0.0, alpha_f);
366  alpha_s = MAX(0.0, alpha_s);
367  alpha_f = sqrt(alpha_f);
368  alpha_s = sqrt(alpha_s);
369  }
370 
371  if (Btmag > 1.e-9) {
372  SELECT( ,
373  beta_y = DSIGN(By); ,
374  beta_y = By/Btmag;
375  beta_z = Bz/Btmag;)
376  } else {
377  SELECT( ,
378  beta_y = 1.0; ,
379  beta_z = beta_y = sqrt_1_2;)
380  }
381 
382  /* -------------------------------------------------------------------
383  6f. Compute non-zero entries of conservative eigenvectors (Rc),
384  wave strength L*dU (=eta) for all 8 (or 7) waves using the
385  expressions given by Eq. [4.18]--[4.21].
386  Fast and slow eigenvectors are multiplied by a^2 while
387  jumps are divided by a^2.
388 
389  Notes:
390  - the expression on the paper has a typo in the very last term
391  of the energy component: it should be + and not - !
392  - with background field splitting: additional terms must be
393  added to the energy component for fast, slow and Alfven waves.
394  To obtain energy element, conservative eigenvector (with
395  total field) must be multiplied by | 0 0 0 0 -B0y -B0z 1 |.
396  Also, H - b2 does not give gas enthalpy. A term b0*btot must
397  be added and eta (wave strength) should contain total field
398  and deviation's delta.
399  ------------------------------------------------------------------- */
400 
401  /* -----------------------
402  Fast wave: u - c_f
403  ----------------------- */
404 
405  k = KFASTM;
406  lambda[k] = u - cf;
407 
408  scrh = alpha_s*cs*sBx;
409  beta_dv = EXPAND(0.0, + beta_y*dV[VXt], + beta_z*dV[VXb]);
410  beta_dB = EXPAND(0.0, + beta_y*dV[BXt], + beta_z*dV[BXb]);
411  beta_v = EXPAND(0.0, + beta_y*v, + beta_z*w);
412 
413  Rc[RHO][k] = alpha_f;
414  EXPAND(Rc[MXn][k] = alpha_f*lambda[k]; ,
415  Rc[MXt][k] = alpha_f*v + scrh*beta_y; ,
416  Rc[MXb][k] = alpha_f*w + scrh*beta_z;)
417  EXPAND( ,
418  Rc[BXt][k] = alpha_s*a*beta_y/sqrt_rho; ,
419  Rc[BXb][k] = alpha_s*a*beta_z/sqrt_rho;)
420 
421  #if EOS == IDEAL
422  Rc[ENG][k] = alpha_f*(Hgas - u*cf) + scrh*beta_v
423  + alpha_s*a*Btmag/sqrt_rho;
424  #if BACKGROUND_FIELD == YES
425  Rc[ENG][k] -= B0y*Rc[BXt][k] + B0z*Rc[BXb][k];
426  #endif
427 
428  eta[k] = alpha_f*(X*dV[RHO] + dV[PRS]) + rho*scrh*beta_dv
429  - rho*alpha_f*cf*dV[VXn] + sqrt_rho*alpha_s*a*beta_dB;
430  #elif EOS == ISOTHERMAL
431  eta[k] = alpha_f*(0.0*X + a2)*dV[RHO] + rho*scrh*beta_dv
432  - rho*alpha_f*cf*dV[VXn] + sqrt_rho*alpha_s*a*beta_dB;
433  #endif
434 
435  eta[k] *= 0.5/a2;
436 
437  /* -----------------------
438  Fast wave: u + c_f
439  ----------------------- */
440 
441  k = KFASTP;
442  lambda[k] = u + cf;
443 
444  Rc[RHO][k] = alpha_f;
445  EXPAND( Rc[MXn][k] = alpha_f*lambda[k]; ,
446  Rc[MXt][k] = alpha_f*v - scrh*beta_y; ,
447  Rc[MXb][k] = alpha_f*w - scrh*beta_z; )
448  EXPAND( ,
449  Rc[BXt][k] = Rc[BXt][KFASTM]; ,
450  Rc[BXb][k] = Rc[BXb][KFASTM]; )
451 
452  #if EOS == IDEAL
453  Rc[ENG][k] = alpha_f*(Hgas + u*cf) - scrh*beta_v
454  + alpha_s*a*Btmag/sqrt_rho;
455 
456  #if BACKGROUND_FIELD == YES
457  Rc[ENG][k] -= B0y*Rc[BXt][k] + B0z*Rc[BXb][k];
458  #endif
459 
460  eta[k] = alpha_f*(X*dV[RHO] + dV[PRS]) - rho*scrh*beta_dv
461  + rho*alpha_f*cf*dV[VXn] + sqrt_rho*alpha_s*a*beta_dB;
462  #elif EOS == ISOTHERMAL
463  eta[k] = alpha_f*(0.*X + a2)*dV[RHO] - rho*scrh*beta_dv
464  + rho*alpha_f*cf*dV[VXn] + sqrt_rho*alpha_s*a*beta_dB;
465  #endif
466 
467  eta[k] *= 0.5/a2;
468 
469  /* -----------------------
470  Entropy wave: u
471  ----------------------- */
472 
473  #if EOS == IDEAL
474  k = KENTRP;
475  lambda[k] = u;
476 
477  Rc[RHO][k] = 1.0;
478  EXPAND( Rc[MXn][k] = u; ,
479  Rc[MXt][k] = v; ,
480  Rc[MXb][k] = w; )
481  Rc[ENG][k] = 0.5*vel2 + (g_gamma - 2.0)/g1*X;
482 
483  eta[k] = ((a2 - X)*dV[RHO] - dV[PRS])/a2;
484  #endif
485 
486  /* -----------------------------------------------------------------
487  div.B wave (u): this wave exists when:
488 
489  1) 8 wave formulation
490  2) CT, since we always have 8 components, but it
491  carries zero jump.
492 
493  With GLM, KDIVB is replaced by KPSI_GLMM, KPSI_GLMP and these
494  two waves should not enter in the Riemann solver (eta = 0.0)
495  since the 2x2 linear system formed by (B,psi) has already
496  been solved.
497  ----------------------------------------------------------------- */
498 
499  #ifdef GLM_MHD
500  lambda[KPSI_GLMP] = glm_ch;
501  lambda[KPSI_GLMM] = -glm_ch;
502  eta[KPSI_GLMP] = eta[KPSI_GLMM] = 0.0;
503  #else
504  k = KDIVB;
505  lambda[k] = u;
506  #if DIVB_CONTROL == EIGHT_WAVES
507  Rc[BXn][k] = 1.0;
508  eta[k] = dU[BXn];
509  #else
510  Rc[BXn][k] = eta[k] = 0.0;
511  #endif
512  #endif
513 
514  #if COMPONENTS > 1
515 
516  /* -----------------------
517  Slow wave: u - c_s
518  ----------------------- */
519 
520  scrh = alpha_f*cf*sBx;
521 
522  k = KSLOWM;
523  lambda[k] = u - cs;
524 
525  Rc[RHO][k] = alpha_s;
526  EXPAND( Rc[MXn][k] = alpha_s*lambda[k]; ,
527  Rc[MXt][k] = alpha_s*v - scrh*beta_y; ,
528  Rc[MXb][k] = alpha_s*w - scrh*beta_z;)
529  EXPAND( ,
530  Rc[BXt][k] = - alpha_f*a*beta_y/sqrt_rho; ,
531  Rc[BXb][k] = - alpha_f*a*beta_z/sqrt_rho; )
532 
533  #if EOS == IDEAL
534  Rc[ENG][k] = alpha_s*(Hgas - u*cs) - scrh*beta_v
535  - alpha_f*a*Btmag/sqrt_rho;
536  #if BACKGROUND_FIELD == YES
537  Rc[ENG][k] -= B0y*Rc[BXt][k] + B0z*Rc[BXb][k];
538  #endif
539 
540  eta[k] = alpha_s*(X*dV[RHO] + dV[PRS]) - rho*scrh*beta_dv
541  - rho*alpha_s*cs*dV[VXn] - sqrt_rho*alpha_f*a*beta_dB;
542  #elif EOS == ISOTHERMAL
543  eta[k] = alpha_s*(0.*X + a2)*dV[RHO] - rho*scrh*beta_dv
544  - rho*alpha_s*cs*dV[VXn] - sqrt_rho*alpha_f*a*beta_dB;
545  #endif
546 
547  eta[k] *= 0.5/a2;
548 
549  /* -----------------------
550  Slow wave: u + c_s
551  ----------------------- */
552 
553  k = KSLOWP;
554  lambda[k] = u + cs;
555 
556  Rc[RHO][k] = alpha_s;
557  EXPAND(Rc[MXn][k] = alpha_s*lambda[k]; ,
558  Rc[MXt][k] = alpha_s*v + scrh*beta_y; ,
559  Rc[MXb][k] = alpha_s*w + scrh*beta_z; )
560  EXPAND( ,
561  Rc[BXt][k] = Rc[BXt][KSLOWM]; ,
562  Rc[BXb][k] = Rc[BXb][KSLOWM];)
563 
564  #if EOS == IDEAL
565  Rc[ENG][k] = alpha_s*(Hgas + u*cs) + scrh*beta_v
566  - alpha_f*a*Btmag/sqrt_rho;
567  #if BACKGROUND_FIELD == YES
568  Rc[ENG][k] -= B0y*Rc[BXt][k] + B0z*Rc[BXb][k];
569  #endif
570 
571  eta[k] = alpha_s*(X*dV[RHO] + dV[PRS]) + rho*scrh*beta_dv
572  + rho*alpha_s*cs*dV[VXn] - sqrt_rho*alpha_f*a*beta_dB;
573  #elif EOS == ISOTHERMAL
574  eta[k] = alpha_s*(0.*X + a2)*dV[RHO] + rho*scrh*beta_dv
575  + rho*alpha_s*cs*dV[VXn] - sqrt_rho*alpha_f*a*beta_dB;
576  #endif
577 
578  eta[k] *= 0.5/a2;
579 
580  #endif
581 
582  #if COMPONENTS == 3
583 
584  /* ------------------------
585  Alfven wave: u - c_a
586  ------------------------ */
587 
588  k = KALFVM;
589  lambda[k] = u - ca;
590 
591  Rc[MXt][k] = - rho*beta_z;
592  Rc[MXb][k] = + rho*beta_y;
593  Rc[BXt][k] = - sBx*sqrt_rho*beta_z;
594  Rc[BXb][k] = sBx*sqrt_rho*beta_y;
595  #if EOS == IDEAL
596  Rc[ENG][k] = - rho*(v*beta_z - w*beta_y);
597  #if BACKGROUND_FIELD == YES
598  Rc[ENG][k] -= B0y*Rc[BXt][k] + B0z*Rc[BXb][k];
599  #endif
600  #endif
601 
602  eta[k] = + beta_y*dV[VXb] - beta_z*dV[VXt]
603  + sBx/sqrt_rho*(beta_y*dV[BXb] - beta_z*dV[BXt]);
604 
605  eta[k] *= 0.5;
606 
607  /* -----------------------
608  Alfven wave: u + c_a
609  ----------------------- */
610 
611  k = KALFVP;
612  lambda[k] = u + ca;
613 
614  Rc[MXt][k] = - Rc[MXt][KALFVM];
615  Rc[MXb][k] = - Rc[MXb][KALFVM];
616  Rc[BXt][k] = Rc[BXt][KALFVM];
617  Rc[BXb][k] = Rc[BXb][KALFVM];
618  #if EOS == IDEAL
619  Rc[ENG][k] = - Rc[ENG][KALFVM];
620  #if BACKGROUND_FIELD == YES
621  Rc[ENG][k] -= B0y*Rc[BXt][k] + B0z*Rc[BXb][k];
622  #endif
623  #endif
624 
625  eta[k] = - beta_y*dV[VXb] + beta_z*dV[VXt]
626  + sBx/sqrt_rho*(beta_y*dV[BXb] - beta_z*dV[BXt]);
627 
628  eta[k] *= 0.5;
629  #endif
630 
631  /* -----------------------------------------
632  6g. Compute maximum signal velocity
633  ----------------------------------------- */
634 
635  cmax[i] = fabs (u) + cf;
636  g_maxMach = MAX (fabs (u / a), g_maxMach);
637  for (k = 0; k < NFLX; k++) alambda[k] = fabs(lambda[k]);
638 
639  /* --------------------------------
640  6h. Entropy Fix
641  -------------------------------- */
642 
643  if (alambda[KFASTM] < 0.5*delta) {
644  alambda[KFASTM] = lambda[KFASTM]*lambda[KFASTM]/delta + 0.25*delta;
645  }
646  if (alambda[KFASTP] < 0.5*delta) {
647  alambda[KFASTP] = lambda[KFASTP]*lambda[KFASTP]/delta + 0.25*delta;
648  }
649  #if COMPONENTS > 1
650  if (alambda[KSLOWM] < 0.5*delta) {
651  alambda[KSLOWM] = lambda[KSLOWM]*lambda[KSLOWM]/delta + 0.25*delta;
652  }
653  if (alambda[KSLOWP] < 0.5*delta) {
654  alambda[KSLOWP] = lambda[KSLOWP]*lambda[KSLOWP]/delta + 0.25*delta;
655  }
656  #endif
657 
658  /* ---------------------------------
659  6i. Compute Roe numerical flux
660  --------------------------------- */
661 
662  for (nv = 0; nv < NFLX; nv++) {
663  scrh = 0.0;
664  for (k = 0; k < NFLX; k++) {
665  scrh += alambda[k]*eta[k]*Rc[nv][k];
666  }
667  state->flux[i][nv] = 0.5*(fL[i][nv] + fR[i][nv] - scrh);
668  }
669  state->press[i] = 0.5*(pL[i] + pR[i]);
670 
671  /* --------------------------------------------------------
672  6j. Check the Roe matrix condition, FR - FL = A*(UR - UL)
673  where A*(UR - UL) = R*lambda*eta.
674  -------------------------------------------------------- */
675 
676  #if CHECK_ROE_MATRIX == YES
677  for (nv = 0; nv < NFLX; nv++){
678  dV[nv] = fR[i][nv] - fL[i][nv];
679  if (nv == MXn) dV[MXn] += pR[i] - pL[i];
680  for (k = 0; k < NFLX; k++){
681  dV[nv] -= Rc[nv][k]*eta[k]*lambda[k];
682  }
683  if (fabs(dV[nv]) > 1.e-4){
684  printf (" ! Roe matrix condition not satisfied, var = %d\n", nv);
685  printf (" ! Err = %12.6e\n",dV[nv]);
686  Show(VL, i);
687  Show(VR, i);
688  exit(1);
689  }
690  }
691  #endif
692 
693  /* -------------------------------------------------------------
694  6k. Save max and min Riemann fan speeds for EMF computation.
695  ------------------------------------------------------------- */
696 
697  SL[i] = lambda[KFASTM];
698  SR[i] = lambda[KFASTP];
699 
700  /* -----------------------------------------------------------------
701  6l. Hybridize with HLL solver: replace occurences of unphysical
702  states (p < 0, rho < 0) with HLL Flux. Reference:
703 
704  "A Positive Conservative Method for MHD based based on HLL
705  and Roe methods", P. Janhunen, JCP (2000), 160, 649
706  ----------------------------------------------------------------- */
707 
708  #if HLL_HYBRIDIZATION == YES
709 
710  if (SL[i] < 0.0 && SR[i] > 0.0){
711  ifail = 0;
712 
713  /* -----------------------
714  check left state
715  ----------------------- */
716 
717  #if EOS == ISOTHERMAL
718  Uv[RHO] = uL[RHO] + (state->flux[i][RHO] - fL[i][RHO])/SL[i];
719  ifail = (Uv[RHO] < 0.0);
720  #else
721  for (nv = NFLX; nv--; ){
722  Uv[nv] = uL[nv] + (state->flux[i][nv] - fL[i][nv])/SL[i];
723  }
724  Uv[MXn] += (state->press[i] - pL[i])/SL[i];
725 
726  vel2 = EXPAND(Uv[MX1]*Uv[MX1], + Uv[MX2]*Uv[MX2], + Uv[MX3]*Uv[MX3]);
727  b2 = EXPAND(Uv[BX1]*Uv[BX1], + Uv[BX2]*Uv[BX2], + Uv[BX3]*Uv[BX3]);
728  pr = Uv[ENG] - 0.5*vel2/Uv[RHO] - 0.5*b2;
729  ifail = (pr < 0.0) || (Uv[RHO] < 0.0);
730  #endif
731 
732  /* -----------------------
733  check right state
734  ----------------------- */
735 
736  #if EOS == ISOTHERMAL
737  Uv[RHO] = uR[RHO] + (state->flux[i][RHO] - fR[i][RHO])/SR[i];
738  ifail = (Uv[RHO] < 0.0);
739  #else
740  for (nv = NFLX; nv--; ){
741  Uv[nv] = uR[nv] + (state->flux[i][nv] - fR[i][nv])/SR[i];
742  }
743  Uv[MXn] += (state->press[i] - pR[i])/SR[i];
744 
745  vel2 = EXPAND(Uv[MX1]*Uv[MX1], + Uv[MX2]*Uv[MX2], + Uv[MX3]*Uv[MX3]);
746  b2 = EXPAND(Uv[BX1]*Uv[BX1], + Uv[BX2]*Uv[BX2], + Uv[BX3]*Uv[BX3]);
747  pr = Uv[ENG] - 0.5*vel2/Uv[RHO] - 0.5*b2;
748  ifail = (pr < 0.0) || (Uv[RHO] < 0.0);
749  #endif
750 
751  /* -------------------------------------------------------
752  Use the HLL flux function if the interface lies
753  within a strong shock. The effect of this switch is
754  visible in the Mach reflection test.
755  ------------------------------------------------------- */
756 
757  #if DIMENSIONS > 1
758  #if EOS == ISOTHERMAL
759  scrh = fabs(vL[RHO] - vR[RHO]);
760  scrh /= MIN(vL[RHO], vR[RHO]);
761  #else
762  scrh = fabs(vL[PRS] - vR[PRS]);
763  scrh /= MIN(vL[PRS], vR[PRS]);
764  #endif
765  if (scrh > 1.0 && (vR[VXn] < vL[VXn])) ifail = 1;
766  #endif
767 
768  if (ifail){
769  scrh = 1.0/(SR[i] - SL[i]);
770  for (nv = 0; nv < NFLX; nv++) {
771  state->flux[i][nv] = SL[i]*SR[i]*(uR[nv] - uL[nv]) +
772  SR[i]*fL[i][nv] - SL[i]*fR[i][nv];
773  state->flux[i][nv] *= scrh;
774  }
775  state->press[i] = (SR[i]*pL[i] - SL[i]*pR[i])*scrh;
776  }
777  }
778  #endif
779  }
780 
781 /* --------------------------------------------------------
782  initialize source term
783  -------------------------------------------------------- */
784 
785  #if DIVB_CONTROL == EIGHT_WAVES
786  Roe_DivBSource (state, beg + 1, end, grid);
787  #endif
788 
789 }
IDEAL
#define IDEAL
Definition: pluto.h:45
MX3
#define MX3
Definition: mod_defs.h:22
EOS
#define EOS
Definition: pluto.h:341
a
static double a
Definition: init.c:135
MAX
#define MAX(a, b)
Definition: macros.h:101
g_gamma
double g_gamma
Definition: globals.h:112
GLM_Solve
void GLM_Solve(const State_1D *state, double **VL, double **VR, int beg, int end, Grid *grid)
Definition: glm.c:24
MX1
#define MX1
Definition: mod_defs.h:20
STATE_1D::flux
double ** flux
upwind flux computed with the Riemann solver
Definition: structs.h:149
print1
void print1(const char *fmt,...)
Definition: amrPluto.cpp:511
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
eta
static double *** eta[3]
Definition: res_functions.c:94
ca
static double ca
Definition: init.c:75
STATE_1D::SR
double * SR
Rightmost velocity in the Riemann fan at i+1/2.
Definition: structs.h:167
func_param::v
double v[NVAR]
Definition: eos.h:106
DSIGN
#define DSIGN(x)
Definition: macros.h:110
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
KPSI_GLMP
Definition: mod_defs.h:75
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
MIN
#define MIN(a, b)
Definition: macros.h:104
MXt
int MXt
Definition: globals.h:74
g_maxMach
double g_maxMach
The maximum Mach number computed during integration.
Definition: globals.h:119
X
#define X
Definition: radiat.modified_by_Bhargav.c:3
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
Roe_DivBSource
void Roe_DivBSource(const State_1D *state, int is, int ie, Grid *grid)
Definition: source.c:39
j
int j
Definition: analysis.c:2
MX2
#define MX2
Definition: mod_defs.h:21
Sph_disk.vel2
list vel2
Definition: Sph_disk.py:39
VXt
int VXt
Definition: globals.h:73
STATE_1D::flag
unsigned char * flag
Definition: structs.h:168
KFASTP
Definition: mod_defs.h:58
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
MXn
int MXn
Definition: globals.h:74
KFASTM
Definition: mod_defs.h:58
sqrt_1_2
#define sqrt_1_2
Definition: roe.c:48
VXn
int VXn
Definition: globals.h:73
STATE_1D::uR
double ** uR
same as vR, in conservative vars
Definition: structs.h:145
BX3
#define BX3
Definition: mod_defs.h:27
ARRAY_1D
#define ARRAY_1D(nx, type)
Definition: prototypes.h:170
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
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
BX1
#define BX1
Definition: mod_defs.h:25
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
BX2
#define BX2
Definition: mod_defs.h:26
glm_ch
double glm_ch
The propagation speed of divergence error.
Definition: glm.c:21
BXt
int BXt
Definition: globals.h:75
KENTRP
Definition: mod_defs.h:59
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
KPSI_GLMM
Definition: mod_defs.h:75
STATE_1D::uL
double ** uL
same as vL, in conservative vars
Definition: structs.h:144

Here is the call graph for this function: