Partial Differential Equations

1 / 78

# Partial Differential Equations - PowerPoint PPT Presentation

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.

## PowerPoint Slideshow about 'Partial Differential Equations' - mairwen

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
• Discretizations and Iterative Solvers, Chenfang Chen
• Parallelization, Dr. Danny Thorne
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
• 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 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
• 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
• 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 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 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
• 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?
• 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?
• Enormous application to computational science, reaching into almost every nook and cranny of the field including, but not limited to: physics, chemistry, etc.
Example
• 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
• 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
• http://www.npac.syr.edu/users/gcf/cps713overI94/
• http://www.cse.uiuc.edu/~rjhartma/pdesong.html
• http://www.maths.soton.ac.uk/teaching/units/ma274/node2.html
• http://www.npac.syr.edu/projects/csep/pde/pde.html
• http://mathworld.wolfram.com/PartialDifferentialEquation.html

### Discretization and Iterative Method for PDEs

Chunfang Chen

Danny Thorne

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)
• Wave equation (Second-order linear)
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 (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 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

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

SOR

Gauss-Seidel

Krylov

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

y

x

Example : Poisson Equation

(-1,1)

(1,1)

(-1,-1)

(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

the only known grid

points

What to solve?

Discretization produces a linear system of

equations.

The A matrix is a

matrix of a standard

form:

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

Form

Matrix

System

PDE

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)

θ = 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.

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

• Easy to program (compared to other advanced methods)
• Fast (theoretical convergence in N steps for an N by N system)
• 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
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
• http://csep1.phy.ornl.gov/CSEP/PDE/PDE.html
• www.mgnet.org/
• 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

### Parallelization

HPC Issues in PDEs, Part 3

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

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
• Consider solving a PDE on this grid.
• Suppose that only half of the nodes fit on a processor.
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
• 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 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.

http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/MeshDecomp-11.html

Decomposition Goals
• 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.

http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/MeshDecomp-13.html

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.

http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/MeshDecomp-14.html

Rod Impact
• Notice that the mesh changes dynamically.

http://www.amtec.com/plotgallery/animation/movies/t-rod.gif

Frisbee
• Need a dynamic mesh for this?

http://www.llnl.gov/CASC/Overture/henshaw/figures/pib.mpg

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.

http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/MeshDecomp-15.html

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.

http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/MeshDecomp-16.html

2D Example

http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/MeshDecomp-17.html

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.

http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/MeshDecomp-18.html

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.

http://www-users.cs.umn.edu/~karypis/publications/Talks/multi-constraint/sld002.htm

Mesh Decomposition
• View mesh decomposition as having two aspects:
• Partitioning a graph,
• Mapping the sub-domains to processors.
• Algorithms may
• 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.

http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/MeshDecomp-22.html

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.

http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/MeshDecomp-24.html

Examples
• 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

http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/MeshDecomp-25.html

Examples
• 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

http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/MeshDecomp-26.html

Space Filling Curves
• Locality preserving.
• Each point lies a unique distance along the curve.
Space Filling Curves
• Subdomain boundaries are sub-optimal.
• Recall: Optimizing load and comm is NP-hard.
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 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 -- http://www.ifi.uio.no/~xingca/PCFD2000/Slides/pcfd2000/
• Numerical Objects Online -- http://www.nobjects.com/Diffpack/
• Related Grid Decomposition Sites -- http://www.erc.msstate.edu/~bilderba/research/other_sites.html
• Domain Decomposition Methods -- http://www.ddm.org/
• People working on Domain Decomposition -- http://www.ddm.org/people.html
• Unstructured Mesh Decomposition -- http://www.epcc.ed.ac.uk/epcc-tec/documents/meshdecomp-slides/
• HPCI Seminar: Parallel Sparse Matrix Solvers -- http://www.epcc.ed.ac.uk/computing/research_activities/HPCI/sparse_seminar.html
• METIS: Family of Multilevel Partitioning Algorithms -- http://www-users.cs.umn.edu/~karypis/metis/index.html
• Multilevel Algorithms for Multi-Constraint Graph Partitioning -- http://www-users.cs.umn.edu/~karypis/publications/Talks/multi-constraint/sld001.htm
• Structured Versus Unstructured Meshes -- http://www.epcc.ed.ac.uk/epcc-tec/courses/images/Meshes.gif
• Domain Decomposition Methods for elliptic PDEs -- http://www.cs.purdue.edu/homes/mav/projects/dom_dec.html
• Multiblock Parti library -- http://www.cs.umd.edu/projects/hpsl/compilers/base_mblock.html
• Chaco: Software for Partitioning Graphs -- http://www.cs.sandia.gov/~bahendr/chaco.html
• A Portable and Efficient Parallel Code for Astrophysical Fluid Dynamics -- http://astro.uchicago.edu/Computing/On_Line/cfd95/camelse.html
• A Multilevel Algorithm for Partitioning Graphs -- http://www.chg.ru/SC95PROC/509_BHEN/SC95.HTM
• Multigrid Workbench: The Algorithm -- http://www.mgnet.org/mgnet/tutorials/xwb/mgframe.html
• Domain DecompositionParallelization of Mesh Based Applications -- http://www.hlrs.de/organization/par/par_prog_ws/online/DomainDecomposition/index.htm
• Parallel ELLPACK Elliptic PDE Solvers -- http://www.cs.purdue.edu/homes/markus/research/pubs/pde-solvers/paper.html
• POOMA -- http://www.acl.lanl.gov/pooma/
• Overture -- http://www.llnl.gov/CASC/Overture/
• PADRE: A Parallel Asynchronous Data Routing Environment -- http://www.c3.lanl.gov:80/~kei/iscope97.1/iscope97.1.html
• PETSc: The Portable, Extensible Toolkit for Scientific Computation -- http://www-fp.mcs.anl.gov/petsc/