interpolation and saving support

Project.toml

@@ -27,7 +27,7 @@ OrdinaryDiffEqTsit5 = "1.1.0, 2"
 PrecompileTools = "1.1"
 RecursiveArrayTools = "3.39.0, 4"
 SafeTestsets = "0.1.0"
-SciMLBase = "2.77.0, 3.1"
+SciMLBase = "3.36"
 SciMLIterators = "1"
 SciMLTesting = "2.1"
 SymbolicIndexingInterface = "0.3.36"

docs/src/devdocs/index.md

@@ -108,3 +108,42 @@ end
     # Done :)
 end
 ```
+
+## Dense output
+
+Saving, `saveat` and continuous callback root-finding all go through a single hook,
+so an algorithm only has to describe how to interpolate *within one of its steps*:
+
+```julia
+function OrdinaryDiffEqOperatorSplitting.splitting_interpolant(
+        integrator, cache::MySimpleFirstOrderCache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}}
+) where {D}
+    # Θ = (t - tprev) / dt is the step-local coordinate, y₀ = u(tprev), y₁ = u(t).
+    # ...
+end
+
+function OrdinaryDiffEqOperatorSplitting.splitting_interpolant!(
+        out, integrator, cache::MySimpleFirstOrderCache, Θ, dt, y₀, y₁, idxs,
+        ::Type{Val{D}}
+) where {D}
+    # In-place variant; must return `out`.
+end
+```
+
+Both fall back to linear interpolation for any `AbstractOperatorSplittingCache`, so
+implementing them is optional -- but with the fallback, saved output is first order
+even for a second-order scheme, and continuous callback event times are the exact
+roots of a straight chord across the step. One method improves both.
+
+The fallback is linear for a structural reason, and it is the same reason saving and
+callbacks live on the outer integrator alone: a splitting step advances its children
+sequentially over staggered subintervals, so a child's own interpolant describes a
+different sub-problem over a different interval, and they do not compose into an
+approximation of the split solution. Only the step endpoints, which the outer
+integrator owns, are states of the full split system -- an inner split is a stage,
+not a step.
+
+When implementing this, note that `Θ` is derived from
+`integrator.t - integrator.tprev`, **not** from `integrator.dt`: once a step is
+accepted, `step_accept_controller!` has already replaced `dt` with the step size
+proposed for the *next* step.

docs/src/usage/index.md

@@ -42,8 +42,9 @@ alg = LieTrotterGodunov(
     (Euler(), Euler())
 )
 
-# Right now OrdinaryDiffEqOperatorSplitting.jl does not implement the SciML solution interface,
-# but we can obtain intermediate solutions via the iterator interface.
+# OrdinaryDiffEqOperatorSplitting.jl implements only part of the SciML solution
+# interface (see "Saving and interpolation" below); intermediate solutions are most
+# directly obtained via the iterator interface.
 integrator = init(prob, alg, dt = 0.1)
 for (u, t) in TimeChoiceIterator(integrator, 0.0:0.5:1.0)
     @show t, u
@@ -65,6 +66,90 @@ for (u, t) in TimeChoiceIterator(integrator, 0.0:0.5:1.0)
 end
 ```
 
+## Saving and interpolation
+
+`solve!` returns a solution whose `t`/`u` are filled as the integration proceeds, so
+the usual saving keywords work and `sol(t)` interpolates between saved points:
+
+```julia
+integrator = init(prob, alg; dt = 0.1, saveat = 0.25)
+sol = solve!(integrator)
+
+sol.t          # [0.0, 0.25, 0.5, 0.75, 1.0]
+sol(0.3)       # interpolated between the saved points
+```
+
+The supported keywords are `saveat`, `save_everystep` (default `false`),
+`save_start`, `save_end` and `save_on`. They apply to the **outer** integrator only:
+the inner splits are stages rather than steps, so their intermediate states do not
+approximate the split solution at any time point.
+
+`saveat` points that fall strictly inside a step are filled from the step's
+interpolant, so requesting output never changes the sequence of steps and therefore
+never changes the splitting error.
+
+!!! note "The interpolant is first order"
+    A splitting step advances its children sequentially over staggered
+    subintervals, so their individual interpolants do not compose into an
+    approximation of the split solution. The dense output available at the outer
+    level is therefore linear between the step endpoints -- exact for the state a
+    `LieTrotterGodunov` step produces, but only first order for the second-order
+    schemes. If you need output that is accurate to the order of the scheme, pass
+    the times as `tstops` so that the integrator lands on them exactly:
+
+    ```julia
+    sol = solve!(init(prob, alg; dt = 0.1, tstops = [0.25, 0.5], saveat = [0.25, 0.5]))
+    ```
+
+    `dense = true` is accepted but has no effect: `sol(t)` already reproduces the
+    integrator's own interpolant from the saved points, so there is no per-step
+    interpolation data left to store. Use `save_everystep` or `saveat` to control
+    how finely that interpolant resolves the trajectory.
+
+## Callbacks and events
+
+The standard SciML callbacks work through the `callback` keyword:
+
+```julia
+# Stop as soon as the first component drops below a threshold.
+cb = ContinuousCallback((u, t, integrator) -> u[1] - 0.5, terminate!)
+sol = solve!(init(prob, alg; dt = 0.1, callback = cb))
+```
+
+The machinery is DiffEqBase's own, so `DiscreteCallback`, `ContinuousCallback`,
+`VectorContinuousCallback`, `CallbackSet` and anything built on top of them (for
+instance DiffEqCallbacks.jl) all behave as they do elsewhere in SciML.
+
+!!! note "Callbacks run on the outer integrator only"
+    A condition is evaluated once per **outer** step, after all the operators of
+    that step have been applied, and never between two inner splits: for the reason
+    given under "Saving and interpolation" above, those intermediate states are
+    stages and approximate the split solution at no time point, so there is nothing
+    meaningful for a condition to test or an `affect!` to modify at that level.
+
+    A consequence worth knowing: an `affect!` that modifies `integrator.u` is
+    propagated into every subintegrator before the next step, so modifying the state
+    from a callback is safe.
+
+### Accuracy of continuous events
+
+Event times are found by root-finding on the step's interpolant, which is linear
+(see "Saving and interpolation" above). Two consequences:
+
+  - The located event time is second order accurate in the step size, and is the
+    *exact* root of the linear interpolant over the step that brackets it. With
+    large steps -- adaptive splittings can grow the step considerably on smooth
+    problems -- the event time degrades accordingly. Cap it with `dtmax` when event
+    accuracy matters.
+  - An event that occurs and reverses **within a single step** cannot be detected,
+    because a linear interpolant has no interior extremum. Raising `interp_points`
+    does not help for the same reason, so setting `interp_points = 0` on the
+    callback avoids a sweep that cannot find anything the endpoints missed.
+
+Once the event time is located, the state there comes from the same interpolant and
+the whole subintegrator tree is re-anchored to it, so integration resumes
+consistently from the event.
+
 ## Configuring individual subintegrators
 
 `init` takes one value per keyword, which is not enough when the operators want

src/OrdinaryDiffEqOperatorSplitting.jl

@@ -42,6 +42,7 @@ abstract type AbstractOperatorSplittingAlgorithm end
 abstract type AbstractOperatorSplittingCache end
 
 @inline SciMLBase.isadaptive(::AbstractOperatorSplittingAlgorithm) = false
+@inline SciMLBase.isdiscrete(::AbstractOperatorSplittingAlgorithm) = false
 @inline isdtchangeable(alg::AbstractOperatorSplittingAlgorithm) = all(isdtchangeable.(alg.inner_algs))
 
 include("function.jl")