DifferentialEquations.jl multithreading

418 views Asked by At

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.

enter image description here

What is happening?

2

There are 2 answers

4
Chris Rackauckas On

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.

0
Leon Chen On
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.