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

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.

"The same operator on every component" is literal, coefficients included: diffusion rejects a non-scalar-eltype κ (src/operators/diffusion.jl:545) and so does scaling(κ::Field) (src/operators/scaling.jl:48), so one scalar κ applies to every component. There is no way today to give each species its own D.

(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 (MFO does not do this yet: componentwise application shares one scalar D across all species, see above);
• 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 "per-species 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 neither exists yet:

Stacked(ops...)              # block-diagonal: op[i] acts on species i  (componentwise gives only op[i] == op[j], one shared κ)
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.

Adjacent precursor: open PR #107 (titled GHScalingOp; in its current diff, scaling(f, args...) building ScalingOp(func, args)) evaluates a pointwise coefficient from several fields on the fly. That is the scalar-diagonal special case of PointwiseBlock — f must return a real scalar, applied identically to every component — and it declares itself self-adjoint unconditionally (isselfadjoint(::ScalingOp) = true, adjoint_operator(S) = S), so it cannot carry the non-symmetric cross-species coupling above (O₂ → CO₂ is not matched by CO₂ → O₂). PointwiseBlock needs a declared per-cell block transpose instead.

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 corrected a claim that was optimistic about componentwise D.
_{claude-opus-5-5[1m] (high) · claude-code 2.1.285 · unreviewed}