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 = ∂ 2 u ∂x 2 + ∂ 2 u ∆ u = f on Ω ∂y 2 where u = g on ∂ Ω Ω = (0 , 1) × (0 , 1)

Poisson’s Equation ∆ u = ∂ 2 u ∂x 2 + ∂ 2 u ∆ u = f on Ω ∂y 2 where u = g on ∂ Ω Ω = (0 , 1) × (0 , 1) Photo by Graham Richardson / CC BY

Poisson’s Equation ∆ u = ∂ 2 u ∂x 2 + ∂ 2 u ∆ u = f on Ω ∂y 2 where u = g on ∂ Ω Ω = (0 , 1) × (0 , 1) 0 -0.5 -1 1 1 0.5 0.5 0 0

Poisson’s Equation ∆ u = ∂ 2 u ∂x 2 + ∂ 2 u ∆ u = f on Ω ∂y 2 where u = g on ∂ Ω Ω = (0 , 1) × (0 , 1) 0 -0.5 -1 1 1 0.5 0.5 0 0

Poisson’s Equation ∆ u = ∂ 2 u ∂x 2 + ∂ 2 u ∆ u = f on Ω ∂y 2 where u = g on ∂ Ω Ω = (0 , 1) × (0 , 1) u : Ω → R y Ω x

Poisson’s Equation ∆ u = ∂ 2 u ∂x 2 + ∂ 2 u ∆ u = f on Ω ∂y 2 where u = g on ∂ Ω Ω = (0 , 1) × (0 , 1) u : Ω → R u h : Ω h → R y y Ω h Ω x x

Multigrid in One Slide ◮ For many simple iterative methods and certain problems we have that the error is 0.2 0.15 0.1 0.05 smooth after a few iterations. 0 -0.05 -0.1 -0.15 -0.2 ◮ Smooth error can be represented on a 1 0.8 1 0.6 0.8 0.4 0.6 0.4 0.2 0.2 0 0 coarser grid. ⇒ Linear system of reduced size. 0.015 0.01 ◮ Use recursion on the coarse problem. 0.005 0 -0.005 -0.01 choose u 1 0.8 0.2 0.4 0.6 0.8 1 0.6 0.4 0.2 0 repeat 0 smooth u 0.015 e ← coarse approx. of e ˜ 0.01 0.005 0 u ← u + ˜ e -0.005 -0.01 1 smooth u 0.8 0.2 0.4 0.6 0.8 1 0.6 0.4 0.2 0 0 until satisfied

Gauß-Seidel Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Bad Parallelization...

Gauß-Seidel Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Bad Parallelization...

Gauß-Seidel Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Bad Parallelization...

Gauß-Seidel Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Bad Parallelization...

Gauß-Seidel Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Bad Parallelization...

Gauß-Seidel Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Bad Parallelization...

Red-Black GS Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 2 Msgs.

Red-Black GS Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 2 Msgs.

Red-Black GS Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 2 Msgs.

Red-Black GS Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 2 Msgs.

Red-Black GS Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 2 Msgs.

Red-Black GS Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 2 Msgs.

Red-Black GS Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 2 Msgs.

Red-Black GS Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 2 Msgs.

Red-Black GS Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 2 Msgs.

Red-Black GS Smoother u x,y = 1 4 ( f x,y − u x,y − 1 − u x − 1 ,y − u x +1 ,y − u x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 2 Msgs.

Jacobi Smoother u n ew x,y = 1 4 ( f x,y − u o ld x,y − 1 − u o ld x − 1 ,y − u o ld x +1 ,y − u o ld x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 1 Msgs.

Jacobi Smoother u n ew x,y = 1 4 ( f x,y − u o ld x,y − 1 − u o ld x − 1 ,y − u o ld x +1 ,y − u o ld x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 1 Msgs.

Jacobi Smoother u n ew x,y = 1 4 ( f x,y − u o ld x,y − 1 − u o ld x − 1 ,y − u o ld x +1 ,y − u o ld x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 1 Msgs.

Jacobi Smoother u n ew x,y = 1 4 ( f x,y − u o ld x,y − 1 − u o ld x − 1 ,y − u o ld x +1 ,y − u o ld x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 1 Msgs.

Jacobi Smoother u n ew x,y = 1 4 ( f x,y − u o ld x,y − 1 − u o ld x − 1 ,y − u o ld x +1 ,y − u o ld x,y +1 ) ◮ Parallelizable ◮ Per Iteration: √ √ ◮ ≈ 4 N/ P Words p. Proc. ◮ 1 Msgs.

Jacobi Smoother u n ew x,y = 1 4 ( f x,y − u o ld x,y − 1 − u o ld x − 1 ,y − u o ld x +1 ,y − u o ld 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 (since lim n →∞ � E n � 1 /n = ρ ( E ) ) 2. for ρ ( E ) . That is E = ( I − ωD − 1 L ) (Jacobi’s Method), E = I − P ℓ L − 1 ℓ +1 P T ℓ L ℓ (Coarse Grid Correction) or ℓ ( I − P ℓ ( I − E γ E ℓ = M ν 2 ℓ +1 )( L ℓ +1 ) − 1 R ℓ L ℓ ) M ν 1 ℓ E ℓ max = 0 (Multigrid Method).

Local Fourier Analysis 1. We seek for � E � for “some” operator E and (since lim n →∞ � E n � 1 /n = ρ ( E ) ) 2. for ρ ( E ) . That is E = ( I − ωD − 1 L ) (Jacobi’s Method), E = I − P ℓ L − 1 ℓ +1 P T ℓ L ℓ (Coarse Grid Correction) or ℓ ( I − P ℓ ( I − E γ E ℓ = M ν 2 ℓ +1 )( L ℓ +1 ) − 1 R ℓ L ℓ ) M ν 1 ℓ E ℓ max = 0 (Multigrid Method).

Stencil Operators ◮ Let Ω h = h 1 Z × · · · × h d Z . ◮ If L is a linear operator on the grid Ω h it can be written in the form � ( Lu )( x ) = s x ( κ ) · u ( x + κ ) κ ∈ h Z d for x ∈ Ω h . ◮ For many operators: s x ( κ ) x ≈ x ′ ⇒ s x ≈ s x ′ . u ( x ) u ( x + κ ) ◮ To analyze the local behavior it is sufficient to deal with constant stencils.

Local Fourier Analysis 1. We seek for � E � for “some” operator E and (since lim n →∞ � E n � 1 /n = ρ ( E ) ) 2. for ρ ( E ) . ◮ Consider the operator given by the constant stencil � ( Lu )( x ) = s ( κ ) · u ( x + κ ) κ ∈ h Z d ◮ 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) F h “Time” Domain Frequency Domain F − 1 h (truncated) time domain frequency domain 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 -10 -5 0 5 10 -3 -2 -1 0 1 2 3 � � vol 1 / 2 vol 1 / 2 u ( ϑ ) e i � ϑ,x � d ϑ f ( x ) e i � ϑ,x � h h ˆ (2 π ) d/ 2 (2 π ) d/ 2 Θ h x ∈ Ω h

Discrete Time Fourier Transform (DTFT) F h “Time” Domain Frequency Domain F − 1 h (truncated) time domain frequency domain 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 -10 -5 0 5 10 -3 -2 -1 0 1 2 3 � � vol 1 / 2 vol 1 / 2 u ( ϑ ) e i � ϑ,x � d ϑ f ( x ) e i � ϑ,x � h h ˆ (2 π ) d/ 2 (2 π ) d/ 2 Θ h x ∈ Ω h

Discrete Time Fourier Transform (DTFT) F h “Time” Domain Frequency Domain F − 1 h (truncated) time domain frequency domain 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 -10 -5 0 5 10 -3 -2 -1 0 1 2 3 � � vol 1 / 2 vol 1 / 2 u ( ϑ ) e i � ϑ,x � d ϑ f ( x ) e i � ϑ,x � h h ˆ (2 π ) d/ 2 (2 π ) d/ 2 Θ h x ∈ Ω h