🟠 3d Swift-Hohenberg equation, Finite differences

• 🟠 3d Swift-Hohenberg equation, Finite differences
• • Plotting with Makie
  • Setting up the problem
  • Choice of linear solver
  • Continuation and bifurcation points
  • Automatic branch switching

We look at the following PDE on a 3d domain, e.g. a cube:

\[-(I+\Delta)^2 u+l\cdot u +\nu u^2-u^3 = 0\tag{E}\]

with Neumann boundary conditions. We use a Sparse Matrix to express the operator $L_1\equiv(I+\Delta)^2$. However, compared to the 2d case (see 2d Swift-Hohenberg equation: snaking), we cannot use directly \ to solve linear systems because the LU décomposition is a bit slow, it uses a lot of memory.

We start by defining the associated functional to encode (E).

using Revise, KrylovKit
using GLMakie # must be imported before BifurcationKit to trigger some imports
using BifurcationKit
using LinearAlgebra, SparseArrays, LinearMaps
const BK = BifurcationKit

function Laplacian3D(Nx, Ny, Nz, lx, ly, lz)
    speye(n) = sparse(I, n, n)
    hx = 2lx/Nx
    hy = 2ly/Ny
    hz = 2lz/Nz
    D2x = spdiagm(0 => -2ones(Nx), 1 => ones(Nx-1), -1 => ones(Nx-1) ) / hx^2
    D2y = spdiagm(0 => -2ones(Ny), 1 => ones(Ny-1), -1 => ones(Ny-1) ) / hy^2
    D2z = spdiagm(0 => -2ones(Nz), 1 => ones(Nz-1), -1 => ones(Nz-1) ) / hz^2

    D2x[1,1] = -1/hx^2
    D2x[end,end] = -1/hx^2

    D2y[1,1] = -1/hy^2
    D2y[end,end] = -1/hy^2

    D2z[1,1] = -1/hz^2
    D2z[end,end] = -1/hz^2

    D2xsp = sparse(D2x)
    D2ysp = sparse(D2y)
    D2zsp = sparse(D2z)

    _A = kron(speye(Ny), D2xsp) + kron(D2ysp, speye(Nx))
    A = kron(speye(Nz), _A) + kron(kron(D2zsp, speye(Ny)), speye(Nx))
    return A, D2x
end

# main functional
function F_sh(u, p)
	(;l, ν, L1) = p
	return -(L1 * u) .+ (l .* u .+ ν .* u.^2 .- u.^3)
end

# differential of the functional
function dF_sh(u, p, du)
	(;l, ν, L1) = p
	return -(L1 * du) .+ (l .+ 2 .* ν .* u .- 3 .* u.^2) .* du
end

# various differentials
d2F_sh(u, p, dx1, dx2) = (2 .* p.ν .* dx2 .- 6 .* dx2 .* u) .* dx1
d3F_sh(u, p, dx1, dx2, dx3) = (-6 .* dx2 .* dx3) .* dx1

# these types are useful to switch to GPU
TY = Float64
AF = Array{TY}

Plotting with Makie

In most tutorials, we have used Plots.jl for the figures. However, it appears that Makie.jl is more convenient for 3d plots. We thus define the following convenience functions to display the solutions of (E).

contour3dMakie(x; k...) = GLMakie.contour(x;  k...)
contour3dMakie(x::AbstractVector; k...) = contour3dMakie(reshape(x,Nx,Ny,Nz); k...)

contour3dMakie(ax, x; k...) = (GLMakie.contour!(ax, x;  k...))
contour3dMakie(ax, x::AbstractVector; k...) = contour3dMakie(ax, reshape(x,Nx,Ny,Nz); k...)

contour3dMakie!(ax, x; k...) = (GLMakie.contour!(ax, x;  k...))
contour3dMakie!(ax, x::AbstractVector; k...) = contour3dMakie!(ax, reshape(x,Nx,Ny,Nz); k...)

Setting up the problem

We provide the parameters defining the PDE:

Nx = Ny = Nz = 22; N = Nx*Ny*Nz
lx = ly = lz = pi

X = -lx .+ 2lx/(Nx) * collect(0:Nx-1)
Y = -ly .+ 2ly/(Ny) * collect(0:Ny-1)
Z = -lz .+ 2lz/(Nz) * collect(0:Nz-1)

# initial guess for newton
sol0 = [(cos(x) .* cos(y )) for x in X, y in Y, z in Z]
sol0 .= sol0 .- minimum(vec(sol0))
sol0 ./= maximum(vec(sol0))
sol0 .*= 1.7

# parameters for PDE
Δ, D2x = Laplacian3D(Nx, Ny, Nz, lx, ly, lz);
L1 = (I + Δ)^2;
par = (l = 0.1, ν = 1.2, L1 = L1);

Choice of linear solver

Let us run a quick benchmark to evaluate the direct linear solvers:

julia> @time cholesky(L1) \ sol_hexa;
  0.152849 seconds (54 allocations: 87.273 MiB)

julia> @time lu(L1) \ sol_hexa;
  0.556157 seconds (87 allocations: 226.210 MiB, 0.49% compilation time)

julia> @time qr(L1) \ sol_hexa;
  1.609175 seconds (8.96 k allocations: 989.285 MiB, 2.67% gc time, 0.67% compilation time)

Hence, cholesky is the big winner but it requires a positive matrix so let's see how to do that.

As said in the introduction, the LU linear solver does not scale well with dimension $N$. Hence, we do something else. We note that the matrix $L_1$ is hermitian positive and use it as a preconditioner. Thus, we pre-factorize it using a Cholesky decomposition:

Pr = cholesky(L1);
using SuiteSparse
# we need this "hack" to be able to use Pr as a preconditioner.
LinearAlgebra.ldiv!(o::Vector, P::SuiteSparse.CHOLMOD.Factor{Float64}, v::Vector) = o .= (P \ v)

# rtol must be small enough to pass the Fold points and to get precise eigenvalues
# we know that the jacobian is symmetric so we tell the solver
ls = GMRESKrylovKit(verbose = 0, rtol = 1e-9, maxiter = 150, Pl = Pr)

Let's try this on a Krylov-Newton computation to refine the guess sol0:

prob = BifurcationProblem(F_sh, AF(vec(sol0)), par, (@optic _.l),
	J = (x, p) -> (dx -> dF_sh(x, p, dx)),
	plot_solution = (ax, x, p; ax1) -> contour3dMakie(ax, x),
	issymmetric = true, # so we dont need to provide the adjoint (normal form)
	record_from_solution = (x, p; k...) -> (n2 = norm(x), n8 = norm(x, 8)))

optnew = NewtonPar(verbose = true, tol = 1e-8, max_iterations = 20, linsolver = ls)
sol_hexa = @time BK.solve(prob, Newton(), optnew)

which gives

┌─────────────────────────────────────────────────────┐
│ Newton step         residual     linear iterations  │
├─────────────┬──────────────────────┬────────────────┤
│       0     │       2.6003e+02     │        0       │
│       1     │       1.5414e+02     │       25       │
│       2     │       2.6040e+02     │       21       │
│       3     │       7.3531e+01     │       21       │
│       4     │       2.0513e+01     │       23       │
│       5     │       6.4608e+00     │       18       │
│       6     │       1.3743e+00     │       18       │
│       7     │       1.7447e-01     │       17       │
│       8     │       4.0924e-03     │       17       │
│       9     │       2.4048e-06     │       17       │
│      10     │       1.8389e-10     │       17       │
└─────────────┴──────────────────────┴────────────────┘
  0.472765 seconds (13.97 k allocations: 134.073 MiB, 1.91% gc time, 0.00% compilation time)

and contour3dMakie(sol_hexa.u) produces

Continuation and bifurcation points

We now switch gears and compute the stability of the trivial solution $u=0$. We will then branch from the detected bifurcation points. However, we wish to show an example of computation of eigenvalues based on a custom preconditioned Shift-Invert strategy.

We thus define our eigensolver based on the previously defined pre-conditioned linear solver ls:

# structure to hold eigensolver
struct SH3dEig{Ts, Tσ} <: BK.AbstractEigenSolver
	# linear solver used for Shift-Invert strategy
	ls::Ts
	# shift of the linear operator
	σ::Tσ
end

# function to extract eigenvectors, used for automatic branch switching
BifurcationKit.geteigenvector(eigsolve::SH3dEig, vecs, n::Union{Int, Array{Int64,1}}) = vecs[n]

# implementation of Shift-invert strategy
function (sheig::SH3dEig)(J, nev::Int; verbosity = 0, kwargs...)
	σ = sheig.σ
	nv = 30
	Jshift = du -> J(du) .- σ .* du
	A = du -> sheig.ls(Jshift, du)[1]
	# we adapt the krylov dimension as function of the requested eigenvalue number
	vals, vec, info = KrylovKit.eigsolve(A, AF(rand(Nx*Ny*Nz)), nev, :LM;
		 tol = 1e-12, maxiter = 20, verbosity = verbosity, ishermitian = true,
		 krylovdim = max(nv, nev + nv))
	vals2 = 1 ./vals .+ σ
	Ind = sortperm(vals2, by = real, rev = true)
	return vals2[Ind], vec[Ind], true, info.numops
end

We can then declare our eigensolver and pass it to the newton parameters

eigSH3d = SH3dEig((@set ls.rtol = 1e-9), 0.1)
@reset optnew.eigsolver = eigSH3d

We are now ready to perform continuation and detection of bifurcation points:

optcont = ContinuationPar(dsmin = 0.0001, dsmax = 0.005, ds= -0.001, p_max = 0.15,
	p_min = -.1, newton_options = setproperties(optnew; tol = 1e-9, max_iterations = 15),
	max_steps = 146, detect_bifurcation = 3, nev = 15, n_inversion = 4, plot_every_step = 1)

br = continuation( re_make(prob, u0 = zeros(N)),
  # we use a particular bordered linear solver to
  # take advantage of our specific linear solver
  PALC(bls = BorderingBLS(solver = optnew.linsolver, check_precision = false)),
  optcont;
  normC = x -> norm(x, Inf),
  plot = true, verbosity = 3)

The following result shows the detected bifurcation points (its takes ~300s)

julia> br
Branch number of points: 34
Branch of Equilibrium
Parameters l from 0.1 to -0.1
Bifurcation points:
 (ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`)
- #  1,    bp at l ≈ +0.06243495 ∈ (+0.06243495, +0.06287689), |δp|=4e-04, [converged], δ = (-1,  0), step =   8, eigenelements in eig[  9], ind_ev =  10
- #  2,    nd at l ≈ +0.06069653 ∈ (+0.06069653, +0.06069826), |δp|=2e-06, [converged], δ = (-6,  0), step =   9, eigenelements in eig[ 10], ind_ev =   9
- #  3,    nd at l ≈ -0.00007046 ∈ (-0.00007046, +0.00015051), |δp|=2e-04, [converged], δ = (-3,  0), step =  18, eigenelements in eig[ 19], ind_ev =   3

We get the following plot during computation:

Automatic branch switching

We can use Branch switching to compute the different branches emanating from the bifurcation points. For example, the following code will perform automatic branch switching from the last bifurcation point of br. Note that this bifurcation point is 3d.

br1 = @time continuation(br, 3, setproperties(optcont; save_sol_every_step = 10,
	detect_bifurcation = 0, p_max = 0.1, plot_every_step = 5, dsmax = 0.02);
	plot = true, verbosity = 3,
	# to set initial point on the branch
	δp = 0.01,
	# remove display of deflated newton iterations
	verbosedeflation = false,
	alg = PALC(tangent = Bordered(), bls = BorderingBLS(solver = optnew.linsolver, check_precision = false)),
	normC = norminf)

We can then plot the branches using BK.plot(br, br1...) where green (resp. red) means stable (resp. unstable) solution.

There are 19 banches that were discovered. You can plot the solutions on the branches using

fig = Figure(resolution = (1200, 900))
	for i=1:length(br1)
		ix = div(i,5)+1; iy = i%5+1
		@show i, ix, iy
		ax = Axis3(fig[ix, iy], title = "$i", aspect = (1, 1, 1))
		hidedecorations!(ax, grid=false)
		contour3dMakie!(ax, br1[i].sol[2].x)
		ax.protrusions = (0, 0, 0, 10)
	end
	display(fig)