Add JWL equation of state and detonation modeling for the five-equation model

Summary

The goal is Jones-Wilkins-Lee (JWL) detonation products in the five-equation Allaire model (model_eqns = 2), selectable per fluid (fluid_pp(i)%eos) and coexisting with an ideal-gas or stiffened-gas ambient, so MFC can simulate explosive products expanding into air or water.

A monolithic implementation was closed in favor of an incremental equation-of-state architecture. This issue tracks that incremental path: a central thermodynamics interface first, then explicit EOS selection, then a generic Mie-Gruneisen backend, and only then JWL as a particular Mie-Gruneisen reference curve with user-supplied parameters.

The two goals of this path are:

MFC gets multi-EOS support in a quick and easy way that is ensured correct and extensible.
MFC carries no lingering export-control risk, because no calibrated explosive catalog ever enters the repository.

Current PR

Add a per-fluid equation-of-state selector #1700 adds fluid_pp(i)%eos, a stable enumeration for a fluid's equation of state, in place of implicit backend selection. Only the two already-supported backends are accepted (stiffened_gas default, ideal_gas_mixture for chemistry builds via Pyrometheus); mie_gruneisen, jwl, and table are reserved in the enum and rejected until their backends exist. The default resolves from the compile-time chemistry flag, so existing cases do not change and the golden suite stays bit-identical.

Plan

The right route is incremental EOS support first, while keeping thermodynamics, reaction kinetics, and material data separate. Pyrometheus remains the gas-chemistry backend; it is not forced to represent every phenomenological explosive burn law.

MFC currently has stiffened-gas thermodynamics for ordinary fluids and a Pyrometheus-generated ideal-gas-mixture path when chemistry = T. What it lacks is an explicit, coherent EOS/backend abstraction.

PR 1: Establish the thermodynamics interface. No new physics, no user-facing EOS choice. Centralize the operations currently scattered through conversion routines, Riemann solvers, post-processing, and coupling code: pressure(state), internal_energy(rho, pressure, composition), temperature(state), sound_speed(state), thermodynamic_derivatives(state), validate_state(state). Implement two adapters, stiffened_gas_thermo and pyrometheus_thermo, the latter wrapping the already-generated m_thermochem and selected automatically when chemistry is enabled. Results stay bit-identical. The goal is to prove all solver paths obtain thermodynamics through one interface.

PR 2: Introduce explicit EOS selection. Add a stable enumeration rather than scattered booleans: EOS_STIFFENED_GAS, EOS_IDEAL_GAS_MIXTURE, EOS_MIE_GRUNEISEN, EOS_JWL, EOS_TABLE. Initially expose only the already-supported choices (for example fluid_pp(1)%eos: "stiffened_gas"); when chemistry = T the reacting gas uses ideal_gas_mixture backed by Pyrometheus. Preserve stiffened gas as the default, reject unsupported combinations explicitly, do not promise arbitrary combinations of EOS families within one cell, and keep the integer representation internal while preferring readable names in case files.

PR 3: Add generic Mie-Gruneisen support. A cleaner second non-chemistry EOS than jumping straight to an explosives-branded path. Provide reference pressure p_ref(rho), reference energy e_ref(rho), the Gruneisen coefficient, pressure/energy inversion, analytic sound speed and derivatives, and manufactured tests. JWL then becomes a particular Mie-Gruneisen reference curve rather than an unrelated special path.

PR 4: Add nonreactive JWL. Add only the JWL reference curve, its derivatives and inversion, user-supplied parameters, and generic products-expansion and shock-tube cases. No named explosive catalog, no reaction source. The repository holds the equation family, not a library of military-relevant materials.

PR 5: Create the external material-file mechanism. Follow the existing Cantera-file model: a material file naming an eos family and its parameters plus a provenance block (citation, release status). Search locations mirror Cantera mechanisms (explicit path, case directory, configured public material directory). Analytic EOS parameters stay runtime values; tabular EOS data stays external; only a custom symbolic EOS would need generated source.

PR 6: Separate chemistry from reaction integration. Formalize three concepts: thermodynamic backend, reaction mechanism, and reaction integrator. For current Pyrometheus chemistry these are the ideal-gas mixture, a mechanism generated from a Cantera file, and an explicit or alpha-QSS integrator. The interface is conceptually reaction_rates(state, mechanism_parameters), advance_reactions(state, dt), reaction_timestep_constraint(state). This lets integration improve without changing EOS selection and lets an EOS be added without implicitly adding reactions.

PR 7: Strengthen the Pyrometheus path. Before coupling it to new EOS families, make it a first-class supported backend: mechanism identity in build hashes, generated-code provenance recorded, mechanism validation before compilation, elemental and mass-conservation diagnostics, constant-volume and constant-pressure reactor tests, CPU/GPU agreement tests, and explicit reporting of the selected thermodynamic and kinetic backends. Keep generated m_thermochem.f90 in the staging/build area, not the source tree.

PR 8: Add generic progress-variable reactions. Distinct from Pyrometheus elementary chemistry. Introduce a small generic mechanism for condensed materials: progress lambda, a supplied rate law d(lambda)/dt, reactant EOS, product EOS, and energy release. Start with a published generic rate law or manufactured source. Test the zero-rate limit, homogeneous-reactor convergence, energy conservation, positivity and bounds 0 <= lambda <= 1, and equilibrium/end-state recovery.

PR 9: Couple reactants and products. Only after the individual EOS implementations work: reactant EOS, product JWL EOS, reaction progress, consistent mixture energy, pressure equilibrium, reaction energy offset, and sound speed/derivatives for intermediate states. This is where the hard thermodynamics lives; it should not be hidden inside "add a reaction source."

PR 10: Add specific reaction families individually. Each family in a separate PR: generic programmed burn, pressure-dependent progress law, ignition-and-growth-style model, and gas-phase afterburn coupled to Pyrometheus where scientifically appropriate. Each PR initially uses synthetic or clearly unrestricted parameters; real material calibration remains external.

Recommended first three PRs

Central thermodynamics interface, no changed answers.
Explicit EOS selector plus existing stiffened-gas and Pyrometheus adapters.
Generic Mie-Gruneisen, followed by JWL as a reference-curve specialization.

That gives MFC a durable foundation. External material files and generic progress-variable reactions can then be added without building explosive-specific assumptions into the solver core.

Where Pyrometheus fits, and where it does not

Pyrometheus is well suited to gas-phase species chemistry, Cantera mechanisms, thermochemical properties, afterburn involving conventional chemical kinetics, and generated portable CPU/GPU reaction kernels. It is not the natural sole representation for program burn, ignition-and-growth progress laws, pressure-triggered empirical condensed-phase reactions, or a reactant-to-JWL-products mixture closure. Those should use the same reaction-integration infrastructure with a different mechanism backend; encoding them all as Cantera chemistry would blur important physics.

Target capability

JWL pressure law for pure products:

$$p = A\left(1-\frac{\omega}{R_1V}\right)e^{-R_1V} + B\left(1-\frac{\omega}{R_2V}\right)e^{-R_2V} + \omega\rho e, \qquad V=\frac{\rho_0}{\rho}.$$

A composition-weighted mixture closure between JWL products and the ambient, with a closed-form pressure-energy inverse and analytic sound speed, degenerating exactly to the pure laws at the mixture endpoints.
Optional reaction sources, off by default: kinematic program burn, products-air afterburn, and JWL++ (Souers 2000) reactive burn.
Coupling with immersed boundaries and Euler-Lagrange bubbles.
Validation: exact Riemann star states, resolved ZND structure against the analytic Hugoniot/Rayleigh/CJ construction, Sedov-Taylor blast scaling, far-field TNT overpressure against Kinney-Graham.

References

Arienti, M., Morano, E., Shepherd, J.E., Shock and detonation modeling with the Mie-Gruneisen equation of state, GALCIT Report FM99-8 (2004).
Souers, P.C. et al., JWL++: A Simple Reactive Flow Code Package for Detonation, Propellants, Explosives, Pyrotechnics 25 (2000).
Garno, J.M. et al., Phys. Rev. Fluids 5, 123201 (2020), reactant energy offset for resolved ZND structure.

Two constraints that surfaced while centralising the EOS expressions (#1762), both of which bound what a JWL backend can cover more than any interface question does.

1. The pressure relaxation needs a closed-form invertible isentrope, which JWL lacks.

s_equilibrate_pressure (m_pressure_relaxation.fpp) already runs a per-cell Newton solve inside a [seq] device routine — so iterative equilibration on device is an existing, already-paid cost, not a new one. But the iteration needs rho(p) and drho/dp per phase in closed form:

rho_K_s(i) = ... * ((pres_relax + ps_inf(i))/(pres_K_init(i) + ps_inf(i)))**(1._wp/gs_min(i))

That is the stiffened-gas isentrope (Saurel et al., JFM 2009). JWL's pressure inversion p(rho, e) is closed-form and ports directly; its isentrope is not analytically invertible. A JWL phase in the six-equation model therefore needs either a nested sub-iteration inside that Newton loop or a reformulated relaxation. Numerical-methods work, not refactoring.

2. Phase change needs a thermodynamically complete EOS; JWL is mechanical.

m_phase_change computes entropy, temperature and per-phase enthalpy (sk, hk, rhok, ek), and is the only consumer of cvs and qvps. Stiffened gas has closed forms for all of them. JWL as normally written is p(rho, e) with no temperature or entropy unless a caloric extension is supplied. This is a physics limit rather than an implementation gap — the likely answer is that JWL with phase change is refused at validation.

Neither blocks a JWL backend for the mechanical relations (energy, pressure, sound speed), which are all closed-form and now go through single operators after #1762. They bound where it can be used.

Also relevant: a per-fluid eos_model dispatch is viable on GPU because the branch is wavefront-uniform — every lane sits at the same fluid index — unlike branching on cell contents. But it must not be implemented with a derived type in the operator signature; see #1714, which cost 20% on gfx90a that way.

Working through what the generic Mie-Gruneisen backend actually requires of the existing operators,
with numerical verification. Short version: the mixture rule needs no change, and the sound speed
needs exactly one extra term.

Mie-Gruneisen maps exactly onto the coefficients MFC already carries

Any Mie-Gruneisen EOS is affine in p at fixed density, so it can be written in MFC's stored form
with density-dependent coefficients:

rho*e = Gamma(rho)*p + Pi(rho)

For JWL in standard form,

p = A(1 - w/(R1 V))exp(-R1 V) + B(1 - w/(R2 V))exp(-R2 V) + w*rho*e

this gives Gamma = 1/w (constant) and Pi(rho) = -g(rho)/w, where g collects the two
exponential terms. Consistency check: an ideal gas has w = gamma - 1, so Gamma = 1/(gamma - 1),
which is MFC's stored form for fluid_pp(i)%gamma.

The five-equation mixture rule survives mixed EOS families unchanged

s_compute_mixture_coefficients computes Gamma_mix = sum(alpha_k Gamma_k) and
Pi_mix = sum(alpha_k Pi_k). That is just pressure equilibrium plus additivity of internal energy,
solved for p. What it requires is that each phase's (rho e)_k be affine in p - not that the
phases share an EOS family. Since Mie-Gruneisen is affine in p, a JWL phase mixes with a stiffened
or ideal phase with no change to the mixture rule at all. JWL products expanding into air or water in
model_eqns = 2 is therefore sound, which is what this issue assumes.

The sound speed needs one extra term, and `c = sqrt(K/rho)` still holds

From c^2 = (dp/drho)_e + (p/rho^2)(dp/de)_rho with rho e = Gamma(rho) p + Pi(rho):

c^2 = [ (Gamma + 1) p + Pi - rho*Pi' ] / (rho*Gamma)  -  p*Gamma'/Gamma

stiffened / ideal gas: Pi' = Gamma' = 0, so this reduces to ((Gamma+1)p + Pi)/(rho*Gamma),
which is exactly today's f_bulk_modulus/rho. Existing cases stay bit-identical by construction.
JWL: Gamma' = 0, so the single new term is - rho*Pi'.
general Mie-Gruneisen: both terms.

So f_bulk_modulus generalizes by one term that is identically zero for the current backends, and
every consumer of it - every sound speed in MFC is K/rho, differing only in how phases are mixed -
is untouched.

Mixture sound speed

Wood's law follows from differentiating volume additivity 1/rho = sum(Y_k/rho_k) at frozen mass
fractions with drho_k = dp/c_k^2, giving 1/(rho c^2) = sum(alpha_k/K_k) with no reference to the
EOS form. It therefore inherits correctness from the per-phase K_k above, and the alt_soundspeed
path works for mixed EOS as soon as f_bulk_modulus carries the new term.

Worth flagging: the default (non-alt_soundspeed) path instead uses Gamma_mix/Pi_mix
directly - the frozen mixture speed. With a density-dependent Pi those two differ, so mixed-EOS
cases should use the Wood path, or the default path needs its own treatment.

Verification

All of the above checked numerically against finite differences of the defining relations:

JWL, c^2 formula vs (dp/drho)_e + (p/rho^2)(dp/de)_rho, 12 states spanning compression and
expansion: agreement to ~1e-10 relative.
A synthetic EOS with both Gamma(rho) and Pi(rho) varying, to exercise the Gamma' term:
agreement to ~1e-10, including a state with c^2 < 0, so the formula is not accidentally correct
only in the physically stable region.
Wood's law with the generalized per-phase K on a genuine stiffened-water + JWL mixture, over
two decades of pressure with alpha_1 from 0.29 to 0.60: agreement to ~1e-9.

Implication for the incremental path

A backend supplies Gamma(rho), Pi(rho) and their density derivatives. Nothing else in the solver
learns that a second EOS exists: no per-operator dispatch, no widened signatures in the hot loops,
and - importantly given #1714 - no derived type anywhere near a [seq] device routine.

#1762 lands PRs 1 and 2 of the plan above, and supersedes #1700 and #1714. A note on what the next
step actually needs, since the shape of PR 3 changed once the operators were written.

What is in place

Every EOS expression now has one definition, in src/common/m_variables_conversion.fpp:

s_compute_mixture_coefficients(alpha_rho, alpha, rho, gamma, pi_inf, qv)
s_compute_energy(pres, alpha_rho, alpha, vel_sum, E)
s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c)      ! _avg for interface averages
f_bulk_modulus(pres, gamma, pi_inf)
f_pressure(e_int, gamma, pi_inf, qv)
f_phase_internal_energy(pres, alpha, alpha_rho, gamma, pi_inf, qv)
f_isentrope_exponent(gamma) / f_isentrope_pressure(pi_inf, gamma) / f_pressure_on_isentrope(...)

Copies collapsed: mixture rule 9, bulk modulus 7, pressure inversion 5, elastic strain energy 4,
per-phase internal energy 3, isentrope parameters 3. fluid_pp(i)%eos is live with
stiffened_gas and ideal_gas, and ideal gas is a real second backend rather than a relabelling -
it does no arithmetic on pi_inf at all. A Python/Fortran enum-agreement test guards the two
tables, because compound registry keys cannot be auto-generated.

PR 3 is smaller than it looks

MFC's stored form is already Mie-Gruneisen. Writing MFC's gamma as $\Gamma$ and pi_inf as $\Pi$,
the solver only ever needs

$$\rho e = \Gamma p + \Pi$$

and any Mie-Gruneisen EOS $p = p_\text{ref}(\rho) + \Gamma_G(\rho)\rho,(e - e_\text{ref}(\rho))$
maps onto it exactly with

$$\Gamma(\rho) = \frac{1}{\Gamma_G(\rho)}, \qquad \Pi(\rho) = \rho,e_\text{ref}(\rho) - \frac{p_\text{ref}(\rho)}{\Gamma_G(\rho)}.$$

For JWL specifically $\Gamma = 1/\omega$ is constant and $\Pi(\rho) = -g(\rho)/\omega$, where $g$ is
the density-only part of the JWL pressure. Checked numerically against the JWL law across a range
of densities and energies; the identity is exact.

So the generalisation is not a new state representation, it is letting $\Gamma$ and $\Pi$ depend on
$\rho$. Two consequences:

The sound speed needs two extra terms. The general form is

$$c^2 = \frac{1}{\Gamma}\left[\frac{(\Gamma+1)p + \Pi}{\rho} - \Pi'(\rho) - p,\Gamma'(\rho)\right]$$

which collapses to the existing ((gamma + 1)*pres + pi_inf)/(gamma*rho) whenever
$\Gamma' = \Pi' = 0$, so stiffened gas and ideal gas stay bit-identical. Verified against a
numerical $(\partial p/\partial\rho)_s$ for JWL to a relative error of 3e-11, and against a
synthetic EOS with both derivatives nonzero.

Mixtures still work. Wood's law is a harmonic mean of bulk moduli, and f_bulk_modulus is
already the single definition of that quantity, so a JWL-in-air or JWL-in-water mixture needs no new
mixture rule - only a per-fluid $K$ that knows about $\Pi'$ and $\Gamma'$. Checked on a
stiffened-gas + JWL mixture.

Concrete shape of the patch

eos_mie_gruneisen in m_constants.fpp and _EOS_NAMES in params/definitions.py (both, or
the enum test fails).
Per-fluid reference-curve parameters on fluid_pp, plus a case_validator rule rejecting a
Mie-Gruneisen fluid where model_eqns or a feature cannot support it.
f_mie_gruneisen_coefficients(rho, i, gamma, pi_inf, dgamma, dpi_inf) returning the four
quantities the sound speed needs.
s_compute_speed_of_sound gains the two derivative terms behind a branch that is dead for the
constant-coefficient backends.
JWL then supplies one reference curve, which is PR 4 and small.

The one structural change: gammas(i)/pi_infs(i) are currently filled once at start-up from
fluid_pp, because a per-cell branch there costs registers in the Riemann kernels. A
state-dependent EOS has to be evaluated per cell instead, so the coefficient formation moves into
s_compute_mixture_coefficients for those fluids only. That is the part worth benchmarking early -
it is in the hot path, and #1714 was abandoned for a 20% regression from exactly this kind of
change.

Two constraints worth knowing before writing device code

Both cost a CI round trip on #1762 and neither is visible from a local build:

CCE OpenACC rejects an array whose bound is a device global (dimension(num_fluids)) passed to a
device routine from inside another GPU_ROUTINE(parallelism='[seq]'); the same call is fine
from a parallel-loop body. It reports ftn-7066 ... Global in accelerator routine without declare
at whatever line it gave up on, and the line walks forward as triggers are removed, so the
reported location is not the cause. Only non-case-optimized lanes fail.
nvfortran 23.11/24.1 segfault in fort2 on a caller that passes a parameter array from
m_thermochem into a declare-target routine. Read those directly in the kernel instead.

Both are now in .claude/rules/common-pitfalls.md.

One more site for the Mie-Gruneisen backend, found while auditing #1762 for inline stiffened-gas
arithmetic and deliberately left alone there.

f_phase_internal_energy(pres, alpha, alpha_rho, gamma, pi_inf, qv) gives
e_k = alpha_k (Gamma_k p + Pi_k) + alpha_rho_k qv_k. It has no inverse in the module, and the one place
that needs it inverts it by hand - src/simulation/m_pressure_relaxation.fpp:139:

pres_K_init(i) = ((q_cons_vf(i + eqn_idx%int_en%beg - 1)%sf(j, k, l) &
                 & - q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)) &
                 & /q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l) - pi_infs(i))/gammas(i)

which is exactly p_k = ((e_k - alpha_rho_k qv_k)/alpha_k - Pi_k)/Gamma_k.

Left as-is in #1762 because it is a single call site, and an f_phase_pressure there would be an
abstraction for one caller. A second EOS creates the second caller: the phasic inversion is EOS-specific,
so this line has to change for JWL, and at that point the forward/inverse pair earns its place beside
f_phase_internal_energy.

Related: the mixture-level inversion is already centralized as f_pressure; this is the per-phase one.
The isentrope inversion that pressure relaxation also needs is tracked separately in #1785.

PR 3 of the plan is split in two. #1805 lands the backend: Mie-Gruneisen as one dispatch over a reference curve, mapped onto MFC's rho e = Gamma p + Pi form so the operators centralized in #1762 are reused unchanged, with stiffened gas as the degenerate member and no solver wiring yet - every existing answer bit-identical. The follow-up wires it into the mixture coefficients and sound speed (frozen per-phase mixing under the five-equation model, Wood's law generalized to N fluids) and adds Hugoniot-recovery and isentropic-release validation cases via the convergence harness. JWL then becomes one more case supplying its reference curve.

Navigation Menu

Repository navigation

Add JWL equation of state and detonation modeling for the five-equation model #1638

Summary

Current PR

Plan

Recommended first three PRs

Where Pyrometheus fits, and where it does not

Target capability

References

Mie-Gruneisen maps exactly onto the coefficients MFC already carries

The five-equation mixture rule survives mixed EOS families unchanged

The sound speed needs one extra term, and `c = sqrt(K/rho)` still holds

Mixture sound speed

Verification

Implication for the incremental path

What is in place

PR 3 is smaller than it looks

Concrete shape of the patch

Two constraints worth knowing before writing device code

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Add JWL equation of state and detonation modeling for the five-equation model #1638

Description

Summary

Current PR

Plan

Recommended first three PRs

Where Pyrometheus fits, and where it does not

Target capability

References

Activity

sbryngelson commented on Aug 26, 2026

sbryngelson commented on Aug 28, 2026

Mie-Gruneisen maps exactly onto the coefficients MFC already carries

The five-equation mixture rule survives mixed EOS families unchanged

The sound speed needs one extra term, and c = sqrt(K/rho) still holds

Mixture sound speed

Verification

Implication for the incremental path

sbryngelson commented on Aug 29, 2026

What is in place

PR 3 is smaller than it looks

Concrete shape of the patch

Two constraints worth knowing before writing device code

sbryngelson commented on Aug 30, 2026

sbryngelson commented on Sep 2, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

The sound speed needs one extra term, and `c = sqrt(K/rho)` still holds