Skip to content

Chemistry split

How a network run's stiff chemistry is chosen and run each step: which driver reacts at all, when a chemistry step fires, which integrator a pipe gets, and where a step is refused. Written 2026-09-29 with the contracts of solver/chemistry_ops.py and the kernels it dispatches to; the standard is Function contracts and flow charts.

What enters. A built network whose pipes may carry a gas-phase law (py_reaction, py_reaction_multi, plus py_reaction_params for the JIT Arrhenius kernel), a washcoat surface law (py_reaction_surface with coverages in surface_state), or a catalyst mechanism (surface_mechanism). Mass fractions are MASS fractions throughout; every rate law is per mass [1/s], heat release is volumetric [W/m^3], and T [K] is recovered from the conserved energy with density and velocity frozen for the step.

Which driver reacts.

  • solver = "strang" runs T(Δt/2) C(Δt) T(Δt/2), or, with chem_every_n > 1 / chem_interval > 0, a Lie super-cycle C(N·Δt) after N transport steps. Each C refreshes γ and is cut into chem_substeps equal sub-steps.
  • An implicit rk runs C(Δt) after each transport step, super-cycled the same way (no γ refresh, no sub-division).
  • solver = "explicit" has no split: gas-phase laws are a source term, a surface law or a catalyst mechanism is refused before the run.
  • solver = "quasi_steady" has its own chemistry step (banded LSODA, linop or fixed_step; see the quasi-steady chart). It shares only the extras reorder and the coverage bound with the kernels below, and it refuses a catalyst mechanism, which it cannot integrate.

Which integrator. chemistry_step_on_pipes sends a pipe with a catalyst mechanism to its own batched LSODA whatever chem_solver says; every other pipe goes to chemistry_step_on_pipe, which reads chem_solver. Two arms fall through to LSODA when the pipe is not theirs: numbalsoda (anything but a compiled surface mechanism) and bs_brinkmeier (a pipe without the Brinkmeier dispatcher; until 2026-09-29 this arm stopped instead, and any other chemistry on that pipe silently did not run). bs_jit refuses a pipe without its parameter pack; auto picks it only for a plain Arrhenius pipe.

Refusals. Every arm commits all-or-nothing. The LSODA arms check before writing; the two Bader-Deuflhard arms write in place and are restored by their wrappers when check_and_raise refuses (also 2026-09-29). A coverage outside [0, 1] by more than _COVERAGE_VIOLATION_TOL is refused; inside it, species are floored at 0 and coverages clipped into [0, 1] at commit on every arm. Density is never checked here: chemistry does not change it, and the driver's next density check refuses a bad cell.

flowchart TD
    A[network step, driver chosen by dispatch_advance] --> B{solver}
    B -- explicit --> BX{surface law or catalyst mechanism?}
    BX -- yes --> XR[ValueError before the run]
    BX -- no --> BE[gas laws as a source term, no split]
    B -- quasi_steady --> QM{catalyst mechanism?}
    QM -- yes --> XR
    QM -- no --> QS[QS chem step: banded LSODA / linop / fixed_step]
    B -- strang --> SG{chem_every_n or chem_interval due?}
    SG -- "N = 1" --> SS["T(dt/2) C(dt) T(dt/2)"]
    SG -- "N > 1" --> SL["T(dt) each step, C(N dt) when due"]
    SS --> SB[_strang_chem_burst: refresh gamma, chem_substeps equal pieces]
    SL --> SB
    B -- implicit rk --> IG["_ImplicitRun.chemistry_supercycle: C(accumulated dt) when due"]
    SB --> P[chemistry_step_on_pipes, per pipe]
    IG --> P
    P -- surface_mechanism --> M[surface_mechanism_step_on_pipe: batched banded LSODA]
    P -- otherwise --> C{chemistry_step_on_pipe: species? chem_solver}
    C -- no species --> N0[untouched]
    C -- auto --> AJ{bs_jit_applicable?}
    AJ -- yes --> J[_chem_step_bs_jit]
    AJ -- no --> L
    C -- bs_jit --> J
    J -- no params --> XJ[RuntimeError]
    J --> JK[bs_chemistry_step_pipe_arrhenius]
    JK --> CR{check_and_raise}
    CR -- a cell not converged --> RB[restore Y and E, ChemistryConvergenceError]
    C -- bs_brinkmeier --> BB{_chem_step_bs_brinkmeier: dispatcher?}
    BB -- yes --> BK[twc_brinkmeier BS step: serial or prange kernel]
    BK --> CR2{check_and_raise}
    CR2 -- a cell not converged --> RB2[restore Y, extras, E and coverages, raise]
    BB -- no --> L
    C -- numbalsoda --> NB{_chem_step_numbalsoda: compiled mechanism?}
    NB -- yes --> NK[numbalsoda_chemistry_step_pipe]
    NK -- backend missing, failed cell, coverage over-draw --> XN[RuntimeError, nothing written]
    NB -- no --> L
    C -- unknown name --> XU[ValueError]
    C -- lsoda --> L{any reaction?}
    L -- none --> N0
    L -- py_reaction_surface --> SF[chemistry_step_surface_lsoda]
    SF --> SV[_surface_chem_solve]
    SV -- failed or short --> XC[ChemistryConvergenceError, nothing written]
    SV --> SA[_assert_surface_coverage]
    SA -- gross over-draw --> XO[RuntimeError, nothing written]
    SA --> SC[_commit_surface_chem, then _snapshot_surface_q_rxn]
    L -- gas law --> G[_gas_lsoda_step_on_pipe: per-cell LSODA, staged]
    G -- a cell failed --> XC
    M -- species order wrong, LSODA failed --> XM[RuntimeError, nothing written]

What comes out. Per pipe: Y, the integrated Y_extras columns, the coverages in surface_state, U[:, 2] (the total energy absorbs the heat release; density and momentum are untouched) and, on the surface paths, q_rxn [W/m^3] at the end state. A BS-Brinkmeier step also leaves its telemetry summary on pipe._chem_telemetry. Pipes earlier in the loop keep their step when a later pipe refuses; the run then aborts.

Functions

  • src/aeolus1d/io/case.py::dispatch_advance
  • src/aeolus1d/io/case.py::_refuse_surface_chemistry_under_explicit
  • src/aeolus1d/io/case.py::_refuse_surface_mechanism_under
  • src/aeolus1d/solver/network_timeloop.py::_strang_symmetric_step
  • src/aeolus1d/solver/network_timeloop.py::_strang_supercycle_step
  • src/aeolus1d/solver/network_timeloop.py::_strang_chem_burst
  • src/aeolus1d/solver/implicit_network.py::_ImplicitRun.chemistry_supercycle
  • src/aeolus1d/solver/chemistry_ops.py::chemistry_step_on_pipes
  • src/aeolus1d/solver/chemistry_ops.py::surface_mechanism_step_on_pipe
  • src/aeolus1d/solver/chemistry_ops.py::chemistry_step_on_pipe
  • src/aeolus1d/solver/chemistry_ops.py::bs_jit_applicable
  • src/aeolus1d/solver/chemistry_ops.py::_chem_step_bs_jit
  • src/aeolus1d/solver/chemistry_ops.py::_chem_step_bs_brinkmeier
  • src/aeolus1d/solver/chemistry_ops.py::_chem_step_numbalsoda
  • src/aeolus1d/solver/chemistry_ops.py::chemistry_step_surface_lsoda
  • src/aeolus1d/solver/chemistry_ops.py::_surface_chem_solve
  • src/aeolus1d/solver/chemistry_ops.py::_assert_surface_coverage
  • src/aeolus1d/solver/chemistry_ops.py::_commit_surface_chem
  • src/aeolus1d/solver/chemistry_ops.py::_snapshot_surface_q_rxn
  • src/aeolus1d/solver/chemistry_ops.py::_gas_lsoda_step_on_pipe
  • src/aeolus1d/solver/chem_bs_jit.py::bs_chemistry_step_pipe_arrhenius
  • src/aeolus1d/solver/chem_bs_brinkmeier.py::bs_chemistry_step_pipe_brinkmeier
  • src/aeolus1d/solver/chem_bs_brinkmeier.py::bs_chemistry_step_pipe_brinkmeier_parallel
  • src/aeolus1d/solver/chem_numbalsoda.py::numbalsoda_chemistry_step_pipe
  • src/aeolus1d/solver/chem_diagnostics.py::check_and_raise