feat(grids): active-cell mask so stencil operators can skip inactive cells

[kylebeggs]

What's missing

MFO has no domain masking. Every mask in src/ is bounds-masking inside adjoint gathers (_maskedget, src/operators/abstract.jl:210) — there is no active-cell indicator, no irregular-geometry / cut-cell support, and no Field of Bool that participates in stencils. BlockForest refinement changes resolution, not participation, so AMR is not a substitute.

Why (downstream)

The microvascular transport model solves a reaction–diffusion system on a hemisphere of tissue embedded in a Cartesian box. Cells outside the hemisphere must be excluded from both the reaction and the diffusion stencil. Today the downstream package carries a per-cell active::Vector that it multiplies into every reaction term by hand, because a masked operator does not exist.

Getting this wrong is not cosmetic: leaving the box corners active created a spurious O₂ sink outside the domain that broke the radial-O₂ comparator entirely. The mask is load-bearing physics, and it currently lives downstream where the operator algebra can't see it.

Sketch

g = CartesianGrid(...; bc = ((Neumann(), Neumann()), ...))
active = scalar_field(Bool, g); set!(active, inside_dome)
L = laplacian(g; active)          # inactive cells contribute no flux and get no update

The interesting constraints:

• Adjoint correctness. Whatever apply! does, apply_adjoint! has to remain its exact discrete transpose — masking a face flux must mask the same face on the way back, or ⟨Lx,y⟩ = ⟨x,Lᵀy⟩ breaks. That's the test that should gate this.
• operator_diagonal. An inactive neighbour looks like a no-flux face, so the diagonal correction is the same (-2 + bc_sign)h⁻² form already handled for BC faces — but it now varies cell to cell, which the current implementation doesn't allow for.
• Singularity. A fully masked-off region turns the operator singular, which matters for the multigrid coarsest-level LU (_build_coarsest, src/multigrid.jl:288) and for coarsening a mask consistently.

Scope note

Uniform CartesianGrid only would already unblock the downstream use. Masked BlockForest / distributed slabs can stay out of scope — same boundary you drew for Diffusion in #54.

---

_{🤖 Filed by Claude, who read the code so Kyle didn't have to scroll.}

[kylebeggs]

Evidence from two Kaboudian/Cherry/Fenton papers

Kyle read Kaboudian, Cherry & Fenton, Sci. Adv. 5:eaav6019 (2019) and the same authors' Computing in Cardiology 50 (2023) paper, "GPU Load Balancing Using Sparse Cartesian Grids". The 2023 one is this issue's problem, solved, with numbers. Five things it changes here.

1. The fill fraction is much worse than "mostly empty"

Across the Strocchi 24-patient four-chamber cohort, tissue occupies 3.42 ± 0.72 % of the ventricular bounding box and 10.39 ± 2.81 % of the atrial one. Their measured wins from compacting it out:

structure	fill	speedup	memory vs. dense
atria (minimal model)	3.24 %	>20×	<4 %
ventricles (41-state OVVR)	11.7 %	>7×	<13 %

Speedup tracks 1/fill, which is the tell: the entire win is never launching work for a dead cell. A predicated mask — the natural reading of the sketch in this issue — gets the correctness fix the microvascular model needs and approximately none of that, because on a SIMT device the warp launches regardless. Worth separating the two goals explicitly in the issue rather than assuming one implies the other.

2. Their compaction is branchless, which is why it survives our invariants

Two integer index arrays: full → compact and compact → full. Every dead cell in the forward map points at one shared sentinel slot, whose reverse entry is the sentinel (-1,-1). A stencil arm that leaves the tissue lands in a designated garbage cell instead of taking a branch.

That detail matters for us specifically: a gather through an index array is still a broadcast (u[nbr]), so it does not break the "leaf bodies stay array-level" invariant — which is exactly the objection #69 raises against per-cell variable stencils.

The adjoint is not free, though. Forward is a gather (many readers, one garbage value — fine). The transpose is a scatter, and every dead neighbour scatters into that same sentinel slot, so a naive KernelAbstractions transpose kernel races on it. The fix is cheap — the sentinel's cotangent is discarded by construction, so mask it or give it a per-thread dump — but it has to be designed in rather than discovered. This issue already names ⟨Lx,y⟩ = ⟨x,Lᵀy⟩ as the gate; this is the concrete thing that will fail it.

3. The per-cell diagonal has a closed form

Both papers impose ∂u/∂n = 0 on the tissue surface by mirroring: if a stencil arm points at a non-tissue cell, substitute the centre value. That is precisely _bc_sign(::Neumann) = +1, the symmetric mirror about the face (src/boundaries.jl:69-71) — for the first ghost layer, reflecting about the face is "take the centre value." Only the locus changes: six box faces → per cell, per face.

So the open question in the body above — "the diagonal now varies cell to cell, which the current implementation doesn't allow for" — has a one-line answer under mirroring:

diag(L)ᵢ = −(number of live neighbours of i) / h²

directly computable from the mask, and exactly what a Jacobi/Chebyshev multigrid smoother needs. It also means _build_coarsest only goes singular for a fully enclosed all-Neumann region, which is the same pure-Neumann singularity we already have, not a new failure mode.

Naming trap: _maskedget (src/operators/abstract.jl:212) already exists and already says "mask," but it is a bounds guard whose out-of-range neighbours contribute zero. A geometry mask at a mirrored face must contribute the centre value. Same word, opposite rule — do not reuse it.

4. The Diffusion leaf is already index-driven, so this is a smaller change than it looks

_apply_raw! (src/operators/diffusion.jl:532-545) is a broadcast over an index set:

yi .= α .* _diff_at.(Ref(x.data), Ref(κ), interior(g), …)

Swapping interior(g) for a vector of live CartesianIndexes — or, for compaction, a full_index array — is a drop-in substitution that leaves the stencil body, the Val{D} recursion in _diff_axes, and the whole AD path untouched. Masking and compaction can share one seam, and the array-level-leaf invariant survives both. This is the strongest architectural argument for doing this at all, and it wasn't in the issue.

Caveat on the GPU half: Diffusion has no KernelAbstractions kernel at all today. It runs through Base broadcast on CartesianGrid and falls back to the per-leaf reference sweep on a packed forest (src/operators/forest_packed.jl covers Laplacian/Derivative/Gradient/Divergence/Advection/Scaling — not Diffusion). The GPU load-balancing win above has no MFO kernel to land in yet; #59 is the prerequisite for any of section 1 and 2 mattering.

5. Pushback: their 20× is against a dense baseline, and we have a forest

For a heart in a box, most 8³ root blocks are entirely dead. Never creating those blocks captures the bulk of the same win with no index indirection and no adjoint scatter hazard — but BlockForest's constructor requires an exact root tiling (src/BlockForest.jl:49-54), so holes aren't expressible today.

Honest ordering:

1. Cell mask on CartesianGrid — correctness; unblocks the microvascular model. Already this issue's scope note.
2. Block-level pruning on BlockForest — allow a forest with holes. Small change, captures most of the memory and launch win.
3. Cell-level compaction — buys only the remaining partially-live boundary blocks.

Against a block-pruned forest, the marginal return on step 3 is far below 20×. Recording that so nobody scopes step 3 expecting the paper's number.

DESIGN.md §9 already made this call once in a different guise — it rejects a cell-based octree because "hanging nodes everywhere break the dense-array leaf stencils," in favour of fixed-size leaf blocks where "every leaf operator runs unchanged per block." Block-level pruning is that same principle applied to participation instead of resolution. Cell-level compaction is the thing that argument was against. Not a veto — the fill fractions are extreme enough that step 3 may still earn its keep — but it should be argued, not assumed.

Pointers for whoever builds this

• Carp2Cartesian — open-source tetrahedral-VTK → uniform Cartesian voxel converter, same group.
• WGLPackedCartesian — their reference compaction implementation.
• The Strocchi cohort is openly published, so a real anatomical test fixture is a download rather than a project. That is a far better acceptance test than a synthetic hemisphere for the performance claims — keep the hemisphere for the correctness ones.

Related: #69 gets the geometric-accuracy half of this (the mask is the φ ∈ {0,1} special case of a phase field), and #59 is the GPU prerequisite.

---

_{🤖 Claude read two cardiac GPU papers so this issue could stop guessing at its own constants.}