Currently, most operational forecasting models use latitude–longitude grids, whose convergence of meridians toward the poles limits parallel scaling. Quasi-uniform grids might avoid this limitation. Thuburn et al. and Ringler et al. have developed a method for arbitrarily structured, orthogonal C grids called TRiSK, which has many of the desirable properties of the C grid on latitude–longitude grids but which works on a variety of quasi-uniform grids. Here, five quasi-uniform, orthogonal grids of the sphere are investigated using TRiSK to solve the shallow-water equations.
Some of the advantages and disadvantages of the hexagonal and triangular icosahedra, a “Voronoi-ized” cubed sphere, a Voronoi-ized skipped latitude–longitude grid, and a grid of kites in comparison to a full latitude–longitude grid are demonstrated. It is shown that the hexagonal icosahedron gives the most accurate results (for least computational cost). All of the grids suffer from spurious computational modes; this is especially true of the kite grid, despite it having exactly twice as many velocity degrees of freedom as height degrees of freedom. However, the computational modes are easiest to control on the hexagonal icosahedron since they consist of vorticity oscillations on the dual grid that can be controlled using a diffusive advection scheme for potential vorticity.
For modeling the global atmosphere, grids with inherent resolution clustering, such as the latitude–longitude grid (lat-lon) and (less severely) the conformal cubed sphere (Rančić et al. 1996; Putman and Lin 2007), are expected to be inefficient on massively parallel computing architectures for which data communication will be a performance bottleneck (e.g., Swinbank and Purser 2006) since resolution clustering leads to time-step restrictions that are circumvented by algorithms that use a larger domain of dependence per time step. In recent years this has reinvigorated research into the use of quasi-uniform grids. However, all proposed quasi-uniform grids have potential problems associated with their grid structure; most notably these include grid imprinting and the existence of computational modes (see below). This paper examines five candidate quasi-uniform spherical grids (section 2) for solving the shallow-water equations (SWEs) and compares them with the lat-lon grid. A common numerical framework is used throughout (section 3) so that the comparison is as clean as possible. The test cases and diagnostics presented (section 4) focus in particular on grid imprinting and computational modes.
No grid of the sphere is perfectly uniform. All grids have some nonuniformity (such as variations in cell size, edge length, cell orientation, or cell shape) and all grids have special points or lines where the structure is locally different from elsewhere (such as the cube vertices and edges on a cubed sphere grid). Numerical truncation errors can be expected to reflect the underlying grid structure, leading to grid imprinting in the numerical solution. The worst cases of grid imprinting can be quite conspicuous (e.g., Tomita et al. 2001). However, the manifestation of grid imprinting may be rather subtle, for example in the form of spurious triggering of physical barotropic or baroclinic instability (e.g., St-Cyr et al. 2008) or a spurious distortion of the energy spectrum. Some grids are heterogeneous at the grid scale; the basic repeating unit is a cluster of cells and a cell’s immediate neighbor may be quite different from the cell itself (e.g., on a triangular grid an “up” triangle is adjacent to a “down” triangle). Such grids may have grid-scale variation in numerical truncation error leading to noisy solutions (e.g., Le Roux et al. 2005; Danilov 2010).
Staniforth and Thuburn (2012) describe computational modes as wave modes supported by the discretization that have no analog among those supported by the continuous equations. They are usually at or near the grid scale. A discretization may support isolated individual computational modes (which are often stationary) or one or more entire branches. Some well-known examples of isolated stationary computational modes are the checkerboard gravity modes on the Arakawa A and B grids, the vertical grid-scale mode supported by the Lorenz-staggered vertical grid, and the Coriolis mode on the Arakawa C grid. Stationary computational modes are associated with a nontrivial null space of one or more discrete operators in the linearized governing equations [this is how Le Roux et al. (2005) define a computational mode]. A serious consequence of the A- and B-grid computational modes is that waves with spatial scale close to the grid scale have group velocity of the wrong sign and so energy propagates in the wrong direction (e.g., Durran 1999).
The C grid has more accurate gravity wave dispersion than the A and B grids but if the Rossby radius is not well resolved, inertia–gravity modes wrongly have decreasing frequency as wavenumber increases (Randall 1994). However, the Rossby radius is well resolved for the deepest modes in the atmosphere and so the C grid on a regular latitude–longitude grid is a common choice among operational forecasting centers (e.g., Davies et al. 2005) and a C grid is used for the results presented in this paper. There is also a stationary, nondivergent Coriolis mode on the C grid. This arises because only an averaged version of the discrete vorticity is seen by the momentum equation.
Computational modes are usually excited if they exist in models of the atmosphere since nonlinear dynamics and forcing from physical parameterizations can project onto them. If they grow then they can overwhelm the well-resolved part of the solution and so must be removed by some form of diffusion. It is therefore desirable to use a numerical method with fewest possible or no computational modes.
Le Roux et al. (2005) and Staniforth and Thuburn (2012) argue that a necessary (but not sufficient) condition for avoiding such spurious branches in 2D is that there should be twice as many velocity degrees of freedom (DOFs) as height. Otherwise, in addition to isolated stationary computational modes, entire branches of computational modes can exist in the discrete dispersion relation. For example, the C-grid discretization on Delaunay triangulations described by Bonaventura and Ringler (2005) was subsequently found to support four rather than two inertio-gravity modes for each horizontal wavenumber (e.g., Danilov 2010). Each triangle has 1 height DOF to just 1½ velocity DOFs. Hence smaller length scale waves in height are represented that cannot interact correctly with velocity because velocity is at a coarser resolution. When the Rossby radius is poorly resolved these computational modes are retarded; the smallest-scale modes have very small frequency and are then easily forced, which results in grid-scale oscillations of height and divergence (e.g., Gassmann 2011) that are coupled with the physical waves by the nonlinear terms.
The hexagonal C grid supports a branch of low-frequency computational Rossby modes in addition to the physical branch (Thuburn 2008). It is associated with the fact that each grid cell has 3 velocity DOFs per height DOF.
Most finite element pairs on triangles also have computational modes, which again may be isolated modes or entire branches (Le Roux et al. 2005). Some are due to collocation of pressure and velocity while those with more than twice as many velocity DOFs as pressure DOFs may have spurious branches of inertial oscillations (e.g., P1DG−P2; Cotter and Ham 2011).
A further issue with the hexagonal (and other polygonal) C grid is the maintenance of geostrophic balance. The hexagonal C grid was analyzed by Purser (1998), Ničković et al. (2002), and Torsvik et al. (2005), who found that it could not support geostrophic modes of zero frequency. Subsequently, Thuburn (2008), Thuburn et al. (2009), and Ringler et al. (2010) showed how the Coriolis term could be discretized to ensure stationary geostrophic modes. (In the rest of this paper their scheme is referred to as TRiSK.) Prior to TRiSK, models using hexagonal icosahedra have either used collocation of pressure and velocity (e.g., Tomita and Satoh 2004; Weller and Weller 2008); collocation of pressure, vorticity, and divergence (a hexagonal Z grid; e.g., Heikes and Randall 1995a); or collocation of both components of velocity (the ZM grid; Ringler and Randall 2002).
TRiSK is a horizontal discretization scheme for the SWEs on arbitrarily structured orthogonal C grids. (“Orthogonal” is defined to mean that the dual grid edges are at right angles to the primal grid edges. The dual grid has vertices at the cell centers of the primal grid and cell centers at the vertices of the primal grid.) TRiSK can maintain discrete geostrophic balance and has a number of other conservation and mimetic properties that will be described in section 3. Thus, these properties, which can be obtained relatively straightforwardly on the lat-lon C grid, can now also be obtained on a range of alternative grid structures. Note, however, that TRiSK does not eliminate the branch of computational Rossby modes on the hexagonal C grid or the branches of computational inertio-gravity modes on the triangular C grid, since it does not alter the numbers of DOF in the velocity and height fields. TRiSK does ensure that the frequency of both physical and computational Rossby modes goes to zero for constant f. Also, TRiSK provides a discrete analog of the potential vorticity (PV) equation with some flexibility in the choice of PV flux; an upwind-biased interpolation of the PV partially controls the computational Rossby modes (Weller 2012; see also section 4 herein).
2. Quasi-uniform orthogonal grids of the sphere
The construction of five quasi-uniform grids of the sphere is described along with some statistics of the resolutions used in the results in section 4. Shallow-water tests will, where possible, be done at very coarse resolution so that the relationship between the errors and the grid can be seen clearly at a global scale and so that fundamental length scales such as the Rossby radius can easily be underresolved, as sometimes happens in 3D models.
Some statistics of the grids, time steps used, and Rossby radii of the test cases are shown in Table 1. For each test case, grids of each type are generated with similar number of DOFs.
a. Skipped latitude–longitude
The skipped lat-lon grid (Purser 1998; see our Figs. 1c,d) has 2:1 coarsening in the longitude direction once the distance between cell centers reaches of the equatorial distance. To produce an orthogonal grid, the edges between fine and coarse cells are moved by first Delaunay triangulating the cell centroids and then taking the Voronoi dual to create a Voronoi grid. The result consists mostly of squares, some pentagons (where the longitudinal resolution reduces toward the pole), and a polygon at each pole.
It is hoped that the skipped lat-lon grid will maintain the accuracy of the lat-lon grid equatorward of the reductions while achieving the properties of TRiSK (see section 1) at the reductions and at the poles.
b. Hexagonal and triangular icosahedra
Hexagonal and triangular icosahedra (Figs. 1e,f) are close to uniform and have the minimum “angular deficiency” (Purser 1998). (The angular deficiency at the vertex of a polyhedron is the difference between the sum of the face angles at that vertex and 2π.) Hence distortions near the edges of the base polyhedron will be minimal. Icosahedral grids can also be constructed to be orthogonal using Delaunay triangulation and Voronoi tessellation.
Results will be presented on a centroidal Voronoi hexagonal icosahedron (e.g., Ringler et al. 2008; see our Fig. 1e) for direct comparison with Ringler et al. (2010); and on the dual grid, the Delaunay triangulation. A centroidal grid has Voronoi generating points collocated with cell centroids, which improves accuracy of the estimate of the cell average value by the cell-center value. But the Voronoi–Delaunay grid crossover points do not converge to the midpoints of the edges with increasing resolution, which reduces the accuracy of the estimate of the edge average value by the edge crossover point value. Alternatively, the Heikes and Randall (1995b) hexagonal icosahedron does not have Voronoi points at the cell centroids but the crossover points do converge to the midpoints of the edges. This version has also been tested; for the results presented here, this makes a little difference (not shown). Ringler (2003) found a similar insensitivity to the grid details.
c. Voronoi-ized cubed sphere
To make a quasi-uniform cubed sphere grid orthogonal and thus suitable for use with TRiSK, it has been “Voronoi-ized” by taking the cell centroids of the equal-angle cubed sphere, Delaunay triangulating them and then taking the Voronoi dual (Fig. 1h). This Voronoi-ized cubed sphere has lost some of the desirable properties of a cubed sphere as the cells are no longer quadrilateral (two short edges have been introduced into each cell) and so there are more than twice as many velocity DOFs as height. However, the cell to cell addressing is still very straightforward, which may allow better code optimization.
The Delaunay triangulation is not always unique on square grids of points; whether a quadrilateral is split into an up-triangle and then a down-triangle or vice versa is subject to rounding error. Therefore the Voronoi-ized cubed sphere will lose the symmetry of the original cubed sphere. (This loss of symmetry is visible in Fig. 7h.)
Grids of quadrilaterals will always have exactly twice the number of velocity DOFs as height DOFs using C-grid staggering, which is a necessary condition for avoiding computational modes. A grid of quadrilaterals can be created by dividing the triangles in any Delaunay triangulation into three kite-shaped quadrilaterals (Figs. 2 and 1g). The decomposition is achieved by lines from the triangle edge midpoints to the triangle circumcenters. The result has an orthogonal dual if the midpoint between the triangle circumcenter and the triangle corner is used as the kite center since the dark-shaded triangle is similar to the light-shaded right-angle triangle (which is a right angle since the Delaunay triangulation is orthogonal). Similar grids were used by Giraldo and Rosmond (2004).
However, the kites lack a number of other desirable grid properties: they are not Voronoi since the kites edges are not midway between the kite centers; they are not centroidal since the center used to make them orthogonal is not the center of mass, so the value at the center is not a second-order accurate approximation of the area average; they are heterogeneous at the grid scale (adjacent kites have different orientations) and they are anisotropic (there are widely varying edge lengths within each cell and widely varying cell center to cell center distances). We will show in this paper that these disadvantages outweigh the advantages of having the correct ratio of DOFs.
e. Cubed sphere
Various versions of cubed sphere grids have been proposed for atmospheric modeling, and have some advantages. For example, they retain a logically rectangular quadrilateral structure, so the full C-grid accuracy of the inertio-gravity wave dispersion may be possible without computational modes. Accurate tracer transport algorithms have been developed for the cubed sphere (e.g., Lauritzen et al. 2010). The SWEs have been solved accurately by Fournier et al. (2004) on the cubed sphere by employing high-order accurate collocated spectral elements, and this model has been extended to three dimensions and coupled with physics parameterizations to produce realistic rainfall simulations (Mishra et al. 2011). Given the popularity of the cubed sphere grids, it would be desirable to include them in the present comparison. Unfortunately, TRiSK, as currently formulated, requires an orthogonal grid. (An extension to nonorthogonal grids is an active area of research.) The conformal cubed sphere (Rančić et al. 1996) is orthogonal but suffers from resolution clustering around the cube corners and so is not quasi-uniform. Various quasi-uniform versions of the cubed sphere have been developed, including the equal-angle cubed sphere (e.g., Rančić et al. 1996), the equal distance cubed sphere, a version modified by spring dynamics, and a version modified based on a variational principle (Putman and Lin 2007). However, these are all nonorthogonal and so unsuitable for TRiSK. None of the proposed schemes for nonorthogonal cubed sphere grids has all of the desirable properties of TRiSK (section 3). We therefore restrict attention in this paper to grids that are suitable for TRiSK in order to maintain a clean comparison.
TRiSK is exactly defined by Thuburn et al. (2009) and Ringler et al. (2010) and here we summarize TRiSK (section 3a) and define some minor alterations relevant for non-Voronoi, orthogonal 2D grids on the sphere (section 3b). (Voronoi implies that the edges between cells are exactly midway between the cell centers.) TRiSK has the following desirable properties:
Exact, local conservation of mass.
Exact, local conservation of potential vorticity. This relies on, among other things, obeying a discrete analog of ∇ × ∇Φ ≡ 0 for any scalar field, Φ.
The evolution of PV is consistent in the sense that constant initial PV remains constant for all time.
The discrete PV is compatible with the momentum equation (i.e., the same result is obtained whether the PV is diagnosed at a particular time step from the prognostic variables of velocity and height or whether it is evolved using its own transport equation and the fluxes from the velocity and height).
Conservation of total energy to within time-truncation error (i.e., neither the Coriolis force nor the pressure gradient force generate energy and the sources and sinks of kinetic and potential energy are exactly equal and opposite).
Ability to maintain discrete geostrophic balance for the linearized equations in the limit of constant f.
a. Concise statement of TRiSK (Thuburn et al. 2009; Ringler et al. 2010)
TRiSK is a mixed finite-volume–finite difference solution method of the nonlinear, rotating SWEs with the continuity equation in flux form and the momentum equation in vector-invariant form:
where h is the fluid depth, u is the horizontal wind, b is the height of the bottom boundary, g is the acceleration due to gravity, is the potential vorticity, u⊥ = k × u, f is the Coriolis parameter, k is the local outward pointing unit normal to the sphere, and is the kinetic energy.
The prognostic variables are the height hi in each cell i and the normal component of the velocity ue at each edge e between cells. The PV qυ is diagnosed at each vertex υ (or equivalently, in each dual cell). The length of cell edges is ℓe, the distance between cell centers is de and the area of each cell is Ai (Fig. 3). Calculation of the lengths and areas for spherical geometry are stated in section 3b. The cell centers are at locations xi, the vertices at xυ, and the crossover points between edges and the lines between cell centers (the edge points) at xe.
1) Discretization of divergence
The divergence of mass flux in the continuity equation (1) and, by analogy, other divergences are discretized using Gauss’s divergence theorem:
where e ∈ i means all the edges e around cell i and is the height at the edge mapped from the height in the cells either side i and j and the sign of ue is switched if necessary so that positive is outwards. In Ringler et al. (2010) midpoint interpolation is used: . On Voronoi grids, this is a second-order approximation of the value at the edge crossover point but only a first-order approximation of the edge-average value since the crossover point is not at the edge center. An alternative interpolation for for non-Voronoi grids will be defined in section 3b below.
2) Discretization of gradient normal to cell edges
The gradient at the cell edge is needed in the direction normal to the cell edge:
where i and j are the cells either side of edge e. This simple two-point formula on an orthogonal grid guarantees that gradients are curl free around vertices. If the grid were nonorthogonal, this formula would still give curl free gradients but the formula would not be consistent (i.e., the truncation error would not go to zero as the grid is refined). This discrete gradient is the dual of the discrete divergence, which is required for energy conservation.
3) Discretization of curl
The vorticity and, by analogy, other curls are discretized using Stoke’s theorem around the dual cells circulating each vertex:
where Aυ is the area of the dual cell surrounding vertex υ and the sign of ue is switched if necessary so that positive implies anticlockwise circulation.
4) Reconstruction of tangential velocities
The tangential velocities do not always need to be calculated in TRiSK as the entire PV flux term qhu⊥ is discretized together as described in section 6 below. The method of reconstructing tangential vector components is reproduced here as it will be needed in subsection 5 below.
Vector quantities along edges are reconstructed from normal components at the edges of the cells either side of the edge in question using weights that guarantee that the divergence of the vector field on the dual grid is a convex combination of the divergence on the primal grid:
where the edges e′ are the edges of the cells i and j that are either side of edge e. In Thuburn et al. (2009), the weights wee′ were derived to be
where the υs are the vertices in a walk between edges e and e′ and Aiυ is the overlapping area between the dual cell around vertex υ with cell i (see Fig. 4). If this walk starts in the positive direction, then the sign is positive if positive ue′ is outwards.
5) Interpolation of PV from vertices to edges
In TRiSK it is necessary to map the PV from the vertices to the edges in the calculation of the acceleration. Ringler et al. (2010) report two possible methods: the arithmetic mean of the vertices at either end of the edge (i.e., midpoint interpolation) or, to enable the removal of potential enstrophy (but not energy) at the smallest scales, the anticipated potential vorticity method (APVM) of Sadourny and Basdevant (1985), which calculates PV at a point upstream of the edge:
where [∇q]e is the gradient of q and δt is the time step. In this new implementation, the vector [∇q]e is reconstructed from the normal gradients of q on the dual grid using the weights described in subsection 4 above.
6) Discretization of the PV flux
The PV flux qhu⊥ of (2) is averaged between target and surrounding edges to guarantee that the Coriolis force conserves energy:
7) Calculation of kinetic energy
The cell kinetic energy is defined only in terms of the normal velocities:
where . This new implementation will use a modified version for non-Voronoi grids, which will be described in section 3b below.
8) Discretization of ∂/∂t
The momentum and continuity equations are advanced using the explicit fourth-order Runge–Kutta scheme for exact comparison with Ringler et al. (2010, 2008) in test cases in sections 4a and 4e. Semi-implicit Crank–Nicolson time stepping (with no off-centering, nonlinear and Coriolis terms treated explicitly, and two outer iterations per time step) is used in sections 4b and 4c to allow longer time steps to facilitate the use of the full lat-lon grid. However, TRiSK is a definition of the spatial discretization rather than temporal and many other temporal schemes could be used instead.
b. Modifications of TRiSK
1) Spherical geometry
When using the vector-invariant form of the equations on a C grid, spherical polar coordinates are not needed. Instead all areas and distances need to be defined on the sphere. The areas Ai, Aiυ, and Aie (Fig. 4) are defined by decomposing each cell into spherical triangles, with each triangle having a parent cell center xi a parent vertex xυ and a parent edge crossing point xe as its corners and great circle lines between the points. The areas Ai, Aiυ, and Aie can all be decomposed if the grid is Pitteway (i.e., if the line between the cell centers crosses the edge between the cells rather than being off the end). The area of each spherical triangle is
where a, b, and c are the angles between points xi, xυ, and xe on the sphere and . This is the formulation that is less sensitive to rounding error in the limit of small triangles (Williams 2011). The angle a between points xi and xυ is
The edge lengths ℓe (Fig. 3) are great circle distances: . The distances between cells, de, are set by so as not to violate conservation of energy.
2) Non-Voronoi orthogonal grids
Ringler et al. (2010) use midpoint interpolation from cell centers to edges: , which gives conservation of energy if the kinetic energy is calculated using (11). Midpoint interpolation is second-order accurate on Voronoi grids as the edges are midway between the cell centers and it produces a conservative mapping between the cells and the edge areas. Instead we use a conservative mapping between cell centers and edges for Voronoi and non-Voronoi grids:
where i and j are the cells either side of edge e, Aie is the area of the triangle between edge e and cell center i (Fig. 4), and Ae = Aie + Aje. With this modification, a conservative transfer of energy between potential and kinetic can only be maintained if the kinetic energy is defined as
c. Implementation in OpenFOAM
The additional TRiSK operators and Runge–Kutta time discretization have been implemented in the OpenFOAM CFD software library (The OpenCFD Foundation 2011) in order to make use of some existing OpenFOAM operators, grid handling, and linear equation solvers. The divergence and normal gradient of TRiSK are the same as those already in OpenFOAM. The reconstruction of tangential velocities, the curl, the PV flux, the calculation of the kinetic energy, and the APVM interpolation have been newly implemented in OpenFOAM.
Selected results are presented that demonstrate strengths and weaknesses of the different grids and that illuminate grid-imprinting and some computational modes. These will start with a direct comparison with Ringler et al. (2010) showing results of the Williamson et al. (1992) test case 5, the flow over a midlatitude mountain (section 4a). Then steady-state, geostrophically balanced solid body rotation of the linearized and nonlinear SWEs (Williamson et al. 1992; test case 2) will be presented (sections 4b and 4c). The normal mode frequencies and some normal modes are presented in section 4d; these help to explain the SWE results. Finally, results of the barotropically unstable jet of Galewsky et al. (2004) will be presented with a subset of the grids described (section 4e).
a. Flow over a mountain (Williamson et al. 1992; test case 5)
TRiSK has been implemented exactly as in Ringler et al. (2010) and so extensive comparisons have been made with the results in that paper in order to validate the implementation, including extensive diagnostics for the Williamson et al. (1992) test case 2 (not shown). Comparisons are presented here for the vorticity of the flow over a midlatitude mountain after 50 days (Fig. 5) using a hexagonal icosahedral grid of 40 962 cells (see Table 1), giving a cell center to cell center distance of around 120 km, a time step of 100 s, and APVM to remove small-scale potential enstrophy.
For Fig. 5, the vorticity was diagnosed at 50 days in two different ways. First, the vorticity on hexagons was calculated using Stoke’s theorem to circulate around each hexagon and given by (6) (middle row of Fig. 5). This compares well with Ringler et al. (2010) (top left), confirming the correct implementation. However, when the vorticity is calculated on the vertices from the normal velocities using Stoke’s theorem (as is done during the simulation) and plotted on the dual grid of triangles (bottom left), then grid scale oscillations are visible despite the Rossby radius (2525 km) being well resolved and despite the APVM implicit diffusion of vorticity. Calculating the vorticity on the primal (hexagonal) grid rather than the dual (triangular) grid smooths the vorticity, thus hiding the grid-scale velocity oscillations. The vorticity oscillations on triangles are the manifestation of the computational Rossby modes of the hexagonal C grid (Thuburn 2008) caused by the excessive number of velocity DOFs in comparison to height DOFs. Using a smoothed version of the vorticity during the simulation would not help. In fact it would make the mode more difficult to control. The mode occurs because a smoothed version of the vorticity is already used in the momentum equation (the vorticity must be averaged from vertices to edges). The momentum equation can only act to reduce grid-scale oscillations if it can see them.
This level of grid-scale noise is considered small for a 50-day run and may not be detrimental in 3D simulations with parameterized forcing and dissipation (T. Ringler 2011, personal communication). However, it may be possible to do better than this with an advection scheme more accurate than APVM for the potential vorticity (e.g., Weller 2012).
b. Solid body rotation of the linearized SWEs
The linearized, rotating SWEs are
where H is the mean fluid depth and h represents the depth perturbations about H. These simplified equations are solved as well as the full nonlinear SWEs since they support the same waves as the full SWEs, the computational modes are the solutions of the linearized equations, and the solution errors of solving the linearized equations help in understanding the errors in solving the full SWEs.
A steady-state analytic solution for solid-body rotation on the sphere is available:
where φ is latitude and f = 2Ω sinφ. We use H = 2000 m and other parameters from test case 2 of Williamson et al. (1992): a = 6.371 22 × 106 m, Ω = 7.292 × 10−5 s−1, g = 9.806 16 m s−2, and u0 = 2πa/(12 day). These parameters give a gravity wave speed of and the minimum Rossby radius is 960 km (at the poles).
The model is initialized by sampling u, υ, and h at the centroids of the edges and cells rather than at the computational points since these values give a second-order approximation to the fluxes. When initializing from a divergence-free velocity field this means that the discrete velocity is discretely divergence-free up to second order in space. This initialization has been found to be closer to discrete geostrophic balance than sampling at the computational points, especially for the noncentroidal grids.
The velocity and height differences from the initial conditions after 5 days for all the grids are shown in Fig. 6. These solutions use midpoint interpolation of f from vertices to edges for calculating the PV flux from (10). All the grids use Crank–Nicolson time-stepping with a time step of 3600 s, apart from the full rotated lat-lon grid, which uses 900 s in order to achieve a similar Courant number near the pole of the grid. The maximum and minimum height errors and the ℓ2 height error norm are shown in the figure captions of Fig. 6. The ℓ2 error norm is defined as , where h0 is the initial condition equal to the sampled analytic solution for this steady-state case. Since TRiSK can maintain geostrophic balance (it has steady geostrophic modes) and conserves energy in the limit of small Δt, the errors in Fig. 6 are due to the difference between the initial conditions and a discretely balanced state. The sampled initial conditions do not contain grid-scale oscillations and so oscillations generated are due to model errors.
The grid that gives the lowest height errors is the hexagonal icosahedron (Fig. 6e) and the errors have large-scale structure with reduced height at the equator and increased height at high latitudes. There is some large-scale grid imprinting based on the fivefold symmetry of the grid, but no grid-scale oscillations in height. The computational mode on hexagons consists of velocity circulating in alternating directions around each vertex (Thuburn 2008) that is present in Fig. 6e.
The unrotated and 45° rotated full lat-lon grids (Figs. 6a,b) also have low height errors, without grid-scale oscillations, although larger height errors are localized near the pole of the 45° rotated lat-lon grid. In an Eulerian operational lat-lon model, these oscillations would be removed by polar filtering. The unrotated skipped lat-lon grid has additional errors at the grid reductions (Fig. 6c), which become grid-scale oscillations when the skipped grid is rotated (Fig. 6d). These oscillations are related to computational modes with grid-scale oscillations on the dual, which will be discussed later.
The sampled initial conditions have triggered large grid-scale oscillations in height on the triangular grid (Fig. 6f), similar in appearance to the spurious inertio-gravity modes on the triangular C grid described by Danilov (2010). The velocity errors on triangles in Fig. 6f are mostly in the tangential component, with velocity circulating in alternating directions around each triangle. This pattern is present in the initial conditions due to the initial TRiSK reconstruction of velocity.
The height and velocity errors on the kite grid are larger than on any other grid (Fig. 6g). In common with other grids, there are negative height errors at the equator but elsewhere the error is dominated by grid-scale oscillations. These are likely to be due to the larger truncation errors of using kites, resulting from the grid inhomogeneities and anisotropies described in section 2d.
The large height errors on the Voronoi-ized cube (Fig. 6h) are mostly large scale rather than grid scale with some large-scale grid imprinting with fourfold symmetry and velocity aligning with the grid rather than along lines of latitude.
More insight can be gained into the errors on different grids by inspecting the relative vorticity after 5 days defined on the dual grid (Fig. 7) (i.e., defined on the vertices of the primary grid). The ℓ2 error norms of PV are shown in the captions where q0 is the initial value of q calculated from the sampled analytic values of u and h.
The unrotated full lat-lon grid has the smallest ℓ2 error norm and so the relative vorticity in Fig. 7a can be regarded as nearly exact. As for the height, the rotated full lat-lon grid has errors around the pole of the grid (Fig. 7b) where small-scale features are localized in the fine resolution. Once grid reductions are included (Fig. 7c) grid-scale oscillations are introduced with high vorticity on the triangles of the dual. This corresponds to cyclonic vorticity on the triangular dual cells that occur at the grid reduction. This is due to the growth of computational Rossby modes resulting from the excessive number of velocity to height at the grid reductions. (The vorticity oscillations must be computational modes since they are stationary and real Rossby modes at this scale would propagate.) These grid-scale oscillations become much worse when the grid is rotated (Fig. 7d) since the mode is now forced by the flow across the grid reductions. The alternating cyclonic and anticyclonic flow forces the grid-scale oscillations in height as seen in Fig. 6d.
Grid-scale oscillations consisting of just two values can occur on grids of squares (checkerboard patterns) and also on triangles but two values cannot create grid-scale oscillations on hexagons on which two adjacent hexagons never share the same value. This implies that 2D grid-scale computational modes cannot be present when values are represented on hexagons. But for the hexagonal C grid, the vorticity lives on the vertices (or dual triangles). Therefore computational modes are present in the vorticity on triangles (Fig. 7e). These do not directly reduce the accuracy of the height field (Fig. 6) since the vorticity is averaged in the Coriolis terms, so the noise is smoothed in its effect on the dynamics. The variable that lives on hexagons cannot have grid-scale oscillations and so the vorticity on the hexagonal dual of the triangular grid (Fig. 7f) is free of grid-scale oscillations. It also has a low ℓ2 error.
The worst grid-scale vorticity oscillations occur for the kite grid (Fig. 7g) whose dual consists of a variety of shapes, with different vorticities on each shape. These grid-scale oscillations are stationary and are therefore computational modes.
The results on the Voronoi-ized cube (Fig. 7h) have similarities with those on the hexagonal grid (Fig. 7e). Both grids are Voronoi with triangular duals and there are grid-scale oscillations on the vorticity on the triangles. These are worse on the Voronoi-ized cube since the grid is less isotropic, with a local mixture of long and short edges. The oscillations change phase around the cube as the grid bends north and south. The dual of the Voronoi-ized cube does not completely retain the symmetry of the cube since, due to the nonuniqueness of the Delaunay triangulation and rounding error, each face of the cube is triangulated differently.
c. Solid body rotation of the nonlinear SWEs
The nonlinear SWEs are solved for Williamson et al. (1992) test case 2 on the same grids and input parameters as section 4b above. The height and velocity errors after 5 days are shown in Fig. 8. These simulations have midpoint interpolation of PV from vertices to edge points (8).
When solving the linearized SWEs, all the grids tended to give negative height errors around the equator and positive toward the poles (Fig. 6). This tendency does not occur for the nonlinear equations (Fig. 8). This could be due to truncation errors of the opposite sign introduced with approximating the nonlinear terms [nonlinear PV flux and kinetic energy gradient in (2)]. Therefore, there could be a cancellation of errors, and the height errors are reduced on the full and skipped unrotated lat-lon grids when solving the nonlinear equations (Figs. 8a and 8c vs Figs. 6a and 6c). However, all the other grids give larger errors. When solving the nonlinear equations, velocity errors that are generated at a grid inhomogeneity are advected away, increasing the errors downstream. In particular, the hexagonal grid shows more large-scale grid imprinting and large errors are generated away from the grid reductions of the rotated skipped grid. The grid-scale errors on the triangular and kite grids retain the same structure but deepen and the errors on the Voronoi-ized cube deepen and display more large-scale grid imprinting.
The vorticity errors that result from solving the nonlinear equations are similar to those that result from solving the linear equations but they grow more quickly due to nonlinear feedbacks onto the computational modes. They are not presented since they do not give further insight.
d. Normal modes
Normal modes and their frequencies are calculated in order to answer the following questions:
Do all the grids have stationary geostrophic modes?
Are these stationary modes physical (in which case they will have zonal symmetry) or computational (consisting of grid-scale oscillations)?
Do grids with the correct ratio DOFs of suffer less from computational modes?
How do the computational modes relate to the errors when simulating steady, geostrophically balanced flow?
How does the grid-scale heterogeneity influence the wave modes?
The normal modes of the linearized SWEs and frequencies can be found for any grid with any numerical algorithm by calculating the eigenvectors and eigenvalues of the matrix which represents the linearized action of the model on any set of initial conditions. The initial conditions are represented as a vector, (h, un)T, where h is the vector of the h values in every cell and un is the vector of the un values normal to every edge. The model matrix is found by running the model for one short time step for each of the set of initial conditions (1, 0, … 0)T, (0, 1, 0 … 0)T, … (0, 0, … , 0, 1)T and the solutions form the columns of . The short time step is to ensure that nonlinear effects are minimal. The eigenvectors v and the eigenvalues λ of are computed and then the vs are the normal modes with frequencies ω calculated from the eigenvalues, λ = αeiωδt. The normal modes are calculated using the parameters given by Thuburn et al. (2009): Earth’s radius, a = 6 371 220 m, constant Coriolis, f = 1.4584 × 10−4 s−1, gH = 105 m2 s−2, and using a time step of 10 s. Coarser versions of the grids shown in Fig. 1 are used due to the computational expense of the eigendecomposition. Grid specifications are given in Table 1. The Rossby radius is 2168 km, which is marginally resolved. The zero frequency modes of the rotating sphere are presented first and then we will look at the remaining modes and compare with the analytic frequencies.
1) Stationary modes of the rotating sphere
The stationary modes are particularly important since they are either the physical, geostrophically balanced states (which are zonally symmetric on the rotating sphere) or they are stationary grid-scale computational modes, which are damaging to the solution. These are the modes associated with nontrivial null spaces in one or more operators of the governing equations. For each of the grids, one of the stationary modes is shown in the first column of Fig. 9 (if it exists) and a stationary or very slow mode is shown in the second column. The second mode is chosen to contrast in structure with the first.
Only the lat-lon grid has stationary geostrophic modes in perfect zonal symmetry as only this grid can represent perfect zonal symmetry (top left of Fig. 9). The other zero-frequency mode shown for the lat-lon grid is a combination of zonally symmetric geostrophic balance and the Coriolis mode (with spurious stationary vortices around each vertex, zero Coriolis force, and zero divergence).
The skipped lat-lon and hexagonal grids have very low-frequency modes in near-zonal symmetry (second column of Fig. 9). The exact zero-frequency modes are contaminated by computational modes: the Coriolis mode on the skipped lat-lon grid and a spurious Rossby mode on the hexagonal grid.
An exact zero-frequency mode with some zonal symmetry has been found for the triangular grid (second column of Fig. 9) but, contrary to the theory of the triangular C grid on a beta plane, there is also a spurious stationary Rossby mode on triangles (first column). However, in spherical geometry with spherically varying Coriolis, different spurious solutions are possible and it has been possible to find a velocity field satisfying ∇ ⋅ u = 0 on the cells and ∇ ⋅ fu = 0 on the vertices.
It was hoped that the kite grid, with the correct ratio of velocity to height DOFs, would not suffer from the spurious Rossby mode branches of the hexagonal grid. However, the kite grid has spurious stationary rotational modes (first column of Fig. 9) as well as spurious very low-frequency rotational modes (second column) and none of the slow modes is zonally symmetric.
The Voronoi-ized cube does not have any exact zero-frequency modes or slow modes in near-zonal symmetry, just the spurious slow Rossby modes of the hexagonal C grid.
2) Nonstationary modes of the f sphere and rotating sphere
Normal modes are calculated on the f sphere as well as the rotating sphere. The f sphere is a purely mathematical concept—a sphere with globally uniform Coriolis. This is used because there is an analytic dispersion relation for the frequency ω for each spherical harmonic with wavenumber n:
The rotating sphere has the usual sinusoidal variation of f with latitude.
The f sphere model’s frequencies are compared with the analytic frequencies from (18) and are shown in Fig. 10. Since the frequencies occur in complex conjugate pairs, only half of them are shown. Ideally the frequencies would be plotted against wavenumber as a dispersion relation. However, grid inhomogeneity causes many mode structures to become localized, making it impossible to define a meaningful wavenumber. Instead, for the f-sphere results, the frequencies are simply sorted into zero values and nonzero values in order to distinguish geostrophic modes from inertio-gravity modes, and the nonzero values are sorted into ascending order. Unfortunately physical modes cannot be distinguished from computational modes on the basis of frequency alone because their frequency spectra overlap.
For the rotating sphere (red frequencies in Fig. 10) the Rossby modes (low frequencies) are separated from the inertio-gravity modes (high frequencies) purely by ordering the frequencies and by assuming that the rotating sphere has the same number of Rossby and geostrophic modes as the f sphere has zero-frequency geostrophic modes. This separation is arbitrary and some modes may appear in the wrong branch. The rotating sphere has less quantization of mode frequencies because of the spatial variation of f.
For the f sphere, the hexagonal grid (Fig. 10c) and Voronoi-ized cube sphere (Fig. 10f) have nearly the same number of zero modes (geostrophic modes) and nonzero modes (the inertio-gravity modes) since, as explained by Thuburn et al. (2009), the number of geostrophic modes is given by the number of vorticity degrees of freedom (the number of vertices minus one) and the number of inertio-gravity modes is the number of mass plus divergence degrees of freedom (twice the number of cells minus one). These Voronoi grids have three cells meeting at one vertex and most cells are hexagonal so the number of each type of mode is about equal. The wave frequencies are similar for these two grids except that there is less quantization of frequencies for the Voronoi-ized cube sphere f sphere than for the hexagonal icosahedron due to the variations in cell sizes and edge lengths for the Voronoi-ized cube sphere. The spurious Rossby modes on the hexagonal f sphere are zero frequency and so are on the x axis.
For each of the grids, the spatial structure of two more of the nonstationary eigenmodes of the rotating, linearized SWEs is plotted in Fig. 9. One of the Rossby modes with large spatial scale and frequency around 10−5 s−1 is plotted for each grid in the third column. The mode with the highest frequency is in the final column.
All of the grids have a similar-looking wavenumber 1 Rossby mode that is symmetric about the equator (third column of Fig. 9). This mode is presented as a sanity check that the analysis reproduces the large-scale wave modes that we expect on the rotating sphere. All grids capture the mode with comparable frequency, except the coarser lat-lon grid.
The final column of Fig. 9 shows the highest-frequency mode which, for most grids is associated with the smallest cells. All of the grids have some variation in cell size and consequently have higher-frequency modes localized in the highest resolution regions. For some grids (lat-lon and kites), the variation in cell size is so extreme that there is an upturn at the end of the frequency graphs (Figs. 10a,b,e) due to gravity waves localized near inhomogeneities of the grid. For the triangular grid, the highest-frequency mode is one of the spurious inertio-gravity modes of the triangular C grid, with largest amplitude at the pole where the Rossby radius is smallest.
For the triangular grid (Fig. 10d) and the kites (Fig. 10e) there is a jump in the frequency part way along the inertio-gravity modes. This change in slope was also recognized for triangles by Thuburn et al. (2009), who attributed this to a switch to a spurious branch of inertio-gravity modes. However, the nature of the jump in frequency for the kite grid is different; beyond this jump, all of the modes are regionally confined (not shown). Because of the inhomogeneous and anisotropic nature of the kite grid, modes become localized in many locations of the grid rather than just close to the two poles of the lat-lon grids. This would be a problem for an operational model: if convection occurred in one of the fine-resolution regions of the kite grid, the resulting gravity waves could not propagate cleanly away.
The Voronoi-ized cube has frequencies that increase approximately linearly (Fig. 10f) for the highest frequencies whereas we would expect the frequencies to flatten out for a uniform, structured C grid with linear differencing (e.g., Randall 1994). So it is likely that this apparently realistic increase in frequency is in fact a computational artifact. This is confirmed by inspecting the highest-frequency mode on this grid in Fig. 9, which has grid-scale oscillations in a limited region of the grid. These high-frequency waves will propagate around the high-resolution regions of the grid rather than propagating out of them.
We can now answer the questions posed at the beginning of this subsection concerning wave modes on the rotation sphere:
All of the grids have stationary geostrophic modes on the rotating sphere, apart from the Voronoi-ized cube.
All of the exactly stationary geostrophic modes on the rotating sphere are spurious computational modes, apart from on the full lat-lon grid, where perfect zonal symmetry is possible.
The triangular, kite, and Voronoi-ized cube grids do not have any near-stationary modes in zonal symmetry. However, the hexagonal grid has a very low-frequency mode in near-zonal symmetry.
The kite grid, with the correct ratio of DOFs, suffers from many computational modes. Merely having the correct ratio of DOFs does not eliminate these.
The grid-scale vorticity errors on hexagons are caused by the stationary computational Rossby mode (first column of Fig. 9).
Heterogeneity and anisotropy at the grid scale leads to wave modes localized near the finest parts of the grids, for all grids. The kite grid has a large proportion of the wave modes localized in different parts of the grid.
e. Barotropically unstable jet (Galewsky et al. 2004)
The simulation of a barotropically unstable jet with an initial perturbation (Galewsky et al. 2004) is a tough test for a low-order model on a non-lat-lon grid since the jet is fast and narrow and numerical truncation errors lead to perturbations that release the instability in a similar manner to the initial perturbation. This test case therefore needs higher resolution than solid-body rotation in order to achieve results qualitatively similar to the high-resolution solution. A high-resolution solution for the relative vorticity after 6 days is shown at the top of Fig. 11 (using a hexagonal icosahedron of 163 782 cells, δx ≈ 60 km, as given in Table 1) and is compared with the results from other the grids at lower resolutions; 70 to 140 km (also in Table 1). Only the icosahedra and skipped lat-lon grids were used since from sections 4b, 4c, and 4d we can see no advantages of the kite and Voronoi-ized cube over the icosahedral grids. The overlaid dual grids are a factor of 8 coarser in both directions for clarity.1
This test case has a Rossby radius of 2525 km based on the midjet values, which is well resolved. The time step is 100 s, which implies that the flow Courant numbers based on all the primal and dual grids are well below one but the Courant number based on the gravity wave speed is close to one. APVM (9) is used for suppressing grid-scale vorticity oscillations and dissipating potential enstrophy.
An attempt is made to find initial conditions in discrete balance so that initial grid-scale divergence and geostrophic imbalance are not larger than the initial perturbations. The initial velocity is found as in Ringler et al. (2011) by sampling ψ at the vertices from (17) then setting ue = k × ∇ψ. This ensures that the initial velocity is divergence free and therefore initially If additionally it is possible to find h such that fu⊥ = −g∇h then there should be no drift from the initial conditions. This h is sought by solving the Poisson equation: −g∇2h = ∇ ⋅ fu⊥ using the TRiSK operators. The Poisson equation can be solved to arbitrary tolerance but this does not guarantee that discrete balance holds, just the divergence of the discrete balance equation; it is only possible to find h satisfying the discretized fu⊥ = −g∇h if the discretized ∇ ⋅ fu = 0 at the vertices. This is because taking the curl of fu⊥ = −g∇h implies that ∇ × fu⊥ = ∇ ⋅ fu = ∇ × ∇h = 0. However, the initial conditions found are closer to discrete balance than sampled initial conditions.
The unrotated, skipped lat-lon grid (Fig. 11b) gives results closest to the high-resolution solution because the flow is aligned with the grid and the jet does not cross the first grid reduction (at 66°N) until day 5 and so the additional truncation errors caused by the grid reduction have not had sufficient time to amplify and contaminate the solution. Spikes in the vorticity can be seen where the jet crosses the first grid reduction in the zoomed region.
When the skipped lat-lon grid is rotated by 90° (Fig. 11c), the poles are now along the equator and the grid reductions are near the equator so the narrow midlatitude jet does not reach the first grid reduction until day 5. Therefore, the grid reduction does not have a big impact on the solution at day 6; the reduced accuracy, phase error of the largest wave, and enhanced spurious release of barotropic instability for the 90° rotated grid are due to the misalignment of the jet with the grid and the change in angle between the grid and the jet around the sphere. It is only reasonable to compare non-lat-lon grids with rotated lat-lon grids since real jets are not aligned with the grid and so the optimal accuracy of the aligned lat-lon grid will never be achieved. So if the quasi-uniform grids can perform as well as the 90° rotated lat-lon grid then they may be sufficiently accurate for operational forecasting. When the skipped grid is rotated by 45° (Fig. 11d) more damage to the solution is done when the jet passes through the grid reductions. APVM imposes an additional time-step restriction and so this case has not been simulated using a full lat-lon grid.
At 120-km resolution, the hexagonal grid does suffer from some spurious release of instability after 6 days (Fig. 11e), similar to the 90° rotated skipped lat-lon grid, but it does not suffer from serious phase error. For TRiSK, which uses the cheapest possible two-point differencing, this is an impressive result that compares well with a collocated model on the same grid which uses much more expensive quadratic differencing (e.g., Fig. 7 of Weller et al. 2009). The spurious checkerboard pattern in vorticity on the dual grid appears for this test case using the hexagonal grid (see zoomed region in Fig. 11e) although APVM has damped this computational mode. This numerical artifact is not visible when viewing the vorticity on the primal grid of hexagons (not shown).
At 70-km resolution, the triangular grid results (Fig. 11f) are worse than the hexagonal grid results at 120 km; larger waves have grown all around the hemisphere and the position of the largest wave has a larger phase error (clear from zoomed region). This is not due to the computational mode of triangles or the grid-scale heterogeneity since there is no grid-scale noise in the height on triangles (not shown). The computational modes on triangles do not appear. Based on numerous other simulations on triangles, this seems to be because these modes are not strongly forced when using balanced initialization. The reduced accuracy of triangles relative to hexagons in this case is therefore a consequence of the lower-order accuracy two-point interpolation on triangles since they are not Voronoi.
For this test case, the hexagonal grid and the lat-lon grid are giving the most accurate results (so long as the jet does not interact with the lat-lon grid reductions). This crucial result implies that the hexagonal grid is as good as the full latitude–longitude grid (since if the jet does not interact with the lat-lon reductions then the results will be as good as the full lat-lon grid). However the hexagonal grid is quasi-uniform and therefore will not have the same parallel scaling problems as the full lat-lon grid. The skipped lat-lon grid leads to large errors when the jet passes through the grid reductions.
TRiSK on all five quasi-uniform grids suffers from computational modes that are usually associated with triangles in the primal or dual grid. The skipped lat-lon grids have triangles in their dual which develop spikes in the vorticity. The dual of the hexagonal icosahedron is triangular and so computational modes consisting of grid-scale vorticity oscillations on the triangles exist [the computational Rossby modes of Thuburn (2008)]. The triangular grid has the computational mode in height on the primal grid. The dual of the Voronoi-ized cube consists of anisotropic triangles and so vorticity oscillations grow on the triangles. The kite grid, despite having the correct ratio of DOFs, suffers from more grid-scale oscillations than any of the other grids. The resolutions of velocity and mass are anisotropic on the kites and so grid-scale noise and high-frequency waves are localized in clusters of cells.
Using TRiSK, computational modes involving vorticity oscillations can be damped using a diffusive advection scheme for PV such as APVM, which leads to a dissipation of potential enstrophy. This does not remove energy or violate any of the other desirable properties of TRiSK, unlike divergence damping, which would be needed to control spurious inertio-gravity modes on triangles. Therefore, computational modes consisting of vorticity oscillations on triangles appear to be easier to deal with. Given this, the hexagonal icosahedron gives the best results on all of the test cases.
The skipped lat-lon grid performs almost perfectly when the flow is aligned with the grid and well when the grid is rotated. However, when the flow interacts with the changes of resolution of the skipped lat-lon grid, much larger errors are generated, which are then advected globally.
The cubed sphere was made perfectly orthogonal while retaining quasi-uniformity by making it into a Voronoi tessellation. Hence the cells are composed of polygons rather than squares. However, this implies that the dual now consists of triangles and so computational Rossby modes exist on the triangles. The triangular dual is less isotropic than the dual of the hexagonal grid and so the errors are larger.
All of these conclusions are relevant to the arbitrarily structured, orthogonal C grid (TRiSK), which is a low-order method with stationary geostrophic modes and conservation of mass, energy, and PV and which is consistent and compatible. No consistent numerical C-grid method currently exists with these properties on nonorthogonal grids. These properties should enable accuracy close to that of the lat-lon C grid. However, the conclusions will be relevant to other methods on these grids. In particular, the computational modes that can be supported when the velocity and mass do not have the same resolution will be problematic for all C-grid numerical methods.
The final conclusion is that TRiSK was developed for hexagonal icosahedra and other Voronoi tessellations that are close to centroidal (Ringler et al. 2011), and indeed TRiSK works best on the hexagonal icosahedra and works very well in comparison to other, more expensive methods. It therefore appears to be a good choice for a quasi-uniform operational weather and climate forecasting model.
H.W. acknowledges support from NCAS Climate and NERC Grants NE/H002774/1, NE/H001166/1, and NE/I022086/1. J.T. acknowledges funding from NERC Grants NE/J005346/1 and NE/H002774/1. C.C. acknowledges funding from NERC Grants NE/I02013X/1 and NE/I016007/1. We also are grateful for discussions with Nigel Wood, Andrew Staniforth, and Terry Davies of the UK Met Office and the insightful comments of three anonymous reviewers.
The hexagonal icosahedral grid at this resolution should have 163 842 cells rather than 163 782. However, during the iterations to make the grid centroidal while retaining uniform resolution, 60 of the cells have vanished. This is unlikely to have a big impact on the solution accuracy.