SciML / MethodOfLines.jl

Automatic Finite Difference PDE solving with Julia SciML
https://docs.sciml.ai/MethodOfLines/stable/
MIT License
162 stars 31 forks source link

Leveraging EnsembleProblem for parameter scans #104

Open thompsonmj opened 2 years ago

thompsonmj commented 2 years ago

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?

ChrisRackauckas commented 2 years ago

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

thompsonmj commented 2 years ago

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
ChrisRackauckas commented 2 years ago

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

thompsonmj commented 2 years ago

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!

ChrisRackauckas commented 2 years ago

haha yes, always need those threads.