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
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)
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!
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:
So I think, writing your functions inplace to remove allocations is a good idea:
function F_sh(du, u, p, t)
@unpack l , L1 = p
@. du = -L1 * u + (l * u) # if L1 and l are numbers ..
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).
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)) + kron(sparse(D2y * Qy), sparse(I, Nx, Nx))
I tried to write F_sh like
function F_sh(du,u, p, t)
@unpack l , d = p
du .= -(I + Δ)^2 * u .+ (l .* u)
Putting d=Δ. But this not seems to be the good way to solve my problem. Maybe I have to writing different the Laplacian function?
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.
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.
Watch out with broadcasting
I though: it doesn’t mean what you think it means.
Ah I missed that. Good point. So here you would avoid
@. and just pay more attention to where to place dots. Thanks!
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.
Yes, and also write (I + Δ)^2 is different from (I + Δ).^2