Example 2: Colored overlapping smoother for graphene

In here we use the ALFA framework to analyze a 4 color overlap smoother for the tight-binding Hamiltonian of graphene.

It corresponds to section 6.1 of [1] Kahl, K., Kintscher, N. Automated local Fourier analysis (ALFA). Bit Numer Math (2020). https://doi.org/10.1007/s10543-019-00797-w.

Definition of the operators

using ALFA
using LinearAlgebra
using Plots

System operator

L = ALFA.gallery.graphene_tight_binding()
plot(L)

Restriction operator of the two-grid method

R = ALFA.gallery.graphene_dirac_restriction(wl=.25, wlh=.25)
plot(R)

Prolongation and coarse grid operator

P = R'
Lc = R*L*P

p1 = plot(P, title="P",  aspect_ratio=:equal)
p2 = plot(Lc, title="Lc",  aspect_ratio=:equal)
(w,h) = p1.attr[:size]
plot(p1, p2, layout=(2,1), size=(w,2h))

Smoother definition

We first rewrite the system operator $L$ with respect to the translational invariance 2A, where $A=$$L.C.L.A$

slat = ALFA.Lattice{2, Float64}(2*L.C.L.A)
Ls = ALFA.wrtLattice(L,slat)
plot(Ls)

We now construct the $4$ operators used in the four color smoother.

idx=[2,3,4,5,6,7] #
shifts = [[i,j] for i in [0,1] for j in [0,1]]
S = []
p = []
for (it,x) in enumerate(shifts)
   se = [y + L.C.L.A*x for y in Ls.C.Domain]
   Ls_tmp = ALFA.ChangeStructureElement(Ls, se, se)
   push!(S, ALFA.CrystalOperatorCopyWithMultipliers(Ls_tmp, idx=idx))
   push!(p, plot(S[it], title="S[$it]"))
end
plot(p..., layout=(2,2), size=(1.5w,1.5h))

... and normalize the operators S for illustrative purposes:

p = []
for (j,s) in enumerate(S)
   S[j] = ALFA.normalize(s)
   push!(p, plot(S[j], title="S[$j]"))
end
plot(p..., layout=(2,2), size=(1.5w,1.5h))

Spectrum of the error propagator of the Smoother

Now we construct the error propagator of the smoother and plot its spectrum

f_s = prod(:(I-0.5*pinv($s)*$L) for s in S)
oc_s = ALFA.OperatorComposition(f_s)
plotspectrum(oc_s, N=41)

Spectrum of the error propagator of the two-grid method

Finally, we analyze the two-grid error propagator

f_cgc = :(I-$P*inv($Lc)*$R*$L)
f_tg = f_s*f_cgc*f_s
oc_tg = ALFA.OperatorComposition(f_tg)
plotspectrum(oc_tg, N=41)