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

.typos.toml

@@ -31,4 +31,4 @@ tru = "tru"           # typo for "true" in "when_tru" - tests dependency keys
 PNGs = "PNGs"
 
 [files]
-extend-exclude = ["docs/documentation/references*", "docs/references.bib", "tests/", "toolchain/cce_simulation_workgroup_256.sh", "build-docs/", "build/", "build_test/", "fp-stability-logs/"]
+extend-exclude = ["docs/documentation/references*", "docs/references.bib", "tests/", "toolchain/cce_simulation_workgroup_256.sh", "build-docs/", "build/", "build_test/", "fp-stability-logs/", "_jwl_verification/", "_jwl_organization/"]

docs/documentation/case.md

@@ -440,6 +440,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: [0] isobaric; [1] Kuhl; [2] p-T equilibrium; [3] Rocflu blend |
 | `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_cj_products/case.py

@@ -0,0 +1,124 @@
+#!/usr/bin/env python3
+"""1D five-equation JWL products with a moving detonation-product state.
+
+The left state is initialized from a simple Rankine-Hugoniot product estimate
+behind a right-running detonation:
+
+    u_p = D * (1 - rho_0 / rho_products)
+    p_products = p_0 + rho_0 * D * u_p
+
+Uses the five-equation model (Allaire et al. JCP 2002) with the isobaric JWL
+mixture closure (jwl_mix_type = 0). The numbers are TNT-like and intended as a
+solver/debugging benchmark, not a calibrated explosive model.
+"""
+
+import json
+
+eps = 1.0e-8
+
+rho0_jwl = 1630.0
+rho_air = 1.225
+p0 = 101325.0
+
+detonation_speed = 6900.0
+compression_ratio = 1.80
+rho_products = compression_ratio * rho0_jwl
+u_products = detonation_speed * (1.0 - rho0_jwl / rho_products)
+p_products = p0 + rho0_jwl * detonation_speed * u_products
+
+jwl = {
+    "A": 3.712e11,
+    "B": 3.231e9,
+    "R1": 4.15,
+    "R2": 0.95,
+    "omega": 0.30,
+    "rho0": rho0_jwl,
+    "E0": 1.0089e10,
+    "air_e0": 2.5575e5,
+    "air_rho0": rho_air,
+    "air_gamma": 0.4,
+}
+
+print(
+    json.dumps(
+        {
+            "run_time_info": "T",
+            "x_domain%beg": 0.0,
+            "x_domain%end": 1.0,
+            "m": 399,
+            "n": 0,
+            "p": 0,
+            "t_step_start": 0,
+            "dt": 1.0e-12,
+            "t_step_stop": 8,
+            "t_step_save": 8,
+            "n_start": 0,
+            "num_patches": 2,
+            "model_eqns": 2,
+            "num_fluids": 2,
+            "jwl_mix_type": 0,
+            "mpp_lim": "T",
+            "mixture_err": "T",
+            "time_stepper": 3,
+            "recon_type": 2,
+            "muscl_order": 2,
+            "muscl_lim": 2,
+            "riemann_solver": 1,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "cons_vars_wrt": "F",
+            "rho_wrt": "T",
+            "pres_wrt": "T",
+            "E_wrt": "T",
+            "c_wrt": "T",
+            "parallel_io": "F",
+            # Ambient air.
+            "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": p0,
+            "patch_icpp(1)%alpha_rho(1)": eps * rho0_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,
+            # Moving detonation products estimated from the jump conditions.
+            "patch_icpp(2)%geometry": 1,
+            "patch_icpp(2)%alter_patch(1)": "T",
+            "patch_icpp(2)%x_centroid": 0.10,
+            "patch_icpp(2)%length_x": 0.20,
+            "patch_icpp(2)%vel(1)": u_products,
+            "patch_icpp(2)%pres": p_products,
+            "patch_icpp(2)%alpha_rho(1)": (1.0 - eps) * rho_products,
+            "patch_icpp(2)%alpha_rho(2)": eps * rho_air,
+            "patch_icpp(2)%alpha(1)": 1.0 - eps,
+            "patch_icpp(2)%alpha(2)": eps,
+            # Fluid 1: fully reacted JWL products.
+            "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": 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_e0": jwl["air_e0"],
+            "fluid_pp(1)%jwl_air_rho0": jwl["air_rho0"],
+            "fluid_pp(1)%jwl_air_gamma": jwl["air_gamma"],
+            # Fluid 2: ideal-gas air.
+            "fluid_pp(2)%eos": 1,
+            "fluid_pp(2)%gamma": 1.0 / 0.4,
+            "fluid_pp(2)%pi_inf": 0.0,
+            "fluid_pp(2)%cv": 717.5,
+        },
+        indent=2,
+    )
+)

examples/1D_jwl_mixture_test/README.md

@@ -0,0 +1,85 @@
+# 1D JWL Mixture Smoke Test
+
+This is the smallest practical test for the mixed JWL/air EOS path. It keeps the geometry simple so that failures are easy to diagnose: a JWL-rich high-pressure driver expands into mostly ideal-gas air.
+
+The case is intentionally modest. It is not a calibrated detonation problem and it is not meant to match explosive test data. It is a quick check that the new EOS path can initialize pressure, convert between primitive and conservative variables, compute sound speed, march a few steps, and post-process output without producing bad thermodynamic states.
+
+## What It Exercises
+
+The case uses:
+
+```text
+model_eqns  2
+num_fluids  2
+fluid 1     JWL explosive products
+fluid 2     ideal-gas air
+```
+
+The JWL fluid is selected with:
+
+```python
+"fluid_pp(1)%eos": 2
+```
+
+Air remains the standard ideal/stiffened-gas path:
+
+```python
+"fluid_pp(2)%eos": 1
+```
+
+The mixed cells use a tiny volume-fraction floor, `eps = 1.0e-8`, to avoid exactly zero volume fraction for either material. That is a numerical safety device for the mixture model, not a physical interface thickness.
+
+The exact EOS used by this case is documented in the root `README-JWL-EOS.md`. In short, cells with `Y_JWL <= 1.0e-2` use the ideal-gas air fallback, partially mixed cells use effective JWL constants, and JWL-rich cells use the standard JWL pressure form.
+
+## Run
+
+From the MFC repository root:
+
+```bash
+./mfc.sh run examples/1D_jwl_mixture_test/case.py --no-build
+```
+
+If the build is stale:
+
+```bash
+./mfc.sh run examples/1D_jwl_mixture_test/case.py
+```
+
+## Expected Result
+
+A successful run completes:
+
+```text
+syscheck
+pre_process
+simulation
+post_process
+```
+
+and exits with code `0`.
+
+Because this is a smoke test, the important result is not a particular shock location. The important result is that pressure, density, energy, and sound speed stay finite and the run completes.
+
+## Parameters
+
+The JWL constants are TNT-like:
+
+```text
+A      3.712e11 Pa
+B      3.231e9 Pa
+R1     4.15
+R2     0.95
+omega  0.30
+rho0   1630 kg/m3
+E0     1.0089e10 J/kg
+```
+
+The mixed JWL/air helper also uses reference air values:
+
+```text
+air_e0      2.5575e5 J/kg
+air_rho0    1.225 kg/m3
+air_gamma   0.4
+```
+
+In this branch, `jwl_air_gamma` is the `gamma - 1` coefficient used by the current mixed JWL/air helper, matching the convention used in the implementation.

examples/1D_jwl_mixture_test/case.py

@@ -0,0 +1,112 @@
+#!/usr/bin/env python3
+import json
+import os
+
+eps = 1.0e-8
+
+# Mixture closure to benchmark: 0 isobaric, 1 Kuhl, 2 p-T equilibrium, 3 Rocflu blend.
+jwl_mix_type = int(os.environ.get("JWL_MIX_TYPE", "0"))
+
+rho_jwl = 1630.0
+rho_air = 1.225
+
+jwl = {
+    "A": 3.712e11,
+    "B": 3.231e9,
+    "R1": 4.15,
+    "R2": 0.95,
+    "omega": 0.30,
+    "rho0": rho_jwl,
+    "E0": 1.0089e10,
+    "air_e0": 2.5575e5,
+    "air_rho0": rho_air,
+    "air_gamma": 0.4,
+}
+
+print(
+    json.dumps(
+        {
+            "run_time_info": "T",
+            # Domain
+            "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": 600,
+            "t_step_save": 150,
+            # Numerics
+            "num_patches": 2,
+            "model_eqns": 2,
+            "num_fluids": 2,
+            "jwl_mix_type": jwl_mix_type,
+            "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,
+            # Output
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "cons_vars_wrt": "F",
+            "rho_wrt": "T",
+            "pres_wrt": "T",
+            "c_wrt": "T",
+            "parallel_io": "F",
+            # Patch 1: mostly air background
+            "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 2: hot JWL detonation-products driver. p must exceed the JWL cold pressure at rho0
+            # (~7.1 GPa here) so the products carry positive internal energy / temperature.
+            "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.30,
+            "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 1: JWL products
+            "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": 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_e0": jwl["air_e0"],
+            "fluid_pp(1)%jwl_air_rho0": jwl["air_rho0"],
+            "fluid_pp(1)%jwl_air_gamma": jwl["air_gamma"],
+            # Fluid 2: ideal-gas air (cv set so p-T equilibrium has R_air = (gamma-1)*cv = 0.4*717.5 = 287 J/kg/K)
+            "fluid_pp(2)%eos": 1,
+            "fluid_pp(2)%gamma": 1.0 / 0.4,
+            "fluid_pp(2)%pi_inf": 0.0,
+            "fluid_pp(2)%cv": 717.5,
+        }
+    )
+)

examples/1D_jwl_single_material_shocktube/README.md

@@ -0,0 +1,23 @@
+# 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
+```
+
+Or run all three 1D verification cases:
+
+```bash
+./scripts/run_1d_verification.sh
+```