From 659523c7664c336184a5f4ff0e870df7aff01acf Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 3 Aug 2026 14:04:10 -0400 Subject: [PATCH 1/4] Add saving, interpolation and callback support to the outer integrator MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Saving and events are implemented on `OperatorSplittingIntegrator` only. The inner splits are stages rather than steps, so their intermediate states do not approximate the split solution at any time point and there is nothing meaningful to save or to test a condition against there. This is enforced structurally rather than by convention: the saving options and `saveiter` are fields of the outer integrator alone, consumed as plain `__init` keywords so they never enter the `ConfigTree` and cannot leak to the leaves, and `save_step!`/`handle_callbacks!` are no-ops on a `SplitSubIntegrator`. Interpolation `(integrator)(t, Val{D}; idxs)` and `(integrator)(val, t, Val{D}; idxs)` now follow the SciML contract, delegating to a cache-dispatched `splitting_interpolant[!]` whose fallback is linear. A higher-order dense output is therefore one method per algorithm, with no changes to saving or callbacks. Note that Θ is derived from `t - tprev`, never from `integrator.dt`: after an accepted step `step_accept_controller!` has already replaced `dt` with the next proposal. `change_t_via_interpolation!` refills `u` from the interpolant and re-anchors the whole subintegrator tree through `rollback_children!`. It decodes `modify_save_endpoint` by dispatch because DiffEqBase passes the literal `Val{:false}` -- a `Val` of the Symbol -- which `if T` would choke on. Saving `SciMLBase.savevalues!` is the single integration point; DiffEqBase's callback code calls exactly that generic. `saveat`, `save_everystep`, `save_start`, `save_end` and `save_on` are supported. `saveat` points inside a step are interpolated, so asking for output never changes the step sequence and therefore never changes the splitting error. `dense = true` is rejected: the interpolant is linear, so `sol(t)` already reproduces it exactly from the saved points. Callbacks `DiscreteCallback`, `ContinuousCallback`, `VectorContinuousCallback` and `CallbackSet` all work, with `save_positions`, `affect_neg!`, `rootfind`, `interp_points`, `terminate!` and `derivative_discontinuity!`. The root finding and `affect!` dispatch are DiffEqBase's; this adds the per-step driver `handle_callbacks!`, the generated `apply_ith_callback!`, `initialize_callbacks!`, and the integrator-side methods DiffEqBase requires. The splitting-specific hazard is child state: every child owns a copy of its slice of `u`, so an `affect!` that modifies `integrator.u` would leave the tree stale, and the palindromic schemes skip the forward sync of their first child so it would not self-correct. `reeval_internals_due_to_modification!` runs `rollback_children!`, which both propagates the new state and restores the `parent.u == child.u` invariant that shortcut depends on. Bug fixes - `init(prob, alg; callback = cb)` threw `FieldError`: DiffEqBase reads `derivative_discontinuity` as a literal field, and this package named it `u_modified` while defining only the setter. Callbacks were unusable. - `postamble!` ran once per step, so callback `finalize` hooks fired every step. SciMLBase's generic `check_error!` finalizes whenever the code is not `Success`, and a mid-solve node legitimately reports `Default`; `check_error!` is now overridden and `postamble!` made idempotent. - `reinit!(integrator, u0)` silently returned wrong results for `StrangMarchuk`. A leaf's `reinit!` restores the `u0` slice captured when it was built, and Strang's skipped first sync meant the stale child state was integrated. `reinit!` now pushes the new state down explicitly, and initializes callbacks after the child reinit so a modifying initializer is not undone. - `reinit!(...; reinit_callbacks = false)` threw `BoundsError` indexing an empty `discrete_callbacks` tuple for a saving callback that never existed. - `handle_tstop!` called `change_t_via_interpolation!`, for which no method existed. - The `saveat` heap included `t0` and duplicated `tf`. - Removed dead `fix_solution_buffer_sizes!`, which referenced fields the integrator does not have. Accuracy limits, documented rather than left to be discovered: interpolated output and located event times are second order accurate, an event is the exact root of the linear interpolant over the bracketing step (so cap `dtmax` when event accuracy matters), and an event that reverses within a single step cannot be detected -- which also means raising `interp_points` cannot help. The SciMLBase lower bound moves to v3, matching the `derivative_discontinuity` field name. Every DiffEqBase 7.x already required SciMLBase 3.x, so the previous `2.77.0` entry was unreachable in practice. Co-Authored-By: Claude Opus 5 --- Project.toml | 2 +- docs/src/devdocs/index.md | 44 ++ docs/src/usage/index.md | 91 +++- src/OrdinaryDiffEqOperatorSplitting.jl | 4 + src/integrator.jl | 561 +++++++++++++++++++++++-- src/utils.jl | 17 +- test/callbacks.jl | 516 +++++++++++++++++++++++ test/qa/qa.jl | 4 + test/saving.jl | 439 +++++++++++++++++++ 9 files changed, 1637 insertions(+), 41 deletions(-) create mode 100644 test/callbacks.jl create mode 100644 test/saving.jl diff --git a/Project.toml b/Project.toml index f0c2129..5602f44 100644 --- a/Project.toml +++ b/Project.toml @@ -27,7 +27,7 @@ OrdinaryDiffEqTsit5 = "1.1.0, 2" PrecompileTools = "1.1" RecursiveArrayTools = "3.39.0, 4" SafeTestsets = "0.1.0" -SciMLBase = "2.77.0, 3.1" +SciMLBase = "3.1" SciMLIterators = "1" SciMLTesting = "2.1" SymbolicIndexingInterface = "0.3.36" diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md index 82708fc..bd63d37 100644 --- a/docs/src/devdocs/index.md +++ b/docs/src/devdocs/index.md @@ -108,3 +108,47 @@ end # Done :) end ``` + +## Dense output + +Saving, `saveat` and continuous callback root-finding all go through a single hook, +so an algorithm only has to describe how to interpolate *within one of its steps*: + +```julia +function OrdinaryDiffEqOperatorSplitting.splitting_interpolant( + integrator, cache::MySimpleFirstOrderCache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}} +) where {D} + # Θ = (t - tprev) / dt is the step-local coordinate, y₀ = u(tprev), y₁ = u(t). + # ... +end + +function OrdinaryDiffEqOperatorSplitting.splitting_interpolant!( + out, integrator, cache::MySimpleFirstOrderCache, Θ, dt, y₀, y₁, idxs, + ::Type{Val{D}} +) where {D} + # In-place variant; must return `out`. +end +``` + +Both fall back to linear interpolation for any `AbstractOperatorSplittingCache`, so +implementing them is optional. Note the accuracy consequences: with the linear +fallback, interpolated output is first order even for a second-order scheme, and +continuous callback event times are the exact roots of a straight chord across the +step. Implementing a higher-order interpolant improves saved output *and* event +location at once. + +The reason the fallback is only linear is structural rather than incidental. A +splitting step advances its children sequentially over staggered subintervals, so at +any time strictly inside the step the children's own interpolants describe different +sub-problems evaluated over different intervals, and they do not compose into an +approximation of the split solution. Only the step endpoints, which the outer +integrator owns, are states of the full split system. + +Two invariants matter when implementing this: + + - `Θ` is derived from `integrator.t - integrator.tprev`, **not** from + `integrator.dt`. Once a step is accepted, `step_accept_controller!` has already + replaced `dt` with the step size proposed for the *next* step. + - Interpolation must not mutate integrator state. `change_t_via_interpolation!` + relies on being able to evaluate the interpolant repeatedly (the callback + root-finder does so many times per step) before committing to a time. diff --git a/docs/src/usage/index.md b/docs/src/usage/index.md index 94e2bcd..00a71c3 100644 --- a/docs/src/usage/index.md +++ b/docs/src/usage/index.md @@ -42,8 +42,9 @@ alg = LieTrotterGodunov( (Euler(), Euler()) ) -# Right now OrdinaryDiffEqOperatorSplitting.jl does not implement the SciML solution interface, -# but we can obtain intermediate solutions via the iterator interface. +# OrdinaryDiffEqOperatorSplitting.jl implements only part of the SciML solution +# interface (see "Saving and interpolation" below); intermediate solutions are most +# directly obtained via the iterator interface. integrator = init(prob, alg, dt = 0.1) for (u, t) in TimeChoiceIterator(integrator, 0.0:0.5:1.0) @show t, u @@ -65,6 +66,92 @@ for (u, t) in TimeChoiceIterator(integrator, 0.0:0.5:1.0) end ``` +## Saving and interpolation + +`solve!` returns a solution whose `t`/`u` are filled as the integration proceeds, so +the usual saving keywords work and `sol(t)` interpolates between saved points: + +```julia +integrator = init(prob, alg; dt = 0.1, saveat = 0.25) +sol = solve!(integrator) + +sol.t # [0.0, 0.25, 0.5, 0.75, 1.0] +sol(0.3) # interpolated between the saved points +``` + +The supported keywords are `saveat`, `save_everystep` (default `false`), +`save_start`, `save_end` and `save_on`. They apply to the **outer** integrator only: +the inner splits are stages rather than steps, so their intermediate states do not +approximate the split solution at any time point. + +`saveat` points that fall strictly inside a step are filled from the step's +interpolant, so requesting output never changes the sequence of steps and therefore +never changes the splitting error. + +!!! note "The interpolant is first order" + A splitting step advances its children sequentially over staggered + subintervals, so their individual interpolants do not compose into an + approximation of the split solution. The dense output available at the outer + level is therefore linear between the step endpoints -- exact for the state a + `LieTrotterGodunov` step produces, but only first order for the second-order + schemes. If you need output that is accurate to the order of the scheme, pass + the times as `tstops` so that the integrator lands on them exactly: + + ```julia + sol = solve!(init(prob, alg; dt = 0.1, tstops = [0.25, 0.5], saveat = [0.25, 0.5])) + ``` + + `dense = true` is rejected: with a linear interpolant, `sol(t)` already + reproduces the integrator's own dense output from the saved points, so storing + per-step interpolation data would buy nothing. + +## Callbacks and events + +The standard SciML callbacks work through the `callback` keyword: + +```julia +# Stop as soon as the first component drops below a threshold. +cb = ContinuousCallback((u, t, integrator) -> u[1] - 0.5, terminate!) +sol = solve!(init(prob, alg; dt = 0.1, callback = cb)) +``` + +`DiscreteCallback`, `ContinuousCallback`, `VectorContinuousCallback` and +`CallbackSet` are all supported, along with `save_positions`, `affect_neg!`, +`rootfind`, `terminate!` and `derivative_discontinuity!`. Because the machinery is +DiffEqBase's own, callbacks built on top of it (for instance those in +DiffEqCallbacks.jl) work too. + +!!! note "Callbacks run on the outer integrator only" + A condition is evaluated once per **outer** step, after all the operators of + that step have been applied. It is never evaluated between two inner splits: a + splitting step applies its operators sequentially over staggered subintervals, + so the intermediate states are stages and do not approximate the solution of the + split system at any time point. There is correspondingly nothing meaningful for + a condition to test or an `affect!` to modify at that level. + + A consequence worth knowing: an `affect!` that modifies `integrator.u` is + propagated into every subintegrator before the next step, so modifying the state + from a callback is safe. + +### Accuracy of continuous events + +Event times are found by root-finding on the step's interpolant, which is linear +(see "Saving and interpolation" above). Two consequences: + + - The located event time is second order accurate in the step size, and is the + *exact* root of the linear interpolant over the step that brackets it. With + large steps -- adaptive splittings can grow the step considerably on smooth + problems -- the event time degrades accordingly. Cap it with `dtmax` when event + accuracy matters. + - An event that occurs and reverses **within a single step** cannot be detected, + because a linear interpolant has no interior extremum. Raising `interp_points` + does not help for the same reason, so setting `interp_points = 0` on the + callback avoids a sweep that cannot find anything the endpoints missed. + +Once the event time is located, the state there comes from the same interpolant and +the whole subintegrator tree is re-anchored to it, so integration resumes +consistently from the event. + ## Configuring individual subintegrators `init` takes one value per keyword, which is not enough when the operators want diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index 1cb9578..765b9d8 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -42,6 +42,10 @@ abstract type AbstractOperatorSplittingAlgorithm end abstract type AbstractOperatorSplittingCache end @inline SciMLBase.isadaptive(::AbstractOperatorSplittingAlgorithm) = false +# `AbstractOperatorSplittingAlgorithm` is not a `SciMLBase.AbstractDEAlgorithm`, so it +# does not pick up the generic trait fallbacks. The continuous callback root finder +# asks for this one. +@inline SciMLBase.isdiscrete(::AbstractOperatorSplittingAlgorithm) = false @inline isdtchangeable(alg::AbstractOperatorSplittingAlgorithm) = all(isdtchangeable.(alg.inner_algs)) include("function.jl") diff --git a/src/integrator.jl b/src/integrator.jl index e6d5499..27e42de 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -22,6 +22,28 @@ Base.@kwdef mutable struct IntegratorOptions{tType, fType, vbType, F3, atolType, end +""" + SaveOptions + +The saving settings of the outermost [`OperatorSplittingIntegrator`](@ref). + +Saving is a property of the outer integrator only: the inner splits are stages +rather than steps, so there is no time point at which their state is a meaningful +approximation of the split solution. These settings are therefore consumed by +`__init` directly and never enter the [`ConfigTree`](@ref), which keeps them from +travelling down to the leaf integrators. + +`saveat` itself is not held here -- it lives in the integrator's `saveat` heap, +because it is consumed as the integration proceeds. +""" +struct SaveOptions + save_on::Bool + save_everystep::Bool + save_start::Bool + save_end::Bool +end + + """ SplitSubIntegratorStatus @@ -103,7 +125,8 @@ mutable struct SplitSubIntegrator{ controller_cache::controllerType force_stepfail::Bool last_step_failed::Bool - u_modified::Bool # TODO we can probably remove this + # Named after the field DiffEqBase's callback machinery reads directly. + derivative_discontinuity::Bool status::SplitSubIntegratorStatus stats::IntegratorStats cache::cacheType @@ -147,6 +170,8 @@ mutable struct OperatorSplittingIntegrator{ tstopsType, saveatType, callbackType, + callbackCacheType, + eventErrType, cacheType, solType, subintTreeType, @@ -172,14 +197,31 @@ mutable struct OperatorSplittingIntegrator{ saveat::heapType _saveat::saveatType callback::callbackType + # Scratch buffers a `VectorContinuousCallback` needs; `nothing` when the callback + # set contains none. + const callback_cache::callbackCacheType + # Continuous-callback bookkeeping. DiffEqBase's root finder reads + # `event_last_time`/`last_event_error` to avoid re-detecting the event it just + # handled, and `handle_callbacks!` is responsible for writing them. + event_last_time::Int + vector_event_last_time::Int + last_event_error::eventErrType advance_to_tstop::Bool last_step_failed::Bool force_stepfail::Bool isout::Bool - u_modified::Bool + # DiffEqBase's callback machinery reads this as a literal field (not through + # `derivative_discontinuity!`), so the name has to match exactly. + derivative_discontinuity::Bool just_hit_tstop::Bool cache::cacheType sol::solType + # Saving state. `saveiter` is the authoritative length of the saved prefix of + # `sol.t`/`sol.u`: `savevalues!` writes index `saveiter` and `postamble!` + # truncates both vectors to it. + const save_opts::SaveOptions + saveiter::Int + postamble_done::Bool # Tuple of SplitSubIntegrator nodes (one per top-level operator). child_subintegrators::subintTreeType child_solution_indices::childSolidxType # Tuple @@ -229,7 +271,11 @@ function SciMLBase.__init( tstops = (), saveat = (), d_discontinuities = (), + save_on = true, save_everystep = false, + save_start = true, + save_end = true, + dense = false, callback = nothing, advance_to_tstop = false, adaptive = nothing, @@ -241,6 +287,19 @@ function SciMLBase.__init( (; u0, p) = prob t0, tf = prob.tspan + # Dense output would need per-step interpolation data on top of the saved + # points. The interpolant of a splitting step is linear (see the interpolation + # section below), so `sol(t)` already reproduces it exactly from `sol.t`/`sol.u` + # via SciMLBase's `LinearInterpolation` and storing `k` would buy nothing. + dense && throw( + ArgumentError( + "dense output is not supported by operator splitting integrators. The \ + interpolant is linear, so `sol(t)` between saved points is already exact \ + for it; use `saveat` or `save_everystep = true` to control the saved points." + ) + ) + save_opts = SaveOptions(save_on, save_everystep, save_start, save_end) + # By default every node adapts exactly if its own algorithm does; a scalar or a # TreeOption overrides the whole tree explicitly. if adaptive === nothing @@ -280,9 +339,27 @@ function SciMLBase.__init( tmp = setup_u(prob, alg, false) uType = typeof(u) - sol = SciMLBase.build_solution(prob, alg, tType[], uType[]) + # `build_solution` defaults `interp` to a `LinearInterpolation` over the very + # vectors passed in here, so `sol(t)` sees every later `savevalues!` push and + # agrees exactly with the integrator's own (linear) interpolant. `stats` has to + # be a real `DEStats`: the callback machinery increments `sol.stats.ncondition`. + sol = SciMLBase.build_solution( + prob, alg, tType[], uType[]; + stats = SciMLBase.DEStats(0), + calculate_error = false, + ) callback = DiffEqBase.CallbackSet(callback) + # Only a VectorContinuousCallback needs the scratch buffers, and they have to be + # wide enough for the widest one in the set. + max_len_cb = DiffEqBase.max_vector_callback_length_int(callback) + eventErrType = real(eltype(u)) + callback_cache = if max_len_cb === nothing + nothing + else + DiffEqBase.CallbackCache(u, max_len_cb, eventErrType, eventErrType) + end + child_subintegrators = build_subintegrators( prob, alg, uprev, u, @@ -315,9 +392,12 @@ function SciMLBase.__init( tstops_internal, tstops, saveat_internal, saveat, callback, + callback_cache, + 0, 1, zero(eventErrType), # event_last_time, vector_event_last_time, last_event_error advance_to_tstop, false, false, false, false, false, cache, sol, + save_opts, 0, false, child_subintegrators, child_solution_indices, child_synchronizers, @@ -330,10 +410,69 @@ function SciMLBase.__init( false, config, ) - DiffEqBase.initialize!(callback, u0, t0, integrator) + # The initial save happens *after* the callbacks are initialized, so that an + # initializer which modifies `u` is reflected in the first saved point. + initialize_callbacks!(integrator) + save_initial_value!(integrator) return integrator end +""" + initialize_callbacks!(integrator::OperatorSplittingIntegrator) + +Run the callbacks' `initialize` hooks, mirroring OrdinaryDiffEqCore's +`initialize_callbacks!`. + +The flag starts out set so that an `initialize` hook can clear it -- the default +initializer does exactly that -- and a hook that did modify `u` has its change +propagated to `uprev` and to the subintegrator tree before the first step. +""" +function initialize_callbacks!(integrator::OperatorSplittingIntegrator) + integrator.derivative_discontinuity = true + # Pass the integrator's own `u`, not the problem's `u0`: with `alias_u0 = false` + # they are different arrays and a hook has to be able to modify the live state. + modified = DiffEqBase.initialize!( + integrator.callback, integrator.u, integrator.t, integrator + ) + if modified + update_uprev!(integrator) + resync_children_after_modification!(integrator) + end + integrator.derivative_discontinuity = false + return nothing +end + +# A callback that changed `u` invalidates every child's copy of its slice. +# `rollback_children!` refills the children from the parent's `u`, re-anchors their +# clocks and marks their state modified so leaf integrators drop FSAL data. It also +# restores the `parent.u == child.u` invariant that lets the palindromic schemes skip +# the forward sync of their first child (`next_sync_is_continuous`). +function resync_children_after_modification!(integrator::OperatorSplittingIntegrator) + rollback_children!(integrator) + return nothing +end + +function SciMLBase.reeval_internals_due_to_modification!( + integrator::OperatorSplittingIntegrator, continuous_modification = true; + callback_initializealg = nothing + ) + resync_children_after_modification!(integrator) + integrator.derivative_discontinuity = false + return nothing +end + +# `save_start` puts u(t0) at index 1, so the first `savevalues!` writes index 2. +function save_initial_value!(integrator::OperatorSplittingIntegrator) + if integrator.save_opts.save_on && integrator.save_opts.save_start + integrator.saveiter = 1 + RecursiveArrayTools.copyat_or_push!(integrator.sol.t, 1, integrator.t) + RecursiveArrayTools.copyat_or_push!(integrator.sol.u, 1, integrator.u) + else + integrator.saveiter = 0 + end + return nothing +end + # --------------------------------------------------------------------------- # reinit! # --------------------------------------------------------------------------- @@ -379,12 +518,9 @@ function DiffEqBase.reinit!( resize!(integrator.sol.t, 0) resize!(integrator.sol.u, 0) end - if reinit_callbacks - DiffEqBase.initialize!(integrator.callback, u0, t0, integrator) - else - saving_callback = integrator.callback.discrete_callbacks[end] - DiffEqBase.initialize!(saving_callback, u0, t0, integrator) - end + integrator.event_last_time = 0 + integrator.vector_event_last_time = 1 + integrator.last_event_error = zero(integrator.last_event_error) if reinit_retcode integrator.sol = SciMLBase.solution_new_retcode( integrator.sol, ReturnCode.Default @@ -400,6 +536,23 @@ function DiffEqBase.reinit!( erase_sol, tstops, saveat, reinit_callbacks, reinit_retcode ) + # A leaf's `reinit!` restores it to the `u0` slice captured when it was built, not + # to the `u0` handed to this call, so the new state has to be pushed down + # explicitly. Relying on the forward sync of the next step is not enough: the + # palindromic schemes skip it for their first child. + resync_children_after_modification!(integrator) + + # After the children, so that an initializer which modifies `u` (and therefore + # resyncs the tree) is not undone by the child reinit, and so that the first saved + # point is the state the initializers left behind. + if reinit_callbacks + initialize_callbacks!(integrator) + end + # Saving is built into the integrator rather than provided by a saving callback, + # so the initial save is redone regardless of `reinit_callbacks`. Stale entries + # past the new `saveiter` are truncated by `postamble!`. + integrator.postamble_done = false + save_initial_value!(integrator) return nothing end @@ -649,7 +802,7 @@ function step_header!(integrator::AnySplitIntegrator) else reject_step!(integrator) end - elseif integrator.u_modified + elseif integrator.derivative_discontinuity update_uprev!(integrator) end increment_iteration(integrator) @@ -682,8 +835,11 @@ is_first_iteration(integrator::AnySplitIntegrator) = integrator.iter == 0 increment_iteration(integrator::AnySplitIntegrator) = integrator.iter += 1 function footer_reset_flags!(integrator) - integrator.u_modified = false + integrator.derivative_discontinuity = false integrator.just_hit_tstop = false + # Re-arm the postamble: stepping on after a `solve!`/`done` (e.g. after pushing a + # further tstop) has to be able to close out the solution again. + integrator.postamble_done = false return end footer_reset_flags!(::SplitSubIntegrator) = nothing @@ -692,14 +848,6 @@ function setup_validity_flags!(integrator, t_next) return end setup_validity_flags!(::SplitSubIntegrator, _) = nothing -function fix_solution_buffer_sizes!(integrator, sol) - resize!(integrator.sol.t, integrator.saveiter) - resize!(integrator.sol.u, integrator.saveiter) - if !(integrator.sol isa SciMLBase.DAESolution) - resize!(integrator.sol.k, integrator.saveiter_dense) - end - 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 @@ -753,6 +901,8 @@ function step_footer!(integrator::AnySplitIntegrator) try_snap_children_to_tstop!.(integrator.child_subintegrators, integrator.t) step_accept_controller!(integrator) validate_time_point(integrator) + # Callbacks run here (outer integrator only) and also perform the saving. + handle_callbacks!(integrator) elseif integrator.force_stepfail # Failure escalation protocol: the failing node's own adaptivity decides. fatal_rc = _fatal_child_retcode(integrator.child_subintegrators) @@ -928,6 +1078,24 @@ function SciMLBase.check_error!(integrator::SplitSubIntegrator) return code end +# SciMLBase's generic `check_error!` finalizes the integrator whenever the code is +# anything but `Success`, but a mid-solve node legitimately reports `Default`: the +# stepping loops call this before every step, so the generic version would run the +# postamble -- callback finalizers, the endpoint save -- once per step. Only a real +# failure finalizes here. +function SciMLBase.check_error!(integrator::OperatorSplittingIntegrator) + code = SciMLBase.check_error(integrator) + # Rebuilding the solution allocates, and this runs before every step, so only + # touch it when the code actually changed. + if integrator.sol.retcode !== code + integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, code) + end + if code ∉ (ReturnCode.Success, ReturnCode.Default) + SciMLBase.postamble!(integrator) + end + return code +end + @unroll function _check_error_children(current_retcode, children::Tuple) @unroll for child in children rc = _child_retcode(child) @@ -949,15 +1117,303 @@ end return (integrator.tmp,) end -function linear_interpolation!(y, t, y1, y2, t1, t2) - return y .= y1 + (t - t1) * (y2 - y1) / (t2 - t1) +# --------------------------------------------------------------------------- +# Interpolation +# +# A splitting step advances its children sequentially over staggered subintervals, +# so the children's own interpolants do not compose into an approximation of the +# split solution: the only dense output well defined at this level is the one built +# from the endpoints the outer integrator owns. The generic fallback is therefore +# linear, which is exact for the state a `LieTrotterGodunov` step produces and first +# order for the second order schemes. An algorithm that can do better implements +# `splitting_interpolant`/`splitting_interpolant!` for its own cache type; nothing +# else in the saving or callback machinery has to change. +# --------------------------------------------------------------------------- + +""" + splitting_interpolant(integrator, cache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}}) + splitting_interpolant!(out, integrator, cache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}}) + +Dense output of a single operator splitting step, evaluated in the step-local +coordinate `Θ = (t - tprev) / dt`, where `y₀`/`y₁` are the state at `tprev`/`t` and +`dt` is the step that was taken. `D` is the requested derivative order and `idxs` +selects components (`nothing` for all of them). + +This is the extension point for higher order dense output: dispatch on the +algorithm's cache type. The fallback for any +[`AbstractOperatorSplittingCache`](@ref) is linear interpolation. +""" +splitting_interpolant, splitting_interpolant! + +splitting_interpolant( + integrator, ::AbstractOperatorSplittingCache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}} +) where {D} = SciMLBase.linear_interpolant(Θ, dt, y₀, y₁, idxs, Val{D}) + +function splitting_interpolant!( + out, integrator, ::AbstractOperatorSplittingCache, Θ, dt, y₀, y₁, idxs, + ::Type{Val{D}} + ) where {D} + SciMLBase.linear_interpolant!(out, Θ, dt, y₀, y₁, idxs, Val{D}) + return out +end + +# The step that was just taken. `integrator.dt` must not be used here: by the time a +# step is finished `step_accept_controller!` has already overwritten it with the +# *next* proposal, so it is not the length of the interval `uprev`/`u` span. +_step_dt(integrator::OperatorSplittingIntegrator) = integrator.t - integrator.tprev + +function _interp_theta(integrator::OperatorSplittingIntegrator, t) + dt = _step_dt(integrator) + # A zero-length interval only happens before the first step, where `uprev == u` + # and every Θ gives the same value. + return iszero(dt) ? zero(t / oneunit(dt)) : (t - integrator.tprev) / dt +end + +function (integrator::OperatorSplittingIntegrator)( + t::Number, ::Type{deriv} = Val{0}; idxs = nothing + ) where {deriv} + return splitting_interpolant( + integrator, integrator.cache, + _interp_theta(integrator, t), _step_dt(integrator), + integrator.uprev, integrator.u, idxs, deriv + ) end -function (integrator::OperatorSplittingIntegrator)(tmp, t) - return linear_interpolation!( - tmp, t, integrator.uprev, integrator.u, integrator.tprev, integrator.t + +function (integrator::OperatorSplittingIntegrator)( + val::AbstractArray, t::Number, ::Type{deriv} = Val{0}; idxs = nothing + ) where {deriv} + return splitting_interpolant!( + val, integrator, integrator.cache, + _interp_theta(integrator, t), _step_dt(integrator), + integrator.uprev, integrator.u, idxs, deriv ) end +# --------------------------------------------------------------------------- +# change_t_via_interpolation! +# --------------------------------------------------------------------------- + +# DiffEqBase's continuous callback path passes the literal `Val{:false}` -- a `Val` of +# the *Symbol*, not of `false` -- so the flag has to be decoded by dispatch. Testing +# `if T` the way OrdinaryDiffEqCore does would throw on a Symbol. +_modify_save_endpoint(::Type{Val{true}}) = true +_modify_save_endpoint(::Type{Val{T}}) where {T} = false + +""" + change_t_via_interpolation!(integrator::OperatorSplittingIntegrator, t, modify_save_endpoint = Val{false}, reinitialize_alg = nothing) + +Move the integrator back to a time `t` inside the step that was just taken. + +`u` is refilled from the step's interpolant and the whole subintegrator tree is +re-anchored to `t` through [`rollback_children!`](@ref), so every child's state and +clock stay consistent with the parent (`validate_time_point` asserts the latter). + +`integrator.dt` is deliberately left alone: it holds the step size proposed for the +next step, and the interpolation coordinate is derived from `t - tprev` rather than +from `dt`. +""" +function SciMLBase.change_t_via_interpolation!( + integrator::OperatorSplittingIntegrator, t, + modify_save_endpoint::Type{Val{T}} = Val{false}, + reinitialize_alg = nothing + ) where {T} + if integrator.tdir * t < integrator.tdir * integrator.tprev + error("Current interpolant only works between tprev and t") + elseif t != integrator.t + integrator(integrator.u, t) + integrator.t = t + # Push the interpolated slices into the children and move their clocks back. + rollback_children!(integrator) + if _modify_save_endpoint(modify_save_endpoint) + solution_endpoint_match_cur_integrator!(integrator) + end + end + return nothing +end + +# --------------------------------------------------------------------------- +# Saving +# +# Only the outermost integrator saves: the inner splits are stages rather than +# steps, so their intermediate states are not approximations of the split solution +# at any time point. +# +# `saveiter` is the authoritative length of the saved prefix of `sol.t`/`sol.u`. +# `sol.interp` is a `LinearInterpolation` aliasing those same vectors, so `sol(t)` +# reproduces this integrator's own interpolant between saved points. +# --------------------------------------------------------------------------- + +_saved_at_current_t(integrator::OperatorSplittingIntegrator) = + integrator.saveiter > 0 && + integrator.sol.t[integrator.saveiter] == integrator.t + +function _save_current!(integrator::OperatorSplittingIntegrator) + integrator.saveiter += 1 + RecursiveArrayTools.copyat_or_push!( + integrator.sol.t, integrator.saveiter, integrator.t + ) + RecursiveArrayTools.copyat_or_push!( + integrator.sol.u, integrator.saveiter, integrator.u + ) + return nothing +end + +""" + savevalues!(integrator::OperatorSplittingIntegrator, force_save = false, reduce_size = true) + +Save the state at the current time point and at any pending `saveat` points that the +step just taken has passed, returning `(saved, savedexactly)`. + +`saveat` points strictly inside the step are filled from the step's interpolant +(see [`splitting_interpolant`](@ref)), so asking for output never changes the +sequence of steps and therefore never changes the splitting error. Because the +generic interpolant is linear, interpolated points are first order accurate even +when the splitting scheme is second order; pass the times as `tstops` instead to +have the integrator land on them exactly. +""" +function SciMLBase.savevalues!( + integrator::OperatorSplittingIntegrator, + force_save = false, reduce_size = true + )::Tuple{Bool, Bool} + saved = false + savedexactly = false + integrator.save_opts.save_on || return saved, savedexactly + (; save_everystep, save_end) = integrator.save_opts + tf = integrator.sol.prob.tspan[2] + + # The heaps store raw times and carry the direction in their ordering, so the + # comparison is scaled by tdir rather than the stored value. + tdir_t = integrator.tdir * integrator.t + while !isempty(integrator.saveat) && + integrator.tdir * first(integrator.saveat) <= tdir_t + curt = pop!(integrator.saveat) + if curt == integrator.t + # `save_end` owns the final point; leave it to the postamble. + (!save_end && curt == tf) && continue + saved = true + savedexactly = true + _save_current!(integrator) + else + saved = true + integrator.saveiter += 1 + RecursiveArrayTools.copyat_or_push!( + integrator.sol.t, integrator.saveiter, curt + ) + # Freshly allocated by the interpolant, so hand over ownership. + RecursiveArrayTools.copyat_or_push!( + integrator.sol.u, integrator.saveiter, integrator(curt), false + ) + end + end + + if force_save || ( + save_everystep && !_saved_at_current_t(integrator) && + (save_end || integrator.t != tf) + ) + saved = true + savedexactly = true + _save_current!(integrator) + end + + return saved, savedexactly +end + +# Make sure the final time point is in the solution. Mirrors OrdinaryDiffEqCore's +# `solution_endpoint_match_cur_integrator!`. +function solution_endpoint_match_cur_integrator!(integrator::OperatorSplittingIntegrator) + integrator.save_opts.save_on || return nothing + integrator.save_opts.save_end || return nothing + _saved_at_current_t(integrator) && return nothing + _save_current!(integrator) + return nothing +end + +# --------------------------------------------------------------------------- +# Callbacks +# +# Callbacks run on the outermost integrator only. A condition evaluated between two +# inner splits would see a state that is not an approximation of the split solution +# at any time point -- the splits are stages, not steps -- so there is nothing +# meaningful for a condition or an `affect!` to act on down there. +# +# The heavy lifting (root finding, `save_positions`, `affect!` dispatch) is +# DiffEqBase's; this is the per-step driver, modelled on OrdinaryDiffEqCore's +# `handle_callbacks!`. It also owns the plain `saveat` save when no callback saved. +# --------------------------------------------------------------------------- + +# Type-stable dispatch into `DiffEqBase.apply_callback!`: the explicit +# `if (return) else (rest)` chain compiles to a switch over a heterogeneous tuple. +@generated function apply_ith_callback!( + integrator, + time, upcrossing, event_idx, cb_idx, + callbacks::NTuple{ + N, + Union{DiffEqBase.ContinuousCallback, DiffEqBase.VectorContinuousCallback}, + } + ) where {N} + ex = quote + throw(BoundsError(callbacks, cb_idx)) + end + for i in 1:N + ex = quote + if (cb_idx == $i) + return DiffEqBase.apply_callback!( + integrator, callbacks[$i], time, upcrossing, event_idx + ) + else + $ex + end + end + end + return ex +end + +handle_callbacks!(::SplitSubIntegrator) = nothing + +function handle_callbacks!(integrator::OperatorSplittingIntegrator) + discrete_callbacks = integrator.callback.discrete_callbacks + continuous_callbacks = integrator.callback.continuous_callbacks + + continuous_modified = false + discrete_modified = false + saved_in_cb = false + + if !(continuous_callbacks isa Tuple{}) + time, upcrossing, event_occurred, event_idx, idx, _counter = + DiffEqBase.find_first_continuous_callback( + integrator, continuous_callbacks... + ) + if event_occurred + integrator.event_last_time = idx + integrator.vector_event_last_time = event_idx + continuous_modified, saved_in_cb = apply_ith_callback!( + integrator, time, upcrossing, event_idx, idx, continuous_callbacks + ) + # The step was cut at the event and the state may have jumped, so the + # controller's error history no longer describes what happens next. + reinit_node_controller!(integrator) + else + # Clearing these is what stops the root finder from nudging `tprev` on a + # step that had no event. + integrator.event_last_time = 0 + integrator.vector_event_last_time = 1 + end + end + + if !integrator.force_stepfail && !(discrete_callbacks isa Tuple{}) + discrete_modified, saved_in_cb = DiffEqBase.apply_discrete_callback!( + integrator, discrete_callbacks... + ) + end + + if !saved_in_cb + SciMLBase.savevalues!(integrator) + end + + integrator.derivative_discontinuity = continuous_modified | discrete_modified + return nothing +end + # --------------------------------------------------------------------------- # Step size control # @@ -1116,7 +1572,17 @@ function SciMLBase.done(integrator::OperatorSplittingIntegrator) end function SciMLBase.postamble!(integrator::OperatorSplittingIntegrator) - return DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator) + # `solve!`, `done` and `check_error!` all call this; only the first one after a + # step may run the finalizers and close out the solution. + integrator.postamble_done && return nothing + integrator.postamble_done = true + DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator) + solution_endpoint_match_cur_integrator!(integrator) + # `saveiter` is authoritative: drop whatever a previous, longer run left behind + # (a `reinit!` without `erase_sol`). + resize!(integrator.sol.t, integrator.saveiter) + resize!(integrator.sol.u, integrator.saveiter) + return nothing end function __step!(integrator::AnySplitIntegrator) @@ -1325,7 +1791,7 @@ function _build_child( 0, 0, # iter, success_iter EEst_val, controller_cache, - false, false, false, # force_stepfail, last_step_failed, u_modified + false, false, false, # force_stepfail, last_step_failed, derivative_discontinuity SplitSubIntegratorStatus(), IntegratorStats(), level_cache, @@ -1412,6 +1878,36 @@ SciMLBase.first_tstop(i::AnySplitIntegrator) = first(i.tstops) SciMLBase.pop_tstop!(i::AnySplitIntegrator) = pop!(i.tstops) DiffEqBase.get_dt(i::AnySplitIntegrator) = i.dt + +# Continuous callbacks relax the step size through this after an event. +# `dtcache` mirrors the standing proposal as a magnitude (see `__init`). +function SciMLBase.set_proposed_dt!(integrator::OperatorSplittingIntegrator, dt) + if integrator.dtcache != abs(dt) + integrator.dtcache = abs(dt) + if !isadaptive(integrator) + integrator.dt = dt + end + end + return nothing +end + +""" + terminate!(integrator::OperatorSplittingIntegrator, retcode = ReturnCode.Terminated) + +Stop the integration from inside a callback's `affect!`. + +Emptying the outermost `tstops` heap is what ends the `solve!` loop and makes +`done` report completion; the children's heaps are irrelevant once the outer loop +exits, and `postamble!` still closes out the solution. +""" +function SciMLBase.terminate!( + integrator::OperatorSplittingIntegrator, retcode = ReturnCode.Terminated + ) + _set_retcode!(integrator, retcode) + empty!(integrator.tstops) + return nothing +end + function set_dt!(i::DEIntegrator, dt) iszero(dt) && error("dt must be nonzero") return i.dt = dt @@ -1434,12 +1930,13 @@ function DiffEqBase.add_saveat!(i::OperatorSplittingIntegrator, t) return nothing end -# SciMLBase v3 renamed `u_modified!` → `derivative_discontinuity!`. +# SciMLBase v3 renamed `u_modified!` → `derivative_discontinuity!`. The flag itself +# is stored in the `derivative_discontinuity` field, which is the name DiffEqBase's +# callback machinery reads directly. @static if isdefined(DiffEqBase, :u_modified!) - DiffEqBase.u_modified!(i::OperatorSplittingIntegrator, bool) = i.u_modified = bool - DiffEqBase.u_modified!(i::SplitSubIntegrator, bool) = i.u_modified = bool + DiffEqBase.u_modified!(i::AnySplitIntegrator, bool) = i.derivative_discontinuity = bool end @static if isdefined(SciMLBase, :derivative_discontinuity!) - SciMLBase.derivative_discontinuity!(i::OperatorSplittingIntegrator, bool) = i.u_modified = bool - SciMLBase.derivative_discontinuity!(i::SplitSubIntegrator, bool) = i.u_modified = bool + SciMLBase.derivative_discontinuity!(i::AnySplitIntegrator, bool) = + i.derivative_discontinuity = bool end diff --git a/src/utils.jl b/src/utils.jl index 642fd79..a635f2a 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -7,17 +7,22 @@ function tstops_and_saveat_heaps(t0, tf, tstops, saveat) tstops = [filter(t -> t0 < t < tf || tf < t < t0, tstops)..., tf] tstops = BinaryHeaps.BinaryHeap{FT, ordering}(tstops) + # `t0` is excluded and `tf` included, matching OrdinaryDiffEqCore's + # `initialize_saveat`: the initial point is owned by `save_start` and the final + # one by `save_end`, so leaving them in the heap would duplicate them. + tdir = tf > t0 ? one(FT) : -one(FT) if isnothing(saveat) - saveat = [t0, tf] + saveat = FT[] elseif saveat isa Number saveat > zero(saveat) || error("saveat value must be positive") - saveat = tf > t0 ? saveat : -saveat - saveat = [t0:saveat:tf..., tf] + step = tdir * saveat + # `tf` is not appended: the range hits it when it divides the span evenly, + # and `save_end` owns the final point otherwise. + saveat = collect((t0 + step):step:tf) else - # We do not need to filter saveat like tstops because the saving - # callback will ignore any times that are not between t0 and tf. - saveat = collect(saveat) + saveat = collect(FT, saveat) end + saveat = filter(t -> tdir * t0 < tdir * t <= tdir * tf, saveat) saveat = BinaryHeaps.BinaryHeap{FT, ordering}(saveat) return tstops, saveat diff --git a/test/callbacks.jl b/test/callbacks.jl new file mode 100644 index 0000000..6db104c --- /dev/null +++ b/test/callbacks.jl @@ -0,0 +1,516 @@ +using OrdinaryDiffEqOperatorSplitting +using Test + +import SciMLBase +import SciMLBase: ReturnCode +import DiffEqBase: DiffEqBase, ODEFunction, DiscreteCallback, ContinuousCallback, + VectorContinuousCallback, CallbackSet +using OrdinaryDiffEqLowOrderRK +using OrdinaryDiffEqTsit5 + +const OS = OrdinaryDiffEqOperatorSplitting + +# --------------------------------------------------------------------------- +# Reference problem. The third component is decoupled and solved by the first +# operator alone, so `u[3](t) = 3exp(-t/10)` exactly -- which gives continuous +# callback tests an analytic event time to compare against. +# --------------------------------------------------------------------------- +const U0 = [1.0, 2.0, 3.0] +const TSPAN = (0.0, 1.0) + +f1 = ODEFunction((du, u, p, t) -> (@. du = -0.1u)) +f2 = ODEFunction( + function (du, u, p, t) + du[1] = -0.01u[2] + du[2] = -0.01u[1] + return nothing + end +) +fsplit = GenericSplitFunction((f1, f2), ([1, 2, 3], [1, 2])) + +make_prob(tspan = TSPAN) = OperatorSplittingProblem(fsplit, copy(U0), tspan) +ltg() = LieTrotterGodunov((Euler(), Euler())) +# Accurate inner solvers, so that an event time's error is the interpolant's alone. +ltg_exact() = LieTrotterGodunov((Tsit5(), Tsit5())) + +# Exact crossing time of u[3] = level. +exact_crossing(level) = -10 * log(level / 3) + +@testset "discrete callbacks" begin + prob = make_prob() + + @testset "condition is evaluated once per outer step" begin + calls = Ref(0) + cb = DiscreteCallback((u, t, integrator) -> (calls[] += 1; false), integrator -> nothing) + DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)) + @test calls[] == 10 + end + + @testset "nested splitting does not run callbacks per inner stage" begin + inner = GenericSplitFunction((f1, f1), ([1, 2], [1, 2])) + outer = GenericSplitFunction((f1, inner), ([1, 2, 3], [1, 2])) + nprob = OperatorSplittingProblem(outer, copy(U0), TSPAN) + alg = LieTrotterGodunov((Euler(), LieTrotterGodunov((Euler(), Euler())))) + calls = Ref(0) + cb = DiscreteCallback((u, t, integrator) -> (calls[] += 1; false), integrator -> nothing) + DiffEqBase.solve!(DiffEqBase.init(nprob, alg; dt = 0.1, callback = cb)) + # Ten outer steps, regardless of how many substeps the inner node takes. + @test calls[] == 10 + end + + @testset "affect! fires and sees the integrator" begin + seen = Float64[] + # Threshold between grid points, compared with isapprox: see the CallbackSet + # testset below for why the nominal grid cannot be trusted exactly. + cb = DiscreteCallback( + (u, t, integrator) -> t >= 0.45, + integrator -> push!(seen, integrator.t) + ) + sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)) + @test seen ≈ [0.5, 0.6, 0.7, 0.8, 0.9, 1.0] + @test sol.retcode == ReturnCode.Success + end + + @testset "save_positions brackets the affect!" begin + cb = DiscreteCallback( + (u, t, integrator) -> isapprox(t, 0.5), + integrator -> (integrator.u[1] += 10.0) + ) + sol = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb, tstops = [0.5]) + ) + # Two points at the event time: the state before and after the affect!. + idxs = findall(t -> isapprox(t, 0.5), sol.t) + @test length(idxs) == 2 + @test sol.u[idxs[2]][1] - sol.u[idxs[1]][1] ≈ 10.0 + + # save_positions = (false, false) records nothing extra. + cb_quiet = DiscreteCallback( + (u, t, integrator) -> isapprox(t, 0.5), + integrator -> (integrator.u[1] += 10.0); + save_positions = (false, false) + ) + sol = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb_quiet, tstops = [0.5]) + ) + @test count(t -> isapprox(t, 0.5), sol.t) == 0 + end + + @testset "a modified u reaches the whole subintegrator tree" begin + # This is the operator-splitting-specific hazard: every child holds its own + # copy of its slice, and the palindromic schemes skip the forward sync of + # their first child on the assumption that it still matches the parent. + for alg in ( + ltg(), + StrangMarchuk((Euler(), Euler())), + PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())), + ) + cb = DiscreteCallback( + (u, t, integrator) -> isapprox(t, 0.5), + integrator -> (integrator.u .= [100.0, 200.0, 300.0]) + ) + integrator = DiffEqBase.init( + prob, alg; dt = 0.1, callback = cb, tstops = [0.5] + ) + while integrator.t < 0.5 + DiffEqBase.step!(integrator) + end + @test integrator.u == [100.0, 200.0, 300.0] + @test integrator.child_subintegrators[1].u == [100.0, 200.0, 300.0] + @test integrator.child_subintegrators[2].u == [100.0, 200.0] + OS.validate_time_point(integrator) + end + end + + @testset "a modified u actually changes the trajectory" begin + plain = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, tstops = [0.5]) + ) + cb = DiscreteCallback( + (u, t, integrator) -> isapprox(t, 0.5), + integrator -> (integrator.u .*= 2) + ) + bumped = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb, tstops = [0.5]) + ) + # Doubling halfway through must roughly double the endpoint of a linear ODE. + @test bumped.u[end] ≈ 2 .* plain.u[end] rtol = 1.0e-8 + end + + @testset "multiple callbacks in a CallbackSet" begin + # Thresholds sit between grid points so that ulp drift in the step times + # (0.1 * 8 is 0.7999999999999999, not 0.8) cannot flip a condition, and the + # recorded times are compared with isapprox for the same reason. + early = Float64[] + late = Float64[] + cb1 = DiscreteCallback((u, t, integrator) -> t >= 0.35, integrator -> push!(early, integrator.t)) + cb2 = DiscreteCallback((u, t, integrator) -> t >= 0.85, integrator -> push!(late, integrator.t)) + sol = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, callback = CallbackSet(cb1, cb2)) + ) + @test early ≈ [0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] + @test late ≈ [0.9, 1.0] + @test sol.retcode == ReturnCode.Success + end + + @testset "sol.stats.ncondition counts condition evaluations" begin + cb = DiscreteCallback((u, t, integrator) -> false, integrator -> nothing) + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb) + sol = DiffEqBase.solve!(integrator) + @test sol.stats.ncondition == 10 + end +end + +@testset "terminate!" begin + prob = make_prob() + + @testset "stops the solve and sets the retcode" begin + cb = DiscreteCallback( + (u, t, integrator) -> t >= 0.5, + integrator -> SciMLBase.terminate!(integrator) + ) + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb) + sol = DiffEqBase.solve!(integrator) + @test sol.retcode == ReturnCode.Terminated + @test sol.t[end] ≈ 0.5 + @test integrator.t ≈ 0.5 + @test SciMLBase.done(integrator) + end + + @testset "accepts an explicit retcode" begin + cb = DiscreteCallback( + (u, t, integrator) -> t >= 0.5, + integrator -> SciMLBase.terminate!(integrator, ReturnCode.Success) + ) + sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)) + @test sol.retcode == ReturnCode.Success + end +end + +@testset "continuous callbacks" begin + prob = make_prob() + level = 2.9 + exact = exact_crossing(level) + + @testset "locates the crossing and fires exactly once" begin + hits = Float64[] + cb = ContinuousCallback( + (u, t, integrator) -> u[3] - level, + integrator -> push!(hits, integrator.t) + ) + sol = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg_exact(); dt = 0.05, callback = cb) + ) + @test length(hits) == 1 + @test hits[1] ≈ exact atol = 1.0e-3 + @test hits[1] in sol.t + @test issorted(sol.t) + @test sol.retcode == ReturnCode.Success + end + + @testset "event time converges with the interpolant" begin + # The interpolant is linear, so the located root is second order accurate. + # (The step grid moves with dt, so the constant is noisy; a factor-of-four + # improvement over a factor-of-four refinement is the honest assertion.) + errs = map((0.1, 0.025)) do dt + hits = Float64[] + cb = ContinuousCallback( + (u, t, integrator) -> u[3] - level, + integrator -> push!(hits, integrator.t) + ) + DiffEqBase.solve!( + DiffEqBase.init(prob, ltg_exact(); dt = dt, callback = cb) + ) + abs(hits[1] - exact) + end + @test errs[2] < errs[1] / 4 + end + + @testset "affect_neg! handles the downcrossing" begin + ups = Ref(0) + downs = Ref(0) + cb = ContinuousCallback( + (u, t, integrator) -> u[3] - level, + integrator -> (ups[] += 1); + affect_neg! = integrator -> (downs[] += 1) + ) + DiffEqBase.solve!(DiffEqBase.init(prob, ltg_exact(); dt = 0.05, callback = cb)) + # u[3] decays through the level, so this is a downcrossing only. + @test ups[] == 0 + @test downs[] == 1 + end + + @testset "the tree is re-anchored to the event time" begin + cb = ContinuousCallback( + (u, t, integrator) -> u[3] - level, + integrator -> nothing + ) + integrator = DiffEqBase.init(prob, ltg_exact(); dt = 0.05, callback = cb) + while integrator.t < exact + DiffEqBase.step!(integrator) + end + # `change_t_via_interpolation!` moved the step back to the event; every child + # clock and buffer must have followed. + @test integrator.t ≈ exact atol = 1.0e-3 + OS.validate_time_point(integrator) + @test integrator.child_subintegrators[1].u ≈ integrator.u + @test integrator.child_subintegrators[2].u ≈ integrator.u[1:2] + end + + @testset "an affect! that modifies u reaches the children" begin + cb = ContinuousCallback( + (u, t, integrator) -> u[3] - level, + integrator -> (integrator.u .= [7.0, 8.0, 9.0]) + ) + integrator = DiffEqBase.init(prob, ltg_exact(); dt = 0.05, callback = cb) + while integrator.t < exact + DiffEqBase.step!(integrator) + end + @test integrator.u == [7.0, 8.0, 9.0] + @test integrator.child_subintegrators[1].u == [7.0, 8.0, 9.0] + @test integrator.child_subintegrators[2].u == [7.0, 8.0] + end + + @testset "NoRootFind fires at the step endpoint" begin + hits = Float64[] + cb = ContinuousCallback( + (u, t, integrator) -> u[3] - level, + integrator -> push!(hits, integrator.t); + rootfind = SciMLBase.NoRootFind + ) + DiffEqBase.solve!(DiffEqBase.init(prob, ltg_exact(); dt = 0.1, callback = cb)) + @test length(hits) == 1 + # No root finding: the event is reported at the end of the bracketing step. + @test hits[1] ≈ 0.4 + end + + @testset "works with an adaptive splitting" begin + hits = Float64[] + cb = ContinuousCallback( + (u, t, integrator) -> u[3] - level, + integrator -> push!(hits, integrator.t) + ) + # `dtmax` matters here: this problem is so nearly linear that the controller + # otherwise grows the step to span the rest of the tspan in one go, and the + # event can only ever be as accurate as the interpolant over the bracketing + # step -- which is linear. + sol = DiffEqBase.solve!( + DiffEqBase.init( + prob, PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())); + dt = 0.05, dtmax = 0.05, callback = cb + ) + ) + @test length(hits) == 1 + @test hits[1] ≈ exact atol = 1.0e-3 + @test sol.retcode == ReturnCode.Success + end + + @testset "a long adaptive step degrades the event time, not the bracket" begin + # Documented consequence of the linear interpolant: the located event is the + # exact root of the interpolant over whatever step brackets it. Pin that + # relationship so the behaviour is a stated property rather than a surprise. + bracket = Ref((0.0, 0.0)) + hits = Float64[] + cb = ContinuousCallback( + (u, t, integrator) -> u[3] - level, + function (integrator) + push!(hits, integrator.t) + return bracket[] = (integrator.tprev, integrator.t) + end + ) + DiffEqBase.solve!( + DiffEqBase.init( + prob, PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())); + dt = 0.05, callback = cb + ) + ) + t0 = bracket[][1] + # u[3] is solved exactly by the first operator alone. + u3(t) = 3 * exp(-0.1t) + # The step originally ended where the event time now sits plus the part of the + # step that was cut off; recover it from the interpolant relation instead. + @test length(hits) == 1 + @test hits[1] > t0 + # The event solves the *linear* interpolant, so g is zero there by + # construction while the true crossing is elsewhere. + @test !isapprox(hits[1], exact; atol = 1.0e-3) + @test u3(hits[1]) < level # the linear chord under-estimates a convex decay + end + + @testset "interp_points = 0 skips the safety sweep" begin + calls = Ref(0) + cb = ContinuousCallback( + (u, t, integrator) -> (calls[] += 1; u[3] - level), + integrator -> nothing; + interp_points = 0 + ) + DiffEqBase.solve!(DiffEqBase.init(prob, ltg_exact(); dt = 0.25, callback = cb)) + sweeping = Ref(0) + cb2 = ContinuousCallback( + (u, t, integrator) -> (sweeping[] += 1; u[3] - level), + integrator -> nothing + ) + DiffEqBase.solve!(DiffEqBase.init(prob, ltg_exact(); dt = 0.25, callback = cb2)) + # The default `interp_points = 11` samples the (collinear) interpolant many + # more times without being able to find anything the endpoints missed. + @test calls[] < sweeping[] + end +end + +@testset "VectorContinuousCallback" begin + prob = make_prob() + levels = (2.9, 2.85) + + # Current DiffEqBase hands `affect!` a mask of simultaneous events rather than a + # single index; accept either shape so this does not depend on that detail. + fired_index(idx::Integer) = Int(idx) + fired_index(mask) = only(findall(!iszero, mask)) + + fired = Tuple{Int, Float64}[] + cb = VectorContinuousCallback( + function (out, u, t, integrator) + out[1] = u[3] - levels[1] + out[2] = u[3] - levels[2] + return nothing + end, + (integrator, idx) -> push!(fired, (fired_index(idx), integrator.t)), + 2 + ) + integrator = DiffEqBase.init(prob, ltg_exact(); dt = 0.05, callback = cb) + @test integrator.callback_cache !== nothing + sol = DiffEqBase.solve!(integrator) + + @test length(fired) == 2 + @test first.(fired) == [1, 2] # the higher level is crossed first + @test last(fired[1]) ≈ exact_crossing(levels[1]) atol = 1.0e-3 + @test last(fired[2]) ≈ exact_crossing(levels[2]) atol = 1.0e-3 + @test sol.retcode == ReturnCode.Success + + @testset "no cache is allocated without a vector callback" begin + scalar_cb = ContinuousCallback((u, t, integrator) -> u[3] - 2.9, integrator -> nothing) + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = scalar_cb) + @test integrator.callback_cache === nothing + end +end + +@testset "callback initialization" begin + prob = make_prob() + + @testset "an initializer can add tstops" begin + # The pattern PresetTimeCallback uses. + hits = Float64[] + cb = DiscreteCallback( + (u, t, integrator) -> t in (0.35, 0.65), + integrator -> push!(hits, integrator.t); + initialize = function (c, u, t, integrator) + DiffEqBase.add_tstop!(integrator, 0.35) + DiffEqBase.add_tstop!(integrator, 0.65) + return SciMLBase.derivative_discontinuity!(integrator, false) + end + ) + sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)) + @test hits == [0.35, 0.65] + @test sol.retcode == ReturnCode.Success + end + + @testset "an initializer that modifies u is reflected everywhere" begin + cb = DiscreteCallback( + (u, t, integrator) -> false, + integrator -> nothing; + initialize = function (c, u, t, integrator) + integrator.u .= [5.0, 6.0, 7.0] + return SciMLBase.derivative_discontinuity!(integrator, true) + end + ) + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb) + @test integrator.u == [5.0, 6.0, 7.0] + @test integrator.uprev == [5.0, 6.0, 7.0] + @test integrator.child_subintegrators[1].u == [5.0, 6.0, 7.0] + @test integrator.child_subintegrators[2].u == [5.0, 6.0] + # The first saved point is the modified state, not the problem's u0. + @test integrator.sol.u[1] == [5.0, 6.0, 7.0] + end + + @testset "finalize runs exactly once" begin + finals = Ref(0) + cb = DiscreteCallback( + (u, t, integrator) -> false, + integrator -> nothing; + finalize = (c, u, t, integrator) -> (finals[] += 1) + ) + DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)) + @test finals[] == 1 + end + + @testset "reinit! keeps an initializer's modification" begin + # The child reinit must not undo what the initializer did, so the callbacks + # have to be initialized after the subintegrator tree is restored. + cb = DiscreteCallback( + (u, t, integrator) -> false, + integrator -> nothing; + initialize = function (c, u, t, integrator) + integrator.u .= [5.0, 6.0, 7.0] + return SciMLBase.derivative_discontinuity!(integrator, true) + end + ) + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb) + DiffEqBase.solve!(integrator) + DiffEqBase.reinit!(integrator) + @test integrator.u == [5.0, 6.0, 7.0] + @test integrator.uprev == [5.0, 6.0, 7.0] + @test integrator.child_subintegrators[1].u == [5.0, 6.0, 7.0] + @test integrator.child_subintegrators[2].u == [5.0, 6.0] + @test integrator.sol.u[1] == [5.0, 6.0, 7.0] + end + + @testset "reinit! re-runs the initializers" begin + inits = Ref(0) + cb = DiscreteCallback( + (u, t, integrator) -> false, + integrator -> nothing; + initialize = function (c, u, t, integrator) + inits[] += 1 + return SciMLBase.derivative_discontinuity!(integrator, false) + end + ) + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb) + @test inits[] == 1 + DiffEqBase.solve!(integrator) + DiffEqBase.reinit!(integrator) + @test inits[] == 2 + DiffEqBase.reinit!(integrator; reinit_callbacks = false) + @test inits[] == 2 + @test DiffEqBase.solve!(integrator).retcode == ReturnCode.Success + end +end + +@testset "callbacks combined with saving" begin + prob = make_prob() + + @testset "saveat and an event coexist" begin + cb = ContinuousCallback( + (u, t, integrator) -> u[3] - 2.9, + integrator -> nothing + ) + sol = DiffEqBase.solve!( + DiffEqBase.init( + prob, ltg_exact(); dt = 0.1, saveat = [0.2, 0.8], callback = cb + ) + ) + @test issorted(sol.t) + @test 0.2 in sol.t + @test 0.8 in sol.t + @test sol.t[1] == 0.0 + @test sol.t[end] == 1.0 + end + + @testset "a discrete callback still leaves a usable interpolation" begin + cb = DiscreteCallback((u, t, integrator) -> false, integrator -> nothing) + sol = DiffEqBase.solve!( + DiffEqBase.init( + prob, ltg(); dt = 0.1, callback = cb, save_everystep = true + ) + ) + @test sol(0.35) ≈ (sol.u[4] .+ sol.u[5]) ./ 2 + end +end diff --git a/test/qa/qa.jl b/test/qa/qa.jl index cd55cef..8de8428 100644 --- a/test/qa/qa.jl +++ b/test/qa/qa.jl @@ -32,6 +32,10 @@ run_qa( :get_EEst, :set_EEst!, :get_current_adaptive_order, # OrdinaryDiffEqCore :gamma_default, :failfactor_default, :AbstractControllerCache, # OrdinaryDiffEqCore :promote_tspan, # SciMLBase + # Shared linear interpolation kernels, reused so that the + # integrator's own dense output and `sol(t)` (which goes through + # SciMLBase's `LinearInterpolation`) cannot drift apart. + :linear_interpolant, :linear_interpolant!, # SciMLBase # Broadcast extension points a TreeOption implements (src/config_tree.jl). :Broadcasted, :broadcastable, :dotview, :materialize!, # Base ), diff --git a/test/saving.jl b/test/saving.jl new file mode 100644 index 0000000..8fc8435 --- /dev/null +++ b/test/saving.jl @@ -0,0 +1,439 @@ +using OrdinaryDiffEqOperatorSplitting +using Test + +import SciMLBase +import SciMLBase: ReturnCode +import DiffEqBase: DiffEqBase, ODEFunction, DiscreteCallback, CallbackSet +using OrdinaryDiffEqLowOrderRK +using OrdinaryDiffEqTsit5 + +const OS = OrdinaryDiffEqOperatorSplitting + +# --------------------------------------------------------------------------- +# Reference problem: two linear operators sharing the first two components. +# --------------------------------------------------------------------------- +const U0 = [1.0, 2.0, 3.0] +const TSPAN = (0.0, 1.0) + +f1 = ODEFunction((du, u, p, t) -> (@. du = -0.1u)) +f2 = ODEFunction( + function (du, u, p, t) + du[1] = -0.01u[2] + du[2] = -0.01u[1] + return nothing + end +) +fsplit = GenericSplitFunction((f1, f2), ([1, 2, 3], [1, 2])) + +make_prob(tspan = TSPAN) = OperatorSplittingProblem(fsplit, copy(U0), tspan) +ltg() = LieTrotterGodunov((Euler(), Euler())) + +@testset "saved time points" begin + prob = make_prob() + + @testset "default saves only the interval endpoints" begin + # `save_everystep` defaults to false, so only `save_start`/`save_end` fire. + sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1)) + @test sol.t == [0.0, 1.0] + @test length(sol.u) == 2 + @test sol.u[1] == U0 + @test sol.retcode == ReturnCode.Success + end + + @testset "save_everystep saves every accepted step" begin + sol = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, save_everystep = true) + ) + @test length(sol.t) == 11 # t0 plus ten steps of dt = 0.1 + @test sol.t[1] == 0.0 + @test sol.t[end] == 1.0 + @test issorted(sol.t) + # No duplicate of the final point from `save_end`. + @test allunique(sol.t) + end + + @testset "saveat vector" begin + sol = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = [0.25, 0.5]) + ) + @test sol.t == [0.0, 0.25, 0.5, 1.0] + end + + @testset "saveat number excludes t0 and does not duplicate tf" begin + sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = 0.25)) + @test sol.t == [0.0, 0.25, 0.5, 0.75, 1.0] + # A step that does not divide the span evenly leaves tf to `save_end`. + sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = 0.3)) + @test sol.t ≈ [0.0, 0.3, 0.6, 0.9, 1.0] + end + + @testset "saveat points outside the tspan are dropped" begin + sol = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = [-1.0, 0.5, 7.0]) + ) + @test sol.t == [0.0, 0.5, 1.0] + end + + @testset "save_start / save_end / save_on" begin + sol = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = [0.5], save_start = false) + ) + @test sol.t == [0.5, 1.0] + + sol = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = [0.5], save_end = false) + ) + @test sol.t == [0.0, 0.5] + + sol = DiffEqBase.solve!( + DiffEqBase.init( + prob, ltg(); dt = 0.1, saveat = [0.5], + save_on = false, save_everystep = true + ) + ) + @test isempty(sol.t) + @test isempty(sol.u) + end + + @testset "backward integration keeps the direction ordering" begin + bprob = make_prob((1.0, 0.0)) + sol = DiffEqBase.solve!( + DiffEqBase.init(bprob, ltg(); dt = 0.1, saveat = [0.25, 0.75]) + ) + @test sol.t == [1.0, 0.75, 0.25, 0.0] + + sol = DiffEqBase.solve!(DiffEqBase.init(bprob, ltg(); dt = 0.1, saveat = 0.25)) + @test sol.t == [1.0, 0.75, 0.5, 0.25, 0.0] + end + + @testset "asking for output does not change the trajectory" begin + # saveat points are interpolated, never stepped onto, so the sequence of + # steps -- and hence the splitting error -- is untouched. + plain = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1)) + dense = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = 0.017, save_everystep = true) + ) + @test dense.u[end] == plain.u[end] + end + + @testset "saveat times passed as tstops are landed on, not interpolated" begin + # The documented way to get output that is accurate to the order of the + # scheme: the integrator steps exactly onto the time, so the saved value is + # the stepped state rather than the (first order) interpolant. + ts = [0.25, 0.55] + landed = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, tstops = ts, saveat = ts) + ) + interpolated = DiffEqBase.solve!( + DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = ts) + ) + @test issorted(landed.t) + @test allunique(landed.t) + for t in ts + @test t in landed.t + @test t in interpolated.t + end + # Same requested output times, two different values: one is the stepped + # state, the other the interpolant of the step that straddles the time. + for t in ts + a = landed.u[findfirst(==(t), landed.t)] + b = interpolated.u[findfirst(==(t), interpolated.t)] + @test a != b + # ... but both approximate the same solution. + @test isapprox(a, b; rtol = 1.0e-3) + end + end + + @testset "add_saveat!" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + DiffEqBase.step!(integrator) + DiffEqBase.add_saveat!(integrator, 0.65) + sol = DiffEqBase.solve!(integrator) + @test sol.t == [0.0, 0.65, 1.0] + # Cannot save behind the current time. + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + DiffEqBase.step!(integrator) + DiffEqBase.step!(integrator) + @test_throws ErrorException DiffEqBase.add_saveat!(integrator, 0.05) + end + + @testset "other algorithms" begin + for alg in ( + StrangMarchuk((Euler(), Euler())), + StrangMarchuk((Tsit5(), Tsit5())), + PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())), + ) + sol = DiffEqBase.solve!( + DiffEqBase.init(prob, alg; dt = 0.1, saveat = [0.5]) + ) + @test sol.t[1] == 0.0 + @test sol.t[end] == 1.0 + @test 0.5 in sol.t + @test sol.retcode == ReturnCode.Success + end + end + + @testset "nested splitting saves only at the outer level" begin + inner = GenericSplitFunction((f1, f1), ([1, 2], [1, 2])) + outer = GenericSplitFunction((f1, inner), ([1, 2, 3], [1, 2])) + nprob = OperatorSplittingProblem(outer, copy(U0), TSPAN) + alg = LieTrotterGodunov((Euler(), LieTrotterGodunov((Euler(), Euler())))) + integrator = DiffEqBase.init(nprob, alg; dt = 0.1, saveat = [0.5]) + sol = DiffEqBase.solve!(integrator) + @test sol.t == [0.0, 0.5, 1.0] + # The inner node owns no solution storage of its own. + @test !hasproperty(integrator.child_subintegrators[2], :saveiter) + end +end + +@testset "interpolation" begin + prob = make_prob() + + @testset "reproduces the step endpoints" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + DiffEqBase.step!(integrator) + @test integrator(integrator.tprev) ≈ integrator.uprev + @test integrator(integrator.t) ≈ integrator.u + end + + @testset "is linear in between" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + DiffEqBase.step!(integrator) + (; tprev, t, uprev, u) = integrator + tmid = tprev + (t - tprev) / 3 + Θ = (tmid - tprev) / (t - tprev) + @test integrator(tmid) ≈ @. (1 - Θ) * uprev + Θ * u + end + + @testset "first derivative is the step slope" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + DiffEqBase.step!(integrator) + (; tprev, t, uprev, u) = integrator + expected = (u .- uprev) ./ (t - tprev) + @test integrator(tprev + (t - tprev) / 2, Val{1}) ≈ expected + # Constant slope: the derivative does not depend on where it is evaluated. + @test integrator(t, Val{1}) ≈ expected + end + + @testset "in-place and idxs" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + DiffEqBase.step!(integrator) + tmid = (integrator.tprev + integrator.t) / 2 + out = similar(integrator.u) + integrator(out, tmid) + @test out ≈ integrator(tmid) + + @test integrator(tmid, Val{0}; idxs = 2) ≈ integrator(tmid)[2] + @test integrator(tmid, Val{0}; idxs = [1, 3]) ≈ integrator(tmid)[[1, 3]] + + sub = similar(integrator.u, 2) + integrator(sub, tmid, Val{0}; idxs = [1, 3]) + @test sub ≈ integrator(tmid)[[1, 3]] + end + + @testset "does not read the stale dt" begin + # After an accepted step the controller has already replaced `dt` with the + # next proposal, so an interpolant keyed off `integrator.dt` would be wrong. + integrator = DiffEqBase.init( + prob, PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())); dt = 0.1 + ) + DiffEqBase.step!(integrator) + DiffEqBase.step!(integrator) + @test integrator.dt != integrator.t - integrator.tprev + @test integrator(integrator.t) ≈ integrator.u + @test integrator(integrator.tprev) ≈ integrator.uprev + end + + @testset "sol(t) matches the integrator interpolant" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, save_everystep = true) + sol = DiffEqBase.solve!(integrator) + # Between two saved step endpoints `sol.interp` is the same linear + # interpolant the integrator uses within a step. + i = 4 + t0, t1 = sol.t[i], sol.t[i + 1] + tmid = (t0 + t1) / 2 + @test sol(tmid) ≈ (sol.u[i] .+ sol.u[i + 1]) ./ 2 + end + + @testset "saveat values are the step interpolant" begin + tmid = 0.05 # strictly inside the first step of length 0.1 + integrator = DiffEqBase.init( + prob, ltg(); dt = 0.1, saveat = [tmid], save_everystep = true + ) + sol = DiffEqBase.solve!(integrator) + idx = findfirst(==(tmid), sol.t) + @test idx !== nothing + # Bracketing saved points are the enclosing step's endpoints. + @test sol.u[idx] ≈ (sol.u[idx - 1] .+ sol.u[idx + 1]) ./ 2 + end +end + +@testset "change_t_via_interpolation!" begin + prob = make_prob() + + @testset "moves the whole tree back to the interpolated state" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + DiffEqBase.step!(integrator) + (; tprev, t) = integrator + tmid = (tprev + t) / 2 + expected = integrator(tmid) + + SciMLBase.change_t_via_interpolation!(integrator, tmid) + @test integrator.t == tmid + @test integrator.u ≈ expected + # Every child clock and buffer follows the parent. + OS.validate_time_point(integrator) + for child in integrator.child_subintegrators + @test child.t == tmid + end + @test integrator.child_subintegrators[1].u ≈ expected[1:3] + @test integrator.child_subintegrators[2].u ≈ expected[1:2] + end + + @testset "refuses to move before tprev" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + DiffEqBase.step!(integrator) + DiffEqBase.step!(integrator) + @test_throws ErrorException SciMLBase.change_t_via_interpolation!( + integrator, integrator.tprev - 0.01 + ) + end + + @testset "moving to the current t is a no-op" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + DiffEqBase.step!(integrator) + u_before = copy(integrator.u) + SciMLBase.change_t_via_interpolation!(integrator, integrator.t) + @test integrator.u == u_before + end + + @testset "accepts the Val{:false} the callback path passes" begin + # DiffEqBase's continuous callback path calls this with `Val{:false}` -- a + # `Val` of the Symbol -- so the flag must be decoded by dispatch. + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + DiffEqBase.step!(integrator) + tmid = (integrator.tprev + integrator.t) / 2 + SciMLBase.change_t_via_interpolation!(integrator, tmid, Val{:false}, nothing) + @test integrator.t == tmid + end + + @testset "modify_save_endpoint rewrites the saved endpoint" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, save_everystep = true) + DiffEqBase.step!(integrator) + nsaved = integrator.saveiter + tmid = (integrator.tprev + integrator.t) / 2 + SciMLBase.change_t_via_interpolation!(integrator, tmid, Val{true}) + @test integrator.saveiter == nsaved + 1 + @test integrator.sol.t[integrator.saveiter] == tmid + end +end + +@testset "solution object plumbing" begin + prob = make_prob() + + @testset "dense output is rejected" begin + @test_throws ArgumentError DiffEqBase.init(prob, ltg(); dt = 0.1, dense = true) + end + + @testset "sol.stats is a real DEStats" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + @test integrator.sol.stats !== nothing + @test integrator.sol.stats.ncondition == 0 + end + + @testset "postamble! is idempotent" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + sol = DiffEqBase.solve!(integrator) + n = length(sol.t) + SciMLBase.postamble!(integrator) + SciMLBase.postamble!(integrator) + @test length(integrator.sol.t) == n + end + + @testset "postamble! does not run once per step" begin + finalized = Ref(0) + cb = DiscreteCallback( + (u, t, integrator) -> false, + integrator -> nothing; + finalize = (c, u, t, integrator) -> (finalized[] += 1) + ) + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb) + DiffEqBase.solve!(integrator) + @test finalized[] == 1 + end +end + +@testset "reinit!" begin + prob = make_prob() + + @testset "restores the saved solution" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = [0.5]) + first_run = copy(DiffEqBase.solve!(integrator).t) + DiffEqBase.reinit!(integrator) + @test integrator.saveiter == 1 + @test integrator.sol.t[1] == 0.0 + DiffEqBase.solve!(integrator) + @test integrator.sol.t == first_run + end + + @testset "reinit_callbacks = false works without a saving callback" begin + # Regression: this used to index an empty `discrete_callbacks` tuple. + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + DiffEqBase.solve!(integrator) + DiffEqBase.reinit!(integrator; reinit_callbacks = false) + @test integrator.t == 0.0 + @test integrator.saveiter == 1 + @test DiffEqBase.solve!(integrator).retcode == ReturnCode.Success + end + + @testset "a new u0 reaches the whole tree" begin + # Regression: a leaf's `reinit!` restores the `u0` slice captured when it was + # built, and the palindromic schemes skip the forward sync of their first + # child -- so `reinit!(integrator, u0)` used to integrate a stale child state + # and silently return the wrong answer. + u0b = [10.0, 20.0, 30.0] + probb = OperatorSplittingProblem(fsplit, copy(u0b), TSPAN) + for alg in ( + ltg(), + StrangMarchuk((Euler(), Euler())), + StrangMarchuk((Tsit5(), Euler())), + PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())), + ) + integrator = DiffEqBase.init(make_prob(), alg; dt = 0.1) + DiffEqBase.solve!(integrator) + DiffEqBase.reinit!(integrator, copy(u0b)) + @test integrator.child_subintegrators[1].u == u0b + @test integrator.child_subintegrators[2].u == u0b[1:2] + got = DiffEqBase.solve!(integrator).u[end] + ref = DiffEqBase.solve!(DiffEqBase.init(probb, alg; dt = 0.1)).u[end] + @test got ≈ ref rtol = 1.0e-12 + end + end + + @testset "erase_sol clears the stored solution" begin + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, save_everystep = true) + DiffEqBase.solve!(integrator) + DiffEqBase.reinit!(integrator; erase_sol = true) + @test integrator.saveiter == 1 + @test length(integrator.sol.t) == 1 + end +end + +@testset "callbacks can be attached" begin + # Regression: DiffEqBase's `initialize!` reads the `derivative_discontinuity` + # field directly, so `init` used to throw as soon as any callback was passed. + prob = make_prob() + calls = Ref(0) + cb = DiscreteCallback((u, t, integrator) -> false, integrator -> (calls[] += 1)) + integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb) + @test integrator.derivative_discontinuity == false + sol = DiffEqBase.solve!(integrator) + @test sol.retcode == ReturnCode.Success + @test calls[] == 0 + + integrator = DiffEqBase.init( + prob, ltg(); dt = 0.1, + callback = CallbackSet(cb, cb) + ) + @test DiffEqBase.solve!(integrator).retcode == ReturnCode.Success +end From e379348e336a289c8f2f23cd8b87754215999e81 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 3 Aug 2026 14:36:22 -0400 Subject: [PATCH 2/4] Simplify saving/callback code and its OrdinaryDiffEq interfaces Cleanup pass over the saving, interpolation and callback support, with no intended behavior change. Drop duplicated upstream code: - `apply_ith_callback!` was a byte-for-byte copy of `OrdinaryDiffEqCore.apply_ith_callback!`, which is generic in `integrator`. Call upstream's and add it to the QA ignore list. - `save_initial_value!` re-inlined `_save_current!` with a hardcoded index. - The two save branches of `savevalues!` now share one `_save_at!` writer. Drop compat shims made dead by the SciMLBase 3.1 bump: the `@static` pairs for `u_modified!` / `derivative_discontinuity!` collapse to single unconditional methods. The `DEVerbosity` shims stay, since they are gated on DiffEqBase. Remove needless indirection: - `resync_children_after_modification!` was a pure alias for `rollback_children!`; its three call sites now go direct. - `_modify_save_endpoint`'s dispatch pair becomes one identity comparison, which tolerates DiffEqBase's `Val{:false}` just as safely. - `_step_dt` and `_interp_theta` merge into `_interp_coords`, so the step length is computed once per interpolation instead of twice. - Rename the local `solution_endpoint_match_cur_integrator!` to `save_endpoint!`: it shadowed an unrelated OrdinaryDiffEqCore function of the same name. Avoid work in the step loop: `savevalues!` returns early when there is nothing to save, and building the `saveat` heap no longer allocates twice. Finally, the "inner splits are stages, not steps" argument was written out eight times across src/ and the docs. It now lives once in the devdocs "Dense output" section, with one-line references from the code. Co-Authored-By: Claude Opus 5 --- docs/src/devdocs/index.md | 5 + docs/src/usage/index.md | 9 +- src/integrator.jl | 230 ++++++++++++-------------------------- src/utils.jl | 28 ++--- test/qa/qa.jl | 1 + 5 files changed, 92 insertions(+), 181 deletions(-) diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md index bd63d37..7e9cd3f 100644 --- a/docs/src/devdocs/index.md +++ b/docs/src/devdocs/index.md @@ -144,6 +144,11 @@ sub-problems evaluated over different intervals, and they do not compose into an approximation of the split solution. Only the step endpoints, which the outer integrator owns, are states of the full split system. +The same structural argument is why saving and callbacks live on the outer integrator +alone: an inner split is a stage, not a step, so there is no time point at which its +state is a meaningful approximation of the split solution for a saved point to record +or for a callback condition to act on. + Two invariants matter when implementing this: - `Θ` is derived from `integrator.t - integrator.tprev`, **not** from diff --git a/docs/src/usage/index.md b/docs/src/usage/index.md index 00a71c3..598821b 100644 --- a/docs/src/usage/index.md +++ b/docs/src/usage/index.md @@ -123,11 +123,10 @@ DiffEqCallbacks.jl) work too. !!! note "Callbacks run on the outer integrator only" A condition is evaluated once per **outer** step, after all the operators of - that step have been applied. It is never evaluated between two inner splits: a - splitting step applies its operators sequentially over staggered subintervals, - so the intermediate states are stages and do not approximate the solution of the - split system at any time point. There is correspondingly nothing meaningful for - a condition to test or an `affect!` to modify at that level. + that step have been applied, and never between two inner splits: for the reason + given under "Saving and interpolation" above, those intermediate states are + stages and approximate the split solution at no time point, so there is nothing + meaningful for a condition to test or an `affect!` to modify at that level. A consequence worth knowing: an `affect!` that modifies `integrator.u` is propagated into every subintegrator before the next step, so modifying the state diff --git a/src/integrator.jl b/src/integrator.jl index 27e42de..e81ed20 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -27,14 +27,10 @@ end The saving settings of the outermost [`OperatorSplittingIntegrator`](@ref). -Saving is a property of the outer integrator only: the inner splits are stages -rather than steps, so there is no time point at which their state is a meaningful -approximation of the split solution. These settings are therefore consumed by -`__init` directly and never enter the [`ConfigTree`](@ref), which keeps them from -travelling down to the leaf integrators. - -`saveat` itself is not held here -- it lives in the integrator's `saveat` heap, -because it is consumed as the integration proceeds. +Saving is a property of the outer integrator only, so these are consumed by `__init` +directly and never enter the [`ConfigTree`](@ref), which is what keeps them from +travelling down to the leaf integrators. `saveat` is not held here: it lives in the +integrator's `saveat` heap, because it is consumed as the integration proceeds. """ struct SaveOptions save_on::Bool @@ -125,7 +121,6 @@ mutable struct SplitSubIntegrator{ controller_cache::controllerType force_stepfail::Bool last_step_failed::Bool - # Named after the field DiffEqBase's callback machinery reads directly. derivative_discontinuity::Bool status::SplitSubIntegratorStatus stats::IntegratorStats @@ -200,9 +195,8 @@ mutable struct OperatorSplittingIntegrator{ # Scratch buffers a `VectorContinuousCallback` needs; `nothing` when the callback # set contains none. const callback_cache::callbackCacheType - # Continuous-callback bookkeeping. DiffEqBase's root finder reads - # `event_last_time`/`last_event_error` to avoid re-detecting the event it just - # handled, and `handle_callbacks!` is responsible for writing them. + # DiffEqBase's root finder reads these to avoid re-detecting the event it just + # handled; `handle_callbacks!` is responsible for writing them. event_last_time::Int vector_event_last_time::Int last_event_error::eventErrType @@ -287,10 +281,6 @@ function SciMLBase.__init( (; u0, p) = prob t0, tf = prob.tspan - # Dense output would need per-step interpolation data on top of the saved - # points. The interpolant of a splitting step is linear (see the interpolation - # section below), so `sol(t)` already reproduces it exactly from `sol.t`/`sol.u` - # via SciMLBase's `LinearInterpolation` and storing `k` would buy nothing. dense && throw( ArgumentError( "dense output is not supported by operator splitting integrators. The \ @@ -436,39 +426,25 @@ function initialize_callbacks!(integrator::OperatorSplittingIntegrator) ) if modified update_uprev!(integrator) - resync_children_after_modification!(integrator) + rollback_children!(integrator) end integrator.derivative_discontinuity = false return nothing end -# A callback that changed `u` invalidates every child's copy of its slice. -# `rollback_children!` refills the children from the parent's `u`, re-anchors their -# clocks and marks their state modified so leaf integrators drop FSAL data. It also -# restores the `parent.u == child.u` invariant that lets the palindromic schemes skip -# the forward sync of their first child (`next_sync_is_continuous`). -function resync_children_after_modification!(integrator::OperatorSplittingIntegrator) - rollback_children!(integrator) - return nothing -end - function SciMLBase.reeval_internals_due_to_modification!( integrator::OperatorSplittingIntegrator, continuous_modification = true; callback_initializealg = nothing ) - resync_children_after_modification!(integrator) + rollback_children!(integrator) integrator.derivative_discontinuity = false return nothing end -# `save_start` puts u(t0) at index 1, so the first `savevalues!` writes index 2. function save_initial_value!(integrator::OperatorSplittingIntegrator) + integrator.saveiter = 0 if integrator.save_opts.save_on && integrator.save_opts.save_start - integrator.saveiter = 1 - RecursiveArrayTools.copyat_or_push!(integrator.sol.t, 1, integrator.t) - RecursiveArrayTools.copyat_or_push!(integrator.sol.u, 1, integrator.u) - else - integrator.saveiter = 0 + _save_current!(integrator) end return nothing end @@ -540,7 +516,7 @@ function DiffEqBase.reinit!( # to the `u0` handed to this call, so the new state has to be pushed down # explicitly. Relying on the forward sync of the next step is not enough: the # palindromic schemes skip it for their first child. - resync_children_after_modification!(integrator) + rollback_children!(integrator) # After the children, so that an initializer which modifies `u` (and therefore # resyncs the tree) is not undone by the child reinit, and so that the first saved @@ -730,10 +706,14 @@ end # 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. +# are moved back to the parent's time. This is also the resync path for a callback +# that modified `u`. +# +# Leaves get their `u` restored here as well, not only by the forward sync before the +# next solve: palindromic algorithms skip that sync for their first child +# (`next_sync_is_continuous`), so a rollback that left the leaf state stale would +# silently resume from the failed attempt. Restoring it here is also what re- +# establishes the `parent.u == child.u` invariant that shortcut relies on. rollback_children!(parent::AnySplitIntegrator) = _rollback_children!( parent.child_subintegrators, parent.child_solution_indices, parent.u, parent.t ) @@ -1085,11 +1065,8 @@ end # failure finalizes here. function SciMLBase.check_error!(integrator::OperatorSplittingIntegrator) code = SciMLBase.check_error(integrator) - # Rebuilding the solution allocates, and this runs before every step, so only - # touch it when the code actually changed. - if integrator.sol.retcode !== code - integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, code) - end + # Rebuilding the solution allocates, and this runs before every step. + integrator.sol.retcode !== code && _set_retcode!(integrator, code) if code ∉ (ReturnCode.Success, ReturnCode.Default) SciMLBase.postamble!(integrator) end @@ -1120,14 +1097,9 @@ end # --------------------------------------------------------------------------- # Interpolation # -# A splitting step advances its children sequentially over staggered subintervals, -# so the children's own interpolants do not compose into an approximation of the -# split solution: the only dense output well defined at this level is the one built -# from the endpoints the outer integrator owns. The generic fallback is therefore -# linear, which is exact for the state a `LieTrotterGodunov` step produces and first -# order for the second order schemes. An algorithm that can do better implements -# `splitting_interpolant`/`splitting_interpolant!` for its own cache type; nothing -# else in the saving or callback machinery has to change. +# See the "Dense output" section of the developer documentation for why the only +# dense output well defined here is the one built from the endpoints the outer +# integrator owns, and hence why the generic fallback is linear. # --------------------------------------------------------------------------- """ @@ -1157,24 +1129,23 @@ function splitting_interpolant!( return out end -# The step that was just taken. `integrator.dt` must not be used here: by the time a -# step is finished `step_accept_controller!` has already overwritten it with the -# *next* proposal, so it is not the length of the interval `uprev`/`u` span. -_step_dt(integrator::OperatorSplittingIntegrator) = integrator.t - integrator.tprev - -function _interp_theta(integrator::OperatorSplittingIntegrator, t) - dt = _step_dt(integrator) - # A zero-length interval only happens before the first step, where `uprev == u` - # and every Θ gives the same value. - return iszero(dt) ? zero(t / oneunit(dt)) : (t - integrator.tprev) / dt +# Step-local coordinate `Θ` and the length of the step that was just taken. +# `integrator.dt` must not be used for the latter: by the time a step is finished +# `step_accept_controller!` has already overwritten it with the *next* proposal, so it +# is not the length of the interval `uprev`/`u` span. A zero-length interval only +# happens before the first step, where `uprev == u` and every Θ gives the same value. +function _interp_coords(integrator::OperatorSplittingIntegrator, t) + dt = integrator.t - integrator.tprev + Θ = iszero(dt) ? zero(t / oneunit(dt)) : (t - integrator.tprev) / dt + return Θ, dt end function (integrator::OperatorSplittingIntegrator)( t::Number, ::Type{deriv} = Val{0}; idxs = nothing ) where {deriv} + Θ, dt = _interp_coords(integrator, t) return splitting_interpolant( - integrator, integrator.cache, - _interp_theta(integrator, t), _step_dt(integrator), + integrator, integrator.cache, Θ, dt, integrator.uprev, integrator.u, idxs, deriv ) end @@ -1182,9 +1153,9 @@ end function (integrator::OperatorSplittingIntegrator)( val::AbstractArray, t::Number, ::Type{deriv} = Val{0}; idxs = nothing ) where {deriv} + Θ, dt = _interp_coords(integrator, t) return splitting_interpolant!( - val, integrator, integrator.cache, - _interp_theta(integrator, t), _step_dt(integrator), + val, integrator, integrator.cache, Θ, dt, integrator.uprev, integrator.u, idxs, deriv ) end @@ -1193,12 +1164,6 @@ end # change_t_via_interpolation! # --------------------------------------------------------------------------- -# DiffEqBase's continuous callback path passes the literal `Val{:false}` -- a `Val` of -# the *Symbol*, not of `false` -- so the flag has to be decoded by dispatch. Testing -# `if T` the way OrdinaryDiffEqCore does would throw on a Symbol. -_modify_save_endpoint(::Type{Val{true}}) = true -_modify_save_endpoint(::Type{Val{T}}) where {T} = false - """ change_t_via_interpolation!(integrator::OperatorSplittingIntegrator, t, modify_save_endpoint = Val{false}, reinitialize_alg = nothing) @@ -1207,10 +1172,8 @@ Move the integrator back to a time `t` inside the step that was just taken. `u` is refilled from the step's interpolant and the whole subintegrator tree is re-anchored to `t` through [`rollback_children!`](@ref), so every child's state and clock stay consistent with the parent (`validate_time_point` asserts the latter). - -`integrator.dt` is deliberately left alone: it holds the step size proposed for the -next step, and the interpolation coordinate is derived from `t - tprev` rather than -from `dt`. +`integrator.dt` is deliberately left alone: it already holds the step size proposed +for the next step, and the interpolation coordinate comes from `t - tprev`. """ function SciMLBase.change_t_via_interpolation!( integrator::OperatorSplittingIntegrator, t, @@ -1222,10 +1185,11 @@ function SciMLBase.change_t_via_interpolation!( elseif t != integrator.t integrator(integrator.u, t) integrator.t = t - # Push the interpolated slices into the children and move their clocks back. rollback_children!(integrator) - if _modify_save_endpoint(modify_save_endpoint) - solution_endpoint_match_cur_integrator!(integrator) + # Compared by identity, not `if T`: DiffEqBase's continuous callback path + # passes `Val{:false}`, a `Val` of the *Symbol*, which `if T` would throw on. + if modify_save_endpoint === Val{true} + save_endpoint!(integrator) end end return nothing @@ -1234,30 +1198,26 @@ end # --------------------------------------------------------------------------- # Saving # -# Only the outermost integrator saves: the inner splits are stages rather than -# steps, so their intermediate states are not approximations of the split solution -# at any time point. -# -# `saveiter` is the authoritative length of the saved prefix of `sol.t`/`sol.u`. -# `sol.interp` is a `LinearInterpolation` aliasing those same vectors, so `sol(t)` -# reproduces this integrator's own interpolant between saved points. +# Outer integrator only, for the reason given in the "Dense output" section of the +# developer documentation. `saveiter` is the authoritative length of the saved prefix +# of `sol.t`/`sol.u`. # --------------------------------------------------------------------------- _saved_at_current_t(integrator::OperatorSplittingIntegrator) = integrator.saveiter > 0 && integrator.sol.t[integrator.saveiter] == integrator.t -function _save_current!(integrator::OperatorSplittingIntegrator) +# `copy_u = false` hands over ownership of a `u` the caller freshly allocated. +function _save_at!(integrator::OperatorSplittingIntegrator, t, u, copy_u = true) integrator.saveiter += 1 - RecursiveArrayTools.copyat_or_push!( - integrator.sol.t, integrator.saveiter, integrator.t - ) - RecursiveArrayTools.copyat_or_push!( - integrator.sol.u, integrator.saveiter, integrator.u - ) + RecursiveArrayTools.copyat_or_push!(integrator.sol.t, integrator.saveiter, t) + RecursiveArrayTools.copyat_or_push!(integrator.sol.u, integrator.saveiter, u, copy_u) return nothing end +_save_current!(integrator::OperatorSplittingIntegrator) = + _save_at!(integrator, integrator.t, integrator.u) + """ savevalues!(integrator::OperatorSplittingIntegrator, force_save = false, reduce_size = true) @@ -1279,6 +1239,8 @@ function SciMLBase.savevalues!( savedexactly = false integrator.save_opts.save_on || return saved, savedexactly (; save_everystep, save_end) = integrator.save_opts + (force_save || save_everystep || !isempty(integrator.saveat)) || + return saved, savedexactly tf = integrator.sol.prob.tspan[2] # The heaps store raw times and carry the direction in their ordering, so the @@ -1290,20 +1252,13 @@ function SciMLBase.savevalues!( if curt == integrator.t # `save_end` owns the final point; leave it to the postamble. (!save_end && curt == tf) && continue - saved = true savedexactly = true _save_current!(integrator) else - saved = true - integrator.saveiter += 1 - RecursiveArrayTools.copyat_or_push!( - integrator.sol.t, integrator.saveiter, curt - ) - # Freshly allocated by the interpolant, so hand over ownership. - RecursiveArrayTools.copyat_or_push!( - integrator.sol.u, integrator.saveiter, integrator(curt), false - ) + # The interpolant allocates a fresh `u`, so hand over ownership. + _save_at!(integrator, curt, integrator(curt), false) end + saved = true end if force_save || ( @@ -1318,9 +1273,8 @@ function SciMLBase.savevalues!( return saved, savedexactly end -# Make sure the final time point is in the solution. Mirrors OrdinaryDiffEqCore's -# `solution_endpoint_match_cur_integrator!`. -function solution_endpoint_match_cur_integrator!(integrator::OperatorSplittingIntegrator) +# Make sure the final time point is in the solution. +function save_endpoint!(integrator::OperatorSplittingIntegrator) integrator.save_opts.save_on || return nothing integrator.save_opts.save_end || return nothing _saved_at_current_t(integrator) && return nothing @@ -1331,43 +1285,13 @@ end # --------------------------------------------------------------------------- # Callbacks # -# Callbacks run on the outermost integrator only. A condition evaluated between two -# inner splits would see a state that is not an approximation of the split solution -# at any time point -- the splits are stages, not steps -- so there is nothing -# meaningful for a condition or an `affect!` to act on down there. -# -# The heavy lifting (root finding, `save_positions`, `affect!` dispatch) is -# DiffEqBase's; this is the per-step driver, modelled on OrdinaryDiffEqCore's -# `handle_callbacks!`. It also owns the plain `saveat` save when no callback saved. +# Outer integrator only, for the reason given in the "Dense output" section of the +# developer documentation. Root finding, `save_positions` handling and `affect!` +# dispatch are DiffEqBase's; this is just the per-step driver, modelled on +# OrdinaryDiffEqCore's `handle_callbacks!`. It also owns the plain `saveat` save +# when no callback saved. # --------------------------------------------------------------------------- -# Type-stable dispatch into `DiffEqBase.apply_callback!`: the explicit -# `if (return) else (rest)` chain compiles to a switch over a heterogeneous tuple. -@generated function apply_ith_callback!( - integrator, - time, upcrossing, event_idx, cb_idx, - callbacks::NTuple{ - N, - Union{DiffEqBase.ContinuousCallback, DiffEqBase.VectorContinuousCallback}, - } - ) where {N} - ex = quote - throw(BoundsError(callbacks, cb_idx)) - end - for i in 1:N - ex = quote - if (cb_idx == $i) - return DiffEqBase.apply_callback!( - integrator, callbacks[$i], time, upcrossing, event_idx - ) - else - $ex - end - end - end - return ex -end - handle_callbacks!(::SplitSubIntegrator) = nothing function handle_callbacks!(integrator::OperatorSplittingIntegrator) @@ -1386,7 +1310,7 @@ function handle_callbacks!(integrator::OperatorSplittingIntegrator) if event_occurred integrator.event_last_time = idx integrator.vector_event_last_time = event_idx - continuous_modified, saved_in_cb = apply_ith_callback!( + continuous_modified, saved_in_cb = OrdinaryDiffEqCore.apply_ith_callback!( integrator, time, upcrossing, event_idx, idx, continuous_callbacks ) # The step was cut at the event and the state may have jumped, so the @@ -1577,7 +1501,7 @@ function SciMLBase.postamble!(integrator::OperatorSplittingIntegrator) integrator.postamble_done && return nothing integrator.postamble_done = true DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator) - solution_endpoint_match_cur_integrator!(integrator) + save_endpoint!(integrator) # `saveiter` is authoritative: drop whatever a previous, longer run left behind # (a `reinit!` without `erase_sol`). resize!(integrator.sol.t, integrator.saveiter) @@ -1880,7 +1804,6 @@ SciMLBase.pop_tstop!(i::AnySplitIntegrator) = pop!(i.tstops) DiffEqBase.get_dt(i::AnySplitIntegrator) = i.dt # Continuous callbacks relax the step size through this after an event. -# `dtcache` mirrors the standing proposal as a magnitude (see `__init`). function SciMLBase.set_proposed_dt!(integrator::OperatorSplittingIntegrator, dt) if integrator.dtcache != abs(dt) integrator.dtcache = abs(dt) @@ -1894,11 +1817,8 @@ end """ terminate!(integrator::OperatorSplittingIntegrator, retcode = ReturnCode.Terminated) -Stop the integration from inside a callback's `affect!`. - -Emptying the outermost `tstops` heap is what ends the `solve!` loop and makes -`done` report completion; the children's heaps are irrelevant once the outer loop -exits, and `postamble!` still closes out the solution. +Stop the integration from inside a callback's `affect!`. Emptying the outermost +`tstops` heap is what ends the `solve!` loop and makes `done` report completion. """ function SciMLBase.terminate!( integrator::OperatorSplittingIntegrator, retcode = ReturnCode.Terminated @@ -1930,13 +1850,5 @@ function DiffEqBase.add_saveat!(i::OperatorSplittingIntegrator, t) return nothing end -# SciMLBase v3 renamed `u_modified!` → `derivative_discontinuity!`. The flag itself -# is stored in the `derivative_discontinuity` field, which is the name DiffEqBase's -# callback machinery reads directly. -@static if isdefined(DiffEqBase, :u_modified!) - DiffEqBase.u_modified!(i::AnySplitIntegrator, bool) = i.derivative_discontinuity = bool -end -@static if isdefined(SciMLBase, :derivative_discontinuity!) - SciMLBase.derivative_discontinuity!(i::AnySplitIntegrator, bool) = - i.derivative_discontinuity = bool -end +SciMLBase.derivative_discontinuity!(i::AnySplitIntegrator, bool) = + i.derivative_discontinuity = bool diff --git a/src/utils.jl b/src/utils.jl index a635f2a..8f4b469 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -7,22 +7,21 @@ function tstops_and_saveat_heaps(t0, tf, tstops, saveat) tstops = [filter(t -> t0 < t < tf || tf < t < t0, tstops)..., tf] tstops = BinaryHeaps.BinaryHeap{FT, ordering}(tstops) - # `t0` is excluded and `tf` included, matching OrdinaryDiffEqCore's - # `initialize_saveat`: the initial point is owned by `save_start` and the final - # one by `save_end`, so leaving them in the heap would duplicate them. + # Keep `t0 < t <= tf` in tdir-space, as OrdinaryDiffEqCore's `initialize_saveat` + # does: the initial point is owned by `save_start` and the final one by + # `save_end`, so leaving either in the heap would duplicate it. tdir = tf > t0 ? one(FT) : -one(FT) - if isnothing(saveat) - saveat = FT[] + saveat = if isnothing(saveat) + FT[] elseif saveat isa Number saveat > zero(saveat) || error("saveat value must be positive") + # The range already satisfies the window, so no filtering is needed: it + # starts one step past `t0` and hits `tf` only if the span divides evenly. step = tdir * saveat - # `tf` is not appended: the range hits it when it divides the span evenly, - # and `save_end` owns the final point otherwise. - saveat = collect((t0 + step):step:tf) + collect((t0 + step):step:tf) else - saveat = collect(FT, saveat) + filter(t -> tdir * t0 < tdir * t <= tdir * tf, collect(FT, saveat)) end - saveat = filter(t -> tdir * t0 < tdir * t <= tdir * tf, saveat) saveat = BinaryHeaps.BinaryHeap{FT, ordering}(saveat) return tstops, saveat @@ -80,14 +79,9 @@ function forward_sync_subintegrator!( end # 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. +# FSAL information. function mark_state_modified!(child::DEIntegrator) - @static if isdefined(SciMLBase, :derivative_discontinuity!) - SciMLBase.derivative_discontinuity!(child, true) - else - SciMLBase.u_modified!(child, true) - end + SciMLBase.derivative_discontinuity!(child, true) return nothing end diff --git a/test/qa/qa.jl b/test/qa/qa.jl index 8de8428..42f9f80 100644 --- a/test/qa/qa.jl +++ b/test/qa/qa.jl @@ -31,6 +31,7 @@ run_qa( :setup_controller_cache, :reinit_controller!, :post_newton_controller!, # OrdinaryDiffEqCore :get_EEst, :set_EEst!, :get_current_adaptive_order, # OrdinaryDiffEqCore :gamma_default, :failfactor_default, :AbstractControllerCache, # OrdinaryDiffEqCore + :apply_ith_callback!, # OrdinaryDiffEqCore :promote_tspan, # SciMLBase # Shared linear interpolation kernels, reused so that the # integrator's own dense output and `sol(t)` (which goes through From b87463303c9b75ded2915df5413efec90eeace42 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 3 Aug 2026 17:18:41 -0400 Subject: [PATCH 3/4] simplify --- docs/src/devdocs/index.md | 38 +++---- docs/src/usage/index.md | 15 ++- src/OrdinaryDiffEqOperatorSplitting.jl | 3 - src/integrator.jl | 129 +++++++++-------------- src/utils.jl | 7 +- test/callbacks.jl | 85 ++++++---------- test/saving.jl | 135 +++++++++---------------- 7 files changed, 153 insertions(+), 259 deletions(-) diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md index 7e9cd3f..30387a5 100644 --- a/docs/src/devdocs/index.md +++ b/docs/src/devdocs/index.md @@ -131,29 +131,19 @@ end ``` Both fall back to linear interpolation for any `AbstractOperatorSplittingCache`, so -implementing them is optional. Note the accuracy consequences: with the linear -fallback, interpolated output is first order even for a second-order scheme, and -continuous callback event times are the exact roots of a straight chord across the -step. Implementing a higher-order interpolant improves saved output *and* event -location at once. - -The reason the fallback is only linear is structural rather than incidental. A -splitting step advances its children sequentially over staggered subintervals, so at -any time strictly inside the step the children's own interpolants describe different -sub-problems evaluated over different intervals, and they do not compose into an +implementing them is optional -- but with the fallback, saved output is first order +even for a second-order scheme, and continuous callback event times are the exact +roots of a straight chord across the step. One method improves both. + +The fallback is linear for a structural reason, and it is the same reason saving and +callbacks live on the outer integrator alone: a splitting step advances its children +sequentially over staggered subintervals, so a child's own interpolant describes a +different sub-problem over a different interval, and they do not compose into an approximation of the split solution. Only the step endpoints, which the outer -integrator owns, are states of the full split system. +integrator owns, are states of the full split system -- an inner split is a stage, +not a step. -The same structural argument is why saving and callbacks live on the outer integrator -alone: an inner split is a stage, not a step, so there is no time point at which its -state is a meaningful approximation of the split solution for a saved point to record -or for a callback condition to act on. - -Two invariants matter when implementing this: - - - `Θ` is derived from `integrator.t - integrator.tprev`, **not** from - `integrator.dt`. Once a step is accepted, `step_accept_controller!` has already - replaced `dt` with the step size proposed for the *next* step. - - Interpolation must not mutate integrator state. `change_t_via_interpolation!` - relies on being able to evaluate the interpolant repeatedly (the callback - root-finder does so many times per step) before committing to a time. +When implementing this, note that `Θ` is derived from +`integrator.t - integrator.tprev`, **not** from `integrator.dt`: once a step is +accepted, `step_accept_controller!` has already replaced `dt` with the step size +proposed for the *next* step. diff --git a/docs/src/usage/index.md b/docs/src/usage/index.md index 598821b..277ed4d 100644 --- a/docs/src/usage/index.md +++ b/docs/src/usage/index.md @@ -101,9 +101,10 @@ never changes the splitting error. sol = solve!(init(prob, alg; dt = 0.1, tstops = [0.25, 0.5], saveat = [0.25, 0.5])) ``` - `dense = true` is rejected: with a linear interpolant, `sol(t)` already - reproduces the integrator's own dense output from the saved points, so storing - per-step interpolation data would buy nothing. + `dense = true` is accepted but has no effect: `sol(t)` already reproduces the + integrator's own interpolant from the saved points, so there is no per-step + interpolation data left to store. Use `save_everystep` or `saveat` to control + how finely that interpolant resolves the trajectory. ## Callbacks and events @@ -115,11 +116,9 @@ cb = ContinuousCallback((u, t, integrator) -> u[1] - 0.5, terminate!) sol = solve!(init(prob, alg; dt = 0.1, callback = cb)) ``` -`DiscreteCallback`, `ContinuousCallback`, `VectorContinuousCallback` and -`CallbackSet` are all supported, along with `save_positions`, `affect_neg!`, -`rootfind`, `terminate!` and `derivative_discontinuity!`. Because the machinery is -DiffEqBase's own, callbacks built on top of it (for instance those in -DiffEqCallbacks.jl) work too. +The machinery is DiffEqBase's own, so `DiscreteCallback`, `ContinuousCallback`, +`VectorContinuousCallback`, `CallbackSet` and anything built on top of them (for +instance DiffEqCallbacks.jl) all behave as they do elsewhere in SciML. !!! note "Callbacks run on the outer integrator only" A condition is evaluated once per **outer** step, after all the operators of diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index 765b9d8..19aaf77 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -42,9 +42,6 @@ abstract type AbstractOperatorSplittingAlgorithm end abstract type AbstractOperatorSplittingCache end @inline SciMLBase.isadaptive(::AbstractOperatorSplittingAlgorithm) = false -# `AbstractOperatorSplittingAlgorithm` is not a `SciMLBase.AbstractDEAlgorithm`, so it -# does not pick up the generic trait fallbacks. The continuous callback root finder -# asks for this one. @inline SciMLBase.isdiscrete(::AbstractOperatorSplittingAlgorithm) = false @inline isdtchangeable(alg::AbstractOperatorSplittingAlgorithm) = all(isdtchangeable.(alg.inner_algs)) diff --git a/src/integrator.jl b/src/integrator.jl index e81ed20..0d68db0 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -27,10 +27,9 @@ end The saving settings of the outermost [`OperatorSplittingIntegrator`](@ref). -Saving is a property of the outer integrator only, so these are consumed by `__init` -directly and never enter the [`ConfigTree`](@ref), which is what keeps them from -travelling down to the leaf integrators. `saveat` is not held here: it lives in the -integrator's `saveat` heap, because it is consumed as the integration proceeds. +Consumed by `__init` directly rather than through the [`ConfigTree`](@ref), which is +what keeps them from travelling down to the leaf integrators. `saveat` lives in the +integrator's `saveat` heap instead, because it is consumed as the integration proceeds. """ struct SaveOptions save_on::Bool @@ -192,11 +191,8 @@ mutable struct OperatorSplittingIntegrator{ saveat::heapType _saveat::saveatType callback::callbackType - # Scratch buffers a `VectorContinuousCallback` needs; `nothing` when the callback - # set contains none. + # Scratch buffers a `VectorContinuousCallback` needs; `nothing` when there is none. const callback_cache::callbackCacheType - # DiffEqBase's root finder reads these to avoid re-detecting the event it just - # handled; `handle_callbacks!` is responsible for writing them. event_last_time::Int vector_event_last_time::Int last_event_error::eventErrType @@ -204,15 +200,12 @@ mutable struct OperatorSplittingIntegrator{ last_step_failed::Bool force_stepfail::Bool isout::Bool - # DiffEqBase's callback machinery reads this as a literal field (not through - # `derivative_discontinuity!`), so the name has to match exactly. derivative_discontinuity::Bool just_hit_tstop::Bool cache::cacheType sol::solType - # Saving state. `saveiter` is the authoritative length of the saved prefix of - # `sol.t`/`sol.u`: `savevalues!` writes index `saveiter` and `postamble!` - # truncates both vectors to it. + # `saveiter`, not `length(sol.t)`, is the authoritative length of the saved + # prefix of `sol.t`/`sol.u`. const save_opts::SaveOptions saveiter::Int postamble_done::Bool @@ -269,6 +262,9 @@ function SciMLBase.__init( save_everystep = false, save_start = true, save_end = true, + # Accepted and ignored: `sol(t)` already interpolates over the saved points + # with the same (linear) interpolant a dense output would store per step. + # Named explicitly so it cannot travel down to the leaf integrators. dense = false, callback = nothing, advance_to_tstop = false, @@ -281,13 +277,6 @@ function SciMLBase.__init( (; u0, p) = prob t0, tf = prob.tspan - dense && throw( - ArgumentError( - "dense output is not supported by operator splitting integrators. The \ - interpolant is linear, so `sol(t)` between saved points is already exact \ - for it; use `saveat` or `save_everystep = true` to control the saved points." - ) - ) save_opts = SaveOptions(save_on, save_everystep, save_start, save_end) # By default every node adapts exactly if its own algorithm does; a scalar or a @@ -331,8 +320,8 @@ function SciMLBase.__init( # `build_solution` defaults `interp` to a `LinearInterpolation` over the very # vectors passed in here, so `sol(t)` sees every later `savevalues!` push and - # agrees exactly with the integrator's own (linear) interpolant. `stats` has to - # be a real `DEStats`: the callback machinery increments `sol.stats.ncondition`. + # agrees with the integrator's own (linear) interpolant. `stats` has to be a real + # `DEStats`: the callback machinery increments `sol.stats.ncondition`. sol = SciMLBase.build_solution( prob, alg, tType[], uType[]; stats = SciMLBase.DEStats(0), @@ -340,8 +329,7 @@ function SciMLBase.__init( ) callback = DiffEqBase.CallbackSet(callback) - # Only a VectorContinuousCallback needs the scratch buffers, and they have to be - # wide enough for the widest one in the set. + # The scratch buffers have to be wide enough for the widest callback in the set. max_len_cb = DiffEqBase.max_vector_callback_length_int(callback) eventErrType = real(eltype(u)) callback_cache = if max_len_cb === nothing @@ -400,27 +388,20 @@ function SciMLBase.__init( false, config, ) - # The initial save happens *after* the callbacks are initialized, so that an - # initializer which modifies `u` is reflected in the first saved point. + # Initialize first, so that an initializer which modifies `u` is reflected in the + # first saved point. initialize_callbacks!(integrator) save_initial_value!(integrator) return integrator end -""" - initialize_callbacks!(integrator::OperatorSplittingIntegrator) - -Run the callbacks' `initialize` hooks, mirroring OrdinaryDiffEqCore's -`initialize_callbacks!`. - -The flag starts out set so that an `initialize` hook can clear it -- the default -initializer does exactly that -- and a hook that did modify `u` has its change -propagated to `uprev` and to the subintegrator tree before the first step. -""" +# Run the callbacks' `initialize` hooks, mirroring OrdinaryDiffEqCore's +# `initialize_callbacks!`. The flag starts out set so that a hook can clear it; the +# default initializer does exactly that. function initialize_callbacks!(integrator::OperatorSplittingIntegrator) integrator.derivative_discontinuity = true - # Pass the integrator's own `u`, not the problem's `u0`: with `alias_u0 = false` - # they are different arrays and a hook has to be able to modify the live state. + # `integrator.u`, not `prob.u0`: with `alias_u0 = false` they are different arrays + # and a hook has to be able to modify the live state. modified = DiffEqBase.initialize!( integrator.callback, integrator.u, integrator.t, integrator ) @@ -514,19 +495,16 @@ function DiffEqBase.reinit!( ) # A leaf's `reinit!` restores it to the `u0` slice captured when it was built, not # to the `u0` handed to this call, so the new state has to be pushed down - # explicitly. Relying on the forward sync of the next step is not enough: the - # palindromic schemes skip it for their first child. + # explicitly. The forward sync of the next step is not enough: the palindromic + # schemes skip it for their first child. rollback_children!(integrator) - # After the children, so that an initializer which modifies `u` (and therefore - # resyncs the tree) is not undone by the child reinit, and so that the first saved - # point is the state the initializers left behind. + # After the children, so the child reinit cannot undo an initializer's changes. if reinit_callbacks initialize_callbacks!(integrator) end - # Saving is built into the integrator rather than provided by a saving callback, - # so the initial save is redone regardless of `reinit_callbacks`. Stale entries - # past the new `saveiter` are truncated by `postamble!`. + # Redone regardless of `reinit_callbacks`, since saving is built into the + # integrator rather than provided by a saving callback. integrator.postamble_done = false save_initial_value!(integrator) return nothing @@ -706,14 +684,13 @@ end # 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. This is also the resync path for a callback -# that modified `u`. +# are moved back to the parent's time. Also the resync path for a callback that +# modified `u`. # # Leaves get their `u` restored here as well, not only by the forward sync before the # next solve: palindromic algorithms skip that sync for their first child # (`next_sync_is_continuous`), so a rollback that left the leaf state stale would -# silently resume from the failed attempt. Restoring it here is also what re- -# establishes the `parent.u == child.u` invariant that shortcut relies on. +# silently resume from the failed attempt. rollback_children!(parent::AnySplitIntegrator) = _rollback_children!( parent.child_subintegrators, parent.child_solution_indices, parent.u, parent.t ) @@ -817,8 +794,7 @@ increment_iteration(integrator::AnySplitIntegrator) = integrator.iter += 1 function footer_reset_flags!(integrator) integrator.derivative_discontinuity = false integrator.just_hit_tstop = false - # Re-arm the postamble: stepping on after a `solve!`/`done` (e.g. after pushing a - # further tstop) has to be able to close out the solution again. + # Re-armed, so that stepping on after a `solve!` can close the solution out again. integrator.postamble_done = false return end @@ -881,8 +857,7 @@ function step_footer!(integrator::AnySplitIntegrator) try_snap_children_to_tstop!.(integrator.child_subintegrators, integrator.t) step_accept_controller!(integrator) validate_time_point(integrator) - # Callbacks run here (outer integrator only) and also perform the saving. - handle_callbacks!(integrator) + handle_callbacks!(integrator) # also does the saving elseif integrator.force_stepfail # Failure escalation protocol: the failing node's own adaptivity decides. fatal_rc = _fatal_child_retcode(integrator.child_subintegrators) @@ -1058,11 +1033,9 @@ function SciMLBase.check_error!(integrator::SplitSubIntegrator) return code end -# SciMLBase's generic `check_error!` finalizes the integrator whenever the code is -# anything but `Success`, but a mid-solve node legitimately reports `Default`: the -# stepping loops call this before every step, so the generic version would run the -# postamble -- callback finalizers, the endpoint save -- once per step. Only a real -# failure finalizes here. +# Only a real failure finalizes here. The generic method finalizes on anything but +# `Success`, but a mid-solve node legitimately reports `Default`, and this is called +# before every step -- so the postamble would run once per step. function SciMLBase.check_error!(integrator::OperatorSplittingIntegrator) code = SciMLBase.check_error(integrator) # Rebuilding the solution allocates, and this runs before every step. @@ -1129,11 +1102,10 @@ function splitting_interpolant!( return out end -# Step-local coordinate `Θ` and the length of the step that was just taken. -# `integrator.dt` must not be used for the latter: by the time a step is finished -# `step_accept_controller!` has already overwritten it with the *next* proposal, so it -# is not the length of the interval `uprev`/`u` span. A zero-length interval only -# happens before the first step, where `uprev == u` and every Θ gives the same value. +# Step-local coordinate `Θ` and the length of the step that was just taken. The latter +# must not come from `integrator.dt`: once a step is accepted `step_accept_controller!` +# has overwritten it with the *next* proposal. A zero-length interval only happens +# before the first step, where `uprev == u` and every Θ gives the same value. function _interp_coords(integrator::OperatorSplittingIntegrator, t) dt = integrator.t - integrator.tprev Θ = iszero(dt) ? zero(t / oneunit(dt)) : (t - integrator.tprev) / dt @@ -1172,8 +1144,7 @@ Move the integrator back to a time `t` inside the step that was just taken. `u` is refilled from the step's interpolant and the whole subintegrator tree is re-anchored to `t` through [`rollback_children!`](@ref), so every child's state and clock stay consistent with the parent (`validate_time_point` asserts the latter). -`integrator.dt` is deliberately left alone: it already holds the step size proposed -for the next step, and the interpolation coordinate comes from `t - tprev`. +`integrator.dt` is deliberately left alone: it already holds the next step's proposal. """ function SciMLBase.change_t_via_interpolation!( integrator::OperatorSplittingIntegrator, t, @@ -1199,8 +1170,7 @@ end # Saving # # Outer integrator only, for the reason given in the "Dense output" section of the -# developer documentation. `saveiter` is the authoritative length of the saved prefix -# of `sol.t`/`sol.u`. +# developer documentation. # --------------------------------------------------------------------------- _saved_at_current_t(integrator::OperatorSplittingIntegrator) = @@ -1225,11 +1195,10 @@ Save the state at the current time point and at any pending `saveat` points that step just taken has passed, returning `(saved, savedexactly)`. `saveat` points strictly inside the step are filled from the step's interpolant -(see [`splitting_interpolant`](@ref)), so asking for output never changes the -sequence of steps and therefore never changes the splitting error. Because the -generic interpolant is linear, interpolated points are first order accurate even -when the splitting scheme is second order; pass the times as `tstops` instead to -have the integrator land on them exactly. +(see [`splitting_interpolant`](@ref)), so asking for output never changes the sequence +of steps and therefore never changes the splitting error. The generic interpolant is +linear, so such points are only first order accurate; pass the times as `tstops` +instead to have the integrator land on them exactly. """ function SciMLBase.savevalues!( integrator::OperatorSplittingIntegrator, @@ -1255,7 +1224,6 @@ function SciMLBase.savevalues!( savedexactly = true _save_current!(integrator) else - # The interpolant allocates a fresh `u`, so hand over ownership. _save_at!(integrator, curt, integrator(curt), false) end saved = true @@ -1273,7 +1241,6 @@ function SciMLBase.savevalues!( return saved, savedexactly end -# Make sure the final time point is in the solution. function save_endpoint!(integrator::OperatorSplittingIntegrator) integrator.save_opts.save_on || return nothing integrator.save_opts.save_end || return nothing @@ -1288,8 +1255,8 @@ end # Outer integrator only, for the reason given in the "Dense output" section of the # developer documentation. Root finding, `save_positions` handling and `affect!` # dispatch are DiffEqBase's; this is just the per-step driver, modelled on -# OrdinaryDiffEqCore's `handle_callbacks!`. It also owns the plain `saveat` save -# when no callback saved. +# OrdinaryDiffEqCore's `handle_callbacks!`, and it also does the saving when no +# callback saved. # --------------------------------------------------------------------------- handle_callbacks!(::SplitSubIntegrator) = nothing @@ -1317,8 +1284,8 @@ function handle_callbacks!(integrator::OperatorSplittingIntegrator) # controller's error history no longer describes what happens next. reinit_node_controller!(integrator) else - # Clearing these is what stops the root finder from nudging `tprev` on a - # step that had no event. + # Clearing these stops the root finder from nudging `tprev` and + # re-detecting the event handled on an earlier step. integrator.event_last_time = 0 integrator.vector_event_last_time = 1 end @@ -1502,8 +1469,7 @@ function SciMLBase.postamble!(integrator::OperatorSplittingIntegrator) integrator.postamble_done = true DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator) save_endpoint!(integrator) - # `saveiter` is authoritative: drop whatever a previous, longer run left behind - # (a `reinit!` without `erase_sol`). + # Drop whatever a previous, longer run left behind (`reinit!` without `erase_sol`). resize!(integrator.sol.t, integrator.saveiter) resize!(integrator.sol.u, integrator.saveiter) return nothing @@ -1803,7 +1769,6 @@ SciMLBase.pop_tstop!(i::AnySplitIntegrator) = pop!(i.tstops) DiffEqBase.get_dt(i::AnySplitIntegrator) = i.dt -# Continuous callbacks relax the step size through this after an event. function SciMLBase.set_proposed_dt!(integrator::OperatorSplittingIntegrator, dt) if integrator.dtcache != abs(dt) integrator.dtcache = abs(dt) diff --git a/src/utils.jl b/src/utils.jl index 8f4b469..ab3e99d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -7,16 +7,13 @@ function tstops_and_saveat_heaps(t0, tf, tstops, saveat) tstops = [filter(t -> t0 < t < tf || tf < t < t0, tstops)..., tf] tstops = BinaryHeaps.BinaryHeap{FT, ordering}(tstops) - # Keep `t0 < t <= tf` in tdir-space, as OrdinaryDiffEqCore's `initialize_saveat` - # does: the initial point is owned by `save_start` and the final one by - # `save_end`, so leaving either in the heap would duplicate it. + # Keep `t0 < t <= tf` in tdir-space: `save_start` owns the initial point and + # `save_end` the final one, so leaving either in the heap would duplicate it. tdir = tf > t0 ? one(FT) : -one(FT) saveat = if isnothing(saveat) FT[] elseif saveat isa Number saveat > zero(saveat) || error("saveat value must be positive") - # The range already satisfies the window, so no filtering is needed: it - # starts one step past `t0` and hits `tf` only if the span divides evenly. step = tdir * saveat collect((t0 + step):step:tf) else diff --git a/test/callbacks.jl b/test/callbacks.jl index 6db104c..3ef3fe9 100644 --- a/test/callbacks.jl +++ b/test/callbacks.jl @@ -54,14 +54,11 @@ exact_crossing(level) = -10 * log(level / 3) calls = Ref(0) cb = DiscreteCallback((u, t, integrator) -> (calls[] += 1; false), integrator -> nothing) DiffEqBase.solve!(DiffEqBase.init(nprob, alg; dt = 0.1, callback = cb)) - # Ten outer steps, regardless of how many substeps the inner node takes. - @test calls[] == 10 + @test calls[] == 10 # outer steps only, whatever the inner node does end @testset "affect! fires and sees the integrator" begin seen = Float64[] - # Threshold between grid points, compared with isapprox: see the CallbackSet - # testset below for why the nominal grid cannot be trusted exactly. cb = DiscreteCallback( (u, t, integrator) -> t >= 0.45, integrator -> push!(seen, integrator.t) @@ -97,9 +94,8 @@ exact_crossing(level) = -10 * log(level / 3) end @testset "a modified u reaches the whole subintegrator tree" begin - # This is the operator-splitting-specific hazard: every child holds its own - # copy of its slice, and the palindromic schemes skip the forward sync of - # their first child on the assumption that it still matches the parent. + # Every child holds its own copy of its slice, and the palindromic schemes + # skip the forward sync of their first child. for alg in ( ltg(), StrangMarchuk((Euler(), Euler())), @@ -138,9 +134,8 @@ exact_crossing(level) = -10 * log(level / 3) end @testset "multiple callbacks in a CallbackSet" begin - # Thresholds sit between grid points so that ulp drift in the step times - # (0.1 * 8 is 0.7999999999999999, not 0.8) cannot flip a condition, and the - # recorded times are compared with isapprox for the same reason. + # Thresholds sit between grid points, and times are compared with isapprox, + # because of ulp drift (0.1 * 8 is 0.7999999999999999, not 0.8). early = Float64[] late = Float64[] cb1 = DiscreteCallback((u, t, integrator) -> t >= 0.35, integrator -> push!(early, integrator.t)) @@ -208,10 +203,10 @@ end @test sol.retcode == ReturnCode.Success end - @testset "event time converges with the interpolant" begin - # The interpolant is linear, so the located root is second order accurate. - # (The step grid moves with dt, so the constant is noisy; a factor-of-four - # improvement over a factor-of-four refinement is the honest assertion.) + @testset "event time converges under refinement" begin + # The located root is at least second order accurate. The step grid moves + # with dt, so the constant is noisy: assert only that a factor-of-four + # refinement buys a factor of four. errs = map((0.1, 0.025)) do dt hits = Float64[] cb = ContinuousCallback( @@ -249,8 +244,6 @@ end while integrator.t < exact DiffEqBase.step!(integrator) end - # `change_t_via_interpolation!` moved the step back to the event; every child - # clock and buffer must have followed. @test integrator.t ≈ exact atol = 1.0e-3 OS.validate_time_point(integrator) @test integrator.child_subintegrators[1].u ≈ integrator.u @@ -290,10 +283,9 @@ end (u, t, integrator) -> u[3] - level, integrator -> push!(hits, integrator.t) ) - # `dtmax` matters here: this problem is so nearly linear that the controller - # otherwise grows the step to span the rest of the tspan in one go, and the - # event can only ever be as accurate as the interpolant over the bracketing - # step -- which is linear. + # `dtmax` matters: this problem is so nearly linear that the controller + # otherwise spans the rest of the tspan in one step, and the event is only + # as accurate as the (linear) interpolant over the bracketing step. sol = DiffEqBase.solve!( DiffEqBase.init( prob, PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())); @@ -305,36 +297,25 @@ end @test sol.retcode == ReturnCode.Success end - @testset "a long adaptive step degrades the event time, not the bracket" begin - # Documented consequence of the linear interpolant: the located event is the - # exact root of the interpolant over whatever step brackets it. Pin that - # relationship so the behaviour is a stated property rather than a surprise. - bracket = Ref((0.0, 0.0)) - hits = Float64[] - cb = ContinuousCallback( - (u, t, integrator) -> u[3] - level, - function (integrator) - push!(hits, integrator.t) - return bracket[] = (integrator.tprev, integrator.t) - end - ) - DiffEqBase.solve!( - DiffEqBase.init( - prob, PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())); - dt = 0.05, callback = cb + @testset "a long adaptive step degrades the event time" begin + # An event is only as accurate as the interpolant over the step that brackets + # it, so capping the step with dtmax improves it. + function locate(; kwargs...) + hits = Float64[] + cb = ContinuousCallback( + (u, t, integrator) -> u[3] - level, + integrator -> push!(hits, integrator.t) ) - ) - t0 = bracket[][1] - # u[3] is solved exactly by the first operator alone. - u3(t) = 3 * exp(-0.1t) - # The step originally ended where the event time now sits plus the part of the - # step that was cut off; recover it from the interpolant relation instead. - @test length(hits) == 1 - @test hits[1] > t0 - # The event solves the *linear* interpolant, so g is zero there by - # construction while the true crossing is elsewhere. - @test !isapprox(hits[1], exact; atol = 1.0e-3) - @test u3(hits[1]) < level # the linear chord under-estimates a convex decay + DiffEqBase.solve!( + DiffEqBase.init( + prob, PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())); + dt = 0.05, callback = cb, kwargs... + ) + ) + return only(hits) + end + + @test abs(locate(dtmax = 0.05) - exact) < abs(locate() - exact) end @testset "interp_points = 0 skips the safety sweep" begin @@ -351,8 +332,8 @@ end integrator -> nothing ) DiffEqBase.solve!(DiffEqBase.init(prob, ltg_exact(); dt = 0.25, callback = cb2)) - # The default `interp_points = 11` samples the (collinear) interpolant many - # more times without being able to find anything the endpoints missed. + # The default `interp_points = 11` sweeps the interpolant for sign changes + # the endpoints missed; with `0` the condition is only seen at the endpoints. @test calls[] < sweeping[] end end @@ -362,7 +343,7 @@ end levels = (2.9, 2.85) # Current DiffEqBase hands `affect!` a mask of simultaneous events rather than a - # single index; accept either shape so this does not depend on that detail. + # single index; accept either shape. fired_index(idx::Integer) = Int(idx) fired_index(mask) = only(findall(!iszero, mask)) diff --git a/test/saving.jl b/test/saving.jl index 8fc8435..bc545f3 100644 --- a/test/saving.jl +++ b/test/saving.jl @@ -28,11 +28,24 @@ fsplit = GenericSplitFunction((f1, f2), ([1, 2, 3], [1, 2])) make_prob(tspan = TSPAN) = OperatorSplittingProblem(fsplit, copy(U0), tspan) ltg() = LieTrotterGodunov((Euler(), Euler())) +# Swept by every test that is not about one particular scheme. +algs() = ( + LieTrotterGodunov((Tsit5(), Tsit5())), + StrangMarchuk((Tsit5(), Tsit5())), + PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())), +) + +function exact(t) + p = 1.5 * exp(-0.11t) # (1, 1) mode + m = -0.5 * exp(-0.09t) # (1, -1) mode + return [p + m, p - m, 3 * exp(-0.1t)] +end +rhs(u) = [-0.1u[1] - 0.01u[2], -0.1u[2] - 0.01u[1], -0.1u[3]] + @testset "saved time points" begin prob = make_prob() @testset "default saves only the interval endpoints" begin - # `save_everystep` defaults to false, so only `save_start`/`save_end` fire. sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1)) @test sol.t == [0.0, 1.0] @test length(sol.u) == 2 @@ -107,8 +120,7 @@ ltg() = LieTrotterGodunov((Euler(), Euler())) end @testset "asking for output does not change the trajectory" begin - # saveat points are interpolated, never stepped onto, so the sequence of - # steps -- and hence the splitting error -- is untouched. + # saveat points are interpolated, never stepped onto. plain = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1)) dense = DiffEqBase.solve!( DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = 0.017, save_everystep = true) @@ -117,9 +129,7 @@ ltg() = LieTrotterGodunov((Euler(), Euler())) end @testset "saveat times passed as tstops are landed on, not interpolated" begin - # The documented way to get output that is accurate to the order of the - # scheme: the integrator steps exactly onto the time, so the saved value is - # the stepped state rather than the (first order) interpolant. + # The documented way to get output accurate to the order of the scheme. ts = [0.25, 0.55] landed = DiffEqBase.solve!( DiffEqBase.init(prob, ltg(); dt = 0.1, tstops = ts, saveat = ts) @@ -133,8 +143,7 @@ ltg() = LieTrotterGodunov((Euler(), Euler())) @test t in landed.t @test t in interpolated.t end - # Same requested output times, two different values: one is the stepped - # state, the other the interpolant of the step that straddles the time. + # Same times, different values: stepped state vs. interpolated one. for t in ts a = landed.u[findfirst(==(t), landed.t)] b = interpolated.u[findfirst(==(t), interpolated.t)] @@ -158,11 +167,7 @@ ltg() = LieTrotterGodunov((Euler(), Euler())) end @testset "other algorithms" begin - for alg in ( - StrangMarchuk((Euler(), Euler())), - StrangMarchuk((Tsit5(), Tsit5())), - PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())), - ) + for alg in algs() sol = DiffEqBase.solve!( DiffEqBase.init(prob, alg; dt = 0.1, saveat = [0.5]) ) @@ -189,82 +194,38 @@ end @testset "interpolation" begin prob = make_prob() - @testset "reproduces the step endpoints" begin - integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) + @testset "$(nameof(typeof(alg)))" for alg in algs() + integrator = DiffEqBase.init(prob, alg; dt = 0.1) DiffEqBase.step!(integrator) - @test integrator(integrator.tprev) ≈ integrator.uprev - @test integrator(integrator.t) ≈ integrator.u - end - - @testset "is linear in between" begin - integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) DiffEqBase.step!(integrator) - (; tprev, t, uprev, u) = integrator - tmid = tprev + (t - tprev) / 3 - Θ = (tmid - tprev) / (t - tprev) - @test integrator(tmid) ≈ @. (1 - Θ) * uprev + Θ * u - end + (; tprev, t) = integrator + tmid = (tprev + t) / 2 - @testset "first derivative is the step slope" begin - integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) - DiffEqBase.step!(integrator) - (; tprev, t, uprev, u) = integrator - expected = (u .- uprev) ./ (t - tprev) - @test integrator(tprev + (t - tprev) / 2, Val{1}) ≈ expected - # Constant slope: the derivative does not depend on where it is evaluated. - @test integrator(t, Val{1}) ≈ expected - end + @test integrator(tprev) ≈ integrator.uprev + @test integrator(t) ≈ integrator.u + @test integrator(tmid) ≈ exact(tmid) rtol = 1.0e-2 + @test integrator(tmid, Val{1}) ≈ rhs(exact(tmid)) rtol = 1.0e-2 + # For an adaptive scheme `dt` is already the next proposal, so the endpoints + # above are only reproduced if the interpolant does not key off it. + SciMLBase.isadaptive(alg) && @test integrator.dt != t - tprev - @testset "in-place and idxs" begin - integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) - DiffEqBase.step!(integrator) - tmid = (integrator.tprev + integrator.t) / 2 out = similar(integrator.u) integrator(out, tmid) @test out ≈ integrator(tmid) - @test integrator(tmid, Val{0}; idxs = 2) ≈ integrator(tmid)[2] - @test integrator(tmid, Val{0}; idxs = [1, 3]) ≈ integrator(tmid)[[1, 3]] - sub = similar(integrator.u, 2) integrator(sub, tmid, Val{0}; idxs = [1, 3]) @test sub ≈ integrator(tmid)[[1, 3]] - end - - @testset "does not read the stale dt" begin - # After an accepted step the controller has already replaced `dt` with the - # next proposal, so an interpolant keyed off `integrator.dt` would be wrong. - integrator = DiffEqBase.init( - prob, PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())); dt = 0.1 - ) - DiffEqBase.step!(integrator) - DiffEqBase.step!(integrator) - @test integrator.dt != integrator.t - integrator.tprev - @test integrator(integrator.t) ≈ integrator.u - @test integrator(integrator.tprev) ≈ integrator.uprev - end - - @testset "sol(t) matches the integrator interpolant" begin - integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, save_everystep = true) - sol = DiffEqBase.solve!(integrator) - # Between two saved step endpoints `sol.interp` is the same linear - # interpolant the integrator uses within a step. - i = 4 - t0, t1 = sol.t[i], sol.t[i + 1] - tmid = (t0 + t1) / 2 - @test sol(tmid) ≈ (sol.u[i] .+ sol.u[i + 1]) ./ 2 - end - @testset "saveat values are the step interpolant" begin - tmid = 0.05 # strictly inside the first step of length 0.1 - integrator = DiffEqBase.init( - prob, ltg(); dt = 0.1, saveat = [tmid], save_everystep = true + sol = DiffEqBase.solve!( + DiffEqBase.init( + prob, alg; dt = 0.1, saveat = [0.05], save_everystep = true + ) ) - sol = DiffEqBase.solve!(integrator) - idx = findfirst(==(tmid), sol.t) - @test idx !== nothing - # Bracketing saved points are the enclosing step's endpoints. - @test sol.u[idx] ≈ (sol.u[idx - 1] .+ sol.u[idx + 1]) ./ 2 + @test sol.u[findfirst(==(0.05), sol.t)] ≈ exact(0.05) rtol = 1.0e-2 + for τ in (0.35, 0.62, 0.97) + @test sol(τ) ≈ exact(τ) rtol = 1.0e-2 + end end end @@ -281,7 +242,6 @@ end SciMLBase.change_t_via_interpolation!(integrator, tmid) @test integrator.t == tmid @test integrator.u ≈ expected - # Every child clock and buffer follows the parent. OS.validate_time_point(integrator) for child in integrator.child_subintegrators @test child.t == tmid @@ -308,8 +268,7 @@ end end @testset "accepts the Val{:false} the callback path passes" begin - # DiffEqBase's continuous callback path calls this with `Val{:false}` -- a - # `Val` of the Symbol -- so the flag must be decoded by dispatch. + # A `Val` of the *Symbol*, so the flag has to be decoded by dispatch. integrator = DiffEqBase.init(prob, ltg(); dt = 0.1) DiffEqBase.step!(integrator) tmid = (integrator.tprev + integrator.t) / 2 @@ -331,8 +290,16 @@ end @testset "solution object plumbing" begin prob = make_prob() - @testset "dense output is rejected" begin - @test_throws ArgumentError DiffEqBase.init(prob, ltg(); dt = 0.1, dense = true) + @testset "dense is accepted and does not reach the leaves" begin + integrator = DiffEqBase.init( + prob, ltg(); dt = 0.1, dense = true, save_everystep = true + ) + sol = DiffEqBase.solve!(integrator) + @test sol.retcode == ReturnCode.Success + @test sol(0.35) ≈ (sol.u[4] .+ sol.u[5]) ./ 2 + # `dense` must be swallowed by `__init`; a leaf handed `dense = true` would + # store interpolation data for stages that are not solution points. + @test integrator.child_subintegrators[1].sol.dense == false end @testset "sol.stats is a real DEStats" begin @@ -389,8 +356,7 @@ end @testset "a new u0 reaches the whole tree" begin # Regression: a leaf's `reinit!` restores the `u0` slice captured when it was # built, and the palindromic schemes skip the forward sync of their first - # child -- so `reinit!(integrator, u0)` used to integrate a stale child state - # and silently return the wrong answer. + # child, so this used to integrate a stale child state. u0b = [10.0, 20.0, 30.0] probb = OperatorSplittingProblem(fsplit, copy(u0b), TSPAN) for alg in ( @@ -420,8 +386,7 @@ end end @testset "callbacks can be attached" begin - # Regression: DiffEqBase's `initialize!` reads the `derivative_discontinuity` - # field directly, so `init` used to throw as soon as any callback was passed. + # Regression: `init` used to throw as soon as any callback was passed. prob = make_prob() calls = Ref(0) cb = DiscreteCallback((u, t, integrator) -> false, integrator -> (calls[] += 1)) From 6ba417d8f0c0588f6b308427b3e1f42cb5e4a26c Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 4 Aug 2026 10:27:45 -0400 Subject: [PATCH 4/4] bump SciMLBase --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 5602f44..a83ae67 100644 --- a/Project.toml +++ b/Project.toml @@ -27,7 +27,7 @@ OrdinaryDiffEqTsit5 = "1.1.0, 2" PrecompileTools = "1.1" RecursiveArrayTools = "3.39.0, 4" SafeTestsets = "0.1.0" -SciMLBase = "3.1" +SciMLBase = "3.36" SciMLIterators = "1" SciMLTesting = "2.1" SymbolicIndexingInterface = "0.3.36"