P573 Scientific Computing Lecture 7 Numeric Solution of Partial Differential Equations

1
P573Scientific ComputingLecture 7 - Numeric
Solution of Partial Differential Equations

• Peter Gottschling
• pgottsch_at_cs.indiana.edu
• www.osl.iu.edu/pgottsch/courses/p573-06

2
Overview

• Partial differential equations in physics models
• Components of simulations
• Grids
• Regular
• Irregular
• Discretizations
• Finite difference method
• Finite element method
• Matrices
• Overview linear solvers

Part of the slides with courtesy of UC
Berkeley www.cs.berkeley.edu/demmel/cs267_Spr05
3
Poisson Equation

• Describes stationary (time-invariant) diffusion
  processes
• Or the final state in case of constant boundaries
• For instance stationary heat distribution
• Or describes concentrations
• Etymology Why is a tool for the Poisson equation
  called fishpack? (www.netlib.org/fishpack/)

Example one-dimensional distribution
x0
x1
4
Heat Equation

• Models heat distribution in time-variant
  processes
• As well as the progress time-variant diffusion
  processes

5
Components of Simulation Software

• Grid generator
• Discretization
• Linear solver
• See next lectures
• Visualization
• Not part of the lecture

6
Overview Grids

• Regular
• Composed
• Semi-regular
• Overlapped
• Irregular
• Adaptively refined

7
Regular Grids

• Regularity refers to the neighborhood of grid
  points
• Geometric distances can vary
• Like in right figure, called curvi-linear
• Different grids can be represented by the same
  matrix
• Regular grids lead to regular matrix structures

8
Composed Grids

• Sub-domains are discretized regularly
• Special has to paid at interfaces
• Complex geometries can be treated in an easy way

9
Semi-Regular Grids

• Based on regular grids
• Rectangular
• Cuboid
• Subset of element can be refined
• Refinement is regular, too
• Adaptive refinement possible (more later)

Example Shock waves in fluid dynamicsTool
SAMRAI http//www.llnl.gov/CASC/SAMRAI/
10
Overlapped Grids

• Goal locally finer-grained resolution with
  regular data structures
• Finer resolution necessary for numeric reasons,
  e.g., eddies

Ref.www.grc.nasa.gov/WWW/wind/valid/nlrflap/nlrfl
ap01/nlrflap01.html
11
Irregular Grids

• To discretize arbitrary domains
• Discretization can be locally finer grained
  (within limits)
• 2D mostly triangles
• 3D mostly tetrahedra
• Other elements are possible but more difficult to
  handle in automatic grid generation
• Example triangulation tool Triangle
  (www-2.cs.cmu.edu/quake/triangle.html)

12
Adaptive Grid Refinement

• Local refinement depending on computation
• Calculation of error estimation / indication
  decide where to refine
• Error estimation quantitative statement on error
  magnitude
• Not necessarily very precise
• Error indication only statement that local error
  is (potentially) relatively large
• After local refinement new computation and check
  of result
• Regions needing refinement can move during
  simulation run
• Therefore some systems provide also local
  coarsening

13
Adaptive Grid Refinement Quadrangles

• Split into 4 quadrangles with half length
• Or 8 cuboids of half size in 3D
• Angle-preserving
• Even for curvi-linear grids
• Small angles lead to numerical problems

14
Closure of Quadrangle Grid

• Certain discretization prohibit the split edges,
  i.e. another edge begins in the middle of it
• Quadrangles with refined neighbors are refined
  into 3 triangles

15
Adaptive Grid Refinement Triangles

• Triangles to be refined are split into 4
  triangles
• Which have the same angles
• Triangles with two split edges are also refined
  this way
• Triangles with one split edge are split into 2
  triangles
• The angle opposite the split edge is cut in half
• To avoid further angle bisections these triangles
  are recomposed and split into 4 smaller triangles
  before any further refinement
• By induction smallest angle in refined graph at
  least half size of smallest angle in starting
  graph

16
Discretization

• Differential equation on continuous domain is
  approximated in grid points
• Represented by a system of linear or non-linear
  equations
• Where variables are associated with grid points
• The systems solution approximates the solution
  of the differential equation

17
Finite Difference Method

• Approximation of partial derivatives with
  difference quotients
• Can be estimated with Taylor series, e.g. for u
• It follows

18
Poisson 1D-Discretization

• Approximation of -u(x)f(x) by
• Applied on all xi ih mit 0ltiltN, with N1/h
  yields Aub
• Matrix entries aij only non-zero if i-jlt1
• Tridiagonal matrix
• Right-hand side bih²f(xi) with xi ih for
  0ltiltN
• Boundary conditions also in b b1h²f(x1) )u(xN)
  and bN-1h²f(xN-1)u(xN)

2 -1 -1 2 -1 -1 2 -1
-1 2 -1 -1 2
Grid
A
i 1 2 3 4 5 6 7
19
Poisson 2D-Diskretization

• Same approximation of -?²u(x,y)/?y²
• Matrix for NN-Domain contains (N-1)(N-1)
  unknowns for inner points
• -(?²u(x,y)/?x²?²u(x,y)/?y²)f(x,y) yields to
  Axb
• With bih²f(xi,yi) in inner points
• For points next to boundary u(xi1,yi1) from
  boundary added

20
Flow Equation Explicit Time Scheme

• Initial value problem (IVP)
• u(x, 0) given for all x
• Boundary value problem (BVP)
• u(0, t) and u(5, t) given
• Computation in discrete time steps
• Spatial derivatives approximated in last time
  step
• u(x, t) depends only on constant values
• Computed in a simple matrix vector product
• Requires very small time steps to be numerically
  stable
• Rarely used these days

t5 t4 t3 t2 t1 t0
u0,0 u1,0 u2,0 u3,0 u4,0
u5,0
21
Instability of Explicit Solution
22
Flow Equation Implicit Solution

• Spatial derivatives approximated in current time
  step
• Similar linear systems as for Poisson equation
• Only elements in diagonal are larger
• The shorter the time steps the larger the
  diagonal elements
• The easier to solve

t5 t4 t3 t2 t1 t0
u0,0 u1,0 u2,0 u3,0 u4,0
u5,0
23
Finite Element Method

• Decomposition of the domain into elements
• Parameterized functions are defined on these
  elements
• A criterion is defined that characterizes how
  well the sum of these functions approximate the
  solution of the PDE
• This is evaluated element-wise and leads to
  nn-Matrices where n3, 4, 6, ...
• Depending on the shape of the elements and the
  functions
• These matrices are assembled into a sparse matrix
• Or not, a matrix vector product can be
  implemented directly with the element matrices
• Theoretical background more solid than for finite
  differences
• Especially if solution and/or boundary are not
  very smooth

24
Sparse Matrices

• All linear systems arising from discretization of
  PDEs contain sparse matrices
• Well, except some really rare exceptions
• Sparse Matrix means
• Number of non-zero elements per row (and column)
  is limited, e.g. lt 5 or 20, and doesnt grow for
  finer discretizations (or only very little)
• Number of rows and columns significantly larger
• Grid refinement only increases the number of rows
  and columns not the number of non-zero per row
• Regular grids yield structured matrices
• All non-zeros have same distance to diagonal
• Coefficients can be equal, e.g. 1D-Poisson -1, 2,
  -1 thus matrix doesnt need to stored (can be
  part of the programs)
• Coefficients can be different if other PDE or
  curvilinear grid
• Irregular grids yield unstructured matrices
• Discretization error usually larger
• Applicability much higher

25
Resume

• Solution of partial differential equations are
  approximated in grid points
• For this purpose we lay a grid over the domain
• Based on grid and PDE a linear system is defined
• In easy cases, like Poisson on regular grid, by
  hand
• Can be already contained in the program
• In more complicated cases, FEM on irregular
  grids, program set up the linear system
• Solution of linear systems in the following
  lectures