writing into shared arrays within a distributed for loop in JULIA

Question

I have an ODE that I need to solver over a wide range of parameters. Previously I have used MATLAB's parfor to divide the parameter ranges between multiple threads. I am new to Julia and need to do the same thing in Julia now. Here is the code that I am using

using DifferentialEquations, SharedArrays, Distributed, Plots


function SingleBubble(du,u,p,t)
    du[1]=@. u[2]
    du[2]=@. ((-0.5*u[2]^2)*(3-u[2]/(p[4]))+(1+(1-3*p[7])*u[2]/p[4])*((p[6]-p[5])/p[2]+2*p[1]/(p[2]*p[8]))*(p[8]/u[1])^(3*p[7])-2*p[1]/(p[2]*u[1])-4*p[3]*u[2]/(p[2]*u[1])-(1+u[2]/p[4])*(p[6]-p[5]+p[10]*sin(2*pi*p[9]*t))/p[2]-p[10]*u[1]*cos(2*pi*p[9]*t)*2*pi*p[9]/(p[2]*p[4]))/((1-u[2]/p[4])*u[1]+4*p[3]/(p[2]*p[4]))
end

R0=2e-6
f=2e6
u0=[R0,0]
LN=1000

RS = SharedArray(zeros(LN))
P = SharedArray(zeros(LN))
bif = SharedArray(zeros(LN,6))

 @distributed     for i= 1:LN
    ps=1e3+i*1e3
    tspan=(0,60/f)
    p=[0.0725,998,1e-3,1481,0,1.01e5,7/5,R0,f,ps]
    prob = ODEProblem(SingleBubble,u0,tspan,p)
    sol=solve(prob,Tsit5(),alg_hints=:stiff,saveat=0.01/f,reltol=1e-8,abstol=1e-8)
    RS[i]=maximum((sol[1,5000:6000])/R0)
    P[i]=ps
    for j=1:6
          nn=5500+(j-1)*100;
          bif[i,j]=(sol[1,nn]/R0);
     end
end


plotly()
scatter(P/1e3,bif,shape=:circle,ms=0.5,label="")#,ma=0.6,mc=:black,mz=1,label="")

When using one worker, the for loop is basically executed as a normal single threaded loop and it works fine. However, when I am using addprocs(n) to add n more workers, nothing gets written into the SharedArrays RS, P and bif. I appreciate any guidance anyone may provide.

Answer 1

These changes are required to make your program work with multiple workers and display the results you need:

1. Whatever packages and functions are used under @distributed loop must be made available in all the processes using @everywhere as explained here. So, in your case it would be DifferentialEquations and SharedArrays packages as well as the SingleBubble() function.
2. You need to draw the plot only after all the workers have finished their tasks. For this, you would need to use @sync along with @distributed.

With these changes, your code would look like:

using Distributed, Plots

@everywhere using DifferentialEquations, SharedArrays

@everywhere function SingleBubble(du,u,p,t)
    du[1]=@. u[2]
    du[2]=@. ((-0.5*u[2]^2)*(3-u[2]/(p[4]))+(1+(1-3*p[7])*u[2]/p[4])*((p[6]-p[5])/p[2]+2*p[1]/(p[2]*p[8]))*(p[8]/u[1])^(3*p[7])-2*p[1]/(p[2]*u[1])-4*p[3]*u[2]/(p[2]*u[1])-(1+u[2]/p[4])*(p[6]-p[5]+p[10]*sin(2*pi*p[9]*t))/p[2]-p[10]*u[1]*cos(2*pi*p[9]*t)*2*pi*p[9]/(p[2]*p[4]))/((1-u[2]/p[4])*u[1]+4*p[3]/(p[2]*p[4]))
end

R0=2e-6
f=2e6
u0=[R0,0]
LN=1000

RS = SharedArray(zeros(LN))
P = SharedArray(zeros(LN))
bif = SharedArray(zeros(LN,6))

@sync @distributed     for i= 1:LN
    ps=1e3+i*1e3
    tspan=(0,60/f)
    p=[0.0725,998,1e-3,1481,0,1.01e5,7/5,R0,f,ps]
    prob = ODEProblem(SingleBubble,u0,tspan,p)
    sol=solve(prob,Tsit5(),alg_hints=:stiff,saveat=0.01/f,reltol=1e-8,abstol=1e-8)
    RS[i]=maximum((sol[1,5000:6000])/R0)
    P[i]=ps
    for j=1:6
          nn=5500+(j-1)*100;
          bif[i,j]=(sol[1,nn]/R0);
     end
end


plotly()
scatter(P/1e3,bif,shape=:circle,ms=0.5,label="")#,ma=0.6,mc=:black,mz=1,label="")

Output using multiple workers: