diff --git a/_typos.toml b/_typos.toml index 625e683..512f4ef 100644 --- a/_typos.toml +++ b/_typos.toml @@ -3,5 +3,13 @@ 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" +# 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/api-reference/index.md b/docs/src/api-reference/index.md index 44b0610..4662502 100644 --- a/docs/src/api-reference/index.md +++ b/docs/src/api-reference/index.md @@ -18,6 +18,9 @@ GenericSplitFunction LieTrotterGodunov StrangMarchuk PalindromicPairLieTrotterGodunov +Ruth3 +Yoshida4 +AdjointPair ``` ## Per-node configuration diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index 3cb2732..6b07def 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -52,3 +52,75 @@ @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} +} + +@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} +} + +@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 518509e..3e37e2c 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,94 @@ 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. + +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 +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..ef63c08 100644 --- a/docs/src/topics/time-integration.md +++ b/docs/src/topics/time-integration.md @@ -165,6 +165,120 @@ $[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 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 +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. + +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. + +### 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 661dc24..593b06b 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -150,11 +150,16 @@ 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, - PalindromicPairLieTrotterGodunov + PalindromicPairLieTrotterGodunov, Ruth3, Yoshida4, AdjointPair export SplitNode, TreeOption include("precompilation.jl") diff --git a/src/integrator.jl b/src/integrator.jl index e5befd1..9e82bbd 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), + integrator.tdir * SciMLBase.pop_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,9 +807,11 @@ _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 = 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) + try_snap_children_to_tstop!.( + integrator.child_subintegrators, tstop, integrator.dt + ) tstop else ttmp @@ -825,16 +820,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 +853,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 @@ -918,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, @@ -939,7 +939,7 @@ end function DiffEqBase.step!(integrator::AnySplitIntegrator) @timeit_debug "step!" if integrator.advance_to_tstop - tstop = SciMLBase.first_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) ∉ ( @@ -973,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 @@ -1235,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 @@ -1484,13 +1481,15 @@ 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 `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 @@ -1547,7 +1546,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 @@ -1616,11 +1615,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 +1626,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 @@ -1865,7 +1876,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 @@ -1873,7 +1884,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/solver.jl b/src/solver.jl deleted file mode 100644 index db76d35..0000000 --- a/src/solver.jl +++ /dev/null @@ -1,288 +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 _ppltg_advance_child!(parent, child, i, dt) - idxs = parent.child_solution_indices[i] - sync = parent.child_synchronizers[i] - - @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync) - @timeit_debug "time solve" advance_solution_by!(parent, child, dt) - if child_failed(child) - parent.force_stepfail = true - return - end - - @timeit_debug "sync <-" backward_sync_subintegrator!(parent, child, idxs, sync) - return -end - -# Forward sequence: A₁(dt) → … → A_N(dt) -@unroll function _ppltg_forward_pass!(parent, children::Tuple, dt) - i = 0 - @unroll for child in children - i += 1 - _ppltg_advance_child!(parent, child, i, dt) - parent.force_stepfail && return - end -end - -# Reverse sequence: A_N(dt) → … → A₁(dt) -@unroll function _ppltg_reverse_pass!(parent, rchildren::Tuple, dt, N) - j = 0 - @unroll for child in rchildren - j += 1 - _ppltg_advance_child!(parent, child, N + 1 - j, dt) - parent.force_stepfail && return - end -end - -function _perform_step!( - parent, - children::Tuple, - cache::PalindromicPairLieTrotterGodunovCache, - dt - ) - (; uforward) = cache - - _ppltg_forward_pass!(parent, children, dt) - parent.force_stepfail && return - uforward .= parent.u - - # Rewind to the initial state of the step; uprev is untouched while stepping. - parent.u .= parent.uprev - rollback_children!(parent) - - _ppltg_reverse_pass!(parent, reverse(children), dt, length(children)) - parent.force_stepfail && return - - # The average of the pair is the second order solution ... - parent.u .= (parent.u .+ uforward) ./ 2 - if parent.controller_cache !== nothing - # ... and half the pair difference the local error of a single sequence. - (; abstol, reltol, internalnorm) = parent.opts - @. uforward = (parent.u - uforward) / - (abstol + max(abs(parent.u), abs(parent.uprev)) * reltol) - OrdinaryDiffEqCore.set_EEst!(parent, internalnorm(uforward, parent.t + dt)) - end - return -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 e7c059e..5fae579 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 @@ -256,10 +266,75 @@ 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 +""" + 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.tdir = -sub.tdir + sub.dt = -sub.dt + sub.dtcache = -sub.dtcache + + _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. + _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 _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..53bab32 --- /dev/null +++ b/test/adjoint_pair.jl @@ -0,0 +1,88 @@ +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()))) + @test_throws ArgumentError AdjointPair(Yoshida4((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..4433980 --- /dev/null +++ b/test/coefficient_schemes.jl @@ -0,0 +1,84 @@ +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 "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 + # 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..2fab596 100644 --- a/test/convergence.jl +++ b/test/convergence.jl @@ -44,21 +44,29 @@ 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), + ("Yoshida4", Yoshida4, 4, 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 +76,9 @@ end end end end + +@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 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)