bdf fix

Author: Shreyas-Ekanathan

lib/OrdinaryDiffEqBDF/src/controllers.jl

@@ -161,7 +161,27 @@ function step_accept_controller!(integrator, cache::Union{QNDFCache, QNDFConstan
             h / q
         end
     end
-    return is_disco ? min((integrator.disco_checkpoint - integrator.t) / 4, new_dt) : new_dt
+    if is_disco
+        return min((integrator.disco_checkpoint - integrator.t) / 4, new_dt)
+    end
+    return new_dt
+end
+
+"""
+    bdf_restart_estimates!(cache)
+
+Restart the order and constant-step estimates for a step that straddles a derivative
+discontinuity, where the history behind it and the solution ahead belong to different
+regimes. We do this because the BDF step-size and order logic is based on the history of the solution, 
+and we have effectively entered a new regime where old estimates no longer apply.
+"""
+function bdf_restart_estimates!(cache)
+    cache.order = 1
+    cache.nconsteps = 0
+    if hasfield(typeof(cache), :qwait)
+        cache.qwait = 3
+    end
+    return nothing
 end
 
 function bdf_step_reject_controller!(integrator, cache, EEst1)
@@ -438,7 +458,11 @@ function step_accept_controller!(
         cache.qwait -= 1 # countdown
     end
     new_dt = integrator.dt / q
-    return is_disco ? min((integrator.disco_checkpoint - integrator.t) / 4, new_dt) : new_dt
+    if is_disco
+        bdf_restart_estimates!(cache)
+        return min((integrator.disco_checkpoint - integrator.t) / 4, new_dt)
+    end
+    return new_dt
 end
 
 function step_reject_controller!(integrator, alg::DFBDF)
@@ -607,5 +631,9 @@ function step_accept_controller!(
         cache.qwait -= 1 # countdown
     end
     new_dt = integrator.dt / q
-    return is_disco ? min((integrator.disco_checkpoint - integrator.t) / 4, new_dt) : new_dt
+    if is_disco
+        bdf_restart_estimates!(cache)
+        return min((integrator.disco_checkpoint - integrator.t) / 4, new_dt)
+    end
+    return new_dt
 end

lib/OrdinaryDiffEqCore/src/disco.jl

@@ -24,6 +24,7 @@ function find_discontinuity(integrator)
     dt = integrator.dt
     u = integrator.u
     uprev = integrator.uprev
+    integrator.is_disco_step = false
     cb === nothing && return one(dt)
     isempty(cb.continuous_callbacks) && return one(dt)
     p = integrator.p

test/Integrators_I/disco_tests.jl

@@ -145,13 +145,6 @@ prob_stiff = ODEProblem(f_stiff_disc!, u0_stiff, tspan_stiff, [500.0])
 cond_stiff(u, t, integrator) = u[1] - 0.5
 cb_stiff = ContinuousCallback(cond_stiff, default_affect!; maybe_discontinuity = true)
 
-sol_disco_radau = solve(prob_stiff, RadauIIA5(); callback=cb_stiff, reltol=1e-9, abstol=1e-11, controller = predictive_disco_controller(RadauIIA5()))
-#  207.500 μs (2285 allocations: 101.72 KiB)
-sol_no_disco_radau = solve(prob_stiff, RadauIIA5(); callback=cb_stiff, reltol=1e-9, abstol=1e-11)
-#  196.458 μs (2054 allocations: 95.58 KiB)
-@test sol_disco_radau.retcode == ReturnCode.Success
-@test sol_disco_radau.stats.nreject <= sol_no_disco_radau.stats.nreject
-
 sol_disco_rosenbrock = solve(prob_stiff, Rodas5P(); callback=cb_stiff, reltol=1e-9, abstol=1e-11, controller = PI_disco_controller(Rodas5P()))
 #  364.833 μs (2030 allocations: 85.34 KiB)
 sol_no_disco_rosenbrock = solve(prob_stiff, Rodas5P(); callback=cb_stiff, reltol=1e-9, abstol=1e-11)