feat(sciml): build the planned SciMLOperators jac_prototype adapter (DESIGN §7) #65

Status

Already scoped, not built. README.md:164-166 and DESIGN.md §7 describe a thin FunctionOperator-style adapter as the jac_prototype hook for implicit OrdinaryDiffEq, in a future ext/…SciMLExt.jl. This issue is to give it a real downstream forcing function rather than to propose it fresh.

The trait prerequisite has landed: #89 made PreparedOperator/PreparedForest declare islinear/isconstant/isselfadjoint/isdiagonal explicitly (src/linalg.jl:75-82) and gave PreparedComposed/PreparedAdjoint their trait methods (src/linalg.jl:308-309, :337-338), so the wrapper can forward traits from the prepared operator without falling through to a default it did not earn.

Why now

The microvascular transport model marches its coupled system with KenCarp47 + KrylovJL_GMRES. To hand OrdinaryDiffEq a matrix-free Jacobian it needs a SciMLOperators-shaped object, and MFO's LinearizedOp / EnzymeLinearizedOp is exactly that object wearing the wrong interface. Without the adapter, the downstream package is pushed into one of two bad options, and has taken both:

• AutoFiniteDiff everywhere, because its preallocated Float64 work buffers reject ForwardDiff.Duals — so it gets finite-difference Jacobians even though MFO can supply exact ones via Enzyme.
• Assemble a sparse W = I − γJ and ILU(0) it, with a hand-built Jacobian sparsity pattern and a Curtis–Powell–Reid colouring, purely to get a preconditioner into LinearSolve. That is a full sparse-matrix detour bolted onto a matrix-free stack, and it exists only because the matrix-free path couldn't reach the integrator.

Both disappear if prepare/LinearizedOp can be presented as a SciML operator.

Minimum useful shape

jac_prototype only, as §7 already scopes it — enough that ODEFunction(f; jac_prototype = sciml(prepare(L))) works and KrylovJL_GMRES calls straight into the zero-alloc flat mul!. Lazy algebra and adjoints through the SciML interface can wait; the integrator hook is what unblocks.

Two pieces of the SciML contract are on that path, though, and are in scope here (moved from #73, which keeps the SplitODEFunction + solver-search layer above this one):

• linearize! exposed as update_coefficients!. An implicit integrator refreshes the Jacobian at the current state through update_coefficients!(J, u, p, t). The package already has exactly those semantics in linearize!(J, u) (src/operators/linearize.jl:110-113 for LinearizedOp, ext/MatrixFreeOperatorsEnzymeExt.jl:66-69 for EnzymeLinearizedOp): it overwrites J.u0.data in place and keeps the existing prepare valid — refresh without re-prepare. The adapter's update_coefficients! should forward to it rather than reinvent it.
• ldiv! on MultigridPreconditioner. It exposes only size/eltype/mul! (src/multigrid.jl:378-395), which is Krylov.jl's preconditioner convention; LinearSolve's Pl/Pr (and therefore the precs argument) want ldiv!. Without it the multigrid V-cycle cannot stand in for the downstream ILU.

Pairs naturally with #62 — a Robin-aware operator_diagonal is what makes MultigridPreconditioner usable as the precs argument, and between the two the downstream ILU path can be retired rather than ported.

—
🤖 Claude gave this issue two more chores from #73.
_{claude-opus-5-5[1m] (high) · claude-code 2.1.285 · unreviewed}