fix floating point issues in tstops

SciML · Jun 9, 2017 · bd0516b · bd0516b
1 parent 5ba1e83
commit bd0516b
Show file tree

Hide file tree

Showing 2 changed files with 14 additions and 5 deletions.
diff --git a/src/callbacks.jl b/src/callbacks.jl
@@ -32,8 +32,10 @@ end
   # Check if the event occured
   if typeof(callback.idxs) <: Void
     previous_condition = callback.condition(integrator.tprev,integrator.uprev,integrator)
-  else
+  elseif typeof(callback.idxs) <: Number
     previous_condition = callback.condition(integrator.tprev,integrator.uprev[callback.idxs],integrator)
+  else
+    previous_condition = callback.condition(integrator.tprev,@view(integrator.uprev[callback.idxs]),integrator)
   end
   if isapprox(previous_condition,0,rtol=callback.reltol,atol=callback.abstol)
     prev_sign = 0.0
@@ -42,9 +44,11 @@ end
   end
   prev_sign_index = 1
   if typeof(callback.idxs) <: Void
-    next_sign = sign(callback.condition(integrator.tprev+integrator.dt,integrator.u,integrator))
+    next_sign = sign(callback.condition(integrator.t,integrator.u,integrator))
+  elseif typeof(callback.idxs) <: Number
+    next_sign = sign(callback.condition(integrator.t,integrator.u[callback.idxs],integrator))
   else
-    next_sign = sign(callback.condition(integrator.tprev+integrator.dt,integrator.u[callback.idxs],integrator))
+    next_sign = sign(callback.condition(integrator.t,@view(integrator.u[callback.idxs]),integrator))
   end
   if ((prev_sign<0 && !(typeof(callback.affect!)<:Void)) || (prev_sign>0 && !(typeof(callback.affect_neg!)<:Void))) && prev_sign*next_sign<0
     event_occurred = true
@@ -172,7 +176,7 @@ function apply_discrete_callback!(integrator::SDEIntegrator,callback::DiscreteCa
   end
 
   integrator.u_modified = true
-  if callback.condition(integrator.tprev+integrator.dt,integrator.u,integrator)
+  if callback.condition(integrator.t,integrator.u,integrator)
     callback.affect!(integrator)
     if callback.save_positions[2]
       savevalues!(integrator)

diff --git a/src/integrators/integrator_utils.jl b/src/integrators/integrator_utils.jl
@@ -122,7 +122,12 @@ end
     integrator.accept_step = (!integrator.isout && integrator.EEst <= 1.0) || (integrator.opts.force_dtmin && integrator.dt <= integrator.opts.dtmin)
     if integrator.accept_step # Accepted
       integrator.tprev = integrator.t
-      integrator.t = ttmp
+      if typeof(integrator.t)<:AbstractFloat && !isempty(integrator.opts.tstops)
+        tstop = top(integrator.opts.tstops)
+        abs(ttmp - tstop) < 10eps(integrator.t) ? (integrator.t = tstop) : (integrator.t = ttmp)
+      else
+        integrator.t = ttmp
+      end
       calc_dt_propose!(integrator)
       handle_callbacks!(integrator)
     end