# PDE simulation with ETD2. For computing quasicrystal

Dear all,

I am trying to solve a pde using the time stepper ETD2 in Fourier space. My goal was to play with the results of the paper

Subramanian, P., A. J. Archer, E. Knobloch, and A. M. Rucklidge. “Spatially Localized Quasicrystals.” New Journal of Physics 20, no. 12 (December 14, 2018): 122002. https://doi.org/10.1088/1367-2630/aaf3bd.

The goal is to solve \frac{d}{dt}u=-\Delta(\mathcal L u+Qu^2-u^3) where \mathcal L is specified by its symbol (or Fourier transform).

However, I am facing the issue that the Fourier transform v of the state diverges… This can be seen using ENV["JULIA_DEBUG"] = Main.

If anyone has an idea, that would be helpful…My guess is that I messed up with the wave vectors.

PS: The following is not a MWE but should be a WE which produces nice (but wrong) figures like


using Revise
using LinearAlgebra, Plots, Setfield, Parameters

TY = Float64
AF = Array{TY}

#################################################################
# to simplify plotting of the solution
plotsol(x; k...) = heatmap(X,Y,reshape(Array(x), Nx, Ny)'; color=:viridis, k...)
plotsol!(x; k...) = heatmap!(X,Y,reshape(Array(x), Nx, Ny)'; color=:viridis, k...)
norminf(x) = maximum(abs,x)

const Nx = 2^8
const Ny = 2^8
lx = 30pi
ly = lx

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

sol0 = [(cos(x) .+ cos(x/2) * cos(sqrt(3) * y/2) ) for x in X, y in Y]

using AbstractFFTs, FFTW
import Base: *, \

@with_kw struct QPLinearOp{Treal, Tcomp, Tk, Tplan, Tiplan}
tmp_real::Treal         # temporary
tmp_complex::Tcomp      # temporary
k1::Tk
k2::Tk
Δ::Treal
Σmul::Treal
Σdiv::Treal
EXP::Treal
#ETD1
ETD0::Treal
#ETD2
ETD1::Treal

fftplan::Tplan
ifftplan::Tiplan
end

A(k,q) = (k^2*(q^2-3)-2q^2+4)*(q^2-k^2)^2*q^4
B(k,q) = (k^2*(3q^2-1)+2q^2-4q^4)*(1-k^2)^2
Σ(k, par) = k^2 * (par.r1 * A(k, par.q) + par.rq * B(k, par.q))/(par.q^4 * (1-par.q^2)^3) +
(par.σ0 * 1)/(par.q^4) * (1-k^2)^2 * (par.q^2-k^2)^2

function QPLinearOp(N, l, par, AF = Array{TY}) where TY
Nx, Ny = N
lx, ly = l
# define the wave vectors
k1 = vcat(collect(0:Nx/2), collect(Nx/2+1:Nx-1) .- Nx)
k2 = vcat(collect(0:Ny/2), collect(Ny/2+1:Ny-1) .- Ny)
# symbol for L (pde)
Σmul = [ Σ( sqrt((pi/lx * kx)^2 + (pi/ly * ky)^2), par) for kx in k1, ky in k2]
# symbol used to invert L (pde)
Σdiv = 1 .+ Σmul
# symbol for laplacian
Δ = [ -(pi/lx * kx)^2 - (pi/ly * ky)^2 for kx in k1, ky in k2]
# symbol for L (time dep pde)
L = -Δ .* Σmul
# factors for ETD algo
EXP = @. exp(L * par.dt)
# ETD0
ETD0 = @. (exp(L * par.dt) - 1) / L
replace!(x->isnan(x) ? par.dt : x, ETD0) # because, ETD0->dt when L->0
# ETD1
ETD1 = @. (exp(L * par.dt) - 1 - L * par.dt) / (L^2 * par.dt)
replace!(x->isnan(x) ? par.dt/2 : x, ETD1) # because, ETD1->dt/2 when L->0
# temp for complex valued data
tmpc = Complex.(AF(zeros(Nx, Ny)))

return QPLinearOp( tmp_real = AF(zeros(Nx, Ny)), tmp_complex = tmpc,
k1 = AF(k1), k2 = AF(k2), Δ = AF(Δ), Σmul = AF(Σmul), Σdiv = AF(Σdiv), EXP = AF(EXP), ETD0 = AF(ETD0), ETD1 = AF(ETD1), fftplan = plan_fft!(tmpc), ifftplan = plan_ifft!(tmpc))
end

function apply(c::QPLinearOp, u, multiplier, op = *)
c.tmp_complex .= Complex.(u)
c.fftplan * c.tmp_complex
c.tmp_complex .= op.(c.tmp_complex, multiplier)
c.ifftplan * c.tmp_complex
c.tmp_real .= real.(c.tmp_complex)
return copy(c.tmp_real)
end

*(c::QPLinearOp, u) = apply(c, u, c.Σmul)
\(c::QPLinearOp, u) = apply(c, u, c.Σdiv, /)

LinearAlgebra.ldiv!(c::QPLinearOp, u) = c \ u

parL = (r1 = 0.5, rq = 1.0, q = 0.5176, σ0 = -10., dt = 0.01)
L = QPLinearOp((Nx, Ny), (lx, ly), parL, AF)
par = (Q = 2., L = L, parL...)

####################################################################################################
# ETD2 method
amplitude(u) = maximum(u) - minimum(u)

function etdQH(u0, par, nmax, ETD2 = false; nplot = 20, rtol = 1e-5)
println("#"^50)
@unpack Q, dt, L, r1, rq = par
v = fft(u0)
N = length(v)
vtmp = copy(v)
u = copy(u0)
nl = copy(v)	#NL term at current time
nlm1 = copy(nl) #NL term at previous time step
for n = 1:nmax
u .= real.(ifft(v))
u .-= sum(u) / N
nl .= fft(Q.*u.^2 .- u.^3) .* (-L.Δ)
vtmp .= @. L.EXP * v + L.ETD0 * nl
if ETD2 #if ETD2, add the missing part
vtmp .+= L.ETD1 .* (nl .- nlm1)
end
if mod(n, nplot) == 0
@debug ETD2 n, norminf(u), norminf(nl), norminf(vtmp)
@debug amplitude(u) norminf(u-real.(ifft(vtmp)))/norminf(u)
plot(plotsol(u),
plotsol(abs.(vtmp) .|> log, title = "t = $(n*dt), r = ($r1, \$rq)")
) |> display
end
if norminf(u-real.(ifft(vtmp)))/norminf(u) < rtol
@show norminf(u-ifft(vtmp))/norminf(u)
u .= real.(ifft(vtmp))
break
end
v .= vtmp
nlm1 .= nl
end
return u
end

# q-hex
@set! par.r1 = -0.6181
@set! par.rq = -0.412
@set par.dt = 0.01
L = QPLinearOp((Nx, Ny), (lx, ly), par, AF)
@set par.L = L
# initial guess
u0 = L \ AF(0.1randn(Nx,Ny))
mask = (X.^2 .+ Y'.^2) .<= 25^2
u0 .= u0 .*  mask # enforce periodicity of the initial guess
u0 .-= sum(u0) / (Nx*Ny)
res = etdQH(u0, par, 30_000, true; nplot = 100, rtol = 1e-7)
# plotsol(res)


Hi! Oof, that looks like a complicated thing you’re doing, you should definitely try to reduce this further and debug simpler things first. Eg start from the heat equation with implicit euler. The line defining k1 and k2 looks reasonable to me, at least for N even. I have the formula ω(i, N) = (i <= div(N,2)+1) ? (i-1) : -(N-i+1) lying around an old code somewhere. The behavior you have is consistent with something going wrong at boundaries, but your definition of X and Y looks fine. Could that be related to your non-periodic initial condition? I would try with a periodic initial condition first.

(also you should probably write the code in a more modular way, decoupling the integrator and the PDE, so you can more easily swap integrators, and eg plug in DiffEq)

Although I cannot claim to fully grasp everything that is going on either in the code or in the paper, I noticed two things:

1. Is there an ! missing in the line @set par.L = L ? It seems to me that otherwise this line will just return a struct with the new L, but not actually change par (which already contains an L that was defined earlier, with different values for r1 and rq).
Similarly for @set par.dt = 0.01, although that is the same value as defined earlier.
2. In the definition of Σ, there seems to be an overall factor of k^2 missing with respect to formula in the paper, but maybe that gets taken into account elsewhere and I couldn’t quickly see it.
1 Like

Thank you a lot for your suggestions!! I put (code updated) a mask to enforce periodic boundary conditions but it did not change the outcome.

Oh YES!! But no change to the end result

Nice dig up!! but they say " The operator L is obtained from Eq. (3) by first dividing by k^2 and then replacing k^2 by −∇ ." in Three-Dimensional Icosahedral Phase Field Quasicrystal 2016. But I agree, that’s an awful way of defining an operator in a paper (you need to read two others for that…).

I would say (not very efficient but…):

using DifferentialEquations, DiffEqOperators
lin(du,u,p,t) = du.=-apply(p.L, u, p.L.Δ .* p.L.Σmul, *)
NL(du,u,p,t) = du.=-apply(p.L, p.Q.*u.^2 .- u.^3, p.L.Δ, *)
prob = SplitODEProblem(lin, NL, u0, (0.0, 10.0), par)
sol = solve(prob, SplitEuler(), dt = 0.1)


Maybe @ChrisRackauckas has a better way.
Also, I coded my own ETD2 because I have an expression for the exponential of the linear part exp(-\Delta\mathcal L t)u. I dont know how to do that with exp(lin) somehow.

I think DiffEq supports this, but then again if you just want the basic 2nd order scheme there’s probably little to gain by using DiffEq (since there’s no timestep adaptivity for these split schemes)

1 Like

Yeah for these the only complex machinery is the connection to ExponentialUtilities.jl, which is you have an analytical expression for the expmv then you might as well use that. It will specialize on some matrix types already though: worth testing how close it naively gets without tweaking settings.