meta(sciml): expose prepared operators to the SciML integrator stack

[kylebeggs]

Umbrella for handing MFO operators to OrdinaryDiffEq / LinearSolve. Nothing SciML exists today: no weak dep, no ext/…SciMLExt.jl, and every example hand-rolls its integrator. README.md currently advertises OrdinaryDiffEq interop that is not there.

Sub-issues, in order

1. fix(linalg): prepared operators declare their traits — PreparedOperator supertype, PreparedComposed/PreparedAdjoint islinear/isconstant #86 — prepared-operator trait hygiene — PreparedOperator has no supertype and PreparedComposed / PreparedAdjoint declare no islinear / isconstant, so a wrapper forwarding traits gets the false default by accident. An invariant hazard independent of SciML; ship first.
2. feat(sciml): build the planned SciMLOperators jac_prototype adapter (DESIGN §7) #65 — the jac_prototype adapter (DESIGN §7). Enough that ODEFunction(f; jac_prototype = sciml(prepare(L))) works and KrylovJL_GMRES calls the zero-alloc flat mul!. The downstream forcing function is the microvascular model's KenCarp47 + sparse-ILU detour.
3. feat(sciml): expose operators as (Split)ODEFunction and search integrators per model #73 — SplitODEFunction for the linear half plus a per-model work-precision search over IMEX / exponential integrators (a benchmark artifact, not a package feature).
4. docs: retract the OrdinaryDiffEq interop claim and document field-level apply! as the explicit-stepping idiom #87 — docs — fix the README claim, and document field-level apply! as the explicit-stepping idiom (the prepared mul! leg pays 15–50% in flat copies that only Krylov needs).

Related, not parented: #63 (stacked multi-species operator) is what the reaction half of a split problem wants to look like; #62 (Robin) is what makes MultigridPreconditioner applicable to the downstream wall-coupled system.

Constraints

• Extension only. SciMLOperators / OrdinaryDiffEq are weak deps, nothing enters [deps] (CLAUDE.md, DESIGN.md §7).
• Do not route around the mul!(y, ::AbstractOperator, x::AbstractVector) guard; a convenience constructor must still go through prepare.
• The affine boundary_rhs term needs a declared home in a split problem.

---

_{🤖 Beep boop — Claude wired the tracker to the integrator; Kyle is not stepping himself forward in time.}

[kylebeggs]

Addendum: are the staging copies actually necessary, or just the shape we ended up with?

#87 / #90 established that the flat mul! path pays 15–50% in flat_to_interior! / interior_to_flat! staging, and steered explicit stepping to field-level apply! to dodge it. That routes around the cost. It does not establish that the cost has to exist. Three questions surfaced reviewing that PR, in increasing order of scope — recording them here because #1 blocks this umbrella directly and #2/#3 are the layer underneath it.

1. What actually stops a SciML explicit stepper from using apply!?

The RHS signature is a closure away, so the blocker is the state type, not the function. Worth establishing from scratch: what does OrdinaryDiffEq genuinely require of a state u — which of similar / copyto! / broadcast / norm / length / indexing, and which of those are needed only by adaptive stepping or by particular integrator families rather than by all of them? Field today defines similar and nothing else (src/Fields.jl:207). So: what is the minimal interface that makes Field a legal state, what does implementing it over a padded array whose interior is a strided view cost, and does an interior-only broadcast/norm sit cleanly with the interior-only-DOF invariant or fight it?

Whatever the answer, this is a gap that belongs on #73's list and currently isn't there.

2. Can the copies be removed from mul! instead of routed around?

Krylov.jl is generic over its vector type and lets you extend its kernel operations, so a Field-valued Krylov workspace is at least expressible. The question is whether it pays. Every workspace vector would become a padded array with a strided interior, so dots, axpys and norms trade contiguity (and GPU coalescing) for the saved copies, and GMRES(k) would carry k ghost rings it never fills. Two O(N) copies per matvec against a penalty on the 5–30 vector ops per iteration is not obviously a win in either direction — it wants a benchmark, not an argument.

The GPU case deserves separate measurement: there the staging is also two extra kernel launches, and on small grids launch overhead may swamp the bandwidth reasoning entirely.

3. Should the layout have been flat from the start?

That is, interior stored contiguously with ghosts held beside it in separate face buffers, rather than inline in a padded array. A solver vector could then alias the interior directly, and filling faces scales as O(surface) rather than two O(volume) passes. It is also plausibly the more natural fit for the #79 distributed work, since MPI wants contiguous face buffers regardless and inline ghosts have to be packed into them.

The cost lands on operator bodies. With ghosts inline, one broadcast over shifted slices is correct at every cell, boundary included — which is what the "leaf bodies stay array-level, device-agnostic, no per-backend code" structural rule rests on. With ghosts beside, the sweep splits into an interior core plus a boundary shell of up to 3ᴺ−1 regions, per leaf, per BC kind, and again for every declared adjoint and AD path. Whether that multiplier is as bad as it sounds is the open question, and it is worth checking against the one existing data point: laplacian_7pt_noflux (DESIGN §5) is already a flat halo-free primitive, deliberately kept as a separate body, pinned to homogeneous no-flux, and locked to the halo-based leaf by a parity test. That precedent is evidence about how far the pattern generalizes — in either direction.

Where the answers should land

• 1 is scoped and belongs on feat(sciml): expose operators as (Split)ODEFunction and search integrators per model #73 as a gap regardless of the verdict.
• 2 is a benchmark; CPU and GPU separately.
• 3 is a DESIGN question, not a SciML one. If the layout is right, it wants an alternatives-considered note so it stops resurfacing; if it is wrong, it wants its own issue rather than a comment buried under a SciML umbrella.

---

_{🤖 Beep boop — Claude filed the questions, not the answers; Kyle gets to keep the fun part.}