Hi everyone

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.

Every help will be welcome, thank you so much!

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.

I’m sorry, here my whole code.

``````#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.

1 Like

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?

That’s just the starting `dt`. I would let the heuristics find a good starting dt.

1 Like

Ok.
Thank you so much for the help!