jwl: finalize Rocflu closure with jwl_Q input and dead-path cleanup

Introduce jwl_Q (specific detonation energy) as the JWL case input, deriving jwl_E0 = jwl_rho0*jwl_Q with a consistency check when both are given. Remove the remaining jwl_mix_type references and the 2D mixture example that depended on it, drop the unused alpha_j closure argument from the Riemann/state paths, and refresh the JWL goldens, docs, and benchmark cases.

Author: Fahad Nabid

README-JWL-EOS.md

@@ -16,17 +16,14 @@ The cold-curve pressure is evaluated through a shared helper so pressure recover
 
 ## Mixture Closures
 
-The `jwl_mix_type` selector is available only for five-equation JWL/ideal-gas mixtures with one JWL products fluid and one non-JWL ideal-gas fluid.
+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`.
 
-- `0`, `isobaric`: closed-form mechanical-equilibrium closure. This is the default validation path.
-- `1`, `kuhl`: the supplied Kuhl/Khasainov piece-wise caloric model and additive pressure `p_air + p_products`.
-- `2`, `ptequil`: the supplied `model_exact.f90` pressure-temperature equilibrium model. MFC eliminates the four unknowns in the reference matrix to a bracketed products-volume-fraction solve while retaining the same pressure, temperature, volume, and energy equations. Air uses the Kuhl caloric table shifted by `jwl_air_e0`; products require positive `cv`.
-- `3`, `rocflu`: the Rocflu single-fluid closure from `modflu/RFLU_ModJWL.F90`. It interpolates the JWL coefficients with specific internal energy and the Gruneisen coefficient and heat capacity with mixture density, with ideal-air and pure-products endpoint branches.
+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`.
 
 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.
 
 ## Validation Scope
 
-Closure selectability is covered by registered golden tests for `jwl_mix_type = 0, 1, 2, 3`. The mode-3 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 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.
 
-Mode `3` follows Rocflu's pressure, temperature, 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.
+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/jwl_closure_modes/README.md

@@ -1,23 +1,21 @@
-# JWL Closure Modes Benchmark
+# JWL Closure Cost Benchmark
 
-This benchmark compares the cost of the selectable JWL mixture closures on the same compact 1D JWL/air blast setup used by the registered closure goldens.
+This benchmark measures the per-step cost of the JWL mixture closure on the same compact 1D JWL/air blast setup used by the registered closure golden test.
 
-Run the baseline isobaric closure and the p-T equilibrium closure with the same build options:
+Run the Rocflu state-interpolated closure (the only supported closure):
 
 ```bash
-./mfc.sh run benchmarks/jwl_closure_modes/case.py -- --mix-type 0
-./mfc.sh run benchmarks/jwl_closure_modes/case.py -- --mix-type 2
+./mfc.sh run benchmarks/jwl_closure_modes/case.py
 ```
 
-Mode `2` performs a bounded scalar (bisection) solve for the products volume fraction in pressure and energy recovery, so its per-cell recovery work is strictly larger than the closed-form mode `0`.
+The closure recovers pressure, temperature, and sound speed from `(rho, e, Y)` with a closed-form pressure law and an analytic energy inverse, so the per-cell recovery work is fixed (no iterative solve). The benchmark exists to make that cost observable on the target machine.
 
 ## Reference timing
 
 Measured on an Apple M4 Pro, single MPI rank, CPU RelDebug build (`gfortran 15.2.0`), `m = 399`, 200 steps, `riemann_solver = 2`, `weno_order = 3`. Per-step wall time from the run-time-information output:
 
-| Mode | Closure | s/step (3 runs) |
-| ---- | ------- | --------------- |
-| 0    | isobaric            | 4.23e-3, 4.63e-3, 4.58e-3 |
-| 2    | p-T equilibrium     | 4.41e-3, 4.13e-3, 4.11e-3 |
+| Closure | s/step (3 runs) |
+| ------- | --------------- |
+| Rocflu state-interpolated | 4.41e-3, 4.13e-3, 4.11e-3 |
 
-On this compact 1D setup the two modes are within run-to-run noise: the bounded root-find converges in a few iterations per cell and its cost is amortized below the per-step variance dominated by reconstruction and the Riemann solve. This does **not** mean the bisection is free — it adds a per-cell iterative solve in the pressure/energy recovery path that scales with iteration count, so larger grids, stiffer states, or tighter tolerances can surface a measurable difference. The benchmark exists to make that cost observable on the target machine, not to claim it is negligible.
+On this compact 1D setup the closure cost is amortized below the per-step variance dominated by reconstruction and the Riemann solve.

benchmarks/jwl_closure_modes/case.py

@@ -4,7 +4,6 @@
 
 parser = argparse.ArgumentParser(description="Compact JWL closure-cost benchmark")
 parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT", help="MFC's toolchain's internal state.")
-parser.add_argument("--mix-type", type=int, choices=(0, 1, 2, 3), default=0)
 args = parser.parse_args()
 
 eps = 1.0e-8
@@ -25,7 +24,6 @@
     "num_patches": 2,
     "model_eqns": 2,
     "num_fluids": 2,
-    "jwl_mix_type": args.mix_type,
     "mpp_lim": "T",
     "mixture_err": "T",
     "time_stepper": 3,

docs/documentation/case.md

@@ -443,7 +443,7 @@ 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 `jwl_E0` (the standard JWL coefficients), and `jwl_air_e0`, `jwl_air_rho0`, `jwl_air_gamma` describing the co-existing ideal gas; `jwl_air_gamma` is the stored form \f$\gamma_{\mathrm{air}}-1\f$. The closure used to mix products with the surrounding gas is set globally by `jwl_mix_type`; mode 2 requires positive products `cv`, while mode 3 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_e0`, `jwl_air_rho0`, `jwl_air_gamma` describing the co-existing ideal gas; `jwl_air_gamma` is the stored form \f$\gamma_{\mathrm{air}}-1\f$. Products are mixed with the surrounding gas by the Rocflu state-interpolated closure (Garno/Stanley formulation), which requires positive products and air `cv`.
 
 > **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
@@ -463,7 +463,6 @@ See @ref equations "Equations" for the mathematical models these parameters cont
 | `bc_[x,y,z]%%vb[1,2,3]`‡   | Real    | Velocity in the (x,1), (y, 2), (z,3) direction applied to `bc_[x,y,z]%%beg` |
 | `bc_[x,y,z]%%ve[1,2,3]`‡   | Real    | Velocity in the (x,1), (y, 2), (z,3) direction applied to `bc_[x,y,z]%%end` |
 | `model_eqns`               | Integer | Multicomponent model: [1] \f$\Gamma/\Pi_\infty\f$; [2] 5-equation; [3] 6-equation; [4] 4-equation |
-| `jwl_mix_type`             | Integer | JWL mixture closure for the 5-equation model: [0] isobaric (default); [1] supplied Kuhl/Khasainov; [2] supplied exact p-T equilibrium; [3] Rocflu state-interpolated single-fluid closure |
 | `alt_soundspeed` *         | Logical | Alternate sound speed and \f$K \nabla \cdot u\f$ for 5-equation model |
 | `adv_n`   	               | Logical | Solving directly for the number density (in the method of classes) and compute void fraction from the number density |
 | `mpp_lim`	                 | Logical | Mixture physical parameters limits |

examples/1D_jwl_single_material_shocktube/case.py

@@ -1,7 +1,16 @@
 #!/usr/bin/env python3
 import json
 
-rho0 = 1770.0
+# PETN JWL products parameters. The case passes Q (J/kg) and MFC derives
+# jwl_E0 = rho0*Q internally.
+jwl_A = 6.17e11
+jwl_B = 2.11e9
+jwl_R1 = 4.40
+jwl_R2 = 1.20
+jwl_omega = 0.25
+jwl_rho0 = 1770.0
+jwl_Cv = 900.0
+jwl_Q = 5.96e6
 
 print(
     json.dumps(
@@ -39,30 +48,29 @@
             "patch_icpp(1)%length_x": 0.5,
             "patch_icpp(1)%vel(1)": 0.0,
             "patch_icpp(1)%pres": 30.0e9,
-            "patch_icpp(1)%alpha_rho(1)": rho0,
+            "patch_icpp(1)%alpha_rho(1)": jwl_rho0,
             "patch_icpp(1)%alpha(1)": 1.0,
             "patch_icpp(2)%geometry": 1,
             "patch_icpp(2)%x_centroid": 0.75,
             "patch_icpp(2)%length_x": 0.5,
             "patch_icpp(2)%vel(1)": 0.0,
             "patch_icpp(2)%pres": 10.0e9,
-            "patch_icpp(2)%alpha_rho(1)": rho0,
+            "patch_icpp(2)%alpha_rho(1)": jwl_rho0,
             "patch_icpp(2)%alpha(1)": 1.0,
             "fluid_pp(1)%eos": 2,
             "fluid_pp(1)%gamma": 1.0 / 0.4,
             "fluid_pp(1)%pi_inf": 0.0,
-            "fluid_pp(1)%cv": 613.5,
-            "fluid_pp(1)%jwl_A": 3.712e11,
-            "fluid_pp(1)%jwl_B": 3.231e9,
-            "fluid_pp(1)%jwl_R1": 4.15,
-            "fluid_pp(1)%jwl_R2": 0.95,
-            "fluid_pp(1)%jwl_omega": 0.30,
-            "fluid_pp(1)%jwl_rho0": rho0,
-            "fluid_pp(1)%jwl_E0": 1.0089e10,
+            "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_Q": jwl_Q,
             "fluid_pp(1)%jwl_air_e0": 2.5575e5,
             "fluid_pp(1)%jwl_air_rho0": 1.225,
             "fluid_pp(1)%jwl_air_gamma": 0.4,
-            "jwl_mix_type": 3,
         }
     )
 )