From eaa54ea1558481bebb0b06697bf9df32fe0f1af8 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas Date: Mon, 29 Dec 2025 07:34:56 -0500 Subject: [PATCH] Add StrangMarchuk second-order symmetric operator splitting MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Implements the classical Strang-Marchuk ABA splitting scheme which achieves second-order accuracy through symmetry. For two operators A and B, the scheme performs: - A(dt/2) -> B(dt) -> A(dt/2) This is an improvement over the first-order Lie-Trotter-Godunov scheme (A(dt) -> B(dt)) and is commonly used for semidiscrete exponential propagator (EP) problems. Closes #10 References: * G. Strang, On the construction and comparison of difference schemes, SIAM Journal on Numerical Analysis, 5 (1968), pp. 506-517 * G. I. Marchuk, On the theory of the splitting-up method, in Numerical Solution of Partial Differential Equations-II, Academic Press, 1971, pp. 469-500 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 --- src/OrdinaryDiffEqOperatorSplitting.jl | 2 +- src/solver.jl | 103 +++++++++++++++++++++++++ test/operator_splitting_api.jl | 4 +- 3 files changed, 106 insertions(+), 3 deletions(-) diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index 94f3152..8d0ce8f 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -24,6 +24,6 @@ include("integrator.jl") include("solver.jl") include("utils.jl") -export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov +export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, StrangMarchuk end diff --git a/src/solver.jl b/src/solver.jl index 16837ca..b285c69 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -51,3 +51,106 @@ end backward_sync_subintegrator!(outer_integrator, subinteg, idxs, synchronizer) end end + +# Strang-Marchuk Splitting Implementation +""" + StrangMarchuk <: AbstractOperatorSplittingAlgorithm + +A second order symmetric operator splitting algorithm attributed to [Str:1968:ccd,Mar:1971:tsm](@cite). +This implements the classical ABA scheme where for operators A and B: +- First half-step with A (dt/2) +- Full step with B (dt) +- Second half-step with A (dt/2) + +For two operators, this gives second-order accuracy in time. + +# References +* G. Strang, On the construction and comparison of difference schemes, SIAM Journal on +Numerical Analysis, 5 (1968), pp. 506–517 +* G. I. Marchuk, On the theory of the splitting-up method, in Numerical Solution of Partial +Differential Equations-II, Academic Press, 1971, pp. 469 – 500 +""" +struct StrangMarchuk{AlgTupleType} <: AbstractOperatorSplittingAlgorithm + inner_algs::AlgTupleType # Tuple of timesteppers for inner problems +end + +struct StrangMarchukCache{uType, uprevType, iiType} <: AbstractOperatorSplittingCache + u::uType + uprev::uprevType + inner_caches::iiType +end + +function init_cache(f::GenericSplitFunction, alg::StrangMarchuk; + uprev::AbstractArray, u::AbstractVector, + inner_caches, + alias_uprev = true, + alias_u = false +) + _uprev = alias_uprev ? uprev : SciMLBase.recursivecopy(uprev) + _u = alias_u ? u : SciMLBase.recursivecopy(u) + StrangMarchukCache(_u, _uprev, inner_caches) +end + +@inline function advance_solution_to!( + outer_integrator::OperatorSplittingIntegrator, + subintegrators::Tuple, solution_indices::Tuple, + synchronizers::Tuple, cache::StrangMarchukCache, tnext) + # Strang-Marchuk ABA splitting scheme + # For two operators A and B: A(dt/2) -> B(dt) -> A(dt/2) + # This achieves second-order accuracy through symmetry + + @unpack inner_caches = cache + tcurr = outer_integrator.t + dt = tnext - tcurr + thalf = tcurr + dt / 2 + + # We require exactly 2 subintegrators for the classical ABA scheme + if length(subintegrators) != 2 + error("StrangMarchuk splitting requires exactly 2 operators") + end + + # First operator A at half step + subinteg_A = subintegrators[1] + synchronizer_A = synchronizers[1] + idxs_A = solution_indices[1] + cache_A = inner_caches[1] + + # Second operator B at full step + subinteg_B = subintegrators[2] + synchronizer_B = synchronizers[2] + idxs_B = solution_indices[2] + cache_B = inner_caches[2] + + # Step 1: A for dt/2 + @timeit_debug "sync ->" forward_sync_subintegrator!(outer_integrator, subinteg_A, idxs_A, synchronizer_A) + @timeit_debug "time solve A (half)" advance_solution_to!( + outer_integrator, subinteg_A, idxs_A, synchronizer_A, cache_A, thalf) + if !(subinteg_A isa Tuple) && + subinteg_A.sol.retcode ∉ + (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success) + return + end + backward_sync_subintegrator!(outer_integrator, subinteg_A, idxs_A, synchronizer_A) + + # Step 2: B for dt + @timeit_debug "sync ->" forward_sync_subintegrator!(outer_integrator, subinteg_B, idxs_B, synchronizer_B) + @timeit_debug "time solve B (full)" advance_solution_to!( + outer_integrator, subinteg_B, idxs_B, synchronizer_B, cache_B, tnext) + if !(subinteg_B isa Tuple) && + subinteg_B.sol.retcode ∉ + (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success) + return + end + backward_sync_subintegrator!(outer_integrator, subinteg_B, idxs_B, synchronizer_B) + + # Step 3: A for dt/2 + @timeit_debug "sync ->" forward_sync_subintegrator!(outer_integrator, subinteg_A, idxs_A, synchronizer_A) + @timeit_debug "time solve A (half)" advance_solution_to!( + outer_integrator, subinteg_A, idxs_A, synchronizer_A, cache_A, tnext) + if !(subinteg_A isa Tuple) && + subinteg_A.sol.retcode ∉ + (SciMLBase.ReturnCode.Default, SciMLBase.ReturnCode.Success) + return + end + backward_sync_subintegrator!(outer_integrator, subinteg_A, idxs_A, synchronizer_A) +end diff --git a/test/operator_splitting_api.jl b/test/operator_splitting_api.jl index 29e7030..5a502a1 100644 --- a/test/operator_splitting_api.jl +++ b/test/operator_splitting_api.jl @@ -68,7 +68,7 @@ f2 = ODEFunction(ode2) fsplit2_outer = GenericSplitFunction((f1, fsplit2_inner), (f1dofs, f2dofs)) prob2 = OperatorSplittingProblem(fsplit2_outer, u0, tspan) - for TimeStepperType in (LieTrotterGodunov,) + for TimeStepperType in (LieTrotterGodunov, StrangMarchuk) @testset "Solver type $TimeStepperType | $tstepper" for (prob, tstepper) in ( (prob1, TimeStepperType((Euler(), Euler()))), (prob1, TimeStepperType((Tsit5(), Euler()))), @@ -137,7 +137,7 @@ f2 = ODEFunction(ode2) fsplit_NaN = GenericSplitFunction((f1, f_NaN), (f1dofs, f_NaN_dofs)) prob_NaN = OperatorSplittingProblem(fsplit_NaN, u0, tspan) - for TimeStepperType in (LieTrotterGodunov,) + for TimeStepperType in (LieTrotterGodunov, StrangMarchuk) @testset "Solver type $TimeStepperType | $tstepper" for (prob, tstepper) in ( (prob1, TimeStepperType((Euler(), Euler()))), (prob1, TimeStepperType((Tsit5(), Euler()))),