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

[fahnab666]

JWL EOS PR Description and Capability Tracker

Use this as the copy-paste source for the MFC feature PR. It is written to keep
the five-equation JWL work review-ready while clearly separating validated
capabilities from stress tests and follow-up work.

PR Description

Implements the Jones-Wilkins-Lee (JWL) equation of state in MFC for
detonation-products modeling and five-equation multi-fluid JWL/ideal-gas
mixtures.

This PR adds:

• Per-fluid EOS selection through fluid_pp(i)%eos.
• JWL pressure recovery, energy inversion, and sound-speed support.
• GPU-resident JWL material parameter tables.
• JWL-aware total-energy reconstruction in HLL/HLLC paths.
• A contact-restoration repair for the five-equation HLLC JWL path.
• Validation artifacts for exact-reference, invariant, and stress-test cases.

Closes #TBD.

Type of Change

• New feature
• Bug fix / numerical correctness fix
• Tests or validation artifacts
• Documentation update

Implemented Capabilities

Capability	Status	Evidence / Notes
Single-material JWL EOS	Implemented	Pure-JWL Riemann cases compare against exact/sampled JWL references.
Five-equation JWL/ideal-gas mixtures	Implemented	Well-posed compressed-air interface compares against an exact two-material Riemann solver.
Per-fluid EOS selection	Implemented	fluid_pp(i)%eos = 1 selects stiffened gas; fluid_pp(i)%eos = 2 selects JWL.
JWL pressure recovery	Implemented	Uses the Mie-Gruneisen form with the JWL cold curve and thermal contribution.
JWL energy inversion	Implemented	Recovers specific internal energy from (rho, p) for JWL reconstruction and initialization.
JWL sound speed	Implemented	Uses the Rocflu/Stanley-style Mie-Gruneisen sound-speed relation.
HLL/HLLC JWL energy reconstruction	Implemented	Replaces the ideal-gas total-energy reconstruction when a JWL material is active.
Five-equation HLLC contact repair	Implemented	Initializes alpha_rho_L(:) and alpha_rho_R(:) before JWL energy reconstruction.
GPU parameter availability	Implemented	JWL parameter arrays are available on the device for GPU-resident EOS evaluation.
Six-equation JWL/air	Not validated	Early-time smoke checks exist; longer final-time exact-interface runs remain under-debug.
TNT-air blast	Stress evidence	Useful positivity/symmetry stress test, not a headline quantitative validation yet.

JWL EOS Form

The pressure is modeled as a Mie-Gruneisen EOS with a two-exponential JWL cold
curve:

p = p_cold(rho) + omega*(rho/rho0)*e

where

p_cold(rho) =
    A*(1 - omega/(R1*V))*exp(-R1*V)
  + B*(1 - omega/(R2*V))*exp(-R2*V)

V = rho0/rho

The sound speed follows the Mie-Gruneisen/Rocflu-style form used by the
reference solver:

c^2 = dp_cold/drho + (omega/rho0)*e + (omega/rho0)*(p/rho)

This avoids explicit temperature evaluation in the hot EOS path.

EOS Selection

Each material selects its EOS with fluid_pp(i)%eos.

Value	EOS	Required parameters
1	Stiffened gas	Existing stiffened-gas material parameters
2	JWL	fluid_pp(i)%jwl_A, jwl_B, jwl_R1, jwl_R2, jwl_omega, jwl_rho0, and mixture parameters when applicable

Below-cold-curve states require careful interpretation. If a requested state
has p < p_cold(rho), exact round-trip identity is not generally expected once
the implementation applies the non-negative-energy or temperature floor. These
states should be described as cold-floor-sensitive rather than as ordinary
positive-temperature JWL states.

Validation Evidence Tracker

Tier	Test family	Reference / metric	Current result	PR verdict
Core	EOS formula and round-trip checks	(rho, p) -> e -> p' residuals	Formula path checked against reference behavior	Pass / core evidence
Core	Pure 1D JWL Riemann	Exact/sampled JWL Riemann solution	MFC profiles overlap reference with expected shock/contact smearing	Pass / core evidence
Core	Five-equation JWL/air exact interface	Exact two-material Riemann solver	Monotone L1 convergence on the well-posed compressed-air case	Pass / core evidence
Core	Smooth JWL entropy wave	Periodic return / consistency check	Useful smooth-code verification with documented order behavior	Pass / core evidence
Core	2D quasi-1D JWL	1D exact slice plus y-invariance	Zero y-deviation; 2D reproduces 1D behavior	Pass / core evidence
Supporting	JWL mixture closure formulas	Analytic closure residuals	Near machine-precision residuals in closure checks	Supporting evidence
Supporting	Shyue water-air	Exact stiffened-gas two-material regression	All checks pass; validates multi-fluid infrastructure, not JWL directly	Supporting evidence
Stress	2D TNT-air blast	Positivity, symmetry, radial profile	Stable circular blast morphology; no direct spherical reference yet	Stress only
Stress	Ambient-pressure JWL/air mixing	Known extreme pressure-ratio behavior	Useful limitation study	Do not use as headline validation
Exploratory	Kamm/Tarver/PBX cases	Literature-style comparisons	Needs documented metrics and pass/fail tables	Follow-up
Under-debug	Six-equation JWL/air	Exact two-material final-time comparison	Longer runs show air-pressure collapse and grid-independent error	Not validated

Reviewer-Facing Summary

The reviewable feature is the five-equation JWL foundation: single-material JWL,
JWL/ideal-gas five-equation mixtures, EOS conversion routines, and JWL-aware
HLL/HLLC reconstruction. The strongest validation stack is exact-reference JWL
Riemann testing, exact two-material JWL/air convergence, smooth EOS
code-verification, and 2D quasi-1D invariance.

The six-equation JWL/air model is intentionally not claimed as complete in this
PR. Existing six-equation results are useful diagnostics, but longer final-time
comparisons still fail and should remain future work.

PR Checklist

• Added JWL EOS capability.
• Added or updated tests/validation artifacts for new behavior.
• Documented EOS selection and JWL material requirements.
• Separated exact-reference validation from qualitative stress tests.
• Replace #TBD with the issue number before opening the PR.
• Confirm whether the orphan jwl_contact_blend parameter should be removed before review.
• Confirm CPU validation numbers from the final branch after source cleanup.
• Confirm GPU results match CPU results on NVIDIA or AMD hardware.
• Trigger AI/code review after the PR is opened:
  • @coderabbitai review for incremental review.
  • @coderabbitai full review for full review from scratch.
  • /review for Qodo review.
  • /improve for Qodo suggestions.
  • @claude full review or label claude-full-review for Claude review.

Known Limitations and Follow-Up

• Remove or consume the jwl_contact_blend parameter if it remains unused in
  solver logic.
• Add explicit pass/fail reports before promoting Kamm/Tarver/PBX cases from
  exploratory/supporting evidence to validation evidence.
• Run an axisymmetric or 3D spherical blast case before claiming direct
  Kingery-Bulmash blast validation.
• Continue six-equation JWL/air debugging separately from this five-equation
  PR foundation.

[sbryngelson]

Maintainer review — full issue list

Thanks for this. The physics is genuinely thoughtful — four mixture closures, exact inverses for the linear ones, fail-fast parameter validation in init, and an honest validation tracker that separates what's verified from what isn't. That part I'm not asking you to change.

The engineering, however, isn't yet at MFC standards. Below is the complete list of what needs to be addressed before this can merge, grouped by severity, with the reason each one matters. Numbering is for ease of reference in replies.

---

Blocking — correctness & maintainability

1. Six near-identical copies of the JWL energy-reconstruction block.
The same ~14-line block (Y_L/Y_R clamp → alpha_jwl_* → s_jwl_mix_energy_pr → E_* = rho*e_jwl + 0.5*rho*vel_rms) appears 4× in src/simulation/m_riemann_solver_hllc.fpp (~lines 245, 632, 908, 1364) and again in m_riemann_solver_hll.fpp and m_riemann_solver_lf.fpp.
Why it matters: any future fix or sign correction will land in one copy and silently miss the other five. This is the single highest-risk maintenance hazard in the PR. Collapse it into one #:def fpp macro or a shared subroutine and call it from all six sites.

2. HLLC and HLL/LF feed the sound-speed routine different inputs for the same physics.
HLLC passes the bounded primitives jwl_Y/jwl_alpha into s_compute_speed_of_sound; HLL, LF, m_data_output.fpp, and m_time_steppers.fpp instead pass alpha_rho_j = q_prim_vf(max(1, jwl_idx))%sf(...).
Why it matters: the mass fraction driving the sound speed comes from a different source depending on which Riemann solver is selected, so the same physical state produces different wave speeds. There is no physical reason for this divergence — it reads like two iterations of the same idea that never got unified. Pick one input path and use it everywhere.

3. max(1, jwl_idx) indexing hack.
Appears in m_riemann_solver_hll.fpp, m_riemann_solver_lf.fpp, m_data_output.fpp, m_time_steppers.fpp:

alpha_rho_j = q_prim_vf(max(1, jwl_idx))%sf(j, k, l)

Why it matters: when no JWL fluid is present (jwl_idx == 0) this silently reads fluid 1 and passes a meaningless value into the sound-speed call on every non-JWL run, relying on present()/NaN guards in the callee to ignore it. That's a band-aid for not guarding at the call site. Guard the call site instead so non-JWL paths never construct the argument.

4. Runtime branch added to the hottest kernel in every build.
Every flux loop in all three Riemann solvers now evaluates if (jwl_idx > 0 .and. model_eqns == model_eqns_5eq). jwl_idx is a runtime device scalar, not a compile-time constant.
Why it matters: this is contrary to MFC's compile-time case-specialization philosophy — the branch is not optimized away, so every simulation (including the 99% that never touch JWL) pays a predicated branch per face in the hottest part of the code. Gate the JWL path behind an fpp/case-optimization flag so non-JWL builds compile it out entirely.

5. Per-face 60-iteration bisection with transcendentals, on GPU.
s_jwl_ptequil_pressure_er / s_jwl_ptequil_energy_pr (m_jwl.fpp ~321, ~376) run a fixed 60-iteration bisection with exp() calls inside, invoked from the energy-reconstruction path that executes on every Riemann face. The convergence test rarely triggers early, so it's effectively always 60 iterations.
Why it matters: on GPU this is a serial 60-iteration transcendental loop per thread with heavy warp divergence — a severe throughput problem for an exascale target. Either provide a convergence/performance benchmark justifying it, or gate jwl_mix_type = 2/3 as host-only/experimental until it's optimized.

6. Silent NaN/floor squashing masks divergence.
Both s_jwl_mix_pressure_er and s_jwl_mix_energy_pr end with:

if (pres /= pres) pres = sgm_eps
if (pres < sgm_eps) pres = sgm_eps

and s_compute_speed_of_sound (m_variables_conversion.fpp ~307) has:

if (c2 /= c2) c2 = jwl_omegas(jwl_idx)*max(pres, 1._wp)/max(rho, sgm_eps)
c2 = max(c2, jwl_omegas(jwl_idx)*max(pres, 1._wp)/max(rho, sgm_eps))

Why it matters: a diverged bisection or a NaN pressure gets converted to a 1e-16 floor and the solver keeps running with garbage instead of aborting — a classic silent-failure mode that makes bugs nearly impossible to localize. Additionally max(pres, 1._wp) hardcodes a 1 Pa dimensional floor, and the unconditional max(c2, ...) raises every sound speed to at least the products-gas value, biasing wave speeds in air-dominated cells. Replace silent squashing with debug-mode assertions; remove the hardcoded 1._wp.

---

Blocking — validation & test integration

7. The unit test exercises a Python re-implementation, not the compiled Fortran.
toolchain/mfc/test/test_jwl_inverse.py (665 lines) is a standalone Python translation of the EOS routines.
Why it matters: it can pass while the Fortran is wrong (or vice versa) — it cannot catch a Fortran/Python divergence, which is exactly the kind of bug it appears to be guarding against. It provides false confidence. Either drive the compiled Fortran, or clearly relabel it as a reference oracle (and wire it into CI discovery so it actually runs).

8. Golden cases not wired into the test registry, generated on a dirty tree.
The four golden case dirs (tests/25B9827F, 52A7542F, DB95E4F7, E3460991) appear with no corresponding change to the test-case generator in the diff, and their golden-metadata.txt records generation from /Users/fahadnabid/.../MFC on jwl-upstream-rebase (dirty).
Why it matters: if they aren't produced by a registered case spec, the suite will never exercise them (or will flag them stale), and goldens generated from an uncommitted working tree are not reproducible. Wire the cases into the generator and regenerate goldens on a clean checkout.

9. Promote validation tiers only with reports. The tracker lists Kamm/Tarver/PBX as exploratory and TNT-air as stress-only.
Why it matters: fine to keep them out of scope, but the PR should not present any of them as validation. Please keep the headline claim scoped to the exact-reference five-equation stack, which is the part that's actually backed.

---

Non-blocking — quality & style

10. Hardcoded air defaults as non-sentinel values, triple-duplicated.
In src/simulation/m_global_parameters.fpp, src/pre_process/m_global_parameters.fpp, src/post_process/m_global_parameters.fpp:

fluid_pp(i)%jwl_air_e0    = 2.5575e5_wp
fluid_pp(i)%jwl_air_rho0  = 1._wp
fluid_pp(i)%jwl_air_gamma = 0.4_wp

Why it matters: every neighboring field defaults to dflt_real so a forgotten input fails fast; these instead silently inject a specific air model. Make them dflt_real with validation (or document them as deliberate), and factor the block into one shared init helper instead of three copies that will drift.

11. Dead/over-broad public API.
s_jwl_pcold, s_jwl_sound_speed_squared, s_jwl_energy_pr, s_jwl_kuhl_pressure_er, s_jwl_kuhl_energy_pr, s_jwl_kuhl_sound_speed_squared are in the module's public list and re-exported again through m_variables_conversion, but are only ever reached internally via the dispatchers.
Why it matters: gratuitous API surface invites accidental external coupling and obscures the real interface. Make them private; only the s_jwl_mix_*, the two sound-speed dispatchers, jwl_idx, and init/finalize need to be public, and they shouldn't be re-exported a second time.

12. ;-packed multi-statement lines.
The p-T-equilibrium routines pack 3–4 statements per line, e.g.:

rj = max(Y_s*rho_s/a_lo, sgm_eps); ra = max((1._wp - Y_s)*rho_s/(1._wp - a_lo), sgm_eps); V = rho0/rj

Why it matters: MFC is one-statement-per-line throughout; this hurts readability and makes git blame/diffs useless. Unpack them.

13. Unnamed magic thresholds.
Guards mix sgm_eps (1e-16), 1.e-12_wp, and a 1.e-4_wp pure-fluid cutoff in the Rocflu routines (s_jwl_rocflu_pressure_er / _energy_pr / _sound_speed_squared).
Why it matters: the 1.e-4_wp is a physical threshold with no name and no rationale in-code. Give it a named constant and a one-line justification.

14. Possibly-negative denominator flooring in the Kuhl closure.
s_jwl_kuhl_* divides by R_mix via max(rho_safe*R_mix, sgm_eps), but R_mix can be negative depending on the air_gamma/omega combination.
Why it matters: max(., sgm_eps) on a negative value clamps it to +1e-16, flipping the sign of the denominator and silently producing wrong (not just floored) results. Check R_mix > 0 explicitly or document why it can't be negative.

15. Energy floor ignores the cold-curve offset.
The inverse routines floor specific internal energy at max(e, 0._wp).
Why it matters: for JWL the physically meaningful floor is the cold-curve energy, not 0, so the round-trip pressure_er(energy_pr(p)) == p breaks near the floor. The Python unit test only samples e ≥ 0 states, so it won't catch this. At minimum document it as a deliberate safety net.

16. LLM-style verbosity vs MFC norms.
m_jwl.fpp carries multi-sentence justification comments inline (e.g. the Wood's-rule block, the reader-addressed "present() guards must be separate if-blocks because Fortran does not short-circuit…") and defensive boilerplate at the end of every dispatcher.
Why it matters: it's far more heavily commented than the rest of the codebase, which makes it read as out-of-place and harder to skim. Trim to MFC's comment density.

17. Leftover dev cruft in .typos.toml.
The PR adds exclusions for _jwl_verification/ and _jwl_organization/, directories that aren't part of the PR.
Why it matters: dangling references to scratch dirs are confusing and will rot. Remove them.

---

Summary

The closures are the strong part — leave them. The asks, in priority order: de-dup the six energy blocks (1), unify the sound-speed input path (2, 3), gate the JWL branch at compile time (4), justify-or-host-gate the bisection (5), replace silent floors with assertions (6), and fix the test/golden story (7, 8). The rest are cleanup. Happy to discuss any of these.

[danieljvickers]

You should change the level of this loop back a layer to make it consistent with the surrounding code. Make the previous line an else if and leave the previous line at the same level as the previous loop.

[danieljvickers]

Same issue here. Don't add unnecessary if-tree nesting if you can prevent it

[danieljvickers]

I unsure why you created the file test_jwl_inverse.py, but this is a significant deviation from how MFC handles testing in general. If you want to add a unit testing framework, we will need a significant PR just for the testing framework itself. It would be unsustainable to add the framework as you have implemented it here, as it would bloat the code significantly.

Also, for your tests added, you add 3 1D cases and a 2D case. It is not inherently bad to add them, but a feature like this really requires a 3D example to be added with comparison to some sort of example. I know that I have seen comparisons of your code to published examples, and those comparisons should make it into the readme if you are going to provide a readme.

The only concern I have about the fortran code outside of the minor comments that I left are with the size of the jwl module. I need to actually take a look at it in an IDE where I have access to grep and other tools, but it looks like some redundant code that we should consider parsing down.