SciML
Unexpected behaviour with continuous_events
I have encountered some behaviour that I don't think is intended when composing systems with callbacks incorporated in them.
Say we use the example of the bouncing ball:
using ModelingToolkit, OrdinaryDiffEq, Plots
function ball_factory(;name)
@variables t x(t)=1 v(t)=0
D = Differential(t)
root_eqs = [x ~ 0] # the event happens at the ground x(t) = 0
affect = [v ~ -v] # the effect is that the velocity changes sign
ball = ODESystem([
D(x) ~ v
D(v) ~ -9.8
], t;name=name, continuous_events = root_eqs => affect) # equation => affect
end
@named ball2 = ball_factory()
@named ball1 = ball_factory()
@named two_balls = compose(ball1, ball2)
tspan = (0.0,5.0)
prob = ODEProblem(two_balls, Pair[], tspan)
sol = solve(prob,Tsit5())
plot(sol)
On my system, the callbacks for the second system (ball2) are never triggered.
But when I then compose with something else again, e.g.:
function decay(;name)
@variables t Y(t)=10.
D = Differential(t)
ODESystem([
D(Y) ~ - Y
], t; name=name)
end
@named decayer = decay()
combined_sys = compose(decayer, two_balls)
combined_prob = ODEProblem(combined_sys, Pair[], (0.0, 5.0))
combined_sol = solve(combined_prob, Tsit5())
plot(sol, vars=[two_balls.ball2.x])
Things seem to be good again as far as callbacks are concerned. So the problem / unexpected behaviour is probably related to the namespacing, I'm guessing.
The problem appears to be that the compose function sidesteps the ODESystem constructor by using
@set! sys.systems = [get_systems(sys); systems]
The ODESystem constructor is responsible for keeping track of the events. This is likely a problem using discrete equations as well, since those are also handled in the constructor. The solution might be to implement the compose in using the system constructor so that we do not have to copy the code that handles events.
A current workaround for your case is to use the ODESystem(..., systems=[ball1, ball2]) form instead of compose
Edit: Actually, scrap that. The callback generation extracts events from all subsystems. I'll look into it
The callbacks are correctly namespaced
julia> ModelingToolkit.continuous_events(two_balls)
2-element Vector{ModelingToolkit.SymbolicContinuousCallback}:
ModelingToolkit.SymbolicContinuousCallback(Equation[x(t) ~ 0], Equation[v(t) ~ -v(t)])
ModelingToolkit.SymbolicContinuousCallback(Equation[ball2₊x(t) ~ 0], Equation[ball2₊v(t) ~ -ball2₊v(t)])
julia> two_balls
Model two_balls with 4 equations
States (4):
x(t) [defaults to 1]
v(t) [defaults to 0]
ball2₊x(t) [defaults to 1]
ball2₊v(t) [defaults to 0]
Parameters (0):
but it appears as if only the first condition ever triggers...
I wonder if the reason is that both events trigger at exactly the same time, and the VectorContinuousCallback only triggers the affect! once, using one index. @ChrisRackauckas if two elements of the condition vector hits 0 at the same time, is affect!(integrator, event_index) guaranteed to be called with both event_index that hit zero?
It does seem to be realted, because if I change the initial condition slightly so that the events occur at different times, I get
prob = ODEProblem(two_balls, Pair[@nonamespace(two_balls.x) => 1.1], tspan)
julia> sol = solve(prob,Tsit5())
ERROR: Double callback crossing floating pointer reducer errored. Report this issue.
Stacktrace:
[1] error(s::String)
@ Base ./error.jl:33
[2] find_callback_time(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Nothing, Float64, Nothing, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Nothing, ODEFunction{true, ModelingToolkit.var"#f#300"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xb0ed6ebf, 0xd944739c, 0xbcf44059, 0x676aefbe, 0x33c6aa61)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xcd20b0e4, 0xd8d63225, 0x59bca120, 0xbaf9c407, 0xbefd8714)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#291#generated_observed#307"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Pairs{Symbol, VectorContinuousCallback{ModelingToolkit.var"#742#752"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x5b99be80, 0x04e6c94c, 0x607fbd44, 0x97508608, 0x94f5fbd2)}}, ModelingToolkit.var"#744#754"{Vector{ModelingToolkit.var"#276#283"{rf_ip, Vector{Int64}} where rf_ip}, Vector{Int64}}, ModelingToolkit.var"#744#754"{Vector{ModelingToolkit.var"#276#283"{rf_ip, Vector{Int64}} where rf_ip}, Vector{Int64}}, typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{VectorContinuousCallback{ModelingToolkit.var"#742#752"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x5b99be80, 0x04e6c94c, 0x607fbd44, 0x97508608, 0x94f5fbd2)}}, ModelingToolkit.var"#744#754"{Vector{ModelingToolkit.var"#276#283"{rf_ip, Vector{Int64}} where rf_ip}, Vector{Int64}}, ModelingToolkit.var"#744#754"{Vector{ModelingToolkit.var"#276#283"{rf_ip, Vector{Int64}} where rf_ip}, Vector{Int64}}, typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}}}}, SciMLBase.StandardODEProblem}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, ModelingToolkit.var"#f#300"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xb0ed6ebf, 0xd944739c, 0xbcf44059, 0x676aefbe, 0x33c6aa61)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xcd20b0e4, 0xd8d63225, 0x59bca120, 0xbaf9c407, 0xbefd8714)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#291#generated_observed#307"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats}, ODEFunction{true, ModelingToolkit.var"#f#300"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xb0ed6ebf, 0xd944739c, 0xbcf44059, 0x676aefbe, 0x33c6aa61)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xcd20b0e4, 0xd8d63225, 0x59bca120, 0xbaf9c407, 0xbefd8714)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#291#generated_observed#307"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(opnorm), Nothing, CallbackSet{Tuple{VectorContinuousCallback{ModelingToolkit.var"#742#752"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x5b99be80, 0x04e6c94c, 0x607fbd44, 0x97508608, 0x94f5fbd2)}}, ModelingToolkit.var"#744#754"{Vector{ModelingToolkit.var"#276#283"{rf_ip, Vector{Int64}} where rf_ip}, Vector{Int64}}, ModelingToolkit.var"#744#754"{Vector{ModelingToolkit.var"#276#283"{rf_ip, Vector{Int64}} where rf_ip}, Vector{Int64}}, typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, DiffEqBase.CallbackCache{Vector{Float64}, Vector{Float64}}, OrdinaryDiffEq.DefaultInit}, callback::VectorContinuousCallback{ModelingToolkit.var"#742#752"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x5b99be80, 0x04e6c94c, 0x607fbd44, 0x97508608, 0x94f5fbd2)}}, ModelingToolkit.var"#744#754"{Vector{ModelingToolkit.var"#276#283"{rf_ip, Vector{Int64}} where rf_ip}, Vector{Int64}}, ModelingToolkit.var"#744#754"{Vector{ModelingToolkit.var"#276#283"{rf_ip, Vector{Int64}} where rf_ip}, Vector{Int64}}, typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT), Float64, Int64, Rational{Int64}, Nothing, Int64}, counter::Int64)
@ DiffEqBase ~/.julia/packages/DiffEqBase/b1nST/src/callbacks.jl:766
[3] find_first_continuous_callback
I tried to test if that was the case by changing the gravitational constant (which changes callback times) for the balls by doing
function ball_factory(;name, g = -9.8)
@variables t x(t)=1 v(t)=0
@parameters g=g
D = Differential(t)
root_eqs = [x ~ 0] # the event happens at the ground x(t) = 0
affect = [v ~ -v] # the effect is that the velocity changes sign
ball = ODESystem([
D(x) ~ v
D(v) ~ g
], t;name=name, continuous_events = root_eqs => affect) # equation => affect
end
@named ball1 = ball_factory(g = -8.)
@named ball2 = ball_factory()
@named two_balls = compose(ball1, ball2)
etc., but this gave an error: Double callback crossing floating pointer reducer errored. Report this issue.
Edit: Seems you found this already :)
That's an open issue IIRC for simultaneous events. The min_event_idx is expected to be an integer, and we need to have something a bit more interesting for that to change:
https://github.com/SciML/DiffEqBase.jl/blob/master/src/callbacks.jl#L747-L806
I think we'd need to extend the API too so that affect! takes the vector of bools for the yes/no of whether it's triggered at that point (if it was a vector of ints then it would allocate every time 😅, unless we do some tricky view to an array that already exists).
