Add StrangMarchuk second-order symmetric operator splitting (#73)

* Add generic Strang-Marchuk second-order symmetric operator splitting
* Add convergence order tests
* Remove stringent errors on subintegration failure

Author: oameye

docs/src/api-reference/index.md

@@ -12,8 +12,9 @@ OperatorSplittingProblem
 GenericSplitFunction
 ```
 
-## Solver
+## Solvers
 
 ```@docs
 LieTrotterGodunov
+StrangMarchuk
 ```

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}
+}

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

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
+```

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")
 

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
 

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

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