diff --git a/Project.toml b/Project.toml index a83ae67..1b9ea66 100644 --- a/Project.toml +++ b/Project.toml @@ -12,6 +12,7 @@ OrdinaryDiffEqLowOrderRK = "1344f307-1e59-4825-a18e-ace9aa3fa4c6" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SciMLLogging = "a6db7da4-7206-11f0-1eab-35f2a5dbe1d1" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" Unrolled = "9602ed7d-8fef-5bc8-8597-8f21381861e8" @@ -28,8 +29,9 @@ PrecompileTools = "1.1" RecursiveArrayTools = "3.39.0, 4" SafeTestsets = "0.1.0" SciMLBase = "3.36" +SciMLLogging = "2" SciMLIterators = "1" -SciMLTesting = "2.1" +SciMLTesting = "2.4" SymbolicIndexingInterface = "0.3.36" Test = "1" TimerOutputs = "0.5.28, 1.0" diff --git a/docs/make.jl b/docs/make.jl index 378cf31..1b61757 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -24,8 +24,7 @@ makedocs( collapselevel = 1 ), sitename = "OrdinaryDiffEqOperatorSplitting.jl", - doctest = false, - warnonly = true, + checkdocs = :exports, draft = false, pages = Any[ "Home" => "index.md", diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md index 30387a5..518509e 100644 --- a/docs/src/devdocs/index.md +++ b/docs/src/devdocs/index.md @@ -1,5 +1,11 @@ # Developer documentation +!!! warning + This page specifies developer extension APIs for packages implementing + operator-splitting solvers. They are versioned for the OrdinaryDiffEq ecosystem, + but are not supported end-user APIs. Application code should use the APIs in the + API Reference instead. + ## Synchronizers API A key part of operator splitting algorithms is the synchronization logic. Parameters of one subproblem might need to be kept in sync with the solution of other subproblems and vice versa. To handle this efficiently OrdinaryDiffEqOperatorSplitting.jl provides a small set of utils. @@ -8,6 +14,8 @@ A key part of operator splitting algorithms is the synchronization logic. Parame OrdinaryDiffEqOperatorSplitting.NoExternalSynchronization OrdinaryDiffEqOperatorSplitting.forward_sync_subintegrator! OrdinaryDiffEqOperatorSplitting.backward_sync_subintegrator! +OrdinaryDiffEqOperatorSplitting.forward_sync_external! +OrdinaryDiffEqOperatorSplitting.backward_sync_external! OrdinaryDiffEqOperatorSplitting.need_sync OrdinaryDiffEqOperatorSplitting.sync_vectors! ``` @@ -37,78 +45,142 @@ tspan = (0.0, 1.0) prob = OperatorSplittingProblem(f, u0, tspan) ``` -## Adding Solvers +## Solver extension API + +```@docs +OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingAlgorithm +OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingCache +OrdinaryDiffEqOperatorSplitting.init_cache +OrdinaryDiffEqOperatorSplitting._perform_step! +OrdinaryDiffEqOperatorSplitting.advance_solution_by! +OrdinaryDiffEqOperatorSplitting.child_failed +OrdinaryDiffEqOperatorSplitting.alg_adaptive_order +OrdinaryDiffEqOperatorSplitting.splitting_interpolant +OrdinaryDiffEqOperatorSplitting.splitting_interpolant! +``` + +## Adding solvers !!! warning - The API is not stable yet and subject to breaking changes. + This is a developer extension API. It may change independently of the end-user + API and should not be used in application code. + +To add a new solver, define two structs -- one describing the algorithm, one for its +cache -- and dispatch the developer extension functions on them: -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: +- `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) + # A failed child must stop the remaining stages of this step. + 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 + +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 :) + # 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. hand the tolerance-scaled error estimate to `OrdinaryDiffEqCore.set_EEst!` at the + end of `_perform_step!`, following the OrdinaryDiffEq convention that a step is + accepted when the estimate is `<= 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) + OrdinaryDiffEqCore.set_EEst!(parent, internalnorm(residual, parent.t + dt)) +end +``` + +`set_EEst!` and its counterpart `OrdinaryDiffEqCore.get_EEst` are the public +OrdinaryDiffEqCore interface for the estimate; go through them rather than touching +the `EEst` field, whose location on the integrator is an implementation detail. + +See [`PalindromicPairLieTrotterGodunov`](@ref) in `src/solver.jl` for a complete +example, and [Adaptive time stepping](@ref) for how the two layers of adaptivity +interact. + ## Dense output Saving, `saveat` and continuous callback root-finding all go through a single hook, diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index 19aaf77..661dc24 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -7,7 +7,7 @@ import Unrolled: @unroll import BinaryHeaps -import SciMLBase, DiffEqBase +import SciMLBase, DiffEqBase, SciMLLogging import SciMLBase: ReturnCode import SciMLBase: DEIntegrator, NullParameters, isadaptive import SymbolicIndexingInterface: variable_symbols @@ -21,8 +21,10 @@ import OrdinaryDiffEqCore: OrdinaryDiffEqCore, isdtchangeable, # In OrdinaryDiffEq v7 / DiffEqBase v7, passing verbose::Bool to inner ODE # integrators is no longer supported. Convert Bool → DEVerbosity when available. @static if isdefined(DiffEqBase, :DEVerbosity) - _inner_verbose(verbose::Bool) = verbose ? DiffEqBase.DEFAULT_VERBOSE : DiffEqBase.DEVerbosity(DiffEqBase.None()) - const DEFAULT_VERBOSITY = DiffEqBase.DEFAULT_VERBOSE + _inner_verbose(verbose::Bool) = verbose ? + DiffEqBase.DEVerbosity(SciMLLogging.Minimal()) : + DiffEqBase.DEVerbosity(SciMLLogging.None()) + const DEFAULT_VERBOSITY = DiffEqBase.DEVerbosity(SciMLLogging.Minimal()) else const DEFAULT_VERBOSITY = false end @@ -37,10 +39,109 @@ _is_verbose(verbose) = true _is_verbose(::DiffEqBase.DEVerbosity{B}) where {B} = B end +""" + AbstractOperatorSplitFunction + +Abstract supertype for functions that define an operator-splitting problem. +Concrete subtypes must provide an operator tree and the local state indices used by +[`OperatorSplittingProblem`](@ref). End users normally construct a +[`GenericSplitFunction`](@ref) rather than implementing this interface directly. +""" abstract type AbstractOperatorSplitFunction <: SciMLBase.AbstractODEFunction{true} end + +""" + AbstractOperatorSplittingAlgorithm + +Developer-only abstract supertype for algorithms that advance an +[`OperatorSplittingProblem`](@ref). + +# Interface requirements +- Store an `inner_algs` tuple whose shape mirrors the associated + [`GenericSplitFunction`](@ref). +- Implement [`init_cache`](@ref) to construct a cache for every node. +- Implement [`_perform_step!`](@ref) to advance that node. + +This interface is versioned for OrdinaryDiffEq solver developers. It is not a +supported end-user API. +""" abstract type AbstractOperatorSplittingAlgorithm end + +""" + AbstractOperatorSplittingCache + +Developer-only abstract supertype for an algorithm's per-node cache. A concrete +subtype holds references to the node's `u` and `uprev` buffers and any additional +temporary storage required by the splitting scheme. Construct it from +[`init_cache`](@ref). + +This interface is versioned for OrdinaryDiffEq solver developers. It is not a +supported end-user API. +""" abstract type AbstractOperatorSplittingCache end +""" + init_cache(f::GenericSplitFunction, alg::AbstractOperatorSplittingAlgorithm; uprev, u) + +Construct the per-node cache for a developer-defined operator-splitting algorithm. + +# Arguments +- `f`: Operator tree at the node being initialized. +- `alg`: Algorithm used at that node. + +# Keyword Arguments +- `uprev`: Mutable state buffer for the preceding accepted state. +- `u`: Mutable state buffer for the state currently being advanced. + +# Returns +A concrete [`AbstractOperatorSplittingCache`](@ref). The cache must retain the +provided buffers by reference; the integrator owns their allocation and restores them +after rejected steps. + +This is a developer extension API, not a supported end-user API. +""" +function init_cache end + +""" + _perform_step!(parent, children::Tuple, cache::AbstractOperatorSplittingCache, dt) + +Advance one node of an operator-splitting tree by `dt`. + +# Arguments +- `parent`: Integrator node that owns the current solution and rollback buffers. +- `children`: Direct child integrators, either `DEIntegrator`s or nested splitting + nodes. +- `cache`: Cache returned by [`init_cache`](@ref) for `parent`'s algorithm. +- `dt`: Signed duration of the splitting step. + +# Interface requirements +- Synchronize a child before and after advancing it. +- Leave `parent.uprev` unchanged; rejection restores the node from that buffer. +- Set `parent.force_stepfail = true` and return immediately when a child fails. +- For adaptive algorithms, pass the tolerance-scaled local error estimate to + `OrdinaryDiffEqCore.set_EEst!` when `parent.controller_cache !== nothing`. + +This is a developer extension API, not a supported end-user API. +""" +function _perform_step! end + +""" + child_failed(child) + +Return whether a direct child integrator has failed while a developer-defined +splitting step is executing. + +# Arguments +- `child`: A leaf `DEIntegrator` or nested operator-splitting integrator. + +# Returns +`true` when the child cannot be used for the remainder of the current splitting +step. In that case, `_perform_step!` must set `parent.force_stepfail = true` and +return without synchronizing the failed state back to its parent. + +This is a developer extension API, not a supported end-user API. +""" +function child_failed end + @inline SciMLBase.isadaptive(::AbstractOperatorSplittingAlgorithm) = false @inline SciMLBase.isdiscrete(::AbstractOperatorSplittingAlgorithm) = false @inline isdtchangeable(alg::AbstractOperatorSplittingAlgorithm) = all(isdtchangeable.(alg.inner_algs)) diff --git a/src/config_tree.jl b/src/config_tree.jl index b6bcad2..2d53472 100644 --- a/src/config_tree.jl +++ b/src/config_tree.jl @@ -210,17 +210,17 @@ Base.setindex!(opt::TreeOption, v, node::SplitNode) = TreeOptionSubtree The target of a broadcast assignment into a [`TreeOption`](@ref), produced by -`Base.dotview`. Assigning into it writes one value to a node and all of its -descendants. +`Base.Broadcast.dotview`. Assigning into it writes one value to a node and all +of its descendants. """ struct TreeOptionSubtree{T} node::TreeOption{T} end -Base.dotview(opt::TreeOption) = TreeOptionSubtree(opt) -Base.dotview(opt::TreeOption, i::Integer, is::Integer...) = +Base.Broadcast.dotview(opt::TreeOption) = TreeOptionSubtree(opt) +Base.Broadcast.dotview(opt::TreeOption, i::Integer, is::Integer...) = TreeOptionSubtree(_option_node(opt, _path(i, is...))) -Base.dotview(opt::TreeOption, node::SplitNode) = +Base.Broadcast.dotview(opt::TreeOption, node::SplitNode) = TreeOptionSubtree(_option_node(opt, node.path)) const _Broadcasted = Base.Broadcast.Broadcasted @@ -228,14 +228,14 @@ const _Broadcasted = Base.Broadcast.Broadcasted # A TreeOption is not a collection as far as broadcasting is concerned. Without this # the default `broadcastable` tries to `collect` it, and `opt .op= x` fails with a # MethodError about iteration instead of the explanation in `_broadcast_value`. -Base.broadcastable(opt::TreeOption) = Ref(opt) +Base.Broadcast.broadcastable(opt::TreeOption) = Ref(opt) # Broadcast assignment into a TreeOption only ever means "give every node of this # subtree the same value". Anything else that lowers to the same call is rejected # rather than silently given a made up meaning -- see `_broadcast_value`. -Base.materialize!(dest::TreeOptionSubtree, bc::_Broadcasted) = +Base.Broadcast.materialize!(dest::TreeOptionSubtree, bc::_Broadcasted) = _fill_subtree!(dest.node, _broadcast_value(bc)) -Base.materialize!(dest::TreeOption, bc::_Broadcasted) = +Base.Broadcast.materialize!(dest::TreeOption, bc::_Broadcasted) = _fill_subtree!(dest, _broadcast_value(bc)) function _broadcast_value(bc::_Broadcasted) diff --git a/src/integrator.jl b/src/integrator.jl index 0d68db0..e5befd1 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -592,9 +592,9 @@ function _subreinit_child!( end # --------------------------------------------------------------------------- -# handle_tstop! +# _handle_tstop! # --------------------------------------------------------------------------- -function OrdinaryDiffEqCore.handle_tstop!(integrator::AnySplitIntegrator) +function _handle_tstop!(integrator::AnySplitIntegrator) if SciMLBase.has_tstop(integrator) # The heaps store raw times; comparisons happen in tdir-space so that # "ahead"/"behind" is direction independent. @@ -724,7 +724,7 @@ end # re-run it -- a transiently failed non-adaptive child would otherwise stay failed # forever and turn every escalated failure fatal. function _reset_child_failure!(child::SplitSubIntegrator) - _child_failed(child) && (child.status.retcode = ReturnCode.Default) + child_failed(child) && (child.status.retcode = ReturnCode.Default) child.last_step_failed = false child.force_stepfail = false return nothing @@ -763,7 +763,7 @@ function step_header!(integrator::AnySplitIntegrator) update_uprev!(integrator) end increment_iteration(integrator) - OrdinaryDiffEqCore.fix_dt_at_bounds!(integrator) + _fix_dt_at_bounds!(integrator) modify_dt_for_tstops!(integrator) integrator.force_stepfail = false return nothing @@ -894,7 +894,7 @@ function step_footer!(integrator::AnySplitIntegrator) end function abort_below_dtmin!(integrator::AnySplitIntegrator) - abs(integrator.dt) > abs(OrdinaryDiffEqCore.timedepentdtmin(integrator)) && return nothing + abs(integrator.dt) > abs(DiffEqBase.timedepentdtmin(integrator)) && return nothing _is_verbose(integrator.opts.verbose) && @warn("dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.") _set_retcode!(integrator, ReturnCode.DtLessThanMin) @@ -928,7 +928,7 @@ function DiffEqBase.solve!(integrator::OperatorSplittingIntegrator) step_footer!(integrator) SciMLBase.has_tstop(integrator) || break end - OrdinaryDiffEqCore.handle_tstop!(integrator) + _handle_tstop!(integrator) end SciMLBase.postamble!(integrator) integrator.sol.retcode != ReturnCode.Default && return integrator.sol @@ -965,7 +965,7 @@ function DiffEqBase.step!(integrator::AnySplitIntegrator) step_footer!(integrator) end end - OrdinaryDiffEqCore.handle_tstop!(integrator) + _handle_tstop!(integrator) return end @@ -991,10 +991,10 @@ function DiffEqBase.step!(integrator::AnySplitIntegrator, dt, stop_at_tdt = fals # root -- and leaving one in the heap once we sit exactly on it makes # the next header compute a zero step-to-tstop gap: the loop would # then spin at fixed `t` forever, growing the child tstop heaps. - OrdinaryDiffEqCore.handle_tstop!(integrator) + _handle_tstop!(integrator) end end - OrdinaryDiffEqCore.handle_tstop!(integrator) + _handle_tstop!(integrator) return nothing end @@ -1006,7 +1006,7 @@ function SciMLBase.check_error(integrator::OperatorSplittingIntegrator) integrator.sol.retcode != ReturnCode.Default return integrator.sol.retcode end - if DiffEqBase.NAN_CHECK(integrator.dtcache) || DiffEqBase.NAN_CHECK(integrator.dt) + if isnan(integrator.dtcache) || isnan(integrator.dt) _is_verbose(integrator.opts.verbose) && @warn("NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.") return ReturnCode.DtNaN @@ -1019,7 +1019,7 @@ function SciMLBase.check_error(integrator::SplitSubIntegrator) integrator.status.retcode != ReturnCode.Default return integrator.status.retcode end - if DiffEqBase.NAN_CHECK(integrator.dtcache) || DiffEqBase.NAN_CHECK(integrator.dt) + if isnan(integrator.dtcache) || isnan(integrator.dt) _is_verbose(integrator.opts.verbose) && @warn("NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.") return ReturnCode.DtNaN @@ -1081,7 +1081,7 @@ end Dense output of a single operator splitting step, evaluated in the step-local coordinate `Θ = (t - tprev) / dt`, where `y₀`/`y₁` are the state at `tprev`/`t` and -`dt` is the step that was taken. `D` is the requested derivative order and `idxs` +`dt` is the step that was taken. `D` is derivative order `0` or `1`, and `idxs` selects components (`nothing` for all of them). This is the extension point for higher order dense output: dispatch on the @@ -1091,14 +1091,37 @@ algorithm's cache type. The fallback for any splitting_interpolant, splitting_interpolant! splitting_interpolant( - integrator, ::AbstractOperatorSplittingCache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}} -) where {D} = SciMLBase.linear_interpolant(Θ, dt, y₀, y₁, idxs, Val{D}) + integrator, ::AbstractOperatorSplittingCache, Θ, dt, y₀, y₁, idxs, deriv +) = _linear_interpolant(Θ, dt, y₀, y₁, idxs, deriv) function splitting_interpolant!( out, integrator, ::AbstractOperatorSplittingCache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}} ) where {D} - SciMLBase.linear_interpolant!(out, Θ, dt, y₀, y₁, idxs, Val{D}) + return _linear_interpolant!(out, Θ, dt, y₀, y₁, idxs, Val{D}) +end + +_linear_interpolant(Θ, dt, y₀, y₁, ::Nothing, ::Type{Val{0}}) = + @. (1 - Θ) * y₀ + Θ * y₁ +_linear_interpolant(Θ, dt, y₀, y₁, idxs, ::Type{Val{0}}) = + @. (1 - Θ) * y₀[idxs] + Θ * y₁[idxs] +_linear_interpolant(Θ, dt, y₀, y₁, ::Nothing, ::Type{Val{1}}) = (y₁ - y₀) / dt +_linear_interpolant(Θ, dt, y₀, y₁, idxs, ::Type{Val{1}}) = (y₁[idxs] - y₀[idxs]) / dt + +function _linear_interpolant!(out, Θ, dt, y₀, y₁, ::Nothing, ::Type{Val{0}}) + @. out = (1 - Θ) * y₀ + Θ * y₁ + return out +end +function _linear_interpolant!(out, Θ, dt, y₀, y₁, idxs, ::Type{Val{0}}) + @. out = (1 - Θ) * y₀[idxs] + Θ * y₁[idxs] + return out +end +function _linear_interpolant!(out, Θ, dt, y₀, y₁, ::Nothing, ::Type{Val{1}}) + @. out = (y₁ - y₀) / dt + return out +end +function _linear_interpolant!(out, Θ, dt, y₀, y₁, idxs, ::Type{Val{1}}) + @. out = (y₁[idxs] - y₀[idxs]) / dt return out end @@ -1261,6 +1284,29 @@ end handle_callbacks!(::SplitSubIntegrator) = nothing +@generated function _apply_ith_continuous_callback!( + integrator, + time, + upcrossing, + event_idx, + callback_index, + callbacks::NTuple{N, Union{SciMLBase.ContinuousCallback, SciMLBase.VectorContinuousCallback}}, + ) where {N} + expression = :(throw(BoundsError(callbacks, callback_index))) + for i in N:-1:1 + expression = quote + if callback_index == $i + return DiffEqBase.apply_callback!( + integrator, callbacks[$i], time, upcrossing, event_idx + ) + else + $expression + end + end + end + return expression +end + function handle_callbacks!(integrator::OperatorSplittingIntegrator) discrete_callbacks = integrator.callback.discrete_callbacks continuous_callbacks = integrator.callback.continuous_callbacks @@ -1277,7 +1323,7 @@ function handle_callbacks!(integrator::OperatorSplittingIntegrator) if event_occurred integrator.event_last_time = idx integrator.vector_event_last_time = event_idx - continuous_modified, saved_in_cb = OrdinaryDiffEqCore.apply_ith_callback!( + continuous_modified, saved_in_cb = _apply_ith_continuous_callback!( integrator, time, upcrossing, event_idx, idx, continuous_callbacks ) # The step was cut at the event and the state may have jumped, so the @@ -1309,7 +1355,7 @@ end # Step size control # # A splitting node runs a controller only if it is adaptive and its algorithm -# provides an error estimate (written to the node's `EEst` by `_perform_step!`); +# provides an error estimate (reported by `_perform_step!` via `set_EEst!`); # its `controller_cache` is `nothing` otherwise. The controllers themselves are # the OrdinaryDiffEqCore ones: `setup_controller_cache` resolves a controller's # knobs against this algorithm and mints the mutable per-solve state that @@ -1480,7 +1526,23 @@ function __step!(integrator::AnySplitIntegrator) return nothing end -# Entry point: dispatch to the algorithm's advance_solution_by! +""" + advance_solution_by!(parent, child, dt) + +Advance one direct child of an operator-splitting node by the signed duration `dt`. + +# Arguments +- `parent`: The node that owns `child` and coordinates the splitting step. +- `child`: A leaf `DEIntegrator` or nested splitting node. +- `dt`: Signed duration to apply to the child. + +# Returns +`nothing`. This function only advances the child; a custom +[`_perform_step!`](@ref) must synchronize the child before and after the call and use +[`child_failed`](@ref) before accepting its state. + +This is a developer extension API, not a supported end-user API. +""" function advance_solution_by!(integrator::AnySplitIntegrator, dt) return advance_solution_by!(integrator, integrator.cache, dt) end @@ -1519,7 +1581,7 @@ end # so it already exhausted its own adaptation); `ReturnCode.Default` if none is. @unroll function _fatal_child_retcode(children::Tuple) @unroll for child in children - if _child_failed(child) && _child_is_adaptive(child) + if child_failed(child) && _child_is_adaptive(child) return _failure_retcode(child) end end @@ -1528,7 +1590,7 @@ end @unroll function _first_failed_child_retcode(children::Tuple) @unroll for child in children - _child_failed(child) && return _failure_retcode(child) + child_failed(child) && return _failure_retcode(child) end return ReturnCode.Failure end diff --git a/src/solver.jl b/src/solver.jl index 156313a..db76d35 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -48,7 +48,7 @@ end @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync) @timeit_debug "time solve" advance_solution_by!(parent, child, dt) - if _child_failed(child) + if child_failed(child) parent.force_stepfail = true return end @@ -111,7 +111,7 @@ end @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync) @timeit_debug "time solve" advance_solution_by!(parent, child, step_dt) - if _child_failed(child) + if child_failed(child) parent.force_stepfail = true return end @@ -132,7 +132,7 @@ end @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync) @timeit_debug "time solve" advance_solution_by!(parent, child, half_dt) - if _child_failed(child) + if child_failed(child) parent.force_stepfail = true return end @@ -227,7 +227,7 @@ function _ppltg_advance_child!(parent, child, i, dt) @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync) @timeit_debug "time solve" advance_solution_by!(parent, child, dt) - if _child_failed(child) + if child_failed(child) parent.force_stepfail = true return end @@ -282,7 +282,7 @@ function _perform_step!( (; abstol, reltol, internalnorm) = parent.opts @. uforward = (parent.u - uforward) / (abstol + max(abs(parent.u), abs(parent.uprev)) * reltol) - parent.EEst = internalnorm(uforward, parent.t + dt) + OrdinaryDiffEqCore.set_EEst!(parent, internalnorm(uforward, parent.t + dt)) end return end diff --git a/src/utils.jl b/src/utils.jl index ab3e99d..e7c059e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -27,9 +27,16 @@ end """ need_sync(a, b) -Determines whether it is necessary to synchronize two objects with any -solution information. A possible reason when no synchronization is necessary -might be that the vectors alias each other in memory. +Return whether copying solution information from `b` into `a` is necessary. + +# Arguments +- `a`: Destination vector or view. +- `b`: Source vector or view. + +# Returns +`false` when both arguments share the same backing storage and copying would be +redundant; `true` otherwise. Extend this function for custom array wrappers whose +aliasing relationship cannot be determined by the built-in vector methods. """ need_sync @@ -41,7 +48,15 @@ need_sync(a::SubArray, b::SubArray) = a.parent !== b.parent """ sync_vectors!(a, b) -Copies the information in `b` into `a` if synchronization is necessary. +Copy solution information from `b` into `a` when [`need_sync`](@ref) determines that +their storage does not alias. + +# Arguments +- `a`: Destination vector or view. +- `b`: Source vector or view. + +# Returns +`nothing`. If no copy is required, `a` is left unchanged. """ function sync_vectors!(a, b) if need_sync(a, b) && a !== b @@ -51,14 +66,25 @@ 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. -""" +Synchronize one child with its parent immediately before advancing that child. + +# Arguments +- `parent_integrator`: Splitting node that owns the current full solution. +- `inner_integrator`: Direct leaf child to receive its local state. +- `solution_indices`: Indices selecting the child state in the parent's solution. +- `sync`: Synchronizer object for any parameter or cross-child synchronization. + +# Interface requirements +The built-in method copies the parent's selected state into `inner_integrator` and +marks its derivative cache stale. To synchronize parameters or coupled child state, +implement `forward_sync_external!(parent_integrator, inner_integrator, sync)` for the +concrete `sync` type. +# Returns +`nothing`. +""" function forward_sync_subintegrator!( parent::AnySplitIntegrator, child::DEIntegrator, @@ -102,14 +128,25 @@ 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. -""" +Synchronize one advanced child back into its parent. + +# Arguments +- `parent_integrator`: Splitting node that owns the full solution. +- `inner_integrator`: Direct leaf child whose state is copied back. +- `solution_indices`: Indices selecting the child state in the parent's solution. +- `sync`: Synchronizer object for parameter or cross-child synchronization. + +# Interface requirements +The built-in method copies the child state into the parent's selected state. To +synchronize parameters or coupled child state, implement +`backward_sync_external!(parent_integrator, inner_integrator, sync)` for the concrete +`sync` type. +# Returns +`nothing`. +""" function backward_sync_subintegrator!( parent::AnySplitIntegrator, child::DEIntegrator, @@ -134,6 +171,50 @@ end # These handle parameter synchronisation via the `sync` object. # --------------------------------------------------------------------------- +""" + forward_sync_external!(parent_integrator, inner_integrator, sync) + +Synchronize external state into a child immediately before that child advances. + +# Arguments +- `parent_integrator`: Splitting parent that owns the full state. +- `inner_integrator`: Direct child receiving the synchronized state. +- `sync`: Synchronizer object selected in the [`GenericSplitFunction`](@ref) tree. + +# Interface requirements +Implement this method for a concrete `sync` type when parameters or coupled child +state must be refreshed before the child solves. The default +[`NoExternalSynchronization`](@ref) implementation is a no-op. + +# Returns +`nothing`. + +This is a developer extension API, not a supported end-user API. +""" +function forward_sync_external! end + +""" + backward_sync_external!(parent_integrator, inner_integrator, sync) + +Synchronize external state after a child has advanced. + +# Arguments +- `parent_integrator`: Splitting parent that owns the full state. +- `inner_integrator`: Direct child whose result may update coupled state. +- `sync`: Synchronizer object selected in the [`GenericSplitFunction`](@ref) tree. + +# Interface requirements +Implement this method for a concrete `sync` type when parameters or coupled child +state must be refreshed after the child solves. The default +[`NoExternalSynchronization`](@ref) implementation is a no-op. + +# Returns +`nothing`. + +This is a developer extension API, not a supported end-user API. +""" +function backward_sync_external! end + # NoExternalSynchronization: no-op for all parent/child combinations forward_sync_external!(parent::DEIntegrator, child::DEIntegrator, ::NoExternalSynchronization) = nothing backward_sync_external!(parent::DEIntegrator, child::DEIntegrator, ::NoExternalSynchronization) = nothing @@ -170,11 +251,11 @@ function synchronize_solution_with_parameters!( end # Time stuff -function OrdinaryDiffEqCore.fix_dt_at_bounds!(integrator::AnySplitIntegrator) +function _fix_dt_at_bounds!(integrator::AnySplitIntegrator) # dtmin/dtmax are magnitudes; clamp |dt| and restore the direction. dtmin wins # over dtmax if the two conflict. dtmax = abs(integrator.opts.dtmax) - dtmin = abs(OrdinaryDiffEqCore.timedepentdtmin(integrator)) + dtmin = abs(DiffEqBase.timedepentdtmin(integrator)) integrator.dt = tdir(integrator) * max(min(abs(integrator.dt), dtmax), dtmin) return nothing end @@ -197,15 +278,15 @@ function validate_time_point(parent, child::DEIntegrator) end # --------------------------------------------------------------------------- -# _child_failed: check whether a child failed +# child_failed: check whether a child failed # # Leaves are checked through `SciMLBase.check_error`, not their stored retcode: # fixed-step leaf integrators complete `step!` with NaN state without flagging it # themselves, and the failure has to be caught *before* the surrounding step is # accepted so that the escalation protocol can retry from a clean `uprev`. # --------------------------------------------------------------------------- -_child_failed(child::DEIntegrator) = +child_failed(child::DEIntegrator) = SciMLBase.check_error(child) ∉ (ReturnCode.Default, ReturnCode.Success) -_child_failed(child::SplitSubIntegrator) = +child_failed(child::SplitSubIntegrator) = child.status.retcode ∉ (ReturnCode.Default, ReturnCode.Success) diff --git a/test/failure_escalation.jl b/test/failure_escalation.jl index 7f48347..4c50d9e 100644 --- a/test/failure_escalation.jl +++ b/test/failure_escalation.jl @@ -125,11 +125,11 @@ PPLTG = PalindromicPairLieTrotterGodunov leaf = integ.child_subintegrators[2] leaf.last_stepfail = true @test SciMLBase.check_error(leaf) == ReturnCode.ConvergenceFailure - @test OS._child_failed(leaf) # the eager detection sees it ... + @test OS.child_failed(leaf) # the eager detection sees it ... OS.reject_step!(integ) # ... and the retry's rollback clears it @test !leaf.last_stepfail - @test !OS._child_failed(leaf) + @test !OS.child_failed(leaf) @test SciMLBase.check_error(leaf) == ReturnCode.Success end diff --git a/test/operator_splitting_api.jl b/test/operator_splitting_api.jl index ba944a9..bc8fa74 100644 --- a/test/operator_splitting_api.jl +++ b/test/operator_splitting_api.jl @@ -8,6 +8,42 @@ using OrdinaryDiffEqTsit5 using ModelingToolkit using SciMLIterators: TimeChoiceIterator, intervals +const OS = OrdinaryDiffEqOperatorSplitting + +struct ExternalFirstOrderAlgorithm{AlgTuple} <: OS.AbstractOperatorSplittingAlgorithm + inner_algs::AlgTuple +end + +struct ExternalFirstOrderCache{U, UPrev} <: OS.AbstractOperatorSplittingCache + u::U + uprev::UPrev +end + +function OS.init_cache( + ::OS.GenericSplitFunction, ::ExternalFirstOrderAlgorithm; + uprev::AbstractArray, u::AbstractVector + ) + return ExternalFirstOrderCache(u, uprev) +end + +function OS._perform_step!( + parent, children::Tuple, ::ExternalFirstOrderCache, dt + ) + for i in eachindex(children) + child = children[i] + indices = parent.child_solution_indices[i] + synchronizer = parent.child_synchronizers[i] + OS.forward_sync_subintegrator!(parent, child, indices, synchronizer) + OS.advance_solution_by!(parent, child, dt) + if OS.child_failed(child) + parent.force_stepfail = true + return nothing + end + OS.backward_sync_subintegrator!(parent, child, indices, synchronizer) + end + return nothing +end + # --------------------------------------------------------------------------- # Reference problem # --------------------------------------------------------------------------- @@ -75,6 +111,33 @@ _sub2_iter_factor(::PalindromicPairLieTrotterGodunov) = 2 # --------------------------------------------------------------------------- # Tests # --------------------------------------------------------------------------- +@testset "generic developer solver interface" begin + f1dofs = [1, 2, 3] + f2dofs = [1, 3] + split_function = GenericSplitFunction((f1, f2), (f1dofs, f2dofs)) + problem = OperatorSplittingProblem(split_function, u0, (0.0, 1.0)) + algorithm = ExternalFirstOrderAlgorithm((Euler(), Euler())) + + integrator = DiffEqBase.init( + problem, algorithm; dt = 0.01, adaptive = false, verbose = false + ) + + @test integrator.cache isa ExternalFirstOrderCache + @test integrator.cache.u === integrator.u + @test integrator.cache.uprev === integrator.uprev + @test all(!OS.child_failed(child) for child in integrator.child_subintegrators) + + solution = DiffEqBase.solve!(integrator) + @test solution.retcode == DiffEqBase.ReturnCode.Success + + reference = DiffEqBase.init( + problem, LieTrotterGodunov((Euler(), Euler())); + dt = 0.01, adaptive = false, verbose = false + ) + DiffEqBase.solve!(reference) + @test integrator.u == reference.u +end + @testset "reinit and convergence" begin dt = 0.01π diff --git a/test/qa/Project.toml b/test/qa/Project.toml index ce90109..6c7332c 100644 --- a/test/qa/Project.toml +++ b/test/qa/Project.toml @@ -8,6 +8,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] Aqua = "0.8" JET = "0.9,0.10,0.11" -SciMLTesting = "2.1" +SciMLTesting = "2.4" Test = "1" julia = "1.10" diff --git a/test/qa/qa.jl b/test/qa/qa.jl index 42f9f80..175052d 100644 --- a/test/qa/qa.jl +++ b/test/qa/qa.jl @@ -5,48 +5,20 @@ using Test run_qa( OrdinaryDiffEqOperatorSplitting; - # target_defined_modules scopes the report to this package's own modules (the - # default target_modules=(pkg,) filter hides via-dependency-driven frames). - jet_kwargs = (; target_defined_modules = true, mode = :basic), ei_kwargs = (; - # Names re-exported through the SciML umbrella chain; accessed via a - # re-exporting dep rather than the owning package. - all_qualified_accesses_via_owners = (; - ignore = ( - :None, # owner SciMLLogging, via DiffEqBase - :timedepentdtmin, # owner DiffEqBase, via OrdinaryDiffEqCore - ), - ), all_qualified_accesses_are_public = (; ignore = ( - :__init, :__solve, :done, :postamble!, :solution_new_retcode, # SciMLBase - :DEFAULT_VERBOSE, :NAN_CHECK, :None, :ODE_DEFAULT_NORM, # DiffEqBase - :fix_dt_at_bounds!, :handle_tstop!, :increment_accept!, # OrdinaryDiffEqCore - :increment_reject!, :initialize_d_discontinuities, :initialize_saveat, # OrdinaryDiffEqCore - :initialize_tstops, :timedepentdtmin, :IController, # OrdinaryDiffEqCore - # Controller-cache protocol names. Several are `public` only in newer - # OrdinaryDiffEqCore (e.g. post_newton_controller! from 4.7); keep them - # all ignored so QA stays valid across the whole [compat] range even - # though the pinned QA manifest resolves a newer version. - :setup_controller_cache, :reinit_controller!, :post_newton_controller!, # OrdinaryDiffEqCore - :get_EEst, :set_EEst!, :get_current_adaptive_order, # OrdinaryDiffEqCore - :gamma_default, :failfactor_default, :AbstractControllerCache, # OrdinaryDiffEqCore - :apply_ith_callback!, # OrdinaryDiffEqCore - :promote_tspan, # SciMLBase - # Shared linear interpolation kernels, reused so that the - # integrator's own dense output and `sol(t)` (which goes through - # SciMLBase's `LinearInterpolation`) cannot drift apart. - :linear_interpolant, :linear_interpolant!, # SciMLBase - # Broadcast extension points a TreeOption implements (src/config_tree.jl). - :Broadcasted, :broadcastable, :dotview, :materialize!, # Base - ), - ), - all_explicit_imports_are_public = (; - ignore = ( - :isdtchangeable, # OrdinaryDiffEqCore - # Public only in newer OrdinaryDiffEqCore; see the note above. - :stepsize_controller!, :step_accept_controller!, # OrdinaryDiffEqCore - :step_reject_controller!, :accept_step_controller, # OrdinaryDiffEqCore + # Broadcast overloading has no public spelling for these two: + # `Base.Broadcast` marks `dotview`, `broadcastable` and + # `BroadcastStyle` public but not `materialize!` or the + # `Broadcasted` type they dispatch on, and neither has an alias. + # src/config_tree.jl needs both to give `opt[...] .= x` its + # subtree-fill meaning. + :Broadcasted, :materialize!, + # https://github.com/SciML/OrdinaryDiffEq.jl/pull/4111 makes this + # public alongside the rest of the per-algorithm controller + # defaults. Drop once a release carrying it is registered. + :failfactor_default, ), ), ),