Add a Mie-Gruneisen equation-of-state backend - #1805
Conversation
There was a problem hiding this comment.
Warning
Copilot couldn't run its full agentic review because it didn't start before the timeout. Make sure your repository has a runner available, or add a copilot-code-review.yml file specifying one with the runs-on attribute. See the docs for more details.
Pull request overview
Adds a new Mie–Grüneisen EOS backend (linear-Hugoniot reference curve mapped to MFC’s rho e = Gamma p + Pi form), along with toolchain parameter support, validation, and math-focused tests, without wiring it into solver paths.
Changes:
- Introduces EOS selector value
mie_gruneisen+ newfluid_pp%mg_*parameters (toolchain + Fortran derived types + initialization). - Adds
s_eos_coefficientsto compute(Gamma, Pi, dPi/drho)for Mie–Grüneisen while keeping stiffened/ideal behavior unchanged. - Adds a pytest suite to pin the reference-curve derivatives, mapping identity, and sound-speed relation.
Reviewed changes
Copilot reviewed 10 out of 10 changed files in this pull request and generated 4 comments.
Show a summary per file
| File | Description |
|---|---|
| toolchain/mfc/test_eos_mie_gruneisen.py | New manufactured-math tests for the MG backend and mapping identity. |
| toolchain/mfc/params/definitions.py | Registers EOS selector value 3 and new mg_* parameters for fluid_pp. |
| toolchain/mfc/case_validator.py | Validates mg_* presence/absence and constraints; forbids qv for MG. |
| src/simulation/m_global_parameters.fpp | Initializes new fluid_pp%mg_* fields to defaults. |
| src/pre_process/m_global_parameters.fpp | Initializes new fluid_pp%mg_* fields to defaults. |
| src/post_process/m_global_parameters.fpp | Initializes new fluid_pp%mg_* fields to defaults. |
| src/common/m_variables_conversion.fpp | Adds device-side storage of EOS selectors/params and introduces s_eos_coefficients. |
| src/common/m_global_parameters_common.fpp | Adds global arrays for EOS selectors/mg_* params + GPU_DECLARE. |
| src/common/m_derived_types.fpp | Extends fluid_pp derived type with mg_* fields. |
| src/common/m_constants.fpp | Adds eos_mie_gruneisen = 3 constant. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| real(wp) :: mu, d, p_ref, e_ref, dp_dmu, de_dmu, G0 | ||
|
|
||
| select case (eoss(i)) | ||
| case (eos_mie_gruneisen) | ||
| ! Linear-Hugoniot reference curve, u_s = c0 + s u_p: p_H = rho0 c0^2 mu (1 + mu)/(1 - (s - 1) mu)^2 on | ||
| ! compression, extended linearly on release, with the Hugoniot energy e_H = p_H mu/(2 rho0 (1 + mu)). | ||
| ! Pole at mu = 1/(s - 1); the validator warns near it. | ||
| mu = rho/mg_rho0s(i) - 1._wp | ||
| if (mu >= 0._wp) then | ||
| d = 1._wp - (mg_ss(i) - 1._wp)*mu |
There was a problem hiding this comment.
You're right that the comment claimed a warning that didn't exist in this PR - that was a deferred item I left documented as done. Fixed in 3ffde53: the comment now states the singularity plainly as a validity limit (mu = 1/(s - 1), the Hugoniot's maximum compression) with no runtime guard, and the validator gains an initial-state check that warns when any patch starts a Mie-Gruneisen fluid within 20% of rho0*s/(s - 1). A runtime clamp is deliberately not added: a state past the pole is outside the EOS, and clamping would return a wrong pressure silently - the same reasoning as #1793. The solver's per-step NaN check is the runtime backstop, and it only sees the Inf once it has become a NaN; that gap is real and now documented rather than implied to be closed.
| eoss(i) = fluid_pp(i)%eos | ||
| mg_rho0s(i) = fluid_pp(i)%mg_rho0 | ||
| mg_c0s(i) = fluid_pp(i)%mg_c0 | ||
| mg_ss(i) = fluid_pp(i)%mg_s | ||
| mg_G0s(i) = fluid_pp(i)%mg_G0 | ||
| if (fluid_pp(i)%eos == eos_mie_gruneisen) any_state_dependent_eos = .true. |
There was a problem hiding this comment.
Fixed in 3ffde53b. Also dropped the initializer, which made the flag an implicit save; it is now set at the top of init.
| real(wp) :: mg_rho0 !< Mie-Gruneisen reference density | ||
| real(wp) :: mg_c0 !< Mie-Gruneisen bulk sound speed at mg_rho0 | ||
| real(wp) :: mg_s !< Mie-Gruneisen linear Hugoniot slope, u_s = c0 + s u_p | ||
| real(wp) :: mg_G0 !< Mie-Gruneisen coefficient |
There was a problem hiding this comment.
Renamed mg_G0 -> mg_gruneisen in 3ffde53b. The stronger reason than Gamma-vs-Gamma_G: fluid_pp%G is already the shear modulus, so mg_G0 sat one line from an unrelated G.
| "title": "Equation of State Selector", | ||
| "category": "Thermodynamic Constraints", | ||
| "math": r"\Pi_\infty = 0 \;\; \text{for an ideal gas}", | ||
| "math": r"\rho e = \Gamma(\rho)\,p + \Pi(\rho), \quad \Gamma = 1/\Gamma_G, \quad \Pi = \rho e_{\mathrm{ref}} - p_{\mathrm{ref}}/\Gamma_G", |
There was a problem hiding this comment.
Fixed in 3ffde53b: Gamma is constant, Pi(rho) carries the density dependence, and the string now says exactly that.
|
Claude Code Review Head SHA: a25b84d Files changed:
Findings:
|
Any Mie-Gruneisen EOS p = p_ref + rho Gamma_G (e - e_ref) is MFC's rho e = Gamma p + Pi with Gamma = 1/Gamma_G and Pi = rho e_ref - p_ref/Gamma_G, so the existing operators need only the coefficients and dPi/drho. s_eos_coefficients is the single dispatch: the Mie-Gruneisen case supplies a linear-Hugoniot reference curve (u_s = c0 + s u_p, linear release) and one shared conversion produces Gamma, Pi and dPi/drho; a second family is one more case. Stiffened and ideal gas return the constants resolved at init, and no caller invokes the routine yet, so every existing answer is bit-identical - 27-case gate. Parameters mg_rho0, mg_c0, mg_s, mg_G0 follow the flat per-fluid style; the validator requires all four under mie_gruneisen, forbids them otherwise, and forbids qv there because e_ref carries the formation energy. Thirty pytest checks pin the maths against finite differences and a numerically integrated isentrope.
|
The pole finding in the Claude Code Review is addressed in 3ffde53b: the source comment no longer claims a guard that did not exist, and the validator now warns when a patch starts a Mie-Gruneisen fluid within 20% of the Hugoniot pole. Details in the reply to the matching Copilot thread. |
a25b84d to
3ffde53
Compare
Lines of Code
|
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1805 +/- ##
==========================================
- Coverage 62.26% 62.24% -0.03%
==========================================
Files 84 84
Lines 21558 21595 +37
Branches 3188 3196 +8
==========================================
+ Hits 13423 13441 +18
- Misses 5937 5956 +19
Partials 2198 2198 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
First step of the Mie-Gruneisen path in #1638: the backend and its maths, with nothing wired into the
solver yet. Every existing answer is bit-identical.
The mapping
Any Mie-Gruneisen equation of state,
p = p_ref(rho) + rho Gamma_G (e - e_ref(rho)), is already MFC'sform
rho e = Gamma p + Piwithso the operators #1762 centralized need only the two coefficients plus
dPi/drhofor the sound speed.Stiffened gas is the degenerate member,
p_ref = -gamma pi_inf,e_ref = 0,Gamma_G = gamma - 1, whichreproduces the stored
gammas(i)andpi_infs(i)exactly.What is added
s_eos_coefficients(rho, i, gamma, pi_inf, dpi)- one dispatch on the fluid's EOS. The Mie-Gruneisencase supplies a linear-Hugoniot reference curve,
u_s = c0 + s u_p, extended linearly on release; a sharedconversion turns any reference curve into
Gamma,Pi,dPi/drho. Adding JWL is one morecase.Stiffened and ideal gas return the constants resolved at init and never enter the per-cell path.
Parameters
fluid_pp(i)%mg_rho0,mg_c0,mg_s,mg_G0, flat likegammaandpi_inf. The validatorrequires all four under
eos = 'mie_gruneisen', forbids them otherwise, requires positivity ands >= 1,and forbids
qvthere becausee_refalready carries the formation energy.Gamma_Gis constant, which is exact for JWL; thep dGamma/drhoterm is kept in the sound-speedexpression so a density-dependent form later is a one-function change.
Verification
toolchain/mfc/test_eos_mie_gruneisen.py, 30 checks, no solver run: reference-curve derivatives againstcentral differences; the
Gamma/Piform inverting to the Mie-Gruneisen pressure to 1e-12; analyticc^2 = [((Gamma+1)p + Pi)/rho - dPi/drho]/Gammaagainst a numerically integrated isentrope to 1e-6;stiffened gas as the degenerate member;
C^1continuity of the reference curve atrho0.bit-identical.
Deliberately not here
Solver wiring - per-cell evaluation in the mixture coefficients, the derivative terms in
f_bulk_modulus, Wood's law for N fluids, frozen mixing for Mie-Gruneisen under the five-equationmodel - and the two validation problems, Hugoniot recovery and isentropic release as
convergencecases.Those follow in a second PR once this backend is reviewed. The reference curve has a pole at
rho/rho0 = s/(s - 1); the validator warning for it lands with the cases, since it needs the initial state.