PLUTO
Macros | Functions
roe.c File Reference

Roe Riemann solver for the HD equations. More...

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

Go to the source code of this file.

Macros

#define ROE_AVERAGE   YES
 
#define CHECK_ROE_MATRIX   NO
 

Functions

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

Detailed Description

Roe Riemann solver for the HD equations.

Solve the Riemann problem for the Euler equations of gasdynamics using the standard Roe solver with local characteristic decomposition. Eigenvectors are identical to the ones given in the book by Toro and can also be derived from the maple script "eigenv.maple" in Src/HD/. The solver can be used for adiabatic and isothermal hydrodynamics.

The macro ROE_AVERAGE specifies how the averaging process is done:

  • ROE_AVERAGE == YES use Roe average (default);
  • ROE_AVERAGE == NO use arithmetic average.

The macro CHECK_ROE_MATRIX can be used to verify that the characteristic decomposition reproduces the Roe matrix.

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:

  • "Riemann Solver and Numerical Methods for Fluid Dynamics" by E.F. Toro (Chapter 11)
Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
Dec 10, 2013

Definition in file roe.c.

Macro Definition Documentation

#define CHECK_ROE_MATRIX   NO

Definition at line 39 of file roe.c.

#define ROE_AVERAGE   YES

Definition at line 38 of file roe.c.

Function Documentation

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

Solve the Riemann problem using the Roe 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 42 of file roe.c.

54 {
55  int nv, i, j, k, nn;
56  double scrh;
57  double um[NFLX], vel2;
58  double a2, a, h;
59  double dv[NFLX], eta[NFLX];
60  double Rc[NFLX][NFLX], alambda[NFLX], lambda[NFLX];
61  double delta, delta_inv, gmm1, gmm1_inv;
62 #if ROE_AVERAGE == YES
63  double s, c, hl, hr;
64 #endif
65  double *ql, *qr, *uL, *uR;
66  static double **fL, **fR, *pL, *pR, *a2L, *a2R;
67 
68  double bmin, bmax, scrh1;
69  double Us[NFLX];
70 
71  delta = 1.e-7;
72  delta_inv = 1.0/delta;
73  #if EOS == IDEAL
74  gmm1 = g_gamma - 1.0;
75  gmm1_inv = 1.0/gmm1;
76  #endif
77 
78  if (fL == NULL){
79  fL = ARRAY_2D(NMAX_POINT, NFLX, double);
80  fR = ARRAY_2D(NMAX_POINT, NFLX, double);
81  pR = ARRAY_1D(NMAX_POINT, double);
82  pL = ARRAY_1D(NMAX_POINT, double);
83  a2R = ARRAY_1D(NMAX_POINT, double);
84  a2L = ARRAY_1D(NMAX_POINT, double);
85  }
86 
87  for (i = NFLX; i--; ) {
88  for (j = NFLX; j--; ) {
89  Rc[i][j] = 0.0;
90  }}
91 
92 /* ----------------------------------------------------
93  compute sound speed & fluxes at zone interfaces
94  ---------------------------------------------------- */
95 
96  SoundSpeed2 (state->vL, a2L, NULL, beg, end, FACE_CENTER, grid);
97  SoundSpeed2 (state->vR, a2R, NULL, beg, end, FACE_CENTER, grid);
98 
99  Flux (state->uL, state->vL, a2L, fL, pL, beg, end);
100  Flux (state->uR, state->vR, a2R, fR, pR, beg, end);
101 
102  for (i = beg; i <= end; i++) {
103 
104  uR = state->uR[i];
105  uL = state->uL[i];
106 
107 #if SHOCK_FLATTENING == MULTID
108 
109  /* ---------------------------------------------
110  HLL switching function as in Quirk (1994).
111  Since the problem is related to multidimensional
112  pathologies, it works in more than 1-D only.
113  Use the HLL flux function if the interface
114  lies within a strong shock.
115  The effect of this switch is visible
116  in the Mach reflection test.
117  --------------------------------------------- */
118 
119  if ((state->flag[i] & FLAG_HLL) || (state->flag[i+1] & FLAG_HLL)){
120  HLL_Speed (state->vL, state->vR, a2L, a2R,
121  &bmin - i, &bmax - i, i, i);
122  a = MAX(fabs(bmin), fabs(bmax));
123  cmax[i] = a;
124  bmin = MIN(0.0, bmin);
125  bmax = MAX(0.0, bmax);
126  scrh = 1.0/(bmax - bmin);
127  for (nv = NFLX; nv--; ){
128  state->flux[i][nv] = bmin*bmax*(uR[nv] - uL[nv])
129  + bmax*fL[i][nv] - bmin*fR[i][nv];
130  state->flux[i][nv] *= scrh;
131  }
132  state->press[i] = (bmax*pL[i] - bmin*pR[i])*scrh;
133  continue;
134  }
135 #endif
136 
137  ql = state->vL[i];
138  qr = state->vR[i];
139 
140  /* ---- Define Wave Jumps ---- */
141 
142  for (nv = NFLX; nv--; ) dv[nv] = qr[nv] - ql[nv];
143 
144  #if ROE_AVERAGE == YES
145  s = sqrt(qr[RHO]/ql[RHO]);
146  um[RHO] = ql[RHO]*s;
147  s = 1.0/(1.0 + s);
148  c = 1.0 - s;
149 
150  EXPAND(um[VX1] = s*ql[VX1] + c*qr[VX1]; ,
151  um[VX2] = s*ql[VX2] + c*qr[VX2]; ,
152  um[VX3] = s*ql[VX3] + c*qr[VX3];)
153 
154  #if EOS == IDEAL
155  vel2 = EXPAND(um[VX1]*um[VX1], + um[VX2]*um[VX2], + um[VX3]*um[VX3]);
156 
157  hl = 0.5*(EXPAND(ql[VX1]*ql[VX1], + ql[VX2]*ql[VX2], + ql[VX3]*ql[VX3]));
158  hl += a2L[i]*gmm1_inv;
159 
160  hr = 0.5*(EXPAND(qr[VX1]*qr[VX1], + qr[VX2]*qr[VX2], + qr[VX3]*qr[VX3]));
161  hr += a2R[i]*gmm1_inv;
162 
163  h = s*hl + c*hr;
164 
165  /* -------------------------------------------------
166  the following should be equivalent to
167 
168  scrh = EXPAND( dv[VX1]*dv[VX1],
169  + dv[VX2]*dv[VX2],
170  + dv[VX3]*dv[VX3]);
171 
172  a2 = s*a2L + c*a2R + 0.5*gmm1*s*c*scrh;
173 
174  and therefore always positive.
175  just work out the coefficiendnts...
176  -------------------------------------------------- */
177 
178  a2 = gmm1*(h - 0.5*vel2);
179  a = sqrt(a2);
180  #endif
181  #else
182  for (nv = NFLX; nv--; ) um[nv] = 0.5*(ql[nv] + qr[nv]);
183  #if EOS == IDEAL
184  a2 = g_gamma*um[PRS]/um[RHO];
185  a = sqrt(a2);
186 
187  vel2 = EXPAND(um[VX1]*um[VX1], + um[VX2]*um[VX2], + um[VX3]*um[VX3]);
188  h = 0.5*vel2 + a2/gmm1;
189  #endif /* EOS == IDEAL */
190  #endif /* ROE_AVERAGE == YES/NO */
191 
192  #if EOS == ISOTHERMAL
193  a2 = 0.5*(a2L[i] + a2R[i]);
194  a = sqrt(a2);
195  #endif
196 
197  /* ----------------------------------------------------------------
198  define non-zero components of conservative eigenvectors Rc,
199  eigenvalues (lambda) and wave strenght eta = L.du
200  ---------------------------------------------------------------- */
201 
202  /* ---- (u - c_s) ---- */
203 
204  nn = 0;
205  lambda[nn] = um[VXn] - a;
206  #if EOS == IDEAL
207  eta[nn] = 0.5/a2*(dv[PRS] - dv[VXn]*um[RHO]*a);
208  #elif EOS == ISOTHERMAL
209  eta[nn] = 0.5*(dv[RHO] - um[RHO]*dv[VXn]/a);
210  #endif
211 
212  Rc[RHO][nn] = 1.0;
213  EXPAND(Rc[MXn][nn] = um[VXn] - a; ,
214  Rc[MXt][nn] = um[VXt]; ,
215  Rc[MXb][nn] = um[VXb];)
216  #if EOS == IDEAL
217  Rc[ENG][nn] = h - um[VXn]*a;
218  #endif
219 
220  /* ---- (u + c_s) ---- */
221 
222  nn = 1;
223  lambda[nn] = um[VXn] + a;
224  #if EOS == IDEAL
225  eta[nn] = 0.5/a2*(dv[PRS] + dv[VXn]*um[RHO]*a);
226  #elif EOS == ISOTHERMAL
227  eta[nn] = 0.5*(dv[RHO] + um[RHO]*dv[VXn]/a);
228  #endif
229 
230  Rc[RHO][nn] = 1.0;
231  EXPAND(Rc[MXn][nn] = um[VXn] + a; ,
232  Rc[MXt][nn] = um[VXt]; ,
233  Rc[MXb][nn] = um[VXb];)
234  #if EOS == IDEAL
235  Rc[ENG][nn] = h + um[VXn]*a;
236  #endif
237 
238  /* ---- (u) ---- */
239 
240  #if EOS == IDEAL
241  nn = 2;
242  lambda[nn] = um[VXn];
243  eta[nn] = dv[RHO] - dv[PRS]/a2;
244  Rc[RHO][nn] = 1.0;
245  EXPAND(Rc[MX1][nn] = um[VX1]; ,
246  Rc[MX2][nn] = um[VX2]; ,
247  Rc[MX3][nn] = um[VX3];)
248  Rc[ENG][nn] = 0.5*vel2;
249  #endif
250 
251  #if COMPONENTS > 1
252 
253  /* ---- (u) ---- */
254 
255  nn++;
256  lambda[nn] = um[VXn];
257  eta[nn] = um[RHO]*dv[VXt];
258  Rc[MXt][nn] = 1.0;
259  #if EOS == IDEAL
260  Rc[ENG][nn] = um[VXt];
261  #endif
262  #endif
263 
264  #if COMPONENTS > 2
265 
266  /* ---- (u) ---- */
267 
268  nn++;
269  lambda[nn] = um[VXn];
270  eta[nn] = um[RHO]*dv[VXb];
271  Rc[MXb][nn] = 1.0;
272  #if EOS == IDEAL
273  Rc[ENG][nn] = um[VXb];
274  #endif
275  #endif
276 
277  /* ---- get max eigenvalue ---- */
278 
279  cmax[i] = fabs(um[VXn]) + a;
280  g_maxMach = MAX(fabs(um[VXn]/a), g_maxMach);
281 
282  #if DIMENSIONS > 1
283 
284  /* ---------------------------------------------
285  use the HLL flux function if the interface
286  lies within a strong shock.
287  The effect of this switch is visible
288  in the Mach reflection test.
289  --------------------------------------------- */
290 
291  #if EOS == IDEAL
292  scrh = fabs(ql[PRS] - qr[PRS]);
293  scrh /= MIN(ql[PRS],qr[PRS]);
294  #elif EOS == ISOTHERMAL
295  scrh = fabs(ql[RHO] - qr[RHO]);
296  scrh /= MIN(ql[RHO],qr[RHO]);
297  scrh *= a*a;
298  #endif
299  if (scrh > 0.5 && (qr[VXn] < ql[VXn])){ /* -- tunable parameter -- */
300  bmin = MIN(0.0, lambda[0]);
301  bmax = MAX(0.0, lambda[1]);
302  scrh1 = 1.0/(bmax - bmin);
303  for (nv = NFLX; nv--; ){
304  state->flux[i][nv] = bmin*bmax*(uR[nv] - uL[nv])
305  + bmax*fL[i][nv] - bmin*fR[i][nv];
306  state->flux[i][nv] *= scrh1;
307  }
308  state->press[i] = (bmax*pL[i] - bmin*pR[i])*scrh1;
309  continue;
310  }
311  #endif
312 
313  #if CHECK_ROE_MATRIX == YES
314  for (nv = 0; nv < NFLX; nv++){
315  um[nv] = 0.0;
316  for (k = 0; k < NFLX; k++){
317  for (j = 0; j < NFLX; j++){
318  um[nv] += Rc[nv][k]*(k==j)*lambda[k]*eta[j];
319  }}
320  }
321  for (nv = 0; nv < NFLX; nv++){
322  scrh = fR[i][nv] - fL[i][nv] - um[nv];
323  if (nv == MXn) scrh += pR[i] - pL[i];
324  if (fabs(scrh) > 1.e-6){
325  print ("! Matrix condition not satisfied %d, %12.6e\n", nv, scrh);
326  Show(state->vL, i);
327  Show(state->vR, i);
328  exit(1);
329  }
330  }
331  #endif
332 
333  /* -----------------------------------------------------------
334  compute Roe flux
335  ----------------------------------------------------------- */
336 
337  for (nv = NFLX; nv--; ) alambda[nv] = fabs(lambda[nv]);
338 
339  /* ---- entropy fix ---- */
340 
341  if (alambda[0] <= delta) {
342  alambda[0] = 0.5*lambda[0]*lambda[0]/delta + 0.5*delta;
343  }
344  if (alambda[1] <= delta) {
345  alambda[1] = 0.5*lambda[1]*lambda[1]/delta + 0.5*delta;
346  }
347 
348  for (nv = NFLX; nv--; ) {
349  state->flux[i][nv] = fL[i][nv] + fR[i][nv];
350  for (k = NFLX; k-- ; ) {
351  state->flux[i][nv] -= alambda[k]*eta[k]*Rc[nv][k];
352  }
353  state->flux[i][nv] *= 0.5;
354  }
355  state->press[i] = 0.5*(pL[i] + pR[i]);
356  }
357 }
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
MX1
#define MX1
Definition: mod_defs.h:20
VX2
#define VX2
Definition: mod_defs.h:29
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
eta
static double *** eta[3]
Definition: res_functions.c:94
STATE_1D::vR
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
VX1
#define VX1
Definition: mod_defs.h:28
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
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
NFLX
#define NFLX
Definition: mod_defs.h:32
FACE_CENTER
#define FACE_CENTER
Definition: pluto.h:206
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
k
int k
Definition: analysis.c:2
print
void print(const char *fmt,...)
Definition: amrPluto.cpp:497
MXn
int MXn
Definition: globals.h:74
menu.c
tuple c
Definition: menu.py:375
VXn
int VXn
Definition: globals.h:73
STATE_1D::uR
double ** uR
same as vR, in conservative vars
Definition: structs.h:145
s
#define s
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
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
VX3
#define VX3
Definition: mod_defs.h:30
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
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: