Transcript Parallelization

Discretization for PDEs
Chunfang Chen
Danny Thorne
Adam Zornes
Classification of PDEs
Different mathematical and physical
behaviors:
 Elliptic Type
 Parabolic Type
 Hyperbolic Type
System of coupled equations for several
variables:
 Time : first-derivative (second-derivative for
wave equation)
 Space: first- and second-derivatives
Classification of PDEs (cont.)
General form of second-order PDEs ( 2
variables)
PDE Model Problems
 Hyperbolic (Propagation)
• Advection equation (First-order linear)
• Wave equation
(Second-order linear )
PDE Model Problems (cont.)
 Parabolic (Time- or space-marching)
• Burger’s equation(Second-order nonlinear)
(Diffusion / dispersion)
• Fourier equation (Second-order linear )
PDE Model Problems (cont.)
 Elliptic (Diffusion, equilibrium problems)
• Laplace/Poisson (second-order linear)
• Helmholtz equation (second-order linear)
PDE Model Problems (cont.)
 System of Coupled PDEs
• Navier-Stokes Equations
Well-Posed Problem
Numerically well-posed
 Discretization equations
 Auxiliary conditions (discretized approximated)
• the computational solution exists (existence)
• the computational solution is unique (uniqueness)
• the computational solution depends continuously on the
approximate auxiliary data
• the algorithm should be well-posed (stable) also
Boundary and Initial
Conditions
Initial conditions: starting
point for propagation problems
Boundary conditions: specified
on domain boundaries to provide
the interior solution in
computational domain
R
R
s
n
Numerical Methods
 Complex geometry
 Complex equations (nonlinear, coupled)
 Complex initial / boundary conditions
 No analytic solutions
 Numerical methods needed !!
Numerical Methods
Objective: Speed, Accuracy at minimum
cost
 Numerical
 Numerical
 Numerical
 Validation
Accuracy (error analysis)
Stability (stability analysis)
Efficiency (minimize cost)
(model/prototype data, field data,
analytic solution, theory, asymptotic solution)
 Reliability and Flexibility (reduce preparation
and debugging time)
 Flow Visualization (graphics and animations)
computational solution
procedures
Governing
Equations
ICS/BCS
Continuous
Solutions
System of
Algebraic
Equations
Equation
(Matrix)
Solver
Discrete
Nodal
Values
Tridiagonal
Ui (x,y,z,t)
ADI
p (x,y,z,t)
Finite-Element
SOR
T (x,y,z,t)
Spectral
Gauss-Seidel
Boundary Element
Krylov
Hybrid
Multigrid
Discretization
Finite-Difference
Finite-Volume
DAE
Approximate
Solution
or
 (,,, )
Discretization
 Time derivatives
almost exclusively by finite-difference methods
 Spatial derivatives
- Finite-difference: Taylor-series expansion
- Finite-element: low-order shape function and
interpolation function, continuous within each
element
- Finite-volume: integral form of PDE in each
control volume
- There are also other methods, e.g. collocation,
spectral method, spectral element, panel
method, boundary element method
Finite Difference
 Taylor series
 Truncation error
 How to reduce truncation errors?
• Reduce grid spacing, use smaller x = x-xo
• Increase order of accuracy, use larger n
Finite Difference Scheme
 Forward difference
 Backward difference
 Central difference
Example : Poisson Equation
(-1,1)
(1,1)
y
x
(-1,-1)
(1,-1)
Example (cont.)
Ui
 1, j
 2Ui , j  Ui
x
2
 1, j

Ui ,
j  1
 2Ui , j  Ui ,
y





  
 (i,j+1)y 
(i-1,j)x (i,j)(i+1,j)
 (i,j-1) 
  
2





j  1
1
Rectangular Grid
After we discretize the Poisson equation on a
rectangular domain, we are left with a finite
number of gird points. The boundary values
of the equation are
0,4 1,4 2,4 3,4 4,4
the only known grid
0,3 1,3 2,3 3,3 4,3
points

0,1
0,0
0,2

1,1
1,0
1,2

2,1
2,0
2,2

3,1
3,0
3,2

4,1
4,0
4,2
What to solve?
Discretization produces a linear system of
equations.
4 1
0
1
0
The A matrix is a
pentadiagonal banded
matrix of a standard
form:
1
4
1
0
1
0
1
4
1
0
1
0
0
1
1
0
A solution method is to be performed for
solving
4 1
1 4
Matrix Storage
We could try and take advantage of the
banded nature of the system, but a
more general solution is the adoption of
a sparse matrix storage strategy.
Limitations of Finite
Differences
 Unfortunately, it is not easy to use
finite differences in complex
geometries.
 While it is possible to formulate
curvilinear finite difference methods,
the resulting equations are usually
pretty nasty.
Finite Element Method
 The finite element method, while more
complicated than finite difference methods,
easily extends to complex geometries.
 A simple (and short) description of the finite
element method is not easy to give.
PDE
Weak
Matrix
Form
System
Finite Element Method
(Variational Formulations)
 Find u in test space H such that a(u,v) = f(v) for all v
in H, where a is a bilinear form and f is a linear
functional.




V(x,y) = j Vjj(x,y), j = 1,…,n
I(V) = .5 j j AijViVj - j biVi, i,j = 1,…,n
Aij =  a( j,  j), i = 1,…,n
Bi =  f j, i = 1,…,n
 The coefficients Vj are computed and the function
V(x,y) is evaluated anyplace that a value is needed.
 The basis functions should have local support (i.e.,
have a limited area where they are nonzero).
Time Stepping Methods
 Standard methods are common:
– Forward Euler (explicit)
– Backward Euler (implicit)
– Crank-Nicolson (implicit)
T jn 1  T jn
t
  LxxT
LxxT j 
n 1
j
 (1   ) LxxT
T j 1  2T j  T j 1
x 2
n
j
θ = 0, Fully-Explicit
θ = 1, Fully-Implicit
θ = ½, Crank-Nicolson
Time Stepping Methods (cont.)
 Variable length time stepping
– Most common in Method of Lines (MOL) codes or
Differential Algebraic Equation (DAE) solvers
Solving the System
 The system may be solved using simple
iterative methods - Jacobi, Gauss-Seidel,
SOR, etc.
•
Some advantages:
•
Some disadvantages
- No explicit storage of the matrix is required
- The methods are fairly robust and reliable
- Really slow (Gauss-Seidel)
- Really slow (Jacobi)
Solving the System
Advanced iterative methods (CG, GMRES)
CG is a much more powerful way to solve the
problem.
• Some advantages:
– Easy to program (compared to other advanced
methods)
– Fast (theoretical convergence in N steps for an N
by N system)
Some disadvantages:
– Explicit representation of the matrix is probably
necessary
– Applies only to SPD matrices
Multigrid Algorithm:
Components
 Residual
compute the error of the approximation
 Iterative method/Smoothing Operator
Gauss-Seidel iteration
 Restriction
obtain a ‘coarse grid’
 Prolongation
from the ‘coarse grid’ back to the
original grid
Residual Vector
 The equation we are to solve is defined as:
 So then the residual is defined to be:
 Where uq is a vector approximation for u
 As the u approximation becomes better,
the components of the residual vector(r) ,
move toward zero
Multigrid Algorithm:
Components
 Residual
compute the error of your approximation
 Iterative method/Smoothing Operator
Gauss-Seidel iteration
 Restriction
obtain a ‘coarse grid’
 Prolongation
from the ‘coarse grid’ back to the original
grid
Multigrid Algorithm:
Components
 Residual
compute the error of your approximation
 Iterative method/Smoothing Operator
Gauss-Seidel iteration
 Restriction
obtain a ‘coarse grid’
 Prolongation
from the ‘coarse grid’ back to the original
grid
The Restriction Operator
 Defined as ‘half-weighted’ restriction.
 Each new point in the courser grid, is
dependent upon it’s neighboring points from
the fine grid
Multigrid Algorithm:
Components
 Residual
compute the error of your approximation
 Iterative method/Smoothing Operator
Gauss-Seidel iteration
 Restriction
obtain a ‘coarse grid’
 Prolongation
from the ‘coarse grid’ back to the original
grid
The Prolongation Operator
 The grid change is exactly the opposite
of restriction
(0,2) 
(1,2) 
(2,2) 
(0,4) (1,4) (2,4) (3,4) (4,4)
(0,1) 
(1,1) 
(2,1) 
(0,2) (1,2) (2,2) (3,2) (4,2)
(0,0) 
(1,0) 
(2,0) 
(0,0) (1,0) (2,0) (3,0) (4,0)
(0,3) (1,3) (2,3) (3,3) (4,3)
(0,1) (1,1) (2,1) (3,1) (4,1)
Prolongation vs. Restriction
 The most efficient multigrid algorithms use
prolongation and restriction operators that
are directly related to each other. In the one
dimensional case, the relation between
prolongation and restriction is as follows:
Full Multigrid Algorithm
1.Smooth initial U vector to receive a new approximation Uq
2. Form residual vector: Rq =b -A Uq
3. Restrict Rq to the next courser grid  Rq-1
4. Smooth Ae= Rq-1 starting with e=0 to obtain eq-1
5.Form a new residual vector using: Rq-1= Rq-1 -A eq-1
6. Restrict R2 (5x? where ?5) down to R1(3x? where ?3)
7. Solve exactly for Ae= R1 to obtain e1
8. Prolongate e1e2 & add to e2 you got from restriction
9. Smooth Ae= R2 using e2 to obtain a new e2
10. Prolongate eq-1 to eq and add to Uq
Reference
www.mgnet.org/
http://csep1.phy.ornl.gov/CSEP/PDE/PDE.html
www.ceprofs.tamu.edu/hchen/
www.cs.cmu.edu/~ph/859B/www/notes/multigrid.pdf
www.cs.ucsd.edu/users/carter/260
www.cs.uh.edu/~chapman/teachpubs/slides04methods.ppt
http://www.ccs.uky.edu/~douglas/Classes/cs521s01/index.html
Homework
Introduce the following processing
by read the book “A Tutorial on Elliptic PDE
Solvers and Their Parallelization”
PDE
Weak
Discrete
Form
Matrix