Improving performance in solving SDEProblem

General Usage Performance
,
1

Hi everyone,
I’m studying part of the 2D Swift Hohenberg model with an additive noise term, but I’m really new to the SDEProblems so I don’t know if my way to build the model is the best way to solve it.
I’ve take a look to the tutorial Stochastic Differential Equations · DifferentialEquations.jl and I’ve defined

function Add_noise(u,p,t)
     
    return fill(0.5,Nx*Ny)
end

As the noise. Then

#Definition of the problem
u0=   rand(Nx,Ny) |> vec
tspan = (0.0,35.0)
prob =@time SDEProblem(F_sh,Add_noise,u0,tspan,p);

Where F_sh is

#Definition of the equation
function F_sh(u, p, t)
    @unpack l , L1 = p
    
    return -L1 * u .+ (l .* u)
end

And L1 is the SH operator.
Then I try to solve the problem using

#Solution of the problem

sol =@time solve(prob,alg_hints=[:additive,:stiff],SRIW1(),  save_everystep=false, save_start=false, saveat=[5,35])

But it seems not working very well because is slow, in particular increasing the resolution of the grid (Nx*Ny).
What kind of error I’m doing?
Thank you so much for the help, every comment is welcome!

2

Allocating large vectors as the output is probably the biggest performance hit (and for vectors of small size you should use StaticArrays). There is a nice tutorial for optimizing DiffEq code:

https://tutorials.sciml.ai/html/introduction/03-optimizing_diffeq_code.html

So I think, writing your functions inplace to remove allocations is a good idea:

function Add_noise(du,u,p,t)
  fill!(du,1//2*one(eltype(u)))
  return nothing
end

function F_sh(du, u, p, t)
    @unpack l , L1 = p

    @. du = -L1 * u + (l * u) #  if L1 and l are numbers ..
    return nothing
end

Since your noise is additive, I think using a specialized solver will result in a speed up. Likely, SOSRA() is a good choice (you’ll find the available solvers at SDE Solvers · DifferentialEquations.jl).

1 Like
3

Thank you so much Frank, but still something goes wrong. I think that the problem is because L1 is not a number but L1=(I+Δ)^2 where Δ is the laplacian operator that I wrote as

#Domain and grid 
Nx = 90
Ny = 90
lx = 30
ly = 30

# the Laplacian operator
function Laplacian2D(Nx, Ny, lx, ly)
	hx = lx/Nx
	hy = ly/Ny
	D2x = CenteredDifference(2, 2, hx, Nx)
	D2y = CenteredDifference(2, 2, hy, Ny)
	Qx = PeriodicBC(hx |> typeof)
	Qy = PeriodicBC(hy |> typeof)
	
	A = kron(sparse(I, Ny, Ny), sparse(D2x * Qx)[1]) + kron(sparse(D2y * Qy)[1], sparse(I, Nx, Nx))
	return A
end

I tried to write F_sh like

function F_sh(du,u, p, t)
    @unpack l , d = p
    
    du .= -(I + Δ)^2 * u .+ (l .* u)
    return nothing
end

Putting d=Δ. But this not seems to be the good way to solve my problem. Maybe I have to writing different the Laplacian function?

4

I see. Δ seems to be defined in global scope now and you don’t want to compute the square in each function call. So, you’ll have to localize L1 within the function by either adding it to the parameter vector (Why did you change that?) or by using a functor.

1 Like
5

Just a remark in passing: I̶f̶ ̶y̶o̶u̶ ̶b̶r̶o̶a̶d̶c̶a̶s̶t̶ ̶w̶i̶t̶h̶ ̶̶@̶.̶̶ ̶a̶t̶ ̶t̶h̶e̶ ̶s̶t̶a̶r̶t̶ ̶o̶f̶ ̶t̶h̶e̶ ̶l̶i̶n̶e̶,̶ ̶y̶o̶u̶ ̶d̶o̶n̶'̶t̶ ̶h̶a̶v̶e̶ ̶t̶o̶ ̶w̶o̶r̶r̶y̶ ̶a̶b̶o̶u̶t̶ ̶w̶h̶e̶r̶e̶ ̶e̶x̶a̶c̶t̶l̶y̶ ̶t̶o̶ ̶p̶l̶a̶c̶e̶ ̶d̶o̶t̶s̶:

@̶.̶ ̶d̶u̶ ̶=̶ ̶-̶(̶I̶ ̶+̶ ̶Δ̶)̶^̶2̶ ̶*̶ ̶u̶ ̶+̶ ̶(̶l̶ ̶*̶ ̶u̶)̶

No, don’t do that! See Chris’s remark below.

6

Watch out with broadcasting I though: it doesn’t mean what you think it means.

7

Ah I missed that. Good point. So here you would avoid @. and just pay more attention to where to place dots. Thanks!

8

Ok thank you, using your hints i reach the solution in

619.053820 seconds (7.22 M allocations: 1.374 TiB, 6.78% gc time)

Then, it’s quite good now.

9

Yes, and also write (I + Δ)^2 is different from (I + Δ).^2

1 Like