Fix Rootfinding in Callbacks

src/callbacks.jl

@@ -536,30 +536,31 @@ function find_callback_time(integrator,callback::ContinuousCallback,counter)
         bottom_θ = typeof(integrator.t)(0)
       end
       if callback.rootfind && !isdiscrete(integrator.alg)
-        zero_func = (Θ) -> begin
-          abst = integrator.tprev+integrator.dt*Θ
+        zero_func = (Θnot) -> begin
+          abst = integrator.tprev+integrator.dt*(1-Θnot)
           return get_condition(integrator, callback, abst)
         end
-        if zero_func(top_Θ) == 0
+        if zero_func(1-top_Θ) == 0
           Θ = top_Θ
         else
           if integrator.event_last_time == counter &&
-            abs(zero_func(bottom_θ)) < 100abs(integrator.last_event_error) &&
+            abs(zero_func(1-bottom_θ)) < 100abs(integrator.last_event_error) &&
             prev_sign_index == 1
 
             # Determined that there is an event by derivative
             # But floating point error may make the end point negative
 
-            sign_top = sign(zero_func(top_Θ))
+            sign_top = sign(zero_func(1-top_Θ))
             bottom_θ += 2eps(typeof(bottom_θ))
             iter = 1
-            while sign(zero_func(bottom_θ)) == sign_top && iter < 12
+            while sign(zero_func(1-bottom_θ)) == sign_top && iter < 12
               bottom_θ *= 5
+              iter += 1
             end
             iter == 12 && error("Double callback crossing floating pointer reducer errored. Report this issue.")
           end
-          Θ = prevfloat(find_zero(zero_func, (bottom_θ,top_Θ), Roots.AlefeldPotraShi(), atol = callback.abstol/100))
-          integrator.last_event_error = ODE_DEFAULT_NORM(zero_func(Θ),integrator.t+integrator.dt*Θ)
+          Θ = 1 - find_zero(zero_func, (1-top_Θ, 1-bottom_θ), Roots.Bisection())
+          integrator.last_event_error = ODE_DEFAULT_NORM(zero_func(1-Θ),integrator.t+integrator.dt*Θ)
         end
         #Θ = prevfloat(...)
         # prevfloat guerentees that the new time is either 1 floating point
@@ -602,33 +603,33 @@ function find_callback_time(integrator,callback::VectorContinuousCallback,counte
         minΘ = nextfloat(top_Θ)
         min_event_idx = -1
         for idx in event_idx
-          zero_func = (Θ) -> begin
-            abst = integrator.tprev+integrator.dt*Θ
+          zero_func = (Θnot) -> begin
+            abst = integrator.tprev+integrator.dt*(1-Θnot)
             return ArrayInterface.allowed_getindex(get_condition(integrator, callback, abst),idx)
           end
-          if zero_func(top_Θ) == 0
+          if zero_func(1-top_Θ) == 0
             Θ = top_Θ
           else
             if integrator.event_last_time == counter &&
               (callback isa VectorContinuousCallback ? integrator.vector_event_last_time == event_idx : true) &&
-              abs(zero_func(bottom_θ)) < 100abs(integrator.last_event_error) &&
+              abs(zero_func(1-bottom_θ)) < 100abs(integrator.last_event_error) &&
               prev_sign_index == 1
 
               # Determined that there is an event by derivative
               # But floating point error may make the end point negative
 
-              sign_top = sign(zero_func(top_Θ))
+              sign_top = sign(zero_func(1-top_Θ))
               bottom_θ += 2eps(typeof(bottom_θ))
               iter = 1
-              while sign(zero_func(bottom_θ)) == sign_top && iter < 12
+              while sign(zero_func(1-bottom_θ)) == sign_top && iter < 12
                 bottom_θ *= 5
                 iter += 1
               end
               iter == 12 && error("Double callback crossing floating pointer reducer errored. Report this issue.")
             end
-            Θ = prevfloat(find_zero(zero_func, (bottom_θ,top_Θ), Roots.AlefeldPotraShi(), atol = callback.abstol/100))
+            Θ = 1 - find_zero(zero_func, (1-top_Θ,1-bottom_θ), Roots.Bisection())
             if Θ < minΘ
-              integrator.last_event_error = ODE_DEFAULT_NORM(zero_func(Θ),integrator.t+integrator.dt*Θ)
+              integrator.last_event_error = ODE_DEFAULT_NORM(zero_func(1-Θ),integrator.t+integrator.dt*Θ)
             end
           end
           if Θ < minΘ

test/downstream/event_detection_tests.jl

@@ -46,6 +46,7 @@ function condition(u,t,integrator) # Event when event_f(u,t) == 0
   u[1]
 end
 function affect!(integrator)
+  @test integrator.u[1] >= 0
   integrator.u[2] = -integrator.u[2]
 end
 cb2 = ContinuousCallback(condition,affect!)
@@ -55,3 +56,74 @@ p = 9.8
 prob = ODEProblem(f,u0,tspan,p)
 sol = solve(prob,Tsit5(),callback=cb2)
 @test minimum(sol') > -40
+
+function vcondition!(out,u,t,integrator)
+  out[1] = u[1]
+  out[2] = u[2]
+end
+
+function vaffect!(integrator, event_idx)
+  @test integrator.u[1] >= 0.0
+  if event_idx == 1
+    integrator.u[2] = -integrator.u[2]
+  else
+    integrator.p = 0.0
+  end
+end
+
+u0 = [50.0,0.0]
+tspan = (0.0,15.0)
+p = 9.8
+prob = ODEProblem(f,u0,tspan,p)
+Vcb = VectorContinuousCallback(vcondition!,vaffect!, 2 , save_positions=(true,true))
+sol = solve(prob,Tsit5(), callback=Vcb)
+
+f = function (du,u,p,t)
+  du[1] = u[2]
+  du[2] = -9.81
+end
+
+function condition(u,t,integrator) # Event when event_f(u,t) == 0
+  u[1]    # Event when height crosses from positive to negative
+end
+
+function affect_neg(integrator)
+  @test integrator.u[1] >= 0.0
+  integrator.u[2] = -0.8*integrator.u[2]
+end
+
+cb = ContinuousCallback(condition,nothing,affect_neg! = affect_neg)
+
+u0 = [1.0,0.0]
+tspan = (0.0, 3.0)
+prob = ODEProblem(f,u0,tspan)
+sol = solve(prob,Tsit5(),saveat=0.01,callback=cb)
+
+f! = function (du,u,p,t)
+  du[1] = u[2]
+  du[2] = -p
+end
+
+function condition!(out,u,t,integrator) 
+  out[1] = u[1]
+  out[2] = u[2]
+end
+
+function affect!(integrator, event_idx)
+  if event_idx == 1
+    integrator.u[2] = -integrator.u[2]
+  else
+    integrator.p = 0.0
+  end
+end
+
+u0 = [50.0,0.0]
+tspan = (0.0,15.0)
+
+begin
+    p = 9.8
+    prob = ODEProblem(f!,u0,tspan,p)
+    Vcb = VectorContinuousCallback(condition!,affect!, 2 , save_positions=(true,true))
+    sol = solve(prob,Tsit5(), callback=Vcb)
+end
+