`PeriodicCallback` is seemingly not thread-safe in `EnsembleProblem` with `EnsembleThreads`

[dpad]

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

[dpad]

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?

[ChrisRackauckas]

We should solve this by making the cache in init