Problem of convergence for the SH model

General Usage
,
1

Hi everyone,

I’m trying to solving the Swift Hohenberg 2D model as a ODE problem, but it works only on a short time range and reducing the resolution of my grid.
I believe that the problem is not only an optimization problem, but what do you think about? Here’s my code:

using DifferentialEquations
using DiffEqOperators, Setfield, Parameters
using BifurcationKit, LinearAlgebra, Plots, SparseArrays
const BK = BifurcationKit


Nx = 128
Ny = 128
lx = 4pi
ly = 4pi

# we use DiffEqOperators to compute the Laplacian operator
function Laplacian2D(Nx, Ny, lx, ly)
	hx = 2*lx/Nx
	hy = 2*ly/Ny
	D2x = CenteredDifference(2, 2, hx, Nx)
	D2y = CenteredDifference(2, 2, hy, Ny)
	Qx = Neumann0BC(hx)
	Qy = Neumann0BC(hy)
	
	A = kron(sparse(I, Ny, Ny), sparse(D2x * Qx)[1]) + kron(sparse(D2y * Qy)[1], sparse(I, Nx, Nx))
	return A
end
function F_sh(u, p, t)
    @unpack l , L1 = p
	return L1 * u .+ (l .* u)
end
X = -lx .+ 2lx/(Nx) * collect(0:Nx-1)
Y = -ly .+ 2ly/(Ny) * collect(0:Ny-1)


# define parameters for the PDE
Δ = Laplacian2D(Nx, Ny, lx, ly)

L1 = (I + Δ)^2;

#l= -0.1

p = (l= -0.1, L1=L1)
u0=   rand(Nx,Ny) |> vec
tspan = (0.0,100.0)
prob = ODEProblem(F_sh,u0,tspan,p)

 sol =solve(prob, alg_hints=[:stiff],save_everystep=false, save_start=false)

And the error is returned:

[CVODES WARNING] CVode
Internal t = 0.0446594 and h = 3.61507e-021 are such that t + h = t on the next step. The solver will continue anyway.
[CVODES WARNING] CVode
The above warning has been issued mxhnil times and will not be issued again for this problem.
[CVODES ERROR] CVode
At t = 0.0446594 and h = 3.61829e-023, the corrector convergence test failed repeatedly or with |h| = hmin.

Thank you so much for the help, every comment is welcome!

2

Try running at a lower tolerance?

3

No, how it is done?

4

I tried now changing the reltol and abstol values, but the error is still returned… maybe the problem is in the definition of the function?

5

You can try splitODE. See:
https://rveltz.github.io/BifurcationKit.jl/dev/tutorialsPD/#Periodic-orbits-from-the-PD-point-(Shooting)-1
for something quite close to your setting

6

It looks like you’re solving just the linear terms of the Swift-Hohenberg equation and have left out the nonlinear damping term (e.g. -u^3). I’d guess the linear equation is unstable for the parameters you’ve chosen. Could it be that the solution is blowing up, causing the time step to shrink in order to meet an absolute tolerance condition?

First thing I’d do is put back in the nonlinear terms and see if the solution converges to something sensible. From https://rveltz.github.io/BifurcationKit.jl/dev/tutorials2/#Snaking-in-the-2d-Swift-Hohenberg-equation-1, that’d be

function F_sh(u, p)
	@unpack l, ν, L1 = p
	return -L1 * u .+ (l .* u .+ ν .* u.^2 .- u.^3)
end

And maybe while you’re at it change the domain size and discretization to match theirs, just to rule out that difference

Nx = 151
Ny = 100
lx = 4*2pi
ly = 2*2pi/sqrt(3)
3 Likes
7

I saw that if I use only L1= Δ the convergence is respected. So the problem begin when I write L1=Δ^2.
There is another way to write the operator L1 successfully?

Thank you very much.

Francesco

8

Yes, for the moment I’m interested to study SH with only the linear term. Here only the solution u=0 is stable, then starting from random initial condition I expect to find the system at u=0. But maybe I’m doing something wrong?

9

there is a sign missing, it should be

-L1 * u .+ l .* u
2 Likes
10

Nice spot. That would explain blowing up on such a short time scale.

11

I’m so sorry for the very stupid mistake.
Thank you so much for the help, probably I’d have never seen it.

1 Like
12

Sometimes it takes someone else’s eyes! That seems to be the nature of software. And personally, I learn things from reading & participating in these discussions, even if the final answer is a typo.

1 Like