Performance history — the Phase Perf series¶

A record of the audit-driven optimization pass that ran 2026-05-07. Captures the seven phases, what changed at each one, and the empirical wall-clock numbers that drove the next decision.

The point of this document is not just what shipped — the commit log already says that. The value is the empirical-vs-speculative pattern: every phase started by profiling (cProfile against scripts/profile_perf1.py) and ended by re-profiling. Several candidate optimizations were rejected after benchmarking even though they looked correct on paper. Those are recorded here too, so the next person who reads "let's parallelize the inner loop" on a similar hot path has the previous evidence to lean on.

For the design rationale behind each Newton-solve path (direct + chord, JFNK + ILU(0), block-Jacobi, pure JFNK) and the preconditioner-shape lesson distilled from the empirical evidence here, see numerics.md.

Canonical workloads¶

Two TOML cases were the calibration targets. Both run under scripts/profile_perf1.py:

Slug	Path	What it stresses
40	`examples/40_implicit_adaptive_pipe.toml`	Single pipe, implicit BDF2, adaptive Δt, porous-Forchheimer plug. Ideal gas. The per-Newton-iter Jacobian/residual hot path with closures.
54	`examples/54_adaptive_picard_turbo.toml`	4-pipe Picard turbo network with adaptive Picard. Ideal gas. The orchestration overhead, BC dispatch, and per-pipe scaling.
58	`examples/58_tabulated_eos_implicit.toml`	Single pipe, implicit BDF2, tabulated CO₂ EOS. Validates the EOS-batched helpers when CoolProp is installed.

Cumulative trajectory¶

Phase	Case 40 wall	Case 54 wall (serial)
Audit baseline	0.323 s	39.66 s
Perf-1 (EOS hot-path batching)	0.323 → 0.323 s*	39.66 → 39.66 s*
Perf-2-quick (D_J broadcast)	0.323 → 0.241 s (−25 %)	39.66 → 27.70 s (−30 %)
Perf-2 (closure-Jac JIT + line-search)	0.241 → 0.110 s (−54 %)	27.70 → 18.09 s (−35 %)
Perf-3 (chord Jacobian reuse)	0.110 → 0.048 s (−56 %)	18.09 → 8.24 s (−54 %)
27c (opt-in parallel pipes)	n/a (single pipe)	8.24 → 5.73 s (−30 % opt-in)
Perf-4 (JIT'd BC dispatchers)	0.048 → 0.020 s (−58 %)	8.24 → 2.67 s (−68 %, serial)
27c-cells (`prange` Jacobian — reverted)	0.020 → 0.041 s (×2 slower)	2.67 → 3.4 s (×1.3 slower)
Final	0.020 s	2.67 s
Cumulative reduction	94 % (16 ×)	93 % (15 ×)

* Case 40 and 54 are ideal-gas; Phase Perf-1 targets the tabulated and cubic-EOS Newton-Jac path which neither workload exercises. Use scripts/bench_perf1_tabulated.py or example 58 for direct Perf-1 measurements.

Phase Perf-1 — EOS / solver hot-path optimization sweep¶

Commit: 338126d. Trigger: initial audit identified redundant bilinear table lookups + FD-fallback derivative computations as the EOS hot-path drag.

What changed

Wired the B23-precomputed ∂p/∂T and ∂p/∂ρ panels directly into eos_dpdT_dpdrho_jit (replaces a 4-call central-difference fallback that the Phase B23 commit had flagged as a future optimization).
Closed-form ∂e/∂T, ∂e/∂ρ for tabulated EOS via _bilinear_grad on the energy panel — same FD-fallback pattern eliminated for the energy derivatives.
Batched _bilinear_two and _bilinear_grad_two helpers share the (T, ρ) cell-index work across two panels in one cell pass.
New fused dispatcher eos_a_and_derivs_jit(T, ρ, kind, params) returns (a, ∂p/∂T, ∂p/∂ρ, ∂a/∂T, ∂a/∂ρ) in a single call — collapses three separate dispatches per face per Newton iter in hllc_jac._eos_primitive_gradients. The cubic-EOS body shares p_v / p_T / p_vv / p_TT / p_vT across all five outputs; the tabulated path takes a single shared bilinear- grad cell pass on the dpdT and dpdρ panels.
Fused tabulated (p, a) recovery in primitives_from_U_eos_jit (was three separate bilinear traversals).
Block-Thomas _inv3x3_into writes through a hoisted scratch buffer, eliminating per-cell np.empty((3, 3)) allocations.

Validation

Bit-exact (rel diff < 1e-12) on Peng-Robinson CO₂ at T = 350 K, ρ = 80 kg/m³ vs the legacy separate-call paths. 82 EOS / HLLC-Jac / real-gas / implicit-real-gas regression tests pass.

Direct wall-clock impact

Cases 40 and 54 don't exercise the analytic-EOS Jacobian path (both are ideal-gas, which has its own fast-path). To measure Perf-1 directly, run scripts/bench_perf1_tabulated.py or profile examples/58_tabulated_eos_implicit.toml.

Phase Perf-2-quick — D_J broadcast¶

Commit: d56a644. Trigger: Post-Perf-1 cProfile flagged _newton_step's self-time at 29 ms (single-pipe) and 15.5 s (network) — most of which was a per-cell Python loop building D_J = I − dt · D_blk.

What changed

# Before:
D_J = np.empty_like(D_blk)
for i in range(pipe.N):
    D_J[i] = _I3 - dt_eff * D_blk[i]

# After:
D_J = -dt_eff * D_blk
D_J[:, _DIAG3, _DIAG3] += 1.0

Wall-clock

case 40: 0.323 → 0.241 s (~25 %)
case 54: 39.66 → 27.70 s (~30 %)
_newton_step self-time: 29 ms → 1 ms (single-pipe); 15.5 s → 3.1 s (network)

The single-broadcast was dramatically more expensive than expected — the per-cell Python interp overhead × N × Newton iters × steps adds up.

Phase Perf-2 — JIT-ify closure-Jacobian dispatcher and Newton line-search¶

Commits: 9ea8bd4 (closure-Jac) + a14d29a (line-search). Trigger: Post-Perf-2-quick cProfile showed _augment_jacobian_with_closures at 25–37 % of single-pipe wall time — a Python for i in range(pipe.N) loop dispatching closure-Jacobian helpers one cell at a time.

What changed

Five real-gas closure helpers were silently plain-Python despite calling @njit kernels internally: _real_gas_T_app_chain, _real_gas_heat_jacobian, _real_gas_heat_jacobian_from_chain, _real_gas_porous_sutherland_jacobian, _real_gas_porous_sutherland_jacobian_from_chain. Decorated @njit(cache=True, fastmath=True, inline="always").
New JIT kernel _augment_closure_jac_analytic_jit runs the per-cell loop body inside @njit, with all dispatch flags passed as integers and missing-closure params passed as length-1 sentinel arrays (numba can't take Python None).
Line-search positivity-guard pieces became _line_search_trial_jit (fused trial-state alloc + ρ>0 / e>0 checks) and _residual_norm_jit (fused F = U_trial − U_hist − dt·R + ‖F‖). Eliminated 5 transient numpy allocs per backtrack.

Wall-clock

case 40: 0.241 → 0.110 s (~54 %)
case 54: 27.70 → 18.09 s (~35 %)

Test note — relaxed test_shared_real_gas_chain_matches_separate_helpers from assert_array_equal to assert_allclose(rtol=1e-12): both helpers now run under fastmath; FMA fusion gives 1 ULP diffs even though the algebra is identical.

Phase Perf-3 — chord (lazy-Jacobian) reuse¶

Commit: 11d32f1. Trigger: Post-Perf-2 the dominant cost on case 54 was the JIT Jacobian assembly itself (5.6 s of 18.1 s = 31 %). The kernel was already optimal; the only remaining lever was to call it less often.

What changed

Hold J = I − dt · ∂R/∂U fixed across multiple inner Newton iterations, accepting linear-with-good-constant convergence in exchange for 30–50 % less Jacobian work per Newton solve.

Two NewtonConfig knobs:

chord_max_age (default 2): hard cap on iterations between refreshes. 0 reproduces the legacy full-Newton path bit-exactly.
chord_rate_threshold (default 0.5): if F_norm_new / F_norm_old > threshold, the chord has stopped matching the local Newton direction — refresh the Jacobian on the next iter. Line-search rejection on a chord-aged iteration also forces a refresh.

The cache is the post-broadcast (L_J, D_J, Up_J) triple, so chord-reuse skips the closure-Jac augmentation, the geometry-Jac augmentation, and the numpy broadcasts too.

Wall-clock

case 40: 0.110 → 0.048 s (~56 %)
case 54: 18.09 → 8.24 s (~54 %)

Side-effect — chord-Jacobian shifted Newton iter counts on each Picard iter, which the adaptive Picard controller's "easy" detection picks up. One regression test (test_adaptive_picard_settles_on_plenum_only_network) was relaxed from == 2 to in (2, 3) to track the shifted but still-valid controller trajectory.

Phase 27c — opt-in parallel per-pipe Newton solves¶

Commit: 4316e08. Trigger: Post-Perf-3 the Python orchestration around _newton_step was still 9 % of network wall time, and the per-pipe Newton solves are theoretically independent during a Picard iteration (junction state is frozen between _fill_component_ghosts calls).

What changed

New NewtonConfig.parallel_pipes flag (default False).
New [case] parallel_pipes TOML field on CaseHeader.
Helper _try_parallel_first_pass runs the per-pipe Newton solves via ThreadPoolExecutor (max 8 workers) and returns (max_n_iters, total_n_iters) on success or None on any pipe's failure.
On success the parallel result accumulates the Newton-iter stats and the sequential per-pipe loop is bypassed. On failure all pipes restore to the start-of-iter snapshot and the legacy sequential dt-halving loop runs unchanged.

Wall-clock

case 54 (parallel, opt-in): 8.24 → 5.73 s (~30 %)

Validation — test_parallel_pipes_phase27c.py: 3 tests (serial-vs-parallel rtol=1e-10, default-off contract, single- pipe bypass).

Why default-off after Perf-4 — see next phase.

Phase Perf-4 — JIT'd BC dispatchers + chord knobs on TOML¶

Commit: 0acfd6a. Trigger: Post-27c parallel-pipes cProfile showed _thread.lock.acquire at 67 % of network wall time. The Python-level BC dispatchers (stagnation_inlet.bc, pressure_outlet.bc) were holding the GIL during the per-pipe apply_ghosts callback path, which serialized the parallel threads.

What changed

make_stagnation_inlet_bc and make_pressure_outlet_bc now delegate to @njit kernels (_stagnation_inlet_apply_jit with a JIT'd _solve_a_b_jit for the characteristic+stagnation quadratic; _pressure_outlet_apply_jit with a Python-side fallback to transmissive_bc on supersonic outflow).
[case] chord_max_age / chord_rate_threshold surface the Phase Perf-3 chord knobs to TOML cases.

Wall-clock

case 40: 0.048 → 0.020 s (~58 %)
case 54 (serial): 8.24 → 2.67 s (~68 %)
case 54 (parallel): 5.73 → 5.10 s (~11 %)

The unexpected finding

JIT-ifying the BCs sped up the serial path more than the parallel one. Once the BC apply is GIL-released and inlined into the @njit region, the per-pipe sequential pipeline is bound by JIT kernel work rather than Python orchestration. Thread-coordination overhead (ThreadPoolExecutor setup per Picard iter, BC closure dict lookups, _thread.lock.acquire) now dominates the parallel path.

Result: serial beats parallel ~2 × on the canonical 4-pipe network case (2.67 s vs 5.10 s). parallel_pipes stays default off — defaulting it on would have been a regression.

Phase 27c-cells — `prange` Jacobian (prototyped, reverted)¶

Commit: b2667ae (negative-result documentation; the prototype lives in git history pre-revert). Trigger: Post-Perf-4, _assemble_jacobian_fo_ideal_uniform_jit is the dominant remaining cost on case 54 (0.85 s of 2.67 s = 32 %). Theoretically parallelizable across the face loop — each face writes to two cells (f-1 and f), so adjacent faces conflict but even/odd-faces colouring gives two disjoint passes.

What was tried

@njit(cache=True, fastmath=True, parallel=True)
def _assemble_jacobian_fo_ideal_uniform_jit(...):
    # Pass 1: faces 0, 2, 4, …
    for k in prange(n_even):
        f = 2 * k
        _accumulate_face_contribution_jit(f, ..., L, D, Up)

    # Pass 2: faces 1, 3, 5, …
    for k in prange(n_odd):
        f = 2 * k + 1
        _accumulate_face_contribution_jit(f, ..., L, D, Up)

The colouring is correct (even faces write to disjoint cell sets, same for odd). Each prange iteration's scratch (dF_dUL, dF_dUR, etc.) was made local to the helper so each parallel iteration gets its own.

Wall-clock — slower than serial

case 40: 0.020 → 0.041 s (~2 × slower)
case 54: 2.67 → 3.4–3.8 s (~1.3–1.4 × slower)

Numba's parallel=True setup overhead per call dominated the parallel work at N = 80 cells per pipe. Break-even is somewhere around N = 200–500 cells per pipe based on rough scaling. Reverted; the colouring scheme is recorded in BACKLOG.md for replay if a long-pipe bench case ever justifies it.

The lesson — kernel-internal prange reads as the obvious optimization (the data dependencies clearly factor), but at small N the framework overhead dominates. Always benchmark.

Phase Perf-5 — odds and ends¶

Commit: b2667ae. What landed

Fixed examples/58_tabulated_eos_implicit.toml (EOS block was under [gas.eos] not [pipe.eos]; missing pipe = field on [[bc]] blocks).
Profile harness learns case 58.
Picard outer-loop audit: profile shows _ghost_snapshot_max_rel_change at 1.7 % of case 54 wall; below the threshold where JIT-ifying pays back. Skipped.
Phase 27c-cells negative result documented.

Phase Perf-6 — chemistry advection JIT¶

Commit: 53f9bab. Trigger: Phase Perf-history-2 profile of examples/41_dpf_regen.toml (DPF + chemistry + thermal cake) showed _species_rate was 56 % of total wall time (1.19 s of 2.13 s). Same Python-loop-around-JIT-kernels pattern as Phase Perf-2's _augment_jacobian_with_closures — the per-cell for i in range(N) body in network_timeloop._species_rate called hllc_flux(...) per cell from Python, with Python attribute access (pipe.U_pad, pipe.Y, pipe.gamma) per iteration.

What landed

Two @njit kernels in aeolus1d/solver/species_advection_njit.py:

species_advection_primary_jit(U_pad, Y_pad, gamma, dx, N, Ng) — primary-species-only path. Returns d(rho*Y)/dt.
species_advection_extras_jit(U_pad, Y_pad, Yx_pad, gamma, dx, N, Ng, n_extra) — primary + Phase-22 extras. Returns (d(rho*Y)/dt, d(rho*Yx)/dt).

hllc_flux is @njit and inlines into the kernel; the two variants exist because numba's handling of zero-column 2-D arrays is awkward and the n_extra=0 / n_extra>0 branches have different return shapes. The Python wrapper in _species_rate keeps:

The species-inflow BC overlay (per-cell Y_pad build with optional inflow schedule lookup) — needs Python-side dict access on pipe.species_inflow_left/right.
The reaction callable invocation (pipe.py_reaction / pipe.py_reaction_multi) — Python callbacks by definition.
The DPF filter-trap overlay via filter_trap_rates (which is itself @njit).

This isolation is intentional: keep the JIT kernel pure-numeric + stable input shapes; let Python handle the per-pipe schema discovery and overlay calls.

Wall-clock

Workload	Baseline	Phase Perf-6	Speedup
`41_dpf_regen` (default `t_end = 0.06`)	4.63 s	1.74 s	2.66 ×

Apples-to-apples wall-clock measurement (best of 3 runs without cProfile, on the same machine, both with pre-warmed numba caches).

cProfile breakdown of _species_rate:

	Baseline	Phase Perf-6
`_species_rate` self-time	1.19 s	0.27 s
`species_advection_extras_jit` (new)	—	0.11 s
End-to-end species-rate cost	~1.19 s	~0.38 s

The Python wrapper still has measurable overhead (Y_pad build + reaction overlay + filter call) — that's ~0.27 s. The JIT kernel itself is ~0.11 s. End-to-end the chemistry hot path is ~3 × faster.

Tests

All 170 chemistry / DPF / SCR / filter integration tests pass.
The Phase Perf-history-2 single-shot regression bench was already in place (tests/integration/test_dpf_*) — bit-exact against the pre-Perf-6 baseline because the algebra is unchanged (we just moved the for-loop from Python to numba-compiled code).

Why no make_strang_chemistry change

The Strang-split path uses _species_rate(pipe, R, include_reaction=False) for the transport half; same kernel fires, same speedup applies. Reaction half is BS_jit / LSODA and was already JIT (BS) or external-C (LSODA).

Phase Perf-6r2 — chemistry / filter Python wrapper hoist¶

Trigger: Re-profiled case 41 post-Perf-6. Wall was 2.55 s down from 4.63 s, but the breakdown told a story — the wrappers that call the JIT kernels still spent more wall in Python attribute-access dispatch than in the kernels themselves:

Wrapper	self-time	% of wall
`_species_rate`	0.386 s	15.1 %
`regen_into`	0.237 s	9.3 %
`cake_into`	0.079 s	3.1 %

Each wrapper was doing pipe.filter.X 10-25 times per call (× ~20 K RHS calls/s = ~250 K attribute walks/s); pipe.n_extra_species is a @property called 3× per filter wrapper; the filter ones did 5+ float()/int() casts on filter parameters inside the call site (so Python evaluated each before dispatching to numba).

False start: First tried to eliminate the Y_pad / Yx_pad allocations in _species_rate by changing the JIT-kernel signature to read pipe.Y / pipe.Y_extras directly with explicit ghost scalars. This regressed 22 % (2.55 s → 3.13 s). Root cause: the new kernel had to insert a per-cell-per-species if right_is_ghost branch into the hot inner loop because the interior and ghost rows now live in different arrays. The old Yx_pad[idx, k] was a single indirect load with no branch; numba SIMD-vectorised that loop. The pad-array allocation cost was tiny compared to the branchless-loop benefit. Reverted.

Actual fix: Hoist every repeated attribute lookup inside the three wrappers into a local variable at the top of the function — f = pipe.filter, then f.kind / f.soot_loading / etc; cache pipe.U, pipe.Y, pipe.Y_extras, pipe.gamma, pipe.n_extra_species; hoist the from aeolus1d.closures.filter import filter_trap_rates out of _species_rate to module scope (was doing an import-system lookup per RHS call even though Python caches modules — the dict lookup itself isn't free).

Result (CORRECTED 2026-05-15)

The original commit message claimed 1.43× wall on case 41 (2.55 → 1.78 s, 7-run median). That was a system-load artefact — measurements drifted between the two benchmark runs because the WSL2 host was under different load.

A proper back-to-back cProfile A/B comparison (PRE → HEAD → PRE → HEAD) shows the real per-function and wall-clock deltas:

Function	PRE-r2	HEAD	Δ
`_species_rate` self	0.275 s	0.226 s	−18 %
`regen_into` self	0.171 s	0.156 s	−9 %
`cake_into` self	0.052 s	0.070 s	+35 %*
`pipe.U` property calls	269 724	249 626	−7 %
Total cProfile wall	2.387 s	2.293 s	−4 %

(*cake_into's cProfile self-time went up because cProfile's per-function overhead is now a larger fraction of a smaller function — actual production wall is faster, just hard to measure with the wrapper's small self-time.)

So the real win is ~4 % wall and 9–18 % on the targeted wrappers, not 43 %. Still positive, still worth keeping, but a much more modest gain than the original claim.

Lessons recorded

When a Python wrapper sits on top of a JIT kernel and the wrapper has many obj.attr walks per call, hoist the walks before reshaping the kernel. The wrapper overhead is often the actual hotspot and the cleanest fix has nothing to do with the kernel.
Wall-clock benchmarks under WSL2 (and any shared hardware) drift with host load. The audit-driven Perf series used cProfile back-to-back deltas on the premise that per-function self-time stays comparable across system load while wall time does not. That premise is false (retracted 2026-10-07; BACKLOG_HISTORY.md): one profile read 82.29 s at -n 4 against 44.64 s alone, because cProfile charges the contended time to whatever function was running. Compare self-time only between runs under the same, recorded contention (the perf records now carry CPU pressure and steal). A bench reporting "1.4×" that's really 1.04× is worse than no bench at all — it burns trust. Always confirm wall-time wins with a tight back-to-back A/B (HEAD → revert → HEAD), not before-and- after measurements separated by minutes of unrelated work.

Phase Perf-6r3 — junction-code Python wrapper hoist¶

Trigger: Re-profiled engine case 50 (2-cyl engine network) to find what's hot. The chemistry wrappers Perf-6r2 targeted are never called on engine cases (no DPF, no chemistry), so that phase had zero effect on engine wall. The actual engine hot path is the junction code:

Function (case 50)	Self-time	Per-call cost
`_packed_residual_fo_ser` (JIT)	0.94 s	36 µs (true compute)
`accumulate_rates`	0.39 s	7.5 µs (Python wrapper)
`_pipe_port_flux`	0.38 s	2.0 µs (Python dispatch)
`_fill_plenum_ghost`	0.30 s	1.5 µs
`_pipe_fuel_key`	0.28 s	0.9 µs
`apply_port_bcs`	0.28 s	5.3 µs

The cumulative junction-Python self-time at HEAD-of-Perf-6r2 was 1.94 s — 27 % of the engine-50 cProfile wall (7.31 s).

Investigation showed two systemic patterns worth fixing in this pass:

accumulate_rates had a per-call import:
```
from aeolus1d.components.species_coupling import (
    iter_pipe_species, read_pipe_boundary_fuel,
)
```
inside the function body. ~52 K calls/case × Python's import-system dict lookup adds up.
_pipe_port_flux had five getattr(pipe, "field", default) calls per invocation (scheme, use_per_cell_gamma, gamma_pad, has_real_gas, muscl_limiter_code). All five are real Pipe dataclass fields with the same defaults the getattr calls supplied — so the getattr machinery was pure overhead. The function also imported boundary_face_flux_muscl_var_gamma_jit inside the muscl-var-γ branch.

Fix

Hoisted both per-call imports to module scope.
Replaced getattr with direct attribute access on the five known-existing Pipe fields.
Cached pipe.U_pad once in _pipe_port_flux (used 4× in the common first-order branch).

Result (back-to-back cProfile A/B, case 50)

Function	PRE	HEAD avg	Δ
`accumulate_rates` self	0.390 s	0.353 s	−10 %
`_pipe_port_flux` self	0.379 s	0.350 s	−8 %
Total cProfile wall	7.31 s	7.23 s	−1 %

Modest because the JIT kernel _packed_residual_fo_ser still owns 13 % of wall (true compute) and the other junction methods (_fill_plenum_ghost, _pipe_fuel_key, apply_port_bcs) weren't touched in this pass — they're candidates for a follow-up. The function-level deltas confirm the hoist hit its targets cleanly.

Phase Perf-6r4 — junction follow-up (remaining wrappers)¶

Trigger: Perf-6r3 hoisted accumulate_rates + _pipe_port_flux. The follow-up profile of case 50 showed three more junction wrappers in the top-8 self-time list with the same systemic problem — per-call imports inside the function body:

_pipe_fuel_key (0.299 s self, 313 K calls) — imported pipe_fuel_key per call.
_pipe_boundary_fuel — imported read_pipe_boundary_fuel per call.
apply_port_bcs (0.286 s self) — imported write_pipe_inflow_fuel once per port, then iter_pipe_species + write_pipe_inflow_fuel again inside the per-port loop body.

Plus _fill_plenum_ghost was reading pipe.U_pad 3× inside its inner write loop.

Fix

Hoisted pipe_fuel_key and write_pipe_inflow_fuel to the module-level import block (iter_pipe_species / read_pipe_boundary_fuel were already there from Perf-6r3).
Cached self.fuel_species_name and self.m_species as locals in the apply_port_bcs per-port loop so the per-iteration self.X walks collapse to single attribute reads at the function top.
Cached Up = pipe.U_pad once in _fill_plenum_ghost.

Result (back-to-back cProfile A/B, case 50)

Function	PRE-r4	HEAD avg	Δ
`_pipe_fuel_key` self	0.299 s	(out of top 8)	eliminated
`apply_port_bcs` self	0.286 s	(out of top 8)	eliminated
Total cProfile wall	7.38 s	7.27 s	−1.5 %

The per-call import cost was the dominant factor for both — once hoisted, they fell off the hot list entirely. New top-of-hot-list entries on case 50 are now cylinder methods (accumulate_rates 0.246 s, apply_port_bcs 0.223 s) which carry the same per-call import pattern; candidates for Perf-6r5.

Phase Perf-6r5 — cylinder-side wrappers¶

Trigger: Post-Perf-6r4 profile of case 50 (2-cyl engine) showed the cylinder analogues of the junction wrappers as the new top-of-hot-list entries:

cylinder.accumulate_rates — 0.351 s self
cylinder.apply_port_bcs — 0.311 s self

Both carry the same per-call-import pattern Perf-6r3/r4 fixed in junction.py: _pipe_fuel_key imports aeolus1d.components.species_coupling.pipe_fuel_key per call, _pipe_boundary_fuel imports read_pipe_boundary_fuel per call, and apply_port_bcs imports write_pipe_inflow_fuel inside the per-port loop (~52 K calls/case-50 each).

Fix

Hoisted all three imports to module scope (extending the import block introduced in Perf-6r3 for junction.py).

Result (back-to-back cProfile A/B, case 50)

Function	PRE	HEAD	Δ
`cylinder.accumulate_rates` self	0.351 s	0.323 s	−8 %
`cylinder.apply_port_bcs` self	0.311 s	0.292 s	−6 %

(Run-to-run cProfile noise was higher than usual on this pass — host system load increased between measurements. Function-level self-time still tracks the import-hoist effect cleanly even when total wall jitters.)

The wins are smaller than Perf-6r4's because case 50 doesn't carry a fuel species column, so the write_pipe_inflow_fuel hot-path inside apply_port_bcs is gated off — only the two key lookups (_pipe_fuel_key / _pipe_boundary_fuel) actually fire per port-step. On a fuel-tracking engine case (case 38, 39, or 52) the savings would compound further.

Phase Perf-6r2-corrections-2 — profile-based regression harness¶

Trigger: The Perf-6r2 commit's wall-claim correction documented the lesson — wall-clock benches drift with system load and a "1.43×" claim that's really 1.04× burns trust. The fix is structural: install a regression gate that uses cProfile per-function self-time instead of wall time, so the gate measures the actual wrapper hotspot at the function level instead of the case at the wall level.

What landed

tests/benchmarks/_timer.py gained a bench_profile() helper that runs a callable once under cProfile and returns a {"basename:line(name)": self_time_seconds} dict, plus an assert_profile_self_time_under() matcher with a list-of- substrings disambiguator (so ["junction.py", "(accumulate_rates)"] pins the junction one, not the cylinder one).
tests/benchmarks/test_profile_regression.py (new) — parametrised gates against case 41 (chemistry wrappers from Perf-6r2) and case 50 (junction wrappers from Perf-6r3/r4). Marker @pytest.mark.profile_regress (registered in pyproject.toml) so the gates only run on pytest -m profile_regress and don't slow the default suite.
The existing bench_regress micro-bench gates stay as-is — they're nanosecond-scale @njit kernel benches where wall-clock median over many reps IS the right tool. cProfile would add 30-50 % overhead per kernel call and dwarf the measurement. The two gates partition the perf-regression surface cleanly:
bench_regress — kernel-level, nanosecond-scale, wall-clock median over 100s of reps. Tolerance ~3×.
profile_regress — full-case, function-level, cProfile self-time. Tolerance ~3× over the post-hoist baseline.

Lesson

The right perf metric depends on what's being measured. For a hot @njit kernel, wall-clock is the only thing that matters and median-of-many-reps filters noise. For a full-case Python wrapper hotspot, per-function self-time under cProfile is the right measure because it's normalised against wall by virtue of being a fraction of the cProfile-overhead-dominated profile call — same function's self-time stays within ~1 % across runs even when wall jitters by tens of percent under host CPU contention.

Phase Perf-6r6 — turbo-side wrapper hoist¶

Trigger: Closes the four-wrapper hoist arc started in Perf-6r3. The pattern was systemic — every pipe-component wrapper that called into species_coupling did so via a per-call import inside the function body. After Perf-6r3/r4 fixed the junction wrappers and Perf-6r5 fixed the cylinder wrappers, the only remaining one was _pass_fuel_through_turbo in aeolus1d.components.turbo, which Compressor and Turbine apply_port_bcs both invoke.

Fix

Hoisted the same three imports (pipe_fuel_key, read_pipe_boundary_fuel, write_pipe_inflow_fuel) to the module-level import block, matching the junction.py and cylinder.py patterns from r3/r4/r5.

Result

Hoist is invisible on case 50 (no turbo) and case 54 (no fuel tracking) — the function only fires when a turbo + fuel-tracking pipes coexist (cases 38, 39). For those cases the savings would match the cylinder wrapper deltas from Perf-6r5 (~6-8 % per function self-time). Mechanical change; no measurable wall regression on the existing canonical bench cases.

This commit closes the four-wrapper-arc:

r3 — junction _pipe_port_flux + accumulate_rates
r4 — junction _pipe_fuel_key + _pipe_boundary_fuel
apply_port_bcs + _fill_plenum_ghost
r5 — cylinder _pipe_fuel_key + _pipe_boundary_fuel
apply_port_bcs
r6 — turbo _pass_fuel_through_turbo

The remaining engine-case hot list is now dominated by the JIT kernels themselves (_packed_residual_fo_ser ~13 % of wall, _assemble_jacobian_* for implicit cases) which are true compute, not Python wrapper overhead.

Phase Perf-8 — getattr-dispatch removal (species_coupling + muscl boundary)¶

Trigger: Post-Perf-6r6 case-50 profile showed {built-in method builtins.getattr} at 0.264 s self-time across 1.49 M calls — explicit getattr(obj, "name", default) dispatch was eating ~1.5× the cost of the JIT residual kernel (true compute). Wrapper-hoist arc was closed; the next layer of Python overhead is the explicit getattr dispatch in 5 hot sites.

Hot sites identified (per-port-per-step call counts × cost):

species_coupling.pipe_fuel_key — 3 getattrs/call, 443 K calls = ~1.3 M getattrs.
species_coupling.iter_pipe_species — 3 getattrs/call.
junction._pipe_port_flux_muscl — 2 getattrs (muscl-only).

All target real Pipe dataclass fields (species_name / species_names) or properties (has_species / has_real_gas / muscl_limiter_code). The getattr(..., default) calls supply the same defaults the dataclass already has, so the dispatch was pure overhead.

Fix

Direct attribute access on every site: pipe.X instead of getattr(pipe, "X", default). Same pattern as Perf-6r3's _pipe_port_flux cleanup but extended to the species-coupling helpers and the muscl boundary path.

Result (back-to-back cProfile A/B, case 50)

Function / metric	PRE	HEAD avg	Δ
`pipe_fuel_key` self	0.236 s	0.20 s	−15 %
`{builtins.getattr}` self	0.264 s	(out of top 12)	eliminated

The 1.49 M getattr calls drop to ~120 K (the legitimate duck-typed sites: the EHR port_u_rates lookup and the apply_fuel_step Component-Protocol fishing). Total wall moves only modestly (the JIT kernel _packed_residual_fo_ser still dominates) but the function-level overhead reduction is unambiguous.

The other getattr sites left intact:

network_timeloop._has_port_u_n and the per-port getattr(p, "u_n", None) snapshot list-comprehension — Port.u_n is now a real Optional Pipe field, but CylinderPort doesn't have it. The getattr is the right pattern for the heterogeneous Component port lists.
Component.apply_fuel_step getattr in advance_network_implicit — same story, only some Component implementations have this method.

Phase Perf-9 — investigate remaining getattr calls (no-op)¶

Trigger: Following the Perf-8 hoist, the {built-in method builtins.getattr} self-time dropped from 0.264 s / 1.49 M calls to 0.031 s / 287 K calls — an 80 % reduction. Phase Perf-9 set out to investigate whether the remaining 287 K calls (the legitimate duck-typed sites in _has_port_u_n, the per-port port.u_n snapshot list-comp, and the Component.apply_fuel_step Protocol fishing) could be replaced with isinstance dispatch for additional savings.

Investigation

The remaining 287 K calls cost 0.031 s self-time on case 50 (0.4 % of wall). Per call: ~0.1 µs — about as fast as getattr can get. The candidate sites:

_has_port_u_n — 3 getattrs per junction per RK stage. Could short-circuit to isinstance(junction, EHRJunction), but EHRJunction lives in aeolus1d.components.junction_ehr which isn't imported by aeolus1d.solver.network_timeloop (and importing it module-level would create a circular-import surface to manage).
port.u_n snapshot — 1 getattr per port per stage. Same EHRPort isinstance check would apply, with the same circular-import cost.
apply_fuel_step Protocol fishing — 1 getattr per junction per step. A clean isinstance against the mixin/Protocol that defines apply_fuel_step would work, but this is exactly the kind of thing Python Protocols exist to avoid (the Component contract is structural, not nominal).

Decision: no change

Estimated savings: ~0.01–0.02 s wall on case 50 — within the profile-noise band. The code-complexity cost (module-level imports of EHRJunction / EHRPort, isinstance dispatches replacing duck-typed Protocol satisfaction) outweighs the perf benefit. Documented as a deliberate no-op so a future contributor doesn't re-investigate.

The Perf-8 → Perf-9 trajectory is the canonical "diminishing returns" pattern: the first wave of cleanup eliminates 80 % of the cost; the long tail of remaining cost isn't worth chasing. Profile-driven discipline says stop and move on; the next perf phase should target a different hot path (engine-case JIT kernels, or a different wrapper layer).

Phase Perf-history-3 — re-profile post-Perf-8 to find next-tier targets¶

Investigation only (no code changes). Re-ran scripts/profile_perf1.py on the three canonical workloads after the wrapper-hoist + getattr arc closed (Perf-6r2..r6 + Perf-8). Documents the new top-of-hot-list so future contributors don't re-investigate the same surfaces and so the next Perf phase has a clear target list.

Case 41 (chemistry / DPF) — 2.32 s cProfile wall:

Function	Self	% wall
`_species_rate` (already hoisted r2)	0.229 s	9.9 %
`regen_into` (already hoisted r2)	0.154 s	6.6 %
`species_advection_extras_jit` (JIT)	0.119 s	5.1 %
`pipe_residual_quasi1d_first_order`	0.117 s	5.0 %
`Pipe.U` property (249 K calls)	0.082 s	3.5 %
`cake_into` (already hoisted r2)	0.049 s	2.1 %
`Pipe.n_extra_species` (148 K calls)	0.037 s	1.6 %
`np.zeros` allocations (100 K calls)	0.028 s	1.2 %

Case 50 (engine network) — 7.05 s cProfile wall:

Function	Self	% wall
`_packed_residual_fo_ser` (JIT)	0.975 s	13.8 %
`junction.accumulate_rates`	0.362 s	5.1 %
`junction._pipe_port_flux`	0.361 s	5.1 %
`junction._fill_plenum_ghost`	0.314 s	4.5 %
`junction._ghost_rho_E`	0.274 s	3.9 %
`cylinder.accumulate_rates`	0.249 s	3.5 %
`cylinder.apply_port_bcs`	0.222 s	3.1 %
`junction.apply_port_bcs`	0.190 s	2.7 %
`hllc.hllc_flux` (JIT, 192 K calls)	0.167 s	2.4 %
`Pipe.has_species` (992 K calls)	0.162 s	2.3 %
`species_coupling.pipe_fuel_key`	0.144 s	2.0 %
`cylinder._kinematics`	0.142 s	2.0 %

Case 52 (knock event) — 1.09 s cProfile wall:

Smaller-scope case; _packed_residual_fo_ser (0.197 s, 18 %) + advance_network (0.185 s, 17 %) own the wall. Cylinder wrappers (accumulate_rates 0.064, apply_port_bcs 0.043) combined ~10 %. Limited optimisation surface.

Future-work candidates (in rough priority order):

Phase Perf-10 — Pipe.has_species property → cached field. 992 K calls on case 50 for a property that just evaluates self.species_name != "". Currently nothing mutates species_name post-construction, so a _has_species: bool computed in __post_init__ would eliminate the dispatch entirely. Realistic win: ~2 % wall on case 50.
Phase Perf-11 — Pipe.U property dispatch reduction. 249 K calls on case 41 + 392 K on case 50. Each does self.U_pad[Ng:Ng + N] (a small numpy slice). Can be cached on Pipe (refreshed when U_pad is reassigned — currently only on __post_init__ and the rare resize-then-reassign path). Realistic win: ~3 % wall on case 41, ~2 % on case 50.
Phase Perf-12 — np.zeros scratch-buffer pre-allocation in regen_into. 100 K allocations per case 41 for the Y_o2 / Y_no2 / Y_no / d_soot_dt / dT_cake_dt arrays. Pre- allocate on FilterState (lazy-initialised on first call, sized to pipe.N). Realistic win: ~1 % wall.
Phase Perf-13 — junction._fill_plenum_ghost per-cell loop fusion. 191 K calls × 3 array writes each (with Ng=1, the loop runs once but does 3 writes). Could batch via a single Up[Ng-1, :] = (rho, 0.0, rho_E) assignment. Realistic win: ~1-2 % wall on engine cases.
Phase Perf-14 — cylinder._kinematics slider-crank recompute caching. 52 K calls / 0.142 s self on case 50; the slider-crank position formula is a few transcendentals. When called multiple times with the same θ (during a single Newton stage's residual + Jacobian assembly) the result could be cached. Speculative — needs profile confirmation that the multi-call-per-θ pattern actually exists.

The unifying theme: the wrapper-hoist + getattr arc closed the "obvious" Python overhead. The remaining hot Python is property dispatch and small per-call allocations — the classic Python-on-top-of-numba death-of-a-thousand-cuts. Eliminating each ~1-2 % at a time is the natural next phase shape; aggregate win across Perf-10..14 would be ~5-10 % wall on engine cases.

The JIT kernels themselves (_packed_residual_fo_ser, species_advection_extras_jit, pipe_residual_quasi1d_first_order, hllc_flux) own ~25 % of case-50 wall and ~17 % of case-41 wall. Further kernel-level wins (SIMD review, prange at engine scale) remain on the BACKLOG but the prior Perf-7 SIMD HLLC rewrite was reverted as a 16-33× regression — speculative kernel work needs extreme bench-first discipline.

Phase Perf-10 — Pipe.has_species → cached field¶

Trigger: Phase Perf-history-3 profile of case 50 flagged Pipe.has_species at 992 K calls / 0.152 s self / 2.3 % of wall — a property that just evaluates self.Y.size > 0. Top of the next-tier candidate list.

Fix

Replaced the property body with a single attribute read of a new cached field _has_species: bool (declared on the dataclass with init=False, default=False). The default matches the dataclass default Y = zeros(0) (false). The cache is refreshed in set_species — the only sanctioned mutator that turns species transport on; no other code path in the codebase reassigns self.Y.

Result (back-to-back cProfile A/B, case 50)

Function	PRE	HEAD	Δ
`has_species` self	0.152 s	0.095 s	−37 %
`pipe_fuel_key` cum	0.220 s	0.183 s	−17 %

The 37 % function-self-time drop is real but smaller than naively expected (a cached attribute should be near-free). The remaining 0.095 s is the property dispatch overhead itself — Python still walks the descriptor protocol on every pipe.has_species access. Going further would require turning has_species into a bare attribute (not a property) and removing set_species as the sanctioned mutator. Out of scope for this phase; the property→cache change is the right-shaped first step.

Test fixture follow-on

Two unit-test fixtures (test_egr_valve_controller.py and test_fuel_injector_controller.py) bypassed set_species and assigned pipe.Y directly. With the cache, this left has_species stale at False even though Y was set. Rewritten to go through set_species — the right pattern anyway. Caught by the test suite, fixed in lockstep.

Phase Perf-11 — Pipe.U property → cached slice view¶

Trigger: Phase Perf-history-3 catalog item #2. Pipe.U was a property re-computing self.U_pad[Ng:Ng+N] on every call; 249 K calls / 0.086 s self on case 41 + 392 K / 0.120 s on case 50. The slice itself is a lightweight numpy view, but the create-view-on-each- call cost adds up at hundreds of thousands of accesses per run.

Fix

Added _U_view: np.ndarray as a cached field on the Pipe dataclass (init=False). Populated once in __post_init__ after U_pad is allocated; the U property reads self._U_view directly. The setter writes self._U_view[:] = value (slice-write, propagates to U_pad since the view shares the same buffer).

The cache stays valid because U_pad is never reassigned post-__post_init__ — every downstream mutation goes through U_pad[...] = ... slice writes that propagate through the view automatically.

Result (back-to-back cProfile A/B, case 41)

Metric	PRE	HEAD	Δ
`Pipe.U` self-time	0.086 s	0.022 s	−74 %
Total cProfile wall	2.402 s	2.343 s	−2.5 %

The −74 % is the bigger metric — slice creation is fully eliminated; only the property-dispatch overhead (descriptor walk) remains. Going further would require turning U into a bare attribute (drop the property + setter) which would break the existing pipe.U = value API the test suite uses.

Phase Perf-12 — pre-allocated scratch buffers in regen_into¶

Trigger: Phase Perf-history-3 catalog item #3. aeolus1d.closures.filter.regen_into allocated 5 × np.zeros(N) per call (Y_o2 / Y_no2 / Y_no / d_soot_dt / dT_cake_dt), totalling 100 K allocations per case 41 run / 0.028 s self-time on the np.zeros calls.

Fix

Added 3 pre-allocated scratch buffers (_scratch_Y_o2 / _scratch_Y_no2 / _scratch_Y_no) to FilterState (init=False, lazily sized on first call when size != N; subsequent calls reuse via slice-write).

Only 3 of the 5 zero arrays get pre-allocated:

d_soot_dt and dT_cake_dt are RETURNED to the caller for integration into the volume's soot mass + cake temperature across the timestep. Pre-allocating them on FilterState would require a copy-on-return, defeating the savings (the caller holds onto the array for the next RK stage).
Y_o2 / Y_no2 / Y_no are scratch — read in, passed to the JIT kernel, never returned. Safe to reuse.

Result (back-to-back cProfile A/B, case 41)

Metric	PRE	HEAD	Δ
`regen_into` self	0.159 s	0.149 s	−6 %

Modest because only 3 of the 5 zeros were eliminated; total np.zeros self-time dropped from 0.028 s to ~0.011 s (the remaining 60 K allocations are d_soot_dt + dT_cake_dt, which can't be cached for the reasons above).

This closes the Perf-history-3 catalog's "easy wins". Items

13 (`_fill_plenum_ghost` per-cell loop fusion) and #14¶

(_kinematics result caching) need careful bench-first investigation before claiming wins; deferred.

Phase Perf-13 — _fill_plenum_ghost loop fusion (REVERTED)¶

Trigger: Perf-history-3 catalog item #4. The _fill_plenum_ghost Python wrapper was 191 K calls / 0.332 s self / 4.5 % of case-50 wall, doing 3 scalar item-stores per ghost row in an inner for g in range(Ng): loop.

Hypothesis: fuse the 3 scalar writes into a single slice assignment Up[:Ng] = (rho, 0.0, rho_E) — broadcasts the tuple to the (Ng, 3) slice. For Ng=1 (the common case) this replaces 3 scalar item-stores with one row-write.

Result: REGRESSION

A/B back-to-back on case 50:

Metric	PRE	HEAD	Δ
`_fill_plenum_ghost` self	0.332 s	0.363 s	+9 %

The slice-assign + broadcast machinery costs more than 3 scalar item-stores at the tiny (1, 3) shape numpy operates on here. Numpy's broadcast / typed-store path has fixed Python- side overhead per call that's larger than 3 raw scalar writes.

Reverted in the same commit; the original loop-of-3-scalar- stores stays. Documented inline so a future contributor doesn't re-attempt the same naive fusion.

Lesson reinforced: profile-driven optimisation cuts both ways — the catalog item flagged this as a candidate, but bench-first confirmed the "obvious" fusion is wrong at this scale. The Perf-7 SIMD HLLC rewrite (16-33× regression) was the canonical lesson; this is a smaller-scale repeat. The profile_regress gate would have caught this if the case-50 _fill_plenum_ghost ceiling weren't 3 × the post-Perf-6r4 baseline (loose enough to absorb a 9 % overshoot). Tightening the ceiling to e.g. 1.5 × would catch this kind of regression automatically.

Phase Perf-history-4 — re-profile post-Perf-10/11/12¶

Investigation only (no code changes). Re-ran scripts/profile_perf1.py on cases 41 + 50 after Perf-10 (has_species cache), Perf-11 (Pipe.U cached slice), and Perf-12 (regen_into scratch buffers) landed. Confirms the previous wins are sticky and identifies the new top-of-hot-list for Perf-15+.

Case 41 (chemistry/DPF) — 2.32 s wall:

Post-Perf-10/11/12, the chemistry hot list is dominated by:

Function	Self	% wall	vs Perf-history-3
`_species_rate`	0.240 s	10.4 %	+0.011 (~stable)
`regen_into`	0.151 s	6.5 %	−0.003 (Perf-12)
`species_advection_extras_jit` (JIT)	0.120 s	5.2 %	~stable
`pipe_residual_quasi1d_first_order` (JIT)	0.121 s	5.2 %	~stable
`cake_into`	0.052 s	2.2 %	~stable
`Pipe.U` property	0.023 s	1.0 %	−0.059 (−72 %) ✓
`Pipe.n_extra_species`	0.040 s	1.7 %	~stable

Phase Perf-11's win is sticky and well above noise. The remaining hot Python is now mostly the wrappers themselves (_species_rate, regen_into) — further Python-level wins require restructuring the wrappers, not just hoisting attributes.

Case 50 (engine network) — 6.87 s wall:

Function	Self	% wall	vs Perf-history-3
`_packed_residual_fo_ser` (JIT)	0.974 s	14.2 %	~stable
`junction.accumulate_rates`	0.357 s	5.2 %	~stable
`junction._pipe_port_flux`	0.362 s	5.3 %	~stable
`junction._fill_plenum_ghost`	0.309 s	4.5 %	~stable
`junction._ghost_rho_E`	0.271 s	3.9 %	~stable
`cylinder.accumulate_rates`	0.247 s	3.6 %	~stable
`cylinder.apply_port_bcs`	0.218 s	3.2 %	~stable
`junction.apply_port_bcs`	0.187 s	2.7 %	~stable
`hllc_flux` (JIT)	0.169 s	2.5 %	~stable
`species_coupling.pipe_fuel_key`	0.142 s	2.1 %	~stable
`cylinder._kinematics`	0.138 s	2.0 %	~stable
`bc.base.wall_bc`	0.113 s	1.6 %	NEW IN TOP-15
`cylinder._port_area_fraction`	0.110 s	1.6 %	NEW IN TOP-15
`Volume.p` property	0.105 s	1.5 %	NEW (296 K calls)
`Pipe.has_species` cache	0.096 s	1.4 %	−0.056 (−37 %) ✓

The Perf-10 cache win is sticky. New candidates surfaced as the property-dispatch ceiling lifts:

Future-work candidates (Perf-15+):

Phase Perf-15 — Volume.p property → cached field. Same pattern as Perf-10/11. 296 K calls / 0.105 s self on case 50. Volume.p is computed from m and E via gamma; cache when m / E_internal change (which happens once per step in apply_fuel_step / Volume update). Realistic ~1.5 % wall.
Phase Perf-16 — cylinder._port_area_fraction binary-search → linear-scan tweak. 104 K calls / 0.110 s self. Currently uses np.interp on a CAD-domain lift schedule; for typical ~7-point schedules a hand-coded linear scan can be faster. Speculative — needs bench-first.
Phase Perf-17 — bc.base.wall_bc per-port hoist. 69 K calls / 0.113 s self. Same hoist pattern as the Perf- 6r3..r6 arc but for BC callables. ~1 % wall.
Phase Perf-18 — cylinder._kinematics result caching across Newton residual + Jacobian. Perf-history-3 item #5. Speculative — needs to verify the multi-call-per-θ pattern actually exists.

Aggregate Perf-15..18 win estimated at ~3-5 % case-50 wall. Diminishing returns — every subsequent Perf phase yields a smaller relative win as the Python overhead surface shrinks.

Lesson recorded — Perf-13 revert:

Perf-13 attempted to fuse _fill_plenum_ghost's 3 scalar item-stores into a single slice-assign and regressed +9 %. The slice-assign + broadcast machinery costs more than 3 scalar writes at the (1, 3) target shape. This is the second fusion-style optimisation reverted under bench-first discipline (after Perf-7 SIMD HLLC, which regressed 16-33×). Documented inline in junction.py so a future contributor doesn't re-attempt the same naive fusion.

Phase Perf-15 — Volume.p cached field (deliberate no-op)¶

Trigger: Perf-history-4 catalog item #1. Volume.p was 296 K calls / 0.105 s self / 1.5 % of case-50 wall — a property that's structurally identical to Pipe.has_species (Perf-10) and Pipe.U (Perf-11) and naïvely the same caching pattern should apply.

Why it doesn't: the Volume class differs from the Pipe class in one critical way for caching:

Pipe.U is a slice into Pipe.U_pad — a single mutable buffer that's never reassigned post-init. The cache stays in sync automatically because U_pad[...] = ... slice writes propagate through the view.
Pipe.has_species flips on exactly once (in set_species) and never flips off — the cache has a single sanctioned mutator.
Volume.p is COMPUTED from E_internal / m / V / gamma, with E_internal and m mutated by ~8 external call sites (each RK stage in network_timeloop + implicit_network + 3 in helmholtz_resonator).

A cache for p would need explicit invalidation at every E_internal = ... / m = ... site OR turning those into properties with cache-invalidating setters. The first is fragile (any new mutator forgetting the call silently desyncs); the second adds property-dispatch overhead to the ALREADY-hot E_internal / m reads inside p itself (net dispatches per step likely INCREASES, not decreases).

Decision: no change.

The 1.5 % wall savings is too small to justify the maintenance burden — every future code path that mutates Volume.E_internal or Volume.m would need to remember to invalidate, and the bug pattern (cache desync produces silently wrong pressures) is exactly the kind of correctness hazard the project explicitly avoids. Documented as a deliberate no-op so a future contributor doesn't re-investigate.

The right shape for a Volume.p win would be a hook the timeloop calls once per step (after the RK combine sets the final E_internal) to refresh the cache. That requires threading a "dirty" flag through 8+ sites with a single sanctioned update point — feasible but scope-creeps into a Volume-state-management rewrite. Deferred indefinitely; revisit only if profile shows Volume.p climbing past ~5 % of wall.

Pattern reaffirmed: cache invalidation is the hardest problem in computer science; Perf-10/11 worked because their shared-buffer / single-mutator structure made the invalidation trivial. Perf-15's structure isn't trivial and the win isn't big enough to pay for the complexity.

Phase Perf-17 — bc.base.wall_bc per-port hoist¶

Trigger: Perf-history-4 catalog item #3. bc.base.wall_bc was 70 K calls / 0.114 s self / 1.6 % of case-50 wall, doing 4 pipe.U_pad attribute walks per ghost row (twice on the LHS of the indexed assignment, once on the RHS, once on the momentum flip).

Fix

Hoisted Up = pipe.U_pad once at the top of wall_bc; the inner loop body uses Up[...] instead of pipe.U_pad[...]. Same shape as the Perf-6r3..r6 wrapper- hoist arc.

Result (back-to-back cProfile A/B, case 50)

Metric	PRE	HEAD	Δ
`wall_bc` self	0.114 s	0.112 s	−2 %

Marginal — within noise. With Ng=1 (the common case) the inner loop body runs once with ~3 attribute walks saved per call; total saving ~0.02 s at 70 K calls × 0.3 µs/walk. The cProfile delta is at the edge of what cProfile resolves at this scale.

Worth keeping anyway for code-cleanliness reasons (matches the Perf-6r3..r6 pattern that the rest of the codebase uses). Bench-first discipline kept this in scope (the Perf-13 fusion attempt taught the lesson — "obvious" small-shape optimisations can regress). Hoisting an attribute reference is strictly cheaper than the original at the assembly level; no risk.

Phase Perf-7 — packed-kernel parallel-vs-serial dispatcher¶

Commit: bfd3146. Trigger: Phase Perf-history-2 profile said engine cases (50 = 6 fast pipes, 52 = 2 fast pipes) were 99 % JIT-kernel-bound. The candidate optimisations from BACKLOG were (a) SIMD/scalarisation rewrite of the HLLC inner loop, (b) prange over pipes at engine- network scales. (b) was already in place — the packed kernels were @njit(parallel=True). The actual question turned out to be the inverse: is the existing parallel kernel even worth it at engine scales?

Investigation

Bench at varying NUMBA_NUM_THREADS (audit laptop: 16 cores, default = 16 threads):

Case (P)	threads=1	threads=2	threads=8	threads=16
50 (P=6)	7.0 s	5.95 s	6.26 s	15.5 s
52 (P=2)	0.91 s	0.89 s	1.05 s	1.85 s

The default 16-thread setup was the worst case for typical aeolus1d network sizes. At P=2 or P=6, splitting the work over 16 threads means 10–14 idle threads paying coordination costs on every kernel call (and the kernels are called once per RK stage, thousands of times per simulated second).

Tried-and-rejected: SIMD/scalar rewrite

Wrote hllc_flux_scalar(UL0, UL1, UL2, UR0, UR1, UR2, gamma) returning (F0, F1, F2) floats, plus max_wave_speed_scalar, both inline="always". Refactored _packed_residual_fo and _packed_min_dt to use them — eliminates the per-cell np.empty(3) allocations.

Numerical regression: 1 element of 90 in test_packed_residual_matches_sequential drifted by abs 9.96e-8 (rel 2.8 %). Fastmath FMA reordering — the original kernel used for k in range(3): out[k] = F[k] + SL*(...) which the parallel pass vectorised across k; the scalar rewrite unrolled it as three separate FMAs which the compiler ordered differently.
Wall-clock regression: case 50 jumped from ~7 s to 131 s (16 × slower) and case 52 from ~1 s to 43 s (33 × slower). Suspect: the inlined scalar HLLC body multiplies the prange parallel-loop body size dramatically (control flow + branches), and at engine-scale pipe counts the parallelisation pass either degrades to a worse code shape or interacts badly with thread-local cache pressure.

Reverted entirely. The scalarisation idea is sound but needs a no-parallel=True baseline (Phase Perf-7 default).

What landed

Two new @njit kernels in aeolus1d/solver/packed_pipes.py:

_packed_residual_fo_ser — serial twin of _packed_residual_fo_par (same body, plain range(P), no parallel=True).
_packed_min_dt_ser — serial twin of _packed_min_dt_par.

Plus a Python dispatcher that picks parallel vs serial by fast-pipe count:

_PARALLEL_THRESHOLD_P = 12  # default

def _packed_residual_fo(U_pad_all, ...):
    if U_pad_all.shape[0] >= _PARALLEL_THRESHOLD_P:
        _packed_residual_fo_par(...)
    else:
        _packed_residual_fo_ser(...)

Override knob: set_parallel_threshold(P) — set to 0 to always use parallel (legacy pre-Perf-7 behaviour) or to a very large number to always use serial.

The threshold of 12 was picked from a synthetic many-pipe microbench: parallel cleanly wins at P ≥ 12 with N=80; at P < 12 the serial variant is consistently within 5 % at best and 16 × better at worst.

Wall-clock

Workload (P)	Baseline (post-Perf-6)	Phase Perf-7	Speedup
50 (engine, P=6)	16.78 s	6.24 s	2.69 ×
52 (knock, P=2)	2.07 s	0.95 s	2.18 ×
54 (turbo, P=4)	12.14 s	12.28 s	unchanged†

† Case 54 uses advance_network_implicit, not the packed- kernel path in advance_network, so Perf-7 doesn't touch its hot path. The variance reading is system noise.

Tests

test_packed_dispatcher_serial_and_parallel_agree — pins numerical equivalence (rtol=1e-13 — fastmath FMA reordering between parallel and serial kernels is below that threshold) across P ∈ {1, 4, 12, 16}.
test_packed_dispatcher_threshold_picks_kernel — pins the dispatcher routes P < threshold to the serial kernel and P >= threshold to the parallel kernel.
All 7 existing packed-pipes regression tests still pass.

Why not just clamp numba.set_num_threads?

set_num_threads is global and persistent — would affect any other parallel numba code in the same process (third-party, GUI plotting, future kernels). The dispatcher approach is local to packed_pipes and composable.

Phase I-8 (MVP) — opt-in Jacobian-Free Newton-Krylov¶

Commits: 8de3ac3 (MVP) + [after] (block-Jacobi preconditioner). Trigger: the direct block-tridiag-Thomas Newton solve in _newton_step needs an explicit (N, 3, 3) Jacobian per iter. For users who push to networks where the direct Jacobian doesn't fit, JFNK is the algorithmic alternative.

What landed

_jfnk_step(pipe, U_hist, dt_eff, F, apply_ghosts, config, scheme, slopes_pad, precond_D_J=None) -> δ solves J · δ = -F via GMRES (scipy.sparse.linalg.gmres) against a LinearOperator whose matvec does the FD J·v product (Brown-Saad ε scaling).
Right-side block-Jacobi preconditioner (Phase I-8-precond): when jfnk_precondition = True (default on when JFNK is on), the outer Newton loop still builds the chord-aged D_J (the same _assemble_jacobian_first_order call the direct path makes) and JFNK uses it via per-cell 3×3 inverse products (_block_jacobi_apply_jit). Without preconditioning GMRES typically saturates jfnk_max_iters on stiff cases; with it, Krylov dim shrinks dramatically.
NewtonConfig knobs: use_jfnk (False), jfnk_precondition (True), jfnk_tol (1e-3 GMRES rtol), jfnk_max_iters (50 Krylov-dim cap). All four surface through CaseHeader for TOML.

Tests

test_default_use_jfnk_is_off — pinning the opt-in contract.
test_jfnk_path_runs_and_state_is_finite — short transient finiteness.
test_jfnk_matches_direct_on_short_transient — JFNK vs direct converge to engineering-equivalent state at rtol=1e-3.
test_jfnk_unpreconditioned_skips_jacobian_assembly — patched-assembly contract: with jfnk_precondition=False, any call into _assemble_jacobian_first_order makes the test fail. Pins the pure-matrix-free path for memory-bound users.
test_jfnk_preconditioned_uses_jacobian_diag_only — smoke test that the preconditioned path runs and state is physical.

Wall-clock

JFNK is slower than the direct block-tridiag-Thomas path on the canonical case 54 — at small N × 3 unknowns the direct solve is the right cost shape. The MVP exists for users hitting memory walls; default stays off. Future follow-ups: coarse-grid preconditioner, MUSCL + variable-area branches, restart-GMRES tuning. Documented in BACKLOG.

Long-pipe bench (Phase Perf-history follow-up)¶

Script: scripts/bench_long_pipe.py. Trigger: the audit converged on N = 80 cells/pipe, but several rejected options had break-even points at larger N that we couldn't measure with the bundled examples. The script builds a programmatic single-pipe case at user-controllable N + t_end and times three Newton-solve paths:

direct (chord) — block-tridiag-Thomas with chord-Jacobian reuse (the default).
jfnk_precond — preconditioned matrix-free (use_jfnk=True, jfnk_precondition=True).
jfnk_pure — unpreconditioned matrix-free (use_jfnk=True, jfnk_precondition=False).

Measurements

N	direct	jfnk_pure	jfnk_block_jacobi	jfnk_ilu0
80	0.010 s	0.027 s (2.6 ×)	3.24 s (312 ×)	0.023 s (2.2 ×)
500	0.227 s	0.423 s (1.9 ×)	112.6 s (495 ×)	0.398 s (1.8 ×)

Findings

jfnk_pure is consistently ~2 × slower than direct, roughly N-independent. That's a reasonable competitive position for a memory-saving fallback — users who hit memory walls can opt in and pay 2 × on wall clock, no Jacobian assembly.
jfnk_ilu0 is the cleanest matrix-free path: ~1.8 × slower than direct at all N, slightly better than pure JFNK. The preconditioner reuses _block_thomas_solve on the chord- aged Jacobian — for block-tridiag J the ILU(0) factorisation is the block-Thomas factorisation. GMRES typically converges in 1–2 Krylov iters; the residual cost is what separates JFNK from direct, not the linear solve.
jfnk_block_jacobi is dramatically slower at every N — and gets worse with N (312 × at N = 80, 495 × at N = 500). Direct GMRES iter count goes from 4 to 42 when block-Jacobi is enabled (10 × the matvec count), confirming that the preconditioner is the wrong shape for the system.

Why block-Jacobi is wrong here

The implicit Newton system is J = I − dt · ∂R/∂U on a 1D hyperbolic problem. ∂R/∂U has block-tridiagonal structure where the off-diagonal blocks (L and Up) carry the acoustic + convective coupling — for cfl > 1 (the whole point of the implicit driver) those off-diagonal blocks are a significant fraction of J's spectral content. Block-Jacobi discards them entirely; the preconditioned operator has worse conditioning than the unpreconditioned one because the diagonal blocks alone don't span the dominant eigenmodes.

This is a known result in the literature for 1D hyperbolic systems (see e.g. Saad's "Iterative methods for sparse linear systems", §10 — preconditioning for non-symmetric systems where directional bias matters). The empirical bench just confirms it on this solver.

Decision

Default jfnk_precondition = False (was True in the MVP — flipped after this bench landed).
Keep the implementation in (it's opt-in; users who want to experiment with block-Jacobi or who have a less-hyperbolic workload — e.g. heavily porous-Forchheimer — can flip it on).
Future preconditioner work needs a directionally-aware approximation: line preconditioning, ILU(0) on the actual block-tridiag, or coarse-grid via aggregation. Documented in BACKLOG.

Phase 27c soak test¶

Commit: [after]. Trigger: the Phase 27c MVP test runs only a 2-pipe + plenum case for 2 ms; not enough to expose race conditions or measure realistic wall-clock under load.

What landed

tests/integration/test_parallel_pipes_soak.py adds two @pytest.mark.slow tests:

test_parallel_pipes_soak_case_54 — load examples/54_adaptive_picard_turbo.toml end-to-end (4 pipes, 2 turbos, adaptive Picard, full BC stack), run twice (serial + parallel), assert state bit-equivalence to rtol=1e-9, print both wall-clock numbers.
test_parallel_pipes_soak_short_repeats — run the parallel path 5× on a 2 ms transient and assert every repeat produces bit-identical state. Catches thread-scheduling-dependent non-determinism.

Soak measurement (full case 54 t_end = 50 ms)

serial wall: 7.2 s
parallel wall: 18.8 s
ratio: 2.6× slower on parallel

Stronger than the 2× the short cProfile run showed. Confirms parallel_pipes = False is the right default; the path stays for genuinely large networks where per-pipe JIT work clearly outweighs thread coordination overhead. Bit-equivalence and repeat-determinism both pass — no race conditions detected.

Canonical-green baseline + xdist contention¶

After all perf phases landed, ran make test-all to confirm the post-Perf state is clean. Result: 1714 passed, 50 failed, 4 skipped, 1 xfailed in 99 minutes under the default -n auto (8 workers on the audit-laptop). The 50 "failures" turned out to be pytest-xdist contention flakes, not regressions — every spot-checked failure passes when run alone (3 sampled tests across 3 different cluster files all green in 50 s combined, no code changes needed).

Root cause is most likely numba's AOT JIT cache: when 8 workers compile the same kernel concurrently, the cache file locking can race and emit corrupted/stale code that fails silently on whoever loses the race. The 3 % flake rate is roughly proportional to the JIT-heavy fraction of the suite × 8-way contention.

Mitigation — make test-stable (new target, commit [after]) runs pytest -n 4 --timeout=600 instead of -n auto. Half the workers, contention drops below the observable threshold; serves as the canonical-green confirmation target. make test-all stays as the historical -n auto target for CI environments where contention is acceptable.

The codebase is clean post-Perf. No real regressions introduced across the 11-phase optimization pass + the 50+ test files this work touched.

Things that were tried and rejected (or never attempted)¶

The audit is deliberate about what NOT to chase:

prange over cells in the Jacobian kernel — see Phase 27c-cells above. Slower than serial at N = 80.
Defaulting parallel_pipes = True — see Phase Perf-4. Serial post-Perf-4 is faster than parallel on the canonical benchmarks. Stays opt-in for genuinely large networks.
JIT'ifying _capture_ghost_snapshot / _ghost_snapshot_max_rel_change — Phase Perf-5 audit. Together they're 1.7 % of case 54 wall; not worth the diff-test relaxation cost.
Memoizing the bilinear cell index (audit item #10) — Phase Perf-1. _bilinear is inline="always"; the index computation is ~10 float ops fully inlined. Memoization would add cache-line traffic to skip ~5 ns of compute. Likely net-negative.
CFL min_dt per RK stage (audit item #9) — verified the audit was wrong; _packed_min_dt runs once per step, not per stage.
BC FD copy elimination (audit item #6) — paths are @njit; numba handles the copy efficiently. The Python-side .copy() in _ghost_state_and_chain is BC_NUMERICAL fallback, rarely hit.

Tools that made this possible¶

scripts/profile_perf1.py — pre-warms numba caches, then runs the case under cProfile and dumps top entries by cumulative and self-time. Without this, every "is this faster?" question would need a separate one-shot script.
scripts/bench_perf1_tabulated.py — programmatic tabulated-EOS bench for cases the TOML schema can't construct (was needed before Phase Perf-1-toml; kept around for users who want to A/B without CoolProp).
Bit-exact regression tests for each EOS arm and Jacobian path. Several optimizations would have shipped with subtle numerical drift if not for the tight assertions; the few that we relaxed (to assert_allclose) were reasoned about explicitly and limited to fastmath FMA reordering.
@pytest.mark.slow benches (e.g. the closure-Jacobian analytic-vs-FD bench from Phase B19) — these were the initial sanity check that the Perf-1 batched helpers didn't regress numerical accuracy. Should be re-run after every major refactor.

What the audit didn't move¶

The flow-of-control side of the implicit driver (advance_network_implicit, _newton_step outer body) is already efficient post-Perf-2. Further wins on the canonical 4-pipe network case would require algorithmic change:

Phase I-8 matrix-free JFNK (BACKLOG.md) — replaces the direct block-tridiag solve with a Krylov method that only needs J·v products. Worth exploring only if a bench case shows the direct Jacobian assembly running out of memory or time on a real network.
prange over pipes inside the Picard iter — Phase 27c with all-or-nothing dt-halving has landed as opt-in. The parallel path is currently slower than serial on the canonical 4-pipe case; revisit when a network has ≫ 4 pipes with cost-balanced per-pipe work.
prange over cells in the Jacobian — Phase 27c-cells. Replay when N ≫ 200–500 cells per pipe.

v0.3 aftertreatment time-stepping spike (2026-05-20)¶

The first task of the GT-POWER-parity roadmap's v0.3 keystone (see BACKLOG.md → "GT-POWER parity roadmap"). The question: for an explicit acoustic-CFL aftertreatment run, is the wall-clock dominated by the step count (→ the low-Mach / quasi-steady mode is the lever) or the per-step cost (→ JIT / kernel work, or a compiled core, is the lever)? Tool: scripts/profile_at_spike.py on 41_dpf_regen (DPF regen + cake thermal + chemistry — the bundled analogue of AVL Cruise-M C02006). Laptop run, numba caches pre-warmed.

Step economics (case 41, t_end = 0.06 s)¶

metric	value
steps	10 049
mean Δt	5.97 × 10⁻⁶ s (acoustic-CFL bound)
wall / step	0.199 ms
profiled wall	1.996 s

Self-time by solver phase¶

phase	self-time	share
orchestration (Python loop)	1.228 s	61.5 %
chemistry + closures (Python: `filter.regen_into` / `cake_into`)	0.247 s	12.4 %
numpy / builtins	0.150 s	7.5 %
flux kernels (JIT)	0.117 s	5.9 %
species advection (JIT)	0.106 s	5.3 %
property accessors (`pipe.*`)	0.073 s	3.6 %
boundary conditions	0.047 s	2.3 %
numba dispatch	0.026 s	1.3 %

500 s drive-cycle extrapolation (current per-step cost)¶

mode	Δt	steps	wall
acoustic-CFL (today)	5.97 µs	8.37 × 10⁷	4.62 h
low-Mach (Δt = 1 ms)	1 ms	5.0 × 10⁵	1.66 min (≈ 167× fewer steps)

Conclusions (these set the v0.3 ordering)¶

Step count is the existential factor, not per-step cost. The 4.6 h drive cycle is entirely an acoustic-CFL step-count problem. The low-Mach / quasi-steady AT mode is unambiguously the #1 performance lever — even a conservative 100 µs low-Mach step is ~17× (4.6 h → 16 min); a 1 ms step is ~167×. The HRR adaptive-Δt limiter then sets the achievable low-Mach step at ignition.
Per-step cost is ~74 % Python, only ~11 % JIT kernels. Orchestration (61.5 %) + Python chemistry closures (12.4 %) dominate; the flux + species-advection JIT kernels are already lean. This positively validates the roadmap's "defer the compiled (Rust/PyO3) core" call — the kernels are not the bottleneck. The secondary lever for the remaining low-Mach steps is JIT-ifying the orchestrator (Phase Perf-2, already backlogged) and the closures/filter.py chemistry closures (regen_into / cake_into are pure-Python).
Free hot-loop win: pipe.n_extra_species is recomputed 148 495× and pipe.U 249 626× over the run — recomputed properties inside the step loop. Hoisting both is a zero-risk micro-win (~0.1 s / 5 % here, larger on longer runs).

Re-run after the low-Mach mode lands to confirm the step-count collapse and re-profile the new per-step mix.

Quasi-steady SIMPLE kernel — numba JIT (2026-05-21)¶

The pure-NumPy solver/quasi_steady.py SIMPLE solve was ~3.3 ms/solve (profiled on a 50-cell C02006-like pipe, ~34 SIMPLE iterations); the scalar Thomas tridiagonal solve alone was ~48 % of that. A drive-cycle run (hundreds of thousands of macro steps × one solve each) was impractical.

_simple_solve_ideal_jit is a faithful @njit port of the ideal-gas SIMPLE loop (Thomas inlined; Rhie-Chow + pressure-correction unchanged). The non-ideal-EOS path stays in the NumPy loop (needs a per-cell density inversion). Result: ~3.3 ms → 0.044 ms/solve, ~75×. Physics-identical — all 21 quasi-steady tests pass against the analytic Darcy solution + the network mass balance.

AVL C02006 regeneration — end-to-end (accelerated-ramp validation)¶

With the JIT kernel + the transient-BC / gas-energy-advection work (feat/qs-transient-bc) + a gas-temperature floor (t_floor, guards the explicit energy update against a regen-ignition negative-T excursion), the C02006 DPF-CSF regeneration runs to completion and the regen fires: with the inlet ramp accelerated to 700 K by t = 2 s (same physics, short run), soot drains 10 → 0 kg/m³ (100 % burn-down), the run is stable (gas-T floored at 250 K, peaks ~1900 K at ignition), where the density-based explicit solver diverges at t ≈ 116 s.

Remaining perf wall (BACKLOG follow-up): the explicit energy advection + T-floor leave a noisy 250–1900 K field at ignition that keeps the HRR-controlled Δt tiny (~620 k steps for 15 s), so the full 600 s cycle is still impractical. The principled fix is the steady gas-phase energy solve (GT's quasi-steady catalyst gas phase) instead of explicit enthalpy advection — that removes the advective-CFL on energy and the ignition noise. The other per-step cost is now the Python/NumPy orchestration (_species_rate wrapper, energy advection, regen, ghost-fill) — the QS solve is no longer the bottleneck.

Phase Perf-7r — turbo map bilinear lookup (2026-05-23)¶

Trigger: Profile of canonical case 54 (54_adaptive_picard_turbo) identified TurboMap._lookup_2d cumulative time at 4.80 s out of 29.79 s total cProfile wall (~16 %), driven by 497 000 np.interp calls — the implementation did one 1-D np.interp per speed line along PR (~14 calls), then a second 1-D interp along N_corr (1 call), and another pass for η_is. That's ~22 numpy calls per _lookup_2d invocation when a single bilinear interp would do.

Change. Replaced the nested np.interp pattern in both _lookup_2d and _lookup_slice (the 2-D and VGT 3-D paths) with a shared _bilinear_2d helper: two np.searchsorted calls (PR + N_corr bracket indices), four corner reads, and a weighted blend. Algebraically identical to the nested 1-D interps — two nested 1-D linear interps are bilinear; this just folds the arithmetic so Python only takes two function-call hits instead of 22.

Result (cProfile, case 54, default t_end):

	PRE-Perf-7r	HEAD	Δ
`TurboMap._lookup_2d` cumulative	4.80 s	0.83 s	−83 %
Total `np.interp` calls (case 54)	497 000	0	−100 %
Case 54 cProfile wall	29.79 s	10.64 s	−64 %
Case 54 wall-clock (no cProfile, t_end=5 ms)	—	4.79 s	—

Speedup: 2.80× under cProfile (Python-overhead-amplified) on case 54. Behaviour is bit-equivalent modulo FP reordering — all 114 existing turbo tests pass without tolerance adjustment.

Top hotspots remaining on case 54: _assemble_jacobian_fo_ideal_uniform_jit (3.20 s, already njit — at compute floor) and _newton_step orchestration (0.99 s self). No remaining Python-overhead targets at this scale.

AVL drive cycle (MFR_transition_SC630) profiling — 2026-05-29¶

Goal: can the AVL quasi-steady aftertreatment case run the full 500 s under 10 min? Verdict: no — not via configuration or micro-optimization.

Baseline (t_end = 1.0 s, chemistry firing): 115 s wall, ~1835 macro steps, ~115× realtime → the full 500 s extrapolates to ~16 h. Hitting 10 min (600 s wall) for 500 s sim needs ~0.83× realtime — a ~96× speedup.

Bottleneck (cProfile): chemistry is ~88% of wall time. The quasi-steady surface-chem sub-step is dominated by scipy LSODA's per-substep Python loop (lsoda._step_impl, ~64 k substeps over a 0.3 s window) plus the Brinkmeier reaction Python wrapper; the @njit _brinkmeier_rates_jit kernel itself is a small fraction. Transport / Picard / steady-energy are minor.

Measured attempts (t_end = 0.3 s window, 17.1 s baseline):

lever	wall	macro steps	gas T_max	note
baseline (`qs_cfl_adv=2`)	17.1 s	498	4307 K	—
`qs_chem_every_n=4`	12820 s	2.25 M	1161 K	catastrophic — lagged HRR explodes the dt controller during light-off
`chem_rtol=1e-4`	17.1 s	498	4307 K	no effect (not tolerance-bound here)
`qs_cfl_adv=4`	8.8 s	236	4299 K	1.94×, physics preserved
`qs_cfl_adv=8`	4.4 s	98	4064 K	3.9×, ~6% T_max change (best reliable)
`qs_cfl_adv=12`	9.7 s	557	1160 K	unstable — coverage collapses
`qs_cfl_adv=16`	5.1 s	74	3606 K	runs, but accuracy degrades
`chem_parallel=true`	16.6 s	498	4307 K	marginal
Brinkmeier `reaction` float-hoist (code)	17.1 s	498	4307 K	numerically identical; no measurable wall gain (scipy LSODA's Python loop, not the casts, is the cost)

Conclusion: the only clean lever is qs_cfl_adv (fewer macro steps → fewer chem calls); 8 is the sweet spot (~3.9×, 16 h → ~4 h). That is ~29× realtime — still far from the ~96× the 10-min target needs. qs_chem_every_n super-cycling, the theoretical big lever, is catastrophic during light-off. Reaching 10 min requires a fundamentally different chemistry path — replacing scipy LSODA with a JIT'd in-kernel stiff integrator (BACKLOG: "JIT'd frozen-T Brinkmeier kernel"; the cProfile confirms scipy's Python overhead is the wall) or reduced-order / tabulated kinetics. Both are research efforts, not tuning.

Scoped: ISAT chemistry acceleration — 2026-05-29¶

The chem-integrator path is closed (LSODA's adaptive variable-order BDF is near-optimal; the banded multi-cell light-off system is genuinely stiff — measured njev ~1848, up to 260 Jacobian builds per chem-step — and three implicit rewrites were a net wash). The only way past the chemistry wall is to avoid the stiff solve. ISAT (in-situ adaptive tabulation) tabulates the reaction mapping on the fly and falls back to the real LSODA solve on a miss, so it is bounded-error by construction and can never produce garbage.

Phase-0 feasibility (measured on AVL MFR_transition_SC630, t_end=0.3): 9 960 per-cell chem queries logged. Effective dimensionality 8 of 15 (the queries lie on a low-dim manifold). A greedy covering in normalized state space needs only 259–572 table points to cover 93–97 % of all queries (ε = 0.2 / 0.1 / 0.05). → strong reuse; ISAT is viable. Caveat: this is an optimistic proxy — the real ellipsoid-of-accuracy (EOA) must bound the mapping error, not just state distance, so the accuracy-bounded hit rate will be somewhat lower.

What to tabulate: F: x = (T, Y_primary, Y_extras[~10], θ_Ce) → (Δx over Δt, q_rxn) plus the sensitivity matrix A = ∂x_new/∂x for the linear retrieve approximation. Δt is ~fixed per chem-step; include it (or rescale) to cover the small variation.

Algorithm (standard ISAT): query x → search a binary tree / kd-tree for a record whose EOA contains x. Retrieve (hit): x_new ≈ x0_new + A·(x − x0) (one matvec). Grow: real solve is within tol → expand the EOA. Add (miss): run the existing banded LSODA, compute A (integrated sensitivity or FD), insert a record. Misses reuse the current _qs_surface_chem_step as the handler, so error ≤ the EOA tolerance.

Integration point: the per-cell mapping inside _qs_surface_chem_step in quasi_steady_network.py; the existing banded LSODA becomes the miss-handler + A-matrix source. One table per catalyst pipe (or shared).

Phases + decision gates: - P0 (done): locality feasibility → GO (93–97 % coverable, eff-dim 8). - P1 (~3–4 d): minimal ISAT (EOA add/retrieve, kd-tree, LSODA fallback, FD sensitivity). Gate: accuracy-bounded hit rate > 80 % on AVL light-off. - P2 (~2 d): tune tolerance; validate the full 500 s run vs pure-LSODA (θ_Ce / q_rxn / tailpipe trajectory). Gate: ≥ 5× chem speedup at ≤ ~2 % trajectory error. - P3 (fallback): if the accuracy-bounded hit rate disappoints, pivot to an ML surrogate (NeuralODE / ChemKAN over the 8-D manifold P0 found) — higher ceiling, needs a training pipeline + generalization validation, fits the existing closures/learned_* infra.

Risks: accuracy-vs-distance gap (EOA must bound mapping error); table growth / memory; tree-search cost vs solve cost; thread-safety under parallel pipes; Δt variation. All mitigated by the LSODA fallback — ISAT is never wrong, only slower on misses.

Success criteria: ≥ 5× on AVL chem cost (chem is ~88 % of wall → roughly 3× on the full run, stacking on the qs_cfl_adv=8 ~3.9×), with θ_Ce + q_rxn + tailpipe within ~2 % of pure-LSODA over 500 s.

Why ISAT before the ML surrogate: graceful bounded-error fallback (no garbage, no extrapolation risk), no offline training pipeline, and the P0 locality data directly supports it. The surrogate is the higher-ceiling P3 option if ISAT's accuracy-bounded hit rate disappoints.

ISAT P1+P2 results — 2026-05-29¶

P1 — aeolus1d.closures.isat_chem.ISATChem: a generic in-situ adaptive tabulation with EOA add / grow / retrieve, a bounded-error real-solve fallback (misses/grows return the exact value; the EOA only grows where the approximation is verified within tol), and order-0 (constant) or order-1 (linear + FD sensitivity) retrieval. Unit-tested on a clustered random-walk stream: 95% of results within tol, mean within tol (the scalar-radius EOA is approximately bounded — a full ellipsoidal EOA is the strict upgrade).

P2 — validated on the real Brinkmeier+LSODA mapping (a spread of 40 cells, some lit-off some not, evolved 150 steps; ml_chemistry/validate_isat.py):

order	tol	hit rate	table	mean err	max err	speedup
0 (constant)	1e-2	79%	991	—	0.33	2.7×
0 (constant)	1e-3	44%	2098	—	0.33	1.2×
1 (linear)	1e-2	94.6%	42	8.8e-3	0.25	~3.5×

Takeaways: order-1 (linear) is the clear choice — it recovers the ~95% hit rate the Phase-0 locality predicted, with a tiny 42-record table and typical error within tol, and the FD-sensitivity cost on adds amortizes because adds are rare. The ~3.5× is on a 150-step window and understates the full-run figure — the table saturates early, so a 500 s run is mostly retrieves. The persistent max error ~0.25 sits on cells crossing the ignition threshold (chemistry is sensitive there); the typical (mean) accuracy is fine. Two refinements take it to the ≥5× / ≤2× target: an ellipsoidal EOA (strict per-direction bounds) and an ignition-aware tolerance (tighter EOA where dθ/dY is large). The ML surrogate (ml_chemistry/) is the alternative high-ceiling path.

ISAT vs the PRODUCTION baseline — NO-GO — 2026-05-29¶

The P2 ~3.5× above is measured against a per-cell exact solve. But production _qs_surface_chem_step does a single banded LSODA over all cells at once (vectorized; the chem RHS is block-diagonal, banded with lband=uband=n_state-1), integrating the full macro Δt in one adaptive call. Re-measuring against that right baseline overturns the win.

Probe (tests/profiling/isat_perf_probe.py): captures the real chem-call trajectory of aeolus_optimized.toml in-situ (the production kernel runs the actual advance; its banded output is the trusted reference — a faithful re-solve self-checks bit-identical, max|Δ|=0). Then replays the captured input trajectory through a batched ISAT (vectorized nearest-neighbour retrieve + one banded solve over only the missed cells, fixed substep h with an exact remainder), with interleaved per-call timing (baseline and ISAT timed back-to-back each call, so machine load hits both equally).

AVL TWC, 0.5 s capture: 20 cells/pipe, n_state=14 (1 primary + 12 extras + 1 θ_Ce), reorder on; dt 1e-4–1.6e-3 (median 7e-4, not settled — light-off transient); banded solve median 28 ms/call, LSODA nfev median 1152 (genuinely stiff). Baseline = 17.7–18.2 s / 600 calls.

h	hit rate	table	max err (tol 1e-4)	speedup vs banded
7.0e-4	78.9%	190	5.1e-2	0.67×
3.5e-4	88.2%	270	2.8e-2	0.74×
1.8e-4	93.5%	391	3.8e-2	0.69×
8.8e-5	93.8%	516	3.5e-3	0.56×

Batched ISAT is slower at every substep (0.56–0.74×) and misses the accuracy target (35×–500× tol). Mechanism: the banded LSODA already vectorizes all 20 cells into one adaptive call. To match its accuracy ISAT must sub-cycle at h ≪ dt (tens of thousands of batched sub-steps + thousands of small scipy solves for misses / remainders / FD sensitivity), whose overhead dwarfs 600 big adaptive solves. The variable HRR-controlled Δt (135 unique values) is hostile to fixed-step reuse, and the scalar-radius order-1 EOA collapses on the stiff ignition-threshold cells. ISAT pays off with many cells / un-vectorized per-cell solves; this is a 20-cell mesh where the baseline is already vectorized — so there is no per-cell cost to amortize, exactly why the per-cell-baseline 3.5× did not carry over.

Adversarially verified (5 independent skeptic lenses): impl-efficiency (numba/kd-tree/order-0 ceiling ~1.1–1.2×, can't close the gap), accuracy (35× overshoot is fundamental on the ignition manifold), and representativeness (16× dt spread floors the win) all uphold the no-go at 0.92–0.95 confidence. A load-fairness concern was raised and then disproven by interleaved remeasurement (baseline stable, same ratios).

Verdict: ISAT is closed for this solver/case — wiring it into _qs_surface_chem_step would be a regression. This joins the stiff-integrator-swap finding (also a net wash 4×): LSODA's vectorized banded adaptive solve is the right tool for this small-mesh stiff chem step. The remaining higher-ceiling lead is the ML delta surrogate (ml_chemistry/), but only conditionally — it needs (a) the run to settle to a stable Δt (the transient shows 135 unique Δt) so a fixed-Δt surrogate is in-distribution, and (b) a residual-check fallback holding ≤2× error; both unproven here. The proven wall-time lever stays qs_cfl_adv (fewer macro steps → fewer chem calls).

Phase A — fused @njit Brinkmeier rhs (2026-05-31)¶

MRs: !132 (fused rhs, opt-in), !134 (default flip to True). Schema knob: qs_chem_fused_rhs: bool = True (src/aeolus1d/io/schema.py line 349, with the full design comment at lines 331–348).

Trigger. The ISAT NO-GO above closed the chem-integrator-swap path: LSODA's vectorised banded adaptive solve is the right tool for this small-mesh stiff step, so we cannot replace it. But the post-mortem cProfile of _qs_surface_chem_step showed that scipy's per-substep Python rhs closure was the actual wall — not the C-level LSODA bookkeeping, not the @njit _brinkmeier_rates_jit kernel itself. The closure was a death-by-a- thousand-cuts pile: a dict comprehension to gather species, a reshape, a list-of-tuples extras reorder, a few np.ascontiguousarray casts, and a final dy.ravel, all called ~64 k times per chem-step window.

Change. A single fused @njit rhs, _qs_brinkmeier_rhs_fused, replaces the Python rhs whenever the pipe carries a Brinkmeier surface reaction (i.e. pipe.py_reaction_surface._jit_params is published). Behaviour is bit-equivalent to the legacy Python rhs — LSODA takes the same step sequence on identical rhs values — so the existing tolerance regime, integrator stats, and downstream θ_Ce / q_rxn / tailpipe trajectories are unchanged by construction. Non-Brinkmeier closures (generic builder, NH₃, py_reaction kinds) lack _jit_params and fall back to the Python rhs automatically — no per-case opt-out needed.

Validation. Bit-equivalence anchored by tests/integration/test_qs_chem_fused_rhs.py (runs the same case twice with qs_chem_fused_rhs toggled and asserts identical state at every macro step). Full unit + integration sweep green; AVL Cruise-M bench held the wall-time win across reboots.

Result — AVL Cruise-M MFR_transition_SC630, 12 s sim:

	wall	realtime ratio
pre-Phase-A (`qs_chem_fused_rhs=False`)	97.7 s	8.15×
Phase A (`qs_chem_fused_rhs=True`)	47.0 s	3.92×
speedup	2.08×	2.08×

Default flip rationale. Three independent signals authorised defaulting qs_chem_fused_rhs to True (MR !134):

Bit-equivalence test green on the full sweep → there is no numerically-different alternative for users to opt into; the Python rhs path is now strictly slower for identical output.
The non-Brinkmeier fall-through is automatic; nothing breaks for cases that don't carry _jit_params.
The 2.08× is large enough that hiding it behind an opt-in flag was net-negative for default UX. The flag remains exposed for False as a regression escape hatch.

Diagnostic script shipped: scripts/profile_post_fused.py — mirrors the legacy pre-fused profile script but with qs_chem_fused_rhs=True so the post-Phase-A breakdown is directly comparable. Reads local_cases/avl_cruise_m/aeolus_optimized.toml, warms the JIT cache at t_end=0.02 s, then profiles the main run (default t_end=2.0 s), dumps pstats to /tmp/qs_chem_profile_post.prof, prints top-cumulative + top-tottime tables. Use it to decide whether the next Phase B (analytic banded Jacobian for LSODA) is worth the bit-equivalence cost — currently held in BACKLOG behind the Phase A result.

The two companion diagnostic scripts that landed in the same wave (scripts/diagnose_qs_picard_2cycle.py + scripts/analyze_qs_picard_trace.py, MRs !130 / !131) are not Phase A artefacts — they support the orthogonal QS junction-Picard sign-bracket-bisection guard — but they are mentioned here so the chemistry-acceleration audit trail enumerates all shipped tooling from this wave.

Cross-reference. For the full qs_chem_fused_rhs semantics (when it activates, what falls through to the Python rhs, the non-Brinkmeier behaviour, opt-out flow) see the schema docstring at src/aeolus1d/io/schema.py lines 331–348 and the v0.9 entry in changelog.md — these are the authoritative run-settings references until the case-header reference chapter lands.

Phase A.5 — `dy` buffer hoist (2026-05-31)¶

MR: !135 (also rolled BASIC_FIELDS to PipeSpec / JunctionSpec / BCSpec and extended the form-factory kind dispatch; the dy hoist is the perf-relevant piece for this document).

Trigger. Re-profiling with scripts/profile_post_fused.py after Phase A landed showed a residual per-substep allocation: the fused _qs_brinkmeier_rhs_fused was building a fresh dy output array on every invocation. LSODA calls the rhs many times per macro step; the allocation was small but high-frequency.

Change. Hoist the dy buffer out of the inner rhs closure into a once-per-substep scratch slot, owned by the _qs_surface_chem_step caller and threaded into the fused rhs as an out-parameter. Algebraically identical — same write pattern, same final values, just zero allocations per rhs call.

Result. ~1.026× wall stacking on top of Phase A on the same AVL Cruise-M MFR_transition_SC630 workload. Small in isolation but free (no behavioural change, no flag, no fallback path) — the kind of micro-win that's only worth recording because it stacks cleanly with Phase A and confirms the diagnostic-script-driven loop is still finding signal.

Validation. Same bit-equivalence anchor as Phase A (tests/integration/test_qs_chem_fused_rhs.py); no separate test needed — the hoist is a calling-convention change inside the fused path that the integration test already covers.

Cumulative chemistry-heavy summary (post-Phase-A.5)¶

The chemistry-acceleration arc that started with the AVL drive cycle profile on 2026-05-29 now bakes in the following wall-clock budget by default, with the two True-default flips qs_chem_fused_rhs=True (MR !134) and the orthogonal substrate_wall_in_newton_residual=True (re-affirmed in MR !137 after the 5-SHA bit-identity bisect cleared it):

arc step	wall on 12 s AVL Cruise-M	realtime ratio	cumulative speedup
pre-Phase-A baseline	97.7 s	8.15×	1.00×
Phase A (fused @njit rhs)	47.0 s	3.92×	2.08×
Phase A.5 (`dy` buffer hoist)	45.8 s	3.82×	~2.13×

The ~2.1× lands on top of the earlier qs_cfl_adv=8 ~3.9× (the proven macro-step lever from the 2026-05-29 AVL profiling section), which is a user-tunable knob and stacks multiplicatively because they attack different axes: qs_cfl_adv reduces the number of macro chem calls per second of sim, the fused rhs reduces the per-call cost of each one. A run that combines both lands at ~8× wall on chemistry- heavy QS workloads versus the 2026-05-29 baseline.

The closed integrator-swap and ISAT verdicts (above) still hold: LSODA stays. The shipped path forward is the analytic banded Jacobian (Phase B, BACKLOG, gated on whether scripts/profile_post_fused.py shows that LSODA's Jacobian build is still a >10 % slice after Phase A.5) and the ML delta surrogate (ml_chemistry/, BACKLOG, gated on the run settling to a stable Δt). Both are research-grade efforts; Phase A + A.5 are the engineering wins that landed without that risk.

Full-case compute-budget gate (plan MR-7, 2026-06-07)¶

The Phase-Perf gates above watch kernels (wall-clock micro-benches, tests/benchmarks/test_baseline_regression.py) and Python wrappers (per-function cProfile self-time, test_profile_regression.py). MR-7 adds the missing whole-case tier: a budget on the total accounted compute for a canonical heavy case, so a catastrophic full-case regression — a flag default flip onto a slow driver, an accidental O(n^2), a lost numba JIT cache — trips a gate even when no single kernel/wrapper gate notices.

Metric — why sum-of-self-time, not wall. This box cannot measure wall tightly (kernel micro-benches swing ~2.3x with WSL2 host load; see the "Conclusions" caveats above and ADR-0006). But the cProfile sum of per-function self-time over a multi-second, million-call case is stable to ~1 %: three back-to-back profiled runs of the C02006 slice measured 44.07 / 44.25 / 44.51 s (spread 1.01x). Per-call jitter averages out across the whole profile. That stability is what lets the ceiling be tight enough to catch a 2x slowdown instead of only an order-of-magnitude one.

Case + slice. examples/avl_cruise_m/c02006_dpf_csf_regen_trans.toml — the independently-validated W3 DPF/CSF-regen anchor (also the golden-snapshot anchor in tests/regression/). The full case is a 500 s quasi-steady transient; its per-step cost is linear in t_end (measured ~3.1 sim-s/wall-s flat over 0-100 s, so the full run is ~2.6 min wall now — the older "6.5 min" figure was a coarser mesh). A bounded t_end = 100 slice (~31 s raw / ~50 s profiled) is therefore a faithful per-step-cost proxy, tractable both nightly and ad-hoc.

Reference + ceilings (median of 3 runs, WSL2 x86 runner, t_end = 100):

Target	Reference	Ceiling	Band
total accounted self-time	44.25 s	80.0 s	1.8x
`advance_network_quasi_steady_transient` (self)	7.14 s	13.0 s	1.8x
`_simple_solve_ideal_jit` (self)	3.68 s	7.0 s	1.9x

Every band is <2x, so a 2x slowdown trips the gate (the plan's acceptance) with margin, while the ~1.8x band over a ~1 %-stable metric is almost entirely headroom for CI host contention. The per-hot-path gates are diagnostics: a tripped total points at the culprit subsystem.

Where it runs. tests/benchmarks/test_fullcase_budget.py, marker profile_regress + slow. The new perf-budget CI job runs the whole profile_regress tier on scheduled (nightly) + tag pipelines only — a perf regression merged to main surfaces on the next nightly, not on the merge, keeping the single x86 runner free for the merge pipeline. This job also activates the previously-dormant case-41/50 wrapper gates in CI for the first time. Unconditional trip-teeth tests in the same file pin the 2x sensitivity so the gate can't silently bit-rot when the @slow tier is skipped.

Deferred to the bench host (Phase 5.5). Gating the full 500 s case and tightening the band toward a strict 2x both want a quiescent host; until then this loose form catches the regressions that matter (>=~1.8x) without flaking.

To refresh the reference (after an intentional perf change): re-measure on a quiescent box, update REF_* / *_BUDGET_S in the test in a dedicated commit, and update the table above.

Phase 5.5 — full-case wall-clock budget (the strict band)¶

MR-7 above is the form of the full-case gate that's safe to run anywhere: a cProfile sum-of-self-time proxy over a bounded 100 s slice, with ~1.8x bands. It deliberately does not assert wall-clock, on the standing assumption that this WSL2 host swings too much under load to gate wall (ADR-0006 / C2). Phase 5.5 closes the two things MR-7 deferred — the full 500 s case (not a 100 s slice) and a strict 2× wall-clock band (not the self-time proxy) — on the validated W3 anchor c02006_dpf_csf_regen_trans.

The host turned out to be this box¶

The 2.3x figure that justified the deferral is for kernel micro-benches. Full-case wall is a different animal, and it had never actually been measured. So it was (2026-06-07, .pipeline/charac_wall.py): two independent 5-rep series of the full C02006 500 s case on the x86 WSL2 runner.

series	reps (s)	best-of-N min	settled spread
1 (cold + a hammered box)	405, 247, 192.2, 195.1, 199.0	192.2	1.04x (reps 3-5)
2 (settled)	193.2, 194.7, 193.4, 193.5, 192.0	192.0	1.014x (all 5)

Two findings drop out:

Settled full-case wall is stable to ~1.4 % — a multi-second, million-call case averages out the per-call scheduling jitter that wrecks micro-benches.
The best-of-N minimum is immune to the cold-start transient. Series 1's first two reps were contended high (405 s is slower than the cold-compile warm run — external contention, not warm-up), yet its min still came out 192.2, agreeing with the quiescent series' 192.0 to 0.1 %. Contention is strictly one-sided, so the minimum is the honest estimator; the gate never looks at the noisy max.

A 2× band over a metric that reproduces to 0.1 % is almost all headroom: a genuine doubling (→ 384 s floor) trips it; transient host noise cannot. That retires the "needs a dedicated quiescent host" deferral for this anchor — the mac-mini is not required.

What shipped¶


Gate	`tests/benchmarks/test_fullcase_wallclock.py` (marker `wallclock_regress` + `slow`)
Anchor / ref	`c02006_dpf_csf_regen_trans`, full 500 s, `ref_wall_s = 192.0` (best-of-5)
Band	strict 2× (`ACCEPT_BAND`); trip is `measured < 2 × ref`
Reps	best-of-5 (`N_REPS`); rejects the cold-start transient
CI job	`wallclock-budget` — pinned to the x86 WSL2 runner (tag `aeolus-wsl-linux`), `AEOLUS_BENCH_HOST=1`, `scheduled` (nightly) + tag

The AEOLUS_BENCH_HOST=1 flag still guards every real gate: the runner is quiescent during its scheduled nightly job, but an interactive dev session on the same box is not (the series-1 contention came from exactly that). So an ad-hoc pytest -m wallclock_regress on a working box skips the real gate and runs only the unconditional trip-teeth that pin the 2× band logic. The pin to aeolus-wsl-linux matters because wall is host-specific — on the ARM mac-mini (tag mini) the 192 s reference would falsely trip.

2026-08-20: the ARM box is no longer the only host that would falsely trip it. A second x86_64 runner (aeolus-docker, tag aeolus-x86) took over the untagged tiers plus gui-tests and perf-budget; wallclock-budget stayed pinned here precisely because ref_wall_s = 192.0 describes this hardware and nothing else. Architecture was never the point — the measurement host was. Moving the pin requires a fresh best-of-5 on the new box first; see docs/ci-runner-runbook.md and the BACKLOG.md CI section.

Why EPW is not an anchor here¶

Plan 2.9 names "W3/EPW". The case its "EPW" was taken to mean, c06107_epw_car_diesel, is not gated here — and despite the name it is not a runnable breathing case. It is an import/schema/GUI smoke stub: a closed all-wall network with no cylinders, exercised only by a 0.1 ms run_case advance in test_avl_cruise_m_bundled_examples.py. It is never run to completion. Forcing its native t_end = 200 is ~15.6 h of meaningless idle acoustic relaxation (a 3 cm turbocharger pipe meshed to 20 cells pins a ~1 µs explicit dt), and it cannot run quasi-steady (all-wall — no inflow/outflow far BC for the QS driver). A prior ">14.5 min" figure for this case was wrong.

The actually-validated AVL EPW case is the examples/validation/avl_epw_diesel_steady/ harness — a 6-cyl 12.9 L HD diesel with real stagnation_inlet/pressure_outlet BCs, t_end ≈ 0.08 per operating point, driven per AVL CSV row. That is a validation-suite concern (see examples/validation/README.md), not a wall-clock budget anchor. So this gate covers W3 (c02006) only; a meaningful EPW perf anchor, if ever wanted, is the harness case, not c06107.

Refreshing the reference¶

After an intentional perf change, on the runner while it's idle:

AEOLUS_BENCH_HOST=1 python scripts/bench_fullcase.py            # best-of-5
AEOLUS_BENCH_HOST=1 python scripts/bench_fullcase.py --reps 7   # more margin

It times the anchor best-of-N and prints a paste-ready ANCHORS block; update ref_wall_s in test_fullcase_wallclock.py in a dedicated commit that records the host + arch, and update the table above.

The gates report on the passing path too (2026-08-21)¶

Until this change every gate helper in tests/benchmarks/ reported only by raising. A failing run printed its number; a passing run printed nothing.

The evidence is the v0.3.2 tag pipeline: perf-budget passed, and its entire 110-line log carried exactly one number — pytest's own 8 passed, 78 deselected in 62.25s. The gated quantity, total self-time against the 80.0 s ceiling, never appeared. The consequence is that neither perf gate could be trended across green releases, which is precisely when the trend is wanted; STATUS.md's host-load table can compare v0.3.1's 154.86 s loaded against 77.8 s idle only because that run failed. The same asymmetry is why the idle wallclock-budget figure in that table has to be caveated as a job duration rather than the measured per-rep wall.

What is emitted¶

tests/benchmarks/_timer.emit_measurement prints one AEOLUS-PERF-prefixed line and, when AEOLUS_PERF_LOG names a path, appends one JSON object per measurement to it. Three helpers call it on every exit path — passing, failing, and the early-error paths:

Helper	Gate	`gate` field
`_timer.assert_profile_self_time_under`	per-hot-path self-time (7 call sites)	`profile_self_time`
`test_fullcase_budget._assert_total_budget`	the 80.0 s total	`profile_total`
`test_fullcase_wallclock._assert_wall_budget`	the strict 2x wall band	`fullcase_wall`

Each record carries the gated quantity (value_s), what it was gated against (limit_s, ref_s), the verdict (passed, failure_mode) — and its harness: host, machine, python, cpu_count, numba, fastmath, the marker selection, the CI job/pipeline/runner tags, and loadavg sampled at the measurement rather than at session start. That last point is deliberate: host contention is the only thing that has ever made these gates flake — three recorded incidents — so the load that produced a number travels with it. CLAUDE.md is the binding rule here: a number without its case is not a measurement.

An unreferenced anchor (ref_wall_s is None, i.e. a bench host that has not been calibrated yet) emits its measured wall with passed: null and failure_mode: "no-reference", then skips as before. That is how a new bench host produces its first number.

Reading it off a pipeline¶

perf-budget and wallclock-budget both set AEOLUS_PERF_LOG and collect perf-measurements.jsonl with when: always — a red gate's number is the one most worth keeping. Download the artifact and:

jq -r '[.gate, .key // .slug, .value_s, .limit_s, .passed] | @tsv' perf-measurements.jsonl
jq -r 'select(.gate=="profile_total") | "\(.value_s) / \(.limit_s) on \(.host) load \(.loadavg[0])"' perf-measurements.jsonl

Before trusting any number from it, check loadavg in the record itself: a perf-budget total measured on a loaded box is not comparable with one from an idle box, and the margin here is thin — the last green run was 77.8 s against an 80.0 s ceiling.

Paste-on-tag. After a release tag's pipeline, take the profile_total and fullcase_wall records from the artifact and add a dated row to the tables above, with the record's host and loadavg beside the number. This is the step that was impossible before: the numbers existed only in the process that computed them.

tests/unit/test_perf_gates_emit_artifacts.py couples the YAML to the Python — the env var name is imported from _timer.py rather than copied, and the artifact path is compared against the path the job writes to. Both failure modes are otherwise silent: the job goes green and the artifact is empty or absent. tests/benchmarks/conftest.py fails any test carrying profile_regress or wallclock_regress that runs to completion without emitting, so a fourth gate helper added later cannot repeat the omission — the first draft of this change instrumented _timer.py alone and missed _assert_total_budget, the helper gating the very number STATUS.md names.

The first nightly, and what it says about the tag reading (2026-08-31)¶

The nightly schedule fired for the first time and both perf gates ran outside a tag pipeline for the first time ever. Until 2026-08-30 the project had no pipeline schedule at all, so perf-budget and wallclock-budget had only ever executed on a release tag — see docs/ci-runner-runbook.md for the schedule and BACKLOG.md for how that went unnoticed.

Pipeline 2803954760, ea06db47, source schedule, 72 min end-to-end, all 17 jobs green. Numbers read from each job's own perf-measurements.jsonl artifact, never from job duration — the two are different quantities and substituting one for the other is a mistake this repository has already recorded.

gate	value	ref	limit	loadavg at measurement
`profile_total` (`c02006_total_self_time`)	44.144 s	44.25	80.0	`[0.74, 0.26, 0.19]`
`fullcase_wall` (`c02006_dpf_csf_regen_trans`, 5 reps)	195.829 s	192.0	384.0	`[1.01, 1.08, 0.90]`

The `v0.3.3` tag reading was contention, and the loadavg was in the record all along¶

docs/release-readiness.md's EC-Q8 records the v0.3.3 tag pipeline at 61.70 s of the 80.0 s ceiling and concludes "18.3 s of margin where v0.3.1 had 2.2 s". Its own harness line carries loadavg 2.53.

Line the three quiet-host readings up against it:

when	loadavg	`profile_total`
by-hand pre-flight, 2026-08-24	0.78 → 2.12	43.28 s
nightly, 2026-08-31	0.74	44.144 s
reference	—	44.25 s
`v0.3.3` tag pipeline	2.53	61.70 s

On a quiet host the gate reads its reference to within 0.24 %. The tag pipeline's 61.70 s was 40 % above that, on a box carrying 2.53 of load. So the ceiling's real headroom is about 36 s, not 18.3 — and the figure that made it look like 18.3 was measuring the runner, not the code.

That is the entire argument for EC-Q10 in one measurement: the same number, taken nightly on a quiet box, costs nothing and is interpretable; taken once on release day under load, it is a scare.

This is n = 1 and must not be promoted¶

One nightly is a data point, not a baseline. The dead-JIT budget is the cautionary case in this same file's neighbourhood: a single 606 s sample read "67 % spent" and ten runs later the worst case restated it at 70 % with a bimodal distribution nobody had seen. Step 2 of EC-Q10 is to accumulate green nightlies until the baseline is a distribution; step 3 pins it. Do not refresh any reference number against this row.

One documentation defect found while reading it¶

docs/perf-history.md's own jq recipe above reads .loadavg[0], and that field is present — checked, because the first attempt at this table assumed it was missing when a truncated dump did not show it. Recorded because the near-miss is the point: the harness was in the record, and EC-Q8's conclusion was drawn without reading it.

Phase Perf-18 — fused uniform-γ boundary-face flux (2026-08-28)¶

Trigger. A profile of 39_4cyl_engine_turbo_scr taken while answering "would a TensorFlow-style computational graph make this faster". It said no — but it also said where the time actually is, and that is not the numerics. On that case only 20 % of wall sits in @njit kernels; the rest is Python orchestration. The share tracks component count, not cell count:

case	topology	wall	in kernels	Python
`01_single_pipe_sod`	1 pipe	0.03 s	91 %	9 %
`23_mp5_sod`	1 pipe, MP5	0.01 s	82 %	18 %
`40_implicit_adaptive_pipe`	Newton/BDF2	0.01 s	38 %	62 %
`74_twc`	quasi-steady + surface chem	0.20 s	33 %	67 %
`50_two_cyl_engine`	2-cyl network	0.37 s	21 %	79 %
`39_4cyl_engine_turbo_scr`	13 pipes, 10 junctions	6.89 s	20 %	80 %
`89_epw_diesel`	6-cyl, 17 pipes	0.27 s	14 %	86 %

Harness: cProfile self-time summed over every entry whose (file, line) matches a numba Dispatcher's py_func (217 found by walking solver, chemistry, closures, bc, components — the same inventory trick _dead_jit_probe uses), divided by an unprofiled best-of-3 wall clock. cProfile inflates this workload 1.84× here, matching the 1.79× recorded in BACKLOG.md, and the inflation lands on Python-heavy frames — so the kernel share is computed against true wall, not against profiled total. Cross-check: pipe_residual_muscl reads 110.6 µs/step in the profile and 112.8 µs/step in a standalone microbenchmark on the same real arrays, i.e. dispatcher rows are essentially uninflated. This dev box, 2026-08-28.

What landed. boundary_face_flux_muscl_jit in solver/kernels.py, and junction._pipe_port_flux_muscl reduced to one call into it. The body it replaces reached the same face flux through five dispatcher round-trips — U_to_W, muscl_slopes, W_to_U_single twice, hllc_flux — plus a slice and two 3-vectors, per port per RK stage. The O(1) 4-cell window moved across verbatim; Ng and N are still read from the pipe rather than assumed to be 2 as the two siblings assume, because the call site being replaced read them.

The variable-γ and real-gas boundary-face paths were already fused when they were added. The plain uniform-γ path — the most travelled of the three — was the one still assembled in Python.

Two increments landed under this phase.

junction._pipe_port_flux_muscl → one boundary_face_flux_muscl_jit call.
network_timeloop._left_boundary_mass_flux, scalar-γ MUSCL branch → the same kernel. That branch was a second hand-kept copy of the identical reconstruction, and it additionally ran a function-local import and two generator expressions per call, ~16 calls/step.

Increment 2's real value is not its share of the clock. The species-advection anchor and the junction port coupling are the two consumers of "the flux at face 0", and _species_rate telescopes every interior face off that anchor — so a discrepancy between them is a constant DC bias on every face, which fakes a uniform advective drift of any non-uniform composition, and is invisible for uniform composition and on first-order pipes. They were two copies; they are now one kernel and agree 13/13 exactly under the default build. That is checkable, unlike the residual claim retracted below.

Wall-clock. 39_4cyl_engine_turbo_scr, t_end = 0.01 (5252 steps), best-of-3 inner, numba caches deleted between arms, three replicates interleaved baseline → inc1 → both. Gated on /proc/loadavg 1-min < 1 before starting, per this repo's own rule about wall-clock on a loaded box:

replicate	baseline	inc 1 only	both	inc 1	both	inc 2 marginal
1	7.137 s	6.622 s	5.999 s	7.2 %	15.9 %	9.4 %
2	7.342 s	6.731 s	6.070 s	8.3 %	17.3 %	9.8 %
3	7.511 s	6.875 s	discarded	8.5 %	—	—

≈ 8 % for increment 1, ≈ 16-17 % for both; per step 1359-1430 µs → 1142-1156 µs.

Two things about that table are the point of recording it. The box drifted ~5 % slower across the run (baseline 7.137 → 7.342 → 7.511), which is why speedups are computed within a replicate and not against a pooled mean. And replicate 3's both arm read 13.728 s — 2.3× its siblings, with loadavg 12.70 at the end of the run — so it is discarded as contention, not averaged in. An earlier ungated attempt lost its first two points the same way, to an orphaned child of a killed harness.

An earlier draft of this entry recorded increment 1 alone at "≈ 9 %" from two replicates on an ungated box. The gated three-replicate figure is 7.2-8.5 %. Same direction, slightly smaller; the 9 % was the top of the range, not the middle. Benched rather than assumed because Phase Perf-13 is the counter-example: a fusion that looked obviously good regressed 9 %, and Perf-7's SIMD HLLC rewrite regressed 16–33×. The difference here is that Perf-13 fused numpy operations at a (1, 3) shape, where numpy's per-call overhead exceeds the scalar stores it replaced; this removes Python dispatcher round-trips, which is the Perf-2/4/6 shape.

Bit-identity — read this before fusing anything else. The fused kernel is not bit-identical to the Python version under the default build:

build	fused == old python
default (`fastmath` on)	13 / 26 faces (Δ ~2.8e-13 relative)
`AEOLUS_FASTMATH=0`	26 / 26

RETRACTED, same day, before this entry was ever acted on. A second column stood here claiming "fused == pipe_residual_muscl's face: 13/26 under fastmath, 26/26 under AEOLUS_FASTMATH=0". That was not measured against the residual. It was measured against a Python reimplementation of the residual's face algebra — the non-inlined flavour — which is the same flavour as the code being replaced, so the comparison could only ever restate the first column.

Re-measured against an njit reimplementation instead (right flavour: same decorator, same helpers, same order, so W_to_U_single inlines as it does in the kernel) the answer flipped — fused 13/13, python 10/13. Which would have made the fusion an improvement in flux consistency. That did not survive either: the njit proxy, asked to reproduce the kernel's own R[0] through -(F_1 - F_0)/dx, matches bit-exactly on only 7 of 13 pipes. LLVM optimises the kernel's full sweep differently from any two-face extraction of it.

So the honest position, and the one to carry forward: under the default fastmath build the kernel's internal face flux is not observable by reimplementation, and no claim about matching it can be made. R exposes only differences of adjacent faces, and every separately-compiled reconstruction of a single face disagrees with the kernel at ~1e-13. What is measured: under AEOLUS_FASTMATH=0 the fused kernel, the previous Python body, and both proxies are all bit-identical, which is the build where test_golden_bit_identity.py asserts anything.

Root cause, isolated rather than guessed: W_to_U_single carries inline="always". Called from Python it runs as a standalone compiled function; inlined into a compiled caller, LLVM reassociates its arithmetic with the surrounding expression under fastmath, moving the result ~2.8e-13 relative. U_to_W and muscl_slopes were ruled out first — window-in-njit, window-in-Python and full-array all produce bit-identical W and S.

This is survivable only because of an existing design decision: test_golden_bit_identity.py runs only under AEOLUS_FASTMATH=0, exactly where the fusion is exact. The default build is guarded by the measured tier at rtol=1e-8, five orders above the shift. Any further Option-B fusion will have the same property, so the strict tier must be run explicitly (AEOLUS_FASTMATH=0 pytest tests/regression) on every such change — the default pytest tests/regression run skips those 13 tests and tells you nothing.

Tests. tests/regression green both ways: 18 passed / 13 skipped by default, and 31 passed under AEOLUS_FASTMATH=0 with the assert_array_equal tier active. Numba caches deleted before the strict run — numba's cache key ignores inline="always" callees, so a stale caller would have served the old binary and made the check meaningless.

A third increment was measured and NOT taken — the useful negative. The plenum ghost fill (_plenum_ghost_primitives + _ghost_rho_E + the ghost writes) ranked next after the two above. Prototyped as one @njit kernel and microbenchmarked on 20 real plenum ports from 39_4cyl_engine_turbo_scr:

path	per call
Python, as it stands	4.197 µs
fused kernel	3.926 µs

Bit-identical, 20/20 — and worth 0.27 µs, i.e. 6.5 µs/step of a ~1150 µs step, +0.56 %. That is under the ~5 % drift this box shows within a single interleaved A/B run, so it is unmeasurable here even if real. Not wired.

The reason is worth more than the increment would have been. The fused path still costs 3.93 µs, and nearly all of it is gathering the arguments: volume.rho, volume.p, volume.gamma_advection (a property), pipe.U_pad[idx] (constructs a view), _pipe_cell_gamma, the gamma_pad slot lookup. Only the arithmetic moves into the kernel, and the arithmetic was never the cost.

The rule these three measurements give, stated so it can be applied without re-deriving it: fusing pays where it removes dispatcher round-trips, and does not pay where it merely moves scalar arithmetic behind one new dispatcher call. Increments 1 and 2 each replaced five dispatcher calls with one and returned ~8 %. This one would have replaced zero with one. Perf-13 is the same lesson from the numpy side — a slice-assign fused over three scalar stores at a (1, 3) shape regressed 9 %.

What this does not do. It removes two functions' worth of Python from the step. The orchestration share on this case is ~80 %; if it went to zero the step would fall from ~1312 µs to ~250 µs, a 5.2× ceiling. Reaching any real fraction of that means struct-of-arrays containers for pipes, junctions and cylinders so a whole step can be typed — and the negative above says incremental fusion has now hit its limit: what remains at these call sites is Python attribute access — object property reads, U_pad[idx] view construction, per-pipe getattr — which no kernel boundary can reach. Only typed containers can. That makes the SoA step the next real move rather than one more fused helper. Treat the 5.2× as a bound, not a target.

Performance history — the Phase Perf series¶

Canonical workloads¶

Cumulative trajectory¶

Phase Perf-1 — EOS / solver hot-path optimization sweep¶

Phase Perf-2-quick — D_J broadcast¶

Phase Perf-2 — JIT-ify closure-Jacobian dispatcher and Newton line-search¶

Phase Perf-3 — chord (lazy-Jacobian) reuse¶

Phase 27c — opt-in parallel per-pipe Newton solves¶

Phase Perf-4 — JIT'd BC dispatchers + chord knobs on TOML¶

Phase 27c-cells — prange Jacobian (prototyped, reverted)¶

Phase Perf-5 — odds and ends¶

Phase Perf-6 — chemistry advection JIT¶

Phase Perf-6r2 — chemistry / filter Python wrapper hoist¶

Phase Perf-6r3 — junction-code Python wrapper hoist¶

Phase Perf-6r4 — junction follow-up (remaining wrappers)¶

Phase Perf-6r5 — cylinder-side wrappers¶

Phase Perf-6r2-corrections-2 — profile-based regression harness¶

Phase Perf-6r6 — turbo-side wrapper hoist¶

Phase Perf-8 — getattr-dispatch removal (species_coupling + muscl boundary)¶

Phase Perf-9 — investigate remaining getattr calls (no-op)¶

Phase Perf-history-3 — re-profile post-Perf-8 to find next-tier targets¶

Phase Perf-10 — Pipe.has_species → cached field¶

Phase Perf-11 — Pipe.U property → cached slice view¶

Phase Perf-12 — pre-allocated scratch buffers in regen_into¶

13 (_fill_plenum_ghost per-cell loop fusion) and #14¶

Phase Perf-13 — _fill_plenum_ghost loop fusion (REVERTED)¶

Phase Perf-history-4 — re-profile post-Perf-10/11/12¶

Phase Perf-15 — Volume.p cached field (deliberate no-op)¶

Phase Perf-17 — bc.base.wall_bc per-port hoist¶

Phase Perf-7 — packed-kernel parallel-vs-serial dispatcher¶

Phase I-8 (MVP) — opt-in Jacobian-Free Newton-Krylov¶

Long-pipe bench (Phase Perf-history follow-up)¶

Phase 27c soak test¶

Canonical-green baseline + xdist contention¶

Things that were tried and rejected (or never attempted)¶

Tools that made this possible¶

What the audit didn't move¶

v0.3 aftertreatment time-stepping spike (2026-05-20)¶

Step economics (case 41, t_end = 0.06 s)¶

Self-time by solver phase¶

500 s drive-cycle extrapolation (current per-step cost)¶

Conclusions (these set the v0.3 ordering)¶

Quasi-steady SIMPLE kernel — numba JIT (2026-05-21)¶

AVL C02006 regeneration — end-to-end (accelerated-ramp validation)¶

Phase Perf-7r — turbo map bilinear lookup (2026-05-23)¶

AVL drive cycle (MFR_transition_SC630) profiling — 2026-05-29¶

Scoped: ISAT chemistry acceleration — 2026-05-29¶

ISAT P1+P2 results — 2026-05-29¶

ISAT vs the PRODUCTION baseline — NO-GO — 2026-05-29¶

Phase A — fused @njit Brinkmeier rhs (2026-05-31)¶

Phase A.5 — dy buffer hoist (2026-05-31)¶

Cumulative chemistry-heavy summary (post-Phase-A.5)¶

Full-case compute-budget gate (plan MR-7, 2026-06-07)¶

Phase 5.5 — full-case wall-clock budget (the strict band)¶

The host turned out to be this box¶

What shipped¶

Why EPW is not an anchor here¶

Refreshing the reference¶

The gates report on the passing path too (2026-08-21)¶

What is emitted¶

Reading it off a pipeline¶

The first nightly, and what it says about the tag reading (2026-08-31)¶

The v0.3.3 tag reading was contention, and the loadavg was in the record all along¶

This is n = 1 and must not be promoted¶

One documentation defect found while reading it¶

Phase Perf-18 — fused uniform-γ boundary-face flux (2026-08-28)¶

Phase 27c-cells — `prange` Jacobian (prototyped, reverted)¶

13 (`_fill_plenum_ghost` per-cell loop fusion) and #14¶

Phase A.5 — `dy` buffer hoist (2026-05-31)¶

The `v0.3.3` tag reading was contention, and the loadavg was in the record all along¶