JWL equation of state: core EOS + linear mixture closure + non-reactive IBM/EL coupling#1659
Draft
fahnab666 wants to merge 3 commits into
Draft
JWL equation of state: core EOS + linear mixture closure + non-reactive IBM/EL coupling#1659fahnab666 wants to merge 3 commits into
fahnab666 wants to merge 3 commits into
Conversation
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.
Contributor
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. |
Comment on lines
349
to
+350
| (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}"), |
Comment on lines
+771
to
+789
| 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") |
Comment on lines
+110
to
+114
| #: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 |
Comment on lines
+4
to
+7
| # 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. |
Comment on lines
+72
to
+73
| # Fluid 1: TNT JWL products (weighted-composition closure). | ||
| "fluid_pp(1)%eos": 2, |
Comment on lines
+10
to
+13
| # 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. |
Comment on lines
+93
to
+94
| # Fluid 1: TNT JWL products (weighted-composition closure). | ||
| "fluid_pp(1)%eos": 2, |
fahnab666
marked this pull request as draft
July 19, 2026 06:49
Member
|
Closing your existing pr #1586 until this one is merged (you can always reopen it from the existing branch on your fork) |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
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).