Assessing the algorithmic scaling behavior of PDE solvers Matthias - - PowerPoint PPT Presentation

Assessing the algorithmic scaling behavior of PDE solvers

Matthias Bolten

University of Wuppertal

2015/09/17

Motivation Stationary problems Time-dependent problems Wrap-up

Outline

Motivation Stationary problems Model problem Discretization Direct solvers Iterative methods Time-dependent problems Model problem Explicit methods Implicit methods Wrap-up

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 1/57

Motivation Stationary problems Time-dependent problems Wrap-up

Outline

Motivation Stationary problems Model problem Discretization Direct solvers Iterative methods Time-dependent problems Model problem Explicit methods Implicit methods Wrap-up

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 2/57

Motivation Stationary problems Time-dependent problems Wrap-up

Motivation

◮ Many problems in scientific computing involve solution of a

PDE

◮ A lot of optimized methods exist, providing good efficiency in

terms of FLOPS

◮ Several scalable methods exist, providing good efficiency in

terms of—well—efficiency

◮ A lot of effort is put on optimizing scientific applications ◮ The underlying method is rarely questioned ◮ But: Many methods’ algorithmic scalability does not scale ◮ Bad algorithmic scaling leads to bad weak scaling behavior

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 3/57

Motivation Stationary problems Time-dependent problems Wrap-up

Motivation: Example

Assume there is an application solving an elliptic partial differential equation in 3d. Now, consider two cases:

• 1. Solution of the PDE using Gaussian elimination at 80% peak
• 2. Solution of the PDE using full multigrid at 0.5% peak

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 4/57

Motivation Stationary problems Time-dependent problems Wrap-up

Motivation: Example

Assume there is an application solving an elliptic partial differential equation in 3d. Now, consider two cases:

• 1. Solution of the PDE using Gaussian elimination at 80% peak
• 2. Solution of the PDE using full multigrid at 0.5% peak

Which is the better approach?

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 4/57

Motivation Stationary problems Time-dependent problems Wrap-up

Motivation: Example

Assume there is an application solving an elliptic partial differential equation in 3d. Now, consider two cases:

• 1. Solution of the PDE using Gaussian elimination at 80% peak
• 2. Solution of the PDE using full multigrid at 0.5% peak

Which is the better approach? Answer: Method 1, as it provides more impressive performance numbers

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 4/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Outline

Motivation Stationary problems Model problem Discretization Direct solvers Iterative methods Time-dependent problems Model problem Explicit methods Implicit methods Wrap-up

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 5/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Model problem: Poisson equation

Let Ω = [0, 1]2. Consider the Poisson equation given by −∆u(x) = f(x), for x ∈ Ω and u(x) = 0 for x ∈ ∂Ω. (1) This is the prototypical elliptic partial differential equation (PDE). A vast amount of other elliptic PDEs exist with

◮ higher dimensionality, ◮ more complex geometries, ◮ varying coefficients, ◮ additional (convective) terms, ◮ . . . .

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 6/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Discretization: Finite differences I

◮ The simplest discretization is obtained by approximating the

differential quotient d dxf(x)|x=x0 = lim

h→0

f(x0 + h) − f(x0) h ≈ f(x0 + h) − f(x0) h

◮ Discretization on regular grid yields stencil in each grid point ◮ For the Laplace operator ∆ we obtain the stencil

1 h2   1 1 −4 1 1  

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 7/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Discretization: Finite differences II

◮ Discretization on grid with N = (n + 1) · (n + 1) grid points

using 5-point scheme yields linear system Lu = f, where − 1 h2 (4ui,j − ui−1,j − ui+1,j − ui,j−1 − ui,j+1) = fi,j, for i, j = 1, . . . , nk, where h = 1/n and ui,j = 0 for i ∈ {0, n + 1} or j ∈ {0, n + 1}.

◮ Eigenvalues and eigenvectors given by

λl,m = 4 − 2 cos(lπh) − 2 cos(mπh), (2) (ϕl,m)i,j = sin(lπih) sin(mπjh), (3)

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 8/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Discretization: Finite elements I

Instead of (5) one can consider its weak formulation, i.e., find u vanishing on ∂Ω s.t. −

• Ω

(∆u)v dx =

• Ω

f vdx for all test functions v. Using Green’s identity we obtain −

• Ω

∇u · ∇v dx =

• Ω

f vdx. The functions u and v have to be chosen from proper function spaces.

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 9/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Discretization: Finite elements II

(Conforming) finite element methods use finite subspaces of these function spaces for both u and v. To obtain this subspace one proceeds as follows: The domain is decomposed into polyhedrons, often triangles are chosen, resulting in a triangulation like the following:

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 10/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Discretization: Finite elements III

The function space under consideration consists of all continuous functions vanishing on the boundary that restricted to the triangles are polynomials of a certain degree p, e.g., p = 1. As basis we choose functions vanishing at all nodes but one, e.g.:

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 11/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Discretization: Finite elements IV

For a regular grid we can choose the following triangulation: Choosing piecewise linear basis functions yields the same stencil as for finite differences: 1 h2   1 1 −4 1 1   .

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 12/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Discretization: General considerations

◮ In either case, problem leads to linear system

Ax = b, (4) where A ∈ Cn×n, x, b ∈ Cn

◮ Even exact solution x is only approximation to continuous

problem, not a sampling of the continuous solution

◮ Error of exact solution is known as discretization error ◮ Error estimates exist ◮ Solution of (4) up to machine precision not required ◮ Additional algebraic error affordable

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 13/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Direct solvers

◮ Usually based on factorization of system matrix A ◮ Well-known methods:

◮ A = LU: LU factorization (Gaussian elimination) ◮ A = LDLT : LDL factorization ◮ A = LL∗: Cholesky factorization (A hermitian, positive

(semi-)definite)

◮ Direct solvers for sparse matrices exist, including heuristics to

minimize fill-in

◮ Parallel implementations of the latter available ◮ Also include special solvers, e.g. using FFTs

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 14/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Direct solvers

◮ Usually based on factorization of system matrix A ◮ Well-known methods:

◮ A = LU: LU factorization (Gaussian elimination) ◮ A = LDLT : LDL factorization ◮ A = LL∗: Cholesky factorization (A hermitian, positive

(semi-)definite)

◮ Direct solvers for sparse matrices exist, including heuristics to

minimize fill-in

◮ Parallel implementations of the latter available ◮ Also include special solvers, e.g. using FFTs ◮ Provide exact solution of (4)

(up to machine precision, subject to round-off)

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 14/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Direct solvers: Algorithmic scaling behavior

◮ Factorizations usually scale like O(N3) ◮ N is the number of unknowns ◮ Doubling resolution in 2d results in 64 times the run time ◮ Doubling resolution in 3d results in 512 times the run time

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 15/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Direct solvers: Algorithmic scaling behavior

◮ Factorizations usually scale like O(N3) ◮ N is the number of unknowns ◮ Doubling resolution in 2d results in 64 times the run time ◮ Doubling resolution in 3d results in 512 times the run time ◮ Specialized solvers like FFT-based methods scale better

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 15/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Iterative solvers

◮ Start with initial guess x(0) ◮ Construct sequence {x(k)}∞ k=0 of approximate solutions ◮ Two classes: stationary and non-stationary methods ◮ Stationary methods characterized by iteration matrix M ◮ Examples for stationary methods:

◮ Jacobi method (M = I − D−1A) ◮ Gauss-Seidel methods (M = I − (D − L)−1A) ◮ Successive over-relaxation (SOR, M(ω) = I − ω(D − ωL)−1A) ◮ Symmetric SOR (SSOR, M(ω) = . . . )

◮ Example for non-stationary methods:

◮ Krylov-subspace methods Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 16/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Iterative solvers

◮ Start with initial guess x(0) ◮ Construct sequence {x(k)}∞ k=0 of approximate solutions ◮ Two classes: stationary and non-stationary methods ◮ Stationary methods characterized by iteration matrix M ◮ Examples for stationary methods:

◮ Jacobi method (M = I − D−1A) ◮ Gauss-Seidel methods (M = I − (D − L)−1A) ◮ Successive over-relaxation (SOR, M(ω) = I − ω(D − ωL)−1A) ◮ Symmetric SOR (SSOR, M(ω) = . . . )

◮ Example for non-stationary methods:

◮ Krylov-subspace methods

◮ Provide approximate solution to (4), only

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 16/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Krylov subspace methods

◮ Non-stationary iterative methods ◮ Approximation of solution in Krylov subspace:

Kk(A, b) = b, Ab, A2b, . . . , Ak−1b ⊆ Cn

◮ Krylov matrix defined by

Kk =

• b
• Ab
• A2b
• . . .
• Ak−1b
• ◮ Kk has reduced QR factorization

Kk = QkRk

◮ Basis Qk of Kk created by Arnoldi iteration

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 17/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Arnoldi iteration I

Suppose, to pass the time while marooned on a desert island, you challenged yourself to devise an algorithm to reduce a nonhermitian matrix to Hessenberg form by or- thogonal similarity transformations, proceeding column by column from a prescribed first column q1. To your surprise you could solve this problem in an hour and still have time to gather coconuts for dinner. The method you would come up with goes by the name of the Arnoldi iteration.

(Lloyd N. Trefethen, Numerical Linear Algebra, SIAM, Philadelphia, 1997)

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 18/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Arnoldi iteration II

◮ Complete orthogonal similarity transform given by

A = QHQ∗ ⇔ AQ = QH

◮ Let Qk ∈ Cn×k consist of the first k columns of Q ◮ Define ˜

Hk as upper-left section of H ˜ Hk =        h1,1 · · · h1,k h2,1 h2,2 ... ... . . . hk,k−1 hk,k hk+1,k       

◮ Then AQk = Qk+1 ˜

Hk (Arnoldi relation)

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 19/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Arnoldi iteration III

◮ k-th column given by (k + 1)-term recurrence:

Aqk = h1,kq1 + · · · + hk,kqk + hk+1,kqk+1 Algorithm (Arnoldi iteration)

b = arbitrary, q1 = b/b for k = 1, 2, . . . do v = Aqk for j = 1, . . . , i do hj,k = q∗

j v

v = v − hj,kqj end for hk+1,k = v qk+1 = v/hk+1,k end for

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 20/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

GMRES

◮ GMRES computes x(k) = Kkc(k) ∈ Kk, c(k) ∈ Ck s.t.

AKkc(k) − b = minimum

◮ Set x(k) = Qky(k) and use Arnoldi relation to reduce linear

system to AQky(k) − b = minimum ⇔ Qk+1 ˜ Hky(k) − b = minimum

◮ Properties of Qk+1 finally yield

˜ Hky − be1 = minimum

◮ Problem is reduced to (k + 1) × k matrix least squares problem ◮ Work further reduced by updating QR factorization

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 21/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Restarted GMRES

◮ Size of least squares problem grows with number of iterations ◮ Size of Qk also grows, increased memory requirement ◮ Approach to limit both grows: Restarting after m iterations ◮ Downside: Method is not guaranteed to converge anymore

Algorithm (restarted GMRES, GMRES(m))

for j = 0, 1, . . . do r(j·m) = b − Ax(j·m) q1 = r(j·m)/r(j·m) for k = 1, . . . , m do Step k of Arnoldi iteration Find y to minimize ˜ Hky − re1 x(j·m+k) = x(j·m) + Qky end for end for

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 22/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Lanzcos iteration I

◮ For hermitian matrices Hessenberg matrix Hk becomes

tridiagonal matrix Tk: Tk =         α1 β1 β1 α2 β2 β2 α3 ... ... ... βk−1 βk−1 αk        

◮ (k + 1)-term recurrence reduces to 3-term recurrence ◮ Arnoldi iteration reduces to Lanczos iteration ◮ GMRES reduces to MINRES

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 23/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Lanzcos iteration II

Algorithm (Lanczos iteration) β0 = 0, q0 = 0, b = arbitrary, q1 = b/b for k = 1, 2, . . . do v = Aqk αk = q∗

kv

v = v − βk−1qk−1 − αkqk βk = v qk+1 = v/βk end for

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 24/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Conjugate gradients (CG)

◮ Can be interpreted as an optimization of the functional

ϕ(x) = 1/2x∗Ax − x∗b

◮ Minimizes A-norm of the error (instead of 2-norm of residual) ◮ Error satisfies ekA/e0A ≤ 2((√κ − 1)/(√κ + 1))k

Algorithm (CG iteration)

x(0) = 0, r(0) = 0, p(0) = 0 for k = 1, 2, . . . do αk = ((r(k−1))∗r(k−1))/((p(k−1))∗Ap(k−1)) x(k) = x(k−1) + αkp(k−1) r(k) = r(k−1) − αkp(k−1) βk = ((r(k))∗r(k))/((r(k−1))∗r(k−1)) p(k) = r(k) + βkp(k−1) end for

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 25/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Biorthogonalization

◮ CG for nonhermitian matrices using normal equation

Ax = b ⇔ A∗Ax = A∗b

◮ Squared condition number leads to slow convergence ◮ Alterative: Tridiagonal biorthogonalization

A = V TV −1 and A∗ = V −∗

• =:W

T ∗(V −∗)−1 = WT ∗W −1 where W ∗V = V ∗W = I

◮ Columns of V and W span K(A, v1) and K(A∗, w1) ◮ Leads to BiCG and variants (including BiCGstab, QMR,. . . )

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 26/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Krylov subspace methods: Algorithmic scaling behavior

◮ Krylov subspace methods only need multiplication with

system matrix A

◮ If A is sparse with a fixed number of unknowns per

row—usually the case for PDEs—the cost for matrix-vector multiplication and thus for each iteration is linear

◮ Number of iterations depends on error reduction ◮ E.g., for CG we have the following bound

ekA/e0A ≤ 2((√κ − 1)/(√κ + 1))k

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 27/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Preconditioning I

◮ Given nonsingular M ∈ Cn×n equation Ax = b equivalent to

M−1Ax = M−1b

◮ Convergence of iterative solver now depends on M−1A ◮ M should be chosen such that convergence rate is improved

(optimal: M = A)

◮ Linear systems with system matrix M should be easy to solve ◮ This is called (left) preconditioning ◮ Right preconditioning: Solve AM−1y = b, then Mx = y ◮ Hermitian preconditioning: M = CC∗, system transformed to

(C−1AC−∗)C∗x = C−1b

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 28/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Preconditioning II

◮ Simple preconditioner: Scale A by M = diag(A) like in Jacobi ◮ Easy to implement ◮ Often very effective ◮ More general: Scaling with M = diag(c) for some vector

c ∈ Cn

◮ Extension: Apply Gauss-Seidel, SOR or SSOR

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 29/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Preconditioning III

◮ Often preconditioners based on factorizations are used ◮ Gaussian elimination produces LU factorization

A = LU

◮ Even for sparse A L and U are usually not as sparse ◮ Incomplete factorization obtained by allowing nonzeros only,

where A was nonzero

◮ Same is possible for Cholesky factorization ◮ Extension possible by introducing a drop tolerance, multiple

levels, . . .

◮ Not as easy to implement ◮ Hard to parallelize

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 30/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Preconditioning: Algorithmic scaling behavior

◮ Cost of one application of the preconditioner varies a lot ◮ Estimating resulting number of iterations requires analysis of

preconditioned system matrix, e.g. of M−1A

◮ E.g. straightforward application of error estimate for CG

results in ekM−1A/e0M−1A ≤ 2((√κ − 1)/(√κ + 1))k, where κ is now the condition of M−1A instead of A

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 31/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Domain Decomposition I

◮ Split the computational domain into subdomains Bi ◮ Solve system (iteratively) on each subdomain B1 B2 B3 B4 ◮ Canonical injection Ij

Ijei = e(Bj)i

◮ Restriction of x onto Bj

xBj = I†

j x ◮ Restriction of A onto Bj

ABj = I†

j AIj

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 32/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Domain decomposition II

1 2 3 4 5 6 7 8 9 10 11 12 13 14

◮ Consider simple restriction returning only indicated nodes with

indices in B = {1, 2, 5} R =   1 · · · 1 · · · 1 · · ·  

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 33/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Domain decomposition III

◮ transpose of restriction RT inserts values at the given nodes ◮ want to correct at given nodes ◮ best correction defined by

w∗ = min

w u∗ − (un + RT w)2 A ◮ leads to

w∗ = RT (RART )−1R(f − Aun)

◮ AB = RART is subblock of A associated with given nodes

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 34/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Domain decomposition IV

Additive Schwarz

for k = 0, 1, . . . do r(k) = b − Ax(k) for j = 1, 2, . . . , nB do x(k+1)

Bj

= x(k)

Bj + A−1 Bj r(k) Bj

end for end for

Multiplicative Schwarz

for k = 0, 1, . . . do for j = 1, 2, . . . , nB do r = b − Ax xBj = xBj + A−1

Bj rBj

end for end for

◮ Block-Jacobi ◮ Embarrassingly Parallel

Schwarz methods in general ⊕ Data parallel ⊕ Computation parallel

◮ Block-Gauss-Seidel ◮ Sequential (→ coloring)

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 35/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Multilevel domain decomposition I

◮ as for other iterative methods, performance deteriorates when

number of processors/subdomains is increased

◮ result of information transfer needed within whole domain for

solution of elliptic PDE

◮ algorithmic scalability not given even for simple elliptic

problems

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 36/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Multilevel domain decomposition II

◮ to overcome this limitation, coarse grid operator is constructed ◮ similar to multigrid, but higher coarsening ratio

(10 − 100 instead of 2 − 4)

◮ if AC is the coarse operator and R is the restriction from the

fine grid, approximate correction can be written as RT A−1

C Rr,

where r is the residual

◮ coarse grid problem can be solved recursively using the same

idea

◮ coarse grid correction is – under certain circumstances –

projection of error onto subspace of solution space

◮ multilevel Schwarz methods have good and desirable

convergence properties

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 37/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Domain decomposition: Algorithmic scaling behavior

◮ A variety of domain decomposition methods exists ◮ Cost of one iteration is the accumulated cost of the solves on

the subdomain plus the coarse level solve

◮ In general single-level domain decomposition shows a

deterioration in convergence rate as the problem size grows or the subdomains shrink

◮ Introduction of a coarse level removes this effect ◮ Coarse level usually gets more expensive when the problem

size increases

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 38/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Multigrid I

Fewer

First Coarse Grid Finest Grid Smooth

The Multigrid V−cycle Restriction Prolongation Dofs

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 39/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Multigrid II

Multigrid cycle xni = MGi(xni, bni) xni ← Sν1

i (xni, bni)

rni ← bni − Aixni rni+1 ← Rirni eni+1 ← 0 if i + 1 = lmax then enlmax ← A−1

lmaxrnlmax

else for j = 1, . . . , γ do eni+1 ← MGi+1(eni+1, rni+1) end for end if eni ← Pienni+1 xni ← xni + eni xni ← ˜ Sν2

i (xni, bni)

coarse grid correction

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 40/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Multigrid methods: Algorithmic scaling behavior

◮ Number of grid points over the whole hierarchy grows linear

as long as the coarsening ratio is strictly smaller than 1

◮ Smoother usually is based on application of system matrix,

resulting in linear complexity

◮ Overall a multigrid iteration’s cost is linear in the number of

unknowns

◮ The convergence rate can be shown to be independent of the

grid size

◮ To compare different solvers, the convergence rate is needed

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 41/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Analyzing convergence rate

◮ In order to assess performance of iterative methods the

convergence rate is needed

◮ Convergence rate allows to predict number of iterations

necessary to achieve prescribed accuracy

◮ Convergence rate usually depends on spectral properties of

the underlying operator

◮ In many cases a local Fourier analysis (LFA) can provide

insight into the convergence behavior

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 42/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

Local Fourier analysis (LFA)

Idea Study effect of linear solvers by

◮ Freezing non-constant coefficients ◮ Neglecting boundary conditions

as a constant coefficient problem on an infinite grid.

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 43/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

LFA: Definitions

◮ d-dimension infinite grid:

Gh = {x = kh := (k1h1, . . . , kdhd), k ∈ Zd}

◮ grid function:

ϕ : Gh → C, x → Ψh(x)

◮ grid basis functions (θ ∈ [−π, π)d):

ϕ(θ, x) = eiθ,x/h

◮ discrete operator Lh with constant coefficients sl, l ∈ V

Lh : Gh × C → Gh × C Ψh(x) → LhΨh(x) =

• l∈V

slΨh(x + lh)

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 44/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

LFA: Relation to Toeplitz matrices

◮ operator Lh described by d-fold tensor product of infinite

Toeplitz matrices: T =      t0 t−1 t−2 · · · t1 t0 t−1 · · · t2 t1 t0 · · · . . . . . . . . . ...     

◮ grid basis functions ϕ(θ, x) (formal) eigenfunctions of Lh ◮ symbol ˜

Lh of Lh: ˜ Lh(θ) =

• l∈V

sleiθ,l

◮ for the (formal) eigenfunctions ϕ(θ, x) follows

Lhϕ(θ, x) = ˜ Lh(θ)ϕ(θ, x), x ∈ Gh

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 45/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

LFA: Analysis of point smoothers I

◮ smoothers have to be effective for modes not damped

efficiently by coarse grid correction

◮ coarse grid with grid spacing 2h given by:

G2h = {x = k2h := (k12h1, . . . , kd2hd), k ∈ Zd}

◮ on coarse grid: 2d − 1 grid functions ϕ(θ, ·) not

distinguishable from ϕ(θ′, ·), θ′ ∈ [−π/2, π/2)d

◮ set of harmonics:

Eθ

h := span{ϕ(θα, ·) : α ∈ {0, 1}d},

where (θα)j =

• θj + αjπ

if θj ≤ 0, θj − αjπ if θj > 0. , αj ∈ {0, 1}, j = 1, . . . , d.

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 46/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

LFA: Analysis of point smoothers II

◮ in splitting method Lhuh = fh can be written as

L+

h ¯

wh + L−

h wh = fh,

where wh is the old iterate and ¯ wh is the new one

◮ for the errors vh and ¯

vh we obtain L+

h ¯

vh + L−

h vh = 0

• r

¯ vh = Shvh

◮ Sh has same (formal) eigenvectors ϕ(θ, ·) and amplification

˜ Sh(θ) = − ˜ L−

h (θ)

˜ L+

h (θ) ◮ smoothing factor µ(Sh) defined as

µ(Sh) = sup{| ˜ Sh(θ)| : θ ∈ [−π, π)d\[−π/2, π/2)d}

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 47/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

LFA: Analysis of multigrid I

◮ ϕ(θ, ·) are formal eigenfunctions of Lh, so Eθ h are invariant

subspaces

◮ considering the coarse level as well results in intermixing of

harmonics, nevertheless Eθ

h are invariant subspaces ◮ iteration matrix of coarse grid correction:

K2h

h = Ih − Ih 2h(L2h)−1I2h h Lh ◮ iteration matrix of twogrid operator:

M2h

h

= Sν2

h K2h h Sν1 h , ◮ LFA analyzes effect of twogrid method on ψ ∈ Eθ h, i.e.

ψ =

• α∈{0,1}d

Aαϕ(θα, ·), Aα ∈ R

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 48/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Discretization Direct solvers Iterative methods

LFA: Analysis of multigrid II

◮ assume that Lh, Ih 2h, L2h and I2h h

are given by stencils

◮ K2h h

• n Eθ

h represented by (2d × 2d)-matrix:

ˆ K2h

h (θ) = ˆ

Ih − ˆ Ih

2h(θ)(ˆ

L2h(2θ))−1 ˆ I2h

h (θ)ˆ

Lh(θ) for each θ ∈ [−π/2, π/2)d, where ˆ Ih ∈ C2d×2d (identity) ˆ Lh(θ) ∈ C2d×2d (fine grid smybol) ˆ Ih

2h ∈ C2d×1

(prolongation symbol) ˆ L2h ∈ C\{0} (coarse grid symbol) ˆ I2h

h ∈ C1×2d

(restriction symbol)

◮ leads to (2d × 2d)-matrix representation of twogrid method ◮ more detailed analysis by 3-grid analysis etc.

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 49/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Explicit methods Implicit methods

Outline

Motivation Stationary problems Model problem Discretization Direct solvers Iterative methods Time-dependent problems Model problem Explicit methods Implicit methods Wrap-up

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 50/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Explicit methods Implicit methods

Model problem: Heat equation

Let Ω = [0, 1]2. Consider the heat equation given by ∂u(x, t) ∂t − α∆u(x, t) = 0, for x ∈ Ω and u(x) = 0 for x ∈ ∂Ω. (5) This is the prototypical parabolic partial differential equation (PDE). Like in the elliptic case many modifications exist. . .

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 51/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Explicit methods Implicit methods

Explicit methods

Discretizing the time derivative using forward finite differences yields the explicit Euler scheme:

u(x,t+τ)−u(x,t) τ

− α∆u(x, t) = 0 ⇔ u(x, t + τ) = u(x, t) + τα∆u(x, t). Given an discretization of ∆ the value in the next time step can be computed without solving a linear system

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 52/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Explicit methods Implicit methods

Explicit methods

Discretizing the time derivative using forward finite differences yields the explicit Euler scheme:

u(x,t+τ)−u(x,t) τ

− α∆u(x, t) = 0 ⇔ u(x, t + τ) = u(x, t) + τα∆u(x, t). Given an discretization of ∆ the value in the next time step can be computed without solving a linear system But: Explicit Euler is only stable for 2ατ ≤ h2.

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 52/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Explicit methods Implicit methods

Explicit methods: Algorithmic scaling behavior

◮ Higher spatial resolution requires much higher temporal

resolution

◮ If this is not required by the application, many more time

steps are necessary

◮ E.g., reducing grid spacing by 1/2 requires 4 times as many

time steps

◮ If this is required by the application the method provides good

scaling behavior

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 53/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Explicit methods Implicit methods

Implicit methods

Discretizing the time derivative using backward finite differences yields the implicit Euler scheme:

u(x,t+τ)−u(x,t) τ

− α∆u(x, t + τ) = 0 ⇔ u(x, t + τ) − α∆u(x, t + τ) = u(x, t), i.e., the solution of a shifted linear system is required to compute the value at the next time step. The solvers discussed before can be used for this purpose. Implicit Euler is unconditionally stable (no restriction on τ).

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 54/57

Motivation Stationary problems Time-dependent problems Wrap-up Model problem Explicit methods Implicit methods

Implict methods: Algorithmic scaling behavior

◮ Inherited from the linear solver used ◮ No restriction on time step

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 55/57

Motivation Stationary problems Time-dependent problems Wrap-up

Outline

Motivation Stationary problems Model problem Discretization Direct solvers Iterative methods Time-dependent problems Model problem Explicit methods Implicit methods Wrap-up

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 56/57

Motivation Stationary problems Time-dependent problems Wrap-up

Wrap up

◮ Solution of linear system arising from PDE discretization does

not have to be 100% accurate

◮ Iterative solver can do this job perfectly ◮ To judge performance of iterative solver, number of iterations

needed

◮ Number of iterations depends on convergence rate ◮ Optimal solvers’ convergence rate does not depend on

discretization parameter (grid spacing)

◮ For many solvers it does... ◮ Time dependent problems allow for explicit solution, no linear

solve involved

◮ Explicit methods pose restrictions on time step

Matthias Bolten, Assessing the algorithmic scaling behavior of PDE solvers 57/57