feat(operators): stacked multi-field operator with cross-component coupling

[kylebeggs]

What's missing

MFO has multi-component fields but not multi-species operators. Field over Array{SVector{M,T}} gives ncomponents, component, and component-fastest flatten — but operators are written generically over element type, so they act componentwise: the same operator on every component, no cross-component block. Only Gradient/Divergence touch components at all, and vector_field(g) only ever builds SVector{N,T} where N is the grid dimension, not a species count.

(Naming trap worth recording: BlockField / PackedBlockField are spatial AMR leaf-blocks, not field blocks. They read like they solve this and they don't.)

Why (downstream)

The microvascular transport model marches a stacked tissue state [O₂; CO₂; H; Lf] on one shared grid. The structure is:

• diffusion — per-species, decoupled, different D per species (this part MFO already does, componentwise);
• reaction — a dense local block coupling all four species in every cell: O₂ consumption drives CO₂ production through the same term, CO₂ hydration produces the proton load H while consuming CO₂, and Lf couples to H/CO₂ via bicarbonate consumption and to O₂ via the PO₂-graded anaerobic source.

So it is "componentwise stencil + dense pointwise block", which is a common enough shape (any reaction–diffusion system, and most of SciML's SplitFunction use cases) that it feels like it belongs in the operator algebra rather than in every downstream RHS.

Today the downstream package hand-writes a fused kernel over a flat 4nc + 5nv state vector and hand-builds the Jacobian sparsity pattern to match, because there is no way to express "this operator is block-diagonal in space and dense in species".

Sketch

Two pieces, and the second is the one that doesn't exist:

Stacked(ops...)              # block-diagonal: op[i] acts on species i  (≈ what componentwise already gives)
PointwiseBlock(f!, nspecies) # dense M×M local coupling at every cell

operator_diagonal for the pair is what makes it worth doing — that is the hook MultigridPreconditioner needs, and getting a species-block diagonal is exactly the information a CPR-style two-stage preconditioner would want too.

Scope note

Sizing: is the species count a type parameter (SVector{M}, unrolled, isbits, GPU-friendly) or dynamic? For the downstream case M = 4 and static is strictly better, and static also keeps Adapt.adapt_structure trivial. Suggest starting static.

---

_{🤖 Claude drafted this. Kyle's opinions are his own; the file:line citations are mine.}

[kylebeggs]

Cross-reference: this is the natural shape of the reaction half of a SplitODEFunction, so it pairs with the SciML tracker #81 (in particular #73, which hands the linear half to the integrator). Left standalone rather than parented there because the operator itself — block-diagonal in space, dense in species, with an operator_diagonal a CPR-style preconditioner can use — is useful with or without SciML.

---

_{🤖 Beep boop — Claude drew the edge; Kyle is not stacking himself.}