Skip to content

Pass an eos_state to s_compute_speed_of_sound - #1714

Open
sbryngelson wants to merge 2 commits into
masterfrom
refactor/eos-state
Open

Pass an eos_state to s_compute_speed_of_sound#1714
sbryngelson wants to merge 2 commits into
masterfrom
refactor/eos-state

Conversation

@sbryngelson

@sbryngelson sbryngelson commented Aug 9, 2026

Copy link
Copy Markdown
Member

Closes the enthalpy half of #1708. Also subsumes #1707.

The problem

s_compute_speed_of_sound took ten positional arguments, and the contract between two of them was invisible at the call site: H must include qv, because the routine subtracts qv/rho internally (m_variables_conversion.fpp, default stiffened-gas branch).

Nothing stated that. Of the five sites that open-coded H:

The change

type(eos_state) in m_derived_types carries the scalars the routine needs, with two ways to build one:

The call becomes (state, adv, c). All 24 call sites across nine files are converted, and the state is added to the private() clause of every GPU parallel loop that builds one.

adv stays a separate argument rather than a component. A derived-type component cannot have a runtime extent, and num_fluids is a parameter only under case optimization, so dimension(num_fluids) does not compile in a general build:

Error: Variable 'num_fluids' cannot appear in the expression at (1)

Padding to num_fluids_max would compile but put ten reals into a per-cell private struct on device.

The safety guarantee is uneven — read this before approving

There are three ways to build a state, and they do not all enforce the invariant:

path sites enforces H includes qv?
s_eos_state 4 yes — derives H, cannot disagree
s_eos_state_roe 4 no — caller supplies H, which is the point for Roe paths
intrinsic eos_state(...) 16 no — every component is caller-supplied

The 16 hot Riemann sites use the intrinsic form for performance (below), so at those sites a wrong H is still writable. What makes them safe in practice is that they all pass an H derived from the conserved energy (H_L = (E_L + pres_L)/rho_L), which carries qv implicitly and cannot drift — the same property that kept them correct before this PR.

So the accurate claim is narrower than "the bug is now impossible": the defect is unrepresentable at the sites that had it (#1707, all diagnostic, all now using s_eos_state), and the sites that never had it keep the construction that made them correct. If you would rather trade the ~5% on HLLC for uniform enforcement, switching those 16 back to s_eos_state_roe is a one-line change per site.

A related hazard worth knowing: the type's component order is rho, pres, ... while the procedure constructors take pres, rho, .... The intrinsic sites are positional and follow the type. Getting that backwards compiles cleanly and silently swaps two doubles; it is caught only by the goldens.

Performance

The first cut cost +5–8% on the Riemann benchmarks. Profiling showed it was not inlining loss — s_compute_speed_of_sound is a real, non-inlined call both before and after. It was the constructor call added on top of it. Object counts for m_riemann_solver_hllc (arm64, GNU 15.2):

build total bl relocs -> speed_of_sound relocs -> eos_state_roe
baseline 179 27 0
constructor call 206 27 27
structure constructor 179 27 0

Nine call sites x 3 relocations = exactly the 27 added calls, one per site. Building the state at the 16 hot Riemann sites with Fortran's intrinsic structure constructor builds it in place and restores the baseline call count exactly. The procedure constructors remain at the diagnostic sites, where deriving H is the point and call cost is irrelevant.

Benchmarks after that fix (idle machine, versus the mean of two baseline runs; measured baseline noise +/-3.5%):

case delta
5eq_rk3_weno3_hll -0.12%
5eq_rk3_weno3_lf -2.12%
hypo_hll -1.00%
ibm -0.57%
igr -0.09%
viscous_weno5_sgb_acoustic -1.26%
5eq_rk3_weno3_hllc +5.63%

Six of seven sit at or below baseline. HLLC is a genuine outstanding cost — consistent across two independent runs and above both baseline samples, though the baseline spread there is wide enough that the magnitude is soft (call it +4-7%). The call count is identical to baseline, so it is not call overhead; the likely cause is materialising the state itself, since HLLC has nine call sites, more than any other solver. Reviewers should weigh that against the interface guarantee.

Relationship to #1707

The three probe sites now use s_eos_state, so the #1707 correction falls out of the interface rather than being a separate patch. That changes the reported probe sound speed where qv /= 0. No golden observes it, because the packer keeps only the last column of probe output (#1711, fixed in #1712).

Verification

build (GNU 15.2, MPI):   clean
full suite:              627 passed, 0 failed
goldens regenerated:     none
precheck (7/7):          pass

A scan confirms no GPU parallel loop uses a state without declaring it private.

Still open in #1708

This covers the enthalpy and sound-speed half. The four open-coded copies of the pressure inversion remain (#1709 fixes the one that was wrong); routing them through a single definition is follow-up work.


Handover notes

Branch and commits

repo    MFlowCode/MFC
branch  refactor/eos-state          (base: master, rebased after #1713 merged)
  936c676  Pass an eos_state to s_compute_speed_of_sound
  102c54c  Build hot-path eos_states with the intrinsic structure constructor

Originally built on top of #1713 and rebased with
git rebase --onto origin/master <removal-sha> refactor/eos-state
once that merged. Re-verified after the rebase: build clean, 627/0.

Environment

./mfc.sh build -j 16
./mfc.sh test  -j 12
./mfc.sh bench --mem 1 -o out.yaml     # ~12 min, 7 cases

GNU 15.2.0, MPI on, no GPU, macOS arm64 (18 cores). This branch has never been built for GPU or with case optimization — both are worth checking in CI, since the type is used inside [seq] device routines and private() clauses.

Where things live

what where
type eos_state src/common/m_derived_types.fpp, next to qbmm_idx_info
s_eos_state, s_eos_state_roe src/common/m_variables_conversion.fpp, immediately before s_compute_speed_of_sound
converted routine same file, s_compute_speed_of_sound(s, adv, c)
call sites 24 across m_riemann_solver_{hll,hllc,hlld,lf}, m_cbc, m_data_output (sim + post), m_time_steppers, post_process/m_start_up

Component order is rho, pres, gamma, pi_inf, qv, vel_sum, H, c_c — note the constructors take pres before rho, the type does not. The structure-constructor call sites are positional and follow the type order; getting this backwards compiles fine and silently swaps two doubles.

Two constructors, and when to use which

  • s_eos_state(...) derives H. Use wherever the exact state enthalpy is wanted. This is what makes bug: simulation probe sound speed omits qv from enthalpy, disagreeing with post-process #1707 unrepresentable.
  • s_eos_state_roe(..., H [, c_c]) takes H. Required for Roe-averaged Riemann paths, the chemistry Roe branch, and the relativistic branch, which pass an H that is deliberately not the exact enthalpy. Do not "simplify" these to the derived form — it changes the scheme.
  • At the 16 hot Riemann sites the state is built with Fortran's intrinsic structure constructor (eos_s_L = eos_state(...)) rather than a call, for performance (below). Keep it that way.

Performance: the whole story

First cut cost +5-8% on the Riemann benchmarks. It was not inlining loss — s_compute_speed_of_sound is a real, non-inlined call both before and after (checked with nm -u on the object). It was the constructor call stacked on top.

Reproduce the object-level check:

O=build/staging/cpu-*/CMakeFiles/simulation.dir/fypp/simulation/m_riemann_solver_hllc.fpp.f90.o
objdump -d "$O" | grep -cE '^\s+[0-9a-f]+:.*\bbl\b'    # total call instructions
objdump -r "$O" | grep -c speed_of_sound               # relocations
objdump -r "$O" | grep -c eos_state_roe
build bl relocs speed_of_sound relocs eos_state_roe
baseline 179 27 0
constructor call 206 27 27
structure constructor 179 27 0

Benchmark methodology matters. Run ./mfc.sh bench in the foreground on an idle machine. I once ran it in the background while the pre-commit precheck was running at -j 12 and got a bogus +50% on igr — a case that executes none of the converted code (run_time_info = F in every benchmark case, and no case sets cfl_dt, so both per-cell loops using the state are switched off). Baseline-vs-baseline noise on this machine is +/-3.5%; run the baseline twice before trusting any delta.

Outstanding: HLLC ~+4-7%

Real, consistent across two independent runs, above both baseline samples. Call count is identical to baseline so it is not call overhead; the likely cause is materialising the state at HLLC's nine call sites, more than any other solver. Not yet chased. Ideas, cheapest first:

  1. More baseline samples on HLLC alone to narrow the error bar — the current estimate rests on two samples per side and the baseline spread is wide (68.82 vs 66.42).
  2. Trim the struct. c_c is read only in the chemistry branch and H only in some branches; a smaller state means fewer stores per site.
  3. Check whether -flto or NVHPC/Cray change the picture; this was measured on gfortran only.

Verification performed

build:                clean
full suite:           627 passed, 0 failed
goldens regenerated:  none
precheck (7/7):       pass
GPU private() audit:  no loop uses a state without declaring it private

The private() audit is worth redoing if you add call sites — a CPU run cannot catch a missing private(), it is a silent device race. Script it by finding $:GPU_PARALLEL_LOOP( directives (excluding END_GPU_PARALLEL_LOOP, which contains the same substring and will produce false positives), joining & continuations, and checking every eos_s* used in the body appears in the directive text.

Related work


Purpose in the series

Closes the enthalpy half of #1708 and gives #1700's selector something to dispatch through.

The underlying problem is that MFC has nine hand-written copies of two EOS expressions — four wrong, one right only by a cancellation no reader can see locally. Adding a second equation of state under #1638 multiplies that: every open-coded site is a place the new backend silently does not reach.

This PR attacks the enthalpy/sound-speed half by making the invisible contract (H must include qv) explicit in a type. The pressure-inversion half remains open — four copies at m_pressure_relaxation.fpp, m_data_output.fpp (x2) and m_bubbles_EL.fpp.

Read the "safety guarantee is uneven" section before approving; the enforcement is not uniform across the 24 sites, deliberately.


Working conventions and hazards (shared across this series)

Collected from the work that produced #1705, #1709, #1712, #1713, #1714, #1716. Every one of these cost real time or produced a wrong result before being caught.

Testing

A regression test that cannot fail is worse than no test. Always verify the negative: revert the fix, rebuild, confirm the case fails, restore. Two ways this silently broke here:

  • git stash push -- <file> has nothing to stash once the fix is committed, so the "reverted" run tests the fixed binary and reports a pass identical to a real one. Use git checkout master -- <file>, rebuild, test, then git checkout HEAD -- <file>.
  • Picking the wrong case configuration. A probe case at model_eqns = 3 cannot detect a sound-speed defect, because the six-equation branch of s_compute_speed_of_sound touches neither H nor qv.

The golden packer discards data. toolchain/mfc/packer/pack.py treated every .dat under D/ as <x> [<y> <z>] <value> and kept only the last column of each row. Probe and integral output are multi-column time series, so most columns were never compared (#1711, fixed in #1712). Before asserting that a golden covers something, check it is actually in golden.txt.

Case labels are load-bearing. The golden UUID is crc32(sha1(str(trace))) — the label chain determines the directory name. Renaming a label renames the golden. Avoid ! in labels (history expansion in interactive bash).

Local suite runs are flaky at high -j. Non-reproducible failures appeared on several unrelated branches at -j 1216 (chemistry cases, probe cases) that passed individually and in clean reruns. Re-run before investigating.

Removing parameters or features

Deregistering a parameter breaks things that are not the source tree. Removing pref/rhoref from the registry broke the entire suite because BASE_CFG in toolchain/mfc/test/case.py set them for every case. Also check fp_stability.py, params_tests/mutation_tests.py, and lint fixtures that use real parameter names as examples.

Grep the generated artifacts, not just the sources. A stale TYPED_DECLS entry naming a deleted type survived removal and did not break the build only because the parameter had also left every target's namelist vars, so it was never emitted. Check generated_decls.fpp, generated_constants.fpp, SIM_GPU_DECL_VARS, and the MPI broadcast generators.

Dead-local tell: after removing a block, a local with exactly one remaining occurrence in its file is almost certainly its own declaration. Two occurrences often means declaration plus a private() entry.

Fortran is case-insensitive. A local pRef shadowed the module global pref in the hardcoded-IC files; the read site was spelled pref and looked like a reference to the global. It is not. Confirm scope before concluding a global is live.

GPU

A CPU test run cannot catch a missing private(). It is a silent device race. Audit by hand or by script when adding per-cell state.

Do not match GPU_PARALLEL_LOOP naivelyEND_GPU_PARALLEL_LOOP contains the same substring and will register as a loop start, producing false positives. Exclude it explicitly.

Derived-type components cannot have runtime extents. dimension(num_fluids) in a type fails to compile outside case-optimized builds, where num_fluids is a parameter.

Benchmarking

Run ./mfc.sh bench in the foreground on an idle machine. Running it in the background while a pre-commit precheck ran at -j 12 produced a bogus +50% regression on a case that executes none of the changed code. Baseline noise here is ±3.5%; run the baseline twice before trusting any delta, and sanity-check that the regressing cases actually execute the modified code.

GitHub mechanics

  • --force-with-lease needs an explicit SHA (--force-with-lease=<branch>:<sha>) when the ref has not been fetched in the current clone; the bare form fails with "stale info".
  • gh run view --log-failed can miss the real output entirely. On the Frontier jobs the failing step carried only a non-zero exit while the actual test output lived in a separate Print Logs step that succeeded. Fetch the full log.
  • gh run rerun <id> --failed refuses while the workflow is still running; retry later.
  • Frontier jobs fail spuriously with some regularity — uniform ~15-minute deaths across heterogeneous jobs, and shared-Lustre install races (file INSTALL cannot set modification time ... No such file or directory, exit 143). Check for a real error before assuming a code fault.

s_compute_speed_of_sound took ten arguments, and the contract between two of
them was invisible at the call site: H must include qv, because the routine
subtracts qv/rho internally. Nothing said so, and four call sites got it wrong
(#1707) while a fifth was right only because two separate omissions cancelled.

Introduces type(eos_state) in m_derived_types carrying the scalars the routine
needs, and two constructors:

  - s_eos_state derives H from the other members, so it cannot disagree with
    qv. The defect in #1707 is unrepresentable through this path.
  - s_eos_state_roe takes H explicitly, for the Roe-averaged Riemann paths,
    the chemistry Roe branch and the relativistic branch, all of which pass an
    H that is deliberately not the exact state enthalpy.

The call becomes (state, adv, c) instead of ten positional arguments. All 24
call sites across nine files are converted, and the state is added to the
private() clause of every GPU parallel loop that builds one.

adv stays a separate argument rather than a component. A derived-type
component cannot have a runtime extent, and num_fluids is a parameter only
under case optimization, so dimension(num_fluids) does not compile in a
general build; padding to num_fluids_max would compile but place ten reals in
a per-cell private struct on device.

The three probe sites in simulation/m_data_output that open-coded H now use
s_eos_state, which incorporates the #1707 correction as a consequence of the
interface rather than as a separate patch. That changes the reported probe
sound speed where qv /= 0; no golden observes it, because the packer keeps
only the last column of probe output (#1711, fixed separately in #1712).

Verified: builds clean; full suite 627 passed, 0 failed, no golden
regenerated. A scan confirms no GPU parallel loop uses a state without
declaring it private.
Profiling the +5-8% regression on the Riemann benchmarks showed it was not
inlining loss: both before and after, s_compute_speed_of_sound is a real
non-inlined call. The cost was the constructor call added on top of it.

Object-code counts for m_riemann_solver_hllc (arm64, GNU 15.2):

    baseline                bl 179, relocs speed_of_sound 27, eos_state 0
    s_eos_state_roe call    bl 206, relocs speed_of_sound 27, eos_state 27
    structure constructor   bl 179, relocs speed_of_sound 27, eos_state 0

Nine call sites x 3 relocations each = the 27 added calls, exactly one per
site. Replacing the constructor call at the 16 Riemann sites with Fortran's
intrinsic structure constructor builds the state in place and restores the
baseline call count.

The derived-type constructors remain for the diagnostic sites, where deriving
H from the other members is the point and the call cost is irrelevant.

Benchmarks (idle machine, mean of two baselines): hll -0.12%, lf -2.12%,
hypo_hll -1.00%, ibm -0.57%, igr -0.09%, viscous -1.26%, all within the
measured +/-3.5% baseline noise band. hllc remains +5.63%, consistent across
two independent runs and above both baseline samples; it has nine call sites,
the most of any solver, so the residual is attributed to materialising the
state rather than to call overhead.

Full suite 627 passed, 0 failed.

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Pull request overview

This PR refactors the equation-of-state sound-speed interface to take a single eos_state derived type rather than many positional scalars, making the H/qv contract explicit and preventing the probe/diagnostic bug where H was open-coded without qv.

Changes:

  • Introduces type(eos_state) in src/common/m_derived_types.fpp and adds two constructors: s_eos_state (exact state; derives H) and s_eos_state_roe (caller-supplied H, optional Roe chemistry term).
  • Replaces s_compute_speed_of_sound(pres, rho, ...) with s_compute_speed_of_sound(state, adv, c) in src/common/m_variables_conversion.fpp.
  • Updates callers across simulation/post-process (including GPU loops) to build an eos_state and pass it to the new routine, adding the state to GPU private(...) lists where needed.

Reviewed changes

Copilot reviewed 11 out of 11 changed files in this pull request and generated 5 comments.

Show a summary per file
File Description
src/common/m_derived_types.fpp Adds the eos_state derived type used to bundle EOS scalars for sound-speed computation.
src/common/m_variables_conversion.fpp Adds state constructors and changes s_compute_speed_of_sound to accept eos_state.
src/simulation/m_time_steppers.fpp Updates timestep CFL sound-speed computation to use eos_state.
src/simulation/m_data_output.fpp Updates stability/probe sound-speed paths to construct an eos_state (fixing the prior qv/H contract issue by construction).
src/simulation/m_cbc.fpp Updates CBC sound-speed computation to the new state-based API.
src/simulation/m_riemann_solver_hll.fpp Updates HLL solver sound-speed calls to the new API (intrinsic state construction).
src/simulation/m_riemann_solver_hllc.fpp Updates HLLC solver sound-speed calls to the new API (intrinsic state construction).
src/simulation/m_riemann_solver_hlld.fpp Updates HLLD solver sound-speed calls to the new API (intrinsic state construction).
src/simulation/m_riemann_solver_lf.fpp Updates LF solver sound-speed calls to the new API (intrinsic state construction).
src/post_process/m_start_up.fpp Updates post-process derived-field sound-speed computation to use eos_state.
src/post_process/m_data_output.fpp Updates post-process energy/probe diagnostics sound-speed computation to use eos_state.
Suppressed comments (1)

src/post_process/m_data_output.fpp:1293

  • H is computed here but no longer used (s_eos_state derives H internally). Removing this assignment avoids redundant per-cell arithmetic in a hot loop.
                    call s_eos_state(eos_s, pres, rho, gamma, pi_inf, qv, 0._wp)
                    call s_compute_speed_of_sound(eos_s, adv, c)

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

Comment on lines +131 to +134
!> `H` is the specific total enthalpy and must include `qv`, because the sound-speed relation
!> subtracts `qv/rho`. Build states with f_eos_state so that invariant holds by construction
!> rather than by convention; f_eos_state_roe exists for the Roe-averaged paths, which supply an
!> `H` that is deliberately not the exact state enthalpy.
real(wp), dimension(num_vels) :: vel
real(wp), dimension(num_fluids) :: adv
integer :: i, j, k, l, s !< looping indices
type(eos_state) :: eos_s

! Compute mixture sound speed
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv)
call s_eos_state_roe(eos_s, pres, rho, gamma, pi_inf, qv, vel_sum, H)
call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l)

call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv)
call s_eos_state_roe(eos_s, pres, rho, gamma, pi_inf, qv, vel_sum, H)
Comment thread src/simulation/m_cbc.fpp

! Compute mixture sound speed
call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c, qv)
call s_eos_state_roe(eos_s, pres, rho, gamma, pi_inf, qv, vel_K_sum, H)
@github-actions

github-actions Bot commented Aug 9, 2026

Copy link
Copy Markdown

Claude Code Review

Head SHA: 102c54c

Files changed:

  • 11
  • src/common/m_derived_types.fpp
  • src/common/m_variables_conversion.fpp
  • src/post_process/m_data_output.fpp
  • src/post_process/m_start_up.fpp
  • src/simulation/m_cbc.fpp
  • src/simulation/m_data_output.fpp
  • src/simulation/m_riemann_solver_hll.fpp
  • src/simulation/m_riemann_solver_hllc.fpp
  • src/simulation/m_riemann_solver_hlld.fpp
  • src/simulation/m_riemann_solver_lf.fpp

Findings:

  • src/common/m_variables_conversion.fpp: the old one-line doc comment !> Compute the speed of sound from thermodynamic state variables, supporting multiple equation-of-state models. was left in place above the new s_eos_state subroutine instead of being moved down to s_compute_speed_of_sound. s_eos_state now carries a doc comment that describes the wrong routine (it builds an EOS state, it does not compute a speed of sound), while s_compute_speed_of_sound below it is left with no doc comment at all.
  • src/common/m_derived_types.fpp: the new eos_state doc comment says "Build states with f_eos_state ... f_eos_state_roe exists for the Roe-averaged paths", but the actual constructors added in this PR are named s_eos_state / s_eos_state_roe (subroutines, per the s_<verb>_<noun> convention) — there is no f_eos_state. This is likely to mislead a future contributor searching for the "sanctioned" constructor referenced by the invariant this type is meant to enforce.
  • src/common/m_variables_conversion.fpp (s_eos_state_roe) and its callers in src/simulation/m_cbc.fpp, src/simulation/m_data_output.fpp, src/simulation/m_time_steppers.fpp, src/post_process/m_start_up.fpp: c_c is declared optional on a routine annotated $:GPU_ROUTINE(..., cray_inline=True) and called from inside GPU_PARALLEL_LOOP/GPU_PARALLEL regions with the argument always omitted, relying on present(c_c) inside the device routine. OPTIONAL dummy arguments combined with present() checks inside GPU-offloaded routines are a known trouble spot on some of the CI-gated backends (this codebase already carries several AMD-flang/Cray-specific device workarounds elsewhere for similar cross-TU/offload quirks), so this pattern is worth explicit verification on the OpenACC/OpenMP-offload and AMD flang builds rather than assumed to behave identically to the host path.

@codecov

codecov Bot commented Aug 9, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 84.14634% with 13 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.28%. Comparing base (b40e08f) to head (102c54c).
⚠️ Report is 1 commits behind head on master.

Files with missing lines Patch % Lines
src/simulation/m_data_output.fpp 50.00% 4 Missing ⚠️
src/simulation/m_riemann_solver_hlld.fpp 0.00% 4 Missing ⚠️
src/common/m_variables_conversion.fpp 91.17% 2 Missing and 1 partial ⚠️
src/post_process/m_data_output.fpp 0.00% 2 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master    #1714      +/-   ##
==========================================
+ Coverage   61.24%   61.28%   +0.04%     
==========================================
  Files          83       83              
  Lines       20700    20741      +41     
  Branches     3072     3073       +1     
==========================================
+ Hits        12677    12712      +35     
- Misses       5969     5975       +6     
  Partials     2054     2054              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Development

Successfully merging this pull request may close these issues.

2 participants