2023-10-27 AMG

Last time

• Multigrid

  • Spectral perspective

  • Factorization perspective

Today

• Algebraic multigrid

• Smoothed aggregation

• Multigrid in PETSc

using Plots
default(linewidth=3)
using LinearAlgebra
using SparseArrays

function my_spy(A)
    cmax = norm(vec(A), Inf)
    s = max(1, ceil(120 / size(A, 1)))
    spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end

function laplace1d(n)
    "1D Laplacion with Dirichlet boundary conditions eliminated"
    h = 2 / (n + 1)
    x = LinRange(-1, 1, n+2)[2:end-1]
    A = diagm(0 => ones(n),
        -1 => -.5*ones(n-1),
        +1 => -.5*ones(n-1))
    x, A # Hermitian(A) xyz
end

function symbol(S, theta)
    if length(S) % 2 != 1
        error("Length of stencil must be odd")
    end
    w = length(S) ÷ 2
    phi = exp.(1im * (-w:w) * theta')
    vec(S * phi) # not! (S * phi)'
end

function plot_symbol(S, deriv=2; plot_ref=true, n_theta=30)
    theta = LinRange(-pi, pi, n_theta)
    sym = symbol(S, theta)
    rsym = real.((-1im)^deriv * sym)
    fig = plot(theta, rsym, marker=:circle, label="stencil")
    if plot_ref
        plot!(fig, th -> th^deriv, label="exact")
    end
    fig
end

function interpolate(m, stride=2)
    s1 = (stride - 1) / stride
    s2 = (stride - 2) / stride
    P = diagm(0 => ones(m),
        -1 => s1*ones(m-1), +1 => s1*ones(m-1),
        -2 => s2*ones(m-2), +2 => s2*ones(m-2))
    P[:, 1:stride:end]
end
n = 50; x, A = laplace1d(n)
P = interpolate(n, 2)
plot(x, P[:, 4:6], marker=:auto, xlims=(-1, 0))

Galerkin approximation of \(A\) in coarse space

\[ A_{2h} u_{2h} = P^T A_h P u_{2h} \]

n = 100
x, A = laplace1d(n)
P = interpolate(n)
@show size(A)
A_2h = P' * A * P
@show size(A_2h)
L_2h = eigvals(A_2h)
plot(L_2h, title="cond $(L_2h[end]/L_2h[1])")

size(A) = (100, 100)
size(A_2h) = (50, 50)

my_spy(A_2h)

Coarse grid correction

Consider the \(A\)-orthogonal projection onto the range of \(P\),

\[ S_c = P A_{2h}^{-1} P^T A \]

Sc = P * (A_2h \ P' * A)
Ls, Xs = eigen(I - Sc)
scatter(real.(Ls))

• This spectrum is typical for a projector. If \(u\) is in the range of \(P\), then \(S_c u = u\). Why?

• For all vectors \(v\) that are \(A\)-orthogonal to the range of \(P\), we know that \(S_c v = 0\). Why?

A two-grid method

w = .67
T = (I - w*A)^1 * (I - Sc) * (I - w*A)^2
Lt = eigvals(T)
scatter(Lt)

• Can analyze these methods in terms of frequency.

• LFAToolkit

Multigrid as factorization

We can interpret factorization as a particular multigrid or domain decomposition method.

We can partition an SPD matrix as

\[\begin{split}A = \begin{bmatrix} F & B^T \\ B & D \end{bmatrix}\end{split}\]

and define the preconditioner by the factorization

\[\begin{split} M = \begin{bmatrix} I & \\ B F^{-1} & I \end{bmatrix} \begin{bmatrix} F & \\ & S \end{bmatrix} \begin{bmatrix} I & F^{-1}B^T \\ & I \end{bmatrix} \end{split}\]

where \(S = D - B F^{-1} B^T\). \(M\) has an inverse

\[\begin{split} \begin{bmatrix} I & -F^{-1}B^T \\ & I \end{bmatrix} \begin{bmatrix} F^{-1} & \\ & S^{-1} \end{bmatrix} \begin{bmatrix} I & \\ -B F^{-1} & I \end{bmatrix} . \end{split}\]

Define the interpolation

\[\begin{split} P_f = \begin{bmatrix} I \\ 0 \end{bmatrix}, \quad P_c = \begin{bmatrix} -F^{-1} B^T \\ I \end{bmatrix} \end{split}\]

and rewrite as

\[\begin{split} M^{-1} = [P_f\ P_c] \begin{bmatrix} F^{-1} P_f^T \\ S^{-1} P_c^T \end{bmatrix} = P_f F^{-1} P_f^T + P_c S^{-1} P_c^T.\end{split}\]

The iteration matrix is thus

\[ I - M^{-1} A = I - P_f F^{-1} P_f^T J - P_c S^{-1} P_c^T A .\]

Permute into C-F split

m = 64
x, A = laplace1d(m)
my_spy(A)

idx = vcat(1:2:m, 2:2:m)
J = A[idx, idx]
my_spy(J)

Coarse basis functions

F = J[1:end÷2, 1:end÷2]
B = J[end÷2+1:end,1:end÷2]
P = [-F\B'; I]
my_spy(P)

idxinv = zeros(Int64, m)
idxinv[idx] = 1:m
Pp = P[idxinv, :]
plot(x, Pp[:, 3:5], marker=:auto, xlims=(-1, -.5))

From factorization to algebraic multigrid

• Factorization as a multigrid (or domain decomposition) method incurs significant cost in multiple dimensions due to lack of sparsity.

  • We can’t choose enough coarse basis functions so that \(F\) is diagonal, thereby making the minimal energy extension \(-F^{-1} B^T\) sparse.

• Algebraic multigrid

  • Use matrix structure to aggregate or define C-points

  • Create an interpolation rule that balances sparsity with minimal energy

Aggregation

Form aggregates from “strongly connected” dofs.

agg = 1 .+ (0:m-1) .÷ 3
mc = maximum(agg)
T = zeros(m, mc)
for (i, j) in enumerate(agg)
    T[i,j] = 1
end
plot(x, T[:, 3:5], marker=:auto, xlims=(-1, -.5))

Sc = T * ((T' * A * T) \ T') * A
w = .67; k = 1
E = (I - w*A)^k * (I - Sc) * (I - w*A)^k
scatter(eigvals(E), ylims=(-1, 1))

• simple and cheap method

• stronger smoothing (bigger k) doesn’t help much; need more accurate coarse grid

Smoothed aggregation

\[ P = (I - w A) T \]

w = 2/3
P = (I - w * A) * T
plot(x, P[:, 3:5], marker=:auto, xlims=(-1,-.5))

Sc = P * ((P' * A * P) \ P') * A
w = 2/3; k = 3
E = (I - w*A)^k * (I - Sc) * (I - w*A)^k
scatter(eigvals(E), ylims=(-1, 1))

• Eigenvalues are closer to zero; stronger smoothing (larger k) helps.

• Smoother can be made stronger using Chebyshev (like varying the damping between iterations in Jacobi)

Multigrid in PETSc

Geometric multigrid

• -pc_type mg

  • needs a grid hierarchy (automatic with DM)

  • PCMGSetLevels()

  • PCMGSetInterpolation()

  • PCMGSetRestriction()

• -pc_mg_galerkin

• -pc_mg_cycle_type [v,w]

• -mg_levels_ksp_type chebyshev

• -mg_levels_pc_type jacobi

• -mg_coarse_pc_type svd

  • -mg_coarse_pc_svd_monitor (report singular values/diagnose singular coarse grids)

Algebraic multigrid

• -pc_type gamg

  • native PETSc implementation of smoothed aggregation (and experimental stuff), all -pc_type mg options apply.

  • MatSetNearNullSpace()

  • MatNullSpaceCreateRigidBody()

  • DMPlexCreateRigidBody()

  • -pc_gamg_threshold .01 drops weak edges in strength graph

• -pc_type ml

  • similar to gamg with different defaults

• -pc_type hypre

  • Classical AMG (based on C-F splitting)

  • Manages its own grid hierarchy

Integration by parts

One dimension

it’s the product rule backwards

(39)\[\begin{align} \int_a^b d(uv) &= \int_a^b u dv + \int_a^b v du \\ (uv)_a^b &= \int_a^b u dv + \int_a^b v du \\ \int_a^b u dv &= (uv)_a^b - \int_a^b v du \end{align}\]

you can move the derivative to the other term; it’ll cost you a minus sign and a boundary term

Multiple dimensions

(40)\[\begin{align} \int_\Omega v \nabla\cdot \mathbf f = -\int_\Omega \nabla v \cdot \mathbf f + \int_{\partial \Omega} v \mathbf f \cdot \mathbf n \end{align}\]

Strong form

\[ -\nabla\cdot(\kappa \nabla u) = 0 \]

Weak form

• multiply by a test function and integrate by parts

(41)\[\begin{align} -\int_\Omega v \nabla\cdot(\kappa \nabla u) = 0, \forall v \\ \int_\Omega \nabla v \cdot \kappa \nabla u - \int_{\partial\Omega} v \underbrace{\kappa \nabla u \cdot \mathbf n}_{\text{boundary condition}} = 0 \end{align}\]