feat(sciml): expose operators as (Split)ODEFunction and search integrators per model

[kylebeggs]

From Dennis Ogiermann, 2026-08-14: "Is there any chance we can funnel this into an ODEFunction or SplitODEFunction and make a search over which solver works best for our individual models?"

Relationship to #65

#65 is the prerequisite and covers the implicit hook that DESIGN.md §7 (:639-651) already scopes: a thin FunctionOperator-style adapter so LinearizedOp can be handed to ODEFunction as jac_prototype. This issue is the layer above it, and it is the part Dennis actually asked for: presenting the operator itself as the linear half of a split problem, and then searching the integrator space per model. Do #65 first.

What exists today

An RHS closure and nothing else. src/linalg.jl:86-88 documents f!(du, u, p, t) = mul!(du, A, u); test/ode_rhs.jl exercises it with a hand-rolled RK4 and loads no ODE package. There is no SciML anywhere in [deps], [weakdeps], [extensions], or ext/.

Meanwhile every time-stepping example rolls its own integrator: explicit Euler in examples/heat_equation.jl:24-28, a hand-written Godunov split in examples/niederer_benchmark.jl:177-186, and examples/monodomain_amr.jl:61-72. None of them can reach an IMEX or exponential integrator, which for a stiff-diffusion / stiff-reaction monodomain problem is where the actual timestep savings live.

Also worth fixing here: README.md:147 advertises "the OrdinaryDiffEq.jl interop" in the documentation, and docs/ contains zero occurrences of ODE, OrdinaryDiffEq, or SciML. That claim is currently false.

The shape

A  = prepare(diffusion(g, κ), u)              # linear, constant, matrix-free
f2 = (du, u, p, t) -> reaction!(du, u, p, t)  # pointwise
prob = ODEProblem(SplitODEFunction(sciml(A), f2), u0, tspan, p)

SplitODEFunction is the right target rather than a plain ODEFunction because it is what unlocks the integrator families that suit this problem class: IMEX (KenCarp*, ARKODE), exponential integrators (ETDRK*, EPIRK*), and the split/LinearExponential methods. The reason the examples all hand-roll Godunov splitting is that there is no way to hand the linear half to a solver as an operator.

Then the search: a work-precision sweep across integrators, per model, comparing wall time at fixed accuracy against a reference solution. That belongs in benchmark/ or examples/, never in the package — same rule that keeps DifferentiationInterface out of [deps] per CLAUDE.md. The output is a work-precision diagram per model plus a recommended default, which is directly what Dennis asked for.

Gaps that block the wrapper

• PreparedOperator has no supertype (src/linalg.jl:24), so islinear(P), isconstant(P), and adjoint(P) are MethodErrors. A SciML operator needs at least the first two.
• PreparedComposed and PreparedAdjoint are <: AbstractOperator (src/linalg.jl:229,248) but declare no trait methods, so islinear on a prepared composed tree silently falls through to the false default at src/operators/abstract.jl:31. Any wrapper that forwards traits from P.op gets a wrong answer that degrades safely by accident rather than by design. This is an invariant hazard — "traits default to the weaker claim" only protects you if the weaker claim is also the declared one — and it should be fixed regardless of whether this issue proceeds.
• No update_coefficients!-shaped hook exists. The closest analogue in the package is linearize!(J, u) (src/operators/linearize.jl:110-113), which overwrites J.u0.data in place and keeps the existing prepare valid — exactly the "refresh without re-prepare" semantics SciML wants, and the natural template.
• isconstant is declared on every leaf and propagated through the algebra, and consumed nowhere in src/. It is the obvious contract for a SciML isconstant / caching decision and is currently inert.
• MultigridPreconditioner speaks the wrong convention. It exposes size/eltype/mul! (src/multigrid.jl:378-395), which is Krylov.jl's preconditioner interface; LinearSolve's Pl/Pr want ldiv!.
• The affine term needs a home. Inhomogeneous BCs live outside the linear map by design and come back through boundary_rhs. In a SplitODEFunction that is a constant additive term, and the wrapper has to either carry an affine slot or document that f2 must include it. Nothing supplies it today. (It also re-allocates on every call — src/linalg.jl:398-413.)
• Don't route around the guard. src/linalg.jl:540-547 deliberately makes mul!(y, ::AbstractOperator, x::AbstractVector) throw with "wrap the operator once with prepare(L, x)". A convenience constructor must not bypass it.
• Fix the README.md:147 claim, or make it true.

Scope note

Extension only — ext/MatrixFreeOperatorsSciMLExt.jl, weak dep, nothing added to [deps], per CLAUDE.md and DESIGN.md:50 ("Drop SciMLOperators as the core … a thin optional adapter"). The solver search is a benchmark artifact, not a package feature.

---

_{🤖 Claude wired the operator to the integrator on paper; Kyle is not quietly reviewing his own tracker.}