Support multiple simultaneous events in VectorContinuousCallback

lib/DiffEqBase/src/callbacks.jl

@@ -196,41 +196,66 @@ end
     bottom_condition = get_condition(integrator, callback, integrator.tprev)
     @. bottom_sign = sign(bottom_condition)
 
+    prev_simultaneous_events = integrator.callback_cache.prev_simultaneous_events
     if integrator.event_last_time == callback_idx
-        nudged_idx = integrator.vector_event_last_time
         # If there was a previous event, nudge tprev on the right
         # side of the root (if necessary) to avoid repeat detection
 
         if callback.interp_points == 0
             addsteps!(integrator)
         end
 
+        # Find the condition value closest to zero across all triggered events
+        min_condition_val = zero(eltype(bottom_condition))
+        min_abs_condition = typemax(eltype(bottom_condition))
+        for idx in 1:callback.len
+            if prev_simultaneous_events[idx]
+                cond_val = ArrayInterface.allowed_getindex(bottom_condition, idx)
+                if abs(cond_val) < min_abs_condition
+                    min_abs_condition = abs(cond_val)
+                    min_condition_val = cond_val
+                end
+            end
+        end
+
         # Evaluate condition slightly in future
-        nudged_t = nudge_tprev(integrator, callback, ArrayInterface.allowed_getindex(bottom_condition, nudged_idx))
+        nudged_t = nudge_tprev(integrator, callback, min_condition_val)
         tmp_condition = get_condition(integrator, callback, nudged_t)
 
-        ArrayInterface.allowed_setindex!(bottom_sign, sign(ArrayInterface.allowed_getindex(tmp_condition, nudged_idx)), nudged_idx)
+        for idx in 1:callback.len
+            if prev_simultaneous_events[idx]
+                ArrayInterface.allowed_setindex!(bottom_sign, sign(ArrayInterface.allowed_getindex(tmp_condition, idx)), idx)
+            end
+        end
     else
-        nudged_idx = -1
         nudged_t = bottom_t
     end
 
     # Check if an event occured
     event_occurred, event_idx, top_t, top_sign =
         check_event_occurence(integrator, callback, bottom_sign)
 
+    # Track simultaneous events
+    (; simultaneous_events) = integrator.callback_cache
+    if event_occurred
+        @. prev_simultaneous_events = !iszero(simultaneous_events)
+        simultaneous_events .= Int8(0)
+    end
+
     # Find callback time if occurence
     if !event_occurred
         callback_t = integrator.t
         min_event_idx = 1
         residual = zero(eltype(bottom_condition))
     elseif isdiscrete(integrator.alg) || callback.rootfind == SciMLBase.NoRootFind
         callback_t = top_t
-        min_event_idx = 1
+        min_event_idx = -1
         for i in 1:length(event_idx)
             if ArrayInterface.allowed_getindex(event_idx, i) == 1
-                min_event_idx = i
-                break
+                if min_event_idx < 0
+                    min_event_idx = i
+                end
+                simultaneous_events[i] = Int8(sign(ArrayInterface.allowed_getindex(bottom_sign, i)))
             end
         end
         residual = zero(eltype(bottom_condition))
@@ -251,16 +276,20 @@ end
                 if iszero(ArrayInterface.allowed_getindex(top_sign, idx))
                     cbi_t = top_t
                 else
-                    if idx == nudged_idx
+                    if integrator.event_last_time == callback_idx && prev_simultaneous_events[idx]
                         cbi_t = find_root(zero_func, (nudged_t, top_t), callback.rootfind)
                     else
                         cbi_t = find_root(zero_func, (bottom_t, top_t), callback.rootfind)
                     end
                 end
                 if integrator.tdir * cbi_t < integrator.tdir * callback_t
+                    simultaneous_events .= Int8(0)
+                end
+                if integrator.tdir * cbi_t <= integrator.tdir * callback_t
                     min_event_idx = idx
                     callback_t = cbi_t
                     residual = zero_func(cbi_t)
+                    simultaneous_events[idx] = Int8(sign(ArrayInterface.allowed_getindex(bottom_sign, idx)))
                 end
             end
         end
@@ -270,7 +299,8 @@ end
         end
     end
 
-    return callback_t, ArrayInterface.allowed_getindex(bottom_sign, min_event_idx),
+    # We still pass around the min_event_idx for now because some stuff in OrdinaryDiffEqCore expects it to be an Int
+    return callback_t, bottom_sign,
         event_occurred::Bool, min_event_idx::Int, residual
 end
 
@@ -444,6 +474,31 @@ function is_event_occurence(prev_sign::Number, next_sign::Number, affect!::F1, a
     ) && prev_sign * next_sign <= 0
 end
 
+"""
+    apply_callback!(integrator, callback, cb_time, prev_sign, event_idx)
+
+Apply a continuous callback at the determined event time.
+
+For `ContinuousCallback`, the `affect!` or `affect_neg!` function is called based on the
+crossing direction (`prev_sign`):
+  - `prev_sign < 0` (upcrossing): `callback.affect!(integrator)` is called
+  - `prev_sign > 0` (downcrossing): `callback.affect_neg!(integrator)` is called
+
+For `VectorContinuousCallback`, `callback.affect!` is called once with the full
+`simultaneous_events::Vector{Int8}` array from the callback cache:
+
+    callback.affect!(integrator, simultaneous_events)
+
+Each element of `simultaneous_events` encodes both whether the event triggered and
+its crossing direction:
+  - `0`: event did not trigger
+  - `-1`: event triggered via upcrossing (condition went from negative to positive)
+  - `+1`: event triggered via downcrossing (condition went from positive to negative)
+
+Multiple events may be nonzero simultaneously when they occur at the same time.
+The `affect_neg!` field is not called for `VectorContinuousCallback`; the user's
+`affect!` function should handle both crossing directions using the sign information.
+"""
 function apply_callback!(
         integrator,
         callback::Union{ContinuousCallback, VectorContinuousCallback},
@@ -474,19 +529,25 @@ function apply_callback!(
 
     integrator.u_modified = true
 
-    if prev_sign < 0
+    if callback isa VectorContinuousCallback
         if callback.affect! === nothing
             integrator.u_modified = false
         else
-            callback isa VectorContinuousCallback ?
-                callback.affect!(integrator, event_idx) : callback.affect!(integrator)
+            callback.affect!(integrator, integrator.callback_cache.simultaneous_events)
         end
-    elseif prev_sign > 0
-        if callback.affect_neg! === nothing
-            integrator.u_modified = false
-        else
-            callback isa VectorContinuousCallback ?
-                callback.affect_neg!(integrator, event_idx) : callback.affect_neg!(integrator)
+    else
+        if prev_sign < 0
+            if callback.affect! === nothing
+                integrator.u_modified = false
+            else
+                callback.affect!(integrator)
+            end
+        elseif prev_sign > 0
+            if callback.affect_neg! === nothing
+                integrator.u_modified = false
+            else
+                callback.affect_neg!(integrator)
+            end
         end
     end
 
@@ -610,6 +671,8 @@ mutable struct CallbackCache{conditionType, signType}
     next_condition::conditionType
     next_sign::signType
     prev_sign::signType
+    simultaneous_events::Vector{Int8}
+    prev_simultaneous_events::Vector{Bool}
 end
 
 function CallbackCache(
@@ -620,7 +683,10 @@ function CallbackCache(
     next_condition = similar(u, conditionType, max_len)
     next_sign = similar(u, signType, max_len)
     prev_sign = similar(u, signType, max_len)
-    return CallbackCache(tmp_condition, next_condition, next_sign, prev_sign)
+    simultaneous_events = zeros(Int8, max_len)
+    prev_simultaneous_events = zeros(Bool, max_len)
+    return CallbackCache(tmp_condition, next_condition, next_sign, prev_sign,
+        simultaneous_events, prev_simultaneous_events)
 end
 
 function CallbackCache(
@@ -631,5 +697,8 @@ function CallbackCache(
     next_condition = zeros(conditionType, max_len)
     next_sign = zeros(signType, max_len)
     prev_sign = zeros(signType, max_len)
-    return CallbackCache(tmp_condition, next_condition, next_sign, prev_sign)
+    simultaneous_events = zeros(Int8, max_len)
+    prev_simultaneous_events = zeros(Bool, max_len)
+    return CallbackCache(tmp_condition, next_condition, next_sign, prev_sign,
+        simultaneous_events, prev_simultaneous_events)
 end

lib/DiffEqBase/test/callbacks.jl

@@ -72,8 +72,11 @@ function DiffEqBase.find_callback_time(
     return 1.0 + counter, 0.9 + counter, true, counter, 0.0
 end
 find_first_integrator = EmptyIntegrator([1.0, 2.0], 1, 0.0)
-vector_affect! = function (integrator, idx)
-    return integrator.u = integrator.u + idx
+vector_affect! = function (integrator, events)
+    for (idx, dir) in enumerate(events)
+        iszero(dir) && continue
+        integrator.u = integrator.u .+ idx
+    end
 end
 
 cond_1(u, t, integrator) = t - 1.0

lib/DiffEqBase/test/downstream/callback_detection.jl

@@ -9,7 +9,12 @@ using OrdinaryDiffEq
 
     vcb = VectorContinuousCallback(
         (out, u, t, integrator) -> out .= (t - 1.0e-8, t - 2.0e-8, t - 2.0e-7),
-        (integrator, event_index) -> push!(record, event_index),
+        (integrator, events) -> begin
+            for (idx, dir) in enumerate(events)
+                iszero(dir) && continue
+                push!(record, idx)
+            end
+        end,
         3;
         abstol = 0.0
     )

lib/DiffEqBase/test/downstream/community_callback_tests.jl

@@ -153,37 +153,40 @@ function condition(out, u, t, integrator)
     return
 end
 
-function collision_affect!(integrator, idx)
-    i = 0
+function collision_affect!(integrator, events)
     u = integrator.u
     n = length(u.nodes)
-    return for k in 1:n
-        for l in (k + 1):n
-            i += 1
-            if idx == i
-                x₁ = u.nodes[k][1:2]
-                v₁ = u.nodes[k][3:4]
-                x₂ = u.nodes[l][1:2]
-                v₂ = u.nodes[l][3:4]
-                # https://stackoverflow.com/a/35212639
-                v₁ = (v₁ - 2 / (1 + 1) * (dot(v₁ - v₂, x₁ - x₂) / sum(abs2, x₁ - x₂) * (x₁ - x₂)))
-                v₂ = -(v₂ - 2 / (1 + 1) * (dot(v₂ - v₁, x₂ - x₁) / sum(abs2, x₂ - x₁) * (x₂ - x₁)))
-
-                println("Collision handled.")
-
-                m = (x₁ + x₂) / 2
-
-                u.nodes[k][3:4] .= v₁
-                u.nodes[l][3:4] .= v₂
-
-                set_u!(integrator, u)
-                println(sqrt(sum(abs2, x₁ .- x₂)) - 100, ":", v₁ ./ v₂)
-                println(
-                    norm(v₁), ":", norm(v₂), ":", integrator.t, ":",
-                    integrator.t - t_last
-                )
-                global t_last = integrator.t
-                break
+    for (event_idx, dir) in enumerate(events)
+        iszero(dir) && continue
+        i = 0
+        for k in 1:n
+            for l in (k + 1):n
+                i += 1
+                if event_idx == i
+                    x₁ = u.nodes[k][1:2]
+                    v₁ = u.nodes[k][3:4]
+                    x₂ = u.nodes[l][1:2]
+                    v₂ = u.nodes[l][3:4]
+                    # https://stackoverflow.com/a/35212639
+                    v₁ = (v₁ - 2 / (1 + 1) * (dot(v₁ - v₂, x₁ - x₂) / sum(abs2, x₁ - x₂) * (x₁ - x₂)))
+                    v₂ = -(v₂ - 2 / (1 + 1) * (dot(v₂ - v₁, x₂ - x₁) / sum(abs2, x₂ - x₁) * (x₂ - x₁)))
+
+                    println("Collision handled.")
+
+                    m = (x₁ + x₂) / 2
+
+                    u.nodes[k][3:4] .= v₁
+                    u.nodes[l][3:4] .= v₂
+
+                    set_u!(integrator, u)
+                    println(sqrt(sum(abs2, x₁ .- x₂)) - 100, ":", v₁ ./ v₂)
+                    println(
+                        norm(v₁), ":", norm(v₂), ":", integrator.t, ":",
+                        integrator.t - t_last
+                    )
+                    global t_last = integrator.t
+                    break
+                end
             end
         end
     end
@@ -217,8 +220,8 @@ function cond!(out, u, t, i)
     out[1] = u[3]
     return nothing
 end
-function terminate_affect!(int, idx)
-    return terminate!(int)
+function terminate_affect!(int, events)
+    any(!iszero, events) && terminate!(int)
 end
 cb = VectorContinuousCallback(cond!, terminate_affect!, nothing, 1)
 

lib/DiffEqBase/test/runtests.jl

@@ -10,7 +10,7 @@ function activate_downstream_env()
     Pkg.develop(
         [
             PackageSpec(path = dirname(@__DIR__)),
-            PackageSpec(path = joinpath(dirname(lib_dir), "..")),
+            PackageSpec(path = dirname(lib_dir)),
             PackageSpec(path = joinpath(lib_dir, "StochasticDiffEq")),
         ]
     )
@@ -29,7 +29,7 @@ function activate_modelingtoolkit_env()
     Pkg.develop(
         [
             PackageSpec(path = dirname(@__DIR__)),
-            PackageSpec(path = joinpath(dirname(lib_dir), "..")),
+            PackageSpec(path = dirname(lib_dir)),
         ]
     )
     return Pkg.instantiate()
@@ -40,7 +40,7 @@ function activate_sundials_env()
     lib_dir = dirname(dirname(@__DIR__))
     Pkg.develop(
         [
-            PackageSpec(path = joinpath(dirname(lib_dir), "..")),
+            PackageSpec(path = dirname(lib_dir)),
         ]
     )
     return Pkg.instantiate()