JWL equation of state: single-material and multi-fluid mixture implementation

README-JWL-EOS.md

@@ -0,0 +1,32 @@
+# JWL EOS Notes
+
+This page documents the Jones-Wilkins-Lee (JWL) equation-of-state support used by the five-equation model in MFC.
+
+## Pure JWL Products
+
+For products density `rho`, specific internal energy `e`, and relative volume `V = rho0/rho`, MFC evaluates
+
+```text
+p = A (1 - omega/(R1 V)) exp(-R1 V)
+  + B (1 - omega/(R2 V)) exp(-R2 V)
+  + omega rho e.
+```
+
+The cold-curve pressure is evaluated through a shared helper so pressure recovery, energy recovery, and sound-speed code use the same expression.
+
+## 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.
+
+- `0`, `isobaric`: closed-form mechanical-equilibrium closure. This is the default validation path.
+- `1`, `kuhl`: Kuhl/Khasainov temperature-form additive closure. It requires positive `cv` on both products and air.
+- `2`, `ptequil`: pressure-temperature equilibrium closure. It solves a bounded scalar root for the products volume fraction and is substantially more expensive than mode `0`.
+- `3`, `rocflu`: Garno/Rocflu-style single-fluid blend. Its sound speed is evaluated by the Rocflu-specific Gruneisen form rather than by phasic volume fractions.
+
+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
+
+The exact-reference validation in this branch is scoped to five-equation JWL cases: a 1D pure-JWL shock tube and a compact 3D quasi-1D repeat of the same published Shyue-style Riemann reference. Closure selectability is covered by registered golden tests for `jwl_mix_type = 0, 1, 2, 3`.
+
+The p-T equilibrium mode remains selectable and physically distinct. Its root-find cost is measured by `benchmarks/jwl_closure_modes`; it is not optimized away by the case-optimization gates.

benchmarks/jwl_closure_modes/README.md

@@ -0,0 +1,23 @@
+# JWL Closure Modes 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.
+
+Run the baseline isobaric closure and the p-T equilibrium closure with the same build options:
+
+```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
+```
+
+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`.
+
+## 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 |
+
+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.

benchmarks/jwl_closure_modes/case.py

@@ -0,0 +1,87 @@
+#!/usr/bin/env python3
+import argparse
+import json
+
+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
+rho_jwl = 1630.0
+rho_air = 1.225
+
+params = {
+    "run_time_info": "T",
+    "x_domain%beg": 0.0,
+    "x_domain%end": 1.0,
+    "m": 399,
+    "n": 0,
+    "p": 0,
+    "dt": 5.0e-8,
+    "t_step_start": 0,
+    "t_step_stop": 200,
+    "t_step_save": 200,
+    "num_patches": 2,
+    "model_eqns": 2,
+    "num_fluids": 2,
+    "jwl_mix_type": args.mix_type,
+    "mpp_lim": "T",
+    "mixture_err": "T",
+    "time_stepper": 3,
+    "recon_type": 1,
+    "weno_order": 3,
+    "weno_eps": 1.0e-16,
+    "mapped_weno": "T",
+    "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",
+    "rho_wrt": "T",
+    "pres_wrt": "T",
+    "c_wrt": "T",
+    "parallel_io": "F",
+    "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)": eps * rho_jwl,
+    "patch_icpp(1)%alpha_rho(2)": (1.0 - eps) * rho_air,
+    "patch_icpp(1)%alpha(1)": eps,
+    "patch_icpp(1)%alpha(2)": 1.0 - eps,
+    "patch_icpp(2)%geometry": 1,
+    "patch_icpp(2)%alter_patch(1)": "T",
+    "patch_icpp(2)%x_centroid": 0.15,
+    "patch_icpp(2)%length_x": 0.3,
+    "patch_icpp(2)%vel(1)": 0.0,
+    "patch_icpp(2)%pres": 1.2e10,
+    "patch_icpp(2)%alpha_rho(1)": (1.0 - eps) * rho_jwl,
+    "patch_icpp(2)%alpha_rho(2)": eps * rho_air,
+    "patch_icpp(2)%alpha(1)": 1.0 - eps,
+    "patch_icpp(2)%alpha(2)": eps,
+    "fluid_pp(1)%eos": 2,
+    "fluid_pp(1)%gamma": 2.5,
+    "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": rho_jwl,
+    "fluid_pp(1)%jwl_E0": 1.0089e10,
+    "fluid_pp(1)%jwl_air_e0": 2.5575e5,
+    "fluid_pp(1)%jwl_air_rho0": rho_air,
+    "fluid_pp(1)%jwl_air_gamma": 0.4,
+    "fluid_pp(2)%eos": 1,
+    "fluid_pp(2)%gamma": 2.5,
+    "fluid_pp(2)%pi_inf": 0.0,
+    "fluid_pp(2)%cv": 717.5,
+}
+
+print(json.dumps(params))

docs/documentation/case.md

@@ -443,6 +443,8 @@ 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`; modes 1 (Kuhl) and 2 (p-T equilibrium) additionally require a positive `cv` on both fluids.
+
 > **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
 > - `pi_inf` = \f$\gamma\,\pi_\infty / (\gamma - 1)\f$, not \f$\pi_\infty\f$ itself
@@ -461,6 +463,7 @@ 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, verified against an exact Riemann reference); [1] Kuhl, [2] p-T equilibrium, [3] Rocflu blend (experimental) |
 | `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 |

docs/module_categories.json

@@ -19,6 +19,7 @@
     {
         "category": "Physics Models",
         "modules": [
+            "m_jwl",
             "m_viscous",
             "m_hb_function",
             "m_surface_tension",

examples/1D_jwl_single_material_shocktube/README.md

@@ -0,0 +1,37 @@
+# 1D JWL Single-Material Shock Tube
+
+This case is a 1D single-fluid JWL shock tube intended for quick verification against a
+reference/exact Riemann workflow.
+
+## Initial conditions
+
+- Left state: `rho = 1770 kg/m^3`, `u = 0 m/s`, `p = 30 GPa`
+- Right state: `rho = 1770 kg/m^3`, `u = 0 m/s`, `p = 10 GPa`
+
+The setup is generated by `case.py` and uses MFC's JWL EOS (`fluid_pp(1)%eos = 2`).
+
+## Run
+
+```bash
+./mfc.sh run examples/1D_jwl_single_material_shocktube/case.py
+```
+
+## Verification
+
+The single-material JWL shock-tube path was verified against an exact Mie-Grüneisen
+JWL Riemann solution under grid refinement (Shyue, Journal of Computational Physics,
+2001). Errors are
+relative L1 norms versus the exact solution; the star state is `p* = 11.53 GPa`,
+`u* = 2026 m/s`.
+
+| Nx   | L1(ρ) [%] | L1(u) [%] | L1(p) [%] | mass drift |
+|------|-----------|-----------|-----------|------------|
+| 100  | 1.885     | 3.044     | 1.287     | 5.5e-10    |
+| 200  | 0.930     | 1.437     | 0.621     | 4.6e-16    |
+| 400  | 0.514     | 0.814     | 0.317     | 0.0        |
+| 800  | 0.251     | 0.369     | 0.158     | 3.4e-16    |
+| 1600 | 0.128     | 0.188     | 0.085     | 3.4e-16    |
+| 3200 | 0.069     | 0.103     | 0.043     | 0.0        |
+
+The global convergence rate is ~1 (shock-limited, as expected for a discontinuous
+solution); mass is conserved to machine precision for Nx ≥ 200.

examples/1D_jwl_single_material_shocktube/case.py

@@ -0,0 +1,67 @@
+#!/usr/bin/env python3
+import json
+
+rho0 = 1770.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-8,
+            "t_step_start": 0,
+            "t_step_stop": 600,
+            "t_step_save": 60,
+            "num_patches": 2,
+            "model_eqns": 2,
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "recon_type": 2,
+            "muscl_order": 2,
+            "muscl_lim": 2,
+            "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_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.25,
+            "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(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(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)%jwl_air_e0": 2.5575e5,
+            "fluid_pp(1)%jwl_air_rho0": 1.225,
+            "fluid_pp(1)%jwl_air_gamma": 0.4,
+        }
+    )
+)