Skip to content

Combustion: condensed-phase reactive burn (5-/6-equation, Arrhenius kinetics) + IBM burning surfaces#1670

Open
sbryngelson wants to merge 26 commits into
MFlowCode:masterfrom
sbryngelson:reactive-burn-ibm-combustion
Open

Combustion: condensed-phase reactive burn (5-/6-equation, Arrhenius kinetics) + IBM burning surfaces#1670
sbryngelson wants to merge 26 commits into
MFlowCode:masterfrom
sbryngelson:reactive-burn-ibm-combustion

Conversation

@sbryngelson

@sbryngelson sbryngelson commented Jul 23, 2026

Copy link
Copy Markdown
Member

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 = 3 the reactant→product qv release is carried through the phasic-pressure relaxation: s_equilibrate_pressure / s_correct_internal_energies now thread the formation energy alpha_rho_k*qv_k, making the phasic and mixture pressures qv-consistent with s_compute_pressure (a nonzero qv otherwise inflated the phasic pressure by rho_k*qv_k/gamma_k). For qv = 0 the 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 rburn derived type (following chem_params):

Parameter Meaning
reactive_burn Enable the burn
rburn%k Rate coefficient [1/s]
rburn%pign Ignition-pressure threshold [Pa]
rburn%pref Reference pressure for the pressure drive [Pa]
rburn%n Pressure-drive exponent
rburn%ta Activation temperature [K]; 0 = pure pressure-driven, >0 multiplies the rate by exp(-rburn%ta/T) with T the reactant phasic temperature

Preconditions are checked in m_checker / case_validator: exactly two fluids sharing the stiffened-gas EOS with an energy-releasing qv gap; model_eqns in {2, 3}; rburn%k > 0, rburn%pign set, rburn%pref > 0, rburn%n >= 0 (each defaults to a sentinel, so an unset value would otherwise silently produce spurious ignition, a NaN via drive**n, or a backward reaction); rburn%ta >= 0, and fluid_pp(1)%cv > 0 when rburn%ta > 0 so 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_species is validated against num_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_IP through s_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

  • New feature
  • Bug fix (IBM+chemistry hang; negative-pressure hardening)

Testing

Built and tested on wingtip-gpu3 (NVHPC 25.11, 4×A100), rebased on current master:

Backend Build (pre_process / simulation / post_process) Affected tests
CPU 20/20 pass
OpenMP target offload (NVIDIA) pass
OpenACC (NVIDIA) pass
  • ./mfc.sh precheck passes (format, spelling, lint, doc/param references, example validation).
  • Reactive-burn goldens: 5-equation, 6-equation (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 the rburn broadcast and pass with it.
  • 6-equation relaxation path: the 13 existing model_eqns = 3 goldens pass unchanged on CPU, OpenMP, and OpenACC, confirming the qv threading is inert for qv = 0.
  • IBM burning: the ibm_burning_grain / ibm_flameholder Example goldens cover surface injection; a new 2D -> IBM -> Vieille Burn Rate golden covers the pressure-coupled burn rate (burn_rate_exp > 0), which both examples otherwise leave off.
  • The reactive goldens were hardened to a small uniform advection velocity so the momentum output is well-resolved rather than near-zero roundoff (the prior golden compared ~0 momentum at 1e-12 and was not portable across recompiles/backends).
  • 1D_propellant_flame is excluded from the Example suite: its finite-rate flame front drifts past the 1e-3 Example 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

  • I added or updated tests for new behavior
  • I updated documentation if user-facing behavior changed
GPU changes (expand if you modified src/simulation/)
  • GPU results match CPU results
  • Tested on NVIDIA GPU or AMD GPU

Copilot AI review requested due to automatic review settings July 23, 2026 23:55
@sbryngelson

Copy link
Copy Markdown
Member Author

@claude full review

@claude

claude Bot commented Jul 23, 2026

Copy link
Copy Markdown

Claude Code is working…

I'll analyze this and get back to you.

View job run

@sbryngelson

Copy link
Copy Markdown
Member Author

@claude full review

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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_burn and 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 chemistry flag rather than present(Ys_IP). That assumes every chemistry call 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

Comment thread src/simulation/m_ibm.fpp
pres_IP = 0._wp
vel_IP = 0._wp

if (chemistry) Ys_IP = 0._wp
Comment thread src/simulation/m_checker.fpp
@claude

claude Bot commented Jul 24, 2026

Copy link
Copy Markdown

Claude Code is working…

I'll analyze this and get back to you.

View job run

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
@sbryngelson

Copy link
Copy Markdown
Member Author

Thanks — both are valid. Addressed in c1c5ada:

  1. rburn_pref divide-by-zero: reactive_burn computes drive = (p - rburn_pign)/rburn_pref, so I added @:PROHIBIT(reactive_burn .and. rburn_pref <= 0._wp, ...) to m_checker.

  2. Ys_IP optional accessed without present(): the branch in s_ibm_correct_state selects the bubbles_euler/qbmm path before the chemistry path, so Ys_IP isn't passed while if (chemistry) inside s_interpolate_image_point is still true — reachable only if chemistry and Euler bubbles are enabled together. Rather than add present() guards to the hot IBM path (the module guards optionals by feature flag, not present()), I disabled the unsupported combination: @:PROHIBIT(chemistry .and. (bubbles_euler .or. qbmm), ...) in m_checker plus the matching constraint in case_validator.py. No existing example uses that combination.

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.
@sbryngelson sbryngelson changed the title Combustion: condensed-phase reactive burn + IBM burning surfaces (fuel injection & Vieille's-law burn rate) Combustion: condensed-phase reactive burn (5-/6-equation, Arrhenius kinetics) + IBM burning surfaces Jul 24, 2026
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

codecov Bot commented Jul 25, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 90.47619% with 8 lines in your changes missing coverage. Please review.
✅ Project coverage is 60.49%. Comparing base (4796d9a) to head (9c8ea1d).
⚠️ Report is 3 commits behind head on master.

Files with missing lines Patch % Lines
src/simulation/m_ibm.fpp 87.09% 0 Missing and 4 partials ⚠️
src/simulation/m_reactive_burn.fpp 89.47% 0 Missing and 2 partials ⚠️
src/simulation/m_checker.fpp 80.00% 1 Missing ⚠️
src/simulation/m_rhs.fpp 75.00% 0 Missing and 1 partial ⚠️
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.
📢 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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Development

Successfully merging this pull request may close these issues.

2 participants