JWL equation of state: core EOS + linear mixture closure + non-reactive IBM/EL coupling#1659
JWL equation of state: core EOS + linear mixture closure + non-reactive IBM/EL coupling#1659fahnab666 wants to merge 4 commits into
Conversation
There was a problem hiding this comment.
Pull request overview
Adds Tier-1 support for a JWL detonation-products equation of state (EOS) to MFC’s five-equation model, including a two-material JWL/ambient mixture closure and plumbing through key solver/diagnostic/coupling paths (Riemann solvers, immersed boundaries, and Euler–Lagrange bubbles). This also extends the Python toolchain parameter system and case validation, and introduces new examples + golden regressions.
Changes:
- Introduces a new JWL EOS/closure module and a per-fluid
fluid_pp(i)%eosselector; wires JWL into conservative↔primitive conversion, Riemann solvers, and sound-speed diagnostics. - Extends non-reactive coupling paths (IBM ghost/fresh-cell rebuild and EL bubble seed pressure / Keller–Miksis sound speed) to use the JWL closure when active.
- Adds toolchain params/validation/help text plus documentation, examples, and new golden tests for JWL scenarios.
Reviewed changes
Copilot reviewed 51 out of 52 changed files in this pull request and generated 8 comments.
Show a summary per file
| File | Description |
|---|---|
| toolchain/mfc/test/cases.py | Adds kernel/golden test cases covering JWL closure and JWL+IBM paths. |
| toolchain/mfc/params/descriptions.py | Adds parameter descriptions for new JWL inputs (but currently missing one for fluid_pp%eos). |
| toolchain/mfc/params/definitions.py | Registers new fluid_pp%eos and JWL parameter definitions in the params source of truth. |
| toolchain/mfc/jwl_products.py | Provides reusable, named JWL products parameter sets for case authors. |
| toolchain/mfc/case.py | Exposes jwl_active to Fypp via case preprocessing for conditional compilation. |
| toolchain/mfc/case_validator.py | Validates EOS/JWL inputs and prohibits incompatible feature combos for eos_jwl. |
| tests/F0754EFF/golden-metadata.txt | Adds golden generation metadata for a new/updated regression. |
| tests/E3460991/golden-metadata.txt | Adds golden generation metadata for new/updated regressions. |
| tests/DC926C1D/golden-metadata.txt | Adds golden generation metadata for new/updated regressions. |
| tests/DB95E4F7/golden-metadata.txt | Adds golden generation metadata for new/updated regressions. |
| tests/C1099A06/golden-metadata.txt | Adds golden generation metadata for new/updated regressions. |
| tests/5F515D10/golden.txt | Adds/updates golden reference outputs for a JWL regression. |
| tests/5F515D10/golden-metadata.txt | Adds golden generation metadata for the 5F515D10 regression. |
| tests/1F54C11E/golden-metadata.txt | Adds golden generation metadata for a new/updated regression. |
| src/simulation/m_riemann_solver_lf.fpp | Hooks JWL closure into Lax–Friedrichs face energy reconstruction and sound speed. |
| src/simulation/m_riemann_solver_hllc.fpp | Hooks JWL closure into HLLC face reconstruction and c/c_avg evaluation paths. |
| src/simulation/m_riemann_solver_hll.fpp | Hooks JWL closure into HLL face energy reconstruction and c/c_avg evaluation paths. |
| src/simulation/m_ibm.fpp | Rebuilds IBM ghost/fresh-cell energy via JWL closure when active; adds marker snapshot for moving bodies. |
| src/simulation/m_global_parameters.fpp | Adds JWL parameter storage arrays and assigns new fluid defaults. |
| src/simulation/m_checker.fpp | Prohibits incompatible runtime feature combinations when a JWL fluid is present. |
| src/simulation/m_bubbles_EL.fpp | Uses JWL closure for bubble seed pressure and liquid sound speed in JWL host cells; rejects seeding inside IB. |
| src/simulation/include/inline_riemann.fpp | Adds an inline macro to reconstruct JWL face energy + sound speed in one fused call. |
| src/pre_process/m_global_parameters.fpp | Initializes new EOS selector and JWL defaults in pre_process. |
| src/pre_process/m_data_output.fpp | Ensures pre_process diagnostic pressure recovery passes JWL products mass fraction when applicable. |
| src/post_process/m_global_parameters.fpp | Initializes new EOS selector and JWL defaults in post_process. |
| src/post_process/m_derived_variables.fpp | Computes sound speed via JWL closure in post_process when JWL is active. |
| src/common/m_variables_conversion.fpp | Adds JWL-aware pressure/energy conversion and a JWL sound-speed helper; initializes/finalizes JWL module. |
| src/common/m_jwl.fpp | New JWL EOS + linear mass-fraction closure implementation with init-time verification scan. |
| src/common/m_global_parameters_common.fpp | Adds common default assignment for new per-fluid JWL parameters across all targets. |
| src/common/m_derived_types.fpp | Extends physical_parameters with EOS selector and JWL parameter fields. |
| src/common/m_constants.fpp | Adds EOS selector constants (eos_stiffened_gas, eos_jwl). |
| src/common/m_checker_common.fpp | Adds common input validation for EOS selector + required JWL parameters. |
| README-JWL-EOS.md | New top-level documentation describing the Tier-1 JWL EOS/closure and case setup. |
| examples/3D_jwl_air_burst/case.py | New 3D inert air-burst example exercising JWL closure in 3D. |
| examples/2D_jwl_ibm_static_cylinder/case.py | New 2D example exercising JWL + static IBM coupling. |
| examples/2D_jwl_ibm_freeflight_cylinder/case.py | New 2D example exercising JWL + moving IBM fresh-cell repopulation. |
| examples/2D_jwl_el_bubbles/input/lag_bubbles.dat | Input bubble seed file for the EL bubbles example. |
| examples/2D_jwl_el_bubbles/case.py | New 2D example exercising JWL + Euler–Lagrange bubble coupling. |
| examples/2D_jwl_detonation/case.py | New 2D inert “detonation products burst” example case. |
| examples/1D_jwl_single_material_shocktube/README.md | New 1D single-material JWL shocktube documentation/verification notes. |
| examples/1D_jwl_single_material_shocktube/case.py | New 1D single-material JWL shocktube example case. |
| examples/1D_jwl_mixture_test/case.py | New 1D JWL/air mixing shocktube example (comments currently use inconsistent closure terminology). |
| examples/1D_jwl_air_interface_advection/case.py | New 1D mechanical-equilibrium interface guard case (comments currently use inconsistent closure terminology). |
| examples/1D_jwl_air_interface_advection/analyze.py | New analysis script for interface-advection metrics (docstring currently outdated w.r.t. t=0 diagnostics). |
| docs/module_categories.json | Registers m_jwl under Physics Models for documentation grouping. |
| docs/documentation/case.md | Documents fluid_pp%eos and the Tier-1 JWL EOS/closure setup and constraints. |
| (r"fluid_pp\((\d+)\)%mu_bulk", "Bulk viscosity (non-Newtonian) for fluid {0}"), | ||
| (r"fluid_pp\((\d+)\)%jwl_A", "JWL A coefficient for fluid {0}"), |
| self.prohibit(num_fluids > 1 and num_fluids - len(jwl_ids) != 1, "the weighted-composition closure requires exactly one non-JWL ambient fluid") | ||
| ji = jwl_ids[0] | ||
| cv_jwl = self.get(f"fluid_pp({ji})%cv") | ||
| self.prohibit(cv_jwl is None or cv_jwl <= 0, f"the weighted-composition closure requires positive fluid_pp({ji})%cv for the JWL fluid") | ||
| air_ids = [i for i in range(1, num_fluids + 1) if self.get(f"fluid_pp({i})%eos") != 2] | ||
| if air_ids: | ||
| cv_air = self.get(f"fluid_pp({air_ids[0]})%cv") | ||
| self.prohibit(cv_air is None or cv_air <= 0, f"the weighted-composition closure requires positive fluid_pp({air_ids[0]})%cv for the non-JWL ambient fluid") | ||
| rho0 = self.get(f"fluid_pp({ji})%jwl_rho0") | ||
| air_rho0 = self.get(f"fluid_pp({ji})%jwl_air_rho0") | ||
| if rho0 is not None and air_rho0 is not None: | ||
| self.prohibit(rho0 <= air_rho0, "the weighted-composition closure requires products reference density above the ambient-gas density") | ||
| E0 = self.get(f"fluid_pp({ji})%jwl_E0") | ||
| Q = self.get(f"fluid_pp({ji})%jwl_Q") | ||
| if E0 is None and Q is not None and rho0 is not None: | ||
| E0 = rho0 * Q | ||
| air_e0 = self.get(f"fluid_pp({ji})%jwl_air_e0") | ||
| if E0 is not None and rho0 is not None and air_e0 is not None: | ||
| self.prohibit(E0 / rho0 <= air_e0, "the weighted-composition closure requires products reference energy above the ambient-gas energy") |
| #:if not MFC_CASE_OPTIMIZATION or jwl_active | ||
| if ((jwl_idx > 0) .and. (.not. mhd) .and. (model_eqns /= model_eqns_4eq) .and. (bubbles_euler .neqv. .true.)) then | ||
| Y_jwl = 1._wp | ||
| if (present(jwl_Y)) Y_jwl = jwl_Y | ||
| eint = energy - dyn_p |
| # 1D air + JWL mixing shock tube exercising the weighted-composition closure. | ||
| # Fluid 1 = TNT JWL products (eos=2); fluid 2 = ambient air (ideal gas, eos=1). | ||
| # A high-pressure products slug (0 <= x <= 0.3 m) drives a shock into ambient air, | ||
| # so the mixed band 0 < Y < 1 forms at the contact and rides the weighted-composition closure. |
| # Fluid 1: TNT JWL products (weighted-composition closure). | ||
| "fluid_pp(1)%eos": 2, |
| # spurious pressure/velocity oscillation generated by the weighted-composition JWL/air mixing | ||
| # closure in the WENO-reconstructed mixed band (0 < Y < 1) at each contact. The | ||
| # weighted-composition closure preserves this state to machine precision (max|p-p0|/p0 ~ 1e-16), | ||
| # so this case is a protective regression guard on mechanical equilibrium. |
| # Fluid 1: TNT JWL products (weighted-composition closure). | ||
| "fluid_pp(1)%eos": 2, |
|
Closing your existing pr #1586 until this one is merged (you can always reopen it from the existing branch on your fork) |
Add a Jones-Wilkins-Lee EOS for detonation products to the five-equation model (model_eqns=2) with a two-material mixture rule that blends products and ambient linearly in the products mass fraction Y. Coefficients depend on Y alone, so pressure stays linear in internal energy: the (rho,p,Y) inverse is closed form and the sound speed is the exact Grueneisen derivative. The blend is endpoint-exact (Y=1 pure JWL, Y=0 ambient) and recovers the ideal closure bit-identically when the ambient pi_inf is zero. Core closure in src/common/m_jwl.fpp is GPU-resident and closed-form, with an init-time envelope scan that aborts on any non-positive sound speed or failed round trip. EOS dispatch is wired through variables conversion, the inline Riemann energy reconstruction, and the HLLC/HLL/LF solvers. Parameters, checkers (sim and toolchain), presets, four basic examples, and their goldens are included. Reactions, the composition-weighted closure, and delta_e ZND structure are deferred to later tiers.
Rebuild the immersed-boundary ghost-cell energy through the JWL closure so near-body pressure and loads match the interior EOS rather than the stiffened-gas relation, and repopulate moving-IB fresh (solid to fluid) cells by inverse-distance extrapolation closed through the JWL energy inverse instead of leaving them at the near-vacuum p=1 Pa solid seed. Both paths are gated on jwl_idx, so base immersed-boundary behavior is bit-identical; the pre-move marker snapshot is allocated only for JWL moving-IB cases. Adds a static-cylinder and a free-flight-cylinder example plus an IBM golden test.
…ample Seed Lagrangian bubble pressure and the Keller-Miksis liquid sound speed through the JWL closure so a bubble placed in a products or mixture cell is initialised and driven with the JWL pressure law rather than a stiffened-gas relation. Both paths are gated on jwl_idx, so base Lagrange-bubble behavior is bit-identical. A bubble seeded inside an immersed body is rejected at input time. Adds a 2D products-blast bubble-cloud example (needs input/lag_bubbles.dat, so example-only) and a 3D inert spherical air-burst example sized for a roughly five-minute eight-core run, the exascale-structure smoke test for the closure and the Riemann/cons-prim path in 3D.
There was a problem hiding this comment.
i'm not adding this file to the repo. i'll come up w/ some other strategy
|
Thanks for the work here. After reviewing, I'm closing this PR in favor of introducing a general EOS interface first. here's a plan that i like:
Much of this work should be reusable.. but we need the general thermochem architecture in place first [see Slack] |
Summary
Works toward #1638. This is a smaller, review-friendly slice of #1586 (which lands the
full detonation stack in one about 10k-line PR). The capability here is one JWL products fluid
plus one ambient gas, the simplest possible mixture rule, and non-reactive immersed-boundary
and Euler-Lagrange bubble coupling. The diff against
masteris 3.8k lines.Why a second PR
#1586 delivers everything at once: the JWL EOS, a composition (heat-capacity) weighted
closure, a reactant/product energy offset, three reaction sources (program burn, afterburn,
JWL++ reactive), and IBM + Euler-Lagrange coupling. That is a lot to review in a single pass.
This PR is the first tier: the pieces can merge and trust on their own, split so each
commit is independently reviewable. Everything reactive and the composition-weighted closure
become follow-up PRs layered on this base.
Core of this PR
model_eqns = 2), selected per fluid withfluid_pp(i)%eos = 2, coexisting with an ideal-gas or stiffened-gas ambient.coefficient is blended linearly in the products mass fraction
Y, so pressure stays linearin internal energy. The
(rho, p, Y) -> einverse is closed form, the sound speed is theexact Grueneisen derivative, and the closure is endpoint-exact (pure JWL at
Y = 1, ambientat
Y = 0) and bit-identical to the ideal closure when the ambientpi_infis zero. Aninit-time envelope scan aborts on any non-physical state.
energy rebuild and the bubble seed-pressure / Keller-Miksis sound-speed paths run through the
JWL closure, all gated on
jwl_idxso base behavior is bit-identical.The base-without-JWL path is bit-identical throughout; incompatible features
(
wave_speeds = 2, CBC boundaries,alt_soundspeed, elasticity,igr,bubbles_euler,mhd,chemistry) are rejected at startup and in./mfc.sh validate.Commits
Deferred to follow-up PRs
Composition (heat-capacity) weighted closure, reactant/product energy offset (ZND structure),
the three reaction sources, and reactive-mode IBM/EL coupling.
Testing
./mfc.sh format,precheck(all checks), andbuild(all three targets) are clean; theJWL golden tests pass; the 2D/3D examples run NaN-free (the 3D air-burst in ~5 min on 8 cores).