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

Activity

1. kylebeggs

   added

   enhancementNew feature or requestNew feature or request

   on Sep 9, 2026
2. kylebeggs

   added sub-issues on Sep 9, 2026

   • fix(linalg): prepared operators declare their traits — PreparedOperator supertype, PreparedComposed/PreparedAdjoint islinear/isconstant #86
   • feat(sciml): expose operators as (Split)ODEFunction and search integrators per model #73
   • feat(sciml): build the planned SciMLOperators jac_prototype adapter (DESIGN §7) #65
   • docs: retract the OrdinaryDiffEq interop claim and document field-level apply! as the explicit-stepping idiom #87
3. kylebeggs

   added the

   Feature

   issue type on Sep 9, 2026
4. kylebeggs

   added 4 commits that reference this issue on Sep 10, 2026

   docs: retract the OrdinaryDiffEq interop claim, step examples via apply!

   ...

   9255717

   docs(design): explicit apply! idiom is a hand-rolled stepper, not ODE…

   ...

   e63982a

   docs: retract the OrdinaryDiffEq interop claim, step examples via apply!

   ...

   d3f12f3

   docs(design): explicit apply! idiom is a hand-rolled stepper, not ODE…

   ...

   a9b8ef9
5. kylebeggs commented on Sep 10, 2026

   

   kylebeggs

   on Sep 10, 2026

   CollaboratorAuthor

   More actions

   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.}