Skip to content

Central thermodynamics interface#1663

Closed
fahnab666 wants to merge 6 commits into
MFlowCode:masterfrom
fahnab666:feature/thermodynamics-interface
Closed

Central thermodynamics interface#1663
fahnab666 wants to merge 6 commits into
MFlowCode:masterfrom
fahnab666:feature/thermodynamics-interface

Conversation

@fahnab666

@fahnab666 fahnab666 commented Jul 21, 2026

Copy link
Copy Markdown

What this PR does

Establishes one central equation-of-state interface (m_thermodynamics) so no solver
path carries its own EOS algebra. No new physics and no user-facing EOS choice yet.
Existing results are bit-identical.

The six operations, all through one interface

Operation Entry point Callers routed
pressure(state) s_compute_pressure conversion cons->prim, data_output probes, start_up
internal_energy(rho, p, composition) s_compute_internal_energy conversion prim->cons
temperature(state) f_compute_temperature chemistry cons->prim
sound_speed(state) s_compute_speed_of_sound Riemann (HLL, HLLC), conversion
thermodynamic_derivatives(state) f_bulk_modulus (rho c^2) alt_soundspeed, RHS K-div-u, post-process derived sound speed
validate_state(state) f_validate_state the two mixture_err sound-speed guards (interface and post-process)

Two adapters

Compile-time chemistry selects the adapter:

  • stiffened_gas adapter: the analytic mixture EOS (non-chemistry branch).
  • pyrometheus adapter: temperature-then-query against the already-generated
    m_thermochem, selected automatically when chemistry is enabled.

GPU device routines cannot use runtime polymorphism, so the adapter is the compile-time
branch, not a per-operation object.

validate_state and the GPU constraint

Device code cannot abort (no s_mpi_abort on the device), so validate_state is a pure
predicate that returns a flag rather than asserting. It centralizes the one runtime
state-validity test that was already scattered: the c^2 >= 0 sound-speed admissibility,
which lived inline in both the interface and the post-process derived path. Each site
keeps its own floor value, so .not. f_validate_state(c) is bit-identical to the previous
c < 0. Later adapters (Mie-Gruneisen, JWL) tighten the admissible domain in this one
predicate.

Coupling code left inline, documented

Of the coupling sites the spec names, the post-process derived sound speed is routed. The
rest stay inline because routing them would change answers, not merely relocate arithmetic:

  • IBM ghost and fresh-cell energy (s_ibm_correct_state): the bubble branch weights by
    (1 - alpha) inside the energy, the 5-equation branch drops the qv term, and its
    pres_IP is real(wp) where the interface argument is real(stp).
  • Lagrangian-bubble host state (s_compute_cson_from_pinf): omits -qv from the liquid
    pressure and -qv/rho from the sound speed, and carries its own mixture_err floor.

These are per-site formula differences, not a common interface call. Routing them is a
behavior change and belongs in a later PR. Some of the qv omissions may be latent bugs
for cases with nonzero qv; they are flagged here but not fixed, to keep this PR
bit-identical.

Testing

  • Full golden suite bit-identical on pre_process, simulation, and post_process. No golden
    file was regenerated.
  • --debug bounds-checked suite (the CI condition that was failing): 615 of 615 passed.
  • The earlier out-of-bounds read in the prim->cons energy call (void-fraction index for
    non-bubble configurations) is fixed and confirmed under bounds-checking.

Related

Part of the incremental EOS ladder tracked in #1638. Ladder so far: PR 1 #1663 (central thermodynamics interface), PR 2 #1664 (explicit EOS selector), PR 3 #1665 (generic Mie-Gruneisen core); JWL follows as a Mie-Gruneisen reference curve. The three open PRs are independent and can be reviewed in parallel.

Introduce m_thermodynamics as the single equation-of-state interface and move s_compute_pressure into it unchanged. m_variables_conversion re-exports the routine, so all ten call sites are untouched. The leaf module uses only global parameters and the thermochemistry surface, so no module cycle forms. The stiffened-gas and Pyrometheus branches are preserved verbatim.

No changed answers: the branch-representative golden subset (stiffened gas, MHD, bubbles, hypoelastic, phase change, chemistry) passes with zero golden regeneration across all three targets.

First increment of PR1. Sound speed and the remaining query sites follow.
Move s_compute_speed_of_sound into m_thermodynamics unchanged and re-export it through m_variables_conversion, matching the pressure op. To keep the interface a leaf module, relocate the pre and post process declaration of the per-fluid property tables (gammas, gs_min, pi_infs, ps_inf, cvs, qvs, qvps) from m_variables_conversion into m_global_parameters_common; the simulation target already keeps them in global parameters. Allocation, population, and deallocation are unchanged.

No changed answers: the full golden suite passes with zero regeneration (616 of 616 runnable cases) across all three targets, including the sound-speed branches (Wood alt_soundspeed, six-equation, bubble, relativistic, chemistry).

Second increment of PR1. Internal energy and temperature query sites follow.
Add s_compute_internal_energy to m_thermodynamics as the inverse of
s_compute_pressure over the same model branches (MHD, five-equation,
bubble, four-equation, chemistry). The primitive-to-conservative energy
build in m_variables_conversion now delegates to it, so the conversion
routine no longer carries its own EOS energy algebra. Arithmetic and the
original branch comments are relocated unchanged; the full golden suite
is bit-identical (616 passed, 0 failed).
…modynamics

Add f_bulk_modulus, the stiffened-gas thermodynamic derivative
((gamma+1)p+pi_inf)/gamma = rho*c^2, as an inlinable pure function in the
central thermodynamics interface. Route the six identical inline sites to it:
the two alt_soundspeed Wood-mixture branches in s_compute_speed_of_sound, the
two five-equation K-div-u blkmod sites in m_rhs, and the two post-process
derived sound-speed blkmod sites. The fused /(gamma*rho) post-process branch is
left inline (different arithmetic, not this derivative). Bit-identical: 616
passed, 0 failed.
@fahnab666 fahnab666 changed the title Central thermodynamics interface (no changed answers) Central thermodynamics interface Jul 21, 2026
@fahnab666
fahnab666 marked this pull request as ready for review July 21, 2026 01:59
@fahnab666
fahnab666 requested a review from sbryngelson as a code owner July 21, 2026 01:59
Copilot AI review requested due to automatic review settings July 21, 2026 01:59

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

Refactors equation-of-state (EOS) and related thermodynamic computations into a new leaf module (src/common/m_thermodynamics.fpp), with the goal of routing “genuine thermodynamic queries” through a single interface while preserving bit-identical results.

Changes:

  • Introduces m_thermodynamics as the central interface for pressure, internal energy, sound speed, and bulk-modulus derivative.
  • Replaces duplicated inline bulk-modulus expressions with f_bulk_modulus in simulation RHS and post-process derived-variable paths.
  • Moves per-fluid property table declarations for non-simulation targets into m_global_parameters_common to support the new leaf-module dependency structure.

Reviewed changes

Copilot reviewed 6 out of 6 changed files in this pull request and generated no comments.

Show a summary per file
File Description
src/simulation/m_rhs.fpp Uses f_bulk_modulus to replace duplicated bulk-modulus algebra in the RHS path.
src/post_process/m_derived_variables.fpp Uses f_bulk_modulus to replace duplicated bulk-modulus algebra in derived sound-speed computation.
src/common/m_variables_conversion.fpp Re-exports thermodynamic routines from m_thermodynamics and routes internal-energy construction through the new interface.
src/common/m_thermodynamics.fpp New central thermodynamics/EOS leaf module providing pressure, energy inversion, sound speed, and bulk modulus derivative.
src/common/m_global_parameters_common.fpp Declares per-fluid property tables for pre/post targets to avoid module cycles and support the new interface.
docs/module_categories.json Registers m_thermodynamics in the documentation module categorization list.

@fahnab666
fahnab666 marked this pull request as draft July 21, 2026 02:05
@fahnab666
fahnab666 marked this pull request as ready for review July 21, 2026 02:54
…nfigs

eqn_idx%alf is 0 outside euler-bubble models (model_eqns=3 in pre/post_process, where the assignment is MFC_SIMULATION-only). PR1 hoisted q_prim_vf(eqn_idx%alf) into an unconditional actual argument, so bounds-checked reldebug CI traps on q_prim_vf(0). The value is read only in the routine's bubble branch, so clamp the index with max(1, ...); results are bit-identical.
Add the two remaining operations so all six live in m_thermodynamics:
f_compute_temperature (ideal-gas T from pressure, routed from the chemistry
cons->prim path) and f_validate_state (sound-speed admissibility predicate).
The two scattered mixture_err sound-speed guards, in the interface and in
post-process, now consult f_validate_state instead of testing c<0 inline, and
the post-process derived sound speed routes through f_bulk_modulus like the
alt_soundspeed branch. Bit-identical: floor values are unchanged and the
temperature routing keeps operand order.

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.

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

@codecov

codecov Bot commented Jul 21, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 70.73171% with 24 lines in your changes missing coverage. Please review.
✅ Project coverage is 59.39%. Comparing base (0857ace) to head (c861a5c).

Files with missing lines Patch % Lines
src/common/m_thermodynamics.fpp 75.36% 13 Missing and 4 partials ⚠️
src/post_process/m_derived_variables.fpp 0.00% 4 Missing ⚠️
src/common/m_variables_conversion.fpp 57.14% 1 Missing and 2 partials ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master    #1663      +/-   ##
==========================================
- Coverage   59.46%   59.39%   -0.07%     
==========================================
  Files          83       84       +1     
  Lines       21170    21172       +2     
  Branches     3135     3135              
==========================================
- Hits        12588    12575      -13     
- Misses       6476     6489      +13     
- Partials     2106     2108       +2     

☔ 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.

@sbryngelson
sbryngelson marked this pull request as draft July 21, 2026 12:32
@sbryngelson

Copy link
Copy Markdown
Member

is it fast

@sbryngelson

Copy link
Copy Markdown
Member

does it scale

@sbryngelson

sbryngelson commented Jul 21, 2026

Copy link
Copy Markdown
Member

i don't trust the performance of this branch or its contents, unfortunately

@fahnab666

fahnab666 commented Jul 21, 2026

Copy link
Copy Markdown
Author

@sbryngelson Fair questions, so I measured with the repository's own harness rather than argue. Summary: on CPU this branch is performance neutral on every benchmark case, and the diff itself shows why it cannot slow the hot path.

Setup. ./mfc.sh bench --mem 1 on the merge base (0857ace) and this branch (c861a5c), identical build locks (Release, MPI on, no fastmath), single node, Apple M4 Pro, gfortran 15.2.0. These are CPU runs. Grind time is ns per grid point per time step, lower is better.

Case base grind PR grind delta
5eq_rk3_weno3_hll 69.32 69.53 +0.31%
5eq_rk3_weno3_hllc 72.01 72.36 +0.48%
5eq_rk3_weno3_lf 67.77 68.05 +0.41%
hypo_hll 60.40 61.29 +1.47%*
ibm 281.43 281.88 +0.16%
igr 142.33 142.35 +0.02%
viscous_weno5_sgb_acoustic 194.82 194.68 -0.07%

*The one case above 1% was re-run four more times alternating branches: base 60.63 and 61.10 against PR 60.58 and 60.93, so the PR side was marginally faster in both pairs and the single-shot delta was machine noise.

Scaling. Same HLLC case at 1, 2, and 4 ranks with fixed memory per rank: deltas of +0.4%, -0.3%, and -0.8%, with parallel efficiency at 4 ranks of 0.91 (base) versus 0.93 (PR). The branch changes no MPI code, so parity is expected; this makes it concrete.

Why the diff cannot touch the hot path.

  • No Riemann solver file is modified. The flux kernels are byte-for-byte the code on master.
  • The per-cell pressure call in the conservative-to-primitive loop already existed as a call on master, with the same Cray NOINLINE directive. This branch only moves the callee into m_thermodynamics; on Cray the generated call is identical by construction.
  • s_compute_speed_of_sound was already a cross-module call for every hot caller (the Riemann solvers) before this branch, so relocating it changes nothing.
  • The only inline algebra converted to calls is f_bulk_modulus in the alt_soundspeed branch of m_rhs and in post_process, both outside the benchmarked hot loops. I have a two-line follow-up marking f_bulk_modulus and f_validate_state cray_inline so even those sites compile to the same inlined code as master on Cray; it builds clean with the golden subset bit-identical, and I will push it here if you reopen.

The one thing I cannot exercise locally is the GPU lanes. If you are open to reopening this PR, I can ran then on Tuolumne.

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.

3 participants