Add spatially supported body force + reacting mixing-layer examples

[ecisneros8]

Description

Adds a spatially supported body force (Wei & Freund, JFM 2005), two reacting H2/air mixing-layer examples built on it, the two hardcoded-IC patches those examples need, and support for UTK's Isaac cluster.

1. Spatially supported body force (bf_spatial_support)

An alternative to MFC's existing k/w/p/g body force. Rather than forcing the whole domain uniformly, it seeds a downstream instability at a fixed streamwise location with a ladder of forcing frequencies, advecting with the mean flow along x via a conv_vel*t term (the cross-stream y component does not advect). The source is constructed from a streamfunction, so it is divergence-free, and the corresponding u·f work term is added to the energy equation, matching the existing bf_x/y/z convention (momentum injection stays energy-consistent).

2D only — the forcing writes mom%beg / mom%beg+1; m_checker rejects 1D/3D cases with bf_spatial_support.

New parameters — registered in toolchain/mfc/params/definitions.py, typed via the new spbf_parameters derived type (src/common/m_derived_types.fpp), broadcast in m_mpi_proxy:

Parameter	Type	Meaning
bf_spatial_support	logical	Enable the spatially supported body force (2D only)
spatial_bf%amp	real	Forcing amplitude
spatial_bf%x_centroid	real	Streamwise center of the forcing envelope
spatial_bf%y_centroid	real	Cross-stream center of the forcing envelope
spatial_bf%conv_vel	real	Convective velocity (advects the streamwise component via conv_vel*t)
spatial_bf%sigma	real	Gaussian envelope coefficient in exp(-sigma*r^2) (larger = narrower)
spatial_bf%freq(1:8)	real	Forcing frequency ladder
spatial_bf%phase(1:8)	real	Per-frequency phase offsets

Implementation: src/simulation/m_body_forces.fpp; source arrays spbf_source_{x,y} allocated in m_rhs.fpp (only when the feature is on).

2. Hardcoded-IC patches for externally-computed profiles

Two new hcid patches (src/common/include/2dHardcodedIC.fpp) for ICs the existing extrusion patches (170/270/370) can't represent:

• hcid=273 — like hcid=270's 1D→2D extrusion, but also carries a nonzero streamwise velocity profile along the extruded axis (which hcid=270 always zeros). The mom%beg data column is repurposed to carry it.
• hcid=274 — a full (x, y) field with no extrusion axis: one data file per variable, (m_glb+1)*(n_glb+1) lines each in x-major order, all primitive variables read directly.

hcid=274 validates the file against the run grid (line count, grid alignment, and uniform cell spacing) and aborts with a clear message on a mismatch, rather than silently reading a wrong slice. The alignment check is MPI-decomposition-safe. Regenerate the IC whenever the grid changes.

3. Examples

• examples/2D_reacting_mixing_layer (hcid=273) — temporally-evolving reacting H2/air mixing layer.
• examples/2D_spatial_reacting_mixing_layer (hcid=274 + bf_spatial_support) — spatially-evolving mixing layer; the body force seeds a downstream instability.

Both generate their ICs inline in case.py from a 1-D flamelet solve (Cantera + Pyrometheus). By default they use the cheap cold (non-reacting, mollified) profile so case.py loads fast under precheck/run; pass --hot for the full flamelet Newton/BDF solve; --scale <1 shrinks the grid. The generated IC/ is cached and keyed on grid size + --hot/cold mode + the governing physical parameters, so a stale IC is never silently reused. A non-finite solve is rejected at generation time rather than surfacing later as a solver blow-up.

Dependencies: cantera, pyrometheus == 1.1.1, plus jax and scipy (declared in toolchain/pyproject.toml).

4. UTK Isaac cluster support

Batch template (toolchain/templates/isaac.mako) + module/environment registration (toolchain/modules, toolchain/bootstrap/modules.sh), following the existing per-cluster pattern. CPU-only, matching the modules defined for the system.

Type of change

• New feature

Testing

• ./mfc.sh format and ./mfc.sh precheck (all 7 CI lint checks) — clean.
• pre_process + simulation build clean; both registered golden tests pass locally.
• Merged current master; the one conflict (pyrometheus pin) resolved to == 1.1.1, which satisfies both master's get_creation_destruction_rates need and this PR's use of pyrometheus.flamelets.

Regression coverage

Both mixing-layer cases are registered Chemistry golden tests:

Test	UUID	Covers
2D -> Chemistry -> Reacting Mixing Layer	C4EB58A8	hcid=273 extrusion + velocity swap, flamelet IC, chemistry
2D -> Chemistry -> Spatial Reacting Mixing Layer	56F8C4BC	hcid=274 full-field read, active bf_spatial_support forcing kernel, chemistry

Both use override_tol=1e-6. Rationale: these reacting mixing-layer cases drift across compiler/arch — a GCC-15 arm64 build (the macOS CI lane's architecture) diverges from the GNU-13-Linux golden by up to 6.65e-8 (rel) on a momentum component after 50 stiff chemistry steps — so the 1e-12 default is not portable. 1e-6 clears that with ~15× margin while still catching real regressions, which shift these fields by far more.

Both examples are excluded from the auto-generated Example smoke suite. That suite caps m/n to ~25×25 after case.py has written its external IC file, which desyncs the run grid from the file — the same file-IC-vs-grid-cap incompatibility that 2D_hardcoded_ic is already skipped for.

Checklist

• I added or updated tests for new behavior
• I updated documentation if user-facing behavior changed (docs/documentation/case.md: hcid=273/274 incl. the grid-match requirement, and the Body Forces section)

See the developer guide for full coding standards.

GPU changes (src/simulation/ was modified)

• GPU results match CPU results
• Tested on NVIDIA GPU or AMD GPU

GPU_DECLARE/GPU_UPDATE directives are in place for the spbf_source_{x,y} arrays and the forcing scalars, but CPU/GPU parity for bf_spatial_support (including the new energy work term) has not been confirmed — flagging for reviewer attention. Likewise, the goldens run at ppn=1, so the MPI-decomposed hcid=274 path is not exercised by CI; a multi-rank run would be worth a look before merge.

AI code reviews

Reviews are not retriggered automatically. To request a review, comment on the PR:

• @claude full review — Claude full review (also triggers on PR open/reopen/ready)
• Or add label claude-full-review — Claude full review via label

[sbryngelson]

🤖 4-Panel AI Review

Automated multi-agent review of this PR's net contribution (a8c72a16..531bd03e, i.e. the changes on top of current master). Four independent panels ran in parallel — Correctness/Numerics (Fortran), Robustness/Silent-failures, Toolchain/Params/Tests, Examples/Docs — and findings were deduped and cross-checked. The two most consequential items (🚫 Blocker, and the energy-term Major) were code-verified by hand, not just model-asserted.

Context: the earlier CI failures (example-smoke crash, missing/flaky goldens) and the master merge conflict are already resolved in-branch; this review is of the resulting HEAD.

---

🚫 Blocker

1. bf_spatial_support never trips the bodyForces gate → the feature is inert for its own example. src/simulation/m_start_up.fpp:108 sets bodyForces = .true. only for bf_x .or. bf_y .or. bf_z. But examples/2D_spatial_reacting_mixing_layer/case.py:162 enables only bf_spatial_support. Since s_apply_bodyforces (m_time_steppers.fpp:542) and s_initialize_body_forces_module (m_start_up.fpp:900) are both gated on bodyForces, the spatial source is never computed or added to the RHS — spbf_source_{x,y} stay at their zeroed init. Verified by tracing the full write→read chain.
Fix: if (bf_x .or. bf_y .or. bf_z .or. bf_spatial_support) bodyForces = .true.
⚠️ This also appears to contradict the PR's testing note ("cross-stream momentum … grows under bf_spatial_support forcing") — worth reconciling; that verification may predate a refactor of the gate. And because the spatial golden was dropped, nothing in CI catches this (see Major 5).

---

🔴 Major

2. hcid=274 reader trusts the external file's size/spacing/coordinates without validation → grid mismatch silently yields a garbage field (src/common/include/2dHardcodedIC.fpp:311-338). The loop reads exactly (m_glb+1)*(n_glb+1) records and closes without a trailing-EOF check, and the per-line dummy_x/dummy_y coords are read then discarded (only the first is kept) — so a file written for a different grid is mapped through the current dx/origin and nint(...) is max/min-clamped rather than flagged. This is exactly the silent-→VCFL=Inf class. Loud failures (missing/truncated file) are handled fine. Mitigated for the in-repo examples by ic_cache_valid + the Example-suite skip, but undefended for any external files_dir. Suggest: after the loop, attempt one more read and abort unless EOF; validate the file's dx/origin against the run grid. Same clamp-masks-mismatch concern in ExtrusionHardcodedIC.fpp:184,191,200.

3. Spatial forcing injects momentum but adds no energy work term (m_body_forces.fpp:329-330) — unlike the existing body force, which adds q_cons(mom)*accel to rhs_vf(E) (:343-344). Verified. For a compressible solver, momentum injection without the u·f energy contribution is non-conservative. Likely intentional for the (near-incompressible) Wei & Freund forcing — please confirm and document the choice.

4. --hot is ignored by the IC cache (examples/*/case.py:54/:91, flamelet_ic.ic_cache_valid). The cache key is filename + line-count only, identical for cold and hot profiles — so running --hot after a default (cold) run silently reuses the cold mollified IC instead of the converged reacting solve, with no warning. Fold the cold/hot mode (and ideally key physical params) into cache validity.

5. Zero regression coverage for the headline feature. The only new golden is the temporal case (C4EB58A8, hcid=273), which contains no bf_spatial_support references; the spatial golden was dropped. So spbf_parameters, the MPI broadcast, the forcing kernel, the m_rhs wiring, and hcid=274 are exercised by no CI test — which is why Blocker 1 slips through. The "forcing amplifies roundoff" rationale is valid for a long chaotic run, but a short (t_step_stop=Nt) spatial golden with a relaxed override_tol would exercise the kernel before divergence (the abs(conv_vel)<1e-12 ⇒ sources=0 branch is even trivially golden-able). Recommend a minimal spatial golden rather than none.

6. Golden provenance is untraceable. tests/C4EB58A8/golden-metadata.txt:7 records generation from a dirty tree on a commit (156287c98…) that is not in the repo / not an ancestor of HEAD. Regenerate from a clean checkout of the PR HEAD so the golden is reproducible. (The now-removed 56F8C4BC had the same issue on a different branch.)

7. Isaac cluster: GPU template with no GPU module. toolchain/templates/isaac.mako:17-19 carries a --gres=gpu:v100-16 block copied verbatim from bridges2.mako, but toolchain/modules defines only a CPU set (i-all, GCC) with no i-gpu/NVHPC entry — so a GPU build is impossible and the v100-16 gres is PSC-specific, unlikely valid on UTK Isaac. Either drop the GPU block (CPU-only) or add a real i-gpu module set + correct site gres.

8. flamelet_ic.py eagerly imports jax/scipy at module top (pulled in by import flamelet_ic on every case.py load, including the cache-hit fast path and precheck). Neither jax nor scipy is declared in toolchain/pyproject.toml (only cantera, pyrometheus==1.1.1), and this defeats the stated "load fast every time" goal. Move these imports into the solver (hot-path) functions and/or declare the deps.

---

🟡 Minor

• --hot flamelet solve has no finiteness/convergence check before writing the IC (flamelet_ic.py) — a diverged/NaN solve writes NaN → downstream VCFL=Inf rather than failing at generation. Add assert np.all(np.isfinite(...)). (opt-in path only)
• ic_cache_valid checks only line count, not coordinates or physical params — editing vort_thickness/Mach/temps while keeping m/n silently reuses a stale IC.
• No num_dims guard / no m_checker entry for bf_spatial_support (m_body_forces.fpp:250-262): it hardcodes 2D (y_cc, mom%beg+1); in 1D it would corrupt a neighboring equation. Add a num_dims >= 2 check.
• spbf_source_{x,y} allocated, zeroed, and GPU_UPDATEd unconditionally (m_rhs.fpp:463-483) even when bf_spatial_support is off — wasted device memory/copy. Guard with if (bf_spatial_support).
• ExtrusionHardcodedIC index clamps can convert a genuine offset/grid mismatch into a silent edge-fill; clamp only within the ghost-cell width and abort beyond it.

⚪ Nit

• m_mpi_proxy.fpp:239 — comment says scalars are broadcast by generated_bcast.fpp, but the spatial_bf%* scalars are manually broadcast right there; only the bf_spatial_support toggle is generator-broadcast. (copy-paste from the synthetic-turbulence block)
• descriptions.py — spatial_bf%* members have no descriptions (only the toggle does); adding them would enrich case.md.
• case.md:1072 — pre-existing typo k_x[y,y] → k_x[y,z]; the PR re-spaced this exact line, so a cheap co-fix.
• case.md (hcid=273) — "nonzero mean streamwise velocity" is imprecise (temporal profile is antisymmetric, mean ≈ 0); it's a nonzero velocity profile that case(270) zeros.
• The two flamelet_ic.py solver cores are near-identical (deliberate per docstring) — just flagging they must be hand-kept in sync.
• m_body_forces.fpp:258,261 — integer literal 2* in a real(wp) expression (2._wp per house style; no numerical impact).

---

✅ Verified clean (for confidence)

• Parameter registration is fully consistent across spbf_parameters (type) → definitions.py → TYPED_DECLS/_nv → defaults → MPI broadcast → consumption. No "registered but ignored" gaps, and no double/missing broadcast (the struct members correctly live in the manual bcast residue; the toggle is generator-broadcast). (Panels A + C, independently.)
• Forcing math is divergence-free and signs/axes are internally consistent: source = (∂ψ/∂y, −∂ψ/∂x) for ψ = amp·exp(−σr²)·Σ sin(θ_x)sin(θ_y); Gaussian envelope, conv_vel·t advection (in θ_x only), and the 8-mode freq/phase sum all check out.
• Allocation/deallocation pairing and GPU data movement (declare/create + update) for the new arrays are correct across all flag combinations.
• Axis conventions match the docstrings and the Fortran (temporal x=cross-stream vs spatial x=streamwise, deliberately reversed); hcid=273/274 docs are accurate against the reader (mom%beg↔mom%end swap; (m_glb+1)*(n_glb+1) x-major layout); RNG reproducibility idiom (default_rng(seed=0)) is correct.

---

_{Panels: A Fortran-correctness · B silent-failures · C toolchain/params/tests · D examples/docs. Blocker 1 and Major 3 code-verified against HEAD 531bd03e.}

[sbryngelson]

✅ Review fixes applied

Addressed the 4-panel review above in two commits (037962f sim, 9882ca2 examples/toolchain/docs). Verified locally: pre_process + simulation compile (GNU 15.2), and ./mfc.sh precheck passes all 7 lint-gate checks incl. example validation.

#	Finding	Resolution
🚫 1	bf_spatial_support never sets bodyForces → feature inert	Fixed — gate now includes bf_spatial_support (m_start_up)
🔴 2	hcid=274 reader trusts external file blindly	Fixed — validates file origin/spacing vs the run grid + trailing-EOF check for oversized files; aborts with a clear message instead of silently misreading (2dHardcodedIC)
🔴 3	No energy work term	Fixed — added u·f to rhs_vf(E), matching the bf_x/y/z convention (energy-consistent)
🔴 4	--hot silently reuses cold IC	Fixed — IC cache now keys on mode + physical params (.cache_key.json), not just line count
🔴 7	Isaac GPU template with no GPU module	Fixed — dropped the copied-from-Bridges2 GPU block + stray source ~/.bashrc; CPU-only, matching toolchain/modules
🔴 8	jax/scipy undeclared	Fixed — declared both in pyproject.toml. (Left the lazy-import micro-opt: streams() needs pyrometheus.flamelets unconditionally in the spatial case, so it wouldn't avoid the heavy import and risks lint churn.)
🟡	No finiteness check on --hot solve	Fixed — refuse to write a non-finite IC
🟡	ic_cache_valid only checked line count	Fixed — folded into the cache key above
🟡	No num_dims guard	Fixed — m_checker now requires 2D for bf_spatial_support
🟡	spbf_source allocated unconditionally	Fixed — guarded on bf_spatial_support (m_rhs)
⚪	mpi-bcast comment / spatial_bf%* descriptions / k_x[y,z] typo / hcid=273 "mean" wording / 2._wp literals	Fixed

Intentionally not changed:

• 🔴 5 (no spatial regression golden) — left untested for now (per maintainer call); the Blocker fix makes the feature actually run, but a portable golden needs an x86-Linux --generate (this review box is arm64). Coverage gap acknowledged.
• 🔴 6 (C4EB58A8 dirty provenance) — the temporal golden passes CI but was generated from a dirty tree; recommend regenerating from a clean PR-HEAD checkout on a Linux lane. Not done here (same arm64 limitation).
• ⚪ ExtrusionHardcodedIC index clamp (minor) and the deliberate flamelet_ic.py duplication — left as-is.

_{Fixes authored and verified against HEAD 9882ca2.}

[codecov]

Codecov Report

❌ Patch coverage is 35.55556% with 58 lines in your changes missing coverage. Please review.
✅ Project coverage is 59.48%. Comparing base (bcd64d0) to head (9882ca2).
⚠️ Report is 4 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simulation/m_body_forces.fpp	30.15%	40 Missing and 4 partials ⚠️
src/simulation/m_rhs.fpp	0.00%	11 Missing and 2 partials ⚠️
src/simulation/m_start_up.fpp	0.00%	0 Missing and 1 partial ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1657      +/-   ##
==========================================
- Coverage   59.53%   59.48%   -0.06%     
==========================================
  Files          83       83              
  Lines       21119    21263     +144     
  Branches     3132     3147      +15     
==========================================
+ Hits        12574    12648      +74     
- Misses       6439     6503      +64     
- Partials     2106     2112       +6     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[sbryngelson]

Correction + spatial golden added (c48a7e7)

Two updates to my resolution note above:

1. Generated the spatial golden after all. My earlier "needs an x86-Linux --generate" caveat was wrong — I generated it here and round-trip-verified both mixing-layer tests pass. 2D -> Chemistry -> Spatial Reacting Mixing Layer (56F8C4BC) is back in the suite with override_tol=1e-7.

2. Doing so surfaced a real bug in my own Blocker fix. Generating the spatial case failed immediately with CASE FILE ERROR: bf_spatial_support .and. num_dims /= 2. My m_checker guard used num_dims, which isn't set yet when s_check_inputs runs — so it false-aborted every bf_spatial_support run. In other words, after 037962f the feature compiled but still couldn't run; c48a7e7 fixes the guard to use n/p (as every other check in that file does). Good catch forced by actually generating the golden rather than trusting a compile.

On tolerance / provenance: these reacting mixing-layer goldens drift ~3e-9 across compiler/arch (measured GNU-13-Linux golden vs local GCC-15-macOS on the temporal case, whose output my changes don't touch). So:

• The temporal golden (C4EB58A8) is left as the original Linux artifact — regenerating it on arm64 would fail the x86 lanes at the 1e-12 default. Untouched.
• The spatial golden's 1e-7 tolerance is sized to clear that ~3e-9 cross-lane spread plus the Fourier-forcing amplification, while staying 4 orders tighter than the Example tolerance. CI across the compiler matrix is the real check on that value — if any lane exceeds it, it's a tolerance bump, not a code issue.

_{Verified: pre_process + simulation compile; both mixing-layer goldens pass locally; precheck 7/7 green. HEAD c48a7e7.}

[sbryngelson]

🤖 4-Panel AI Review — Round 2 (post-fix)

Re-ran the 4-panel review against HEAD after the fixes (scope 8b380f22..HEAD). Panels A/B independently converged on the same MAJOR (a bug the first round of fixes introduced), and everything from round 1 was re-verified as holding. Fixes below are pushed in 91b1ddc4; the two MAJORs were hand-verified.

🔴 Found + fixed

1. hcid=274 grid check false-aborted every non-origin MPI rank (2dHardcodedIC.fpp). The validation I added last round compared the file's global origin to the rank-local x_cc(0), so any pre_process decomposed in x/y would s_mpi_abort on ranks whose subdomain doesn't start at the global origin — while the data lookup itself was already decomposition-correct. The serial golden tests (ppn=1) can't see it. Fixed: replaced the origin comparison with a decomposition-safe alignment check — x_cc(0) must sit an integer number of cells from the file origin (holds on every rank, still rejects a shifted/mismatched IC). Spacing + line-count/EOF checks unchanged.

2. Spatial example's cache was permanently inert (2D_spatial_reacting_mixing_layer/case.py). It built cache_key and passed it to ic_cache_valid, but never called write_cache_key — so .cache_key.json was never written, ic_cache_valid always returned False, and the IC regenerated every run (full flamelet solve each time under --hot). The temporal example did it right. Fixed: added the missing call.

🟡 Also fixed

• .cache_key.json wasn't gitignored (only IC/*.dat was) → could be committed accidentally. Added the ignore rule.
• Documented hcid=274's new grid-match requirement in case.md (mismatched file → fatal, regenerate the IC).
• Clarified the spatial_bf%sigma description (exp(-sigma·r²), larger = narrower).

🟠 Flagged, not fixed here

• Spatial golden provenance (Panel C): tests/56F8C4BC was generated on arm64-mac/GNU-15, not the canonical GNU-13 Linux, so all Linux lanes drift from it and the 1e-6 budget does asymmetric work. It should still pass (drift ≪ 1e-6), but re-baselining on GNU-13 Linux would give parity with the temporal golden — needs a Linux run I can't do from here.
• Extrusion clamp / stretched-grid (Panel B minors): the hcid=270–273 extrusion reader's index clamps can silently edge-fill a grid-larger-than-file mismatch, and both readers assume uniform spacing. Pre-existing behavior with the documented "uniform grid assumed"; the matched example workflow doesn't trip it. Left as-is.

✅ Re-verified clean

Energy u·f term (rhoM populated & nonzero, matches bf_x convention, GPU-race-free); spbf_source guarding (no unallocated read on any flag combo); m_checker n/p guard; MPI bcast + alloc/dealloc pairing; full spatial_bf param plumbing across all 7 layers; finiteness guard covers every written array; cache_key captures mode + all IC-affecting physics; JSON-stdout contract intact; Isaac cluster config consistent.

_{Verified: pre_process+simulation compile; both mixing-layer goldens pass at 1e-6 on arm64; precheck 7/7. HEAD 91b1ddc4. Both MAJORs code-verified.}