Using EnsembleThreads in DifferentialEquations.jl fails for SDEs with trajectories > nthreads

Specific Domains Modelling & Simulations
,
1

I’m trying to run some stochastic simulations for a research problem using DifferentialEquations.jl. I’ve had no problem setting up my differential equation, and solving it is fine until I try to setup an ensembleproblem, but only if the trajectories option is set higher than the number of threads I have available. If I run with EnsembleSerial() there are no issues, but I would really like to take advantage of multithreading if I can. Here is the code I’m running (Julia V 1.5.3, DifferentialEquations.jl V 6.16.0)

 using DifferentialEquations, Parameters
@with_kw struct pars
    β::Array{Float64,2} = [-0.05 0.15 -0.2
                           -0.01 -0.0267 0.05
                           0.1 -0.1 -0.0148]
    r::Vector{Float64} = [6.0,4.0,2.0]
    σ::Float64 = 0.1
end

function LotkaVolterra!(du,u,p,t)
    @unpack r,β = p
    du .= u.*(r.+β*u)
end

#Noise function
function σ_LV!(du,u,p,t)
    @unpack σ = p
    du .= σ*u
end

u0 = [500.0, 400.0, 275.0]
tspan = (0.0,2.0)
dt = 0.005
prob = SDEProblem(LotkaVolterra!,σ_LV!,u0,tspan,pars(σ=0.1),
                           noise = RealWienerProcess(0.0,zeros(size(u0)),zeros(size(u0))))
ensemble_prob = EnsembleProblem(prob)
sol = solve(ensemble_prob, SRIW1(), saveat=dt, EnsembleThreads(), trajectories=100)

I have 8 threads available. I get an assortment of errors when running this with trajectories > 7 , ranging from the solver aborting in every instance because dt < dtmin (which doesn’t happen for trajectories <=7), to errors about trying to access an array at index 0.

2

Not a sulation but it works without the noise part

3

Do a prob_func where you create a new noise process each time?

4

Yes, defining a prob_func resetting the noise process each time works. Seems like a bit of a work around though.