Operator-split α-QSS sub-stepped chemistry integrator with per-rank adaptive substepping#1647
Open
sbryngelson wants to merge 13 commits into
Open
Operator-split α-QSS sub-stepped chemistry integrator with per-rank adaptive substepping#1647sbryngelson wants to merge 13 commits into
sbryngelson wants to merge 13 commits into
Conversation
… with per-rank adaptive substepping
…a-QSS substepping
sbryngelson
marked this pull request as ready for review
July 16, 2026 20:04
Contributor
There was a problem hiding this comment.
Pull request overview
Note
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-setup-steps.yml file specifying one with the runs-on attribute. See the docs for more details.
Adds an operator-split α‑QSS chemistry integrator with fixed sub-stepping and optional per-rank adaptive sub-step selection, wiring it into the simulation timestep and exposing new configuration parameters and examples.
Changes:
- Introduces operator-split reaction integration (
s_chemistry_reaction_substep) and disables reaction RHS coupling when sub-stepping is enabled. - Adds new chemistry parameters (
reaction_substeps,adap_substeps,reaction_substeps_max) across definitions, defaults, and input checking. - Expands tests/examples coverage and updates dependency pinning and goldens.
Reviewed changes
Copilot reviewed 18 out of 18 changed files in this pull request and generated 7 comments.
Show a summary per file
| File | Description |
|---|---|
| toolchain/pyproject.toml | Switches pyrometheus dependency to a pinned release version. |
| toolchain/mfc/test/cases.py | Adds new example cases and a chemistry sub-stepping test case. |
| toolchain/mfc/params/definitions.py | Registers new chemistry parameters for the toolchain. |
| tests/5D59DDF0/golden.txt | Adds/updates golden outputs for a test case. |
| tests/5D59DDF0/golden-metadata.txt | Adds/updates metadata for the golden generation run. |
| src/simulation/m_time_steppers.fpp | Calls operator-split chemistry sub-stepper after the flow update. |
| src/simulation/m_rhs.fpp | Skips reaction RHS contribution when operator-split sub-stepping is enabled. |
| src/simulation/m_global_parameters.fpp | Adds defaults for new chemistry parameters (simulation). |
| src/simulation/m_checker.fpp | Validates new chemistry parameters and combinations. |
| src/pre_process/m_global_parameters.fpp | Adds defaults for new chemistry parameters (pre-process). |
| src/post_process/m_global_parameters.fpp | Adds defaults for new chemistry parameters (post-process). |
| src/common/m_derived_types.fpp | Extends chemistry_parameters with new fields and docstrings. |
| src/common/m_chemistry.fpp | Implements α‑QSS operator-split sub-stepped reactor with adaptive per-rank sub-stepping. |
| examples/2D_reactive_shock_bubble/h2o2_xe.yaml | Adds a new Cantera mechanism YAML for an example case. |
| examples/2D_reactive_shock_bubble/case.py | Adds a new reactive shock-bubble example using sub-stepped chemistry. |
| examples/2D_detonation_cell/case.py | Adds a cellular detonation example using adaptive sub-stepping. |
| examples/1D_reactive_shocktube_adaptive/case.py | Adds a 1D adaptive sub-stepping shocktube example. |
| docs/documentation/case.md | Documents reaction_substeps (but currently describes the wrong algorithm). |
Comment on lines
+171
to
+175
| #:if not MFC_CASE_OPTIMIZATION and USING_AMD | ||
| real(wp), dimension(10) :: Ys, cdot, ddot, y0, prod0, Lloss, alp | ||
| #:else | ||
| real(wp), dimension(num_species) :: Ys, cdot, ddot, y0, prod0, Lloss, alp | ||
| #:endif |
Comment on lines
+251
to
+266
| Ysum = 0._wp | ||
| $:GPU_LOOP(parallelism='[seq]') | ||
| do eqn = 1, num_species | ||
| wr = molecular_weights(eqn)/rho | ||
| prod_p = wr*cdot(eqn) | ||
| loss_p = wr*ddot(eqn) | ||
| Lbar = 0.5_wp*(Lloss(eqn) + loss_p/max(Ys(eqn), y_floor)) | ||
| pbar = alp(eqn)*prod_p + (1._wp - alp(eqn))*prod0(eqn) | ||
| Ys(eqn) = y0(eqn) + dt_sub*(pbar - Lbar*y0(eqn))/(1._wp + alp(eqn)*dt_sub*Lbar) | ||
| if (Ys(eqn) < 0._wp) Ys(eqn) = 0._wp | ||
| Ysum = Ysum + Ys(eqn) | ||
| end do | ||
| $:GPU_LOOP(parallelism='[seq]') | ||
| do eqn = 1, num_species | ||
| Ys(eqn) = Ys(eqn)/Ysum | ||
| end do |
| | `chem_params%%reactions` | Logical | Enable chemical reactions | | ||
| | `chem_params%%gamma_method` | Integer | Methodology for calculating the heat capacity ratio | | ||
| | `chem_params%%transport_model` | Integer | Methodology for calculating the diffusion coefficients | | ||
| | `chem_params%%reaction_substeps` | Integer | Sub-steps for operator-split reaction integration (0 = off) | |
| | `cantera_file` | String | Cantera-format mechanism file (e.g., .yaml) | | ||
|
|
||
| - `chem_params%%transport_model` specifies the methodology for calculating diffusion coefficients and other transport properties, `1` for mixture-average, `2` for Unity-Lewis | ||
| - `chem_params%%reaction_substeps` controls how the reaction source is integrated. With `0` (default) the net production rates are added to the flow right-hand side and advanced by the flow time stepper (fine for hydrogen). With a value `> 0`, the reaction is instead integrated by operator splitting after each flow update: every cell's constant-density, constant-internal-energy reactor is advanced over the timestep with that many forward-Euler sub-steps. This decouples the (often much faster) chemical timescale from the flow timestep and is required for stiff mechanisms — e.g. hydrocarbons such as GRI-Mech methane, which otherwise diverge on the first step |
| "cantera>=3.1.0", | ||
| #"pyrometheus == 1.0.5", | ||
| "pyrometheus @ git+https://github.com/pyrometheus/pyrometheus.git", | ||
| "pyrometheus == 1.1.1", |
|
|
||
| Invocation: test --generate --only 5D59DDF0 -j 8 | ||
| Lock: mpi=Yes & gpu=Acc & debug=Yes & reldebug=No & gcov=No & unified=No & single=No & mixed=No & fastmath=No | ||
| Git: 30a411e98a697d4b7c12a342d9e26f3b46d852ca on chemistry (dirty) |
Comment on lines
+389
to
+391
| chem_params%reaction_substeps = 0 | ||
| chem_params%adap_substeps = .false. | ||
| chem_params%reaction_substeps_max = 0 |
…lpha-QSS + adaptive controls
…istic provenance)
Member
Author
|
Thanks for the review — addressed in Fixed
Respectfully keeping as-is (with rationale)
|
…er exponent, document stiffness metric
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #1647 +/- ##
==========================================
- Coverage 60.95% 59.61% -1.35%
==========================================
Files 83 83
Lines 20003 21194 +1191
Branches 2983 3137 +154
==========================================
+ Hits 12193 12634 +441
- Misses 5782 6452 +670
- Partials 2028 2108 +80 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…d reactor seed q_T_sf holds the pre-update RK-stage temperature, so a just-shock-heated cell was probed at its stale (cool) temperature, under-sizing the adaptive nsub exactly at an ignition front. Recompute T from the fresh post-advection internal energy before the stiffness probe and the first predictor sub-step. Reactor golden 5D59DDF0 unchanged (at-rest reactor: same energy -> same T).
Address review nits: the reactor header now notes the adaptive nsub range; the example comments say one nsub per rank (sized by the rank's stiffest cell) instead of implying spatially-local sub-cycling; stiff_target=0.5 is documented as an uncalibrated accuracy/cost knob (never stability); the corrector's alpha reuse is noted as an intentional CHEMEQ2 simplification; and the header records the first-order Lie-Trotter split and the sum-to-one (not element-conserving) renorm.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Adds an operator-split integrator for the chemistry reaction source, run as a fixed number of sub-steps per flow step (
chem_params%reaction_substeps). The sub-step uses an α-QSS (Mott/CHEMEQ2 predictor–corrector) update on a constant-(ρ,e) reactor — matrix-free, Jacobian-free, and GPU-friendly (deterministic per-cell work, no thread divergence). This keeps stiff mechanisms stable where an explicit source in the flow RHS diverges (e.g. GRI-Mech methane).On top of that, an optional per-rank adaptive sub-step count (
chem_params%adap_substeps): each rank sizes itsnsubfrom a local max-reduction of a stiffness metric over its own cells (floorreaction_substeps, ceilingreaction_substeps_max). Zero MPI collectives — bottleneck-equivalent to a global reduction since the slowest rank dominates anyway, but safe at scale. Default off, so existing behavior is unchanged.Dependency
The integrator calls
get_creation_destruction_rates(a single-pass routine filling creation + destruction in one concentration/rate evaluation, halving the rate-of-progress work per sub-step), added in pyrometheus/pyrometheus#116 and released in pyrometheus 1.1.1.toolchain/pyproject.tomlpinspyrometheus == 1.1.1— verified end-to-end that MFC installs it from PyPI and builds against it.Validation
5D59DDF0) shifts (forward-Euler → α-QSS), max ~0.5% and confined to trace radicals; allreactions=F/substeps=0goldens unchanged. Golden regenerated for α-QSS.nsubverified to track stiffness across [floor, ceiling] (floor in quiescent gas, ceiling only at the reaction front).substeps=0stability limit);substeps=0at that dt NaNs.Files
src/common/m_chemistry.fpp— α-QSS integrator + adaptive probem_derived_types.fpp,toolchain/mfc/params/definitions.py,m_global_parameters.fpp(×3),src/simulation/m_checker.fppsrc/simulation/m_rhs.fpp,m_time_steppers.fpp— operator-split call site2D_reactive_shock_bubble,1D_reactive_shocktube_adaptive,2D_detonation_cell(the latter two demonstrate the automatedcfl_const_dt+ adaptive pipeline)