From d84627c8c03d2b37134a4216b0d6232649678aa6 Mon Sep 17 00:00:00 2001
From: termi-official <9196588+termi-official@users.noreply.github.com>
Date: Mon, 6 Oct 2025 19:10:19 +0200
Subject: [PATCH 01/11] Prototype for using NonlinearSolve in IDSolve as the
 first step towards generic adaptivity

---
 lib/ImplicitDiscreteSolve/Project.toml        | 44 ++++++-------
 .../src/ImplicitDiscreteSolve.jl              | 14 ++--
 lib/ImplicitDiscreteSolve/src/algorithms.jl   | 21 ++++++
 lib/ImplicitDiscreteSolve/src/cache.jl        | 53 +++++++++++----
 lib/ImplicitDiscreteSolve/src/solve.jl        | 64 +++++++++----------
 lib/ImplicitDiscreteSolve/test/runtests.jl    | 36 +++++------
 6 files changed, 137 insertions(+), 95 deletions(-)
 create mode 100644 lib/ImplicitDiscreteSolve/src/algorithms.jl

diff --git a/lib/ImplicitDiscreteSolve/Project.toml b/lib/ImplicitDiscreteSolve/Project.toml
index d25e9f75d1..ef8116612f 100644
--- a/lib/ImplicitDiscreteSolve/Project.toml
+++ b/lib/ImplicitDiscreteSolve/Project.toml
@@ -4,38 +4,38 @@ authors = ["vyudu <vincent.duyuan@gmail.com>"]
 version = "1.2.0"
 
 [deps]
-SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
-SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
-[extras]
-JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
-OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
-AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a"
+[sources]
+OrdinaryDiffEqCore = {path = "../OrdinaryDiffEqCore"}
 
 [compat]
-Test = "1.10.0"
+AllocCheck = "0.2"
+Aqua = "0.8.11"
+DiffEqBase = "6.176"
+JET = "0.9.18, 0.10.4"
+NonlinearSolveFirstOrder = "1.9.0"
+OrdinaryDiffEqCore = "1.29.0"
 OrdinaryDiffEqSDIRK = "1.6.0"
+Reexport = "1.2"
 SciMLBase = "2.99"
-SimpleNonlinearSolve = "2.7"
-OrdinaryDiffEqCore = "1.29.0"
-Aqua = "0.8.11"
 SymbolicIndexingInterface = "0.3.38"
-julia = "1.10"
-JET = "0.9.18, 0.10.4"
+Test = "1.10.0"
 UnPack = "1.0.2"
-AllocCheck = "0.2"
-DiffEqBase = "6.176"
-Reexport = "1.2"
+julia = "1.10"
+
+[extras]
+AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a"
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["OrdinaryDiffEqSDIRK", "Test", "JET", "Aqua", "AllocCheck"]
-
-[sources.OrdinaryDiffEqCore]
-path = "../OrdinaryDiffEqCore"
diff --git a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
index b52829f5ed..0146afcee6 100644
--- a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
+++ b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
@@ -1,25 +1,21 @@
 module ImplicitDiscreteSolve
 
 using SciMLBase
-using SimpleNonlinearSolve
+using NonlinearSolveFirstOrder
 using UnPack
 using SymbolicIndexingInterface: parameter_symbols
 import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache,
                            OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal,
                            initialize!, perform_step!, isdiscretecache, isdiscretealg,
                            alg_order, beta2_default, beta1_default, dt_required,
-                           _initialize_dae!, DefaultInit, BrownFullBasicInit, OverrideInit
+                           _initialize_dae!, DefaultInit, BrownFullBasicInit, OverrideInit,
+                           OrdinaryDiffEqNewtonAdaptiveAlgorithm, @muladd, @..,
+                           AutoForwardDiff, _process_AD_choice, _unwrap_val
 
 using Reexport
 @reexport using SciMLBase
 
-"""
-    IDSolve()
-
-Simple solver for `ImplicitDiscreteSystems`. Uses `SimpleNewtonRaphson` to solve for the next state at every timestep.
-"""
-struct IDSolve <: OrdinaryDiffEqAlgorithm end
-
+include("algorithms.jl")
 include("cache.jl")
 include("solve.jl")
 include("alg_utils.jl")
diff --git a/lib/ImplicitDiscreteSolve/src/algorithms.jl b/lib/ImplicitDiscreteSolve/src/algorithms.jl
new file mode 100644
index 0000000000..d3523ae98f
--- /dev/null
+++ b/lib/ImplicitDiscreteSolve/src/algorithms.jl
@@ -0,0 +1,21 @@
+"""
+    IDSolve()
+
+First order solver for `ImplicitDiscreteSystems`.
+"""
+# struct IDSolve{CS, AD, NLS, FDT, ST, CJ} <:
+struct IDSolve{NLS} <:
+       OrdinaryDiffEqAlgorithm
+    nlsolve::NLS
+    extrapolant::Symbol
+    controller::Symbol
+end
+
+function IDSolve(; 
+        nlsolve = NewtonRaphson(), #NLNewton(),
+        extrapolant = :constant,
+        controller = :PI,
+    )
+
+    IDSolve{typeof(nlsolve)}(nlsolve, extrapolant, controller)
+end
diff --git a/lib/ImplicitDiscreteSolve/src/cache.jl b/lib/ImplicitDiscreteSolve/src/cache.jl
index f9c7835d82..0f7352a7d6 100644
--- a/lib/ImplicitDiscreteSolve/src/cache.jl
+++ b/lib/ImplicitDiscreteSolve/src/cache.jl
@@ -1,36 +1,65 @@
-mutable struct ImplicitDiscreteState{uType, pType, tType}
+struct ImplicitDiscreteState{uType, pType, tType}
     u::uType
     p::pType
-    t_next::tType
+    t::tType
 end
 
-mutable struct IDSolveCache{uType} <: OrdinaryDiffEqMutableCache
+mutable struct IDSolveCache{uType, cType} <: OrdinaryDiffEqMutableCache
     u::uType
     uprev::uType
-    state::ImplicitDiscreteState
-    prob::Union{Nothing, SciMLBase.AbstractNonlinearProblem}
+    z::uType
+    nlcache::cType
 end
 
 function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
         dt, reltol, p, calck,
         ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    state = ImplicitDiscreteState(isnothing(u) ? nothing : zero(u), p, t)
-    IDSolveCache(u, uprev, state, nothing)
+    state = ImplicitDiscreteState(zero(u), p, t)
+    f_nl = (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t)
+
+    u_len = length(u)
+    nlls = !isnothing(f.resid_prototype) && (length(f.resid_prototype) != u_len)
+    prob = if nlls
+        NonlinearLeastSquaresProblem{isinplace(f)}(
+            NonlinearFunction(f_nl; resid_prototype = f.resid_prototype),
+            u, state)
+    else
+        NonlinearProblem{isinplace(f)}(f_nl, u, state)
+    end
+
+    nlcache = init(prob, alg.nlsolve)
+
+    IDSolveCache(u, uprev, state.u, nlcache)
 end
 
 isdiscretecache(cache::IDSolveCache) = true
 
-struct IDSolveConstantCache <: OrdinaryDiffEqConstantCache
-    prob::Union{Nothing, SciMLBase.AbstractNonlinearProblem}
-end
+# struct IDSolveConstantCache <: OrdinaryDiffEqConstantCache
+#     prob::Union{Nothing, SciMLBase.AbstractNonlinearProblem}
+# end
 
 function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
         dt, reltol, p, calck,
         ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    state = ImplicitDiscreteState(isnothing(u) ? nothing : zero(u), p, t)
-    IDSolveCache(u, uprev, state, nothing)
+    state = ImplicitDiscreteState(zero(u), p, t)
+    f_nl = (u_next, p) -> f(u_next, p.u, p.p, p.t)
+
+    u_len = length(u)
+    nlls = !isnothing(f.resid_prototype) && (length(f.resid_prototype) != u_len)
+    prob = if nlls
+        NonlinearLeastSquaresProblem{isinplace(f)}(
+            NonlinearFunction(f_nl; resid_prototype = f.resid_prototype),
+            u, state)
+    else
+        NonlinearProblem{isinplace(f)}(f_nl, u, state)
+    end
+
+    nlcache = init(prob, alg.nlsolve)
+
+    # FIXME Use IDSolveConstantCache
+    IDSolveCache(u, uprev, state.u, nlcache)
 end
 
 get_fsalfirstlast(cache::IDSolveCache, rate_prototype) = (nothing, nothing)
diff --git a/lib/ImplicitDiscreteSolve/src/solve.jl b/lib/ImplicitDiscreteSolve/src/solve.jl
index e0580799d8..cd1dc38dd6 100644
--- a/lib/ImplicitDiscreteSolve/src/solve.jl
+++ b/lib/ImplicitDiscreteSolve/src/solve.jl
@@ -1,38 +1,33 @@
-# Remake the nonlinear problem, then update
 function perform_step!(integrator, cache::IDSolveCache, repeat_step = false)
-    @unpack alg, u, uprev, dt, t, f, p = integrator
-    @unpack state, prob = cache
-    state.u .= uprev
-    state.t_next = t
-    prob = remake(prob, p = state)
+    (; alg, u, uprev, dt, t, f, p) = integrator
 
-    u = solve(prob, SimpleNewtonRaphson())
-    integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, u.retcode)
-    integrator.u = u
-end
-
-function initialize!(integrator, cache::IDSolveCache)
-    integrator.u isa AbstractVector && (cache.state.u .= integrator.u)
-    cache.state.p = integrator.p
-    cache.state.t_next = integrator.t
-    f = integrator.f
-
-    _f = if isinplace(f)
-        (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t_next)
-    else
-        (u_next, p) -> f(u_next, p.u, p.p, p.t_next)
+    # initial guess
+    if alg.extrapolant == :linear
+        @.. broadcast=false cache.z=integrator.uprev + dt * (integrator.uprev - integrator.uprev2)
+    else # :constant
+        cache.z .= integrator.u
     end
-    u_len = isnothing(integrator.u) ? 0 : length(integrator.u)
-    nlls = !isnothing(f.resid_prototype) && (length(f.resid_prototype) != u_len)
+    state = ImplicitDiscreteState(cache.z, p, t)
 
-    prob = if nlls
-        NonlinearLeastSquaresProblem{isinplace(f)}(
-            NonlinearFunction(_f; resid_prototype = f.resid_prototype),
-            cache.state.u, cache.state)
-    else
-        NonlinearProblem{isinplace(f)}(_f, cache.state.u, cache.state)
+    # nonlinear solve step
+    SciMLBase.reinit!(cache.nlcache, p=state)
+    # TODO compute convergence rate estimate
+    # for i in 1:10
+    #     step!(cache.nlcache)
+    #     # ...
+    # end
+    solve!(cache.nlcache)
+    if cache.nlcache.retcode != ReturnCode.Success
+        integrator.force_stepfail = true
+        return
     end
-    cache.prob = prob
+
+    # Accept step
+    u .= cache.nlcache.u
+end
+
+function initialize!(integrator, cache::IDSolveCache)
+    integrator.u isa AbstractVector && (cache.z .= integrator.u)
 end
 
 function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem,
@@ -43,13 +38,14 @@ function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem,
         _initialize_dae!(integrator, prob,
             OverrideInit(atol), x)
     else
-        @unpack u, p, t, f = integrator
+        (; u, p, t, f) = integrator
+
         initstate = ImplicitDiscreteState(u, p, t)
 
         _f = if isinplace(f)
-            (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t_next)
+            (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t)
         else
-            (u_next, p) -> f(u_next, p.u, p.p, p.t_next)
+            (u_next, p) -> f(u_next, p.u, p.p, p.t)
         end
 
         nlls = !isnothing(f.resid_prototype) &&
@@ -60,7 +56,7 @@ function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem,
         else
             NonlinearProblem{isinplace(f)}(_f, u, initstate)
         end
-        sol = solve(prob, SimpleNewtonRaphson())
+        sol = solve(prob, integrator.alg.nlsolve)
         if sol.retcode == ReturnCode.Success
             integrator.u = sol
         else
diff --git a/lib/ImplicitDiscreteSolve/test/runtests.jl b/lib/ImplicitDiscreteSolve/test/runtests.jl
index 4c94621302..1410f2b446 100644
--- a/lib/ImplicitDiscreteSolve/test/runtests.jl
+++ b/lib/ImplicitDiscreteSolve/test/runtests.jl
@@ -4,7 +4,7 @@ using ImplicitDiscreteSolve
 using OrdinaryDiffEqCore
 using OrdinaryDiffEqSDIRK
 using SciMLBase
-using JET
+# using JET
 
 # Test implicit Euler using ImplicitDiscreteProblem
 @testset "Implicit Euler" begin
@@ -71,16 +71,16 @@ end
     end
 end
 
-@testset "Handle nothing in u0" begin
-    function empty(u_next, u, p, t)
-        nothing
-    end
+# @testset "Handle nothing in u0" begin
+#     function empty(u_next, u, p, t)
+#         nothing
+#     end
 
-    tsteps = 5
-    u0 = nothing
-    idprob = ImplicitDiscreteProblem(empty, u0, (0, tsteps), [])
-    @test_nowarn integ = init(idprob, IDSolve())
-end
+#     tsteps = 5
+#     u0 = nothing
+#     idprob = ImplicitDiscreteProblem(empty, u0, (0, tsteps), [])
+#     @test_nowarn integ = init(idprob, IDSolve())
+# end
 
 @testset "Create NonlinearLeastSquaresProblem" begin
     function over(u_next, u, p, t)
@@ -92,7 +92,7 @@ end
     idprob = ImplicitDiscreteProblem(
         ImplicitDiscreteFunction(over, resid_prototype = zeros(3)), u0, (0, tsteps), [])
     integ = init(idprob, IDSolve())
-    @test integ.cache.prob isa NonlinearLeastSquaresProblem
+    @test integ.cache.nlcache.prob isa NonlinearLeastSquaresProblem
 
     function under(u_next, u, p, t)
         [u_next[1] - u_next[2] - 1]
@@ -100,7 +100,7 @@ end
     idprob = ImplicitDiscreteProblem(
         ImplicitDiscreteFunction(under; resid_prototype = zeros(1)), u0, (0, tsteps), [])
     integ = init(idprob, IDSolve())
-    @test integ.cache.prob isa NonlinearLeastSquaresProblem
+    @test integ.cache.nlcache.prob isa NonlinearLeastSquaresProblem
 
     function full(u_next, u, p, t)
         [u_next[1]^2 - 3, u_next[2] - u[1]]
@@ -108,7 +108,7 @@ end
     idprob = ImplicitDiscreteProblem(
         ImplicitDiscreteFunction(full; resid_prototype = zeros(2)), u0, (0, tsteps), [])
     integ = init(idprob, IDSolve())
-    @test integ.cache.prob isa NonlinearProblem
+    @test integ.cache.nlcache.prob isa NonlinearProblem
 end
 
 @testset "InitialFailure thrown" begin
@@ -125,9 +125,9 @@ end
     @test !SciMLBase.successful_retcode(sol)
 end
 
-@testset "JET Tests" begin
-    test_package(
-        ImplicitDiscreteSolve, target_defined_modules = true, mode = :typo)
-end
+# @testset "JET Tests" begin
+#     test_package(
+#         ImplicitDiscreteSolve, target_defined_modules = true, mode = :typo)
+# end
 
-include("qa.jl")
+# include("qa.jl")

From f02492f26be861b2850c7da62073a76e566f9b1c Mon Sep 17 00:00:00 2001
From: termi-official <9196588+termi-official@users.noreply.github.com>
Date: Mon, 6 Oct 2025 19:15:27 +0200
Subject: [PATCH 02/11] Reenable JET tests

---
 lib/ImplicitDiscreteSolve/test/runtests.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/lib/ImplicitDiscreteSolve/test/runtests.jl b/lib/ImplicitDiscreteSolve/test/runtests.jl
index 1410f2b446..7b94d8e743 100644
--- a/lib/ImplicitDiscreteSolve/test/runtests.jl
+++ b/lib/ImplicitDiscreteSolve/test/runtests.jl
@@ -4,7 +4,7 @@ using ImplicitDiscreteSolve
 using OrdinaryDiffEqCore
 using OrdinaryDiffEqSDIRK
 using SciMLBase
-# using JET
+using JET
 
 # Test implicit Euler using ImplicitDiscreteProblem
 @testset "Implicit Euler" begin
@@ -125,9 +125,9 @@ end
     @test !SciMLBase.successful_retcode(sol)
 end
 
-# @testset "JET Tests" begin
-#     test_package(
-#         ImplicitDiscreteSolve, target_defined_modules = true, mode = :typo)
-# end
+@testset "JET Tests" begin
+    test_package(
+        ImplicitDiscreteSolve, target_defined_modules = true, mode = :typo)
+end
 
-# include("qa.jl")
+include("qa.jl")

From c96581c885f66f0fc49dbf4f32fde2e68fc77a43 Mon Sep 17 00:00:00 2001
From: termi-official <9196588+termi-official@users.noreply.github.com>
Date: Mon, 6 Oct 2025 21:50:42 +0200
Subject: [PATCH 03/11] Add poor mans adaptivity to get the ball rolling

---
 .../src/ImplicitDiscreteSolve.jl              |  4 +-
 lib/ImplicitDiscreteSolve/src/alg_utils.jl    |  4 +-
 lib/ImplicitDiscreteSolve/src/algorithms.jl   |  3 +-
 lib/ImplicitDiscreteSolve/src/controller.jl   | 80 +++++++++++++++++++
 lib/ImplicitDiscreteSolve/src/solve.jl        |  6 +-
 lib/ImplicitDiscreteSolve/test/runtests.jl    | 18 ++++-
 6 files changed, 104 insertions(+), 11 deletions(-)
 create mode 100644 lib/ImplicitDiscreteSolve/src/controller.jl

diff --git a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
index 0146afcee6..5f86796166 100644
--- a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
+++ b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
@@ -9,8 +9,7 @@ import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMut
                            initialize!, perform_step!, isdiscretecache, isdiscretealg,
                            alg_order, beta2_default, beta1_default, dt_required,
                            _initialize_dae!, DefaultInit, BrownFullBasicInit, OverrideInit,
-                           OrdinaryDiffEqNewtonAdaptiveAlgorithm, @muladd, @..,
-                           AutoForwardDiff, _process_AD_choice, _unwrap_val
+                           @muladd, @..,  _unwrap_val, OrdinaryDiffEqCore, isadaptive
 
 using Reexport
 @reexport using SciMLBase
@@ -19,6 +18,7 @@ include("algorithms.jl")
 include("cache.jl")
 include("solve.jl")
 include("alg_utils.jl")
+include("controller.jl")
 
 export IDSolve
 
diff --git a/lib/ImplicitDiscreteSolve/src/alg_utils.jl b/lib/ImplicitDiscreteSolve/src/alg_utils.jl
index 5c25c2c217..77470c5bc1 100644
--- a/lib/ImplicitDiscreteSolve/src/alg_utils.jl
+++ b/lib/ImplicitDiscreteSolve/src/alg_utils.jl
@@ -11,9 +11,11 @@ end
 SciMLBase.isdiscrete(alg::IDSolve) = true
 
 isfsal(alg::IDSolve) = false
-alg_order(alg::IDSolve) = 0
+alg_order(alg::IDSolve) = 1
 beta2_default(alg::IDSolve) = 0
 beta1_default(alg::IDSolve, beta2) = 0
 
 dt_required(alg::IDSolve) = false
 isdiscretealg(alg::IDSolve) = true
+
+isadaptive(alg::IDSolve) = true
diff --git a/lib/ImplicitDiscreteSolve/src/algorithms.jl b/lib/ImplicitDiscreteSolve/src/algorithms.jl
index d3523ae98f..70fc865cd4 100644
--- a/lib/ImplicitDiscreteSolve/src/algorithms.jl
+++ b/lib/ImplicitDiscreteSolve/src/algorithms.jl
@@ -3,7 +3,6 @@
 
 First order solver for `ImplicitDiscreteSystems`.
 """
-# struct IDSolve{CS, AD, NLS, FDT, ST, CJ} <:
 struct IDSolve{NLS} <:
        OrdinaryDiffEqAlgorithm
     nlsolve::NLS
@@ -12,7 +11,7 @@ struct IDSolve{NLS} <:
 end
 
 function IDSolve(; 
-        nlsolve = NewtonRaphson(), #NLNewton(),
+        nlsolve = NewtonRaphson(), 
         extrapolant = :constant,
         controller = :PI,
     )
diff --git a/lib/ImplicitDiscreteSolve/src/controller.jl b/lib/ImplicitDiscreteSolve/src/controller.jl
new file mode 100644
index 0000000000..816c013d2a
--- /dev/null
+++ b/lib/ImplicitDiscreteSolve/src/controller.jl
@@ -0,0 +1,80 @@
+# function should_accept_step(integrator::ThunderboltTimeIntegrator, cache::HomotopyPathSolverCache, controller::Deuflhard2004_B_DiscreteContinuationControllerVariant)
+#     (; Θks) = cache.inner_solver_cache
+#     (; Θreject) = controller
+#     if cache.inner_solver_cache.parameters.enforce_monotonic_convergence
+#         result = all(Θks .≤ Θreject)
+#         return result
+#     else
+#         return all(isfinite.(Θks))
+#     end
+# end
+# function reject_step!(integrator::ThunderboltTimeIntegrator, cache::HomotopyPathSolverCache, controller::Deuflhard2004_B_DiscreteContinuationControllerVariant)
+#     # Reset solution
+#     integrator.u .= integrator.uprev
+
+#     @inline g(x) = √(1+4x) - 1
+
+#     # Shorten dt according to (Eq. 5.24)
+#     (; Θks) = cache.inner_solver_cache
+#     (; Θbar, Θreject, γ, Θmin, qmin, qmax, p) = controller
+#     for Θk in Θks
+#         if Θk > Θreject
+#             q = clamp(γ * (g(Θbar)/g(Θk))^(1/p), qmin, qmax)
+#             integrator.dt = q * integrator.dt
+#             return
+#         end
+#     end
+# end
+
+# function adapt_dt!(integrator::ThunderboltTimeIntegrator, cache::HomotopyPathSolverCache, controller::Deuflhard2004_B_DiscreteContinuationControllerVariant)
+#     @inline g(x) = √(1+4x) - 1
+
+#     # Adapt dt with a priori estimate (Eq. 5.24)
+#     (; Θks) = cache.inner_solver_cache
+#     (; Θbar, γ, Θmin, qmin, qmax, p) = controller
+
+#     Θ₀ = length(Θks) > 0 ? max(first(Θks), Θmin) : Θmin
+#     q = clamp(γ * (g(Θbar)/(g(Θ₀)))^(1/p), qmin, qmax)
+#     integrator.dt = q * integrator.dt
+# end
+
+Base.@kwdef struct KantorovichTypeController <: OrdinaryDiffEqCore.AbstractController
+    Θmin::Float64
+    p::Int64
+    Θreject::Float64 = 0.95
+    Θbar::Float64 = 0.5
+    γ::Float64    = 0.95
+    qmin::Float64 = 1/5
+    qmax::Float64 = 5.0
+end
+
+function OrdinaryDiffEqCore.default_controller(alg::IDSolve, cache::IDSolveCache, _1, _2, _3)
+    return KantorovichTypeController(;Θmin=1//8, p=1)
+end
+
+function OrdinaryDiffEqCore.stepsize_controller!(
+    integrator, alg::IDSolve
+)
+    # @inline g(x) = √(1+4x) - 1
+
+    # # Adapt dt with a priori estimate (Eq. 5.24)
+    # (; Θks) = cache.inner_solver_cache
+    # (; Θbar, γ, Θmin, qmin, qmax, p) = controller
+
+    # Θ₀ = length(Θks) > 0 ? max(first(Θks), Θmin) : Θmin
+    # q = clamp(γ * (g(Θbar)/(g(Θ₀)))^(1/p), qmin, qmax)
+    return 1.1
+end
+
+function OrdinaryDiffEqCore.step_accept_controller!(
+    integrator, controller::KantorovichTypeController, alg::IDSolve, q
+)
+    integrator.dtpropose = q * integrator.dt
+end
+
+function OrdinaryDiffEqCore.step_reject_controller!(
+    integrator, controller::KantorovichTypeController, alg::IDSolve
+)
+    @info "Reject"
+    integrator.dt *= typeof(integrator.dt)(1//2)
+end
diff --git a/lib/ImplicitDiscreteSolve/src/solve.jl b/lib/ImplicitDiscreteSolve/src/solve.jl
index cd1dc38dd6..8ca3b157c6 100644
--- a/lib/ImplicitDiscreteSolve/src/solve.jl
+++ b/lib/ImplicitDiscreteSolve/src/solve.jl
@@ -1,13 +1,13 @@
 function perform_step!(integrator, cache::IDSolveCache, repeat_step = false)
-    (; alg, u, uprev, dt, t, f, p) = integrator
+    (; alg, u, uprev, dt, t, tprev, f, p) = integrator
 
     # initial guess
     if alg.extrapolant == :linear
-        @.. broadcast=false cache.z=integrator.uprev + dt * (integrator.uprev - integrator.uprev2)
+        @.. broadcast=false cache.z=integrator.u + dt * (integrator.u - integrator.uprev)
     else # :constant
         cache.z .= integrator.u
     end
-    state = ImplicitDiscreteState(cache.z, p, t)
+    state = ImplicitDiscreteState(cache.z, p, t+dt)
 
     # nonlinear solve step
     SciMLBase.reinit!(cache.nlcache, p=state)
diff --git a/lib/ImplicitDiscreteSolve/test/runtests.jl b/lib/ImplicitDiscreteSolve/test/runtests.jl
index 7b94d8e743..7ddc02808e 100644
--- a/lib/ImplicitDiscreteSolve/test/runtests.jl
+++ b/lib/ImplicitDiscreteSolve/test/runtests.jl
@@ -4,7 +4,7 @@ using ImplicitDiscreteSolve
 using OrdinaryDiffEqCore
 using OrdinaryDiffEqSDIRK
 using SciMLBase
-using JET
+# using JET
 
 # Test implicit Euler using ImplicitDiscreteProblem
 @testset "Implicit Euler" begin
@@ -22,7 +22,7 @@ using JET
     tspan = (0.0, 0.5)
 
     idprob = ImplicitDiscreteProblem(f!, u0, tspan, []; dt = 0.01)
-    idsol = solve(idprob, IDSolve())
+    idsol = solve(idprob, IDSolve(); adaptive=false)
 
     oprob = ODEProblem(lotkavolterra, u0, tspan)
     osol = solve(oprob, ImplicitEuler())
@@ -45,7 +45,7 @@ using JET
     tspan = (0, 0.2)
 
     idprob = ImplicitDiscreteProblem(g!, u0, tspan, []; dt = 0.01)
-    idsol = solve(idprob, IDSolve())
+    idsol = solve(idprob, IDSolve(); adaptive=false)
 
     oprob = ODEProblem(ff, u0, tspan)
     osol = solve(oprob, ImplicitEuler())
@@ -71,6 +71,18 @@ end
     end
 end
 
+@testset "Hard problem" begin
+    function hard!(resid, u, u_prev, p, t)
+        resid[1] = tanh((u[1]-10t)^2)/2
+    end
+
+    u0 = [0.0]
+    idprob = ImplicitDiscreteProblem(hard!, u0, (0.0, 1.0), [])
+    integrator = init(idprob, IDSolve())
+    idsol = solve!(integrator)
+    @test idsol.retcode == ReturnCode.Success
+end
+
 # @testset "Handle nothing in u0" begin
 #     function empty(u_next, u, p, t)
 #         nothing

From e54472a4527e8493613d7a10a8b02f1f945754e0 Mon Sep 17 00:00:00 2001
From: termi-official <9196588+termi-official@users.noreply.github.com>
Date: Mon, 6 Oct 2025 23:02:10 +0200
Subject: [PATCH 04/11] Add a modified variant of Deuflhard's controller using
 Kantorovich estimates.

---
 lib/ImplicitDiscreteSolve/Project.toml        |  6 +-
 .../src/ImplicitDiscreteSolve.jl              |  3 +-
 lib/ImplicitDiscreteSolve/src/alg_utils.jl    | 24 ++++++++
 lib/ImplicitDiscreteSolve/src/cache.jl        |  9 +--
 lib/ImplicitDiscreteSolve/src/controller.jl   | 15 ++---
 lib/ImplicitDiscreteSolve/src/solve.jl        | 56 ++++++++++++++++---
 6 files changed, 90 insertions(+), 23 deletions(-)

diff --git a/lib/ImplicitDiscreteSolve/Project.toml b/lib/ImplicitDiscreteSolve/Project.toml
index ef8116612f..f609fb5bac 100644
--- a/lib/ImplicitDiscreteSolve/Project.toml
+++ b/lib/ImplicitDiscreteSolve/Project.toml
@@ -4,13 +4,14 @@ authors = ["vyudu <vincent.duyuan@gmail.com>"]
 version = "1.2.0"
 
 [deps]
+ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [sources]
 OrdinaryDiffEqCore = {path = "../OrdinaryDiffEqCore"}
@@ -18,8 +19,10 @@ OrdinaryDiffEqCore = {path = "../OrdinaryDiffEqCore"}
 [compat]
 AllocCheck = "0.2"
 Aqua = "0.8.11"
+ConcreteStructs = "0.2.3"
 DiffEqBase = "6.176"
 JET = "0.9.18, 0.10.4"
+NonlinearSolveBase = "2.0.0"
 NonlinearSolveFirstOrder = "1.9.0"
 OrdinaryDiffEqCore = "1.29.0"
 OrdinaryDiffEqSDIRK = "1.6.0"
@@ -27,7 +30,6 @@ Reexport = "1.2"
 SciMLBase = "2.99"
 SymbolicIndexingInterface = "0.3.38"
 Test = "1.10.0"
-UnPack = "1.0.2"
 julia = "1.10"
 
 [extras]
diff --git a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
index 5f86796166..ab809cad4b 100644
--- a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
+++ b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
@@ -1,8 +1,9 @@
 module ImplicitDiscreteSolve
 
 using SciMLBase
+using NonlinearSolveBase
 using NonlinearSolveFirstOrder
-using UnPack
+using ConcreteStructs
 using SymbolicIndexingInterface: parameter_symbols
 import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache,
                            OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal,
diff --git a/lib/ImplicitDiscreteSolve/src/alg_utils.jl b/lib/ImplicitDiscreteSolve/src/alg_utils.jl
index 77470c5bc1..55d23c555f 100644
--- a/lib/ImplicitDiscreteSolve/src/alg_utils.jl
+++ b/lib/ImplicitDiscreteSolve/src/alg_utils.jl
@@ -19,3 +19,27 @@ dt_required(alg::IDSolve) = false
 isdiscretealg(alg::IDSolve) = true
 
 isadaptive(alg::IDSolve) = true
+
+# @concrete struct ConvergenceRateTracing
+#     inner_tracing
+# end
+
+# @concrete struct ConvergenceRateTraceTrick
+#     incrementL2norms
+#     residualL2norms
+#     trace_wrapper
+# end
+
+# function NonlinearSolveBase.init_nonlinearsolve_trace(
+#         prob, alg::IDSolve, u, fu, J, δu;
+#         trace_level::ConvergenceRateTracing, kwargs... # This kind of dispatch does not work. Need to figure out a different way.
+# )
+#     inner_trace = NonlinearSolveBase.init_nonlinearsolve_trace(
+#         prob, alg, u, fu, J, δu;
+#         trace_level.inner_tracing, kwargs...
+#     )
+
+#     return ConvergenceRateTraceTrick(eltype(δu)[], eltype(fu)[], inner_trace)
+# end
+
+
diff --git a/lib/ImplicitDiscreteSolve/src/cache.jl b/lib/ImplicitDiscreteSolve/src/cache.jl
index 0f7352a7d6..f3ef6ac751 100644
--- a/lib/ImplicitDiscreteSolve/src/cache.jl
+++ b/lib/ImplicitDiscreteSolve/src/cache.jl
@@ -4,11 +4,12 @@ struct ImplicitDiscreteState{uType, pType, tType}
     t::tType
 end
 
-mutable struct IDSolveCache{uType, cType} <: OrdinaryDiffEqMutableCache
+mutable struct IDSolveCache{uType, cType, thetaType} <: OrdinaryDiffEqMutableCache
     u::uType
     uprev::uType
     z::uType
     nlcache::cType
+    Θks::thetaType
 end
 
 function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -30,7 +31,7 @@ function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits},
 
     nlcache = init(prob, alg.nlsolve)
 
-    IDSolveCache(u, uprev, state.u, nlcache)
+    IDSolveCache(u, uprev, state.u, nlcache, uBottomEltypeNoUnits[])
 end
 
 isdiscretecache(cache::IDSolveCache) = true
@@ -58,8 +59,8 @@ function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits},
 
     nlcache = init(prob, alg.nlsolve)
 
-    # FIXME Use IDSolveConstantCache
-    IDSolveCache(u, uprev, state.u, nlcache)
+    # FIXME Use IDSolveConstantCache?
+    IDSolveCache(u, uprev, state.u, nlcache, uBottomEltypeNoUnits[])
 end
 
 get_fsalfirstlast(cache::IDSolveCache, rate_prototype) = (nothing, nothing)
diff --git a/lib/ImplicitDiscreteSolve/src/controller.jl b/lib/ImplicitDiscreteSolve/src/controller.jl
index 816c013d2a..c55afa22a1 100644
--- a/lib/ImplicitDiscreteSolve/src/controller.jl
+++ b/lib/ImplicitDiscreteSolve/src/controller.jl
@@ -55,15 +55,16 @@ end
 function OrdinaryDiffEqCore.stepsize_controller!(
     integrator, alg::IDSolve
 )
-    # @inline g(x) = √(1+4x) - 1
+    @inline g(x) = √(1+4x) - 1
 
-    # # Adapt dt with a priori estimate (Eq. 5.24)
-    # (; Θks) = cache.inner_solver_cache
-    # (; Θbar, γ, Θmin, qmin, qmax, p) = controller
+    # Adapt dt with a priori estimate (Eq. 5.24)
+    (; Θks) = integrator.cache
+    (; Θbar, γ, Θmin, qmin, qmax, p) = integrator.opts.controller
 
-    # Θ₀ = length(Θks) > 0 ? max(first(Θks), Θmin) : Θmin
-    # q = clamp(γ * (g(Θbar)/(g(Θ₀)))^(1/p), qmin, qmax)
-    return 1.1
+    Θ₀ = length(Θks) > 0 ? max(first(Θks), Θmin) : Θmin
+    q = clamp(γ * (g(Θbar)/(g(Θ₀)))^(1/p), qmin, qmax)
+
+    return q
 end
 
 function OrdinaryDiffEqCore.step_accept_controller!(
diff --git a/lib/ImplicitDiscreteSolve/src/solve.jl b/lib/ImplicitDiscreteSolve/src/solve.jl
index 8ca3b157c6..830e93571f 100644
--- a/lib/ImplicitDiscreteSolve/src/solve.jl
+++ b/lib/ImplicitDiscreteSolve/src/solve.jl
@@ -1,5 +1,6 @@
 function perform_step!(integrator, cache::IDSolveCache, repeat_step = false)
     (; alg, u, uprev, dt, t, tprev, f, p) = integrator
+    (; nlcache, Θks) = cache
 
     # initial guess
     if alg.extrapolant == :linear
@@ -10,20 +11,57 @@ function perform_step!(integrator, cache::IDSolveCache, repeat_step = false)
     state = ImplicitDiscreteState(cache.z, p, t+dt)
 
     # nonlinear solve step
-    SciMLBase.reinit!(cache.nlcache, p=state)
-    # TODO compute convergence rate estimate
-    # for i in 1:10
-    #     step!(cache.nlcache)
-    #     # ...
-    # end
-    solve!(cache.nlcache)
-    if cache.nlcache.retcode != ReturnCode.Success
+    SciMLBase.reinit!(nlcache, p=state)
+
+    # solve!(nlcache)
+    # The solve here is simply unrolled by hand to quiery the convergence rate estimates "manually" for now
+    if nlcache.retcode == ReturnCode.InitialFailure
+        integrator.force_stepfail = true
+        return
+    end
+
+    resize!(Θks, 0)
+    residualnormprev = zero(eltype(u))
+    while NonlinearSolveBase.not_terminated(nlcache)
+        step!(nlcache)
+        residualnorm = NonlinearSolveBase.L2_NORM(nlcache.fu)
+        if nlcache.nsteps > 1
+            # Θk = min(residualnorm/residualnormprev, incrementnorm/incrementnormprev)
+            Θk = residualnorm/residualnormprev
+            if residualnormprev ≈ 0.0 #|| incrementnormprev ≈ 0.0
+                push!(Θks, 0.0)
+            else
+                push!(Θks, Θk)
+            end
+            # if nlcache.parameters.enforce_monotonic_convergence && Θk ≥ 1.0
+            #     @debug "Newton-Raphson diverged. Aborting. ||r|| = $residualnorm" _group=:nlsolve
+            #     return false
+            # end
+        end
+        residualnormprev = residualnorm
+    end
+
+    # The solver might have set a different `retcode`
+    if nlcache.retcode == ReturnCode.Default
+        nlcache.retcode = ifelse(
+            nlcache.nsteps ≥ nlcache.maxiters, ReturnCode.MaxIters, ReturnCode.Success
+        )
+    end
+
+    NonlinearSolveBase.update_from_termination_cache!(nlcache.termination_cache, nlcache)
+
+    NonlinearSolveBase.update_trace!(
+        nlcache.trace, nlcache.nsteps, NonlinearSolveBase.get_u(nlcache), NonlinearSolveBase.get_fu(nlcache), nothing, nothing, nothing;
+        last = Val(true)
+    )
+
+    if nlcache.retcode != ReturnCode.Success
         integrator.force_stepfail = true
         return
     end
 
     # Accept step
-    u .= cache.nlcache.u
+    u .= nlcache.u
 end
 
 function initialize!(integrator, cache::IDSolveCache)

From cd0fbf679be7e67f4b59bd4ef983d3f31300b4e6 Mon Sep 17 00:00:00 2001
From: termi-official <9196588+termi-official@users.noreply.github.com>
Date: Mon, 6 Oct 2025 23:20:25 +0200
Subject: [PATCH 05/11] Cleanup debris and rejection implementation

---
 lib/ImplicitDiscreteSolve/src/controller.jl | 61 ++++++---------------
 1 file changed, 16 insertions(+), 45 deletions(-)

diff --git a/lib/ImplicitDiscreteSolve/src/controller.jl b/lib/ImplicitDiscreteSolve/src/controller.jl
index c55afa22a1..c6cfe6d84d 100644
--- a/lib/ImplicitDiscreteSolve/src/controller.jl
+++ b/lib/ImplicitDiscreteSolve/src/controller.jl
@@ -1,43 +1,3 @@
-# function should_accept_step(integrator::ThunderboltTimeIntegrator, cache::HomotopyPathSolverCache, controller::Deuflhard2004_B_DiscreteContinuationControllerVariant)
-#     (; Θks) = cache.inner_solver_cache
-#     (; Θreject) = controller
-#     if cache.inner_solver_cache.parameters.enforce_monotonic_convergence
-#         result = all(Θks .≤ Θreject)
-#         return result
-#     else
-#         return all(isfinite.(Θks))
-#     end
-# end
-# function reject_step!(integrator::ThunderboltTimeIntegrator, cache::HomotopyPathSolverCache, controller::Deuflhard2004_B_DiscreteContinuationControllerVariant)
-#     # Reset solution
-#     integrator.u .= integrator.uprev
-
-#     @inline g(x) = √(1+4x) - 1
-
-#     # Shorten dt according to (Eq. 5.24)
-#     (; Θks) = cache.inner_solver_cache
-#     (; Θbar, Θreject, γ, Θmin, qmin, qmax, p) = controller
-#     for Θk in Θks
-#         if Θk > Θreject
-#             q = clamp(γ * (g(Θbar)/g(Θk))^(1/p), qmin, qmax)
-#             integrator.dt = q * integrator.dt
-#             return
-#         end
-#     end
-# end
-
-# function adapt_dt!(integrator::ThunderboltTimeIntegrator, cache::HomotopyPathSolverCache, controller::Deuflhard2004_B_DiscreteContinuationControllerVariant)
-#     @inline g(x) = √(1+4x) - 1
-
-#     # Adapt dt with a priori estimate (Eq. 5.24)
-#     (; Θks) = cache.inner_solver_cache
-#     (; Θbar, γ, Θmin, qmin, qmax, p) = controller
-
-#     Θ₀ = length(Θks) > 0 ? max(first(Θks), Θmin) : Θmin
-#     q = clamp(γ * (g(Θbar)/(g(Θ₀)))^(1/p), qmin, qmax)
-#     integrator.dt = q * integrator.dt
-# end
-
 Base.@kwdef struct KantorovichTypeController <: OrdinaryDiffEqCore.AbstractController
     Θmin::Float64
     p::Int64
@@ -53,13 +13,13 @@ function OrdinaryDiffEqCore.default_controller(alg::IDSolve, cache::IDSolveCache
 end
 
 function OrdinaryDiffEqCore.stepsize_controller!(
-    integrator, alg::IDSolve
+    integrator, controller::KantorovichTypeController, alg::IDSolve
 )
     @inline g(x) = √(1+4x) - 1
 
     # Adapt dt with a priori estimate (Eq. 5.24)
     (; Θks) = integrator.cache
-    (; Θbar, γ, Θmin, qmin, qmax, p) = integrator.opts.controller
+    (; Θbar, γ, Θmin, qmin, qmax, p) = controller
 
     Θ₀ = length(Θks) > 0 ? max(first(Θks), Θmin) : Θmin
     q = clamp(γ * (g(Θbar)/(g(Θ₀)))^(1/p), qmin, qmax)
@@ -70,12 +30,23 @@ end
 function OrdinaryDiffEqCore.step_accept_controller!(
     integrator, controller::KantorovichTypeController, alg::IDSolve, q
 )
-    integrator.dtpropose = q * integrator.dt
+    @info integrator.dt, q
+    return q * integrator.dt
 end
 
 function OrdinaryDiffEqCore.step_reject_controller!(
     integrator, controller::KantorovichTypeController, alg::IDSolve
 )
-    @info "Reject"
-    integrator.dt *= typeof(integrator.dt)(1//2)
+    @inline g(x) = √(1+4x) - 1
+
+    # Shorten dt according to (Eq. 5.24)
+    (; Θks) = cache.inner_solver_cache
+    (; Θbar, Θreject, γ, Θmin, qmin, qmax, p) = controller
+    for Θk in Θks
+        if Θk > Θreject
+            q = clamp(γ * (g(Θbar)/g(Θk))^(1/p), qmin, qmax)
+            integrator.dt = q * integrator.dt
+            return
+        end
+    end
 end

From c25238c7be608d6c78c9205bf4e047be818b99fa Mon Sep 17 00:00:00 2001
From: termi-official <9196588+termi-official@users.noreply.github.com>
Date: Mon, 6 Oct 2025 23:30:38 +0200
Subject: [PATCH 06/11] :)

---
 lib/ImplicitDiscreteSolve/test/runtests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/ImplicitDiscreteSolve/test/runtests.jl b/lib/ImplicitDiscreteSolve/test/runtests.jl
index 7ddc02808e..6207fc004a 100644
--- a/lib/ImplicitDiscreteSolve/test/runtests.jl
+++ b/lib/ImplicitDiscreteSolve/test/runtests.jl
@@ -4,7 +4,7 @@ using ImplicitDiscreteSolve
 using OrdinaryDiffEqCore
 using OrdinaryDiffEqSDIRK
 using SciMLBase
-# using JET
+using JET
 
 # Test implicit Euler using ImplicitDiscreteProblem
 @testset "Implicit Euler" begin

From e9a58fad76912e07158518bec788893279018b67 Mon Sep 17 00:00:00 2001
From: termi-official <9196588+termi-official@users.noreply.github.com>
Date: Tue, 7 Oct 2025 16:29:49 +0200
Subject: [PATCH 07/11] Add partial handling of empty u

---
 lib/ImplicitDiscreteSolve/src/cache.jl      | 24 +++++++++----------
 lib/ImplicitDiscreteSolve/src/controller.jl |  2 +-
 lib/ImplicitDiscreteSolve/src/solve.jl      |  2 +-
 lib/ImplicitDiscreteSolve/test/runtests.jl  | 26 +++++++++++++--------
 4 files changed, 30 insertions(+), 24 deletions(-)

diff --git a/lib/ImplicitDiscreteSolve/src/cache.jl b/lib/ImplicitDiscreteSolve/src/cache.jl
index f3ef6ac751..f3200762d4 100644
--- a/lib/ImplicitDiscreteSolve/src/cache.jl
+++ b/lib/ImplicitDiscreteSolve/src/cache.jl
@@ -16,17 +16,18 @@ function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
         dt, reltol, p, calck,
         ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    state = ImplicitDiscreteState(zero(u), p, t)
+    state = ImplicitDiscreteState(isnothing(u) ? nothing : zero(u), p, t)
     f_nl = (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t)
 
-    u_len = length(u)
+    u_len = isnothing(u) ? 0 : length(u)
     nlls = !isnothing(f.resid_prototype) && (length(f.resid_prototype) != u_len)
+    unl = isnothing(u) ? Float64[] : u # FIXME nonlinear solve cannot handle nothing for u
     prob = if nlls
         NonlinearLeastSquaresProblem{isinplace(f)}(
             NonlinearFunction(f_nl; resid_prototype = f.resid_prototype),
-            u, state)
+            unl, state)
     else
-        NonlinearProblem{isinplace(f)}(f_nl, u, state)
+        NonlinearProblem{isinplace(f)}(f_nl, unl, state)
     end
 
     nlcache = init(prob, alg.nlsolve)
@@ -36,25 +37,24 @@ end
 
 isdiscretecache(cache::IDSolveCache) = true
 
-# struct IDSolveConstantCache <: OrdinaryDiffEqConstantCache
-#     prob::Union{Nothing, SciMLBase.AbstractNonlinearProblem}
-# end
-
 function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
         dt, reltol, p, calck,
         ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    state = ImplicitDiscreteState(zero(u), p, t)
+    @assert !isnothing(u) "Empty u not supported with out of place functions yet."
+
+    state = ImplicitDiscreteState(isnothing(u) ? nothing : zero(u), p, t)
     f_nl = (u_next, p) -> f(u_next, p.u, p.p, p.t)
 
-    u_len = length(u)
+    u_len = isnothing(u) ? 0 : length(u)
     nlls = !isnothing(f.resid_prototype) && (length(f.resid_prototype) != u_len)
+    unl = isnothing(u) ? Float64[] : u # FIXME nonlinear solve cannot handle nothing for u
     prob = if nlls
         NonlinearLeastSquaresProblem{isinplace(f)}(
             NonlinearFunction(f_nl; resid_prototype = f.resid_prototype),
-            u, state)
+            unl, state)
     else
-        NonlinearProblem{isinplace(f)}(f_nl, u, state)
+        NonlinearProblem{isinplace(f)}(f_nl, unl, state)
     end
 
     nlcache = init(prob, alg.nlsolve)
diff --git a/lib/ImplicitDiscreteSolve/src/controller.jl b/lib/ImplicitDiscreteSolve/src/controller.jl
index c6cfe6d84d..cf4bf722a2 100644
--- a/lib/ImplicitDiscreteSolve/src/controller.jl
+++ b/lib/ImplicitDiscreteSolve/src/controller.jl
@@ -40,7 +40,7 @@ function OrdinaryDiffEqCore.step_reject_controller!(
     @inline g(x) = √(1+4x) - 1
 
     # Shorten dt according to (Eq. 5.24)
-    (; Θks) = cache.inner_solver_cache
+    (; Θks) = integrator.cache.nlcache
     (; Θbar, Θreject, γ, Θmin, qmin, qmax, p) = controller
     for Θk in Θks
         if Θk > Θreject
diff --git a/lib/ImplicitDiscreteSolve/src/solve.jl b/lib/ImplicitDiscreteSolve/src/solve.jl
index 830e93571f..ffb3520010 100644
--- a/lib/ImplicitDiscreteSolve/src/solve.jl
+++ b/lib/ImplicitDiscreteSolve/src/solve.jl
@@ -14,7 +14,7 @@ function perform_step!(integrator, cache::IDSolveCache, repeat_step = false)
     SciMLBase.reinit!(nlcache, p=state)
 
     # solve!(nlcache)
-    # The solve here is simply unrolled by hand to quiery the convergence rate estimates "manually" for now
+    # The solve here is simply unrolled by hand to query the convergence rate estimates "manually" for now
     if nlcache.retcode == ReturnCode.InitialFailure
         integrator.force_stepfail = true
         return
diff --git a/lib/ImplicitDiscreteSolve/test/runtests.jl b/lib/ImplicitDiscreteSolve/test/runtests.jl
index 6207fc004a..e6f3849369 100644
--- a/lib/ImplicitDiscreteSolve/test/runtests.jl
+++ b/lib/ImplicitDiscreteSolve/test/runtests.jl
@@ -83,16 +83,22 @@ end
     @test idsol.retcode == ReturnCode.Success
 end
 
-# @testset "Handle nothing in u0" begin
-#     function empty(u_next, u, p, t)
-#         nothing
-#     end
-
-#     tsteps = 5
-#     u0 = nothing
-#     idprob = ImplicitDiscreteProblem(empty, u0, (0, tsteps), [])
-#     @test_nowarn integ = init(idprob, IDSolve())
-# end
+@testset "Handle nothing in u0" begin
+    function emptyiip(residual, u_next, u, p, t) # TODO OOP variant does not work yet
+        nothing
+    end
+    function emptyoop(u_next, u, p, t) # TODO OOP variant does not work yet
+        nothing
+    end
+
+    tsteps = 5
+    u0 = nothing
+    idprob = ImplicitDiscreteProblem(emptyiip, u0, (0, tsteps), [])
+    @test_nowarn integ = init(idprob, IDSolve())
+
+    idprob2 = ImplicitDiscreteProblem(emptyoop, u0, (0, tsteps), [])
+    @test_throws AssertionError("Empty u not supported with out of place functions yet.") integ = init(idprob2, IDSolve())
+end
 
 @testset "Create NonlinearLeastSquaresProblem" begin
     function over(u_next, u, p, t)

From bdfad629d940b1d963490c81b5e413e34c710092 Mon Sep 17 00:00:00 2001
From: termi-official <9196588+termi-official@users.noreply.github.com>
Date: Tue, 7 Oct 2025 16:30:30 +0200
Subject: [PATCH 08/11] Remove debug message

---
 lib/ImplicitDiscreteSolve/src/controller.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/lib/ImplicitDiscreteSolve/src/controller.jl b/lib/ImplicitDiscreteSolve/src/controller.jl
index cf4bf722a2..724772d3c2 100644
--- a/lib/ImplicitDiscreteSolve/src/controller.jl
+++ b/lib/ImplicitDiscreteSolve/src/controller.jl
@@ -30,7 +30,6 @@ end
 function OrdinaryDiffEqCore.step_accept_controller!(
     integrator, controller::KantorovichTypeController, alg::IDSolve, q
 )
-    @info integrator.dt, q
     return q * integrator.dt
 end
 

From 4b597fb82a9ea3469049ec41c348d83581c36640 Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 14 Oct 2025 17:50:51 +0200
Subject: [PATCH 09/11] Format

---
 .../src/ImplicitDiscreteSolve.jl              |  2 +-
 lib/ImplicitDiscreteSolve/src/alg_utils.jl    |  2 --
 lib/ImplicitDiscreteSolve/src/algorithms.jl   |  9 ++++----
 lib/ImplicitDiscreteSolve/src/controller.jl   | 23 ++++++++++---------
 lib/ImplicitDiscreteSolve/src/solve.jl        | 12 ++++++----
 lib/ImplicitDiscreteSolve/test/runtests.jl    | 21 +++++++++--------
 6 files changed, 35 insertions(+), 34 deletions(-)

diff --git a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
index ab809cad4b..60d05501aa 100644
--- a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
+++ b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
@@ -10,7 +10,7 @@ import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMut
                            initialize!, perform_step!, isdiscretecache, isdiscretealg,
                            alg_order, beta2_default, beta1_default, dt_required,
                            _initialize_dae!, DefaultInit, BrownFullBasicInit, OverrideInit,
-                           @muladd, @..,  _unwrap_val, OrdinaryDiffEqCore, isadaptive
+                           @muladd, @.., _unwrap_val, OrdinaryDiffEqCore, isadaptive
 
 using Reexport
 @reexport using SciMLBase
diff --git a/lib/ImplicitDiscreteSolve/src/alg_utils.jl b/lib/ImplicitDiscreteSolve/src/alg_utils.jl
index 55d23c555f..a2d6e6ecb5 100644
--- a/lib/ImplicitDiscreteSolve/src/alg_utils.jl
+++ b/lib/ImplicitDiscreteSolve/src/alg_utils.jl
@@ -41,5 +41,3 @@ isadaptive(alg::IDSolve) = true
 
 #     return ConvergenceRateTraceTrick(eltype(δu)[], eltype(fu)[], inner_trace)
 # end
-
-
diff --git a/lib/ImplicitDiscreteSolve/src/algorithms.jl b/lib/ImplicitDiscreteSolve/src/algorithms.jl
index 70fc865cd4..aa2c567598 100644
--- a/lib/ImplicitDiscreteSolve/src/algorithms.jl
+++ b/lib/ImplicitDiscreteSolve/src/algorithms.jl
@@ -10,11 +10,10 @@ struct IDSolve{NLS} <:
     controller::Symbol
 end
 
-function IDSolve(; 
-        nlsolve = NewtonRaphson(), 
+function IDSolve(;
+        nlsolve = NewtonRaphson(),
         extrapolant = :constant,
-        controller = :PI,
-    )
-
+        controller = :PI
+)
     IDSolve{typeof(nlsolve)}(nlsolve, extrapolant, controller)
 end
diff --git a/lib/ImplicitDiscreteSolve/src/controller.jl b/lib/ImplicitDiscreteSolve/src/controller.jl
index 724772d3c2..63a27c9f84 100644
--- a/lib/ImplicitDiscreteSolve/src/controller.jl
+++ b/lib/ImplicitDiscreteSolve/src/controller.jl
@@ -3,47 +3,48 @@ Base.@kwdef struct KantorovichTypeController <: OrdinaryDiffEqCore.AbstractContr
     p::Int64
     Θreject::Float64 = 0.95
     Θbar::Float64 = 0.5
-    γ::Float64    = 0.95
-    qmin::Float64 = 1/5
+    γ::Float64 = 0.95
+    qmin::Float64 = 1 / 5
     qmax::Float64 = 5.0
 end
 
-function OrdinaryDiffEqCore.default_controller(alg::IDSolve, cache::IDSolveCache, _1, _2, _3)
-    return KantorovichTypeController(;Θmin=1//8, p=1)
+function OrdinaryDiffEqCore.default_controller(
+        alg::IDSolve, cache::IDSolveCache, _1, _2, _3)
+    return KantorovichTypeController(; Θmin = 1 // 8, p = 1)
 end
 
 function OrdinaryDiffEqCore.stepsize_controller!(
-    integrator, controller::KantorovichTypeController, alg::IDSolve
+        integrator, controller::KantorovichTypeController, alg::IDSolve
 )
-    @inline g(x) = √(1+4x) - 1
+    @inline g(x) = √(1 + 4x) - 1
 
     # Adapt dt with a priori estimate (Eq. 5.24)
     (; Θks) = integrator.cache
     (; Θbar, γ, Θmin, qmin, qmax, p) = controller
 
     Θ₀ = length(Θks) > 0 ? max(first(Θks), Θmin) : Θmin
-    q = clamp(γ * (g(Θbar)/(g(Θ₀)))^(1/p), qmin, qmax)
+    q = clamp(γ * (g(Θbar) / (g(Θ₀)))^(1 / p), qmin, qmax)
 
     return q
 end
 
 function OrdinaryDiffEqCore.step_accept_controller!(
-    integrator, controller::KantorovichTypeController, alg::IDSolve, q
+        integrator, controller::KantorovichTypeController, alg::IDSolve, q
 )
     return q * integrator.dt
 end
 
 function OrdinaryDiffEqCore.step_reject_controller!(
-    integrator, controller::KantorovichTypeController, alg::IDSolve
+        integrator, controller::KantorovichTypeController, alg::IDSolve
 )
-    @inline g(x) = √(1+4x) - 1
+    @inline g(x) = √(1 + 4x) - 1
 
     # Shorten dt according to (Eq. 5.24)
     (; Θks) = integrator.cache.nlcache
     (; Θbar, Θreject, γ, Θmin, qmin, qmax, p) = controller
     for Θk in Θks
         if Θk > Θreject
-            q = clamp(γ * (g(Θbar)/g(Θk))^(1/p), qmin, qmax)
+            q = clamp(γ * (g(Θbar) / g(Θk))^(1 / p), qmin, qmax)
             integrator.dt = q * integrator.dt
             return
         end
diff --git a/lib/ImplicitDiscreteSolve/src/solve.jl b/lib/ImplicitDiscreteSolve/src/solve.jl
index ffb3520010..d5a511fbab 100644
--- a/lib/ImplicitDiscreteSolve/src/solve.jl
+++ b/lib/ImplicitDiscreteSolve/src/solve.jl
@@ -8,10 +8,10 @@ function perform_step!(integrator, cache::IDSolveCache, repeat_step = false)
     else # :constant
         cache.z .= integrator.u
     end
-    state = ImplicitDiscreteState(cache.z, p, t+dt)
+    state = ImplicitDiscreteState(cache.z, p, t + dt)
 
     # nonlinear solve step
-    SciMLBase.reinit!(nlcache, p=state)
+    SciMLBase.reinit!(nlcache, p = state)
 
     # solve!(nlcache)
     # The solve here is simply unrolled by hand to query the convergence rate estimates "manually" for now
@@ -27,7 +27,7 @@ function perform_step!(integrator, cache::IDSolveCache, repeat_step = false)
         residualnorm = NonlinearSolveBase.L2_NORM(nlcache.fu)
         if nlcache.nsteps > 1
             # Θk = min(residualnorm/residualnormprev, incrementnorm/incrementnormprev)
-            Θk = residualnorm/residualnormprev
+            Θk = residualnorm / residualnormprev
             if residualnormprev ≈ 0.0 #|| incrementnormprev ≈ 0.0
                 push!(Θks, 0.0)
             else
@@ -51,7 +51,8 @@ function perform_step!(integrator, cache::IDSolveCache, repeat_step = false)
     NonlinearSolveBase.update_from_termination_cache!(nlcache.termination_cache, nlcache)
 
     NonlinearSolveBase.update_trace!(
-        nlcache.trace, nlcache.nsteps, NonlinearSolveBase.get_u(nlcache), NonlinearSolveBase.get_fu(nlcache), nothing, nothing, nothing;
+        nlcache.trace, nlcache.nsteps, NonlinearSolveBase.get_u(nlcache),
+        NonlinearSolveBase.get_fu(nlcache), nothing, nothing, nothing;
         last = Val(true)
     )
 
@@ -98,7 +99,8 @@ function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem,
         if sol.retcode == ReturnCode.Success
             integrator.u = sol
         else
-            integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, ReturnCode.InitialFailure)
+            integrator.sol = SciMLBase.solution_new_retcode(
+                integrator.sol, ReturnCode.InitialFailure)
         end
     end
 end
diff --git a/lib/ImplicitDiscreteSolve/test/runtests.jl b/lib/ImplicitDiscreteSolve/test/runtests.jl
index e6f3849369..7837874636 100644
--- a/lib/ImplicitDiscreteSolve/test/runtests.jl
+++ b/lib/ImplicitDiscreteSolve/test/runtests.jl
@@ -9,20 +9,20 @@ using JET
 # Test implicit Euler using ImplicitDiscreteProblem
 @testset "Implicit Euler" begin
     function lotkavolterra(u, p, t)
-        [1.5*u[1] - u[1]*u[2], -3.0*u[2] + u[1]*u[2]]
+        [1.5 * u[1] - u[1] * u[2], -3.0 * u[2] + u[1] * u[2]]
     end
 
     function f!(resid, u_next, u, p, t)
         lv = lotkavolterra(u_next, p, t)
-        resid[1] = u_next[1] - u[1] - 0.01*lv[1]
-        resid[2] = u_next[2] - u[2] - 0.01*lv[2]
+        resid[1] = u_next[1] - u[1] - 0.01 * lv[1]
+        resid[2] = u_next[2] - u[2] - 0.01 * lv[2]
         nothing
     end
     u0 = [1.0, 1.0]
     tspan = (0.0, 0.5)
 
     idprob = ImplicitDiscreteProblem(f!, u0, tspan, []; dt = 0.01)
-    idsol = solve(idprob, IDSolve(); adaptive=false)
+    idsol = solve(idprob, IDSolve(); adaptive = false)
 
     oprob = ODEProblem(lotkavolterra, u0, tspan)
     osol = solve(oprob, ImplicitEuler())
@@ -37,15 +37,15 @@ using JET
 
     function g!(resid, u_next, u, p, t)
         f = ff(u_next, p, t)
-        resid[1] = u_next[1] - u[1] - 0.01*f[1]
-        resid[2] = u_next[2] - u[2] - 0.01*f[2]
+        resid[1] = u_next[1] - u[1] - 0.01 * f[1]
+        resid[2] = u_next[2] - u[2] - 0.01 * f[2]
         nothing
     end
     u0 = [10.0, 0.0]
     tspan = (0, 0.2)
 
     idprob = ImplicitDiscreteProblem(g!, u0, tspan, []; dt = 0.01)
-    idsol = solve(idprob, IDSolve(); adaptive=false)
+    idsol = solve(idprob, IDSolve(); adaptive = false)
 
     oprob = ODEProblem(ff, u0, tspan)
     osol = solve(oprob, ImplicitEuler())
@@ -55,7 +55,7 @@ end
 
 @testset "Solver initializes" begin
     function periodic!(resid, u_next, u, p, t)
-        resid[1] = u_next[1] - u[1] - sin(t*π/4)
+        resid[1] = u_next[1] - u[1] - sin(t * π / 4)
         resid[2] = 16 - u_next[2]^2 - u_next[1]^2
     end
 
@@ -73,7 +73,7 @@ end
 
 @testset "Hard problem" begin
     function hard!(resid, u, u_prev, p, t)
-        resid[1] = tanh((u[1]-10t)^2)/2
+        resid[1] = tanh((u[1] - 10t)^2) / 2
     end
 
     u0 = [0.0]
@@ -97,7 +97,8 @@ end
     @test_nowarn integ = init(idprob, IDSolve())
 
     idprob2 = ImplicitDiscreteProblem(emptyoop, u0, (0, tsteps), [])
-    @test_throws AssertionError("Empty u not supported with out of place functions yet.") integ = init(idprob2, IDSolve())
+    @test_throws AssertionError("Empty u not supported with out of place functions yet.") integ=init(
+        idprob2, IDSolve())
 end
 
 @testset "Create NonlinearLeastSquaresProblem" begin

From 43962fbc803e68f1cb1c37e4eef5a0a8547fb561 Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 28 Oct 2025 15:48:14 +0100
Subject: [PATCH 10/11] Add docs and rejection logic

---
 lib/ImplicitDiscreteSolve/src/algorithms.jl |  4 +--
 lib/ImplicitDiscreteSolve/src/controller.jl | 39 ++++++++++++++++++++-
 2 files changed, 39 insertions(+), 4 deletions(-)

diff --git a/lib/ImplicitDiscreteSolve/src/algorithms.jl b/lib/ImplicitDiscreteSolve/src/algorithms.jl
index aa2c567598..a11d41f12a 100644
--- a/lib/ImplicitDiscreteSolve/src/algorithms.jl
+++ b/lib/ImplicitDiscreteSolve/src/algorithms.jl
@@ -7,13 +7,11 @@ struct IDSolve{NLS} <:
        OrdinaryDiffEqAlgorithm
     nlsolve::NLS
     extrapolant::Symbol
-    controller::Symbol
 end
 
 function IDSolve(;
         nlsolve = NewtonRaphson(),
         extrapolant = :constant,
-        controller = :PI
 )
-    IDSolve{typeof(nlsolve)}(nlsolve, extrapolant, controller)
+    IDSolve{typeof(nlsolve)}(nlsolve, extrapolant)
 end
diff --git a/lib/ImplicitDiscreteSolve/src/controller.jl b/lib/ImplicitDiscreteSolve/src/controller.jl
index 63a27c9f84..95550d89eb 100644
--- a/lib/ImplicitDiscreteSolve/src/controller.jl
+++ b/lib/ImplicitDiscreteSolve/src/controller.jl
@@ -1,3 +1,30 @@
+"""
+    KantorovichTypeController
+
+Default controller for implicit discrete solvers. Assuming a Newton method is used
+to solve the nonlinear problem, this controller uses convergence rate estimates to
+adapt the step size based on an posteriori estimate of the change in the Newton
+convergence radius between steps as described below.
+
+Given the convergence rate estimate Θ₀ for the first iteration, the step size controller 
+adapts the time step as `dtₙ₊₁ = γ * (g(Θbar) / (g(Θ₀)))^(1 / p) dtₙ`
+with `g(x) = √(1 + 4x) - 1`. `p` denotes the order of the solver -- i.e. the order of
+the extrapolation algorithm to compute the initial guess for the solve at `tₙ` given a 
+solution at `tₙ₋₁` -- and `Θbar` denotes the  target convergence rate. `γ` is a safety factor.
+
+A factor `Θreject` controls the rejection of a time step if any `Θₖ > Θreject`.
+In this case the first Θₖ violating this criterion is taken and the time step is adapted
+such that `dtₙ₊₁ = γ * (g(Θbar) / (g(Θk)))^(1 / p) dtₙ`. This behavior can be changed
+by setting `strict = false`. In this case the step is accepted whenever the Newton
+solver converges.
+
+The controller furthermore limits the growth and shrinkage of the time step by a factor
+between `qmin` and `qmax`.
+
+The baseline algorithm has been derived in Peter Deuflhard's book "Newton Methods for
+Nonlinear Problems" in Section 5.1.3 (Adaptive pathfollowing algorithms). Please note
+that some implementation details deviate from the original algorithm.
+"""
 Base.@kwdef struct KantorovichTypeController <: OrdinaryDiffEqCore.AbstractController
     Θmin::Float64
     p::Int64
@@ -6,6 +33,7 @@ Base.@kwdef struct KantorovichTypeController <: OrdinaryDiffEqCore.AbstractContr
     γ::Float64 = 0.95
     qmin::Float64 = 1 / 5
     qmax::Float64 = 5.0
+    strict::Bool = true
 end
 
 function OrdinaryDiffEqCore.default_controller(
@@ -40,7 +68,7 @@ function OrdinaryDiffEqCore.step_reject_controller!(
     @inline g(x) = √(1 + 4x) - 1
 
     # Shorten dt according to (Eq. 5.24)
-    (; Θks) = integrator.cache.nlcache
+    (; Θks) = integrator.cache
     (; Θbar, Θreject, γ, Θmin, qmin, qmax, p) = controller
     for Θk in Θks
         if Θk > Θreject
@@ -50,3 +78,12 @@ function OrdinaryDiffEqCore.step_reject_controller!(
         end
     end
 end
+
+function OrdinaryDiffEqCore.accept_step_controller(integrator, controller::KantorovichTypeController)
+    (; Θks) = integrator.cache
+    if controller.strict
+        return all(controller.Θreject .< Θks)
+    else
+        return true
+    end
+end

From 5793442d6a0f56456d7f70fffbe5aeb9d59c29ed Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 28 Oct 2025 15:50:54 +0100
Subject: [PATCH 11/11] Remove linear extrapolant

---
 lib/ImplicitDiscreteSolve/src/solve.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/lib/ImplicitDiscreteSolve/src/solve.jl b/lib/ImplicitDiscreteSolve/src/solve.jl
index d5a511fbab..f6b58de08b 100644
--- a/lib/ImplicitDiscreteSolve/src/solve.jl
+++ b/lib/ImplicitDiscreteSolve/src/solve.jl
@@ -3,10 +3,10 @@ function perform_step!(integrator, cache::IDSolveCache, repeat_step = false)
     (; nlcache, Θks) = cache
 
     # initial guess
-    if alg.extrapolant == :linear
-        @.. broadcast=false cache.z=integrator.u + dt * (integrator.u - integrator.uprev)
-    else # :constant
+    if alg.extrapolant == :constant
         cache.z .= integrator.u
+    else
+        error("Unknown extrapolant $(alg.extrapolant).")
     end
     state = ImplicitDiscreteState(cache.z, p, t + dt)