Rebase Relaxation-Controller

Project.toml

@@ -121,8 +121,8 @@ OrdinaryDiffEqQPRK = "1"
 OrdinaryDiffEqRKN = "1"
 OrdinaryDiffEqRosenbrock = "1"
 OrdinaryDiffEqSDIRK = "1"
-OrdinaryDiffEqStabilizedIRK = "1"
 OrdinaryDiffEqSSPRK = "1"
+OrdinaryDiffEqStabilizedIRK = "1"
 OrdinaryDiffEqStabilizedRK = "1"
 OrdinaryDiffEqSymplecticRK = "1"
 OrdinaryDiffEqTsit5 = "1"

lib/OrdinaryDiffEqCore/src/integrators/controllers.jl

@@ -5,6 +5,11 @@ using OrdinaryDiffEqCore
     stepsize_controller!(integrator, integrator.opts.controller, alg)
 end
 
+# checks whether the controller should accept a step based on the error estimate
+@inline function accept_step_controller(integrator, controller::AbstractController)
+    return integrator.EEst <= 1
+end
+
 @inline function step_accept_controller!(integrator, alg, q)
     step_accept_controller!(integrator, integrator.opts.controller, alg, q)
 end
@@ -23,6 +28,8 @@ reset_alg_dependent_opts!(controller::AbstractController, alg1, alg2) = nothing
 
 DiffEqBase.reinit!(integrator::ODEIntegrator, controller::AbstractController) = nothing
 
+@inline next_time_controller(::ODEIntegrator, ::AbstractController, ttmp, dt) = ttmp
+
 # Standard integral (I) step size controller
 """
     IController()
@@ -60,6 +67,30 @@ end
 struct DummyController <: AbstractController
 end
 
+
+"""
+    NonAdaptiveController()
+This Controller exists to match the interface when one does not want to use a controller,
+basically if you want to keep a fixed time step.
+"""
+struct NonAdaptiveController <: AbstractController
+end
+
+@inline function stepsize_controller!(integrator, ::NonAdaptiveController, alg)
+    nothing
+end
+
+@inline function accept_step_controller(integrator, ::NonAdaptiveController)
+    return true
+end
+
+function step_accept_controller!(integrator, ::NonAdaptiveController, alg, q)
+    integrator.dt
+end
+
+function step_reject_controller!(integrator, ::NonAdaptiveController, alg)
+end
+
 @inline function stepsize_controller!(integrator, controller::IController, alg)
     @unpack qmin, qmax, gamma = integrator.opts
     EEst = DiffEqBase.value(integrator.EEst)
@@ -311,11 +342,6 @@ end
     return dt_factor
 end
 
-# checks whether the controller should accept a step based on the error estimate
-@inline function accept_step_controller(integrator, controller::AbstractController)
-    return integrator.EEst <= 1
-end
-
 @inline function accept_step_controller(integrator, controller::PIDController)
     return integrator.qold >= controller.accept_safety
 end
@@ -400,3 +426,85 @@ function post_newton_controller!(integrator, alg)
     integrator.dt = integrator.dt / integrator.opts.failfactor
     nothing
 end
+
+# Relaxation step size controller
+"""
+    RelaxationController(controller, T)
+
+Controller to perform a relaxation on a step of a Runge-Kuttas method.
+
+## References
+ - Sebastian Bleecke, Hendrik Ranocha (2023)
+    Step size control for explicit relaxation Runge-Kutta methods preserving invariants
+    [DOI: 10.1145/641876.641877](https://doi.org/10.1145/641876.641877)
+"""
+mutable struct RelaxationController{CON, T} <: AbstractController
+    controller::CON
+    gamma::T
+    function RelaxationController(controller::AbstractController, T=Float64)
+        new{typeof(controller), T}(controller, one(T))
+    end
+end
+
+mutable struct Relaxation{INV, T}
+    invariant::INV
+    gamma_min::T
+    gamma_max::T
+    function Relaxation(invariant, gamma_min = 4/5, gamma_max = 6/5)
+        new{typeof(invariant), typeof(gamma_min)}(invariant, gamma_min, gamma_max)
+    end
+end
+
+
+function _relaxation_nlsolve(gamma, p)
+    u, uprev, invariant = p
+    invariant(@.. gamma * u + (1-gamma) * uprev) .- invariant(uprev)
+end
+
+@muladd function (r::Relaxation)(integrator)
+    @unpack t, dt, uprev, u = integrator
+    @unpack invariant, gamma_min, gamma_max = integrator.opts.relaxation
+    gamma = one(t)
+    p = (u, uprev, invariant)
+    if _relaxation_nlsolve(gamma_min, p) * _relaxation_nlsolve(gamma_max, p) ≤ 0
+        gamma = solve(IntervalNonlinearProblem{false}(_relaxation_nlsolve, (gamma_min, gamma_max), p),
+            DiffEqBase.InternalITP()).left
+    end
+    gamma
+end
+
+function Base.show(io::IO, controller::RelaxationController)
+    print(io, controller.controller)
+    print(io, "\n Relaxation parameters = ", controller.gamma)
+end
+
+@inline function next_time_controller(integrator::ODEIntegrator, controller::RelaxationController, ttmp, dt)
+    gamma = integrator.opts.relaxation(integrator)
+    integrator.dt *= oftype(dt, gamma)
+    vdt = integrator.dt
+    modify_dt_for_tstops!(integrator)
+    if integrator.dt != vdt
+        gamma = integrator.dt/dt
+    end
+    @. integrator.u =  integrator.uprev + gamma*(integrator.u .- integrator.uprev)
+    @. integrator.fsallast = integrator.fsalfirst + gamma*(integrator.fsallast - integrator.fsalfirst)
+    controller.gamma = gamma
+    ttmp + integrator.dt - dt
+end
+
+@inline function stepsize_controller!(integrator, controller::RelaxationController, alg)
+    stepsize_controller!(integrator, controller.controller, alg)
+end
+
+@inline function accept_step_controller(integrator, controller::RelaxationController)
+    accept_step_controller(integrator, controller.controller)
+end
+
+function step_accept_controller!(integrator, controller::RelaxationController, alg, dt_factor)
+    step_accept_controller!(integrator, controller.controller, alg, dt_factor)
+end
+
+function step_reject_controller!(integrator, controller::RelaxationController, alg)
+    integrator.dt = integrator.dt * integrator.qold / controller.gamma
+end
+

lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl

@@ -238,6 +238,7 @@ function _loopfooter!(integrator)
                 q)) *
                     oneunit(integrator.dt)
             integrator.tprev = integrator.t
+            ttmp = next_time_controller(integrator, integrator.opts.controller, ttmp, integrator.dt)
             integrator.t = fixed_t_for_floatingpoint_error!(integrator, ttmp)
             calc_dt_propose!(integrator, dtnew)
             handle_callbacks!(integrator)

lib/OrdinaryDiffEqCore/src/integrators/type.jl

@@ -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
+    relaxation
 end
 
 """

lib/OrdinaryDiffEqCore/src/solve.jl

@@ -72,6 +72,7 @@ function DiffEqBase.__init(
         alias_u0 = false,
         alias_du0 = false,
         initializealg = DefaultInit(),
+        relaxation = nothing,
         kwargs...) where {recompile_flag}
     if prob isa DiffEqBase.AbstractDAEProblem && alg isa OrdinaryDiffEqAlgorithm
         error("You cannot use an ODE Algorithm with a DAEProblem")
@@ -366,6 +367,8 @@ function DiffEqBase.__init(
         controller = default_controller(_alg, cache, qoldinit, beta1, beta2)
     end
 
+    controller = relaxation !== nothing ? RelaxationController(controller, eltype(u)) : controller
+
     save_end_user = save_end
     save_end = save_end === nothing ?
                save_everystep || isempty(saveat) || saveat isa Number ||
@@ -414,7 +417,8 @@ function DiffEqBase.__init(
         unstable_check,
         verbose, calck, force_dtmin,
         advance_to_tstop,
-        stop_at_next_tstop)
+        stop_at_next_tstop,
+        relaxation)
 
     stats = SciMLBase.DEStats(0)
     differential_vars = prob isa DAEProblem ? prob.differential_vars :

lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl

@@ -11,7 +11,7 @@ import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, qsteady_max_default,
                            calculate_residuals, calculate_residuals!,
                            get_current_adaptive_order, get_fsalfirstlast,
                            ode_interpolant, ode_interpolant!, trivial_limiter!,
-                           generic_solver_docstring
+                           generic_solver_docstring, DummyController
 using MuladdMacro, FastBroadcast, RecursiveArrayTools
 import LinearAlgebra: rmul!
 import Static: False

lib/OrdinaryDiffEqNordsieck/src/controllers.jl

@@ -1,6 +1,3 @@
-struct DummyController <: AbstractController
-end
-
 # JVODE
 function stepsize_controller!(integrator, alg::JVODE)
     if iszero(integrator.EEst)

src/OrdinaryDiffEq.jl

@@ -58,7 +58,8 @@ import OrdinaryDiffEqCore: trivial_limiter!, CompositeAlgorithm, alg_order,
                            _change_t_via_interpolation!, ODEIntegrator, _ode_interpolant!,
                            current_interpolant, resize_nlsolver!, _ode_interpolant,
                            handle_tstop!, _postamble!, update_uprev!, resize_J_W!,
-                           DAEAlgorithm, get_fsalfirstlast, strip_cache
+                           DAEAlgorithm, get_fsalfirstlast, strip_cache, Relaxation,
+                           RelaxationController, NonAdaptiveController
 
 export CompositeAlgorithm, ShampineCollocationInit, BrownFullBasicInit, NoInit
 AutoSwitch
@@ -228,5 +229,5 @@ export AutoSwitch, AutoTsit5, AutoDP5,
        AutoVern6, AutoVern7, AutoVern8, AutoVern9
 
 import OrdinaryDiffEqCore: IController, PIController, PIDController
-export IController, PIController, PIDController
+export IController, PIController, PIDController, NonAdaptiveController, Relaxation, RelaxationController
 end # module

test/relaxation/harmonic_oscillator.jl

@@ -0,0 +1,23 @@
+using OrdinaryDiffEq, DiffEqDevTools, LinearAlgebra, Test
+
+@testset "Harmonic oscillator" begin
+
+dts = (1 / 2) .^ (6:-1:4)
+
+f = (u, p, t) -> [-u[2],u[1]]
+prob = ODEProblem(
+    ODEFunction(f; analytic = (u0, p, t) -> [cos(t), sin(t)]),
+    [1.0, 0.0],
+    (0.0, 1.0))
+
+# Convergence with the old method Tsit5()
+sim = test_convergence(dts, prob, Tsit5(), adaptive = true)
+@test sim.𝒪est[:final] ≈ 5 atol=0.2
+
+# Convergence with relaxation with FSAL-R, i.e  f(uᵧ,ₙ₊₁) ≈ f(uᵧ,ₙ) + γ ( f(uₙ₊₁) - f(uᵧ,ₙ))
+relaxation = Relaxation(norm)
+controller = RelaxationController(NonAdaptiveController())
+sim_relax = test_convergence(dts, prob, Tsit5(); relaxation, controller)
+@test sim_relax.𝒪est[:final] ≈ 5 atol=0.2
+
+end

test/relaxation/non_linear_oscillator.jl

@@ -0,0 +1,24 @@
+using OrdinaryDiffEq, DiffEqDevTools, LinearAlgebra, Test
+
+@testset "Non linear harmonic oscillator" begin
+
+dts = (1 / 2) .^ (6:-1:4)
+
+f = (u, p, t) -> [-u[2]/(u[1]^2 + u[2]^2),u[1]/(u[1]^2 + u[2]^2)]
+prob = ODEProblem(
+    ODEFunction(f; analytic = (u0, p, t) -> [cos(t), sin(t)]),
+    [1.0, 0.0],
+    (0.0, 1.0))
+
+
+# Convergence with the method Tsit5()
+sim = test_convergence(dts, prob, Tsit5())
+@test sim.𝒪est[:final] ≈ 5.4 atol=0.2
+
+# Convergence with relaxation with FSAL-R, i.e  f(uᵧ,ₙ₊₁) ≈ f(uᵧ,ₙ) + γ ( f(uₙ₊₁) - f(uᵧ,ₙ))
+relaxation = Relaxation(norm)
+controller = RelaxationController(NonAdaptiveController())
+sim_relax = test_convergence(dts, prob, Tsit5(); relaxation, controller, adaptive=true)
+@test sim.𝒪est[:final] ≈ 5.4 atol=0.2
+
+end

test/relaxation/non_linear_pendulum.jl

@@ -0,0 +1,24 @@
+using OrdinaryDiffEq, DiffEqDevTools
+
+printstyled("Non Linear Pendulum\n"; bold = true)
+
+dts = (1 / 2) .^ (6:-1:4)
+
+f = (u, p, t) -> [-sin(u[2]), u[1]]
+prob = ODEProblem(
+    f,
+    [1.0, 0.0],
+    (0.0, 1.0))
+
+invariant(x) = x[1]^2/2 -  cos(x[2])
+
+test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
+
+# Convergence with the method Tsit5()
+sim = analyticless_test_convergence(dts, prob, Tsit5(), test_setup)
+println("order of convergence of older perform_step! : "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with FSAL-R, i.e  f(uᵧ,ₙ₊₁) ≈ f(uᵧ,ₙ) + γ ( f(uₙ₊₁) - f(uᵧ,ₙ)) 
+r = Relaxation(invariant)
+sim = analyticless_test_convergence(dts, prob, Tsit5(), test_setup; relaxation = r)
+println("order with relaxation with FSAL-R modification: "*string(sim.𝒪est[:final]))

test/relaxation/test.jl

@@ -0,0 +1,13 @@
+using OrdinaryDiffEq, DiffEqDevTools,  Test, BenchmarkTools, LinearAlgebra
+
+f = (u, p, t) -> [-u[2],u[1]]
+prob = ODEProblem(
+    ODEFunction(f; analytic = (u0, p, t) -> [cos(t), sin(t)]),
+                [1.0, 0.0],
+                (0.0, 1.0))
+invariant(x) = norm(x)
+
+r = Relaxation(invariant)
+
+sol_relax = solve(prob, Tsit5(); relaxation = r)
+sol = solve(prob, Tsit5())

test/relaxation/time_dependant_harmonic_oscillator.jl

@@ -0,0 +1,25 @@
+using OrdinaryDiffEq, DiffEqDevTools, LinearAlgebra, Test
+
+
+@testset "Time-dependent harmonic oscillator" begin
+
+dts = (1 / 2) .^ (6:-1:4)
+
+f = (u, p, t) -> [-(1+0.5 * sin(t))*u[2], (1+0.5 * sin(t))*u[1]]
+prob = ODEProblem(
+    ODEFunction(f;
+                analytic = (u0, p, t)->[cos(0.5)*cos(t-0.5*cos(t))-sin(0.5)*sin(t-0.5*cos(t)),
+                                        sin(0.5)*cos(t-0.5*cos(t))+cos(0.5)*sin(t-0.5*cos(t))]),
+    [1.0, 0.0],
+    (0.0, 1.0))
+
+# Convergence with the method Tsit5()
+sim = test_convergence(dts, prob, Tsit5())
+@test sim.𝒪est[:final] ≈ 5.2 atol=0.2
+
+# Convergence with relaxation with FSAL-R, i.e  f(uᵧ,ₙ₊₁) ≈ f(uᵧ,ₙ) + γ ( f(uₙ₊₁) - f(uᵧ,ₙ))
+relaxation = Relaxation(norm)
+controller = RelaxationController(NonAdaptiveController())
+sim_relax = test_convergence(dts, prob, Tsit5(); relaxation, controller, adaptive=true)
+@test sim.𝒪est[:final] ≈ 5.2 atol=0.2
+end