DifferentialEquations.jl multithreading

Question

The code below is solving an ODE with different variables, and PAs corresponding to diffrent Δs the only data I'm interested in after each loop. Since the loops are independent I assume I can multithreading the loops so I add Threads.@threads in front of the for loop:

#example
using DifferentialEquations
using Plots;
gr(grid = false, ms = 3, legend = false, msw = 1.5 ,mc = :white)
using Findpeaks

#parameters
P = 1.5
κ = 0.5
Γ = 0.0005

#ODE settings
fs = 400
T = 2*π
y0 = [1, 1]
t0 = [0, 100 *T]
t1 = [0, 50 *T]
function route(dx, x, p, t)
    dx[1] =  1im * (Δ * x[1] - 2 * real(x[2]) * x[1] - 0.5) - κ * x[1]
    dx[2] = -1im * (0.5 * P * abs(x[1])^2 + x[2]) - Γ * x[2]
end
Plots.scatter(1,1) #create new figure

Threads.@threads for i =  -1.3: 0.1 :-0.4
    global Δ = i
    global p =[P, Δ, κ, Γ]
    global prob = ODEProblem(route, collect(Complex{Float64}, y0), t0, p)
    global y1 = solve(prob)
    global y10 = last(y1);#omit unstable parts

    global prob2 = ODEProblem(route, collect(Complex{Float64}, y10), t1, p)
    global y2 = solve(prob2) #calculate steady state data

    global abss = abs.(y2)
    global a = abss[1,:].^2 

    global PA = begin
            global peaka = findpeaks(real(a), y2.t)
            global peaka = real(a)[peaka] #peak values
            global pa = Vector{Float64}(undef, length(peaka))
            pa[:] .= Δ #peak locations corresponding to Δ
            pa .+ 1im*peaka
    end
    global fig = scatter!(real(PA), imag(PA))
    display(fig)
end

But it will cause julia terminal to crash(using vscode), seems related to GR.

What is happening?

score 2 · Answer 1 · answered Apr 18 '22 at 11:58

2

It looks like things thrown from GR. I'm not sure GR is multithreading safe, I would double check that instead of the differential equation solve.

answered Apr 18 '22 at 11:58

Chris Rackauckas

18,645
3
50
81

I have aviod this problem by moving the plot part out of the loop. Still didn't figure out why that happen though. – BenXylona Apr 19 '22 at 05:55
GR's runtime is not thread safe. – Chris Rackauckas Apr 20 '22 at 15:26
So what backend is thread safe? – BenXylona Apr 21 '22 at 07:10
1

I'm not sure any of the Plots.jl backends are thread safe. That's something that would need to be explored. The issue is generally whether the C library itself is, but many of them probably didn't consider this. I'd check whether Makie.jl is thread safe, that might have a higher chance of working out because it doesn't use a bunch of global information like Plots. – Chris Rackauckas Apr 21 '22 at 11:02

score 0 · Answer 2 · answered Aug 16 '23 at 04:09

c = ReentrantLock()
Threads.@threads for i =  -1.3: 0.1 :-0.4
    ...
    lock(c) do
        global fig = scatter!(real(PA), imag(PA))
        display(fig)
    end
end

Just use a thread safe lock to wrap your plotting code like above. Since plotting is not the CPU sensitive part of your program (compared to ODE equation), it should not affect the performance very much.

I met similar issue with CairoMakie too. I guess plotting uses internal list/queue to track the current drawing panel, data race condition might occur. Let's just use locks to make it happy.

DifferentialEquations.jl multithreading

2 Answers2