Project topics and materials
Get Complete Project Material File(s) Now! »
Mesh generation and descriptions
This section describes the generation of the different meshes used in this study. In what follows, the following definitions are needed:
Definition 2.1.1 (Spherical polygon). It is a part of the spherical domain bounded by spherical arcs. A spherical arc, in turn, is the spherical equivalent of a straight segment, i.e. the shortest-length curve joining two points on the sphere.
Definition 2.1.2 (Mesh cell). Each polygon within a partition of the spherical domain into spherical polygons, not containing smaller polygons of the same mesh, is referred to as a cell of the mesh.
Definition 2.1.3 (Neighbouring cells). Cells adjacent to each other by sharing a spherical arc are called neighbouring cells. The common spherical arc is an edge of both cells.
Quasi-uniform Voronoi-mesh generation
The following paragraph describes the construction of the quasi-uniform hexagonal-pentagonal Voronoi mesh used in this study. This mesh is inspired, in particular, by the mesh of Sadourny et al. [1968], followed by the optimization procedure of Du et al. [1999].
Definition 2.1.4 (Voronoi mesh). A Voronoi mesh is one generated from any distribution of points, each of which is the generator of a cell. The latter is composed of all the points of the domain which are closer to the considered generator than to any other generator.
One result of this definition is that the cells of a Voronoi mesh are polygons. In practice, we start with a regular icosahedron, which is a polyhedron with 20 flat triangular faces. The size of this icosahedron does not matter, as long as the triangular faces are equilateral, and equal to one another. The icosahedron is then projected, radially, onto a sphere. The result of this step is a sphere whose surface is divided into 20 equal spherical triangles, with curved edges and surfaces (Fig. 2.1). The resulting surface is referred to as an icosahedral mesh.
The following step consists of refining this icosahedral mesh. To do so, we follow the method of Sadourny et al. [1968]. Each spherical triangle of the spherical icosahedral mesh is considered. Two of its edges are divided into n arcs of equal length. The end of each arc is numbered from 1 to n increasingly from the point common to the two edges, and until the end of the considered edge. Each end of arc is joined with the corresponding end of arc on the second edge of the considered triangle. The lines joining these points together are in turn arcs. Each of the latter is divided into sub-arcs in turn, with an increasing number as we move away from the common triangle vertex. The first arc is kept as 1 segment, the second is divided into 2, and so forth until we divide the nth arc into n sub-arcs. This final arc is, in fact, the third edge of the initial triangle. The ends of these segments are joined together to result in n2 sub-triangles within the initial triangle. This step’s procedure falls into the category of triangular tessellation [Du et al., 1999]. The edges of the resultant spherical sub-triangles, which are no longer equal, have edges which vary to within 18% of each other’s lengths [Wang and Lee, 2011]. This resultant mesh is a Delaunay triangulated mesh.
The refined triangular mesh of the above paragraph is referred to, in this study, as the dual mesh. It is only used in some sub-calculations detailed in the sections below. The cells of this mesh are referred to as dual cells. The next step in mesh generation is forming what we call our primal mesh, where the main computations of the numerical schemes are made.
To describe, completely, the primal mesh cells, we first obtain their vertices—i.e. the points where the cell edges intersect. By design of the Voronoi mesh, each vertex is the circumcenter of a triangle of the dual mesh. Another con-sequence of the definition of a Voronoi mesh is that the Voronoi cell edges and dual edges are orthogonal [Augenbaum and Peskin, 1985]. The shapes of the cells of this mesh are primarily hexagons, with 12 pentagons, one surrounding each vertex of the initial 12-vertex icosahedron. Each primal cell also has a centroid, which, in general, does not coin-cide with the generator. Following Lloyd’s iterative algorithm as in Du et al. [1999], the generators of the Voronoi mesh are made to coincide with the centroids of the Voronoi cells by an iterative process. This procedure produces a more uniform mesh.
Variable-resolution Voronoi-mesh generation
Similarly to the quasi-uniform Voronoi meshes, a variable resolution Voronoi mesh is built upon a Delaunay triangulation. However, to account for the variable resolution, the starting point cannot, in this case, be a regular icosahedron. Instead, a set of points is randomly distributed on the surface of a unit sphere, with a higher probability in regions where higher resolution is sought. A Delaunay triangulation is formed using these points, and then the Voronoi dual is taken in the same way as in the previous case. Also, the target resolution is taken into account during the Lloyd iteration process [Du et al., 1999]. It is important to note that, with variable resolution, the Voronoi mesh is no longer hexagonal-pentagonal: at least some cells have seven edges or more.
Voronoi and Delaunay mesh elements
The following section describes the different parts of the primal and dual meshes, which are defined and annotated for use in the numerical schemes to come. The definitions below apply to the primal and dual meshes simultaneously. The indices used in this section are broadly inspired by the notations of Thuburn et al. [2009].
Mesh cells
Each cell is referred to using a fixed cell index, k. These indices are unstructured: there is no particular direction for the index increment (Fig. 2.2). A neighbouring cell may or may not have an index with a unit difference from the current cell index. The total number of indices is equal to the total number of cells, with no two cells sharing the same index. Dual mesh cells are referred to by an index, v.
Cell areas
The cell of index k has an area Ak. A dual cell with index v has an area Av.
Primal-cell centers
Each primal cell has a geometric center, or centroid, referred to by its position, xk. The number of cell centers is—quite evidently—equal to the number of cells, and so these centers are assigned the same index, k, as their respective cells. The position of each cell center, in spherical coordinates, is found using the arithmetic mean position of all the points belonging to the cell: xk = RVk x d A (2.1) where Vk represents the primal cell itself, regarded as a subset of the sphere.
Cell edges
In the Voronoi mesh, each edge of the mesh is referred to using an unrepeated index, e. Each primal cell, k, has a set of edges, Ek. A dual cell , v, has the set of edges Ev.
Primal-cell vertices
Similarly to the edges, each cell of the Voronoi mesh has a number of vertices, Nk, equal to its number of edges. The indexing is analogous to that of the edges, with a unique index v for each vertex, ordered arbitrarily.
Primal-edge lengths
The values of the primal edge lengths are designated by le, where e refers to the considered edge. le has a value corresponding to the arc length of the edge. Similarly, dual edges have lengths de.
Edge normal
Each edge has a normal vector, ne, where e refers to the considered edge. As each edge is common to two cells, this vector points in an arbitrary direction in or out of the considered cell k. A factor nk;e = 1 is defined in such a way that vector nk;ene points outwards from cell k. This vector is used in the gradient computation detailed in Section 2.2.4.
Discrete representation of fields
Discrete degrees of freedom
The flux-form advection equation (1.2) is an initial-value problem (IVP) for the density field = . This field, (x), defines infinitely many real numbers, one for each value of the position x. The discretization of fields is the process of reducing this information to a finite set of numbers, also called discrete degrees of freedom (DOFs), which can then be manipulated by computers. Of course, this procedure incurs a loss of information, and subsequent errors called discretization errors. One important property of discretization schemes is their order of accuracy, which describes how fast the discretization error decreases as the number of discrete DOFs increases. If spatial discretization errors are O( x)p, then the scheme is of order p. For a given problem, the error of higher-order schemes decreases faster with increasing resolution than that of lower-order schemes.
Reducing a field to a finite set of real numbers can be done following multiple approaches. Taking a finite-difference approach, values of (x) are chosen at meaningful points, (xk). Following a finite-volume approach, on the other
36 CHAPTER 2. DISCRETIZATION
hand, the field (x) is represented through its averages over a finite set of control volumes. In this study, we adopt a finite-volume approach, where the density of a species is represented by the cell-averaged value of its continuous distribution. Thus, a discrete value of the field is attached to each cell and referred to with the same index: k = RVk Ak (2.2)
Similarly, the discretized density, k, is defined as a cell average: k = RVkAkd : (2.3)
Initial value specified from a closed-form expression
In the practical sense, if we have an analytical expression for the density, ( ; ; t), then its value is known at any point of the sphere. Using quadrature formulae, its average over primal cells can be computed almost exactly. However, a simple evaluation at one point, the cell centroid, yields a relatively accurate estimate of the true average. Indeed, this can be seen by considering the integral: I =(x)d2x: (2.4)
Adding and subtracting xk, the position of the centroid Ck of the cell, the expression becomes: I = ZVk (xk + (x xk))d2x
Taking a second-order Taylor expansion about the centroid, xk, we obtain: I = Vk » (xk) + (x xk) @ + O kx xkk2 # d2x
which is equivalent to: I (xk) ZVk d2x = @x xk ZVk (x xk) d2x + ZVk O kx xkk2 d2x
Now, the definition of the cell barycenter is precisely such that the first term on the right-hand-side vanishes. Dividing by the area of a cell, Ak = RVk d2x, we obtain: k(xk) = O x2 : (2.8)
From Eq. 2.8, it can be seen that the cell-averaged value can be approximated to a point-wise value at the cell center, with second-order accuracy. This property is used to compute initial values for k given a closed-form expression for (x). Conversely, during the course of the computation, whenever point-wise values of (x) are sought, they can be estimated—cheaply—at xk from the known cell-averaged values. The same can be done for the density, . As a result, the mixing ratio, , can be estimated at xk with second-order accuracy as the ratio: (xk) = k + O x2 where k = k= k: (2.9)
These cell-center values are used for the reconstruction of edge values, needed for the advection schemes (detailed in Section 2.3.2).
Edge values of the mixing ratio
The mixing ratio at the midpoints of the cell edges is denoted by e. This edge value is needed in the numerical schemes detailed in Section 2.3. Its evaluation depends on the chosen numerical scheme, starting from the value at the respective cell center. This step of the numerical formulation of schemes is referred to as reconstruction [Dubey et al., 2015]. To transport the tracer field from one cell to one of its neighbouring cells, we reconstruct the tracer mixing ratio at the edge common to the two cells.
Discretized velocity field
The velocity field, u, on the other hand, is a vector field, and a different approach is used to discretize it. To each edge, of index e, we attach a real number, ue, which we define as the velocity vector—at the edge midpoint—projected onto the normal unit vector of the edge.
Discretized gradient
This quantity is useful for the second order numerical schemes requiring a discrete gradient (Section 2.3.3). This is the discrete counterpart of the continuous gradient of the continuous mixing ratio field. In our study, we go by the work of Dubey et al. [2015], who demonstrate that a first order approximation of the gradient is sufficient to maintain second order scheme accuracy. Following their first order approximation, a first-order estimate of the gradient, rv is obtained at a each dual cell (with index v) as the dual-cell-averaged gradient: rv = r dA (2.10) where Av is the area of the dual cell v, and Vv is the dual cell itself (i.e. a triangular subset of the sphere). The vertices of these dual cells correspond to the centers of the the primal cells, where the mixing ratio can be estimated with second-order accuracy as discussed above.
Using the Green-Gauss theorem, the above area integral is converted into a line integral: ZVv r dA = Z nd (2.11) where is the perimeter of the vth cell, and n is the normal relative to the infinitesimal edge, d , directed outwards.
Discretizing the above expression using second-order-accurate estimates of the mixing ratio at primal-cell centers, we obtain the discrete gradient of the continuous scalar field at the dual cell v: rv = Av e2Ev ( v;e v) mv;env;ede where nv;e is the vector normal to the dual edge e, and mv;e = 1 is the factor ensuring the vector is directed out of dual cell v (analogous to the nk;e of primal cells). de is the length of this dual edge, v;e is the half sum of estimated at the extremities of the dual edge (which are cell centers) and v is the equal-weight average of at dual-cell vertices, which are also cell centers. We apply the Green-Gauss formula (2.11) to v, which has the same gradient as . This ensures that rv is zero whenever the mixing ratio is uniform. In a second step, rk is computed for each cell k by performing a dual-cell-area-weighted average of rv for all dual cells v intersecting the primal cell k.
Numerical solution of the advection equation
Time discretization
The solution of the advection equation cannot be computed, in practice, at infinitely many values of time t. Instead, we estimate the solution at integer multiples of a time step, t, such that: t = T (2.13) where T is the experiment length, and n is the total number of time steps. In practice, as the problem is an IVP, given (kn), the discrete representation of (x) at t = n t, one then estimates (kn+1) corresponding to time t+ t = (n+1) t. (kn) is then set aside to proceeds further in time. This repeated procedure is the time iteration, or time loop.
The numerical schemes marching the approximate solution of the advection equation from time t to time t + t are classified depending on their separate or coupled treatment of space and time. In the first category is the method of lines, consisting of regarding time as continuous at first, while discretizing the right-hand side of the evolution equation of . This step involves only spatial derivatives and is called the spatial discretization. It transforms the original PDE into a set of many coupled ordinary differential equations (ODEs), one for each cell index k. In the second step of the method of lines, time integration schemes for ODEs are used to solve numerically the previously obtained set of ODEs.
In the second category, the difference (kn+1) (kn) is deduced from the advection equation and the finite-volume definition of the discrete degrees of freedom k, and approximated. Such schemes are called space-time coupled. Space-time coupled schemes make explicit use of the knowledge of the problem being solved, while this knowledge is not used during the time integration step of the method of lines.
In what follows, the space and time decoupled schemes are detailed, followed by an explanation of the used coupled scheme, in the Voronoi mesh framework. The Cartesian schemes are sketched in Section 2.5.3.
Method of lines: time-integration schemes
In this section, we describe the time-integration schemes used as part of the method of lines. Hence, it is assumed that the spatial discretization, described in the next subsection, is known: there is a well-defined procedure to compute, for each cell index k, @ k=@t, given the values of k for all k, and other known quantities, such as the discretized wind and air density fields.
Below, we describe the first- and second-order Euler and Heun schemes, respectively. For the simplicity of notations, we omit the subscript k. Thus, we consider the model problem with a single degree of freedom: @t = f( (t)) (2.14)
Euler method
The simplest time integration scheme, in terms of computational complexity, is the first-order Euler scheme. To begin solving the initial-value problem following the Euler method, a first-order finite-difference approximation is invoked: (t + t) (t) (2.15) which yields an explicit expression: (t + t) = (t) + tf( (t)) + O(( t)2): (2.16)
From the above expression, an iterative discretized form can be stated: (n+1) = (n) + tf( (n)) (2.17) where O(( t)2) is the error term of a first-order finite-difference approximation with step size t.
Heun method
This time scheme is a second-order evolution from the Euler scheme [Hairer, 1996]. It takes the prediction step of first-order accuracy, and adds a correction step which brings the scheme to second order.
Again, the starting point is the model ODE (2.14), from which an advance-in-time algorithm is developed. An Euler ~(n+1) of at time t+ t. The time derivative of at t+ t is then calculated step is first applied, leading to an estimate ~ (n+1) and used to estimate (n+1) (n) by the trapezoidal rule, with second order accuracy: from ~(n+1) = (n) + f( (t)) t; (n+1) = (n) + 1 t hf( (n)) + f(~(n+1)) i :
Method of lines: spatial schemes
Having dealt with the time-integration problem, we are left with the purely spatial problem of estimating the right-hand-side of the flux-form advection equation (1.2), or, more precisely, its discrete representation.
Upwind scheme
The upwind, or Godunov, scheme [Godunov, 1959] is a 1st-order spatial scheme based on the direction of wind flow. The cell-averaged values of tracer density are the main input of the scheme. These are obtained from an initial condition or from the previous time step. In addition to k, the scheme takes as input the cell-averaged air density, , and the discretized wind field, or, more precisely, the discretized air-flux: Ue = u ne: (2.20)
From these values, as explained before, one can immediatelty estimate the mixing ratio at cell barycenters as k = k= k. Then, the edge values of are approximated according to the following conditional equation following the notations in Fig. 2.3: e = 8 i; if Ue < 0: (2.21) < : j; if Ue > 0: where i and j are the concentration ratios at the adjacent cells indexed by i and j respectively, and e is the value at the common edge. Between these two cells, the upwind one is the cell from which the wind field arrives at the edge under consideration.
Finally, to complete the method of lines, an advance in time is made using one of the aforementioned time integration schemes. For consistency in scheme order, the upwind spatial scheme is usually paired with the Euler method (Section 2.3.2) : (kn+1) = (kn) 1 nk;efe (2.22) where Ek is the set of edges of cell k, and the tracer flux is: fe = Uele e(n) t: (2.23)
Green-Gauss gradient scheme
The Green-Gauss gradient scheme is an upwind-biased 2nd-order spatial scheme, based on the upwind scheme, but with a linear reconstruction of the tracer mixing ratio—instead of a piecewise-constant reconstruction—to improve the convergence order. The exact scheme used is the one detailed by Dubey et al. [2015].
The Green-Gauss gradient is first calculated in each primal cell. Then, edge values are computed using the linear Taylor expansion : e = k + rk (xe xk) (2.24) where k is the index of the cell upwind of edge e, and rk is approximated as detailed in Section 2.2.4.
Following the method of lines, this scheme is usually paired with the 2nd-order Heun time scheme to maintain its order of convergence (Section 2.3.2).
Coupled time and space scheme
We also consider a space-time coupled scheme, i.e. one whose iterations cover both time and space dimensions in a coupled manner. The discretized equations of such a scheme account for time and space variations together.
The Semi-Lagrangian Finite-Volume (SLFV) scheme is a second-order space-time coupled scheme for the integration of the advection equation formulated for a spherical unstructured mesh by Dubey et al. [2015]. Following this finite volume method, where the cell-averaged is predicted, the flux form of the advection equation (1.2) is taken as a starting point for the discretization.
Table of contents :
1 Introduction
1.1 Context
1.1.1 Sharp plumes
1.1.2 The advection equation
1.2 Methods for the numerical resolution of the advection equation
1.2.1 Cartesian grid formulations on the sphere
1.2.2 Triangular and polygonal meshes
1.2.3 Numerical schemes
1.2.4 Applications
1.3 Research questions
1.4 General outline
2 Discretization
2.1 Mesh generation and descriptions
2.1.1 Quasi-uniform Voronoi-mesh generation
2.1.2 Variable-resolution Voronoi-mesh generation
2.1.3 Voronoi and Delaunay mesh elements
2.2 Discrete representation of fields
2.2.1 Discrete degrees of freedom
2.2.2 Edge values of the mixing ratio
2.2.3 Discretized velocity field
2.2.4 Discretized gradient
2.3 Numerical solution of the advection equation
2.3.1 Time discretization
2.3.2 Method of lines: time-integration schemes
2.3.3 Method of lines: spatial schemes
2.3.4 Coupled time and space scheme
2.4 Slope limiting
2.5 Cartesian meshes
2.5.1 Cartesian-mesh generation and elements
2.5.2 Operator splitting
2.5.3 One-dimensional transport schemes
3 Two-Dimensional Test Cases
3.1 Tracer-concentration distributions
3.1.1 Uniform distribution
3.1.2 Double cosine-bell distribution
3.1.3 Single cosine-bell distribution
3.2 Wind fields
3.2.1 Solid-body rotation
3.2.2 Tilted solid-body-type periodic oscillation
3.2.3 Dual vortex with solid-body rotation (NL2010)
3.3 Limiting to the Cartesian band
3.4 Metrics
3.4.1 Stability
3.4.2 Monotonicity
3.4.3 Convergence
3.4.4 Numerical diffusion
4 Findings on the Numerical Mesh-Scheme Pairs in 2D
4.1 Stability
4.2 Monotonicity
4.3 Numerical diffusion
4.3.1 Entropy
4.3.2 Preservation of non-linear relationships between tracers
4.4 Convergence
4.4.1 Zonal solid-body-rotation convergence
4.4.2 Tilted periodic solid-body convergence
4.4.3 NL2010 low-shear convergence
4.4.4 NL2010 intermediate-shear convergence
4.4.5 NL2010 default-shear convergence
4.5 Partial conclusions
5 Puyehue-Cord´on Caulle 3D Simulation Results
5.1 Puyehue-Cord´on Caulle eruption
5.2 3D simulations
5.3 Trajectory of the plume
5.3.1 SO2 concentration at 330 hPa
5.3.2 Column-integrated SO2 concentration
5.3.3 Plume trajectory
5.4 Summary
6 Conclusion
6.1 Main results
6.2 Perspectives for chemistry-transport modeling
7 Sommaire en franc¸ais
7.1 R´ esultats principaux
7.2 Perspectives pour la mod´ elisation chimie-transport
