PLUTO
Functions
hllc.c File Reference

Implement the HLLC Riemann solver for relativistic HD. 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, double *cmax, Grid *grid)
 

Detailed Description

Implement the HLLC Riemann solver for relativistic HD.

Solve the Riemann problem for the relativistic hydro (RHD) equations using the HLLC solver of Mignone & Bodo (2005).

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:

  • "An HLLC Riemann solver for relativistic flows - I. Hydrodynamics", Mignone and Bodo, MNRAS (2005) 364,1126.
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 hllc.c.

Function Documentation

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

Solve the RHD 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;
43 
44  double usl[NFLX], usr[NFLX], vm[NFLX], us, ps;
45  double AL, BL, AR, BR, a,b,c;
46  double vxl, vxr;
47 
48  double *vL, *vR, *uL, *uR;
49  static double *SL, *SR;
50  static double **Uhll, **Fhll;
51  static double **fL, **fR;
52  static double *pR, *pL;
53  static double *a2L, *a2R, *hL, *hR;
54 
55  if (fL == NULL){
56  fL = ARRAY_2D(NMAX_POINT, NFLX, double);
57  fR = ARRAY_2D(NMAX_POINT, NFLX, double);
58 
59  Uhll = ARRAY_2D(NMAX_POINT, NFLX, double);
60  Fhll = ARRAY_2D(NMAX_POINT, NFLX, double);
61 
62  pR = ARRAY_1D(NMAX_POINT, double);
63  pL = ARRAY_1D(NMAX_POINT, double);
64  SR = ARRAY_1D(NMAX_POINT, double);
65  SL = ARRAY_1D(NMAX_POINT, double);
66 
67  a2L = ARRAY_1D(NMAX_POINT, double);
68  a2R = ARRAY_1D(NMAX_POINT, double);
69 
70  hL = ARRAY_1D(NMAX_POINT, double);
71  hR = ARRAY_1D(NMAX_POINT, double);
72  }
73 
74 /* ----------------------------------------------------
75  compute sound speed & fluxes at zone interfaces
76  ---------------------------------------------------- */
77 
78  SoundSpeed2 (state->vL, a2L, hL, beg, end, FACE_CENTER, grid);
79  SoundSpeed2 (state->vR, a2R, hR, beg, end, FACE_CENTER, grid);
80 
81  Flux (state->uL, state->vL, a2L, fL, pL, beg, end);
82  Flux (state->uR, state->vR, a2R, fR, pR, beg, end);
83 
84  HLL_Speed (state->vL, state->vR, a2L, a2R, SL, SR, beg, end);
85  for (i = beg; i <= end; i++) {
86 
87  scrh = MAX(fabs(SL[i]), fabs(SR[i]));
88  cmax[i] = scrh;
89 
90 /* --------------------------------------------------
91  compute HLL state and flux
92  -------------------------------------------------- */
93 /*
94  scrh = 1.0/(Sr - Sl);
95  for (nv = NFLX; nv--; ){
96  Uhll[i][nv] = Sr*ur[i][nv] - Sl*ul[i][nv]
97  + fl[i][nv] - fr[i][nv];
98  Uhll[i][nv] *= scrh;
99 
100  Fhll[i][nv] = Sl*Sr*(ur[i][nv] - ul[i][nv])
101  + Sr*fl[i][nv] - Sl*fr[i][nv];
102  Fhll[i][nv] *= scrh;
103  }
104  Uhll[i][MXn] += (pl[i] - pr[i])*scrh;
105  Fhll[i][MXn] += (Sr*pl[i] - Sl*pr[i])*scrh;
106 */
107 /* --------------------------------------------------
108  compute HLLC flux
109  -------------------------------------------------- */
110 
111  if (SL[i] >= 0.0){
112 
113  for (nv = NFLX; nv--; ) state->flux[i][nv] = fL[i][nv];
114  state->press[i] = pL[i];
115 
116  }else if (SR[i] <= 0.0){
117 
118  for (nv = NFLX; nv--; ) state->flux[i][nv] = fR[i][nv];
119  state->press[i] = pR[i];
120 
121  }else{
122 
123  vL = state->vL[i]; uL = state->uL[i];
124  vR = state->vR[i]; uR = state->uR[i];
125 
126 #if SHOCK_FLATTENING == MULTID
127  if ((state->flag[i] & FLAG_HLL) || (state->flag[i+1] & FLAG_HLL)){
128  scrh = 1.0/(SR[i] - SL[i]);
129  for (nv = NFLX; nv--; ){
130  state->flux[i][nv] = SL[i]*SR[i]*(uR[nv] - uL[nv])
131  + SR[i]*fL[i][nv] - SL[i]*fR[i][nv];
132  state->flux[i][nv] *= scrh;
133  }
134  state->press[i] = (SR[i]*pL[i] - SL[i]*pR[i])*scrh;
135  continue;
136  }
137 #endif
138 
139  vxl = vL[VXn];
140  vxr = vR[VXn];
141 
142 /* ---------------------------------------
143  get u*
144  --------------------------------------- */
145 
146  AL = SL[i]*uL[ENG] - fL[i][ENG];
147  AR = SR[i]*uR[ENG] - fR[i][ENG];
148 
149  BL = SL[i]*uL[MXn] - fL[i][MXn] - pL[i];
150  BR = SR[i]*uR[MXn] - fR[i][MXn] - pR[i];
151 
152  a = AR*SL[i] - AL*SR[i];
153  b = AL + BL*SR[i] - AR - BR*SL[i];
154  c = BR - BL;
155 /*
156  if (fabs(a) > 1.e-9){
157  usp = 0.5*(- b + sqrt(b*b - 4.0*a*c))/a;
158  usm = 0.5*(- b - sqrt(b*b - 4.0*a*c))/a;
159  }else{
160  usp = usm = -c/b;
161  }
162 */
163  scrh = -0.5*(b + DSIGN(b)*sqrt(b*b - 4.0*a*c));
164  us = c/scrh;
165 
166  ps = (AL*us - BL)/(1.0 - us*SL[i]);
167 
168  usl[RHO] = uL[RHO]*(SL[i] - vxl)/(SL[i] - us);
169  usr[RHO] = uR[RHO]*(SR[i] - vxr)/(SR[i] - us);
170  EXPAND(usl[MXn] = (SL[i]*(uL[ENG] + ps) - uL[MXn])*us/(SL[i] - us);
171  usr[MXn] = (SR[i]*(uR[ENG] + ps) - uR[MXn])*us/(SR[i] - us); ,
172  usl[MXt] = uL[MXt]*(SL[i] - vxl)/(SL[i] - us);
173  usr[MXt] = uR[MXt]*(SR[i] - vxr)/(SR[i] - us); ,
174  usl[MXb] = uL[MXb]*(SL[i] - vxl)/(SL[i] - us);
175  usr[MXb] = uR[MXb]*(SR[i] - vxr)/(SR[i] - us);)
176 
177  usl[ENG] = uL[ENG] + (usl[MXn] - uL[MXn])/SL[i];
178  usr[ENG] = uR[ENG] + (usr[MXn] - uR[MXn])/SR[i];
179 
180  /* ---- Compute HLLC flux ---- */
181 
182  if (us >= 0.0) {
183  for (nv = NFLX; nv--; ) {
184  state->flux[i][nv] = fL[i][nv] + SL[i]*(usl[nv] - uL[nv]);
185  }
186  state->press[i] = pL[i];
187  }else {
188  for (nv = NFLX; nv--; ) {
189  state->flux[i][nv] = fR[i][nv] + SR[i]*(usr[nv] - uR[nv]);
190  }
191  state->press[i] = pR[i];
192  }
193  } /* -- end loop on speed signs -- */
194  } /* -- end loop on points -- */
195 }
a
static double a
Definition: init.c:135
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
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
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
STATE_1D::flag
unsigned char * flag
Definition: structs.h:168
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
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: