I’m studying a pde (swift-hohenberg model) converting it in the ODE form.
I’ve discretized the laplacian:
function Laplacian2D(Nx, lx)
hx = lx/Nx
D2x = CenteredDifference(2, 2, hx, Nx)
Qx = PeriodicBC(hx |> typeof)
A = kron(sparse(I, Nx, Nx), sparse(D2x * Qx)[1]) + kron(sparse(D2x * Qx)[1], sparse(I, Nx, Nx))
return A
end
After that I have built the function for the ODEProblem:
function F_sh2!(du,u,p,t)
@unpack l , d, d2 = p
du .= - u .- d*u .- d2*u .+ (l .* u .- u.^3)
return nothing
end
Where l is a float, d is the laplacian and d2 is laplacian^2.
For no large discretizations (64x64 grid) it’s all ok, but if I try to increase the grid, for instance 128x128, the solver returns error: ReadOnlyMemoryError()
I’m wondering if it’s only due to my computational limits or if I’m missing some errors.
Please share the whole script. I cannot necessarily even know what solver was chosen from what you showed. My guess is that it switched to a method for stiff ODEs and a full Jacobian because sparsity was not specified.
#Domain and grid
Nx = 128
lx = 16*pi;
# the Laplacian operator
function Laplacian2D(Nx, lx)
hx = lx/Nx
D2x = CenteredDifference(2, 2, hx, Nx)
Qx = PeriodicBC(hx |> typeof)
A = kron(sparse(I, Nx, Nx), sparse(D2x * Qx)[1]) + kron(sparse(D2x * Qx)[1], sparse(I, Nx, Nx))
return A
end
function F_sh2!(du,u,p,t)
@unpack l , d, d2 = p
du .= - u .- d*u .- d2*u .+ (l .* u .- u.^3)
return nothing
end
# definition of the parameter for the PDE
del = Laplacian2D(Nx, lx)
D=2 .*del
D2=del^2;
p=@time (l=0.1,d=D,d2=D2 );
#initial conditions and tspan
tspan = (0.0,110)
u0= rand(MersenneTwister(1234),Nx*Nx)
#Definition of the problem
prob =@time ODEProblem(F_sh2!,u0,tspan,p);
#Solution of the problem
sol =@time solve(prob,lsoda(),abstol=1e-5,reltol=1e-3,save_everystep=false,save_start=false,dt=0.1)
Yes, LSODA will be a bad idea for this code given it uses a dense Jacobian. That will make it fail with a memory issue if the problem is big enough. I would recommend just not picking an algorithm at all and letting it use the defaults.
Thank you so much! Now it works. I had only to put dt not more fixed because it gave me instabilities (also if it should satisfy the CFL condition). Do you reccomend to not fix the dt? Or maybe I have to do another choice?