Contents

2021-10-25 Domain Decomposition

Last time

  • Experiments with preconditioners and diagnostics in PETSc

  • Basic preconditioners and incomplete factorization

Today

  • Domain decomposition as a framework

  • Asymptotics for one-level methods

    • constants, exact and inexact subdomain solves

  • Multigrid as a domain decomposition method

using Plots
using LinearAlgebra
using SparseArrays

default(linewidth=4)

function advdiff_matrix(n; kappa=1, wind=[0, 0])
    "Advection-diffusion with Dirichlet boundary conditions eliminated"
    h = 2 / (n + 1)
    rows = Vector{Int64}()
    cols = Vector{Int64}()
    vals = Vector{Float64}()
    idx((i, j),) = (i-1)*n + j
    in_domain((i, j),) = 1 <= i <= n && 1 <= j <= n
    stencil_advect = [-wind[1], -wind[2], 0, wind[1], wind[2]] / h
    stencil_diffuse = [-1, -1, 4, -1, -1] * kappa / h^2
    stencil = stencil_advect + stencil_diffuse
    for i in 1:n
        for j in 1:n
            neighbors = [(i-1, j), (i, j-1), (i, j), (i+1, j), (i, j+1)]
            mask = in_domain.(neighbors)
            append!(rows, idx.(repeat([(i,j)], 5))[mask])
            append!(cols, idx.(neighbors)[mask])
            append!(vals, stencil[mask])
        end
    end
    sparse(rows, cols, vals)
end

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

my_spy (generic function with 1 method)

Experiments in PETSc

  • Using src/snes/tutorials/ex15.c

    • Newton-Krlov for “p-Bratu”: combines p-Laplacian with Bratu nonlinearity

      \[ -\nabla\cdot\big(\kappa(\nabla u) \nabla u\big) - \lambda e^u = f \]

-pc_type jacobi

\[ M_{\text{Jacobi}}^{-1} = D^{-1} \]

  • -pc_jacobi_type [diagonal|rowmax|rowsum]

-pc_type sor

\[ M_{GS}^{-1} = (L+D)^{-1} \]

or the symmetric form

\[ M_{SGS}^{-1} = (L+U)^{-1} A \Big( I - (L+D)^{-1} \Big) . \]

  • -pc_sor_omega 1. is default (Gauss-Seidel)

-pc_type ilu

  • -pc_factor_levels 2

-pc_type lu

  • -pc_factor_mat_solver_type [petsc|umfpack|mumps|superlu]

  • -pc_factor_mat_ordering_type [natural|nd|rcm]

Krylov diagnostics

  • -ksp_monitor

  • -ksp_monitor_true_residual

  • -ksp_monitor_singular_value

  • -ksp_view_eigenvalues [draw]

Domain decomposition

Alternating Schwarz method

bc_circ = guess()
while not converged:
    u_circ = solve(A_circ, bc_circ)
    bc_rect = eval(u_circ, rect)
    u_rect = solve(A_rect, bc_rect)
    bc_circ = eval(u_rect, circ)

This method was proposed in 1870 by Hermann Schwarz as a theoretical tool, and is now called a multiplicative Schwarz method because the solves depend on each other. We can see it as a generalization of Gauss-Seidel in which we solve on subdomains instead of at individual grid points. As with Gauss-Seidel, it is difficult to expose parallelism.

Additive Schwarz

The additive Schwarz method is more comparable to Jacobi, with each domain solved in parallel. Our fundamental operation will be an embedding of each subdomain \(V_i\) into the global domain \(V\), which we call the prolongation

\[ P_i : V_i \to V \]

The transpose of prolongation, \(P_i^T\), will sometimes be called restriction. Let’s work an example.

n = 21
x = LinRange(-1, 1, n)
overlap = 0
domains = [(1, n÷3), (n÷3+1, 2*n÷3), (2*n÷3+1, n)]
Id = diagm(ones(n))
P = []
for (i, (left, right)) in enumerate(domains)
    push!(P, Id[:, left:right])
end

u = 1 .+ cos.(3*x)
plot(x, u)
u_2 = P[2]' * u
scatter!(P[2]' * x, u_2)
plot!(x, P[2] * u_2, marker=:auto)
../_images/2021-10-25-domain-decomposition_10_0.svg

The algorithm

Define the subdomain operator

\[A_i = P_i^T A P_i\]
The additive Schwarz preconditioner is
\[ M^{-1} = \sum_i P_i A_i^{-1} P_i^T \]
and it has error iteration matrix
\[ I - \sum_i P_i A_i^{-1} P_i^T A \]

A = diagm(0 => 2*ones(n), -1 => -ones(n-1), 1 => -ones(n-1))
b = zeros(n); b[1] = 1
u = zero(b) # initial guess
us = [u]
for _ in 1:5
    r = b - A * u # residual
    u_next = copy(u);
    for Pi in P
        Ai = Pi' * A * Pi
        ui = Ai \ (Pi' * r)
        u_next += Pi * ui
    end
    u = u_next; push!(us, u)
end
plot(x, us, marker=:auto)
../_images/2021-10-25-domain-decomposition_13_0.svg

Theory

Given a linear operator \(A : V \to V\), suppose we have a collection of prolongation operators \(P_i : V_i \to V\). The columns of \(P_i\) are “basis functions” for the subspace \(V_i\). The Galerkin operator \(A_i = P_i^T A P_i\) is the action of the original operator \(A\) in the subspace.

Define the subspace projection

\[ S_i = P_i A_i^{-1} P_i^T A . \]

  • \(S_i\) is a projector: \(S_i^2 = S_i\)

  • If \(A\) is SPD, \(S_i\) is SPD with respect to the \(A\) inner product \(x^T A y\)

  • \(I - S_i\) is \(A\)-orthogonal to the range of \(P_i\)

A2 = P[2]' * A * P[2]
S2 = P[2] * inv(A2) * P[2]' * A
norm(S2^2 - S2)

1.2036131040197588e-15

norm(P[2]' * A * (I - S2))

1.1831170202396054e-15

Note: The concept of \(A\)-orthogonality is meaningful only when \(A\) is SPD. Does the mathematical expression \( P_i^T A (I - S_i) = 0 \) hold even when \(A\) is nonsymmetric?

Convergence theory

The formal convergence is beyond the scope of this course, but the following estimates are useful. We let \(h\) be the element diameter, \(H\) be the subdomain diameter, and \(\delta\) be the overlap, each normalized such that the global domain diameter is 1. We express the convergence in terms of the condition number \(\kappa\) for the preconditioned operator.

  • (Block) Jacobi: \(\delta=0\), \(\kappa \sim H^{-2} H/h = (Hh)^{-1}\)

  • Overlapping Schwarz: \(\kappa \sim H^{-2} H/\delta = (H \delta)^{-1}\)

  • 2-level overlapping Schwarz: \(\kappa \sim H/\delta\)

  • Additive \( I - \sum_{i=0}^n S_i, \)

  • Multiplicative \( \prod_{i=0}^n (I - S_i), \)

  • Hybrid \( (I - S_0) (I - \sum_{i=1}^n S_i) . \)

In each case above, the action is expressed in terms of the error iteration operator.

PETSc experiments

  • Compare domain decomposition and multigrid preconditioning

  • -pc_type asm (Additive Schwarz)

  • -pc_asm_type basic (symmetric, versus restrict)

  • -pc_asm_overlap 2 (increase overlap)

  • Effect of direct subdomain solver: -sub_pc_type lu

  • Symmetric example: src/snes/examples/tutorials/ex5.c

  • Nonsymmetric example: src/snes/examples/tutorials/ex19.c

  • Compare preconditioned versus unpreconditioned norms.

  • Compare BiCG versus GMRES

  • -pc_type mg (Geometric Multigrid)

  • Use monitors:

  • -ksp_monitor_true_residual

  • -ksp_monitor_singular_value

  • -ksp_converged_reason

  • Explain methods: -snes_view

  • Performance info: -log_view