E aStencils Overview Task: Solve a PDE (efficiently). Multigrid - - PowerPoint PPT Presentation

The Mathematics of ExaStencils

Hannah Rittich

joint work with: Matthias Bolten and Karsten Kahl

Bergische Universit¨ at Wuppertal

April 15, 2015

E aStencils

Overview

◮ Task: Solve a PDE (efficiently). ◮ Multigrid methods is a framework for constucting algorithms. ◮ These methods are build from a set of different ingredients.

⇒ Which ingredients to choose?

◮ How do these ingredients influence (the computational and)

numerical performance of method?

◮ Estimate the (numerical) performance is needed.

⇒ Local Fourier Analysis. (What is this?)

Poisson’s Equation

∆u = f

• n

Ω u = g

• n

∂Ω where ∆u = ∂2u

∂x2 + ∂2u ∂y2

Ω = (0, 1) × (0, 1)

Poisson’s Equation

∆u = f

• n

Ω u = g

• n

∂Ω where ∆u = ∂2u

∂x2 + ∂2u ∂y2

Ω = (0, 1) × (0, 1)

Photo by Graham Richardson / CC BY

Poisson’s Equation

∆u = f

• n

Ω u = g

• n

∂Ω where ∆u = ∂2u

∂x2 + ∂2u ∂y2

Ω = (0, 1) × (0, 1)

0.5 1

• 1
• 0.5

0.5 1

Poisson’s Equation

∆u = f

• n

Ω u = g

• n

∂Ω where ∆u = ∂2u

∂x2 + ∂2u ∂y2

Ω = (0, 1) × (0, 1)

• 1
• 0.5

1 0.5 1 0.5

Poisson’s Equation

∆u = f

• n

Ω u = g

• n

∂Ω where ∆u = ∂2u

∂x2 + ∂2u ∂y2

Ω = (0, 1) × (0, 1) x y Ω u : Ω → R

Poisson’s Equation

∆u = f

• n

Ω u = g

• n

∂Ω where ∆u = ∂2u

∂x2 + ∂2u ∂y2

Ω = (0, 1) × (0, 1) x y Ω u : Ω → R x y Ωh uh : Ωh → R

Multigrid in One Slide

• 0.2
• 0.15
• 0.1
• 0.05

0.05 0.1 0.15 0.2 0.8 1 0.6 0.2 0.4 0.2 0.4 0.6 0.8 1

• 0.01
• 0.005

0.005 0.01 0.015 0.8 1 0.4 0.6 0.2 0.2 0.4 0.6 0.8 1

• 0.01
• 0.005

0.005 0.01 0.015 0.4 0.6 0.8 1 0.2 0.2 0.4 0.6 0.8 1

◮ For many simple iterative methods and

certain problems we have that the error is smooth after a few iterations.

◮ Smooth error can be represented on a

coarser grid. ⇒ Linear system of reduced size.

◮ Use recursion on the coarse problem.

choose u repeat smooth u ˜ e ← coarse approx. of e u ← u + ˜ e smooth u until satisfied

Gauß-Seidel Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Bad Parallelization...

Gauß-Seidel Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Bad Parallelization...

Gauß-Seidel Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Bad Parallelization...

Gauß-Seidel Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Bad Parallelization...

Gauß-Seidel Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Bad Parallelization...

Gauß-Seidel Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Bad Parallelization...

Red-Black GS Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 2 Msgs.

Red-Black GS Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 2 Msgs.

Red-Black GS Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 2 Msgs.

Red-Black GS Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 2 Msgs.

Red-Black GS Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 2 Msgs.

Red-Black GS Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 2 Msgs.

Red-Black GS Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 2 Msgs.

Red-Black GS Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 2 Msgs.

Red-Black GS Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 2 Msgs.

Red-Black GS Smoother

ux,y = 1

4(fx,y − ux,y−1 − ux−1,y − ux+1,y − ux,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 2 Msgs.

Jacobi Smoother

unew

x,y = 1 4(fx,y − uold x,y−1 − uold x−1,y − uold x+1,y − uold x,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 1 Msgs.

Jacobi Smoother

unew

x,y = 1 4(fx,y − uold x,y−1 − uold x−1,y − uold x+1,y − uold x,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 1 Msgs.

Jacobi Smoother

unew

x,y = 1 4(fx,y − uold x,y−1 − uold x−1,y − uold x+1,y − uold x,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 1 Msgs.

Jacobi Smoother

unew

x,y = 1 4(fx,y − uold x,y−1 − uold x−1,y − uold x+1,y − uold x,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 1 Msgs.

Jacobi Smoother

unew

x,y = 1 4(fx,y − uold x,y−1 − uold x−1,y − uold x+1,y − uold x,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 1 Msgs.

Jacobi Smoother

unew

x,y = 1 4(fx,y − uold x,y−1 − uold x−1,y − uold x+1,y − uold x,y+1) ◮ Parallelizable ◮ Per Iteration:

◮ ≈ 4

√ N/ √ P Words

• p. Proc.

◮ 1 Msgs.

Smoother Summary

GS RB-GS Jacobi Words — 4 √ N/ √ P 4 √ N/ √ P Msgs — 2 1 Sm. 0.499817 0.25 1 2-lv 0.192336 0.0737621 1 3-lv 0.198407 0.104202 1

Local Fourier Analysis

• 1. We seek for E for “some” operator E and
• 2. for ρ(E).

(since limn→∞ En1/n = ρ(E)) That is E = (I − ωD−1L) (Jacobi’s Method), E = I − PℓL−1

ℓ+1P T ℓ Lℓ

(Coarse Grid Correction) or Eℓ = Mν2

ℓ (I − Pℓ(I − Eγ ℓ+1)(Lℓ+1)−1RℓLℓ)Mν1 ℓ

Eℓmax = 0 (Multigrid Method).

Local Fourier Analysis

• 1. We seek for E for “some” operator E and
• 2. for ρ(E).

(since limn→∞ En1/n = ρ(E)) That is E = (I − ωD−1L) (Jacobi’s Method), E = I − PℓL−1

ℓ+1P T ℓ Lℓ

(Coarse Grid Correction) or Eℓ = Mν2

ℓ (I − Pℓ(I − Eγ ℓ+1)(Lℓ+1)−1RℓLℓ)Mν1 ℓ

Eℓmax = 0 (Multigrid Method).

Stencil Operators

◮ Let Ωh = h1Z × · · · × hdZ. ◮ If L is a linear operator on the grid Ωh it

can be written in the form (Lu)(x) =

• κ∈hZd

sx(κ) · u(x + κ) for x ∈ Ωh.

◮ For many operators:

x ≈ x′ ⇒ sx ≈ sx′.

◮ To analyze the local behavior it is

sufficient to deal with constant stencils.

u(x) u(x + κ) sx(κ)

Local Fourier Analysis

• 1. We seek for E for “some” operator E and
• 2. for ρ(E).

(since limn→∞ En1/n = ρ(E))

◮ Consider the operator given by the constant stencil

(Lu)(x) =

• κ∈hZd

s(κ) · u(x + κ)

◮ Grid size independent analysis.

◮ “Move boundary conditions to infinity.”

⇒ Consider operators on an infinite grid.

◮ Exploits simple structure in the frequency Domain. ◮ Quantitative estimates.

Discrete Time Fourier Transform (DTFT)

“Time” Domain Frequency Domain Fh F −1

h

• 10
• 5

5 10 0.2 0.4 0.6 0.8 1 (truncated) time domain

• 3
• 2
• 1

1 2 3 0.2 0.4 0.6 0.8 1 frequency domain vol1/2

h

(2π)d/2

• x∈Ωh

f(x)eiϑ,x

vol1/2

h

(2π)d/2

• Θh

ˆ u(ϑ)eiϑ,x dϑ

Discrete Time Fourier Transform (DTFT)

“Time” Domain Frequency Domain Fh F −1

h

• 10
• 5

5 10 0.2 0.4 0.6 0.8 1 (truncated) time domain

• 3
• 2
• 1

1 2 3 0.2 0.4 0.6 0.8 1 frequency domain vol1/2

h

(2π)d/2

• x∈Ωh

f(x)eiϑ,x

vol1/2

h

(2π)d/2

• Θh

ˆ u(ϑ)eiϑ,x dϑ

Discrete Time Fourier Transform (DTFT)

“Time” Domain Frequency Domain Fh F −1

h

• 10
• 5

5 10 0.2 0.4 0.6 0.8 1 (truncated) time domain

• 3
• 2
• 1

1 2 3 0.2 0.4 0.6 0.8 1 frequency domain vol1/2

h

(2π)d/2

• x∈Ωh

f(x)eiϑ,x

vol1/2

h

(2π)d/2

• Θh

ˆ u(ϑ)eiϑ,x dϑ

Discrete Time Fourier Transform (DTFT)

“Time” Domain Frequency Domain Fh F −1

h

• 10
• 5

5 10 0.2 0.4 0.6 0.8 1 (truncated) time domain

• 3
• 2
• 1

1 2 3 0.2 0.4 0.6 0.8 1 frequency domain vol1/2

h

(2π)d/2

• x∈Ωh

f(x)eiϑ,x

vol1/2

h

(2π)d/2

• Θh

ˆ u(ϑ)eiϑ,x dϑ

Discrete Time Fourier Transform (DTFT)

“Time” Domain Frequency Domain Fh F −1

h

• 10
• 5

5 10 0.2 0.4 0.6 0.8 1 (truncated) time domain

• 3
• 2
• 1

1 2 3 0.2 0.4 0.6 0.8 1 frequency domain vol1/2

h

(2π)d/2

• x∈Ωh

f(x)eiϑ,x

vol1/2

h

(2π)d/2

• Θh

ˆ u(ϑ)eiϑ,x dϑ

Local Fourier Analysis

Stencil operators

• Lu
• (x) =
• κ∈hZd

s(κ) · u(x + κ) . exhibit a simple representation under the DTFT. u f ˆ u ˆ f L ˆ f = L · ˆ u Fh Fh Theorem Let L ∈ L∞(Θh) be the symbol of L. Then Lℓ2 = ess sup

ϑ

| L(ϑ)| ρ(L) = ess sup

ϑ

ρ( L(ϑ))

1D Poisson Smoothing

• E(ϑ) = cos(hϑ)
• 3
• 2
• 1

1 2 3 0.2 0.4 0.6 0.8 1

1D Poi. Eq.

∂2u ∂x2 = f

Sm. GS 0.50 RB-GS 0.25 Jac. 1.00

1D Poisson Smoothing

• E(ϑ) = exp(−ihϑ)

2−exp(ihϑ)

• 3
• 2
• 1

1 2 3 0.2 0.4 0.6 0.8 1

1D Poi. Eq.

∂2u ∂x2 = f

Sm. GS 0.50 RB-GS 0.25 Jac. 1.00

Harmonic Frequencies

◮ Red-Black Gauß-Seidel is not shift-invariant. ◮ Low and High Frequencies interact.

u1 u2 u3 u4 u5 u6 u1 u2 u3 u4 u5 u0 ϑlow ∈ [− π

2h, π 2h)

ϑhigh = [ϑlow + 2π/h]h

• 4
• 2

2 4

• 1
• 0.5

0.5 1

1D Poisson Smoothing

1 2

• A(ϑ)+1 −A(ϑ+π)+1

−A(ϑ)+1 A(ϑ+π)+1

• · 1

2

• A(ϑ)+1 A(ϑ+π)−1

A(ϑ)−1 A(ϑ+π)+1

• 1.5
• 1
• 0.5

0.5 1 1.5 0.2 0.4 0.6 0.8 1 low → low

• 1.5
• 1
• 0.5

0.5 1 1.5 0.2 0.4 0.6 0.8 1 high → low

• 1.5
• 1
• 0.5

0.5 1 1.5 0.2 0.4 0.6 0.8 1 low → high

• 1.5
• 1
• 0.5

0.5 1 1.5 0.2 0.4 0.6 0.8 1 high → high

1D Poi. Eq.

∂2u ∂x2 = f

Sm. GS 0.50 RB-GS 0.25 Jac. 1.00

Harmonic Freq. / Block Toeplitz Operators

Rˆ u = (R1ˆ u, . . . , Rnˆ u)T = (ˆ u (ϑ + σ1) , . . . , ˆ u (ϑ + σn))T {σ1, . . . , σn} = { 2π

hS z | 0 ≤ zj < S, z ∈ Zd}

ˆ u

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 1 2 3 4 5 6

Rˆ u

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

Harmonic Freq. / Block Toeplitz Operators

Rˆ u = (R1ˆ u, . . . , Rnˆ u)T = (ˆ u (ϑ + σ1) , . . . , ˆ u (ϑ + σn))T {σ1, . . . , σn} = { 2π

hS z | 0 ≤ zj < S, z ∈ Zd}

u f ˆ u ˆ f

• ˆ

u

• ˆ

f L

• ˆ

f = L · ˆ u Fh R Fh R Lℓ2 = ess sup

ϑ

• L(ϑ)2

ρ(L) = ess sup

ϑ

ρ( L(ϑ)) ˆ u

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 1 2 3 4 5 6

Rˆ u

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

Block Stencil Operators

Block Stencil Operators

Block Stencil Operators

Block Stencil Operators

• L(ϑ) = F ∗(G(ϑ) ◦ F)

where Gkj = Rj L(k) Fkj =

1 √neiσj,τk

• is the Haramard product

u f ˆ u ˆ f

• ˆ

u

• ˆ

f L

• ˆ

f = L · ˆ u Fh R Fh R

◮ Symbol of a block stencil

• perator can be

computed from the symbols of its stencils.

Analysis of General Operators

◮ How to analyze more general operators like

Jacobi’s Method E = (I − ωD−1L) ?

◮ Use linearity

• A + λB =

A + λ B and

◮ the composition of operators

BA = B A .

◮ And in case of different block sizes:

ΠT  

• L(ϑ+µ1)

...

• L(ϑ+µm)

  Π′

◮ How about a whole Multigrid method?

Intergrid Transfer Operators

Restriction R = RinjRst with (Rinju)(x) = u(x)

• Rinj =

1 √n(1, . . . , 1)

Interpolation P = PstPinj with (Pinju)(x) =

• u(x)

if x ∈ ΩhS

• therwise
• Pinj =

1 √n(1, . . . , 1)T

ˆ u

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 1 2 3 4 5 6

Rˆ u

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

• 2
• 1.5
• 1
• 0.5

0.5 1 1.5 2 0.75 1.25

A Jumping Coefficient Problem

−div(a grad(u)) = f

• n Ω

u = g

• n ∂Ω

◮ a : Ω → R ◮ Finite Volume Discretization ◮ a discontinuous at interfaces

a = 106 a = 1

u(x) 5 10 15 20 25 30 35 X 5 10 15 20 25 30 35 Y 2e-05 4e-05 6e-05 8e-05 0.0001 0.00012 u(x) 5e-05 1e-05 1e-06 1e-07 1e-09 1e-11 5 10 15 20 25 30 35 X 5 10 15 20 25 30 35

Multigrid Solver

−div(a grad(u)) = f

• n Ω

u = g

• n ∂Ω

◮ a discontinuous at interfaces

a = 106 a = 1

◮ Pointwise smoother (ω-Jacobi, Gauß-Seidel).

Mℓ = (I − B−1

ℓ Lℓ) ◮ Galerkin coarse grid approximation.

Lℓ+1 = P T

ℓ LℓPℓ ◮ Operator dependent interpolation (see [ABDP81]). ◮ Full Coarsening (Two-Level).

Analysis of the JCP

◮ LFA analyzes local behavior.

⇒ Capture the essentials of the operator by a local (and repeated) block stencil operator.

◮ Repeated jumps yield the same behavior as the original

• perator.

. . . · · · · · · . . .

Analysis of the JCP

−div(a grad(u)) = f

• n Ω

u = g

• n ∂Ω

◮ a discontinuous at interfaces

a = 106 a = 1

Jacobi’s Method E = (I − ωD−1L)

Jac (0.8), Poisson 5 10 15 20 25 30 5 10 15 20 25 30 0.2 0.4 0.6 0.8 1 Jac (0.8), JCP (4x4) 5 10 15 20 25 30 5 10 15 20 25 30 0.2 0.4 0.6 0.8 1

Analysis of the JCP

−div(a grad(u)) = f

• n Ω

u = g

• n ∂Ω

◮ a discontinuous at interfaces

a = 106 a = 1

Coarse Grid Correction E = I − PℓL−1

ℓ+1P T ℓ Lℓ

GCA (Bilinear), JCP (4x4) 5 10 15 20 25 30 5 10 15 20 25 30 0.5 1 1.5 2 2.5 GCA (Adaptive), JCP (4x4) 5 10 15 20 25 30 5 10 15 20 25 30 0.5 1 1.5 2 2.5

Analysis of the JCP

−div(a grad(u)) = f

• n Ω

u = g

• n ∂Ω

◮ a discontinuous at interfaces

a = 106 a = 1

Two-Grid Method E = Mℓ(I − PℓL−1

ℓ+1P T ℓ Lℓ)Mℓ

Jac (0.8), Dep., GCA, JCP (4x4) 5 10 15 20 25 30 5 10 15 20 25 30 0.2 0.4 0.6 0.8 1

Summary and Outlook

◮ To maximize the usage of your hardware choose from different

algorithms. ⇒ Need for (automatic) analysis of the convergence behavior.

◮ LFA can be used to analyze complex problems.

◮ Write all your operators in block stencil form. ◮ Use the formula

L = F ∗(G ◦ F) to compute its frequency representation.

◮ Combine the operators (linearity and composition). ◮ Easily compute the norm or spectral radius.

◮ A Software tool is in preparation for easy application of the

analysis.

Bibliography

• R. Alcouffe, A. Brandt, J. Dendy, Jr., and J. Painter.

The multi-grid method for the diffusion equation with strongly discontinuous coefficients. SIAM Journal on Scientific and Statistical Computing, 2(4):430–454, 1981.