JWL equation of state: core EOS + linear mixture closure + non-reactive IBM/EL coupling

README-JWL-EOS.md

@@ -0,0 +1,76 @@
+# JWL equation of state (simple form)
+
+This build adds a Jones-Wilkins-Lee (JWL) equation of state for detonation
+products to MFC's five-equation model (`model_eqns = 2`), together with a simple
+two-material mixture closure. It is the first tier of a staged contribution: one
+JWL products fluid plus one ambient gas, a linear mass-fraction mixture rule, and
+no reactions. The reactive and composition-weighted layers are listed under
+Roadmap below and land as later pull requests.
+
+## Model
+
+Pure JWL products follow the Mie-Grueneisen form referenced to an isentrope,
+with `V = rho0/rho`:
+
+```
+p = A (1 - w/(R1 V)) exp(-R1 V) + B (1 - w/(R2 V)) exp(-R2 V) + w rho e
+```
+
+`A`, `B` (Pa), `R1`, `R2`, and `w = omega` are cylinder-test fits. A cell holding
+a products mass fraction `Y = alpha_rho_products / rho` mixes products and ambient
+by blending the coefficients linearly in `Y`:
+
+```
+An = Y A,   Bn = Y B,   omega = air_gamma + Y (omega0 - air_gamma),
+cv  = Y cv_products + (1 - Y) cv_air
+```
+
+The blend is exact at the endpoints: `Y = 1` recovers the pure JWL law and
+`Y = 0` recovers the ambient law. A stiffened-gas ambient (for example water)
+adds `pi_hat = (1 - Y) pi_inf`; an ideal-gas ambient (`pi_inf = 0`) recovers the
+ideal closure bit-identically. Because every coefficient depends on `Y` alone and
+never on `rho` or `e`, the pressure stays linear in `e`, so the `(rho, p, Y) -> e`
+inverse is closed form and the sound speed is the exact Grueneisen derivative.
+At start-up `s_jwl_verify_closure` sweeps the `(rho, e, Y)` envelope and aborts on
+any non-positive sound speed or failed pressure/energy round trip. The full
+derivation lives in `src/common/m_jwl.fpp`.
+
+## Case setup
+
+Select the EOS per fluid with `fluid_pp(i)%eos`: `1` stiffened gas (default), `2`
+JWL. At most one fluid may be JWL, and it requires `model_eqns = 2`. A JWL fluid
+is defined by:
+
+| Parameter | Meaning | Requirement |
+| :--- | :--- | :--- |
+| `jwl_A`, `jwl_B`, `jwl_R1`, `jwl_R2`, `jwl_omega` | JWL products coefficients | required |
+| `jwl_rho0` | products reference density | required |
+| `jwl_Q` or `jwl_E0` | detonation energy, specific (J/kg) or volumetric (J/m^3); given `jwl_Q`, MFC sets `jwl_E0 = jwl_rho0 * jwl_Q` | one of the two |
+| `jwl_air_rho0` | density of the co-existing ambient gas | required |
+| `jwl_air_e0` or `jwl_air_p0` | ambient specific internal energy or pressure | one of the two |
+| `cv` | positive heat capacity, on both the JWL fluid and the ambient fluid | required |
+
+The ambient Grueneisen coefficient `omega0_air = 1/fluid_pp(air)%gamma` comes from
+the non-JWL fluid (or the JWL fluid's own `gamma` when it is alone). Recall that
+`fluid_pp%gamma` stores `1/(gamma - 1)`, so air is `2.5`, not `1.4`.
+
+JWL is prohibited with `wave_speeds = 2`, characteristic (CBC) boundaries,
+`alt_soundspeed`, elasticity, `igr`, `bubbles_euler`, `mhd`, and `chemistry`;
+each of those paths evaluates a stiffened-gas pressure that would bypass the JWL
+closure. Immersed boundaries (`ib`) and Lagrangian bubbles (`bubbles_lagrange`)
+are supported.
+
+## Examples
+
+- `examples/1D_jwl_single_material_shocktube` — single JWL fluid, two-pressure shock tube.
+- `examples/1D_jwl_mixture_test` — products and air across the linear closure.
+- `examples/1D_jwl_air_interface_advection` — advected products/air interface.
+- `examples/2D_jwl_detonation` — cylindrical products charge bursting into air.
+
+## Roadmap (Tier 2, later pull requests)
+
+- Composition (heat-capacity) weighted closure `w = Y cv_j / (Y cv_j + (1 - Y) cv_a)`
+  in place of the linear `w = Y` blend.
+- Reactant/product energy offset (`jwl_delta_e`) for resolved ZND detonation structure.
+- Reaction sources: program burn, products-air afterburn, and JWL++ pressure-driven reactive burn.
+- Reactive-mode immersed-boundary and Lagrangian-bubble coupling.

docs/documentation/case.md

@@ -438,6 +438,7 @@ A particle cloud is a compact specification of a bed of identical circular (2D)
 | `qvp`  ** | Real   | Stiffened-gas parameter $q'$ of fluid.         |
 | `sigma`   | Real   | Surface tension coefficient                    |
 | `G`       | Real   | Shear modulus of solid.                        |
+| `eos` *   | Integer| Equation of state: `1` stiffened gas, `2` JWL. |
 
 Fluid material's parameters. All parameters except for sigma should be prepended with `fluid_pp(i)` where $i$ is the fluid index.
 
@@ -459,6 +460,25 @@ 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), or `2` Jones-Wilkins-Lee (JWL) for detonation products. JWL setup is described in [JWL equation of state](#sec-jwl-eos) below.
+
+#### JWL equation of state {#sec-jwl-eos}
+
+The JWL equation of state models detonation products. It is supported only with `model_eqns = 2`, and at most one fluid may set `fluid_pp(i)%%eos = 2`. A JWL fluid is defined by these `fluid_pp(i)%%` parameters:
+
+| Parameter | Meaning | Requirement |
+| :--- | :--- | :--- |
+| `jwl_A`, `jwl_B`, `jwl_R1`, `jwl_R2`, `jwl_omega` | JWL products EOS coefficients | required |
+| `jwl_rho0` | products reference density \f$\rho_0\f$ | required |
+| `jwl_Q` or `jwl_E0` | detonation energy: specific (J/kg) or volumetric (J/m³). Given `jwl_Q`, MFC sets `jwl_E0 = jwl_rho0 * jwl_Q` | one of the two |
+| `jwl_air_rho0` | density of the co-existing ambient gas | required |
+| `jwl_air_e0` or `jwl_air_p0` | ambient specific internal energy or pressure | one of the two |
+| `cv` | positive heat capacity, on both the JWL fluid and the ambient fluid | required |
+
+The ambient Grüneisen coefficient is \f$\gamma - 1\f$, obtained from the ambient gas fluid's own stored `gamma` (which holds \f$1/(\gamma - 1)\f$, so its reciprocal recovers \f$\gamma - 1\f$). With a single JWL fluid and no separate ambient fluid, the JWL fluid's own `gamma` is used.
+
+Products mix with the surrounding gas through a linear mass-fraction closure: every coefficient is blended linearly in the products mass fraction \f$Y\f$, so the pressure stays linear in internal energy and the closure recovers pressure, temperature, and sound speed from \f$(\rho, e, Y)\f$ in closed form. It degenerates exactly to the pure-JWL law at \f$Y = 1\f$ and to the ambient law at \f$Y = 0\f$, and recovers the ideal-gas closure bit-identically when the ambient `pi_inf` is zero. A stiffened-gas ambient (e.g. water) is supported by setting the non-JWL fluid's `pi_inf`. JWL is not supported with `wave_speeds = 2`, characteristic (CBC) boundaries, `alt_soundspeed`, elasticity, `igr`, `bubbles_euler`, `mhd`, or `chemistry`; immersed boundaries (`ib`) and Lagrangian bubbles (`bubbles_lagrange`) are supported. See `src/common/m_jwl.fpp` for the full closure derivation.
+
 > **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

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_air_interface_advection/analyze.py

@@ -0,0 +1,72 @@
+#!/usr/bin/env python3
+"""Report the Abgrall interface-advection metrics for the JWL/air closure.
+
+The case initializes one uniform pressure p0 and one uniform velocity u0 across
+a JWL-products slab embedded in air, so the exact solution holds p == p0 and
+u == u0 for all time while the two material contacts merely translate. This
+script reads the primitive-variable output and reports, for each saved step,
+
+    max|p - p0| / p0     and     max|u - u0|
+
+split into the pure regions (products volume fraction alpha < 0.01 or > 0.99)
+and the WENO-reconstructed mixed interface band (0.01 <= alpha <= 0.99). p0 and
+u0 are the uniform reference, taken from the median of the pure-region field at
+each step (robust to the few transition cells).
+
+Both pre_process (t=0) and the simulation (t>0) recover pressure with the cell's
+own products mass fraction, so every saved step holds p == p0 to machine
+precision, including the t=0 diagnostic.
+"""
+
+import glob
+import os
+import sys
+
+import numpy as np
+
+D = os.path.join(os.path.dirname(os.path.abspath(__file__)), "D")
+
+
+def series(var):
+    files = sorted(glob.glob(os.path.join(D, f"prim.{var}.*.dat")))
+    if not files:
+        sys.exit(f"no prim.{var}.*.dat files found in {D}; run the case first")
+    return files
+
+
+def col2(path):
+    return np.loadtxt(path)[:, 1]
+
+
+def main():
+    pres_files = series(4)
+    vel_files = series(3)
+    alpha_files = series(5)
+
+    print(f"{'step':>10} {'region':>7} {'p0 (Pa)':>14} " f"{'max|dp|/p0':>12} {'max|du| (m/s)':>14}")
+
+    worst_pure_dp = worst_pure_du = 0.0
+    for pf, uf, af in zip(pres_files, vel_files, alpha_files):
+        p, u, a = col2(pf), col2(uf), col2(af)
+        mixed = (a >= 0.01) & (a <= 0.99)
+        pure = ~mixed
+        p0 = float(np.median(p[pure]))
+        u0 = float(np.median(u[pure]))
+        step = os.path.basename(pf).split(".")[3]
+
+        for label, mask in (("pure", pure), ("mixed", mixed)):
+            if not np.any(mask):
+                continue
+            dp = float(np.max(np.abs(p[mask] - p0))) / abs(p0)
+            du = float(np.max(np.abs(u[mask] - u0)))
+            print(f"{step:>10} {label:>7} {p0:>14.6e} {dp:>12.3e} {du:>14.3e}")
+            if label == "pure":
+                worst_pure_dp = max(worst_pure_dp, dp)
+                worst_pure_du = max(worst_pure_du, du)
+
+    print()
+    print(f"pure-region worst: max|dp|/p0 = {worst_pure_dp:.3e}, " f"max|du| = {worst_pure_du:.3e} m/s")
+
+
+if __name__ == "__main__":
+    main()

examples/1D_jwl_air_interface_advection/case.py

@@ -0,0 +1,114 @@
+#!/usr/bin/env python3
+import json
+
+# 1D JWL/air mechanical-equilibrium interface guard (mixture-closure isolation).
+#
+# A JWL-products slab (0.25 <= x <= 0.75) sits in ambient air. The whole domain
+# is initialized at ONE uniform pressure p0 and a QUIESCENT field (u0 = 0), so
+# the exact solution holds p == p0 and u == 0 for all time (the two material
+# contacts are stationary). Any deviation max|p-p0|/p0 or max|u-u0| is a
+# spurious pressure/velocity oscillation generated by the linear mass-fraction JWL/air mixing
+# closure in the WENO-reconstructed mixed band (0 < Y < 1) at each contact. The
+# linear mass-fraction closure preserves this state to machine precision (max|p-p0|/p0 ~ 1e-16),
+# so this case is a protective regression guard on mechanical equilibrium.
+#
+# KNOWN BEHAVIOR (Abgrall moving-contact oscillation): setting u0 != 0 advects
+# the contacts and produces a decaying spurious pressure dip riding each moving
+# interface (~15% at u0=1000). This is the classic Abgrall signature of the 5-eq
+# Allaire model with a general (energy/density-dependent) EOS -- it is a
+# model-level limitation, NOT specific to the JWL closure, and is documented here
+# rather than guarded against. The stationary (u0=0) case is the shippable guard.
+#
+# Fluid 1 = TNT JWL products (eos=2); fluid 2 = ambient air (ideal gas, eos=1).
+# Periodic BCs keep both contacts in the domain.
+u0 = 0.0
+p0 = 1.0e8
+rho_products = 400.0  # expanded products: p_cold(rho) < p0, so specific e > 0
+rho_air = 50.0
+
+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
+
+print(
+    json.dumps(
+        {
+            "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": 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": -1,
+            "bc_x%end": -1,
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            # Patch 1: ambient air over the whole domain at (p0, u0).
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%length_x": 1.0,
+            "patch_icpp(1)%vel(1)": u0,
+            "patch_icpp(1)%pres": p0,
+            "patch_icpp(1)%alpha_rho(1)": 1.0e-8 * rho_products,
+            "patch_icpp(1)%alpha_rho(2)": 0.99999999 * rho_air,
+            "patch_icpp(1)%alpha(1)": 1.0e-8,
+            "patch_icpp(1)%alpha(2)": 0.99999999,
+            # Patch 2: JWL-products slab (0.25 <= x <= 0.75) at the SAME (p0, u0).
+            "patch_icpp(2)%geometry": 1,
+            "patch_icpp(2)%alter_patch(1)": "T",
+            "patch_icpp(2)%x_centroid": 0.5,
+            "patch_icpp(2)%length_x": 0.5,
+            "patch_icpp(2)%vel(1)": u0,
+            "patch_icpp(2)%pres": p0,
+            "patch_icpp(2)%alpha_rho(1)": 0.99999999 * rho_products,
+            "patch_icpp(2)%alpha_rho(2)": 1.0e-8 * rho_air,
+            "patch_icpp(2)%alpha(1)": 0.99999999,
+            "patch_icpp(2)%alpha(2)": 1.0e-8,
+            # Fluid 1: TNT JWL products (linear mass-fraction closure).
+            "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_e0": 2.5575e5,
+            "fluid_pp(1)%jwl_air_rho0": 1.225,
+            # Fluid 2: ambient air (ideal gas).
+            "fluid_pp(2)%eos": 1,
+            "fluid_pp(2)%gamma": 2.5,
+            "fluid_pp(2)%pi_inf": 0.0,
+            "fluid_pp(2)%cv": 717.5,
+        }
+    )
+)