jwl: harden Rocflu closure and correct pre_process t=0 pressure diagnostic

- m_jwl: derive Gamma_air from the ideal-gas fluid, smoothstep omega(rho) and the Y blend, use ej_rho_ref, mass-weighted cv, and a c^2 floor; scan (rho,e,Y) at init.
- HLLC/HLL: use the closure-consistent averaged-state pressure.
- pre_process serial writer: pass the cell's own jwl_Y into the t=0 primitive-pressure recovery so pure-air cells use Gamma_air instead of the products omega (was writing 0.75x for air; conserved energy was always correct).
- Add 1D JWL/air interface-advection mechanical-equilibrium guard (quiescent, max|dp|/p0 ~ 1e-16) and a JWL/air mixture test; regenerate the affected golden files.

Author: Fahad Nabid

README-JWL-EOS.md

@@ -16,14 +16,16 @@ 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 interpolates the JWL `A` and `B` coefficients with specific internal energy and interpolates `omega` and heat capacity with mixture density. The closure uses adapted mixture coefficients through products mass fraction `Y <= 0.99` and pure-products JWL coefficients for `Y > 0.99`.
+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.
 
-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, 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. 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_rho0 > jwl_air_e0`.
+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, energy, and sound-speed floors are applied only after explicit finite checks. NaNs are intentionally preserved so bad states are visible during debugging instead of being converted into plausible-looking floor values.
+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.
 
 ## 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.
 
+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.
+
 The closure follows Rocflu's pressure, temperature, inverse-energy, and sound-speed formulas, but replaces its case-specific hard-coded air values and explosive energy divisor with the corresponding MFC material inputs. Its inverse energy selects the exact low-, blended-, or high-energy branch of that pressure law; this removes the legacy fallback's pressure/energy round-trip mismatch.

benchmarks/3D_jwl_spherical_tnt_free_air_validation/README.md

@@ -0,0 +1,244 @@
+# 3D JWL Spherical TNT Free-Air Benchmark
+
+This benchmark is a 3D JWL/TNT free-air blast candidate for MFC. A spherical
+region of TNT detonation products expands into ambient air, and native point
+probes record pressure-time histories for arrival time, peak incident
+overpressure, and positive-phase impulse extraction.
+
+Reference values are intentionally not fabricated. The MFC probe extraction is
+implemented and the case runs, but the published reference values still need to
+be filled from the paper/CONWEP table before this should be described as
+completed validation.
+
+## Citation
+
+Giam, Toh, and Tan, "Numerical Review of Jones-Wilkins-Lee Parameters for
+Trinitrotoluene Explosive in Free-Air Blast," Journal of Applied Mechanics,
+2020. DOI: `10.1115/1.4046243`.
+
+The case is motivated by that free-air TNT/JWL review, but this repository does
+not currently include an accessible table of the paper's exact benchmark
+values. Until those values, Kingery-Bulmash/CONWEP values, or another trusted
+numeric free-air TNT reference are added, this is a benchmark candidate rather
+than completed validation.
+
+## Physics Scope
+
+This benchmark initializes TNT detonation products directly and validates 3D JWL product-air blast expansion. It does not model detonation initiation, reaction-zone structure, afterburn, or structural coupling.
+
+The benchmark exercises:
+
+- 3D JWL product-air expansion.
+- Radial free-air blast propagation in a Cartesian grid.
+- Native MFC pressure probes and post-processing.
+- A comparison framework for published or standard free-air blast quantities.
+
+The benchmark does not validate detonation initiation, reaction-zone physics,
+afterburn, structural coupling, confined blast, tunnel blast, or ground
+reflection.
+
+## Why This Is Genuinely 3D
+
+The default case resolves an octant of a spherical products region in a 3D
+Cartesian domain. The shock expands in x, y, and z, with active transverse
+reconstruction and Riemann fluxes. This is not a 1D shock tube extruded through
+the 3D solver path.
+
+Symmetry planes at `x = 0`, `y = 0`, and `z = 0` recover the full spherical
+solution while reducing the cell count by a factor of about eight. The probes
+lie near the positive x radial direction, at the nearest cell centers to the
+target radii. Keeping probes on cell centers avoids boundary and
+MPI-decomposition interpolation artifacts on the coarse default grid.
+
+## Geometry
+
+| Quantity | Value |
+|---|---:|
+| Charge center | `(0, 0, 0) m` in the octant representation |
+| Charge radius | `0.05 m` |
+| TNT/product reference density | `1630 kg/m^3` |
+| Full-sphere TNT mass | `0.853466 kg` |
+| Domain | `x, y, z in [0, 0.5] m` |
+| Final time | `2.0e-4 s` |
+
+The simulated octant contains one eighth of the full sphere. The symmetry
+boundaries make the resolved field equivalent to the full charge, so scaled
+distance uses the full TNT mass:
+
+```text
+W = (4/3) pi (0.05 m)^3 (1630 kg/m^3) = 0.853466 kg
+Z = r / W^(1/3)
+```
+
+## Grid
+
+Default local grid:
+
+| Quantity | Value |
+|---|---:|
+| `m = n = p` | `63` |
+| Cells in octant | `64^3` |
+| Cell size | `0.0078125 m` |
+| Reconstruction | mapped WENO3 |
+| Riemann solver | HLLC |
+| Time stepper | RK3 |
+| CFL target | `0.3` |
+
+Lower-resolution smoke run:
+
+```bash
+./mfc.sh run benchmarks/3D_jwl_spherical_tnt_free_air_validation/case.py -n 4 -- --grid 31
+python3 benchmarks/3D_jwl_spherical_tnt_free_air_validation/gauges.py --grid 31
+```
+
+The `64^3` octant default gives a visibly cleaner spherical shock than the
+older `32^3` smoke grid, but it takes several minutes and writes larger output
+files.
+
+## Boundary Conditions
+
+| Boundary | MFC value | Meaning |
+|---|---:|---|
+| `bc_x%beg`, `bc_y%beg`, `bc_z%beg` | `-2` | symmetry planes |
+| `bc_x%end`, `bc_y%end`, `bc_z%end` | `-3` | non-reflecting/open boundaries |
+
+## EOS And Initial Conditions
+
+Ambient air fills the full octant before the products sphere is overlaid:
+
+| Parameter | Value |
+|---|---:|
+| Air EOS | ideal gas/stiffened gas, `eos = 1` |
+| Air pressure | `101325 Pa` |
+| Air density | `1.225 kg/m^3` |
+| Air physical gamma | `1.4` |
+| MFC air `fluid_pp(2)%gamma` | `2.5` |
+| Air sound speed | `340.3 m/s` |
+| Air `cv` | `717.5 J/(kg K)` |
+
+TNT products use the JWL constants already present in the MFC JWL examples and
+benchmarks. They are documented here as the repo-local parameter source, not as
+a claim that these are the exact Giam et al. tabulated values:
+
+| Parameter | Value |
+|---|---:|
+| Products EOS | JWL, `eos = 2` |
+| Products density/reference density `rho0` | `1630 kg/m^3` |
+| Products volume-fraction floor | `1.0e-6` |
+| JWL `A` | `3.712e11 Pa` |
+| JWL `B` | `3.231e9 Pa` |
+| JWL `R1` | `4.15` |
+| JWL `R2` | `0.95` |
+| JWL `omega` | `0.30` |
+| Initial specific energy equivalent | `E0/rho0 = 6.1908e6 J/kg` |
+| Initial internal energy density `E0` | `1.0089e10 J/m^3` |
+| JWL cold pressure at `rho = rho0` | `6.2837e9 Pa` |
+| Initial products pressure | `9.3104e9 Pa` |
+| Products `cv` | `613.5 J/(kg K)` |
+| `jwl_air_e0` | `2.5575e5 J/kg` |
+| `jwl_air_rho0` | `1.225 kg/m^3` |
+| air Grüneisen (derived from `fluid_pp(2)%gamma`) | `1/2.5 = 0.4` |
+
+The initialized products pressure is computed in `case.py` from:
+
+```text
+p = A (1 - omega/(R1 V)) exp(-R1 V)
+  + B (1 - omega/(R2 V)) exp(-R2 V)
+  + omega E0,
+V = rho0 / rho = 1.
+```
+
+## Probe Locations
+
+The default `--grid 63` run uses the nearest cell centers to target x positions
+`0.15, 0.25, 0.35, 0.45 m`:
+
+| Gauge | x (m) | y (m) | z (m) | r (m) | Z (m/kg^(1/3)) |
+|---:|---:|---:|---:|---:|---:|
+| 1 | 0.152344 | 0.003906 | 0.003906 | 0.152444 | 0.160712 |
+| 2 | 0.253906 | 0.003906 | 0.003906 | 0.253966 | 0.267741 |
+| 3 | 0.347656 | 0.003906 | 0.003906 | 0.347700 | 0.366558 |
+| 4 | 0.449219 | 0.003906 | 0.003906 | 0.449253 | 0.473618 |
+
+MFC's native probe interpolation can write a local initial pressure that differs
+from the analytic ambient state on this coarse octant grid. The reducer therefore
+uses each probe's first sample as its local baseline and reports pressure rise
+relative to that baseline.
+
+Arrival time is the first time the pressure rise exceeds 5 percent of ambient:
+
+```text
+p(t) - p(t=0) > 0.05 p_ambient
+```
+
+Peak incident overpressure is estimated from the same baseline-corrected trace:
+
+```text
+dp_peak = max(p(t) - p(t=0))
+```
+
+Positive-phase impulse is integrated over the baseline-corrected history:
+
+```text
+I = integral max(p(t) - p(t=0), 0) dt
+```
+
+## Run Commands
+
+From the repository root:
+
+```bash
+./mfc.sh run benchmarks/3D_jwl_spherical_tnt_free_air_validation/case.py -n 4
+python3 benchmarks/3D_jwl_spherical_tnt_free_air_validation/gauges.py
+./mfc.sh precheck
+```
+
+The probe reducer reads:
+
+```text
+benchmarks/3D_jwl_spherical_tnt_free_air_validation/D/probe<i>_prim.dat
+```
+
+and writes:
+
+```text
+benchmarks/3D_jwl_spherical_tnt_free_air_validation/gauge_results.csv
+```
+
+## Comparison Table
+
+| Gauge | r (m) | Z (m/kg^(1/3)) | MFC arrival (s) | Ref arrival (s) | MFC peak dp (Pa) | Ref peak dp (Pa) | MFC impulse (Pa s) | Ref impulse (Pa s) | Error |
+|---|---:|---:|---:|---:|---:|---:|---:|---:|---:|
+| 1 | 0.152444 | 0.160712 | 8.000000e-06 | pending reference | 9.189414e+06 | pending reference | 1.142617e+03 | pending reference | pending |
+| 2 | 0.253966 | 0.267741 | 3.700000e-05 | pending reference | 7.494431e+06 | pending reference | 3.127588e+02 | pending reference | pending |
+| 3 | 0.347700 | 0.366558 | 7.200000e-05 | pending reference | 5.750937e+06 | pending reference | 2.453974e+02 | pending reference | pending |
+| 4 | 0.449253 | 0.473618 | 1.150000e-04 | pending reference | 4.023872e+06 | pending reference | 3.029903e+02 | pending reference | pending |
+
+Reference values may be filled from Giam et al. if numeric values are available,
+from Kingery-Bulmash/CONWEP free-air TNT data if used for comparison, or from
+another clearly cited spherical free-air TNT reference. If values are digitized
+from a figure, label them as digitized estimates.
+
+## Pass/Fail Expectations
+
+This candidate benchmark passes its local sanity checks when:
+
+- The case completes without NaNs.
+- Density, pressure, and internal energy remain positive.
+- Shock arrival time increases monotonically with radius.
+- Peak incident overpressure decreases monotonically with radius.
+- The pressure histories show a radially reasonable outward blast.
+
+Quantitative error should be reported only after trusted reference arrival,
+peak incident overpressure, and impulse values are added.
+
+## Limitations And Remaining Work
+
+- This is not completed validation until the reference columns are populated.
+- A `32^3` octant run is available as a faster smoke test, but the default is now `64^3`.
+- The products sphere initialization does not model a detonation wave or
+  reaction zone.
+- No afterburn, structural coupling, confinement, or ground reflection is
+  represented.
+- A higher-resolution run should be used before drawing quantitative
+  conclusions.