Example 3: Half-hybrid smoother for the curl-curl equation
In here we use this framework to analyze the half hybrid smoother for the curl-curl equation as described in [2] https://epubs.siam.org/doi/abs/10.1137/S0036142997326203, and analyzed in [3] https://epubs.siam.org/doi/abs/10.1137/070679119.
It corresponds to section 6.2 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
using SparseArrays
System & restriction operator
L = ALFA.gallery.curlcurl(.01)
p1 = plot(L, title="L")
R = ALFA.gallery.curlcurl_restriction()
p2 = plot(R, title="R")
P = R'
Lc = R*L*P
(w,h) = p1.attr[:size]
plot(p1, p2, layout=(2,1), size=(w,2h))
Hybrid smoother
The Discrete gradient operator
A = [1 0; 0 1]
Domain = [[.5, 0],[0, .5]]
Codomain = [[0,0]]
C = ALFA.Crystal{2,Float64}(A, Domain, Codomain)
Rs = ALFA.CrystalOperator{2,Float64}(C)
push!(Rs, ALFA.Multiplier([0 0], [-1 -1]))
push!(Rs, ALFA.Multiplier([-1 0], [1 0]))
push!(Rs, ALFA.Multiplier([0 -1], [0 1]))
Ps = Rs'
Ls = Rs*L*Ps
p1 = plot(Rs, title="Rs")
p2 = plot(Ps, title="Ps")
p3 = plot(Ls, title="Ls")
plot(p1, p2, p3, layout=(3,1), size=(w,3h))
Gauss-Seidel on vertices
# We have to change to lexicographic ordering in order to obtain the same operator as described in [3]
S1 = ALFA.CrystalOperatorCopyLowerTriangle(Ls,perm=[2,1])
plot(S1)
scalar Gauss-Seidel on edges
# change the structure element in order to obtain the same operator as described in [3]
s0_lex = [[.5,0],[1,-.5]]
L_lex = ALFA.ChangeStructureElement(L, s0_lex, s0_lex)
S0 = ALFA.CrystalOperatorCopyLowerTriangle(L_lex, perm=[2,1])
#Now, S0 describes a block-Gauss-Seidel smoother. We need to use the lower triangle of the central multiplier.
m = ALFA.find_multiplier(S0, [0,0])
m.mat = tril(m.mat)
plot(S0)
Spectrum of the smoother
f_edge = :(I-pinv($S0)*$L)
f_node = :(I-$Ps*pinv($S1)*$Rs*$L)
oc_s = ALFA.OperatorComposition(f_edge*f_node)
plotspectrum(oc_s, N=30)
Spectrum of the twogrid method
f_cgc = :(I-$P*inv($Lc)*$R*$L)
oc_tg = ALFA.OperatorComposition(f_node*f_edge*f_cgc*f_node*f_edge)
plotspectrum(oc_tg, N=30)
Double check result with a twogrid implementation
wrtL = ALFA.Lattice{2, Float64}(30*[1.0 0; 0 1.0])
Lm = SparseMatrixCSC{Float64,Int}(ALFA.construct_matrix(L,wrtL))
Rm = SparseMatrixCSC{Float64,Int}(ALFA.construct_matrix(R,wrtL))
Pm = SparseMatrixCSC{Float64,Int}(ALFA.construct_matrix(P,wrtL))
Lcm = SparseMatrixCSC{Float64,Int}(ALFA.construct_matrix(Lc,wrtL))
S0m = SparseMatrixCSC{Float64,Int}(ALFA.construct_matrix(S0,wrtL))
S1m = SparseMatrixCSC{Float64,Int}(ALFA.construct_matrix(S1,wrtL))
Rsm = SparseMatrixCSC{Float64,Int}(ALFA.construct_matrix(Rs,wrtL))
Psm = SparseMatrixCSC{Float64,Int}(ALFA.construct_matrix(Ps,wrtL))
Definition of the application of the smoother and coarse grid correction
S0mLU = lu(S0m)
S1mLU = lu(S1m)
LcmLU = lu(Lcm)
function smooth_edge(b, x)
r = b-Lm*x
x += (S0mLU\r)
return x
end
function smooth_node(b, x)
r = b-Lm*x
rc = Rsm*r
xc = S1mLU\rc
x += Psm*xc
return x
end
function cgc(b, x)
r = b-Lm*x
rc = Rm*r
xc = LcmLU\rc
x += Pm*xc
return x
end
Testrun of the twogrid method.
# initialize rhs and initial guess
global x, b, casym_vec, resnorm_vec, num_iter
n = size(Lm,1);
x = Lm*rand(n);
x = x/norm(x)
b = 0*Lm*rand(n);
casym_vec = [1.0]
resnorm_vec = [1.0]
num_iter = 0
while resnorm_vec[end] > 1e-200 && num_iter < 175
global x, b, casym_vec, resnorm_vec, num_iter
num_iter = num_iter + 1
x = smooth_edge(b,x)
x = smooth_node(b,x)
x = cgc(b,x)
x = smooth_edge(b,x)
x = smooth_node(b,x)
push!(resnorm_vec, norm(b - Lm*x))
push!(casym_vec, resnorm_vec[end]/resnorm_vec[end-1])
end
# plot convergence behavior.
plot(resnorm_vec, yaxis=:log, xlabel="iteration", ylabel="residual norm", title="measured asymptotic convrate="*string(casym_vec[end])*"\n conv.rate from analysis: "*string(abs(ALFA.eigvals(oc_tg,N=30)[end])))