JWL equation of state and detonation modeling for the five-equation model

[fahnab666]

Summary

Closes #1638.

This PR adds a Jones-Wilkins-Lee (JWL) equation of state and a detonation modeling
capability to the five-equation Allaire model (model_eqns = 2). A JWL fluid is selected
per material with fluid_pp(i)%eos = 2 and coexists with an ideal-gas or stiffened-gas
ambient. On top of the EOS, the PR adds three optional reaction sources (program burn,
afterburn, and self-propagating reactive detonation) and couples the whole capability with
immersed boundaries and Euler-Lagrange bubbles.

Mixture closure

For a pure products cell the JWL pressure law is

$$p = A\left(1-\frac{\omega}{R_1V}\right)e^{-R_1V} + B\left(1-\frac{\omega}{R_2V}\right)e^{-R_2V} + \omega\rho e, \qquad V=\frac{\rho_0}{\rho}.$$

Products and the surrounding gas are mixed with a composition (heat-capacity) weighted
closure. With $Y$ the JWL mass fraction, the products' heat-capacity share is

$$ w = \frac{Y,c_{v,\text{prod}}}{Y,c_{v,\text{prod}} + (1-Y),c_{v,\text{amb}}}, $$

and the mixture coefficients are

$$ A_n = w A, \qquad B_n = w B, \qquad \omega = \Gamma_{\text{amb}} + w,(\omega_0 - \Gamma_{\text{amb}}), \qquad \hat{\pi} = (1-w),\pi_\infty , $$

where $\hat{\pi}$ is a cold offset that is active only for a stiffened ambient
($\pi_\infty > 0$). Every coefficient depends on $Y$ alone, never on $\rho$ or $e$, which is
what makes the closure cheap and exact at the endpoints:

• the pressure-to-energy inverse is a single closed-form expression (no iterative solve),
• the sound speed is the exact analytic Grüneisen derivative, and
• the closure degenerates exactly to pure JWL at $Y=1$ and to the ambient law at $Y=0$,
  so no handover smoothing is needed.

An initialization-time scan verifies the $e \to p \to e$ round-trip over the reachable state
envelope and aborts on any non-physical state.

Reaction sources

All three sources default off, and each carries a $(1-x)/\Delta t$ rate clamp so the SSP-RK
stages cannot overshoot the $Y,q$ energy budget.

Source	Model	Notes
prog_burn	Kinematic program burn: a front expands from the detonation point at $D_{cj}$ and deposits $Q$ over a band of width pb_width	Front position is exact by construction
jwl_afterburn	Products-air afterburn with an advected progress variable; mixing-rate or Arrhenius law	Ideal-gas ambient; requires HLLC
jwl_reactive	JWL++ (Souers 2000) pressure-driven burn, self-propagating from a hot spot	Requires HLLC; exclusive with prog_burn

The JWL++ reaction progress evolves as

$$ \frac{d\lambda}{dt} = G,p^{,b},(1-\lambda). $$

For jwl_reactive, the optional fluid_pp(i)%jwl_delta_e ($\le 0$, default $0$ = off)
applies the Garno et al. (2020) reactant energy offset: the thermal term uses

$$ e_{\text{eff}} = e + Y,(1-\lambda),\Delta e , $$

putting unreacted explosive on a stiffer Hugoniot so a resolved detonation shows genuine ZND
structure (a von Neumann spike decaying to the Chapman-Jouguet state) instead of a monotonic
energy source. The $Y$ factor confines the offset to the explosive, and the inverse and sound
speed stay closed-form. The two progress variables ride the HLLC contact flux like the color
function.

Capabilities added

Capability	Detail
JWL EOS	Composition-weighted mixture of JWL products with an ideal-gas or stiffened-gas ambient, model_eqns = 2
Program burn	Kinematic front from a detonation point (pb_* parameters)
Afterburn	Advected-progress products-air energy release (mixing-rate or Arrhenius)
Reactive burn	Pressure-driven JWL++ detonation with optional ZND offset
Immersed boundaries	Static and moving ib bodies rebuild ghost energy and progress scalars through the closure
Lagrangian bubbles	bubbles_lagrange seed pressure and Keller-Miksis sound speed go through the closure
Diagnostics	pre_process and simulation reach the same closure leaf routine for the $t=0$ pressure

Prohibited with eos_jwl (their pressure paths bypass the closure): wave_speeds = 2,
CBC boundaries, alt_soundspeed, elasticity, igr, bubbles_euler, mhd, and chemistry.
Constraints are enforced both at startup (m_checker*) and in ./mfc.sh validate
(case_validator.py).

Validation

• Initialization self-scan and energy round-trip at the floating-point floor.
• Exact Riemann star states for products/air and products/water interfaces.
• ZND detonation against the analytic Hugoniot, Rayleigh, and CJ construction
  (examples/1D_jwl_znd_detonation).
• Sedov-Taylor blast scaling (examples/2D_sedov_taylor_blast) and far-field TNT
  overpressure against Kinney-Graham
  (benchmarks/3D_jwl_spherical_tnt_free_air_validation).
• Golden regression tests cover single-material and mixture shocktubes, interface advection,
  program burn, afterburn, reactive burn, and the immersed-boundary and bubble couplings.

Tasks

• Composition-weighted mixture closure with closed-form inverse and analytic sound speed
• Program burn, afterburn, and JWL++ reactive sources with the $(1-x)/\Delta t$ clamp
• jwl_delta_e reactant offset for resolved ZND structure
• Immersed-boundary coupling (ghost energy and progress scalars through the closure)
• Euler-Lagrange bubble coupling
• Feature-compatibility checks in m_checker* and case_validator.py
• Worked examples, validation benchmarks, and golden regression tests
• docs/documentation/case.md parameters and capabilities overview updated
• ./mfc.sh format, precheck (7/7), and a full three-target build pass
• JWL golden regression suite passes
• GPU-hardware run of the new JWL paths (constructs have in-repo GPU precedent; no NVIDIA hardware available locally)
• Maintainer review

AI review

• Request an automated AI code review (e.g. @copilot review) before human review
• AI review comments addressed

Notes for reviewers

The GPU device paths use the GPU_* Fypp macros only, and the new sources and progress-variable
advection are feature-gated and GPU-wrapped. The Mie-Grüneisen $N$-constituent
pressure-equilibrium closure (for products, water, and air together) is the sanctioned next
step and is described in README-JWL-IMPLEMENTATION.md; it is not part of this PR.

[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]

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.

[sbryngelson]

Automated re-review (run on behalf of @sbryngelson)

🤖 This is an automated re-review against the current head. It re-checks the 17 items from the prior maintainer review by reading the actual diff (not commit messages) and looks for regressions from the revision churn (22 commits, two master merges). CI is fully green across all NVHPC CPU/GPU compilers, Intel, GNU, Convergence, and Verrou FP-stability.

Nice work on the revision — most of the prior review is genuinely addressed in the code, not just the commit log. The energy-reconstruction macro de-duplicates correctly, the sound-speed inputs are unified, the indexing hack and Python re-implementation test are gone, the defaults are sentinel-ized through a shared helper, the public API is tight, and eos_idxs is removed with no dangling references. The qv double-count fix and the Rocflu omega-term fix are present and look correct.

The following items are still open and should each be resolved before merge.

Resolved

For the record, these are confirmed fixed in the current diff: de-duplicated JWL_RECONSTRUCT_ENERGY macro (#1); unified (jwl_Y, jwl_alpha) sound-speed inputs across hllc/hll/lf (#2); removal of max(1, jwl_idx) (#3); deleted test_jwl_inverse.py (#7); dflt_real air defaults via s_assign_jwl_fluid_defaults (#10); narrowed public list with no re-export (#11); named jwl_pure_cutoff constant (#13); documented energy floor (#15); .typos.toml cleanup (#17); plus the validation-tier wording (#9) and most of the style pass (#12, #16).

Still open

A. Golden-case provenance and registration.
tests/25B9827F and tests/E3460991 carry golden-metadata.txt recording generation on a dirty tree (on jwl-upstream-rebase (dirty)), the personal branch name, and four hardcoded absolute paths (/Users/fahadnabid/Documents/default/MFC/...). There is also no change to the test-case registry (tests/cases.py) in the diff, so it cannot be confirmed that these cases are actually generated and exercised by the suite — they may be orphan golden files. Regenerate the goldens from a clean checkout via the case generator, and include the registry entry that produces these hashes.

B. Hot-path JWL branch is a runtime test, not compile-time gated.
The live five-equation path is still else if (jwl_idx > 0 .and. model_eqns == model_eqns_5eq) evaluated per face in all three Riemann solvers (and else if (jwl_idx > 0) in s_compute_speed_of_sound, m_data_output, m_time_steppers). The statically-dead 6eq/4eq/bubbles copies were correctly removed, but the live branch was not moved behind an fpp / case-optimization parameter, so the commit description "compile-time gated" is inaccurate. Either gate it at compile time consistent with MFC's case-specialization model, or adjust the description to match what the code does.

C. jwl_mix_type = 2 per-cell root-find cost.
The p-T-equilibrium closure still contains two do it = 1, 60 bisection loops with exp() calls inside the loop body (m_jwl.fpp, ~lines 320–338 and ~392–407). The default isobaric closure is correctly closed-form and jwl_pure_cutoff skips near-pure cells, but mixed cells under jwl_mix_type = 2 pay the full 60-iteration transcendental cost per face on GPU. Green CI does not measure this (the regression cases use the cheap closures). Please either provide a throughput measurement for mode 2 or document it as an experimental/host-preferred path; the revision should not be described as having reduced this cost, since the loops are unchanged.

D. Residual silent mask with a misleading comment.
s_compute_speed_of_sound now has c = sqrt(max(c2, sgm_eps)) with a comment stating a genuinely bad c2 (e.g. NaN) propagates so the failure surfaces. max(NaN, sgm_eps) is compiler/flag-dependent and frequently returns the finite argument, so a NaN sound speed can be silently floored to sqrt(sgm_eps) rather than surfacing. Either make the NaN handling explicit (self-compare check) or correct the comment so it does not claim a guarantee the code does not provide.

E. Unconditionally-passed Y_jwl can be uninitialized.
In m_data_output.fpp and m_time_steppers.fpp, Y_jwl is assigned only inside if (jwl_idx > 0) but is passed by keyword on every call:

if (jwl_idx > 0) then
    Y_jwl = min(max(q_prim_vf(jwl_idx)%sf(j, k, l)/max(rho, sgm_eps), 0._wp), 1._wp)
end if
call s_compute_speed_of_sound(..., jwl_Y=Y_jwl)

When jwl_idx == 0 this passes an uninitialized GPU-private variable. The callee only reads it under its own jwl_idx > 0 guard so behavior is currently correct, but this is undefined-behavior-flavored and -Wuninitialized/sanitizers will flag it. Initialize Y_jwl = 1._wp unconditionally.

F. Kuhl R_mix guard is a magnitude floor, not a sign check.
The Kuhl closure uses max(rho_safe*R_mix, sgm_eps), which clamps a negative R_mix to +sgm_eps and flips the denominator sign rather than failing. Inputs are validated positive in init so this cannot trigger in practice, but the prior request was an explicit R_mix > 0 check; please either add the sign check or document why R_mix is provably positive.

G. Remaining ;-packed statements.
Two multi-statement lines remain in init code (m_jwl.fpp ~700, ~704): jwl_cv_prod = 0._wp; jwl_cv_air = 0._wp and jwl_cv_air = fluid_pp(i)%cv; exit. Unpack to one statement per line.

H. 3D example with published comparison.
Per @danieljvickers' note, a feature of this scope warrants a 3D case compared against a published reference, and the code-vs-literature comparisons that exist should be included in the README alongside the cases. This is still outstanding.

---

Once A–H are addressed, this looks ready for a final human pass. The physics and the structural refactor are in good shape.

🤖 Generated by Claude Code at @sbryngelson's request. A human maintainer will make the final call.

[sbryngelson]

Automated re-review — additional code-quality items (run on behalf of @sbryngelson)

🤖 Follow-up to the previous automated pass. Those items (A–H) still stand. This pass covers MFC-specific code-quality expectations not yet raised, each verified against the current diff. Items continue the same lettering.

I. Closure test coverage.
The PR ships four mixture closures (jwl_mix_type 0/1/2/3 = isobaric / Kuhl / p-T-equilibrium / Rocflu), but only the isobaric default is exercised. Neither 1D_jwl_single_material_shocktube/case.py nor 2D_jwl_mixture_test/case.py sets jwl_mix_type, so both run mode 0, and the two golden cases (tests/25B9827F, tests/E3460991) inherit that. The Kuhl, p-T-equilibrium, and Rocflu closures have no automated regression coverage — the only evidence they work is the hand-written convergence table in 2D_jwl_mixture_test/README.md. Please add a golden case per closure, or gate the uncovered closures off as experimental until they have coverage. Shipping selectable GPU code paths that CI never executes is the largest remaining quality gap.

J. Single-precision (wp = real32) portability.
The p-T-equilibrium bisection convergence tests use fixed relative tolerances below single-precision epsilon (~1.2e-7):

! m_jwl.fpp:331
if (abs(f_m) <= 1.e-12_wp*max(abs(pres), sgm_eps)) exit
! m_jwl.fpp:400
if (abs(g_m) <= 1.e-12_wp*p_s) exit

In a single-precision build these can never be satisfied, so both do it = 1, 60 loops always run the full 60 iterations on the hot path. Relatedly, the pure-cell shortcut Y >= 1._wp - 1.e-12_wp collapses to Y >= 1._wp in single precision, silently shifting the branch threshold. Scale these tolerances with epsilon(1._wp) (and add a convergence-failure diagnostic if the bracket never tightens).

K. EOS selectors are bare magic numbers.
fluid_pp%eos is compared against literal 1/2, and jwl_mix_type against literal 1/2/3, throughout (m_jwl.fpp:667,703,710, m_variables_conversion.fpp), including in the abort message text ("use 1 for stiffened gas or 2 for JWL"). MFC's surrounding hot-path code uses named constants (model_eqns_5eq, etc. from m_constants). Define named parameters (e.g. eos_stiffened_gas, eos_jwl, jwl_mix_isobaric, …) and use them in the comparisons and messages.

L. Validation lives outside the checker convention.
All JWL parameter validation is in s_initialize_jwl_module; nothing was added to any m_checker* module, which is where MFC validates case inputs. The aborts themselves are correct (standard s_mpi_abort with informative messages — good), but they're in the wrong stage. Compounding this, the comment added in m_global_parameters_common.fpp states the air parameters "are validated in m_checker_common," which is not true — no such check exists. Move the validation to the checker stage, or at minimum correct the stale comment.

M. Missing artifacts for a feature of this scope.

• No benchmark case under benchmarks/, despite the change touching the Riemann hot path and adding a 60-iteration root-find.
• No theory/equations docs page — only the case.md parameter table. 2D_jwl_mixture_test/README.md references a root README-JWL-EOS.md that is not committed (dangling reference). Either commit that file or drop the reference.

N. Incomplete/inconsistent Doxygen.
New routines have !> summary lines but no !! @param/!! @return blocks. The new jwl_Y/jwl_alpha args on s_compute_pressure are added with no !< documentation, while the same args on s_compute_speed_of_sound do have !< field comments — inconsistent. Bring the new public-interface docs up to the codebase's Doxygen level.

O. Cold-curve duplication across closures.
The cold-curve expression A*(1 - omega*rho1/(R1*rho0))*exp(...) + B*... is recomputed verbatim in s_jwl_pressure_er, s_jwl_energy_pr, and the Kuhl/p-T/Rocflu variants, even though an s_jwl_pcold helper exists and is used elsewhere. Routing all callers through the helper would remove the duplication and shrink the module (this is the concrete form of the module-size concern raised earlier in the thread).

P. present() of optionals inside GPU [seq] routines.
s_compute_speed_of_sound relies on present() of optional arguments inside an offloaded [seq] routine, with NaN-self-compare guards. The non-short-circuit caveat is documented in comments, but present() on optionals inside device routines has historically been fragile across NVHPC versions. CI is green today; please confirm this is intentional and considered portable, since it's an easy source of a future silent device/host divergence.

---

For completeness, these were checked and are fine: fluid_pp%eos defaults to 1 so existing cases are unaffected; toolchain registration is complete via params/definitions.py and descriptions.py (this layout has no case_dicts.py); no TODO/FIXME left in the new source; and the PR correctly fixes several neighboring 5.e-1 → 5.e-1_wp suffixes on lines it touches.

🤖 Generated by Claude Code at @sbryngelson's request. A human maintainer will make the final call.

[fahnab666]

I am still working on this PR in more detail. I am currently changing the implementation format to make the JWL mixture-closure logic cleaner, more modular, and easier to review. The goal is to move toward a better closure-model structure where each approach is more clearly separated, documented, and easier to validate.

So far, I have been focusing on the mixture EOS closure workflow, including the connection between the 5-equation model, the explosive-products mass fraction, the caloric EOS, and the thermal EOS. I am also comparing different closure directions, including p–T equilibrium, Kuhl/Khasainov-style interpolation, linear/quadratic fitting approaches, and possible reduced-fit models for mixture pressure, temperature, and sound speed.

I will keep refining the implementation and keep up here in the comment. The code is far for ready.