Parallel Numerical Algorithms Chapter 7 Differential Equations - PowerPoint PPT Presentation

Domain Decomposition Computation with Grids Scalability Parallel Numerical Algorithms Chapter 7 – Differential Equations Section 7.2 – Partial Differential Equations Michael T. Heath Department of Computer Science University of Illinois at Urbana-Champaign CS 554 / CSE 512 Michael T. Heath Parallel Numerical Algorithms 1 / 49

Domain Decomposition Computation with Grids Scalability Outline Domain Decomposition 1 Overlapping Subdomains Non-Overlapping Subdomains Computation with Grids 2 Parallel Computation with Grids Multigrid Scalability 3 Scalability Michael T. Heath Parallel Numerical Algorithms 2 / 49

Domain Decomposition Computation with Grids Scalability Outline Domain Decomposition 1 Overlapping Subdomains Non-Overlapping Subdomains Computation with Grids 2 Parallel Computation with Grids Multigrid Scalability 3 Scalability Michael T. Heath Parallel Numerical Algorithms 3 / 49

Domain Decomposition Computation with Grids Scalability Numerical Methods for Partial Differential Equations Partial differential equations are typically solved numerically by finite difference, finite element, finite volume, or spectral discretization Such discretization yields system of linear or nonlinear algebraic equations whose solution gives approximate solution to PDE Solving linear or nonlinear algebraic system is one major source of parallelism in solving PDEs numerically We will consider domain decomposition methods that exploit natural parallelism in PDE and its discretization Michael T. Heath Parallel Numerical Algorithms 4 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Alternating Schwarz Method Consider elliptic PDE L u = f on domain Ω = Ω 1 ∪ Ω 2 , with boundary condition u = g on ∂ Ω Γ 1 Ω 1 Ω 2 Γ 2 Michael T. Heath Parallel Numerical Algorithms 5 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Alternating Schwarz Method Given initial guess u (0) on Ω 2 , for k = 0 , 1 , . . . 2 On Ω 1 , solve L u ( k +1) = f 1 with boundary conditions u ( k +1) = g on ∂ Ω 1 \ Γ 1 1 u ( k +1) u ( k ) = on Γ 1 1 2 On Ω 2 , solve L u ( k +1) = f 2 with boundary conditions u ( k +1) = g on ∂ Ω 2 \ Γ 2 2 u ( k +1) u ( k +1) = on Γ 2 2 1 Michael T. Heath Parallel Numerical Algorithms 6 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Alternating Schwarz Method Alternating iterations continue until convergence to solution u on entire domain Ω Schwarz proposed this method in 1870 to deal with regions for which analytical solutions are not known Today it is of interest, in discretized form, for suggesting one of two major paradigms for solving PDEs numerically by domain decomposition Overlapping subdomains (Schwarz) Non-overlapping subdomains (Schur) Michael T. Heath Parallel Numerical Algorithms 7 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Discretized Schwarz Method Discretization yields n × n symmetric positive definite linear algebraic system Ax = b For i = 1 , 2 , let S i be set of n i indices of grid points in interior of Ω i , where n i = | S i | Because subdomains overlap, S 1 ∩ S 2 � = ∅ and n 1 + n 2 > n For i = 1 , 2 , let R i be n i × n Boolean restriction matrix such that for any n -vector v , v i = R i v contains precisely those components of v corresponding to indices in S i (i.e., nodes in Ω i ) Michael T. Heath Parallel Numerical Algorithms 8 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Discretized Schwarz Method Conversely, n × n i extension matrix R T i expands n i -vector v i to n -vector v whose components corresponding to indices in S i are same as those of v i , and whose remaining components are all zero Principal submatrices of A of order n 1 and n 2 corresponding to two subdomains are given by A 1 O R 1 AR T = A 1 1 R 2 AR T A 2 = 2 O A 2 Michael T. Heath Parallel Numerical Algorithms 9 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Discretized Schwarz Method For discretized problem, alternating Schwarz iteration takes form x k + R T 1 A − 1 x k + 1 = 1 R 1 ( b − Ax k ) 2 2 + R T 2 A − 1 x k +1 = x k + 1 2 R 2 ( b − Ax k + 1 2 ) This method is analogous to block Gauss-Seidel, but with overlapping blocks Overall error e k = x − x k updated as e k +1 = Ge k where G = ( I − R T 2 A − 1 2 R 2 A )( I − R T 1 A − 1 1 R 1 A ) so this is known as multiplicative Schwarz method Michael T. Heath Parallel Numerical Algorithms 10 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Discretized Schwarz Method We have as yet achieved no parallelism, since two subproblems must be solved sequentially for each iteration, but instead of Gauss-Seidel, we can use the block Jacobi approach 1 A − 1 x k + R T x k + 1 = 1 R 1 ( b − Ax k ) 2 2 A − 1 2 + R T x k +1 = x k + 1 2 R 2 ( b − Ax k ) whose subproblems can be solved simultaneously With either Gauss-Seidel or Jacobi version, it can be shown that iteration converges at rate independent of mesh size, provided overlap area between subdomains is sufficiently large (and mesh is refined uniformly) Michael T. Heath Parallel Numerical Algorithms 11 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Discretized Schwarz Method Eliminating x k + 1 2 in Jacobi version, we obtain x k +1 = x k + ( R T 1 A − 1 1 R 1 + R T 2 A − 1 2 R 2 )( b − Ax k ) which is just Richardson iteration with additive Schwarz preconditioner M − 1 = R T 1 A − 1 2 A − 1 1 R 1 + R T 2 R 2 Symmetry of preconditioned system means it can be used in conjunction with conjugate gradient method to accelerate convergence Michael T. Heath Parallel Numerical Algorithms 12 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Discretized Schwarz Method Multiplicative Schwarz iteration matrix is not symmetric, but can be made symmetric by additional step with A − 1 each 1 iteration x k + R T 1 A − 1 = 1 R 1 ( b − Ax k ) x k + 1 3 3 + R T 2 A − 1 x k + 2 = x k + 1 2 R 2 ( b − Ax k + 1 3 ) 3 3 + R T 1 A − 1 = 1 R 1 ( b − Ax k + 2 3 ) x k +1 x k + 2 which yields symmetric preconditioner that can be used in conjunction with conjugate gradient method to accelerate convergence Michael T. Heath Parallel Numerical Algorithms 13 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Many Overlapping Subdomains To achieve higher degree of parallelism with Schwarz method, we can apply two-domain algorithm recursively or use many subdomains If there are p overlapping subdomains, then define matrices R i and A i as before, i = 1 , . . . , p Additive Schwarz preconditioner then takes form p � R T i A − 1 i R i i =1 Michael T. Heath Parallel Numerical Algorithms 14 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Many Overlapping Subdomains Resulting generalization of block-Jacobi iteration is highly parallel, but not algorithmically scalable because convergence rate degrades as p grows Convergence rate can be restored by using coarse grid correction to provide global coupling If R 0 and R T 0 are restriction and interpolation matrices between coarse and fine grids, and A 0 = R 0 AR T 0 , then additive Schwarz preconditioner becomes p � R T i A − 1 i R i i =0 Michael T. Heath Parallel Numerical Algorithms 15 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Many Overlapping Subdomains Multiplicative Schwarz iteration for p domains is defined analogously As with classical Gauss-Seidel vs. Jacobi, multiplicative Schwarz has faster convergence rate than corresponding additive Schwarz (though it still requires coarse grid correction to remain scalable) But unfortunately, multiplicative Schwarz appears to provide no parallelism, as p subproblems per iteration must be solved sequentially As with classical Gauss-Seidel, parallelism can be introduced by coloring subdomains to identify independent subproblems that can be solved simultaneously Michael T. Heath Parallel Numerical Algorithms 16 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Many Overlapping Subdomains Ω 1 Ω 3 Ω 2 Ω 4 Ω 9 Ω 11 Ω 10 Ω 12 Ω 5 Ω Ω 6 Ω 8 7 Ω Ω Ω Ω 13 15 14 16 Michael T. Heath Parallel Numerical Algorithms 17 / 49

Domain Decomposition Overlapping Subdomains Computation with Grids Non-Overlapping Subdomains Scalability Non-Overlapping Subdomains We now consider adjacent subdomains whose only points in common are along their mutual boundary Ω 1 Ω 2 Γ We partition indices of unknowns in corresponding discrete linear system into three sets, S 1 and S 2 corresponding to interior nodes in Ω 1 and Ω 2 , respectively, and S 3 corresponding to interface nodes in Γ Michael T. Heath Parallel Numerical Algorithms 18 / 49