Unstructured Meshes and Finite Elements in Space Plasma Modelling:

Advanced Methods for Space Simulations, edited by H. Usui and Y. Omura, pp. 111–143.
c TERRAPUB, Tokyo, 2007.
Unstructured Meshes and Finite Elements in Space Plasma Modelling:
Principles and Applications
R. Marchand, J. Y. Lu, K. Kabin, and R. Rankin
Department of Physics, University of Alberta, Edmonton Alberta, Canada
This chapter describes the general principles of unstructured mesh generation in
two and three dimensions, together with the use of finite elements in modelling space
plasmas. Particular attention is given to the requirement of generating meshes that are
aligned along magnetic field lines, in order to properly capture the strong anisotropy
of magnetised plasmas. Example results are given in 2D and 3D, for the dynamics of
shear Alfvén waves in Earth’s magnetosphere.
1
Introduction
Computational space plasma physics has evolved into a mature area of research,
where it can be used as a reliable tool for interpreting observations, or for conducting
numerical experiments to study physical effects that are not readily observable. This
success is due, on one hand, to constant research and development of evermore efficient algorithms and numerical methods and, on the other hand, to the phenomenal
advances in computer hardware. As a result, there now exist several well-developed
numerical approaches for solving a wide range of problems encountered in space,
extending from global to microscopic scales.
Most operational models of space plasmas are based on the following components:
1. a set of equations that is deemed to capture the essential physics of the problem,
2. a discretisation of the equations that is suitable for numerical analysis,
3. a practical set of solution techniques,
4. diagnostics capable of generating model predictions that can be directly compared with observations.
In this chapter, we are primarily concerned with a set of models in which the
plasma is described in the fluid, or macroscopic approximation, and its dynamics is
determined by a set of conservation equations. Important instances of this approach
are based on the ideal MHD approximation [Cabannes, 1970]. The techniques that
will be presented, however, are not limited to this model, and can be equally applied to other fluid descriptions, such as double adiabatic [Chew et al., 1956], Hall
MHD [Huba, 1995], multi-fluid [Freidberg, 1987] or multiple moment [Grad, 1957]
approaches.
111
112
R. Marchand et al.
The discretisation mentioned in component 2 can be obtained in one of two broad
ways:
• The first, and most common, relies on a discretisation of space in terms of
points and cells at which the solution is calculated. The collection of these
points and their connectivity constitutes the discretisation mesh. This is the
basis of the most familiar numerical approaches, such as finite difference, finite
volume or finite element discretisation.
• The second approach relies on a representation of the solution as a superposition of basis functions in a finite dimensional Hilbert space. This approach
leads to the general class of solutions referred to as spectral methods. In this
approach, space is not discretised, but the vector space of possible solutions is.
It is worth noting that finite elements have something in common with both of
these general approaches. As with spectral methods, they rely on a representation
of the unknown functions as a superposition of known basis functions. Contrary
to standard spectral methods, however, which generally use high order interpolating
functions that are non zero on the entire domain, the finite element discretisation
makes use of interpolation functions with a compact support. These are functions that
are non-vanishing over only a small subset of the simulation domain, as prescribed
by the cells on a given mesh.
In what follows, we limit our attention to solution techniques that rely on a discretisation of space. More specifically, we consider the use of unstructured meshes
to solve sets of coupled partial differential equations describing space plasmas in the
fluid approximation. In the remainder of this chapter, we explain the distinction between structured and unstructured meshes. We then introduce some basic concepts
and techniques used in generating a particular type of unstructured mesh; that is,
those based on simplices. We then review finite element discretisation on an unstructured mesh, and point out the importance of mesh alignment along the magnetic field.
Finally, example applications are given in two and three dimensional geometries.
2
The Discretisation of Space: Structured vs. Unstructured Meshes
With the exception of spectral methods, which as stated above, do not require a
discretisation of space, every numerical approach used to model space plasmas relies
on a discretisation of space. This, in turn, is obtained from a mesh consisting of
a distribution of points in space, together with their connectivity. These points and
their connectivity are then assumed to form cells, which cover a prescribed simulation
domain. Before going further, let us introduce the following definitions:
Domain Region of space delimited by an outer boundary and, possibly, by inner
boundaries. In the following, we only consider connected domains; that is,
domains in which it is possible to connect any two interior points with a continuous line entirely contained in the domain. We will also assume that the
domains are classical Euclidian domains and their boundaries are of integer
Unstructured Meshes and Finite Elements in Space Plasma Modelling
113
dimensionality (with the dimension of the boundary being that of the domain
minus one). Fractal domains or boundaries are therefore excluded in the following discussion.
Set of points This will refer to a set of a finite number of points N in a space. For
simplicity these points will be assumed to be contained in a bounded domain.
In particular, in 3D, if xi , yi , z i represent the coordinates of point i, they will
be contained entirely in a rectangular prism xmin ≤ xi ≤ xmax , ymin ≤ yi ≤
ymax , z min ≤ z i ≤ z max , for 1 ≤ i ≤ N . Finally, we will assume that all points
are distinct; that is, xi = x j , yi = y j , z i = z j implies i = j.
Convex envelope to a set of points This is the convex boundary of the domain with
the smallest volume that contains all the points (points contained in the boundary are contained in the domain). We recall that a convex boundary is such that
none of its tangent (line in 2D or plane in 3D) intersects the boundary at points
other than those where it is tangent.
Simplex For a space of dimension D, this is a domain delimited by the convex envelope of D+1 points not contained in the same hyperplane. In 2D, this is a
triangle, and in 3D, it is a tetrahedron.
Cells, or elements We refer to the most elementary components of a mesh as cells or
elements. In the following discussion, we concentrate our attention on meshes
in which cells are made of simplices. We will also use the words cell and
element interchangeably.
Mesh A mesh consists of a set of points, together with their connectivity. In principle, that connectivity need not define cells in the simulation domain but, in
practice, it almost universally does. Here, we assume that the connectivity
makes it possible to define a finite number of disjoint cells, and that every point
in the mesh belongs to one or more cells.
Coverage A mesh consisting of set of N connected points Pi and M cells C j is said to
cover a simulation
domain D if the union of all the cells in the mesh coincides
M
with the domain;
i=1 Ci = D. We recall that our cells are assumed to be
disjoint; i.e., Ci C j = φ whenever i = j and their volume is assumed to be
larger than zero.
Connectivity Given a set of points, a connectivity is a graph in which every point is
connected to a number of neighbouring points, so as to form cells or elements.
There are two generic types of meshes that can be used: structured and unstructured.
2.1 Structured meshes
The simplest examples of structured meshes are rectangular Cartesian coordinates, cylindrical or spherical coordinates in 3D, or polar coordinates in 2D. Structured meshes can also be curvilinear, and they may be orthogonal or not. For example, Fig. 1 illustrates a typical non-orthogonal curvilinear structured mesh in 2D.
114
R. Marchand et al.
1.4
1.2
1
0.8
0.6
0.4
0.2
0
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Fig. 1. Example of a structured mesh constructed on curvilinear coordinates.
The essential feature of structured meshes is that there is a simple and implicit relation between the vertices and their neighbours, or equivalently between the cells and
their neighbours. In the case of a Cartesian 3D mesh, for example, the direct neighbours of a point with indices (i, j, k) will have indices (i ± 1, j, k), (i, j ± 1, k) and
(i, j, k ± 1). In that case, cells which are rectangular prisms, can also be labelled with
three indices, and adjacency between cells can be recognised with similar operations
on the indices.
2.2 Unstructured meshes
Unstructured meshes are all those that don’t fall in the category of structured
meshes. In particular, there is no simple and implicit relation between vertices in
the mesh and the points that they are connected to. In the remainder of this chapter,
we limit our attention on unstructured meshes only and, in particular, on simplexbased meshes; that is, on meshes constructed with triangles in 2D and tetrahedra in
3D. For example, Fig. 2 illustrates a simple triangular mesh in two dimensions. It is
worth noting that this is not the only type of unstructured mesh used in modelling.
For example, an alternative approach is based on a so-called quadtree (2D) or octree
(3D) decomposition of an existing domain [DeZeeuw, 1993]. In this approach, the
elements of an existing mesh may be subdivided into smaller elements by introducing
vertices halfway along their sides. This technique can equally be applied to simplex
meshes or to meshes based on rectangular cells. As an example, the code BATS-RUS [Powell et al., 1999] uses an octree decomposition of a rectangular grid to achieve
local adaptivity. Figure 3 illustrates how this approach can be used to generate a mesh
with variable spatial resolution in a rectangular domain.
2.3 Pros and cons of unstructured meshes
The main advantages of structured meshes are their simplicity and the fact that
they lend themselves to high order difference operators. Their main disadvantage is
that they are more difficult to apply in geometries with multiple connectivity and/or
to problems with irregular boundaries. Structured meshes are also not well suited
Unstructured Meshes and Finite Elements in Space Plasma Modelling
115
10
8
6
4
2
0
0
2
4
6
8
10
12
14
16
Fig. 2. Example of a triangular unstructured mesh in 2D.
1
0.8
0.6
0.4
0.2
0
0
0.2
0.4
0.6
0.8
1
Fig. 3. Example of an adaptive mesh based on a quadtree decomposion .
for local mesh refinement; that is, to refinement of the mesh local to a particular
region of space, without affecting the mesh resolution away from that region. Unstructured meshes, on the other hand, are naturally adapted to represent domains of
arbitrary connectivity, or with boundaries of arbitrary shapes. They can also be refined locally. Their principal disadvantage is their relative complexity and the added
data structures that are required in the simulation codes that use them. Although
high order interpolation functions can be defined on unstructured meshes [Šolín et
al., 2003], the interpolation functions used in most codes are generally limited to C 0
continuity; that is, the interpolated solutions are continuous at the boundary between
elements, but their derivatives there are in general discontinuous. It should be noted
that interpolation functions with a higher level of continuity can be constructed on
triangular meshes. For example, several schemes have been proposed to construct
functions with C 1 continuity on triangular elements [Braess, 2002]. These, however,
116
R. Marchand et al.
are somewhat more complicated to work with, and their use in plasma modelling
is only starting to appear in the literature [Jardin, 2005; Ferraro and Jardin, 2006].
Finally, with certain problems characterised by a strong transport anisotropy, it is essential to align the mesh (structured or not) along the direction in which transport is
the strongest. This point is discussed in more detail in Sec. 5, after the introduction of
the method of finite elements. Although mesh alignment can be achieved with structured meshes for complex topologies, this alignment may become cumbersome and
impractical. Simplex unstructured meshes, on the other hand, can always be aligned
locally. Alignment of meshes based on a quad/octree decomposition of rectangular
elements is in general not practical. In principle, quad/octree decomposition could
be applied to an initially structured and aligned mesh. The resulting refined (and
unstructured) mesh would then be approximately aligned. This approach, however,
would be somewhat complex and to the authors’ knowledge, it is not used in any
existing space plasma physics model.
3
Basic Concepts of Simplex-based Unstructured Mesh Generation
Much development has gone into the generation and use of unstructured mesh
generation. This activity has initially been driven by engineering applications. It was
made possible by the advent of affordable fast computers with extensive memory.
Unstructured grids are now also being used more generally in science, and in applications where it is important to model complex systems while accurately accounting
for their detailed geometry. Excellent reviews and textbooks have been written on the
subject [Frey and George, 1999; De Berg et al., 2000; Gruau and Coupez, 2005; Yang
et al., 2005], and the following section, unfortunately, does not do justice to this very
active field of research and development. It merely attempts to summarise the essential concepts and techniques used by the authors in the context of magnetospheric
plasma modelling.
The problem of mesh generation can be divided into the following five steps.
1. Defining a domain.
2. Defining a metric.
3. Distributing mesh points along the boundaries.
4. Distributing mesh points inside the domain.
5. Connecting these points so as to construct elements with desired properties and
shapes.
These steps, of course, are interdependent, and a good mesh generation strategy is
one that accounts for the ultimate mesh quality at each step of the generation process.
Let us briefly explain each of these steps separately.
Unstructured Meshes and Finite Elements in Space Plasma Modelling
117
3.1
Definition of a domain
Simulation domains are defined by specifying an outer boundary and, possibly, a number of interior boundaries. These boundaries, in turn, are conveniently
parametrised in terms of closed curves (often polylines; that is a sequence of connected straight segments defined in terms of a string of vertices) in 2D, or closed surfaces (often defined in terms of connected triangles) in 3D. These boundaries have
the property that they delimit the simulation domain, and the mesh used to discretise
it. They are therefore domain delimiting. Generally, in a computer model, all domain
delimiting boundaries will be used to impose some form of boundary condition. They
will therefore also be boundary condition bearing. It is worth noting that, in some
applications, it is useful to introduce other boundary-like objects in a mesh that are
not used to delimit the domain. We will loosely refer to such structures as boundaries,
as well. These may, however, be used to specify certain boundary conditions within
the domain. Such boundaries may therefore be boundary condition bearing or not.
A possible use of such interior boundaries, on which no boundary condition is specified, is to guide, or align the mesh along certain prescribed directions. This is the
approach used in the 2D mesh generator Aranea [Marchand et al., 2001], in the creation of meshes aligned along the magnetic field lines. In Aranea, these non domain
delimiting boundaries on which no boundary conditions are specified, are referred to
as guiding boundaries.
3.2 Metric
In mesh generation, the metric is a suitable prescription of size, orientation and
aspect ratio of the elements that constitute the mesh. This can be achieved analytically
or numerically, by using a background mesh on which the metric is discretised, and
from which it can be interpolated. This is referred to in what follows as the domain
metric. In our approach to mesh generation, the metric is assumed to be isotropic;
i.e., without a preferred orientation or stretching of the elements. When modelling
anisotropic magnetised plasmas, however, it is necessary to construct meshes with elements aligned along magnetic field lines. This is achieved by using guiding boundaries, together with a prescribed (isotropic) metric.
3.3 Distribution of mesh points along boundaries
With the domain and the metric being defined, the next step is to distribute mesh
points along boundaries. This distribution is determined from the domain metric as
introduced above, as well as from specific characteristics of the boundaries, such as
curvature, corners or edges in three dimensions. Let us, for now, limit our attention
to mesh distribution in two dimensions. The process then goes as follows:
1. Corners are identified and mesh points are placed on each of them. A corner is
a point where the direction of a line boundary changes suddenly. These points
usually represent locations where important changes in the geometry or in the
physics occur. For open boundaries, the end points are formally treated as corners, and they are assigned mesh points. It is possible for closed boundaries
not to have any corner (e.g., circles or closed smooth curves). In that case, a
starting mesh point may be assigned anywhere along the boundary. A possibil-
118
R. Marchand et al.
ity, then, is to assign a point where the curvature is maximum, or simply, to the
first point used to parametrise the boundary.
2. Mesh points are distributed between corners as prescribed by the metric, and
by the characteristics of the boundary. For example, in Aranea, the number of
intervals created by inserting mesh points between two successive corners A
and B, is given by the closest integer N to
B
1
N =
dl
(1)
δ
A
where δ is the desired separation between adjacent mesh points along that segment of the boundary. This is determined from the domain metric, as well as
from particular properties of the boundary, such as its curvature, or proximity
to other boundaries. Mesh points are then distributed along the segment AB
according to the following prescription,
Pi
N
1
(2)
dl = i.
N A
δ
In 3D, boundary characteristics include corners, edges and curvature. In 2D,
boundaries are conveniently parametrised in terms of polylines, while in 3D, they are
parametrised with triangulated surfaces. In 3D, edges are curves where the unit normal vector to the boundary changes in an abrupt and appreciable manner. Corners are
points where the orientation of the surrounding boundary changes in a discontinuous
and appreciable way. Corners can occur at the intersections between edges, or they
can occur independently of any edge as, for example, at the apex of a cone. Edges and
corners typically correspond to important features in the domain, and it is required
that mesh points be placed at or along them in order to accurately reproduce these
features. When meshing a conical boundary, for example, one would make sure that
mesh points be placed at the summit, as well as along the elliptic edge at the base.
A boundary partition is a subdivision of that boundary into a number of elements
that are disjoint and that cover the entire boundary (their union produces the entire
initial boundary). Partitions are dictated by the problem at hand. In particular, partitions are required whenever different boundary conditions need to be specified on
different portions of a given boundary. In that case, it is important to distribute mesh
points at the junctions of adjacent boundary segments (in 2D) or around closed loops
delimiting partition patches (in 3D) in order to create a mesh with elements that do
not overlap multiple partition segments. That is, every element with a face on the
boundary should have its face entirely in one of the partition segments, and only one.
This is illustrated in Fig. 4, where a valid and invalid distribution of mesh points is
shown for a partitioned boundary.
Once mesh points have been assigned for every corner and around partition segments (patches in 3D), additional mesh points must be distributed between these
points, in order to complete the boundary mesh. The spacing between mesh points
Unstructured Meshes and Finite Elements in Space Plasma Modelling
119
a
10
8
6
4
2
0
0
2
4
6
8
10
12
14
16
10
12
14
16
b
10
8
6
4
2
0
0
2
4
6
8
Fig. 4. Example of a mesh conforming (a) and not conforming (b) to a boundary partition, for a triangular
mesh in two dimensions. The boundary is partitioned so as to accommodate different types of boundary
conditions above the rectangle located betwen x = 5 and x = 9 at the bottom of the figure.
is specified by the domain metric, and also by the local curvature of the boundary.
In 2D, a convenient parametrisation of the curvature is obtained in terms of κ, the
reciprocal of the radius of curvature on the curve. This quantity can be calculated
at each of the parametrising points and interpolated linearly along the curve. In 3D,
the curvature of a surface is parametrised with two principal radii. If the mesh points
distributed on the surface are used to generate an isotropic triangular mesh, then one
120
R. Marchand et al.
would characterise the curvature at each point by the smallest of the two principal
radii. With a representative reciprocal curvature κ, the associated metric would then
be δκ−1 = ακ where δκ represents the desired spacing between adjacent points in a
curve, due to its curvature only. For example, in order to distribute 36 points around
a circular boundary, one would choose α = 5.73.
3.4 Distribution of mesh points in the domain
Given boundaries and a distribution of mesh points along them, there are a number of standard approaches for distributing mesh points within a domain. The main
ones are based on advancing front approaches, and the reader is referred to the specialised literature for more details [Frey and George, 1999]. The approach that we
follow in 2D is that proposed by Rebay [Rebay, 1993], later implemented in a Java
program with full GUI [Marchand et al., 2001]. The strategy that we chose for generating meshes in 3D is somewhat more specialised, and it is dictated by the very
special conditions encountered in space plasmas. Indeed, in order to capture the high
anisotropy that characterises certain processes in space plasmas, it is necessary to
have a mesh that is oriented along the direction of the anisotropy. As explained in
Sec. 5, in order to model the propagation of resonant shear Alfvén waves, or thermal
transport in a magnetised plasma, it is necessary to use a mesh with elements that are
aligned along the magnetic field lines. We use such meshes in two types of simulation
domains, referred to as a) flux tubes and b) toroidal domain.
a) A flux tube is delimited by a simply connected (typically a circle) patch at the
surface of Earth, and its boundary is determined by the extent of the magnetic
field lines that cross that patch boundary. Flux tubes start with a circular cross
section in one hemisphere and, for typical solar wind conditions, map into an
area of oval shape in the other hemisphere. Open flux tubes can, in principle, be
considered as well, but they would have to be terminated at some appropriate
distance from Earth. Flux tubes are useful for simulating phenomena that are
localised to a given set of field lines, such as the propagation of shear wave
impulses. They can also be used to calculate derivatives and metric properties
in the vicinity of given field lines [Rankin et al., 2004, 2005].
b) A toroidal domain is constructed as flux tubes, except that its projections on
Earth surface are entire annuli surrounding the magnetic poles. As flux tubes,
their extent is determined by the magnetic field lines that cross the annular
boundaries.
3.5
Connecting a given set of points to form simplices
The connectivity criterion adopted here is that of Delaunay [Delaunay, 1934].
With that connectivity, each simplex element of the mesh is such that its circumscribing ball (circle in 2D, sphere in 3D) contains no other points in the mesh than the
vertices of the simplex on which it is constructed. Given a set of points, degeneracies
notwithstanding, a mesh satisfying this connectivity criterion is known to exist and to
be unique. Degeneracies occur when more than the number of points needed to define
a simplex appear on the boundary of the same ball. In that case, the connection of
Unstructured Meshes and Finite Elements in Space Plasma Modelling
121
these points is not unique. For most of the discussion that follows, we assume that the
distribution of mesh points is not degenerate. Given a set of points, this connectivity
criterion tends to produce meshes with elements that are “as equilateral as possible”.
For many practical problems, this is the mesh with the best properties; that is, it leads
to the smallest discretisation errors.
The underlying concept of Delaunay connectivity have been exposed many years
ago by Dirichlet [Dirichlet, 1850], Voronoi [Voronoi, 1908] and Delaunay [Delaunay,
1934]. In their seminal articles, Dirichlet, and later, Voronoi, considered the properties of cells around each point Pi with the property that: Given a set of points Pi ,
Ti is the region of space made of all the points P from where the closest mesh point
is Pi . In two (three) dimensions, this defines a set of cells which are delimited by,
generally, irregular polygonal (polyhedral) boundaries. This is illustrated in Fig. 5
for a distribution of points in two dimensions. The partition of a domain in terms of
such cells is referred to as a Dirichlet or Voronoi tessellation, and the cells themselves
are referred to as Voronoi cells. It was later noted by Delaunay that elements with
minimal circumscribing balls (i.e., containing no mesh point other than the vertices of
the elements themselves) could be constructed by connecting all mesh points within
Voronoi cells that share a common boundary. A mesh generation algorithm based on
the prior construction of a Dirichlet tessellation, however, would be numerically very
cumbersome, and not useful in practice. The first efficient algorithm for constructing
Delaunay cells from a given distribution of points was proposed independently by
Bowyer [Bowyer, 1981] and Watson [Watson, 1981], and it is often referred to as the
Bowyer-Watson algorithm. It consists of recognising that a Delaunay mesh can be
constructed incrementally by successively adding mesh points to an existing mesh,
and successively updating the connectivity. The starting point of this algorithm consists of a primitive mesh composed of a single super-element enclosing all the mesh
points in the domain. The algorithm then relies on the uniqueness of the Delaunay
connectivity (ignoring degeneracy); hence, on the independence of the order in which
the points are added to the mesh. Another important factor, affecting numerical efficiency is the fact that the corrections to the mesh, when a new point is added, are
generally local to that point. That is, the number of elements that need to be modified
is roughly constant, independent of the number of elements in the mesh. This results
in a computation time that scales linearly with the number N of mesh points to be
inserted. The reader is encouraged to refer to Boyer and Whatson’s original articles
for more details. This algorithm has been implemented in a number of efficient computer programs. An article by Sloan [Sloan, 1987], for example, provides a detailed
explanation of the method, together with an efficient FORTRAN implementation for
an arbitrary distribution of points in a plane.
The task of incrementally constructing a Delaunay mesh in three dimensions is
conceptually the same as in two dimensions, despite some technical differences. For
example, when adding a point to the interior of a convex triangular mesh in two dimensions, there always results an increase in the number of elements by two. Specifically, when a point is added to an existing triangle, that triangle is deleted and three
new triangles are created, as illustrated in Fig. 6. The connectivity of the mesh sur-
122
R. Marchand et al.
12
10
8
6
4
2
0
0
2
4
6
8
10
12
14
16
Fig. 5. Illustration of Voronoi cells (dotted) and their dual Delaunay elements (solid) in two dimensions.
rounding the new triangles is then updated by recursively checking the Delaunay
criterion with adjacent triangles that are opposite the added point. These triangles,
together with the modified triangles, form quadrilaterals, and updating the mesh connectivity amounts to swapping the diagonals, or not, as prescribed by the Delaunay criterion. When a diagonal is swapped, the two original triangles that made the
quadrilateral are deleted, and two new ones are created. The total number of triangles, however, remains constant. This process is repeated recursively until all the
elements adjacent to the new triangles satisfy the Delaunay criterion. This procedure
is described in more detail by Sloan [Sloan, 1987].
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
Fig. 6. When a point is added to a two dimensional mesh, the triangle containing that point is deleted and
three new triangles are created. This results in a net increase in the number of elements in the mesh by
2.
Unstructured Meshes and Finite Elements in Space Plasma Modelling
123
In three dimensions, when a point is added to a convex Delaunay mesh, the result
is not a fixed increase in the number of tetrahedra. The element containing the added
point is deleted, as in the 2D problem, but a cavity surrounding that point has to be
constructed by removing all the elements with a circumsphere containing the added
point. The point is then connected to the faces of that cavity to form new tetrahedra.
In general the number of new tetrahedra formed, when a new mesh point is added, is
not fixed as in the two dimensional case, and it will vary from point to point. In three
dimensions, it is not possible to design an algorithm for checking and correcting the
Delaunay connectivity, that is equivalent to a simple diagonal swap in two dimensions. Finally, in the process of constructing the cavity, care must be taken in order to
avoid the formation of so-called sliver elements; that is, elements with nearly vanishing volumes. These elements tend to occur with (nearly) degenerate distributions of
mesh points (where more than N+1 points are (nearly) on the same circumsphere).
Slivers typically occur because of roundoff errors in floating point arithmetic, and because of the resulting inconsistencies in assessing the circumsphere criterion between
neighbouring elements. A detailed discussion of 3D mesh generation strategies is
beyond the scope of this chapter, and the interested reader is referred to the excellent
review by Frey and George [Frey and George, 1999] and references therein.
3.6
Data structure
As mentioned earlier, one difficulty in working with unstructured meshes, has to
do with the added overhead related to a more complex data structure. With simplex
cells, the mesh is specified by two lists. One gives the coordinates of each vertex in
the mesh, and another one gives the ordered vertices that every element is made of.
For example, the simple two dimensional mesh illustrated in Fig. 7 would be specified
by the following data structure:
Fig. 7. Simple representation of a triangular domain in terms of four triangular elements. Vertices are
numbered from 1 to 6 (small numbers) and elements are numbered at cell centres, from 1 to 4 (larger
numbers).
124
R. Marchand et al.
Coordinates
0.0
1.0
0.5
2.0
1.5
1.0
Elements
1
2
2
3
nv=6
0.0
0.0
0.5
0.0
0.5
1.0
ne=4
2
4
5
5
3
5
3
6
adjacency
−1
−1
2
3
3
−1
4
−1
−2
3
1
−2
This represents a triangular domain subdivided into four triangular elements. The
first list gives the x, y coordinates of the vertices in the mesh (6 in this case). The
second list gives, for each element, the indices (referring to the first list) of the vertices of which the element is made. By convention, in plane mesh generation, vertices
are ordered counterclockwise. For example, in this mesh, element 4 is made of vertices of indices 3, 5 and 6 in the first list. The adjacency list indicates the elements
or boundaries to which a given element is adjacent. In our programs, we identify
boundaries with a negative index. Thus, in this example, there are two boundaries on
which, possibly different boundary conditions could be imposed. Those are boundary
1 (the polyline defined from vertices 1-2-4-5-6) and boundary 2 (polyline 6-3-1). The
adjacency list refers to sides 1, 2 and 3 of a given element where, by convention, sides
1, 2 and 3 are defined by vertices 1-2, 2-3 and 3-1 respectively. In our example, side
1 of element 4 (defined by vertices of indices 3-5) is adjacent to element 3, side 3 is
adjacent to boundary 2, etc. The output generated by some mesh generators may or
may not include the adjacency list. If it does not, it is usually necessary to reconstruct
this list in the finite element code.
The data structure used in a three dimensional mesh is similar. It provides a list of
vertices defined in terms of their x, y, z coordinates, followed by a list of elements,
and possibly of adjacency for each element. In 3D, it is obviously not possible to
order vertices counterclockwise. It is nonetheless important to have a specific order
that will serve to determine whether a point is interior or exterior to an element. In
our work, we order vertices so that, as seen from vertex 1 of an element, a rotation
from vertices 2, 3 and 4 is clockwise. Adjacency also has to be defined differently in
3D. While referring to sides in a triangular mesh, the adjacency list refers to faces
of a given element (tetrahedron) in 3D. In our mesh generation and finite element
codes, we associate the four faces of a tetrahedron with the vertices to which they are
opposite. Thus, for example, if an element were defined by
12
27
3
900
-4
3
64
2019
its face opposite to the first vertex (with index 12); that is, the face made of vertices
with indices 27, 3 and 900, would be adjacent to boundary 4. Similarly, the face
Unstructured Meshes and Finite Elements in Space Plasma Modelling
125
opposite the second vertex (with index 27) would be adjacent to element 3.
We have not yet discussed the solution of systems of partial differential equations with
finite elements, but it is appropriate at this point, to mention how their solutions may
be represented. Given an unstructured simplex mesh as described above, a solution
simply consists of values of unknown functions (possibly several) at each mesh point.
This is conveniently provided by an ordered list of these values for each dependent
variable. For example, in a solution of a two-dimensional problem involving the
temperature and density as dependent variables, the solution would be presented as
follows:
Dependent variable: density
1.0000
1.1000
...
Dependent variable: temperature
44000.00
47000.00
...
Such solution fields may be appended to the file containing the definition of the mesh,
or they may be contained in separate files. The data structure presented here for the
mesh and for the fields of dependent variables is sufficient for completely defining a
solution, and for visualising it. The exact format of these fields will vary depending on
the specific visualisation tool being used and, in some cases, the choice of the optimal
output format may be dictated by the choice of a specific visualisation software.
3.7 Performing miscellaneous tasks
The construction of unstructured meshes requires performing several tasks such
as calculating volume elements, interpolating functions defined at grid points, and
finding the location of a point in the grid. These operations are generally straightforward, and can be understood from simple principles of geometry. As an example, let
us see how the index of an element containing a given point in space can be found. To
start with, let us determine whether a point is interior to an element, on its boundary
or exterior to it. For that purpose, the element t0 being considered is parametrised by
three vector positions r0 , r1 and r2 . These vectors define the positions of the vertices
0, 1 and 2 of the triangle. As stated previously, we follow the common convention in
plane triangular mesh generation, and orient element vertices counterclockwise. Let
r be the position of point P in our system of coordinates. Referring to Fig. 8, it is
clear that P is interior to t0 if and only if the z components of the vector products
(
r1 − r0 ) × (
r − r0 ), (
r2 − r1 ) × (
r − r1 ) and (
r0 − r2 ) × (
r − r2 ) are all positive. If P is
on one of the boundaries but not on one of the vertices, however, one of these vector
products will vanish and the other two will be positive. If P happens to coincide exactly on one of the vertices, then two of the vector products will vanish and the other
will be positive. Finally, if the point is outside t0 (P’ in Fig. 8, for example), then one
of the vector products will have a negative z component.
126
R. Marchand et al.
Fig. 8. Illustration of the strategy used to find the index of an element containing a point P in the plane.
An efficient approach for finding the element containing a point consists of using
the conditions described above, to ”march” in approximately a straight line, from
an arbitrary element, toward the point until it is found to be interior to an element.
This can be achieved by successively computing the vector products given above
and, whenever a cross product is found to be negative, by moving to the element
adjacent to the corresponding side. This is repeated until the point is found to be
interior to the element or on a boundary. For example, in Fig. 8, if point P’ were
considered instead of P, the z component of the cross product (
r2 − r1 ) × (
r − r1 )
would be negative. The next element to be tested would then be t1, which is adjacent
to side 1-2 of t0, and so on. One point worth noting here is that this strategy will
only work if the mesh has a convex boundary and if it is simply connected. This
remark applies equally well to unstructured meshes in three or more dimensions.
If the simulation domain considered is concave or if it is multiply connected, the
marching scheme described above could lead to elements on the boundary with no
neighbouring element in the direction of the destination point. For that reason, while
constructing a mesh, the working mesh has to be simply connected and convex. If
the final mesh to be constructed is multiply connected or non convex, the exterior
elements are removed from the mesh at the very end of the construction process.
4
Discretisation of Governing Equations
In this section we discuss how partial differential equations can be discretised on
an unstructured grid. In doing so, we concentrate on the method of finite elements
because this is the approach used in our models, and because we feel that it is a relatively straightforward one to apply with unstructured grids. Finite elements have
also been studied extensively over many years, and their mathematical properties are
now well established. Before proceeding, however, it is useful to mention the main
Unstructured Meshes and Finite Elements in Space Plasma Modelling
127
discretisation methods that can be used with unstructured grids, and point out some
of their features.
1. The most important alternative to finite elements is, in our view, based on a
finite volume discretisation. This form of discretisation is particularly useful
when working with governing equations written in the conservative form; that
is,
∂U
+ ∇ · F = 0,
(3)
∂t
where U may represent a multi-component field, and F is the associated (multicomponent) flux vector. The finite volume discretisation then relies on an integration of Eq. (3) over cell volumes, on the use of the divergence theorem,
and on accounting for the detailed fluxes in and out of every cell at their boundaries. This approach is conceptually simple and it has the advantage of ensuring
detailed conservation laws on a cell-to-cell basis. It can be used in conjunction with various methods for calculating fluxes at cell interfaces, such as ones
based on high order interpolations, or a combination of high and low order interpolations so as to preserve positivity and monotonicity [Dendy, 2002]. When
finite volume discretisation is applied to unstructured simplex-based meshes,
the basic cells being considered are the dual cells or Voronoi cells, and not the
simplices (triangles or tetrahedra) themselves [Whitaker, 1988; Michev, 1996;
Frink, and Pirzadeh 1999].
2. Another method relies on a finite difference approximation of the partial differential operators on a given grid [Erlebacher, 1985]. While conceptually simple
with structured orthogonal grids, this approach requires more care when working with simplex unstructured meshes. The resulting finite difference operators
and the approximate finite integral operators then need to be constructed carefully in order to ensure consistency, conservation laws and, where applicable,
positivity and monotonicity. Finite differences are being used to solve a number of problems in physics and engineering, but they are seldom encountered
in complex problems formulated with unstructured grids.
3. With finite elements, it is the solution space that is discretised, and not the differential operators themselves, as with finite differences. This is achieved by
expressing unknowns as a superposition of finite basis functions with a compact support. The finite element approach does not lead to detailed conservation schemes on a cell-by-cell basis, as with finite volumes. It does, however,
preserve conservation equations in the integral sense. In particular, because finite elements involve spatial integrations over simplex elements (not over their
duals as with finite volumes), they lead to time variations of physical quantities
integrated over elements, that are consistent with fluxes across their boundaries
through conservation laws [Lomax, 1977; Kuzmin, 2003].
128
R. Marchand et al.
Each of these methods has its advantages and disadvantages. In modern computational physics, both the finite volume and the finite element approaches are essentially
on par with respect to the types of problems that they can be applied to, as well as
for the quality of the solutions that they produce. In the remainder of this section,
we limit our attention to finite elements, with which we have more experience. No
attempt here is made to give a comprehensive description of this method. We rather
concentrate on a simple example diffusion problem, which is sufficient for explaining
the basic workings of the method, and for pointing out some technical considerations
related to the use of unstructured meshes in strongly anisotropic media.
Let us consider the problem of modelling anisotropic diffusion for a scalar field u
as described by
∂u
(4)
= ∇ · (κ − κ⊥ )b̂b̂ + κ⊥ I · ∇u ,
∂t
where b̂ is a unit vector along a prescribed vector field (in the examples considered
below, that will be the magnetic field), κ and κ⊥ are diffusion transport coefficients
in the direction parallel and perpendicular to b̂ respectively. This model equation
provides a good approximation to electron thermal transport in magnetised laboratory experiments such as tokamak fusion experiments. In these devices, electron
thermal transport is primarily determined by diffusion along magnetic field lines,
and the parallel termal diffusion coefficient κ is typically many orders of magnitude
larger than the diffusion coefficient in the direction perpendicular to the field lines
κ⊥ [Marchand and Simard, 1997]. This problem is also relevant to modelling shear
wave propagation (Alfvén or whistler) in a magnetised plasma, in which propagation
is predominantly along magnetic field lines. The solution procedure is similar to that
followed with standard spectral (or more generally, weighted residual) methods. The
unknown solution u(
r , t) is expressed as a superposition of known basis functions
Ni , i = 1, ..., M as
M
u(
r , t) =
a j (t)N j (
r ).
(5)
j=1
This expression for u is substituted in the transport equation (4) and the resulting
equation is then successively multiplied by M independent weight functions i (
r ),
and integrated over the simulation domain to yield a set of M coupled ordinary differential equations for the coefficients ai (t). In the Galerkin formulation, the weight
functions are the same as the interpolation functions; that is, i (
r ) = Ni (
r ). After
integrating by parts, this becomes
M ∂a j (t)
j=1
∂t
·∇ N j ) − a j
dωNi N j + a j
dω((κ − κ⊥ )b̂ · ∇ Ni b̂ · ∇ N j + κ⊥ ∇ Ni
= 0.
dγ b̂ · n̂ Ni (κ − κ⊥ )b̂ · ∇ N j + Ni κ⊥ n̂ · ∇ N j
(6)
Unstructured Meshes and Finite Elements in Space Plasma Modelling
129
The last integral is carried over the boundary of the domain, and n̂ represents
the unit vector perpendicular to , and pointing outward. This integral is used to
impose boundary conditions of the Neumann, or mixed type. In finite elements, the
first integral appearing in Eq. (6) is referred to as the mass matrix, while the second is
a generalisation of the stress matrix. The time derivatives of the coefficients a j may
in turn be discretised with finite differences. For example, with
− a kj
a k+1
∂a j (t)
j
,
∂t
t
(7)
where the superscripts indicate the time step of the coefficients, Eq. (6) then leads to
a system of M equations for the M unknown coefficients a j of the form
k+1
M
a j − a kj
j=1
t
·∇ N j ) −
a hj
dωNi N j + a hj
dω((κ − κ⊥ )b̂ · ∇ Ni b̂ · ∇ N j + κ⊥ ∇ Ni
dγ b̂ · n̂ Ni (κ − κ⊥ )b̂ · ∇ N j + Ni κ⊥ n̂ · ∇ N j
= 0. (8)
In Eq. (8), if h = k, the scheme is said to be explicit. If h = k + 1/2, or h =
k + 1, on the other hand, it is said to be centred (or semi-implicit), or fully implicit,
respectively. We note that, in the special case where κ and κ⊥ are independent of
the unknown function u, then, this system of equations is linear. In the more general
case (and this usually applies when modelling plasmas) where these coefficients are
functions of u, then the equations are nonlinear and they generally require suitable
iterative solution techniques.
In Eq. (6), if the interpolation functions Ni were chosen to be analytic functions,
such as sine and cosine, or a set of independent polynomial functions defined over
the entire domain, this approach would reduce to a spectral method. In that case, the
system of equations (8) would generally be described by a full matrix A. That is, the
system of equations would be of the form
AX = B,
(9)
where the elements of the matrix A would mostly be non zero. The distinctive advantage of the finite element method over standard spectral methods comes from the use
of interpolating functions Ni with compact support. Specifically, Ni is non zero in the
elements that are connected to mesh point i, and it vanishes in every other element.
As an example, for piecewise linear interpolation functions, Ni is equal to 1 at mesh
point i, and it decreases linearly toward other grid points to which i is connected. It
can be seen that such a superposition of interpolating functions will lead to a representation of the solution with C 0 continuity; that is, the represented solution will be
continuous, but its derivative will, in general, be discontinuous at cell boundaries. An
obvious consequence of the interpolation functions having a compact support is that
they lead to global matrices (matrix A in Eq. (9)) that are sparse. Indeed, the only
130
R. Marchand et al.
elements ai j of A that don’t vanish are those corresponding to nodes i and j that are
directly connected in the mesh. Unfortunately, the main strength of finite element discretisation is also at the root of one of some of its difficulties. Using compact support
on a triangular mesh implies that it is generally not practical to construct high order
interpolating functions, with more than C 0 or C 1 continuity. It is however relatively
straightforward to define high order interpolation functions within a single element.
This can be achieved by introducing additional mesh points within cells. The resulting interpolation and its derivatives are then continuous within each of the elements.
The derivatives of these functions, however, typically remains discontinuous across
elements. The interested reader is referred to a recent book by Šolín et al. [Šolín et
al., 2003] for more information. This can be a limitation of the approach in problems
for which high order of interpolation functions are required. In such problems, it
may be necessary to define intermediate derivatives in terms of additional dependent
variables. As mentioned earlier, simulation codes are now being developed with C 1
continuous interpolation functions to model magnetized laboratory plasmas [Jardin,
2005; Ferraro and Jardin, 2006], and the same approach should be applicable to space
plasmas. These techniques are of interest because they can be used directly to solve
systems of partial differential equations with differential operators of higher orders
(up to order four when combined with integrations by parts) without the need of introducing additional dependent variables. These methods, however, are somewhat more
complex to implement and, at present, they are seldom used in practice. It is also the
authors’ understanding that their use in magnetized plasma modelling has, thus far,
been limited to triangular meshes that are effectively structured and rectangular.
Before closing this section, it must be mentioned that the particular formulation
of finite elements considered here, that based on a Galerkin discretisation, is by no
means unique. As pointed out in the previous paragraphs, there are many possible
choices of basis functions Ni , and weight functions j . Using the Galerkin discretisation with linear interpolation functions led to a mass matrix with non zero
off-diagonal elements. Moreover, the formulation of diffusion problems as in Eq. (4)
typical requires a fully implicit formulation (h = k + 1 in Eq. (8)) in order to ensure numerically stability for acceptable timesteps t. This implies that the solution
of time dependent problems require the solution of a matrix equation. For large 3D
problems, such solutions are considerably more computer intensive than explicit sok
lutions, in which each unknown a k+1
j can be written explicitly in terms of a j without
resorting to any matrix inversion. Such explicit solutions schemes can readily be obtained from the finite volume discretisation scheme. They can also be obtained with
finite elements, albeit, with different choices of basis Ni and weight j functions.
5
The Importance of Mesh Alignment
In the previous section we introduced the rudiments of the finite element method
with a simple anisotropic diffusion equation. As in this example, the finite element
discretisation of a set of partial differential equations generally results in a sparse
system of equations. One crucial aspect of plasma modelling, which is not usually
encountered in other fields, is related to the very strong anisotropy that characterises
Unstructured Meshes and Finite Elements in Space Plasma Modelling
131
magnetised plasmas. In laboratory fusion experiments, for example, the electron heat
thermal conductivity in the direction parallel to the magnetic field is typically much
larger than that in the perpendicular direction. In space plasmas, the propagation of
shear Alfvén waves is also highly anisotropic, as it takes place along magnetic field
lines. In such cases, it is important to have a mesh that is aligned in the direction
along which transport is the strongest, as having a mesh that is not aligned would
lead to unphysical transport in the perpendicular direction. This is best seen from the
example considered above, if we set κ⊥ = 0 and keep κ finite. Ignoring curvature
for simplicity, it can be seen from Eq. (6) that the diffusion operator will couple two
vertices of a given element only if they are along the same field line. This follows
from the definition of the linear interpolation functions Ni , and the fact that, for each
element, ∇ Ni is perpendicular to the side (or to the face in 3D) opposite mesh point
i. Specifically, the second integral in Eq. (6) vanishes if the segment joining i and j
is not parallel to b̂, and it does not vanish otherwise. The situation is more complex
if curvature of the field b̂ is taken into account. In this case, the second domain
integration in Eq. (6) may not vanish when i and j are on different field lines, even
if the mesh is perfectly aligned along the field. In such situations, when κ κ⊥ ,
it is therefore necessary to take special precautions in the construction of the global
matrix, so as to prevent any of the strong parallel diffusion terms from spuriously
contributing to the (much weaker) perpendicular direction. This can be achieved
by artificially setting to zero, the terms involving κ , that would otherwise couple
mesh points located on different field lines. Alternatively, the magnetic field (or more
generally, the direction of strong anisotropy) can be treated as piecewise linear in
every element; with its direction coinciding with that of the edge that is known to be
field aligned.
At this point, it worth noting that aligned meshes can readily be constructed with
simplex-based elements. This is not the case, however, with other types of unstructured meshes. For example, meshes based on quadtree (2D) or octree (3D) domain
decomposition are generally not aligned along the magnetic field, or along a general direction of strong transport anisotropy. As a result, while such meshes can be
refined locally, and used efficiently to model processes with very different spatial
scale lengths, their use would lead to unacceptable numerical diffusion in problems
characterised by strong transport anisotropies.
6
Numerical Solution
Let us now turn to some general remarks concerning the numerical solutions of
the equations that result from the finite element discretisation presented above. In
their simplest form, these equations result in a large sparse systems of linear equations. Moreover, as a result of the use of unstructured meshes in the model, the sparse
matrix itself has no regular structure. For example, Fig. 9 illustrates the locations of
non zero elements of the matrix that would result from discretising Laplace’s equation
with a five-point stencil on a 4 × 4 cartesian grid.
The matrix is seen to be pentadiagonal, with two non zero diagonals four rows
below and above the main diagonal. This regular structure, together with the prop-
132
R. Marchand et al.
16
14
12
10
8
6
4
2
2
4
6
8
10
12
14
16
Fig. 9. Illustration of the structure of a banded matrix resulting from the discretisation of Laplace’s
equation on a four by four cartesian mesh.
erties of the matrix itself, can the be exploited to construct various efficient solution
algorithms. If Laplace’s equation is discretised on an unstructured grid, however, the
non zero elements of the resulting matrix would look typically as shown in Fig. 10.
The lack of regularity in this matrix has some consequences on the types of solution techniques that can be used, and on their numerical efficiency. The construction
of efficient solution algorithms for unstructured matrices such as the one illustrated in
Fig. 11 is still very much an area of research, and it is beyond the scope of this short
review. The interested reader is referred, for example, to an excellent monograph by
Saad [Saad, 2003] and to numerous recent articles in specialised periodicals for more
information. There are basically two types of approaches for solving such systems,
referred to as direct and iterative.
Direct solvers often rely on performing a LU decomposition of the system of
equations [Li and Demmel, 2003; Gupta, 2002, 2006; Al-Kurdi and Kincaid, 2006].
They have the advantage of being robust and, assuming that the matrix is not (nearly)
singular, of always producing a solution. Unfortunately, they tend to require a considerable amount of computer memory and, for large systems of equations, they require
too much computing resources (both memory and time) to be practical. It should be
noted, in passing, that the amount of memory required to perform a complete LU decomposition depends critically on the way mesh points are numbered. This follows
from the fact that, in the global matrix, non zero off-diagonal elements correspond to
connected mesh points in the grid. If the mesh points are numbered in such a way that
connected points differ in their indices by only a small number, then connected points
will correspond to non zero elements close to the diagonal in the global matrix. If, on
Unstructured Meshes and Finite Elements in Space Plasma Modelling
133
70
60
50
40
30
20
10
0
0
10
20
30
40
50
60
70
Fig. 10. Illustration of the structure of a matrix that results from the discretisation of Laplace’s equation
on the unstructured mesh shown in Fig. 2. In this case, the matrix has a profile of length 1207.
the contrary, connected mesh points have indices that are very different, this will result in non zero elements being far from the main diagonal. The numerical efficiency
of a direct LU decomposition depends on the proximity of non zero elements to the
main diagonal. This follows from the fact that, in carrying this decomposition, it is
necessary to allocate storage for many elements of the global matrix A that are initially zero. For example, with the so-called skyline storage, allocation must be made
for all entries up to the last non zero element above every diagonal element, as well as
up to the last non zero element on the left. This typically includes many zero elements
in the original matrix. These elements, however, are generally non zero after the LU
decomposition has been performed. Thus, for matrices with non zero elements that
are scattered far from the main diagonal, it is necessary to allocate more storage than
for those in which non zero elements are tightly located near the diagonal. Various
algorithms have been designed to minimise the storage required when working with
unstructured meshes [Sloan, 1989; Löhner, 1998; Scott, 1999]. In our two dimensional simulations, we typically use a direct LU solver with skyline storage. For this,
Sloan’s profiler typically produces matrices with the smallest profiles (the profile is
the length of the array required to store the elements of the LU decomposition of
the matrix). We note that some solvers come with their own integrated profiler, or
reordering algorithms. This is the case, for example, with WSMP [Gupta, 2006] and
with most optimised solvers provided by high performance computer vendors. When
using such solvers, users need not be concerned with the specifics of vertex ordering.
Iterative solvers have the advantage that they require considerably less storage
134
R. Marchand et al.
70
60
50
40
30
20
10
0
0
10
20
30
40
50
60
70
Fig. 11. Illustration of the structure of a matrix that results from the discretisation of Laplace’s equation
on the same unstructured mesh as used in Fig. 10, but for which mesh points have been renumbered so
as to minimise the matrix profile. The renumbering led to a reduction in the profile from 1207 to 490.
For reference, full matrix storage, in this case, would require memory alloation for 68 × 68 = 4624
matrix elements.
and computer time than their direct counterpart. This is particularly true for large
(3D) problems. They are also less sensitive to the ordering of mesh points. Their
main disadvantage is in their, sometimes, slow convergence, or lack thereof. The
design of efficient and robust iterative solvers for general unstructured matrices is
still an area of active research, and and their use tends to rely on adjusting various
parameters that are problem dependent, and require some trial and error. The reader
is invited to refer to the specialised literature for more details [Li, 2003; Saad, 2003;
Chi et al., 2005; Garcia, 2005]. It should be noted that major computer vendors
provide scientific libraries that are optimised for their computing architecture. These
typically include a variety of direct and iterative solvers for sparse linear systems of
equations. It is our experience, however, that for the larger problems, these may not be
adequate. In practice, in our larger 2D and 3D simulations, we successfully used the
GMRES iteration scheme with Saad’s incomplete LU preconditioning [Saad, 2004]
where generic vendor solvers had failed. We note, in closing, that other types of
iterative solutions are possible, such as those based on a multigrid approach [Briggs,
2000; Šolín et al., 2003]. These methods, on the other hand, are not optimally suited
for all types of meshes and problems, and no attempt is made here to explain this
potentially important class of solution methods.
Unstructured Meshes and Finite Elements in Space Plasma Modelling
135
7
Example Applications in 2D
Our 2D model has already been used to investigate shear Alfvén wave resonances
in the magnetosphere [Lu et al., 2003, 2005a, 2005b]. These calculations were done
by solving the so-called reduced-MHD equations for low frequency plasma. Here, we
present a nonlinear MHD calculation of the interaction of shear Alfvén field line resonances and compressional modes in the inhomogeneous magnetospheric plasma. For
simplicity, we neglect the effect of wave dispersion and assume an isotropic plasma
pressure. Under these assumptions, the interaction of shear Alfvén waves and compressional modes is governed by the following set of ideal MHD equations:
∂B
= −∇ × E,
∂t
∂E
µ0 j + ε0
= ∇ × B,
∂t
E = − (v × B) ,
(10)
(11)
(12)
∂ρ
+ ρ∇ · v + ∇ρ · v = 0,
(13)
∂t
∂v
ρ
+ ρv · ∇v = −∇ P + j × B,
(14)
∂t
∂P
+ γ P∇ · v + v · ∇ P = 0.
(15)
∂t
It should be noted that these equations are not linearised. Also, the displacement
current is included in the Eq. (11). The inclusion of this term from Maxwell’s equations becomes important in regions where the Alfvén velocity approaches the speed
of light.
To illustrate the procedure, we consider a dipole magnetic field and assume a
stationary background plasma. An extension of this analysis to a more complicated
geometry with dispersive effects is straightforward [Lu et al., 2003; Rankin et al.,
2004]. In this example, a finite element discretisation is made between magnetic
shells starting at L = 5 and ending at L = 11. We use an unstructured, triangular and
orthogonal mesh aligned with the magnetic flux surfaces to avoid spurious perpendicular diffusive transport. A constant background density n 0 = 106 m−3 is used for
simplicity. In our simulation, the shear Alfvén wave is excited by driving a magnetic
perturbation over a finite area located near the equator, with a frequency matching the
eigenfrequency on the shell L = 7.7.
Initially, shear Alfvén waves are excited over a broad region of space around the
resonant shell L = 7.7. After a few periods, however, only waves in the immediate
vicinity of the resonant shell grow to a large amplitude, which results in a tightly
localised mode structure near L = 7.7. The localisation of this structure is characterised by a scale length which monotonically decreases with time.
Figures 12 and 13 show respectively the parallel current density and perpendicular electric field after driving the mode for five periods. We see from Figs. (12) and
R. Marchand et al.
Z (Re)
136
X (Re)
Fig. 12. Parallel current density calculated in the vicinity of the ionosphere, associated with the driven
shear Alfvén wave, after five periods.
Fig. 13. Perpendicular electric field associated with the driven shear Alfvén wave, after five periods.
(13) that the mode has already developed a very narrowly localised structure in the
vicinity of the ionosphere. This strong localisation is a result of phase mixing and of
the preferential growth of the mode in the vicinity of the resonant field lines. This
strong spatial localisation is tied to the strongly anisotropic wave propagation along
magnetic field lines. We stress that it was possible to obtain these features in our
Unstructured Meshes and Finite Elements in Space Plasma Modelling
137
model because our unstructured mesh is aligned along the field lines. Had we used a
mesh that was not aligned along the field lines (a mesh based on quadtree decomposition, for example), it would have been impossible to obtain such a level of spatial
localisation near the resonant surface.
Fig. 14. Perturbed density resulting from the ponderomotive force of the driven shear Alfvén wave, after
five periods.
The density variation in the magnetosphere can affect the localisation and structure of the resonance. The nonlinear coupling between shear and compressional
modes, together with the ponderomotive force in Eq. (14), cause the plasma to move
along the magnetic field lines, from high latitudes, down to the equator. This produces
density cavities near the ionosphere, and density bumps near the equator. Figure 14
shows the resulting density perturbation after driving the mode for five periods. This
density redistribution along a field line causes the wave eigenfrequency to decrease,
and it was used to explain features observed in discrete auroral arcs [Lu et al., 2003;
Rankin et al., 2004].
8
Example Applications in 3D
The extension of unstructured meshes and the finite element technique to three
dimensions is straightforward. In that case, elements are made of tetrahedra, and the
procedure for discretising partial differential equations is otherwise identical to what
has already been described in the previous sections. In what follows, we illustrate the
use of unstructured meshes with two example applications. The first one consists of
using unstructured meshes as a convenient means of parametrising a magnetic field
bundle, along which one dimensional eigenmode equations can be solved. In the
second example, we show results from a full three dimensional simulation of mode
coversion of compressional waves into shear Alfvén waves in a numerically generated
138
R. Marchand et al.
$1
8.4762
10.707
12.938
15.169
17.4
19.631
Fig. 15. Illustration of a flux tube surrounding a reference field line, for the purpose of calculating the
frequency of a shear Alfvén wave. The colour scale refers to the logarithm of the reciprocal of the
volume of the mesh elements.
(non dipolar) magnetospheric geometry.
8.1 Representation of a magnetic flux tube
Shear Alfvén waves are known for being a singular solution of the linearised
MHD equations when the time dependence of the solution is assumed to be of the
form exp(−iωt). In practice, this implies that, in the ideal MHD approximation,
non zero solutions of the wave equations can exist on single field lines or surfaces,
while vanishing everywhere else. This leads to the analysis of shear Alfvén waves
in the one dimensional limit, where the wave frequency is obtained from an eigenvalue problem along selected magnetic field lines. Several such analyses have been
described for this mode, in the past. Recently, Rankin et al. [Rankin et al., 2005]
formulated the problem in terms of metric coefficients along magnetic field lines represented in terms of Euler potential, and local derivatives along those lines. In this
approach, magnetic field lines need to be traced near a reference field line, along
which the required derivatives and metric coefficients are calculated. A straightforward calculation of these derivatives may prove to be numerically very intensive, and
time consuming, as it would require numerically tracing magnetic field lines many
times near the reference line. An alternative would be to construct a mesh enclosing
a small flux tube around the reference field line, and use it to calculate the necessary metric coefficients at the mesh points. Interpolations could then be conveniently
calculated at arbitrary locations along the field line.
An example of such a small flux tube is illustrated in Fig. 15. In this case, the
reference field line is traced from a starting point A in the northern hemisphere, extending to the southern hemisphere. A number of field lines are then traced similarly,
Unstructured Meshes and Finite Elements in Space Plasma Modelling
139
from points distributed on a small circle around A. Mesh points are then distributed
along these field lines, and a tetrahedral mesh is constructed. It is then straightforward
to calculate all the derivatives needed for the analysis of a one dimensional (singular)
shear Alfvén wave at each mesh point, and interpolate their values along the line.
8.2 Three dimensional simulation of mode conversion
Let us now illustrate the full power of our numerical approach by showing modelling results of the coupling between a compressional wave and a shear Alfvén wave,
in a numerically generated magnetospheric geometry. In this example, the global
MHD code BATS-R-US is used to calculate the state of the magnetosphere for an
inclination of Earth magnetic dipole by 35◦ and representative solar wind conditions.
These are the same conditions as considered in a recent article on the topology of
magnetic fields [Watanabe et al., 2005]. This solution is then used to trace a number
of magnetic field lines and, as explained in an earlier section, generate an unstructured mesh that is aligned along magnetic field lines. For simplicity, we solve for
a set of linearised MHD equations, and assume a cold plasma (i.e., pressure terms
are neglected). Specifically, we solve for the linearised form of Eqs. (10), (11), (12)
and (14). In the absence of pressure terms in the plasma momentum equation, and
in the linear approximation, it is not necessary to solve for Eqs. (13) and (15). The
boundary of the mesh used in this calculation is illustrated in Fig. 16. In this simuZ
Y
VU
X
Z
Y
X
Fig. 16. Illustration of the boundary of a three dimensional unstructured mesh constructed with a magnetospheric magnetic field computed with BATS-R-US.
lation, compressional waves are driven at the boundary of the domain, by effectively
“pushing” and “pulling” plasma in and out of the simulation domain at the dawn and
dusk sides, with a prescribed frequency, corresponding to a period of 200 s. The
140
R. Marchand et al.
induced compressional wave then propagates inward and mode converts into shear
waves. For simplicity, this simulation assumed a perfectly conducting ionosphere.
Results obtained after three oscillation periods, are shown in Fig. 17. It can be seen
Z
YX
VU
1
Z
3
2
3
Y
2
X
1
vy
-3.7533e+05
-2.1373e+05
-52130
1.0947e+05
2.7107e+05
4.3267e+0
Fig. 17. Illustration of the component of a velocity field associated the shear mode driven by mode
conversion of a compressional wave imposed at the simulation boundary. This result was obtained after
driving the compressional mode for three period.
that, even though the mode has not yet had time to phase mix into a strongly localised
mode, there is a clear signature of a shear wave in the evening and day sides of the
magnetosphere.
9
Summary and Conclusion
We have shown how simplex-based unstructured meshes can be instrumental in
representing simulation domains in space plasmas. Such meshes can be used, in
conjunction with a finite element discretisation of the appropriate fluid (MHD, multifluid, or other) equations to simulate the dynamics of the magnetosphere of space
plasmas in general. The approach described here is somewhat unconventional, as
it differs in several significant ways, from the most common numerical approaches
used in space plasmas. It does rely, however, on sound mathematical grounds and on
proven techniques used commonly in other fields of research. The main advantages
of using simplex-based unstructured meshes to represent space plasmas resides with
their local adaptivity, and the possibility of orienting these meshes along directions
of strong transport anisotropy. Another advantage is that they can naturally be used
to represent domains with complex boundaries. Their disadvantage comes from their
added complexity, and the fact that a Galerkin finite element discretisation of partial
differential equations on such grids leads to a sparse system of equations with matrices that do not have a simple structure. In general, the solution of these equations
Unstructured Meshes and Finite Elements in Space Plasma Modelling
141
therefore requires advanced iterative solution techniques that tend to be computationally more intensive than those resulting from a discretisation on structured grids.
Simplex-based unstructured grids and finite elements are therefore not presented here
as a replacement of existing numerical techniques for simulating space plasmas in
general. They are rather proposed as a valuable complementary tool for modelling
plasmas when other techniques would fail or would not be adequate. This is the case,
for example, when modelling plasmas with strong anisotropy in magnetic fields with
complex topologies.
Acknowledgments. This work was supported by the Natural Science and Engineering Research Council of Canada, the Canadian Space Agency, and the University of Alberta. Some
of our calculations were made using the computing infrastructure provided by WestGrid.
References
Al-Kurdi, A. and D. R. Kincaid, LU-decomposition with iterative refinement for solving sparse linear
systems, J. Comput. Appl. Math., 185, 391–403, 2006.
Bowyer, A., Computing Dirichlet tessallations, Comput. J., 24, 162–166, 1981.
Braess, D., Finite elements: Theory, Fast Solvers, and Applications in Solid Mechanics, Cambridge University Press, 2nd Ed., 2002.
Briggs, W. L., V. E. Henson, and S. F. McCormick, A Multigrid Tutorial, SIAM, 2000.
Cabannes, H., Theoretical Magnetofluiddynamics, Academic Press, New York, 1970.
Chew, G. F., M. L. Goldberger, and F. E. Low, The Boltzmann Equation and the One-Fluid Hydromagnetic
Equations in the Absence of Particle Collisions, Proc. Roy. Soc. (London), A236, 112–118, 1956.
Chi, S., Z. Jun, and K. W. Kai, Distributed block independent set algorithms and parallel multilevel ILU
preconditioners, Journal of Parallel and Distributed Computing 65, 331–46, 2005.
De Berg, M., M. van Kreveld, M. Overmars, and O. Schwartzkopf, Computational Geometry: Algorithms
and Applications, Springer, 2000.
Delaunay, B., Sur la sphère vide, Bull. Acad. Science USSR VII: Class. Sci. Mat. Nat., 793–800, 1934.
Dendy, E. D., A general-purpose finite-volume advection scheme for continuous and discontinuous fields
on unstructured grids, J. Comp. Phys., 180, 559–583, 2002.
DeZeeuw, D. W. and K. G. Powell, An Adaptively-refined Cartesian Mesh Solver for the Euler Equations.,
J. Comp. Phys., 104, 55–68, 1993.
Dirichlet, G. L., Über die Reduction der positiven quadratichen Formen mit drei unbestimmten ganzen
Zohlen, J. reine angewan. Math., 40, 209–227, 1850.
Erlebacher, G., Finite-Difference Operators on Unstructured Triangular Meshes, Springer-Verlag, Lecture
Notes in Physics, 238, 22–53, 1985.
Ferraro, N. M. and S. C. Jardin, Finite element implementation of Braginskii’s gyroviscous stress with
application to the gravitational instability, Phys. Plasmas, 13, 092101-1–12, 2006.
Freidberg, J. P., Ideal Magnetohydrodynamics, Plenum Press, 1987.
Frey, P. J. and P. L. George, Maillages: applications aux éléments finis, Hermes Sciences Publicat. 1999.
Frink, N. T. and S. S. Pirzadeh, Tetrahedral finite-volume solutions to the Navier-Stokes equations on
complex configurations, International Journal for Numerical Methods in Fluids, 31, 175–187, 1999.
Garcia, M. D. A. Suarez, L. Gonzalez, and G. Montero, New implementation of QMR-type algorithms,
Computers and Structures, 83, 2414–2422, 2005.
Grad, H., Principles of kinetic theory of gases, Handbuch der Physik, 12, Springer-Verlag OHG, Berlin,
1957.
Gruau, C. and T. Coupez, 3D tetrahedral, unstructured and anisotropic mesh generation with adaptation
to natural and multidomain metric, Computer Methods in Applied Mechanics and Engineering, 194,
4951–4976, 2005.
142
R. Marchand et al.
Gupta, A., Recent advances in direct methods for solving unsymmetric sparse systems of linear equations,
ACM Transactions on Mathematical Software, 28, 301–324, 2002.
Gupta, A., A shared- and distributed-memory parallel sparse direct solver, Lecture Notes in Computer
Science 3732, 778–787, 2006.
Huba, J. D., Hall magnetohydrodynamics in space and laboratory plasmas, Plasma Phys., 2, 2504–2513,
1995.
Jardin, S. C., Implicit solution of the four-field extended-magnetohydrodynamic equations using highorder high-continuity finite elements, Phys. Plasmas, 12, 56101-1–10, 2005.
Kuzmin, D, M. Moller, and S. Turek, Multidimensional FEM-FCT schemes for arbitrary time stepping,
International Journal for Numerical Methods in Fluids, 42, 265–295, 2003.
Li, N., Y. Saad, and E. Chow, Grout versions of ILU for general sparse matrices, SIAM Journal of Scientific
Computing, 25, 716–728, 2003.
Li, X. S. and J. W. Demmel, SuperLU DIST: A scalable distributed-memory sparse direct solver for unsymmetric linear systems, ACM Transactions on Mathematical Software, 29, 110–140, 2003.
Löhner, R., Renumbering strategies for unstructured-grid solvers operating on shared-memory, cachebased parallel machines, Computer Methods in Applied Mechanics and Engineering, 163, 95–109, 1998.
Lomax, R. J., Preservation of the conservation properties of the finite element method under local mesh
refinement, Computer Methods in Applied Mechanics and Engineering, 12, 309–314, 1977.
Lu, J. Y., R. Rankin, R. Marchand, V. T. Tikhonchuk, and J. Wanliss, Finite element modelling of nonlinear
dispersive field line resonances, J. Geophys. Res., 108, 1394, doi:10.1029/2003JA010035, 2003.
Lu, J. Y., R. Rankin, R. Marchand, and V. T. Tikhonchuk, Nonlinear electron heating by resonant shear
Alfvén waves in the ionosphere, Geophys. Res. Lett., 32, L01106, 10.1029/2004GL021830, 2005a.
Lu, J. Y., R. Rankin, R. Marchand, and V. T. Tikhonchuk, Reply to comment on “Nonlinear Electron
Heating by Resonant Shear Alfvén Waves in the Ionosphere” by J.-P. St.-Maurice, Geophys. Res. Lett.,
32, L13103, 10.1029/2005GL023149, 2005b.
Marchand, R. and M. Simard, Finite element modelling of TdeV edge and divertor with E × B drifts, Nucl.
Fusion, 37, 1629–1639, 1997.
Marchand, R., M. Charbonneau-Lefort, M. Dumberry, and B. Pronovost, ARANEA, a program for generating unstructured triangular meshes with a JAVA Graphics User Interface, Comp. Phys. Comm., 139,
172–185, 2001.
Michev, I. D., Finite volume and finite volume element methods for nonsymmetric problems, Ph.D. dissertation, Texas A&M University, 1996.
Powell, K. G., L. Roe, T. J. Linde, T. I. Gombosi, and K. L. DeZeeuw, J. A Solution-Adaptive Upwind
Scheme for Ideal Magnetohydrodynamics, J. Comp. Phys., 154, 284–309, 1999.
Rankin, R., J. Y. Lu, R. Marchand, and V. T. Tikhonchuk, Spatiotemporal characteristics of ultra-low
frequency dispersive scale shear Alfvén waves in Earth’s magnetosphere, Phys. Plasmas., 11, 1268–
1276, 2004.
Rankin, R., K. Kabin, and R. Marchand, Alfvénic field line resonances in arbitrary magnetic field topology,
Adv. Space Res., 38, 1720–1729, 2006.
Rebay, S., Efficient unstructured mesh generation by means of Delaunay triangulations and BowyerWatson algorithm, J. Comput. Phys., 106, 125–138, 1993.
Saad, Y., Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, 2003.
Scott, J. A., On ordering elements for a frontal solver, Commun. Numer. Meth. Engng, 15, 309–323, 1999.
Sloan, S. W., A fast algorithm for constructing Delaunay triangulations in the plane, Adv. Eng. Software,
9, 34–55, 1987.
Sloan, S. W., A FORTRAN program for profile and wavefront reduction, Int. J. Meth. Eng., 28, 2651–2679,
1989.
Šolín, P., K. Segeth, and I. Doležel, Higher Order Finite Element Methods, in Studies in Advanced Mathemathics, Chapman, Hall/CRC, 2003.
Voronoï, M. G., Nouvelles applications des paramètres continus à la théorie des formes quadratiques, J.
reine angewan. Math., 134, 198–287, 1908.
Unstructured Meshes and Finite Elements in Space Plasma Modelling
143
Watanabe, M., K. Kabin, G. J. Sofko, R. Rankin, T. I. Gombosi, A. J. Ridley, and C. R. Clauer, Internal
reconnection for northward interplanetary magnetic field, J. Geophys. Res., 110, A06210, doi:10.1029/
2004JA010832, 2005.
Watson, D. F., Computing the n-dimensional Delaunay tessellation with application to Voronoi polytopes,
Comput. J. 24, 167–172, 1981.
Whitaker, D. L., Two-dimensional Euler computations on a triangular mesh using an upwind, finite volume
scheme Ph.D. dissertation, Virginia Polytechnic Institute and State University, 1988.
Yang, Y. J., J. H. Yong, and J. G. Sun, An algorithm for tetrahedral mesh generation based on conforming
constrained Delaunay tetrahedralization Computers and Graphics (Pergamon) 29, 606–615, 2005.