DiffEqCallbacks.jl
DiffEqCallbacks.jl copied to clipboard
`PeriodicCallback` is seemingly not thread-safe in `EnsembleProblem` with `EnsembleThreads`
See the minimum example below, which loops until resulting in the shown crash.
I did some simple debugging (looking at the thread IDs) and found that PeriodicCallback seems to have a race issue with the index[] value getting reset (in the initialisation function) by different threads, leading to another thread trying to add a new tstop value in the wrong place (i.e. where tdir_tnew < integrator.t).
I can fix this issue by changing tnew = t0[] + (index[] + 1) * Δt (in PeriodicCallback function) to tnew = integrator.t + Δt -- however, I think the issue probably still remains in the condition function. So it probably needs a better thought out solution.
Could it be simply that EnsembleProblem is sharing callbacks between the threads, rather than making copies? I didn't look into this. [UPDATE: I did. See comments below. This particular error only occurs on OrbitalTrajectories v0.1.10.)
Status `Project.toml`
[459566f4] DiffEqCallbacks v2.16.1
[0c46a032] DifferentialEquations v6.16.0
[31c24e10] Distributions v0.24.18
[2b613a20] OrbitalTrajectories v0.1.10
[1dea7af3] OrdinaryDiffEq v5.52.4
[0bca4576] SciMLBase v1.13.0
using OrbitalTrajectories
using DifferentialEquations
using Distributions
# Nominal initial spacecraft orbital state
system = CR3BP(:earth, :moon)
u0 = [0.7671053516112541, 0.0, 0.0, 0.0, 0.4859512921111289, 0.0]
state = State(system, u0, (0.0, 4.381287609345167))
# (OI) Orbit insertion errors
SD_OI = [1.2226878842036528e-6, 1.2226878842036528e-6, 3.303858750933275e-6, 8.512720156555774e-6, 1.0469667318982388e-5, 1.761252446183953e-5]
N_OI = Normal.(u0, SD_OI)
σ_OI = truncated.(N_OI, u0 - 2*SD_OI, u0 + 2*SD_OI)
# Function to generate a stochastic sample initial spacecraft state
prob_func(prob,i,repeat) = remake(prob; u0=rand.(σ_OI))
# Periodic callback that does nothing
Nt = 10
dt = (state.tspan[end] - state.tspan[begin]) / Nt
callback = PeriodicCallback(integrator -> nothing, dt)
# Simulate a stochastic ensemble of orbits
ensemble_prob = EnsembleProblem(state; prob_func)
while true # Loop until we crash due to race condition
# NOTE: trajectories=1 does not lead to a crash. trajectories>1 does.
ensemble_orbits = solve(ensemble_prob, Vern7(), EnsembleThreads(); trajectories=2, callback)
end
ERROR: TaskFailedException
Stacktrace:
[1] wait
@ .\task.jl:322 [inlined]
[2] threading_run(func::Function)
@ Base.Threads .\threadingconstructs.jl:34
[3] macro expansion
@ .\threadingconstructs.jl:93 [inlined]
[4] tmap(f::Function, args::UnitRange{Int64})
@ SciMLBase C:\Users\Dan\.julia\packages\SciMLBase\9EjAY\src\ensemble\basic_ensemble_solve.jl:220
[5] solve_batch(prob::EnsembleProblem{State{CR3BP, SynodicFrame{true}, ODEProblem{StaticArrays.MVector{6, Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::Vern7, ensemblealg::EnsembleThreads, II::UnitRange{Int64}, pmap_batch_size::Int64; kwargs::Base.Iterators.Pairs{Symbol, DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}}}})
@ SciMLBase C:\Users\Dan\.julia\packages\SciMLBase\9EjAY\src\ensemble\basic_ensemble_solve.jl:212
[6] macro expansion
@ .\timing.jl:287 [inlined]
[7] __solve(prob::EnsembleProblem{State{CR3BP, SynodicFrame{true}, ODEProblem{StaticArrays.MVector{6, Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::Vern7, ensemblealg::EnsembleThreads; trajectories::Int64, batch_size::Int64, pmap_batch_size::Int64, kwargs::Base.Iterators.Pairs{Symbol, DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}}}})
@ SciMLBase C:\Users\Dan\.julia\packages\SciMLBase\9EjAY\src\ensemble\basic_ensemble_solve.jl:108
[8] #solve#59
@ C:\Users\Dan\.julia\packages\DiffEqBase\jhLIm\src\solve.jl:96 [inlined]
[9] top-level scope
@ .\REPL[20]:4
nested task error: Tried to add a tstop that is behind the current time. This is strictly forbidden
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:33
[2] add_tstop!
@ C:\Users\Dan\.julia\packages\OrdinaryDiffEq\vxMSM\src\integrators\integrator_interface.jl:100 [inlined]
[3] (::DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}})(integrator::OrdinaryDiffEq.ODEIntegrator{Vern7, true, Vector{Float64}, Nothing, Float64, Vector{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, Vector{Float64}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Vern7, OrdinaryDiffEq.InterpolationData{ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Vern7Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Vern7Tableau{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Vern7Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Vern7Tableau{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Dict{Any, Any}, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
@ DiffEqCallbacks C:\Users\Dan\.julia\dev\DiffEqCallbacks\src\iterative_and_periodic.jl:93
[4] apply_discrete_callback!
@ C:\Users\Dan\.julia\packages\DiffEqBase\jhLIm\src\callbacks.jl:859 [inlined]
[5] handle_callbacks!(integrator::OrdinaryDiffEq.ODEIntegrator{Vern7, true, Vector{Float64}, Nothing, Float64, Vector{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, Vector{Float64}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Vern7, OrdinaryDiffEq.InterpolationData{ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Vern7Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Vern7Tableau{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Vern7Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Vern7Tableau{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Dict{Any, Any}, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
@ OrdinaryDiffEq C:\Users\Dan\.julia\packages\OrdinaryDiffEq\vxMSM\src\integrators\integrator_utils.jl:261
[6] _loopfooter!(integrator::OrdinaryDiffEq.ODEIntegrator{Vern7, true, Vector{Float64}, Nothing, Float64, Vector{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, Vector{Float64}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Vern7, OrdinaryDiffEq.InterpolationData{ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Vern7Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Vern7Tableau{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Vern7Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Vern7Tableau{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Dict{Any, Any}, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
@ OrdinaryDiffEq C:\Users\Dan\.julia\packages\OrdinaryDiffEq\vxMSM\src\integrators\integrator_utils.jl:204
[7] loopfooter!
@ C:\Users\Dan\.julia\packages\OrdinaryDiffEq\vxMSM\src\integrators\integrator_utils.jl:168 [inlined]
[8] solve!(integrator::OrdinaryDiffEq.ODEIntegrator{Vern7, true, Vector{Float64}, Nothing, Float64, Vector{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, Vector{Float64}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, Vern7, OrdinaryDiffEq.InterpolationData{ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Vern7Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Vern7Tableau{Float64, Float64}}}, DiffEqBase.DEStats}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, OrdinaryDiffEq.Vern7Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Vern7Tableau{Float64, Float64}}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64,
typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryMinHeap{Float64}, DataStructures.BinaryMinHeap{Float64}, Dict{Any, Any}, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit})
@ OrdinaryDiffEq C:\Users\Dan\.julia\packages\OrdinaryDiffEq\vxMSM\src\solve.jl:456
[9] #__solve#403
@ C:\Users\Dan\.julia\packages\OrdinaryDiffEq\vxMSM\src\solve.jl:5 [inlined]
[10] #solve_call#56
@ C:\Users\Dan\.julia\packages\DiffEqBase\jhLIm\src\solve.jl:61 [inlined]
[11] solve_up(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing,
Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, sensealg::Nothing, u0::Vector{Float64}, p::Vector{Float64}, args::Vern7; kwargs::Base.Iterators.Pairs{Symbol, Any, NTuple{4, Symbol}, NamedTuple{(:userdata, :reltol, :abstol, :callback), Tuple{Dict{Any, Any}, Float64, Float64, DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}}}})
@ DiffEqBase C:\Users\Dan\.julia\packages\DiffEqBase\jhLIm\src\solve.jl:82
[12] #solve#57
@ C:\Users\Dan\.julia\packages\DiffEqBase\jhLIm\src\solve.jl:70 [inlined]
[13] #__solve#5
@ c:\Users\Dan\Documents\github\OrbitalTrajectories.jl\src\dynamics\orbital_trajectories.jl:117 [inlined]
[14] #solve#3
@ c:\Users\Dan\Documents\github\OrbitalTrajectories.jl\src\dynamics\orbital_trajectories.jl:101 [inlined]
[15] batch_func(i::Int64, prob::EnsembleProblem{State{CR3BP, SynodicFrame{true}, ODEProblem{StaticArrays.MVector{6, Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, alg::Vern7; kwargs::Base.Iterators.Pairs{Symbol, DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}}}})
@ SciMLBase C:\Users\Dan\.julia\packages\SciMLBase\9EjAY\src\ensemble\basic_ensemble_solve.jl:143
[16] #461
@ C:\Users\Dan\.julia\packages\SciMLBase\9EjAY\src\ensemble\basic_ensemble_solve.jl:213 [inlined]
[17] macro expansion
@ C:\Users\Dan\.julia\packages\SciMLBase\9EjAY\src\ensemble\basic_ensemble_solve.jl:221 [inlined]
[18] (::SciMLBase.var"#447#threadsfor_fun#464"{SciMLBase.var"#461#463"{Base.Iterators.Pairs{Symbol, DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}}}}, EnsembleProblem{State{CR3BP, SynodicFrame{true}, ODEProblem{StaticArrays.MVector{6, Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, Vern7}, Tuple{UnitRange{Int64}}, Vector{Trajectory}, UnitRange{Int64}})(onethread::Bool)
@ SciMLBase .\threadingconstructs.jl:81
[19] (::SciMLBase.var"#447#threadsfor_fun#464"{SciMLBase.var"#461#463"{Base.Iterators.Pairs{Symbol, DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{DiscreteCallback{DiffEqCallbacks.var"#44#49"{Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, DiffEqCallbacks.var"#46#51"{Bool, DiffEqCallbacks.var"#47#52"{Bool}, Float64, DiffEqCallbacks.var"#45#50"{var"#5#6", Float64, Base.RefValue{Int64}, Base.RefValue{Float64}}, Base.RefValue{Int64}, Base.RefValue{Float64}}, typeof(DiffEqBase.FINALIZE_DEFAULT)}}}}, EnsembleProblem{State{CR3BP, SynodicFrame{true}, ODEProblem{StaticArrays.MVector{6, Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ODEFunction{true, CR3BP, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}}, typeof(prob_func), typeof(SciMLBase.DEFAULT_OUTPUT_FUNC), typeof(SciMLBase.DEFAULT_REDUCTION), Nothing}, Vern7}, Tuple{UnitRange{Int64}}, Vector{Trajectory}, UnitRange{Int64}})()
@ SciMLBase .\threadingconstructs.jl:48
Thinking through this a bit more, it does seem like it's more of an issue of not making copies of the callbacks when solving an EnsembleProblem. I've injected a simple deepcopy(callback) in my own __solve() dispatch (which passes down to the underlying __solve(ODEProblem)) and it seems to have solved this issue (NOTE: the above issue is therefore no longer present in OrbitalTrajectories >= 0.1.11).
I wonder if this should be default behaviour, or if there should be some way to trigger it automatically or otherwise warn very clearly that this can happen?
We should solve this by making the cache in init