Discretization for PDEs

1
Discretization for PDEs

• Chunfang Chen
• Danny Thorne
• Adam Zornes

2
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

3
Classification of PDEs (cont.)

• General form of second-order PDEs ( 2 variables)

4
PDE Model Problems

• Hyperbolic (Propagation)
• Advection equation (First-order linear)
• Wave equation (Second-order linear )

5
PDE Model Problems (cont.)

• Parabolic (Time- or space-marching)
• Burgers equation(Second-order nonlinear)
• Fourier equation (Second-order linear )

(Diffusion / dispersion)
6
PDE Model Problems (cont.)

• Elliptic (Diffusion, equilibrium problems)
• Laplace/Poisson (second-order linear)
• Helmholtz equation (second-order linear)

7
PDE Model Problems (cont.)

• System of Coupled PDEs
• Navier-Stokes Equations

8
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

9
Boundary and InitialConditions
?R

• Initial conditions starting point for
  propagation problems
• Boundary conditions specified on domain
  boundaries to provide the interior solution in
  computational domain

R
s
n
10
Numerical Methods

• Complex geometry
• Complex equations (nonlinear, coupled)
• Complex initial / boundary conditions
• No analytic solutions
• Numerical methods needed !!

11
Numerical Methods

• Objective Speed, Accuracy at minimum cost
• Numerical Accuracy (error analysis)
• Numerical Stability (stability analysis)
• Numerical Efficiency (minimize cost)
• Validation (model/prototype data, field data,
  analytic solution, theory, asymptotic solution)
• Reliability and Flexibility (reduce preparation
  and debugging time)
• Flow Visualization (graphics and animations)

12
computational solution procedures
Governing Equations ICS/BCS
System of Algebraic Equations
Equation (Matrix) Solver
ApproximateSolution
Discretization
Ui (x,y,z,t) p (x,y,z,t) T (x,y,z,t) or ?
(?,?,?,? )
Continuous Solutions
Finite-Difference Finite-Volume Finite-Element
Spectral Boundary Element Hybrid
Discrete Nodal Values
Tridiagonal ADI SOR Gauss-Seidel
Krylov Multigrid DAE
13
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

14
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

15
Finite Difference Scheme

• Forward difference
• Backward difference
• Central difference

16
Example Poisson Equation
(-1,1)
(1,1)
(-1,-1)
(1,-1)
17
Example (cont.)
18
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
• the only known grid
• points

19
What to solve?

• Discretization produces a linear system of
• equations.
• The A matrix is a
• pentadiagonal banded
• matrix of a standard
• form
• A solution method is to be performed for
• solving

20
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.

21
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.

22
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.

Weak Form
Matrix System
PDE
23
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).

24
Time Stepping Methods

• Standard methods are common
• Forward Euler (explicit)
• Backward Euler (implicit)
• Crank-Nicolson (implicit)

? 0, Fully-Explicit ? 1, Fully-Implicit ?
½, Crank-Nicolson
25
Time Stepping Methods (cont.)

• Variable length time stepping
• Most common in Method of Lines (MOL) codes or
  Differential Algebraic Equation (DAE) solvers

26
Solving the System

• The system may be solved using simple iterative
  methods - Jacobi, Gauss-Seidel, SOR, etc.
• Some advantages
• - No explicit storage of the matrix is required
• - The methods are fairly robust and reliable
• Some disadvantages
• - Really slow (Gauss-Seidel)
• - Really slow (Jacobi)

27
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

28
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

29
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

30
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

31
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

32
The Restriction Operator

• Defined as half-weighted restriction.
• Each new point in the courser grid, is dependent
  upon its neighboring points from the fine grid

33
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

34
The Prolongation Operator

• The grid change is exactly the opposite of
  restriction

35
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

36
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 e0 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
37
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/slides04-methods.
  ppt
• http//www.ccs.uky.edu/douglas/Classes/cs521-s01/
  index.html

38
Homework

• Introduce the following processing
• by read the book A Tutorial on Elliptic PDE
  Solvers and Their Parallelization

Weak Form
Discrete Matrix
PDE