Combustion: condensed-phase reactive burn (5-/6-equation, Arrhenius kinetics) + IBM burning surfaces#1670
Conversation
…gy at immersed-boundary points
…al-ballistics feedback
…ng cylinder, wake flame)
…nst negative pressure
…nerated on master base)
|
@claude full review |
|
I'll analyze this and get back to you. |
|
@claude full review |
There was a problem hiding this comment.
Pull request overview
This PR adds two new combustion capabilities to MFC’s simulation target—(1) a condensed-phase programmed reactive-burn source for the multi-fluid model and (2) burning/transpiring immersed-boundary (IBM) surfaces with optional pressure-coupled burn-rate scaling—along with new example cases, parameter plumbing, documentation updates, and regression goldens.
Changes:
- Introduces
m_reactive_burnand integrates it into the RHS path, with new runtime parameters (reactive_burn,rburn_*) and documentation. - Extends IBM patches to support fuel injection (
v_blow,inj_species) and optional Vieille’s-law scaling (burn_rate_exp,burn_rate_pref), plus chemistry-related ghost-state handling updates. - Adds new toolchain test cases/examples and accompanying regression golden files/metadata; improves parameter validation error suggestions.
Reviewed changes
Copilot reviewed 32 out of 37 changed files in this pull request and generated 2 comments.
Show a summary per file
| File | Description |
|---|---|
| toolchain/mfc/test/cases.py | Adds a new regression test case definition for condensed-phase reactive burn. |
| toolchain/mfc/params/validate.py | Improves unknown family-attribute errors by ordering valid attributes using fuzzy suggestions. |
| toolchain/mfc/params/descriptions.py | Adds user-facing descriptions for reactive_burn and rburn_* parameters. |
| toolchain/mfc/params/definitions.py | Registers reactive-burn params and adds new IBM patch attributes (inj_species, v_blow, burn_rate_*). |
| tests/F200F862/golden-metadata.txt | Adds generated golden metadata for a new/updated regression test. |
| tests/F08B0D0E/golden.txt | Adds generated golden output for a new/updated regression test. |
| tests/F08B0D0E/golden-metadata.txt | Adds generated golden metadata for a new/updated regression test. |
| tests/E5B66084/golden-metadata.txt | Adds generated golden metadata for a new/updated regression test. |
| tests/CF4FE4C9/golden-metadata.txt | Adds generated golden metadata for a new/updated regression test. |
| tests/85405397/golden-metadata.txt | Adds generated golden metadata for a new/updated regression test. |
| tests/8430E4EA/golden-metadata.txt | Adds generated golden metadata for a new/updated regression test. |
| tests/64F05A05/golden-metadata.txt | Adds generated golden metadata for a new/updated regression test. |
| src/simulation/m_rhs.fpp | Hooks reactive-burn source application into the simulation RHS. |
| src/simulation/m_reactive_burn.fpp | New module implementing the programmed pressure-driven reactive-burn source term. |
| src/simulation/m_mpi_proxy.fpp | Broadcasts newly added IBM patch parameters in simulation. |
| src/simulation/m_ibm.fpp | Adds IBM surface injection/burning behavior and chemistry-consistent ghost state handling. |
| src/simulation/m_global_parameters.fpp | Initializes defaults for new IBM patch parameters in simulation. |
| src/simulation/m_checker.fpp | Adds reactive-burn precondition checks (num fluids, EOS consistency, qv ordering). |
| src/pre_process/m_mpi_proxy.fpp | Broadcasts newly added IBM patch parameters in pre_process. |
| src/pre_process/m_global_parameters.fpp | Initializes defaults for new IBM patch parameters in pre_process. |
| src/common/m_global_parameters_common.fpp | Initializes defaults for new reactive-burn parameters. |
| src/common/m_derived_types.fpp | Extends patch_ib derived type with injection/burn-rate parameters. |
| src/common/include/2dHardcodedIC.fpp | Adds hardcoded IC (HCID 275) for reactive shock–flame initialization. |
| examples/2D_shock_flame/case.py | New example case for shock–flame (RM) interaction with chemistry enabled. |
| examples/2D_ibm_flameholder/case.py | New example case exercising IBM injection + chemistry. |
| examples/2D_ibm_burning_grain/case.py | New example case for IBM “burning grain” with optional pressure-coupled injection scaling. |
| examples/2D_hybrid_slab/case.py | New example case for wall-injected hybrid-rocket-like diffusion flame. |
| examples/1D_reactive_burn/case.py | New 1D example for condensed-phase reactive burn detonation. |
| examples/1D_propellant_flame/case.py | New 1D example for a premixed propellant-flame surrogate with chemistry. |
| docs/module_categories.json | Registers m_reactive_burn for documentation/module categorization. |
| docs/documentation/case.md | Documents new reactive_burn/rburn_* parameters and their meaning. |
Comments suppressed due to low confidence (1)
src/simulation/m_ibm.fpp:889
- Ys_IP is optional, but this block is guarded by the global
chemistryflag rather thanpresent(Ys_IP). That assumes everychemistrycall path always passes the optional array, which is fragile and can lead to invalid memory access.
if (chemistry) then
$:GPU_LOOP(parallelism='[seq]')
do q = 1, num_species
Ys_IP(q) = Ys_IP(q) + coeff*q_prim_vf(eqn_idx%species%beg + q - 1)%sf(i, j, k)
end do
| pres_IP = 0._wp | ||
| vel_IP = 0._wp | ||
|
|
||
| if (chemistry) Ys_IP = 0._wp |
|
I'll analyze this and get back to you. |
Addresses PR review (Copilot): - reactive_burn divides the pressure drive by rburn_pref; prohibit rburn_pref <= 0 - s_interpolate_image_point selects the bubbles/QBMM branch before chemistry, so the Ys_IP species state is not carried when both are on; forbid chemistry + bubbles_euler/qbmm (Fortran m_checker + Python case_validator) rather than touch the hot IBM path
|
Thanks — both are valid. Addressed in c1c5ada:
|
Addresses independent review findings on MFlowCode#1670: - reactive_burn's exact volume-fraction swap and single-pressure read assume the 5-equation pressure-equilibrium model, but num_fluids=2 alone also admits model_eqns=3; prohibit model_eqns /= 2 (m_checker + case_validator.check_reactive_burn + PHYSICS_DOCS) - Vieille's-law v_blow ~ p^n diverges for a negative exponent as surface pressure -> 0; prohibit patch_ib%burn_rate_exp < 0 (case_validator.check_ibm)
… kinetics Extends the condensed-phase reactive burn beyond the 5-equation model to the 6-equation (Saurel/Pelanti) multi-fluid model, and adds an optional temperature-driven ignition term. 6-equation support: the reactant->product qv release is carried through the phasic-pressure relaxation. m_pressure_relaxation.fpp threaded the formation energy (alpha_rho_k*qv_k) through s_equilibrate_pressure and s_correct_internal_energies so the phasic pressures and the recovered mixture pressure are qv-consistent with s_compute_pressure -- otherwise a nonzero qv inflated the phasic pressure by rho_k*qv_k/gamma_k. The existing 5-eq burn source then works unchanged on model_eqns=3. For qv=0 the new terms are exactly zero, so all 13 existing 6-eq relaxation goldens are bit-identical (verified CPU/OMP/ACC). Arrhenius kinetics: new rburn_ta (activation temperature [K]); default 0 keeps the pure pressure-driven rate bit-identical. When rburn_ta>0 the rate gains a exp(-rburn_ta/T_r) factor, with T_r=(p+pi_inf)/((gamma-1)*cv*rho) the reactant phasic temperature. Guarded (checker + case_validator) so rburn_ta>0 with cv=0 is rejected instead of silently disabling the factor. Test: adds a model_eqns=3 reactive-burn regression golden (nested off the 5-eq test), and hardens both reactive goldens by giving the uniform case a small advection velocity so the momentum output is well-resolved rather than near-zero roundoff (the old golden compared ~0 momentum at 1e-12 and failed across recompiles/backends). Validated on CPU, OpenMP-offload, and OpenACC.
Nested off the reactive-burn detonation test with rburn_ta=500 and cv=1500, chosen so the exp(-rburn_ta/T_r) factor is O(0.4) at the IC temperature -- exercising the temperature branch (all other reactive goldens use rburn_ta=0). Verified on CPU, OpenMP-offload, and OpenACC.
…ess review panel + CI Post-review hardening from a 5-reviewer panel, plus the CI-surfaced flame-example fragility. rburn derived type: the five rburn_* scalars (k, pign, pref, n, ta) become a rburn type instance, matching the sibling chem_params pattern (type in m_derived_types, TYPED_DECLS registration, member _r/_nv, descriptions regex, auto GPU_DECLARE + GPU_UPDATE(device=[rburn]) in m_global_parameters and m_start_up so the type is device-resident). Numerically inert (~1 ULP recompile roundoff); reactive goldens regenerated for the namelist-key change. Verified CPU/OMP/ACC. Completeness guards (review): rburn%k (>0), rburn%pign (set), rburn%pref (>0), rburn%n (>=0) are now validated in both m_checker and case_validator. Each defaults to the dflt_real sentinel, so an unset value silently produced spurious ignition (pign), NaN via drive**n (n), or a backward reaction (k). IB inj_species bound (review): m_checker now requires patch_ib%inj_species <= num_species when chemistry + ib; it indexes Ys_IP(1:num_species) at image points, so an out-of-range value was an out-of-bounds write. Docs (review): reactive_burn requires num_fluids = 2 (checker aborts otherwise), corrected from >= 2 in case.md and the module docstring. CI: skip the 1D_propellant_flame example (golden removed). Its finite-rate flame front drifts past the 1e-3 Example tolerance across compilers (nvhpc 25.11 golden vs 24.3/GNU/Intel/CCE/AMD) -- same non-portability class as the flame examples already skipped. The reactive_burn detonation goldens remain portable (they passed every CI lane).
The rburn refactor moved the five rburn_* scalars (auto-broadcast) into a derived type, but derived-type members are not auto-broadcast -- they are hand-emitted per type in fortran_gen.py (as for chem_params/lag_params). rburn had no emitter, so on a multi-rank reactive_burn run every non-root rank kept rburn at the dflt_real sentinel: a silent wrong-physics bug the single-rank golden tests cannot catch. Adds _emit_rburn (all members REAL -> mpi_p, under the reactive_burn guard, sim target), mirroring _emit_chem_params. Verified in the regenerated generated_bcast.fpp.
…roadcast All other reactive goldens run single-rank, so a broken rburn broadcast is invisible to them. This ppn=2 variant of the detonation case burns rank 1's half of the domain with the broadcast rburn; verified to FAIL when the rburn broadcast emitter is removed and PASS with it, so it locks the fix in CI.
The ibm_burning_grain/ibm_flameholder Example goldens already cover surface injection (v_blow + the inj_species composition swap), but both leave burn_rate_exp at 0, so the pressure-coupled scaling v_blow*(p/burn_rate_pref)^n shipped with no regression coverage. Registers the same example with --burn_exp 0.5 at the Example suite's 25x25/50-step size, so it differs from the existing golden only by the coupling: turning it on moves the solution by ~5e-3 relative (1.3 absolute in energy), well above the 1e-3 tolerance, so a broken multiplier fails the test. Tolerance pinned to 1e-3 to match the sibling Example golden -- compute_tolerance would otherwise pick 1e-10 via the ib branch, which a 50-step IBM + finite-rate-chemistry run cannot hold across compilers.
…pilers) The diffusion flame over the fuel slab is positioned by accumulated stiff-kinetics roundoff, so by step 50 the transverse momentum drifts past the 1e-3 Example tolerance: nvhpc lanes pass, Intel and Frontier CCE disagree with the golden by ~2e-3 absolute (up to 44% relative on a small component). Same non-portability class as 1D_propellant_flame and the other flame examples already skipped there, so no single golden passes every lane. The case remains a runnable example; the reactive_burn, shock_flame, ibm_burning_grain and ibm_flameholder goldens are unaffected and pass on all lanes.
…patches particle_cloud_ibs is allocated without default initialization, and s_add_cloud_particle set only the geometry/motion members. s_reduce_ib_patch_array then copies the whole struct into patch_ib, overwriting the zeroed defaults from s_assign_default_values_to_user_inputs -- so v_blow, inj_species and burn_rate_exp/pref reached the solver as uninitialized memory. That was latent until this branch added the first read of those members in m_ibm: a garbage v_blow > 0 superimposes a garbage wall-normal velocity on the ghost points, which NaNs the field. It reproduced only where the allocation is not already zero-filled -- Frontier AMD failed 2D_mibm_particle_cloud with 'ICFL is NaN' at step 2 while every NVHPC lane (and local CPU/OMP/ACC) passed. Sets the four fields for the inert particle surfaces; audited that these are exactly the patch_ib members m_ibm reads (airfoil_id/model_id/length_x remain unset but are never read for circle/sphere particles).
…mmon-pitfalls Both bit this branch, both fail silently, and neither is visible to local testing: a derived-type param regrouping dropped its MPI broadcast (non-root ranks kept the sentinel, invisible to single-rank goldens), and an uninitialized patch_ib member reached the solver through the whole-struct copy from particle_cloud_ibs (garbage v_blow -> ICFL NaN, reproducing only on Frontier AMD where the allocation is not zero-filled). Documents the required steps and the platform-only-NaN signature.
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #1670 +/- ##
==========================================
+ Coverage 59.62% 60.49% +0.87%
==========================================
Files 83 84 +1
Lines 21246 21189 -57
Branches 3147 3134 -13
==========================================
+ Hits 12668 12819 +151
+ Misses 6467 6303 -164
+ Partials 2111 2067 -44 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|

Description
Adds two combustion capabilities to the solver, along with example cases and regression tests. The work targets solid-fuel / condensed-phase and immersed-boundary combustion (internal ballistics, hybrid rockets, bluff-body flame holding, shock–flame interaction).
1. Condensed-phase reactive burn (
src/simulation/m_reactive_burn.fpp,src/simulation/m_pressure_relaxation.fpp)A programmed, pressure-driven burn on the multi-fluid model that converts a reactant fluid to a product fluid, releasing chemical energy — a lightweight route to detonation/deflagration without finite-rate kinetics.
It runs on both the 5-equation (pressure-equilibrium) and 6-equation (Saurel/Pelanti phasic-energy) models. For
model_eqns = 3the reactant→productqvrelease is carried through the phasic-pressure relaxation:s_equilibrate_pressure/s_correct_internal_energiesnow thread the formation energyalpha_rho_k*qv_k, making the phasic and mixture pressures qv-consistent withs_compute_pressure(a nonzeroqvotherwise inflated the phasic pressure byrho_k*qv_k/gamma_k). Forqv = 0the new terms are exactly zero, so every existing 6-equation golden is bit-identical. An optional Arrhenius term gives temperature-driven ignition instead of a pure pressure switch.Parameters are grouped in a
rburnderived type (followingchem_params):reactive_burnrburn%krburn%pignrburn%prefrburn%nrburn%ta0= pure pressure-driven,>0multiplies the rate byexp(-rburn%ta/T)withTthe reactant phasic temperaturePreconditions are checked in
m_checker/case_validator: exactly two fluids sharing the stiffened-gas EOS with an energy-releasingqvgap;model_eqnsin {2, 3};rburn%k > 0,rburn%pignset,rburn%pref > 0,rburn%n >= 0(each defaults to a sentinel, so an unset value would otherwise silently produce spurious ignition, a NaN viadrive**n, or a backward reaction);rburn%ta >= 0, andfluid_pp(1)%cv > 0whenrburn%ta > 0so the Arrhenius factor is not silently disabled. The source is hardened against negative pressure.2. IBM burning surfaces (
src/simulation/m_ibm.fpp)Immersed boundaries can now act as burning/transpiring surfaces: wall-normal fuel injection (
patch_ib%v_blow,patch_ib%inj_species) with an optional pressure-coupled burn rate following Vieille's law,v_blow*(p/p_ref)^n(patch_ib%burn_rate_exp,patch_ib%burn_rate_pref) for internal-ballistics feedback (e.g. a regressing solid-propellant grain).inj_speciesis validated againstnum_species(it indexes the image-point mass fractions, so an out-of-range value would be an out-of-bounds write).3. IBM + chemistry hang fix
Interpolating consistent species mass fractions and energy at immersed-boundary image points (threads
Ys_IPthroughs_interpolate_image_point), which fixes a hang when IBM and chemistry were enabled together.4. Reactive shock-flame hardcoded IC (
hcid 275, 2D)A sinusoidal burned|fresh flame interface for Richtmyer–Meshkov shock–flame roll-up studies.
Type of change
Testing
Built and tested on
wingtip-gpu3(NVHPC 25.11, 4×A100), rebased on currentmaster:./mfc.sh precheckpasses (format, spelling, lint, doc/param references, example validation).model_eqns = 3), Arrhenius (rburn%ta > 0), and a 2-rank MPI variant. The MPI case exists because every other reactive golden runs single-rank, which cannot catch a parameter that fails to reach non-root ranks; it was verified to fail without therburnbroadcast and pass with it.model_eqns = 3goldens pass unchanged on CPU, OpenMP, and OpenACC, confirming theqvthreading is inert forqv = 0.ibm_burning_grain/ibm_flameholderExample goldens cover surface injection; a new2D -> IBM -> Vieille Burn Rategolden covers the pressure-coupled burn rate (burn_rate_exp > 0), which both examples otherwise leave off.1e-12and was not portable across recompiles/backends).1D_propellant_flameis excluded from the Example suite: its finite-rate flame front drifts past the1e-3Example tolerance across compilers (nvhpc 25.11 vs 24.3/GNU/Intel/CCE/AMD), the same non-portability class as the flame examples already skipped there. The case itself remains as a runnable example.New example cases:
1D_reactive_burn,1D_propellant_flame,2D_hybrid_slab,2D_ibm_burning_grain,2D_ibm_flameholder,2D_shock_flame.Checklist
GPU changes (expand if you modified
src/simulation/)