meta(amr): forest performance — the coarse–fine exchange, not regridding, is the AMR tax

Start with examples/niederer_benchmark.jl and run the documented full or NIEDERER_SMOKE benchmark to reproduce the phase table and inspect examples/niederer_cost.png. Read the linked candidate issues #83, #42, #84, #85, #78, and #72 before choosing a focused task; this parent tracker is complete when the selected optimization is measured without changing the nine-point correctness results.

Status (2026-09-30): #42, #84 and #85 have landed; the phase tables below predate them and need a re-run.

The new Niederer 2011 benchmark example (examples/niederer_benchmark.jl) is the first 3D BlockForest in this repo and the first workload that times the adaptive path against a like-for-like uniform baseline. It came out 2.34× faster overall — but the phase breakdown says the overhead is not where I expected, so this issue records the numbers and the optimization candidates they point at.

What the benchmark is

Monodomain cardiac EP on a 20 × 7 × 3 mm slab, ten Tusscher–Panfilov 2006 epicardial kinetics (18 gating/concentration states per cell, carried through regrid! as one SVector{18}-valued BlockField), anisotropic diffusion D_T ∇² + (D_L − D_T) ∂ₓₓ, Godunov splitting with explicit Euler on both halves, dt = 0.01 ms to t = 60 ms. Base 48×16×8, blocksize = (8,8,8), maxlevel = 2, so the finest spacing is ≈0.1 mm — Niederer's reference resolution — and the equal-resolution uniform grid is 192 × 64 × 32 = 393,216 cells.

Correctness first, because the perf numbers only matter if the answer is right. The scored metric is activation time at nine fixed points. Against an independent finite-difference code at dx = 0.1 mm the mean |Δ| is 0.40 ms (max 0.80). Against a uniform forest at the same finest spacing, all nine points agree to 0.00 ms. Every 3D coarse–fine path — Martin–Cartwright interpolation, flux-matching restriction, the balance cascade, the regrid! transfer of both a scalar and an 18-component field — was previously unexercised, and it is exact here.

The cost breakdown

5 threads, macOS/aarch64. The example now times both runs and splits the wall clock by phase.

phase	adaptive	uniform	speedup
diffusion	19.1 s	8.5 s	0.44×
reaction	12.1 s	68.9 s	5.70×
regrid	1.9 s	—	—
total	33.2 s	77.5 s	2.34×

The forest carries 18.8% of the uniform cell count on average, so the identity is 5.3× fewer cell-updates ÷ 2.28× dearer per update = 2.34× net (75 ns vs 33 ns per cell-update). Slightly more than half the cell saving is handed back.

Two things fall out:

• Reaction scales as hoped. It is pointwise, tracks the cell count, and at 5.70× marginally beats the 5.3× cell ratio because the smaller working set is kinder to cache.
• Regridding is cheap. 120 regrids — marking, topology edit, transfer of two fields, schedule rebuild — cost 1.9 s, under 6% of the run. This is the thing I assumed would dominate, and it does not.

Diffusion is the problem: 2.2× slower on 5.3× fewer cells, a ~12× per-cell penalty.

Where the diffusion time goes

Timing the pieces directly on a representative adapted forest (166 leaves, levels 0–2, refined in a slab-crossing band the way the wavefront refines it) against the uniform one:

	adapted, 166 leaves	uniform, 768 leaves
halo_update!	0.94 ms	0.22 ms
apply_bc!	~0 ms	0.02 ms
one laplacian apply	1.01 ms	0.53 ms
the anisotropic apply	2.06 ms	1.07 ms
exchange share of an apply	92%	45%
per cell-update	24.2 ns	2.7 ns

On an adapted forest the halo exchange is the operator — 92% of an apply is ghost-filling, because every coarse–fine face runs quadratic interpolation and flux-matching restriction instead of a slab copy. Per leaf that is ~20× the cost of a same-level exchange.

Optimization candidates

This issue is the parent tracker for AMR forest performance. The numbers above that ranked these candidates are the pre-#84/#85 phase tables — the 92% exchange share and the 0.44× diffusion figure both predate those fusions — so re-run examples/niederer_benchmark.jl and post the updated phase table before prioritizing the open items.

• #42 — make GhostFill isbits with a static term count (the prerequisite for #84).
• #84 — _run_fills! made one broadcast pass per term: 19 terms per interpolation fill in 3D, 16 fills per coarse–fine face, ~300 broadcasts per face. Fused into one gather per fill.
• #85 — Added exchanged once per operand on a forest: two identical exchanges per step for a two-leaf Added, which with the exchange at 92% of an apply roughly doubled diffusion cost (~9 s of the 33 s run, estimated).
• #94 — the adjoint half of #84: halo_update_adjoint! still runs one scatter-add broadcast per CSR term, and fusing it bit-identically needs a source-centric transposed CSR built at schedule time.
• #78 — batch the coarse–fine flux rewrite into one device launch for packed Diffusion.
• #72 — regrid cadence, buffer zone, and prepare reuse across timesteps (where the per-step BlockField allocation in examples/monodomain_amr.jl lives).
• Measure packed-on-CPU: PackedBlockField for a regridding workload. The packed path's reported 30–95× wins are GPU numbers. This workload regrids every 50 steps, and the documented contract is that a packed field cannot be regridded — you regrid the reference field and re-pack. Nobody has measured whether 50 packed steps amortize a re-pack on CPU. Given #83 threads the same reference path, worth deciding before optimizing it twice.

Cross-reference — threading, tracked under #51: #83. The reference forest path is single-threaded end to end: grep -rn "Threads\.\|@spawn\|@threads" src/ returns nothing; the leaf sweeps and the exchange descriptor loops are serial for loops. Likely the single biggest lever; it shares the executor decision with #51, which is why it is parented there rather than here.

Reproducing

julia -t auto --project=examples examples/niederer_benchmark.jl     # both runs, ~2 min
NIEDERER_SMOKE=true ... # seconds, reduced levels/duration
NIEDERER_SKIP_UNIFORM=true ... # adaptive only

The example prints the phase table and writes examples/niederer_cost.png (cells over time, and wall time by phase). Full write-up in docs/pages/benchmark.qmd.

One caveat on all of the above: these are single-machine, 5-thread, CPU, BlockField numbers at one problem size. The relative weight of exchange versus stencil will move with block size, level count and thread count, and #83 in particular could reorder the rest.

—
🤖 Claude ticked the boxes here. The stopwatch still points at the exchange.
_{claude-opus-5-5[1m] (high) · claude-code 2.1.285 · unreviewed}