PLUTO
Functions
hllc.c File Reference

HLLC Riemann solver for MHD. More...

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

Go to the source code of this file.

Functions

void HLLC_Solver (const State_1D *state, int beg, int end, real *cmax, Grid *grid)
 

Detailed Description

HLLC Riemann solver for MHD.

Solve the Riemann problem for the HD equations using the two-state HLLC solver by Toro.

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:

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

Definition in file hllc.c.

Function Documentation

void HLLC_Solver ( const State_1D *  state,
int  beg,
int  end,
real *  cmax,
Grid *  grid 
)

Solve Riemann problem using the HLLC 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 28 of file hllc.c.

40 {
41  int nv, i;
42  double scrh, vxr, vxl;
43  double usL[NFLX], usR[NFLX], vs;
44  double qL, qR, wL, wR;
45  double *vL, *vR, *uL, *uR;
46  #if EOS == ISOTHERMAL
47  double rho, mx;
48  #endif
49  static double *pL, *pR, *SL, *SR, *a2L, *a2R;
50  static double **fL, **fR;
51 
52 /* -- Allocate memory -- */
53 
54  if (fL == NULL){
55  fL = ARRAY_2D(NMAX_POINT, NFLX, double);
56  fR = ARRAY_2D(NMAX_POINT, NFLX, double);
57 
58  pR = ARRAY_1D(NMAX_POINT, double);
59  pL = ARRAY_1D(NMAX_POINT, double);
60  SR = ARRAY_1D(NMAX_POINT, double);
61  SL = ARRAY_1D(NMAX_POINT, double);
62 
63  a2R = ARRAY_1D(NMAX_POINT, double);
64  a2L = ARRAY_1D(NMAX_POINT, double);
65  }
66 
67 /* ----------------------------------------------------
68  Compute sound speed & fluxes at zone interfaces
69  ---------------------------------------------------- */
70 
71  SoundSpeed2 (state->vL, a2L, NULL, beg, end, FACE_CENTER, grid);
72  SoundSpeed2 (state->vR, a2R, NULL, beg, end, FACE_CENTER, grid);
73 
74  Flux (state->uL, state->vL, a2L, fL, pL, beg, end);
75  Flux (state->uR, state->vR, a2R, fR, pR, beg, end);
76 
77  HLL_Speed (state->vL, state->vR, a2L, a2R, SL, SR, beg, end);
78 
79  for (i = beg; i <= end; i++) {
80 
81  scrh = MAX(fabs(SL[i]), fabs(SR[i]));
82  cmax[i] = scrh;
83 
84  if (SL[i] > 0.0){
85 
86  for (nv = NFLX; nv--; ) state->flux[i][nv] = fL[i][nv];
87  state->press[i] = pL[i];
88 
89  }else if (SR[i] < 0.0){
90 
91  for (nv = NFLX; nv--; ) state->flux[i][nv] = fR[i][nv];
92  state->press[i] = pR[i];
93 
94  }else{
95 
96  vR = state->vR[i]; uR = state->uR[i];
97  vL = state->vL[i]; uL = state->uL[i];
98 
99  vxr = vR[VXn];
100  vxl = vL[VXn];
101 
102 #if SHOCK_FLATTENING == MULTID
103  if ((state->flag[i] & FLAG_HLL) || (state->flag[i+1] & FLAG_HLL)){
104  scrh = 1.0/(SR[i] - SL[i]);
105  for (nv = NFLX; nv--; ){
106  state->flux[i][nv] = SL[i]*SR[i]*(uR[nv] - uL[nv])
107  + SR[i]*fL[i][nv] - SL[i]*fR[i][nv];
108  state->flux[i][nv] *= scrh;
109  }
110  state->press[i] = (SR[i]*pL[i] - SL[i]*pR[i])*scrh;
111  continue;
112  }
113 #endif
114 
115  /* ---------------------------------------
116  get u*
117  --------------------------------------- */
118 
119  #if HAVE_ENERGY
120  qL = vL[PRS] + uL[MXn]*(vL[VXn] - SL[i]);
121  qR = vR[PRS] + uR[MXn]*(vR[VXn] - SR[i]);
122 
123  wL = vL[RHO]*(vL[VXn] - SL[i]);
124  wR = vR[RHO]*(vR[VXn] - SR[i]);
125 
126  vs = (qR - qL)/(wR - wL); /* wR - wL > 0 since SL < 0, SR > 0 */
127 /*
128  vs = vR[PRS] - vL[PRS] + uL[MXn]*(SL[i] - vxl)
129  - uR[MXn]*(SR[i] - vxr);
130  vs /= vL[RHO]*(SL[i] - vxl) - vR[RHO]*(SR[i] - vxr);
131 */
132 
133  usL[RHO] = uL[RHO]*(SL[i] - vxl)/(SL[i] - vs);
134  usR[RHO] = uR[RHO]*(SR[i] - vxr)/(SR[i] - vs);
135  EXPAND(usL[MXn] = usL[RHO]*vs; usR[MXn] = usR[RHO]*vs; ,
136  usL[MXt] = usL[RHO]*vL[VXt]; usR[MXt] = usR[RHO]*vR[VXt]; ,
137  usL[MXb] = usL[RHO]*vL[VXb]; usR[MXb] = usR[RHO]*vR[VXb];)
138 
139  usL[ENG] = uL[ENG]/vL[RHO]
140  + (vs - vxl)*(vs + vL[PRS]/(vL[RHO]*(SL[i] - vxl)));
141  usR[ENG] = uR[ENG]/vR[RHO]
142  + (vs - vxr)*(vs + vR[PRS]/(vR[RHO]*(SR[i] - vxr)));
143 
144  usL[ENG] *= usL[RHO];
145  usR[ENG] *= usR[RHO];
146  #elif EOS == ISOTHERMAL
147  scrh = 1.0/(SR[i] - SL[i]);
148  rho = (SR[i]*uR[RHO] - SL[i]*uL[RHO] - fR[i][RHO] + fL[i][RHO])*scrh;
149  mx = (SR[i]*uR[MXn] - SL[i]*uL[MXn] - fR[i][MXn] + fL[i][MXn])*scrh;
150 
151  usL[RHO] = usR[RHO] = rho;
152  usL[MXn] = usR[MXn] = mx;
153  vs = ( SR[i]*fL[i][RHO] - SL[i]*fR[i][RHO]
154  + SR[i]*SL[i]*(uR[RHO] - uL[RHO]));
155  vs *= scrh;
156  vs /= rho;
157  EXPAND( ,
158  usL[MXt] = rho*vL[VXt]; usR[MXt] = rho*vR[VXt]; ,
159  usL[MXb] = rho*vL[VXb]; usR[MXb] = rho*vR[VXb];)
160  #endif
161 
162 
163  /* ---- verify consistency condition ------ */
164  /*
165  for (nv = 0; nv < NFLX; nv++) {
166  scrh = (vs - SL[i])*usL[nv] + (SR[i] - vs)*usR[nv];
167  scrh -= SR[i]*uR[i][nv] - SL[i]*uL[i][nv] +
168  fl[i][nv] - fr[i][nv];
169  if (fabs(scrh) > 1.e-9){
170  printf (" Consistency condition violated\n");
171  printf ("%d %d %12.6e\n",i,nv, scrh);
172  }
173  }
174  */
175 /* ---- Compute HLLC flux ---- */
176 
177  if (vs >= 0.0){
178  for (nv = NFLX; nv--; ) {
179  state->flux[i][nv] = fL[i][nv] + SL[i]*(usL[nv] - uL[nv]);
180  }
181  state->press[i] = pL[i];
182  } else {
183  for (nv = NFLX; nv--; ) {
184  state->flux[i][nv] = fR[i][nv] + SR[i]*(usR[nv] - uR[nv]);
185  }
186  state->press[i] = pR[i];
187  }
188  }
189  } /* end loops on points */
190 }
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
RHO
#define RHO
Definition: mod_defs.h:19
configure.scrh
tuple scrh
Definition: configure.py:200
STATE_1D::vR
double ** vR
Primitive variables to the right of the interface, .
Definition: structs.h:139
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
MXt
int MXt
Definition: globals.h:74
NFLX
#define NFLX
Definition: mod_defs.h:32
FACE_CENTER
#define FACE_CENTER
Definition: pluto.h:206
VXt
int VXt
Definition: globals.h:73
STATE_1D::flag
unsigned char * flag
Definition: structs.h:168
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
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
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: