From b771a1ffedbcdf1d625730a9e14c3f9adf3aef4b Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 3 Aug 2026 11:26:00 -0400 Subject: [PATCH] Update devdocs solver extension API and fix sync docstrings The "Adding Solvers" example documented `advance_solution_to!`, which no longer exists, along with an `init_cache` signature taking `inner_caches` and alias flags. Rewrite it against the current API: `init_cache(f, alg; uprev, u)` plus `_perform_step!(parent, children, cache, dt)`, which reads the solution indices and synchronizers off `parent` and signals failure via `force_stepfail` instead of inspecting child retcodes inline. Add a section on what an adaptive splitting algorithm has to provide (`isadaptive`, `alg_adaptive_order`, writing `EEst`). Also attach the `forward_sync_subintegrator!` and `backward_sync_subintegrator!` docstrings: a blank line between each docstring and its function left both undocumented, so the `@docs` blocks on the devdocs page rendered empty. Their signature lines said `OperatorSplittingIntegrator` where the methods take `AnySplitIntegrator`. Co-Authored-By: Claude Opus 5 --- docs/src/devdocs/index.md | 146 ++++++++++++++++++++++++++------------ src/utils.jl | 6 +- 2 files changed, 101 insertions(+), 51 deletions(-) diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md index 82708fc..af5d07f 100644 --- a/docs/src/devdocs/index.md +++ b/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 +``` diff --git a/src/utils.jl b/src/utils.jl index 642fd79..6e4d05f 100644 --- a/src/utils.jl +++ b/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,