PLUTO
Functions
rusanov-dw.c File Reference
#include "pluto.h"
Include dependency graph for rusanov-dw.c:

Go to the source code of this file.

Functions

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

Function Documentation

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

Solve Riemann problem using the Lax-Friedrichs Rusanov solver:

F(i+1/2) = [FL + FR - c*(UR - UL)]*0.5

where c = max_speed[(UR + UL)*0.5]

LAST_MODIFIED

April 4th 2006, by Andrea Mignone (migno.nosp@m.ne@t.nosp@m.o.ast.nosp@m.ro.i.nosp@m.t)

Definition at line 4 of file rusanov-dw.c.

18 {
19  int nv, i;
20  static double **fL, **fR, **vRL;
21  static double *cRL_min, *cRL_max, *pL, *pR, *a2L, *a2R;
22  double *uR, *uL, *flux;
23 
24 double num, den, lambda, dU, dF, vn;
25 
26  if (fR == NULL){
27  fR = ARRAY_2D(NMAX_POINT, NFLX, double);
28  fL = ARRAY_2D(NMAX_POINT, NFLX, double);
29  vRL = ARRAY_2D(NMAX_POINT, NFLX, double);
30  cRL_min = ARRAY_1D(NMAX_POINT, double);
31  cRL_max = ARRAY_1D(NMAX_POINT, double);
32  pR = ARRAY_1D(NMAX_POINT, double);
33  pL = ARRAY_1D(NMAX_POINT, double);
34  a2R = ARRAY_1D(NMAX_POINT, double);
35  a2L = ARRAY_1D(NMAX_POINT, double);
36  }
37 
38 /* ----------------------------------------------------
39  compute sound speed & fluxes at zone interfaces
40  ---------------------------------------------------- */
41 
42  SoundSpeed2 (state->vL, a2L, NULL, beg, end, FACE_CENTER, grid);
43  SoundSpeed2 (state->vR, a2R, NULL, beg, end, FACE_CENTER, grid);
44 
45  Flux (state->uL, state->vL, a2L, fL, pL, beg, end);
46  Flux (state->uR, state->vR, a2R, fR, pR, beg, end);
47 
48 /* ---------------------------------------------------------------------
49  Compute average state in order to get the local max characteristic
50  velocity.
51  In order to avoid underestimating this speed at strong reflecting
52  boundaries (or when vxL = -vxR in general) , we average the
53  absolute value of the normal velocity.
54  ------------------------------------------------------------------- */
55 
56  for (i = beg; i <= end; i++) {
57  for (nv = NFLX; nv--; ) {
58  vRL[i][nv] = 0.5*(state->vL[i][nv] + state->vR[i][nv]);
59  }
60  vRL[i][VXn] = 0.5*(fabs(state->vL[i][VXn]) + fabs(state->vR[i][VXn]));
61  }
62 
63 /* ---------------------------------------------------------------------
64  Compute sound speed, max and min signal speeds
65  --------------------------------------------------------------------- */
66 
67  SoundSpeed2 (vRL, a2R, NULL, beg, end, FACE_CENTER, grid);
68  MaxSignalSpeed (vRL, a2R, cRL_min, cRL_max, beg, end);
69 
70  for (i = beg; i <= end; i++) {
71  uR = state->uR[i];
72  uL = state->uL[i];
73  flux = state->flux[i];
74 
75  /* -- compute max eigenvalue -- */
76 
77  cmax[i] = MAX(fabs(cRL_max[i]), fabs(cRL_min[i]));
78  g_maxMach = MAX(g_maxMach, fabs(vRL[i][VXn])/sqrt(a2R[i]));
79 
80  /* -- compute fluxes -- */
81 num = den = 0.0;
82 for (nv = 0; nv < NFLX; nv++) {
83  dU = uR[nv] - uL[nv];
84  dF = fR[i][nv] - fL[i][nv] + (nv == MXn ? (pR[i] - pL[i]):0.0);
85 
86  num += dU*dF;
87  den += dU*dU;
88 }
89 lambda = fabs(num)/(fabs(den) + 1.e-12);
90 
91 if (lambda > cmax[i]) lambda = cmax[i];
92 
93 vn = 0.5*( fabs(state->vR[i][VXn]) + fabs(state->vL[i][VXn]) );
94 vn = MAX( fabs(state->vR[i][VXn]), fabs(state->vL[i][VXn]) );
95 if (lambda < vn) lambda = vn;
96 
97 /*
98 cL = MIN(fabs(cmax_L[i]), fabs(cmin_L[i]));
99 cR = MIN(fabs(cmax_R[i]), fabs(cmin_R[i]));
100 */
101 
102  for (nv = NFLX; nv--; ) {
103  flux[nv] = 0.5*(fL[i][nv] + fR[i][nv] - lambda*(uR[nv] - uL[nv]));
104  }
105  state->press[i] = 0.5*(pL[i] + pR[i]);
106  }
107 }
MAX
#define MAX(a, b)
Definition: macros.h:101
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
STATE_1D::vR
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
SoundSpeed2
void SoundSpeed2(double **v, double *cs2, double *h, int beg, int end, int pos, Grid *grid)
Definition: eos.c:16
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
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
i
int i
Definition: analysis.c:2
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
MaxSignalSpeed
void MaxSignalSpeed(double **v, double *cs2, double *cmin, double *cmax, int beg, int end)
Definition: eigenv.c:34
ARRAY_2D
#define ARRAY_2D(nx, ny, type)
Definition: prototypes.h:171
STATE_1D::uL
double ** uL
same as vL, in conservative vars
Definition: structs.h:144

Here is the call graph for this function: