From 7a034609631c392edd93b373b9a5ae7eda631dc4 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 4 Aug 2026 16:30:32 -0400 Subject: [PATCH 1/6] Add order 3+ splitting schemes with adjoint-pair adaptivity Implements the practical parts of Auzinger, Hofstaetter, Ketcheson & Koch, BIT 57:55-74 (2017): coefficient-table splitting schemes and the section 3 adjoint-pair local error estimator. Backward substeps. Every real splitting scheme of order three or above has negative coefficients, so some substeps run against the tree's direction -- previously unsupported, because a child's direction is fixed at construction and the stepping loop re-signs dt back before perform_step! sees it. reverse_direction! flips a leaf's tdir, dt, dtcache, dtpropose and dtmax and rebuilds its tdir-keyed heaps. A splitting node needs different handling: its heap stores raw times with the direction in the heap's type, so keys cannot be re-signed and only behind-times are dropped, and the reversal recurses into the subtree because add_tstop! propagates to descendants eagerly. Coefficient tables. SplittingCoefficients holds an s-by-N table as nested NTuples, checking each operator's consistency condition at construction. Ruth3 is the first table (exactly rational, order 3). One generic _perform_step! walks stages by operators and skips zero coefficients. LieTrotterGodunov and StrangMarchuk keep their hand-written steps. AdjointPair(base) averages a base scheme of odd order p with its adjoint for a solution of order p+1, and takes half their difference as the error estimate driving the controller. Even-order bases are rejected: their leading error terms are equal rather than opposite, so averaging would raise no order and the difference would not be an error estimate. The adjoint reverses the base's whole flat sequence of flows, reusing the same table. AdjointPair(LieTrotterGodunov(...)) is PalindromicPairLieTrotterGodunov, which a test asserts to the bit. Also fixes a pre-existing bug the new schemes exposed: the snap window absorbing child clock drift was scaled by the child's own dt, which after a substep can be arbitrarily small and collapse the window below the drift it must absorb. It now uses the outer step size, the scale the drift accumulated against. Co-Authored-By: Claude Opus 5 --- docs/src/api-reference/index.md | 2 + src/OrdinaryDiffEqOperatorSplitting.jl | 2 +- src/integrator.jl | 48 ++-- src/solver.jl | 294 ++++++++++++++++++++++++- src/utils.jl | 80 +++++++ test/adjoint_pair.jl | 87 ++++++++ test/callbacks.jl | 1 + test/coefficient_schemes.jl | 72 ++++++ test/convergence.jl | 32 ++- test/saving.jl | 2 + 10 files changed, 592 insertions(+), 28 deletions(-) create mode 100644 test/adjoint_pair.jl create mode 100644 test/coefficient_schemes.jl diff --git a/docs/src/api-reference/index.md b/docs/src/api-reference/index.md index 44b0610..5dfdc4d 100644 --- a/docs/src/api-reference/index.md +++ b/docs/src/api-reference/index.md @@ -18,6 +18,8 @@ GenericSplitFunction LieTrotterGodunov StrangMarchuk PalindromicPairLieTrotterGodunov +Ruth3 +AdjointPair ``` ## Per-node configuration diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index 661dc24..4af8b30 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -154,7 +154,7 @@ include("solver.jl") include("utils.jl") export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, StrangMarchuk, - PalindromicPairLieTrotterGodunov + PalindromicPairLieTrotterGodunov, Ruth3, AdjointPair export SplitNode, TreeOption include("precompilation.jl") diff --git a/src/integrator.jl b/src/integrator.jl index e5befd1..04f36ea 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -816,7 +816,9 @@ function fixed_t_for_floatingpoint_error!(integrator::AnySplitIntegrator, ttmp) return if DiffEqBase.has_tstop(integrator) tstop = DiffEqBase.first_tstop(integrator) if abs(ttmp - tstop) < _snap_window(integrator.t, tstop, integrator.dt) - try_snap_children_to_tstop!.(integrator.child_subintegrators, tstop) + try_snap_children_to_tstop!.( + integrator.child_subintegrators, tstop, integrator.dt + ) tstop else ttmp @@ -825,16 +827,20 @@ function fixed_t_for_floatingpoint_error!(integrator::AnySplitIntegrator, ttmp) ttmp end end -function try_snap_children_to_tstop!(integrator::SplitSubIntegrator, tstop) - if abs(tstop - integrator.t) < _snap_window(integrator.t, tstop, integrator.dt) +# `scale` is the *outer* step size, and stays the outer one all the way down the tree. +# All the drift being absorbed here was accumulated by subdividing that one step, so it +# is the right yardstick; a child's own `dt` is whatever fraction it last landed on and +# can be arbitrarily smaller, collapsing the window below the drift it has to absorb. +function try_snap_children_to_tstop!(integrator::SplitSubIntegrator, tstop, scale) + if abs(tstop - integrator.t) < _snap_window(integrator.t, tstop, scale) integrator.t = tstop else @warn "Failed to snap timestep for integrator $(integrator.t) with parent integrator hitting the tstop $(tstop)." end - return try_snap_children_to_tstop!.(integrator.child_subintegrators, tstop) + return try_snap_children_to_tstop!.(integrator.child_subintegrators, tstop, scale) end -function try_snap_children_to_tstop!(integrator::DEIntegrator, tstop) - return if abs(tstop - integrator.t) < _snap_window(integrator.t, tstop, integrator.dt) +function try_snap_children_to_tstop!(integrator::DEIntegrator, tstop, scale) + return if abs(tstop - integrator.t) < _snap_window(integrator.t, tstop, scale) integrator.t = tstop else @warn "Failed to snap timestep for integrator $(integrator.t) with parent integrator hitting the tstop $(tstop)." @@ -854,7 +860,9 @@ function step_footer!(integrator::AnySplitIntegrator) # halves) accumulate ulp-level drift from the parent's exact `t`. # Re-anchor children to the parent's canonical time here so the # drift cannot accumulate across outer steps. - try_snap_children_to_tstop!.(integrator.child_subintegrators, integrator.t) + try_snap_children_to_tstop!.( + integrator.child_subintegrators, integrator.t, integrator.dt + ) step_accept_controller!(integrator) validate_time_point(integrator) handle_callbacks!(integrator) # also does the saving @@ -1616,11 +1624,10 @@ function advance_solution_by!( return end -# `dt` stays signed through the splitting tree, following the SciML `step!` -# convention: the sign has to match the child's own integration direction, which -# equals the tree's. Advancing a child *against* its direction (negative substeps -# of higher order compositions) is not supported yet: leaf ODEIntegrators cannot -# step against the tdir their tspan fixed at construction. +# `dt` stays signed through the splitting tree, following the SciML `step!` convention: +# the sign has to match the child's own integration direction. Schemes of order three +# and above have negative coefficients, so a substep can oppose the tree's direction, +# and a child's direction is fixed at construction -- hence `reverse_direction!`. # Recursion dispatch function advance_solution_by!( @@ -1628,13 +1635,26 @@ function advance_solution_by!( sub::SplitSubIntegrator, dt ) - SciMLBase.step!(sub, dt, true) + _step_signed!(sub, dt) return nothing end # Leaf dispatch function advance_solution_by!(outer::AnySplitIntegrator, child::DEIntegrator, dt) - SciMLBase.step!(child, dt, true) + _step_signed!(child, dt) + return nothing +end + +# Reversing first is also what satisfies `step!`'s own direction guard: once `tdir` is +# flipped, a negative `dt` agrees with it. +function _step_signed!(child, dt) + if child.tdir * dt < zero(dt) + reverse_direction!(child) + SciMLBase.step!(child, dt, true) + reverse_direction!(child) + else + SciMLBase.step!(child, dt, true) + end return nothing end diff --git a/src/solver.jl b/src/solver.jl index db76d35..6389a25 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -221,7 +221,7 @@ function init_cache( return PalindromicPairLieTrotterGodunovCache(u, uprev, similar(u)) end -function _ppltg_advance_child!(parent, child, i, dt) +function _advance_child!(parent, child, i, dt) idxs = parent.child_solution_indices[i] sync = parent.child_synchronizers[i] @@ -241,7 +241,7 @@ end i = 0 @unroll for child in children i += 1 - _ppltg_advance_child!(parent, child, i, dt) + _advance_child!(parent, child, i, dt) parent.force_stepfail && return end end @@ -251,7 +251,7 @@ end j = 0 @unroll for child in rchildren j += 1 - _ppltg_advance_child!(parent, child, N + 1 - j, dt) + _advance_child!(parent, child, N + 1 - j, dt) parent.force_stepfail && return end end @@ -286,3 +286,291 @@ function _perform_step!( end return end + +# --------------------------------------------------------------------------- +# Coefficient-table splitting schemes +# --------------------------------------------------------------------------- +""" + SplittingCoefficients(stages::NTuple{N, T}...) + +Coefficients of an `S`-stage splitting scheme over `N` operators, one tuple per stage: +stage `j` advances operator `i` by `stages[j][i] * dt`. + +This is the generalization to `N` operators of the two-operator (`AB`) and +three-operator (`ABC`) coefficient tables of +[AuzHofKetKoc:2017:psm](@cite); their tables are the `N = 2` and `N = 3` cases. + +Each operator's coefficients must sum to one, the consistency condition, and that is +checked here. The remaining order conditions are not, so a table that constructs +successfully can still fail to attain the order it claims. +""" +struct SplittingCoefficients{S, N, T} + a::NTuple{S, NTuple{N, T}} + + # Do not simplify this to `NTuple{S, NTuple{N, T}}`: the empty tuple matches it for + # any element type, leaving parameters unbound for `S == 0` (and `T` unbound for + # `N == 0`). A leading element plus a counted `Vararg` rules both out. + function SplittingCoefficients( + stage1::Tuple{T, Vararg{T, K}}, + rest::Tuple{T, Vararg{T, K}}... + ) where {T, K} + a = (stage1, rest...) + N = K + 1 + S = length(a) + for i in 1:N + total = sum(a[j][i] for j in 1:S) + total ≈ one(T) || throw( + ArgumentError( + "operator $i's coefficients sum to $total rather than 1, so the \ + scheme is not consistent." + ) + ) + end + return new{S, N, T}(a) + end +end + +""" + coefficients(alg) + +The [`SplittingCoefficients`](@ref) table of a coefficient-driven splitting algorithm. +""" +function coefficients end + +""" + order(alg) + +Order of consistency of a splitting algorithm, counting the splitting error only. +""" +function order end + +""" + Ruth3 <: AbstractOperatorSplittingAlgorithm + +Third-order splitting scheme of [Rut:1983:cim](@cite), in three stages. + +Its coefficients are exactly rational, and -- as is unavoidable for any real +splitting scheme of order three or above -- some of them are negative, so parts of +the step run backward in time. + +As for every splitting scheme here the order statement covers the *splitting* error +only, and presumes the inner solvers resolve their subproblems accurately relative +to it. +""" +struct Ruth3{AlgTupleType <: Tuple} <: AbstractOperatorSplittingAlgorithm + inner_algs::AlgTupleType + + function Ruth3(inner_algs::Tuple) + n = length(inner_algs) + n == 2 || throw( + ArgumentError( + "Ruth3 is a two-operator (AB) table but got $n operators. Group the \ + operators into a nested GenericSplitFunction to use it with more." + ) + ) + return new{typeof(inner_algs)}(inner_algs) + end +end + +function Base.show(io::IO, alg::Ruth3) + print(io, "Ruth3 (") + for inner_alg in alg.inner_algs[1:(end - 1)] + Base.show(io, inner_alg) + print(io, " -> ") + end + length(alg.inner_algs) > 0 && Base.show(io, alg.inner_algs[end]) + return print(io, ")") +end + +const RUTH3_COEFFICIENTS = SplittingCoefficients( + (7 // 24, 2 // 3), (3 // 4, -2 // 3), (-1 // 24, 1 // 1) +) + +coefficients(::Ruth3) = RUTH3_COEFFICIENTS +order(::Ruth3) = 3 + +order(::StrangMarchuk) = 2 +order(::PalindromicPairLieTrotterGodunov) = 2 + +# Lie-Trotter keeps its hand-written step; the table exists only so it can serve as an +# `AdjointPair` base, which is what makes `AdjointPair(LieTrotterGodunov(...))` and +# `PalindromicPairLieTrotterGodunov` the same scheme. +order(::LieTrotterGodunov) = 1 +coefficients(alg::LieTrotterGodunov) = + SplittingCoefficients(ntuple(_ -> 1 // 1, length(alg.inner_algs))) + +struct SplittingCoefficientsCache{uType, uprevType, coeffType} <: AbstractOperatorSplittingCache + u::uType + uprev::uprevType + coeffs::coeffType +end + +function init_cache( + f::GenericSplitFunction, alg::Ruth3; + uprev::AbstractArray, u::AbstractVector, + ) + return SplittingCoefficientsCache(u, uprev, coefficients(alg)) +end + +function _perform_step!( + parent, + children::Tuple, + cache::SplittingCoefficientsCache, + dt + ) + # Deliberately no `mark_next_sync_continuous`: that shortcut needs the previous step + # to have left `parent.u` equal to the buffer of the child solved first, which holds + # for StrangMarchuk only because its reverse pass ends on operator 1. A general table + # ends on operator N, so operator 1's buffer is stale and skipping its forward sync + # resumes from stale state -- the first step stays exact and every later one is wrong. + _table_stages!(parent, children, cache.coeffs.a, dt) + return +end + +@unroll function _table_stages!(parent, children, stages::Tuple, dt) + @unroll for stage in stages + _table_stage!(parent, children, stage, dt) + parent.force_stepfail && return + end +end + +@unroll function _table_stage!(parent, children::Tuple, stage, dt) + i = 0 + @unroll for child in children + i += 1 + coefficient = stage[i] + # A zero coefficient is the identity flow, so it is skipped sync and all. + if !iszero(coefficient) + _advance_child!(parent, child, i, coefficient * dt) + parent.force_stepfail && return + end + end +end + +# --------------------------------------------------------------------------- +# Adjoint pairs +# --------------------------------------------------------------------------- +""" + AdjointPair(base) <: AbstractOperatorSplittingAlgorithm + +Adaptive splitting scheme of order `p+1` built from a `base` scheme of odd order `p`, +following [AuzHofKetKoc:2017:psm](@cite), eq. (3.2). + +One step runs `base` and its adjoint ``\\mathcal{S}^*`` from the same initial value. +Their leading error terms are ``C h^{p+1}`` and ``(-1)^p C h^{p+1}``, so for odd `p` +the signs oppose: the average is a solution of order `p+1`, and half the difference is +an asymptotically correct estimate of the base scheme's local error, which drives the +step size controller. A step therefore costs twice the base scheme. + +The base scheme's order must be **odd**. For even `p` the two leading terms are +*equal* rather than opposite, so averaging cancels nothing and the difference stops +being an error estimate. + +``\\mathcal{S}^*(h, u) = \\mathcal{S}^{-1}(-h, u)`` is the base scheme's entire flat +sequence of flows reversed, every coefficient keeping its sign and its operator, so it +reuses the same table and needs no extra coefficients. + +[`PalindromicPairLieTrotterGodunov`](@ref) is this construction at `p = 1`. + +As everywhere here, the order and the estimate cover the *splitting* error only and +presume the inner solvers resolve their subproblems accurately relative to it. +""" +struct AdjointPair{BaseType, AlgTupleType <: Tuple} <: AbstractOperatorSplittingAlgorithm + base::BaseType + inner_algs::AlgTupleType # aliases `base.inner_algs`, so the tree machinery works unchanged + + function AdjointPair(base::AbstractOperatorSplittingAlgorithm) + p = order(base) + isodd(p) || throw( + ArgumentError( + "AdjointPair needs a base scheme of odd order, got order $p. For even \ + orders a scheme and its adjoint share the same leading error term, so \ + averaging them raises no order and their difference is not an error \ + estimate." + ) + ) + return new{typeof(base), typeof(base.inner_algs)}(base, base.inner_algs) + end +end + +function Base.show(io::IO, alg::AdjointPair) + print(io, "AdjointPair (") + Base.show(io, alg.base) + return print(io, ")") +end + +coefficients(alg::AdjointPair) = coefficients(alg.base) +order(alg::AdjointPair) = order(alg.base) + 1 + +@inline SciMLBase.isadaptive(::AdjointPair) = true +# The estimate measures the *base* scheme's leading error, not the pair's. +alg_adaptive_order(alg::AdjointPair) = order(alg.base) + +struct AdjointPairCache{uType, uprevType, uforwardType, coeffType} <: AbstractOperatorSplittingCache + u::uType + uprev::uprevType + uforward::uforwardType # end state of the base sequence; reused as the residual buffer + coeffs::coeffType +end + +function init_cache( + f::GenericSplitFunction, alg::AdjointPair; + uprev::AbstractArray, u::AbstractVector, + ) + return AdjointPairCache(u, uprev, similar(u), coefficients(alg)) +end + +function _perform_step!( + parent, + children::Tuple, + cache::AdjointPairCache, + dt + ) + (; uforward, coeffs) = cache + + _table_stages!(parent, children, coeffs.a, dt) + parent.force_stepfail && return + uforward .= parent.u + + # Rewind to the initial state of the step; uprev is untouched while stepping. + parent.u .= parent.uprev + rollback_children!(parent) + + _table_stages_adjoint!( + parent, reverse(children), reverse(coeffs.a), dt, length(children) + ) + parent.force_stepfail && return + + # The average of the pair is the order p+1 solution ... + parent.u .= (parent.u .+ uforward) ./ 2 + if parent.controller_cache !== nothing + # ... and half the pair difference the local error of the base scheme. + (; abstol, reltol, internalnorm) = parent.opts + @. uforward = (parent.u - uforward) / + (abstol + max(abs(parent.u), abs(parent.uprev)) * reltol) + parent.EEst = internalnorm(uforward, parent.t + dt) + end + return +end + +# The adjoint reverses the whole flat sequence of flows, so both the stage order and +# the operator order within each stage are reversed. +@unroll function _table_stages_adjoint!(parent, rchildren, rstages::Tuple, dt, N) + @unroll for stage in rstages + _table_stage_adjoint!(parent, rchildren, stage, dt, N) + parent.force_stepfail && return + end +end + +@unroll function _table_stage_adjoint!(parent, rchildren::Tuple, stage, dt, N) + j = 0 + @unroll for child in rchildren + j += 1 + i = N + 1 - j + coefficient = stage[i] + if !iszero(coefficient) + _advance_child!(parent, child, i, coefficient * dt) + parent.force_stepfail && return + end + end +end diff --git a/src/utils.jl b/src/utils.jl index e7c059e..63df700 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -260,6 +260,86 @@ function _fix_dt_at_bounds!(integrator::AnySplitIntegrator) return nothing end +""" + reverse_direction!(integrator) + +Flip an already-initialized integrator's direction of integration in place. + +Negating `dt` alone does not work: `fix_dt_at_bounds!` re-signs it to `tdir` and would +replace it with `+dtmin` before `perform_step!` saw it. + +The `tstops`/`saveat`/`d_discontinuities` heaps store `tdir`-scaled times, so negating +every key re-expresses the same raw times under the new direction -- which inverts the +heap order, hence the rebuild. Times now *behind* are dropped, because `handle_tstop!` +errors on an unconsumed stop that `t` has passed. +""" +function reverse_direction!(integrator::DEIntegrator) + integrator.tdir = -integrator.tdir + integrator.dt = -integrator.dt + integrator.dtcache = -integrator.dtcache + integrator.dtpropose = -integrator.dtpropose + + opts = integrator.opts + opts.dtmax = -opts.dtmax + + threshold = integrator.tdir * integrator.t + _reverse_time_heap!(opts.tstops, threshold) + _reverse_time_heap!(opts.saveat, threshold) + _reverse_time_heap!(opts.d_discontinuities, threshold) + + return integrator +end + +function reverse_direction!(sub::SplitSubIntegrator) + sub.dt = -sub.dt + sub.dtcache = -sub.dtcache + sub.tdir = -sub.tdir + + # This heap stores *raw* times and carries the direction in its ordering, which is + # part of its type, so unlike a leaf's its keys cannot be re-signed. + _drop_times_behind!(sub.tstops, sub.tdir, sub.t) + + # `add_tstop!` propagates eagerly to every descendant, so a child left facing the + # old direction would reject the reversed node's next tstop as behind it. + _reverse_children!(sub.child_subintegrators) + + return sub +end + +@unroll function _reverse_children!(children::Tuple) + @unroll for child in children + reverse_direction!(child) + end +end + +function _drop_times_behind!(heap, tdir, t) + isempty(heap) && return heap + old = [pop!(heap)] + while !isempty(heap) + push!(old, pop!(heap)) + end + threshold = tdir * t + for key in old + tdir * key > threshold && push!(heap, key) + end + return heap +end + +function _reverse_time_heap!(heap, threshold) + isempty(heap) && return heap + # Drain first: the new keys are the negated old ones, so pushing them back into + # the heap being read from would corrupt the ordering mid-traversal. + old = [pop!(heap)] + while !isempty(heap) + push!(old, pop!(heap)) + end + for key in old + reversed = -key + reversed > threshold && push!(heap, reversed) + end + return heap +end + # Check time-step information consistency validate_time_point(integrator::AnySplitIntegrator) = validate_time_point(integrator, integrator.child_subintegrators) function validate_time_point(parent, child::SplitSubIntegrator) diff --git a/test/adjoint_pair.jl b/test/adjoint_pair.jl new file mode 100644 index 0000000..979a1fc --- /dev/null +++ b/test/adjoint_pair.jl @@ -0,0 +1,87 @@ +using OrdinaryDiffEqOperatorSplitting +using Test + +import DiffEqBase: DiffEqBase, ODEFunction +import SciMLBase +import SciMLBase: ReturnCode +using OrdinaryDiffEqTsit5 + +using OrdinaryDiffEqOperatorSplitting: order, coefficients + +M = [-1.0 0.5; 0.5 -2.0] +odeA(du, u, p, t) = (du[1] = -u[1]; du[2] = -2 * u[2]; nothing) +odeB(du, u, p, t) = (du[1] = 0.5 * u[2]; du[2] = 0.5 * u[1]; nothing) +dofs = [1, 2] +u0 = [1.0, 1.0] +fsplit = GenericSplitFunction((ODEFunction(odeA), ODEFunction(odeB)), (dofs, dofs)) +uexact = exp(M) * u0 + +function splitting_error(alg, dt) + prob = OperatorSplittingProblem(fsplit, copy(u0), (0.0, 1.0)) + integ = DiffEqBase.init( + prob, alg; dt = dt, adaptive = false, abstol = 1.0e-14, reltol = 1.0e-14 + ) + DiffEqBase.solve!(integ) + @assert integ.sol.retcode == ReturnCode.Success + return maximum(abs, integ.u .- uexact) +end + +@testset "Adjoint pairs" begin + @testset "the base scheme must be of odd order" begin + # Odd p is what makes the averaging work: the adjoint's leading error is + # (-1)^p C h^(p+1), so the signs oppose and cancel. For even p they are equal, + # the average stays order p, and the difference is no longer an error estimate. + @test AdjointPair(Ruth3((Tsit5(), Tsit5()))) isa AdjointPair + @test_throws ArgumentError AdjointPair(StrangMarchuk((Tsit5(), Tsit5()))) + end + + @testset "the pair raises the order by one" begin + alg = AdjointPair(Ruth3((Tsit5(), Tsit5()))) + @test order(alg) == 4 + # The estimate measures the *base* scheme's leading error, so that is the + # order the controller sees. + @test OrdinaryDiffEqOperatorSplitting.alg_adaptive_order(alg) == 3 + @test SciMLBase.isadaptive(alg) + end + + @testset "AdjointPair(LieTrotterGodunov) reproduces PalindromicPairLieTrotterGodunov" begin + # PPLTG is exactly this construction at p = 1, hand-written. Agreement to + # machine precision tests the whole generic path -- table execution, the + # flat-sequence adjoint, averaging, rollback -- against trusted code. + prob_pair = OperatorSplittingProblem(fsplit, copy(u0), (0.0, 1.0)) + pair = DiffEqBase.init( + prob_pair, AdjointPair(LieTrotterGodunov((Tsit5(), Tsit5()))); + dt = 0.1, adaptive = false, abstol = 1.0e-14, reltol = 1.0e-14 + ) + DiffEqBase.solve!(pair) + + prob_ppltg = OperatorSplittingProblem(fsplit, copy(u0), (0.0, 1.0)) + ppltg = DiffEqBase.init( + prob_ppltg, PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())); + dt = 0.1, adaptive = false, abstol = 1.0e-14, reltol = 1.0e-14 + ) + DiffEqBase.solve!(ppltg) + + @test pair.u == ppltg.u + end + + @testset "the pair is more accurate than its base scheme" begin + # The order itself is checked in test/convergence.jl. What matters here is that + # averaging buys accuracy rather than merely costing two sequences. + @test splitting_error(AdjointPair(Ruth3((Tsit5(), Tsit5()))), 0.1) < + splitting_error(Ruth3((Tsit5(), Tsit5())), 0.1) + end + + @testset "the pair difference drives adaptive stepping" begin + prob = OperatorSplittingProblem(fsplit, copy(u0), (0.0, 1.0)) + integ = DiffEqBase.init( + prob, AdjointPair(Ruth3((Tsit5(), Tsit5()))); + dt = 0.5, abstol = 1.0e-10, reltol = 1.0e-8 + ) + DiffEqBase.solve!(integ) + @test integ.sol.retcode == ReturnCode.Success + @test integ.t ≈ 1.0 + @test isfinite(integ.EEst) + @test maximum(abs, integ.u .- uexact) < 1.0e-6 + end +end diff --git a/test/callbacks.jl b/test/callbacks.jl index 3ef3fe9..2645e76 100644 --- a/test/callbacks.jl +++ b/test/callbacks.jl @@ -100,6 +100,7 @@ exact_crossing(level) = -10 * log(level / 3) ltg(), StrangMarchuk((Euler(), Euler())), PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())), + AdjointPair(Ruth3((Tsit5(), Tsit5()))), ) cb = DiscreteCallback( (u, t, integrator) -> isapprox(t, 0.5), diff --git a/test/coefficient_schemes.jl b/test/coefficient_schemes.jl new file mode 100644 index 0000000..e8a3053 --- /dev/null +++ b/test/coefficient_schemes.jl @@ -0,0 +1,72 @@ +using OrdinaryDiffEqOperatorSplitting +using Test + +import DiffEqBase: DiffEqBase, ODEFunction +import SciMLBase +import SciMLBase: ReturnCode +using OrdinaryDiffEqTsit5 + +using OrdinaryDiffEqOperatorSplitting: SplittingCoefficients, coefficients, order + +# Non-commuting pair with a known exact solution, as in test/backward.jl. +M = [-1.0 0.5; 0.5 -2.0] +odeA(du, u, p, t) = (du[1] = -u[1]; du[2] = -2 * u[2]; nothing) +odeB(du, u, p, t) = (du[1] = 0.5 * u[2]; du[2] = 0.5 * u[1]; nothing) +dofs = [1, 2] +u0 = [1.0, 1.0] +fsplit = GenericSplitFunction((ODEFunction(odeA), ODEFunction(odeB)), (dofs, dofs)) + +# Convergence orders live in test/convergence.jl, which runs every algorithm over two +# and three operators, nested, forward and backward. What is left here is what is +# specific to building a scheme out of a coefficient table. +@testset "Coefficient-table splitting schemes" begin + @testset "each operator's coefficients must sum to one" begin + # Ruth's third order table (Ruth, 1983): both columns sum to 1. + @test SplittingCoefficients((2 // 3, 7 // 24), (-2 // 3, 3 // 4), (1 // 1, -1 // 24)) isa + SplittingCoefficients + # Operator 1 sums to 0.9, so the scheme is not even consistent. + @test_throws ArgumentError SplittingCoefficients((0.5, 0.5), (0.4, 0.5)) + end + + @testset "Ruth3 exposes its table and order" begin + alg = Ruth3((Tsit5(), Tsit5())) + @test order(alg) == 3 + @test coefficients(alg) isa SplittingCoefficients + end + + @testset "tables stay compile-time constants" begin + # The stage count is a type parameter derived from `length`, and `@unroll` can + # only unroll the stage loop if that survives inference. Lie-Trotter's table is + # the one actually built at run time, from the length of `inner_algs`. + @test (@inferred coefficients(Ruth3((Tsit5(), Tsit5())))) isa + SplittingCoefficients{3, 2} + @test (@inferred coefficients(LieTrotterGodunov((Tsit5(), Tsit5())))) isa + SplittingCoefficients{1, 2} + end + + @testset "the table is re-applied faithfully on every step" begin + # A table ends on the last operator, so at the next step's start operator 1's + # buffer is stale and the `mark_next_sync_continuous` shortcut is invalid for + # it. That bug leaves the *first* step exact and corrupts every later one, so + # a single-step check cannot see it: compare several steps against the + # scheme's own composition of exact linear flows. + A = [-1.0 0.0; 0.0 -2.0] + B = [0.0 0.5; 0.5 0.0] + a = coefficients(Ruth3((Tsit5(), Tsit5()))).a + + h = 0.1 + reference = copy(u0) + for _ in 1:4, stage in a + reference = exp(B * (float(stage[2]) * h)) * (exp(A * (float(stage[1]) * h)) * reference) + end + + prob = OperatorSplittingProblem(fsplit, copy(u0), (0.0, 4h)) + integ = DiffEqBase.init( + prob, Ruth3((Tsit5(), Tsit5())); dt = h, abstol = 1.0e-14, reltol = 1.0e-14 + ) + DiffEqBase.solve!(integ) + + @test integ.iter == 4 + @test integ.u ≈ reference atol = 1.0e-12 + end +end diff --git a/test/convergence.jl b/test/convergence.jl index bb836d8..e59a707 100644 --- a/test/convergence.jl +++ b/test/convergence.jl @@ -44,21 +44,28 @@ function convergence_rates(f, alg, tspan, ustart, utarget; dts = (0.1, 0.05, 0.0 end @testset "Convergence order" begin - for (AlgType, expected_order) in ( - (LieTrotterGodunov, 1), - (StrangMarchuk, 2), - (PalindromicPairLieTrotterGodunov, 2), - ) - cases = ( - ("two operators", f_two, AlgType((Tsit5(), Tsit5()))), - ("three operators", f_three, AlgType((Tsit5(), Tsit5(), Tsit5()))), - ("nested", f_nested, AlgType((Tsit5(), AlgType((Tsit5(), Tsit5()))))), + # `build` takes the tuple of inner algorithms, so schemes that wrap a base rather + # than taking `inner_algs` directly (AdjointPair) fit the same loop. Ruth3 is a + # two-operator table, so it sits out the three-operator case. + for (name, build, expected_order, supports_three) in ( + ("LieTrotterGodunov", LieTrotterGodunov, 1, true), + ("StrangMarchuk", StrangMarchuk, 2, true), + ("PalindromicPairLieTrotterGodunov", PalindromicPairLieTrotterGodunov, 2, true), + ("Ruth3", Ruth3, 3, false), + ("AdjointPair(Ruth3)", inner -> AdjointPair(Ruth3(inner)), 4, false), ) + two = ("two operators", f_two, build((Tsit5(), Tsit5()))) + nested = ("nested", f_nested, build((Tsit5(), build((Tsit5(), Tsit5()))))) + cases = if supports_three + (two, ("three operators", f_three, build((Tsit5(), Tsit5(), Tsit5()))), nested) + else + (two, nested) + end directions = ( ("forward", (0.0, 1.0), u0, uT), ("backward", (1.0, 0.0), uT, u0), ) - @testset "$(nameof(AlgType)) | $case | $dir (order $expected_order)" for + @testset "$name | $case | $dir (order $expected_order)" for (case, f, alg) in cases, (dir, tspan, ustart, utarget) in directions rates = convergence_rates(f, alg, tspan, ustart, utarget) @@ -68,3 +75,8 @@ end end end end + +@testset "Ruth3 rejects operator counts its table cannot address" begin + @test_throws ArgumentError Ruth3((Tsit5(), Tsit5(), Tsit5())) + @test_throws ArgumentError Ruth3((Tsit5(),)) +end diff --git a/test/saving.jl b/test/saving.jl index bc545f3..2220e1c 100644 --- a/test/saving.jl +++ b/test/saving.jl @@ -33,6 +33,8 @@ algs() = ( LieTrotterGodunov((Tsit5(), Tsit5())), StrangMarchuk((Tsit5(), Tsit5())), PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())), + Ruth3((Tsit5(), Tsit5())), + AdjointPair(Ruth3((Tsit5(), Tsit5()))), ) function exact(t) From 0594c49402a8d378d25f361c674aaa56fa3a61bd Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 4 Aug 2026 16:44:46 -0400 Subject: [PATCH 2/6] Add Yoshida4 and the missing bibliography entries Yoshida's triple jump, derived in this package's stage convention rather than transcribed: composing three Strang steps of lengths w1*h, w0*h and w1*h and merging the adjacent operator-1 flows gives four stages and eight flow evaluations rather than nine. The merge zeroes the last stage's second coefficient, so this is also the first table to exercise the zero-coefficient skip in the generic step. Measured order 4.02, 4.00, 4.00 over dt 0.2 to 0.025. w0 is negative and large: the second operator's cumulative time reaches 1.35*h before returning through -0.35*h, so this leans hard on backward substeps. Being of even order, Yoshida4 cannot serve as an AdjointPair base, which is now asserted alongside the StrangMarchuk case. Ruth3's two-operator restriction is factored into a shared helper, since Yoshida4's table is also AB only. Also adds the three bibliography entries the docstrings had been citing without defining -- Rut:1983:cim, Yos:1990:cho and AuzHofKetKoc:2017:psm -- which would have failed the docs build. Co-Authored-By: Claude Opus 5 --- docs/src/api-reference/index.md | 1 + docs/src/assets/references.bib | 36 ++++++++++++ src/OrdinaryDiffEqOperatorSplitting.jl | 2 +- src/solver.jl | 77 +++++++++++++++++++++++--- test/adjoint_pair.jl | 1 + test/coefficient_schemes.jl | 12 ++++ test/convergence.jl | 4 +- 7 files changed, 124 insertions(+), 9 deletions(-) diff --git a/docs/src/api-reference/index.md b/docs/src/api-reference/index.md index 5dfdc4d..4662502 100644 --- a/docs/src/api-reference/index.md +++ b/docs/src/api-reference/index.md @@ -19,6 +19,7 @@ LieTrotterGodunov StrangMarchuk PalindromicPairLieTrotterGodunov Ruth3 +Yoshida4 AdjointPair ``` diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index 3cb2732..de7763c 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -52,3 +52,39 @@ @incollection{Mar:1971:tsm year={1971}, publisher={Academic Press} } + +@article{Rut:1983:cim, + title={A canonical integration technique}, + author={Ruth, Ronald D.}, + journal={IEEE Transactions on Nuclear Science}, + volume={NS-30}, + number={4}, + pages={2669--2671}, + year={1983}, + doi={10.1109/TNS.1983.4332919}, + publisher={IEEE} +} + +@article{Yos:1990:cho, + title={Construction of higher order symplectic integrators}, + author={Yoshida, Haruo}, + journal={Physics Letters A}, + volume={150}, + number={5--7}, + pages={262--268}, + year={1990}, + doi={10.1016/0375-9601(90)90092-3}, + publisher={Elsevier} +} + +@article{AuzHofKetKoc:2017:psm, + title={Practical splitting methods for the adaptive integration of nonlinear evolution equations. {P}art {I}: {C}onstruction of optimized schemes and pairs of schemes}, + author={Auzinger, Winfried and Hofst{\"a}tter, Harald and Ketcheson, David and Koch, Othmar}, + journal={BIT Numerical Mathematics}, + volume={57}, + number={1}, + pages={55--74}, + year={2017}, + doi={10.1007/s10543-016-0626-9}, + publisher={Springer} +} diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index 4af8b30..60e1cad 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -154,7 +154,7 @@ include("solver.jl") include("utils.jl") export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, StrangMarchuk, - PalindromicPairLieTrotterGodunov, Ruth3, AdjointPair + PalindromicPairLieTrotterGodunov, Ruth3, Yoshida4, AdjointPair export SplitNode, TreeOption include("precompilation.jl") diff --git a/src/solver.jl b/src/solver.jl index 6389a25..137bb80 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -290,6 +290,17 @@ end # --------------------------------------------------------------------------- # Coefficient-table splitting schemes # --------------------------------------------------------------------------- +function _require_two_operators(scheme, inner_algs) + n = length(inner_algs) + n == 2 || throw( + ArgumentError( + "$scheme is a two-operator (AB) table but got $n operators. Group the \ + operators into a nested GenericSplitFunction to use it with more." + ) + ) + return nothing +end + """ SplittingCoefficients(stages::NTuple{N, T}...) @@ -361,13 +372,7 @@ struct Ruth3{AlgTupleType <: Tuple} <: AbstractOperatorSplittingAlgorithm inner_algs::AlgTupleType function Ruth3(inner_algs::Tuple) - n = length(inner_algs) - n == 2 || throw( - ArgumentError( - "Ruth3 is a two-operator (AB) table but got $n operators. Group the \ - operators into a nested GenericSplitFunction to use it with more." - ) - ) + _require_two_operators("Ruth3", inner_algs) return new{typeof(inner_algs)}(inner_algs) end end @@ -389,6 +394,64 @@ const RUTH3_COEFFICIENTS = SplittingCoefficients( coefficients(::Ruth3) = RUTH3_COEFFICIENTS order(::Ruth3) = 3 +""" + Yoshida4 <: AbstractOperatorSplittingAlgorithm + +Fourth-order splitting scheme of [Yos:1990:cho](@cite), the "triple jump". + +Built by composing three Strang steps of lengths ``w_1 h``, ``w_0 h`` and ``w_1 h`` +with ``w_1 = 1/(2 - 2^{1/3})`` and ``w_0 = -2^{1/3} w_1``, then merging the adjacent +flows the composition leaves next to each other. That merging is what makes it eight +flow evaluations rather than nine, and it leaves the last stage's second coefficient +zero. + +``w_0`` is negative, so a substantial part of each step runs backward in time -- the +second operator's cumulative time reaches ``1.35\\,h`` before returning through +``-0.35\\,h``. + +As for every splitting scheme here the order statement covers the *splitting* error +only, and presumes the inner solvers resolve their subproblems accurately relative +to it. +""" +struct Yoshida4{AlgTupleType <: Tuple} <: AbstractOperatorSplittingAlgorithm + inner_algs::AlgTupleType + + function Yoshida4(inner_algs::Tuple) + _require_two_operators("Yoshida4", inner_algs) + return new{typeof(inner_algs)}(inner_algs) + end +end + +function Base.show(io::IO, alg::Yoshida4) + print(io, "Yoshida4 (") + for inner_alg in alg.inner_algs[1:(end - 1)] + Base.show(io, inner_alg) + print(io, " -> ") + end + length(alg.inner_algs) > 0 && Base.show(io, alg.inner_algs[end]) + return print(io, ")") +end + +const YOSHIDA4_W1 = 1 / (2 - cbrt(2)) +const YOSHIDA4_W0 = -cbrt(2) * YOSHIDA4_W1 + +const YOSHIDA4_COEFFICIENTS = SplittingCoefficients( + (YOSHIDA4_W1 / 2, YOSHIDA4_W1), + ((YOSHIDA4_W1 + YOSHIDA4_W0) / 2, YOSHIDA4_W0), + ((YOSHIDA4_W1 + YOSHIDA4_W0) / 2, YOSHIDA4_W1), + (YOSHIDA4_W1 / 2, 0.0), +) + +coefficients(::Yoshida4) = YOSHIDA4_COEFFICIENTS +order(::Yoshida4) = 4 + +function init_cache( + f::GenericSplitFunction, alg::Yoshida4; + uprev::AbstractArray, u::AbstractVector, + ) + return SplittingCoefficientsCache(u, uprev, coefficients(alg)) +end + order(::StrangMarchuk) = 2 order(::PalindromicPairLieTrotterGodunov) = 2 diff --git a/test/adjoint_pair.jl b/test/adjoint_pair.jl index 979a1fc..53bab32 100644 --- a/test/adjoint_pair.jl +++ b/test/adjoint_pair.jl @@ -33,6 +33,7 @@ end # the average stays order p, and the difference is no longer an error estimate. @test AdjointPair(Ruth3((Tsit5(), Tsit5()))) isa AdjointPair @test_throws ArgumentError AdjointPair(StrangMarchuk((Tsit5(), Tsit5()))) + @test_throws ArgumentError AdjointPair(Yoshida4((Tsit5(), Tsit5()))) end @testset "the pair raises the order by one" begin diff --git a/test/coefficient_schemes.jl b/test/coefficient_schemes.jl index e8a3053..4433980 100644 --- a/test/coefficient_schemes.jl +++ b/test/coefficient_schemes.jl @@ -34,6 +34,18 @@ fsplit = GenericSplitFunction((ODEFunction(odeA), ODEFunction(odeB)), (dofs, dof @test coefficients(alg) isa SplittingCoefficients end + @testset "Yoshida4's merged table carries a zero coefficient" begin + # Merging the adjacent flows the triple jump leaves next to each other is what + # brings it down to eight evaluations, and it zeroes the last stage's second + # coefficient. Its convergence in test/convergence.jl is therefore also the + # coverage for the zero-coefficient skip path. + table = coefficients(Yoshida4((Tsit5(), Tsit5()))).a + @test order(Yoshida4((Tsit5(), Tsit5()))) == 4 + @test length(table) == 4 + @test count(iszero, Iterators.flatten(table)) == 1 + @test iszero(table[end][end]) + end + @testset "tables stay compile-time constants" begin # The stage count is a type parameter derived from `length`, and `@unroll` can # only unroll the stage loop if that survives inference. Lie-Trotter's table is diff --git a/test/convergence.jl b/test/convergence.jl index e59a707..2fab596 100644 --- a/test/convergence.jl +++ b/test/convergence.jl @@ -52,6 +52,7 @@ end ("StrangMarchuk", StrangMarchuk, 2, true), ("PalindromicPairLieTrotterGodunov", PalindromicPairLieTrotterGodunov, 2, true), ("Ruth3", Ruth3, 3, false), + ("Yoshida4", Yoshida4, 4, false), ("AdjointPair(Ruth3)", inner -> AdjointPair(Ruth3(inner)), 4, false), ) two = ("two operators", f_two, build((Tsit5(), Tsit5()))) @@ -76,7 +77,8 @@ end end end -@testset "Ruth3 rejects operator counts its table cannot address" begin +@testset "two-operator tables reject other operator counts" begin @test_throws ArgumentError Ruth3((Tsit5(), Tsit5(), Tsit5())) @test_throws ArgumentError Ruth3((Tsit5(),)) + @test_throws ArgumentError Yoshida4((Tsit5(), Tsit5(), Tsit5())) end From adcbcd57bb19b67d407fc0d414f23442d14977b7 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Tue, 4 Aug 2026 16:49:25 -0400 Subject: [PATCH 3/6] Allow "Ket" in the spellchecker for the Ketcheson citation key The citation key AuzHofKetKoc:2017:psm splits into Auz/Hof/Ket/Koc, and typos reads "Ket" as a misspelling of "Kept" or "Key". Added alongside the existing "Tro" exception for the Trotter key. Co-Authored-By: Claude Opus 5 --- _typos.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/_typos.toml b/_typos.toml index 625e683..e3bfffd 100644 --- a/_typos.toml +++ b/_typos.toml @@ -3,5 +3,7 @@ Strang = "Strang" # Citation key for Trotter (Tro:1959:psg) used in docstrings Tro = "Tro" +# David Ketcheson, in the citation key AuzHofKetKoc:2017:psm +Ket = "Ket" # Splitting variable name BA = "BA" From 5615f479d6ca4740096d4f9e682942c344a0da34 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Sun, 9 Aug 2026 00:11:22 -0400 Subject: [PATCH 4/6] respond to review --- src/solver.jl | 639 -------------------------------------------------- 1 file changed, 639 deletions(-) delete mode 100644 src/solver.jl diff --git a/src/solver.jl b/src/solver.jl deleted file mode 100644 index 137bb80..0000000 --- a/src/solver.jl +++ /dev/null @@ -1,639 +0,0 @@ -# --------------------------------------------------------------------------- -# Lie-Trotter-Godunov operator splitting -# --------------------------------------------------------------------------- -""" - LieTrotterGodunov <: AbstractOperatorSplittingAlgorithm - -First-order sequential operator splitting algorithm attributed to -[Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite). -""" -struct LieTrotterGodunov{AlgTupleType} <: AbstractOperatorSplittingAlgorithm - inner_algs::AlgTupleType # Tuple of timesteppers for inner problems -end - -function Base.show(io::IO, alg::LieTrotterGodunov) - print(io, "LTG (") - for inner_alg in alg.inner_algs[1:(end - 1)] - Base.show(io, inner_alg) - print(io, " -> ") - end - length(alg.inner_algs) > 0 && Base.show(io, alg.inner_algs[end]) - return print(io, ")") -end - -struct LieTrotterGodunovCache{uType, uprevType} <: AbstractOperatorSplittingCache - u::uType - uprev::uprevType -end - -function init_cache( - f::GenericSplitFunction, alg::LieTrotterGodunov; - uprev::AbstractArray, u::AbstractVector, - ) - return LieTrotterGodunovCache(u, uprev) -end - -@unroll function _perform_step!( - parent, - children::Tuple, - cache::LieTrotterGodunovCache, - dt - ) - i = 0 - @unroll for child in children - i += 1 - - idxs = parent.child_solution_indices[i] - sync = parent.child_synchronizers[i] - - @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync) - @timeit_debug "time solve" advance_solution_by!(parent, child, dt) - if child_failed(child) - parent.force_stepfail = true - return - end - - @timeit_debug "sync <-" backward_sync_subintegrator!(parent, child, idxs, sync) - end -end - -# --------------------------------------------------------------------------- -# Strang-Marchuk operator splitting -# --------------------------------------------------------------------------- -""" - StrangMarchuk <: AbstractOperatorSplittingAlgorithm - -Second-order symmetric (palindromic) operator splitting algorithm attributed to -[Str:1968:ccd,Mar:1971:tsm](@cite). - -For ``N`` operators the scheme performs - -``A_1(\\Delta t/2) \\to \\cdots \\to A_{N-1}(\\Delta t/2) \\to A_N(\\Delta t) \\to A_{N-1}(\\Delta t/2) \\to \\cdots \\to A_1(\\Delta t/2)`` - -achieving second-order accuracy through symmetry. -""" -struct StrangMarchuk{AlgTupleType} <: AbstractOperatorSplittingAlgorithm - inner_algs::AlgTupleType # Tuple of timesteppers for inner problems -end - -function Base.show(io::IO, alg::StrangMarchuk) - print(io, "SM (") - for inner_alg in alg.inner_algs[1:(end - 1)] - Base.show(io, inner_alg) - print(io, " -> ") - end - length(alg.inner_algs) > 0 && Base.show(io, alg.inner_algs[end]) - return print(io, ")") -end - -struct StrangMarchukCache{uType, uprevType} <: AbstractOperatorSplittingCache - u::uType - uprev::uprevType -end - -function init_cache( - f::GenericSplitFunction, alg::StrangMarchuk; - uprev::AbstractArray, u::AbstractVector, - ) - return StrangMarchukCache(u, uprev) -end - -# Forward pass: A₁(dt/2) → … → Aₙ₋₁(dt/2) → Aₙ(dt) -@unroll function _sm_forward_pass!(parent, children::Tuple, half_dt, dt) - N = length(children) - i = 0 - @unroll for child in children - i += 1 - step_dt = i < N ? half_dt : dt - - idxs = parent.child_solution_indices[i] - sync = parent.child_synchronizers[i] - - @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync) - @timeit_debug "time solve" advance_solution_by!(parent, child, step_dt) - if child_failed(child) - parent.force_stepfail = true - return - end - - @timeit_debug "sync <-" backward_sync_subintegrator!(parent, child, idxs, sync) - end -end - -# Reverse pass: Aₙ₋₁(dt/2) → … → A₁(dt/2) -@unroll function _sm_reverse_pass!(parent, rev_front::Tuple, half_dt, N) - j = 0 - @unroll for child in rev_front - j += 1 - i = N - j - - idxs = parent.child_solution_indices[i] - sync = parent.child_synchronizers[i] - - @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync) - @timeit_debug "time solve" advance_solution_by!(parent, child, half_dt) - if child_failed(child) - parent.force_stepfail = true - return - end - - @timeit_debug "sync <-" backward_sync_subintegrator!(parent, child, idxs, sync) - end -end - -function _perform_step!( - parent, - children::Tuple, - cache::StrangMarchukCache, - dt - ) - half_dt = dt / 2 - - # Skip sync of for first solve, because it is already in sync - mark_next_sync_continuous(parent) - - _sm_forward_pass!(parent, children, half_dt, dt) - parent.force_stepfail && return - - _sm_reverse_pass!(parent, reverse(children[1:(end - 1)]), half_dt, length(children)) - parent.force_stepfail && return - - return -end - -# --------------------------------------------------------------------------- -# Palindromic pair of Lie-Trotter-Godunov splittings -# --------------------------------------------------------------------------- -""" - PalindromicPairLieTrotterGodunov <: AbstractOperatorSplittingAlgorithm - -Second-order sequential operator splitting algorithm. - -One step solves the palindromic pair of [`LieTrotterGodunov`](@ref) sequences - -``A_1(\\Delta t) \\to \\cdots \\to A_N(\\Delta t)`` and -``A_N(\\Delta t) \\to \\cdots \\to A_1(\\Delta t)`` - -from the same initial value. The leading splitting error of a Lie-Trotter sequence -is ``\\frac{\\Delta t^2}{2}\\sum_{i ") - end - length(alg.inner_algs) > 0 && Base.show(io, alg.inner_algs[end]) - return print(io, ")") -end - -@inline SciMLBase.isadaptive(::PalindromicPairLieTrotterGodunov) = true -# The pair difference estimates the O(dt²) leading error term of a first order -# sequence, so the controller sees a first order error estimator. -alg_adaptive_order(::PalindromicPairLieTrotterGodunov) = 1 - -struct PalindromicPairLieTrotterGodunovCache{uType, uprevType, uforwardType} <: AbstractOperatorSplittingCache - u::uType - uprev::uprevType - uforward::uforwardType # end state of the A₁ → A₂ sequence; reused as the residual buffer -end - -function init_cache( - f::GenericSplitFunction, alg::PalindromicPairLieTrotterGodunov; - uprev::AbstractArray, u::AbstractVector, - ) - return PalindromicPairLieTrotterGodunovCache(u, uprev, similar(u)) -end - -function _advance_child!(parent, child, i, dt) - idxs = parent.child_solution_indices[i] - sync = parent.child_synchronizers[i] - - @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync) - @timeit_debug "time solve" advance_solution_by!(parent, child, dt) - if child_failed(child) - parent.force_stepfail = true - return - end - - @timeit_debug "sync <-" backward_sync_subintegrator!(parent, child, idxs, sync) - return -end - -# Forward sequence: A₁(dt) → … → A_N(dt) -@unroll function _ppltg_forward_pass!(parent, children::Tuple, dt) - i = 0 - @unroll for child in children - i += 1 - _advance_child!(parent, child, i, dt) - parent.force_stepfail && return - end -end - -# Reverse sequence: A_N(dt) → … → A₁(dt) -@unroll function _ppltg_reverse_pass!(parent, rchildren::Tuple, dt, N) - j = 0 - @unroll for child in rchildren - j += 1 - _advance_child!(parent, child, N + 1 - j, dt) - parent.force_stepfail && return - end -end - -function _perform_step!( - parent, - children::Tuple, - cache::PalindromicPairLieTrotterGodunovCache, - dt - ) - (; uforward) = cache - - _ppltg_forward_pass!(parent, children, dt) - parent.force_stepfail && return - uforward .= parent.u - - # Rewind to the initial state of the step; uprev is untouched while stepping. - parent.u .= parent.uprev - rollback_children!(parent) - - _ppltg_reverse_pass!(parent, reverse(children), dt, length(children)) - parent.force_stepfail && return - - # The average of the pair is the second order solution ... - parent.u .= (parent.u .+ uforward) ./ 2 - if parent.controller_cache !== nothing - # ... and half the pair difference the local error of a single sequence. - (; abstol, reltol, internalnorm) = parent.opts - @. uforward = (parent.u - uforward) / - (abstol + max(abs(parent.u), abs(parent.uprev)) * reltol) - OrdinaryDiffEqCore.set_EEst!(parent, internalnorm(uforward, parent.t + dt)) - end - return -end - -# --------------------------------------------------------------------------- -# Coefficient-table splitting schemes -# --------------------------------------------------------------------------- -function _require_two_operators(scheme, inner_algs) - n = length(inner_algs) - n == 2 || throw( - ArgumentError( - "$scheme is a two-operator (AB) table but got $n operators. Group the \ - operators into a nested GenericSplitFunction to use it with more." - ) - ) - return nothing -end - -""" - SplittingCoefficients(stages::NTuple{N, T}...) - -Coefficients of an `S`-stage splitting scheme over `N` operators, one tuple per stage: -stage `j` advances operator `i` by `stages[j][i] * dt`. - -This is the generalization to `N` operators of the two-operator (`AB`) and -three-operator (`ABC`) coefficient tables of -[AuzHofKetKoc:2017:psm](@cite); their tables are the `N = 2` and `N = 3` cases. - -Each operator's coefficients must sum to one, the consistency condition, and that is -checked here. The remaining order conditions are not, so a table that constructs -successfully can still fail to attain the order it claims. -""" -struct SplittingCoefficients{S, N, T} - a::NTuple{S, NTuple{N, T}} - - # Do not simplify this to `NTuple{S, NTuple{N, T}}`: the empty tuple matches it for - # any element type, leaving parameters unbound for `S == 0` (and `T` unbound for - # `N == 0`). A leading element plus a counted `Vararg` rules both out. - function SplittingCoefficients( - stage1::Tuple{T, Vararg{T, K}}, - rest::Tuple{T, Vararg{T, K}}... - ) where {T, K} - a = (stage1, rest...) - N = K + 1 - S = length(a) - for i in 1:N - total = sum(a[j][i] for j in 1:S) - total ≈ one(T) || throw( - ArgumentError( - "operator $i's coefficients sum to $total rather than 1, so the \ - scheme is not consistent." - ) - ) - end - return new{S, N, T}(a) - end -end - -""" - coefficients(alg) - -The [`SplittingCoefficients`](@ref) table of a coefficient-driven splitting algorithm. -""" -function coefficients end - -""" - order(alg) - -Order of consistency of a splitting algorithm, counting the splitting error only. -""" -function order end - -""" - Ruth3 <: AbstractOperatorSplittingAlgorithm - -Third-order splitting scheme of [Rut:1983:cim](@cite), in three stages. - -Its coefficients are exactly rational, and -- as is unavoidable for any real -splitting scheme of order three or above -- some of them are negative, so parts of -the step run backward in time. - -As for every splitting scheme here the order statement covers the *splitting* error -only, and presumes the inner solvers resolve their subproblems accurately relative -to it. -""" -struct Ruth3{AlgTupleType <: Tuple} <: AbstractOperatorSplittingAlgorithm - inner_algs::AlgTupleType - - function Ruth3(inner_algs::Tuple) - _require_two_operators("Ruth3", inner_algs) - return new{typeof(inner_algs)}(inner_algs) - end -end - -function Base.show(io::IO, alg::Ruth3) - print(io, "Ruth3 (") - for inner_alg in alg.inner_algs[1:(end - 1)] - Base.show(io, inner_alg) - print(io, " -> ") - end - length(alg.inner_algs) > 0 && Base.show(io, alg.inner_algs[end]) - return print(io, ")") -end - -const RUTH3_COEFFICIENTS = SplittingCoefficients( - (7 // 24, 2 // 3), (3 // 4, -2 // 3), (-1 // 24, 1 // 1) -) - -coefficients(::Ruth3) = RUTH3_COEFFICIENTS -order(::Ruth3) = 3 - -""" - Yoshida4 <: AbstractOperatorSplittingAlgorithm - -Fourth-order splitting scheme of [Yos:1990:cho](@cite), the "triple jump". - -Built by composing three Strang steps of lengths ``w_1 h``, ``w_0 h`` and ``w_1 h`` -with ``w_1 = 1/(2 - 2^{1/3})`` and ``w_0 = -2^{1/3} w_1``, then merging the adjacent -flows the composition leaves next to each other. That merging is what makes it eight -flow evaluations rather than nine, and it leaves the last stage's second coefficient -zero. - -``w_0`` is negative, so a substantial part of each step runs backward in time -- the -second operator's cumulative time reaches ``1.35\\,h`` before returning through -``-0.35\\,h``. - -As for every splitting scheme here the order statement covers the *splitting* error -only, and presumes the inner solvers resolve their subproblems accurately relative -to it. -""" -struct Yoshida4{AlgTupleType <: Tuple} <: AbstractOperatorSplittingAlgorithm - inner_algs::AlgTupleType - - function Yoshida4(inner_algs::Tuple) - _require_two_operators("Yoshida4", inner_algs) - return new{typeof(inner_algs)}(inner_algs) - end -end - -function Base.show(io::IO, alg::Yoshida4) - print(io, "Yoshida4 (") - for inner_alg in alg.inner_algs[1:(end - 1)] - Base.show(io, inner_alg) - print(io, " -> ") - end - length(alg.inner_algs) > 0 && Base.show(io, alg.inner_algs[end]) - return print(io, ")") -end - -const YOSHIDA4_W1 = 1 / (2 - cbrt(2)) -const YOSHIDA4_W0 = -cbrt(2) * YOSHIDA4_W1 - -const YOSHIDA4_COEFFICIENTS = SplittingCoefficients( - (YOSHIDA4_W1 / 2, YOSHIDA4_W1), - ((YOSHIDA4_W1 + YOSHIDA4_W0) / 2, YOSHIDA4_W0), - ((YOSHIDA4_W1 + YOSHIDA4_W0) / 2, YOSHIDA4_W1), - (YOSHIDA4_W1 / 2, 0.0), -) - -coefficients(::Yoshida4) = YOSHIDA4_COEFFICIENTS -order(::Yoshida4) = 4 - -function init_cache( - f::GenericSplitFunction, alg::Yoshida4; - uprev::AbstractArray, u::AbstractVector, - ) - return SplittingCoefficientsCache(u, uprev, coefficients(alg)) -end - -order(::StrangMarchuk) = 2 -order(::PalindromicPairLieTrotterGodunov) = 2 - -# Lie-Trotter keeps its hand-written step; the table exists only so it can serve as an -# `AdjointPair` base, which is what makes `AdjointPair(LieTrotterGodunov(...))` and -# `PalindromicPairLieTrotterGodunov` the same scheme. -order(::LieTrotterGodunov) = 1 -coefficients(alg::LieTrotterGodunov) = - SplittingCoefficients(ntuple(_ -> 1 // 1, length(alg.inner_algs))) - -struct SplittingCoefficientsCache{uType, uprevType, coeffType} <: AbstractOperatorSplittingCache - u::uType - uprev::uprevType - coeffs::coeffType -end - -function init_cache( - f::GenericSplitFunction, alg::Ruth3; - uprev::AbstractArray, u::AbstractVector, - ) - return SplittingCoefficientsCache(u, uprev, coefficients(alg)) -end - -function _perform_step!( - parent, - children::Tuple, - cache::SplittingCoefficientsCache, - dt - ) - # Deliberately no `mark_next_sync_continuous`: that shortcut needs the previous step - # to have left `parent.u` equal to the buffer of the child solved first, which holds - # for StrangMarchuk only because its reverse pass ends on operator 1. A general table - # ends on operator N, so operator 1's buffer is stale and skipping its forward sync - # resumes from stale state -- the first step stays exact and every later one is wrong. - _table_stages!(parent, children, cache.coeffs.a, dt) - return -end - -@unroll function _table_stages!(parent, children, stages::Tuple, dt) - @unroll for stage in stages - _table_stage!(parent, children, stage, dt) - parent.force_stepfail && return - end -end - -@unroll function _table_stage!(parent, children::Tuple, stage, dt) - i = 0 - @unroll for child in children - i += 1 - coefficient = stage[i] - # A zero coefficient is the identity flow, so it is skipped sync and all. - if !iszero(coefficient) - _advance_child!(parent, child, i, coefficient * dt) - parent.force_stepfail && return - end - end -end - -# --------------------------------------------------------------------------- -# Adjoint pairs -# --------------------------------------------------------------------------- -""" - AdjointPair(base) <: AbstractOperatorSplittingAlgorithm - -Adaptive splitting scheme of order `p+1` built from a `base` scheme of odd order `p`, -following [AuzHofKetKoc:2017:psm](@cite), eq. (3.2). - -One step runs `base` and its adjoint ``\\mathcal{S}^*`` from the same initial value. -Their leading error terms are ``C h^{p+1}`` and ``(-1)^p C h^{p+1}``, so for odd `p` -the signs oppose: the average is a solution of order `p+1`, and half the difference is -an asymptotically correct estimate of the base scheme's local error, which drives the -step size controller. A step therefore costs twice the base scheme. - -The base scheme's order must be **odd**. For even `p` the two leading terms are -*equal* rather than opposite, so averaging cancels nothing and the difference stops -being an error estimate. - -``\\mathcal{S}^*(h, u) = \\mathcal{S}^{-1}(-h, u)`` is the base scheme's entire flat -sequence of flows reversed, every coefficient keeping its sign and its operator, so it -reuses the same table and needs no extra coefficients. - -[`PalindromicPairLieTrotterGodunov`](@ref) is this construction at `p = 1`. - -As everywhere here, the order and the estimate cover the *splitting* error only and -presume the inner solvers resolve their subproblems accurately relative to it. -""" -struct AdjointPair{BaseType, AlgTupleType <: Tuple} <: AbstractOperatorSplittingAlgorithm - base::BaseType - inner_algs::AlgTupleType # aliases `base.inner_algs`, so the tree machinery works unchanged - - function AdjointPair(base::AbstractOperatorSplittingAlgorithm) - p = order(base) - isodd(p) || throw( - ArgumentError( - "AdjointPair needs a base scheme of odd order, got order $p. For even \ - orders a scheme and its adjoint share the same leading error term, so \ - averaging them raises no order and their difference is not an error \ - estimate." - ) - ) - return new{typeof(base), typeof(base.inner_algs)}(base, base.inner_algs) - end -end - -function Base.show(io::IO, alg::AdjointPair) - print(io, "AdjointPair (") - Base.show(io, alg.base) - return print(io, ")") -end - -coefficients(alg::AdjointPair) = coefficients(alg.base) -order(alg::AdjointPair) = order(alg.base) + 1 - -@inline SciMLBase.isadaptive(::AdjointPair) = true -# The estimate measures the *base* scheme's leading error, not the pair's. -alg_adaptive_order(alg::AdjointPair) = order(alg.base) - -struct AdjointPairCache{uType, uprevType, uforwardType, coeffType} <: AbstractOperatorSplittingCache - u::uType - uprev::uprevType - uforward::uforwardType # end state of the base sequence; reused as the residual buffer - coeffs::coeffType -end - -function init_cache( - f::GenericSplitFunction, alg::AdjointPair; - uprev::AbstractArray, u::AbstractVector, - ) - return AdjointPairCache(u, uprev, similar(u), coefficients(alg)) -end - -function _perform_step!( - parent, - children::Tuple, - cache::AdjointPairCache, - dt - ) - (; uforward, coeffs) = cache - - _table_stages!(parent, children, coeffs.a, dt) - parent.force_stepfail && return - uforward .= parent.u - - # Rewind to the initial state of the step; uprev is untouched while stepping. - parent.u .= parent.uprev - rollback_children!(parent) - - _table_stages_adjoint!( - parent, reverse(children), reverse(coeffs.a), dt, length(children) - ) - parent.force_stepfail && return - - # The average of the pair is the order p+1 solution ... - parent.u .= (parent.u .+ uforward) ./ 2 - if parent.controller_cache !== nothing - # ... and half the pair difference the local error of the base scheme. - (; abstol, reltol, internalnorm) = parent.opts - @. uforward = (parent.u - uforward) / - (abstol + max(abs(parent.u), abs(parent.uprev)) * reltol) - parent.EEst = internalnorm(uforward, parent.t + dt) - end - return -end - -# The adjoint reverses the whole flat sequence of flows, so both the stage order and -# the operator order within each stage are reversed. -@unroll function _table_stages_adjoint!(parent, rchildren, rstages::Tuple, dt, N) - @unroll for stage in rstages - _table_stage_adjoint!(parent, rchildren, stage, dt, N) - parent.force_stepfail && return - end -end - -@unroll function _table_stage_adjoint!(parent, rchildren::Tuple, stage, dt, N) - j = 0 - @unroll for child in rchildren - j += 1 - i = N + 1 - j - coefficient = stage[i] - if !iszero(coefficient) - _advance_child!(parent, child, i, coefficient * dt) - parent.force_stepfail && return - end - end -end From aa18c94192cf4bb95c624c916a91c25b0600f2ee Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Sun, 9 Aug 2026 15:13:05 -0400 Subject: [PATCH 5/6] actually commit solvers --- _typos.toml | 3 + docs/src/assets/references.bib | 24 +++ docs/src/devdocs/index.md | 85 ++++++++- docs/src/topics/adaptivity.md | 12 +- docs/src/topics/time-integration.md | 108 ++++++++++++ src/OrdinaryDiffEqOperatorSplitting.jl | 7 +- src/integrator.jl | 52 +++--- src/solvers/adjoint_pair.jl | 228 +++++++++++++++++++++++++ src/solvers/coefficients.jl | 125 ++++++++++++++ src/solvers/common.jl | 48 ++++++ src/solvers/lie_trotter_godunov.jl | 47 +++++ src/solvers/strang_marchuk.jl | 75 ++++++++ src/solvers/tables.jl | 93 ++++++++++ src/utils.jl | 43 +++-- 14 files changed, 889 insertions(+), 61 deletions(-) create mode 100644 src/solvers/adjoint_pair.jl create mode 100644 src/solvers/coefficients.jl create mode 100644 src/solvers/common.jl create mode 100644 src/solvers/lie_trotter_godunov.jl create mode 100644 src/solvers/strang_marchuk.jl create mode 100644 src/solvers/tables.jl diff --git a/_typos.toml b/_typos.toml index e3bfffd..e73ba1a 100644 --- a/_typos.toml +++ b/_typos.toml @@ -5,5 +5,8 @@ Strang = "Strang" Tro = "Tro" # David Ketcheson, in the citation key AuzHofKetKoc:2017:psm Ket = "Ket" +# Citation keys for the order barrier: Sheng (She:1989:slp), Suzuki (Suz:1991:gtf) +She = "She" +Suz = "Suz" # Splitting variable name BA = "BA" diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index de7763c..c4654ab 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -88,3 +88,27 @@ @article{AuzHofKetKoc:2017:psm doi={10.1007/s10543-016-0626-9}, publisher={Springer} } + +@article{She:1989:slp, + title={Solving linear partial differential equations by exponential splitting}, + author={Sheng, Qin}, + journal={IMA Journal of Numerical Analysis}, + volume={9}, + number={2}, + pages={199--212}, + year={1989}, + doi={10.1093/imanum/9.2.199}, + publisher={Oxford University Press} +} + +@article{Suz:1991:gtf, + title={General theory of fractal path integrals with applications to many-body theories and statistical physics}, + author={Suzuki, Masuo}, + journal={Journal of Mathematical Physics}, + volume={32}, + number={2}, + pages={400--407}, + year={1991}, + doi={10.1063/1.529425}, + publisher={AIP} +} diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md index 518509e..78a675d 100644 --- a/docs/src/devdocs/index.md +++ b/docs/src/devdocs/index.md @@ -146,8 +146,9 @@ end This example is written for exactly two operators. A scheme that works for any number of them loops over `children` with `Unrolled.@unroll`, as the built-in algorithms in -`src/solver.jl` do; a plain `for` loop over the heterogeneously typed tuple would be -type unstable. +`src/solvers/` do; a plain `for` loop over the heterogeneously typed tuple would be +type unstable. `advance_one_child!` above is `src/solvers/common.jl`'s +`_advance_child!`, which the built-in schemes share. ### Adaptive algorithms @@ -177,9 +178,83 @@ end OrdinaryDiffEqCore interface for the estimate; go through them rather than touching the `EEst` field, whose location on the integrator is an implementation detail. -See [`PalindromicPairLieTrotterGodunov`](@ref) in `src/solver.jl` for a complete -example, and [Adaptive time stepping](@ref) for how the two layers of adaptivity -interact. +See [`PalindromicPairLieTrotterGodunov`](@ref) in `src/solvers/adjoint_pair.jl` for a +complete example, and [Adaptive time stepping](@ref) for how the two layers of +adaptivity interact. + +## Coefficient tables + +Schemes of order three and above are not written out as passes over the children. +They are described by a table of coefficients and stepped by one shared traversal, +which lives in `src/solvers/coefficients.jl`; the schemes themselves +(`src/solvers/tables.jl`) are then little more than the table plus two interface +methods. See [Higher order splittings](@ref theory_higher-order) for where the +tables come from. + +```@docs +OrdinaryDiffEqOperatorSplitting.SplittingCoefficients +OrdinaryDiffEqOperatorSplitting.coefficients +OrdinaryDiffEqOperatorSplitting.order +``` + +A table-driven scheme is added by giving it a struct, a table, `coefficients`, +`order`, and an `init_cache` returning the shared +`SplittingCoefficientsCache` — no stepping code of its own: + +```julia +struct MyThirdOrder{AlgTupleType <: Tuple} <: + OrdinaryDiffEqOperatorSplitting.AbstractOperatorSplittingAlgorithm + inner_algs::AlgTupleType +end + +const MY_COEFFICIENTS = OrdinaryDiffEqOperatorSplitting.SplittingCoefficients( + (7 // 24, 2 // 3), (3 // 4, -2 // 3), (-1 // 24, 1 // 1) +) + +OrdinaryDiffEqOperatorSplitting.coefficients(::MyThirdOrder) = MY_COEFFICIENTS +OrdinaryDiffEqOperatorSplitting.order(::MyThirdOrder) = 3 + +function OrdinaryDiffEqOperatorSplitting.init_cache( + f::GenericSplitFunction, alg::MyThirdOrder; + uprev::AbstractArray, u::AbstractVector, +) + return OrdinaryDiffEqOperatorSplitting.SplittingCoefficientsCache( + u, uprev, OrdinaryDiffEqOperatorSplitting.coefficients(alg)) +end +``` + +Only the consistency condition is checked when the table is built, so a table that +constructs successfully can still fail to reach the order it claims; a convergence +test is the only real check. + +Two properties of the traversal are worth knowing when reading or extending it: + +- **A zero coefficient is skipped entirely**, synchronization included, since the + flow is the identity. `Yoshida4`'s last stage relies on this. +- **The first flow of a step is always synchronized.** Strang-Marchuk skips that + sync because its reverse pass ends on operator 1, leaving that child's buffer + current; a general table ends on operator `N`, so operator 1's buffer is stale and + the same shortcut would silently corrupt every step after the first. + +Because a table with negative coefficients steps some children backward, a +table-driven scheme also depends on the direction reversal described under +[Backward sub-steps](@ref devdocs_backward-substeps). + +If the scheme's order is odd, wrapping it in `AdjointPair` makes it adaptive for +free: the adjoint is the same table traversed backwards, so nothing further is +needed from the scheme. + +## [Backward sub-steps](@id devdocs_backward-substeps) + +An inner integrator fixes its direction of integration when it is constructed, and +a scheme with negative coefficients has to step it the other way. Rather than build +two integrators per child, a sub-step against the child's direction reverses the +child in place, steps, and reverses it back: + +```@docs +OrdinaryDiffEqOperatorSplitting.reverse_direction! +OrdinaryDiffEqOperatorSplitting.tstops_and_saveat_heaps +``` ## Dense output diff --git a/docs/src/topics/adaptivity.md b/docs/src/topics/adaptivity.md index c2f28c8..f0a14db 100644 --- a/docs/src/topics/adaptivity.md +++ b/docs/src/topics/adaptivity.md @@ -6,10 +6,14 @@ Two independent layers of a splitting tree can adapt their step sizes: they would outside of this package. Their behavior is governed by the tolerances and options they receive. - **Splitting nodes** adapt the splitting step itself. This requires an algorithm - that produces an error estimate of the splitting error; - [`PalindromicPairLieTrotterGodunov`](@ref) is currently the only one. It advances - with the average of the two mutually reversed Lie-Trotter sequences and uses half - their difference as the local error estimate. + that produces an error estimate of the splitting error, which here means running a + scheme together with its adjoint: the average of the two is the solution and half + their difference is the local error estimate (see + [Adjoint pairs](@ref theory_higher-order)). + [`PalindromicPairLieTrotterGodunov`](@ref) is that construction for the + Lie-Trotter sequence, and [`AdjointPair`](@ref) builds it from any odd-order + scheme — `AdjointPair(Ruth3((Tsit5(), Tsit5())))` is a fourth-order adaptive + splitting. Both cost two sweeps per step. By default every node adapts exactly if its own algorithm can: for diff --git a/docs/src/topics/time-integration.md b/docs/src/topics/time-integration.md index 9ecb7f5..69e62ea 100644 --- a/docs/src/topics/time-integration.md +++ b/docs/src/topics/time-integration.md @@ -165,6 +165,114 @@ $[L_1, L_2] = L_1 L_2 - L_2 L_1$ to cancel, leaving a local truncation error of $O(t^3)$ and hence second-order global accuracy. The same argument extends to the general $N$-operator palindromic scheme. +## [Higher order splittings](@id theory_higher-order) + +Lie-Trotter-Godunov and Strang-Marchuk are both instances of a more general +construction. A splitting scheme applies the sub-flows in a fixed sequence, each for +a fixed fraction of the step, so for $N$ operators it is completely described by an +$S \times N$ table of coefficients $a_{ji}$: stage $j$ advances operator $i$ by +$a_{ji} \Delta t$. In the linear two-operator case, + +```math +\mathcal{S}(\Delta t) = \prod_{j=1}^{S} e^{a_{jN} \Delta t L_N} \cdots e^{a_{j1} \Delta t L_1} \, . +``` + +Lie-Trotter-Godunov is the one-stage table $a = (1, 1)$ and Strang-Marchuk is the +two-stage table $a = \bigl((\tfrac{1}{2}, 1), (\tfrac{1}{2}, 0)\bigr)$. + +Each operator's coefficients must sum to one, + +```math +\sum_{j=1}^{S} a_{ji} = 1 \quad \text{for every } i \, , +``` + +since otherwise the scheme does not even advance every sub-problem by $\Delta t$. +This is the consistency condition, and it is the only one this package checks when a +table is constructed. Attaining order $p$ imposes further conditions, one per +independent commutator up to order $p$, obtained by matching the +Baker-Campbell-Hausdorff expansion of the product above against that of +$e^{\Delta t (L_1 + L_2)}$ — the same computation as in the two analyses above, +carried further. + +### Composition: the triple jump + +Solving the order conditions directly gets unpleasant quickly. A cheaper route is +*composition*: build a higher-order scheme out of a symmetric one of lower order. +If $\mathcal{S}_2$ is any symmetric second-order scheme, then + +```math +\mathcal{S}_4(\Delta t) = \mathcal{S}_2(w_1 \Delta t) \, \mathcal{S}_2(w_0 \Delta t) \, \mathcal{S}_2(w_1 \Delta t) +``` + +is symmetric for any weights, and hence of even order. It is of order four as soon +as the weights satisfy + +```math +2 w_1 + w_0 = 1 \, , \qquad 2 w_1^3 + w_0^3 = 0 \, , +``` + +the first being consistency and the second the cancellation of the third-order term. +The real solution is $w_1 = 1/(2 - 2^{1/3})$ and $w_0 = -2^{1/3} w_1$, giving +Yoshida's "triple jump" [Yos:1990:cho](@cite), implemented here as +[`Yoshida4`](@ref). Writing the three Strang steps out as a flat sequence of flows +and merging the adjacent flows of the same operator that the composition leaves next +to each other collapses nine flows to eight — which is exactly the four-stage table +`Yoshida4` carries, its last stage having a zero second coefficient. + +### The order barrier and negative coefficients + +Note that $w_0 < 0$ above. This is not an artifact of the construction: no splitting +scheme of order greater than two has all coefficients positive +[She:1989:slp,Suz:1991:gtf](@cite). Any third- or higher-order splitting therefore +integrates some sub-problem *backward in time* during part of every step, which has +two practical consequences. + +First, the sub-problems must admit a backward flow. For a parabolic sub-problem — +diffusion, say — the backward evolution is ill-posed and the negative sub-steps are +violently unstable, so on a reaction-diffusion system the higher-order schemes here +are not usable on the diffusion operator, however attractive their order. This is +the reason Strang-Marchuk remains the workhorse despite being only second order. + +Second, the implementation has to actually run its sub-integrators backwards. An +inner integrator fixes its direction of integration at construction, so a negative +sub-step temporarily reverses it; see the developer documentation for the details. + +### Adjoint pairs + +The *adjoint* of a scheme is + +```math +\mathcal{S}^*(\Delta t) = \mathcal{S}(-\Delta t)^{-1} \, , +``` + +which for a splitting scheme is simply its whole sequence of flows run in reverse +order, every coefficient keeping its sign and its operator. A scheme is symmetric +exactly when $\mathcal{S}^* = \mathcal{S}$, which is why Strang-Marchuk — a +palindrome — gains an order over Lie-Trotter-Godunov. + +If $\mathcal{S}$ has order $p$ with leading local error $C \Delta t^{p+1}$, then +$\mathcal{S}^*$ has the same order with leading error $(-1)^p C \Delta t^{p+1}$. For +**odd** $p$ the two signs oppose, so running the pair from the same initial value +gives, at twice the cost of one scheme, + +```math +\frac{\mathcal{S} + \mathcal{S}^*}{2} \quad \text{of order } p+1 \, , +\qquad +\frac{\mathcal{S} - \mathcal{S}^*}{2} \quad \text{an estimate of the local error of } \mathcal{S} \, , +``` + +the latter being asymptotically correct as $\Delta t \to 0$ +[AuzHofKetKoc:2017:psm](@cite). This is the Milne device applied to a scheme and its +adjoint, and it is what makes the splitting error itself estimable and hence the +splitting step adaptive — see [Adaptive time stepping](@ref). The construction is +[`AdjointPair`](@ref); at $p = 1$, with Lie-Trotter-Godunov as the base, it is the +pair of mutually reversed sequences implemented directly as +[`PalindromicPairLieTrotterGodunov`](@ref). + +For even $p$ the two leading terms are *equal* rather than opposite: averaging +cancels nothing and the difference is not an error estimate, which is why +[`AdjointPair`](@ref) rejects an even-order base scheme. + ## References ```@bibliography diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index 60e1cad..593b06b 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -150,7 +150,12 @@ include("function.jl") include("config_tree.jl") include("problem.jl") include("integrator.jl") -include("solver.jl") +include("solvers/common.jl") +include("solvers/coefficients.jl") +include("solvers/lie_trotter_godunov.jl") +include("solvers/strang_marchuk.jl") +include("solvers/tables.jl") +include("solvers/adjoint_pair.jl") include("utils.jl") export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, StrangMarchuk, diff --git a/src/integrator.jl b/src/integrator.jl index 04f36ea..59a6e62 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -134,8 +134,6 @@ end # --- SplitSubIntegrator interface --- -tdir(integrator::SplitSubIntegrator) = sign(integrator.dt) - # proposed-dt interface (mirrors ODEIntegrator) function SciMLBase.set_proposed_dt!(sub::SplitSubIntegrator, dt) if sub.dtcache != dt # only touch if actually changing @@ -306,9 +304,6 @@ function SciMLBase.__init( tstops = () end - # Heaps store raw times and carry the integration direction in their ordering. - # (OrdinaryDiffEqCore's initialize_tstops stores tdir-scaled times instead, which - # is incompatible with the heaps reinit! rebuilds and breaks backward tspans.) tstops_internal, saveat_internal = tstops_and_saveat_heaps( t0, tf, (tstops..., d_discontinuities...), saveat ) @@ -447,8 +442,8 @@ function DiffEqBase.reinit!( reinit_callbacks = true, reinit_retcode = true ) - # The heap ordering types and every node's tdir are fixed at init, so a reinit! - # cannot flip the integration direction. + # Every node's tdir is fixed at init, and the heaps it hands to its children are + # scaled by it, so a reinit! cannot flip the integration direction. (tf > t0) == (integrator.tdir > 0) || error("reinit! cannot change the direction of integration. Build a new integrator instead.") @@ -596,22 +591,20 @@ end # --------------------------------------------------------------------------- function _handle_tstop!(integrator::AnySplitIntegrator) if SciMLBase.has_tstop(integrator) - # The heaps store raw times; comparisons happen in tdir-space so that - # "ahead"/"behind" is direction independent. - tdir_t = tdir(integrator) * integrator.t - tdir_tstop = tdir(integrator) * SciMLBase.first_tstop(integrator) + tdir_t = integrator.tdir * integrator.t + tdir_tstop = SciMLBase.first_tstop(integrator) if tdir_t == tdir_tstop while tdir_t == tdir_tstop SciMLBase.pop_tstop!(integrator) SciMLBase.has_tstop(integrator) ? - (tdir_tstop = tdir(integrator) * SciMLBase.first_tstop(integrator)) : break + (tdir_tstop = SciMLBase.first_tstop(integrator)) : break end notify_integrator_hit_tstop!(integrator) elseif tdir_t > tdir_tstop if !integrator.dtchangeable SciMLBase.change_t_via_interpolation!( integrator, - SciMLBase.pop_tstop!(integrator), + pop_next_tstop!(integrator), Val{true} ) notify_integrator_hit_tstop!(integrator) @@ -772,7 +765,7 @@ end function modify_dt_for_tstops!(integrator) if SciMLBase.has_tstop(integrator) tdir_t = integrator.tdir * integrator.t - tdir_tstop = integrator.tdir * SciMLBase.first_tstop(integrator) + tdir_tstop = SciMLBase.first_tstop(integrator) if integrator.opts.adaptive integrator.dt = integrator.tdir * min(abs(integrator.dt), abs(tdir_tstop - tdir_t)) # step! to the end @@ -814,7 +807,7 @@ _snap_window(t, tstop, dt) = function fixed_t_for_floatingpoint_error!(integrator::AnySplitIntegrator, ttmp) return if DiffEqBase.has_tstop(integrator) - tstop = DiffEqBase.first_tstop(integrator) + tstop = next_tstop(integrator) if abs(ttmp - tstop) < _snap_window(integrator.t, tstop, integrator.dt) try_snap_children_to_tstop!.( integrator.child_subintegrators, tstop, integrator.dt @@ -926,8 +919,7 @@ end function DiffEqBase.solve!(integrator::OperatorSplittingIntegrator) while !isempty(integrator.tstops) - while tdir(integrator) * integrator.t < - tdir(integrator) * SciMLBase.first_tstop(integrator) + while integrator.tdir * integrator.t < SciMLBase.first_tstop(integrator) step_header!(integrator) @timeit_debug "check_error" SciMLBase.check_error!(integrator) ∉ ( ReturnCode.Success, ReturnCode.Default, @@ -947,7 +939,7 @@ end function DiffEqBase.step!(integrator::AnySplitIntegrator) @timeit_debug "step!" if integrator.advance_to_tstop - tstop = SciMLBase.first_tstop(integrator) + tstop = next_tstop(integrator) while !reached_tstop(integrator, tstop) step_header!(integrator) @timeit_debug "check_error" SciMLBase.check_error!(integrator) ∉ ( @@ -981,7 +973,7 @@ end # integration. function DiffEqBase.step!(integrator::AnySplitIntegrator, dt, stop_at_tdt = false) @timeit_debug "step!" begin - tdir(integrator) * dt < zero(dt) && error("Cannot step backward.") + integrator.tdir * dt < zero(dt) && error("Cannot step backward.") stop_at_tdt && !integrator.dtchangeable && error("Cannot stop at t + dt if dtchangeable is false") tnext = integrator.t + dt @@ -1243,12 +1235,9 @@ function SciMLBase.savevalues!( return saved, savedexactly 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) + while !isempty(integrator.saveat) && first(integrator.saveat) <= tdir_t + curt = integrator.tdir * pop!(integrator.saveat) if curt == integrator.t # `save_end` owns the final point; leave it to the postamble. (!save_end && curt == tf) && continue @@ -1492,8 +1481,13 @@ end # Time helpers -tdir(integrator) = - integrator.tstops.ordering isa BinaryHeaps.FasterForward ? 1 : -1 +# +# `tstops`/`saveat` keys are tdir-scaled (see `tstops_and_saveat_heaps`), so a raw time +# is compared against them as `tdir * t`, and a key is turned back into a time by the +# same multiplication. +tdir(integrator::AnySplitIntegrator) = integrator.tdir +next_tstop(integrator) = integrator.tdir * SciMLBase.first_tstop(integrator) +pop_next_tstop!(integrator) = integrator.tdir * SciMLBase.pop_tstop!(integrator) is_past_t(integrator, t) = tdir(integrator) * (t - integrator.t) ≤ zero(integrator.t) function reached_tstop(integrator, tstop, stop_at_tstop = integrator.dtchangeable) @@ -1555,7 +1549,7 @@ function advance_solution_by!(integrator::AnySplitIntegrator, dt) return advance_solution_by!(integrator, integrator.cache, dt) end -# Algorithm-level dispatch (implemented in solver.jl per algorithm) +# Algorithm-level dispatch (implemented per algorithm under src/solvers/) function advance_solution_by!( integrator::AnySplitIntegrator, cache::AbstractOperatorSplittingCache, dt @@ -1885,7 +1879,7 @@ function DiffEqBase.add_tstop!(i::AnySplitIntegrator, t) error("Cannot add a tstop at $t because that is behind the current \ integrator time $(i.t)") DiffEqBase.add_tstop!.(i.child_subintegrators, t) - push!(i.tstops, t) + push!(i.tstops, i.tdir * t) return nothing end @@ -1893,7 +1887,7 @@ function DiffEqBase.add_saveat!(i::OperatorSplittingIntegrator, t) is_past_t(i, t) && error("Cannot add a saveat point at $t because that is behind the \ current integrator time $(i.t)") - push!(i.saveat, t) + push!(i.saveat, i.tdir * t) return nothing end diff --git a/src/solvers/adjoint_pair.jl b/src/solvers/adjoint_pair.jl new file mode 100644 index 0000000..fe3164b --- /dev/null +++ b/src/solvers/adjoint_pair.jl @@ -0,0 +1,228 @@ +# --------------------------------------------------------------------------- +# Pairs of a scheme and its adjoint +# +# Running a scheme together with its adjoint from the same initial value gives both a +# solution one order higher and an error estimate, which is what makes these the +# adaptive splitting schemes of the package. See [AuzHofKetKoc:2017:psm](@cite), §3. +# --------------------------------------------------------------------------- + +""" + _pair_average_and_estimate!(parent, uforward, dt) + +Combine the two members of a scheme/adjoint pair, `parent.u` and `uforward`, into the +higher-order solution and, when the node is adaptive, the local error estimate. + +This is the Milne device of [AuzHofKetKoc:2017:psm](@cite), §3: for a base scheme of +odd order `p` the two members have leading error terms `±C h^{p+1}`, so their average +is of order `p+1` and half their difference estimates the base scheme's local error. +`uforward` is overwritten with the scaled residual. +""" +function _pair_average_and_estimate!(parent, uforward, dt) + parent.u .= (parent.u .+ uforward) ./ 2 + if parent.controller_cache !== nothing + (; abstol, reltol, internalnorm) = parent.opts + @. uforward = (parent.u - uforward) / + (abstol + max(abs(parent.u), abs(parent.uprev)) * reltol) + OrdinaryDiffEqCore.set_EEst!(parent, internalnorm(uforward, parent.t + dt)) + end + return +end + +# --------------------------------------------------------------------------- +# Palindromic pair of Lie-Trotter-Godunov splittings +# --------------------------------------------------------------------------- +""" + PalindromicPairLieTrotterGodunov <: AbstractOperatorSplittingAlgorithm + +Second-order sequential operator splitting algorithm. + +One step solves the palindromic pair of [`LieTrotterGodunov`](@ref) sequences + +``A_1(\\Delta t) \\to \\cdots \\to A_N(\\Delta t)`` and +``A_N(\\Delta t) \\to \\cdots \\to A_1(\\Delta t)`` + +from the same initial value. The leading splitting error of a Lie-Trotter sequence +is ``\\frac{\\Delta t^2}{2}\\sum_{i ") + +order(::PalindromicPairLieTrotterGodunov) = 2 + +@inline SciMLBase.isadaptive(::PalindromicPairLieTrotterGodunov) = true +# The pair difference estimates the O(dt²) leading error term of a first order +# sequence, so the controller sees a first order error estimator. +alg_adaptive_order(::PalindromicPairLieTrotterGodunov) = 1 + +struct PalindromicPairLieTrotterGodunovCache{uType, uprevType, uforwardType} <: AbstractOperatorSplittingCache + u::uType + uprev::uprevType + uforward::uforwardType # end state of the A₁ → A₂ sequence; reused as the residual buffer +end + +function init_cache( + f::GenericSplitFunction, alg::PalindromicPairLieTrotterGodunov; + uprev::AbstractArray, u::AbstractVector, + ) + return PalindromicPairLieTrotterGodunovCache(u, uprev, similar(u)) +end + +# Forward sequence: A₁(dt) → … → A_N(dt) +@unroll function _ppltg_forward_pass!(parent, children::Tuple, dt) + i = 0 + @unroll for child in children + i += 1 + _advance_child!(parent, child, i, dt) + parent.force_stepfail && return + end +end + +# Reverse sequence: A_N(dt) → … → A₁(dt) +@unroll function _ppltg_reverse_pass!(parent, rchildren::Tuple, dt, N) + j = 0 + @unroll for child in rchildren + j += 1 + _advance_child!(parent, child, N + 1 - j, dt) + parent.force_stepfail && return + end +end + +function _perform_step!( + parent, + children::Tuple, + cache::PalindromicPairLieTrotterGodunovCache, + dt + ) + (; uforward) = cache + + _ppltg_forward_pass!(parent, children, dt) + parent.force_stepfail && return + uforward .= parent.u + + # Rewind to the initial state of the step; uprev is untouched while stepping. + parent.u .= parent.uprev + rollback_children!(parent) + + _ppltg_reverse_pass!(parent, reverse(children), dt, length(children)) + parent.force_stepfail && return + + _pair_average_and_estimate!(parent, uforward, dt) + return +end + +# --------------------------------------------------------------------------- +# Adjoint pairs of a table-driven scheme +# --------------------------------------------------------------------------- +""" + AdjointPair(base) <: AbstractOperatorSplittingAlgorithm + +Adaptive splitting scheme of order `p+1` built from a `base` scheme of odd order `p`, +following [AuzHofKetKoc:2017:psm](@cite), eq. (3.2). + +One step runs `base` and its adjoint ``\\mathcal{S}^*`` from the same initial value. +Their leading error terms are ``C h^{p+1}`` and ``(-1)^p C h^{p+1}``, so for odd `p` +the signs oppose: the average is a solution of order `p+1`, and half the difference is +an asymptotically correct estimate of the base scheme's local error, which drives the +step size controller. A step therefore costs twice the base scheme. + +The base scheme's order must be **odd**. For even `p` the two leading terms are +*equal* rather than opposite, so averaging cancels nothing and the difference stops +being an error estimate. + +``\\mathcal{S}^*(h, u) = \\mathcal{S}^{-1}(-h, u)`` is the base scheme's entire flat +sequence of flows reversed, every coefficient keeping its sign and its operator, so it +reuses the same table and needs no extra coefficients. + +[`PalindromicPairLieTrotterGodunov`](@ref) is this construction at `p = 1`. + +As everywhere here, the order and the estimate cover the *splitting* error only and +presume the inner solvers resolve their subproblems accurately relative to it. +""" +struct AdjointPair{BaseType, AlgTupleType <: Tuple} <: AbstractOperatorSplittingAlgorithm + base::BaseType + inner_algs::AlgTupleType # aliases `base.inner_algs`, so the tree machinery works unchanged + + function AdjointPair(base::AbstractOperatorSplittingAlgorithm) + p = order(base) + isodd(p) || throw( + ArgumentError( + "AdjointPair needs a base scheme of odd order, got order $p. For even \ + orders a scheme and its adjoint share the same leading error term, so \ + averaging them raises no order and their difference is not an error \ + estimate." + ) + ) + return new{typeof(base), typeof(base.inner_algs)}(base, base.inner_algs) + end +end + +function Base.show(io::IO, alg::AdjointPair) + print(io, "AdjointPair (") + Base.show(io, alg.base) + return print(io, ")") +end + +coefficients(alg::AdjointPair) = coefficients(alg.base) +order(alg::AdjointPair) = order(alg.base) + 1 + +@inline SciMLBase.isadaptive(::AdjointPair) = true +# The estimate measures the *base* scheme's leading error, not the pair's. +alg_adaptive_order(alg::AdjointPair) = order(alg.base) + +struct AdjointPairCache{uType, uprevType, uforwardType, coeffType} <: AbstractOperatorSplittingCache + u::uType + uprev::uprevType + uforward::uforwardType # end state of the base sequence; reused as the residual buffer + coeffs::coeffType +end + +function init_cache( + f::GenericSplitFunction, alg::AdjointPair; + uprev::AbstractArray, u::AbstractVector, + ) + return AdjointPairCache(u, uprev, similar(u), coefficients(alg)) +end + +function _perform_step!( + parent, + children::Tuple, + cache::AdjointPairCache, + dt + ) + (; uforward, coeffs) = cache + + _table_stages!(parent, children, coeffs.a, dt) + parent.force_stepfail && return + uforward .= parent.u + + # Rewind to the initial state of the step; uprev is untouched while stepping. + parent.u .= parent.uprev + rollback_children!(parent) + + _table_stages_adjoint!( + parent, reverse(children), reverse(coeffs.a), dt, length(children) + ) + parent.force_stepfail && return + + _pair_average_and_estimate!(parent, uforward, dt) + return +end diff --git a/src/solvers/coefficients.jl b/src/solvers/coefficients.jl new file mode 100644 index 0000000..48e6f63 --- /dev/null +++ b/src/solvers/coefficients.jl @@ -0,0 +1,125 @@ +# --------------------------------------------------------------------------- +# Coefficient tables +# +# A splitting scheme of order three or above is specified by a table of coefficients +# rather than by hand-written passes. This file holds the table type, the interface +# every table-driven scheme implements (`coefficients` and `order`), and the two +# traversals of a table: the scheme itself and its adjoint. +# --------------------------------------------------------------------------- + +""" + SplittingCoefficients(stages::NTuple{N, T}...) + +Coefficients of an `S`-stage splitting scheme over `N` operators, one tuple per stage: +stage `j` advances operator `i` by `stages[j][i] * dt`. + +This is the generalization to `N` operators of the two-operator (`AB`) and +three-operator (`ABC`) coefficient tables of +[AuzHofKetKoc:2017:psm](@cite); their tables are the `N = 2` and `N = 3` cases. + +Each operator's coefficients must sum to one, the consistency condition, and that is +checked here. The remaining order conditions are not, so a table that constructs +successfully can still fail to attain the order it claims. +""" +struct SplittingCoefficients{S, N, T} + a::NTuple{S, NTuple{N, T}} + + # Do not simplify this to `NTuple{S, NTuple{N, T}}`: the empty tuple matches it for + # any element type, leaving parameters unbound for `S == 0` (and `T` unbound for + # `N == 0`). A leading element plus a counted `Vararg` rules both out. + function SplittingCoefficients( + stage1::Tuple{T, Vararg{T, K}}, + rest::Tuple{T, Vararg{T, K}}... + ) where {T, K} + a = (stage1, rest...) + N = K + 1 + S = length(a) + for i in 1:N + total = sum(a[j][i] for j in 1:S) + total ≈ one(T) || throw( + ArgumentError( + "operator $i's coefficients sum to $total rather than 1, so the \ + scheme is not consistent." + ) + ) + end + return new{S, N, T}(a) + end +end + +""" + coefficients(alg) + +The [`SplittingCoefficients`](@ref) table of a coefficient-driven splitting algorithm. +""" +function coefficients end + +""" + order(alg) + +Order of consistency of a splitting algorithm, counting the splitting error only. +""" +function order end + +struct SplittingCoefficientsCache{uType, uprevType, coeffType} <: AbstractOperatorSplittingCache + u::uType + uprev::uprevType + coeffs::coeffType +end + +function _perform_step!( + parent, + children::Tuple, + cache::SplittingCoefficientsCache, + dt + ) + # Deliberately no `mark_next_sync_continuous`: that shortcut needs the previous step + # to have left `parent.u` equal to the buffer of the child solved first, which holds + # for StrangMarchuk only because its reverse pass ends on operator 1. A general table + # ends on operator N, so operator 1's buffer is stale and skipping its forward sync + # resumes from stale state -- the first step stays exact and every later one is wrong. + _table_stages!(parent, children, cache.coeffs.a, dt) + return +end + +@unroll function _table_stages!(parent, children, stages::Tuple, dt) + @unroll for stage in stages + _table_stage!(parent, children, stage, dt) + parent.force_stepfail && return + end +end + +@unroll function _table_stage!(parent, children::Tuple, stage, dt) + i = 0 + @unroll for child in children + i += 1 + coefficient = stage[i] + # A zero coefficient is the identity flow, so it is skipped sync and all. + if !iszero(coefficient) + _advance_child!(parent, child, i, coefficient * dt) + parent.force_stepfail && return + end + end +end + +# The adjoint reverses the whole flat sequence of flows, so both the stage order and +# the operator order within each stage are reversed. +@unroll function _table_stages_adjoint!(parent, rchildren, rstages::Tuple, dt, N) + @unroll for stage in rstages + _table_stage_adjoint!(parent, rchildren, stage, dt, N) + parent.force_stepfail && return + end +end + +@unroll function _table_stage_adjoint!(parent, rchildren::Tuple, stage, dt, N) + j = 0 + @unroll for child in rchildren + j += 1 + i = N + 1 - j + coefficient = stage[i] + if !iszero(coefficient) + _advance_child!(parent, child, i, coefficient * dt) + parent.force_stepfail && return + end + end +end diff --git a/src/solvers/common.jl b/src/solvers/common.jl new file mode 100644 index 0000000..0f01fca --- /dev/null +++ b/src/solvers/common.jl @@ -0,0 +1,48 @@ +# --------------------------------------------------------------------------- +# Pieces shared by every splitting scheme +# --------------------------------------------------------------------------- + +""" + _advance_child!(parent, child, i, dt) + +Advance the `i`-th operator of `parent` by `dt`, syncing state into the child before +and out of it afterwards. A failed child is reported through `parent.force_stepfail`, +which every scheme checks between flows. +""" +function _advance_child!(parent, child, i, dt) + idxs = parent.child_solution_indices[i] + sync = parent.child_synchronizers[i] + + @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync) + @timeit_debug "time solve" advance_solution_by!(parent, child, dt) + if child_failed(child) + parent.force_stepfail = true + return + end + + @timeit_debug "sync <-" backward_sync_subintegrator!(parent, child, idxs, sync) + return +end + +function _require_two_operators(scheme, inner_algs) + n = length(inner_algs) + n == 2 || throw( + ArgumentError( + "$scheme is a two-operator (AB) table but got $n operators. Group the \ + operators into a nested GenericSplitFunction to use it with more." + ) + ) + return nothing +end + +# Every scheme prints as `NAME (inner -> inner)`; `separator` is what sits between the +# inner algorithms. +function _show_scheme(io::IO, name, inner_algs, separator = " -> ") + print(io, name, " (") + for inner_alg in inner_algs[1:(end - 1)] + Base.show(io, inner_alg) + print(io, separator) + end + length(inner_algs) > 0 && Base.show(io, inner_algs[end]) + return print(io, ")") +end diff --git a/src/solvers/lie_trotter_godunov.jl b/src/solvers/lie_trotter_godunov.jl new file mode 100644 index 0000000..080f89b --- /dev/null +++ b/src/solvers/lie_trotter_godunov.jl @@ -0,0 +1,47 @@ +# --------------------------------------------------------------------------- +# Lie-Trotter-Godunov operator splitting +# --------------------------------------------------------------------------- +""" + LieTrotterGodunov <: AbstractOperatorSplittingAlgorithm + +First-order sequential operator splitting algorithm attributed to +[Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite). +""" +struct LieTrotterGodunov{AlgTupleType} <: AbstractOperatorSplittingAlgorithm + inner_algs::AlgTupleType # Tuple of timesteppers for inner problems +end + +Base.show(io::IO, alg::LieTrotterGodunov) = _show_scheme(io, "LTG", alg.inner_algs) + +struct LieTrotterGodunovCache{uType, uprevType} <: AbstractOperatorSplittingCache + u::uType + uprev::uprevType +end + +function init_cache( + f::GenericSplitFunction, alg::LieTrotterGodunov; + uprev::AbstractArray, u::AbstractVector, + ) + return LieTrotterGodunovCache(u, uprev) +end + +@unroll function _perform_step!( + parent, + children::Tuple, + cache::LieTrotterGodunovCache, + dt + ) + i = 0 + @unroll for child in children + i += 1 + _advance_child!(parent, child, i, dt) + parent.force_stepfail && return + end +end + +# Lie-Trotter keeps its hand-written step; the table exists only so it can serve as an +# `AdjointPair` base, which is what makes `AdjointPair(LieTrotterGodunov(...))` and +# `PalindromicPairLieTrotterGodunov` the same scheme. +order(::LieTrotterGodunov) = 1 +coefficients(alg::LieTrotterGodunov) = + SplittingCoefficients(ntuple(_ -> 1 // 1, length(alg.inner_algs))) diff --git a/src/solvers/strang_marchuk.jl b/src/solvers/strang_marchuk.jl new file mode 100644 index 0000000..30b12ed --- /dev/null +++ b/src/solvers/strang_marchuk.jl @@ -0,0 +1,75 @@ +# --------------------------------------------------------------------------- +# Strang-Marchuk operator splitting +# --------------------------------------------------------------------------- +""" + StrangMarchuk <: AbstractOperatorSplittingAlgorithm + +Second-order symmetric (palindromic) operator splitting algorithm attributed to +[Str:1968:ccd,Mar:1971:tsm](@cite). + +For ``N`` operators the scheme performs + +``A_1(\\Delta t/2) \\to \\cdots \\to A_{N-1}(\\Delta t/2) \\to A_N(\\Delta t) \\to A_{N-1}(\\Delta t/2) \\to \\cdots \\to A_1(\\Delta t/2)`` + +achieving second-order accuracy through symmetry. +""" +struct StrangMarchuk{AlgTupleType} <: AbstractOperatorSplittingAlgorithm + inner_algs::AlgTupleType # Tuple of timesteppers for inner problems +end + +Base.show(io::IO, alg::StrangMarchuk) = _show_scheme(io, "SM", alg.inner_algs) + +order(::StrangMarchuk) = 2 + +struct StrangMarchukCache{uType, uprevType} <: AbstractOperatorSplittingCache + u::uType + uprev::uprevType +end + +function init_cache( + f::GenericSplitFunction, alg::StrangMarchuk; + uprev::AbstractArray, u::AbstractVector, + ) + return StrangMarchukCache(u, uprev) +end + +# Forward pass: A₁(dt/2) → … → Aₙ₋₁(dt/2) → Aₙ(dt) +@unroll function _sm_forward_pass!(parent, children::Tuple, half_dt, dt) + N = length(children) + i = 0 + @unroll for child in children + i += 1 + _advance_child!(parent, child, i, i < N ? half_dt : dt) + parent.force_stepfail && return + end +end + +# Reverse pass: Aₙ₋₁(dt/2) → … → A₁(dt/2) +@unroll function _sm_reverse_pass!(parent, rev_front::Tuple, half_dt, N) + j = 0 + @unroll for child in rev_front + j += 1 + _advance_child!(parent, child, N - j, half_dt) + parent.force_stepfail && return + end +end + +function _perform_step!( + parent, + children::Tuple, + cache::StrangMarchukCache, + dt + ) + half_dt = dt / 2 + + # Skip sync of for first solve, because it is already in sync + mark_next_sync_continuous(parent) + + _sm_forward_pass!(parent, children, half_dt, dt) + parent.force_stepfail && return + + _sm_reverse_pass!(parent, reverse(children[1:(end - 1)]), half_dt, length(children)) + parent.force_stepfail && return + + return +end diff --git a/src/solvers/tables.jl b/src/solvers/tables.jl new file mode 100644 index 0000000..68a76c0 --- /dev/null +++ b/src/solvers/tables.jl @@ -0,0 +1,93 @@ +# --------------------------------------------------------------------------- +# Schemes defined by a coefficient table +# +# Each of these is a struct, a table, and the two interface methods; the stepping +# itself is the generic table traversal in coefficients.jl. +# --------------------------------------------------------------------------- +""" + Ruth3 <: AbstractOperatorSplittingAlgorithm + +Third-order splitting scheme of [Rut:1983:cim](@cite), in three stages. + +Its coefficients are exactly rational, and -- as is unavoidable for any real +splitting scheme of order three or above -- some of them are negative, so parts of +the step run backward in time. + +As for every splitting scheme here the order statement covers the *splitting* error +only, and presumes the inner solvers resolve their subproblems accurately relative +to it. +""" +struct Ruth3{AlgTupleType <: Tuple} <: AbstractOperatorSplittingAlgorithm + inner_algs::AlgTupleType + + function Ruth3(inner_algs::Tuple) + _require_two_operators("Ruth3", inner_algs) + return new{typeof(inner_algs)}(inner_algs) + end +end + +Base.show(io::IO, alg::Ruth3) = _show_scheme(io, "Ruth3", alg.inner_algs) + +const RUTH3_COEFFICIENTS = SplittingCoefficients( + (7 // 24, 2 // 3), (3 // 4, -2 // 3), (-1 // 24, 1 // 1) +) + +coefficients(::Ruth3) = RUTH3_COEFFICIENTS +order(::Ruth3) = 3 + +function init_cache( + f::GenericSplitFunction, alg::Ruth3; + uprev::AbstractArray, u::AbstractVector, + ) + return SplittingCoefficientsCache(u, uprev, coefficients(alg)) +end + +""" + Yoshida4 <: AbstractOperatorSplittingAlgorithm + +Fourth-order splitting scheme of [Yos:1990:cho](@cite), the "triple jump". + +Built by composing three Strang steps of lengths ``w_1 h``, ``w_0 h`` and ``w_1 h`` +with ``w_1 = 1/(2 - 2^{1/3})`` and ``w_0 = -2^{1/3} w_1``, then merging the adjacent +flows the composition leaves next to each other. That merging is what makes it eight +flow evaluations rather than nine, and it leaves the last stage's second coefficient +zero. + +``w_0`` is negative, so a substantial part of each step runs backward in time -- the +second operator's cumulative time reaches ``1.35\\,h`` before returning through +``-0.35\\,h``. + +As for every splitting scheme here the order statement covers the *splitting* error +only, and presumes the inner solvers resolve their subproblems accurately relative +to it. +""" +struct Yoshida4{AlgTupleType <: Tuple} <: AbstractOperatorSplittingAlgorithm + inner_algs::AlgTupleType + + function Yoshida4(inner_algs::Tuple) + _require_two_operators("Yoshida4", inner_algs) + return new{typeof(inner_algs)}(inner_algs) + end +end + +Base.show(io::IO, alg::Yoshida4) = _show_scheme(io, "Yoshida4", alg.inner_algs) + +const YOSHIDA4_W1 = 1 / (2 - cbrt(2)) +const YOSHIDA4_W0 = -cbrt(2) * YOSHIDA4_W1 + +const YOSHIDA4_COEFFICIENTS = SplittingCoefficients( + (YOSHIDA4_W1 / 2, YOSHIDA4_W1), + ((YOSHIDA4_W1 + YOSHIDA4_W0) / 2, YOSHIDA4_W0), + ((YOSHIDA4_W1 + YOSHIDA4_W0) / 2, YOSHIDA4_W1), + (YOSHIDA4_W1 / 2, 0.0), +) + +coefficients(::Yoshida4) = YOSHIDA4_COEFFICIENTS +order(::Yoshida4) = 4 + +function init_cache( + f::GenericSplitFunction, alg::Yoshida4; + uprev::AbstractArray, u::AbstractVector, + ) + return SplittingCoefficientsCache(u, uprev, coefficients(alg)) +end diff --git a/src/utils.jl b/src/utils.jl index 63df700..e167ad5 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,15 +1,25 @@ -# helper function for setting up min/max heaps for tstops and saveat +""" + tstops_and_saveat_heaps(t0, tf, tstops, saveat) + +Build the `tstops` and `saveat` heaps of a node covering `t0` to `tf`. + +Both store **`tdir`-scaled** times, the convention OrdinaryDiffEqCore's +`initialize_tstops` follows, so that "the next one" is always the heap minimum and +"ahead of `t`" is always `key > tdir * t`, whichever way the node integrates. It is +also what lets [`reverse_direction!`](@ref) flip a node by re-signing its keys: the +alternative, storing raw times and encoding the direction in the heap's ordering, +puts the direction in the heap's *type*, where it cannot be changed in place. +""" function tstops_and_saveat_heaps(t0, tf, tstops, saveat) FT = typeof(tf) - ordering = tf > t0 ? BinaryHeaps.FasterForward : BinaryHeaps.FasterReverse + tdir = tf > t0 ? one(FT) : -one(FT) # ensure that tstops includes tf and only has values ahead of t0 tstops = [filter(t -> t0 < t < tf || tf < t < t0, tstops)..., tf] - tstops = BinaryHeaps.BinaryHeap{FT, ordering}(tstops) + tstops = BinaryHeaps.BinaryHeap{FT, BinaryHeaps.FasterForward}(tdir .* tstops) # 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 @@ -19,7 +29,7 @@ function tstops_and_saveat_heaps(t0, tf, tstops, saveat) else filter(t -> tdir * t0 < tdir * t <= tdir * tf, collect(FT, saveat)) end - saveat = BinaryHeaps.BinaryHeap{FT, ordering}(saveat) + saveat = BinaryHeaps.BinaryHeap{FT, BinaryHeaps.FasterForward}(tdir .* saveat) return tstops, saveat end @@ -272,6 +282,10 @@ The `tstops`/`saveat`/`d_discontinuities` heaps store `tdir`-scaled times, so ne every key re-expresses the same raw times under the new direction -- which inverts the heap order, hence the rebuild. Times now *behind* are dropped, because `handle_tstop!` errors on an unconsumed stop that `t` has passed. + +A `SplitSubIntegrator` follows the same convention (see +[`tstops_and_saveat_heaps`](@ref)) and is reversed the same way; it just has fewer +things to reverse. """ function reverse_direction!(integrator::DEIntegrator) integrator.tdir = -integrator.tdir @@ -291,13 +305,11 @@ function reverse_direction!(integrator::DEIntegrator) end function reverse_direction!(sub::SplitSubIntegrator) + sub.tdir = -sub.tdir sub.dt = -sub.dt sub.dtcache = -sub.dtcache - sub.tdir = -sub.tdir - # This heap stores *raw* times and carries the direction in its ordering, which is - # part of its type, so unlike a leaf's its keys cannot be re-signed. - _drop_times_behind!(sub.tstops, sub.tdir, sub.t) + _reverse_time_heap!(sub.tstops, sub.tdir * sub.t) # `add_tstop!` propagates eagerly to every descendant, so a child left facing the # old direction would reject the reversed node's next tstop as behind it. @@ -312,19 +324,6 @@ end end end -function _drop_times_behind!(heap, tdir, t) - isempty(heap) && return heap - old = [pop!(heap)] - while !isempty(heap) - push!(old, pop!(heap)) - end - threshold = tdir * t - for key in old - tdir * key > threshold && push!(heap, key) - end - return heap -end - function _reverse_time_heap!(heap, threshold) isempty(heap) && return heap # Drain first: the new keys are the negated old ones, so pushing them back into From 8e837e1b2fd5aaddef84eec572cbd656e133e34a Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Tue, 11 Aug 2026 12:37:39 +0200 Subject: [PATCH 6/6] address review --- _typos.toml | 3 +++ docs/src/assets/references.bib | 12 ++++++++++++ docs/src/devdocs/index.md | 11 +++++++++++ docs/src/topics/time-integration.md | 14 ++++++++++---- src/integrator.jl | 17 +++++++---------- src/utils.jl | 6 +----- 6 files changed, 44 insertions(+), 19 deletions(-) diff --git a/_typos.toml b/_typos.toml index e73ba1a..512f4ef 100644 --- a/_typos.toml +++ b/_typos.toml @@ -8,5 +8,8 @@ Ket = "Ket" # Citation keys for the order barrier: Sheng (She:1989:slp), Suzuki (Suz:1991:gtf) She = "She" Suz = "Suz" +# Hansen and Ostermann, in the citation key HanOst:2009:hos +Han = "Han" +Ost = "Ost" # Splitting variable name BA = "BA" diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index c4654ab..6b07def 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -112,3 +112,15 @@ @article{Suz:1991:gtf doi={10.1063/1.529425}, publisher={AIP} } + +@article{HanOst:2009:hos, + title={High order splitting methods for analytic semigroups exist}, + author={Hansen, Eskil and Ostermann, Alexander}, + journal={BIT Numerical Mathematics}, + volume={49}, + number={3}, + pages={527--542}, + year={2009}, + doi={10.1007/s10543-009-0236-x}, + publisher={Springer} +} diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md index 78a675d..3e37e2c 100644 --- a/docs/src/devdocs/index.md +++ b/docs/src/devdocs/index.md @@ -244,6 +244,17 @@ If the scheme's order is odd, wrapping it in `AdjointPair` makes it adaptive for free: the adjoint is the same table traversed backwards, so nothing further is needed from the scheme. +Any published two- or three-operator table with real coefficients — the `AB`/`ABC` +tables of [AuzHofKetKoc:2017:psm](@cite) and the collections derived from them — can +be added this way, needing nothing but the coefficients themselves. + +The [complex-coefficient schemes](@ref theory_higher-order) are the exception. The +traversal would not change for them, but everything downstream of `coefficient * dt` +would: the sub-problems, their states and their inner integrators would all have to +be complex, as would the error norms and step size controllers. That is a property of +the problem being split rather than of the table, so it is not something a table +alone can opt into. + ## [Backward sub-steps](@id devdocs_backward-substeps) An inner integrator fixes its direction of integration when it is constructed, and diff --git a/docs/src/topics/time-integration.md b/docs/src/topics/time-integration.md index 69e62ea..ef63c08 100644 --- a/docs/src/topics/time-integration.md +++ b/docs/src/topics/time-integration.md @@ -222,10 +222,10 @@ to each other collapses nine flows to eight — which is exactly the four-stage ### The order barrier and negative coefficients Note that $w_0 < 0$ above. This is not an artifact of the construction: no splitting -scheme of order greater than two has all coefficients positive -[She:1989:slp,Suz:1991:gtf](@cite). Any third- or higher-order splitting therefore -integrates some sub-problem *backward in time* during part of every step, which has -two practical consequences. +scheme with *real* coefficients of order greater than two has all of them positive +[She:1989:slp,Suz:1991:gtf](@cite). Any third- or higher-order real splitting +therefore integrates some sub-problem *backward in time* during part of every step, +which has two practical consequences. First, the sub-problems must admit a backward flow. For a parabolic sub-problem — diffusion, say — the backward evolution is ill-posed and the negative sub-steps are @@ -233,6 +233,12 @@ violently unstable, so on a reaction-diffusion system the higher-order schemes h are not usable on the diffusion operator, however attractive their order. This is the reason Strang-Marchuk remains the workhorse despite being only second order. +The barrier is a statement about real coefficients only. Allowing *complex* +coefficients with positive real part, high-order splittings do exist for analytic +semigroups [HanOst:2009:hos](@cite): the sub-steps then move along rays into the +complex time plane rather than backward along the real axis, which keeps a parabolic +sub-flow well posed. + Second, the implementation has to actually run its sub-integrators backwards. An inner integrator fixes its direction of integration at construction, so a negative sub-step temporarily reverses it; see the developer documentation for the details. diff --git a/src/integrator.jl b/src/integrator.jl index 59a6e62..9e82bbd 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -604,7 +604,7 @@ function _handle_tstop!(integrator::AnySplitIntegrator) if !integrator.dtchangeable SciMLBase.change_t_via_interpolation!( integrator, - pop_next_tstop!(integrator), + integrator.tdir * SciMLBase.pop_tstop!(integrator), Val{true} ) notify_integrator_hit_tstop!(integrator) @@ -807,7 +807,7 @@ _snap_window(t, tstop, dt) = function fixed_t_for_floatingpoint_error!(integrator::AnySplitIntegrator, ttmp) return if DiffEqBase.has_tstop(integrator) - tstop = next_tstop(integrator) + tstop = integrator.tdir * SciMLBase.first_tstop(integrator) if abs(ttmp - tstop) < _snap_window(integrator.t, tstop, integrator.dt) try_snap_children_to_tstop!.( integrator.child_subintegrators, tstop, integrator.dt @@ -939,7 +939,7 @@ end function DiffEqBase.step!(integrator::AnySplitIntegrator) @timeit_debug "step!" if integrator.advance_to_tstop - tstop = next_tstop(integrator) + tstop = integrator.tdir * SciMLBase.first_tstop(integrator) while !reached_tstop(integrator, tstop) step_header!(integrator) @timeit_debug "check_error" SciMLBase.check_error!(integrator) ∉ ( @@ -1483,16 +1483,13 @@ end # Time helpers # # `tstops`/`saveat` keys are tdir-scaled (see `tstops_and_saveat_heaps`), so a raw time -# is compared against them as `tdir * t`, and a key is turned back into a time by the -# same multiplication. -tdir(integrator::AnySplitIntegrator) = integrator.tdir -next_tstop(integrator) = integrator.tdir * SciMLBase.first_tstop(integrator) -pop_next_tstop!(integrator) = integrator.tdir * SciMLBase.pop_tstop!(integrator) +# is compared against them as `integrator.tdir * t`, and a key is turned back into a +# time by the same multiplication. is_past_t(integrator, t) = - tdir(integrator) * (t - integrator.t) ≤ zero(integrator.t) + integrator.tdir * (t - integrator.t) ≤ zero(integrator.t) function reached_tstop(integrator, tstop, stop_at_tstop = integrator.dtchangeable) if stop_at_tstop - tdir(integrator) * (integrator.t - tstop) > zero(integrator.t) && + integrator.tdir * (integrator.t - tstop) > zero(integrator.t) && error("Integrator missed stop at $tstop (current time=$(integrator.t)). Aborting.") return integrator.t ≈ tstop else diff --git a/src/utils.jl b/src/utils.jl index e167ad5..5fae579 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -266,7 +266,7 @@ function _fix_dt_at_bounds!(integrator::AnySplitIntegrator) # over dtmax if the two conflict. dtmax = abs(integrator.opts.dtmax) dtmin = abs(DiffEqBase.timedepentdtmin(integrator)) - integrator.dt = tdir(integrator) * max(min(abs(integrator.dt), dtmax), dtmin) + integrator.dt = integrator.tdir * max(min(abs(integrator.dt), dtmax), dtmin) return nothing end @@ -282,10 +282,6 @@ The `tstops`/`saveat`/`d_discontinuities` heaps store `tdir`-scaled times, so ne every key re-expresses the same raw times under the new direction -- which inverts the heap order, hence the rebuild. Times now *behind* are dropped, because `handle_tstop!` errors on an unconsumed stop that `t` has passed. - -A `SplitSubIntegrator` follows the same convention (see -[`tstops_and_saveat_heaps`](@ref)) and is reversed the same way; it just has fewer -things to reverse. """ function reverse_direction!(integrator::DEIntegrator) integrator.tdir = -integrator.tdir