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

[fahnab666]

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:

1. MFC gets multi-EOS support in a quick and easy way that is ensured correct and extensible.
2. 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

1. Central thermodynamics interface, no changed answers.
2. Explicit EOS selector plus existing stiffened-gas and Pyrometheus adapters.
3. 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.