Operator-split integration for the condensed-phase reactive burn - #1831
Operator-split integration for the condensed-phase reactive burn#1831sbryngelson wants to merge 2 commits into
Conversation
…ubsteps > 0 Added to the flow RHS, the burn ties the reaction time scale to the acoustic CFL: a fast burn forces a smaller step for the whole simulation. With rburn%substeps > 0 the flow is frozen and the ODE is integrated over the step in equal sub-steps, re-evaluating the mixture pressure from the frozen internal energy each one so the rate feels the coefficients moving as reactant becomes product. rburn%substeps = 0, the default, keeps the source in the RHS unchanged. Completing a sub-step hands over the reactant's remaining mass along with the last of its volume fraction; capping the two separately strands mass at zero volume, which the EOS then divides by. The rate law is now stated once and called by both integrators.
There was a problem hiding this comment.
Warning
Copilot couldn't run its full agentic review because it didn't start before the timeout. Make sure your repository has a runner available, or add a copilot-code-review.yml file specifying one with the runs-on attribute. See the docs for more details.
Pull request overview
Adds an operator-split path for condensed-phase reactive burn by introducing rburn%substeps, allowing the burn ODE to be integrated in post-flow substeps so fast chemistry doesn’t globally constrain the acoustic CFL time step.
Changes:
- Register/validate new parameter
rburn%substepsand ensure its MPI broadcast uses the correct MPI datatype (integer vs real). - Integrate reactive burn via operator splitting after the flow update when
rburn%substeps > 0, leaving the existing RHS source path forsubsteps == 0. - Add a new 2-rank golden test case to exercise the new substepping path and catch broadcast-kind issues.
Reviewed changes
Copilot reviewed 11 out of 11 changed files in this pull request and generated 3 comments.
Show a summary per file
| File | Description |
|---|---|
| toolchain/mfc/test/cases.py | Adds a new reactive-burn test case with rburn%substeps=10 on 2 MPI ranks. |
| toolchain/mfc/params/generators/fortran_gen.py | Updates generated MPI_BCAST code to pick MPI type per rburn member (supports new INT member). |
| toolchain/mfc/params/definitions.py | Registers rburn%substeps as an INT reactive-burn parameter. |
| toolchain/mfc/case_validator.py | Validates rburn%substeps >= 0 with an explanatory error message. |
| src/simulation/m_time_steppers.fpp | Calls s_reactive_burn_substep after the flow update when substeps > 0. |
| src/simulation/m_rhs.fpp | Skips adding burn source to RHS when substeps > 0 (operator-split path instead). |
| src/simulation/m_reactive_burn.fpp | Factors out s_burn_rate and implements operator-split substep integrator. |
| src/common/m_global_parameters_common.fpp | Initializes rburn%substeps default to 0. |
| src/common/m_derived_types.fpp | Extends reactive burn parameter type with integer substeps. |
| tests/86893F55/golden.txt | Adds new golden outputs for the new substepping test case. |
| tests/86893F55/golden-metadata.txt | Adds metadata for the new golden generation. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| ! Optional Arrhenius dependence on the reactant phasic temperature from the stiffened-gas | ||
| ! EOS. rburn%ta = 0, the default, leaves the pure pressure-driven rate unchanged. | ||
| if (rburn%ta > 0._wp) then | ||
| rate = rate*exp(-rburn%ta/f_sg_thermal(pres, alpha_rho_react/alpha_react, isentrope_n(1), isentrope_B(1), cvs(1))) |
There was a problem hiding this comment.
Fixed in ed51f33 — the divisor is now max(alpha_react, sgm_eps), matching the form the Riemann solvers already use for phasic densities.
Worth noting it is reachable: completing a sub-step sets the reactant volume fraction to exactly zero, and the lambda < 1 test only covers the same iteration while the volume fractions still sum to one, which the six-equation model does not guarantee.
| alpha_rho(i) = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(x, y, z) | ||
| alpha(i) = q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(x, y, z) | ||
| end do | ||
| rho = alpha_rho(1) + alpha_rho(2) |
There was a problem hiding this comment.
The second alternative already holds: case_validator.py:1932 enforces reactive_burn requires num_fluids = 2 (reactant then product), so fluids 1 and 2 are the whole mixture here. s_compute_reactive_burn computes rho the same way for the same reason. Leaving as is.
| rsub = self.get("rburn%substeps") | ||
| self.prohibit( | ||
| self._is_numeric(rsub) and rsub < 0, | ||
| "reactive_burn requires rburn%substeps >= 0 (operator-split sub-steps per time step; 0 adds the source to the flow RHS)", | ||
| ) |
There was a problem hiding this comment.
Leaving this lower-bounded only, for consistency with the sibling parameter: chem_params%reaction_substeps drives the identical operator-split pattern for chemistry and is likewise validated as >= 0 with no ceiling. Adding a bound to one and not the other would be the surprising choice, and a defensible upper limit is case-dependent rather than universal.
Completing a sub-step drives the reactant volume fraction to exactly zero, and the lambda >= 1 test only covers the same iteration if the volume fractions still sum to one, which the six-equation model does not guarantee. Use the max(alpha, sgm_eps) form the Riemann solvers already use for this.
Lines of Code
|
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1831 +/- ##
==========================================
+ Coverage 62.30% 62.35% +0.04%
==========================================
Files 84 84
Lines 21583 21622 +39
Branches 3195 3199 +4
==========================================
+ Hits 13448 13483 +35
- Misses 5937 5939 +2
- Partials 2198 2200 +2 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Added to the flow RHS, the condensed-phase burn ties the reaction time scale to the acoustic CFL: a fast burn forces a smaller step for the entire simulation, even where nothing is reacting.
This adds
rburn%substeps. Above zero, the flow is frozen and the burn ODE is integrated across the step in equal sub-steps, with the mixture pressure re-evaluated from the frozen internal energy each sub-step so the rate feels its own feedback as reactant becomes product.rburn%substeps = 0is the default and leaves the source in the flow RHS exactly as before.Two details worth review:
s_burn_ratenow states the pressure drive, the(1 - lambda)factor and the optional Arrhenius term once; both integrators call it. Previously the new path would have carried a second copy that had to stay in step forever.Testing
1D -> Reactive Burn -> Condensed Programmed Detonation -> substeps(86893F55),substeps = 10on 2 ranks. Two ranks deliberately:substepsis the only integer among therburnmembers, and the broadcast kind is now registry-driven, so a wrong kind leaves rank 1 sub-stepping a garbage count — invisible to a single-rank golden.