diff --git a/Project.toml b/Project.toml index f0c2129..0829855 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" @@ -29,7 +30,8 @@ RecursiveArrayTools = "3.39.0, 4" SafeTestsets = "0.1.0" SciMLBase = "2.77.0, 3.1" SciMLIterators = "1" -SciMLTesting = "2.1" +SciMLLogging = "2" +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..7686b9f 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -24,8 +24,8 @@ makedocs( collapselevel = 1 ), sitename = "OrdinaryDiffEqOperatorSplitting.jl", - doctest = false, - warnonly = true, + doctest = 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 82708fc..63a7b80 100644 --- a/docs/src/devdocs/index.md +++ b/docs/src/devdocs/index.md @@ -1,110 +1,87 @@ -# Developer documentation +# Developer Documentation -## Synchronizers API +!!! warning "Developer extension interface v1" + These interfaces are for packages implementing operator-splitting solvers + or synchronizers. They are not end-user APIs. Version 1 may change only in + a breaking release of OrdinaryDiffEqOperatorSplitting.jl. -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. +## Solver Extensions -```@docs -OrdinaryDiffEqOperatorSplitting.NoExternalSynchronization -OrdinaryDiffEqOperatorSplitting.forward_sync_subintegrator! -OrdinaryDiffEqOperatorSplitting.backward_sync_subintegrator! -OrdinaryDiffEqOperatorSplitting.need_sync -OrdinaryDiffEqOperatorSplitting.sync_vectors! -``` - -## Adding Synchronizers - -!!! warning - - The API is not stable yet and subject to breaking changes. +A solver extension defines an algorithm type with an `inner_algs` tuple, a +concrete cache, `init_cache`, and `_perform_step!`. The tuple shape must mirror +the corresponding `GenericSplitFunction`. A step implementation synchronizes +each child before and after advancing it, and sets `parent.force_stepfail` if a +child reports failure. -You need to provide dispatches for - -```@docs; canonical=false -OrdinaryDiffEqOperatorSplitting.forward_sync_subintegrator! -OrdinaryDiffEqOperatorSplitting.backward_sync_subintegrator! +```@docs +OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingAlgorithm +OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingCache +OrdinaryDiffEqOperatorSplitting.init_cache +OrdinaryDiffEqOperatorSplitting._perform_step! ``` -with your custom synchronizer object and add it to the split function construction as follows: +The following sweeps the operators in reverse order, which is a complete +first-order splitting in its own right (it coincides with `LieTrotterGodunov` +exactly when the operators commute): ```julia -f1, f2 = generate_individual_functions() # assuming 3 unknowns each -i1, i2 = generate_solution_indices() # e.g. ([1,2,3], Int[]) -synchronizer_tree = generate_my_synchronizer_tree() # e.g. (MySynchronizer([1,2,3]), NoExternalSynchronization()) -f = GenericSplitFunction((f1, f2), (i1, i2), synchronizer_tree) -u0 = [-1.0, 1.0, 0.0] -tspan = (0.0, 1.0) -prob = OperatorSplittingProblem(f, u0, tspan) -``` - -## Adding Solvers - -!!! warning - - 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: +import OrdinaryDiffEqOperatorSplitting as OS -```julia -using SciMLBase, OrdinaryDiffEqOperatorSplitting -struct MySimpleFirstOrderAlgorithm{InnerAlgorithmTypes} <: - OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingAlgorithm - inner_algs::InnerAlgorithmTypes # Tuple of solver for the problem sequence +struct ReverseLieTrotterGodunov{AlgTupleType} <: OS.AbstractOperatorSplittingAlgorithm + inner_algs::AlgTupleType end -struct MySimpleFirstOrderCache{uType, uprevType, iiType} <: - OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingCache +struct ReverseLieTrotterGodunovCache{uType, uprevType} <: OS.AbstractOperatorSplittingCache u::uType uprev::uprevType - inner_caches::iiType end -function OrdinaryDiffEqOperatorSplitting.init_cache( - f::GenericSplitFunction, alg::MySimpleFirstOrderAlgorithm; +function OS.init_cache( + f::GenericSplitFunction, alg::ReverseLieTrotterGodunov; 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) + ) + return ReverseLieTrotterGodunovCache(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) - return +function OS._perform_step!( + parent, children::Tuple, cache::ReverseLieTrotterGodunovCache, dt + ) + for i in reverse(eachindex(children)) + child = children[i] + 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) + parent.force_stepfail = true + return nothing + end + OS.backward_sync_subintegrator!(parent, child, idxs, sync) end - OrdinaryDiffEqOperatorSplitting.backward_sync_subintegrator!( - outer_integrator, subintegrators[2], solution_indices[2], synchronizers[2]) - - # Done :) + return nothing end ``` + +The algorithm is then used like any built-in one, with one inner algorithm per +leaf operator: + +```julia +alg = ReverseLieTrotterGodunov((Euler(), Euler())) +sol = solve(OperatorSplittingProblem(f, u0, tspan), alg; dt = 0.01) +``` + +## Synchronizer Extensions + +Synchronizers communicate parameters and solution buffers between parent and +child integrators. Define the external synchronization methods for the custom +synchronizer type, then supply one synchronizer for every operator in a +`GenericSplitFunction`. + +```@docs +OrdinaryDiffEqOperatorSplitting.NoExternalSynchronization +OrdinaryDiffEqOperatorSplitting.forward_sync_subintegrator! +OrdinaryDiffEqOperatorSplitting.backward_sync_subintegrator! +OrdinaryDiffEqOperatorSplitting.need_sync +OrdinaryDiffEqOperatorSplitting.sync_vectors! +``` diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index 1cb9578..7a3179a 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,12 @@ 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 + # `DEVerbosity()` is what DiffEqBase itself defaults `verbose` to; spelling it + # out keeps this off the non-public `DEFAULT_VERBOSE` binding. + _inner_verbose(verbose::Bool) = verbose ? + DiffEqBase.DEVerbosity() : + DiffEqBase.DEVerbosity(SciMLLogging.None()) + const DEFAULT_VERBOSITY = DiffEqBase.DEVerbosity() else const DEFAULT_VERBOSITY = false end @@ -38,9 +42,63 @@ _is_verbose(verbose) = true end abstract type AbstractOperatorSplitFunction <: SciMLBase.AbstractODEFunction{true} end + +""" + AbstractOperatorSplittingAlgorithm + +Abstract supertype for operator-splitting algorithms. + +## Extension interface v1 + +Solver extensions must provide an `inner_algs` tuple matching the associated +[`GenericSplitFunction`](@ref), an [`init_cache`](@ref) method, and an +[`_perform_step!`](@ref) method. This is a developer extension interface, not +an end-user API. It may change in a breaking release of this package. +""" abstract type AbstractOperatorSplittingAlgorithm end + +""" + AbstractOperatorSplittingCache + +Abstract supertype for caches used by operator-splitting algorithms. + +## Extension interface v1 + +Define a concrete cache subtype and return it from [`init_cache`](@ref) when +implementing an [`AbstractOperatorSplittingAlgorithm`](@ref). This is a +developer extension interface, not an end-user API. It may change in a +breaking release of this package. +""" abstract type AbstractOperatorSplittingCache end +""" + init_cache(f::GenericSplitFunction, alg::AbstractOperatorSplittingAlgorithm; uprev, u) + +Create the cache for an operator-splitting algorithm. + +## Extension interface v1 + +Extensions must dispatch on their concrete algorithm type and return a concrete +[`AbstractOperatorSplittingCache`](@ref). `u` and `uprev` are the mutable local +solution buffers for the current node. This is a developer extension interface, +not an end-user API. It may change in a breaking release of this package. +""" +function init_cache end + +""" + _perform_step!(parent, children::Tuple, cache::AbstractOperatorSplittingCache, dt) + +Advance one operator-splitting step for a developer-defined algorithm. + +## Extension interface v1 + +Extensions must synchronize every child before and after advancing it and set +`parent.force_stepfail` when a child fails. This is a developer extension +interface, not an end-user API. It may change in a breaking release of this +package. +""" +function _perform_step! end + @inline SciMLBase.isadaptive(::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..7dcce4b 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 +`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/function.jl b/src/function.jl index b06b4b5..a4d9010 100644 --- a/src/function.jl +++ b/src/function.jl @@ -2,7 +2,18 @@ GenericSplitFunction(functions::Tuple, solution_indices::Tuple) GenericSplitFunction(functions::Tuple, solution_indices::Tuple, synchronizers::Tuple) -This type of function describes a set of connected inner functions in mass-matrix form, as usually found in operator splitting procedures. +Construct an operator tree for an [`OperatorSplittingProblem`](@ref). + +`functions` contains ODE operators or nested `GenericSplitFunction`s. +`solution_indices[i]` gives the entries of the parent solution used by +`functions[i]`. The three-argument form accepts one synchronizer per operator; +the two-argument form uses [`NoExternalSynchronization`](@ref) for every +operator. + +# Arguments +- `functions`: Tuple of ODE operators or nested split functions. +- `solution_indices`: Tuple of parent-solution index collections. +- `synchronizers`: Tuple of synchronization objects, one for each operator. """ struct GenericSplitFunction{fSetType <: Tuple, idxSetType <: Tuple, sSetType <: Tuple} <: AbstractOperatorSplitFunction @@ -36,7 +47,13 @@ num_operators(f::GenericSplitFunction) = length(f.functions) """ NoExternalSynchronization() -Indicator that no synchronization between parameters and solution vectors is necessary. +Marker indicating that a split operator requires no external parameter +synchronization. + +## Extension interface v1 + +This marker is for synchronizer implementations. It is not an end-user API and +may change in a breaking release of this package. """ struct NoExternalSynchronization end diff --git a/src/integrator.jl b/src/integrator.jl index e6d5499..b9e23ab 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -789,7 +789,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) @@ -901,7 +901,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 @@ -914,7 +914,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 diff --git a/src/problem.jl b/src/problem.jl index 35365db..4743e0e 100644 --- a/src/problem.jl +++ b/src/problem.jl @@ -1,5 +1,19 @@ """ OperatorSplittingProblem(f::AbstractOperatorSplitFunction, u0, tspan, p::Tuple) + +Define an ODE problem whose right-hand side is an operator tree. Use it with +[`LieTrotterGodunov`](@ref) or [`StrangMarchuk`](@ref) through the SciML +`init`/`solve` interface. + +# Arguments +- `f`: [`GenericSplitFunction`](@ref) describing the operator tree. +- `u0`: Initial state of the full system. +- `tspan`: Integration interval. +- `p`: Optional tuple of operator parameters. Null parameters are created when + it is omitted. + +# Keywords +Additional keywords are stored with the problem and forwarded by the solver. """ mutable struct OperatorSplittingProblem{ fType <: AbstractOperatorSplitFunction, uType, tType, pType <: Tuple, K, diff --git a/src/solver.jl b/src/solver.jl index 156313a..aac786c 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -2,10 +2,14 @@ # Lie-Trotter-Godunov operator splitting # --------------------------------------------------------------------------- """ - LieTrotterGodunov <: AbstractOperatorSplittingAlgorithm + LieTrotterGodunov(inner_algs) -First-order sequential operator splitting algorithm attributed to -[Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite). +Construct a first-order sequential operator-splitting algorithm attributed to +[Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite). Each step advances the +operators in tuple order. + +# Arguments +- `inner_algs`: Tuple containing one SciML algorithm for each leaf operator. """ struct LieTrotterGodunov{AlgTupleType} <: AbstractOperatorSplittingAlgorithm inner_algs::AlgTupleType # Tuple of timesteppers for inner problems @@ -61,7 +65,7 @@ end # Strang-Marchuk operator splitting # --------------------------------------------------------------------------- """ - StrangMarchuk <: AbstractOperatorSplittingAlgorithm + StrangMarchuk(inner_algs) Second-order symmetric (palindromic) operator splitting algorithm attributed to [Str:1968:ccd,Mar:1971:tsm](@cite). @@ -71,6 +75,9 @@ For ``N`` operators the scheme performs ``A_1(\\Delta t/2) \\to \\cdots \\to A_{N-1}(\\Delta t/2) \\to A_N(\\Delta t) \\to A_{N-1}(\\Delta t/2) \\to \\cdots \\to A_1(\\Delta t/2)`` achieving second-order accuracy through symmetry. + +# Arguments +- `inner_algs`: Tuple containing one SciML algorithm for each leaf operator. """ struct StrangMarchuk{AlgTupleType} <: AbstractOperatorSplittingAlgorithm inner_algs::AlgTupleType # Tuple of timesteppers for inner problems diff --git a/src/utils.jl b/src/utils.jl index 642fd79..3e3ed3c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -24,11 +24,15 @@ function tstops_and_saveat_heaps(t0, tf, tstops, saveat) end """ - need_sync(a, b) + need_sync(a::AbstractVector, b::AbstractVector) -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. +Determine whether two solution buffers require synchronization. + +## Extension interface v1 + +Synchronizer implementations may specialize this function for their storage +types. It is a developer extension interface, not an end-user API, and may +change in a breaking release of this package. """ need_sync @@ -40,7 +44,12 @@ 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 `b` into `a` when [`need_sync`](@ref) reports distinct storage. + +## Extension interface v1 + +This is a developer extension interface, not an end-user API. It may change in +a breaking release of this package. """ function sync_vectors!(a, b) if need_sync(a, b) && a !== b @@ -50,14 +59,17 @@ function sync_vectors!(a, b) end """ - forward_sync_subintegrator!(parent_integrator::OperatorSplittingIntegrator, inner_integrator::DEIntegrator, solution_indices, sync) + forward_sync_subintegrator!(parent::AnySplitIntegrator, child::DEIntegrator, solution_indices, synchronizer) -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 a child integrator from its parent before a substep. + +Custom synchronizers should define `forward_sync_external!` for their type. +## Extension interface v1 + +This is a developer extension interface, not an end-user API. It may change in +a breaking release of this package. +""" function forward_sync_subintegrator!( parent::AnySplitIntegrator, child::DEIntegrator, @@ -106,14 +118,17 @@ end """ - backward_sync_subintegrator!(parent_integrator::OperatorSplittingIntegrator, inner_integrator::DEIntegrator, solution_indices, sync) + backward_sync_subintegrator!(parent::AnySplitIntegrator, child::DEIntegrator, solution_indices, synchronizer) -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 a child integrator back into its parent after a substep. +Custom synchronizers should define `backward_sync_external!` for their type. + +## Extension interface v1 + +This is a developer extension interface, not an end-user API. It may change in +a breaking release of this package. +""" function backward_sync_subintegrator!( parent::AnySplitIntegrator, child::DEIntegrator, @@ -178,7 +193,7 @@ function OrdinaryDiffEqCore.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 diff --git a/test/operator_splitting_api.jl b/test/operator_splitting_api.jl index ba944a9..5902c64 100644 --- a/test/operator_splitting_api.jl +++ b/test/operator_splitting_api.jl @@ -1,4 +1,5 @@ using OrdinaryDiffEqOperatorSplitting +import OrdinaryDiffEqOperatorSplitting as OS using Test import SciMLBase: ReturnCode @@ -71,6 +72,84 @@ _sub2_iter_factor(::LieTrotterGodunov) = 1 _sub2_iter_factor(::StrangMarchuk) = 1 _sub2_iter_factor(::PalindromicPairLieTrotterGodunov) = 2 +# Mock solver extension exercising the documented developer extension interface v1: +# an algorithm carrying `inner_algs`, a concrete cache, `init_cache`, and +# `_perform_step!`. Sweeping the operators in reverse order is still a complete +# first-order splitting. It does not give a different answer than +# LieTrotterGodunov on the reference problem -- `trueA` is a multiple of the +# identity, so the two operators commute and the sweep order cancels out. +struct ReverseLieTrotterGodunov{AlgTupleType} <: OS.AbstractOperatorSplittingAlgorithm + inner_algs::AlgTupleType +end + +struct ReverseLieTrotterGodunovCache{uType, uprevType} <: OS.AbstractOperatorSplittingCache + u::uType + uprev::uprevType +end + +function OS.init_cache( + f::GenericSplitFunction, alg::ReverseLieTrotterGodunov; + uprev::AbstractArray, u::AbstractVector, + ) + return ReverseLieTrotterGodunovCache(u, uprev) +end + +function OS._perform_step!( + parent, children::Tuple, cache::ReverseLieTrotterGodunovCache, dt + ) + for i in reverse(eachindex(children)) + child = children[i] + 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) + parent.force_stepfail = true + return nothing + end + OS.backward_sync_subintegrator!(parent, child, idxs, sync) + end + return nothing +end + +@testset "solver extension interface v1" begin + split_f = GenericSplitFunction((f1, f2), ([1, 2, 3], [1, 3])) + prob = OperatorSplittingProblem(split_f, u0, tspan) + alg = ReverseLieTrotterGodunov((Euler(), Euler())) + + integrator = DiffEqBase.init(prob, alg; dt = 0.01, verbose = false) + @test integrator.cache isa ReverseLieTrotterGodunovCache + DiffEqBase.solve!(integrator) + @test integrator.sol.retcode == ReturnCode.Success + + function err(dt) + integ = DiffEqBase.init(prob, alg; dt = dt, verbose = false) + DiffEqBase.solve!(integ) + return maximum(abs, integ.u - trueu) + end + @test err(0.1) / err(0.01) ≈ 10 rtol = 0.1 +end + +@testset "nested rollback restores child buffers" begin + f1dofs = [1, 2, 3] + f2dofs = [1, 3] + f3dofs = [1, 2] + nested_f = GenericSplitFunction((f3, f3), (f3dofs, f3dofs)) + split_f = GenericSplitFunction((f1, nested_f), (f1dofs, f2dofs)) + prob = OperatorSplittingProblem(split_f, u0, (0.0, 1.0)) + alg = LieTrotterGodunov((Euler(), LieTrotterGodunov((Euler(), Euler())))) + integrator = DiffEqBase.init(prob, alg; dt = 0.1, adaptive = false, verbose = false) + + nested_child = integrator.child_subintegrators[2] + nested_child.u .= 0 + nested_child.uprev .= 0 + OS.rollback_children!(integrator) + + @test nested_child.u == integrator.u[nested_child.solution_indices] + @test nested_child.uprev == nested_child.u +end + # --------------------------------------------------------------------------- # Tests 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 cd55cef..6e191b1 100644 --- a/test/qa/qa.jl +++ b/test/qa/qa.jl @@ -1,47 +1,24 @@ using SciMLTesting using OrdinaryDiffEqOperatorSplitting -using JET -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 - :promote_tspan, # 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: `Base.Broadcast` marks + # `dotview`, `broadcastable` and `BroadcastStyle` public but not + # `materialize!` or the `Broadcasted` type they dispatch on, and there + # is no alias for either. src/config_tree.jl needs both to give + # `opt[...] .= x` its subtree-fill meaning. + :Broadcasted, :materialize!, + # Public from OrdinaryDiffEqCore 4.13; the [compat] floor is 4.4 and + # 4.12 is the newest registered 4.x, so the check still resolves a + # version without them. Drop once the floor moves past the release. + :fix_dt_at_bounds!, :handle_tstop!, + # https://github.com/SciML/OrdinaryDiffEq.jl/pull/4111 makes this + # public alongside the rest of the per-algorithm controller defaults. + :failfactor_default, ), ), ),