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))
../_images/d5d66e2d9f9b88017e9db1e0b445623095f4d0d38a7fbef40c576f8823ef04dc.svg

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

  • -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}\]