Centralise the equation-of-state expressions in the Riemann solvers - #1762
Open
sbryngelson wants to merge 6 commits into
Open
Centralise the equation-of-state expressions in the Riemann solvers#1762sbryngelson wants to merge 6 commits into
sbryngelson wants to merge 6 commits into
Conversation
s_compute_speed_of_sound took H, |u|^2 and qv and then undid them: for a real state H = ((Gamma+1)p + Pi + qv)/rho + |u|^2/2, so c^2 = (H - |u|^2/2 - qv/rho)/Gamma reduces to ((Gamma+1)p + Pi)/(Gamma rho). The three arguments cancel. Callers therefore no longer supply them, and the invisible 'H must include qv' contract cannot be stated. That contract is what produced MFlowCode#1707: five sites open-coded H and three dropped the qv this routine went on to subtract. Those three are fixed here by construction. An average of two states is not a state - its enthalpy is a free input - so the four interface-averaged sites use s_compute_speed_of_sound_avg, whose arithmetic is unchanged. Falls out as dead: qv_sf (a full-domain array whose only reader was one of the open-coded enthalpies), the H dummy threaded through s_save_data and p_main, and a per-cell qv accumulation in post-process. Claude-Session: https://claude.ai/code/session_011BUQDw64EtzvDzTvtWTNj4
|
Claude Code Review Head SHA: 9974b68 Files changed:
Findings:
|
sbryngelson
marked this pull request as ready for review
August 25, 2026 15:56
sbryngelson
requested
a lite review from Copilot
and removed request for
Copilot
August 25, 2026 15:56
Contributor
There was a problem hiding this comment.
Pull request overview
Derives real-state sound speed directly from thermodynamic state while preserving enthalpy-based calculations for interface averages.
Changes:
- Simplifies real-state sound-speed calls across simulation and post-processing.
- Adds a dedicated interface-average sound-speed routine.
- Removes obsolete enthalpy threading and
qv_sfstorage.
Reviewed changes
Copilot reviewed 12 out of 12 changed files in this pull request and generated 1 comment.
Show a summary per file
| File | Description |
|---|---|
src/common/m_variables_conversion.fpp |
Implements split sound-speed APIs and removes qv_sf. |
src/simulation/m_time_steppers.fpp |
Uses state-derived sound speed. |
src/simulation/m_riemann_solver_lf.fpp |
Simplifies LF sound-speed calls. |
src/simulation/m_riemann_solver_hypo_hlld.fpp |
Simplifies hypoelastic HLLD calls. |
src/simulation/m_riemann_solver_hlld.fpp |
Simplifies HLLD calls. |
src/simulation/m_riemann_solver_hllc.fpp |
Separates real and averaged-state calculations. |
src/simulation/m_riemann_solver_hll.fpp |
Separates real and averaged-state calculations. |
src/simulation/m_data_output.fpp |
Corrects diagnostic sound-speed computation. |
src/simulation/m_cbc.fpp |
Removes unnecessary enthalpy calculation. |
src/post_process/p_main.fpp |
Removes obsolete enthalpy argument. |
src/post_process/m_start_up.fpp |
Simplifies sound-speed output generation. |
src/post_process/m_data_output.fpp |
Removes redundant qv and enthalpy accumulation. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| ! Compute mixture sound Speed | ||
| call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, ((gamma + 1._wp)*pres + pi_inf)/rho, alpha, 0._wp, & | ||
| & 0._wp, c, qv) | ||
| call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, alpha, c) |
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #1762 +/- ##
==========================================
+ Coverage 61.67% 61.71% +0.04%
==========================================
Files 84 84
Lines 21619 21596 -23
Branches 3196 3197 +1
==========================================
- Hits 13334 13329 -5
+ Misses 6093 6074 -19
- Partials 2192 2193 +1 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
s_accumulate_mixture_properties and the else branch of s_convert_species_to_mixture_variables_kernel contained the same four-line accumulation loop. Move the routine from m_riemann_state.fpp into m_variables_conversion.fpp and have the kernel call it, below the mpp_lim clipping so it still sees the clipped volume fractions. The two routines are not merged. s_accumulate_mixture_properties is a parameterised subset accumulator - callers pass num_fluids or num_fluids - 1, and raw or limited volume fractions - while the kernel always spans num_fluids, clips in place, special-cases num_fluids == 1 with bubbles_euler, and optionally emits Re_K and G_K. Merging would need a clipping flag and optional dummies on a [seq] device routine, which is not portable across the offload backends. This is the single place per-fluid EOS dispatch will enter the stiffened-gas mixture path. Claude-Session: https://claude.ai/code/session_011BUQDw64EtzvDzTvtWTNj4
Adds s_compute_energy to m_variables_conversion: the stiffened-gas total energy, thermodynamic terms only. Magnetic and elastic energies stay at their call sites because they are not equation-of-state terms, and the chemistry and relativistic forms are left open-coded because they are not this relation. Converts the ten plain sites in hll, hllc, hypo_hlld and lf. The four MHD sites and the two hllc bubbles sites that omit qv are deliberately left for follow-up commits; the latter changes behaviour. Also deletes a dead energy assignment in m_cbc: E is read only inside the chemistry branch of the flux update, so the non-chemistry assignment was never used. Measured on MI210 / amdflang: all 60 Riemann kernel resource profiles identical (scratch, VGPR, AGPR); 5eq_rk3_weno3_hllc 2.290 / 2.385 / 2.407 against a base of 2.374 / 2.422 / 2.361. Claude-Session: https://claude.ai/code/session_011BUQDw64EtzvDzTvtWTNj4
hll, hllc, hlld, hypo_hlld and lf each derived the stiffened-gas mixture coefficients their own way: two called a shared accumulator, two hand-rolled the loop inline, and hllc carried a three-way branch. They now all call s_compute_mixture_coefficients, which owns the rule including its one special case. That special case is required, not incidental. Under bubbles_euler with num_fluids == 1 the sole advection slot aliases the void fraction (eqn_idx%alf == eqn_idx%adv%end), so alpha is not a composition there and the mixture rule does not apply; the coefficients are the liquid's. Clipping stays with the callers because it is genuinely not uniform - hll clips local arrays while hllc clips the shared reconstruction buffers in place - and it can never coincide with the special case, since case_validator prohibits mpp_lim with num_fluids == 1. Behaviour fix, not a refactor: hll and lf previously lacked the special case and computed gamma = alf*gammas(1) under bubbles_euler, which is physically meaningless. No test or example exercises that combination, but it runs today, so results change for anyone using it. hllc bubbles at num_fluids == 2 likewise moves from pure-liquid to standard accumulation; also unexercised. s_accumulate_mixture_properties is merged in and deleted: after hllc's three-way branch collapsed, its num_fluids - 1 call site was gone and it had exactly one caller. The gpuParallelization.md example that used it is updated. s_compute_energy now takes composition rather than coefficients, deriving them through the same rule, so no call site holds equation-of-state knowledge. Measured on MI210 / amdflang: hypo_hlld drops 24 VGPR and 24 AGPR because the caller no longer keeps gamma, pi_inf and qv live across the call; hllc and lf gain 8 VGPR at 118-142, far from any occupancy limit. 5eq_rk3_weno3_hllc 2.343 / 2.319 / 2.330 against 2.304 / 2.283 / 2.366 before. Claude-Session: https://claude.ai/code/session_011BUQDw64EtzvDzTvtWTNj4
sbryngelson
changed the title
The operator returns the thermodynamic energy; pres_mag is added at the call site, because magnetic energy is not an equation-of-state term and a second EOS backend must not have to know about it. Also trims the operator doc comments to the facts a reader cannot derive from the code.
The bubbles_euler branch computed E without qv while every other energy site in the file includes it, and s_compute_pressure's bubbles branch subtracts qv when inverting. Forward and inverse transforms disagreed, so E was inconsistent with the pressure it came from whenever fluid_pp(1)%qv is non-zero. E feeds the energy flux directly, so the error was not confined to diagnostics. Routing the site through s_compute_energy restores qv. No test or example combines bubbles_euler with a non-zero qv, so the suite neither demonstrates the bug nor the fix; the argument is the forward/inverse inconsistency above. Measured on MI210 / amdflang: HLLC drops 226 to 214 VGPR; 5eq_rk3_weno3_hllc 2.339 / 2.343 / 2.283.
sbryngelson
changed the title
Lines of Code
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Centralizes three equation-of-state expression families so that adding a second EOS (#1638, JWL first)
means adding a branch in one routine rather than editing every solver. Supersedes the approach in
#1714, which severely regressed GPU performance. Closes the enthalpy/sound-speed half of #1708 and
subsumes #1707.
Each family now has exactly one implementation, and no Riemann solver holds equation-of-state
knowledge: the energy operator takes composition rather than coefficients.
Hsupplied by the callers_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c); averaged states use_avgs_compute_mixture_coefficientss_compute_energy(pres, alpha_rho, adv, vel_sum, E)This PR contains a behavior fix, not only refactoring.
hllandlflacked thebubbles_eulerspecial case and computed
gamma = alf*gammas(1), which is physically meaningless -alfis a voidfraction, because under
bubbles_eulerwithnum_fluids == 1the sole advection slot aliases it.No test or example exercises that combination, but it runs today, so results change for anyone using
it.
hllcbubbles atnum_fluids == 2likewise moves from pure-liquid to standard accumulation.Supersedes the approach in #1714, which severely degrades GPU performance. Closes the enthalpy/sound-speed half of #1708 and subsumes #1707.
The observation
s_compute_speed_of_soundtookH,|u|^2andqv, and then undid them. For a real stateso the stiffened-gas relation
c^2 = (H - |u|^2/2 - qv/rho)/Gammareduces toEnthalpy, velocity, and heat of formation cancel identically. The diagnostic call sites make this vivid: they passed
vel_sum = 0and aqv, and the routine immediately subtracted it back off.Those three arguments were never needed, and asking for them produced #1707 — five sites open-coded
H, and three dropped theqv. Here there is nothing to open-code, so the defect cannot be written. That is strictly stronger than a type: #1714'seos_stateenforced nothing at 16 of its 24 sites, by its own admission.The change
s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)for real states — 22 of 26 call sites.s_compute_speed_of_sound_avg(...)keeps the enthalpy for the 4 interface-averaged sites, where it genuinely is a free input (a Roe or arithmetic average of two states is not a state). Its arithmetic is unchanged from master.1 + (Gamma + 1)p/rhomoves inside the routine, so relativistic sites need no enthalpy either.Falls out as dead:
qv_sf(a full-domain array with halos whose only reader was one of the open-coded enthalpies), theHdummy threaded throughs_save_dataandp_main, and a per-cellqvaccumulation in post-process.Why #1714 regresses, measured
One MI210 (gfx90a), AFAR amdflang, OpenMP offload,
5eq_rk3_weno3_hllc, same node/compiler/case/rank count throughout.valuePutting a derived type in a
declare target [seq]signature forces a memory ABI on an address-taken aggregate, and amdflang's register promotion collapses across the whole enclosing loop body — roughly two-thirds of occupancy. It is not the constructor call (#1714's second commit removed that and recovered nothing) and not the 192 bytes of struct copies (the effect is ~25x larger). Declaring the dummyvaluedoes not help: the code object is byte-identical.Kernel numbers are read from the AMDGPU code-object metadata, so they are properties of the binary, not the board.
Verification
chemistryxc_cxrelativityxalt_soundspeedx6eqx5eq+bubblesselect master's branch in both new routines; the only arithmetic that changes anywhere is the stiffened-gas branch at real-state sites../mfc.sh formatclean,./mfc.sh precheck7/7.The closed form is also more accurate: the old expression is a cancellation of nearly equal numbers. Over 200k randomized states, it agrees to ~1e-16 in ordinary regimes, but at kinetic/internal ~5.7e7 the old form loses ~7e-9 relative accuracy in
c^2while the closed form is exact.The measurements above predate the rebase onto current master, which added two
m_riemann_solver_hypo_hlldcall sites; the suite result and the build are from the rebased commit.Second commit: a single home for the mixture accumulation
s_accumulate_mixture_properties(inm_riemann_state.fpp, used by the Riemann solvers) and theelsebranch ofs_convert_species_to_mixture_variables_kernel(inm_variables_conversion.fpp,used by IBM, CBC, the CFL calculation, Lagrangian bubbles and the conservative-to-primitive
conversion) contained the same four-line accumulation of the stiffened-gas mixture rule. The routine
moves into
m_variables_conversion.fppand the kernel calls it, below thempp_limclipping so itstill sees the clipped volume fractions.
The two routines are not merged, deliberately:
s_accumulate_mixture_propertiesis aparameterised subset accumulator (callers pass
num_fluidsornum_fluids - 1, and raw or limitedvolume fractions), while the kernel always spans
num_fluids, clips in place, special-casesnum_fluids == 1 .and. bubbles_euler, and optionally emitsRe_K/G_K. Merging would need aclipping flag and optional dummies on a
[seq]device routine, which is not portable across theoffload backends.
Measured on the same commit base, MI210 / amdflang / OpenMP offload:
s_convert_conservative_to_primitive_variables— 540 B scratch, 118 VGPR, unchanged5eq_rk3_weno3_hllc, three runs each: base 2.453 / 2.279 / 2.377 (mean 2.370), with this commit2.374 / 2.422 / 2.361 (mean 2.386). +0.68%, against a base spread of 7.3% — noise.
Not in scope
s_compute_enthalpynow has no consumer for its enthalpy at any of its three callers.c_avgis consumed only underwave_speeds = pressure; in the default configuration, every solver computes it per face and discards it.gamma/pi_inf, which are stiffened-gas parameters. That needs the same question answered for energy, pressure inversion, and enthalpy at once, and is scoped as follow-up work.