CS 267 Applications of Parallel Computers Lecture 24: Solving Linear Systems arising from PDEs - IOutlinePoisson’s equation arises in many modelsRelation of Poisson’s equation to Gravity, ElectrostaticsPoisson’s equation in 1D2D Poisson’s equationAlgorithms for 2D Poisson Equation with N unknownsShort explanations of algorithms on previous slideComments on practical meshesComposite mesh from a mechanical structureConverting the mesh to a matrixIrregular mesh: NASA Airfoil in 2D (direct solution)Irregular mesh: Tapered Tube (multigrid)Adaptive Mesh Refinement (AMR)Jacobi’s MethodParallelizing Jacobi’s MethodSuccessive Overrelaxation (SOR)Gauss-SeidelSlide 19Conjugate Gradient (CG) for solving A*x = bSummary of Jacobi, SOR and CGSolving the Poisson equation with the FFTSerial FFTUsing the 1D FFT for filteringUsing the 2D FFT for image compressionRelated TransformsSerial Algorithm for the FFTDivide and Conquer FFTDivide-and-Conquer FFTAn Iterative AlgorithmParallel 1D FFTBlock Layout of 1D FFTCyclic Layout of 1D FFTParallel ComplexityFFT With “Transpose”Why is the Communication Step Called a Transpose?Complexity of the FFT with TransposeComment on the 1D Parallel FFTHigher Dimension FFTsCS267 L24 Solving PDEs.1Demmel Sp 1999CS 267 Applications of Parallel ComputersLecture 24: Solving Linear Systems arising from PDEs - IJames Demmelhttp://www.cs.berkeley.edu/~demmel/cs267_Spr99CS267 L24 Solving PDEs.2Demmel Sp 1999Outline°Review Poisson equation°Overview of Methods for Poisson Equation°Jacobi’s method°Red-Black SOR method°Conjugate Gradients°FFT°Multigrid (next lecture)Reduce to sparse-matrix-vector multiplyNeed them to understand MultigridCS267 L24 Solving PDEs.3Demmel Sp 1999Poisson’s equation arises in many models°Heat flow: Temperature(position, time)°Diffusion: Concentration(position, time)°Electrostatic or Gravitational Potential:Potential(position)°Fluid flow: Velocity,Pressure,Density(position,time)°Quantum mechanics: Wave-function(position,time)°Elasticity: Stress,Strain(position,time)CS267 L24 Solving PDEs.4Demmel Sp 1999Relation of Poisson’s equation to Gravity, Electrostatics°Force on particle at (x,y,z) due to particle at 0 is -(x,y,z)/r^3, where r = sqrt(x +y +z )°Force is also gradient of potential V = -1/r = -(d/dx V, d/dy V, d/dz V) = -grad V°V satisfies Poisson’s equation (try it!)2 2 2CS267 L24 Solving PDEs.5Demmel Sp 1999Poisson’s equation in 1D2 -1 -1 2 -1 -1 2 -1 -1 2 -1 -1 2T =2-1 -1Graph and “stencil”CS267 L24 Solving PDEs.6Demmel Sp 19992D Poisson’s equation°Similar to the 1D case, but the matrix T is now°3D is analogous4 -1 -1-1 4 -1 -1 -1 4 -1 -1 4 -1 -1 -1 -1 4 -1 -1 -1 -1 4 -1 -1 4 -1 -1 -1 4 -1 -1 -1 4T =4-1-1-1-1Graph and “stencil”CS267 L24 Solving PDEs.7Demmel Sp 1999Algorithms for 2D Poisson Equation with N unknownsAlgorithm Serial PRAM Memory #Procs°Dense LU N3N N2N2°Band LU N2N N3/2N°Jacobi N2 N N N°Explicit Inv. N log N N N°Conj.Grad. N 3/2N 1/2 *log N N N°RB SOR N 3/2N 1/2N N°Sparse LU N 3/2N 1/2N*log N N°FFT N*log N log N N N°Multigrid N log2 N N N°Lower bound N log N NPRAM is an idealized parallel model with zero cost communication2 22CS267 L24 Solving PDEs.8Demmel Sp 1999Short explanations of algorithms on previous slide°Sorted in two orders (roughly):•from slowest to fastest on sequential machines•from most general (works on any matrix) to most specialized (works on matrices “like” Poisson)°Dense LU: Gaussian elimination; works on any N-by-N matrix°Band LU: exploit fact that T is nonzero only on sqrt(N) diagonals nearest main diagonal, so faster°Jacobi: essentially does matrix-vector multiply by T in inner loop of iterative algorithm°Explicit Inverse: assume we want to solve many systems with T, so we can precompute and store inv(T) “for free”, and just multiply by it•It’s still expensive!°Conjugate Gradients: uses matrix-vector multiplication, like Jacobi, but exploits mathematical properies of T that Jacobi does not°Red-Black SOR (Successive Overrelaxation): Variation of Jacobi that exploits yet different mathematical properties of T•Used in Multigrid°Sparse LU: Gaussian elimination exploiting particular zero structure of T°FFT (Fast Fourier Transform): works only on matrices very like T°Multigrid: also works on matrices like T, that come from elliptic PDEs°Lower Bound: serial (time to print answer); parallel (time to combine N inputs)°Details in class notes and www.cs.berkeley.edu/~demmel/ma221CS267 L24 Solving PDEs.9Demmel Sp 1999Comments on practical meshes°Regular 1D, 2D, 3D meshes•Important as building blocks for more complicated meshes•We will discuss these first°Practical meshes are often irregular•Composite meshes, consisting of multiple “bent” regular meshes joined at edges•Unstructured meshes, with arbitrary mesh points and connectivities•Adaptive meshes, which change resolution during solution process to put computational effort where needed°In later lectures we will talk about some methods on unstructured meshes; lots of open problemsCS267 L24 Solving PDEs.10Demmel Sp 1999Composite mesh from a mechanical structureCS267 L24 Solving PDEs.11Demmel Sp 1999Converting the mesh to a matrixCS267 L24 Solving PDEs.12Demmel Sp 1999Irregular mesh: NASA Airfoil in 2D (direct solution)CS267 L24 Solving PDEs.13Demmel Sp 1999Irregular mesh: Tapered Tube (multigrid)CS267 L24 Solving PDEs.14Demmel Sp 1999Adaptive Mesh Refinement (AMR)°Adaptive mesh around an explosion°John Bell and Phil Colella at LBL (see class web page for URL)°Goal of Titanium is to make these algorithms easier to implementin parallelCS267 L24 Solving PDEs.15Demmel Sp 1999Jacobi’s Method°To derive Jacobi’s method, write Poisson as: u(i,j) = (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) + b(i,j))/4°Let u(i,j,m) be approximation for u(i,j) after m steps u(i,j,m+1) = (u(i-1,j,m) + u(i+1,j,m) + u(i,j-1,m) + u(i,j+1,m) + b(i,j)) / 4°I.e., u(i,j,m+1) is a weighted average of neighbors°Motivation: u(i,j,m+1) chosen to exactly satisfy equation at (i,j)°Convergence is proportional to
View Full Document
Unlocking...