From b6aaeec9b347ff518a0bad7d333230da2b745603 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 13 Aug 2024 19:56:09 +0200 Subject: [PATCH 01/37] init --- examples/conduction-velocity-benchmark.jl | 2 +- src/solver/operator_splitting.jl | 1 + src/solver/operator_splitting/adaptivity.jl | 14 +++++++++++++ src/solver/operator_splitting/integrator.jl | 22 +++++++++++++++------ src/solver/operator_splitting/solver.jl | 13 +++++++----- src/solver/time/time_integrator.jl | 8 ++++++++ 6 files changed, 48 insertions(+), 12 deletions(-) create mode 100644 src/solver/operator_splitting/adaptivity.jl diff --git a/examples/conduction-velocity-benchmark.jl b/examples/conduction-velocity-benchmark.jl index 9c3db53b..3d41d61d 100644 --- a/examples/conduction-velocity-benchmark.jl +++ b/examples/conduction-velocity-benchmark.jl @@ -74,7 +74,7 @@ timestepper = OS.LieTrotterGodunov(( solution_vector_type=Vector{Float32}, reaction_threshold=0.1f0, ), -)) +), OS.ReactionTangentController(0.5, 1.0, (0.01, 0.3))) problem = OS.OperatorSplittingProblem(odeform, u₀, tspan) diff --git a/src/solver/operator_splitting.jl b/src/solver/operator_splitting.jl index 5a71fae6..beac1104 100644 --- a/src/solver/operator_splitting.jl +++ b/src/solver/operator_splitting.jl @@ -12,6 +12,7 @@ abstract type AbstractOperatorSplittingAlgorithm end abstract type AbstractOperatorSplittingCache end include("operator_splitting/utils.jl") +include("operator_splitting/adaptivity.jl") include("operator_splitting/function.jl") include("operator_splitting/problem.jl") include("operator_splitting/integrator.jl") diff --git a/src/solver/operator_splitting/adaptivity.jl b/src/solver/operator_splitting/adaptivity.jl new file mode 100644 index 00000000..57cdee2c --- /dev/null +++ b/src/solver/operator_splitting/adaptivity.jl @@ -0,0 +1,14 @@ +abstract type AbstractTimeAdaptionAlgorithm end + +struct NoTimeAdaption <: AbstractTimeAdaptionAlgorithm end + +struct ReactionTangentController{T <: Real} <: AbstractTimeAdaptionAlgorithm + σ_s::T + σ_c::T + Δt_bounds::NTuple{2,T} +end + +function get_next_dt(R, controller::ReactionTangentController ) + @unpack σ_s, σ_c, Δt_bounds = controller + (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] +end diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index a0a191ae..3b7528dd 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -45,7 +45,7 @@ end # called by DiffEqBase.init and DiffEqBase.solve function DiffEqBase.__init( prob::OperatorSplittingProblem, - alg::AbstractOperatorSplittingAlgorithm, + alg::LieTrotterGodunov{<:Any, AdaptivityController}, args...; dt, tstops = (), @@ -56,7 +56,7 @@ function DiffEqBase.__init( save_func = (u, t) -> copy(u), # custom kwarg dtchangeable = true, # custom kwarg kwargs..., -) +) where AdaptivityController (; u0, p) = prob t0, tf = prob.tspan @@ -94,6 +94,7 @@ function DiffEqBase.__init( callback, advance_to_tstop, false, + # AdaptivityController isa NoTimeAdaption ? false : true, cache, sol, subintegrators, @@ -217,21 +218,25 @@ end function DiffEqBase.postamble!(integrator::OperatorSplittingIntegrator) DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator) end - -function __step!(integrator) +@inline function get_reaction_tangent(_) + return 0.0 +end +function __step!(integrator::OperatorSplittingIntegrator{<:Any, <:LieTrotterGodunov{<:Any, AdaptivityController}}) where AdaptivityController (; dtchangeable, tstops) = integrator _dt = DiffEqBase.get_dt(integrator) # update dt before incrementing u; if dt is changeable and there is # a tstop within dt, reduce dt to tstop - t integrator.dt = - !isempty(tstops) && dtchangeable ? tdir(integrator) * min(_dt, abs(first(tstops) - integrator.t)) : - tdir(integrator) * _dt + !isempty(tstops) && dtchangeable ? tdir(integrator) * min(integrator.dt, abs(first(tstops) - integrator.t)) : + tdir(integrator) * integrator.dt tnext = integrator.t + integrator.dt + R = sum(get_reaction_tangent.(integrator.subintegrators)) # Solve inner problems advance_solution_to!(integrator, tnext) + R = max(R, sum(get_reaction_tangent.(integrator.subintegrators))) # Update integrator # increment t by dt, rounding to the first tstop if that is roughly @@ -242,6 +247,11 @@ function __step!(integrator) integrator.tprev = integrator.t integrator.t = !isempty(tstops) && abs(first(tstops) - tnext) < max_t_error ? first(tstops) : tnext + controller = integrator.alg.time_adaption_alg + @info "before" integrator.dt + integrator.dt = get_next_dt(R, controller) + @info "after" integrator.dt + # Propagate information down into the subintegrators synchronize_subintegrators!(integrator) diff --git a/src/solver/operator_splitting/solver.jl b/src/solver/operator_splitting/solver.jl index d6294c6b..1ff6bf4f 100644 --- a/src/solver/operator_splitting/solver.jl +++ b/src/solver/operator_splitting/solver.jl @@ -4,10 +4,14 @@ LieTrotterGodunov <: AbstractOperatorSplittingAlgorithm A first order operator splitting algorithm attributed to [Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite). """ -struct LieTrotterGodunov{AlgTupleType} <: AbstractOperatorSplittingAlgorithm +struct LieTrotterGodunov{AlgTupleType, TimeAdaptionAlgorithm <: AbstractTimeAdaptionAlgorithm} <: AbstractOperatorSplittingAlgorithm inner_algs::AlgTupleType # Tuple of timesteppers for inner problems + time_adaption_alg::TimeAdaptionAlgorithm # transfer_algs::TransferTupleType # Tuple of transfer algorithms from the master solution into the individual ones end +function LieTrotterGodunov(inner_algs) + LieTrotterGodunov(inner_algs, NoTimeAdaption()) +end struct LieTrotterGodunovCache{uType, tmpType, iiType} <: AbstractOperatorSplittingCache u::uType @@ -18,9 +22,9 @@ struct LieTrotterGodunovCache{uType, tmpType, iiType} <: AbstractOperatorSplitti end # Dispatch for outer construction -function init_cache(prob::OperatorSplittingProblem, alg::LieTrotterGodunov; dt, kwargs...) # TODO +function init_cache(prob::OperatorSplittingProblem{<:GenericSplitFunction}, alg::LieTrotterGodunov; dt, kwargs...) # TODO @unpack f = prob - @assert f isa GenericSplitFunction + # @assert f isa GenericSplitFunction u = copy(prob.u0) uprev = copy(prob.u0) @@ -54,5 +58,4 @@ end advance_solution_to!(subinteg, inner_caches[i], tnext) finalize_local_step!(subinteg) end -end - +end diff --git a/src/solver/time/time_integrator.jl b/src/solver/time/time_integrator.jl index c4909852..8b556fd2 100644 --- a/src/solver/time/time_integrator.jl +++ b/src/solver/time/time_integrator.jl @@ -217,3 +217,11 @@ function perform_step!(integ::ThunderboltTimeIntegrator, cache::AbstractTimeSolv end return true end + +@inline function OS.get_reaction_tangent(_::ThunderboltTimeIntegrator) + return 0.0 +end +@inline function OS.get_reaction_tangent(integrator::ThunderboltTimeIntegrator{<:PointwiseODEFunction}) + φₘidx = transmembranepotential_index(integrator.f.ode) + return maximum(@view integrator.cache.dumat[:, φₘidx]) +end From 83d4b6b11b24b987a3bee10a5ace3e8eb3ce90b5 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 14 Aug 2024 10:55:54 +0200 Subject: [PATCH 02/37] now it stops --- src/solver/operator_splitting.jl | 2 +- src/solver/operator_splitting/integrator.jl | 18 ++++++++---------- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/src/solver/operator_splitting.jl b/src/solver/operator_splitting.jl index beac1104..4f02f713 100644 --- a/src/solver/operator_splitting.jl +++ b/src/solver/operator_splitting.jl @@ -15,8 +15,8 @@ include("operator_splitting/utils.jl") include("operator_splitting/adaptivity.jl") include("operator_splitting/function.jl") include("operator_splitting/problem.jl") -include("operator_splitting/integrator.jl") include("operator_splitting/solver.jl") +include("operator_splitting/integrator.jl") export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, DiffEqBase, init, TimeChoiceIterator, diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index 3b7528dd..d60186ef 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -54,7 +54,7 @@ function DiffEqBase.__init( callback = nothing, advance_to_tstop = false, save_func = (u, t) -> copy(u), # custom kwarg - dtchangeable = true, # custom kwarg + dtchangeable = AdaptivityController isa NoTimeAdaption ? false : true, # custom kwarg kwargs..., ) where AdaptivityController (; u0, p) = prob @@ -94,7 +94,6 @@ function DiffEqBase.__init( callback, advance_to_tstop, false, - # AdaptivityController isa NoTimeAdaption ? false : true, cache, sol, subintegrators, @@ -228,8 +227,12 @@ function __step!(integrator::OperatorSplittingIntegrator{<:Any, <:LieTrotterGodu # update dt before incrementing u; if dt is changeable and there is # a tstop within dt, reduce dt to tstop - t integrator.dt = - !isempty(tstops) && dtchangeable ? tdir(integrator) * min(integrator.dt, abs(first(tstops) - integrator.t)) : - tdir(integrator) * integrator.dt + !isempty(tstops) && dtchangeable ? tdir(integrator) * min(_dt, abs(first(tstops) - integrator.t)) : + tdir(integrator) * _dt + + + # Propagate information down into the subintegrators + synchronize_subintegrators!(integrator) tnext = integrator.t + integrator.dt R = sum(get_reaction_tangent.(integrator.subintegrators)) @@ -248,12 +251,7 @@ function __step!(integrator::OperatorSplittingIntegrator{<:Any, <:LieTrotterGodu integrator.t = !isempty(tstops) && abs(first(tstops) - tnext) < max_t_error ? first(tstops) : tnext controller = integrator.alg.time_adaption_alg - @info "before" integrator.dt - integrator.dt = get_next_dt(R, controller) - @info "after" integrator.dt - - # Propagate information down into the subintegrators - synchronize_subintegrators!(integrator) + integrator._dt = get_next_dt(R, controller) # remove tstops that were just reached while !isempty(tstops) && reached_tstop(integrator, first(tstops)) From 93d1eba20f5319d4abbb8d6793d1bca0eb0f479c Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 14 Aug 2024 12:09:56 +0200 Subject: [PATCH 03/37] . --- src/solver/operator_splitting/integrator.jl | 6 +++--- src/utils.jl | 2 +- test/test_integrators.jl | 3 +++ 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index d60186ef..122e1afc 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -235,11 +235,11 @@ function __step!(integrator::OperatorSplittingIntegrator{<:Any, <:LieTrotterGodu synchronize_subintegrators!(integrator) tnext = integrator.t + integrator.dt - R = sum(get_reaction_tangent.(integrator.subintegrators)) + (AdaptivityController == NoTimeAdaption) || (R = sum(get_reaction_tangent.(integrator.subintegrators))) # Solve inner problems advance_solution_to!(integrator, tnext) - R = max(R, sum(get_reaction_tangent.(integrator.subintegrators))) + (AdaptivityController == NoTimeAdaption) || (R = max(R, sum(get_reaction_tangent.(integrator.subintegrators)))) # Update integrator # increment t by dt, rounding to the first tstop if that is roughly @@ -251,7 +251,7 @@ function __step!(integrator::OperatorSplittingIntegrator{<:Any, <:LieTrotterGodu integrator.t = !isempty(tstops) && abs(first(tstops) - tnext) < max_t_error ? first(tstops) : tnext controller = integrator.alg.time_adaption_alg - integrator._dt = get_next_dt(R, controller) + (AdaptivityController == NoTimeAdaption) || (integrator._dt = get_next_dt(R, controller)) # remove tstops that were just reached while !isempty(tstops) && reached_tstop(integrator, first(tstops)) diff --git a/src/utils.jl b/src/utils.jl index 81420870..4a32278c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -377,7 +377,7 @@ end # Transfer the element data into a vector function ea_collapse!(b::Vector, bes::EAVector) ndofs = size(b, 1) - @batch minbatch=ndofs÷Threads.nthreads() for dof ∈ 1:ndofs + @batch minbatch=max(1, ndofs÷Threads.nthreads()) for dof ∈ 1:ndofs _ea_collapse_kernel!(b, dof, bes) end end diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 0e1b2d4f..50851b52 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -73,6 +73,9 @@ timestepper = LieTrotterGodunov( (DummyForwardEuler(), DummyForwardEuler()) ) +timestepper_adaptive = LieTrotterGodunov( + (DummyForwardEuler(), DummyForwardEuler()) + , OS.ReactionTangentController(0.5, 1.0, (0.01, 0.3))) # The remaining code works as usual. integrator = DiffEqBase.init(prob, timestepper, dt=0.01, verbose=true) DiffEqBase.solve!(integrator) From ec3bd11c773a2a2a04fe5fcfe9e8b589c10b0721 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 14 Aug 2024 12:32:37 +0200 Subject: [PATCH 04/37] quick quick --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index 4a32278c..19ddfa1c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -160,7 +160,7 @@ function mul!(y::AbstractVector{<:Number}, A_::ThreadedSparseMatrixCSR, x::Abstr A.n == size(x, 1) || throw(DimensionMismatch()) A.m == size(y, 1) || throw(DimensionMismatch()) - @batch minbatch = size(y, 1) ÷ Threads.nthreads() for row in 1:size(y, 1) + @batch minbatch = max(1, size(y, 1) ÷ Threads.nthreads()) for row in 1:size(y, 1) @inbounds begin v = zero(eltype(y)) for nz in nzrange(A, row) From 85bc9ba6b087e54421bc9f777fe6aea445555d25 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 14 Aug 2024 15:56:35 +0200 Subject: [PATCH 05/37] restructuring --- examples/conduction-velocity-benchmark.jl | 7 +- src/Thunderbolt.jl | 2 + src/solver/adaptivity.jl | 55 +++++++++ src/solver/operator_splitting.jl | 1 - src/solver/operator_splitting/adaptivity.jl | 14 --- src/solver/operator_splitting/integrator.jl | 23 ++-- src/solver/operator_splitting/solver.jl | 12 +- src/solver/time/time_integrator.jl | 8 -- test/test_integrators.jl | 120 ++++++++++---------- 9 files changed, 138 insertions(+), 104 deletions(-) create mode 100644 src/solver/adaptivity.jl delete mode 100644 src/solver/operator_splitting/adaptivity.jl diff --git a/examples/conduction-velocity-benchmark.jl b/examples/conduction-velocity-benchmark.jl index 3d41d61d..723cacdb 100644 --- a/examples/conduction-velocity-benchmark.jl +++ b/examples/conduction-velocity-benchmark.jl @@ -64,7 +64,7 @@ steady_state_initializer!(u₀, odeform) # io = ParaViewWriter("spiral-wave-test") -timestepper = OS.LieTrotterGodunov(( +_timestepper = OS.LieTrotterGodunov(( BackwardEulerSolver( solution_vector_type=Vector{Float32}, system_matrix_type=Thunderbolt.ThreadedSparseMatrixCSR{Float32, Int32}, @@ -74,7 +74,10 @@ timestepper = OS.LieTrotterGodunov(( solution_vector_type=Vector{Float32}, reaction_threshold=0.1f0, ), -), OS.ReactionTangentController(0.5, 1.0, (0.01, 0.3))) +)) +controller = Thunderbolt.ReactionTangentController(0.5, 1.0, (0.01, 0.3), 0.0, 0.0) + +timestepper = Thunderbolt.AdaptiveOperatorSplittingAlgorithm(_timestepper, controller) problem = OS.OperatorSplittingProblem(odeform, u₀, tspan) diff --git a/src/Thunderbolt.jl b/src/Thunderbolt.jl index 513b274b..d61c2e9e 100644 --- a/src/Thunderbolt.jl +++ b/src/Thunderbolt.jl @@ -71,6 +71,8 @@ include("solver/interface.jl") include("solver/linear.jl") include("solver/nonlinear.jl") include("solver/time_integration.jl") +include("solver/adaptivity.jl") + include("processing/ecg.jl") diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl new file mode 100644 index 00000000..c77bde9b --- /dev/null +++ b/src/solver/adaptivity.jl @@ -0,0 +1,55 @@ +abstract type AbstractTimeAdaptionAlgorithm end + +struct NoTimeAdaption <: AbstractTimeAdaptionAlgorithm end + +mutable struct ReactionTangentController{T <: Real} <: AbstractTimeAdaptionAlgorithm + const σ_s::T + const σ_c::T + const Δt_bounds::NTuple{2,T} + Rₙ₊₁::T + Rₙ::T +end + +function get_next_dt(controller::ReactionTangentController) + @unpack σ_s, σ_c, Δt_bounds, Rₙ₊₁, Rₙ = controller + R = max(Rₙ, Rₙ₊₁) + (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] +end + +struct AdaptiveOperatorSplittingAlgorithm{TOperatorSplittingAlg <: OS.AbstractOperatorSplittingAlgorithm, TTimeAdaptionAlgorithm <: AbstractTimeAdaptionAlgorithm} <: OS.AbstractOperatorSplittingAlgorithm + operator_splitting_algorithm::TOperatorSplittingAlg + timestep_controller::TTimeAdaptionAlgorithm +end + +is_adaptive(::AdaptiveOperatorSplittingAlgorithm) = true + + +@inline function get_reaction_tangent(integrator::OS.OperatorSplittingIntegrator) + for subintegrator in integrator.subintegrators + subintegrator.f isa PointwiseODEFunction || continue + φₘidx = transmembranepotential_index(subintegrator.f.ode) + return maximum(@view subintegrator.cache.dumat[:, φₘidx]) + end +end + + +@inline function OS.update_controller!(integrator::OS.OperatorSplittingIntegrator{<:Any, <:AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov, <:ReactionTangentController}}) + controller = integrator.alg.timestep_controller + controller.Rₙ = controller.Rₙ₊₁ + controller.Rₙ₊₁ = get_reaction_tangent(integrator) +end + +@inline function OS.update_dt!(integrator::OS.OperatorSplittingIntegrator{<:Any, <:AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov, <:ReactionTangentController}}) + controller = integrator.alg.timestep_controller + integrator._dt = get_next_dt(controller) +end + +# Dispatch for outer construction +function OS.init_cache(prob::OS.OperatorSplittingProblem, alg::AdaptiveOperatorSplittingAlgorithm; dt, kwargs...) # TODO + OS.init_cache(prob, alg.operator_splitting_algorithm;dt = dt, kwargs...) +end + +# Dispatch for recursive construction +function OS.construct_inner_cache(f::OS.AbstractOperatorSplitFunction, alg::AdaptiveOperatorSplittingAlgorithm, u::AbstractArray, uprev::AbstractArray) + OS.construct_inner_cache(f, alg.operator_splitting_algorithm, u, uprev) +end \ No newline at end of file diff --git a/src/solver/operator_splitting.jl b/src/solver/operator_splitting.jl index 4f02f713..ff550adc 100644 --- a/src/solver/operator_splitting.jl +++ b/src/solver/operator_splitting.jl @@ -12,7 +12,6 @@ abstract type AbstractOperatorSplittingAlgorithm end abstract type AbstractOperatorSplittingCache end include("operator_splitting/utils.jl") -include("operator_splitting/adaptivity.jl") include("operator_splitting/function.jl") include("operator_splitting/problem.jl") include("operator_splitting/solver.jl") diff --git a/src/solver/operator_splitting/adaptivity.jl b/src/solver/operator_splitting/adaptivity.jl deleted file mode 100644 index 57cdee2c..00000000 --- a/src/solver/operator_splitting/adaptivity.jl +++ /dev/null @@ -1,14 +0,0 @@ -abstract type AbstractTimeAdaptionAlgorithm end - -struct NoTimeAdaption <: AbstractTimeAdaptionAlgorithm end - -struct ReactionTangentController{T <: Real} <: AbstractTimeAdaptionAlgorithm - σ_s::T - σ_c::T - Δt_bounds::NTuple{2,T} -end - -function get_next_dt(R, controller::ReactionTangentController ) - @unpack σ_s, σ_c, Δt_bounds = controller - (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] -end diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index 122e1afc..9a588aa8 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -45,7 +45,7 @@ end # called by DiffEqBase.init and DiffEqBase.solve function DiffEqBase.__init( prob::OperatorSplittingProblem, - alg::LieTrotterGodunov{<:Any, AdaptivityController}, + alg::AbstractOperatorSplittingAlgorithm, args...; dt, tstops = (), @@ -54,9 +54,9 @@ function DiffEqBase.__init( callback = nothing, advance_to_tstop = false, save_func = (u, t) -> copy(u), # custom kwarg - dtchangeable = AdaptivityController isa NoTimeAdaption ? false : true, # custom kwarg + dtchangeable = is_adaptive(alg), # custom kwarg kwargs..., -) where AdaptivityController +) (; u0, p) = prob t0, tf = prob.tspan @@ -184,7 +184,8 @@ function (integrator::OperatorSplittingIntegrator)(tmp, t) linear_interpolation!(tmp, t, integrator.uprev, integrator.u, integrator.tprev, integrator.t) end - +@inline update_controller!(::OperatorSplittingIntegrator) = nothing +@inline update_dt!(::OperatorSplittingIntegrator) = nothing # helper functions for dealing with time-reversed integrators in the same way # that OrdinaryDiffEq.jl does @@ -217,10 +218,8 @@ end function DiffEqBase.postamble!(integrator::OperatorSplittingIntegrator) DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator) end -@inline function get_reaction_tangent(_) - return 0.0 -end -function __step!(integrator::OperatorSplittingIntegrator{<:Any, <:LieTrotterGodunov{<:Any, AdaptivityController}}) where AdaptivityController + +function __step!(integrator) (; dtchangeable, tstops) = integrator _dt = DiffEqBase.get_dt(integrator) @@ -230,16 +229,13 @@ function __step!(integrator::OperatorSplittingIntegrator{<:Any, <:LieTrotterGodu !isempty(tstops) && dtchangeable ? tdir(integrator) * min(_dt, abs(first(tstops) - integrator.t)) : tdir(integrator) * _dt - # Propagate information down into the subintegrators synchronize_subintegrators!(integrator) - tnext = integrator.t + integrator.dt - (AdaptivityController == NoTimeAdaption) || (R = sum(get_reaction_tangent.(integrator.subintegrators))) # Solve inner problems advance_solution_to!(integrator, tnext) - (AdaptivityController == NoTimeAdaption) || (R = max(R, sum(get_reaction_tangent.(integrator.subintegrators)))) + update_controller!(integrator) # Update integrator # increment t by dt, rounding to the first tstop if that is roughly @@ -250,8 +246,7 @@ function __step!(integrator::OperatorSplittingIntegrator{<:Any, <:LieTrotterGodu integrator.tprev = integrator.t integrator.t = !isempty(tstops) && abs(first(tstops) - tnext) < max_t_error ? first(tstops) : tnext - controller = integrator.alg.time_adaption_alg - (AdaptivityController == NoTimeAdaption) || (integrator._dt = get_next_dt(R, controller)) + update_dt!(integrator) # remove tstops that were just reached while !isempty(tstops) && reached_tstop(integrator, first(tstops)) diff --git a/src/solver/operator_splitting/solver.jl b/src/solver/operator_splitting/solver.jl index 1ff6bf4f..65a8a789 100644 --- a/src/solver/operator_splitting/solver.jl +++ b/src/solver/operator_splitting/solver.jl @@ -4,14 +4,12 @@ LieTrotterGodunov <: AbstractOperatorSplittingAlgorithm A first order operator splitting algorithm attributed to [Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite). """ -struct LieTrotterGodunov{AlgTupleType, TimeAdaptionAlgorithm <: AbstractTimeAdaptionAlgorithm} <: AbstractOperatorSplittingAlgorithm +struct LieTrotterGodunov{AlgTupleType} <: AbstractOperatorSplittingAlgorithm inner_algs::AlgTupleType # Tuple of timesteppers for inner problems - time_adaption_alg::TimeAdaptionAlgorithm # transfer_algs::TransferTupleType # Tuple of transfer algorithms from the master solution into the individual ones end -function LieTrotterGodunov(inner_algs) - LieTrotterGodunov(inner_algs, NoTimeAdaption()) -end + +@inline is_adaptive(::AbstractOperatorSplittingAlgorithm) = false struct LieTrotterGodunovCache{uType, tmpType, iiType} <: AbstractOperatorSplittingCache u::uType @@ -22,9 +20,9 @@ struct LieTrotterGodunovCache{uType, tmpType, iiType} <: AbstractOperatorSplitti end # Dispatch for outer construction -function init_cache(prob::OperatorSplittingProblem{<:GenericSplitFunction}, alg::LieTrotterGodunov; dt, kwargs...) # TODO +function init_cache(prob::OperatorSplittingProblem, alg::LieTrotterGodunov; dt, kwargs...) # TODO @unpack f = prob - # @assert f isa GenericSplitFunction + @assert f isa GenericSplitFunction u = copy(prob.u0) uprev = copy(prob.u0) diff --git a/src/solver/time/time_integrator.jl b/src/solver/time/time_integrator.jl index 8b556fd2..c4909852 100644 --- a/src/solver/time/time_integrator.jl +++ b/src/solver/time/time_integrator.jl @@ -217,11 +217,3 @@ function perform_step!(integ::ThunderboltTimeIntegrator, cache::AbstractTimeSolv end return true end - -@inline function OS.get_reaction_tangent(_::ThunderboltTimeIntegrator) - return 0.0 -end -@inline function OS.get_reaction_tangent(integrator::ThunderboltTimeIntegrator{<:PointwiseODEFunction}) - φₘidx = transmembranepotential_index(integrator.f.ode) - return maximum(@view integrator.cache.dumat[:, φₘidx]) -end diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 50851b52..73eb534e 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -69,66 +69,70 @@ tspan = (0.0,100.0) prob = OperatorSplittingProblem(fsplit1, u0, tspan) # The time stepper carries the individual solver information. -timestepper = LieTrotterGodunov( - (DummyForwardEuler(), DummyForwardEuler()) -) - -timestepper_adaptive = LieTrotterGodunov( - (DummyForwardEuler(), DummyForwardEuler()) - , OS.ReactionTangentController(0.5, 1.0, (0.01, 0.3))) -# The remaining code works as usual. -integrator = DiffEqBase.init(prob, timestepper, dt=0.01, verbose=true) -DiffEqBase.solve!(integrator) -ufinal = copy(integrator.u) -@test ufinal ≉ u0 # Make sure the solve did something - -DiffEqBase.reinit!(integrator, u0; tspan) -for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) -end -@test isapprox(ufinal, integrator.u, atol=1e-8) -DiffEqBase.reinit!(integrator, u0; tspan) -for (uprev, tprev, u, t) in DiffEqBase.intervals(integrator) +@testset "OperatorSplitting" begin + for TimeStepperType in (LieTrotterGodunov,) + for controller in (OS.NoTimeAdaption(), + OS.ReactionTangentController(0.5, 1.0, (0.01, 0.3)),) + timestepper = TimeStepperType( + (DummyForwardEuler(), DummyForwardEuler()) + ) + # The remaining code works as usual. + integrator = DiffEqBase.init(prob, timestepper, dt=0.01, verbose=true) + DiffEqBase.solve!(integrator) + ufinal = copy(integrator.u) + @test ufinal ≉ u0 # Make sure the solve did something + + DiffEqBase.reinit!(integrator, u0; tspan) + for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) + end + @test isapprox(ufinal, integrator.u, atol=1e-8) + + DiffEqBase.reinit!(integrator, u0; tspan) + for (uprev, tprev, u, t) in DiffEqBase.intervals(integrator) + end + @test isapprox(ufinal, integrator.u, atol=1e-8) + + DiffEqBase.reinit!(integrator, u0; tspan) + DiffEqBase.solve!(integrator) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success + + # Now some recursive splitting + function ode3(du, u, p, t) + du[1] = 0.005u[2] + du[2] = 0.005u[1] + end + f3 = ODEFunction(ode3) + + # Note that we define the dof indices w.r.t the parent function. + # Hence the indices for `fsplit2_inner` are. + f1dofs = [1,2,3] + f2dofs = [1,3] + f3dofs = [1,3] + fsplit2_inner = GenericSplitFunction((f3,f3), (f3dofs, f3dofs)) + fsplit2_outer = GenericSplitFunction((f1,fsplit2_inner), (f1dofs, f2dofs)) + + timestepper_inner = TimeStepperType( + (DummyForwardEuler(), DummyForwardEuler()) + ) + timestepper2 = TimeStepperType( + (DummyForwardEuler(), timestepper_inner) + ) + prob2 = OperatorSplittingProblem(fsplit2_outer, u0, tspan) + integrator2 = DiffEqBase.init(prob2, timestepper2, dt=0.01, verbose=true) + DiffEqBase.solve!(integrator2) + + DiffEqBase.reinit!(integrator2, u0; tspan) + for (u, t) in DiffEqBase.TimeChoiceIterator(integrator2, 0.0:5.0:100.0) + end + @test isapprox(ufinal, integrator2.u, atol=1e-8) + + DiffEqBase.reinit!(integrator2, u0; tspan) + DiffEqBase.solve!(integrator2) + @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Success + end + end end -@test isapprox(ufinal, integrator.u, atol=1e-8) - -DiffEqBase.reinit!(integrator, u0; tspan) -DiffEqBase.solve!(integrator) -@test integrator.sol.retcode == DiffEqBase.ReturnCode.Success - -# Now some recursive splitting -function ode3(du, u, p, t) - du[1] = 0.005u[2] - du[2] = 0.005u[1] -end -f3 = ODEFunction(ode3) - -# Note that we define the dof indices w.r.t the parent function. -# Hence the indices for `fsplit2_inner` are. -f1dofs = [1,2,3] -f2dofs = [1,3] -f3dofs = [1,3] -fsplit2_inner = GenericSplitFunction((f3,f3), (f3dofs, f3dofs)) -fsplit2_outer = GenericSplitFunction((f1,fsplit2_inner), (f1dofs, f2dofs)) - -timestepper_inner = LieTrotterGodunov( - (DummyForwardEuler(), DummyForwardEuler()) -) -timestepper2 = LieTrotterGodunov( - (DummyForwardEuler(), timestepper_inner) -) -prob2 = OperatorSplittingProblem(fsplit2_outer, u0, tspan) -integrator2 = DiffEqBase.init(prob2, timestepper2, dt=0.01, verbose=true) -DiffEqBase.solve!(integrator2) - -DiffEqBase.reinit!(integrator2, u0; tspan) -for (u, t) in DiffEqBase.TimeChoiceIterator(integrator2, 0.0:5.0:100.0) -end -@test isapprox(ufinal, integrator2.u, atol=1e-8) - -DiffEqBase.reinit!(integrator2, u0; tspan) -DiffEqBase.solve!(integrator2) -@test integrator2.sol.retcode == DiffEqBase.ReturnCode.Success # tnext = tspan[1]+0.01 # @btime OS.advance_solution_to!($integrator, $tnext) setup=(DiffEqBase.reinit!(integrator, u0; tspan)) From a5ccab2f700ab209e0d826a42374a30b96f8ca68 Mon Sep 17 00:00:00 2001 From: PotatoTuna Date: Wed, 14 Aug 2024 20:13:59 +0200 Subject: [PATCH 06/37] test? --- test/test_integrators.jl | 104 +++++++++++++++++++++------------------ 1 file changed, 56 insertions(+), 48 deletions(-) diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 73eb534e..680d9f7d 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -68,68 +68,76 @@ u0 = [0.7611944793397108 tspan = (0.0,100.0) prob = OperatorSplittingProblem(fsplit1, u0, tspan) +# Now some recursive splitting +function ode3(du, u, p, t) + du[1] = 0.005u[2] + du[2] = 0.005u[1] +end +f3 = ODEFunction(ode3) # The time stepper carries the individual solver information. @testset "OperatorSplitting" begin + + + # Note that we define the dof indices w.r.t the parent function. + # Hence the indices for `fsplit2_inner` are. + f1dofs = [1,2,3] + f2dofs = [1,3] + f3dofs = [1,3] + fsplit2_inner = GenericSplitFunction((f3,f3), (f3dofs, f3dofs)) + fsplit2_outer = GenericSplitFunction((f1,fsplit2_inner), (f1dofs, f2dofs)) for TimeStepperType in (LieTrotterGodunov,) - for controller in (OS.NoTimeAdaption(), - OS.ReactionTangentController(0.5, 1.0, (0.01, 0.3)),) + for controller in (Thunderbolt.ReactionTangentController(0.5, 1.0, (0.01, 0.3), 0.0, 0.0),) timestepper = TimeStepperType( (DummyForwardEuler(), DummyForwardEuler()) ) - # The remaining code works as usual. - integrator = DiffEqBase.init(prob, timestepper, dt=0.01, verbose=true) - DiffEqBase.solve!(integrator) - ufinal = copy(integrator.u) - @test ufinal ≉ u0 # Make sure the solve did something - - DiffEqBase.reinit!(integrator, u0; tspan) - for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) - end - @test isapprox(ufinal, integrator.u, atol=1e-8) - - DiffEqBase.reinit!(integrator, u0; tspan) - for (uprev, tprev, u, t) in DiffEqBase.intervals(integrator) - end - @test isapprox(ufinal, integrator.u, atol=1e-8) - - DiffEqBase.reinit!(integrator, u0; tspan) - DiffEqBase.solve!(integrator) - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success - - # Now some recursive splitting - function ode3(du, u, p, t) - du[1] = 0.005u[2] - du[2] = 0.005u[1] - end - f3 = ODEFunction(ode3) - - # Note that we define the dof indices w.r.t the parent function. - # Hence the indices for `fsplit2_inner` are. - f1dofs = [1,2,3] - f2dofs = [1,3] - f3dofs = [1,3] - fsplit2_inner = GenericSplitFunction((f3,f3), (f3dofs, f3dofs)) - fsplit2_outer = GenericSplitFunction((f1,fsplit2_inner), (f1dofs, f2dofs)) - + timestepper_adaptive = Thunderbolt.AdaptiveOperatorSplittingAlgorithm(timestepper, controller) timestepper_inner = TimeStepperType( (DummyForwardEuler(), DummyForwardEuler()) ) + timestepper_inner_adaptive = Thunderbolt.AdaptiveOperatorSplittingAlgorithm(timestepper_inner, controller) timestepper2 = TimeStepperType( (DummyForwardEuler(), timestepper_inner) ) - prob2 = OperatorSplittingProblem(fsplit2_outer, u0, tspan) - integrator2 = DiffEqBase.init(prob2, timestepper2, dt=0.01, verbose=true) - DiffEqBase.solve!(integrator2) - - DiffEqBase.reinit!(integrator2, u0; tspan) - for (u, t) in DiffEqBase.TimeChoiceIterator(integrator2, 0.0:5.0:100.0) + timestepper2_adaptive = Thunderbolt.AdaptiveOperatorSplittingAlgorithm(timestepper2, controller) + + for (tstepper1, tstepper_inner, tstepper2) in ( + (timestepper, timestepper_inner, timestepper2), + (timestepper_adaptive, timestepper_inner_adaptive, timestepper2_adaptive) + ) + # The remaining code works as usual. + integrator = DiffEqBase.init(prob, timestepper, dt=0.01, verbose=true) + DiffEqBase.solve!(integrator) + ufinal = copy(integrator.u) + @test ufinal ≉ u0 # Make sure the solve did something + + DiffEqBase.reinit!(integrator, u0; tspan) + for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) + end + @test isapprox(ufinal, integrator.u, atol=1e-8) + + DiffEqBase.reinit!(integrator, u0; tspan) + for (uprev, tprev, u, t) in DiffEqBase.intervals(integrator) + end + @test isapprox(ufinal, integrator.u, atol=1e-8) + + DiffEqBase.reinit!(integrator, u0; tspan) + DiffEqBase.solve!(integrator) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success + + prob2 = OperatorSplittingProblem(fsplit2_outer, u0, tspan) + integrator2 = DiffEqBase.init(prob2, timestepper2, dt=0.01, verbose=true) + DiffEqBase.solve!(integrator2) + + DiffEqBase.reinit!(integrator2, u0; tspan) + for (u, t) in DiffEqBase.TimeChoiceIterator(integrator2, 0.0:5.0:100.0) + end + @test isapprox(ufinal, integrator2.u, atol=1e-8) + + DiffEqBase.reinit!(integrator2, u0; tspan) + DiffEqBase.solve!(integrator2) + @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Success end - @test isapprox(ufinal, integrator2.u, atol=1e-8) - - DiffEqBase.reinit!(integrator2, u0; tspan) - DiffEqBase.solve!(integrator2) - @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Success end end end From 0bf9c9ed21af604bec4f79ac29d50b0d48a5c125 Mon Sep 17 00:00:00 2001 From: PotatoTuna Date: Wed, 14 Aug 2024 20:54:58 +0200 Subject: [PATCH 07/37] error and a couple of returns --- src/solver/adaptivity.jl | 3 +++ src/solver/operator_splitting.jl | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index c77bde9b..acbb6daa 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -30,6 +30,7 @@ is_adaptive(::AdaptiveOperatorSplittingAlgorithm) = true φₘidx = transmembranepotential_index(subintegrator.f.ode) return maximum(@view subintegrator.cache.dumat[:, φₘidx]) end + @error "PointwiseODEFunction not found" end @@ -37,11 +38,13 @@ end controller = integrator.alg.timestep_controller controller.Rₙ = controller.Rₙ₊₁ controller.Rₙ₊₁ = get_reaction_tangent(integrator) + return nothing end @inline function OS.update_dt!(integrator::OS.OperatorSplittingIntegrator{<:Any, <:AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov, <:ReactionTangentController}}) controller = integrator.alg.timestep_controller integrator._dt = get_next_dt(controller) + return nothing end # Dispatch for outer construction diff --git a/src/solver/operator_splitting.jl b/src/solver/operator_splitting.jl index ff550adc..5a71fae6 100644 --- a/src/solver/operator_splitting.jl +++ b/src/solver/operator_splitting.jl @@ -14,8 +14,8 @@ abstract type AbstractOperatorSplittingCache end include("operator_splitting/utils.jl") include("operator_splitting/function.jl") include("operator_splitting/problem.jl") -include("operator_splitting/solver.jl") include("operator_splitting/integrator.jl") +include("operator_splitting/solver.jl") export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, DiffEqBase, init, TimeChoiceIterator, From aaf7eec049451528f60f97aa5a485542344063c8 Mon Sep 17 00:00:00 2001 From: PotatoTuna Date: Wed, 14 Aug 2024 21:20:10 +0200 Subject: [PATCH 08/37] some docs --- docs/src/api-reference/solver.md | 12 ++++++++ docs/src/assets/references.bib | 10 ++++++ src/solver/adaptivity.jl | 34 ++++++++++++++++++--- src/solver/operator_splitting/integrator.jl | 9 ++++++ src/solver/operator_splitting/solver.jl | 4 +++ 5 files changed, 65 insertions(+), 4 deletions(-) diff --git a/docs/src/api-reference/solver.md b/docs/src/api-reference/solver.md index 93d886e7..f7f0e364 100644 --- a/docs/src/api-reference/solver.md +++ b/docs/src/api-reference/solver.md @@ -34,3 +34,15 @@ Thunderbolt.OS.LieTrotterGodunov Thunderbolt.OS.GenericSplitFunction Thunderbolt.OS.OperatorSplittingIntegrator ``` + +## Operator Splitting Adaptivity + +```@docs +Thunderbolt.AdaptiveOperatorSplittingAlgorithm +Thunderbolt.OS.is_adaptive +Thunderbolt.ReactionTangentController +Thunderbolt.get_reaction_tangent +Thunderbolt.OS.update_controller! +Thunderbolt.OS.update_dt! +Thunderbolt.get_next_dt +``` diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index d22e951e..e905f0ad 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -278,3 +278,13 @@ @article{PotDubRicVinGul:2006:cmb pages={2425-2435}, doi={10.1109/TBME.2006.880875} } +@article{OgiBalPer:2023:seats, +author = {Ogiermann, Dennis and Perotti, Luigi E. and Balzani, Daniel}, +title = {A simple and efficient adaptive time stepping technique for low-order operator splitting schemes applied to cardiac electrophysiology}, +journal = {International Journal for Numerical Methods in Biomedical Engineering}, +volume = {39}, +number = {2}, +pages = {e3670}, +url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/cnm.3670}, +year = {2023} +} diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index acbb6daa..4888ff44 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -1,7 +1,16 @@ abstract type AbstractTimeAdaptionAlgorithm end -struct NoTimeAdaption <: AbstractTimeAdaptionAlgorithm end - +""" + ReactionTangentController{T <: Real} <: AbstractTimeAdaptionAlgorithm +A timestep length controller for [`LieTrotterGodunov`](@ref) [Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite) +operator splitting using the reaction tangent as proposed in [OgiBalPer:2023:seats](@cite) +# Fields +- `σ_s::T`: steepness +- `σ_c::T`: offset in R axis +- `Δt_bounds::NTuple{2,T}`: lower and upper timestep length bounds +- `Rₙ₊₁::T`: updated maximal reaction magnitude +- `Rₙ::T`: previous reaction magnitude +""" mutable struct ReactionTangentController{T <: Real} <: AbstractTimeAdaptionAlgorithm const σ_s::T const σ_c::T @@ -10,20 +19,37 @@ mutable struct ReactionTangentController{T <: Real} <: AbstractTimeAdaptionAlgor Rₙ::T end +""" + get_next_dt(controller::ReactionTangentController) +Returns the next timestep length using [`ReactionTangentController`](@ref) calculated as +```math +\\sigma\\left(R_{\\max }\\right):=\\left(1.0-\\frac{1}{1+\\exp \\left(\\left(\\sigma_{\\mathrm{c}}-R_{\\max }\\right) \\cdot \\sigma_{\\mathrm{s}}\\right)}\\right) \\cdot\\left(\\Delta t_{\\max }-\\Delta t_{\\min }\\right)+\\Delta t_{\\min } +``` +""" function get_next_dt(controller::ReactionTangentController) @unpack σ_s, σ_c, Δt_bounds, Rₙ₊₁, Rₙ = controller R = max(Rₙ, Rₙ₊₁) (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] end +""" + AdaptiveOperatorSplittingAlgorithm{TOperatorSplittingAlg <: OS.AbstractOperatorSplittingAlgorithm, TTimeAdaptionAlgorithm <: AbstractTimeAdaptionAlgorithm} <: OS.AbstractOperatorSplittingAlgorithm +A generic operator splitting algorithm `operator_splitting_algorithm` with adaptive timestepping using the controller `timestep_controller`. +# Fields +- `operator_splitting_algorithm::TOperatorSplittingAlg`: steepness +- `timestep_controller::TTimeAdaptionAlgorithm`: offset in R axis +""" struct AdaptiveOperatorSplittingAlgorithm{TOperatorSplittingAlg <: OS.AbstractOperatorSplittingAlgorithm, TTimeAdaptionAlgorithm <: AbstractTimeAdaptionAlgorithm} <: OS.AbstractOperatorSplittingAlgorithm operator_splitting_algorithm::TOperatorSplittingAlg timestep_controller::TTimeAdaptionAlgorithm end -is_adaptive(::AdaptiveOperatorSplittingAlgorithm) = true - +@inline OS.is_adaptive(::AdaptiveOperatorSplittingAlgorithm) = true +""" + get_reaction_tangent(integrator::OS.OperatorSplittingIntegrator) +Returns the maximal reaction magnitude using the [`PointwiseODEFunction`](@ref) of an operator splitting integrator that uses [`LieTrotterGodunov`](@ref) [Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite). +""" @inline function get_reaction_tangent(integrator::OS.OperatorSplittingIntegrator) for subintegrator in integrator.subintegrators subintegrator.f isa PointwiseODEFunction || continue diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index 9a588aa8..a27fd08f 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -184,7 +184,16 @@ function (integrator::OperatorSplittingIntegrator)(tmp, t) linear_interpolation!(tmp, t, integrator.uprev, integrator.u, integrator.tprev, integrator.t) end +""" + update_controller!(::OperatorSplittingIntegrator) +Updates the controller using the current state of the integrator if the operator splitting algorithm is adaptive. +""" @inline update_controller!(::OperatorSplittingIntegrator) = nothing + +""" + update_dt!(::OperatorSplittingIntegrator) +Updates `_dt` of the integrator if the operator splitting algorithm is adaptive. +""" @inline update_dt!(::OperatorSplittingIntegrator) = nothing # helper functions for dealing with time-reversed integrators in the same way diff --git a/src/solver/operator_splitting/solver.jl b/src/solver/operator_splitting/solver.jl index 65a8a789..e8f358cc 100644 --- a/src/solver/operator_splitting/solver.jl +++ b/src/solver/operator_splitting/solver.jl @@ -9,6 +9,10 @@ struct LieTrotterGodunov{AlgTupleType} <: AbstractOperatorSplittingAlgorithm # transfer_algs::TransferTupleType # Tuple of transfer algorithms from the master solution into the individual ones end +""" + is_adaptive(::AbstractOperatorSplittingAlgorithm) +Returns whether the operator splitting algorithm is adaptive or not. +""" @inline is_adaptive(::AbstractOperatorSplittingAlgorithm) = false struct LieTrotterGodunovCache{uType, tmpType, iiType} <: AbstractOperatorSplittingCache From 03827ed2e5c4a5ecd03c1f9e09f1e49d5946f8a8 Mon Sep 17 00:00:00 2001 From: PotatoTuna Date: Wed, 14 Aug 2024 21:40:15 +0200 Subject: [PATCH 09/37] ugh zero default and one more test and commented tests because we need a way to find which operator in the split is B --- src/solver/adaptivity.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index 4888ff44..96e08141 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -19,6 +19,10 @@ mutable struct ReactionTangentController{T <: Real} <: AbstractTimeAdaptionAlgor Rₙ::T end +function ReactionTangentController(σ_s::T, σ_c::T, Δt_bounds::NTuple{2,T}) where {T <: Real} + return ReactionTangentController(σ_s, σ_c, Δt_bounds, 0.0, 0.0) +end + """ get_next_dt(controller::ReactionTangentController) Returns the next timestep length using [`ReactionTangentController`](@ref) calculated as From 645103050443af798e5a277372e95346669fd719 Mon Sep 17 00:00:00 2001 From: PotatoTuna Date: Wed, 14 Aug 2024 21:40:42 +0200 Subject: [PATCH 10/37] . --- test/test_integrators.jl | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 680d9f7d..75f65033 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -87,7 +87,7 @@ f3 = ODEFunction(ode3) fsplit2_inner = GenericSplitFunction((f3,f3), (f3dofs, f3dofs)) fsplit2_outer = GenericSplitFunction((f1,fsplit2_inner), (f1dofs, f2dofs)) for TimeStepperType in (LieTrotterGodunov,) - for controller in (Thunderbolt.ReactionTangentController(0.5, 1.0, (0.01, 0.3), 0.0, 0.0),) + for controller in (Thunderbolt.ReactionTangentController(0.5, 1.0, (0.01, 0.3)),) timestepper = TimeStepperType( (DummyForwardEuler(), DummyForwardEuler()) ) @@ -95,7 +95,7 @@ f3 = ODEFunction(ode3) timestepper_inner = TimeStepperType( (DummyForwardEuler(), DummyForwardEuler()) ) - timestepper_inner_adaptive = Thunderbolt.AdaptiveOperatorSplittingAlgorithm(timestepper_inner, controller) + timestepper_inner_adaptive = Thunderbolt.AdaptiveOperatorSplittingAlgorithm(timestepper_inner, controller) #TODO: Copy the controller instead timestepper2 = TimeStepperType( (DummyForwardEuler(), timestepper_inner) ) @@ -103,10 +103,10 @@ f3 = ODEFunction(ode3) for (tstepper1, tstepper_inner, tstepper2) in ( (timestepper, timestepper_inner, timestepper2), - (timestepper_adaptive, timestepper_inner_adaptive, timestepper2_adaptive) + # (timestepper_adaptive, timestepper_inner_adaptive, timestepper2_adaptive) ) # The remaining code works as usual. - integrator = DiffEqBase.init(prob, timestepper, dt=0.01, verbose=true) + integrator = DiffEqBase.init(prob, tstepper1, dt=0.01, verbose=true) DiffEqBase.solve!(integrator) ufinal = copy(integrator.u) @test ufinal ≉ u0 # Make sure the solve did something @@ -126,7 +126,7 @@ f3 = ODEFunction(ode3) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success prob2 = OperatorSplittingProblem(fsplit2_outer, u0, tspan) - integrator2 = DiffEqBase.init(prob2, timestepper2, dt=0.01, verbose=true) + integrator2 = DiffEqBase.init(prob2, tstepper2, dt=0.01, verbose=true) DiffEqBase.solve!(integrator2) DiffEqBase.reinit!(integrator2, u0; tspan) @@ -138,6 +138,11 @@ f3 = ODEFunction(ode3) DiffEqBase.solve!(integrator2) @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Success end + # integrator = DiffEqBase.init(prob, timestepper, dt=0.01, verbose=true) + # for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) end + # integrator_adaptive = DiffEqBase.init(prob, timestepper_adaptive, dt=0.01, verbose=true) + # for (u, t) in DiffEqBase.TimeChoiceIterator(integrator_adaptive, 0.0:5.0:100.0) end + # @test isapprox(integrator_adaptive.u, integrator.u, atol=1e-5) end end end From 7a39637b87e623a76f1ab2f7280eb1f5e99b76e3 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Thu, 15 Aug 2024 11:56:25 +0200 Subject: [PATCH 11/37] Now it doesn't allocate --- src/Thunderbolt.jl | 2 +- src/solver/adaptivity.jl | 12 ++++++++---- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/src/Thunderbolt.jl b/src/Thunderbolt.jl index d61c2e9e..b7fa909e 100644 --- a/src/Thunderbolt.jl +++ b/src/Thunderbolt.jl @@ -2,7 +2,7 @@ module Thunderbolt using TimerOutputs -import Unrolled: @unroll +import Unrolled: @unroll, unrolled_filter using Reexport, UnPack import LinearAlgebra: mul! diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index 96e08141..fc6d8354 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -55,15 +55,19 @@ end Returns the maximal reaction magnitude using the [`PointwiseODEFunction`](@ref) of an operator splitting integrator that uses [`LieTrotterGodunov`](@ref) [Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite). """ @inline function get_reaction_tangent(integrator::OS.OperatorSplittingIntegrator) - for subintegrator in integrator.subintegrators - subintegrator.f isa PointwiseODEFunction || continue + reaction_subintegrators = unrolled_filter(subintegrator -> subintegrator.f isa PointwiseODEFunction, integrator.subintegrators) + isempty(reaction_subintegrators) && @error "PointwiseODEFunction not found" + _get_reaction_tangent(reaction_subintegrators) +end + +@unroll function _get_reaction_tangent(subintegrators) + @unroll for subintegrator in subintegrators + #TODO: It should be either all the the same type or just one subintegrator after filtering, don't unroll? φₘidx = transmembranepotential_index(subintegrator.f.ode) return maximum(@view subintegrator.cache.dumat[:, φₘidx]) end - @error "PointwiseODEFunction not found" end - @inline function OS.update_controller!(integrator::OS.OperatorSplittingIntegrator{<:Any, <:AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov, <:ReactionTangentController}}) controller = integrator.alg.timestep_controller controller.Rₙ = controller.Rₙ₊₁ From 5d78f9162b55ecb0f00e2b409a1a5b3f9b053e87 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Thu, 15 Aug 2024 18:57:56 +0200 Subject: [PATCH 12/37] some formatting --- docs/src/api-reference/solver.md | 3 +- src/solver/adaptivity.jl | 49 ++++++++++++--------- src/solver/operator_splitting/integrator.jl | 36 ++++++++++----- src/solver/operator_splitting/solver.jl | 7 +-- 4 files changed, 55 insertions(+), 40 deletions(-) diff --git a/docs/src/api-reference/solver.md b/docs/src/api-reference/solver.md index f7f0e364..1174c9a4 100644 --- a/docs/src/api-reference/solver.md +++ b/docs/src/api-reference/solver.md @@ -39,10 +39,9 @@ Thunderbolt.OS.OperatorSplittingIntegrator ```@docs Thunderbolt.AdaptiveOperatorSplittingAlgorithm -Thunderbolt.OS.is_adaptive Thunderbolt.ReactionTangentController Thunderbolt.get_reaction_tangent -Thunderbolt.OS.update_controller! +Thunderbolt.OS.stepsize_controller! Thunderbolt.OS.update_dt! Thunderbolt.get_next_dt ``` diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index fc6d8354..3f170c89 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -23,32 +23,20 @@ function ReactionTangentController(σ_s::T, σ_c::T, Δt_bounds::NTuple{2,T}) wh return ReactionTangentController(σ_s, σ_c, Δt_bounds, 0.0, 0.0) end -""" - get_next_dt(controller::ReactionTangentController) -Returns the next timestep length using [`ReactionTangentController`](@ref) calculated as -```math -\\sigma\\left(R_{\\max }\\right):=\\left(1.0-\\frac{1}{1+\\exp \\left(\\left(\\sigma_{\\mathrm{c}}-R_{\\max }\\right) \\cdot \\sigma_{\\mathrm{s}}\\right)}\\right) \\cdot\\left(\\Delta t_{\\max }-\\Delta t_{\\min }\\right)+\\Delta t_{\\min } -``` -""" -function get_next_dt(controller::ReactionTangentController) - @unpack σ_s, σ_c, Δt_bounds, Rₙ₊₁, Rₙ = controller - R = max(Rₙ, Rₙ₊₁) - (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] -end """ AdaptiveOperatorSplittingAlgorithm{TOperatorSplittingAlg <: OS.AbstractOperatorSplittingAlgorithm, TTimeAdaptionAlgorithm <: AbstractTimeAdaptionAlgorithm} <: OS.AbstractOperatorSplittingAlgorithm -A generic operator splitting algorithm `operator_splitting_algorithm` with adaptive timestepping using the controller `timestep_controller`. +A generic operator splitting algorithm `operator_splitting_algorithm` with adaptive timestepping using the controller `controller`. # Fields - `operator_splitting_algorithm::TOperatorSplittingAlg`: steepness -- `timestep_controller::TTimeAdaptionAlgorithm`: offset in R axis +- `controller::TTimeAdaptionAlgorithm`: offset in R axis """ struct AdaptiveOperatorSplittingAlgorithm{TOperatorSplittingAlg <: OS.AbstractOperatorSplittingAlgorithm, TTimeAdaptionAlgorithm <: AbstractTimeAdaptionAlgorithm} <: OS.AbstractOperatorSplittingAlgorithm operator_splitting_algorithm::TOperatorSplittingAlg - timestep_controller::TTimeAdaptionAlgorithm + controller::TTimeAdaptionAlgorithm end -@inline OS.is_adaptive(::AdaptiveOperatorSplittingAlgorithm) = true +@inline DiffEqBase.isadaptive(::AdaptiveOperatorSplittingAlgorithm) = true """ get_reaction_tangent(integrator::OS.OperatorSplittingIntegrator) @@ -60,7 +48,7 @@ Returns the maximal reaction magnitude using the [`PointwiseODEFunction`](@ref) _get_reaction_tangent(reaction_subintegrators) end -@unroll function _get_reaction_tangent(subintegrators) +@inline @unroll function _get_reaction_tangent(subintegrators) @unroll for subintegrator in subintegrators #TODO: It should be either all the the same type or just one subintegrator after filtering, don't unroll? φₘidx = transmembranepotential_index(subintegrator.f.ode) @@ -68,19 +56,36 @@ end end end -@inline function OS.update_controller!(integrator::OS.OperatorSplittingIntegrator{<:Any, <:AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov, <:ReactionTangentController}}) - controller = integrator.alg.timestep_controller +@inline function OS.stepsize_controller!(integrator::OS.OperatorSplittingIntegrator, alg::AdaptiveOperatorSplittingAlgorithm) + OS.stepsize_controller!(integrator, alg.controller, alg) +end + +@inline function OS.step_accept_controller!(integrator::OS.OperatorSplittingIntegrator, alg::AdaptiveOperatorSplittingAlgorithm, q) + OS.step_accept_controller!(integrator, alg.controller, alg, q) +end + +@inline function OS.step_reject_controller!(integrator::OS.OperatorSplittingIntegrator, alg::AdaptiveOperatorSplittingAlgorithm, q) + OS.step_reject_controller!(integrator, alg.controller, alg, q) +end + +@inline function OS.stepsize_controller!(integrator::OS.OperatorSplittingIntegrator, controller::ReactionTangentController, alg::AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov}) + @unpack σ_s, σ_c, Δt_bounds, Rₙ₊₁, Rₙ = controller controller.Rₙ = controller.Rₙ₊₁ controller.Rₙ₊₁ = get_reaction_tangent(integrator) return nothing end -@inline function OS.update_dt!(integrator::OS.OperatorSplittingIntegrator{<:Any, <:AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov, <:ReactionTangentController}}) - controller = integrator.alg.timestep_controller - integrator._dt = get_next_dt(controller) +@inline function OS.step_accept_controller!(integrator::OS.OperatorSplittingIntegrator, controller::ReactionTangentController, alg::AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov}, q) + @unpack σ_s, σ_c, Δt_bounds, Rₙ₊₁, Rₙ = controller + R = max(Rₙ, Rₙ₊₁) + integrator._dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] return nothing end +@inline function OS.step_reject_controller!(integrator::OS.OperatorSplittingIntegrator, controller::ReactionTangentController, alg::AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov}, q) + return nothing # Do nothing +end + # Dispatch for outer construction function OS.init_cache(prob::OS.OperatorSplittingProblem, alg::AdaptiveOperatorSplittingAlgorithm; dt, kwargs...) # TODO OS.init_cache(prob, alg.operator_splitting_algorithm;dt = dt, kwargs...) diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index a27fd08f..9755d9c0 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -54,7 +54,7 @@ function DiffEqBase.__init( callback = nothing, advance_to_tstop = false, save_func = (u, t) -> copy(u), # custom kwarg - dtchangeable = is_adaptive(alg), # custom kwarg + dtchangeable = DiffEqBase.isadaptive(alg), # custom kwarg kwargs..., ) (; u0, p) = prob @@ -166,8 +166,6 @@ function DiffEqBase.step!(integrator::OperatorSplittingIntegrator, dt, stop_at_t end end - - # TimeChoiceIterator API @inline function DiffEqBase.get_tmp_cache(integrator::OperatorSplittingIntegrator) DiffEqBase.get_tmp_cache(integrator, integrator.alg, integrator.cache) @@ -185,16 +183,34 @@ function (integrator::OperatorSplittingIntegrator)(tmp, t) end """ - update_controller!(::OperatorSplittingIntegrator) + stepsize_controller!(::OperatorSplittingIntegrator) Updates the controller using the current state of the integrator if the operator splitting algorithm is adaptive. """ -@inline update_controller!(::OperatorSplittingIntegrator) = nothing +@inline function stepsize_controller!(integrator::OperatorSplittingIntegrator) + algorithm = integrator.alg + DiffEqBase.isadaptive(algorithm) || return nothing + stepsize_controller!(integrator, algorithm) +end + +""" + step_accept_controller!(::OperatorSplittingIntegrator) +Updates `_dt` of the integrator if the step is accepted and the operator splitting algorithm is adaptive. +""" +@inline function step_accept_controller!(integrator::OperatorSplittingIntegrator) + algorithm = integrator.alg + DiffEqBase.isadaptive(algorithm) || return nothing + step_accept_controller!(integrator, algorithm, nothing) +end """ - update_dt!(::OperatorSplittingIntegrator) -Updates `_dt` of the integrator if the operator splitting algorithm is adaptive. + step_reject_controller!(::OperatorSplittingIntegrator) +Updates `_dt` of the integrator if the step is rejected and the the operator splitting algorithm is adaptive. """ -@inline update_dt!(::OperatorSplittingIntegrator) = nothing +@inline function step_reject_controller!(integrator::OperatorSplittingIntegrator) + algorithm = integrator.alg + DiffEqBase.isadaptive(algorithm) || return nothing + step_reject_controller!(integrator, algorithm, nothing) +end # helper functions for dealing with time-reversed integrators in the same way # that OrdinaryDiffEq.jl does @@ -244,7 +260,7 @@ function __step!(integrator) # Solve inner problems advance_solution_to!(integrator, tnext) - update_controller!(integrator) + stepsize_controller!(integrator) # Update integrator # increment t by dt, rounding to the first tstop if that is roughly @@ -255,7 +271,7 @@ function __step!(integrator) integrator.tprev = integrator.t integrator.t = !isempty(tstops) && abs(first(tstops) - tnext) < max_t_error ? first(tstops) : tnext - update_dt!(integrator) + step_accept_controller!(integrator) # remove tstops that were just reached while !isempty(tstops) && reached_tstop(integrator, first(tstops)) diff --git a/src/solver/operator_splitting/solver.jl b/src/solver/operator_splitting/solver.jl index e8f358cc..69d9640e 100644 --- a/src/solver/operator_splitting/solver.jl +++ b/src/solver/operator_splitting/solver.jl @@ -1,4 +1,3 @@ - # Lie-Trotter-Godunov Splitting Implementation """ LieTrotterGodunov <: AbstractOperatorSplittingAlgorithm @@ -9,11 +8,7 @@ struct LieTrotterGodunov{AlgTupleType} <: AbstractOperatorSplittingAlgorithm # transfer_algs::TransferTupleType # Tuple of transfer algorithms from the master solution into the individual ones end -""" - is_adaptive(::AbstractOperatorSplittingAlgorithm) -Returns whether the operator splitting algorithm is adaptive or not. -""" -@inline is_adaptive(::AbstractOperatorSplittingAlgorithm) = false +@inline DiffEqBase.isadaptive(::AbstractOperatorSplittingAlgorithm) = false struct LieTrotterGodunovCache{uType, tmpType, iiType} <: AbstractOperatorSplittingCache u::uType From f909309644b005cf93be069ae2a1c6c83f795f3e Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Thu, 15 Aug 2024 19:00:41 +0200 Subject: [PATCH 13/37] example formatting --- examples/conduction-velocity-benchmark.jl | 31 ++++++++++++----------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/examples/conduction-velocity-benchmark.jl b/examples/conduction-velocity-benchmark.jl index 723cacdb..41384c09 100644 --- a/examples/conduction-velocity-benchmark.jl +++ b/examples/conduction-velocity-benchmark.jl @@ -63,21 +63,22 @@ steady_state_initializer!(u₀, odeform) # io = ParaViewWriter("spiral-wave-test") - -_timestepper = OS.LieTrotterGodunov(( - BackwardEulerSolver( - solution_vector_type=Vector{Float32}, - system_matrix_type=Thunderbolt.ThreadedSparseMatrixCSR{Float32, Int32}, - inner_solver=LinearSolve.KrylovJL_CG(atol=1.0f-6, rtol=1.0f-5), - ), - AdaptiveForwardEulerSubstepper( - solution_vector_type=Vector{Float32}, - reaction_threshold=0.1f0, - ), -)) -controller = Thunderbolt.ReactionTangentController(0.5, 1.0, (0.01, 0.3), 0.0, 0.0) - -timestepper = Thunderbolt.AdaptiveOperatorSplittingAlgorithm(_timestepper, controller) +timestepper = Thunderbolt.AdaptiveOperatorSplittingAlgorithm( + OS.LieTrotterGodunov(( + BackwardEulerSolver( + solution_vector_type=Vector{Float32}, + system_matrix_type=Thunderbolt.ThreadedSparseMatrixCSR{Float32, Int32}, + inner_solver=LinearSolve.KrylovJL_CG(atol=1.0f-6, rtol=1.0f-5), + ), + AdaptiveForwardEulerSubstepper( + solution_vector_type=Vector{Float32}, + reaction_threshold=0.1f0, + ) + )), + Thunderbolt.ReactionTangentController( + 0.5, 1.0, (0.01, 0.3) + ) +) problem = OS.OperatorSplittingProblem(odeform, u₀, tspan) From ca0cd7c5a97e094a51d7291a88e9f2e2e1473c74 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Thu, 15 Aug 2024 19:03:24 +0200 Subject: [PATCH 14/37] retcode --- test/test_integrators.jl | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 75f65033..082ead09 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -107,34 +107,46 @@ f3 = ODEFunction(ode3) ) # The remaining code works as usual. integrator = DiffEqBase.init(prob, tstepper1, dt=0.01, verbose=true) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default DiffEqBase.solve!(integrator) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success ufinal = copy(integrator.u) @test ufinal ≉ u0 # Make sure the solve did something DiffEqBase.reinit!(integrator, u0; tspan) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) end + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success @test isapprox(ufinal, integrator.u, atol=1e-8) DiffEqBase.reinit!(integrator, u0; tspan) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default for (uprev, tprev, u, t) in DiffEqBase.intervals(integrator) end + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success @test isapprox(ufinal, integrator.u, atol=1e-8) DiffEqBase.reinit!(integrator, u0; tspan) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default DiffEqBase.solve!(integrator) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success prob2 = OperatorSplittingProblem(fsplit2_outer, u0, tspan) integrator2 = DiffEqBase.init(prob2, tstepper2, dt=0.01, verbose=true) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default DiffEqBase.solve!(integrator2) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success DiffEqBase.reinit!(integrator2, u0; tspan) + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default for (u, t) in DiffEqBase.TimeChoiceIterator(integrator2, 0.0:5.0:100.0) end + @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success @test isapprox(ufinal, integrator2.u, atol=1e-8) DiffEqBase.reinit!(integrator2, u0; tspan) + @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Default DiffEqBase.solve!(integrator2) @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Success end From e5dbc21a183a5439e4a89cacb514df8d603baeb3 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Thu, 15 Aug 2024 19:53:31 +0200 Subject: [PATCH 15/37] add \pi so it doesn't align with tstops hehe and retcode oopsie --- src/solver/adaptivity.jl | 2 +- src/solver/operator_splitting/integrator.jl | 6 +++++- test/test_integrators.jl | 15 +++++++-------- 3 files changed, 13 insertions(+), 10 deletions(-) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index 3f170c89..6acfbae5 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -78,7 +78,7 @@ end @inline function OS.step_accept_controller!(integrator::OS.OperatorSplittingIntegrator, controller::ReactionTangentController, alg::AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov}, q) @unpack σ_s, σ_c, Δt_bounds, Rₙ₊₁, Rₙ = controller R = max(Rₙ, Rₙ₊₁) - integrator._dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] + integrator.dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] return nothing end diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index 9755d9c0..ae2e7e6d 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -111,6 +111,7 @@ function DiffEqBase.reinit!( tstops = integrator._tstops, saveat = integrator._saveat, reinit_callbacks = true, + reinit_retcode = true ) (t0,tf) = tspan integrator.u .= u0 @@ -126,6 +127,9 @@ function DiffEqBase.reinit!( saving_callback = integrator.callback.discrete_callbacks[end] DiffEqBase.initialize!(saving_callback, u0, t0, integrator) end + if reinit_retcode + integrator.sol = DiffEqBase.solution_new_retcode(integrator.sol, DiffEqBase.ReturnCode.Default) + end end # called by DiffEqBase.solve @@ -246,7 +250,7 @@ end function __step!(integrator) (; dtchangeable, tstops) = integrator - _dt = DiffEqBase.get_dt(integrator) + _dt = DiffEqBase.isadaptive(integrator.alg) ? DiffEqBase.get_dt(integrator) : integrator.dt # update dt before incrementing u; if dt is changeable and there is # a tstop within dt, reduce dt to tstop - t diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 082ead09..593e280b 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -106,7 +106,7 @@ f3 = ODEFunction(ode3) # (timestepper_adaptive, timestepper_inner_adaptive, timestepper2_adaptive) ) # The remaining code works as usual. - integrator = DiffEqBase.init(prob, tstepper1, dt=0.01, verbose=true) + integrator = DiffEqBase.init(prob, tstepper1, dt=0.01π, verbose=true) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default DiffEqBase.solve!(integrator) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success @@ -117,14 +117,12 @@ f3 = ODEFunction(ode3) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) end - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success @test isapprox(ufinal, integrator.u, atol=1e-8) DiffEqBase.reinit!(integrator, u0; tspan) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default for (uprev, tprev, u, t) in DiffEqBase.intervals(integrator) end - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success @test isapprox(ufinal, integrator.u, atol=1e-8) DiffEqBase.reinit!(integrator, u0; tspan) @@ -133,17 +131,18 @@ f3 = ODEFunction(ode3) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success prob2 = OperatorSplittingProblem(fsplit2_outer, u0, tspan) - integrator2 = DiffEqBase.init(prob2, tstepper2, dt=0.01, verbose=true) - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default + integrator2 = DiffEqBase.init(prob2, tstepper2, dt=0.02π, verbose=true) + @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Default DiffEqBase.solve!(integrator2) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success + ufinal2 = copy(integrator2.u) + @test ufinal2 ≉ u0 # Make sure the solve did something DiffEqBase.reinit!(integrator2, u0; tspan) - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default + @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Default for (u, t) in DiffEqBase.TimeChoiceIterator(integrator2, 0.0:5.0:100.0) end - @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success - @test isapprox(ufinal, integrator2.u, atol=1e-8) + @test isapprox(ufinal2, integrator2.u, atol=1e-8) DiffEqBase.reinit!(integrator2, u0; tspan) @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Default From d2af0f94e9047941fc93de177d21cb288736430a Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Thu, 15 Aug 2024 20:17:12 +0200 Subject: [PATCH 16/37] nans but commented out --- test/test_integrators.jl | 32 +++++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 593e280b..56c1265f 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -76,16 +76,23 @@ end f3 = ODEFunction(ode3) # The time stepper carries the individual solver information. +# Note that we define the dof indices w.r.t the parent function. +# Hence the indices for `fsplit2_inner` are. +f1dofs = [1,2,3] +f2dofs = [1,3] +f3dofs = [1,3] +fsplit2_inner = GenericSplitFunction((f3,f3), (f3dofs, f3dofs)) +fsplit2_outer = GenericSplitFunction((f1,fsplit2_inner), (f1dofs, f2dofs)) + +function ode_NaN(du, u, p, t) + du[1] = NaN + du[2] = 0.01u[1] +end +f_NaN = ODEFunction(ode_NaN) +f_NaN_dofs = f3dofs +fsplit_NaN = GenericSplitFunction((f1,f_NaN), (f1dofs, f_NaN_dofs)) + @testset "OperatorSplitting" begin - - - # Note that we define the dof indices w.r.t the parent function. - # Hence the indices for `fsplit2_inner` are. - f1dofs = [1,2,3] - f2dofs = [1,3] - f3dofs = [1,3] - fsplit2_inner = GenericSplitFunction((f3,f3), (f3dofs, f3dofs)) - fsplit2_outer = GenericSplitFunction((f1,fsplit2_inner), (f1dofs, f2dofs)) for TimeStepperType in (LieTrotterGodunov,) for controller in (Thunderbolt.ReactionTangentController(0.5, 1.0, (0.01, 0.3)),) timestepper = TimeStepperType( @@ -148,6 +155,13 @@ f3 = ODEFunction(ode3) @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Default DiffEqBase.solve!(integrator2) @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Success + @testset "NaNs" begin + # prob_NaN = OperatorSplittingProblem(fsplit_NaN, u0, tspan) + # integrator_NaN = DiffEqBase.init(prob, tstepper1, dt=0.01π, verbose=true) + # @test integrator_NaN.sol.retcode == DiffEqBase.ReturnCode.Default + # DiffEqBase.solve!(integrator_NaN) + # @test integrator_NaN.sol.retcode == DiffEqBase.ReturnCode.Failure + end end # integrator = DiffEqBase.init(prob, timestepper, dt=0.01, verbose=true) # for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) end From 071ecd4d9c05e9d0512522f1be11fabf9cf32e54 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Thu, 15 Aug 2024 20:24:03 +0200 Subject: [PATCH 17/37] assume only one problem --- src/solver/adaptivity.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index 6acfbae5..8dcc4266 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -41,6 +41,7 @@ end """ get_reaction_tangent(integrator::OS.OperatorSplittingIntegrator) Returns the maximal reaction magnitude using the [`PointwiseODEFunction`](@ref) of an operator splitting integrator that uses [`LieTrotterGodunov`](@ref) [Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite). +It is assumed that the problem containing the reaction tangent is a [`PointwiseODEFunction`](@ref). """ @inline function get_reaction_tangent(integrator::OS.OperatorSplittingIntegrator) reaction_subintegrators = unrolled_filter(subintegrator -> subintegrator.f isa PointwiseODEFunction, integrator.subintegrators) @@ -48,12 +49,13 @@ Returns the maximal reaction magnitude using the [`PointwiseODEFunction`](@ref) _get_reaction_tangent(reaction_subintegrators) end -@inline @unroll function _get_reaction_tangent(subintegrators) - @unroll for subintegrator in subintegrators +@inline #=@unroll=# function _get_reaction_tangent(subintegrators) + subintegrator = subintegrators[1] + # @unroll for subintegrator in subintegrators #TODO: It should be either all the the same type or just one subintegrator after filtering, don't unroll? - φₘidx = transmembranepotential_index(subintegrator.f.ode) - return maximum(@view subintegrator.cache.dumat[:, φₘidx]) - end + φₘidx = transmembranepotential_index(subintegrator.f.ode) + return maximum(@view subintegrator.cache.dumat[:, φₘidx]) + # end end @inline function OS.stepsize_controller!(integrator::OS.OperatorSplittingIntegrator, alg::AdaptiveOperatorSplittingAlgorithm) From 921ab740ba9d4868da3b580879ce9a543c0504fd Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Thu, 15 Aug 2024 20:24:40 +0200 Subject: [PATCH 18/37] docs oopsie --- docs/src/api-reference/solver.md | 3 --- 1 file changed, 3 deletions(-) diff --git a/docs/src/api-reference/solver.md b/docs/src/api-reference/solver.md index 1174c9a4..a586bd80 100644 --- a/docs/src/api-reference/solver.md +++ b/docs/src/api-reference/solver.md @@ -41,7 +41,4 @@ Thunderbolt.OS.OperatorSplittingIntegrator Thunderbolt.AdaptiveOperatorSplittingAlgorithm Thunderbolt.ReactionTangentController Thunderbolt.get_reaction_tangent -Thunderbolt.OS.stepsize_controller! -Thunderbolt.OS.update_dt! -Thunderbolt.get_next_dt ``` From c98c348a597dc319403f7b9aaa4c19c03e784120 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 16 Aug 2024 10:36:36 +0200 Subject: [PATCH 19/37] adaptive test UwU --- test/integration/test_electrophysiology.jl | 47 +++++++++++++++++----- 1 file changed, 36 insertions(+), 11 deletions(-) diff --git a/test/integration/test_electrophysiology.jl b/test/integration/test_electrophysiology.jl index 57f590a6..264d6d98 100644 --- a/test/integration/test_electrophysiology.jl +++ b/test/integration/test_electrophysiology.jl @@ -24,7 +24,7 @@ using Thunderbolt end end - function solve_waveprop(mesh, coeff, subdomains) + function solve_waveprop(mesh, coeff, subdomains, isadaptive = false) cs = CoordinateSystemCoefficient(CartesianCoordinateSystem(mesh)) model = MonodomainModel( ConstantCoefficient(1.0), @@ -45,10 +45,18 @@ using Thunderbolt mesh ) - timestepper = LieTrotterGodunov(( + _timestepper = LieTrotterGodunov(( BackwardEulerSolver(), ForwardEulerCellSolver() )) + if isadaptive + timestepper = Thunderbolt.AdaptiveOperatorSplittingAlgorithm( + _timestepper, + Thunderbolt.ReactionTangentController(0.5, 1.0, (0.01, 0.3)) + ) + else + timestepper = _timestepper + end u₀ = zeros(Float64, OS.function_size(odeform)) simple_initializer!(u₀, odeform) @@ -60,26 +68,43 @@ using Thunderbolt DiffEqBase.solve!(integrator) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success @test integrator.u ≉ u₀ + return integrator.u end mesh = generate_mesh(Hexahedron, (4, 4, 4), Vec{3}((0.0,0.0,0.0)), Vec{3}((1.0,1.0,1.0))) coeff = ConstantCoefficient(SymmetricTensor{2,3,Float64}((4.5e-5, 0, 0, 2.0e-5, 0, 1.0e-5))) - solve_waveprop(mesh, coeff, [""]) + u = solve_waveprop(mesh, coeff, [""]) + u_adaptive = solve_waveprop(mesh, coeff, [""], true) + @test u ≈ u_adaptive mesh = generate_ideal_lv_mesh(4,1,1) coeff = ConstantCoefficient(SymmetricTensor{2,3,Float64}((4.5e-5, 0, 0, 2.0e-5, 0, 1.0e-5))) - solve_waveprop(mesh, coeff, [""]) - + u = solve_waveprop(mesh, coeff, [""]) + u_adaptive = solve_waveprop(mesh, coeff, [""], true) + @test u ≈ u_adaptive + mesh = to_mesh(generate_mixed_grid_2D()) coeff = ConstantCoefficient(SymmetricTensor{2,2,Float64}((4.5e-5, 0, 2.0e-5))) - solve_waveprop(mesh, coeff, ["Pacemaker", "Myocardium"]) - solve_waveprop(mesh, coeff, ["Pacemaker"]) - solve_waveprop(mesh, coeff, ["Myocardium"]) + u = solve_waveprop(mesh, coeff, ["Pacemaker", "Myocardium"]) + u_adaptive = solve_waveprop(mesh, coeff, ["Pacemaker", "Myocardium"], true) + @test u ≈ u_adaptive + u = solve_waveprop(mesh, coeff, ["Pacemaker"]) + u_adaptive = solve_waveprop(mesh, coeff, ["Pacemaker"], true) + @test u ≈ u_adaptive + u = solve_waveprop(mesh, coeff, ["Myocardium"]) + u_adaptive = solve_waveprop(mesh, coeff, ["Myocardium"], true) + @test u ≈ u_adaptive mesh = to_mesh(generate_mixed_dimensional_grid_3D()) coeff = ConstantCoefficient(SymmetricTensor{2,3,Float64}((4.5e-5, 0, 0, 2.0e-5, 0, 1.0e-5))) - solve_waveprop(mesh, coeff, ["Ventricle"]) + u = solve_waveprop(mesh, coeff, ["Ventricle"]) + u_adaptive = solve_waveprop(mesh, coeff, ["Ventricle"]) + @test u ≈ u_adaptive coeff = ConstantCoefficient(SymmetricTensor{2,3,Float64}((5e-5, 0, 0, 5e-5, 0, 5e-5))) - solve_waveprop(mesh, coeff, ["Purkinje"]) - solve_waveprop(mesh, coeff, ["Ventricle", "Purkinje"]) + u = solve_waveprop(mesh, coeff, ["Purkinje"]) + u_adaptive = solve_waveprop(mesh, coeff, ["Purkinje"]) + @test u ≈ u_adaptive + u = solve_waveprop(mesh, coeff, ["Ventricle", "Purkinje"]) + u_adaptive = solve_waveprop(mesh, coeff, ["Ventricle", "Purkinje"]) + @test u ≈ u_adaptive end From b5b41957212fbf81aa8f86ac0f96075bc1578ec1 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 19 Aug 2024 15:43:59 +0200 Subject: [PATCH 20/37] Better? --- examples/conduction-velocity-benchmark.jl | 6 +- src/solver/adaptivity.jl | 86 ++++++++++----------- src/solver/operator_splitting/integrator.jl | 11 ++- src/solver/operator_splitting/solver.jl | 3 + test/integration/test_electrophysiology.jl | 72 ++++++++--------- test/test_integrators.jl | 8 +- 6 files changed, 89 insertions(+), 97 deletions(-) diff --git a/examples/conduction-velocity-benchmark.jl b/examples/conduction-velocity-benchmark.jl index 41384c09..f3a9a619 100644 --- a/examples/conduction-velocity-benchmark.jl +++ b/examples/conduction-velocity-benchmark.jl @@ -63,7 +63,7 @@ steady_state_initializer!(u₀, odeform) # io = ParaViewWriter("spiral-wave-test") -timestepper = Thunderbolt.AdaptiveOperatorSplittingAlgorithm( +timestepper = Thunderbolt.ReactionTangentController( OS.LieTrotterGodunov(( BackwardEulerSolver( solution_vector_type=Vector{Float32}, @@ -75,9 +75,7 @@ timestepper = Thunderbolt.AdaptiveOperatorSplittingAlgorithm( reaction_threshold=0.1f0, ) )), - Thunderbolt.ReactionTangentController( - 0.5, 1.0, (0.01, 0.3) - ) + 0.5, 1.0, (0.01, 0.3) ) problem = OS.OperatorSplittingProblem(odeform, u₀, tspan) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index 8dcc4266..f0469c58 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -1,42 +1,36 @@ abstract type AbstractTimeAdaptionAlgorithm end """ - ReactionTangentController{T <: Real} <: AbstractTimeAdaptionAlgorithm + ReactionTangentController{T <: Real} <: OS.AbstractOperatorSplittingAlgorithm A timestep length controller for [`LieTrotterGodunov`](@ref) [Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite) operator splitting using the reaction tangent as proposed in [OgiBalPer:2023:seats](@cite) # Fields - `σ_s::T`: steepness - `σ_c::T`: offset in R axis - `Δt_bounds::NTuple{2,T}`: lower and upper timestep length bounds -- `Rₙ₊₁::T`: updated maximal reaction magnitude -- `Rₙ::T`: previous reaction magnitude """ -mutable struct ReactionTangentController{T <: Real} <: AbstractTimeAdaptionAlgorithm - const σ_s::T - const σ_c::T - const Δt_bounds::NTuple{2,T} - Rₙ₊₁::T - Rₙ::T +struct ReactionTangentController{LTG <: OS.LieTrotterGodunov, T <: Real} <: OS.AbstractOperatorSplittingAlgorithm + ltg::LTG + σ_s::T + σ_c::T + Δt_bounds::NTuple{2,T} end -function ReactionTangentController(σ_s::T, σ_c::T, Δt_bounds::NTuple{2,T}) where {T <: Real} - return ReactionTangentController(σ_s, σ_c, Δt_bounds, 0.0, 0.0) +mutable struct ReactionTangentControllerCache{T <: Real, LTGCache <: OS.LieTrotterGodunovCache} <: OS.AbstractOperatorSplittingCache + const ltg_cache::LTGCache #It has Arrays so it can be const? + Rₙ₊₁::T + Rₙ::T end +@inline OS.get_u(cache::ReactionTangentControllerCache) = OS.get_u(cache.ltg_cache) +@inline OS.get_uprev(cache::ReactionTangentControllerCache) = OS.get_uprev(cache.ltg_cache) +@inline DiffEqBase.get_tmp_cache(integrator::OS.OperatorSplittingIntegrator, alg::OS.AbstractOperatorSplittingAlgorithm, cache::ReactionTangentControllerCache) = DiffEqBase.get_tmp_cache(integrator, alg, cache.ltg_cache) -""" - AdaptiveOperatorSplittingAlgorithm{TOperatorSplittingAlg <: OS.AbstractOperatorSplittingAlgorithm, TTimeAdaptionAlgorithm <: AbstractTimeAdaptionAlgorithm} <: OS.AbstractOperatorSplittingAlgorithm -A generic operator splitting algorithm `operator_splitting_algorithm` with adaptive timestepping using the controller `controller`. -# Fields -- `operator_splitting_algorithm::TOperatorSplittingAlg`: steepness -- `controller::TTimeAdaptionAlgorithm`: offset in R axis -""" -struct AdaptiveOperatorSplittingAlgorithm{TOperatorSplittingAlg <: OS.AbstractOperatorSplittingAlgorithm, TTimeAdaptionAlgorithm <: AbstractTimeAdaptionAlgorithm} <: OS.AbstractOperatorSplittingAlgorithm - operator_splitting_algorithm::TOperatorSplittingAlg - controller::TTimeAdaptionAlgorithm +@inline function OS.advance_solution_to!(subintegrators::Tuple, cache::ReactionTangentControllerCache, tnext) + OS.advance_solution_to!(subintegrators, cache.ltg_cache, tnext) end -@inline DiffEqBase.isadaptive(::AdaptiveOperatorSplittingAlgorithm) = true +@inline DiffEqBase.isadaptive(::ReactionTangentController) = true """ get_reaction_tangent(integrator::OS.OperatorSplittingIntegrator) @@ -58,42 +52,42 @@ end # end end -@inline function OS.stepsize_controller!(integrator::OS.OperatorSplittingIntegrator, alg::AdaptiveOperatorSplittingAlgorithm) - OS.stepsize_controller!(integrator, alg.controller, alg) -end - -@inline function OS.step_accept_controller!(integrator::OS.OperatorSplittingIntegrator, alg::AdaptiveOperatorSplittingAlgorithm, q) - OS.step_accept_controller!(integrator, alg.controller, alg, q) -end - -@inline function OS.step_reject_controller!(integrator::OS.OperatorSplittingIntegrator, alg::AdaptiveOperatorSplittingAlgorithm, q) - OS.step_reject_controller!(integrator, alg.controller, alg, q) -end - -@inline function OS.stepsize_controller!(integrator::OS.OperatorSplittingIntegrator, controller::ReactionTangentController, alg::AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov}) - @unpack σ_s, σ_c, Δt_bounds, Rₙ₊₁, Rₙ = controller - controller.Rₙ = controller.Rₙ₊₁ - controller.Rₙ₊₁ = get_reaction_tangent(integrator) +@inline function OS.stepsize_controller!(integrator::OS.OperatorSplittingIntegrator, alg::ReactionTangentController) + integrator.cache.Rₙ = integrator.cache.Rₙ₊₁ + integrator.cache.Rₙ₊₁ = get_reaction_tangent(integrator) return nothing end -@inline function OS.step_accept_controller!(integrator::OS.OperatorSplittingIntegrator, controller::ReactionTangentController, alg::AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov}, q) - @unpack σ_s, σ_c, Δt_bounds, Rₙ₊₁, Rₙ = controller +@inline function OS.step_accept_controller!(integrator::OS.OperatorSplittingIntegrator, alg::ReactionTangentController, q) + @unpack Rₙ₊₁, Rₙ = integrator.cache + @unpack σ_s, σ_c, Δt_bounds = alg R = max(Rₙ, Rₙ₊₁) - integrator.dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] + integrator._dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] return nothing end -@inline function OS.step_reject_controller!(integrator::OS.OperatorSplittingIntegrator, controller::ReactionTangentController, alg::AdaptiveOperatorSplittingAlgorithm{<:OS.LieTrotterGodunov}, q) +@inline function OS.step_reject_controller!(integrator::OS.OperatorSplittingIntegrator, alg::ReactionTangentController, q) return nothing # Do nothing end # Dispatch for outer construction -function OS.init_cache(prob::OS.OperatorSplittingProblem, alg::AdaptiveOperatorSplittingAlgorithm; dt, kwargs...) # TODO - OS.init_cache(prob, alg.operator_splitting_algorithm;dt = dt, kwargs...) +function OS.init_cache(prob::OS.OperatorSplittingProblem, alg::ReactionTangentController; dt, kwargs...) # TODO + @unpack f = prob + @assert f isa GenericSplitFunction + + u = copy(prob.u0) + uprev = copy(prob.u0) + + # Build inner integrator + return OS.construct_inner_cache(f, alg, u, uprev) end # Dispatch for recursive construction -function OS.construct_inner_cache(f::OS.AbstractOperatorSplitFunction, alg::AdaptiveOperatorSplittingAlgorithm, u::AbstractArray, uprev::AbstractArray) - OS.construct_inner_cache(f, alg.operator_splitting_algorithm, u, uprev) +function OS.construct_inner_cache(f::OS.AbstractOperatorSplitFunction, alg::ReactionTangentController, u::AbstractArray{T}, uprev::AbstractArray) where T <: Number + ltg_cache = OS.construct_inner_cache(f, alg.ltg, u, uprev) + return ReactionTangentControllerCache(ltg_cache, zero(T), zero(T)) +end + +function OS.build_subintegrators_recursive(f::GenericSplitFunction, synchronizers::Tuple, p::Tuple, cache::ReactionTangentControllerCache, u::AbstractArray, uprev::AbstractArray, t, dt, dof_range, uparent, tstops, _tstops, saveat, _saveat) + OS.build_subintegrators_recursive(f, synchronizers, p, cache.ltg_cache, u, uprev, t, dt, dof_range, uparent, tstops, _tstops, saveat, _saveat) end \ No newline at end of file diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index ae2e7e6d..8b5cfc70 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -73,14 +73,17 @@ function DiffEqBase.__init( callback = DiffEqBase.CallbackSet(callback) cache = init_cache(prob, alg; dt, kwargs...) + + u = get_u(cache) + uprev = get_uprev(cache) - subintegrators = build_subintegrators_recursive(prob.f, prob.f.synchronizers, p, cache, cache.u, cache.uprev, t0, dt, 1:length(u0), cache.u, tstops, _tstops, saveat, _saveat) + subintegrators = build_subintegrators_recursive(prob.f, prob.f.synchronizers, p, cache, u, uprev, t0, dt, 1:length(u0), u, tstops, _tstops, saveat, _saveat) integrator = OperatorSplittingIntegrator( prob.f, alg, - cache.u, - cache.uprev, + u, + uprev, p, t0, copy(t0), @@ -250,7 +253,7 @@ end function __step!(integrator) (; dtchangeable, tstops) = integrator - _dt = DiffEqBase.isadaptive(integrator.alg) ? DiffEqBase.get_dt(integrator) : integrator.dt + _dt = DiffEqBase.get_dt(integrator) # update dt before incrementing u; if dt is changeable and there is # a tstop within dt, reduce dt to tstop - t diff --git a/src/solver/operator_splitting/solver.jl b/src/solver/operator_splitting/solver.jl index 69d9640e..529fcaf3 100644 --- a/src/solver/operator_splitting/solver.jl +++ b/src/solver/operator_splitting/solver.jl @@ -18,6 +18,9 @@ struct LieTrotterGodunovCache{uType, tmpType, iiType} <: AbstractOperatorSplitti inner_caches::iiType end +get_u(cache#=::LieTrotterGodunovCache=#) = cache.u +get_uprev(cache#=::LieTrotterGodunovCache=#) = cache.uprev + # Dispatch for outer construction function init_cache(prob::OperatorSplittingProblem, alg::LieTrotterGodunov; dt, kwargs...) # TODO @unpack f = prob diff --git a/test/integration/test_electrophysiology.jl b/test/integration/test_electrophysiology.jl index 264d6d98..ac5e8163 100644 --- a/test/integration/test_electrophysiology.jl +++ b/test/integration/test_electrophysiology.jl @@ -24,7 +24,7 @@ using Thunderbolt end end - function solve_waveprop(mesh, coeff, subdomains, isadaptive = false) + function solve_waveprop(mesh, coeff, subdomains, timestepper) cs = CoordinateSystemCoefficient(CartesianCoordinateSystem(mesh)) model = MonodomainModel( ConstantCoefficient(1.0), @@ -45,19 +45,6 @@ using Thunderbolt mesh ) - _timestepper = LieTrotterGodunov(( - BackwardEulerSolver(), - ForwardEulerCellSolver() - )) - if isadaptive - timestepper = Thunderbolt.AdaptiveOperatorSplittingAlgorithm( - _timestepper, - Thunderbolt.ReactionTangentController(0.5, 1.0, (0.01, 0.3)) - ) - else - timestepper = _timestepper - end - u₀ = zeros(Float64, OS.function_size(odeform)) simple_initializer!(u₀, odeform) @@ -71,40 +58,49 @@ using Thunderbolt return integrator.u end + timestepper = LieTrotterGodunov(( + BackwardEulerSolver(), + ForwardEulerCellSolver() + )) + timestepper_adaptive = Thunderbolt.ReactionTangentController( + timestepper, + 0.5, 1.0, (0.98, 1.02) + ) + mesh = generate_mesh(Hexahedron, (4, 4, 4), Vec{3}((0.0,0.0,0.0)), Vec{3}((1.0,1.0,1.0))) coeff = ConstantCoefficient(SymmetricTensor{2,3,Float64}((4.5e-5, 0, 0, 2.0e-5, 0, 1.0e-5))) - u = solve_waveprop(mesh, coeff, [""]) - u_adaptive = solve_waveprop(mesh, coeff, [""], true) - @test u ≈ u_adaptive + u = solve_waveprop(mesh, coeff, [""], timestepper) + u_adaptive = solve_waveprop(mesh, coeff, [""], timestepper_adaptive) + @test u ≈ u_adaptive rtol = 1e-4 mesh = generate_ideal_lv_mesh(4,1,1) coeff = ConstantCoefficient(SymmetricTensor{2,3,Float64}((4.5e-5, 0, 0, 2.0e-5, 0, 1.0e-5))) - u = solve_waveprop(mesh, coeff, [""]) - u_adaptive = solve_waveprop(mesh, coeff, [""], true) - @test u ≈ u_adaptive + u = solve_waveprop(mesh, coeff, [""], timestepper) + u_adaptive = solve_waveprop(mesh, coeff, [""], timestepper_adaptive) + @test u ≈ u_adaptive rtol = 1e-4 mesh = to_mesh(generate_mixed_grid_2D()) coeff = ConstantCoefficient(SymmetricTensor{2,2,Float64}((4.5e-5, 0, 2.0e-5))) - u = solve_waveprop(mesh, coeff, ["Pacemaker", "Myocardium"]) - u_adaptive = solve_waveprop(mesh, coeff, ["Pacemaker", "Myocardium"], true) - @test u ≈ u_adaptive - u = solve_waveprop(mesh, coeff, ["Pacemaker"]) - u_adaptive = solve_waveprop(mesh, coeff, ["Pacemaker"], true) - @test u ≈ u_adaptive - u = solve_waveprop(mesh, coeff, ["Myocardium"]) - u_adaptive = solve_waveprop(mesh, coeff, ["Myocardium"], true) - @test u ≈ u_adaptive + u = solve_waveprop(mesh, coeff, ["Pacemaker", "Myocardium"], timestepper) + u_adaptive = solve_waveprop(mesh, coeff, ["Pacemaker", "Myocardium"], timestepper_adaptive) + @test u ≈ u_adaptive rtol = 1e-4 + u = solve_waveprop(mesh, coeff, ["Pacemaker"], timestepper) + u_adaptive = solve_waveprop(mesh, coeff, ["Pacemaker"], timestepper_adaptive) + @test u ≈ u_adaptive rtol = 1e-4 + u = solve_waveprop(mesh, coeff, ["Myocardium"], timestepper) + u_adaptive = solve_waveprop(mesh, coeff, ["Myocardium"], timestepper_adaptive) + @test u ≈ u_adaptive rtol = 1e-4 mesh = to_mesh(generate_mixed_dimensional_grid_3D()) coeff = ConstantCoefficient(SymmetricTensor{2,3,Float64}((4.5e-5, 0, 0, 2.0e-5, 0, 1.0e-5))) - u = solve_waveprop(mesh, coeff, ["Ventricle"]) - u_adaptive = solve_waveprop(mesh, coeff, ["Ventricle"]) - @test u ≈ u_adaptive + u = solve_waveprop(mesh, coeff, ["Ventricle"], timestepper) + u_adaptive = solve_waveprop(mesh, coeff, ["Ventricle"], timestepper_adaptive) + @test u ≈ u_adaptive rtol = 1e-4 coeff = ConstantCoefficient(SymmetricTensor{2,3,Float64}((5e-5, 0, 0, 5e-5, 0, 5e-5))) - u = solve_waveprop(mesh, coeff, ["Purkinje"]) - u_adaptive = solve_waveprop(mesh, coeff, ["Purkinje"]) - @test u ≈ u_adaptive - u = solve_waveprop(mesh, coeff, ["Ventricle", "Purkinje"]) - u_adaptive = solve_waveprop(mesh, coeff, ["Ventricle", "Purkinje"]) - @test u ≈ u_adaptive + u = solve_waveprop(mesh, coeff, ["Purkinje"], timestepper) + u_adaptive = solve_waveprop(mesh, coeff, ["Purkinje"], timestepper_adaptive) + @test u ≈ u_adaptive rtol = 1e-4 + u = solve_waveprop(mesh, coeff, ["Ventricle", "Purkinje"], timestepper) + u_adaptive = solve_waveprop(mesh, coeff, ["Ventricle", "Purkinje"], timestepper_adaptive) + @test u ≈ u_adaptive rtol = 1e-4 end diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 56c1265f..58ff2209 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -94,19 +94,18 @@ fsplit_NaN = GenericSplitFunction((f1,f_NaN), (f1dofs, f_NaN_dofs)) @testset "OperatorSplitting" begin for TimeStepperType in (LieTrotterGodunov,) - for controller in (Thunderbolt.ReactionTangentController(0.5, 1.0, (0.01, 0.3)),) timestepper = TimeStepperType( (DummyForwardEuler(), DummyForwardEuler()) ) - timestepper_adaptive = Thunderbolt.AdaptiveOperatorSplittingAlgorithm(timestepper, controller) + timestepper_adaptive = Thunderbolt.ReactionTangentController(timestepper, 0.5, 1.0, (0.01, 0.3)) timestepper_inner = TimeStepperType( (DummyForwardEuler(), DummyForwardEuler()) ) - timestepper_inner_adaptive = Thunderbolt.AdaptiveOperatorSplittingAlgorithm(timestepper_inner, controller) #TODO: Copy the controller instead + timestepper_inner_adaptive = Thunderbolt.ReactionTangentController(timestepper_inner, 0.5, 1.0, (0.01, 0.3)) #TODO: Copy the controller instead timestepper2 = TimeStepperType( (DummyForwardEuler(), timestepper_inner) ) - timestepper2_adaptive = Thunderbolt.AdaptiveOperatorSplittingAlgorithm(timestepper2, controller) + timestepper2_adaptive = Thunderbolt.ReactionTangentController(timestepper2, 0.5, 1.0, (0.01, 0.3)) for (tstepper1, tstepper_inner, tstepper2) in ( (timestepper, timestepper_inner, timestepper2), @@ -162,7 +161,6 @@ fsplit_NaN = GenericSplitFunction((f1,f_NaN), (f1dofs, f_NaN_dofs)) # DiffEqBase.solve!(integrator_NaN) # @test integrator_NaN.sol.retcode == DiffEqBase.ReturnCode.Failure end - end # integrator = DiffEqBase.init(prob, timestepper, dt=0.01, verbose=true) # for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) end # integrator_adaptive = DiffEqBase.init(prob, timestepper_adaptive, dt=0.01, verbose=true) From 0ccbd0252570c2bd651ba52670c2120b0df05da3 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 19 Aug 2024 15:50:23 +0200 Subject: [PATCH 21/37] Docs --- docs/src/api-reference/solver.md | 2 -- docs/src/assets/references.bib | 6 +++--- src/solver/adaptivity.jl | 4 ++++ 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/docs/src/api-reference/solver.md b/docs/src/api-reference/solver.md index a586bd80..0e4b2c5b 100644 --- a/docs/src/api-reference/solver.md +++ b/docs/src/api-reference/solver.md @@ -38,7 +38,5 @@ Thunderbolt.OS.OperatorSplittingIntegrator ## Operator Splitting Adaptivity ```@docs -Thunderbolt.AdaptiveOperatorSplittingAlgorithm Thunderbolt.ReactionTangentController -Thunderbolt.get_reaction_tangent ``` diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index e905f0ad..fe429478 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -280,11 +280,11 @@ @article{PotDubRicVinGul:2006:cmb } @article{OgiBalPer:2023:seats, author = {Ogiermann, Dennis and Perotti, Luigi E. and Balzani, Daniel}, -title = {A simple and efficient adaptive time stepping technique for low-order operator splitting schemes applied to cardiac electrophysiology}, journal = {International Journal for Numerical Methods in Biomedical Engineering}, +title = {A simple and efficient adaptive time stepping technique for low-order operator splitting schemes applied to cardiac electrophysiology}, +year = {2023} volume = {39}, number = {2}, pages = {e3670}, -url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/cnm.3670}, -year = {2023} +doi = {https://doi.org/10.1002/cnm.3670}, } diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index f0469c58..62ce3fa2 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -4,6 +4,10 @@ abstract type AbstractTimeAdaptionAlgorithm end ReactionTangentController{T <: Real} <: OS.AbstractOperatorSplittingAlgorithm A timestep length controller for [`LieTrotterGodunov`](@ref) [Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite) operator splitting using the reaction tangent as proposed in [OgiBalPer:2023:seats](@cite) +The next timestep length is calculated as +```math +\\sigma\\left(R_{\\max }\\right):=\\left(1.0-\\frac{1}{1+\\exp \\left(\\left(\\sigma_{\\mathrm{c}}-R_{\\max }\\right) \\cdot \\sigma_{\\mathrm{s}}\\right)}\\right) \\cdot\\left(\\Delta t_{\\max }-\\Delta t_{\\min }\\right)+\\Delta t_{\\min } +``` # Fields - `σ_s::T`: steepness - `σ_c::T`: offset in R axis From b505899c3865a78bf9cd356ad5f3c728bbc57f9a Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 19 Aug 2024 15:56:07 +0200 Subject: [PATCH 22/37] oopsie --- src/solver/adaptivity.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index 62ce3fa2..7d2b704e 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -1,7 +1,7 @@ abstract type AbstractTimeAdaptionAlgorithm end """ - ReactionTangentController{T <: Real} <: OS.AbstractOperatorSplittingAlgorithm + ReactionTangentController{LTG <: OS.LieTrotterGodunov, T <: Real} <: OS.AbstractOperatorSplittingAlgorithm A timestep length controller for [`LieTrotterGodunov`](@ref) [Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite) operator splitting using the reaction tangent as proposed in [OgiBalPer:2023:seats](@cite) The next timestep length is calculated as @@ -9,6 +9,7 @@ The next timestep length is calculated as \\sigma\\left(R_{\\max }\\right):=\\left(1.0-\\frac{1}{1+\\exp \\left(\\left(\\sigma_{\\mathrm{c}}-R_{\\max }\\right) \\cdot \\sigma_{\\mathrm{s}}\\right)}\\right) \\cdot\\left(\\Delta t_{\\max }-\\Delta t_{\\min }\\right)+\\Delta t_{\\min } ``` # Fields +- `ltg`::`LTG`: `LieTrotterGodunov` algorithm - `σ_s::T`: steepness - `σ_c::T`: offset in R axis - `Δt_bounds::NTuple{2,T}`: lower and upper timestep length bounds From f3ad8d3b89fb5f57fde74c23136be1b15e263b49 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 20 Aug 2024 13:16:01 +0200 Subject: [PATCH 23/37] Next: Error checking for NaNs --- src/solver/adaptivity.jl | 22 ++++--- src/solver/operator_splitting/integrator.jl | 1 + test/test_integrators.jl | 72 ++++++++++++++------- 3 files changed, 63 insertions(+), 32 deletions(-) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index 7d2b704e..7cf1aa8f 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -43,18 +43,20 @@ Returns the maximal reaction magnitude using the [`PointwiseODEFunction`](@ref) It is assumed that the problem containing the reaction tangent is a [`PointwiseODEFunction`](@ref). """ @inline function get_reaction_tangent(integrator::OS.OperatorSplittingIntegrator) - reaction_subintegrators = unrolled_filter(subintegrator -> subintegrator.f isa PointwiseODEFunction, integrator.subintegrators) - isempty(reaction_subintegrators) && @error "PointwiseODEFunction not found" - _get_reaction_tangent(reaction_subintegrators) + _get_reaction_tangent(integrator.subintegrators) end -@inline #=@unroll=# function _get_reaction_tangent(subintegrators) - subintegrator = subintegrators[1] - # @unroll for subintegrator in subintegrators - #TODO: It should be either all the the same type or just one subintegrator after filtering, don't unroll? - φₘidx = transmembranepotential_index(subintegrator.f.ode) - return maximum(@view subintegrator.cache.dumat[:, φₘidx]) - # end +@inline @unroll function _get_reaction_tangent(subintegrators) + @unroll for subintegrator in subintegrators + if subintegrator isa Tuple + temp = _get_reaction_tangent(subintegrator) + isnan(temp) || return temp + elseif subintegrator.f isa PointwiseODEFunction + φₘidx = transmembranepotential_index(subintegrator.f.ode) + return maximum(@view subintegrator.cache.dumat[:, φₘidx]) + end + end + return NaN end @inline function OS.stepsize_controller!(integrator::OS.OperatorSplittingIntegrator, alg::ReactionTangentController) diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index 8b5cfc70..3ed90b00 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -248,6 +248,7 @@ function DiffEqBase.SciMLBase.done(integrator::OperatorSplittingIntegrator) end function DiffEqBase.postamble!(integrator::OperatorSplittingIntegrator) + # DiffEqBase.check_error!(integrator) != ReturnCode.Success DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator) end diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 58ff2209..a3f6adbc 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -7,18 +7,33 @@ using UnPack ODEFunction = Thunderbolt.DiffEqBase.ODEFunction # For testing purposes -struct DummyForwardEuler +struct DummyForwardEuler <: Thunderbolt.AbstractSolver end -mutable struct DummyForwardEulerCache{duType, uType} <: Thunderbolt.AbstractTimeSolverCache +mutable struct DummyForwardEulerCache{duType, uType, duMatType} <: Thunderbolt.AbstractTimeSolverCache du::duType + dumat::duMatType uₙ::uType uₙ₋₁::uType end # Dispatch for leaf construction function OS.construct_inner_cache(f::ODEFunction, alg::DummyForwardEuler, u::AbstractArray, uprev::AbstractArray) - DummyForwardEulerCache(copy(uprev), copy(uprev), copy(uprev)) + dumat = reshape(uprev, (:,1)) + DummyForwardEulerCache(copy(uprev), dumat, copy(uprev), copy(uprev)) +end +Thunderbolt.num_states(::Any) = 2 +Thunderbolt.transmembranepotential_index(::Any) = 1 + +function Thunderbolt.setup_solver_cache(f::PointwiseODEFunction, solver::DummyForwardEuler, t₀) + @unpack npoints, ode = f + + du = Thunderbolt.create_system_vector(Vector{Float64}, f) + uₙ = Thunderbolt.create_system_vector(Vector{Float64}, f) + uₙ₋₁ = Thunderbolt.create_system_vector(Vector{Float64}, f) + dumat = reshape(du, (:,1)) + + return DummyForwardEulerCache(du, dumat, uₙ, uₙ₋₁) end # Dispatch innermost solve @@ -26,7 +41,7 @@ function OS.advance_solution_to!(integ::ThunderboltTimeIntegrator, cache::DummyF @unpack f, dt, u, p, t = integ @unpack du = cache - f(du, u, p, t) + f isa Thunderbolt.PointwiseODEFunction ? f.ode(du, u, p, t) : f(du, u, p, t) @. u += dt * du end @@ -58,7 +73,13 @@ f2 = ODEFunction(ode2) # ode_true and ode1/ode2 side by side to see how they connect. f1dofs = [1,2,3] f2dofs = [1,3] -fsplit1 = GenericSplitFunction((f1,f2), (f1dofs, f2dofs)) +fpw = PointwiseODEFunction( + 1, + f2, + [0.0] +) + +fsplit1 = GenericSplitFunction((f1,fpw), (f1dofs, f2dofs)) # Now the usual setup just with our new problem type. # u0 = rand(3) @@ -81,38 +102,45 @@ f3 = ODEFunction(ode3) f1dofs = [1,2,3] f2dofs = [1,3] f3dofs = [1,3] -fsplit2_inner = GenericSplitFunction((f3,f3), (f3dofs, f3dofs)) +fsplit2_inner = GenericSplitFunction((fpw,f3), (f3dofs, f3dofs)) fsplit2_outer = GenericSplitFunction((f1,fsplit2_inner), (f1dofs, f2dofs)) function ode_NaN(du, u, p, t) du[1] = NaN du[2] = 0.01u[1] end + f_NaN = ODEFunction(ode_NaN) +fpw_NaN = PointwiseODEFunction( + 1, + f_NaN, + [0.0] +) f_NaN_dofs = f3dofs -fsplit_NaN = GenericSplitFunction((f1,f_NaN), (f1dofs, f_NaN_dofs)) - +fsplit_NaN = GenericSplitFunction((f1,fpw_NaN), (f1dofs, f_NaN_dofs)) +dt = 0.01π +adaptive_tstep_range = (dt * 1, dt * 5) @testset "OperatorSplitting" begin for TimeStepperType in (LieTrotterGodunov,) timestepper = TimeStepperType( (DummyForwardEuler(), DummyForwardEuler()) ) - timestepper_adaptive = Thunderbolt.ReactionTangentController(timestepper, 0.5, 1.0, (0.01, 0.3)) + timestepper_adaptive = Thunderbolt.ReactionTangentController(timestepper, 0.5, 1.0, adaptive_tstep_range) timestepper_inner = TimeStepperType( (DummyForwardEuler(), DummyForwardEuler()) ) - timestepper_inner_adaptive = Thunderbolt.ReactionTangentController(timestepper_inner, 0.5, 1.0, (0.01, 0.3)) #TODO: Copy the controller instead + timestepper_inner_adaptive = Thunderbolt.ReactionTangentController(timestepper_inner, 0.5, 1.0, adaptive_tstep_range) #TODO: Copy the controller instead timestepper2 = TimeStepperType( (DummyForwardEuler(), timestepper_inner) ) - timestepper2_adaptive = Thunderbolt.ReactionTangentController(timestepper2, 0.5, 1.0, (0.01, 0.3)) + timestepper2_adaptive = Thunderbolt.ReactionTangentController(timestepper2, 0.5, 1.0, adaptive_tstep_range) for (tstepper1, tstepper_inner, tstepper2) in ( (timestepper, timestepper_inner, timestepper2), - # (timestepper_adaptive, timestepper_inner_adaptive, timestepper2_adaptive) + (timestepper_adaptive, timestepper_inner_adaptive, timestepper2_adaptive) ) # The remaining code works as usual. - integrator = DiffEqBase.init(prob, tstepper1, dt=0.01π, verbose=true) + integrator = DiffEqBase.init(prob, tstepper1, dt=dt, verbose=true) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default DiffEqBase.solve!(integrator) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success @@ -137,7 +165,7 @@ fsplit_NaN = GenericSplitFunction((f1,f_NaN), (f1dofs, f_NaN_dofs)) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success prob2 = OperatorSplittingProblem(fsplit2_outer, u0, tspan) - integrator2 = DiffEqBase.init(prob2, tstepper2, dt=0.02π, verbose=true) + integrator2 = DiffEqBase.init(prob2, tstepper2, dt=dt, verbose=true) @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Default DiffEqBase.solve!(integrator2) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success @@ -155,17 +183,17 @@ fsplit_NaN = GenericSplitFunction((f1,f_NaN), (f1dofs, f_NaN_dofs)) DiffEqBase.solve!(integrator2) @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Success @testset "NaNs" begin - # prob_NaN = OperatorSplittingProblem(fsplit_NaN, u0, tspan) - # integrator_NaN = DiffEqBase.init(prob, tstepper1, dt=0.01π, verbose=true) - # @test integrator_NaN.sol.retcode == DiffEqBase.ReturnCode.Default + prob_NaN = OperatorSplittingProblem(fsplit_NaN, u0, tspan) + integrator_NaN = DiffEqBase.init(prob, tstepper1, dt=dt, verbose=true) + @test integrator_NaN.sol.retcode == DiffEqBase.ReturnCode.Default # DiffEqBase.solve!(integrator_NaN) # @test integrator_NaN.sol.retcode == DiffEqBase.ReturnCode.Failure end - # integrator = DiffEqBase.init(prob, timestepper, dt=0.01, verbose=true) - # for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) end - # integrator_adaptive = DiffEqBase.init(prob, timestepper_adaptive, dt=0.01, verbose=true) - # for (u, t) in DiffEqBase.TimeChoiceIterator(integrator_adaptive, 0.0:5.0:100.0) end - # @test isapprox(integrator_adaptive.u, integrator.u, atol=1e-5) + integrator = DiffEqBase.init(prob, timestepper, dt=dt, verbose=true) + for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) end + integrator_adaptive = DiffEqBase.init(prob, timestepper_adaptive, dt=dt, verbose=true) + for (u, t) in DiffEqBase.TimeChoiceIterator(integrator_adaptive, 0.0:5.0:100.0) end + @test isapprox(integrator_adaptive.u, integrator.u, atol=1e-5) end end end From c685bf36968091b854df68a9c8801647f89ce0db Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Tue, 20 Aug 2024 15:51:10 +0200 Subject: [PATCH 24/37] Try mabe CI works now? --- src/solver/adaptivity.jl | 6 +- src/solver/operator_splitting/integrator.jl | 7 +- test/test_integrators.jl | 221 ++++++++++---------- 3 files changed, 121 insertions(+), 113 deletions(-) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index 7cf1aa8f..39015137 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -69,7 +69,11 @@ end @unpack Rₙ₊₁, Rₙ = integrator.cache @unpack σ_s, σ_c, Δt_bounds = alg R = max(Rₙ, Rₙ₊₁) - integrator._dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] + if DiffEqBase.NAN_CHECK(R) + # TODO: Maybe throw a warning? + else + integrator._dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] + end return nothing end diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index 3ed90b00..440ea670 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -147,7 +147,11 @@ function DiffEqBase.solve!(integrator::OperatorSplittingIntegrator) __step!(integrator) end DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator) - integrator.sol = DiffEqBase.solution_new_retcode(integrator.sol, DiffEqBase.ReturnCode.Success) + if DiffEqBase.NAN_CHECK(integrator.u) + integrator.sol = DiffEqBase.solution_new_retcode(integrator.sol, DiffEqBase.ReturnCode.Failure) + else + integrator.sol = DiffEqBase.solution_new_retcode(integrator.sol, DiffEqBase.ReturnCode.Success) + end return integrator.sol end @@ -248,7 +252,6 @@ function DiffEqBase.SciMLBase.done(integrator::OperatorSplittingIntegrator) end function DiffEqBase.postamble!(integrator::OperatorSplittingIntegrator) - # DiffEqBase.check_error!(integrator) != ReturnCode.Success DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator) end diff --git a/test/test_integrators.jl b/test/test_integrators.jl index a3f6adbc..0272f5c2 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -4,124 +4,126 @@ using UnPack @testset "Operator Splitting API" begin -ODEFunction = Thunderbolt.DiffEqBase.ODEFunction + ODEFunction = Thunderbolt.DiffEqBase.ODEFunction -# For testing purposes -struct DummyForwardEuler <: Thunderbolt.AbstractSolver -end + # For testing purposes + struct DummyForwardEuler <: Thunderbolt.AbstractSolver + end -mutable struct DummyForwardEulerCache{duType, uType, duMatType} <: Thunderbolt.AbstractTimeSolverCache - du::duType - dumat::duMatType - uₙ::uType - uₙ₋₁::uType -end + mutable struct DummyForwardEulerCache{duType, uType, duMatType} <: Thunderbolt.AbstractTimeSolverCache + du::duType + dumat::duMatType + uₙ::uType + uₙ₋₁::uType + end -# Dispatch for leaf construction -function OS.construct_inner_cache(f::ODEFunction, alg::DummyForwardEuler, u::AbstractArray, uprev::AbstractArray) - dumat = reshape(uprev, (:,1)) - DummyForwardEulerCache(copy(uprev), dumat, copy(uprev), copy(uprev)) -end -Thunderbolt.num_states(::Any) = 2 -Thunderbolt.transmembranepotential_index(::Any) = 1 + # Dispatch for leaf construction + function OS.construct_inner_cache(f::ODEFunction, alg::DummyForwardEuler, u::AbstractArray, uprev::AbstractArray) + dumat = reshape(uprev, (:,1)) + DummyForwardEulerCache(copy(uprev), dumat, copy(uprev), copy(uprev)) + end + Thunderbolt.num_states(::Any) = 2 + Thunderbolt.transmembranepotential_index(::Any) = 1 -function Thunderbolt.setup_solver_cache(f::PointwiseODEFunction, solver::DummyForwardEuler, t₀) - @unpack npoints, ode = f + function Thunderbolt.setup_solver_cache(f::PointwiseODEFunction, solver::DummyForwardEuler, t₀) + @unpack npoints, ode = f - du = Thunderbolt.create_system_vector(Vector{Float64}, f) - uₙ = Thunderbolt.create_system_vector(Vector{Float64}, f) - uₙ₋₁ = Thunderbolt.create_system_vector(Vector{Float64}, f) - dumat = reshape(du, (:,1)) + du = Thunderbolt.create_system_vector(Vector{Float64}, f) + uₙ = Thunderbolt.create_system_vector(Vector{Float64}, f) + uₙ₋₁ = Thunderbolt.create_system_vector(Vector{Float64}, f) + dumat = reshape(du, (:,1)) - return DummyForwardEulerCache(du, dumat, uₙ, uₙ₋₁) -end + return DummyForwardEulerCache(du, dumat, uₙ, uₙ₋₁) + end -# Dispatch innermost solve -function OS.advance_solution_to!(integ::ThunderboltTimeIntegrator, cache::DummyForwardEulerCache, tend) - @unpack f, dt, u, p, t = integ - @unpack du = cache + # Dispatch innermost solve + function OS.advance_solution_to!(integ::ThunderboltTimeIntegrator, cache::DummyForwardEulerCache, tend) + @unpack f, dt, u, p, t = integ + @unpack du = cache - f isa Thunderbolt.PointwiseODEFunction ? f.ode(du, u, p, t) : f(du, u, p, t) - @. u += dt * du -end + f isa Thunderbolt.PointwiseODEFunction ? f.ode(du, u, p, t) : f(du, u, p, t) + @. u += dt * du + cache.dumat[:,1] .= du + end -# Operator splitting + # Operator splitting -# Reference -function ode_true(du, u, p, t) - du .= -0.1u - du[1] += 0.01u[3] - du[3] += 0.01u[1] -end + # Reference + function ode_true(du, u, p, t) + du .= -0.1u + du[1] += 0.01u[3] + du[3] += 0.01u[1] + end -# Setup individual functions -# Diagonal components -function ode1(du, u, p, t) - @. du = -0.1u -end -# Offdiagonal components -function ode2(du, u, p, t) - du[1] = 0.01u[2] - du[2] = 0.01u[1] -end + # Setup individual functions + # Diagonal components + function ode1(du, u, p, t) + @. du = -0.1u + end + # Offdiagonal components + function ode2(du, u, p, t) + du[1] = 0.01u[2] + du[2] = 0.01u[1] + end -f1 = ODEFunction(ode1) -f2 = ODEFunction(ode2) - -# Here we describe index sets f1dofs and f2dofs that map the -# local indices in f1 and f2 into the global problem. Just put -# ode_true and ode1/ode2 side by side to see how they connect. -f1dofs = [1,2,3] -f2dofs = [1,3] -fpw = PointwiseODEFunction( - 1, - f2, - [0.0] -) - -fsplit1 = GenericSplitFunction((f1,fpw), (f1dofs, f2dofs)) - -# Now the usual setup just with our new problem type. -# u0 = rand(3) -u0 = [0.7611944793397108 - 0.9059606424982555 - 0.5755174199139956] -tspan = (0.0,100.0) -prob = OperatorSplittingProblem(fsplit1, u0, tspan) - -# Now some recursive splitting -function ode3(du, u, p, t) - du[1] = 0.005u[2] - du[2] = 0.005u[1] -end -f3 = ODEFunction(ode3) -# The time stepper carries the individual solver information. - -# Note that we define the dof indices w.r.t the parent function. -# Hence the indices for `fsplit2_inner` are. -f1dofs = [1,2,3] -f2dofs = [1,3] -f3dofs = [1,3] -fsplit2_inner = GenericSplitFunction((fpw,f3), (f3dofs, f3dofs)) -fsplit2_outer = GenericSplitFunction((f1,fsplit2_inner), (f1dofs, f2dofs)) - -function ode_NaN(du, u, p, t) - du[1] = NaN - du[2] = 0.01u[1] -end + f1 = ODEFunction(ode1) + f2 = ODEFunction(ode2) + + # Here we describe index sets f1dofs and f2dofs that map the + # local indices in f1 and f2 into the global problem. Just put + # ode_true and ode1/ode2 side by side to see how they connect. + f1dofs = [1,2,3] + f2dofs = [1,3] + fpw = PointwiseODEFunction( + 1, + f2, + [0.0] + ) + + fsplit1 = GenericSplitFunction((f1,fpw), (f1dofs, f2dofs)) + + # Now the usual setup just with our new problem type. + # u0 = rand(3) + u0 = [0.7611944793397108 + 0.9059606424982555 + 0.5755174199139956] + tspan = (0.0,100.0) + prob = OperatorSplittingProblem(fsplit1, u0, tspan) + + # Now some recursive splitting + function ode3(du, u, p, t) + du[1] = 0.005u[2] + du[2] = 0.005u[1] + end + f3 = ODEFunction(ode3) + # The time stepper carries the individual solver information. + + # Note that we define the dof indices w.r.t the parent function. + # Hence the indices for `fsplit2_inner` are. + f1dofs = [1,2,3] + f2dofs = [1,3] + f3dofs = [1,3] + fsplit2_inner = GenericSplitFunction((fpw,f3), (f3dofs, f3dofs)) + fsplit2_outer = GenericSplitFunction((f1,fsplit2_inner), (f1dofs, f2dofs)) + + function ode_NaN(du, u, p, t) + du[1] = NaN + du[2] = 0.01u[1] + end -f_NaN = ODEFunction(ode_NaN) -fpw_NaN = PointwiseODEFunction( - 1, - f_NaN, - [0.0] -) -f_NaN_dofs = f3dofs -fsplit_NaN = GenericSplitFunction((f1,fpw_NaN), (f1dofs, f_NaN_dofs)) -dt = 0.01π -adaptive_tstep_range = (dt * 1, dt * 5) -@testset "OperatorSplitting" begin - for TimeStepperType in (LieTrotterGodunov,) + f_NaN = ODEFunction(ode_NaN) + fpw_NaN = PointwiseODEFunction( + 1, + f_NaN, + [0.0] + ) + f_NaN_dofs = f3dofs + fsplit_NaN = GenericSplitFunction((f1,fpw_NaN), (f1dofs, f_NaN_dofs)) + prob_NaN = OperatorSplittingProblem(fsplit_NaN, u0, tspan) + dt = 0.01π + adaptive_tstep_range = (dt * 1, dt * 5) + @testset "OperatorSplitting" begin + for TimeStepperType in (LieTrotterGodunov,) timestepper = TimeStepperType( (DummyForwardEuler(), DummyForwardEuler()) ) @@ -183,12 +185,12 @@ adaptive_tstep_range = (dt * 1, dt * 5) DiffEqBase.solve!(integrator2) @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Success @testset "NaNs" begin - prob_NaN = OperatorSplittingProblem(fsplit_NaN, u0, tspan) - integrator_NaN = DiffEqBase.init(prob, tstepper1, dt=dt, verbose=true) + integrator_NaN = DiffEqBase.init(prob_NaN, tstepper1, dt=dt, verbose=true) @test integrator_NaN.sol.retcode == DiffEqBase.ReturnCode.Default - # DiffEqBase.solve!(integrator_NaN) - # @test integrator_NaN.sol.retcode == DiffEqBase.ReturnCode.Failure + DiffEqBase.solve!(integrator_NaN) + @test integrator_NaN.sol.retcode == DiffEqBase.ReturnCode.Failure end + end integrator = DiffEqBase.init(prob, timestepper, dt=dt, verbose=true) for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, 0.0:5.0:100.0) end integrator_adaptive = DiffEqBase.init(prob, timestepper_adaptive, dt=dt, verbose=true) @@ -196,7 +198,6 @@ adaptive_tstep_range = (dt * 1, dt * 5) @test isapprox(integrator_adaptive.u, integrator.u, atol=1e-5) end end -end # tnext = tspan[1]+0.01 # @btime OS.advance_solution_to!($integrator, $tnext) setup=(DiffEqBase.reinit!(integrator, u0; tspan)) From 06a541b017bba89bfec3d5bb879b6fe0ba6f43b0 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Wed, 21 Aug 2024 16:16:45 +0200 Subject: [PATCH 25/37] Don't merge, homeoffice push --- src/solver/adaptivity.jl | 7 +++++-- src/solver/operator_splitting/solver.jl | 4 ++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index 39015137..2287a2f8 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -47,13 +47,16 @@ It is assumed that the problem containing the reaction tangent is a [`PointwiseO end @inline @unroll function _get_reaction_tangent(subintegrators) + has_more_than_one_reaction_tangent = false + R = NaN @unroll for subintegrator in subintegrators if subintegrator isa Tuple temp = _get_reaction_tangent(subintegrator) isnan(temp) || return temp elseif subintegrator.f isa PointwiseODEFunction + isnan(R) || @error "More than one reaction tangent were found" φₘidx = transmembranepotential_index(subintegrator.f.ode) - return maximum(@view subintegrator.cache.dumat[:, φₘidx]) + R = maximum(@view subintegrator.cache.dumat[:, φₘidx]) end end return NaN @@ -70,7 +73,7 @@ end @unpack σ_s, σ_c, Δt_bounds = alg R = max(Rₙ, Rₙ₊₁) if DiffEqBase.NAN_CHECK(R) - # TODO: Maybe throw a warning? + @error "reaction tangent cannot be NaN" else integrator._dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] end diff --git a/src/solver/operator_splitting/solver.jl b/src/solver/operator_splitting/solver.jl index 529fcaf3..e824a2f8 100644 --- a/src/solver/operator_splitting/solver.jl +++ b/src/solver/operator_splitting/solver.jl @@ -18,8 +18,8 @@ struct LieTrotterGodunovCache{uType, tmpType, iiType} <: AbstractOperatorSplitti inner_caches::iiType end -get_u(cache#=::LieTrotterGodunovCache=#) = cache.u -get_uprev(cache#=::LieTrotterGodunovCache=#) = cache.uprev +get_u(cache::LieTrotterGodunovCache) = cache.u +get_uprev(cache::LieTrotterGodunovCache) = cache.uprev # Dispatch for outer construction function init_cache(prob::OperatorSplittingProblem, alg::LieTrotterGodunov; dt, kwargs...) # TODO From ffe212845cceb7ae3c69d849f5d1f84bf6d3de5e Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 26 Aug 2024 15:18:03 +0200 Subject: [PATCH 26/37] nans test --- src/solver/adaptivity.jl | 26 +++++++++------------ src/solver/operator_splitting/integrator.jl | 11 +++++++++ 2 files changed, 22 insertions(+), 15 deletions(-) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index 2287a2f8..6f4dfdf9 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -43,23 +43,23 @@ Returns the maximal reaction magnitude using the [`PointwiseODEFunction`](@ref) It is assumed that the problem containing the reaction tangent is a [`PointwiseODEFunction`](@ref). """ @inline function get_reaction_tangent(integrator::OS.OperatorSplittingIntegrator) - _get_reaction_tangent(integrator.subintegrators) + R, _ = _get_reaction_tangent(integrator.subintegrators) + return R end -@inline @unroll function _get_reaction_tangent(subintegrators) - has_more_than_one_reaction_tangent = false - R = NaN +@inline @unroll function _get_reaction_tangent(subintegrators, n_reaction_tangents::Int = 0) + R = 0.0 @unroll for subintegrator in subintegrators if subintegrator isa Tuple - temp = _get_reaction_tangent(subintegrator) - isnan(temp) || return temp + R, n_reaction_tangents = _get_reaction_tangent(subintegrator, n_reaction_tangents) elseif subintegrator.f isa PointwiseODEFunction - isnan(R) || @error "More than one reaction tangent were found" + n_reaction_tangents += 1 φₘidx = transmembranepotential_index(subintegrator.f.ode) - R = maximum(@view subintegrator.cache.dumat[:, φₘidx]) + R = max(R, maximum(@view subintegrator.cache.dumat[:, φₘidx])) end end - return NaN + @assert n_reaction_tangents == 1 "No or multiple integrators using PointwiseODEFunction found" + return (R, n_reaction_tangents) end @inline function OS.stepsize_controller!(integrator::OS.OperatorSplittingIntegrator, alg::ReactionTangentController) @@ -72,11 +72,7 @@ end @unpack Rₙ₊₁, Rₙ = integrator.cache @unpack σ_s, σ_c, Δt_bounds = alg R = max(Rₙ, Rₙ₊₁) - if DiffEqBase.NAN_CHECK(R) - @error "reaction tangent cannot be NaN" - else - integrator._dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] - end + integrator._dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] return nothing end @@ -104,4 +100,4 @@ end function OS.build_subintegrators_recursive(f::GenericSplitFunction, synchronizers::Tuple, p::Tuple, cache::ReactionTangentControllerCache, u::AbstractArray, uprev::AbstractArray, t, dt, dof_range, uparent, tstops, _tstops, saveat, _saveat) OS.build_subintegrators_recursive(f, synchronizers, p, cache.ltg_cache, u, uprev, t, dt, dof_range, uparent, tstops, _tstops, saveat, _saveat) -end \ No newline at end of file +end diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index 440ea670..4fb13690 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -144,6 +144,7 @@ end # either called directly (after init), or by DiffEqBase.solve (via __solve) function DiffEqBase.solve!(integrator::OperatorSplittingIntegrator) while !isempty(integrator.tstops) + DiffEqBase.check_error!(integrator) ∉ (DiffEqBase.ReturnCode.Success, DiffEqBase.ReturnCode.Default) && return __step!(integrator) end DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator) @@ -159,13 +160,22 @@ function DiffEqBase.step!(integrator::OperatorSplittingIntegrator) if integrator.advance_to_tstop tstop = first(integrator.tstops) while !reached_tstop(integrator, tstop) + DiffEqBase.check_error!(integrator) ∉ (DiffEqBase.ReturnCode.Success, DiffEqBase.ReturnCode.Default) && return __step!(integrator) end else + DiffEqBase.check_error!(integrator) ∉ (DiffEqBase.ReturnCode.Success, DiffEqBase.ReturnCode.Default) && return __step!(integrator) end end +function DiffEqBase.check_error!(integrator::OperatorSplittingIntegrator) + if DiffEqBase.NAN_CHECK(integrator._dt) # replace with https://github.com/SciML/OrdinaryDiffEq.jl/blob/373a8eec8024ef1acc6c5f0c87f479aa0cf128c3/lib/OrdinaryDiffEqCore/src/iterator_interface.jl#L5-L6 after moving to sciml integrators + integrator.sol = DiffEqBase.solution_new_retcode(integrator.sol, DiffEqBase.ReturnCode.Failure) + end + return integrator.sol.retcode +end + function DiffEqBase.step!(integrator::OperatorSplittingIntegrator, dt, stop_at_tdt = false) # OridinaryDiffEq lets dt be negative if tdir is -1, but that's inconsistent dt <= zero(dt) && error("dt must be positive") @@ -173,6 +183,7 @@ function DiffEqBase.step!(integrator::OperatorSplittingIntegrator, dt, stop_at_t tnext = integrator.t + tdir(integrator) * dt stop_at_tdt && DiffEqBase.add_tstop!(integrator, tnext) while !reached_tstop(integrator, tnext, stop_at_tdt) + DiffEqBase.check_error!(integrator) ∉ (DiffEqBase.ReturnCode.Success, DiffEqBase.ReturnCode.Default) && return __step!(integrator) end end From 61ecbc9c4c6a6ef83afdae7ebc8f0d865fcc7443 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 26 Aug 2024 15:28:01 +0200 Subject: [PATCH 27/37] multiple pwodef --- test/test_integrators.jl | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 0272f5c2..7d2842f1 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -106,6 +106,8 @@ using UnPack fsplit2_inner = GenericSplitFunction((fpw,f3), (f3dofs, f3dofs)) fsplit2_outer = GenericSplitFunction((f1,fsplit2_inner), (f1dofs, f2dofs)) + prob2 = OperatorSplittingProblem(fsplit2_outer, u0, tspan) + function ode_NaN(du, u, p, t) du[1] = NaN du[2] = 0.01u[1] @@ -120,6 +122,11 @@ using UnPack f_NaN_dofs = f3dofs fsplit_NaN = GenericSplitFunction((f1,fpw_NaN), (f1dofs, f_NaN_dofs)) prob_NaN = OperatorSplittingProblem(fsplit_NaN, u0, tspan) + + fsplit_multiple_pwode_outer = GenericSplitFunction((fpw,fsplit2_inner), (f3dofs, f2dofs)) + + prob_multiple_pwode = OperatorSplittingProblem(fsplit_multiple_pwode_outer, u0, tspan) + dt = 0.01π adaptive_tstep_range = (dt * 1, dt * 5) @testset "OperatorSplitting" begin @@ -166,7 +173,6 @@ using UnPack DiffEqBase.solve!(integrator) @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success - prob2 = OperatorSplittingProblem(fsplit2_outer, u0, tspan) integrator2 = DiffEqBase.init(prob2, tstepper2, dt=dt, verbose=true) @test integrator2.sol.retcode == DiffEqBase.ReturnCode.Default DiffEqBase.solve!(integrator2) @@ -196,6 +202,10 @@ using UnPack integrator_adaptive = DiffEqBase.init(prob, timestepper_adaptive, dt=dt, verbose=true) for (u, t) in DiffEqBase.TimeChoiceIterator(integrator_adaptive, 0.0:5.0:100.0) end @test isapprox(integrator_adaptive.u, integrator.u, atol=1e-5) + @testset "Multiple `PointwiseODEFunction`s" begin + integrator_multiple_pwode = DiffEqBase.init(prob_multiple_pwode, timestepper2_adaptive, dt=dt, verbose=true) + @test_throws AssertionError("No or multiple integrators using PointwiseODEFunction found") DiffEqBase.solve!(integrator_multiple_pwode) + end end end From fbeed4092799e1f713d6677fb1119dc414fb0446 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 26 Aug 2024 16:18:28 +0200 Subject: [PATCH 28/37] Inf --- src/solver/adaptivity.jl | 4 ++++ test/test_integrators.jl | 19 +++++++++++++++++++ 2 files changed, 23 insertions(+) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index 6f4dfdf9..c8615564 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -72,6 +72,10 @@ end @unpack Rₙ₊₁, Rₙ = integrator.cache @unpack σ_s, σ_c, Δt_bounds = alg R = max(Rₙ, Rₙ₊₁) + if isinf(σ_s) && R == σ_c + # Handle it? + throw(ErrorException("Δt is undefined for R = σ_c where σ_s = ∞")) + end integrator._dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] return nothing end diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 7d2842f1..595a9fa6 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -127,6 +127,20 @@ using UnPack prob_multiple_pwode = OperatorSplittingProblem(fsplit_multiple_pwode_outer, u0, tspan) + function ode2_force_half(du, u, p, t) + du[1] = 0.5 + du[2] = 0.5 + end + + fpw_force_half = PointwiseODEFunction( + 1, + ODEFunction(ode2_force_half), + [0.0] + ) + + fsplit_force_half = GenericSplitFunction((f1,fpw_force_half), (f1dofs, f2dofs)) + prob_force_half = OperatorSplittingProblem(fsplit_force_half, u0, tspan) + dt = 0.01π adaptive_tstep_range = (dt * 1, dt * 5) @testset "OperatorSplitting" begin @@ -206,6 +220,11 @@ using UnPack integrator_multiple_pwode = DiffEqBase.init(prob_multiple_pwode, timestepper2_adaptive, dt=dt, verbose=true) @test_throws AssertionError("No or multiple integrators using PointwiseODEFunction found") DiffEqBase.solve!(integrator_multiple_pwode) end + @testset "σ_s = Inf, R = σ_c" begin + timestepper_stepfunc_adaptive = Thunderbolt.ReactionTangentController(timestepper, Inf, 0.5, adaptive_tstep_range) + integrator_stepfunc_adaptive = DiffEqBase.init(prob_force_half, timestepper_stepfunc_adaptive, dt=dt, verbose=true) + @test_throws ErrorException("Δt is undefined for R = σ_c where σ_s = ∞") DiffEqBase.solve!(integrator_stepfunc_adaptive) + end end end From e39a77357f0efb2bff15e2b05032a6746c261dd4 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 26 Aug 2024 16:19:19 +0200 Subject: [PATCH 29/37] indentation --- docs/src/assets/references.bib | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib index fe429478..0bf31467 100644 --- a/docs/src/assets/references.bib +++ b/docs/src/assets/references.bib @@ -279,12 +279,12 @@ @article{PotDubRicVinGul:2006:cmb doi={10.1109/TBME.2006.880875} } @article{OgiBalPer:2023:seats, -author = {Ogiermann, Dennis and Perotti, Luigi E. and Balzani, Daniel}, -journal = {International Journal for Numerical Methods in Biomedical Engineering}, -title = {A simple and efficient adaptive time stepping technique for low-order operator splitting schemes applied to cardiac electrophysiology}, -year = {2023} -volume = {39}, -number = {2}, -pages = {e3670}, -doi = {https://doi.org/10.1002/cnm.3670}, + author = {Ogiermann, Dennis and Perotti, Luigi E. and Balzani, Daniel}, + journal = {International Journal for Numerical Methods in Biomedical Engineering}, + title = {A simple and efficient adaptive time stepping technique for low-order operator splitting schemes applied to cardiac electrophysiology}, + year = {2023} + volume = {39}, + number = {2}, + pages = {e3670}, + doi = {https://doi.org/10.1002/cnm.3670}, } From e19fd50ee294149d156f6f488331e676c31bfa92 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Mon, 26 Aug 2024 16:32:13 +0200 Subject: [PATCH 30/37] access u from caches directly --- src/solver/adaptivity.jl | 10 +++++++--- src/solver/operator_splitting/integrator.jl | 4 ++-- src/solver/operator_splitting/solver.jl | 3 --- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/solver/adaptivity.jl b/src/solver/adaptivity.jl index c8615564..b298c28e 100644 --- a/src/solver/adaptivity.jl +++ b/src/solver/adaptivity.jl @@ -21,14 +21,18 @@ struct ReactionTangentController{LTG <: OS.LieTrotterGodunov, T <: Real} <: OS.A Δt_bounds::NTuple{2,T} end -mutable struct ReactionTangentControllerCache{T <: Real, LTGCache <: OS.LieTrotterGodunovCache} <: OS.AbstractOperatorSplittingCache +mutable struct ReactionTangentControllerCache{T <: Real, LTGCache <: OS.LieTrotterGodunovCache, uType} <: OS.AbstractOperatorSplittingCache const ltg_cache::LTGCache #It has Arrays so it can be const? + u::uType + uprev::uType # True previous solution Rₙ₊₁::T Rₙ::T + function ReactionTangentControllerCache(ltg_cache::LTGCache, Rₙ₊₁::T, Rₙ::T) where {T, LTGCache <: OS.LieTrotterGodunovCache} + uType = typeof(ltg_cache.u) + return new{T, LTGCache, uType}(ltg_cache, ltg_cache.u, ltg_cache.uprev, Rₙ₊₁, Rₙ) + end end -@inline OS.get_u(cache::ReactionTangentControllerCache) = OS.get_u(cache.ltg_cache) -@inline OS.get_uprev(cache::ReactionTangentControllerCache) = OS.get_uprev(cache.ltg_cache) @inline DiffEqBase.get_tmp_cache(integrator::OS.OperatorSplittingIntegrator, alg::OS.AbstractOperatorSplittingAlgorithm, cache::ReactionTangentControllerCache) = DiffEqBase.get_tmp_cache(integrator, alg, cache.ltg_cache) @inline function OS.advance_solution_to!(subintegrators::Tuple, cache::ReactionTangentControllerCache, tnext) diff --git a/src/solver/operator_splitting/integrator.jl b/src/solver/operator_splitting/integrator.jl index 4fb13690..7044467c 100644 --- a/src/solver/operator_splitting/integrator.jl +++ b/src/solver/operator_splitting/integrator.jl @@ -74,8 +74,8 @@ function DiffEqBase.__init( cache = init_cache(prob, alg; dt, kwargs...) - u = get_u(cache) - uprev = get_uprev(cache) + u = cache.u + uprev = cache.uprev subintegrators = build_subintegrators_recursive(prob.f, prob.f.synchronizers, p, cache, u, uprev, t0, dt, 1:length(u0), u, tstops, _tstops, saveat, _saveat) diff --git a/src/solver/operator_splitting/solver.jl b/src/solver/operator_splitting/solver.jl index e824a2f8..69d9640e 100644 --- a/src/solver/operator_splitting/solver.jl +++ b/src/solver/operator_splitting/solver.jl @@ -18,9 +18,6 @@ struct LieTrotterGodunovCache{uType, tmpType, iiType} <: AbstractOperatorSplitti inner_caches::iiType end -get_u(cache::LieTrotterGodunovCache) = cache.u -get_uprev(cache::LieTrotterGodunovCache) = cache.uprev - # Dispatch for outer construction function init_cache(prob::OperatorSplittingProblem, alg::LieTrotterGodunov; dt, kwargs...) # TODO @unpack f = prob From 68bd74c620392bf9979728844275f7fd934a9c42 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 30 Aug 2024 12:13:28 +0200 Subject: [PATCH 31/37] remove unroll filter --- src/Thunderbolt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Thunderbolt.jl b/src/Thunderbolt.jl index 5bdf7bed..4a4b1607 100644 --- a/src/Thunderbolt.jl +++ b/src/Thunderbolt.jl @@ -2,7 +2,7 @@ module Thunderbolt using TimerOutputs -import Unrolled: @unroll, unrolled_filter +import Unrolled: @unroll using Reexport, UnPack import LinearAlgebra: mul! From b7f6773ddf056a72e769b42a04f480abdb36e3aa Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 30 Aug 2024 12:14:42 +0200 Subject: [PATCH 32/37] adaptivity -> time/rtc --- src/Thunderbolt.jl | 1 - src/solver/{adaptivity.jl => time/rtc.jl} | 0 src/solver/time_integration.jl | 2 ++ 3 files changed, 2 insertions(+), 1 deletion(-) rename src/solver/{adaptivity.jl => time/rtc.jl} (100%) diff --git a/src/Thunderbolt.jl b/src/Thunderbolt.jl index 4a4b1607..b46e3e21 100644 --- a/src/Thunderbolt.jl +++ b/src/Thunderbolt.jl @@ -72,7 +72,6 @@ include("solver/interface.jl") include("solver/linear.jl") include("solver/nonlinear.jl") include("solver/time_integration.jl") -include("solver/adaptivity.jl") include("processing/ecg.jl") diff --git a/src/solver/adaptivity.jl b/src/solver/time/rtc.jl similarity index 100% rename from src/solver/adaptivity.jl rename to src/solver/time/rtc.jl diff --git a/src/solver/time_integration.jl b/src/solver/time_integration.jl index 0477f3b0..b77ff99e 100644 --- a/src/solver/time_integration.jl +++ b/src/solver/time_integration.jl @@ -2,3 +2,5 @@ include("time/time_integrator.jl") include("time/euler.jl") include("time/load_stepping.jl") include("time/partitioned_solver.jl") +include("time/rtc.jl") + From 934a4d30a2340c42a1535b64158fbed5b92226f5 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 30 Aug 2024 12:15:11 +0200 Subject: [PATCH 33/37] remove unused type --- src/solver/time/rtc.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/solver/time/rtc.jl b/src/solver/time/rtc.jl index b298c28e..0353e521 100644 --- a/src/solver/time/rtc.jl +++ b/src/solver/time/rtc.jl @@ -1,5 +1,3 @@ -abstract type AbstractTimeAdaptionAlgorithm end - """ ReactionTangentController{LTG <: OS.LieTrotterGodunov, T <: Real} <: OS.AbstractOperatorSplittingAlgorithm A timestep length controller for [`LieTrotterGodunov`](@ref) [Lie:1880:tti,Tro:1959:psg,God:1959:dmn](@cite) From 084302099312ca27bc2a8857b6aa71d3d57888c7 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 30 Aug 2024 12:15:44 +0200 Subject: [PATCH 34/37] yes --- src/solver/time/rtc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solver/time/rtc.jl b/src/solver/time/rtc.jl index 0353e521..e463e540 100644 --- a/src/solver/time/rtc.jl +++ b/src/solver/time/rtc.jl @@ -20,7 +20,7 @@ struct ReactionTangentController{LTG <: OS.LieTrotterGodunov, T <: Real} <: OS.A end mutable struct ReactionTangentControllerCache{T <: Real, LTGCache <: OS.LieTrotterGodunovCache, uType} <: OS.AbstractOperatorSplittingCache - const ltg_cache::LTGCache #It has Arrays so it can be const? + const ltg_cache::LTGCache u::uType uprev::uType # True previous solution Rₙ₊₁::T From c912e1a0ac515eeb3d4efc126b8381034e6558cd Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 30 Aug 2024 12:16:24 +0200 Subject: [PATCH 35/37] or not to do --- src/solver/time/rtc.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solver/time/rtc.jl b/src/solver/time/rtc.jl index e463e540..4e5f8405 100644 --- a/src/solver/time/rtc.jl +++ b/src/solver/time/rtc.jl @@ -87,7 +87,7 @@ end end # Dispatch for outer construction -function OS.init_cache(prob::OS.OperatorSplittingProblem, alg::ReactionTangentController; dt, kwargs...) # TODO +function OS.init_cache(prob::OS.OperatorSplittingProblem, alg::ReactionTangentController; dt, kwargs...) @unpack f = prob @assert f isa GenericSplitFunction From 3a64dcc7f76c3e6242c2e39d60ab05bcfe614166 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 30 Aug 2024 12:32:22 +0200 Subject: [PATCH 36/37] use only current R --- src/solver/time/rtc.jl | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/solver/time/rtc.jl b/src/solver/time/rtc.jl index 4e5f8405..6120e2fc 100644 --- a/src/solver/time/rtc.jl +++ b/src/solver/time/rtc.jl @@ -23,11 +23,10 @@ mutable struct ReactionTangentControllerCache{T <: Real, LTGCache <: OS.LieTrott const ltg_cache::LTGCache u::uType uprev::uType # True previous solution - Rₙ₊₁::T - Rₙ::T - function ReactionTangentControllerCache(ltg_cache::LTGCache, Rₙ₊₁::T, Rₙ::T) where {T, LTGCache <: OS.LieTrotterGodunovCache} + R::T + function ReactionTangentControllerCache(ltg_cache::LTGCache, R::T) where {T, LTGCache <: OS.LieTrotterGodunovCache} uType = typeof(ltg_cache.u) - return new{T, LTGCache, uType}(ltg_cache, ltg_cache.u, ltg_cache.uprev, Rₙ₊₁, Rₙ) + return new{T, LTGCache, uType}(ltg_cache, ltg_cache.u, ltg_cache.uprev, R) end end @@ -65,15 +64,14 @@ end end @inline function OS.stepsize_controller!(integrator::OS.OperatorSplittingIntegrator, alg::ReactionTangentController) - integrator.cache.Rₙ = integrator.cache.Rₙ₊₁ - integrator.cache.Rₙ₊₁ = get_reaction_tangent(integrator) + integrator.cache.R = get_reaction_tangent(integrator) return nothing end @inline function OS.step_accept_controller!(integrator::OS.OperatorSplittingIntegrator, alg::ReactionTangentController, q) - @unpack Rₙ₊₁, Rₙ = integrator.cache + @unpack R = integrator.cache @unpack σ_s, σ_c, Δt_bounds = alg - R = max(Rₙ, Rₙ₊₁) + if isinf(σ_s) && R == σ_c # Handle it? throw(ErrorException("Δt is undefined for R = σ_c where σ_s = ∞")) @@ -101,7 +99,7 @@ end # Dispatch for recursive construction function OS.construct_inner_cache(f::OS.AbstractOperatorSplitFunction, alg::ReactionTangentController, u::AbstractArray{T}, uprev::AbstractArray) where T <: Number ltg_cache = OS.construct_inner_cache(f, alg.ltg, u, uprev) - return ReactionTangentControllerCache(ltg_cache, zero(T), zero(T)) + return ReactionTangentControllerCache(ltg_cache, zero(T)) end function OS.build_subintegrators_recursive(f::GenericSplitFunction, synchronizers::Tuple, p::Tuple, cache::ReactionTangentControllerCache, u::AbstractArray, uprev::AbstractArray, t, dt, dof_range, uparent, tstops, _tstops, saveat, _saveat) From 5023e3870f227314500b234ff983fa53ab8ac993 Mon Sep 17 00:00:00 2001 From: Abdulaziz Date: Fri, 30 Aug 2024 12:52:14 +0200 Subject: [PATCH 37/37] inf --- src/solver/time/rtc.jl | 8 ++++---- test/test_integrators.jl | 3 ++- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/solver/time/rtc.jl b/src/solver/time/rtc.jl index 6120e2fc..fdb37850 100644 --- a/src/solver/time/rtc.jl +++ b/src/solver/time/rtc.jl @@ -72,11 +72,11 @@ end @unpack R = integrator.cache @unpack σ_s, σ_c, Δt_bounds = alg - if isinf(σ_s) && R == σ_c - # Handle it? - throw(ErrorException("Δt is undefined for R = σ_c where σ_s = ∞")) + if isinf(σ_s) + integrator._dt = R > σ_c ? Δt_bounds[1] : Δt_bounds[2] + else + integrator._dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] end - integrator._dt = (1 - 1/(1+exp((σ_c - R)*σ_s)))*(Δt_bounds[2] - Δt_bounds[1]) + Δt_bounds[1] return nothing end diff --git a/test/test_integrators.jl b/test/test_integrators.jl index 595a9fa6..be21344c 100644 --- a/test/test_integrators.jl +++ b/test/test_integrators.jl @@ -223,7 +223,8 @@ using UnPack @testset "σ_s = Inf, R = σ_c" begin timestepper_stepfunc_adaptive = Thunderbolt.ReactionTangentController(timestepper, Inf, 0.5, adaptive_tstep_range) integrator_stepfunc_adaptive = DiffEqBase.init(prob_force_half, timestepper_stepfunc_adaptive, dt=dt, verbose=true) - @test_throws ErrorException("Δt is undefined for R = σ_c where σ_s = ∞") DiffEqBase.solve!(integrator_stepfunc_adaptive) + DiffEqBase.solve!(integrator_stepfunc_adaptive) + @test integrator_stepfunc_adaptive._dt == timestepper_stepfunc_adaptive.Δt_bounds[2] end end end