Skip to content

Low-Mach / quasi-steady aftertreatment solver — phased plan

Status: HISTORICAL / COMPLETE (kept as a design record). The quasi-steady solver this plan scoped has shipped and is the validated aftertreatment path: solver = "quasi_steady" runs the AVL C02006 keystone (500 s drive cycle in ~6.5 min) with a CI validation gate. For current usage read run-settings.md (the qs_* knobs) and the QS operating limits in scope-and-limitations.md; the v1 whole-network-QS scope restriction below still holds (mixed engine + AT multi-rate networks are tracked in BACKLOG.md).

User request 2026-05-20 (GT-POWER parity roadmap, v0.3 keystone): add an opt-in pressure-projection (SIMPLE/PISO) low-Mach solver mode for aftertreatment pipe networks, so a drive-cycle EAT run is no longer bound by the acoustic CFL. This document scopes the work into reviewable phases and pins the validation gate each phase must pass before the next starts.

The motivation is measured, not theoretical. The v0.3 profiling spike (scripts/profile_at_spike.py, recorded in docs/perf-history.md) showed that on 41_dpf_regen the explicit solver takes ~6 µs acoustic-CFL steps; a 500 s aftertreatment drive cycle is ~83.7 M steps ≈ 4.6 h wall, of which a low-Mach Δt (set by the ~1 ms thermal timescale) would remove ~167× → ~1.7 min. The spike also showed per-step cost is ~74 % Python orchestration and only ~11 % JIT kernels, so the lever is fewer steps (remove the acoustic floor), not faster kernels or a compiled core.

This is a multi-week solver-architecture project, delivered as one MR with a commit per phase (per the project's batching + commit-per-phase conventions).

Scope (v1)

  • Whole-network quasi-steady only. v1 supports networks where every pipe is quasi_steady. This matches the target case — the AVL Cruise-M C02006 EAT chain is a standalone DOC → DPF → SCR brick train with a prescribed inlet ṁ and outlet p, no engine breathing. A network mixing compressible engine pipes with quasi-steady AT pipes (true multi-rate sub-cycling) is an explicit non-goal / follow-up — engine pipes would force the small acoustic Δt back and erase the speedup.
  • Opt-in, compressible core untouched. The projection path is a new advance branch gated behind a per-pipe quasi_steady flag. Non-AT pipes never enter it, so the existing density-based explicit / implicit solvers are bit-identical. This is an acceptance gate, verified by a byte-for-byte regression run of the engine-breathing examples.
  • Composes with, does not replace, hllc_low_mach. The existing hllc_lm.py is a flux-dissipation tweak that is still acoustic-CFL bound; it is orthogonal to this mode (which changes the time-advance strategy, not the flux).

The method

For an AT pipe at Mach ≪ 1 with prescribed inlet ṁ and outlet p, acoustics are physically irrelevant. We march a segregated pressure-projection step on the existing collocated finite-volume grid, then advance the slow energy / species / soot / solid-temperature equations at the thermal Δt. Per macro step Δt_slow:

  1. Predictor (momentum). With the current pressure field and densities, solve the discrete momentum balance for a provisional velocity / face mass flux u*. In 1D quasi-steady low-Mach this is near-algebraic: pressure gradient balanced by convective flux + friction / porous / cake drag (the same Darcy-Forchheimer + cake_resistance·μ·u·soot terms already in kernels.py's quasi-1D residual). u* does not yet satisfy continuity.

  2. Pressure correction (Poisson). Form the continuity defect from u* and solve for a pressure correction p' whose discrete Laplacian ∇·((1/A_p)∇p') equals the mass imbalance (A_p = momentum diagonal coefficient). In 1D this is a scalar tridiagonal system → reuse solver/implicit_kernels.py:_block_thomas_solve at 1×1 block size. Poisson BCs: fixed-flux (Neumann on p') at the ṁ inlet, fixed-pressure (Dirichlet p'=0) at the p outlet. That inlet/outlet pairing is exactly what makes the system non-singular.

  3. Correct. u = u* − (1/A_p)∇p'; p ← p + α_p·p' (SIMPLE under-relaxation α_p; PISO adds a second corrector for less under-relaxation). Update face mass fluxes consistently.

  4. Density from EOS. ρ from the thermodynamic pressure + temperature (closures/eos.py:eos_pressure, ideal ρ=p/(R·T)). The thermodynamic pressure (≈ imposed outlet p plus the small frictional/cake drop) is decoupled from the dynamic pressure (∼ρu²) — standard low-Mach splitting.

  5. Advance slow equations. Energy (with regen / catalyst / chemistry heat release), species transport + reaction, soot loading, T_solid, T_cake — using the corrected mass flux for advection, at the thermal/chemical-limited Δt. Reuse the existing _into-style closure overlays and the post-RK operator-split tail verbatim.

Rhie-Chow (the load-bearing numerical decision)

aeolus1d is collocated FV (U_pad is cell-centered). A naive collocated pressure-correction decouples odd/even cells (checkerboarding). Rhie-Chow momentum interpolation is mandatory: face mass fluxes are built from a momentum-weighted face velocity that retains the face pressure gradient (not the averaged cell gradient), giving the Poisson a compact 3-point stencil. The alternative — a staggered velocity grid — was rejected: it would fork the entire U_pad / HLLC / junction-flux machinery. Rhie-Chow on the existing collocated layout is the right call, and it is built in from phase (a) because it is structural.

What this reuses (do not fork advance_network)

The shared post-RK tail in network_timeloop.py (regen mass balance, urea spray, fuel inventory, coolant, the dt_hrr_fraction controller) is factored into helpers callable from both the existing driver and the new one. The new advance_network_quasi_steady reuses: fill_network_ghosts, junction accumulate_rates, every _into closure overlay, the species/soot/T_cake/T_solid update blocks, and the post-RK tail. It replaces only the dt selection (drop the acoustic pipe_min_dt; use the thermal / HRR / fixed Δt) and the momentum-energy transport (predictor → Poisson → correct instead of explicit HLLC RK). The BCs it needs already exist: bc/mass_flow_inlet.py and bc/pressure_outlet.py. The Picard junction-coupling pattern is mirrored from solver/implicit_network.py.

With acoustics gone, dt_hrr_fraction (shipped in MR !41) becomes the primary Δt limiter — it now needs a sane floor/ceiling since there is no acoustic min() backstop.

Phased plan (one MR, commit per phase)

Each phase ends green (tests + ruff + mypy-strict) before the next starts.

  • (a) Single-pipe quasi-steady prototype. New solver/quasi_steady.py: predictor + scalar pressure-Poisson + corrector with Rhie-Chow, ideal-gas, uniform area, ṁ-in/p-out, no closures. pipe.quasi_steady flag (default False). Gate: steady ρu=const + analytic Darcy Δp; smooth p-field (the checkerboard test).
  • (b) SIMPLE/PISO outer loop. Under-relaxation + continuity-residual convergence (+ optional PISO 2nd corrector); map the tridiagonal solve onto _block_thomas_solve (1×1); EOS density; quasi-1D area variation. Gate: monotone residual history; nozzle steady; ideal-vs-real-gas ρ.
  • (c) Network / junction coupling. advance_network_quasi_steady with a junction Picard sweep; v1 all-QS, plenum-ghost junctions only (raise on neck/LODI/py, as advance_network_implicit does). Gate: 2-pipe + plenum global mass balance to machine tolerance.
  • (d) Closures + HRR Δt. Slot the slow-equation advance (species, cake, regen, catalyst, soot, T_cake, T_solid) + the post-RK tail + the HRR controller as the primary limiter. Gate: 41_dpf_regen in QS mode matches the explicit run's physics (regen timing, soot drain, O₂/CO₂ balance) within tolerance, at ≫ 1 ms steps.
  • (e) Schema + dispatch wiring. Full surface chain (io/schema.py → io/parsers.py → io/build.py → dispatch_advance → driver), a validate_case_pre_run gate rejecting unsupported combos (neck junctions, mixed compressible/QS network, AT pipe without a pressure anchor), GUI writer round-trip + test. Follows the aeolus1d-schema-extend surface chain to avoid the silent-parser-drop bug.
  • (f) Validation. C02006 drive cycle (step-count ~167× + physics + stable through the t≈116 s regen ignition) and a 41_dpf_regen QS variant. Bit-identical compressible-core regression gate. changelog + BACKLOG close. Then open the MR.

Hard parts & risks

  1. Rhie-Chow correctness — phase (a), highest risk; the checkerboard test is the gate.
  2. Bit-identical compressible core — gate everything behind quasi_steady; verify with a byte-for-byte engine-example regression.
  3. Poisson well-posedness — every QS pipe network needs a pressure anchor (the p-outlet Dirichlet); validate_case_pre_run rejects networks without one.
  4. Operator-split ordering — projection → energy/species → regen mass balance must preserve the existing sign/mass conventions.
  5. HRR Δt as sole limiter — needs a floor/ceiling now that the acoustic min() is gone.
  6. Multi-rate is a non-goal in v1 — document the boundary loudly so a mixed engine+AT case fails fast with a clear message rather than running slow or wrong.

Key files

New: src/aeolus1d/solver/quasi_steady.py. Edited: src/aeolus1d/solver/network_timeloop.py (new driver + factored tail helpers), src/aeolus1d/components/pipe.py (quasi_steady flag), src/aeolus1d/io/{schema,parsers,build,case}.py (surface chain), aeolus1d_gui/io/case_io.py (writer round-trip). Reused: src/aeolus1d/solver/implicit_kernels.py:_block_thomas_solve, src/aeolus1d/closures/{gas,eos}.py, src/aeolus1d/bc/{mass_flow_inlet, pressure_outlet}.py. Prior art mirrored: src/aeolus1d/solver/implicit_network.py, src/aeolus1d/components/coolant_pipe.py.