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.

Thank you for your help.

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!! :heart_eyes:

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.