From 667a6647bf536e794765fb0b8f90620afe09d2e3 Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Sun, 10 May 2026 17:52:35 +0200 Subject: [PATCH 1/2] Move toward reusing OrdinaryDiffEqCore infrastructure (refs #23) Replace custom IntegratorStats with SciMLBase.DEStats and move tstops/saveat/callback/advance_to_tstop off the top-level integrator struct into a DEOptions (OrdinaryDiffEqCore v4) via a new build_split_deoptions helper. This aligns OperatorSplittingIntegrator more closely with the broader OrdinaryDiffEq integrator interface. SplitSubIntegrator retains its own lightweight IntegratorOptions. --- Project.toml | 15 +- src/OrdinaryDiffEqOperatorSplitting.jl | 18 +- src/integrator.jl | 219 ++++++++++++++++--------- test/consistency.jl | 16 +- test/operator_splitting_api.jl | 91 ++++++++-- 5 files changed, 249 insertions(+), 110 deletions(-) diff --git a/Project.toml b/Project.toml index 8e0295d..892a798 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ authors = ["Dennis Ogiermann and contr [deps] CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +Reuse OrdinaryDiffEqCore infrastructure = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" OrdinaryDiffEqLowOrderRK = "1344f307-1e59-4825-a18e-ace9aa3fa4c6" @@ -18,16 +19,16 @@ Unrolled = "9602ed7d-8fef-5bc8-8597-8f21381861e8" [compat] CommonSolve = "0.2.4" DataStructures = "0.18.22, 0.19" -DiffEqBase = "6.165.1, 7" +DiffEqBase = "7" ExplicitImports = "1" -ModelingToolkit = "10" -OrdinaryDiffEqCore = "1.19.0, 2, 3.1, 4" -OrdinaryDiffEqLowOrderRK = "1.7, 2" -OrdinaryDiffEqTsit5 = "1.1.0, 2" +ModelingToolkit = "10.31.2, 11" +OrdinaryDiffEqCore = "4" +OrdinaryDiffEqLowOrderRK = "2" +OrdinaryDiffEqTsit5 = "2" PrecompileTools = "1.0" -RecursiveArrayTools = "3.39.0, 4" +RecursiveArrayTools = "4" SafeTestsets = "0.1.0" -SciMLBase = "2.77.0, 3.1" +SciMLBase = "3.1" TimerOutputs = "0.5.28" Unrolled = "0.1.5" julia = "1.10" diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index d752233..d7a534d 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -6,21 +6,19 @@ timeit_debug_enabled() = false import Unrolled: @unroll import SciMLBase, DiffEqBase, DataStructures -import SciMLBase: ReturnCode +import SciMLBase: ReturnCode, DEStats import SciMLBase: DEIntegrator, NullParameters, isadaptive import RecursiveArrayTools +import LinearAlgebra import OrdinaryDiffEqCore: OrdinaryDiffEqCore, isdtchangeable, - stepsize_controller!, step_accept_controller!, step_reject_controller! - -# In OrdinaryDiffEq v7 / DiffEqBase v7, passing verbose::Bool to inner ODE -# integrators is no longer supported. Convert Bool → DEVerbosity when available. -@static if isdefined(DiffEqBase, :DEVerbosity) - _inner_verbose(verbose::Bool) = verbose ? DiffEqBase.DEFAULT_VERBOSE : DiffEqBase.DEVerbosity(DiffEqBase.None()) -else - _inner_verbose(verbose::Bool) = verbose -end + stepsize_controller!, step_accept_controller!, step_reject_controller!, + DEOptions + +# DiffEqBase v7 no longer accepts verbose::Bool for inner ODE integrators; convert to DEVerbosity. +_inner_verbose(verbose::Bool) = verbose ? DiffEqBase.DEFAULT_VERBOSE : DiffEqBase.DEVerbosity(DiffEqBase.None()) +_inner_verbose(verbose::DiffEqBase.DEVerbosity) = verbose abstract type AbstractOperatorSplitFunction <: SciMLBase.AbstractODEFunction{true} end abstract type AbstractOperatorSplittingAlgorithm end diff --git a/src/integrator.jl b/src/integrator.jl index 01c5e23..5c5a8c0 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -1,11 +1,3 @@ -mutable struct IntegratorStats - naccept::Int64 - nreject::Int64 - # TODO inner solver stats -end - -IntegratorStats() = IntegratorStats(0, 0) - Base.@kwdef mutable struct IntegratorOptions{tType, fType, F3} adaptive::Bool dtmin::tType = eps(Float64) @@ -94,7 +86,7 @@ mutable struct SplitSubIntegrator{ last_step_failed::Bool u_modified::Bool # TODO we can probably remove this status::SplitSubIntegratorStatus - stats::IntegratorStats + stats::DEStats cache::cacheType child_subintegrators::childSubintType # Tuple solution_indices::solidxType @@ -132,17 +124,13 @@ mutable struct OperatorSplittingIntegrator{ uType, tType, pType, - heapType, - tstopsType, - saveatType, - callbackType, cacheType, solType, subintTreeType, childSolidxType, childSyncType, controllerType, - optionsType, + optsType, } <: SciMLBase.AbstractODEIntegrator{algType, true, uType, tType} const f::fType const alg::algType @@ -155,31 +143,82 @@ mutable struct OperatorSplittingIntegrator{ dt::tType # Time step length used during time marching dtcache::tType # Proposed time step length const dtchangeable::Bool - tstops::heapType - _tstops::tstopsType - saveat::heapType - _saveat::saveatType - callback::callbackType - advance_to_tstop::Bool last_step_failed::Bool force_stepfail::Bool isout::Bool u_modified::Bool + just_hit_tstop::Bool cache::cacheType sol::solType - # Tuple of SplitSubIntegrator nodes (one per top-level operator). - child_subintegrators::subintTreeType + child_subintegrators::subintTreeType # Tuple of SplitSubIntegrator nodes (one per top-level operator) child_solution_indices::childSolidxType # Tuple child_synchronizers::childSyncType # Tuple iter::Int controller::controllerType - opts::optionsType - stats::IntegratorStats + opts::optsType # DEOptions + stats::DEStats tdir::tType end const AnySplitIntegrator = Union{SplitSubIntegrator, OperatorSplittingIntegrator} +# --------------------------------------------------------------------------- +# build_split_deoptions +# --------------------------------------------------------------------------- +function build_split_deoptions(tType; + tstops, + saveat, + d_discontinuities, + tstops_cache, + saveat_cache, + d_discontinuities_cache, + callback, + adaptive = false, + dtmin = tType(eps(Float64)), + dtmax = tType(Inf), + failfactor = tType(4), + verbose = true, + isoutofdomain = DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN, + advance_to_tstop = false, + save_everystep = false, + save_on = true, + save_start = true, + save_end = true, + ) + QT = tType + return DEOptions( + 1_000_000, # maxiters + save_everystep, + adaptive, + nothing, # abstol + nothing, # reltol + QT(failfactor), + tType(dtmax), + tType(dtmin), + DiffEqBase.ODE_DEFAULT_NORM, + LinearAlgebra.opnorm, + nothing, # save_idxs + tstops, saveat, d_discontinuities, + tstops_cache, saveat_cache, d_discontinuities_cache, + nothing, # userdata + false, 0, "ODE", # progress, progress_steps, progress_name + DiffEqBase.ODE_DEFAULT_PROG_MESSAGE, + :ode, # progress_id + true, false, # timeseries_errors, dense_errors + nothing, false, # delta, dense + save_on, save_start, save_end, + false, false, # save_noise, save_discretes + nothing, # save_end_user + callback, + isoutofdomain, + DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK, + verbose, + false, false, # calck, force_dtmin + advance_to_tstop, + false, # stop_at_next_tstop + ) +end + # --------------------------------------------------------------------------- # __init # --------------------------------------------------------------------------- @@ -192,11 +231,13 @@ function SciMLBase.__init( saveat = (), d_discontinuities = (), save_everystep = false, + save_on = true, + save_start = true, + save_end = true, callback = nothing, advance_to_tstop = false, adaptive = isadaptive(alg), - controller = nothing, - # controller = OrdinaryDiffEqCore.PIController(0.14, 0.08), + controller = nothing, # e.g. OrdinaryDiffEqCore.PIController(0.14, 0.08) alias_u0 = false, verbose = true, kwargs... @@ -209,15 +250,13 @@ function SciMLBase.__init( dt = tf > t0 ? dt : -dt tType = typeof(dt) - (!isadaptive(alg) && adaptive && verbose) && + (!isadaptive(alg) && adaptive && (verbose isa Bool ? verbose : true)) && @warn("The algorithm $alg is not adaptive.") dtchangeable = isdtchangeable(alg) - if tstops isa AbstractArray || tstops isa Tuple || tstops isa Number - _tstops = nothing - else - _tstops = tstops + tstops_cache = tstops + if !(tstops isa AbstractArray || tstops isa Tuple || tstops isa Number) tstops = () end @@ -255,28 +294,43 @@ function SciMLBase.__init( child_solution_indices = ntuple(i -> prob.f.solution_indices[i], length(prob.f.functions)) child_synchronizers = ntuple(i -> prob.f.synchronizers[i], length(prob.f.functions)) + tdir_val = tType(tstops_internal.ordering isa DataStructures.FasterForward ? 1 : -1) + + opts = build_split_deoptions(tType; + tstops = tstops_internal, + saveat = saveat_internal, + d_discontinuities = d_discontinuities_internal, + tstops_cache, + saveat_cache = saveat, + d_discontinuities_cache = d_discontinuities, + callback, + adaptive, + verbose, + advance_to_tstop, + save_everystep, + save_on, + save_start, + save_end, + ) + integrator = OperatorSplittingIntegrator( prob.f, alg, u, uprev, tmp, p, - t0, copy(dt), + t0, t0, dt, dtcache, dtchangeable, - tstops_internal, tstops, - saveat_internal, saveat, - callback, - advance_to_tstop, - false, false, false, false, + false, false, false, false, false, # last_step_failed, force_stepfail, isout, u_modified, just_hit_tstop cache, sol, child_subintegrators, child_solution_indices, child_synchronizers, 0, controller, - IntegratorOptions(; verbose, adaptive), - IntegratorStats(), - tType(tstops_internal.ordering isa DataStructures.FasterForward ? 1 : -1) + opts, + DEStats(0), + tdir_val, ) DiffEqBase.initialize!(callback, u0, t0, integrator) return integrator @@ -294,8 +348,8 @@ function DiffEqBase.reinit!( tf = integrator.sol.prob.tspan[2], dt = isadaptive(integrator) ? nothing : integrator.dtcache, erase_sol = false, - tstops = integrator._tstops, - saveat = integrator._saveat, + tstops = integrator.opts.tstops_cache, + saveat = integrator.opts.saveat_cache, reinit_callbacks = true, reinit_retcode = true ) @@ -305,8 +359,9 @@ function DiffEqBase.reinit!( integrator.tprev = t0 if dt !== nothing integrator.dt = dt + integrator.dtcache = dt end - integrator.tstops, integrator.saveat = + integrator.opts.tstops, integrator.opts.saveat = tstops_and_saveat_heaps(t0, tf, tstops, saveat) integrator.iter = 0 if erase_sol @@ -314,9 +369,9 @@ function DiffEqBase.reinit!( resize!(integrator.sol.u, 0) end if reinit_callbacks - DiffEqBase.initialize!(integrator.callback, u0, t0, integrator) + DiffEqBase.initialize!(integrator.opts.callback, u0, t0, integrator) else - saving_callback = integrator.callback.discrete_callbacks[end] + saving_callback = integrator.opts.callback.discrete_callbacks[end] DiffEqBase.initialize!(saving_callback, u0, t0, integrator) end if reinit_retcode @@ -431,7 +486,11 @@ function OrdinaryDiffEqCore.handle_tstop!(integrator::AnySplitIntegrator) return nothing end -notify_integrator_hit_tstop!(integrator::AnySplitIntegrator) = nothing +notify_integrator_hit_tstop!(integrator::SplitSubIntegrator) = nothing +function notify_integrator_hit_tstop!(integrator::OperatorSplittingIntegrator) + integrator.just_hit_tstop = true + return nothing +end # --------------------------------------------------------------------------- @@ -495,7 +554,7 @@ end function rollback_child!(child::SplitSubIntegrator, u_master) child.u .= @view u_master[child.solution_indices] RecursiveArrayTools.recursivecopy!(child.uprev, child.u) - _rollback_children!(child.child_subintegrators, u_master) + rollback_children!(child.child_subintegrators, u_master) return nothing end function rollback_child!(child::DEIntegrator, u_master) @@ -545,8 +604,9 @@ end is_first_iteration(integrator::AnySplitIntegrator) = integrator.iter == 0 increment_iteration(integrator::AnySplitIntegrator) = integrator.iter += 1 -function footer_reset_flags!(integrator) +function footer_reset_flags!(integrator::OperatorSplittingIntegrator) integrator.u_modified = false + integrator.just_hit_tstop = false return end footer_reset_flags!(::SplitSubIntegrator) = nothing @@ -555,15 +615,6 @@ function setup_validity_flags!(integrator, t_next) return end setup_validity_flags!(::SplitSubIntegrator, _) = nothing -function fix_solution_buffer_sizes!(integrator, sol) - resize!(integrator.sol.t, integrator.saveiter) - resize!(integrator.sol.u, integrator.saveiter) - if !(integrator.sol isa SciMLBase.DAESolution) - resize!(integrator.sol.k, integrator.saveiter_dense) - end - return -end - function fixed_t_for_floatingpoint_error!(integrator::AnySplitIntegrator, ttmp) return if DiffEqBase.has_tstop(integrator) tstop = integrator.tdir * DiffEqBase.first_tstop(integrator) @@ -633,7 +684,7 @@ function SciMLBase.__solve( end function DiffEqBase.solve!(integrator::OperatorSplittingIntegrator) - while !isempty(integrator.tstops) + while !isempty(integrator.opts.tstops) while tdir(integrator) * integrator.t < SciMLBase.first_tstop(integrator) step_header!(integrator) @timeit_debug "check_error" SciMLBase.check_error!(integrator) ∉ ( @@ -652,8 +703,8 @@ function DiffEqBase.solve!(integrator::OperatorSplittingIntegrator) ) end -function DiffEqBase.step!(integrator::AnySplitIntegrator) - @timeit_debug "step!" if integrator.advance_to_tstop +function DiffEqBase.step!(integrator::OperatorSplittingIntegrator) + @timeit_debug "step!" if integrator.opts.advance_to_tstop tstop = SciMLBase.first_tstop(integrator) while !reached_tstop(integrator, tstop) step_header!(integrator) @@ -750,7 +801,7 @@ end end _child_retcode(child::DEIntegrator) = SciMLBase.check_error(child) -_child_retcode(child::SplitSubIntegrator) = child.status.retcode +_child_retcode(child::SplitSubIntegrator) = SciMLBase.check_error(child) function setup_u(prob::OperatorSplittingProblem, solver, alias_u0) return alias_u0 ? prob.u0 : RecursiveArrayTools.recursivecopy(prob.u0) @@ -803,8 +854,7 @@ end # Time helpers -tdir(integrator) = - integrator.tstops.ordering isa DataStructures.FasterForward ? 1 : -1 +tdir(integrator::OperatorSplittingIntegrator) = integrator.tdir is_past_t(integrator, t) = tdir(integrator) * (t - integrator.t) ≤ zero(integrator.t) function reached_tstop(integrator, tstop, stop_at_tstop = integrator.dtchangeable) @@ -820,7 +870,7 @@ end # SciMLBase integrator interface function SciMLBase.done(integrator::OperatorSplittingIntegrator) integrator.sol.retcode ∉ (ReturnCode.Default, ReturnCode.Success) && return true - if isempty(integrator.tstops) + if isempty(integrator.opts.tstops) SciMLBase.postamble!(integrator) return true end @@ -828,7 +878,7 @@ function SciMLBase.done(integrator::OperatorSplittingIntegrator) end function SciMLBase.postamble!(integrator::OperatorSplittingIntegrator) - return DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator) + return DiffEqBase.finalize!(integrator.opts.callback, integrator.u, integrator.t, integrator) end function __step!(integrator::AnySplitIntegrator) @@ -912,7 +962,7 @@ function advance_solution_by!( 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.") + $(sub.status.retcode) at t=$(sub.t). Aborting.") end return nothing end @@ -1033,13 +1083,13 @@ function _build_child( controller, false, false, false, # force_stepfail, last_step_failed, u_modified SplitSubIntegratorStatus(), - IntegratorStats(), + DEStats(0), level_cache, child_subintegrators, solution_indices, child_solution_indices, child_synchronizers, - IntegratorOptions(; verbose, adaptive), + IntegratorOptions(; verbose = verbose isa Bool ? verbose : true, adaptive), one(tType), ) @@ -1062,10 +1112,16 @@ function _build_child( controller = nothing ) where {S, T, P, F} u = uouter[solution_indices] + # MTK v11 compiled systems require a symbolic map for u0; plain SciMLFunctions accept arrays. + u0 = if f isa SciMLBase.AbstractSciMLFunction + u + else + SciMLBase.variable_symbols(f) .=> u + end prob2 = if p isa NullParameters - SciMLBase.ODEProblem(f, u, (t0, tf)) + SciMLBase.ODEProblem(f, u0, (t0, tf)) else - SciMLBase.ODEProblem(f, u, (t0, tf), p) + SciMLBase.ODEProblem(f, u0, (t0, tf), p) end integrator = SciMLBase.__init( @@ -1087,9 +1143,13 @@ end # --------------------------------------------------------------------------- SciMLBase.has_stats(::AnySplitIntegrator) = true -SciMLBase.has_tstop(i::AnySplitIntegrator) = !isempty(i.tstops) -SciMLBase.first_tstop(i::AnySplitIntegrator) = first(i.tstops) -SciMLBase.pop_tstop!(i::AnySplitIntegrator) = pop!(i.tstops) +SciMLBase.has_tstop(i::OperatorSplittingIntegrator) = !isempty(i.opts.tstops) +SciMLBase.first_tstop(i::OperatorSplittingIntegrator) = first(i.opts.tstops) +SciMLBase.pop_tstop!(i::OperatorSplittingIntegrator) = pop!(i.opts.tstops) + +SciMLBase.has_tstop(i::SplitSubIntegrator) = !isempty(i.tstops) +SciMLBase.first_tstop(i::SplitSubIntegrator) = first(i.tstops) +SciMLBase.pop_tstop!(i::SplitSubIntegrator) = pop!(i.tstops) DiffEqBase.get_dt(i::AnySplitIntegrator) = i.dt function set_dt!(i::DiffEqBase.DEIntegrator, dt) @@ -1097,7 +1157,16 @@ function set_dt!(i::DiffEqBase.DEIntegrator, dt) return i.dt = dt end -function DiffEqBase.add_tstop!(i::AnySplitIntegrator, t) +function DiffEqBase.add_tstop!(i::OperatorSplittingIntegrator, t) + is_past_t(i, 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.opts.tstops, t) + return nothing +end + +function DiffEqBase.add_tstop!(i::SplitSubIntegrator, t) is_past_t(i, t) && error("Cannot add a tstop at $t because that is behind the current \ integrator time $(i.t)") @@ -1110,7 +1179,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.opts.saveat, t) return nothing end diff --git a/test/consistency.jl b/test/consistency.jl index da60f10..e08729e 100644 --- a/test/consistency.jl +++ b/test/consistency.jl @@ -15,8 +15,18 @@ splitting_solver = LieTrotterGodunov((Euler(),)) integrator1 = init(prob1, splitting_solver; dt = dt) integrator2 = init(prob2, Euler(); dt = dt) -for ((u1, t1), (u2, t2)) in zip(TimeChoiceIterator(integrator1, tspan[1]:(2dt):tspan[2]), TimeChoiceIterator(integrator2, tspan[1]:(2dt):tspan[2])) - @test u1 ≈ u2 - @test t1 ≈ t2 + +# Compare solutions at every 2*dt time point. +# Original test used TimeChoiceIterator (removed in DiffEqBase v7). +nsteps = Int(round((tspan[2] - tspan[1]) / dt)) +@test integrator1.u ≈ integrator2.u +@test integrator1.t ≈ integrator2.t +for _ in 1:(nsteps ÷ 2) + step!(integrator1) + step!(integrator1) + step!(integrator2) + step!(integrator2) + @test integrator1.u ≈ integrator2.u + @test integrator1.t ≈ integrator2.t end @test integrator1.iter == integrator2.iter diff --git a/test/operator_splitting_api.jl b/test/operator_splitting_api.jl index 2f11a30..2f9fae3 100644 --- a/test/operator_splitting_api.jl +++ b/test/operator_splitting_api.jl @@ -54,17 +54,12 @@ f3 = ODEFunction(ode3) @independent_variables time Dt = Differential(time) -@mtkmodel TestModelODE2 begin - @variables begin - u1(time) - u2(time) - end - @equations begin - Dt(u1) ~ -0.01u2 - Dt(u2) ~ -0.01u1 - end -end -@named testmodel2 = TestModelODE2() +@variables u1(time) u2(time) +eqs = [ + Dt(u1) ~ -0.01u2, + Dt(u2) ~ -0.01u1, +] +@named testmodel2 = System(eqs, time) testsys2 = mtkcompile(testmodel2; sort_eqs = false) # Test whether adaptive code path works in principle @@ -202,7 +197,8 @@ end DiffEqBase.reinit!(integrator; dt = dt) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default - for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, tspan[1]:5.0:tspan[2]) + while !SciMLBase.done(integrator) + DiffEqBase.step!(integrator) end @test isapprox(ufinal, integrator.u, atol = 1.0e-12) @test integrator.t ≈ tspan[2] @@ -211,7 +207,8 @@ end DiffEqBase.reinit!(integrator; dt = dt) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default - for (uprev, tprev, u, t) in DiffEqBase.intervals(integrator) + while !SciMLBase.done(integrator) + DiffEqBase.step!(integrator) end @test isapprox(ufinal, integrator.u, atol = 1.0e-12) @test integrator.t ≈ tspan[2] @@ -255,7 +252,8 @@ end @test integrator.dt == dt @test integrator.dt == integrator.dtcache @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default - for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, tspan[1]:5.0:tspan[2]) + while !SciMLBase.done(integrator) + DiffEqBase.step!(integrator) end @test isapprox(ufinal, integrator.u, atol = 1.0e-12) @test integrator.t ≈ tspan[2] @@ -264,7 +262,8 @@ end DiffEqBase.reinit!(integrator; dt = dt) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default - for (uprev, tprev, u, t) in DiffEqBase.intervals(integrator) + while !SciMLBase.done(integrator) + DiffEqBase.step!(integrator) end @test isapprox(ufinal, integrator.u, atol = 1.0e-12) @test integrator.t ≈ tspan[2] @@ -313,4 +312,66 @@ end end end end + + @testset "reinit! resets dtcache when dt changes" begin + dt_coarse = 0.1 + dt_fine = dt_coarse / 2 + + fsplit_rc = GenericSplitFunction((f1, f2), ([1, 2, 3], [1, 3])) + prob_rc = OperatorSplittingProblem(fsplit_rc, u0, tspan) + tstepper_rc = LieTrotterGodunov((Euler(), Euler())) + + integrator_rc = DiffEqBase.init(prob_rc, tstepper_rc; dt = dt_coarse, alias_u0 = false) + DiffEqBase.solve!(integrator_rc) + @test integrator_rc.sol.retcode == DiffEqBase.ReturnCode.Success + @test integrator_rc.dtcache ≈ dt_coarse + + DiffEqBase.reinit!(integrator_rc; dt = dt_fine) + @test integrator_rc.dt ≈ dt_fine + @test integrator_rc.dtcache ≈ dt_fine + DiffEqBase.solve!(integrator_rc) + @test integrator_rc.sol.retcode == DiffEqBase.ReturnCode.Success + @test integrator_rc.dtcache ≈ dt_fine + end + + @testset "Nested instability propagation" begin + dt = 0.01π + + function ode_nan_nested(du, u, p, t) + du[1] = NaN + du[2] = 0.01u[1] + end + f_nan_nested = ODEFunction(ode_nan_nested) + + # Inner split: one leg produces NaN, nested inside an outer split. + # f3dofs_n = [1,2] indexes into the 2-element view selected by f2dofs = [1,3]. + fsplit_inner_n = GenericSplitFunction((f3, f_nan_nested), ([1, 2], [1, 2])) + fsplit_outer_n = GenericSplitFunction((f1, fsplit_inner_n), ([1, 2, 3], [1, 3])) + prob_nested = OperatorSplittingProblem(fsplit_outer_n, u0, tspan) + + tstepper_n = LieTrotterGodunov((Euler(), LieTrotterGodunov((Euler(), Euler())))) + integrator_n = DiffEqBase.init(prob_nested, tstepper_n; dt = dt, alias_u0 = false) + DiffEqBase.solve!(integrator_n) + @test integrator_n.sol.retcode ∈ (DiffEqBase.ReturnCode.Unstable, DiffEqBase.ReturnCode.DtNaN) + end + + @testset "verbose=false suppresses warnings without affecting retcode" begin + dt = 0.01π + + function ode_nan_quiet(du, u, p, t) + du[1] = NaN + du[2] = 0.01u[1] + end + f_nan_quiet = ODEFunction(ode_nan_quiet) + + fsplit_q = GenericSplitFunction((f1, f_nan_quiet), ([1, 2, 3], [1, 3])) + prob_q = OperatorSplittingProblem(fsplit_q, u0, tspan) + + integrator_q = DiffEqBase.init( + prob_q, LieTrotterGodunov((Euler(), Euler())); + dt = dt, verbose = false, alias_u0 = false + ) + DiffEqBase.solve!(integrator_q) + @test integrator_q.sol.retcode ∈ (DiffEqBase.ReturnCode.Unstable, DiffEqBase.ReturnCode.DtNaN) + end end From 9243e06177b854e6f7c9276a2ddd33373d719d8c Mon Sep 17 00:00:00 2001 From: Orjan Ameye Date: Fri, 15 May 2026 11:26:44 +0200 Subject: [PATCH 2/2] move to IntegratorStats again --- Project.toml | 8 +++++--- src/OrdinaryDiffEqOperatorSplitting.jl | 2 +- src/integrator.jl | 16 ++++++++++++---- test/consistency.jl | 17 ++++------------- 4 files changed, 22 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 892a798..43c9a25 100644 --- a/Project.toml +++ b/Project.toml @@ -6,8 +6,8 @@ authors = ["Dennis Ogiermann and contr [deps] CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -Reuse OrdinaryDiffEqCore infrastructure = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8" OrdinaryDiffEqLowOrderRK = "1344f307-1e59-4825-a18e-ace9aa3fa4c6" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" @@ -21,7 +21,7 @@ CommonSolve = "0.2.4" DataStructures = "0.18.22, 0.19" DiffEqBase = "7" ExplicitImports = "1" -ModelingToolkit = "10.31.2, 11" +ModelingToolkit = "11" OrdinaryDiffEqCore = "4" OrdinaryDiffEqLowOrderRK = "2" OrdinaryDiffEqTsit5 = "2" @@ -29,6 +29,7 @@ PrecompileTools = "1.0" RecursiveArrayTools = "4" SafeTestsets = "0.1.0" SciMLBase = "3.1" +SciMLIterators = "1" TimerOutputs = "0.5.28" Unrolled = "0.1.5" julia = "1.10" @@ -38,7 +39,8 @@ ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +SciMLIterators = "efe4dcdd-3aed-4391-894d-a9dbd16a2f14" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["ExplicitImports", "ModelingToolkit", "OrdinaryDiffEqTsit5", "SafeTestsets", "Test"] +test = ["ExplicitImports", "ModelingToolkit", "OrdinaryDiffEqTsit5", "SafeTestsets", "SciMLIterators", "Test"] diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl index d7a534d..11471ee 100644 --- a/src/OrdinaryDiffEqOperatorSplitting.jl +++ b/src/OrdinaryDiffEqOperatorSplitting.jl @@ -6,7 +6,7 @@ timeit_debug_enabled() = false import Unrolled: @unroll import SciMLBase, DiffEqBase, DataStructures -import SciMLBase: ReturnCode, DEStats +import SciMLBase: ReturnCode import SciMLBase: DEIntegrator, NullParameters, isadaptive import RecursiveArrayTools diff --git a/src/integrator.jl b/src/integrator.jl index 5c5a8c0..19b9c4a 100644 --- a/src/integrator.jl +++ b/src/integrator.jl @@ -1,3 +1,11 @@ +mutable struct IntegratorStats + naccept::Int64 + nreject::Int64 + # TODO inner solver stats +end + +IntegratorStats() = IntegratorStats(0, 0) + Base.@kwdef mutable struct IntegratorOptions{tType, fType, F3} adaptive::Bool dtmin::tType = eps(Float64) @@ -86,7 +94,7 @@ mutable struct SplitSubIntegrator{ last_step_failed::Bool u_modified::Bool # TODO we can probably remove this status::SplitSubIntegratorStatus - stats::DEStats + stats::IntegratorStats cache::cacheType child_subintegrators::childSubintType # Tuple solution_indices::solidxType @@ -156,7 +164,7 @@ mutable struct OperatorSplittingIntegrator{ iter::Int controller::controllerType opts::optsType # DEOptions - stats::DEStats + stats::IntegratorStats tdir::tType end @@ -329,7 +337,7 @@ function SciMLBase.__init( 0, controller, opts, - DEStats(0), + IntegratorStats(), tdir_val, ) DiffEqBase.initialize!(callback, u0, t0, integrator) @@ -1083,7 +1091,7 @@ function _build_child( controller, false, false, false, # force_stepfail, last_step_failed, u_modified SplitSubIntegratorStatus(), - DEStats(0), + IntegratorStats(), level_cache, child_subintegrators, solution_indices, diff --git a/test/consistency.jl b/test/consistency.jl index e08729e..228df1d 100644 --- a/test/consistency.jl +++ b/test/consistency.jl @@ -1,5 +1,6 @@ using OrdinaryDiffEqLowOrderRK using OrdinaryDiffEqOperatorSplitting +using SciMLIterators: TimeChoiceIterator using Test f(du, u, p, t) = @. du = -u @@ -15,18 +16,8 @@ splitting_solver = LieTrotterGodunov((Euler(),)) integrator1 = init(prob1, splitting_solver; dt = dt) integrator2 = init(prob2, Euler(); dt = dt) - -# Compare solutions at every 2*dt time point. -# Original test used TimeChoiceIterator (removed in DiffEqBase v7). -nsteps = Int(round((tspan[2] - tspan[1]) / dt)) -@test integrator1.u ≈ integrator2.u -@test integrator1.t ≈ integrator2.t -for _ in 1:(nsteps ÷ 2) - step!(integrator1) - step!(integrator1) - step!(integrator2) - step!(integrator2) - @test integrator1.u ≈ integrator2.u - @test integrator1.t ≈ integrator2.t +for ((u1, t1), (u2, t2)) in zip(TimeChoiceIterator(integrator1, tspan[1]:(2dt):tspan[2]), TimeChoiceIterator(integrator2, tspan[1]:(2dt):tspan[2])) + @test u1 ≈ u2 + @test t1 ≈ t2 end @test integrator1.iter == integrator2.iter