High Performance Parallel Programming

1
High Performance Parallel Programming

• Dirk van der Knijff
• Advanced Research Computing
• Information Division

2
High Performance Parallel Programming

• Lecture 18 Mesh Generation

3
Review

• So far weve looked at regular grids
• suitable for infinite problems (with wrap-around)
• many simulations can use regular grids with no
  loss of generality
• easy to program
• simple difference schemes
• efficient solvers
• Not always suitable
• shape of boundaries may be important
• may want to use different accuracies

4
Meshes

• Structured regular connectivity
• Unstructured irregular connectivity

5
Mesh characteristics

• Density high density gives more accuracy but
  computation takes longer (adaptive grids)
• Smoothness large variations in grid
  density/shape can cause numerical
  diffusion/anti-diffusion and dispersion/refraction
  of waves
• Shape boundary layers in fluid flow, and
  convergence requirements in Finite Element
  Methods (later)
• Complexity of domain structured grid generation
  is user-intensive for complex domains,
  unstructured grids are automatic and fast (but
  more computation)

6
Structured meshes

• Rectangular meshes are trivial
• Non-rectangular meshes are meshed using some
  boundary-fitting coordinates
• Overset / chimera grids
• fully boundary-fitted coordinates
• multiblock

7
Overset or chimera grids

• Use a boundary-fitting grid near each boundary,
  and simple rectangular grid in the interior, and
  interpolate between them
• Easy to setup since each grid is local
• Good for moving boundaries but interpolation
  changes as grids move with boundaries

8
Boundary-fitted grids

• Fit a contiguous set of rectangular computational
  domains to a physical domain with curved
  boundaries
• Good accuracy at boundaries
• No interpolation required

9
Boundary fitting grids, single block

• Real space curvilinear coordinates
• Computational space rectangle or cuboid
• Need a mapping between these
• interpolate from block boundary curves / surfaces
• can have numerical errors
• fast - good for 3D
• not smoothing - not ideal for fluid flow
• PDE methods
• use elliptic equations to generate grid lines
• solve PDE in terms of elliptic coordinates
• good smoothing

10
Multiblock grids

• In theory, complex geometries can be mapped to
  one rectangular region (perhaps not
  simply-connected)
• Generally, the physical region is broken into
  regions that have a simple mapping (relatively)
  from the computational grid
• Extensions of single block C, O, H grids
• Composite grids multiblock

11
C, O, H grids

• Computational region need not be a rectangle
• Corners dont need to correspond
• (Need care with difference formulation at these
  points)

12
C, O, H (cont.)

• The region can be multiply connected
• But cuts can be made to transform them to
  simply-connected domains, mappable to a rectangle

13
C, O, H (cont.)

• O - connect the ends together..

14
C, O, H (cont.)

• C - one side is partially joined

15
C, O, H (cont.)

• H - either the multiply connected region shown
  before, or with the interior reduced to a slit

16
C, O, H (cont.)

• C blocks can have coordinate directions change at
  cut
• may need to take account of this
• can be avoided by using a halo
• same difference formula can be used
• haloes must be copied at each iteration
• All of the interior boundaries are mapped to
  exterior boundaries in the computational space

17
Multiblock

• For more complex configurations, use composite
  grids
• Segment the physical region into sub-regions
  bounded by 4 curved sides in 2D, 6 curved
  surfaces in 3D

18
Multiblock (cont.)

• Within each sub-region, generate a coordinate
  system
• The overall coordinate system is formed by
  joining these sub-systems together
• the degree of continuity of grid lines at
  boundaries of sub-regions can be anything from
  complete to none
• if algebraic methods are used, the grid positions
  are fixed on the boundaries and the slopes can be
  controlled
• if elliptic methods are used, the equations can
  be solved in the whole domain, and haloes used
  for complete continuity
• The locations of the boundaries can be fixed or
  left to the grid generator (e.g. elliptic method)

19
Example
20
Example (cont.)
21
Finite Element Methods

• Before we go on to unstructured meshes we will
  look at Finite Element Methods for a few slides
• Consider the problem of finding a function u(x)
  to represent a one dimensional field

22
FEM (cont.)

• Could use a polynomial
• very convenient as polynomials can be
  differentiated and integrated easily
• as the degree is increased the data points are
  fitted with increasing accuracy but high-order
  polynomials can oscillate

23
FEM (cont.)

• To avoid this we divide into sub-regions and use
  low-order polynomials over each sub-region
  (element)
• e.g. for linear polynomials

24
FEM (cont.)

• Continuity
• let u1 and u2 be the values at the end points of
  the element
• define a linear variation between these two
  values by
• where is a normalized
  measure of distance
• define
• such that

25
Basis Functions

• These functions are referred to as the basis
  functions associated with the nodal parameters.
• The basis functions are
  straight lines varying between 0 and 1 as below

j2(x)
1
1
0
x
26
Quadratic Basis Functions

• A quadratic variation of u over an element
  requires three nodal parameters, u1, u2 and u3

j1(x)
j2(x)
j3(x)
1
1
1
1
1
1
0
0
0
x
x
x
0.5
0.5
0.5
27
2D Basis functions

• 4 bilinear basis functions

28
Triangular basis functions

• Based on area (usually)
• consider the ratio of the area of triangle P23 to
  triangle 123

where
, and
29
Triangular element
30
Triangular basis functions

• Similarly
• where
• now
• and we can use
• where

31
FEM mesh requirements

• The mesh must have the following properties
• The union of all elements is an approximation of
  the domain
• Two distinct elements can intersect only as a
  FULL edge or vertex
• Corners and other singular points of the
  continuous domain must be vertices of the
  approximate domain
• Vertices of the approximate boundary must be on
  the continuous boundary
• Note that the continuous domain and the
  approximate domain are equal if and only if the
  first is polygonal.

32
Unstructured Meshes

• Triangulation / tetrahedralization can fit
  irregular boundaries and allow a progressive
  change of element size without distortion
• Linear tetrahedra are not that good for FEM (too
  stiff)
• need lots to give acceptable results
• quadratic hexahedra are much better, but it can
  be difficult to make and all-hexahedral mesh
• quadratic tetrahedra can be as good as quadratic
  hexahedra

33
Mesh generation methods

• Decomposition and mapping
• Earliest method
• Decompose domain (c.f. multiblock) in CAD or by
  hand
• map simple grids onto blocks
• Regular grid or quadtree/octree grid
• Advancing Front
• Delaunay
• Other / combination

34
Grid-based methods

• Triangular / rectangular grid overlaid
• Cropped to domain
• Edge points moved to boundary
• Characteristics
• fast
• good grid in interior, poor grid at boundaries

35
Smoothing

• Laplacian smoothing
• Moves each point to the barycentre of the points
  connected to it
• After 2-5 iterations the process converges
• may need to unsmooth negative areas / volumese.g.

36
Advancing front

• Discretize boundary to create an initial front
• Add triangles / tetrahedra into the domain with
  (at least) one edge / face on the front
• Generally need the domain to be bounded, but can
  advance front into an unbounded domain, stopping
  at some distance (e.g. from the aeroplane you are
  modelling)
• Use mesh parameters defined on some background
  grid to control triangle / tetrahedra generation

37
Advancing front
38
Mesh parameters

• The mesh parameters are the average node spacing
  d, and any stretch parameters
• The spatial distribution of mesh parameters is
  given by a linear interpolation from values
  specified at the nodes of a background mesh of
  triangles
• The background grid can be one triangle enclosing
  the whole domain (for constant mesh size), or a
  coarse grid created interactively, or a grid
  already used for computation (as in adaptive
  meshing using the inverse of the local error to
  specify new local grid size)

39
Boundary discretization

• Boundary is made up of closed loops of curved
  segments
• Compute the length l of each segment
• Interpolate mesh parameters onto each segment
• Calculate the number of sides Ns required per
  segment as
• Node positions li are where the following takes
  integer values
• Effectively working in parameter space of each
  segment

40
Boundary discretization

• Boundary discretization quality is important,
  especially in 3D

41
Element generation example - Peraire

• Choose side from front for base of new element,
  usually smallest side, call it AB
• Interpolate mesh parameters to midpoint of AB
• Find position of ideal point C1
• Inequalities ensure no excessive distortion -
  empirical

42
Element generation (cont.)

• Check for nodes on currentfront within radius
  d1for possible use
• Check sides do notcut other sides oncurrent
  front - if souse one of C2...

43
Performance

• Since this is a greedy triangulation
• storage requirements in 2D are
• time complexity is at most
• Several points in the procedure require searching
  and various data structures can be used to speed
  this up.
• Store background grid as some form of tree
• Store the front as a list ordered by segment
  length
• use a quadtree to find nearby nodes
• etc
• Using these reduces time complexity to

44
3D advancing front

• Much more sensitive to geometric tolerances than
  2D
• The quality of the initial front (surface
  triangulation) is critical - ill-shaped triangles
  produce poor tetrahedra
• Performance could be but
  usually with because trees become
  unbalanced
• The algorithm is basically same as 2D
• Clusters of ill-shaped segments can fail - stall
  / retry
• Modifications to improve performance are possible
• e.g. insert well-spaced points first then use them

45
Delaunay triangulation

• Given a set of points in a plane, a Voronoi
  polygon about point P is the set of points closer
  to, or as close to, P as any other point.

Delaunay triangulation
Voronoi tesselation
circumcircle
46
Delaunay triangulation (cont.)

• The Delaunay triangulation is the dual of the
  Voronoi tesselation, and has the following
  properties
• No point is contained in the circumcircle of any
  triangle
• Maximises the minimum angle for all triangular
  elements(note we would like to minimise the
  maximum...)
• The delaunay triangulation is unique except for
  degenerate distributions of points
• 2D 4 points on a circle

47
Delaunay triangulation (cont.)

• 3D 5 points on a sphere
• 3D 6 points (octahedron) can give an invalid mesh

48
Boyer-Watson algorithm

• Adds points sequentially into an existing
  triangulation
• Add a point
• Find all existing triangles whose circumcircle is
  intersected by the new point
• Delete these triangles to leave a convex (always)
  cavity
• Join the new point to all the vertices of the
  cavity

49
Point generation

• The Delaunay triangulation assumes that the
  points to be triangulated are already known
• points can be pre-generated from the vertices of
  overlapping structured grids - may need filtering
  and smoothing
• points can be generated simultaneously with the
  triangles to improve quality of triangles
• Grids are less smooth than advancing front but
  generation is much quicker
• There can be robustness problems, especially in
  the initial phase when triangles may be highly
  distorted

50
Boundaries

• The Delaunay construction triangulates a set of
  points, and does not necessarily conform to
  imposed boundaries

51
Constraining / conforming

• In 2D, the constrained Delaunay triangulation is
  well defined
• In 3D, no constrained DT is defined, but it is
  possible to recover the boundary edges and faces
  by using face-edge swapping
• Alternatively, it is possible to make the DT
  conform to a boundary curve / surface by adding
  enough points on the boundary to ensure that the
  DT will include the edges on the surface.

52
Face-edge swapping
53
Other methods

• Paving advancing layers of quadrilaterals /
  hexahedra
• need collision rules when fronts merge
• good elements at surfaces
• Whisker weaving and the spatial twist continuum
• Structured near surfaces, unstructured elsewhere
• hybrid methods
• advancing layers (for advancing front)
• advancing normals (for Dalaunay point insertion)
• Recursive decomposition
• etc.

54
High Performance Parallel Programming

• Tomorrow - mesh decomposition