Merge pull request SciML#2515 from SciML/predictive

Refactor predictive controller
Shreyas-Ekanathan · Nov 10, 2024 · 6f18857 · 6f18857
2 parents f65bcca + 19312aa
commit 6f18857
Show file tree

Hide file tree

Showing 2 changed files with 53 additions and 54 deletions.
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
@@ -400,3 +400,56 @@ function post_newton_controller!(integrator, alg)
     integrator.dt = integrator.dt / integrator.opts.failfactor
     nothing
 end
+
+@inline function stepsize_controller!(integrator, controller::PredictiveController, alg)
+    @unpack qmin, qmax, gamma = integrator.opts
+    EEst = DiffEqBase.value(integrator.EEst)
+    if iszero(EEst)
+        q = inv(qmax)
+    else
+        if fac_default_gamma(alg)
+            fac = gamma
+        else
+            if isfirk(alg)
+                @unpack iter = integrator.cache
+                @unpack maxiters = alg
+            else
+                @unpack iter, maxiters = integrator.cache.nlsolver
+            end
+            fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
+        end
+        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
+        qtmp = FastPower.fastpower(EEst, expo) / fac
+        @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
+        integrator.qold = q
+    end
+    q
+end
+
+function step_accept_controller!(integrator, controller::PredictiveController, alg, q)
+    @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
+
+    EEst = DiffEqBase.value(integrator.EEst)
+
+    if integrator.success_iter > 0
+        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
+        qgus = (integrator.dtacc / integrator.dt) *
+               FastPower.fastpower((EEst^2) / integrator.erracc, expo)
+        qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
+        qacc = max(q, qgus)
+    else
+        qacc = q
+    end
+    if qsteady_min <= qacc <= qsteady_max
+        qacc = one(qacc)
+    end
+    integrator.dtacc = integrator.dt
+    integrator.erracc = max(1e-2, EEst)
+
+    return integrator.dt / qacc
+end
+
+function step_reject_controller!(integrator, controller::PredictiveController, alg)
+    @unpack dt, success_iter, qold = integrator
+    integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
+end
diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
@@ -1,52 +1,3 @@
-@inline function stepsize_controller!(integrator, controller::PredictiveController, alg)
-    @unpack qmin, qmax, gamma = integrator.opts
-    EEst = DiffEqBase.value(integrator.EEst)
-    if iszero(EEst)
-        q = inv(qmax)
-    else
-        if fac_default_gamma(alg)
-            fac = gamma
-        else
-            if isfirk(alg)
-                @unpack iter = integrator.cache
-                @unpack maxiters = alg
-            else
-                @unpack iter, maxiters = integrator.cache.nlsolver
-            end
-            fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
-        end
-        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
-        qtmp = FastPower.fastpower(EEst, expo) / fac
-        @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
-        integrator.qold = q
-    end
-    q
-end
-
-function step_accept_controller!(integrator, controller::PredictiveController, alg, q)
-    @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
-
-    EEst = DiffEqBase.value(integrator.EEst)
-
-    if integrator.success_iter > 0
-        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
-        qgus = (integrator.dtacc / integrator.dt) *
-               FastPower.fastpower((EEst^2) / integrator.erracc, expo)
-        qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
-        qacc = max(q, qgus)
-    else
-        qacc = q
-    end
-    if qsteady_min <= qacc <= qsteady_max
-        qacc = one(qacc)
-    end
-    integrator.dtacc = integrator.dt
-    integrator.erracc = max(1e-2, EEst)
-
-    return integrator.dt / qacc
-end
-
-
 function step_accept_controller!(integrator, controller::PredictiveController, alg::AdaptiveRadau, q)
     @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = integrator.opts
     @unpack cache = integrator
@@ -85,11 +36,6 @@ function step_accept_controller!(integrator, controller::PredictiveController, a
     return integrator.dt / qacc
 end
 
-function step_reject_controller!(integrator, controller::PredictiveController, alg)
-    @unpack dt, success_iter, qold = integrator
-    integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
-end
-
 function step_reject_controller!(integrator, controller::PredictiveController, alg::AdaptiveRadau)
     @unpack dt, success_iter, qold = integrator
     @unpack cache = integrator