MethodOfLines.jl icon indicating copy to clipboard operation
MethodOfLines.jl copied to clipboard

Leveraging EnsembleProblem for parameter scans

Open thompsonmj opened this issue 2 years ago • 5 comments

In the MOL docs, it discusses how to Remake with different parameter values, but can this be accelerated with EnsembleProblem interface and use parameter design matrices?

thompsonmj avatar May 25 '22 19:05 thompsonmj

Yes, if you use remake inside of an EnsembleProblem then you can use this all to for example multithread over multiple solves.

ChrisRackauckas avatar May 25 '22 21:05 ChrisRackauckas

Huh, it's working, but the ensemble is much slower ...

Using the tutorial problem,

@parameters t x
@parameters Dn, Dp
@variables u(..) v(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

eqs  = [Dt(u(t,x)) ~ Dn * Dxx(u(t,x)) + u(t,x)*v(t,x), 
        Dt(v(t,x)) ~ Dp * Dxx(v(t,x)) - u(t,x)*v(t,x)]
bcs = [u(0,x) ~ sin(pi*x/2),
       v(0,x) ~ sin(pi*x/2),
       u(t,0) ~ 0.0, Dx(u(t,1)) ~ 0.0,
       v(t,0) ~ 0.0, Dx(v(t,1)) ~ 0.0]

domains = [t ∈ Interval(0.0,1.0),
           x ∈ Interval(0.0,1.0)]

@named pdesys = PDESystem(eqs,bcs,domains,[t,x],[u(t,x),v(t,x)],[Dn=>0.5, Dp=>2])
discretization = MOLFiniteDifference([x=>0.1],t)
prob = discretize(pdesys,discretization) # This gives an ODEProblem since it's time-dependent

N = 10000
lb = [0.01,0.01]
ub = [2.0,2.0]
par = QuasiMonteCarlo.sample(N,lb,ub,LatinHypercubeSample())
sols1 = []
@time begin
for (Dnval, Dpval) in zip(par[1,:], par[2,:])
    newprob = remake(prob, p=[Dnval, Dpval])
    push!(sols1, solve(newprob, Tsit5()))
end
end
# 5.7 s

function prob_func(prob,i,repeat)
    remake(prob,p=par[:,i])
end

ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
@time begin
sols2 = solve(ensemble_prob, Tsit5(), trajectories=N)
end
# 91.0 s

thompsonmj avatar May 25 '22 21:05 thompsonmj

Your function solves too fast for parallelism (multithreading) to be beneficial. We may be able to batch that differently.

ChrisRackauckas avatar May 25 '22 21:05 ChrisRackauckas

Ah, well it would help too if I started Julia with more than 1 thread.

For my real problem, the Ensemble speeds it up nicely when I start it up right. Could be a helpful tutorial to include!

thompsonmj avatar May 26 '22 19:05 thompsonmj

haha yes, always need those threads.

ChrisRackauckas avatar May 26 '22 20:05 ChrisRackauckas