From 6a12462c8f49c3e7889ad2c85eb73758843cc7f6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Th=C3=A9o=20Duez?=
 <42682941+Theozeud@users.noreply.github.com>
Date: Tue, 9 Jul 2024 16:58:06 +0200
Subject: [PATCH 1/4] Relaxation implementation

---
 Project.toml                                  |  1 +
 .../src/integrators/controllers.jl            | 83 +++++++++++++++++++
 .../src/integrators/integrator_utils.jl       |  5 +-
 .../src/integrators/type.jl                   |  1 +
 lib/OrdinaryDiffEqCore/src/solve.jl           |  6 +-
 test/relaxation/harmonic_oscillator.jl        | 23 +++++
 test/relaxation/non_linear_oscillator.jl      | 22 +++++
 test/relaxation/non_linear_pendulum.jl        | 24 ++++++
 test/relaxation/test.jl                       | 13 +++
 .../time_dependant_harmonic_oscillator.jl     | 24 ++++++
 10 files changed, 199 insertions(+), 3 deletions(-)
 create mode 100644 test/relaxation/harmonic_oscillator.jl
 create mode 100644 test/relaxation/non_linear_oscillator.jl
 create mode 100644 test/relaxation/non_linear_pendulum.jl
 create mode 100644 test/relaxation/test.jl
 create mode 100644 test/relaxation/time_dependant_harmonic_oscillator.jl

diff --git a/Project.toml b/Project.toml
index 2410fd74ca..5a90dbbd31 100644
--- a/Project.toml
+++ b/Project.toml
@@ -62,6 +62,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
index a13c3daf88..cad6fcd82b 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
+++ b/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()
@@ -400,3 +407,79 @@ 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)
+        new{typeof(controller), T}(controller, one(T))
+    end
+end
+
+mutable struct Relaxation{INV, OPT, T}
+    invariant::INV
+    opt::OPT
+    gamma_min::T
+    gamma_max::T
+    function Relaxation(invariant, opt = AlefeldPotraShi, gamma_min = 4//5, gamma_max = 6//5)
+        new{typeof(invariant), typeof(opt), typeof(gamma_min)}(invariant, opt, gamma_min, gamma_max)
+    end
+end
+
+@muladd function (r::Relaxation)(integrator)
+    @unpack t, dt, uprev, u = integrator
+    @unpack invariant, opt, gamma_min, gamma_max = integrator.opts.relaxation
+    gamma = one(t)
+    S_u = u .- uprev
+    if (invariant(gamma_min * S_u .+ uprev) .- invariant(uprev)) * (invariant(gamma_max * S_u .+ uprev) .- invariant(uprev)) ≤ 0 
+        gamma = find_zero(gamma -> invariant(gamma*S_u .+ uprev) .- invariant(uprev), (gamma_min, gamma_max), opt())
+    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
+
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index fc56ed8d50..848646a461 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -91,7 +91,7 @@ function _savevalues!(integrator, force_save, reduce_size)::Tuple{Bool, Bool}
                 copyat_or_push!(integrator.sol.u, integrator.saveiter,
                     integrator.u[integrator.opts.save_idxs], false)
             end
-            if isdiscretealg(integrator.alg) || integrator.opts.dense
+            if integrator.alg isa FunctionMap || integrator.opts.dense
                 integrator.saveiter_dense += 1
                 if integrator.opts.dense
                     if integrator.opts.save_idxs === nothing
@@ -124,7 +124,7 @@ function _savevalues!(integrator, force_save, reduce_size)::Tuple{Bool, Bool}
                 integrator.u[integrator.opts.save_idxs], false)
         end
         copyat_or_push!(integrator.sol.t, integrator.saveiter, integrator.t)
-        if isdiscretealg(integrator.alg) || integrator.opts.dense
+        if integrator.alg isa FunctionMap || integrator.opts.dense
             integrator.saveiter_dense += 1
             if integrator.opts.dense
                 if integrator.opts.save_idxs === nothing
@@ -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)
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/type.jl b/lib/OrdinaryDiffEqCore/src/integrators/type.jl
index 1f9bef031d..17752dd0e9 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/type.jl
+++ b/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
 
 """
diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl
index a5efe16d43..1ba0052e8e 100644
--- a/lib/OrdinaryDiffEqCore/src/solve.jl
+++ b/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 :
diff --git a/test/relaxation/harmonic_oscillator.jl b/test/relaxation/harmonic_oscillator.jl
new file mode 100644
index 0000000000..74f7c2f6ba
--- /dev/null
+++ b/test/relaxation/harmonic_oscillator.jl
@@ -0,0 +1,23 @@
+using OrdinaryDiffEq, DiffEqDevTools, LinearAlgebra
+
+printstyled("Harmonic Oscillator\n"; bold = true)
+
+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))
+
+invariant(x) = norm(x)
+
+# Convergence with the old method Tsit5()
+sim = test_convergence(dts, prob, Tsit5(), adaptive = true)
+println("order of convergence of Tsit5 without relaxation : "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with FSAL-R, i.e  f(uᵧ,ₙ₊₁) ≈ f(uᵧ,ₙ) + γ ( f(uₙ₊₁) - f(uᵧ,ₙ)) 
+r = Relaxation(invariant)
+sim_relax = test_convergence(dts, prob, Tsit5(); relaxation = r, adaptive = true)
+println("order with relaxation with FSAL-R modification: "*string(sim_relax.𝒪est[:final]))
+
diff --git a/test/relaxation/non_linear_oscillator.jl b/test/relaxation/non_linear_oscillator.jl
new file mode 100644
index 0000000000..7a156e40b2
--- /dev/null
+++ b/test/relaxation/non_linear_oscillator.jl
@@ -0,0 +1,22 @@
+using OrdinaryDiffEq, DiffEqDevTools
+
+printstyled("Non linear harmonic oscillator\n"; bold = true)
+
+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))
+
+invariant(x) = norm(x)
+
+# Convergence with the method Tsit5()
+sim = test_convergence(dts, prob, Tsit5())
+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 = test_convergence(dts, prob, Tsit5(); relaxation = r)
+println("order with relaxation with FSAL-R modification: "*string(sim.𝒪est[:final]))
diff --git a/test/relaxation/non_linear_pendulum.jl b/test/relaxation/non_linear_pendulum.jl
new file mode 100644
index 0000000000..2f759ae113
--- /dev/null
+++ b/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]))
\ No newline at end of file
diff --git a/test/relaxation/test.jl b/test/relaxation/test.jl
new file mode 100644
index 0000000000..c8ad871f3c
--- /dev/null
+++ b/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())
\ No newline at end of file
diff --git a/test/relaxation/time_dependant_harmonic_oscillator.jl b/test/relaxation/time_dependant_harmonic_oscillator.jl
new file mode 100644
index 0000000000..370c65bf83
--- /dev/null
+++ b/test/relaxation/time_dependant_harmonic_oscillator.jl
@@ -0,0 +1,24 @@
+using OrdinaryDiffEq, DiffEqDevTools
+
+printstyled("Time-dependent harmonic oscillator\n"; bold = true)
+
+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))
+
+invariant(x) = norm(x)
+
+# Convergence with the method Tsit5()
+sim = test_convergence(dts, prob, Tsit5(), adaptative = true)
+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 = test_convergence(dts, prob, Tsit5(); relaxation = r, adaptative = true)
+println("order with relaxation with FSAL-R modification: "*string(sim.𝒪est[:final]))
\ No newline at end of file

From b3c8e19548dda330eb7a7741b4ee74600dc662ed Mon Sep 17 00:00:00 2001
From: oscarddssmith <oscar.smith@juliacomputing.com>
Date: Tue, 3 Sep 2024 18:09:42 -0400
Subject: [PATCH 2/4] fixes

---
 Project.toml                                  |  3 +-
 .../src/integrators/controllers.jl            | 29 ++++++++++---------
 .../src/integrators/integrator_utils.jl       |  4 +--
 src/OrdinaryDiffEq.jl                         |  4 +--
 4 files changed, 20 insertions(+), 20 deletions(-)

diff --git a/Project.toml b/Project.toml
index 5a90dbbd31..abfe81bbc6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -62,7 +62,6 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
-Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226"
@@ -122,8 +121,8 @@ OrdinaryDiffEqQPRK = "1"
 OrdinaryDiffEqRKN = "1"
 OrdinaryDiffEqRosenbrock = "1"
 OrdinaryDiffEqSDIRK = "1"
-OrdinaryDiffEqStabilizedIRK = "1"
 OrdinaryDiffEqSSPRK = "1"
+OrdinaryDiffEqStabilizedIRK = "1"
 OrdinaryDiffEqStabilizedRK = "1"
 OrdinaryDiffEqSymplecticRK = "1"
 OrdinaryDiffEqTsit5 = "1"
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
index cad6fcd82b..63966d4b6b 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
@@ -318,11 +318,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
@@ -412,7 +407,7 @@ end
 """
     RelaxationController(controller, T)
 
-Controller to perform a relaxation on a step of a Runge-Kuttas method.  
+Controller to perform a relaxation on a step of a Runge-Kuttas method.
 
 ## References
  - Sebastian Bleecke, Hendrik Ranocha (2023)
@@ -427,23 +422,29 @@ mutable struct RelaxationController{CON, T} <: AbstractController
     end
 end
 
-mutable struct Relaxation{INV, OPT, T}
+mutable struct Relaxation{INV, T}
     invariant::INV
-    opt::OPT
     gamma_min::T
     gamma_max::T
-    function Relaxation(invariant, opt = AlefeldPotraShi, gamma_min = 4//5, gamma_max = 6//5)
-        new{typeof(invariant), typeof(opt), typeof(gamma_min)}(invariant, opt, gamma_min, gamma_max)
+    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, opt, gamma_min, gamma_max = integrator.opts.relaxation
+    @unpack invariant, gamma_min, gamma_max = integrator.opts.relaxation
     gamma = one(t)
-    S_u = u .- uprev
-    if (invariant(gamma_min * S_u .+ uprev) .- invariant(uprev)) * (invariant(gamma_max * S_u .+ uprev) .- invariant(uprev)) ≤ 0 
-        gamma = find_zero(gamma -> invariant(gamma*S_u .+ uprev) .- invariant(uprev), (gamma_min, gamma_max), opt())
+    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
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index 848646a461..58f599c34e 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -91,7 +91,7 @@ function _savevalues!(integrator, force_save, reduce_size)::Tuple{Bool, Bool}
                 copyat_or_push!(integrator.sol.u, integrator.saveiter,
                     integrator.u[integrator.opts.save_idxs], false)
             end
-            if integrator.alg isa FunctionMap || integrator.opts.dense
+            if isdiscretealg(integrator.alg) || integrator.opts.dense
                 integrator.saveiter_dense += 1
                 if integrator.opts.dense
                     if integrator.opts.save_idxs === nothing
@@ -124,7 +124,7 @@ function _savevalues!(integrator, force_save, reduce_size)::Tuple{Bool, Bool}
                 integrator.u[integrator.opts.save_idxs], false)
         end
         copyat_or_push!(integrator.sol.t, integrator.saveiter, integrator.t)
-        if integrator.alg isa FunctionMap || integrator.opts.dense
+        if isdiscretealg(integrator.alg) || integrator.opts.dense
             integrator.saveiter_dense += 1
             if integrator.opts.dense
                 if integrator.opts.save_idxs === nothing
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
index 503484ac33..01d8301d78 100644
--- a/src/OrdinaryDiffEq.jl
+++ b/src/OrdinaryDiffEq.jl
@@ -58,7 +58,7 @@ 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
 
 export CompositeAlgorithm, ShampineCollocationInit, BrownFullBasicInit, NoInit
 AutoSwitch
@@ -228,5 +228,5 @@ export AutoSwitch, AutoTsit5, AutoDP5,
        AutoVern6, AutoVern7, AutoVern8, AutoVern9
 
 import OrdinaryDiffEqCore: IController, PIController, PIDController
-export IController, PIController, PIDController
+export IController, PIController, PIDController, Relaxation
 end # module

From b7f895c2634525c22021c220760d3a3d4c4c2c43 Mon Sep 17 00:00:00 2001
From: oscarddssmith <oscar.smith@juliacomputing.com>
Date: Tue, 3 Sep 2024 18:51:11 -0400
Subject: [PATCH 3/4] add NonAdaptiveController

---
 .../src/integrators/controllers.jl            | 24 +++++++++++++++++++
 .../src/OrdinaryDiffEqNordsieck.jl            |  2 +-
 .../src/controllers.jl                        |  3 ---
 src/OrdinaryDiffEq.jl                         |  5 ++--
 4 files changed, 28 insertions(+), 6 deletions(-)

diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
index 63966d4b6b..8cd99c394f 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
@@ -67,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)
diff --git a/lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl b/lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl
index fa8b873e42..b346fd1795 100644
--- a/lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl
+++ b/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
diff --git a/lib/OrdinaryDiffEqNordsieck/src/controllers.jl b/lib/OrdinaryDiffEqNordsieck/src/controllers.jl
index 216eb36168..c74c50f736 100644
--- a/lib/OrdinaryDiffEqNordsieck/src/controllers.jl
+++ b/lib/OrdinaryDiffEqNordsieck/src/controllers.jl
@@ -1,6 +1,3 @@
-struct DummyController <: AbstractController
-end
-
 # JVODE
 function stepsize_controller!(integrator, alg::JVODE)
     if iszero(integrator.EEst)
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
index 01d8301d78..e989a5e8f9 100644
--- a/src/OrdinaryDiffEq.jl
+++ b/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, Relaxation
+                           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, Relaxation
+export IController, PIController, PIDController, NonAdaptiveController, Relaxation, RelaxationController
 end # module

From 49c80575975a60f6b32116f285ee4d5c3a415022 Mon Sep 17 00:00:00 2001
From: oscarddssmith <oscar.smith@juliacomputing.com>
Date: Tue, 3 Sep 2024 19:10:21 -0400
Subject: [PATCH 4/4] fix tests

---
 .../src/integrators/controllers.jl            |  2 +-
 test/relaxation/harmonic_oscillator.jl        | 18 ++++++-------
 test/relaxation/non_linear_oscillator.jl      | 18 +++++++------
 .../time_dependant_harmonic_oscillator.jl     | 25 ++++++++++---------
 4 files changed, 33 insertions(+), 30 deletions(-)

diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
index 8cd99c394f..5445bb5d7a 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
@@ -441,7 +441,7 @@ Controller to perform a relaxation on a step of a Runge-Kuttas method.
 mutable struct RelaxationController{CON, T} <: AbstractController
     controller::CON
     gamma::T
-    function RelaxationController(controller::AbstractController, T)
+    function RelaxationController(controller::AbstractController, T=Float64)
         new{typeof(controller), T}(controller, one(T))
     end
 end
diff --git a/test/relaxation/harmonic_oscillator.jl b/test/relaxation/harmonic_oscillator.jl
index 74f7c2f6ba..623e6fec9d 100644
--- a/test/relaxation/harmonic_oscillator.jl
+++ b/test/relaxation/harmonic_oscillator.jl
@@ -1,6 +1,6 @@
-using OrdinaryDiffEq, DiffEqDevTools, LinearAlgebra
+using OrdinaryDiffEq, DiffEqDevTools, LinearAlgebra, Test
 
-printstyled("Harmonic Oscillator\n"; bold = true)
+@testset "Harmonic oscillator" begin
 
 dts = (1 / 2) .^ (6:-1:4)
 
@@ -10,14 +10,14 @@ prob = ODEProblem(
     [1.0, 0.0],
     (0.0, 1.0))
 
-invariant(x) = norm(x)
-
 # Convergence with the old method Tsit5()
 sim = test_convergence(dts, prob, Tsit5(), adaptive = true)
-println("order of convergence of Tsit5 without relaxation : "*string(sim.𝒪est[:final]))
+@test sim.𝒪est[:final] ≈ 5 atol=0.2
 
-# Convergence with relaxation with FSAL-R, i.e  f(uᵧ,ₙ₊₁) ≈ f(uᵧ,ₙ) + γ ( f(uₙ₊₁) - f(uᵧ,ₙ)) 
-r = Relaxation(invariant)
-sim_relax = test_convergence(dts, prob, Tsit5(); relaxation = r, adaptive = true)
-println("order with relaxation with FSAL-R modification: "*string(sim_relax.𝒪est[:final]))
+# 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
diff --git a/test/relaxation/non_linear_oscillator.jl b/test/relaxation/non_linear_oscillator.jl
index 7a156e40b2..e7530f8856 100644
--- a/test/relaxation/non_linear_oscillator.jl
+++ b/test/relaxation/non_linear_oscillator.jl
@@ -1,6 +1,6 @@
-using OrdinaryDiffEq, DiffEqDevTools
+using OrdinaryDiffEq, DiffEqDevTools, LinearAlgebra, Test
 
-printstyled("Non linear harmonic oscillator\n"; bold = true)
+@testset "Non linear harmonic oscillator" begin
 
 dts = (1 / 2) .^ (6:-1:4)
 
@@ -10,13 +10,15 @@ prob = ODEProblem(
     [1.0, 0.0],
     (0.0, 1.0))
 
-invariant(x) = norm(x)
 
 # Convergence with the method Tsit5()
 sim = test_convergence(dts, prob, Tsit5())
-println("order of convergence of older perform_step! : "*string(sim.𝒪est[:final]))
+@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ᵧ,ₙ)) 
-r = Relaxation(invariant)
-sim = test_convergence(dts, prob, Tsit5(); relaxation = r)
-println("order with relaxation with FSAL-R modification: "*string(sim.𝒪est[:final]))
+# 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
diff --git a/test/relaxation/time_dependant_harmonic_oscillator.jl b/test/relaxation/time_dependant_harmonic_oscillator.jl
index 370c65bf83..4d25a07bdd 100644
--- a/test/relaxation/time_dependant_harmonic_oscillator.jl
+++ b/test/relaxation/time_dependant_harmonic_oscillator.jl
@@ -1,24 +1,25 @@
-using OrdinaryDiffEq, DiffEqDevTools
+using OrdinaryDiffEq, DiffEqDevTools, LinearAlgebra, Test
 
-printstyled("Time-dependent harmonic oscillator\n"; bold = true)
+
+@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)), 
+    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))
 
-invariant(x) = norm(x)
-
 # Convergence with the method Tsit5()
-sim = test_convergence(dts, prob, Tsit5(), adaptative = true)
-println("order of convergence of older perform_step! : "*string(sim.𝒪est[:final]))
+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ᵧ,ₙ)) 
-r = Relaxation(invariant)
-sim = test_convergence(dts, prob, Tsit5(); relaxation = r, adaptative = true)
-println("order with relaxation with FSAL-R modification: "*string(sim.𝒪est[:final]))
\ No newline at end of file
+# 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