[WIP] Relaxation Runge-Kutta

Project.toml

@@ -27,6 +27,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -43,6 +44,7 @@ SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
 ADTypes = "0.2, 1"

src/OrdinaryDiffEq.jl

@@ -77,6 +77,8 @@ import DiffEqBase: resize!, deleteat!, addat!, full_cache, user_cache, u_cache,
                    add_saveat!, set_reltol!,
                    set_abstol!, postamble!, last_step_failed,
                    isautodifferentiable
+
+export change_u!, change_dt!, apriori_bounds_dt
 
 using DiffEqBase: check_error!, @def, _vec, _reshape
 
@@ -128,6 +130,10 @@ import Preferences
 
 DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
 
+export AbstractPerformStepCallback, NoPerformStepCallback, PerformStepCallback
+
+include("performstep_callback.jl")
+
 include("doc_utils.jl")
 include("misc_utils.jl")
 
@@ -218,6 +224,7 @@ include("perform_step/prk_perform_step.jl")
 include("perform_step/pdirk_perform_step.jl")
 include("perform_step/dae_perform_step.jl")
 include("perform_step/qprk_perform_step.jl")
+include("perform_step/test_for_relaxation_perform_step.jl")
 
 include("dense/generic_dense.jl")
 include("dense/interpolants.jl")
@@ -377,7 +384,7 @@ export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3,
        OwrenZen5,
        BS3, BS5, RK46NL, DP5, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8,
        Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65, PFRK87,
-       RKM, MSRK5, MSRK6, Stepanov5, SIR54, QPRK98, PSRK4p7q6, PSRK3p6q5, PSRK3p5q4
+       RKM, MSRK5, MSRK6, Stepanov5, SIR54, QPRK98, PSRK4p7q6, PSRK3p6q5, PSRK3p5q4, Tsit5_for_relaxation
 
 export SSPRK22, SSPRK33, KYKSSPRK42, SSPRK53, SSPRK53_2N1, SSPRK53_2N2, SSPRK53_H, SSPRK63,
        SSPRK73, SSPRK83, SSPRK43, SSPRK432,

src/alg_utils.jl

@@ -588,6 +588,7 @@ alg_order(alg::OwrenZen4) = 4
 alg_order(alg::OwrenZen5) = 5
 alg_order(alg::DP5) = 5
 alg_order(alg::Tsit5) = 5
+alg_order(alg::Tsit5_for_relaxation) = 5
 alg_order(alg::DP8) = 8
 alg_order(alg::Vern6) = 6
 alg_order(alg::Vern7) = 7

src/algorithms/explicit_rk.jl

@@ -361,6 +361,19 @@ function Tsit5(stage_limiter!, step_limiter! = trivial_limiter!)
     Tsit5(stage_limiter!, step_limiter!, False())
 end
 
+Base.@kwdef struct Tsit5_for_relaxation{StageLimiter, StepLimiter, Thread} <:
+    OrdinaryDiffEqAdaptiveAlgorithm
+stage_limiter!::StageLimiter = trivial_limiter!
+step_limiter!::StepLimiter = trivial_limiter!
+thread::Thread = False()
+end
+
+# for backwards compatibility
+function Tsit5_for_relaxation(stage_limiter!, step_limiter! = trivial_limiter!)
+Tsit5_for_relaxation(stage_limiter!, step_limiter!, False())
+end
+
+
 @doc explicit_rk_docstring(
     "Hairer's 8/5/3 adaption of the Dormand-Prince Runge-Kutta method. (7th order interpolant).",
     "DP8",

src/caches/low_order_rk_caches.jl

@@ -448,6 +448,44 @@ if isdefined(Base, :Experimental) && isdefined(Base.Experimental, :silence!)
     Base.Experimental.silence!(Tsit5Cache)
 end
 
+@cache struct Tsit5Cache_for_relaxation{uType, rateType, uNoUnitsType, StageLimiter, StepLimiter,
+    Thread} <: OrdinaryDiffEqMutableCache
+    u::uType
+    uprev::uType
+    k1::rateType
+    k2::rateType
+    k3::rateType
+    k4::rateType
+    k5::rateType
+    k6::rateType
+    k7::rateType
+    utilde::uType
+    tmp::uType
+    atmp::uNoUnitsType
+    stage_limiter!::StageLimiter
+    step_limiter!::StepLimiter
+    thread::Thread
+end
+
+function alg_cache(alg::Tsit5_for_relaxation, u, rate_prototype, ::Type{uEltypeNoUnits},
+    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+    dt, reltol, p, calck,
+    ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    k1 = zero(rate_prototype)
+    k2 = zero(rate_prototype)
+    k3 = zero(rate_prototype)
+    k4 = zero(rate_prototype)
+    k5 = zero(rate_prototype)
+    k6 = zero(rate_prototype)
+    k7 = zero(rate_prototype)
+    utilde = zero(u)
+    atmp = similar(u, uEltypeNoUnits)
+    recursivefill!(atmp, false)
+    tmp = zero(u)
+    Tsit5Cache_for_relaxation(u, uprev, k1, k2, k3, k4, k5, k6, k7, utilde, tmp, atmp,
+        alg.stage_limiter!, alg.step_limiter!, alg.thread)
+end
+
 @cache struct RK46NLCache{uType, rateType, TabType, StageLimiter, StepLimiter, Thread} <:
               OrdinaryDiffEqMutableCache
     u::uType
@@ -547,6 +585,13 @@ function alg_cache(alg::Tsit5, u, rate_prototype, ::Type{uEltypeNoUnits},
     Tsit5ConstantCache()
 end
 
+function alg_cache(alg::Tsit5_for_relaxation, u, rate_prototype, ::Type{uEltypeNoUnits},
+    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+    dt, reltol, p, calck,
+    ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+Tsit5ConstantCache_for_relaxation()
+end
+
 @cache struct DP5Cache{uType, rateType, uNoUnitsType, StageLimiter, StepLimiter,
     Thread} <: OrdinaryDiffEqMutableCache
     u::uType

src/integrators/integrator_interface.jl

@@ -506,3 +506,22 @@ function DiffEqBase.set_u!(integrator::ODEIntegrator, u)
 end
 
 DiffEqBase.has_stats(i::ODEIntegrator) = true
+
+
+function change_dt!(integrator::ODEIntegrator, dt)
+    integrator.dt_has_changed = true
+    integrator.dt_changed = dt
+end
+
+function change_u!(integrator::ODEIntegrator, u)
+    integrator.u = u
+end
+
+function change_fsallast!(integrator::ODEIntegrator, fsallast)
+    integrator.fsallast = fsallast
+end
+
+has_poststep_callback(integrator::ODEIntegrator) = has_poststep_callback(integrator.opts.performstepcallback)
+has_postfsal_callback(integrator::ODEIntegrator) = has_postfsal_callback(integrator.opts.performstepcallback)
+has_postEEst_callback(integrator::ODEIntegrator) = has_postEEst_callback(integrator.opts.performstepcallback)
+

src/integrators/type.jl

@@ -1,6 +1,6 @@
 mutable struct DEOptions{absType, relType, QT, tType, Controller, F1, F2, F3, F4, F5, F6,
     F7, tstopsType, discType, ECType, SType, MI, tcache, savecache,
-    disccache}
+    disccache, MT}
     maxiters::MI
     save_everystep::Bool
     adaptive::Bool
@@ -46,6 +46,7 @@ mutable struct DEOptions{absType, relType, QT, tType, Controller, F1, F2, F3, F4
     force_dtmin::Bool
     advance_to_tstop::Bool
     stop_at_next_tstop::Bool
+    performstepcallback::AbstractPerformStepCallback
 end
 
 TruncatedStacktraces.@truncate_stacktrace DEOptions
@@ -134,6 +135,10 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori
     stats::SciMLBase.DEStats
     initializealg::IA
     differential_vars::DV
+
+    dt_has_changed::Bool
+    dt_changed::tType
+
     fsalfirst::FSALType
     fsallast::FSALType
 
@@ -155,7 +160,7 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori
             accept_step, isout, reeval_fsal, u_modified,
             reinitialize, isdae,
             opts, stats,
-            initializealg, differential_vars) where {algType, IIP, uType,
+            initializealg, differential_vars, dt_has_changed, dt_changed) where {algType, IIP, uType,
             duType, tType, pType,
             eigenType, EEstT,
             tTypeNoUnits, tdirType,
@@ -177,7 +182,7 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori
             do_error_check,
             event_last_time, vector_event_last_time, last_event_error,
             accept_step, isout, reeval_fsal, u_modified, reinitialize, isdae,
-            opts, stats, initializealg, differential_vars) # Leave off fsalfirst and last
+            opts, stats, initializealg, differential_vars, dt_has_changed, dt_changed) # Leave off fsalfirst and last
     end
 end