PLUTO
Macros | Functions
roe.bckfld.c File Reference
#include "pluto.h"
Include dependency graph for roe.bckfld.c:

Go to the source code of this file.

Macros

#define sqrt_1_2   (0.70710678118654752440)
 

Functions

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

Macro Definition Documentation

#define sqrt_1_2   (0.70710678118654752440)

Definition at line 3 of file roe.bckfld.c.

Function Documentation

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

Definition at line 6 of file roe.bckfld.c.

13  {
14  int nv, i, j, k;
15  double gmm1;
16  double scrh0, scrh1, scrh2, scrh3, scrh4;
17  double rho, u, v, w, vel2, bx, by, bz, pr;
18  double b1x, b1y, b1z;
19  double a2, a, ca2, cf2, cs2;
20  double cs, ca, cf, b2, A2, At2;
21  double S;
22  double alpha_f, alpha_s, beta_y, beta_z;
23  double dw[NVAR], *vL, *vR, *uL, *uR;
24  double *SL, *SR;
25  double eta[NVAR], lambda[NVAR];
26 
27  double tau, sqrt_rho;
28  double delta, delta_inv, gmm1_inv;
29 
30  static double *a2L, *a2R, *pR, *pL;
31  static double **fL, **fR, **Rp, **Rc;
32  static double **VL, **VR, **UL, **UR;
33 
34  double **bgf;
35 
36  gmm1 = g_gamma - 1.0;
37  delta = 1.e-7;
38  delta_inv = 1.0 / delta;
39  gmm1_inv = 1.0 / gmm1;
40 
41  if (fL == NULL){
42  fL = ARRAY_2D(NMAX_POINT, NVAR, double);
43  fR = ARRAY_2D(NMAX_POINT, NVAR, double);
44 
45  a2R = ARRAY_1D(NMAX_POINT, double);
46  a2L = ARRAY_1D(NMAX_POINT, double);
47 
48  pR = ARRAY_1D(NMAX_POINT, double);
49  pL = ARRAY_1D(NMAX_POINT, double);
50  #ifdef GLM_MHD
51  VL = ARRAY_2D(NMAX_POINT, NVAR, double);
52  VR = ARRAY_2D(NMAX_POINT, NVAR, double);
53  UL = ARRAY_2D(NMAX_POINT, NVAR, double);
54  UR = ARRAY_2D(NMAX_POINT, NVAR, double);
55  #endif
56 
57  Rc = ARRAY_2D(NMAX_POINT, NVAR, double);
58  Rp = ARRAY_2D(NMAX_POINT, NVAR, double);
59 
60  }
61 
62  #if BACKGROUND_FIELD == YES
63  bgf = GetBackgroundField (beg, end, FACE_CENTER, grid);
64  #endif
65 
66  #ifdef GLM_MHD
67  GLM_Solve (state, VL, VR, beg, end, grid);
68  PrimToCons (VL, UL, beg, end);
69  PrimToCons (VR, UR, beg, end);
70  #else
71  VL = state->vL; UL = state->uL;
72  VR = state->vR; UR = state->uR;
73  #endif
74 
75  SoundSpeed2 (VL, a2L, NULL, beg, end, FACE_CENTER, grid);
76  SoundSpeed2 (VR, a2R, NULL, beg, end, FACE_CENTER, grid);
77 
78  Flux (UL, VL, a2L, bgf, fL, pL, beg, end);
79  Flux (UR, VR, a2R, bgf, fR, pR, beg, end);
80 
81  SL = state->SL; SR = state->SR;
82 
83 /* Some eigenvector components will always be zero so set
84  R = 0 initially */
85 
86  for (k = NFLX; k--; ) {
87  for (j = NFLX; j--; ) {
88  Rp[k][j] = Rc[k][j] = 0.0;
89  }}
90 
91 
92  for (i = beg; i <= end; i++) {
93 
94  vL = VL[i]; vR = VR[i];
95  uL = UL[i]; uR = UR[i];
96 
97  #if SHOCK_FLATTENING == MULTID
98 
99  /* -- revert to HLL in proximity of strong shock -- */
100 
101  if (CheckZone(i, FLAG_HLL) || CheckZone(i+1,FLAG_HLL)){
102  HLL_SPEED (VL, VR, NULL, SL, SR, i, i);
103 
104  scrh0 = MAX(fabs(SL[i]), fabs(SR[i]));
105  cmax[i] = scrh0;
106 
107  if (SL[i] > 0.0) {
108  for (nv = NFLX; nv--; ) state->flux[i][nv] = fL[i][nv];
109  state->press[i] = pL[i];
110  } else if (SR[i] < 0.0) {
111  for (nv = NFLX; nv--; ) state->flux[i][nv] = fR[i][nv];
112  state->press[i] = pR[i];
113  }else{
114  scrh0 = 1.0/(SR[i] - SL[i]);
115  for (nv = NFLX; nv--; ){
116  state->flux[i][nv] = SR[i]*SL[i]*(uR[nv] - uL[nv])
117  + SR[i]*fL[i][nv] - SL[i]*fR[i][nv];
118  state->flux[i][nv] *= scrh0;
119  }
120  state->press[i] = (SR[i]*pL[i] - SL[i]*pR[i])*scrh0;
121  }
122  continue;
123  }
124  #endif
125 
126  rho = 0.5*(vL[RHO] + vR[RHO]);
127  EXPAND (u = 0.5*(vL[VXn] + vR[VXn]); ,
128  v = 0.5*(vL[VXt] + vR[VXt]); ,
129  w = 0.5*(vL[VXb] + vR[VXb]); )
130  EXPAND (bx = b1x = 0.5*(vL[BXn] + vR[BXn]); ,
131  by = b1y = 0.5*(vL[BXt] + vR[BXt]); ,
132  bz = b1z = 0.5*(vL[BXb] + vR[BXb]);)
133 
134  #if BACKGROUND_FIELD == YES
135  EXPAND (bx += bgf[i][BXn]; ,
136  by += bgf[i][BXt]; ,
137  bz += bgf[i][BXb];)
138  #endif
139 
140  pr = 0.5*(vL[PRS] + vR[PRS]);
141 
142  for (nv = 0; nv < NFLX; nv++) dw[nv] = vR[nv] - vL[nv];
143 
144  tau = 1.0/rho;
145  sqrt_rho = sqrt(rho);
146 
147  vel2 = EXPAND(u*u, +v*v, +w*w);
148  a2 = g_gamma*pr*tau;
149 
150  scrh2 = bx*bx; /* > 0 */
151  scrh3 = SELECT(0.0, by*by, by*by + bz*bz); /* > 0 */
152 
153  b2 = scrh2 + scrh3; /* Magnetic field module, >0 */
154  ca2 = scrh2*tau; /* >0 if tau >0 */
155  A2 = b2*tau; /* >0 '' '' */
156  At2 = scrh3*tau; /* >0 */
157 
158  scrh1 = a2 - A2;
159  scrh0 = sqrt(scrh1*scrh1 + 4.0*a2*At2); /* >0 */
160 
161 /* Now get fast and slow speeds */
162 
163  cf2 = 0.5*(a2 + A2 + scrh0);
164  cs2 = a2*ca2/cf2;
165 
166  cf = sqrt(cf2); /* > 0 */
167  cs = sqrt(cs2); /* > 0 */
168  ca = sqrt(ca2); /* > 0 */
169  a = sqrt(a2); /* > 0 */
170 
171  scrh0 = 1.0/scrh0;
172  alpha_f = (a2 - cs2)*scrh0;
173  alpha_s = (cf2 - a2)*scrh0;
174 
175  alpha_f = MAX(0.0, alpha_f);
176  alpha_s = MAX(0.0, alpha_s);
177 
178  alpha_f = sqrt(alpha_f);
179  alpha_s = sqrt(alpha_s);
180 
181  scrh0 = sqrt(scrh3);
182 
183  if (scrh0 > 1.e-8) {
184  SELECT (, beta_y = DSIGN(by); ,
185  beta_y = by / scrh0; beta_z = bz / scrh0;)
186  } else {
187  SELECT (, beta_y = 1.0; ,
188  beta_z = beta_y = sqrt_1_2;)
189  }
190 
191  S = (bx >= 0.0 ? 1.0 : -1.0);
192 
193  /* --------------------------------------------------------
194  Define non-zero entries of primitive right
195  eigenvectors (Rp), wave strength l*dw (=eta) for all
196  8 (or 7) waves.
197  left eigenvectors for fast & slow waves can be defined
198  in terms of right eigenvectors (see page 296)
199  -------------------------------------------------------- */
200 
201  /* ---- FAST WAVE (u - c_f) ---- */
202 
203  k = KFASTM;
204  lambda[k] = u - cf;
205 
206  scrh0 = alpha_s*cs*S;
207  scrh1 = alpha_s*sqrt_rho*a;
208  scrh2 = 0.5 / a2;
209  scrh3 = scrh2*tau;
210  scrh4 = alpha_f*g_gamma*pr;
211 
212  Rp[RHO][k] = rho*alpha_f; /* right eigenvectors */
213  EXPAND(Rp[VXn][k] = -cf*alpha_f; ,
214  Rp[VXt][k] = scrh0*beta_y; ,
215  Rp[VXb][k] = scrh0*beta_z;)
216  EXPAND( ,
217  Rp[BXt][k] = scrh1*beta_y; ,
218  Rp[BXb][k] = scrh1*beta_z;)
219  Rp[PRS][k] = scrh4;
220 
221  eta[k] = (EXPAND( Rp[VXn][k]*dw[VXn],
222  + Rp[VXt][k]*dw[VXt],
223  + Rp[VXb][k]*dw[VXb]))*scrh2;
224  eta[k] += (EXPAND( 0.0 ,
225  + Rp[BXt][k]*dw[BXt] ,
226  + Rp[BXb][k]*dw[BXb])
227  + alpha_f*dw[PRS])*scrh3;
228 
229  /* ---- FAST WAVE (u + c_f) ---- */
230 
231  k = KFASTP;
232  lambda[k] = u + cf;
233 
234  Rp[RHO][k] = rho*alpha_f;
235  EXPAND(Rp[VXn][k] = cf*alpha_f; ,
236  Rp[VXt][k] = -scrh0*beta_y; ,
237  Rp[VXb][k] = -scrh0*beta_z;)
238  EXPAND(, Rp[BXt][k] = scrh1*beta_y; ,
239  Rp[BXb][k] = scrh1*beta_z;)
240  Rp[PRS][k] = scrh4;
241 
242  eta[k] = (EXPAND( Rp[VXn][k]*dw[VXn],
243  + Rp[VXt][k]*dw[VXt],
244  + Rp[VXb][k]*dw[VXb]))*scrh2;
245  eta[k] += (EXPAND( 0.0,
246  + Rp[BXt][k]*dw[BXt],
247  + Rp[BXb][k]*dw[BXb])
248  + alpha_f*dw[PRS])*scrh3;
249 
250  /* ---- ENTROPY WAVE (u) ---- */
251 
252  k = KENTRP;
253  lambda[k] = u;
254 
255  Rp[RHO][k] = 1.0;
256 
257  eta[k] = dw[RHO] - dw[PRS]/a2;
258 
259  /* ---- MAGNETIC FLUX, for 8-wave formulation only (u) ---- */
260 
261  #ifdef GLM_MHD
262  lambda[KPSI_GLMP] = glm_ch;
263  lambda[KPSI_GLMM] = -glm_ch;
264  eta[KPSI_GLMP] = eta[KPSI_GLMM] = 0.0;
265  #else
266  k = KDIVB;
267  lambda[k] = u;
268  #if DIVB_CONTROL == EIGHT_WAVES
269  Rc[BXn][k] = 1.0;
270  eta[k] = dU[BXn];
271  #else
272  Rc[BXn][k] = eta[k] = 0.0;
273  #endif
274  #endif
275 
276  #if COMPONENTS > 1 /* ---- SLOW WAVE (u - c_s) ---- */
277  k = KSLOWM;
278  lambda[k] = u - cs;
279 
280  scrh0 = alpha_f*cf*S;
281  scrh1 = alpha_f*sqrt_rho*a;
282  scrh4 = alpha_s*g_gamma*pr;
283 
284  Rp[RHO][k] = rho*alpha_s;
285  EXPAND(Rp[VXn][k] = -cs*alpha_s; ,
286  Rp[VXt][k] = -scrh0*beta_y; ,
287  Rp[VXb][k] = -scrh0*beta_z;)
288  EXPAND(, Rp[BXt][k] = -scrh1*beta_y; ,
289  Rp[BXb][k] = -scrh1*beta_z;)
290  Rp[PRS][k] = scrh4;
291 
292  eta[k] = (EXPAND( Rp[VXn][k]*dw[VXn],
293  + Rp[VXt][k]*dw[VXt],
294  + Rp[VXb][k]*dw[VXb]))*scrh2;
295  eta[k] +=(EXPAND( 0.0 ,
296  + Rp[BXt][k]*dw[BXt],
297  + Rp[BXb][k]*dw[BXb])
298  +alpha_s*dw[PRS])*scrh3;
299 
300  /* ---- SLOW WAVE (u + c_s) ---- */
301 
302  k = KSLOWP;
303  lambda[k] = u + cs;
304 
305  Rp[RHO][k] = rho*alpha_s;
306  EXPAND(Rp[VXn][k] = cs*alpha_s; ,
307  Rp[VXt][k] = scrh0*beta_y; ,
308  Rp[VXb][k] = scrh0*beta_z;)
309 
310  EXPAND ( , Rp[BXt][k] = -scrh1*beta_y; ,
311  Rp[BXb][k] = -scrh1*beta_z;)
312  Rp[PRS][k] = scrh4;
313 
314  eta[k] = (EXPAND( Rp[VXn][k]*dw[VXn],
315  + Rp[VXt][k]*dw[VXt],
316  + Rp[VXb][k]*dw[VXb]))*scrh2;
317  eta[k] +=(EXPAND( 0.0 ,
318  + Rp[BXt][k]*dw[BXt],
319  + Rp[BXb][k]*dw[BXb])
320  +alpha_s*dw[PRS])*scrh3;
321  #endif
322 
323  #if COMPONENTS == 3
324 
325  /* ---- ALFVEN WAVE (u-c_a) ---- */
326 
327  k = KALFVM;
328  lambda[k] = u - ca;
329 
330  scrh2 = beta_y*sqrt_1_2;
331  scrh3 = beta_z*sqrt_1_2;
332  Rp[VXt][k] = -scrh3;
333  Rp[VXb][k] = scrh2;
334  Rp[BXt][k] = -scrh3*sqrt_rho*S;
335  Rp[BXb][k] = scrh2*sqrt_rho*S;
336 
337  eta[k] = Rp[VXt][k]*dw[VXt] + Rp[VXb][k]*dw[VXb]
338  +(Rp[BXt][k]*dw[BXt] + Rp[BXb][k]*dw[BXb])*tau;
339 
340  /* ---- ALFVEN WAVE (u+c_a) ---- */
341 
342  k = KALFVP;
343  lambda[k] = u + ca;
344  Rp[VXt][k] = -scrh3;
345  Rp[VXb][k] = scrh2;
346  Rp[BXt][k] = scrh3*sqrt_rho*S;
347  Rp[BXb][k] = -scrh2*sqrt_rho*S;
348 
349  eta[k] = Rp[VXt][k]*dw[VXt] + Rp[VXb][k]*dw[VXb]
350  +(Rp[BXt][k]*dw[BXt] + Rp[BXb][k]*dw[BXb])*tau;
351  #endif
352 
353  /* -------------------------------------------------------------------
354  Find conservative eigenvectors; this is done by vector
355  multiplication, as shown on the reference paper ("A solution
356  adaptive upwind scheme for MHD", JCP 154, 284 (1999)), on page
357  297:
358 
359  Rc = dU/dW * Rp
360 
361  Primitive left eigenvectors are not necessary, since the jump
362  is invariant, i.e.:
363 
364  Lp*dW = Lc*dU
365  ------------------------------------------------------------------- */
366 
367  for (k = 0; k < NFLX; k++) {
368  Rc[RHO][k] = Rp[RHO][k];
369  EXPAND (Rc[MXn][k] = u*Rp[RHO][k] + rho*Rp[VXn][k]; ,
370  Rc[MXt][k] = v*Rp[RHO][k] + rho*Rp[VXt][k]; ,
371  Rc[MXb][k] = w*Rp[RHO][k] + rho*Rp[VXb][k];)
372  EXPAND (Rc[BXn][k] = Rp[BXn][k]; ,
373  Rc[BXt][k] = Rp[BXt][k]; ,
374  Rc[BXb][k] = Rp[BXb][k];)
375 
376  scrh0 = EXPAND(u*Rp[VXn][k], + v*Rp[VXt][k], + w*Rp[VXb][k]);
377  scrh1 = EXPAND(b1x*Rp[BXn][k], + b1y*Rp[BXt][k], + b1z*Rp[BXb][k]);
378  Rc[ENG][k] = 0.5*vel2*Rp[RHO][k] + rho*scrh0
379  + scrh1 + Rp[PRS][k]*gmm1_inv;
380  }
381 
382  /* ---- COMPUTE MAX EIGENVALUE ---- */
383 
384  cmax[i] = fabs (u) + cf;
385  g_maxMach = MAX (fabs (u / a), g_maxMach);
386 
387  SL[i] = lambda[KFASTM];
388  SR[i] = lambda[KFASTP];
389 
390  /* ---- ADD 'VISCOUS' FLUX FIRST ---- */
391 
392  for (nv = 0; nv < NFLX; nv++) {
393  scrh0 = 0.0;
394  for (k = 0; k < NFLX; k++) {
395  scrh1 = fabs(lambda[k]);
396  if ((k == KFASTM || k == KFASTP || k == KSLOWM || k == KSLOWP)
397  && scrh1 < 0.5*delta){ /* entropy fix */
398  scrh1 = scrh1*scrh1/delta + 0.25*delta;
399  }
400  scrh0 += scrh1*eta[k]*Rc[nv][k];
401  }
402  state->flux[i][nv] = 0.5*(fL[i][nv] + fR[i][nv] - scrh0);
403  }
404  state->press[i] = 0.5*(pL[i] + pR[i]);
405  }
406 
407  #if DIVB_CONTROL == EIGHT_WAVES
408  ROE_DivBSource (state, grid, beg, end);
409  #endif
410 }
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
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
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
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
STATE_1D::SL
double * SL
Leftmost velocity in the Riemann fan at i+1/2.
Definition: structs.h:166
sqrt_1_2
#define sqrt_1_2
Definition: roe.bckfld.c:3
j
int j
Definition: analysis.c:2
Sph_disk.vel2
list vel2
Definition: Sph_disk.py:39
VXt
int VXt
Definition: globals.h:73
KFASTP
Definition: mod_defs.h:58
k
int k
Definition: analysis.c:2
MXn
int MXn
Definition: globals.h:74
KFASTM
Definition: mod_defs.h:58
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
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
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
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: