Merge pull request #35 from SciML/dw/discrete

Fix #34

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "StochasticDelayDiffEq"
 uuid = "29a0d76e-afc8-11e9-03a4-eda52ae4b960"
 authors = ["Henrik Sykora <[email protected]>"]
-version = "0.2.5"
+version = "0.2.6"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"

src/integrators/interface.jl

@@ -221,8 +221,28 @@ end
 
 DiffEqBase.change_t_via_interpolation!(integrator::SDDEIntegrator,t,modify_save_endpoint::Type{Val{T}}=Val{false}) where T = StochasticDiffEq.change_t_via_interpolation!(integrator, t, modify_save_endpoint)
 
-function u_modified!(integrator::SDDEIntegrator,bool::Bool)
-  integrator.u_modified = bool
+# update integrator when u is modified by callbacks
+function StochasticDiffEq.handle_callback_modifiers!(integrator::SDDEIntegrator)
+    # copied from StochasticDiffEq
+    if integrator.P !== nothing && integrator.opts.adaptive
+        if integrator.cache isa StochasticDiffEqMutableCache
+            oldrate = integrator.P.cache.currate
+            integrator.P.cache.rate(oldrate, integrator.u, integrator.p, integrator.t)
+        else
+            integrator.P.cache.currate = integrator.P.cache.rate(
+                integrator.u, integrator.p, integrator.t
+            )
+        end
+    end
+
+    # update heap of discontinuities
+    # discontinuity is assumed to be of order 0, i.e. solution x is discontinuous
+    push!(integrator.opts.d_discontinuities, Discontinuity(integrator.tdir * integrator.t, 0))
+end
+
+function DiffEqBase.u_modified!(integrator::SDDEIntegrator, bool::Bool)
+    integrator.u_modified = bool
+    nothing
 end
 
 get_proposed_dt(integrator::SDDEIntegrator) = integrator.dtpropose

test/events.jl

@@ -0,0 +1,21 @@
+using StochasticDelayDiffEq
+using Test
+
+@testset "discrete callback (#34)" begin
+    f(u, h, p, t) = 0.3 * h(p, t - p)
+    g(u, h, p, t) = 0.1
+    h(p, t) = 0.0
+    prob = SDDEProblem(f, g, h, (0.0, 4.0), 0.5; constant_lags=(0.5,))
+
+    # event at `t = 3`
+    cb = DiscreteCallback(
+        (u, t, integrator) -> t == 3,
+        integrator -> (integrator.u = -integrator.u)
+    )
+
+    sol = solve(prob, RKMil(), tstops=(3,), callback=cb)
+    ts = findall(x -> x == 3, sol.t)
+    @test length(ts) == 2
+    @test sol.u[ts[1]] == -sol.u[ts[2]]
+    @test sol(3.0; continuity=:right) == -sol(3.0; continuity=:left)
+end

test/runtests.jl

@@ -2,3 +2,4 @@ using SafeTestsets
 
 @safetestset "SDDEProblem, solve" begin include("test_prob_sol.jl") end
 @safetestset "Analyticless Convergence Tests" begin include("analyticless_convergence_tests.jl") end
+@safetestset "Event handling" begin include("events.jl") end