SOTA: integrating large, stiff soot/PAH mechanisms efficiently

Reference notes from the CRECK “cold-band” CVODE investigation (PR #435, closed unmerged; its standalone Cantera bench scripts/creck_cold_band_bench.py lives on that branch). Captures what we found mechanistically, what’s actionable now, and what the literature offers for going further. Not a plan — a map of options, so the next person doesn’t have to re-derive or re-search this.

Status (2026-09-20)

• The opt-in rtol / atol TubeFurnace YAML properties described below shipped with PR #437.

• The first real model to need them is CH5RVP PH1 (PR #423): 300 K unheated legs on both sides of the heated zone. At the 1e-12 default the nominal case was not reproducible (crash locally; two CI runs of the same commit differing by 2x in outlet C2H4). atol: 1.0e-14 fixed that and cut the CI solve from 243 s to 42 s; soot yield and H2 then converge to <0.3 % across rtol 1e-6..1e-8 and a dense-Newton check, while outlet C2H2/C2H4 still scatter ~15 % between settings (decided in the exit-ramp quench).

• The full explanation of the cold-leg cost, and the options to skip the frozen chemistry there, are in issue #442. Even at atol 1e-14 the two legs are 55 % of the CI march and 85 % of a local one. The team is leaning towards dynamic adaptive chemistry (below).

What we actually found

Marching a closed parcel from torch-outlet conditions (~1777 K) through a wall-cooled tube toward 300 K on the full CRECK 452-species soot mechanism, CVODE’s cost is dominated by a Newton-failure storm that gets worse, not better, as the gas keeps cooling — not confined to a narrow band as originally hypothesized. In the severe case it doesn’t just get slow, it crashes: CVODE eventually accepts a state with negative moles for a few steps before a fatal, unrecoverable error.

The mechanism, traced to the solver’s own math: CVODE weighs convergence per state component as

w_i = 1 / (rtol·|y_i| + atol_i)

(SUNDIALS CVODE Mathematical Considerations). When |y_i| is small, atol alone sets the weight — rtol does nothing for trace species. That’s already the textbook reason atol matters here. What makes it worse for IdealGasConstPressureMoleReactor specifically: its species state is absolute moles (kmol), an extensive quantity that scales with however much total gas happens to be in the reactor — confirmed directly:

r.component_index('CH4')            # -> 14
r.get_state()[14]                   # -> 2.317e-06  (matches X(CH4) * n_total, not X(CH4) itself)

IdealGasConstPressureReactor’s species state is mass fraction Y_k — an intensive, self-normalized [0,1] quantity, independent of reactor size. So the same numerical atol=1e-12 means two very different things for the two reactor families: “resolve to within 1e-12 kmol absolute” (which, for a trace species whose entire absolute amount in a small reactor is already ~1e-12 to 1e-13 kmol, is essentially no resolution) versus “resolve to within 1e-12 of total mass” (a scale that never drifts). This is very likely why IdealGasConstPressureReactor showed 48x fewer Newton failures and ran 6.5x faster on the same severe case in the bench’s test 7, with zero tuning.

SUNDIALS’ C API supports a per-component atol vector (SVtolerances); Cantera’s Python ReactorNet.atol does not — it rejects a list (TypeError: must be real number, not list, checked directly against Cantera 3.2.0). So true per-species scaling isn’t reachable from Python today without a Cantera-side change.

Near-term lever: scale atol to the reactor’s own total moles

Actionable now, cheap to test, no new dependencies. Since the mole reactor’s atol is absolute-extensive, a fixed constant is only well-scaled for whatever total mole count it happened to be tuned against. At reactor-build time (build_closed_carrier_net already has density, volume, mean_molecular_weight), compute

n_total = density * volume / mean_molecular_weight
atol = n_total * relative_floor        # e.g. relative_floor ~ 1e-10

instead of hardcoding 1e-12. This makes the tolerance self-scaling across geometries and mass flows instead of requiring per-model manual tuning (which is what the opt-in atol/rtol TubeFurnace YAML properties, now in PR #437, require the user to do by hand — PH1 sets atol: 1.0e-14). Not yet built or tested — needs the same evidence bar as everything else in this investigation: checked against both the bench’s severe case and a real, non-pathological case (tube_furnace/BASELINE) before it’s believed.

SOTA landscape for large-scale stiff mechanisms

None of these validate literally switching between a small, hand-picked set of mechanisms at fixed temperature/composition thresholds — that risks discontinuities and mass-conservation artifacts right at the switch boundary, the same class of problem this investigation has been chasing, just relocated. The literature’s answers are smoother and more automatic versions of the same idea.

Dynamic Adaptive Chemistry (DAC / TDAC)

The closest real match to “one mechanism per regime” — except continuous and automatic. At every computational state, a graph-based reduction (DRG / DRGEP) builds a local skeletal mechanism sized to what’s actually active there. Reported 15-70x speedups on 578-1099 species mechanisms (larger than CRECK’s 452). TDAC adds ISAT tabulation on top for more. Substantial build: nothing in Cantera does this automatically; it needs a DRGEP reduction layer wrapping every chemistry call.

Sketch for the Lagrangian parcel march (issue #442). The natural granularity is the checkpoint: the march already stops 200 times per tube to record the state and refresh the wall terms. A “light DAC” would, at each checkpoint:

1. Evaluate the net rates of progress of all reactions at the current state (Kinetics.net_rates_of_progress) and build the DRGEP graph from them: for each species pair, the fraction of A’s production/consumption that runs through reactions involving B. Path-flux (PFA) is the usual refinement if two-step couplings matter.

2. Starting from a target set (H2, C2H2, the BIN soot classes, CO — the KPI species), keep every species whose DRGEP coefficient exceeds a threshold (typically 1e-3 to 1e-2) and every reaction whose participants are all kept.

3. Zero the multipliers of the discarded reactions (Kinetics.set_multiplier(0, i)), restore the kept ones to 1, and ReactorNet.reinitialize() before the next advance. The species state vector keeps its full length; discarded species simply stop changing, so nothing is lost at a switch and elements are conserved exactly.

Expected effect: in a 300 K leg the local skeleton collapses to a few dozen radical-recombination reactions and the Newton iteration has nothing left to fail on; in the 1673 K plateau most of the mechanism stays active. The threshold, the target set and the checkpoint count are per-model, opt-in YAML properties, never a runtime heuristic in the shared engine (see the reverted attempts in PR #435 for why).

Acceptance bar: tube_furnace/BASELINE and every model that does not opt in byte-identical; for a model that opts in, KPI species within the tolerance scatter documented for it (PH1: soot yield and H2 within 0.3 %, C2H2/C2H4 within ~15 %); the leg vs heated-zone time split reported before and after. Cantera provides the multipliers and the rates; the DRGEP graph (a few hundred lines) has to be written, or borrowed from an offline skeletal-reduction tool such as pyMARS and applied per checkpoint.

• The use of dynamic adaptive chemistry in combustion simulation of gasoline surrogate fuels

• A dynamic adaptive chemistry scheme for reactive flow computations

• Assessment of On-the-Fly Chemistry Reduction and Tabulation Approaches (TDAC)

Quasi-Steady-State Approximation (QSSA)

Instead of switching mechanisms, algebraically collapse the fastest species (typically short-lived radicals) into zero-derivative constraints, removing exactly the stiffest eigenvalues from the Jacobian rather than working around them with tolerances. Well-precedented; tooling exists (Cantera- compatible QSS reduction scripts, CHEMKIN’s QSS preprocessor). A mechanism-file-level transformation of CRECK, not a solver setting.

• Quasi-Steady-State Assumptions (QSSA) in Kinetics

• Autoignition Problem: GQL versus QSSA Combined with DRG

In Situ Adaptive Tabulation (ISAT)

Cache/interpolate the chemistry mapping in composition space instead of re-integrating every time; up to ~1000x reported in turbulent CFD. Less obviously useful for a single 1D Lagrangian march that doesn’t revisit similar states the way a 3D CFD grid with millions of cells does — the benefit comes from reuse, which one march has less of.

• In situ adaptive tabulation — Wikipedia

• Pope, Computationally efficient implementation of combustion chemistry using ISAT (1997)

Flamelet Generated Manifolds (FGM) with soot progress variables

Built specifically for a “fast gas chemistry + slow soot” structure: dual-scale FGM uses two progress variables, one for major species and a separate one for PAH/soot evolution, precomputed offline into a lookup table. Strongest conceptual match to “mixture of mechanisms per regime” for this exact combustion+soot case, but it’s a paradigm shift — table lookup replacing finite-rate integration entirely, not a tweak to the current approach.

• Systematic definition of progress variables for Flamelet Generated Manifolds

• Soot formation in turbulent non-premixed jet flames based on a dual-scale FGM

Operator splitting: soot population balance vs. gas-phase chemistry

The most incremental-feeling option: decouple the (slow) soot population balance from the (fast) deterministic gas-phase chemistry solve via Strang-type splitting, each on its own appropriate timescale/solver. Directly operationalizes “leverage the fast/slow split” without replacing the modeling paradigm. Specific literature exists on exactly this pairing.

• Coupling a stochastic soot population balance to gas-phase chemistry using operator splitting

Domain-specific prior art: OpenSMOKE++

Built by the same group that publishes CRECK (Cuoci, Frassoldati, Faravelli, Ranzi) — a purpose-built, object-oriented C++ framework specifically for large detailed mechanisms with PAH/soot chemistry. The mechanism’s own authors evidently found generic tools insufficient for exactly this problem class, which is a strong signal this is a known, real limitation and not specific to how Bloc uses Cantera.

What’s actually in it:

• Kinetic preprocessor — CHEMKIN-format compatible, so CRECK mechanisms load natively.

• Ideal reactors — batch, plug-flow, CSTR/jet-stirred, shock-tube, rapid compression machine.

• 1D laminar flames — freely-propagating, burner-stabilized, counterflow diffusion; plus a flamelet/lookup-table generator (their own FGM-style tool, see above).

• Soot module — PAHs with ≥20 carbon atoms plus an aerosol/particle population tracked via a sectional method, i.e. soot is not just more finite-rate kinetic species piled onto the same stiff ODE the way CRECK-in-Cantera currently treats it — it’s a separate representation designed for the slow end of the timescale split this whole investigation has been circling.

• Own stiff ODE/DAE numerics — a dedicated library (BzzMath, developed and maintained by the same authors), plus, per the suite’s own docs, “coupling to a wide range of external ODE, DAE, and NLS solvers” and dense/sparse direct and iterative linear solvers — not simply “Cantera’s CVODE by another name.”

• OpenSMOKE++4OpenFOAM / laminarSMOKE — couples the same chemistry core into OpenFOAM for full 2D/3D reacting-flow CFD, not just 0D/1D.

• GPP (Graphical Post-Processor) — sensitivity analysis, rate-of- production analysis, reaction-path diagrams via GraphViz.

Practical caveat before treating this as a real option: the project’s own site states it is “completely free for Academic use” with registration required. That phrasing does not extend to commercial use, and nothing found states commercial terms — this would need to be checked directly with the CRECK Modeling Lab / Politecnico di Milano before assuming Spark could adopt it. Beyond licensing, this is an architecture decision, not a small change: a different language ecosystem (C++, not the Cantera/Python stack bloc is built on), different reactor/flame APIs, and a soot representation (sectional method) that isn’t a drop-in replacement for CRECK’s finite-rate PAH species — it would mean re-deriving soot output against the sectional model’s own moments, not just swapping a solver call.

• OpenSMOKE++: An object-oriented framework for reactive systems with detailed kinetic mechanisms

• Kinetic Modeling of Soot Formation in Turbulent Nonpremixed Flames (Cuoci, Frassoldati, Faravelli, Ranzi)

• OpenSMOKE++ Suite (components, CHEMKIN compatibility, GPP)

• OpenSMOKE homepage (academic-use license statement)

• laminarSMOKE — OpenSMOKE++4OpenFOAM on GitHub

• BzzOde: a C++ class for stiff/non-stiff ODE systems

Complementary: analytical/sparse Jacobians (pyJac)

A performance angle rather than a stiffness-avoidance strategy: generate fully analytical, sparse Jacobians instead of finite-differencing them. Cantera’s own sparse AdaptivePreconditioner already captures much of this benefit for us — the bench’s test 4 confirmed removing it is 47x worse, not better — so this is lower incremental value for Bloc specifically than for a from-scratch solver.

• pyJac: Analytical Jacobian generator for chemical kinetics

Effort vs. payoff

Option	Effort	Notes
Auto-scaled atol (this doc, above)	Cheap, same-day testable	No new dependencies; extends PR #435’s opt-in knob
QSSA-reduced CRECK variant	Medium — mechanism-level project	Removes stiffness at the source; tooling exists
Operator-split soot vs. gas-phase	Medium-large — numerical-methods project	Directly matches the fast/slow structure asked about
DAC / TDAC	Large — new infrastructure	Best-documented speedups on comparable mechanism sizes
FGM / tabulated soot chemistry	Large — paradigm shift	Closest conceptual match to “regime-based mechanisms”
Switch to OpenSMOKE++	Large — architecture decision	Prior art from CRECK’s own authors

None of these should be started speculatively. Each needs the same evidence bar this investigation has held throughout: checked against both a severe reproduction case and a real, currently-working case (e.g. tube_furnace/BASELINE) before anything lands. Earlier in this same investigation, a plausible-sounding runtime fix that wasn’t checked against a working case cost real regressions (non-deterministic outputs, and a 3.7x-11x slowdown on BASELINE) before it was reverted — see the PR #435 discussion for the full account.

Full reference list

• SUNDIALS CVODE Mathematical Considerations (WRMS norm, atol/rtol weights): https://sundials.readthedocs.io/en/latest/cvode/Mathematics_link.html

• Dynamic adaptive chemistry, gasoline surrogates: https://www.sciencedirect.com/science/article/abs/pii/S0010218009000819

• Dynamic adaptive chemistry scheme for reactive flow computations: https://sciencedirect.com/science/article/pii/S1540748908000941

• Tabulated dynamic adaptive chemistry (TDAC) for MILD combustion: https://pubs.acs.org/doi/10.1021/acs.energyfuels.8b01001

• In situ adaptive tabulation (ISAT) — Wikipedia: https://en.wikipedia.org/wiki/In_situ_adaptive_tabulation

• Pope (1997), ISAT original paper: https://tcg.mae.cornell.edu/pubs/Pope_CTM_97.pdf

• QSSA overview: https://www.emergentmind.com/topics/quasi-steady-state-assumptions-qssa

• QSSA combined with DRG for autoignition: https://doi.org/10.3390/modelling4040027

• Progress variables for Flamelet Generated Manifolds: https://www.sciencedirect.com/science/article/abs/pii/S0010218011003087

• Dual-scale FGM for soot in turbulent jet flames: https://pubs.aip.org/aip/pof/article-abstract/38/7/075122/3397552/Soot-formation-in-turbulent-non-premixed-jet

• Operator splitting for soot population balance + gas-phase chemistry: https://www.researchgate.net/publication/222526279_Coupling_a_stochastic_soot_population_balance_to_gas-phase_chemistry_using_operator_splitting

• OpenSMOKE++ framework: https://www.sciencedirect.com/science/article/abs/pii/S0010465515000715

• Kinetic modeling of soot formation in turbulent nonpremixed flames (CRECK group): https://journals.sagepub.com/doi/abs/10.1089/ees.2007.0193

• pyJac analytical Jacobian generator: https://arxiv.org/abs/1605.03262