Update devdocs solver extension API and fix sync docstrings

docs/src/devdocs/index.md

@@ -43,68 +43,120 @@ prob = OperatorSplittingProblem(f, u0, tspan)
 
     The API is not stable yet and subject to breaking changes.
 
-To add a new solver just define two new structs, one for the algorithm description and one for the algorithm cache and dispatch internal functions, as follows:
+To add a new solver, define two structs -- one describing the algorithm, one for its
+cache -- and dispatch three internal functions on them:
+
+- `init_cache(f, alg; uprev, u)` builds the cache for one node of the splitting tree.
+- `_perform_step!(parent, children, cache, dt)` advances that node by `dt`.
+- `alg_adaptive_order(alg)`, only if the algorithm is adaptive (see below).
+
+The algorithm struct has to carry the inner algorithms of the problem sequence in a
+field named `inner_algs`, because that is how the tree of integrators is built
+alongside the tree of split functions. The cache is where a scheme keeps the buffers
+it needs beyond `u`/`uprev`.
 
 ```julia
 using SciMLBase, OrdinaryDiffEqOperatorSplitting
+import OrdinaryDiffEqOperatorSplitting as OS
+
 struct MySimpleFirstOrderAlgorithm{InnerAlgorithmTypes} <:
-       OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingAlgorithm
-    inner_algs::InnerAlgorithmTypes # Tuple of solver for the problem sequence
+    OS.AbstractOperatorSplittingAlgorithm
+    inner_algs::InnerAlgorithmTypes # Tuple of solvers for the problem sequence
 end
 
-struct MySimpleFirstOrderCache{uType, uprevType, iiType} <:
-       OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingCache
+struct MySimpleFirstOrderCache{uType, uprevType} <: OS.AbstractOperatorSplittingCache
     u::uType
     uprev::uprevType
-    inner_caches::iiType
 end
 
-function OrdinaryDiffEqOperatorSplitting.init_cache(
+function OS.init_cache(
         f::GenericSplitFunction, alg::MySimpleFirstOrderAlgorithm;
-        uprev::AbstractArray, u::AbstractVector,
-        inner_caches,
-        alias_uprev = true,
-        alias_u = false
-)
-    @assert length(inner_caches) == 2
-    _uprev = alias_uprev ? uprev : SciMLBase.recursivecopy(uprev)
-    _u = alias_u ? u : SciMLBase.recursivecopy(u)
-    return MySimpleFirstOrderAlgorithmCache(_u, _uprev, inner_caches)
+        uprev::AbstractArray, u::AbstractVector
+    )
+    # `u` and `uprev` are the buffers of *this* node; the integrator owns them and
+    # the cache only keeps references. Allocate additional buffers with
+    # `similar(u)` if the scheme needs them.
+    return MySimpleFirstOrderCache(u, uprev)
 end
+```
 
-@inline function OrdinaryDiffEqOperatorSplitting.advance_solution_to!(
-        outer_integrator::OperatorSplittingIntegrator, subintegrators::Tuple,
-        solution_indices::Tuple, synchronizers::Tuple,
-        cache::MySimpleFirstOrderAlgorithmCache, tnext)
-    # We assume that the integrators are already synced
-    (;inner_caches) = cache
-
-    # Advance first subproblem
-    OrdinaryDiffEqOperatorSplitting.forward_sync_subintegrator!(
-        outer_integrator, subintegrators[1], solution_indices[1], synchronizers[1])
-    OrdinaryDiffEqOperatorSplitting.advance_solution_to!(
-        outer_integrator, subintegrators[1], solution_indices[1],
-        synchronizers[1], inner_caches[1], tnext)
-    if subintegrators[1].sol.retcode ∉
-       (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success)
-        return
-    end
-    OrdinaryDiffEqOperatorSplitting.backward_sync_subintegrator!(
-        outer_integrator, subintegrators[1], solution_indices[1], synchronizers[1])
-
-    # Advance second subproblem
-    OrdinaryDiffEqOperatorSplitting.forward_sync_subintegrator!(
-        outer_integrator, subintegrators[2], solution_indices[2], synchronizers[2])
-    OrdinaryDiffEqOperatorSplitting.advance_solution_to!(
-        outer_integrator, subintegrators[2], solution_indices[2],
-        synchronizers[2], inner_caches[2], tnext)
-    if subintegrators[2].sol.retcode ∉
-       (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success)
+The stepping function receives the node whose step is being performed (`parent`) and
+the tuple of its child integrators, which may be leaf `DEIntegrator`s or nested
+`SplitSubIntegrator`s -- the same code handles both. Everything else a step needs is
+reachable from `parent`: `parent.child_solution_indices[i]` are the indices of the
+`i`-th child in this node's solution vector, and `parent.child_synchronizers[i]` is its
+synchronizer.
+
+Advancing one child means synchronizing into it, stepping it, and synchronizing back:
+
+```julia
+function advance_one_child!(parent, child, i, dt)
+    idxs = parent.child_solution_indices[i]
+    sync = parent.child_synchronizers[i]
+
+    OS.forward_sync_subintegrator!(parent, child, idxs, sync)
+    OS.advance_solution_by!(parent, child, dt)
+    if OS._child_failed(child)
+        # Signal the failure and return immediately. `step_footer!` runs the failure
+        # escalation protocol from here: retry at this node if it is adaptive,
+        # otherwise escalate to the parent. Never keep stepping after a failure --
+        # the state of the failed child is meaningless.
+        parent.force_stepfail = true
         return
     end
-    OrdinaryDiffEqOperatorSplitting.backward_sync_subintegrator!(
-        outer_integrator, subintegrators[2], solution_indices[2], synchronizers[2])
+    OS.backward_sync_subintegrator!(parent, child, idxs, sync)
+    return
+end
 
-    # Done :)
+function OS._perform_step!(
+        parent, children::Tuple, cache::MySimpleFirstOrderCache, dt
+    )
+    advance_one_child!(parent, children[1], 1, dt)
+    parent.force_stepfail && return
+
+    advance_one_child!(parent, children[2], 2, dt)
+    parent.force_stepfail && return
+
+    # Done :) The solution of the step is in `parent.u`; `parent.uprev` holds the
+    # state at the beginning of the step and must be left untouched, as rollback
+    # after a rejected step restores from it.
+    return
 end
 ```
+
+This example is written for exactly two operators. A scheme that works for any number
+of them loops over `children` with `Unrolled.@unroll`, as the built-in algorithms in
+`src/solver.jl` do; a plain `for` loop over the heterogeneously typed tuple would be
+type unstable.
+
+### Adaptive algorithms
+
+A splitting node runs a step size controller only if its algorithm both declares
+itself adaptive and produces an error estimate. Such an algorithm has to
+
+1. define `SciMLBase.isadaptive(::MyAlgorithm) = true`,
+2. define [`OrdinaryDiffEqOperatorSplitting.alg_adaptive_order`](@ref), the order of
+   its error estimator, and
+3. write the tolerance-scaled error estimate to `parent.EEst` at the end of
+    `_perform_step!`, following the OrdinaryDiffEq convention that a step is accepted
+    when `EEst <= 1`.
+
+Only do the last step when `parent.controller_cache !== nothing`: the node may have
+been configured non-adaptive, in which case no controller consumes the estimate. The
+tolerances and the norm to scale with live in `parent.opts`:
+
+```julia
+if parent.controller_cache !== nothing
+    (; abstol, reltol, internalnorm) = parent.opts
+    @. residual = error_of_the_step / (abstol + max(abs(parent.u), abs(parent.uprev)) * reltol)
+    parent.EEst = internalnorm(residual, parent.t + dt)
+end
+```
+
+See [`PalindromicPairLieTrotterGodunov`](@ref) in `src/solver.jl` for a complete
+example, and [Adaptive time stepping](@ref) for how the two layers of adaptivity
+interact.
+
+```@docs
+OrdinaryDiffEqOperatorSplitting.alg_adaptive_order
+```

src/utils.jl

@@ -50,14 +50,13 @@ function sync_vectors!(a, b)
 end
 
 """
-    forward_sync_subintegrator!(parent_integrator::OperatorSplittingIntegrator, inner_integrator::DEIntegrator, solution_indices, sync)
+    forward_sync_subintegrator!(parent_integrator::AnySplitIntegrator, inner_integrator::DEIntegrator, solution_indices, sync)
 
 This function is responsible of copying the solution and parameters of the parent integrator and the synchronized subintegrators with the information given into the inner integrator.
 If the inner integrator is synchronized with other inner integrators using `sync`, the function `forward_sync_external!` shall be dispatched for `sync`.
 The `sync` object is passed from the outside and is the main entry point to dispatch custom types on for parameter synchronization.
 The `solution_indices` are indices into the parent integrators solution vectors.
 """
-
 function forward_sync_subintegrator!(
         parent::AnySplitIntegrator,
         child::DEIntegrator,
@@ -106,14 +105,13 @@ end
 
 
 """
-    backward_sync_subintegrator!(parent_integrator::OperatorSplittingIntegrator, inner_integrator::DEIntegrator, solution_indices, sync)
+    backward_sync_subintegrator!(parent_integrator::AnySplitIntegrator, inner_integrator::DEIntegrator, solution_indices, sync)
 
 This function is responsible of copying the solution of the inner integrator back into parent integrator and the synchronized subintegrators.
 If the inner integrator is synchronized with other inner integrators using `sync`, the function `backward_sync_external!` shall be dispatched for `sync`.
 The `sync` object is passed from the outside and is the main entry point to dispatch custom types on for parameter synchronization.
 The `solution_indices` are indices in the parent integrators solution vectors.
 """
-
 function backward_sync_subintegrator!(
         parent::AnySplitIntegrator,
         child::DEIntegrator,