diff --git a/Project.toml b/Project.toml index fe3296b..f0c2129 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "OrdinaryDiffEqOperatorSplitting" uuid = "760fc936-9fa0-4281-9c1a-468957eedc89" -version = "0.3.6" +version = "0.4.0" authors = ["Dennis Ogiermann and contributors"] [deps] @@ -19,7 +19,7 @@ Unrolled = "9602ed7d-8fef-5bc8-8597-8f21381861e8" [compat] BinaryHeaps = "1.0.0" CommonSolve = "0.2.4" -DiffEqBase = "6.165.1, 7" +DiffEqBase = "7.5" ModelingToolkit = "11" OrdinaryDiffEqCore = "4.4" OrdinaryDiffEqLowOrderRK = "1.7, 2" diff --git a/_typos.toml b/_typos.toml index fbc9267..625e683 100644 --- a/_typos.toml +++ b/_typos.toml @@ -3,3 +3,5 @@ Strang = "Strang" # Citation key for Trotter (Tro:1959:psg) used in docstrings Tro = "Tro" +# Splitting variable name +BA = "BA" diff --git a/docs/Project.toml b/docs/Project.toml index e458a97..1039bf9 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -6,4 +6,4 @@ OrdinaryDiffEqOperatorSplitting = "760fc936-9fa0-4281-9c1a-468957eedc89" [compat] Documenter = "1.16.1" DocumenterCitations = "1.4.1" -OrdinaryDiffEqOperatorSplitting = "0.2.3, 0.3" +OrdinaryDiffEqOperatorSplitting = "0.2.3, 0.3, 0.4" diff --git a/docs/make.jl b/docs/make.jl index 0f968bf..378cf31 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -31,6 +31,7 @@ makedocs( "Home" => "index.md", "usage/index.md", "Theory Manual" => "topics/time-integration.md", + "Adaptive time stepping" => "topics/adaptivity.md", "api-reference/index.md", "devdocs/index.md", "references.md", diff --git a/docs/src/api-reference/index.md b/docs/src/api-reference/index.md index bb2d674..44b0610 100644 --- a/docs/src/api-reference/index.md +++ b/docs/src/api-reference/index.md @@ -17,6 +17,7 @@ GenericSplitFunction ```@docs LieTrotterGodunov StrangMarchuk +PalindromicPairLieTrotterGodunov ``` ## Per-node configuration diff --git a/docs/src/topics/adaptivity.md b/docs/src/topics/adaptivity.md new file mode 100644 index 0000000..c2f28c8 --- /dev/null +++ b/docs/src/topics/adaptivity.md @@ -0,0 +1,96 @@ +# Adaptive time stepping + +Two independent layers of a splitting tree can adapt their step sizes: + +- **Leaf solvers** adapt their own internal steps within each sub-solve, exactly as + they would outside of this package. Their behavior is governed by the tolerances + and options they receive. +- **Splitting nodes** adapt the splitting step itself. This requires an algorithm + that produces an error estimate of the splitting error; + [`PalindromicPairLieTrotterGodunov`](@ref) is currently the only one. It advances + with the average of the two mutually reversed Lie-Trotter sequences and uses half + their difference as the local error estimate. + +By default every node adapts exactly if its own algorithm can: for + +```julia +alg = PalindromicPairLieTrotterGodunov((Tsit5(), Euler())) +integrator = init(prob, alg; dt = 0.1) +``` + +the splitting node runs its step size controller, the `Tsit5` leaf adapts its inner +steps, and the `Euler` leaf steps fixed. Passing `adaptive = false` (or a +per-node [`TreeOption`](@ref)) overrides this. + +## The two tolerance layers + +`abstol` and `reltol` are understood at *every* node: + +- At a **splitting node** they scale the splitting error estimate: the node accepts + a step when `‖(u₁ - u₂)/2 ./ (abstol .+ max.(|u|) .* reltol)‖ ≤ 1`, where + `u₁, u₂` are the results of the two sequences of the palindromic pair. +- At a **leaf** they control the accuracy of the inner solves. + +A scalar passed to `init` configures the whole tree with the same value. A +[`TreeOption`](@ref) configures each node individually: + +```julia +reltol = TreeOption(f, 1.0e-8) # splitting node(s): tight +reltol[1] = 1.0e-10 # leaf 1: tighter still +reltol[2] = 1.0e-10 # leaf 2: tighter still +integrator = init(prob, alg; dt, reltol) +``` + +## Rule: inner tolerances must be tighter than the splitting tolerances + +The palindromic error estimator measures the **difference of two inner-solved +sequences**. Whatever error the inner solvers commit enters that difference as +noise the estimator cannot distinguish from splitting error. Two consequences: + +1. **The estimator is blind to inner error.** With coarse fixed-step inner solvers + (say `Euler()` stepping at the splitting step size) the overall method degrades + to the inner order while the controller happily reports small `EEst`. +2. **Loose inner tolerances choke the controller.** When a sub-problem has fast + internal dynamics (the typical reason for splitting it off), its solver commits + error on the order of its tolerance in *every* sub-solve, independent of the + splitting step size. If that exceeds the splitting tolerance — e.g. an inner + `abstol` 10× *larger* than the splitting `abstol` — then `EEst > 1` no matter + how small the splitting `dt` becomes. The controller keeps rejecting and + shrinking `dt` without the estimate improving, until the solve aborts with + `ReturnCode.DtLessThanMin` (or, with a tiny `dtmin`, grinds down to absurdly + small splitting steps). + +As a rule of thumb, keep the leaf tolerances at least **one to two orders of +magnitude tighter** than the splitting tolerances (adaptive leaves), or the fixed +inner steps well below the splitting step. The safe default is the scalar spread — +identical tolerances everywhere are already borderline; never configure leaves +*looser* than their splitting node. + +## Symptoms and causes + +| Symptom | Likely cause | +|---|---| +| `ReturnCode.DtLessThanMin`, `dt` collapsed, solution up to that point looks fine | Inner tolerances looser than the splitting tolerances (noise floor in the estimator), or genuinely unreachable tolerances | +| Splitting `dt` grows to the full interval immediately | The operators (nearly) commute, the splitting error is ≈ 0; harmless — the inner solvers carry the accuracy | +| Result visibly less accurate than the splitting tolerances suggest | Inner solves under-resolved: the estimator cannot see inner error (see rule above) | +| A few rejections right after `init` or `reinit!` | Initial `dt` too large for the tolerance; harmless, the controller recovers | +| Immediate abort with `DtLessThanMin` although tolerances look consistent | A scalar `dtmin` travels to the leaves too and may forbid the sub-steps they need. Restrict it to the splitting node with a `TreeOption` (`dtmin = TreeOption(f, 0.0); dtmin[] = 1e-3`) | + +## Controllers + +An adaptive splitting node runs an `OrdinaryDiffEqCore` step size controller +(default: `IController`). The standard knobs (`qmin`, `qmax`, `gamma`, +`qsteady_min`, `qsteady_max`, `failfactor`) can be passed to `init` — they are +folded into the default controller — or a controller object can be passed +explicitly via `controller` (per node via a `TreeOption`), e.g. a `PIController` +whose memory smooths the step size sequence. + +## Failure handling + +A failing *adaptive* node (leaf or splitting node) is fatal: it already exhausted +its own step size adaptation, and its return code propagates to the root. A failing +*non-adaptive* node escalates the failure to the nearest adaptive ancestor, which +rolls the whole subtree back and retries with a `failfactor`-shrunken step — +shrinking the effective step of every non-adaptive descendant — until it either +succeeds or falls below `dtmin`. Without any adaptive ancestor the integration +stops with the escalated return code. diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index 5857b97..1cb9578 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -15,7 +15,8 @@ import SymbolicIndexingInterface: variable_symbols import RecursiveArrayTools import OrdinaryDiffEqCore: OrdinaryDiffEqCore, isdtchangeable, - stepsize_controller!, step_accept_controller!, step_reject_controller! + stepsize_controller!, step_accept_controller!, step_reject_controller!, + accept_step_controller # In OrdinaryDiffEq v7 / DiffEqBase v7, passing verbose::Bool to inner ODE # integrators is no longer supported. Convert Bool → DEVerbosity when available. @@ -50,7 +51,8 @@ include("integrator.jl") include("solver.jl") include("utils.jl") -export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, StrangMarchuk +export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, StrangMarchuk, + PalindromicPairLieTrotterGodunov export SplitNode, TreeOption include("precompilation.jl") diff --git a/src/config_tree.jl b/src/config_tree.jl index d118373..b6bcad2 100644 --- a/src/config_tree.jl +++ b/src/config_tree.jl @@ -367,6 +367,30 @@ signed_dt_tree(config::ConfigTree, tdir, ::Type{tType}) where {tType} = ConfigTr map(child -> signed_dt_tree(child, tdir, tType), config.children) ) +""" + default_adaptive_option(f, alg) + +The default `adaptive` setting of `init`: a [`TreeOption`](@ref) in which every node +adapts exactly if its own algorithm is adaptive. Passing a scalar `adaptive` instead +configures the whole tree, including leaves whose algorithm cannot comply. +""" +function default_adaptive_option( + f::GenericSplitFunction, alg::AbstractOperatorSplittingAlgorithm + ) + opt = TreeOption(f, SciMLBase.isadaptive(alg)) + _default_adaptive!(opt, alg) + return opt +end + +function _default_adaptive!(opt::TreeOption, alg) + opt.value = SciMLBase.isadaptive(alg) + inner = alg isa AbstractOperatorSplittingAlgorithm ? alg.inner_algs : () + for (child, inner_alg) in zip(opt.children, inner) + _default_adaptive!(child, inner_alg) + end + return +end + """ warn_non_adaptive(alg, config) @@ -393,8 +417,14 @@ const NODE_OPTION_KEYS = (:dt, :adaptive, :verbose, :controller) inner_values(values::NamedTuple) = NamedTuple{filter(!in(NODE_OPTION_KEYS), keys(values))}(values) -# ... of which a splitting node understands these. -const SPLIT_OPTION_KEYS = (:dtmin, :dtmax, :failfactor, :isoutofdomain) +# ... of which a splitting node understands these. The step-size controller knobs +# (qmin, qmax, gamma, qsteady_min, qsteady_max) are not integrator options: they are +# folded into the default controller of an adaptive node (`default_controller`) and +# travel to the leaves like any other inner option. +const SPLIT_OPTION_KEYS = ( + :dtmin, :dtmax, :failfactor, :isoutofdomain, + :abstol, :reltol, :internalnorm, +) function split_integrator_options(values::NamedTuple) inner = inner_values(values) diff --git a/src/integrator.jl b/src/integrator.jl index 856474b..e6d5499 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -6,13 +6,19 @@ end IntegratorStats() = IntegratorStats(0, 0) -Base.@kwdef mutable struct IntegratorOptions{tType, fType, vbType, F3} +Base.@kwdef mutable struct IntegratorOptions{tType, fType, vbType, F3, atolType, rtolType, normType} adaptive::Bool dtmin::tType = eps(Float64) dtmax::tType = Inf failfactor::fType = 4.0 verbose::vbType = DEFAULT_VERBOSITY isoutofdomain::F3 = DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN + # Error control, used when the node's algorithm provides an error estimate. The + # step size controller knobs (qmin, qmax, gamma, ...) are not held here: they + # live in the node's controller cache, resolved by OrdinaryDiffEqCore. + abstol::atolType = 1.0e-6 + reltol::rtolType = 1.0e-3 + internalnorm::normType = DiffEqBase.ODE_DEFAULT_NORM end @@ -49,9 +55,12 @@ its children's synchronizers, solution indices, and sub-integrators. It does - `t`, `dt`, `dtcache` — time tracking `dtchangeable`, `stops` - `iter` — step counter at this level +- `success_iter` — accepted-step counter at this level - `EEst` — error estimate (`NaN` for non-adaptive, `1.0` default for adaptive) -- `controller` — step-size controller (or `nothing` for non-adaptive) +- `controller_cache` — `OrdinaryDiffEqCore.AbstractControllerCache` holding the + step-size controller and its state (or `nothing` for + non-adaptive) - `force_stepfail` — flag: current step must be retried - `last_step_failed` — flag: previous step failed (double-failure detection) - `status` — [`SplitSubIntegratorStatus`](@ref) for retcode communication @@ -59,8 +68,9 @@ its children's synchronizers, solution indices, and sub-integrators. It does this level - `child_subintegrators` — tuple of direct children (`SplitSubIntegrator` or `DEIntegrator`) -- `solution_indices` — global indices (into parent `u`) **owned by this node** -- `child_solution_indices` — tuple of per-child global solution indices +- `solution_indices` — indices into the *parent's* solution vector **owned by + this node** (indices are parent-relative at every level) +- `child_solution_indices` — tuple of per-child indices into *this node's* `u` - `child_synchronizers` — tuple of per-child synchronizer objects """ mutable struct SplitSubIntegrator{ @@ -88,8 +98,9 @@ mutable struct SplitSubIntegrator{ const dtchangeable::Bool tstops::tstopsType iter::Int + success_iter::Int EEst::EEstType - controller::controllerType + controller_cache::controllerType force_stepfail::Bool last_step_failed::Bool u_modified::Bool # TODO we can probably remove this @@ -174,7 +185,10 @@ mutable struct OperatorSplittingIntegrator{ child_solution_indices::childSolidxType # Tuple child_synchronizers::childSyncType # Tuple iter::Int - controller::controllerType + success_iter::Int + # Step-size controller plus its state, or `nothing` when this node is not adaptive. + controller_cache::controllerType + EEst::tType # error estimate of the last attempted step (NaN when no controller runs) opts::optionsType stats::IntegratorStats tdir::tType @@ -218,9 +232,8 @@ function SciMLBase.__init( save_everystep = false, callback = nothing, advance_to_tstop = false, - adaptive = isadaptive(alg), + adaptive = nothing, controller = nothing, - # controller = OrdinaryDiffEqCore.PIController(0.14, 0.08), alias_u0 = false, verbose = true, kwargs... @@ -228,6 +241,12 @@ function SciMLBase.__init( (; u0, p) = prob t0, tf = prob.tspan + # By default every node adapts exactly if its own algorithm does; a scalar or a + # TreeOption overrides the whole tree explicitly. + if adaptive === nothing + adaptive = default_adaptive_option(prob.f, alg) + end + # Every setting is either one value for the whole tree or a TreeOption carrying a # value per node. Beyond the four this integrator handles itself, whatever the # caller passes travels down to the leaf integrators. @@ -249,12 +268,11 @@ function SciMLBase.__init( tstops = () end - tstops_internal = OrdinaryDiffEqCore.initialize_tstops( - tType, tstops, d_discontinuities, prob.tspan - ) - saveat_internal = OrdinaryDiffEqCore.initialize_saveat(tType, saveat, prob.tspan) - d_discontinuities_internal = OrdinaryDiffEqCore.initialize_d_discontinuities( - tType, d_discontinuities, prob.tspan + # Heaps store raw times and carry the integration direction in their ordering. + # (OrdinaryDiffEqCore's initialize_tstops stores tdir-scaled times instead, which + # is incompatible with the heaps reinit! rebuilds and breaks backward tspans.) + tstops_internal, saveat_internal = tstops_and_saveat_heaps( + t0, tf, (tstops..., d_discontinuities...), saveat ) u = setup_u(prob, alg, alias_u0) @@ -283,6 +301,9 @@ function SciMLBase.__init( child_solution_indices = ntuple(i -> prob.f.solution_indices[i], length(prob.f.functions)) child_synchronizers = ntuple(i -> prob.f.synchronizers[i], length(prob.f.functions)) + root_controller_cache = _node_controller_cache(alg, cache, config.values, tType) + EEst = root_controller_cache === nothing ? tType(NaN) : one(tType) + integrator = OperatorSplittingIntegrator( prob.f, alg, @@ -300,11 +321,12 @@ function SciMLBase.__init( child_subintegrators, child_solution_indices, child_synchronizers, - 0, - config.values.controller, + 0, 0, # iter, success_iter + root_controller_cache, + EEst, split_integrator_options(config.values), IntegratorStats(), - tType(tstops_internal.ordering isa BinaryHeaps.FasterForward ? 1 : -1), + tType(tf > t0 ? 1 : -1), false, config, ) @@ -329,6 +351,11 @@ function DiffEqBase.reinit!( reinit_callbacks = true, reinit_retcode = true ) + # The heap ordering types and every node's tdir are fixed at init, so a reinit! + # cannot flip the integration direction. + (tf > t0) == (integrator.tdir > 0) || + error("reinit! cannot change the direction of integration. Build a new integrator instead.") + # Without a `dt` every node is restored to the step size it was configured with, # so a multi-rate setup survives a reinit!. Passing one is equivalent to passing # it to `init`: a scalar reconfigures the whole tree, a TreeOption node by node. @@ -346,6 +373,8 @@ function DiffEqBase.reinit!( integrator.tstops, integrator.saveat = tstops_and_saveat_heaps(t0, tf, tstops, saveat) integrator.iter = 0 + integrator.success_iter = 0 + reinit_node_controller!(integrator) if erase_sol resize!(integrator.sol.t, 0) resize!(integrator.sol.u, 0) @@ -439,13 +468,11 @@ function _subreinit_child!( SciMLBase.set_proposed_dt!(sub, config.values.dt) set_dt!(sub, config.values.dt) sub.iter = 0 + sub.success_iter = 0 sub.force_stepfail = false sub.last_step_failed = false sub.status = SplitSubIntegratorStatus(ReturnCode.Default) - # Reset EEst to its appropriate default - if isadaptive(sub) - sub.EEst = one(sub.EEst) - end + reinit_node_controller!(sub) # Recurse into this node's children _subreinit_tuple!( f_child, @@ -462,20 +489,22 @@ end # --------------------------------------------------------------------------- function OrdinaryDiffEqCore.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. tdir_t = tdir(integrator) * integrator.t - tdir_tstop = SciMLBase.first_tstop(integrator) + tdir_tstop = tdir(integrator) * SciMLBase.first_tstop(integrator) if tdir_t == tdir_tstop while tdir_t == tdir_tstop SciMLBase.pop_tstop!(integrator) SciMLBase.has_tstop(integrator) ? - (tdir_tstop = SciMLBase.first_tstop(integrator)) : break + (tdir_tstop = tdir(integrator) * SciMLBase.first_tstop(integrator)) : break end notify_integrator_hit_tstop!(integrator) elseif tdir_t > tdir_tstop if !integrator.dtchangeable SciMLBase.change_t_via_interpolation!( integrator, - tdir(integrator) * SciMLBase.pop_tstop!(integrator), + SciMLBase.pop_tstop!(integrator), Val{true} ) notify_integrator_hit_tstop!(integrator) @@ -499,37 +528,40 @@ end # --------------------------------------------------------------------------- function reject_step!(integrator::AnySplitIntegrator) OrdinaryDiffEqCore.increment_reject!(integrator.stats) - return reject_step!(integrator, integrator.cache, integrator.controller) -end -function reject_step!(integrator::AnySplitIntegrator, cache, controller) - integrator.u .= integrator.uprev - rollback_children!(integrator) - return nothing -end -function reject_step!(integrator::AnySplitIntegrator, cache, ::Nothing) if length(integrator.uprev) == 0 error("Cannot roll back integrator. Aborting time integration step at $(integrator.t).") end + integrator.u .= integrator.uprev + rollback_children!(integrator) return nothing end function should_accept_step(integrator::OperatorSplittingIntegrator) - integrator.force_stepfail || integrator.isout && return false - return should_accept_step(integrator, integrator.cache, integrator.controller) + (integrator.force_stepfail || integrator.isout) && return false + return should_accept_step(integrator, integrator.cache, integrator.controller_cache) end function should_accept_step(integrator::SplitSubIntegrator) integrator.force_stepfail && return false - return should_accept_step(integrator, integrator.cache, integrator.controller) + return should_accept_step(integrator, integrator.cache, integrator.controller_cache) end function should_accept_step(integrator::AnySplitIntegrator, cache, ::Nothing) return !(integrator.force_stepfail) end +# An active controller additionally requires the error estimate to pass. +function should_accept_step( + integrator::AnySplitIntegrator, cache, + controller_cache::OrdinaryDiffEqCore.AbstractControllerCache + ) + return accept_step_controller(integrator, controller_cache, integrator.alg) +end +# `stats.naccept` is counted in `step_footer!` (which sees every accepted attempt, +# including the final one); this header-side bookkeeping only prepares the next step. function accept_step!(integrator::AnySplitIntegrator) - OrdinaryDiffEqCore.increment_accept!(integrator.stats) - return accept_step!(integrator, integrator.cache, integrator.controller) + integrator.success_iter += 1 + return accept_step!(integrator, integrator.cache, integrator.controller_cache) end -function accept_step!(integrator::AnySplitIntegrator, cache, controller) +function accept_step!(integrator::AnySplitIntegrator, cache, controller_cache) return store_previous_info!(integrator) end function store_previous_info!(integrator::AnySplitIntegrator) @@ -543,23 +575,67 @@ function update_uprev!(integrator::AnySplitIntegrator) return nothing end -# Roll back each child's local buffer to match master u. -# For DEIntegrators the leaf will be re-synced via forward_sync before the -# next attempt, so there is nothing to do here. -rollback_children!(integrator::OperatorSplittingIntegrator) = rollback_children!(integrator.child_subintegrators, integrator.u) -@unroll function rollback_children!(children::Tuple, u_master) +# Roll a node's children back to the state of the node itself: the local solution +# buffers are refilled from the parent's (already restored) `u` and the child clocks +# are moved back to the parent's time. Leaves get their `u` restored here as well, +# not only by the forward sync before the next solve: palindromic algorithms skip +# that sync for the first child (`next_sync_is_continuous`), so a rollback that left +# the leaf state stale would silently resume from the failed attempt. +rollback_children!(parent::AnySplitIntegrator) = _rollback_children!( + parent.child_subintegrators, parent.child_solution_indices, parent.u, parent.t +) +@unroll function _rollback_children!(children::Tuple, solution_indices::Tuple, u_parent, t) + i = 0 @unroll for child in children - rollback_child!(child, u_master) + i += 1 + rollback_child!(child, u_parent, solution_indices[i], t) end end -function rollback_child!(child::SplitSubIntegrator, u_master) - child.u .= @view u_master[child.solution_indices] +function rollback_child!(child::SplitSubIntegrator, u_parent, solution_indices, t) + child.u .= @view u_parent[solution_indices] + RecursiveArrayTools.recursivecopy!(child.uprev, child.u) + child.t = t + child.tprev = t + _reset_child_failure!(child) + rollback_children!(child) + return nothing +end +function rollback_child!(child::DEIntegrator, u_parent, solution_indices, t) + child.u .= @view u_parent[solution_indices] RecursiveArrayTools.recursivecopy!(child.uprev, child.u) - _rollback_children!(child.child_subintegrators, u_master) + mark_state_modified!(child) + child.t = t + child.tprev = t + _reset_child_failure!(child) return nothing end -function rollback_child!(child::DEIntegrator, u_master) - # forward_sync before the next sub-step will restore this correctly. + +# Clear a failed child's retcode so the retry of the nearest adaptive ancestor can +# 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.last_step_failed = false + child.force_stepfail = false + return nothing +end +function _reset_child_failure!(child::DEIntegrator) + # Based on the *stored* retcode (the sticky part), not on `check_error`: the + # rollback just restored this child's state, so a state-based check would + # already come back clean while the stored retcode still blocks re-stepping. + if child.sol.retcode ∉ (ReturnCode.Default, ReturnCode.Success) + child.sol = SciMLBase.solution_new_retcode(child.sol, ReturnCode.Default) + end + # Resetting the retcode is not sufficient: SciMLBase's generic check_error + # re-derives ConvergenceFailure from a sticky `last_stepfail` on non-adaptive + # leaves (e.g. after an inner Newton failure), which would turn the retry the + # escalation protocol just set up into an immediate failure again. + if hasfield(typeof(child), :last_stepfail) + child.last_stepfail = false + end + if hasfield(typeof(child), :force_stepfail) + child.force_stepfail = false + end return nothing end @@ -586,7 +662,7 @@ end function modify_dt_for_tstops!(integrator) if SciMLBase.has_tstop(integrator) tdir_t = integrator.tdir * integrator.t - tdir_tstop = SciMLBase.first_tstop(integrator) + tdir_tstop = integrator.tdir * SciMLBase.first_tstop(integrator) if integrator.opts.adaptive integrator.dt = integrator.tdir * min(abs(integrator.dt), abs(tdir_tstop - tdir_t)) # step! to the end @@ -625,12 +701,17 @@ function fix_solution_buffer_sizes!(integrator, sol) return end +# Window for absorbing floating point drift when landing on a time point. Scaled +# by the local time scale *and* the step size: near t = 0 (e.g. integrating +# backward to zero) a purely value-relative window collapses below the ulp drift +# the subdivided child steps accumulate. +_snap_window(t, tstop, dt) = + 100 * eps(float(max(abs(t), abs(tstop), abs(dt)) / oneunit(t))) * oneunit(t) + function fixed_t_for_floatingpoint_error!(integrator::AnySplitIntegrator, ttmp) return if DiffEqBase.has_tstop(integrator) - tstop = integrator.tdir * DiffEqBase.first_tstop(integrator) - if abs(ttmp - tstop) < - 100eps(float(max(integrator.t, tstop) / oneunit(integrator.t))) * - oneunit(integrator.t) + tstop = DiffEqBase.first_tstop(integrator) + if abs(ttmp - tstop) < _snap_window(integrator.t, tstop, integrator.dt) try_snap_children_to_tstop!.(integrator.child_subintegrators, tstop) tstop else @@ -641,8 +722,7 @@ function fixed_t_for_floatingpoint_error!(integrator::AnySplitIntegrator, ttmp) end end function try_snap_children_to_tstop!(integrator::SplitSubIntegrator, tstop) - if abs(tstop - integrator.t) < - 100eps(float(max(integrator.t, tstop) / oneunit(integrator.t))) * oneunit(integrator.t) + if abs(tstop - integrator.t) < _snap_window(integrator.t, tstop, integrator.dt) integrator.t = tstop else @warn "Failed to snap timestep for integrator $(integrator.t) with parent integrator hitting the tstop $(tstop)." @@ -650,8 +730,7 @@ function try_snap_children_to_tstop!(integrator::SplitSubIntegrator, tstop) return try_snap_children_to_tstop!.(integrator.child_subintegrators, tstop) end function try_snap_children_to_tstop!(integrator::DEIntegrator, tstop) - return if abs(tstop - integrator.t) < - 100eps(float(max(integrator.t, tstop) / oneunit(integrator.t))) * oneunit(integrator.t) + return if abs(tstop - integrator.t) < _snap_window(integrator.t, tstop, integrator.dt) integrator.t = tstop else @warn "Failed to snap timestep for integrator $(integrator.t) with parent integrator hitting the tstop $(tstop)." @@ -659,10 +738,11 @@ function try_snap_children_to_tstop!(integrator::DEIntegrator, tstop) end function step_footer!(integrator::AnySplitIntegrator) - ttmp = integrator.t + tdir(integrator) * integrator.dt + ttmp = integrator.t + integrator.dt # dt is signed by the integration direction footer_reset_flags!(integrator) setup_validity_flags!(integrator, ttmp) if should_accept_step(integrator) + OrdinaryDiffEqCore.increment_accept!(integrator.stats) integrator.last_step_failed = false integrator.tprev = integrator.t integrator.t = fixed_t_for_floatingpoint_error!(integrator, ttmp) @@ -674,19 +754,52 @@ function step_footer!(integrator::AnySplitIntegrator) step_accept_controller!(integrator) validate_time_point(integrator) elseif integrator.force_stepfail - if isadaptive(integrator) - step_reject_controller!(integrator) + # Failure escalation protocol: the failing node's own adaptivity decides. + fatal_rc = _fatal_child_retcode(integrator.child_subintegrators) + if fatal_rc !== ReturnCode.Default + # An *adaptive* child failed: it already exhausted its own step size + # adaptation, so retrying on a smaller interval cannot help. Propagate + # its diagnosis and stop. + _set_retcode!(integrator, fatal_rc) + elseif integrator.controller_cache !== nothing + # A non-adaptive descendant failed and this is the nearest adaptive + # ancestor: retry the step on a failfactor-shrunken interval, which + # shrinks the effective dt of the whole subtree. The failed inner solve + # tells us nothing about the error, so no step size law here. The + # header-side reject_step! resets the subtree, including retcodes. OrdinaryDiffEqCore.post_newton_controller!(integrator, integrator.alg) - elseif integrator.dtchangeable - integrator.dt /= integrator.opts.failfactor - elseif integrator.last_step_failed - return + integrator.dtcache = abs(integrator.dt) + abort_below_dtmin!(integrator) + else + # Non-adaptive node: escalate the failure to this node's parent. At a + # non-adaptive root this stops the time integration. + _set_retcode!( + integrator, + _first_failed_child_retcode(integrator.child_subintegrators) + ) end integrator.last_step_failed = true + else + # The controller rejected the step (EEst > 1): shrink dt and retry. + step_reject_controller!(integrator) + abort_below_dtmin!(integrator) + integrator.last_step_failed = true end return nothing end +function abort_below_dtmin!(integrator::AnySplitIntegrator) + abs(integrator.dt) > abs(OrdinaryDiffEqCore.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) + return nothing +end + +_set_retcode!(integrator::OperatorSplittingIntegrator, code) = + integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, code) +_set_retcode!(integrator::SplitSubIntegrator, code) = integrator.status.retcode = code + # --------------------------------------------------------------------------- # __solve / solve! / step! # --------------------------------------------------------------------------- @@ -700,7 +813,8 @@ end function DiffEqBase.solve!(integrator::OperatorSplittingIntegrator) while !isempty(integrator.tstops) - while tdir(integrator) * integrator.t < SciMLBase.first_tstop(integrator) + while tdir(integrator) * integrator.t < + tdir(integrator) * SciMLBase.first_tstop(integrator) step_header!(integrator) @timeit_debug "check_error" SciMLBase.check_error!(integrator) ∉ ( ReturnCode.Success, ReturnCode.Default, @@ -750,12 +864,14 @@ function DiffEqBase.step!(integrator::AnySplitIntegrator) return end +# SciML convention: `dt` is signed and its sign has to match the direction of +# integration. function DiffEqBase.step!(integrator::AnySplitIntegrator, dt, stop_at_tdt = false) @timeit_debug "step!" begin - dt <= zero(dt) && error("dt must be positive") + tdir(integrator) * dt < zero(dt) && error("Cannot step backward.") stop_at_tdt && !integrator.dtchangeable && error("Cannot stop at t + dt if dtchangeable is false") - tnext = integrator.t + tdir(integrator) * dt + tnext = integrator.t + dt stop_at_tdt && DiffEqBase.add_tstop!(integrator, tnext) while !reached_tstop(integrator, tnext, stop_at_tdt) step_header!(integrator) @@ -764,6 +880,13 @@ function DiffEqBase.step!(integrator::AnySplitIntegrator, dt, stop_at_tdt = fals ) && return __step!(integrator) step_footer!(integrator) + # Pop every tstop as soon as it is reached, exactly like the solve! + # loop does. Intermediate stops before `tnext` do occur -- the stale + # `tnext` of an earlier failed attempt, or stops pushed down from the + # 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) end end OrdinaryDiffEqCore.handle_tstop!(integrator) @@ -835,35 +958,134 @@ function (integrator::OperatorSplittingIntegrator)(tmp, t) ) end -# Stepsize controller hooks -@inline function stepsize_controller!(integrator::AnySplitIntegrator) - isadaptive(integrator.alg) || return nothing - stepsize_controller!(integrator, integrator.alg) - return nothing -end -@inline function stepsize_controller!(integrator::AnySplitIntegrator, alg::AbstractOperatorSplittingAlgorithm) - isadaptive(alg) || return nothing - #stepsize_controller!(integrator, integrator.controller) - return nothing -end -@inline function step_accept_controller!(integrator::AnySplitIntegrator) - isadaptive(integrator.alg) || return nothing - step_accept_controller!(integrator, integrator.alg) +# --------------------------------------------------------------------------- +# 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!`); +# 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 +# `stepsize_controller!`/`step_accept_controller!`/`step_reject_controller!` +# thread between the steps (e.g. `errold` of a PIController). Following the +# OrdinaryDiffEq conventions the step is accepted if `EEst <= 1`. +# --------------------------------------------------------------------------- + +const CONTROLLER_KNOB_KEYS = (:qmin, :qmax, :gamma, :qsteady_min, :qsteady_max, :failfactor) + +""" + default_controller(alg::AbstractOperatorSplittingAlgorithm, values::NamedTuple) + +The step size controller an adaptive splitting node runs with when `init` is not +given an explicit `controller`. The controller knobs the caller passed to `init` +(`qmin`, `qmax`, `gamma`, `qsteady_min`, `qsteady_max`, `failfactor`) ride along +as overrides; unset ones resolve to the algorithm defaults in +`setup_controller_cache`. +""" +default_controller(::AbstractOperatorSplittingAlgorithm, values::NamedTuple) = + OrdinaryDiffEqCore.IController( + NamedTuple{filter(in(CONTROLLER_KNOB_KEYS), keys(values))}(values) +) + +function _node_controller(alg::AbstractOperatorSplittingAlgorithm, values::NamedTuple) + (values.adaptive && SciMLBase.isadaptive(alg)) || return nothing + values.controller === nothing || return values.controller + return default_controller(alg, values) +end + +function _node_controller_cache( + alg::AbstractOperatorSplittingAlgorithm, level_cache, + values::NamedTuple, ::Type{tType} + ) where {tType} + controller = _node_controller(alg, values) + controller === nothing && return nothing + if _wants_discontinuity_detection(controller) + throw( + ArgumentError( + "discontinuity_detection is not supported by operator splitting nodes. \ + Construct the controller for this node without it." + ) + ) + end + return OrdinaryDiffEqCore.setup_controller_cache(alg, level_cache, controller, tType) +end + +# The discontinuity handling of OrdinaryDiffEqCore's controllers needs integrator +# state (callbacks, checkpoints) a splitting node does not have, so refuse it +# up front instead of failing deep inside a rejected step. +_wants_discontinuity_detection(controller) = + hasfield(typeof(controller), :basic) && _basic_discontinuity_detection(controller.basic) +_basic_discontinuity_detection(basic::NamedTuple) = + get(basic, :discontinuity_detection, false) === true +_basic_discontinuity_detection(basic) = + hasfield(typeof(basic), :discontinuity_detection) && basic.discontinuity_detection + +reinit_node_controller!(integrator::AnySplitIntegrator) = + reinit_node_controller!(integrator, integrator.controller_cache) +function reinit_node_controller!(integrator::AnySplitIntegrator, ::Nothing) + integrator.EEst = oftype(integrator.EEst, NaN) return nothing end -@inline function step_accept_controller!(integrator::AnySplitIntegrator, alg::AbstractOperatorSplittingAlgorithm) - isadaptive(alg) || return nothing - #step_accept_controller!(integrator, integrator.controller) +function reinit_node_controller!( + integrator::AnySplitIntegrator, + controller_cache::OrdinaryDiffEqCore.AbstractControllerCache + ) + integrator.EEst = one(integrator.EEst) + OrdinaryDiffEqCore.reinit_controller!(integrator, controller_cache) return nothing end -@inline function step_reject_controller!(integrator::AnySplitIntegrator) - isadaptive(integrator.alg) || return nothing - step_reject_controller!(integrator, integrator.alg) + +""" + alg_adaptive_order(alg::AbstractOperatorSplittingAlgorithm) + +Order of the error estimator of an adaptive operator splitting algorithm; every +algorithm with `SciMLBase.isadaptive(alg) == true` has to implement this. +""" +function alg_adaptive_order end + +# The controller caches consume the error estimate and the estimator order through +# these OrdinaryDiffEqCore hooks. Our nodes keep `EEst` on the integrator itself. +@inline OrdinaryDiffEqCore.get_EEst(integrator::AnySplitIntegrator) = integrator.EEst +@inline OrdinaryDiffEqCore.set_EEst!(integrator::AnySplitIntegrator, val) = + integrator.EEst = oftype(integrator.EEst, val) +OrdinaryDiffEqCore.get_current_adaptive_order( + alg::AbstractOperatorSplittingAlgorithm, cache +) = alg_adaptive_order(alg) +# The generic fallbacks of these two are meant for inner solvers with their own +# tuning (gamma_default falls back to 0, which would ruin the step size law). +OrdinaryDiffEqCore.gamma_default(::AbstractOperatorSplittingAlgorithm) = 9 // 10 +OrdinaryDiffEqCore.failfactor_default(::AbstractOperatorSplittingAlgorithm) = 4 + +@inline step_accept_controller!(integrator::AnySplitIntegrator) = + step_accept_controller!(integrator, integrator.controller_cache) +step_accept_controller!(integrator::AnySplitIntegrator, ::Nothing) = nothing +function step_accept_controller!( + integrator::AnySplitIntegrator, + controller_cache::OrdinaryDiffEqCore.AbstractControllerCache + ) + q = stepsize_controller!(integrator, controller_cache, integrator.alg) + dtnew = step_accept_controller!(integrator, controller_cache, integrator.alg, q) + # The proposal derives from the step actually taken -- which + # `modify_dt_for_tstops!` may have clipped to a tstop. After such a step the + # proposal is rebased and regrows at up to qmax per step; that matches what + # the error-based law knows (EEst belongs to the clipped step) and mirrors + # OrdinaryDiffEq. `dtcache` mirrors the standing proposal for `reinit!` and + # introspection. + integrator.dt = dtnew + integrator.dtcache = abs(dtnew) return nothing end -@inline function step_reject_controller!(integrator::AnySplitIntegrator, alg::AbstractOperatorSplittingAlgorithm) - isadaptive(integrator.alg) || return nothing - # step_reject_controller!(integrator, integrator.controller) + +@inline step_reject_controller!(integrator::AnySplitIntegrator) = + step_reject_controller!(integrator, integrator.controller_cache) +step_reject_controller!(integrator::AnySplitIntegrator, ::Nothing) = nothing +function step_reject_controller!( + integrator::AnySplitIntegrator, + controller_cache::OrdinaryDiffEqCore.AbstractControllerCache + ) + stepsize_controller!(integrator, controller_cache, integrator.alg) + step_reject_controller!(integrator, controller_cache, integrator.alg) # sets dt + integrator.dtcache = abs(integrator.dt) return nothing end @@ -875,7 +1097,7 @@ is_past_t(integrator, t) = tdir(integrator) * (t - integrator.t) ≤ zero(integrator.t) function reached_tstop(integrator, tstop, stop_at_tstop = integrator.dtchangeable) if stop_at_tstop - integrator.t > tstop && + tdir(integrator) * (integrator.t - tstop) > zero(integrator.t) && error("Integrator missed stop at $tstop (current time=$(integrator.t)). Aborting.") return integrator.t ≈ tstop else @@ -899,7 +1121,6 @@ end function __step!(integrator::AnySplitIntegrator) advance_solution_by!(integrator, integrator.dt) - stepsize_controller!(integrator) # FIXME this should go into the footer return nothing end @@ -932,20 +1153,33 @@ function advance_solution_by!( cache::AbstractOperatorSplittingCache, dt ) + # Success and failure are both handled in step_footer! via the failure + # escalation protocol; nothing to decide here. _perform_step!(outer, children, cache, dt) + return +end - if outer.force_stepfail && all(isadaptive.(children)) - # We do not know recover at this point, as an decrease in the solve - # interval is unlikely to help here. - outer.sol = SciMLBase.solution_new_retcode( - outer.sol, - ReturnCode.Failure - ) - return +# Retcode of the first failed child whose failure is fatal (the child is adaptive, +# 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) + return _failure_retcode(child) + end end + return ReturnCode.Default +end - return +@unroll function _first_failed_child_retcode(children::Tuple) + @unroll for child in children + _child_failed(child) && return _failure_retcode(child) + end + return ReturnCode.Failure end +_failure_retcode(child::DEIntegrator) = SciMLBase.check_error(child) +_failure_retcode(child::SplitSubIntegrator) = child.status.retcode +_child_is_adaptive(child::DEIntegrator) = child.opts.adaptive +_child_is_adaptive(child::SplitSubIntegrator) = child.controller_cache !== nothing function advance_solution_by!( outer::SplitSubIntegrator, @@ -955,17 +1189,21 @@ function advance_solution_by!( ) _perform_step!(outer, children, cache, dt) - if outer.force_stepfail - outer.status = SplitSubIntegratorStatus(ReturnCode.Failure) - return + # On force_stepfail the status is left clean: step_footer! either retries at + # this level (adaptive) or escalates by writing the failure into the status. + if !outer.force_stepfail + outer.status.retcode = ReturnCode.Success end - # All children succeeded: advance this node's time and counter - outer.status = SplitSubIntegratorStatus(ReturnCode.Success) - return end +# `dt` stays signed through the splitting tree, following the SciML `step!` +# convention: the sign has to match the child's own integration direction, which +# equals the tree's. Advancing a child *against* its direction (negative substeps +# of higher order compositions) is not supported yet: leaf ODEIntegrators cannot +# step against the tdir their tspan fixed at construction. + # Recursion dispatch function advance_solution_by!( outer::AnySplitIntegrator, @@ -1038,14 +1276,21 @@ function _build_child( dt = config.values.dt tType = typeof(dt) - # Recurse: build each consecutive child + u_sub = RecursiveArrayTools.recursivecopy(uouter[solution_indices]) + uprev_sub = RecursiveArrayTools.recursivecopy(uprevouter[solution_indices]) + + # Recurse: build each consecutive child. Solution indices are relative to the + # *parent* at every level, so the children have to address this node's buffers, + # not the outer ones: a child that wires itself as a view into the handed + # buffer (instead of copying, as the stock leaves do) would otherwise alias + # the wrong slots of the root vector. child_subintegrators = ntuple( i -> _build_child( prob, alg.inner_algs[i], get_operator(f, i), p[i], - uprevouter, uouter, u_master, + uprev_sub, u_sub, u_master, f.solution_indices[i], t0, tf, tstops, saveat, d_discontinuities, callback, @@ -1057,11 +1302,8 @@ function _build_child( child_solution_indices = ntuple(i -> f.solution_indices[i], length(f.functions)) child_synchronizers = ntuple(i -> f.synchronizers[i], length(f.functions)) - u_sub = RecursiveArrayTools.recursivecopy(uouter[solution_indices]) - uprev_sub = RecursiveArrayTools.recursivecopy(uprevouter[solution_indices]) - - tstops_internal = OrdinaryDiffEqCore.initialize_tstops( - tType, tstops, d_discontinuities, prob.tspan + tstops_internal, _ = tstops_and_saveat_heaps( + t0, tf, (tstops..., d_discontinuities...), () ) level_cache = init_cache( @@ -1069,7 +1311,8 @@ function _build_child( uprev = uprev_sub, u = u_sub, ) - EEst_val = isadaptive(alg) ? one(tType) : tType(NaN) + controller_cache = _node_controller_cache(alg, level_cache, config.values, tType) + EEst_val = controller_cache === nothing ? tType(NaN) : one(tType) sub = SplitSubIntegrator( alg, @@ -1079,9 +1322,9 @@ function _build_child( t0, t0, dt, dt, # t, tprev, dt, dtcache isdtchangeable(alg), tstops_internal, - 0, # iter + 0, 0, # iter, success_iter EEst_val, - config.values.controller, + controller_cache, false, false, false, # force_stepfail, last_step_failed, u_modified SplitSubIntegratorStatus(), IntegratorStats(), @@ -1091,7 +1334,7 @@ function _build_child( child_solution_indices, child_synchronizers, split_integrator_options(config.values), - one(tType), + sign(dt), # dt was signed by signed_dt_tree, so its sign is the direction ) return sub @@ -1170,7 +1413,7 @@ SciMLBase.pop_tstop!(i::AnySplitIntegrator) = pop!(i.tstops) DiffEqBase.get_dt(i::AnySplitIntegrator) = i.dt function set_dt!(i::DEIntegrator, dt) - dt <= zero(dt) && error("dt must be positive") + iszero(dt) && error("dt must be nonzero") return i.dt = dt end diff --git a/src/precompilation.jl b/src/precompilation.jl index 4b1072e..a966ca1 100644 --- a/src/precompilation.jl +++ b/src/precompilation.jl @@ -48,4 +48,10 @@ end integrator_sm = DiffEqBase.init(prob_sm, tstepper_sm, dt = 0.01, verbose = false) step!(integrator_sm) solve!(integrator_sm) + + # Precompile PalindromicPairLieTrotterGodunov (runs its controller by default) + tstepper_pp = PalindromicPairLieTrotterGodunov((Euler(), Euler())) + integrator_pp = DiffEqBase.init(prob_sm, tstepper_pp, dt = 0.01, verbose = false) + step!(integrator_pp) + solve!(integrator_pp) end diff --git a/src/solver.jl b/src/solver.jl index 2718b40..156313a 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -160,3 +160,129 @@ function _perform_step!( return end + +# --------------------------------------------------------------------------- +# Palindromic pair of Lie-Trotter-Godunov splittings +# --------------------------------------------------------------------------- +""" + PalindromicPairLieTrotterGodunov <: AbstractOperatorSplittingAlgorithm + +Second-order sequential operator splitting algorithm. + +One step solves the palindromic pair of [`LieTrotterGodunov`](@ref) sequences + +``A_1(\\Delta t) \\to \\cdots \\to A_N(\\Delta t)`` and +``A_N(\\Delta t) \\to \\cdots \\to A_1(\\Delta t)`` + +from the same initial value. The leading splitting error of a Lie-Trotter sequence +is ``\\frac{\\Delta t^2}{2}\\sum_{i ") + end + length(alg.inner_algs) > 0 && Base.show(io, alg.inner_algs[end]) + return print(io, ")") +end + +@inline SciMLBase.isadaptive(::PalindromicPairLieTrotterGodunov) = true +# The pair difference estimates the O(dt²) leading error term of a first order +# sequence, so the controller sees a first order error estimator. +alg_adaptive_order(::PalindromicPairLieTrotterGodunov) = 1 + +struct PalindromicPairLieTrotterGodunovCache{uType, uprevType, uforwardType} <: AbstractOperatorSplittingCache + u::uType + uprev::uprevType + uforward::uforwardType # end state of the A₁ → A₂ sequence; reused as the residual buffer +end + +function init_cache( + f::GenericSplitFunction, alg::PalindromicPairLieTrotterGodunov; + uprev::AbstractArray, u::AbstractVector, + ) + return PalindromicPairLieTrotterGodunovCache(u, uprev, similar(u)) +end + +function _ppltg_advance_child!(parent, child, i, dt) + idxs = parent.child_solution_indices[i] + sync = parent.child_synchronizers[i] + + @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync) + @timeit_debug "time solve" advance_solution_by!(parent, child, dt) + if _child_failed(child) + parent.force_stepfail = true + return + end + + @timeit_debug "sync <-" backward_sync_subintegrator!(parent, child, idxs, sync) + return +end + +# Forward sequence: A₁(dt) → … → A_N(dt) +@unroll function _ppltg_forward_pass!(parent, children::Tuple, dt) + i = 0 + @unroll for child in children + i += 1 + _ppltg_advance_child!(parent, child, i, dt) + parent.force_stepfail && return + end +end + +# Reverse sequence: A_N(dt) → … → A₁(dt) +@unroll function _ppltg_reverse_pass!(parent, rchildren::Tuple, dt, N) + j = 0 + @unroll for child in rchildren + j += 1 + _ppltg_advance_child!(parent, child, N + 1 - j, dt) + parent.force_stepfail && return + end +end + +function _perform_step!( + parent, + children::Tuple, + cache::PalindromicPairLieTrotterGodunovCache, + dt + ) + (; uforward) = cache + + _ppltg_forward_pass!(parent, children, dt) + parent.force_stepfail && return + uforward .= parent.u + + # Rewind to the initial state of the step; uprev is untouched while stepping. + parent.u .= parent.uprev + rollback_children!(parent) + + _ppltg_reverse_pass!(parent, reverse(children), dt, length(children)) + parent.force_stepfail && return + + # The average of the pair is the second order solution ... + parent.u .= (parent.u .+ uforward) ./ 2 + if parent.controller_cache !== nothing + # ... and half the pair difference the local error of a single sequence. + (; 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) + end + return +end diff --git a/src/utils.jl b/src/utils.jl index cd78bff..642fd79 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -69,20 +69,15 @@ function forward_sync_subintegrator!( reset_next_sync_continuous(parent) return nothing end - forward_sync_internal!(parent.u, child, solution_indices) + forward_sync_internal!(parent.u, parent.uprev, child, solution_indices) @timeit_debug "external sync" forward_sync_external!(parent, child, sync) return nothing end -# Shared internal helper: copy master u slice → leaf DEIntegrator u/uprev -function forward_sync_internal!(u_source, child::DEIntegrator, solution_indices) - @views usrc = u_source[solution_indices] - @timeit_debug "sync vectors" begin - sync_vectors!(child.u, usrc) - sync_vectors!(child.uprev, child.u) - end - # SciMLBase v3 renamed this to `derivative_discontinuity!`; call the - # appropriate name based on which SciMLBase is loaded. +# Tell a leaf integrator that its state was changed from the outside so it discards +# FSAL information. SciMLBase v3 renamed `u_modified!` → `derivative_discontinuity!`; +# call the appropriate name based on which SciMLBase is loaded. +function mark_state_modified!(child::DEIntegrator) @static if isdefined(SciMLBase, :derivative_discontinuity!) SciMLBase.derivative_discontinuity!(child, true) else @@ -91,6 +86,24 @@ function forward_sync_internal!(u_source, child::DEIntegrator, solution_indices) return nothing end +# Shared internal helper: copy the parent u slice → child DEIntegrator u/uprev. +# `uprev_parent` is the parent's rollback buffer: the refresh of the child's own +# rollback anchor must never write through a child `uprev` that aliases it, because +# that buffer has to survive the whole step untouched (rejection and the palindromic +# mid-step rewind restore from it). When they alias, the slice already holds the +# interval start state and no copy is needed. +function forward_sync_internal!(u_source, uprev_parent, child::DEIntegrator, solution_indices) + @views usrc = u_source[solution_indices] + @timeit_debug "sync vectors" begin + sync_vectors!(child.u, usrc) + if need_sync(child.uprev, uprev_parent) + sync_vectors!(child.uprev, child.u) + end + end + mark_state_modified!(child) + return nothing +end + """ backward_sync_subintegrator!(parent_integrator::OperatorSplittingIntegrator, inner_integrator::DEIntegrator, solution_indices, sync) @@ -162,17 +175,11 @@ end # Time stuff function OrdinaryDiffEqCore.fix_dt_at_bounds!(integrator::AnySplitIntegrator) - if tdir(integrator) > 0 - integrator.dt = min(integrator.opts.dtmax, integrator.dt) - else - integrator.dt = max(integrator.opts.dtmax, integrator.dt) - end - dtmin = OrdinaryDiffEqCore.timedepentdtmin(integrator) - if tdir(integrator) > 0 - integrator.dt = max(integrator.dt, dtmin) - else - integrator.dt = min(integrator.dt, dtmin) - end + # 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)) + integrator.dt = tdir(integrator) * max(min(abs(integrator.dt), dtmax), dtmin) return nothing end @@ -194,10 +201,15 @@ function validate_time_point(parent, child::DEIntegrator) end # --------------------------------------------------------------------------- -# _child_failed: check whether a child reported a failure +# _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.sol.retcode ∉ (ReturnCode.Default, ReturnCode.Success) + SciMLBase.check_error(child) ∉ (ReturnCode.Default, ReturnCode.Success) _child_failed(child::SplitSubIntegrator) = child.status.retcode ∉ (ReturnCode.Default, ReturnCode.Success) diff --git a/test/adaptivity.jl b/test/adaptivity.jl new file mode 100644 index 0000000..4183046 --- /dev/null +++ b/test/adaptivity.jl @@ -0,0 +1,365 @@ +using OrdinaryDiffEqOperatorSplitting +import OrdinaryDiffEqOperatorSplitting as OS +using Test + +import DiffEqBase: DiffEqBase, ODEFunction +import SciMLBase +import SciMLBase: ReturnCode +import OrdinaryDiffEqCore +using OrdinaryDiffEqLowOrderRK +using OrdinaryDiffEqTsit5 + +# Non-commuting pair of linear operators, [A, B] ≠ 0, so the splitting error is +# nonzero and the palindromic pair's error estimator has something real to measure. +A = [-1.0 0.0; 0.0 -2.0] +B = [0.0 0.5; 0.5 0.0] +odeA(du, u, p, t) = (du[1] = -u[1]; du[2] = -2 * u[2]; nothing) +odeB(du, u, p, t) = (du[1] = 0.5 * u[2]; du[2] = 0.5 * u[1]; nothing) +fA = ODEFunction(odeA) +fB = ODEFunction(odeB) + +dofs = [1, 2] +u0 = [1.0, 1.0] +tspan = (0.0, 1.0) +fsplit = GenericSplitFunction((fA, fB), (dofs, dofs)) +prob = OperatorSplittingProblem(fsplit, u0, tspan) +trueu = exp(tspan[2] * (A + B)) * u0 + +PPLTG = PalindromicPairLieTrotterGodunov + +@testset "PPLTG adaptivity" begin + @testset "Error estimate of a single step" begin + # With Euler sub-solvers taking a single step per pass the whole pair is + # reproducible by hand: forward sequence uA → uAB, backward uB → uBA, the + # solution is the pair average, and EEst the tolerance-scaled half difference. + dt = 0.05 + abstol, reltol = 1.0e-6, 1.0e-3 + integ = DiffEqBase.init( + prob, PPLTG((Euler(), Euler())); dt, abstol, reltol + ) + + uA = u0 .+ dt .* (A * u0) + uAB = uA .+ dt .* (B * uA) + uB = u0 .+ dt .* (B * u0) + uBA = uB .+ dt .* (A * uB) + u_expected = (uAB .+ uBA) ./ 2 + resid = (uBA .- uAB) ./ 2 ./ + (abstol .+ max.(abs.(u_expected), abs.(u0)) .* reltol) + EEst_expected = sqrt(sum(abs2, resid) / length(resid)) # ODE_DEFAULT_NORM + + DiffEqBase.step!(integ) + @test integ.u ≈ u_expected + @test integ.EEst ≈ EEst_expected + @test integ.t ≈ dt + + # The accepted step drives the I controller with exponent + # 1/(alg_adaptive_order + 1) = 1/2 and gamma = 9/10; the first accepted step + # may grow by up to qmax_first_step = 10^4. The rtol absorbs the fastpower + # approximation upstream while still discriminating a wrong exponent (~10%). + q_expected = clamp(sqrt(EEst_expected) / (9 / 10), 1 / 10^4, 5) + @test integ.dt ≈ dt / q_expected rtol = 1.0e-3 + @test integ.dtcache ≈ dt / q_expected rtol = 1.0e-3 + end + + @testset "Error estimate scales with dt²" begin + # For single Euler passes over linear operators the pair difference is + # exactly dt²[A,B]u₀ and u₀ dominates the residual scaling, so halving dt + # divides EEst by exactly four. + EEsts = map((0.02, 0.01)) do dt + integ = DiffEqBase.init(prob, PPLTG((Euler(), Euler())); dt) + DiffEqBase.step!(integ) + integ.EEst + end + @test EEsts[1] / EEsts[2] ≈ 4 rtol = 1.0e-6 + end + + @testset "Controller expands the step size under loose tolerances" begin + integ = DiffEqBase.init( + prob, PPLTG((Tsit5(), Tsit5())); + dt = 1.0e-3, reltol = 1.0e-2, abstol = 1.0e-4 + ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + # A fixed-step run would need 1000 steps; the controller settles near the + # equilibrium step size within a handful of them. + @test integ.stats.naccept < 100 + @test integ.dtcache > 1.0e-2 + @test maximum(abs, integ.u .- trueu) < 0.05 + end + + @testset "Commuting operators: steps grow to the tolerance-free maximum" begin + # With commuting operators the splitting is exact, EEst ≈ 0, and the + # controller must grow dt by qmax every step until the tstop truncates it. + odeC1(du, u, p, t) = (du .= -u; nothing) + odeC2(du, u, p, t) = (du .= -0.5 .* u; nothing) + fsplit_c = GenericSplitFunction( + (ODEFunction(odeC1), ODEFunction(odeC2)), (dofs, dofs) + ) + prob_c = OperatorSplittingProblem(fsplit_c, u0, tspan) + integ = DiffEqBase.init(prob_c, PPLTG((Tsit5(), Tsit5())); dt = 1.0e-3) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.stats.naccept ≤ 6 + @test maximum(abs, integ.u .- exp(-1.5) .* u0) < 1.0e-3 + end + + @testset "Rejected steps roll back cleanly" begin + integ = DiffEqBase.init( + prob, PPLTG((Tsit5(), Tsit5())); + dt = 0.5, reltol = 1.0e-6, abstol = 1.0e-8 + ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.stats.nreject ≥ 1 + @test maximum(abs, integ.u .- trueu) < 1.0e-3 + end + + @testset "Rejections with time-dependent operators" begin + # If a rejected step failed to rewind the child integrator clocks, the retry + # would integrate the wrong time interval, which a time-dependent operator + # turns into a visible error against the reference solve. + odeAt(du, u, p, t) = (c = 1 + sin(2π * t); du[1] = -c * u[1]; du[2] = -2c * u[2]; nothing) + fsplit_t = GenericSplitFunction((ODEFunction(odeAt), fB), (dofs, dofs)) + prob_t = OperatorSplittingProblem(fsplit_t, u0, tspan) + + ref = DiffEqBase.init( + prob_t, PPLTG((Tsit5(), Tsit5())); dt = 1.0e-3, adaptive = false + ) + DiffEqBase.solve!(ref) + + integ = DiffEqBase.init( + prob_t, PPLTG((Tsit5(), Tsit5())); + dt = 0.4, reltol = 1.0e-5, abstol = 1.0e-8 + ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.stats.nreject ≥ 1 + @test maximum(abs, integ.u .- ref.u) < 1.0e-3 + end + + @testset "Tighter tolerances give more steps and smaller errors" begin + results = map((1.0e-2, 1.0e-4, 1.0e-6)) do reltol + integ = DiffEqBase.init( + prob, PPLTG((Tsit5(), Tsit5())); + dt = 0.1, reltol, abstol = reltol * 1.0e-2 + ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + (err = maximum(abs, integ.u .- trueu), naccept = integ.stats.naccept) + end + @test issorted(collect(r.naccept for r in results)) + @test issorted(collect(r.err for r in results); rev = true) + end + + @testset "Loose inner tolerances choke the estimator (documented)" begin + # docs/src/topics/adaptivity.md: the pair difference contains the inner + # solver error, which for sub-problems with fast internal dynamics tracks + # the *inner tolerance* rather than the splitting dt. Inner tolerances + # looser than the splitting tolerances then put a dt-independent floor + # into EEst and the controller shrinks dt until DtLessThanMin. The fast + # dynamics matter: on a slow smooth problem a high-order leaf undershoots + # a loose tolerance by orders of magnitude and no floor appears. + odeAslow(du, u, p, t) = (du[1] = -u[1]; du[2] = -1.01 * u[2]; nothing) + odeBfast(du, u, p, t) = (du[1] = 1000.0 * u[2]; du[2] = -1000.0 * u[1]; nothing) + f_fast = GenericSplitFunction( + (ODEFunction(odeAslow), ODEFunction(odeBfast)), (dofs, dofs) + ) + prob_fast = OperatorSplittingProblem(f_fast, copy(u0), tspan) + trueu_fast = exp([-1.0 1000.0; -1000.0 -1.01]) * u0 + + function solve_with_leaf_tols(atol_leaf, rtol_leaf; kwargs...) + abstol = TreeOption(f_fast, 1.0e-6) # splitting node target + reltol = TreeOption(f_fast, 1.0e-4) + abstol[1] = atol_leaf + abstol[2] = atol_leaf + reltol[1] = rtol_leaf + reltol[2] = rtol_leaf + dtmin = TreeOption(f_fast, 0.0) + dtmin[] = 1.0e-3 # splitting node only! + integ = DiffEqBase.init( + prob_fast, PPLTG((Tsit5(), Tsit5())); + dt = 0.1, abstol, reltol, dtmin, verbose = false, kwargs... + ) + DiffEqBase.solve!(integ) + return integ + end + + # Leaves 100x looser than the splitting tolerances: EEst floor, abort. + choked = solve_with_leaf_tols(1.0e-4, 1.0e-2) + @test choked.sol.retcode == ReturnCode.DtLessThanMin + + # Same splitting tolerances and dtmin, leaves 10^4x tighter: fine. + # (Global error ~ splitting reltol accumulated over the ~60 steps.) + healthy = solve_with_leaf_tols(1.0e-10, 1.0e-8) + @test healthy.sol.retcode == ReturnCode.Success + @test maximum(abs, healthy.u .- trueu_fast) < 5.0e-3 + end + + @testset "Unreachable tolerances abort with DtLessThanMin" begin + integ = DiffEqBase.init( + prob, PPLTG((Tsit5(), Tsit5())); + dt = 0.1, reltol = 1.0e-12, abstol = 1.0e-14, dtmin = 0.01, + verbose = false + ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.DtLessThanMin + end + + @testset "Non-adaptive PPLTG runs fixed steps" begin + integ = DiffEqBase.init( + prob, PPLTG((Euler(), Euler())); dt = 0.05, adaptive = false + ) + @test integ.controller_cache === nothing + @test isnan(integ.EEst) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.iter == 20 + @test integ.dtcache ≈ 0.05 + @test isnan(integ.EEst) + end + + @testset "reject_step! restores the whole subtree, including leaf u" begin + # StrangMarchuk skips the forward sync of its first child on a retry + # (`next_sync_is_continuous`), so the rollback itself has to restore leaf + # states -- rewinding only the clocks would silently resume the retry from + # the failed attempt's state. + integ = DiffEqBase.init(prob, StrangMarchuk((Euler(), Euler())); dt = 0.05) + DiffEqBase.step!(integ) + for child in integ.child_subintegrators + child.u .+= 1.0 # pollute, as if a failed attempt had advanced the child + child.t += 0.5 + end + OS.reject_step!(integ) + @test integ.u == integ.uprev + for child in integ.child_subintegrators + @test child.u == integ.u[dofs] + @test child.t == integ.t + @test child.tprev == integ.t + end + end + + @testset "Per-node adaptive defaults" begin + # Without an `adaptive` keyword every node adapts exactly if its own + # algorithm can: the splitting node stays fixed-step while a Tsit5 leaf + # adapts and an Euler leaf does not. + integ = DiffEqBase.init(prob, LieTrotterGodunov((Tsit5(), Euler())); dt = 0.1) + @test integ.opts.adaptive == false + @test integ.controller_cache === nothing + @test integ.child_subintegrators[1].opts.adaptive == true + @test integ.child_subintegrators[2].opts.adaptive == false + end + + @testset "discontinuity_detection is refused up front" begin + @test_throws ArgumentError DiffEqBase.init( + prob, PPLTG((Tsit5(), Tsit5())); + dt = 0.1, + controller = OrdinaryDiffEqCore.IController(discontinuity_detection = true) + ) + end + + @testset "Adaptive root with a nested non-adaptive splitting node" begin + # B split once more into two identical halves handled by an inner LTG node. + # The inner node stays non-adaptive by default, and a rejection at the root + # has to roll the whole subtree back. + odeBhalf(du, u, p, t) = (du[1] = 0.25 * u[2]; du[2] = 0.25 * u[1]; nothing) + fBh = ODEFunction(odeBhalf) + f_nested = GenericSplitFunction( + (fA, GenericSplitFunction((fBh, fBh), (dofs, dofs))), (dofs, dofs) + ) + prob_nested = OperatorSplittingProblem(f_nested, u0, tspan) + alg_nested = PPLTG((Tsit5(), LieTrotterGodunov((Tsit5(), Tsit5())))) + + integ = DiffEqBase.init( + prob_nested, alg_nested; dt = 0.4, reltol = 1.0e-5, abstol = 1.0e-8 + ) + @test integ.controller_cache !== nothing + @test integ.child_subintegrators[2].controller_cache === nothing + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.stats.nreject ≥ 1 + @test maximum(abs, integ.u .- trueu) < 1.0e-3 + end + + @testset "Three-operator PPLTG adapts" begin + # A plus the strictly upper and lower triangles of B: pairwise non-commuting, + # so the N-ary pair difference measures a genuine splitting error. + odeB1(du, u, p, t) = (du[1] = 0.5 * u[2]; du[2] = 0.0; nothing) + odeB2(du, u, p, t) = (du[1] = 0.0; du[2] = 0.5 * u[1]; nothing) + fsplit3 = GenericSplitFunction( + (fA, ODEFunction(odeB1), ODEFunction(odeB2)), (dofs, dofs, dofs) + ) + prob3 = OperatorSplittingProblem(fsplit3, u0, tspan) + integ = DiffEqBase.init( + prob3, PPLTG((Tsit5(), Tsit5(), Tsit5())); + dt = 0.5, reltol = 1.0e-6, abstol = 1.0e-8 + ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.stats.nreject ≥ 1 + @test maximum(abs, integ.u .- trueu) < 1.0e-3 + end + + @testset "Nested PPLTG: the inner node adapts and rejects on its own" begin + # B split into its strictly lower and upper triangle, which do not commute, + # handled by an inner adaptive PPLTG. The inner node gets a much tighter + # splitting tolerance than the root, so it has to reject and subcycle + # within the intervals the root hands it. + odeBu(du, u, p, t) = (du[1] = 0.5 * u[2]; du[2] = 0.0; nothing) + odeBl(du, u, p, t) = (du[1] = 0.0; du[2] = 0.5 * u[1]; nothing) + f_nested = GenericSplitFunction( + (fA, GenericSplitFunction((ODEFunction(odeBu), ODEFunction(odeBl)), (dofs, dofs))), + (dofs, dofs) + ) + prob_nested = OperatorSplittingProblem(f_nested, u0, tspan) + alg_nested = PPLTG((Tsit5(), PPLTG((Tsit5(), Tsit5())))) + + reltol = TreeOption(f_nested, 1.0e-3) + reltol[2] = 1.0e-8 + integ = DiffEqBase.init( + prob_nested, alg_nested; dt = 0.5, reltol, abstol = 1.0e-10 + ) + sub = integ.child_subintegrators[2] + @test sub.controller_cache !== nothing + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + # The inner node ran its own controller: it rejected at least once and + # accepted far more (sub-cycled) steps than the root. + @test sub.stats.nreject ≥ 1 + @test sub.stats.naccept > integ.stats.naccept + @test sub.EEst ≤ 1 + @test maximum(abs, integ.u .- trueu) < 1.0e-2 + end + + @testset "PIController threads its state between steps" begin + controller = OrdinaryDiffEqCore.PIController(0.35, 0.2) + integ = DiffEqBase.init( + prob, PPLTG((Tsit5(), Tsit5())); + dt = 0.1, controller, reltol = 1.0e-6, abstol = 1.0e-8 + ) + @test integ.controller_cache isa OrdinaryDiffEqCore.PIControllerCache + errold0 = integ.controller_cache.errold + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.controller_cache.errold != errold0 # the error history was carried + @test maximum(abs, integ.u .- trueu) < 1.0e-3 + ufinal = copy(integ.u) + niters = integ.iter + # reinit! resets the controller memory, so a rerun is bit-identical. + DiffEqBase.reinit!(integ) + @test integ.controller_cache.errold == errold0 + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.u == ufinal + @test integ.iter == niters + end + + @testset "Explicitly passed controller is used" begin + integ = DiffEqBase.init( + prob, PPLTG((Tsit5(), Tsit5())); + dt = 0.1, controller = OrdinaryDiffEqCore.IController() + ) + @test integ.controller_cache isa OrdinaryDiffEqCore.IControllerCache + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + end +end diff --git a/test/backward.jl b/test/backward.jl new file mode 100644 index 0000000..ebc99d8 --- /dev/null +++ b/test/backward.jl @@ -0,0 +1,99 @@ +using OrdinaryDiffEqOperatorSplitting +using Test + +import DiffEqBase: DiffEqBase, ODEFunction +import SciMLBase +import SciMLBase: ReturnCode +using OrdinaryDiffEqTsit5 + +# Same non-commuting pair as in test/adaptivity.jl. Integrating the exact solution +# at t = 1 backward to t = 0 has to recover u0. +M = [-1.0 0.5; 0.5 -2.0] +odeA(du, u, p, t) = (du[1] = -u[1]; du[2] = -2 * u[2]; nothing) +odeB(du, u, p, t) = (du[1] = 0.5 * u[2]; du[2] = 0.5 * u[1]; nothing) +fA = ODEFunction(odeA) +fB = ODEFunction(odeB) + +dofs = [1, 2] +u0 = [1.0, 1.0] +uT = exp(M) * u0 +fsplit = GenericSplitFunction((fA, fB), (dofs, dofs)) + +PPLTG = PalindromicPairLieTrotterGodunov + +@testset "Backward-in-time integration" begin + @testset "Fixed-step consistency | $(nameof(typeof(alg)))" for (alg, atol) in ( + (LieTrotterGodunov((Tsit5(), Tsit5())), 1.0e-2), + (StrangMarchuk((Tsit5(), Tsit5())), 1.0e-4), + (PPLTG((Tsit5(), Tsit5())), 1.0e-4), + ) + prob = OperatorSplittingProblem(fsplit, copy(uT), (1.0, 0.0)) + integ = DiffEqBase.init(prob, alg; dt = 0.01, adaptive = false) + @test integ.dt < 0 + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.t ≈ 0.0 + @test integ.iter == 100 + @test maximum(abs, integ.u .- u0) < atol + + # A reinitialized backward solve is deterministic. + ufinal = copy(integ.u) + DiffEqBase.reinit!(integ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.u == ufinal + end + + @testset "Forward-backward roundtrip" begin + prob_fwd = OperatorSplittingProblem(fsplit, copy(u0), (0.0, 1.0)) + fwd = DiffEqBase.init(prob_fwd, PPLTG((Tsit5(), Tsit5())); dt = 0.01, adaptive = false) + DiffEqBase.solve!(fwd) + @test fwd.sol.retcode == ReturnCode.Success + + prob_bwd = OperatorSplittingProblem(fsplit, copy(fwd.u), (1.0, 0.0)) + bwd = DiffEqBase.init(prob_bwd, PPLTG((Tsit5(), Tsit5())); dt = 0.01, adaptive = false) + DiffEqBase.solve!(bwd) + @test bwd.sol.retcode == ReturnCode.Success + @test maximum(abs, bwd.u .- u0) < 1.0e-4 + end + + @testset "Adaptivity works backward" begin + prob = OperatorSplittingProblem(fsplit, copy(uT), (1.0, 0.0)) + integ = DiffEqBase.init( + prob, PPLTG((Tsit5(), Tsit5())); + dt = 0.5, reltol = 1.0e-6, abstol = 1.0e-8 + ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.t ≈ 0.0 + # The controller both rejected oversized steps and kept dt pointing backward. + @test integ.stats.nreject ≥ 1 + @test integ.dt < 0 + @test maximum(abs, integ.u .- u0) < 1.0e-5 + end + + @testset "Nested splitting works backward" begin + odeB1(du, u, p, t) = (du[1] = 0.5 * u[2]; du[2] = 0.0; nothing) + odeB2(du, u, p, t) = (du[1] = 0.0; du[2] = 0.5 * u[1]; nothing) + f_nested = GenericSplitFunction( + (fA, GenericSplitFunction((ODEFunction(odeB1), ODEFunction(odeB2)), (dofs, dofs))), + (dofs, dofs) + ) + prob = OperatorSplittingProblem(f_nested, copy(uT), (1.0, 0.0)) + alg = PPLTG((Tsit5(), PPLTG((Tsit5(), Tsit5())))) + integ = DiffEqBase.init(prob, alg; dt = 0.1, reltol = 1.0e-6, abstol = 1.0e-8) + sub = integ.child_subintegrators[2] + @test sub.tdir < 0 + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.t ≈ 0.0 + @test sub.t ≈ 0.0 + @test maximum(abs, integ.u .- u0) < 1.0e-5 + end + + @testset "reinit! cannot flip the direction" begin + prob = OperatorSplittingProblem(fsplit, copy(uT), (1.0, 0.0)) + integ = DiffEqBase.init(prob, PPLTG((Tsit5(), Tsit5())); dt = 0.01, adaptive = false) + @test_throws ErrorException DiffEqBase.reinit!(integ; t0 = 0.0, tf = 1.0) + end +end diff --git a/test/convergence.jl b/test/convergence.jl new file mode 100644 index 0000000..bb836d8 --- /dev/null +++ b/test/convergence.jl @@ -0,0 +1,70 @@ +using OrdinaryDiffEqOperatorSplitting +using Test + +import DiffEqBase: DiffEqBase, ODEFunction +import SciMLBase +using OrdinaryDiffEqTsit5 + +# Convergence orders of every splitting algorithm, in both integration directions, +# on pairwise non-commuting operators so the splitting error dominates. The inner +# solves use adaptive Tsit5, which resolves the sub-problems well below the +# splitting error at these step sizes. +# +# A = diag(-1, -2), and B = 0.5·offdiag split once more into its strictly upper +# and lower triangles Bu, Bl: [A, B], [A, Bu], [A, Bl], [Bu, Bl] are all nonzero. +M = [-1.0 0.5; 0.5 -2.0] +odeA(du, u, p, t) = (du[1] = -u[1]; du[2] = -2 * u[2]; nothing) +odeB(du, u, p, t) = (du[1] = 0.5 * u[2]; du[2] = 0.5 * u[1]; nothing) +odeBu(du, u, p, t) = (du[1] = 0.5 * u[2]; du[2] = 0.0; nothing) +odeBl(du, u, p, t) = (du[1] = 0.0; du[2] = 0.5 * u[1]; nothing) +fA = ODEFunction(odeA) +fB = ODEFunction(odeB) +fBu = ODEFunction(odeBu) +fBl = ODEFunction(odeBl) + +dofs = [1, 2] +u0 = [1.0, 1.0] +uT = exp(M) * u0 + +f_two = GenericSplitFunction((fA, fB), (dofs, dofs)) +f_three = GenericSplitFunction((fA, fBu, fBl), (dofs, dofs, dofs)) +f_nested = GenericSplitFunction( + (fA, GenericSplitFunction((fBu, fBl), (dofs, dofs))), (dofs, dofs) +) + +function convergence_rates(f, alg, tspan, ustart, utarget; dts = (0.1, 0.05, 0.025)) + errs = map(dts) do dt + prob = OperatorSplittingProblem(f, copy(ustart), tspan) + integ = DiffEqBase.init(prob, alg; dt, adaptive = false, verbose = false) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == SciMLBase.ReturnCode.Success + maximum(abs, integ.u .- utarget) + end + return [log2(errs[i] / errs[i + 1]) for i in 1:(length(errs) - 1)] +end + +@testset "Convergence order" begin + for (AlgType, expected_order) in ( + (LieTrotterGodunov, 1), + (StrangMarchuk, 2), + (PalindromicPairLieTrotterGodunov, 2), + ) + cases = ( + ("two operators", f_two, AlgType((Tsit5(), Tsit5()))), + ("three operators", f_three, AlgType((Tsit5(), Tsit5(), Tsit5()))), + ("nested", f_nested, AlgType((Tsit5(), AlgType((Tsit5(), Tsit5()))))), + ) + directions = ( + ("forward", (0.0, 1.0), u0, uT), + ("backward", (1.0, 0.0), uT, u0), + ) + @testset "$(nameof(AlgType)) | $case | $dir (order $expected_order)" for + (case, f, alg) in cases, (dir, tspan, ustart, utarget) in directions + + rates = convergence_rates(f, alg, tspan, ustart, utarget) + for rate in rates + @test rate ≈ expected_order atol = 0.3 + end + end + end +end diff --git a/test/failure_escalation.jl b/test/failure_escalation.jl new file mode 100644 index 0000000..7f48347 --- /dev/null +++ b/test/failure_escalation.jl @@ -0,0 +1,145 @@ +using OrdinaryDiffEqOperatorSplitting +import OrdinaryDiffEqOperatorSplitting as OS +using Test + +import DiffEqBase: DiffEqBase, ODEFunction +import SciMLBase +import SciMLBase: ReturnCode +using OrdinaryDiffEqLowOrderRK +using OrdinaryDiffEqTsit5 + +# Failure protocol: a failing *adaptive* node is fatal (it exhausted its own +# adaptation); a failing non-adaptive node escalates to the nearest adaptive +# ancestor, which retries on a failfactor-shrunken interval; a non-adaptive root +# stops with the escalated retcode. + +M = [-1.0 0.5; 0.5 -2.0] +odeA(du, u, p, t) = (du[1] = -u[1]; du[2] = -2 * u[2]; nothing) +fA = ODEFunction(odeA) +dofs = [1, 2] +u0 = [1.0, 1.0] +trueu = exp(M) * u0 + +# The coupling operator B (or a part of it), but producing NaNs for a few +# evaluations once `t` first crosses `tfail`: a transient inner failure in the +# middle of the solve. Time-triggered (leaf integrators evaluate `f` at init, so a +# purely count-based trigger would be consumed before any stepping), and emitting +# several NaNs (a single one lands in the FSAL cache, which the forward sync's +# u_modified re-evaluation heals without the step ever failing). +mutable struct TransientFailure{F} + countdown::Int + triggered::Bool + tfail::Float64 + f::F +end +TransientFailure(tfail, f) = TransientFailure(3, false, tfail, f) +function (tf::TransientFailure)(du, u, p, t) + tf.triggered |= t >= tf.tfail + if tf.triggered && tf.countdown > 0 + tf.countdown -= 1 + du .= NaN + else + tf.f(du, u, p, t) + end + return nothing +end +odeB(du, u, p, t) = (du[1] = 0.5 * u[2]; du[2] = 0.5 * u[1]; nothing) +odeB1(du, u, p, t) = (du[1] = 0.5 * u[2]; du[2] = 0.0; nothing) +odeB2(du, u, p, t) = (du[1] = 0.0; du[2] = 0.5 * u[1]; nothing) +ode_nan(du, u, p, t) = (du .= NaN; nothing) + +PPLTG = PalindromicPairLieTrotterGodunov + +@testset "Failure escalation" begin + @testset "Transient non-adaptive failure recovers via the adaptive root" begin + ffail = TransientFailure(0.35, odeB) + f = GenericSplitFunction((fA, ODEFunction(ffail)), (dofs, dofs)) + prob = OperatorSplittingProblem(f, copy(u0), (0.0, 1.0)) + # Per-node defaults: adaptive PPLTG root, non-adaptive Euler leaf for B. + integ = DiffEqBase.init(prob, PPLTG((Tsit5(), Euler())); dt = 0.1) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.stats.nreject ≥ 1 # the failed attempt was retried + @test ffail.triggered # ... and the failure was actually hit + @test maximum(abs, integ.u .- trueu) < 0.05 + end + + @testset "Adaptive child failure is fatal, without retries" begin + f = GenericSplitFunction((fA, ODEFunction(ode_nan)), (dofs, dofs)) + prob = OperatorSplittingProblem(f, copy(u0), (0.0, 1.0)) + integ = DiffEqBase.init( + prob, PPLTG((Tsit5(), Tsit5())); dt = 0.1, verbose = false + ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode ∈ + (ReturnCode.Unstable, ReturnCode.DtNaN, ReturnCode.DtLessThanMin) + # The fatal branch stops immediately instead of shrinking dt to dtmin. + @test integ.iter ≤ 2 + end + + @testset "Failure escalates past a non-adaptive intermediate node" begin + ffail = TransientFailure(0.35, odeB2) + f_nested = GenericSplitFunction( + (fA, GenericSplitFunction((ODEFunction(odeB1), ODEFunction(ffail)), (dofs, dofs))), + (dofs, dofs) + ) + prob = OperatorSplittingProblem(f_nested, copy(u0), (0.0, 1.0)) + # The inner LieTrotterGodunov node cannot adapt, so the failure of its + # Euler leaf has to bubble up to the PPLTG root. + alg = PPLTG((Tsit5(), LieTrotterGodunov((Euler(), Euler())))) + integ = DiffEqBase.init(prob, alg; dt = 0.1) + sub = integ.child_subintegrators[2] + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.stats.nreject ≥ 1 + @test ffail.triggered + @test SciMLBase.successful_retcode(sub.status.retcode) + @test maximum(abs, integ.u .- trueu) < 0.05 + end + + @testset "A transient failure under a non-adaptive root still stops" begin + ffail = TransientFailure(0.35, odeB) + f = GenericSplitFunction((fA, ODEFunction(ffail)), (dofs, dofs)) + prob = OperatorSplittingProblem(f, copy(u0), (0.0, 1.0)) + integ = DiffEqBase.init( + prob, LieTrotterGodunov((Euler(), Euler())); dt = 0.1, adaptive = false + ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Unstable + @test ffail.triggered + @test integ.t < 1.0 # aborted mid-solve, nobody could retry + end + + @testset "Rollback clears sticky Newton-style leaf failure flags" begin + # SciMLBase's generic check_error re-derives ConvergenceFailure from + # `last_stepfail` on non-adaptive leaves even after the retcode is reset, + # so the escalation retry has to clear the flag during rollback or it dies + # again immediately. An implicit inner solver sets the flag organically on + # a Newton failure; it is injected here because the test dependencies are + # explicit solvers only. + f = GenericSplitFunction((fA, ODEFunction(odeB)), (dofs, dofs)) + prob = OperatorSplittingProblem(f, copy(u0), (0.0, 1.0)) + integ = DiffEqBase.init(prob, PPLTG((Tsit5(), Euler())); dt = 0.1, verbose = false) + DiffEqBase.step!(integ) + + 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 ... + + OS.reject_step!(integ) # ... and the retry's rollback clears it + @test !leaf.last_stepfail + @test !OS._child_failed(leaf) + @test SciMLBase.check_error(leaf) == ReturnCode.Success + end + + @testset "Persistent non-adaptive failure exhausts the adaptive root" begin + f = GenericSplitFunction((fA, ODEFunction(ode_nan)), (dofs, dofs)) + prob = OperatorSplittingProblem(f, copy(u0), (0.0, 1.0)) + integ = DiffEqBase.init( + prob, PPLTG((Tsit5(), Euler())); dt = 0.1, verbose = false + ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.DtLessThanMin + end +end diff --git a/test/operator_splitting_api.jl b/test/operator_splitting_api.jl index 99d786b..ba944a9 100644 --- a/test/operator_splitting_api.jl +++ b/test/operator_splitting_api.jl @@ -1,5 +1,4 @@ using OrdinaryDiffEqOperatorSplitting -import OrdinaryDiffEqOperatorSplitting: OrdinaryDiffEqOperatorSplitting as OS using Test import SciMLBase: ReturnCode @@ -63,77 +62,14 @@ eqs = [ @named testmodel2 = System(eqs, time) testsys2 = mtkcompile(testmodel2; sort_eqs = false) -# Test whether adaptive code path works in principle -struct FakeAdaptiveAlgorithm{T, T2} <: OS.AbstractOperatorSplittingAlgorithm - alg::T - inner_algs::T2 # delegate inner_algs to the wrapped algorithm -end -FakeAdaptiveAlgorithm(alg) = FakeAdaptiveAlgorithm(alg, alg.inner_algs) - -struct FakeAdaptiveAlgorithmCache{T} <: OS.AbstractOperatorSplittingCache - cache::T -end - -@inline DiffEqBase.isadaptive(::FakeAdaptiveAlgorithm) = true - -@inline function OS.stepsize_controller!( - integrator::OS.OperatorSplittingIntegrator, alg::FakeAdaptiveAlgorithm - ) - return nothing -end - -@inline function OS.step_accept_controller!( - integrator::OS.OperatorSplittingIntegrator, alg::FakeAdaptiveAlgorithm, q - ) - integrator.dt = integrator.dtcache - return nothing -end -@inline function OS.step_reject_controller!( - integrator::OS.OperatorSplittingIntegrator, alg::FakeAdaptiveAlgorithm, q - ) - error("The tests should never run into this scenario!") - return nothing -end - -# Override init_cache to wrap the inner cache in FakeAdaptiveAlgorithmCache -function OS.init_cache( - f::GenericSplitFunction, alg::FakeAdaptiveAlgorithm; - kwargs... - ) - inner_cache = OS.init_cache(f, alg.alg; kwargs...) - return FakeAdaptiveAlgorithmCache(inner_cache) -end - -@inline DiffEqBase.get_tmp_cache( - integrator::OS.OperatorSplittingIntegrator, - alg::OS.AbstractOperatorSplittingAlgorithm, - cache::FakeAdaptiveAlgorithmCache -) = DiffEqBase.get_tmp_cache(integrator, alg, cache.cache) - -@inline function OS._perform_step!( - outer_integrator, - subintegrators::Tuple, - cache::FakeAdaptiveAlgorithmCache, - dt - ) - return OS._perform_step!( - outer_integrator, subintegrators, cache.cache, dt - ) -end - -FakeAdaptiveLTG(inner) = FakeAdaptiveAlgorithm(LieTrotterGodunov(inner)) -FakeAdaptiveSM(inner) = FakeAdaptiveAlgorithm(StrangMarchuk(inner)) - -function Base.show(io::IO, alg::FakeAdaptiveAlgorithm) - print(io, "FAKE (") - Base.show(io, alg.alg) - return print(io, ")") -end - -# StrangMarchuk steps child 1 twice per outer step (two half-steps). +# Steps a child takes per outer step: StrangMarchuk steps child 1 twice (two +# half-steps), the palindromic pair steps every child twice (once per sequence). _sub1_iter_factor(::LieTrotterGodunov) = 1 _sub1_iter_factor(::StrangMarchuk) = 2 -_sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg) +_sub1_iter_factor(::PalindromicPairLieTrotterGodunov) = 2 +_sub2_iter_factor(::LieTrotterGodunov) = 1 +_sub2_iter_factor(::StrangMarchuk) = 1 +_sub2_iter_factor(::PalindromicPairLieTrotterGodunov) = 2 # --------------------------------------------------------------------------- @@ -163,7 +99,7 @@ _sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg) nsteps = ceil(Int, (tspan[2] - tspan[1]) / dt) - for TimeStepperType in (LieTrotterGodunov, FakeAdaptiveLTG, StrangMarchuk, FakeAdaptiveSM) + for TimeStepperType in (LieTrotterGodunov, StrangMarchuk, PalindromicPairLieTrotterGodunov) @testset "$tstepper" for (prob, tstepper) in ( (prob1a, TimeStepperType((Euler(), Euler()))), (prob1a, TimeStepperType((Tsit5(), Euler()))), @@ -188,6 +124,7 @@ _sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg) sub1 = integrator.child_subintegrators[1] sub2 = integrator.child_subintegrators[2] expected_sub1_iters = _sub1_iter_factor(tstepper) * nsteps + expected_sub2_iters = _sub2_iter_factor(tstepper) * nsteps DiffEqBase.solve!(integrator) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success @@ -201,7 +138,7 @@ _sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg) @test sub1.iter == expected_sub1_iters @test sub2.t ≈ tspan[2] - @test sub2.iter == nsteps + @test sub2.iter == expected_sub2_iters DiffEqBase.reinit!(integrator; dt = dt) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default @@ -233,58 +170,60 @@ _sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg) @test sub1.iter == expected_sub1_iters @test sub2.t ≈ tspan[2] - @test sub2.iter == nsteps + @test sub2.iter == expected_sub2_iters end end - for TimeStepperType in (FakeAdaptiveLTG, FakeAdaptiveSM) - @testset "Adaptive solver type $TimeStepperType | $tstepper" for (prob, tstepper) in ( - (prob1a, TimeStepperType((Tsit5(), Tsit5()))), - (prob2, TimeStepperType((Tsit5(), TimeStepperType((Tsit5(), Tsit5()))))), - ) - integrator = DiffEqBase.init( - prob, tstepper, dt = dt, verbose = true, alias_u0 = false, adaptive = true - ) - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default - DiffEqBase.solve!(integrator) - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success - ufinal = copy(integrator.u) - @test isapprox(ufinal, trueu, atol = 1.0e-6) - @test integrator.t ≈ tspan[2] - @test integrator.dtcache ≈ dt - @test integrator.iter == nsteps - - DiffEqBase.reinit!(integrator; dt = dt) - @test integrator.dt == dt - @test integrator.dt == integrator.dtcache - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default - for (u, t) in TimeChoiceIterator(integrator, tspan[1]:5.0:tspan[2]) - end - @test isapprox(ufinal, integrator.u, atol = 1.0e-12) - @test integrator.t ≈ tspan[2] - @test integrator.dtcache ≈ dt - @test integrator.iter == nsteps - - DiffEqBase.reinit!(integrator; dt = dt) - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default - for (uprev, tprev, u, t) in intervals(integrator) - end - @test isapprox(ufinal, integrator.u, atol = 1.0e-12) - @test integrator.t ≈ tspan[2] - @test integrator.dtcache ≈ dt - @test integrator.iter == nsteps + @testset "Adaptive splitting | $tstepper" for (prob, tstepper) in ( + (prob1a, PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5()))), + ( + prob2, PalindromicPairLieTrotterGodunov( + (Tsit5(), PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5()))) + ), + ), + ) + # PPLTG is adaptive by default; the integrator interface has to keep working + # while the controller reshapes the step sequence. + integrator = DiffEqBase.init( + prob, tstepper, dt = dt, verbose = true, alias_u0 = false + ) + @test integrator.opts.adaptive + @test integrator.controller_cache !== nothing + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default + DiffEqBase.solve!(integrator) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success + ufinal = copy(integrator.u) + @test isapprox(ufinal, trueu, atol = 1.0e-6) + @test integrator.t ≈ tspan[2] + niters = integrator.iter + + # A reinitialized adaptive solve is deterministic. + DiffEqBase.reinit!(integrator; dt = dt) + @test integrator.dt == dt + @test integrator.dt == integrator.dtcache + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default + DiffEqBase.solve!(integrator) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success + @test isapprox(ufinal, integrator.u, atol = 1.0e-12) + @test integrator.t ≈ tspan[2] + @test integrator.iter == niters + + # The iteration protocols land on the same time points despite the + # controller-driven step sequence in between. + DiffEqBase.reinit!(integrator; dt = dt) + for (u, t) in TimeChoiceIterator(integrator, tspan[1]:5.0:tspan[2]) + end + @test isapprox(integrator.u, trueu, atol = 1.0e-6) + @test integrator.t ≈ tspan[2] - DiffEqBase.reinit!(integrator; dt = dt) - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default - DiffEqBase.solve!(integrator) - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success - @test integrator.t ≈ tspan[2] - @test integrator.dtcache ≈ dt - @test integrator.iter == nsteps + DiffEqBase.reinit!(integrator; dt = dt) + for (uprev, tprev, u, t) in intervals(integrator) end + @test isapprox(integrator.u, trueu, atol = 1.0e-6) + @test integrator.t ≈ tspan[2] end - @testset "StrangMarchuk with 3 operators" begin + @testset "Three operators" begin dt = 0.01π # f1 + f3 + f3 = f1 + f2, so the reference solution is the same trueu. f1dofs = [1, 2, 3] @@ -293,10 +232,18 @@ _sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg) prob3 = OperatorSplittingProblem(fsplit3, u0, tspan) nsteps = ceil(Int, (tspan[2] - tspan[1]) / dt) + # StrangMarchuk solves all but the last operator twice (two half-steps); + # the palindromic pair solves every operator once per sequence. + sub_iter_factors(::StrangMarchuk) = (2, 2, 1) + sub_iter_factors(::PalindromicPairLieTrotterGodunov) = (2, 2, 2) + @testset "$tstepper" for tstepper in ( StrangMarchuk((Euler(), Euler(), Euler())), StrangMarchuk((Tsit5(), Euler(), Tsit5())), StrangMarchuk((Tsit5(), Tsit5(), Tsit5())), + PalindromicPairLieTrotterGodunov((Euler(), Euler(), Euler())), + PalindromicPairLieTrotterGodunov((Tsit5(), Euler(), Tsit5())), + PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5(), Tsit5())), ) integrator = DiffEqBase.init( prob3, tstepper, dt = dt, verbose = true, alias_u0 = false, adaptive = false @@ -307,61 +254,15 @@ _sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg) @test integrator.t ≈ tspan[2] @test integrator.iter == nsteps - sub1 = integrator.child_subintegrators[1] - sub2 = integrator.child_subintegrators[2] - sub3 = integrator.child_subintegrators[3] - # Palindromic: children 1 & 2 get two half-steps, child 3 gets one full step - @test sub1.iter == 2 * nsteps - @test sub2.iter == 2 * nsteps - @test sub3.iter == nsteps - end - end - - @testset "Convergence order" begin - # Use non-commuting operators so splitting error is non-zero. - # A = diag(-1,-2), B = [0 0.5; 0.5 0] have [A,B] ≠ 0. - function ode_conv_A(du, u, p, t) - du[1] = -u[1] - return du[2] = -2 * u[2] - end - function ode_conv_B(du, u, p, t) - du[1] = 0.5 * u[2] - return du[2] = 0.5 * u[1] - end - fA = ODEFunction(ode_conv_A) - fB = ODEFunction(ode_conv_B) - - conv_tspan = (0.0, 1.0) - conv_u0 = [1.0, 1.0] - conv_trueu = exp(conv_tspan[2] * [-1.0 0.5; 0.5 -2.0]) * conv_u0 - - conv_dofs = [1, 2] - fsplit_conv = GenericSplitFunction((fA, fB), (conv_dofs, conv_dofs)) - prob_conv = OperatorSplittingProblem(fsplit_conv, conv_u0, conv_tspan) - - dts = [0.1, 0.05, 0.025] - for (TimeStepperType, expected_order) in ( - (LieTrotterGodunov, 1), - (StrangMarchuk, 2), - ) - @testset "$TimeStepperType (order $expected_order)" begin - errors = map(dts) do dt_i - tstepper = TimeStepperType((Tsit5(), Tsit5())) - integrator = DiffEqBase.init( - prob_conv, tstepper, dt = dt_i, verbose = false, - alias_u0 = false, adaptive = false - ) - DiffEqBase.solve!(integrator) - maximum(abs, integrator.u .- conv_trueu) - end - for i in 1:(length(errors) - 1) - rate = log2(errors[i] / errors[i + 1]) - @test rate ≈ expected_order atol = 0.3 - end + for (i, factor) in pairs(sub_iter_factors(tstepper)) + @test integrator.child_subintegrators[i].t ≈ tspan[2] + @test integrator.child_subintegrators[i].iter == factor * nsteps end end end + # Convergence orders are covered systematically in test/convergence.jl. + @testset "Instability detection" begin dt = 0.01π @@ -377,7 +278,18 @@ _sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg) fsplit_NaN = GenericSplitFunction((f1, f_NaN), (f1dofs, f3dofs)) prob_NaN = OperatorSplittingProblem(fsplit_NaN, u0, tspan) - for TimeStepperType in (LieTrotterGodunov, StrangMarchuk) + for TimeStepperType in (LieTrotterGodunov, StrangMarchuk, PalindromicPairLieTrotterGodunov) + # An adaptive root (PPLTG by default) retries escalated non-adaptive + # failures until its dt reaches dtmin, so it may also end in + # DtLessThanMin instead of surfacing the child diagnosis directly. + expected_retcodes = if TimeStepperType === PalindromicPairLieTrotterGodunov + ( + DiffEqBase.ReturnCode.Unstable, DiffEqBase.ReturnCode.DtNaN, + DiffEqBase.ReturnCode.DtLessThanMin, + ) + else + (DiffEqBase.ReturnCode.Unstable, DiffEqBase.ReturnCode.DtNaN) + end @testset "Solver type $TimeStepperType | $tstepper" for tstepper in ( TimeStepperType((Euler(), Euler())), TimeStepperType((Tsit5(), Euler())), @@ -385,12 +297,11 @@ _sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg) TimeStepperType((Tsit5(), Tsit5())), ) integrator_NaN = DiffEqBase.init( - prob_NaN, tstepper, dt = dt, verbose = true, alias_u0 = false + prob_NaN, tstepper, dt = dt, verbose = false, alias_u0 = false ) @test integrator_NaN.sol.retcode == DiffEqBase.ReturnCode.Default DiffEqBase.solve!(integrator_NaN) - @test integrator_NaN.sol.retcode ∈ - (DiffEqBase.ReturnCode.Unstable, DiffEqBase.ReturnCode.DtNaN) + @test integrator_NaN.sol.retcode ∈ expected_retcodes end end end diff --git a/test/qa/qa.jl b/test/qa/qa.jl index 901b7c4..cd55cef 100644 --- a/test/qa/qa.jl +++ b/test/qa/qa.jl @@ -5,13 +5,8 @@ using Test run_qa( OrdinaryDiffEqOperatorSplitting; - # JET reports 2 genuine errors in src/integrator.jl on the SplitSubIntegrator - # rollback path: `rollback_children!(::SplitSubIntegrator)` has no matching - # method and `_rollback_children!` is called (line 514) but never defined. - # Tracked in https://github.com/SciML/OrdinaryDiffEqOperatorSplitting.jl/issues/87 # target_defined_modules scopes the report to this package's own modules (the - # default target_modules=(pkg,) filter hides these via-dependency-driven frames). - jet_broken = true, + # 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 @@ -25,10 +20,17 @@ run_qa( all_qualified_accesses_are_public = (; ignore = ( :__init, :__solve, :done, :postamble!, :solution_new_retcode, # SciMLBase - :DEFAULT_VERBOSE, :NAN_CHECK, :None, # DiffEqBase + :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, :post_newton_controller!, :timedepentdtmin, # 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 @@ -37,6 +39,9 @@ run_qa( 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 ), ), ), diff --git a/test/sync.jl b/test/sync.jl new file mode 100644 index 0000000..13244c2 --- /dev/null +++ b/test/sync.jl @@ -0,0 +1,74 @@ +using OrdinaryDiffEqOperatorSplitting +import OrdinaryDiffEqOperatorSplitting as OS +using Test + +import SciMLBase +import DiffEqBase: DiffEqBase, ODEFunction +using OrdinaryDiffEqLowOrderRK + +# Minimal stand-in for a leaf integrator whose buffers downstream packages may alias +# into the parent's buffers (e.g. views into the master solution for GPU setups). +mutable struct MockLeaf{U, UP} <: SciMLBase.AbstractODEIntegrator{Nothing, true, U, Float64} + u::U + uprev::UP + u_modified::Bool +end + +@static if isdefined(SciMLBase, :derivative_discontinuity!) + SciMLBase.derivative_discontinuity!(m::MockLeaf, b) = m.u_modified = b +end +@static if isdefined(DiffEqBase, :u_modified!) + DiffEqBase.u_modified!(m::MockLeaf, b) = m.u_modified = b +end + +@testset "Nested children address their parent's buffers" begin + # dof indices are parent-relative at every level, so a grandchild's initial + # state is parent-slice-of-parent-slice -- NOT the root vector indexed with + # parent-relative indices. Stock (copying) leaves hide a mix-up because the + # first forward sync overwrites their state; leaves wired as views into the + # handed buffer would alias the wrong root slots permanently. + ode1(du, u, p, t) = (du .= -0.1 .* u; nothing) + ode2(du, u, p, t) = (du[1] = -0.01 * u[2]; du[2] = -0.01 * u[1]; nothing) + u0 = [10.0, 20.0, 30.0] + f1dofs = [1, 2, 3] + f2dofs = [1, 3] + f3dofs = [1, 2] # relative to the inner node's [1, 3] slice → root slots 1 and 3 + inner = GenericSplitFunction((ODEFunction(ode2), ODEFunction(ode2)), (f3dofs, f3dofs)) + outer = GenericSplitFunction((ODEFunction(ode1), inner), (f1dofs, f2dofs)) + prob = OperatorSplittingProblem(outer, u0, (0.0, 1.0)) + alg = LieTrotterGodunov((Euler(), LieTrotterGodunov((Euler(), Euler())))) + integ = DiffEqBase.init(prob, alg; dt = 0.1, adaptive = false) + + sub = integ.child_subintegrators[2] + @test sub.u == u0[f2dofs] # [10, 30] + @test sub.child_subintegrators[1].u == u0[f2dofs][f3dofs] # [10, 30], not [10, 20] + @test sub.child_subintegrators[2].u == u0[f2dofs][f3dofs] +end + +@testset "forward_sync_internal!" begin + @testset "independent child buffers are synced" begin + u_parent = [1.0, 2.0, 3.0] + uprev_parent = [10.0, 20.0, 30.0] + idxs = [1, 2] + child = MockLeaf([0.0, 0.0], [0.0, 0.0], false) + + OS.forward_sync_internal!(u_parent, uprev_parent, child, idxs) + @test child.u == u_parent[idxs] + @test child.uprev == u_parent[idxs] + @test child.u_modified + end + + @testset "a child uprev aliasing the parent rollback buffer survives" begin + # A child whose uprev is a view into the parent's uprev: the sync must not + # scribble on the parent's rollback anchor mid-step. + u_parent = [1.0, 2.0, 3.0] + uprev_parent = [10.0, 20.0, 30.0] + idxs = [1, 2] + child = MockLeaf([0.0, 0.0], view(uprev_parent, idxs), false) + + OS.forward_sync_internal!(u_parent, uprev_parent, child, idxs) + @test child.u == u_parent[idxs] + @test uprev_parent == [10.0, 20.0, 30.0] # untouched + @test child.u_modified + end +end