CS 294-73 Software Engineering for Scientific Computing - - PowerPoint PPT Presentation

CS 294-73 
 Software Engineering for Scientific Computing
 
 Lecture 16: Multigrid (structured grids revisited)
 


Guest lecturer: Hans Johansen, hjohansen@lbl.gov

10/17/2017 CS294-73 – Lecture 16

Laplacian, Poisson’s equation, Heat Equation

Laplacian: Poisson’s equation: Heat equation:

2

Blue:Ω Black: ∂Ω

10/17/2017 CS294-73 – Lecture 16

Discretizing on a rectangle (or a union of rectangles)

Defined only on interior (blue)

• nodes. Values at black nodes

are set to zero if we are using the boundary conditions given above.

3

10/17/2017 CS294-73 – Lecture 16

Solving Poisson’s equation

Want to solve linear system of equations

4

For ϕ defined on the interior (blue) grid points i . This makes sense, since the stencil for the operator requires only nearest neighbors. and we have values defined by the boundary conditions to be set to zero (black points).

10/17/2017 CS294-73 – Lecture 16

Point Jacobi

As in the unstructured grid case, we can use point Jacobi to solve these equations.

5

where is the relaxation parameter for our iterative scheme. Each iteration corrects eigenmode of at different rate:

λ ∆hφh

0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

N=16, k=1, 16 iterations … operations to get this far …

O(N 2)

(ref: Briggs, et al, A Multigrid Tutorial)

kth

10/17/2017 CS294-73 – Lecture 16

Smoothing properties of point Jacobi

Define the error

6

We can also see how the error behaves under point Jacobi Even though we don’t know the error, we can compute the residual to provide a measure for the error: e.g. convergence if .

10/17/2017 CS294-73 – Lecture 16

Smoothing properties of point Jacobi

For example, choose . Then

7

10/17/2017 CS294-73 – Lecture 16

Smoothing properties of point Jacobi

8

• The value at the new iteration is an average of the values at the
• ld iteration: weighted sum, positive weights that sum to 1.
• The max and min values are always strictly decreasing.

Smoothes the local error very rapidly.

• Error equation looks like forward Euler for the heat equation –

smooths higher wavenumbers faster.

10/17/2017 CS294-73 – Lecture 16

Point Jacobi - smoothing example

9

0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

N=128, k=1, for 128 iterations: 64x the work, but worse error! What about other modes?

0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

Point Jacobi for Poisson requires work for a specified level of error! (and leaves “smooth” error) Why? Discrete Fourier analysis …

O(N D+2)

k=2,4,8 for 128 iterations: each closer to final solution.

10/17/2017 CS294-73 – Lecture 16

Discrete Fourier Analysis of Point Jacobi

10

Assume periodic bc’s for Poisson, and recall discrete Fourier modes: And finite difference operators are diagonalized:

10/17/2017 CS294-73 – Lecture 16

Discrete Fourier Analysis of Point Jacobi

11

So what does Point Jacobi do in Fourier space?

δl+1 = δl + λ∆hδl ⇒ Fk(δl+1) = Fk(δl) + λ Λ(k)Fk(δl) = (1 + λ Λ(k)) Fk(δl) = ω(k) Fk(δl)

Where each error mode’s damping factor is:

ω(k) = 1 + λΛ(k) = 1 + λ h2

• zk − 2 + z−k

= 1 + 2λ h2 (cos(2πkh) − 1) ω(k) = 1

10/17/2017 CS294-73 – Lecture 16

Discrete Fourier Analysis of Point Jacobi

12

What does this look like for (real spectrum) ?

β = 2πkh , k ∈ [0, N 2 ]

0.5 1 1.5 2 2.5 3 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

λ = 1/2

+/- modes undamped!

ω(β) = 1 + 2λ (cos(β) − 1) λ = 1/4

highest k modes completely damped No choice helps lowest k modes

λ = 1/8

meh

λ = 1/3

most high-k modes damped

10/17/2017 CS294-73 – Lecture 16

Point Jacobi Fourier Analysis: Conclusions

13

• A few iterations of point Jacobi can be tuned to reduce the

“high frequency” error with very just a few sweeps.

• These modes can be fixed locally – across a small number of

points – while the “smooth” error persists globally.

• As resolution decreases, N increases, and smooth errors ( j=1)

are less and less responsive to point Jacobi:

ω(1) = 1 + 2λ h2 (cos(2πh) − 1) max(λ) = h2 2 ≈ 1 − Ch2

• Ideally, we’d use PJ for those modes separately with a coarse

and with good convergence rate and is stable on coarse grid …

• Then somehow blend the answers?

Multigrid!

λc Nc, hc

10/17/2017 CS294-73 – Lecture 16

Multigrid

14

• Do a few iterations of point Jacobi on

your grid to obtain

• Compute the residual .
• Average the residual down onto a grid

coarsened by a factor of 2 :

• Apply point Jacobi to solve the residual

equation for :

• Interpolate correction back onto fine grid

solution:

• Smooth again using point Jacobi.

10/17/2017 CS294-73 – Lecture 16

Multigrid

15

• If the number of grid points is 2M+1, can

apply this recursively.

multigrid(phi,f,h,level) phi:= phi + lambda*(L(phi) – f); (p times) if (level > 0) R = L(phi) – f; Rcoarse = Av(R); delta = 0.; call multigrid(delta,Rcoarse,2h, level-1); phi += Int(delta); endif; phi:= phi + lambda*(L(phi) – f); (p times)

10/17/2017 CS294-73 – Lecture 16

Averaging and Interpolation

16

• Conservation of total charge.
• Adjoint & order conditions, other considerations

(see Trottenberg, et al, Multigrid).

• For our second-order accurate discretization on a

nodal-point grid, we can use the trapezoidal rule (and bilinear interpolation). Even if the grid is not a power of 2, can apply a direct solve at the bottom. 3 levels in 3D leads to a reduction by 512 in the number of unknowns.

¹ ¹ ¹ ¹ ² ² ² ²

⁴

¹⁄₁₆ ×

1

½ ½ ½ ½ ¼ ¼ ¼ ¼

A I

10/17/2017 CS294-73 – Lecture 16

Multigrid Convergence: 1D example

17

0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

N=128, exact solution

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −30 −20 −10 10 20 30

N=128, RHS step functions

10/17/2017 CS294-73 – Lecture 16

Multigrid Convergence: 1D example

18

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

N=128, Point Jacobi only

Iterations Max Error N ~0.9 8 * N ~0.45 N^2 ~2e-6

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

Coarse solution only

Depth Max Error 1 ~9e-4 2 ~4e-3 3 ~2e-2 4 ~6e-2 5 ~0.25

10/17/2017 CS294-73 – Lecture 16

Multigrid Convergence: 1D example

19

What’s the error from a coarse grid solution? Look at Fourier modes …

• Solution has lots of high

wave-number content

10 20 30 40 50 60 10

−2

10

−1

10 10

1

10

2

10/17/2017 CS294-73 – Lecture 16

Multigrid Convergence: 1D example

20

What’s the error from a coarse grid solution? Look at Fourier modes …

• Solution has high-k

(wave number) content

• Solution on 1st coarse

level matches low-k well, but has high-k error (piecewise linear interp.)

10 20 30 40 50 60 10

−2

10

−1

10 10

1

10

2

10/17/2017 CS294-73 – Lecture 16

Multigrid Convergence: 1D example

21

What’s the error from a coarse grid solution? Look at Fourier modes …

10 20 30 40 50 60 10

−2

10

−1

10 10

1

10

2

• Solution has high-k

(wave number) content

• Solution on 1st coarse

level matches low-k well, but has high-k error (piecewise linear interp.)

• Solution on 2nd coarse

level has mid-k error too (PJ less effective)

10/17/2017 CS294-73 – Lecture 16

Multigrid Convergence: 1D example

22

What’s the error from a coarse grid solution? Look at Fourier modes …

• Solution has high-k

(wave number) content

• Solution on 1st coarse

level matches low-k well, but has high-k error (piecewise linear interp.)

• Solution on 2nd coarse

level has mid-k error too (PJ less effective)

• Solution on 3rd coarse

level has errors in most k (solving wrong problem)

10 20 30 40 50 60 10

−2

10

−1

10 10

1

10

2

10/17/2017 CS294-73 – Lecture 16

10 20 30 40 50 60 10

−2

10

−1

10 10

1

10

2

Multigrid Convergence: 1D example

23

What’s the error from a coarse grid solution? Look at Fourier modes …

• Solution has high-k

(wave number) content

• Solution on 1st coarse

level matches low-k well, but has high-k error (piecewise linear interp.)

• Solution on 2nd coarse

level has mid-k error too (PJ less effective)

• Solution on 3rd coarse

level has errors in most k (solving wrong problem)

• MG V-cycle gets all-k

10/17/2017 CS294-73 – Lecture 16

Multigrid Convergence: 1D example

24

1 2 3 4 5 6 7 8 9 10 10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10

How much work was it?

• O(N) PJ on each level, O(N) cost to average/interpolate à O(N)
• Each coarser level is 1/2 the work à O(N/2) + O(N/4) + O(N/8) …

à Overall O(N) – compared to O(N3) for point Jacobi alone!

• NOTE: ideal convergence for model problem, real world is harder!

Error vs. MG iteration

N=128 *and* 2048

• N=4 coarsest grid
• 2 Point Jacobi before

(to reduce aliasing)

• 2 Point Jacobi after

(to reduce high-k error) à Each MG V-cycle reduces error by ~12x

A I

10/17/2017 CS294-73 – Lecture 16

Multigrid Analysis: Conclusions for 1D example

25

• Coarse grid correction is close to fine solution for low-k, but has

high-k content from interpolation operator.

• Coarsest point relaxation reduces low-k wave number errors.
• The recursive algorithm reduces all wave numbers in O(N) ops:

(relax – coarsen – (do coarse) – interpolate – relax) MG approach is problem-dependent … 1.2-20x reductions typical:

• Choice of point relaxation: some more effective than others,

depending on the operator (discontinuous coefficients, etc.).

• Average and interpolation operators should represent smooth

components of fine residuals/coarse solutions well (accuracy), with errors in the “near null space” of point relaxation scheme.

• See Trottenberg, et al, Multigrid for theory, unfriendly examples.

10/17/2017 CS294-73 – Lecture 16

Results for real world examples

26

10/17/2017 CS294-73 – Lecture 16

Finite-Volume Methods on Structured Grids

Using the Divergence Theorem, this leads to a natural discretization of the Laplacian.

This discretization satisfies a discrete analogue of the divergence theorem.

We write the Laplacian as the divergence of fluxes.

27

10/17/2017 CS294-73 – Lecture 16

Finite-Volume Methods on Structured Grids

We use centered-differences and the midpoint rule to approximate fluxes: Away from boundaries, this leads to the same finite-difference approximation to the Laplacian. At boundaries, it is different. It also interacts differently with multigrid.

28

10/17/2017 CS294-73 – Lecture 16

Finite-Volume Methods on Structured Grids

Coarsening is done by aggregating volumes, rather than by sampling gridpoints. Boundary conditions are expressed in terms of computing fluxes at boundary faces.

29

10/17/2017 CS294-73 – Lecture 16

Finite-Volume Methods on Structured Grids

Averaging and interpolation: “Charge conservation” motivated by PDE Piecewise constant interpolation!

30

10/17/2017 CS294-73 – Lecture 16

Finite-Volume Methods on Structured Grids

Why would this work? 1D example …

31

hφciic−1 hφciic+1 hφciic ⌦ φf↵

i

⌦ φf↵

i−1

⌦ φf↵

i+1

Piecewise constant is linear with small +/- mode, which PJ removes

10/17/2017 CS294-73 – Lecture 16

Discrete Fourier Analysis: 2D Poisson example

N=322 initial error Fourier modes 1 MG V-cycle 3 MG V-cycle

• O(N) PJ + average/interpolate on each level
• Coarser levels O(N/4) + O(N/16) + O(N/64) + …

work, better in 2D, even better in 3D

• Choosing , , and relaxation scheme are critical
• Note: need to shift FFT by for FV

A I

32

e2π h

2

10/17/2017 CS294-73 – Lecture 16

Results (Fourth-Order Finite-Volume)

33

10/17/2017 CS294-73 – Lecture 16

The Heat Equation

Explicit discretization in time leads to time step restrictions of the form So instead, we use implicit discretizations, e.g. Backward Euler: More generally, any implicit method for the heat equation will involve solving

34

L(φ) = f, L = (I − µ∆h), µ > 0

10/17/2017 CS294-73 – Lecture 16

The Heat Equation – Point Jacobi

Can do the same residual / error analysis as in Poisson to obtain

35

φl+1 = φl + λ

• f − (I − µ∆)φl

δl

i ≡ φl i − φi

(Lδ)l

i = Rl i , Rl = Lφl − f

δl+1 = δl − λ Lδl

Discrete Fourier analysis leads to:

ω(k) = 1 − λ + 2µλ h2 (cos(2πkh) − 1)

For large (big or strong diffusion) looks like Poisson, but for small this looks more like identity operator (converges much faster for low-k):

µ ∆t µ ω(k) ≈ 1 − λ

10/17/2017 CS294-73 – Lecture 16

The Heat Equation – Point Jacobi

Weighted sum, positive weights, but sum to less than 1, more so as gets smaller … heat equation converges faster than Poisson, sometimes very fast.

36

Relaxation parameter that eliminates highest-k wave numbers: Note as h à 0 this becomes smaller (more Poisson-like), but ~1 for More generally for any D:

µ ω(N/2) = 1 − λ + 2µλ h2 (cos(π) − 1) ⇒ λ = 1 1 + 4µ

h2

µ ⌧ h2 4 λ = 1 1 + 4Dµ

h2

δl+1

i

= λ ✓2Dµ h2 ◆ δl

i + 1

2D X

d

δl

i−ed + δl i+ed

!