COSC 6374 Parallel Computation Partial Differential Equations Edgar Gabriel Fall 2011 Edgar Gabriel Numerical differentiation forward difference formula y From the definition of derivatives f x f x h f x f x lim h 0 h f xn one can derive an approximation for the 1st derivative f x f x h f x h xn h xn h x The same formula can be obtained from the Taylor 2 series e g f x h f x hf x h f Parallel Computation Edgar Gabriel 2 f x h hf x h f x f h 2 1 Center Difference Formula A better formula is derived if looking at the following two terms y 2 3 2 1 2 2 f x h f x hf x f x h f x hf x h h f x f 1 2 3 2 f x h h f x f 2 2 3 Subtracting equation 2 2 from 2 1 leads to f x f xn 3 hh 2 1 h f x h f x h 2h 12 xn h xn h x h is quadratic in the error term Parallel Computation Edgar Gabriel Center Difference Formula for 2nd Derivatives Extend 2 1 and 2 2 by an additional term f x h f x hf x h2 h3 h 4 4 f x f x f 1 2 3 4 f x h f x hf x h2 h3 h 4 4 f x f x f 2 2 3 4 Adding both equations leads to f x 1 h4 f x h 2 f x f x h 2 h 12 Parallel Computation Edgar Gabriel 2 Numerical differentiation summary Forward difference formula f x f x h f x h Center difference formula for the 1st derivative f x 1 f x h f x h 2h Center difference formula for the 2nd derivative f x 1 f x h 2 f x f x h h2 Parallel Computation Edgar Gabriel Differential equations terminology Differential equations equations containing the derivative of a function as a variable An ordinary differential equation ODE only contains functions of one independent variable A partial differential equation PDE contains functions of multiple independent variables and their partial derivatives The order of a differential equation is that of the highest derivative that it contains The goal is to find a function y t whose derivatives fulfill the given differential equations e g y n t f t y y y y n 1 Parallel Computation Edgar Gabriel 3 Finite Differences Approach for Solving Differential Equations If the analytic solution of the DE can not be determined calculate an approximate solution in discrete locations Replace the derivatives in the DE by an according approximation formula y t 1 y t h y t h 2h y t 1 y t h 2 y t y t h h2 Parallel Computation Edgar Gabriel Example I Solve the following two point boundary value problem using the finite difference method d2y dy 2 10 x 0 dx 2 dx y 0 1 y 1 2 0 x 1 For simplicity lets assume the points of interest are equally spaced b a xi a ih 0 i n 1 h n 1 e g for h 0 2 the mesh points are x0 0 x1 0 2 x2 0 4 x3 0 6 x4 0 8 x5 1 0 Due to the boundary values y0 y x0 1 y5 y x5 2 y1 y4 are unknown Parallel Computation Edgar Gabriel 4 Example II Discrete version of the ODE using central differences 1 1 yi 1 2 yi yi 1 2 yi 1 yi 1 10 xi 0 2 h 2h 1 1 yi 1 yi 1 10 xi 0 yi 1 2 yi yi 1 2 2 0 2 0 4 25 yi 1 2 yi yi 1 5 yi 1 yi 1 10 xi 0 20 yi 1 50 yi 30 yi 1 10 xi Parallel Computation Edgar Gabriel Example III i 1 20 y0 50 y1 30 y2 10 x1 20 50 y1 30 y2 10 0 2 i 1 i 2 i 3 i 4 or 22 50 y1 30 y2 20 y1 50 y2 30 y3 4 20 y 2 50 y3 30 y4 6 20 y3 50 y 4 68 50 30 y1 22 20 50 30 y 4 2 20 50 30 y3 6 20 50 y4 68 A y b Parallel Computation Edgar Gabriel 5 Solving Ay b using an iterative solver Given A b and an initial guess y0 r0 b Ay0 r such that r 0T r0 0 0 0 1 v0 p0 0 e g the Bi CGSTAB algorithm Given for i 1 2 i r 0T ri 1 Matrix vector multiplication i i 1 i 1 pi ri 1 pi 1 i 1vi 1 vi Api Scalar product i T r 0 v i s ri 1 vi t As i tT s tTt yi yi 1 pi i s ri s i t Parallel Computation Edgar Gabriel Scalar product in parallel Process with rank 0 Scalar product a 0 N N 1 s a i b i 2 1 b 0 N 2 Process with rank 1 1 a N N 2 b N N 2 i 0 Parallel algorithm s N 2 1 a i b i i 0 N 2 1 N 1 a i b i i N 2 N 2 1 a i blocal i alocal i blocal i i 0 i 0 1 44424443 1 44424443 local rank 0 rank 1 requires communication between the processes Parallel Computation Edgar Gabriel 6 Matrix vector product in parallel 50 30 x1 rhs1 20 50 30 x rhs 2 2 20 50 30 x3 rhs3 20 50 x4 rhs 4 50 x1 30 x 2 Process 0 Process 1 rhs1 rhs 2 20 x1 50 x2 30 x3 20 x2 50 x3 30 x4 rhs3 Process 0 needs x3 Process 1 needs x2 20 x3 50 x4 rhs4 Parallel Computation Edgar Gabriel Matrix vector product in parallel II Introduction of ghost cells Process zero x1 x2 x3 Process one x2 x3 x4 Looking at the source code e g pi ri 1 pi 1 i 1vi 1 vi Api since the vector used in the matrix vector multiplication changes every iteration you always have to update the ghost cells before doing the calculation Parallel Computation Edgar Gabriel 7 Matrix vector product in parallel III so the parallel algorithm for the same area is pi ri 1 pi 1 i 1vi 1 Update the ghost cells of p e g Process 0 sends p 2 to Process 1 Process 1 sends p 3 to Process 0 vi Api Parallel Computation Edgar Gabriel 2D Example Laplace equation I 2 D Laplace equation 2 2 u x y 2 u x y 0 x 2 y Central discretization leads to ui 1 j 2ui j ui 1 j h2 ui j 1 2ui j ui j 1 h2 0 i j 1 i 1 j i j i 1 j i j 1 Parallel Computation …