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) 


Laplacian, Poisson’s equation, Heat Equation Laplacian: Poisson’s equation: Heat equation: Blue: Ω Black: ∂Ω 2 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 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). 5 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 Even though we don’t know the error, we can compute the residual to provide a measure for the error: e.g. convergence if . We can also see how the error behaves under point Jacobi 7 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 • The value at the new iteration is an average of the values at the old 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. 9 10/24/2019 CS294-73 – Lecture 16

Discrete Fourier Analysis of Point Jacobi ∆ h W ( k ) = 1 h 2 ( − 4 + 2 cos (2 π k x h ) + 2 cos (2 π k y h )) W ( k ) ≡ σ k W ( k ) − ( ∆ h ) − 1 W ( k ) = − 1 W ( k ) Let f = W (k) . Then . σ k Point Jacobi approximates this by φ l = (1 + λ ( σ k + 1)) φ l − 1 δ l = (1 + λσ k ) δ l − 1 = (1 + λσ k ) l δ 0 For the error to converge to zero for any Fourier mode, we require that 0 ≤ 1 + λσ k < 1 ⇣ N 2 , N : σ k = − 8 /h 2 , 1 + λσ k = 0 ⌘ k = 2 | k | h << 1 : λσ k = O (1) , 1 + λσ k = 1 + O ( h 2 ) Generally, high wave numbers in he error are damped efficiently, but not others. But this is relative to the grid spacing. 10 10/24/2019 CS294-73 – Lecture 16

Multigrid 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: ∆ h δ = R h , R h = ∆ h ˜ φ − f, φ = ˜ φ + δ ⇓ ∆ h φ = f captures the smooth remnants of the error. • Have to iterate, since this isn’t exact. 11 Deep mathematics: point Jacobi smooths <-> special property of Laplacian. 10/24/2019 CS294-73 – Lecture 16

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

Multigrid • If the number of grid points is 2 M +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) 13 10/24/2019 CS294-73 – Lecture 16

Averaging and Interpolation • 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. 14 10/24/2019 CS294-73 – Lecture 16

Results 15 10/24/2019 CS294-73 – Lecture 16

Finite-Volume Methods on Structured Grids We write the Laplacian as the divergence of fluxes. Using the Divergence Theorem, this leads to a natural discretization of the Laplacian. This discretization satisfies a discrete analogue of the divergence theorem. 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 Boundary conditions are expressed in terms of computing fluxes at boundary faces. Coarsening is done by aggregating volumes, rather than by sampling gridpoints. 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 vcycle ( φ , ρ ) { φ := φ + λ ( L ( φ ) − ρ ) p times if ( level > 0) { R = ρ − L ( φ ) R c = A ( R ) δ : B c → R , δ = 0 vcycle ( δ , R c ) φ := φ + I ( δ ) φ := φ + λ ∗ ( L ( φ ) − ρ ) p times } else { φ := φ + λ ∗ ( L ( φ ) − ρ ) p B times { } At the top level, iterate until residual is reduced by some factor. 23 10/24/2019 CS294-73 – Lecture 16

Implementing Multigrid The algorithm is naturally recursive – implement as nested set of calls to vcycle. 24 10/24/2019 CS294-73 – Lecture 16

main ... 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; } 25 10/24/2019 CS294-73 – Lecture 16

Multigrid.H ... 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. ... 26 10/24/2019 CS294-73 – Lecture 16

Multigrid::define ... 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); ... 27 10/24/2019 CS294-73 – Lecture 16

Multigrid Pseudocode vcycle ( φ , ρ ) { φ := φ + λ ( L ( φ ) − ρ ) p times if ( level > 0) { R = ρ − L ( φ ) R c = A ( R ) δ : B c → R , δ = 0 vcycle ( δ , R c ) φ := φ + I ( δ ) φ := φ + λ ∗ ( L ( φ ) − ρ ) p times } else { φ := φ + λ ∗ ( L ( φ ) − ρ ) p B times { } 28 10/24/2019 CS294-73 – Lecture 16