Discretization for PDEs by Sn69ct4Y

VIEWS: 17 PAGES: 38

									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         R
point for propagation problems
                                  R
Boundary conditions: specified
on domain boundaries to provide
the interior solution in          s
computational domain
                                      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   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)
             computational solution
                 procedures
Governing                        System of      Equation    Approximate
Equations      Discretization    Algebraic      (Matrix)     Solution
ICS/BCS                          Equations       Solver

Continuous   Finite-Difference   Discrete    Tridiagonal    Ui (x,y,z,t)
 Solutions                        Nodal
               Finite-Volume                    ADI         p (x,y,z,t)
                                 Values
              Finite-Element                    SOR
                                                            T (x,y,z,t)
                 Spectral                    Gauss-Seidel
                                                                or
             Boundary Element                  Krylov
                  Hybrid                      Multigrid      (,,, )
                                                DAE
             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     1, j       Ui ,   j  1    2Ui , j  Ui ,   j  1
                                                                                      1
                     2                                             2
                x                                            y
                                                            
                                      (i,j+1)y             
                              (i-1,j)x (i,j)(i+1,j)        
                                      (i,j-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,2   1,2   2,2   3,2   4,2
                     0,1   1,1   2,1   3,1   4,1

                     0,0   1,0   2,0   3,0   4,0
             What to solve?
Discretization produces a linear system of
equations.             4 1       0    1   0
The A matrix is a          1    4    1    0    1
pentadiagonal banded       0    1    4    1    0
matrix of a standard
                           1     0    1    4   1
form:
                           0    1     0    1    4
A solution method is to be performed for
solving
        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.


                     Weak                Matrix
 PDE
                     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                                             θ = 0, Fully-Explicit
                   LxxT    n 1
                                     (1   ) LxxT    n

     t
                            j                         j
                                                           θ = 1, Fully-Implicit
                    T j 1  2T j  T j 1
        LxxT j                                            θ = ½, Crank-Nicolson
                             x 2
  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:
     - No explicit storage of the matrix is required
     - The methods are fairly robust and reliable


•   Some disadvantages
     - 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,3) (1,3) (2,3) (3,3) (4,3)

 (0,1)    (1,1)    (2,1)    (0,2) (1,2) (2,2) (3,2) (4,2)
                               (0,1) (1,1) (2,1) (3,1) (4,1)
 (0,0)    (1,0)    (2,0)    (0,0) (1,0) (2,0) (3,0) (4,0)
  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/slides04-
methods.ppt
http://www.ccs.uky.edu/~douglas/Classes/cs521-
s01/index.html
              Homework
   Introduce the following processing
by read the book “A Tutorial on Elliptic PDE
     Solvers and Their Parallelization”


                   Weak             Discrete
 PDE
                   Form              Matrix

								
To top