Skip to content
Snippets Groups Projects
Commit 1e44cacd authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Fixed rebase merge conflict in HLLC Riemann solver.

parent 6466e76c
Branches
Tags
1 merge request!588Gizmo mfm clean
Hide whitespace changes
Inline Side-by-side
src/riemann/riemann_hllc.h
+ 109
− 110
View file @ 1e44cacd
...@@ -43,19 +43,21 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux( ...@@ -43,19 +43,21 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
/* Handle pure vacuum */ /* Handle pure vacuum */
if (!WL[0] && !WR[0]) { if (!WL[0] && !WR[0]) {
totflux[0] = 0.f; totflux[0] = 0.0f;
totflux[1] = 0.f; totflux[1] = 0.0f;
totflux[2] = 0.f; totflux[2] = 0.0f;
totflux[3] = 0.f; totflux[3] = 0.0f;
totflux[4] = 0.f; totflux[4] = 0.0f;
return; return;
} }
/* STEP 0: obtain velocity in interface frame */ /* STEP 0: obtain velocity in interface frame */
const float uL = WL[1] * n[0] + WL[2] * n[1] + WL[3] * n[2]; const float uL = WL[1] * n[0] + WL[2] * n[1] + WL[3] * n[2];
const float uR = WR[1] * n[0] + WR[2] * n[1] + WR[3] * n[2]; const float uR = WR[1] * n[0] + WR[2] * n[1] + WR[3] * n[2];
const float aL = sqrtf(hydro_gamma * WL[4] / WL[0]); const float rhoLinv = 1.0f / WL[0];
const float aR = sqrtf(hydro_gamma * WR[4] / WR[0]); const float rhoRinv = 1.0f / WR[0];
const float aL = sqrtf(hydro_gamma * WL[4] * rhoLinv);
const float aR = sqrtf(hydro_gamma * WR[4] * rhoRinv);
/* Handle vacuum: vacuum does not require iteration and is always exact */ /* Handle vacuum: vacuum does not require iteration and is always exact */
if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) { if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) {
...@@ -64,98 +66,92 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux( ...@@ -64,98 +66,92 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
} }
/* STEP 1: pressure estimate */ /* STEP 1: pressure estimate */
const float rhobar = 0.5f * (WL[0] + WR[0]); const float rhobar = WL[0] + WR[0];
const float abar = 0.5f * (aL + aR); const float abar = aL + aR;
const float pPVRS = 0.5f * (WL[4] + WR[4]) - 0.5f * (uR - uL) * rhobar * abar; const float pPVRS =
const float pstar = max(0.f, pPVRS); 0.5f * ((WL[4] + WR[4]) - 0.25f * (uR - uL) * rhobar * abar);
const float pstar = max(0.0f, pPVRS);
/* STEP 2: wave speed estimates /* STEP 2: wave speed estimates
all these speeds are along the interface normal, since uL and uR are */ all these speeds are along the interface normal, since uL and uR are */
float qL = 1.f; float qL = 1.0f;
if (pstar > WL[4] && WL[4] > 0.f) { if (pstar > WL[4] && WL[4] > 0.0f) {
qL = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma * qL = sqrtf(1.0f +
(pstar / WL[4] - 1.f)); 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
(pstar / WL[4] - 1.0f));
} }
float qR = 1.f; float qR = 1.0f;
if (pstar > WR[4] && WR[4] > 0.f) { if (pstar > WR[4] && WR[4] > 0.0f) {
qR = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma * qR = sqrtf(1.0f +
(pstar / WR[4] - 1.f)); 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
(pstar / WR[4] - 1.0f));
} }
const float SL = uL - aL * qL; const float SLmuL = -aL * qL;
const float SR = uR + aR * qR; const float SRmuR = aR * qR;
const float Sstar = const float Sstar =
(WR[4] - WL[4] + WL[0] * uL * (SL - uL) - WR[0] * uR * (SR - uR)) / (WR[4] - WL[4] + WL[0] * uL * SLmuL - WR[0] * uR * SRmuR) /
(WL[0] * (SL - uL) - WR[0] * (SR - uR)); (WL[0] * SLmuL - WR[0] * SRmuR);
/* STEP 3: HLLC flux in a frame moving with the interface velocity */ /* STEP 3: HLLC flux in a frame moving with the interface velocity */
if (Sstar >= 0.f) { if (Sstar >= 0.0f) {
const float rhoLuL = WL[0] * uL;
const float v2 = WL[1] * WL[1] + WL[2] * WL[2] + WL[3] * WL[3];
const float eL =
WL[4] * rhoLinv * hydro_one_over_gamma_minus_one + 0.5f * v2;
const float SL = SLmuL + uL;
/* flux FL */ /* flux FL */
totflux[0] = WL[0] * uL; totflux[0] = rhoLuL;
/* these are the actual correct fluxes in the boosted lab frame /* these are the actual correct fluxes in the boosted lab frame
(not rotated to interface frame) */ (not rotated to interface frame) */
totflux[1] = WL[0] * WL[1] * uL + WL[4] * n[0]; totflux[1] = rhoLuL * WL[1] + WL[4] * n[0];
totflux[2] = WL[0] * WL[2] * uL + WL[4] * n[1]; totflux[2] = rhoLuL * WL[2] + WL[4] * n[1];
totflux[3] = WL[0] * WL[3] * uL + WL[4] * n[2]; totflux[3] = rhoLuL * WL[3] + WL[4] * n[2];
const float v2 = WL[1] * WL[1] + WL[2] * WL[2] + WL[3] * WL[3]; totflux[4] = rhoLuL * eL + WL[4] * uL;
const float eL = WL[4] * hydro_one_over_gamma_minus_one / WL[0] + 0.5f * v2;
totflux[4] = WL[0] * eL * uL + WL[4] * uL; if (SL < 0.0f) {
if (SL < 0.f) {
const float starfac = SLmuL / (SL - Sstar) - 1.0f;
float UstarL[5]; const float rhoLSL = WL[0] * SL;
const float SstarmuL = Sstar - uL;
/* add flux FstarL */ const float rhoLSLstarfac = rhoLSL * starfac;
UstarL[0] = 1.f; const float rhoLSLSstarmuL = rhoLSL * SstarmuL;
/* we need UstarL in the lab frame:
* subtract the velocity in the interface frame from the lab frame totflux[0] += rhoLSLstarfac;
* velocity and then add Sstar in interface frame */ totflux[1] += rhoLSLstarfac * WL[1] + rhoLSLSstarmuL * n[0];
UstarL[1] = WL[1] + (Sstar - uL) * n[0]; totflux[2] += rhoLSLstarfac * WL[2] + rhoLSLSstarmuL * n[1];
UstarL[2] = WL[2] + (Sstar - uL) * n[1]; totflux[3] += rhoLSLstarfac * WL[3] + rhoLSLSstarmuL * n[2];
UstarL[3] = WL[3] + (Sstar - uL) * n[2]; totflux[4] += rhoLSLstarfac * eL +
UstarL[4] = eL + (Sstar - uL) * (Sstar + WL[4] / (WL[0] * (SL - uL))); rhoLSLSstarmuL * (Sstar + WL[4] / (WL[0] * SLmuL));
UstarL[0] *= WL[0] * (SL - uL) / (SL - Sstar);
UstarL[1] *= WL[0] * (SL - uL) / (SL - Sstar);
UstarL[2] *= WL[0] * (SL - uL) / (SL - Sstar);
UstarL[3] *= WL[0] * (SL - uL) / (SL - Sstar);
UstarL[4] *= WL[0] * (SL - uL) / (SL - Sstar);
totflux[0] += SL * (UstarL[0] - WL[0]);
totflux[1] += SL * (UstarL[1] - WL[0] * WL[1]);
totflux[2] += SL * (UstarL[2] - WL[0] * WL[2]);
totflux[3] += SL * (UstarL[3] - WL[0] * WL[3]);
totflux[4] += SL * (UstarL[4] - WL[0] * eL);
} }
} else { } else {
/* flux FR */ const float rhoRuR = WR[0] * uR;
totflux[0] = WR[0] * uR;
totflux[1] = WR[0] * WR[1] * uR + WR[4] * n[0];
totflux[2] = WR[0] * WR[2] * uR + WR[4] * n[1];
totflux[3] = WR[0] * WR[3] * uR + WR[4] * n[2];
const float v2 = WR[1] * WR[1] + WR[2] * WR[2] + WR[3] * WR[3]; const float v2 = WR[1] * WR[1] + WR[2] * WR[2] + WR[3] * WR[3];
const float eR = WR[4] * hydro_one_over_gamma_minus_one / WR[0] + 0.5f * v2; const float eR =
totflux[4] = WR[0] * eR * uR + WR[4] * uR; WR[4] * rhoRinv * hydro_one_over_gamma_minus_one + 0.5f * v2;
if (SR > 0.f) { const float SR = SRmuR + uR;
float UstarR[5]; /* flux FR */
totflux[0] = rhoRuR;
/* add flux FstarR */ totflux[1] = rhoRuR * WR[1] + WR[4] * n[0];
UstarR[0] = 1.f; totflux[2] = rhoRuR * WR[2] + WR[4] * n[1];
totflux[3] = rhoRuR * WR[3] + WR[4] * n[2];
/* we need UstarR in the lab frame: totflux[4] = rhoRuR * eR + WR[4] * uR;
* subtract the velocity in the interface frame from the lab frame
* velocity and then add Sstar in interface frame */ if (SR > 0.0f) {
UstarR[1] = WR[1] + (Sstar - uR) * n[0];
UstarR[2] = WR[2] + (Sstar - uR) * n[1]; const float starfac = SRmuR / (SR - Sstar) - 1.0f;
UstarR[3] = WR[3] + (Sstar - uR) * n[2]; const float rhoRSR = WR[0] * SR;
UstarR[4] = eR + (Sstar - uR) * (Sstar + WR[4] / (WR[0] * (SR - uR))); const float SstarmuR = Sstar - uR;
UstarR[0] *= WR[0] * (SR - uR) / (SR - Sstar); const float rhoRSRstarfac = rhoRSR * starfac;
UstarR[1] *= WR[0] * (SR - uR) / (SR - Sstar); const float rhoRSRSstarmuR = rhoRSR * SstarmuR;
UstarR[2] *= WR[0] * (SR - uR) / (SR - Sstar);
UstarR[3] *= WR[0] * (SR - uR) / (SR - Sstar); totflux[0] += rhoRSRstarfac;
UstarR[4] *= WR[0] * (SR - uR) / (SR - Sstar); totflux[1] += rhoRSRstarfac * WR[1] + rhoRSRSstarmuR * n[0];
totflux[0] += SR * (UstarR[0] - WR[0]); totflux[2] += rhoRSRstarfac * WR[2] + rhoRSRSstarmuR * n[1];
totflux[1] += SR * (UstarR[1] - WR[0] * WR[1]); totflux[3] += rhoRSRstarfac * WR[3] + rhoRSRSstarmuR * n[2];
totflux[2] += SR * (UstarR[2] - WR[0] * WR[2]); totflux[4] += rhoRSRstarfac * eR +
totflux[3] += SR * (UstarR[3] - WR[0] * WR[3]); rhoRSRSstarmuR * (Sstar + WR[4] / (WR[0] * SRmuR));
totflux[4] += SR * (UstarR[4] - WR[0] * eR);
} }
} }
...@@ -189,11 +185,11 @@ riemann_solve_for_middle_state_flux(const float *WL, const float *WR, ...@@ -189,11 +185,11 @@ riemann_solve_for_middle_state_flux(const float *WL, const float *WR,
/* Handle pure vacuum */ /* Handle pure vacuum */
if (!WL[0] && !WR[0]) { if (!WL[0] && !WR[0]) {
totflux[0] = 0.f; totflux[0] = 0.0f;
totflux[1] = 0.f; totflux[1] = 0.0f;
totflux[2] = 0.f; totflux[2] = 0.0f;
totflux[3] = 0.f; totflux[3] = 0.0f;
totflux[4] = 0.f; totflux[4] = 0.0f;
return; return;
} }
...@@ -205,37 +201,40 @@ riemann_solve_for_middle_state_flux(const float *WL, const float *WR, ...@@ -205,37 +201,40 @@ riemann_solve_for_middle_state_flux(const float *WL, const float *WR,
/* Handle vacuum: vacuum does not require iteration and is always exact */ /* Handle vacuum: vacuum does not require iteration and is always exact */
if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) { if (riemann_is_vacuum(WL, WR, uL, uR, aL, aR)) {
totflux[0] = 0.f; totflux[0] = 0.0f;
totflux[1] = 0.f; totflux[1] = 0.0f;
totflux[2] = 0.f; totflux[2] = 0.0f;
totflux[3] = 0.f; totflux[3] = 0.0f;
totflux[4] = 0.f; totflux[4] = 0.0f;
return; return;
} }
/* STEP 1: pressure estimate */ /* STEP 1: pressure estimate */
const float rhobar = 0.5f * (WL[0] + WR[0]); const float rhobar = WL[0] + WR[0];
const float abar = 0.5f * (aL + aR); const float abar = aL + aR;
const float pPVRS = 0.5f * (WL[4] + WR[4]) - 0.5f * (uR - uL) * rhobar * abar; const float pPVRS =
0.5f * ((WL[4] + WR[4]) - 0.25f * (uR - uL) * rhobar * abar);
const float pstar = max(0.f, pPVRS); const float pstar = max(0.f, pPVRS);
/* STEP 2: wave speed estimates /* STEP 2: wave speed estimates
all these speeds are along the interface normal, since uL and uR are */ all these speeds are along the interface normal, since uL and uR are */
float qL = 1.f; float qL = 1.0f;
if (pstar > WL[4] && WL[4] > 0.f) { if (pstar > WL[4] && WL[4] > 0.0f) {
qL = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma * qL = sqrtf(1.0f +
(pstar / WL[4] - 1.f)); 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
(pstar / WL[4] - 1.0f));
} }
float qR = 1.f; float qR = 1.0f;
if (pstar > WR[4] && WR[4] > 0.f) { if (pstar > WR[4] && WR[4] > 0.0f) {
qR = sqrtf(1.f + 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma * qR = sqrtf(1.0f +
(pstar / WR[4] - 1.f)); 0.5f * hydro_gamma_plus_one * hydro_one_over_gamma *
(pstar / WR[4] - 1.0f));
} }
const float SL = uL - aL * qL; const float SLmuL = -aL * qL;
const float SR = uR + aR * qR; const float SRmuR = aR * qR;
const float Sstar = const float Sstar =
(WR[4] - WL[4] + WL[0] * uL * (SL - uL) - WR[0] * uR * (SR - uR)) / (WR[4] - WL[4] + WL[0] * uL * SLmuL - WR[0] * uR * SRmuR) /
(WL[0] * (SL - uL) - WR[0] * (SR - uR)); (WL[0] * SLmuL - WR[0] * SRmuR);
totflux[0] = 0.0f; totflux[0] = 0.0f;
totflux[1] = pstar * n[0]; totflux[1] = pstar * n[0];
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment