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. Spatially localized quasicrystalline structures - IOPscience.
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)
