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)
 


10/24/2019 CS294-73 – Lecture 16

Laplacian, Poisson’s equation, Heat Equation

Laplacian: Poisson’s equation: Heat equation:

2

Blue:Ω Black: ∂Ω

10/24/2019 CS294-73 – Lecture 16

FFT1DW problem

void FFT1DW::forwardFFTCC(vector<complex<double> > & a_fHat, const vector<complex<double> >& a_f) const { for (int i = 0; i < a_f.size(); i++) { m_in[i] = a_f[i]; } fftw_execute(m_forward); ... } Compiler error. Fixed by declaring in .H file ... mutable vector<complex<double > m_in;

3

10/24/2019 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.

4

10/24/2019 CS294-73 – Lecture 16

Solving Poisson’s equation

Want to solve linear system of equations

5

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/24/2019 CS294-73 – Lecture 16

Point Jacobi

We have used point Jacobi to solve these equations.

6

10/24/2019 CS294-73 – Lecture 16

Smoothing properties of point Jacobi

Define the error

7

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/24/2019 CS294-73 – Lecture 16

Smoothing properties of point Jacobi

For example, choose . Then

8

10/24/2019 CS294-73 – Lecture 16

Smoothing properties of point Jacobi

9

• 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.

Smooths the local error very rapidly.

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

smooths higher wavenumbers faster.

10/24/2019 CS294-73 – Lecture 16

10

Let f = W(k) . Then . Point Jacobi approximates this by For the error to converge to zero for any Fourier mode, we require that Generally, high wave numbers in he error are damped efficiently, but not others. But this is relative to the grid spacing.

0 ≤ 1 + λσk < 1

Discrete Fourier Analysis of Point Jacobi

∆hW (k) = 1 h2 (−4 + 2cos(2πkxh) + 2cos(2πkyh))W (k) ≡ σkW (k) −(∆h)−1W (k) = − 1 σk W (k) φl = (1 + λ(σk + 1))φl−1 δl = (1 + λσk)δl−1 = (1 + λσk)lδ0 k = ⇣N 2 , N 2 ⌘ :σk = −8/h2 , 1 + λσk = 0 |k|h << 1 :λσk = O(1) , 1 + λσk = 1 + O(h2)

10/24/2019 CS294-73 – Lecture 16

11

Idea:

• 1. Apply point Jacobi to the current approximation to the solution.
• 2. Solve a coarse-grid version of the problem on a coarsened grid.
• 3. Correct the fine-grid solution using the coarse-grid solution.

Why should this work ? Point Jacobi eliminates the high-wavenumber components of the solution so that what remains is representable on the coarse grid. How does this work ?

• Residual-correction form:

captures the smooth remnants of the error.

• Have to iterate, since this isn’t exact.

Deep mathematics: point Jacobi smooths <-> special property of Laplacian.

Multigrid

∆hδ = Rh , Rh = ∆h ˜ φ − f, φ = ˜ φ + δ ⇓ ∆hφ = f

10/24/2019 CS294-73 – Lecture 16

Multigrid

12

• 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/24/2019 CS294-73 – Lecture 16

Multigrid

13

• 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/24/2019 CS294-73 – Lecture 16

Averaging and Interpolation

14

• Conservation of total charge.
• Adjoint conditions, order conditions (see

Trottenberg). 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.

¹ ¹ ¹ ¹ ² ² ² ² ⁴ ¹⁄₁₆× ₁ ½ ½ ½ ½ ¼ ¼ ¼ ¼

10/24/2019 CS294-73 – Lecture 16

Results

15

10/24/2019 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.

10/24/2019 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.

10/24/2019 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.

10/24/2019 CS294-73 – Lecture 16

Finite-Volume Methods on Structured Grids

Averaging and interpolation:

10/24/2019 CS294-73 – Lecture 16

Results (Fourth-Order Finite-Volume)

10/24/2019 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

10/24/2019 CS294-73 – Lecture 16

The Heat Equation – Point Jacobi

Can do the same residual / error analysis as in Poisson to obtain Weighted sum, positive weights, but sum to less than 1, more so as h gets larger. Fourier analysis works the same way as Poisson.

10/24/2019 CS294-73 – Lecture 16

Multigrid Pseudocode

23

At the top level, iterate until residual is reduced by some factor.

vcycle(φ, ρ) { φ := φ + λ(L(φ) − ρ) p times if (level > 0) { R = ρ − L(φ) Rc = A(R) δ : Bc → R , δ = 0 vcycle(δ, Rc) φ := φ + I(δ) φ := φ + λ ∗ (L(φ) − ρ) p times } else { φ := φ + λ ∗ (L(φ) − ρ) pB times { }

10/24/2019 CS294-73 – Lecture 16

Implementing Multigrid

24

The algorithm is naturally recursive – implement as nested set of calls to

vcycle.

10/24/2019 CS294-73 – Lecture 16

main

25

... Multigrid mg(domain,dx,numLevels); ... int maxiter; double tol; cin >> maxiter >> tol; double resnorm0 = mg.resnorm(phi,rho); cout << "initial residual = " << resnorm0 << endl; for (int iter = 0; iter < maxiter; iter++) { mg.vCycle(phi,rho); double resnorm = mg.resnorm(phi,rho); cout << "iter = " << iter << ", resnorm = " << resnorm << endl; if (resnorm < tol*resnorm0) break; }

10/24/2019 CS294-73 – Lecture 16

Multigrid.H

26

... private: BoxData<double > m_res; //Residual at current level. BoxData<double > m_resc; // Average of residual. BoxData<double > m_delta; // Coarse-level correction. Box m_box; Box m_bdry[2*DIM]; int m_domainSize; shared_ptr<Multigrid> m_coarsePtr; double m_dx; double m_lambda; int m_level; // What level ? (m_level = 0 is the bottom). int m_preRelax = 2*DIM; int m_postRelax = 2*DIM; int m_bottomRelax = 10; int m_flops = 0; // I’m going to use this to count flops. ...

10/24/2019 CS294-73 – Lecture 16

Multigrid::define

27

... Multigrid::define(const Box& a_box, double a_dx, int a_level){ m_box = a_box; m_level = a_level; m_res.define(m_box); m_dx = a_dx; m_domainSize = m_box.low()[0] + 1; if (m_level > 0) { Box bCoarse = m_box.coarsen(2); m_resc.define(bCoarse); m_delta.define(bCoarse.grow(1)); m_coarsePtr = shared_ptr<Multigrid> 
 (new Multigrid(bCoarse,2*m_dx,m_level-1); } m_lambda = m_dx*m_dx/(4*DIM); ...

10/24/2019 CS294-73 – Lecture 16

Multigrid Pseudocode

28

vcycle(φ, ρ) { φ := φ + λ(L(φ) − ρ) p times if (level > 0) { R = ρ − L(φ) Rc = A(R) δ : Bc → R , δ = 0 vcycle(δ, Rc) φ := φ + I(δ) φ := φ + λ ∗ (L(φ) − ρ) p times } else { φ := φ + λ ∗ (L(φ) − ρ) pB times { }

10/24/2019 CS294-73 – Lecture 16

Multigrid::vcycle

29

Multigrid::vCycle( BoxData<double >& a_phi, const BoxData<double >& a_rhs ){ pointRelax(a_phi,a_rhs,m_preRelax); if (m_level > 0) { residual(m_res,a_phi,a_rhs); avgDown(m_resc,m_res); m_delta.setVal(0.); m_coarsePtr->vCycle(m_delta,m_resc); fineInterp(a_phi,m_delta); pointRelax(a_phi,a_rhs,m_postRelax); } else{ pointRelax(a_phi,a_rhs,m_bottomRelax); }

10/24/2019 CS294-73 – Lecture 16

Multigrid::residual

30

Multigrid::residual( BoxData<double >& a_res, BoxData<double >& a_phi, const BoxData<double >& a_rhs ){ getGhost(a_phi); res.setVal(0.); double dxsqi = 1./(m_dx*m_dx); for (auto it = m_box.begin();it.done();++it){ Point pt = *it; a_res[pt] = -2*DIM*a_phi[pt]; for (int dir = 0; dir < DIM ; dir++) { Point edir = getUnitv(dir); a_res[pt] += (a_phi[pt + edir] + a_phi[pt-edir]); }; a_res[pt] = -a_res[pt]*dxsqi + a_rhs[pt]; } };

10/24/2019 CS294-73 – Lecture 16

Multigrid::pointRelax

31

Multigrid::pointRelax(BoxData<double >& a_phi, const BoxData<double >& a_rhs, int a_numIter){ for (int iter = 0; iter < a_numiter; iter++){ for (auto it = m_box.begin();it.done();++it){ Point = *it; res[pt] = .5*a_phi[pt]; for (int dir = 0; dir < DIM ; dir++){ Point edir = Basis(dir); res[pt] += (a_phi[pt + edir] + a_phi[pt-edir])/(4*DIM); } } for (auto it = m_box.begin();it.done();++it { Point pt =*it; a_phi[pt] = res[pt] - m_lambda*a_rhs[pt]; } } }

10/24/2019 CS294-73 – Lecture 16

Multigrid::getGhost

32

void Multigrid::getGhost(BoxData<double >& a_phi ){ ... for (int k = 0; k < 2*DIM; k++) { Box bx = m_bdry[k]; for (auto it = m_box.begin();it.done();++it){ Point pt = *it; int image[DIM]; for (int dir = 0; dir < DIM; dir++) { image[dir] = (pt[dir] + m_domainSize)%m_domainSize; } Point ptimage(image); a_phi[pt] = a_phi[ptimage]; } } ... };