PLUTO
Functions
tvdlf.c File Reference

Lax-Friedrechs (Rusanov) Riemann solver for MHD. More...

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

Go to the source code of this file.

Functions

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

Detailed Description

Lax-Friedrechs (Rusanov) Riemann solver for MHD.

Solve the Riemann problem for the adiabatic and isothermal MHD equations using the Lax-Friedrichs Rusanov Riemann solver with local maximum characteristic speed:

\[ \hat{F}_{i+\HALF} = \frac{F^L_{i+\HALF} + F^R_{i+\HALF}}{2} - c^{\rm max}_{i+\HALF}\frac{U^R_{i+\HALF} - U^L_{i-\HALF}}{2} \qquad{\rm where}\quad c^{\rm max} = \lambda_{\rm max}\left( \frac{V^R_{i+\HALF} + V^L_{i+\HALF}}{2}\right) \]

where $\lambda_{\rm \max}$ is a function of the arithmetic average between the left and the right states.

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:

  • "Comparison of Some Flux Correced Transport and Total Variation Diminishing Numerical Schemes for Hydrodynamics and Magnetohydrodynamics Problems", Toth and Odstrcil, JCP (1996), 128,82
Authors
A. Mignone (migno.nosp@m.ne@p.nosp@m.h.uni.nosp@m.to.i.nosp@m.t)
Date
April 7, 2014

Definition in file tvdlf.c.

Function Documentation

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

Solve Riemann problem for the adiabatic MHD equations using the Lax-Friedrichs (Rusanov) Riemann 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 38 of file tvdlf.c.

51 {
52  int nv, i;
53  double cRL;
54  double *uR, *uL;
55  static double **fL, **fR, **vRL;
56  static double *pR, *pL, *cmin_RL, *cmax_RL, *a2L, *a2R;
57  static double **VL, **VR, **UL, **UR;
58  double **bgf;
59 
60  if (fR == NULL){
61  fR = ARRAY_2D(NMAX_POINT, NFLX, double);
62  fL = ARRAY_2D(NMAX_POINT, NFLX, double);
63  vRL = ARRAY_2D(NMAX_POINT, NFLX, double);
64  cmin_RL = ARRAY_1D(NMAX_POINT, double);
65  cmax_RL = ARRAY_1D(NMAX_POINT, double);
66  pR = ARRAY_1D(NMAX_POINT, double);
67  pL = ARRAY_1D(NMAX_POINT, double);
68 
69  a2R = ARRAY_1D(NMAX_POINT, double);
70  a2L = ARRAY_1D(NMAX_POINT, double);
71  #ifdef GLM_MHD
72  VL = ARRAY_2D(NMAX_POINT, NVAR, double);
73  VR = ARRAY_2D(NMAX_POINT, NVAR, double);
74  UL = ARRAY_2D(NMAX_POINT, NVAR, double);
75  UR = ARRAY_2D(NMAX_POINT, NVAR, double);
76  #endif
77  }
78 
79  #if BACKGROUND_FIELD == YES
80  bgf = GetBackgroundField (beg, end, FACE_CENTER, grid);
81  #endif
82 
83  #ifdef GLM_MHD
84  GLM_Solve (state, VL, VR, beg, end, grid);
85  PrimToCons (VL, UL, beg, end);
86  PrimToCons (VR, UR, beg, end);
87  #else
88  VL = state->vL; UL = state->uL;
89  VR = state->vR; UR = state->uR;
90  #endif
91 
92 /* ----------------------------------------------------
93  compute sound speed & fluxes at zone interfaces
94  ---------------------------------------------------- */
95 
96  SoundSpeed2 (VL, a2L, NULL, beg, end, FACE_CENTER, grid);
97  SoundSpeed2 (VR, a2R, NULL, beg, end, FACE_CENTER, grid);
98 
99  Flux (UL, VL, a2L, bgf, fL, pL, beg, end);
100  Flux (UR, VR, a2R, bgf, fR, pR, beg, end);
101 
102 /* ------------------------------------------------------
103  Compute max eigenvalue and fluxes
104  ------------------------------------------------------ */
105 
106  for (i = beg; i <= end; i++) {
107  VAR_LOOP(nv) vRL[i][nv] = 0.5*(VL[i][nv] + VR[i][nv]);
108  #if EOS == IDEAL
109  g_maxMach = MAX(g_maxMach,
110  fabs(vRL[i][VXn])/sqrt(g_gamma*vRL[i][PRS]/vRL[i][RHO]));
111  #elif EOS == BAROTROPIC
112  cRL = sqrt(g_gamma*BAROTROPIC_PR(vRL[i][RHO])/vRL[i][RHO]);
113  g_maxMach = MAX(g_maxMach, fabs(vRL[i][VXn])/cRL);
114  #elif EOS == ISOTHERMAL
115  g_maxMach = MAX(g_maxMach, fabs(vRL[i][VXn])/g_isoSoundSpeed);
116  #endif
117  }
118 
119  SoundSpeed2 (vRL, a2R, NULL, beg, end, FACE_CENTER, grid);
120  MaxSignalSpeed (vRL, a2R, cmin_RL, cmax_RL, bgf, beg, end);
121 
122  for (i = beg; i <= end; i++) {
123  cRL = MAX(fabs(cmin_RL[i]), fabs(cmax_RL[i]));
124  state->SL[i] = -cRL;
125  state->SR[i] = cRL;
126  cmax[i] = cRL;
127  uL = UL[i];
128  uR = UR[i];
129  for (nv = NFLX; nv--; ) {
130  state->flux[i][nv] = 0.5*(fL[i][nv] + fR[i][nv] - cRL*(uR[nv] - uL[nv]));
131  }
132  state->press[i] = 0.5*(pL[i] + pR[i]);
133  }
134 
135 /* --------------------------------------------------------
136  initialize source term
137  -------------------------------------------------------- */
138 
139  #if DIVB_CONTROL == EIGHT_WAVES
140  Roe_DivBSource (state, beg + 1, end, grid);
141  #endif
142 }
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
STATE_1D::SR
double * SR
Rightmost velocity in the Riemann fan at i+1/2.
Definition: structs.h:167
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
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
VAR_LOOP
#define VAR_LOOP(n)
Definition: macros.h:226
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
Roe_DivBSource
void Roe_DivBSource(const State_1D *state, int is, int ie, Grid *grid)
Definition: source.c:39
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
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
NVAR
#define NVAR
Definition: pluto.h:609
STATE_1D::uL
double ** uL
same as vL, in conservative vars
Definition: structs.h:144

Here is the call graph for this function: