ADR-0015 — Struct-of-arrays network state, so a whole step can be compiled¶

Status: Proposed — phases 1-4 and 6a implemented and measured on feat/adr-0015-soa-network-state; see "What the implementation measured" below before accepting.

Context¶

Profiling 39_4cyl_engine_turbo_scr on 2026-08-28 established that the solver's cost is not its numerics. Only 20 % of wall time is inside @njit kernels; the rest is Python orchestration. The share tracks component count, not mesh size:

case	topology	wall	in kernels	Python
`01_single_pipe_sod`	1 pipe	0.03 s	91 %	9 %
`74_twc`	1 pipe, 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 %

A single pipe is already a compiled program. A 17-pipe engine is Python with numba subroutines. Phase Perf-18 then took two increments of the obvious remedy — fuse a Python helper into one kernel call — and got ~8 % each. A third was measured and rejected at +0.56 %, which is where the useful part starts.

Why more fusing will not work¶

The rejected increment (plenum ghost fill) is bit-identical and still costs 3.93 µs against the Python path's 4.20 µs, because the fused version still has to gather its arguments in Python. Measured cost of a single access on this box (200 K reps, best-of-7, net of an empty-lambda baseline of 34.2 ns):

access	net cost
`pipe.Ng` (plain field)	15 ns
`pipe.U_pad` (array field)	10 ns
`pipe.U_pad[idx]` (row view)	78 ns
`pipe.U_pad[idx, 0]` (scalar)	75 ns
`volume.rho` (property)	79 ns
`volume.p` (computed property)	363 ns
`volume.gamma_advection` (property)	55 ns
`_pipe_cell_gamma(pipe, idx)` (helper)	101 ns

The rule those measurements give, and the reason this ADR exists: fusing pays where it removes dispatcher round-trips and does not pay where it merely moves scalar arithmetic behind one new dispatcher call. What is left at these call sites is Python attribute access, and no kernel boundary reaches it.

The trap this ADR must not fall into¶

Struct-of-arrays is not itself the speedup. Measured on the same box:

arr[idx, 0]        (plain local ndarray)  70 ns
pipe.U_pad[idx, 0] (through the object)   75 ns

Flattening state into arrays and then indexing it from Python buys 5 ns. Numpy scalar indexing costs ~70 ns whoever owns the array. SoA is worth doing only as the enabler for compiling the loop that does the indexing — it makes the state typeable so the sweep can move inside @njit, where an index is an instruction rather than an object protocol. Any SoA change that does not end with a compiled sweep has bought nothing, and this is the failure mode to review for.

Proposal¶

1. Layout — flat with offsets, not padded¶

Pipes are ragged: case 39 runs 30–80 cells, case 89 runs 8–20. Padding to Nmax costs 2.05× the memory on case 39 (25.6 kB vs 12.5 kB), measured. It does not cost compute in a loop-over-pipes kernel, which skips the pad via N[i] — the two layouts benchmark identically (below). It would cost compute in a fully vectorised tensor formulation, which evaluates every lane. Flat is chosen for the memory and because no lane ever holds a value nobody wrote.

U_flat      (sum(N_i + 2*Ng_i), 3)  float64   — conserved state, all pipes
p_off       (P+1,)                  int64     — CSR-style start offsets
p_N, p_Ng   (P,)                    int64
p_dx, p_gamma (P,)                  float64
p_scheme, p_limiter, p_flags (P,)   int64     — flags is a bitfield
p_friction_kind (P,), p_friction_params (P, k)
... one parallel array per scalar field the sweep reads

Volumes and ports are the connectivity, as an edge list — which is the SoA form of a graph:

v_m, v_E, v_V, v_gamma, v_p   (V,)   — v_p refreshed once per step
port_pipe, port_vol (Pt,) int64
port_end            (Pt,) int64      — 0 = LEFT, 1 = RIGHT
port_A, port_Cd     (Pt,) float64

2. Migration — views, not copies¶

The step that makes this incremental and verifiable: Pipe.U_pad becomes a view into U_flat, not a copy of it. Both representations are then the same memory, coherent by construction, with no synchronisation to get wrong. Every existing Python path keeps working unchanged while sweeps move across one at a time, each A/B'd bit-for-bit against the path it replaces.

This constrains one thing: U_pad may no longer be reassigned after build (resize, set_species growth). Those sites must rebuild the flat buffer and re-view, or reserve up front.

3. Fast path and fallback¶

Not every pipe can be typed: Python closures (py_reaction, py_friction, the UDF seams), real-gas EOS, catalyst and filter objects are Python by construction. Measured eligibility — pipes carrying none of those:

case	pipes	typed-eligible
`50_two_cyl_engine`	6	6 (100 %)
`89_epw_diesel`	17	16 (94 %)
`39_4cyl_engine_turbo_scr`	13	12 (92 %)
`41_dpf_regen`	1	0 (0 %)

The split lands well: the engine networks, which are 80–86 % Python and have the most to gain, are 92–100 % eligible. The aftertreatment cases are not eligible at all — and are also not what this fixes, since their overhead is per-pipe chemistry wrappers rather than network orchestration. A pipe that is ineligible keeps today's Python sweep; the two must produce identical results, which is the acceptance test.

4. What comes free¶

Volume.p is a computed property costing 363 ns, read ~64 times per step (~2 % of wall). Phase Perf-15 declined to cache it: E_internal and m are mutated at 8+ sites, so a cache needs invalidation everywhere and "cache desync produces silently wrong pressures". Perf-15 names the correct fix and defers it — "a hook the timeloop calls once per step … scope-creeps into a Volume-state-management rewrite". That rewrite is this ADR. Under SoA v_p is a column refreshed once per step at a single sanctioned point, and the invalidation problem does not exist to be got wrong.

What has been tested, and what has not¶

This ADR is Proposed, so the distinction matters more than the proposal.

Tested — the layout claims hold.

The view migration works. Every pipe's U_pad rebound to a slice of one flat buffer (_U_view refreshed with it), then 39_4cyl_engine_turbo_scr run to t_end = 0.01: 13/13 pipes bit-identical, worst relative difference exactly 0.0, t_final equal. All 13 views report C_CONTIGUOUS, so numba compiles the same specialisation and no kernel sees a different array type. There is exactly one U_pad assignment site in the tree (pipe.py:588, __post_init__), so nothing reassigns it out from under a view.

The flat+offsets layout costs nothing over padded. One compiled sweep over the real 13-pipe geometry, best-of-9:

	µs	vs today
per-pipe loop (today)	58.21	1.00×
padded batched, serial	16.05	3.63×
padded batched, prange	8.45	6.89×
flat+offsets, serial	16.23	3.59×
flat+offsets, prange	8.50	6.85×

Flat matches padded to within noise and halves the memory. Correct to 7.6e-13 against the reference kernel (fastmath reassociation, the same order as every other fusion in Phase Perf-18).

The acceptance test this ADR set itself — a compiled sweep over the port edge list — has now been built and passes.

Prototype: port_ghost_sweep(U_flat, p_off, p_Ng, p_N, p_gamma, v_rho, v_p, v_gamma, pt_pipe, pt_vol, pt_right) — one call filling the ghost cells of every plenum port in the network, same algebra as _plenum_ghost_primitives + _ghost_rho_E + the writes. Case 39: 13 pipes, 10 volumes, 20 plenum ports. Bit-identical to the Python path, worst relative difference exactly 0.0.

	total	per port
Python port loop (today)	52.43 µs	2621 ns
compiled sweep, volume state still gathered in Python	8.88 µs	444 ns
— of which the Python gather	7.32 µs (82 %)
— kernel alone	1.10 µs	55 ns

5.9× as prototyped; 47.7× with the gather removed, which a full SoA does by construction. End-to-end on this subsystem: 7.6 % now, 8.9 % ceiling.

The comparison that matters is against Phase Perf-18's rejected third increment. That fused the same arithmetic behind a per-port dispatcher call and returned +0.56 %. Moving the loop inside the kernel returns +7.6 % from identical algebra — a 13× difference produced entirely by where the loop lives. This is the thesis of the ADR, measured rather than argued.

The 82 % figure also independently confirms the design's v_p decision: the gather is 30 property reads (10 volumes × rho, p, gamma_advection) at ~244 ns each, dominated by the 363 ns Volume.p that Phase Perf-15 declined to cache. Under SoA it is a column and the cost is not paid at all.

Still not tested. The prototype covers the scalar-γ ideal plenum port variant (20 of 20 ports in this case). The py_junction_bc, lodi and neck variants are excluded and need either their own kernels or a Python fallback. And apply_port_bcs does more than fill ghosts: the fuel and non-fuel species writes that follow are dict-keyed (m_species.get(name), iter_pipe_species) and do not type. They were ~28 % of that function's cumulative cost in the profile and stay in Python under this design, so the port loop cannot go to zero — only its ghost-fill half can.

Per-subsystem gains do not compound into 5.2× on their own. Each sweep converted returns 7-9 % because each subsystem is 8-10 % of the step. Reaching the ceiling means converting all of them; stopping after two or three leaves most of the 80 % in place.

Species addressing — global IDs, sparse storage¶

The species path is the part of the port loop that does not type today: m_species.get(name), iter_pipe_species, species_inflow_*[key] are all name-keyed dicts. The obvious remedy is "give every species a fixed global index". Measured, at a realistic 5 columns per port:

approach	ns/port
dict by name (today), in Python	310
integer scalar index, in Python	975 — 3× slower
integer fancy index, in Python	731 — 2.4× slower
`@njit` kernel called once per port	256
`@njit` sweep over 20 ports	21 — 15×

Converting to integer indices and staying in Python is a regression, not an optimisation. Python dict lookup is ~50 ns; numpy scalar indexing is 70-200 ns and fancy indexing carries ~1 µs of fixed overhead. The dicts are not slow — they are simply unreachable from a kernel. So global species IDs are worth adopting for exactly one reason: they make the species path compilable, and the win arrives only with the sweep. Half-done, it is 3× worse than today. This is the same shape as the port-sweep result above and as Phase Perf-18's rejected increment; it is now confirmed for the species path specifically.

Storage stays sparse. A census of the bundled corpus: 28 pipes carry extras, n_extra ranges 1-12 with a mean of 3.2, against a canonical SPECIES_NAMES of 12 — 30 % mean occupancy. Storing Y_extras dense by global ID would inflate the per-cell arrays ~3.3× and advect zeros. The proposal is therefore:

a global species ID registry, resolved once at build time, extensible per case — it cannot be the hardcoded canonical 12, because cases carry species outside it (41_dpf_regen advects PM1);
per-pipe sparse storage unchanged — Y_extras stays (N, n_extra);
a per-pipe int64[] col_sid mapping column → global ID, so kernels address species by integer and never by name;
volume composition dense by global ID (v_Y[V, S]), which is cheap because S is small and V is not per-cell.

This generalises ADR-0003 rather than replacing it. That ADR's invariant — "at the kernel boundary, species are addressed by canonical name→slot, never by raw column index" — stays true; what changes is that the name→slot resolution moves from per-call pad/unpack to build time, and the slot space becomes the case's registry rather than one chemistry kernel's fixed 12.

What the implementation measured (2026-08-29)¶

Phases 1-4 and 6a are built. Phase 5 (global species IDs) and the rest of phase 6 (accumulate_rates, the BC fill) are not.

case	adopts	min Δ	median Δ
`89_epw_diesel`	15 ports, 15 pipes	+10.8 %	+14.1 %
`50_two_cyl_engine`	10 ports	+10.7 %	+9.2 %
`39_4cyl_engine_turbo_scr`	20 ports, 12 pipes	+9.7 %	+7.8 %
`38`, `01`, `41`, `74`	nothing	within noise	sign disagrees

Seven replicates, identical step counts, CPU time rather than wall — another session was saturating the box for part of this work and the same comparison swung -3.4 % to +15.9 % on wall clock. The four non-adopting cases disagree in sign between minima and medians, which is what "adopts nothing" should look like.

Bit-identical to the Python path under AEOLUS_FASTMATH=0 (13/13 and 17/17 pipes on the two engine cases); ~2.6e-15 under the default fastmath build. tests/unit and tests/integration exit 0, tests/regression exit 0 both ways with the assert_array_equal tier active.

Four things the ADR did not predict¶

The argument count is a first-order cost. Numba charges roughly 160 ns per array argument to unbox and type-check. The first draft passed 13 arrays, so a zero-port call cost 2066 ns — on a 2-port network that was the entire result, and the sweep measured a 5 % regression. Packing the columns into three 2-D blocks took it to four arrays and 414 ns. Any future sweep should be designed around this from the start.

Both sweeps need a minimum size, and the crossovers are unresolved. The residual sweep copies into a shared buffer where the per-pipe path returned the kernel's own array: -1.9 % at 1 pipe and -7.6 % at 3, against +9.4 % at 12. The ghost sweep is -1.0 % at 2 ports even though the arithmetic says it should already pay. Thresholds are 8 pipes and 4 ports; the corpus jumps 3 -> 12 and 2 -> 10 with nothing between, so neither number is a measured boundary.

A convenience property on a per-step path cost more than a rejected optimisation. n_ports as a @property, read twice per step, was a consistent ~1.5 % loss on a case that adopts nothing — larger than the +0.56 % increment Phase Perf-18 rejected as too small to bother with.

A persistent skip flag is the wrong shape. Marking ports as "the sweep owns this" and honouring it unconditionally broke every driver that fills ghosts without sweeping — the implicit driver, analysis.linearize, two bench harnesses. Fifteen tests caught it. The skip must be scoped to the pass that actually swept, so correctness does not depend on having enumerated the callers right.

What is left, and what it is worth¶

Each converted sweep returns roughly what its subsystem is worth, and the two built here return ~10 % together on engine networks. The 5.2 × bound needs the species path, accumulate_rates, the BC fill and the cylinder ports as well. This banks about a tenth of it. Whether the rest is worth the same again in effort is the decision this ADR still needs.

Consequences¶

Expected gain. If orchestration went to zero the step would fall from ~1150 µs to ~250 µs, a 5.2× ceiling. It will not go to zero. Treat 5.2× as a bound, not a target, and require a load-gated interleaved A/B per increment — this box drifts ~5 % within a single run and has produced 2.3× contaminated outliers.

Bit-identity. Every fusion so far shifts the default fastmath build by ~1e-13, because inline="always" helpers reassociate once their caller is compiled. This is survivable only because test_golden_bit_identity.py runs only under AEOLUS_FASTMATH=0, where the fusions are exact. Changing the memory layout can additionally change vectorisation, so the strict tier must be run explicitly (AEOLUS_FASTMATH=0 pytest tests/regression) on every increment — the default run skips those 13 tests — and numba's *.nbc / *.nbi must be deleted first, since its cache key ignores inline="always" callees.

Blast radius. Pipe, Volume, Junction, Port, build.py, network_timeloop.py, and every consumer that reads a pipe attribute. This is the largest structural change proposed for the solver.

Alternatives considered¶

More helper fusion. Rejected on measurement: the third increment returned +0.56 %, below what this box resolves. The remaining cost is attribute access.
numba.experimental.jitclass / typed.List. Keeps the object shape but carries its own boxing and reflection costs, and would freeze the class layout against a codebase that still adds fields. Plain parallel arrays impose no such constraint.
A TensorFlow/XLA-style graph. Measured and rejected separately: a 480-cell network is far below the size where a GPU's ~150 µs launch floor pays, time-stepping is strictly sequential, and torch.compile on CPU lost to numba on every shape but one. GPU remains the right answer for ensembles (DoE/UQ, ≥512 independent cases, 8–15× at fp32) — a separate, additive axis.
Do nothing. Defensible. The solver is fast enough for present cases, and this is a large change with a bit-level risk surface. The counter-argument is that component count is the growth direction, and the Python share grows with it: 79 % at 6 pipes, 86 % at 17.