jwl: stiffened-gas ambient, reaction sources, and detonation validation

Extend the Rocflu JWL closure and add reactive-burn capability, all gated
behind default-off flags so non-JWL cases are byte-for-byte unchanged.

EOS closure:
- Support a stiffened-gas ambient (e.g. water) alongside ideal-gas air by
  carrying a cold-stiffness offset that is independent of density and energy,
  preserving the exact pressure<->energy inverse and the Grueneisen sound
  speed. Blend in mass fraction for a stiffened ambient to keep c^2 convex.
- Fuse the Riemann-path sound speed into a single closure evaluation
  (inverse + forward state) to roughly halve the per-face JWL EOS cost.
- Correct the sound-speed floor labelling, widen the start-up self-scan to
  reflected-shock states, and make the temperature floor unit-agnostic.

Reaction energy sources (m_jwl_sources):
- Program burn: kinematic lighting-time front from a detonation point.
- Afterburn: advected progress variable with mixing-rate or Arrhenius law,
  budget-capped by the afterburn energy.
- JWL++ reactive burn: pressure-driven, self-propagating detonation.

Diagnostics and guards:
- JWL-aware sound speed in post_process c_wrt output.
- Pass the cell products mass fraction to probe pressure output.
- Reject configurations that assume stiffened-gas mixture relations for JWL
  cells (pressure-based wave speeds, characteristic BCs, alt_soundspeed,
  elasticity) in both the Fortran checkers and the case validator.

Validation:
- Registered golden tests for the stiffened closure and all reaction models.
- Underwater products/water shock tube matching the exact Riemann star state.
- Free-air TNT blast checked against the Kinney-Graham overpressure curve.
- Documented the theory and citations in m_jwl and README-JWL-EOS.

Author: Fahad Nabid

README-JWL-EOS.md

@@ -16,15 +16,27 @@ The cold-curve pressure is evaluated through a shared helper so pressure recover
 
 ## Mixture Closures
 
-The JWL mixture closure applies to five-equation JWL/ideal-gas mixtures with one JWL products fluid and one non-JWL ideal-gas fluid. Whenever a JWL fluid is present, MFC applies the Rocflu single-fluid closure from `modflu/RFLU_ModJWL.F90`. It ramps the JWL `A` and `B` coefficients linearly in specific internal energy between the ambient-gas energy `jwl_air_e0` and the products reference energy `e_j = jwl_E0/jwl_ej_rho_ref`, ramps `omega` with mixture density via a smoothstep between `jwl_air_rho0` and `jwl_rho0`, and blends the heat capacity as the mass-weighted `cv = Y*cv_products + (1 - Y)*cv_air`. The transition from these mixture coefficients to pure-products JWL coefficients is a smoothstep in products mass fraction `Y` over `[0.95, 0.999]`, replacing the earlier hard cutoff at `Y = 0.99`; because the smoothstep weight depends only on `Y`, the analytic pressure-to-energy inverse remains exact under the blend.
+The JWL mixture closure applies to five-equation mixtures with one JWL products fluid and one non-JWL ambient fluid, which may be an ideal gas (air) or a stiffened gas (e.g. water; see "Stiffened-Gas Ambient" below). Whenever a JWL fluid is present, MFC applies the Rocflu single-fluid closure from `modflu/RFLU_ModJWL.F90`. It ramps the JWL `A` and `B` coefficients linearly in specific internal energy between the ambient-gas energy `jwl_air_e0` and the products reference energy `e_j = jwl_E0/jwl_ej_rho_ref`, ramps `omega` with mixture density via a smoothstep between `jwl_air_rho0` and `jwl_rho0`, and blends the heat capacity as the mass-weighted `cv = Y*cv_products + (1 - Y)*cv_air`. The transition from these mixture coefficients to pure-products JWL coefficients is a smoothstep in products mass fraction `Y` over `[0.95, 0.999]`, replacing the earlier hard cutoff at `Y = 0.99`; because the smoothstep weight depends only on `Y`, the analytic pressure-to-energy inverse remains exact under the blend.
 
 This closure is not TNT-specific. MFC reads `jwl_A`, `jwl_B`, `jwl_R1`, `jwl_R2`, `jwl_omega`, `jwl_rho0`, `cv`, and either `jwl_Q` or `jwl_E0` from the JWL fluid, plus `jwl_air_rho0` and either `jwl_air_e0` or `jwl_air_p0` for the co-existing gas, so any explosive products model with a valid JWL parameter set can be used. `jwl_Q` is the specific detonation energy in J/kg; internally MFC derives `jwl_E0 = jwl_rho0*jwl_Q`. If both are provided, they must be consistent. The ambient-gas Grüneisen coefficient is taken from the ideal-gas fluid's own `gamma` (`Gamma_air = 1/gamma`; with a single JWL fluid the JWL fluid's own `gamma` is used), and the optional `jwl_ej_rho_ref` sets the products-energy reference density (default `jwl_rho0`). Unlike Rocflupicl's case-specific implementation, MFC does not hard-code TNT density, TNT energy scaling, ambient density, ambient energy, or air heat capacity inside the EOS. The reference parameters must satisfy `jwl_rho0 > jwl_air_rho0` and `jwl_E0/jwl_ej_rho_ref > jwl_air_e0`.
 
-Finite pressure, temperature, and energy floors are applied only after explicit finite checks. The state routine returns the raw squared sound speed; the public wrappers bound it below by the ideal-gas value `Gamma_air*p/rho`. NaNs are otherwise intentionally preserved so bad states are visible during debugging instead of being converted into plausible-looking floor values.
+## Stiffened-Gas Ambient (Underwater/Condensed)
+
+When the non-JWL fluid has `pi_inf > 0`, MFC treats it as a stiffened gas and switches the blend variable: `omega`, the `A`/`B` energy ramp, and the cold stiffness `pi_c = (Gamma + 1)*pi_hat` all blend in products mass fraction `Y` (smoothstep over [0, 1]) instead of density, because a density ramp of a stiff liquid's large Grueneisen coefficient makes the mixture sound speed non-convex. The pure-ambient limit then recovers `p = Gamma*rho*e - (Gamma + 1)*pi_inf` exactly at any density, so pure-water shocks follow the exact stiffened-gas law. The true stiffness is derived from the ambient fluid's own inputs (`pi_inf_true = fluid_pp%pi_inf/(fluid_pp%gamma + 1)`); no new case parameters are needed, and `jwl_air_p0` resolves the ambient energy through the stiffened form `e0 = (p0*gamma + pi_inf)/rho0` (MFC-convention `gamma`/`pi_inf`). Because `pi_c` is independent of both `e` and `rho`, the analytic pressure-to-energy inverse remains exact and the Grueneisen sound-speed identity is unchanged. The pressure floor doubles as the standard cavitation (pressure-cutoff) model, and the sound-speed safety floor becomes `min(Gamma, omega)*(p + pi_hat)/rho`. The ideal-gas path (`pi_inf = 0`) is bit-identical to the closure described above. Note the closure's temperature remains `T = e/cv` in the pure-ambient limit, which for a stiffened liquid counts the cold (stiffness) energy as thermal; `T` is a diagnostic on this path (pressure and sound speed are `cv`-free) and is not a physical liquid temperature. A 1D products/water shock tube (`examples/1D_jwl_underwater_shocktube`) reproduces the exact two-material Riemann star state (JWL isentrope vs. water shock Hugoniot) to under 1%.
+
+Finite pressure, temperature, and energy floors are applied only after explicit finite checks. The state routine returns the raw squared sound speed; the public wrappers bound it below by `Gamma_air*p/rho`, a safety floor that sits below any physical mixture sound speed (pure air gives `(Gamma_air + 1)*p/rho`, pure products `(omega + 1)*p/rho` plus cold-curve terms) so it engages only on unphysical states. NaNs are otherwise intentionally preserved so bad states are visible during debugging instead of being converted into plausible-looking floor values.
+
+## Reaction Sources
+
+Three optional detonation-energy sources sit on top of the closure (all default off, so a plain products-expansion run is unaffected):
+
+- **Program burn** (`prog_burn`): a Rocflu-style kinematic front expands from the detonation point (`pb_x_det`/`pb_y_det`/`pb_z_det`, time `pb_t_det`) at speed `pb_D_cj`, depositing the detonation energy `jwl_Q` over a reaction zone of width `pb_width`. No extra field is transported; each swept cell receives exactly `jwl_Q` per unit explosive mass.
+- **Afterburn** (`jwl_afterburn`): products-air combustion energy `jwl_q_ab` released through an advected progress variable, either mixing-rate (`jwl_ab_model = 1`, time scale `jwl_ab_tau`) or Arrhenius (`jwl_ab_model = 2`, default; `jwl_ab_A`, `jwl_ab_theta`, `jwl_ab_n`). Requires an ideal-gas ambient. Since `jwl_q_ab` adds to `jwl_Q`, use a detonation-only JWL fit to avoid double counting.
+- **Reactive burn** (`jwl_reactive`, JWL++ after Souers 2000): a self-propagating detonation driven by the local pressure, `dλ/dt = jwl_G·p^jwl_b_exp·(1−λ)`, releasing `jwl_Q` as the explosive reacts. Mutually exclusive with `prog_burn`; `jwl_G`/`jwl_b_exp` are calibrated per explosive to reproduce its CJ velocity.
 
 ## Validation Scope
 
-The production closure is covered by a registered golden test. The regression includes a homogeneous 50/50 products-air slab so the Rocflu density and energy interpolation is exercised in addition to its endpoint branches.
+The production closure is covered by registered golden tests: the Rocflu products-air closure (with a homogeneous 50/50 slab exercising the density/energy interpolation), a stiffened-gas (water) ambient, the program-burn plus afterburn sources, and the JWL++ reactive burn. Beyond the golden regression, verified physics includes: the underwater products/water shock tube against the exact two-material Riemann star state (<1%); the afterburn and JWL++ energy budgets closing to their `Y·q` bounds (<1%); the program-burn front speed matching `pb_D_cj`; and the 3D free-air TNT blast tracking the Kinney-Graham overpressure curve in the far field (`benchmarks/3D_jwl_spherical_tnt_free_air_validation/validate_kingery.py`).
 
 At initialization MFC scans the closure over the configured material's `(rho, e, Y)` envelope and aborts if any state yields a non-positive or non-finite squared sound speed, or if the pressure-to-energy round-trip fails to recover the input energy. An inconsistent parameter set therefore fails fast at startup rather than producing silently wrong states during the run.
 

benchmarks/3D_jwl_spherical_tnt_free_air_validation/validate_kingery.py

@@ -0,0 +1,78 @@
+#!/usr/bin/env python3
+"""Physics validation of free-air JWL/TNT blast against a reference air-blast curve.
+
+Reads the reduced gauge metrics (gauge_results.csv: peak overpressure vs scaled
+distance Z = R / W^(1/3)) and compares MFC's peak incident overpressure to the
+Kinney & Graham (1985) closed-form correlation for a spherical TNT charge in
+free air. This checks the physics (are the JWL coefficients and closure giving a
+correct blast), not just implementation self-consistency.
+
+Kinney-Graham peak overpressure ratio (Z in m/kg^(1/3), valid ~0.05 <= Z <= 40):
+
+    Po/Pa = 808 [1 + (Z/4.5)^2]
+            / sqrt(1 + (Z/0.048)^2) / sqrt(1 + (Z/0.32)^2) / sqrt(1 + (Z/1.35)^2)
+"""
+
+import argparse
+import csv
+import math
+from pathlib import Path
+
+P_AMB = 101325.0
+
+
+def kinney_graham_overpressure(z):
+    """Peak incident overpressure (Pa) at scaled distance z (m/kg^(1/3))."""
+    num = 808.0 * (1.0 + (z / 4.5) ** 2)
+    den = math.sqrt(1.0 + (z / 0.048) ** 2) * math.sqrt(1.0 + (z / 0.32) ** 2) * math.sqrt(1.0 + (z / 1.35) ** 2)
+    return P_AMB * num / den
+
+
+def main():
+    script_dir = Path(__file__).resolve().parent
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--csv", type=Path, default=script_dir / "gauge_results.csv")
+    parser.add_argument(
+        "--tol",
+        type=float,
+        default=0.20,
+        help="max acceptable relative error in the reliable far-field band",
+    )
+    parser.add_argument(
+        "--z-farfield",
+        type=float,
+        default=0.30,
+        help="scaled distance above which Kinney-Graham and the grid resolution are both reliable",
+    )
+    args = parser.parse_args()
+
+    with args.csv.open() as stream:
+        rows = [r for r in csv.DictReader(stream) if r.get("status") == "ok"]
+    if not rows:
+        raise SystemExit(f"no usable gauge rows in {args.csv}")
+
+    # The near-field (small Z) is expected to deviate: the charge is burst-initialized
+    # rather than detonated, the steep gradient is under-resolved on a coarse grid, and
+    # Kinney-Graham itself has larger uncertainty below Z~0.2. The physically meaningful
+    # test is the far-field band, together with the convergence trend toward the curve.
+    print(f"{'Z (m/kg^1/3)':>13} {'MFC (kPa)':>12} {'K-G (kPa)':>12} {'rel err':>9}  band")
+    worst_far = 0.0
+    for r in rows:
+        z = float(r["scaled_distance"])
+        mfc = float(r["peak_overpressure_pa"])
+        ref = kinney_graham_overpressure(z)
+        rel = abs(mfc - ref) / ref
+        band = "far" if z >= args.z_farfield else "near"
+        if band == "far":
+            worst_far = max(worst_far, rel)
+        print(f"{z:13.4f} {mfc / 1e3:12.1f} {ref / 1e3:12.1f} {rel:9.2%}  {band}")
+
+    print(f"\nworst far-field (Z >= {args.z_farfield}) error: {worst_far:.2%} (tolerance {args.tol:.0%})")
+    if worst_far > args.tol:
+        raise SystemExit(f"FAIL: far-field overpressure deviates from Kinney-Graham by {worst_far:.1%}")
+    print("PASS: far-field blast overpressures track the Kinney-Graham TNT curve;")
+    print("      near-field converges toward it as Z increases (expected for burst init + coarse grid)")
+
+
+if __name__ == "__main__":
+    main()

docs/documentation/case.md

@@ -443,7 +443,9 @@ Details of implementation of viscosity in MFC can be found in \cite Coralic15.
 
 - `fluid_pp(i)%%G` is required for `hypoelasticity`.
 
-- `fluid_pp(i)%%eos` selects the equation of state for the $i$-th fluid: [1] stiffened gas (default); [2] Jones-Wilkins-Lee (JWL) for detonation products, supported with `model_eqns = 2` and at most one JWL fluid. A JWL fluid requires `jwl_A`, `jwl_B`, `jwl_R1`, `jwl_R2`, `jwl_omega`, `jwl_rho0`, and either `jwl_Q` (specific detonation energy, J/kg) or `jwl_E0` (volumetric detonation energy, J/m³; MFC derives `jwl_E0 = jwl_rho0*jwl_Q`), plus `jwl_air_rho0` and either `jwl_air_e0` (specific internal energy) or `jwl_air_p0` (pressure) describing the co-existing ideal gas. The ideal gas Grüneisen coefficient \f$\Gamma_{\mathrm{air}} = 1/\gamma\f$ is taken from that gas fluid's own `gamma`; with a single JWL fluid and no separate ideal-gas fluid, the JWL fluid's own `gamma` is used. The optional `jwl_ej_rho_ref` sets the products-energy reference density (default `jwl_rho0`, so \f$e_j = E_0/\rho_0\f$). Products are mixed with the surrounding gas by the Rocflu state-interpolated closure (Garno/Stanley formulation), which requires positive products and air `cv`.
+- `fluid_pp(i)%%eos` selects the equation of state for the $i$-th fluid: [1] stiffened gas (default); [2] Jones-Wilkins-Lee (JWL) for detonation products, supported with `model_eqns = 2` and at most one JWL fluid. A JWL fluid requires `jwl_A`, `jwl_B`, `jwl_R1`, `jwl_R2`, `jwl_omega`, `jwl_rho0`, and either `jwl_Q` (specific detonation energy, J/kg) or `jwl_E0` (volumetric detonation energy, J/m³; MFC derives `jwl_E0 = jwl_rho0*jwl_Q`), plus `jwl_air_rho0` and either `jwl_air_e0` (specific internal energy) or `jwl_air_p0` (pressure) describing the co-existing ideal gas. The ideal gas Grüneisen coefficient \f$\Gamma_{\mathrm{air}} = 1/\gamma\f$ is taken from that gas fluid's own `gamma`; with a single JWL fluid and no separate ideal-gas fluid, the JWL fluid's own `gamma` is used. The optional `jwl_ej_rho_ref` sets the products-energy reference density (default `jwl_rho0`, so \f$e_j = E_0/\rho_0\f$). Products are mixed with the surrounding gas by the Rocflu state-interpolated closure (Garno/Stanley formulation), which requires positive products and air `cv`. A stiffened-gas ambient (e.g. water) is supported by setting the non-JWL fluid's `pi_inf`; see README-JWL-EOS.md.
+
+- JWL reaction sources (simulation stage, require a JWL fluid): `prog_burn` enables a Rocflu-style kinematic program burn — a front expands from the detonation point (`pb_x_det`, `pb_y_det`, `pb_z_det`) starting at time `pb_t_det` with speed `pb_D_cj`, releasing the detonation energy `jwl_Q` over a reaction zone of width `pb_width`. `jwl_afterburn` enables products-air afterburn energy release with an advected progress variable (adds one equation): `jwl_ab_model` selects the rate law ([1] mixing-rate with time scale `jwl_ab_tau`; [2] Arrhenius, the default, with prefactor `jwl_ab_A`, activation temperature `jwl_ab_theta`, and pressure exponent `jwl_ab_n`), and `jwl_q_ab` is the specific afterburn energy (J/kg products, in addition to `jwl_Q` — use a detonation-only JWL fit to avoid double counting). Afterburn requires `riemann_solver = 2` and an ideal-gas ambient. `jwl_reactive` enables a JWL++ (Souers 2000) pressure-driven reactive burn (adds one equation): the reaction progress evolves as `dλ/dt = jwl_G * p^jwl_b_exp * (1 - λ)` and releases `jwl_Q` as the explosive reacts, so a detonation self-propagates from a high-pressure hot spot rather than a prescribed front (mutually exclusive with `prog_burn`; requires `riemann_solver = 2`).
 
 > **Stored-form parameters:** The values `gamma`, `pi_inf`, and `Re(1)`/`Re(2)` are **not** the raw physical quantities. MFC expects transformed stored forms:
 > - `gamma` = \f$1/(\gamma-1)\f$, not \f$\gamma\f$ itself

docs/module_categories.json

@@ -20,6 +20,7 @@
         "category": "Physics Models",
         "modules": [
             "m_jwl",
+            "m_jwl_sources",
             "m_viscous",
             "m_hb_function",
             "m_surface_tension",
@@ -100,4 +101,4 @@
             "m_boundary_conditions"
         ]
     }
-]
+]

examples/1D_jwl_underwater_shocktube/case.py

@@ -0,0 +1,101 @@
+#!/usr/bin/env python3
+import json
+
+# 1D underwater explosion shock tube exercising the stiffened-gas ambient branch of
+# the Rocflu JWL closure. Fluid 1 = TNT JWL products (eos=2); fluid 2 = water as a
+# stiffened gas (gamma = 6.12, pi_inf = 3.43e8 Pa). A constant-volume-burst
+# products slug (rho0, e = E0/rho0 -> ~9.3 GPa; 0 <= x <= 0.1 m) drives a
+# shock into water at atmospheric pressure.
+jwl_A = 3.712e11
+jwl_B = 3.231e9
+jwl_R1 = 4.15
+jwl_R2 = 0.95
+jwl_omega = 0.30
+jwl_rho0 = 1630.0
+jwl_E0 = 1.0089e10
+jwl_Cv = 613.5
+
+# Water stiffened gas: gamma = 6.12, pi_inf = 3.43e8 Pa. MFC convention:
+# fluid_pp%gamma = 1/(gamma-1), fluid_pp%pi_inf = gamma*pi_inf/(gamma-1).
+water_gamma = 6.12
+water_pi_inf = 3.43e8
+water_gamma_mfc = 1.0 / (water_gamma - 1.0)
+water_pi_inf_mfc = water_gamma * water_pi_inf / (water_gamma - 1.0)
+
+print(
+    json.dumps(
+        {
+            "run_time_info": "T",
+            "x_domain%beg": 0.0,
+            "x_domain%end": 1.0,
+            "m": 399,
+            "n": 0,
+            "p": 0,
+            "dt": 1.0e-7,
+            "t_step_start": 0,
+            "t_step_stop": 400,
+            "t_step_save": 100,
+            "num_patches": 2,
+            "model_eqns": 2,
+            "num_fluids": 2,
+            "mpp_lim": "T",
+            "mixture_err": "T",
+            "alt_soundspeed": "F",
+            "time_stepper": 3,
+            "weno_order": 3,
+            "weno_eps": 1.0e-16,
+            "mapped_weno": "T",
+            "null_weights": "F",
+            "mp_weno": "F",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -3,
+            "bc_x%end": -3,
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            # Patch 1: ambient water background over the whole domain.
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%length_x": 1.0,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%pres": 101325.0,
+            "patch_icpp(1)%alpha_rho(1)": 1.63e-5,
+            "patch_icpp(1)%alpha_rho(2)": 999.99999,
+            "patch_icpp(1)%alpha(1)": 1.0e-8,
+            "patch_icpp(1)%alpha(2)": 0.99999999,
+            # Patch 2: high-pressure JWL products slug (0 <= x <= 0.1 m).
+            "patch_icpp(2)%geometry": 1,
+            "patch_icpp(2)%alter_patch(1)": "T",
+            "patch_icpp(2)%x_centroid": 0.05,
+            "patch_icpp(2)%length_x": 0.1,
+            "patch_icpp(2)%vel(1)": 0.0,
+            "patch_icpp(2)%pres": 9.311e9,
+            "patch_icpp(2)%alpha_rho(1)": 1629.9999837,
+            "patch_icpp(2)%alpha_rho(2)": 1.0e-5,
+            "patch_icpp(2)%alpha(1)": 0.99999999,
+            "patch_icpp(2)%alpha(2)": 1.0e-8,
+            # Fluid 1: TNT JWL products (Rocflu closure, stiffened ambient branch).
+            "fluid_pp(1)%eos": 2,
+            "fluid_pp(1)%gamma": 2.5,
+            "fluid_pp(1)%pi_inf": 0.0,
+            "fluid_pp(1)%cv": jwl_Cv,
+            "fluid_pp(1)%jwl_A": jwl_A,
+            "fluid_pp(1)%jwl_B": jwl_B,
+            "fluid_pp(1)%jwl_R1": jwl_R1,
+            "fluid_pp(1)%jwl_R2": jwl_R2,
+            "fluid_pp(1)%jwl_omega": jwl_omega,
+            "fluid_pp(1)%jwl_rho0": jwl_rho0,
+            "fluid_pp(1)%jwl_E0": jwl_E0,
+            "fluid_pp(1)%jwl_air_p0": 101325.0,
+            "fluid_pp(1)%jwl_air_rho0": 1000.0,
+            # Fluid 2: water (stiffened gas).
+            "fluid_pp(2)%eos": 1,
+            "fluid_pp(2)%gamma": water_gamma_mfc,
+            "fluid_pp(2)%pi_inf": water_pi_inf_mfc,
+            "fluid_pp(2)%cv": 4186.0,
+        }
+    )
+)