Massively Parallel ODE Solving with Event Handling and Callbacks

using DiffEqGPU, StaticArrays, OrdinaryDiffEq, CUDA
function f(u, p, t)
    du1 = -u[1]
    return SVector{1}(du1)
end

u0 = @SVector [10.0f0]
prob = ODEProblem{false}(f, u0, (0.0f0, 10.0f0))
prob_func = (prob, i, repeat) -> remake(prob, p = prob.p)
monteprob = EnsembleProblem(prob, safetycopy = false)

condition(u, t, integrator) = t == 4.0f0
affect!(integrator) = integrator.u += @SVector[10.0f0]

gpu_cb = DiscreteCallback(condition, affect!; save_positions = (false, false))

sol = solve(monteprob, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()),
    trajectories = 10,
    adaptive = false, dt = 0.01f0, callback = gpu_cb, merge_callbacks = true,
    tstops = [4.0f0])

EnsembleSolution Solution of length 10 with uType:
SciMLBase.ODESolution{Float32, 2, SubArray{StaticArraysCore.SVector{1, Float32}, 1, Matrix{StaticArraysCore.SVector{1, Float32}}, Tuple{UnitRange{Int64}, Int64}, true}, Nothing, Nothing, SubArray{Float32, 1, Matrix{Float32}, Tuple{UnitRange{Int64}, Int64}, true}, Nothing, DiffEqGPU.ImmutableODEProblem{StaticArraysCore.SVector{1, Float32}, Tuple{Float32, Float32}, false, SciMLBase.NullParameters, SciMLBase.ODEFunction{false, SciMLBase.AutoSpecialize, typeof(Main.f), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, GPUTsit5, SciMLBase.LinearInterpolation{SubArray{Float32, 1, Matrix{Float32}, Tuple{UnitRange{Int64}, Int64}, true}, SubArray{StaticArraysCore.SVector{1, Float32}, 1, Matrix{StaticArraysCore.SVector{1, Float32}}, Tuple{UnitRange{Int64}, Int64}, true}}, Nothing, Nothing}