diff --git a/docs/src/api-reference/index.md b/docs/src/api-reference/index.md index 1e84973..7c630e8 100644 --- a/docs/src/api-reference/index.md +++ b/docs/src/api-reference/index.md @@ -12,8 +12,9 @@ OperatorSplittingProblem GenericSplitFunction ``` -## Solver +## Solvers ```@docs LieTrotterGodunov +StrangMarchuk ``` diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index 7184cca..3cb2732 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -32,3 +32,23 @@ @article{God:1959:dmn year={1959}, publisher={Russian Academy of Sciences, Steklov Mathematical Institute} } + +@article{Str:1968:ccd, + title={On the construction and comparison of difference schemes}, + author={Strang, Gilbert}, + journal={SIAM Journal on Numerical Analysis}, + volume={5}, + number={3}, + pages={506--517}, + year={1968}, + publisher={SIAM} +} + +@incollection{Mar:1971:tsm, + title={On the theory of the splitting-up method}, + author={Marchuk, Guri Ivanovich}, + booktitle={Numerical Solution of Partial Differential Equations-{II}}, + pages={469--500}, + year={1971}, + publisher={Academic Press} +} diff --git a/docs/src/topics/time-integration.md b/docs/src/topics/time-integration.md index f6fccf1..9ecb7f5 100644 --- a/docs/src/topics/time-integration.md +++ b/docs/src/topics/time-integration.md @@ -116,6 +116,55 @@ $n \in \mathbb{N}$ the following bound which implies stability of the scheme. +### Strang-Marchuk Splitting + +A natural way to improve the accuracy of operator splitting is to symmetrize the +scheme. The Strang-Marchuk splitting [Str:1968:ccd,Mar:1971:tsm](@cite) achieves +second-order accuracy for $N$ operators by performing a palindromic sweep + +```math +F_1(\Delta t/2) \to \cdots \to F_{N-1}(\Delta t/2) \to F_N(\Delta t) \to F_{N-1}(\Delta t/2) \to \cdots \to F_1(\Delta t/2) +``` + +More formally, for the simplest case of two operators $F_1$ and $F_2$ + +```math +\begin{aligned} + \text{Solve} \quad d_t u^1(t) &= F_1(u^1(t), p, t) & & \quad \text{on} \; [t_0, t_0 + \Delta t/2] \; \text{with} \; u^1(t_0) = u_0 \\ + \text{Solve} \quad d_t u^2(t) &= F_2(u^2(t), p, t) & & \quad \text{on} \; [t_0, t_0 + \Delta t] \; \text{with} \; u^2(t_0) = u^1(t_0 + \Delta t/2) \\ + \text{Solve} \quad d_t u^3(t) &= F_1(u^3(t), p, t) & & \quad \text{on} \; [t_0 + \Delta t/2, t_0 + \Delta t] \; \text{with} \; u^3(t_0 + \Delta t/2) = u^2(t_0 + \Delta t) +\end{aligned} +``` + +yielding $u(t_0 + \Delta t) \approx u^3(t_0 + \Delta t)$. + +### Analysis of Strang-Marchuk + +We show the second-order accuracy for two bounded linear operators $L_1$ and +$L_2$. The Strang-Marchuk approximation reads + +```math +\tilde{u}(t) = e^{L_1 t/2} \, e^{L_2 t} \, e^{L_1 t/2} \, u_0 \, . +``` + +Expanding the exponentials: + +```math +\begin{aligned} +e^{L_1 t/2} \, e^{L_2 t} \, e^{L_1 t/2} +&= \bigl(I + \tfrac{t}{2}L_1 + \tfrac{t^2}{8}L_1^2 + \cdots\bigr) + \bigl(I + t L_2 + \tfrac{t^2}{2}L_2^2 + \cdots\bigr) + \bigl(I + \tfrac{t}{2}L_1 + \tfrac{t^2}{8}L_1^2 + \cdots\bigr) \\ +&= I + t(L_1 + L_2) + \tfrac{t^2}{2}(L_1 + L_2)^2 + O(t^3) +\end{aligned} +``` + +which matches the Taylor expansion of $e^{(L_1+L_2)t}$ through the $t^2$ term. +The symmetry of the scheme causes the first-order commutator term +$[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. + ## References ```@bibliography diff --git a/docs/src/usage/index.md b/docs/src/usage/index.md index 9dc456c..4e3d0e4 100644 --- a/docs/src/usage/index.md +++ b/docs/src/usage/index.md @@ -49,3 +49,18 @@ for (u, t) in TimeChoiceIterator(integrator, 0.0:0.5:1.0) @show t, u end ``` + +For second-order accuracy, use the `StrangMarchuk` algorithm instead. +It performs the symmetric palindromic splitting +A₁(Δt/2) → … → Aₙ(Δt) → … → A₁(Δt/2): + +```julia +alg = StrangMarchuk( + (Euler(), Euler()) +) + +integrator = init(prob, alg, dt = 0.1) +for (u, t) in TimeChoiceIterator(integrator, 0.0:0.5:1.0) + @show t, u +end +``` diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index d752233..8c69eb5 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -35,7 +35,7 @@ include("integrator.jl") include("solver.jl") include("utils.jl") -export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov +export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, StrangMarchuk include("precompilation.jl") diff --git a/src/integrator.jl b/src/integrator.jl index 01c5e23..2f9a288 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -605,7 +605,13 @@ function step_footer!(integrator::AnySplitIntegrator) integrator.last_step_failed = false integrator.tprev = integrator.t integrator.t = fixed_t_for_floatingpoint_error!(integrator, ttmp) + # Children that step with subdivided dt (e.g. StrangMarchuk's `dt/2` + # 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) step_accept_controller!(integrator) + validate_time_point(integrator) elseif integrator.force_stepfail if isadaptive(integrator) step_reject_controller!(integrator) @@ -617,7 +623,6 @@ function step_footer!(integrator::AnySplitIntegrator) end integrator.last_step_failed = true end - validate_time_point(integrator) return nothing end @@ -907,26 +912,12 @@ function advance_solution_by!( dt ) SciMLBase.step!(sub, dt, true) - - # Unrecoverable failure: error immediately regardless of adaptive/non-adaptive - if !SciMLBase.successful_retcode(sub.status.retcode) && - sub.status.retcode != ReturnCode.Default - error("Inner integrator failed unrecoverably with retcode \ - $(sub.status.retcode) at t=$(child.t). Aborting.") - end return nothing end # Leaf disptach function advance_solution_by!(outer::AnySplitIntegrator, child::DEIntegrator, dt) SciMLBase.step!(child, dt, true) - - # Unrecoverable failure: error immediately regardless of adaptive/non-adaptive - if !SciMLBase.successful_retcode(child.sol.retcode) && - child.sol.retcode != ReturnCode.Default - error("Inner integrator failed unrecoverably with retcode \ - $(child.sol.retcode) at t=$(child.t). Aborting.") - end return nothing end diff --git a/src/precompilation.jl b/src/precompilation.jl index f638e74..df96482 100644 --- a/src/precompilation.jl +++ b/src/precompilation.jl @@ -34,10 +34,18 @@ end fsplit = GenericSplitFunction((f1, fsplitinner), (f1dofs, [1, 2, 3])) prob = OperatorSplittingProblem(fsplit, u0, tspan) - tstepper = LieTrotterGodunov((Euler(), LieTrotterGodunov((Euler(), Euler())))) - # Precompile init and a few steps - integrator = DiffEqBase.init(prob, tstepper, dt = 0.01, verbose = false) + # Precompile LieTrotterGodunov + tstepper_ltg = LieTrotterGodunov((Euler(), LieTrotterGodunov((Euler(), Euler())))) + integrator = DiffEqBase.init(prob, tstepper_ltg, dt = 0.01, verbose = false) step!(integrator) solve!(integrator) + + # Precompile StrangMarchuk + fsplit_sm = GenericSplitFunction((f1, f2), (f1dofs, f2dofs)) + prob_sm = OperatorSplittingProblem(fsplit_sm, u0, tspan) + tstepper_sm = StrangMarchuk((Euler(), Euler())) + integrator_sm = DiffEqBase.init(prob_sm, tstepper_sm, dt = 0.01, verbose = false) + step!(integrator_sm) + solve!(integrator_sm) end diff --git a/src/solver.jl b/src/solver.jl index 3163474..a72c62f 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -56,3 +56,104 @@ end 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 + + 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 + + backward_sync_subintegrator!(parent, child, idxs, sync) + end +end + +function _perform_step!( + parent, + children::Tuple, + cache::StrangMarchukCache, + dt + ) + half_dt = dt / 2 + + _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/test/operator_splitting_api.jl b/test/operator_splitting_api.jl index 2f11a30..7cabb90 100644 --- a/test/operator_splitting_api.jl +++ b/test/operator_splitting_api.jl @@ -126,6 +126,7 @@ end end FakeAdaptiveLTG(inner) = FakeAdaptiveAlgorithm(LieTrotterGodunov(inner)) +FakeAdaptiveSM(inner) = FakeAdaptiveAlgorithm(StrangMarchuk(inner)) function Base.show(io::IO, alg::FakeAdaptiveAlgorithm) print(io, "FAKE (") @@ -133,6 +134,11 @@ function Base.show(io::IO, alg::FakeAdaptiveAlgorithm) return print(io, ")") end +# StrangMarchuk steps child 1 twice per outer step (two half-steps). +_sub1_iter_factor(::LieTrotterGodunov) = 1 +_sub1_iter_factor(::StrangMarchuk) = 2 +_sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg) + # --------------------------------------------------------------------------- # Tests @@ -161,7 +167,7 @@ end nsteps = ceil(Int, (tspan[2] - tspan[1]) / dt) - for TimeStepperType in (LieTrotterGodunov, FakeAdaptiveLTG) + for TimeStepperType in (LieTrotterGodunov, FakeAdaptiveLTG, StrangMarchuk, FakeAdaptiveSM) @testset "$tstepper" for (prob, tstepper) in ( (prob1a, TimeStepperType((Euler(), Euler()))), (prob1a, TimeStepperType((Tsit5(), Euler()))), @@ -185,6 +191,7 @@ end sub1 = integrator.child_subintegrators[1] sub2 = integrator.child_subintegrators[2] + expected_sub1_iters = _sub1_iter_factor(tstepper) * nsteps DiffEqBase.solve!(integrator) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success @@ -195,7 +202,7 @@ end @test integrator.iter == nsteps @test sub1.t ≈ tspan[2] - @test sub1.iter == nsteps + @test sub1.iter == expected_sub1_iters @test sub2.t ≈ tspan[2] @test sub2.iter == nsteps @@ -227,14 +234,14 @@ end @test integrator.iter == nsteps @test sub1.t ≈ tspan[2] - @test sub1.iter == nsteps + @test sub1.iter == expected_sub1_iters @test sub2.t ≈ tspan[2] @test sub2.iter == nsteps end end - for TimeStepperType in (FakeAdaptiveLTG,) + for TimeStepperType in (FakeAdaptiveLTG, FakeAdaptiveSM) @testset "Adaptive solver type $TimeStepperType | $tstepper" for (prob, tstepper) in ( (prob1a, TimeStepperType((Tsit5(), Tsit5()))), (prob2, TimeStepperType((Tsit5(), TimeStepperType((Tsit5(), Tsit5()))))), @@ -281,6 +288,84 @@ end end end + @testset "StrangMarchuk with 3 operators" begin + dt = 0.01π + # f1 + f3 + f3 = f1 + f2, so the reference solution is the same trueu. + f1dofs = [1, 2, 3] + f3dofs = [1, 3] + fsplit3 = GenericSplitFunction((f1, f3, f3), (f1dofs, f3dofs, f3dofs)) + prob3 = OperatorSplittingProblem(fsplit3, u0, tspan) + nsteps = ceil(Int, (tspan[2] - tspan[1]) / dt) + + @testset "$tstepper" for tstepper in ( + StrangMarchuk((Euler(), Euler(), Euler())), + StrangMarchuk((Tsit5(), Euler(), Tsit5())), + StrangMarchuk((Tsit5(), Tsit5(), Tsit5())), + ) + integrator = DiffEqBase.init( + prob3, tstepper, dt = dt, verbose = true, alias_u0 = false, adaptive = false + ) + DiffEqBase.solve!(integrator) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success + @test isapprox(integrator.u, trueu, atol = 1.0e-6) + @test integrator.t ≈ tspan[2] + @test integrator.iter == nsteps + + sub1 = integrator.child_subintegrators[1] + sub2 = integrator.child_subintegrators[2] + sub3 = integrator.child_subintegrators[3] + # Palindromic: children 1 & 2 get two half-steps, child 3 gets one full step + @test sub1.iter == 2 * nsteps + @test sub2.iter == 2 * nsteps + @test sub3.iter == nsteps + end + end + + @testset "Convergence order" begin + # Use non-commuting operators so splitting error is non-zero. + # A = diag(-1,-2), B = [0 0.5; 0.5 0] have [A,B] ≠ 0. + function ode_conv_A(du, u, p, t) + du[1] = -u[1] + return du[2] = -2 * u[2] + end + function ode_conv_B(du, u, p, t) + du[1] = 0.5 * u[2] + return du[2] = 0.5 * u[1] + end + fA = ODEFunction(ode_conv_A) + fB = ODEFunction(ode_conv_B) + + conv_tspan = (0.0, 1.0) + conv_u0 = [1.0, 1.0] + conv_trueu = exp(conv_tspan[2] * [-1.0 0.5; 0.5 -2.0]) * conv_u0 + + conv_dofs = [1, 2] + fsplit_conv = GenericSplitFunction((fA, fB), (conv_dofs, conv_dofs)) + prob_conv = OperatorSplittingProblem(fsplit_conv, conv_u0, conv_tspan) + + dts = [0.1, 0.05, 0.025] + for (TimeStepperType, expected_order) in ( + (LieTrotterGodunov, 1), + (StrangMarchuk, 2), + ) + @testset "$TimeStepperType (order $expected_order)" begin + errors = map(dts) do dt_i + tstepper = TimeStepperType((Tsit5(), Tsit5())) + integrator = DiffEqBase.init( + prob_conv, tstepper, dt = dt_i, verbose = false, + alias_u0 = false, adaptive = false + ) + DiffEqBase.solve!(integrator) + maximum(abs, integrator.u .- conv_trueu) + end + for i in 1:(length(errors) - 1) + rate = log2(errors[i] / errors[i + 1]) + @test rate ≈ expected_order atol = 0.3 + end + end + end + end + @testset "Instability detection" begin dt = 0.01π @@ -296,7 +381,7 @@ end fsplit_NaN = GenericSplitFunction((f1, f_NaN), (f1dofs, f3dofs)) prob_NaN = OperatorSplittingProblem(fsplit_NaN, u0, tspan) - for TimeStepperType in (LieTrotterGodunov,) + for TimeStepperType in (LieTrotterGodunov, StrangMarchuk) @testset "Solver type $TimeStepperType | $tstepper" for tstepper in ( TimeStepperType((Euler(), Euler())), TimeStepperType((Tsit5(), Euler())),