# 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, -wind, 0, wind, wind] / 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)]
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.

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' * u
scatter!(P' * x, u_2)
plot!(x, P * u_2, marker=:auto) # 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
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) # 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' * A * P
S2 = P * inv(A2) * P' * A
norm(S2^2 - S2)

1.2036131040197588e-15

norm(P' * 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