partial differential equations n.
Skip this Video
Loading SlideShow in 5 Seconds..
Partial Differential Equations PowerPoint Presentation
Download Presentation
Partial Differential Equations

Loading in 2 Seconds...

play fullscreen
1 / 78

Partial Differential Equations - PowerPoint PPT Presentation

  • Uploaded on

Partial Differential Equations. Introduction, Adam Zornes Discretizations and Iterative Solvers, Chenfang Chen Parallelization, Dr. Danny Thorne. What do You Stand For?. A PDE is a P artial D ifferential E quation This is an equation with derivatives of at least two variables in it.

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
Download Presentation

PowerPoint Slideshow about 'Partial Differential Equations' - mairwen

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
partial differential equations
Partial Differential Equations
  • Introduction, Adam Zornes
  • Discretizations and Iterative Solvers, Chenfang Chen
  • Parallelization, Dr. Danny Thorne
what do you stand for
What do You Stand For?
  • A PDE is a Partial Differential Equation
  • This is an equation with derivatives of at least two variables in it.
  • In general, partial differential equations are much more difficult to solve analytically than are ordinary differential equations
what does a pde look like
What Does a PDE Look Like
  • Let u be a function of x and y. There are several ways to write a PDE, e.g.,
    • ux + uy = 0
    • du/dx + du/dy = 0
the baskin robin s esq characterization of pde s
The Baskin Robin’s esq Characterization of PDE’s
  • The order is determined by the maximum number of derivatives of any term.
  • Linear/Nonlinear
    • A nonlinear PDE has the solution times a partial derivative or a partial derivative raised to some power in it
  • Elliptic/Parabolic/Hyperbolic
six one way
Six One Way
  • Say we have the following: Auxx + 2Buxy + Cuyy + Dux + Euy + F = 0.
  • Look at B2 - AC
    • < 0 elliptic
    • = 0 parabolic
    • > 0 hyperbolic
or half a dozen another
Or Half a Dozen Another
  • A general linear PDE of order 2:
  • Assume symmetry in coefficients so that A = [aij] is symmetric. Eig(A) are real. Let P and Z denote the number of positive and zero eigenvalues of A.
    • Elliptic: Z = 0 and P = n or Z = 0 and P = 0..
    • Parabolic: Z > 0 (det(A) = 0).
    • Hyperbolic: Z=0 and P = 1 or Z = 0 and P = n-1.
    • Ultra hyperbolic: Z = 0 and 1 < P < n-1.
elliptic not just for exercise anymore
Elliptic, Not Just For Exercise Anymore
  • Elliptic partial differential equations have applications in almost all areas of mathematics, from harmonic analysis to geometry to Lie theory, as well as numerous applications in physics.
  • The basic example of an elliptic partial differential equation is Laplace’s Equation
    • uxx - uyy = 0
the others
The Others
  • The heat equation is the basic Hyperbolic
    • ut - uxx - uyy = 0
  • The wave equations are the basic Parabolic
    • ut - ux - uy = 0
    • utt - uxx - uyy = 0
  • Theoretically, all problems can be mapped to one of these
what happens where you can t tell what will happen
What Happens Where You Can’t Tell What Will Happen
  • Types of boundary conditions
    • Dirichlet: specify the value of the function on a surface
    • Neumann: specify the normal derivative of the function on a surface
    • Robin: a linear combination of both
  • Initial Conditions
is it worth the effort
Is It Worth the Effort?
  • Basically, is it well-posed?
    • A solution to the problem exists.
    • The solution is unique.
    • The solution depends continuously on the problem data.
  • In practice, this usually involves correctly specifying the boundary conditions
why should you stay awake for the remainder of the talk
Why Should You Stay Awake for the Remainder of the Talk?
  • Enormous application to computational science, reaching into almost every nook and cranny of the field including, but not limited to: physics, chemistry, etc.
  • Laplace’s equation involves a steady state in systems of electric or magnetic fields in a vacuum or the steady flow of incompressible non-viscous fluids
  • Poisson’s equation is a variation of Laplace when an outside force is applied to the system
computational fluid dynamics
Computational Fluid Dynamics
  • CFD can be defined narrowly as confined to aerodynamic flow around vehicles but it can be generalized to include as well such areas as weather and climate simulation, flow of pollutants in the earth, and flow of liquids in oil fields (reservoir modelling).
  • Involve Huge PDE’s
  • Computational Science only Realistic Solution
discretization and iterative method for pdes

Discretization and Iterative Method for PDEs

Chunfang Chen

Danny Thorne

Adam Zornes

classification of pdes
Classification of PDEs

Different mathematical and physical


  • Elliptic Type
  • Parabolic Type
  • Hyperbolic Type

System of coupled equations for several


  • Time : first-derivative (second-derivative for wave equation)
  • Space: first- and second-derivatives
classification of pdes cont
Classification of PDEs (cont.)

General form of second-order PDEs ( 2 variables)

pde model problems
PDE Model Problems
  • Hyperbolic (Propagation)
  • Advection equation (First-order linear)
  • Wave equation (Second-order linear)
pde model problems cont
PDE Model Problems (cont.)
  • Parabolic (Time- or space-marching)
  • Burger’s equation(Second-order nonlinear)
  • Fourier equation (Second-order linear)

(Diffusion / dispersion)

pde model problems cont1
PDE Model Problems (cont.)
  • Elliptic (Diffusion, equilibrium problems)
  • Laplace/Poisson (second-order linear)
  • Helmholtz equation (second-order linear)
pde model problems cont2
PDE Model Problems (cont.)
  • System of Coupled PDEs
  • Navier-Stokes Equations
well posed problem
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
Boundary and InitialConditions


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




numerical methods
Numerical Methods
  • Complex geometry
  • Complex equations (nonlinear, coupled)
  • Complex initial / boundary conditions
  • No analytic solutions
  • Numerical methods needed !!
numerical methods1
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
computational solution procedures

Governing Equations ICS/BCS

System of Algebraic Equations

Equation (Matrix) Solver



Ui (x,y,z,t)

p (x,y,z,t)

T (x,y,z,t)


 (,,, )

Continuous Solutions





Boundary Element


Discrete Nodal Values








  • 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


- 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
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
Finite Difference Scheme
  • Forward difference
  • Backward difference
  • Central difference
example poisson equation



Example : Poisson Equation





rectangular grid
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


what to solve
What to solve?

Discretization produces a linear system of


The A matrix is a

pentadiagonal banded

matrix of a standard


A solution method is to be performed for


matrix storage
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
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
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.






finite element method variational formulations
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
Time Stepping Methods
  • Standard methods are common:
    • Forward Euler (explicit)
    • Backward Euler (implicit)
    • Crank-Nicolson (implicit)

θ = 0, Fully-Explicit

θ = 1, Fully-Implicit

θ = ½, Crank-Nicolson

time stepping methods cont
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
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 system1
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
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
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 components1
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 components2
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
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 components3
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 Prolongation Operator
  • The grid change is exactly the opposite of restriction
prolongation vs restriction
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
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



HPC Issues in PDEs, Part 3

Danny Thorne, CS521, University of Kentucky, April 9, 2002

parallel computation
Parallel Computation
  • Serious calculations today are mostly done on a parallel computer.
  • The domain is partitioned into subdomains that may or may not overlap slightly.
  • Goal is to calculate as many things in parallel as possible even if some things have to be calculated on several processors in order to avoid communication.
  • ”Communication is the Darth Vader of parallel computing.”
example original mesh
Example: Original Mesh
  • Consider solving a PDE on this grid.
  • Suppose that only half of the nodes fit on a processor.
example mesh on two processors
Example: Mesh on Two Processors
  • Divide into two connected subsets and renumber within the subdomains.
  • Communication occurs between neighbors that cross the processor boundary.
  • Ghost points (or overlap) can reduce communication sometimes at the expense of extra computation.
  • Computation is communication per word.
mesh decomposition
Mesh Decomposition
  • Goals are to maximize interior while minimizing connections between subdomains. That is, minimize communication.
  • Such decomposition problems have been studied in load balancing for parallel computation.
  • Lots of choices:
    • METIS, ParMETIS -- University of Minnesota.
    • PARTI -- University of Maryland,
    • CHACO -- Sandia National Laboratories,
    • JOSTLE -- University of Greenwich,
    • PARTY -- University of Paderborn,
    • SCOTCH -- Université Bordeaux,
    • TOP/DOMDEC -- NAS at NASA Ames Research Center.
mesh decomposition1
Mesh Decomposition
  • Quality of mesh decomposition has a highly significant effect on performance.
    • Arriving at a “good” decomposition is a complex task.
    • “Good” may be problem/architecture dependent.
  • A wide variety of well-established methods exist.
    • Several packages/libraries provide implementations of many of these methods.
    • Major practical difficulty is differences in file formats.

decomposition goals
Decomposition Goals
  • Load balance.
    • Distribute elements evenly across processors.
    • Each processor should have equal share of work.
  • Communication costs should be minimized.
    • Minimize sub-domain boundary elements.
    • Minimize number of neighboring domains.
  • Distribution should reflect machine architecture.
    • Communication versus calculation.
    • Bandwidth versus latency.
  • Note that optimizing load balance and communication cost simultaneously is an NP-hard problem.

static mesh versus dynamic mesh
Static Mesh versus Dynamic Mesh
  • Static mesh
    • Decomposition need only be carried out once.
    • Static decomposition may therefore be carried out as a preprocessing step, often done in serial.
  • Dynamic mesh
    • Decomposition must be adapted as underlying mesh changes.
    • Dynamic decomposition therefore becomes part of the calculation itself and cannot be carried out solely as a pre-processing step.

rod impact
Rod Impact
  • Notice that the mesh changes dynamically.

  • Need a dynamic mesh for this?

dynamic decomposition
Dynamic Decomposition
  • To redistribute a mesh
    • A new decomposition must be found and required changes made (shuffling of data between processors).
    • The time taken may outweigh the benefit gained so an assessment of whether it is worthwhile is needed.
  • The decomposition must run in parallel and be fast.
    • Take into account the old decomposition.
    • Make minimal changes.

dual graphs
Dual Graphs
  • To make a decomposition we need
    • A representation of the basic entities being distributed (e.g. the elements, in the case of a FE mesh),
    • An idea of how communication takes place between them.
  • A dual graph, based on the mesh, fills this role.
    • Vertices in the graph represent the entities (elements).
    • Edges in the graph represent communication.
  • Graph depends on how data is transferred.
    • For finite elements it could be via nodes, edges or faces.
    • So a single mesh can have more than one dual graph.
  • Edges (communications) or vertices (calculations) may be weighted.

2d example
2D Example

graph partitioning
Graph Partitioning
  • This problem occurs in several applications.
    • VLSI circuit layout design,
    • Efficient orderings for parallel sparse matrix factorization,
    • Mesh decomposition.
  • This problem has been extensively studied.
    • A large body of literature already exists.
    • Much work is still being produced in the field.

graph partitioning1
Graph Partitioning
  • Given a graph we want to partition it into parts such that:
    • Each part has roughly the same number of vertices,
    • Number of edges that cross partitions is minimized.

mesh decomposition2
Mesh Decomposition
  • View mesh decomposition as having two aspects:
    • Partitioning a graph,
    • Mapping the sub-domains to processors.
  • Algorithms may
    • Address both issues,
    • Or simplify the problem by ignoring the latter.
  • Whether this is an issue depends on target machine.
    • With very fast communications, e.g. the T3D, it's probably less important.
    • On a cluster of workstations it may be absolutely critical.

decomposition algorithms
Decomposition Algorithms
  • Global methods
    • Direct k-way partitioning.
    • Recursive application of some -way technique, e.g. bisection.
  • Local refinement techniques
    • Incrementally improve quality of an existing decomposition.
  • Multi-level variants
    • Work with coarse approximation of graph to increase performance.
  • Hybrid techniques
    • Using various combinations of above.

  • Simple, direct k-way algorithms
    • Random, Scattered & Linear Bandwidth Reduction
    • The Greedy Algorithm
  • Optimization algorithms
    • Simulated Annealing
    • Chained Local Optimization
    • Genetic Algorithms
  • Geometry based recursive algorithms
    • Coordinate Partitioning
    • Inertial Partitioning

  • Graph based recursive algorithms
    • Layered Partitioning
    • Spectral Partitioning
  • Local refinement algorithms
    • Kernighan and Lin
    • Jostle
  • Multi-level and hybrid variants
    • Multi-level Kernighan and Lin Partitioning
    • Multi-level Spectral RQL / Symmlq Partitioning

space filling curves
Space Filling Curves
  • Locality preserving.
  • Each point lies a unique distance along the curve.
space filling curves1
Space Filling Curves
  • Optimal load balance.
  • Subdomain boundaries are sub-optimal.
  • Recall: Optimizing load and comm is NP-hard.
space filling curves2
Space Filling Curves
  • The main advantages of this partition method are:
    • It is fast compared to graph partitioning heuristics,
    • It runs in parallel,
    • It requires no administration and no storage of processor neighborhoods.
    • The knowledge of the separators is enough to compute where to find a node and which processor to ask for it.
parallel pde tools
Parallel PDE Tools
  • Parallel ELLPACK -- Purdue University
  • POOMA -- Los Alamos National Laboratory
  • Overture -- Lawrence Livermore National Laboratory
  • PETSc -- Argonne National Laboratory
  • Unstructured Grids -- University of Heidelberg
  • A Software Framework for Easy Parallelization of PDE Solvers --
  • Numerical Objects Online --
  • Xing Cai's home page at University of Oslo --
  • Related Grid Decomposition Sites --
  • Domain Decomposition Methods --
  • People working on Domain Decomposition --
  • Unstructured Mesh Decomposition --
  • HPCI Seminar: Parallel Sparse Matrix Solvers --
  • METIS: Family of Multilevel Partitioning Algorithms --
  • Multilevel Algorithms for Multi-Constraint Graph Partitioning --
  • Structured Versus Unstructured Meshes --
  • Domain Decomposition Methods for elliptic PDEs --
  • Multiblock Parti library --
  • Chaco: Software for Partitioning Graphs --
  • A Portable and Efficient Parallel Code for Astrophysical Fluid Dynamics --
  • Links for Graph Partitioning --
  • A Multilevel Algorithm for Partitioning Graphs --
  • Multigrid Workbench: The Algorithm --
  • Domain DecompositionParallelization of Mesh Based Applications --
  • Parallel ELLPACK Elliptic PDE Solvers --
  • POOMA --
  • Overture --
  • PADRE: A Parallel Asynchronous Data Routing Environment --
  • PETSc: The Portable, Extensible Toolkit for Scientific Computation --
  • UG Home Page --
  • Space-filling Curves --
  • Separate IMAGE for Basic foil 30 Two Space Filling Curves --
  • Using Space-filling Curves for Multi-dimensional Indexing --
  • The origins of fractals --
  • Adaptive Parallel Multigrid Methods --