Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Use Roe-averaged state in HLLE to compute wave speeds. #93

Merged
merged 3 commits into from
Jan 30, 2024
Merged

Use Roe-averaged state in HLLE to compute wave speeds. #93

Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
fd9538c
HLLE to use wave speeds from L/R/Roe
Nov 30, 2023
8f41d31
cpp-py-formatter
Jan 27, 2024
7574933
Merge branch 'main' into forrestglines/hlle-roe-wavespeeds
forrestglines Jan 30, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 32 additions & 29 deletions src/hydro/rsolvers/hydro_hlle.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ struct Riemann<Fluid::euler, RiemannSolver::hlle> {
Real gm1 = gamma - 1.0;
Real igm1 = 1.0 / gm1;
parthenon::par_for_inner(member, il, iu, [&](const int i) {
Real wli[(NHYDRO)], wri[(NHYDRO)];
Real wli[(NHYDRO)], wri[(NHYDRO)], wroe[(NHYDRO)];
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We're only using three out of five variable in wroe[(NHYDRO)] so from a performance point of view (kernel size) might be worth to just use three plain variables.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll make a note of this for the flux optimizations I'm doing. Register pressure might be a big bottle neck on the flux kernel performance. This fix would be useful, although I suspect the compile is already optimizing it out.

Real fl[(NHYDRO)], fr[(NHYDRO)], flxi[(NHYDRO)];
//--- Step 1. Load L/R states into local variables
wli[IDN] = wl(IDN, i);
Expand All @@ -65,34 +65,37 @@ struct Riemann<Fluid::euler, RiemannSolver::hlle> {
wri[IV3] = wr(ivz, i);
wri[IPR] = wr(IPR, i);

//--- Step 2. Compute middle state estimates with PVRS (Toro 10.5.2)
Real al, ar, el, er;
Real cl = eos.SoundSpeed(wli);
Real cr = eos.SoundSpeed(wri);
el = wli[IPR] * igm1 +
0.5 * wli[IDN] * (SQR(wli[IV1]) + SQR(wli[IV2]) + SQR(wli[IV3]));
er = wri[IPR] * igm1 +
0.5 * wri[IDN] * (SQR(wri[IV1]) + SQR(wri[IV2]) + SQR(wri[IV3]));
Real rhoa = .5 * (wli[IDN] + wri[IDN]); // average density
Real ca = .5 * (cl + cr); // average sound speed
Real pmid = .5 * (wli[IPR] + wri[IPR] + (wli[IV1] - wri[IV1]) * rhoa * ca);

//--- Step 3. Compute sound speed in L,R
Real ql, qr;
ql = (pmid <= wli[IPR])
? 1.0
: (1.0 + (gamma + 1) / std::sqrt(2 * gamma) * (pmid / wli[IPR] - 1.0));
qr = (pmid <= wri[IPR])
? 1.0
: (1.0 + (gamma + 1) / std::sqrt(2 * gamma) * (pmid / wri[IPR] - 1.0));

//--- Step 4. Compute the max/min wave speeds based on L/R states

al = wli[IV1] - cl * ql;
ar = wri[IV1] + cr * qr;

Real bp = ar > 0.0 ? ar : 0.0;
Real bm = al < 0.0 ? al : 0.0;
//--- Step 2. Compute Roe-averaged state
Real sqrtdl = std::sqrt(wli[IDN]);
Real sqrtdr = std::sqrt(wri[IDN]);
Real isdlpdr = 1.0 / (sqrtdl + sqrtdr);

wroe[IDN] = sqrtdl * sqrtdr;
wroe[IV1] = (sqrtdl * wli[IV1] + sqrtdr * wri[IV1]) * isdlpdr;
wroe[IV2] = (sqrtdl * wli[IV2] + sqrtdr * wri[IV2]) * isdlpdr;
wroe[IV3] = (sqrtdl * wli[IV3] + sqrtdr * wri[IV3]) * isdlpdr;

// Following Roe(1981), the enthalpy H=(E+P)/d is averaged for adiabatic flows,
// rather than E or P directly. sqrtdl*hl = sqrtdl*(el+pl)/dl = (el+pl)/sqrtdl
const Real el = wli[IPR] * igm1 +
0.5 * wli[IDN] * (SQR(wli[IV1]) + SQR(wli[IV2]) + SQR(wli[IV3]));
const Real er = wri[IPR] * igm1 +
0.5 * wri[IDN] * (SQR(wri[IV1]) + SQR(wri[IV2]) + SQR(wri[IV3]));
const Real hroe = ((el + wli[IPR]) / sqrtdl + (er + wri[IPR]) / sqrtdr) * isdlpdr;

//--- Step 3. Compute sound speed in L,R, and Roe-averaged states

const Real cl = eos.SoundSpeed(wli);
const Real cr = eos.SoundSpeed(wri);
Real q = hroe - 0.5 * (SQR(wroe[IV1]) + SQR(wroe[IV2]) + SQR(wroe[IV3]));
const Real a = (q < 0.0) ? 0.0 : std::sqrt(gm1 * q);

//--- Step 4. Compute the max/min wave speeds based on L/R and Roe-averaged values
const Real al = std::min((wroe[IV1] - a), (wli[IV1] - cl));
const Real ar = std::max((wroe[IV1] + a), (wri[IV1] + cr));

Real bp = ar > 0.0 ? ar : TINY_NUMBER;
Real bm = al < 0.0 ? al : TINY_NUMBER;

//-- Step 5. Compute L/R fluxes along lines bm/bp: F_L - (S_L)U_L; F_R - (S_R)U_R
Real vxl = wli[IV1] - bm;
Expand Down