GHOST : Geoscientific Hollow Sphere Tesselation

I present in this work the GHOST (Geoscientific HOllow Sphere Tesselation) software which allows for the fast generation of computational meshes in hollow sphere geometries counting up to a hundred millions of cells. Each mesh is composed of concentric spherical shells which are built out of quadrilaterals or triangles. I focus here on three commonly used meshes used in geodynamics/geophysics and demonstrate the accuracy of shell surfaces and mesh volume measurements as a function of resolution. I further benchmark the built-in gravity and gravitational potential procedures in the simple case of a 5 constant density geometry and finally show how the produced meshes can be used to visualise the S40RTS mantle tomography model. The code is open source and is available on the GitHub sharing platform.


Introduction
In the last 40 years, numerical mantle convection studies have improved our understanding of mantle dynamics as a whole (Schubert et al., 2001).While early studies looked at aspects of fluid dynamics (Busse, 1975;Christensen and Harder, 1991), more recent studies have been exploring a wide variety of topics -for example, mantle mixing (van Keken et al., 2002), melting (Tackley, 2012;van Heck et al., 2016;Dannberg and Heister, 2016), the effect of plate motion history on the longevity of deep mantle heterogeneities (Bull et al., 2014) or assimilating lithosphere and slab history in 4-D Earth models (Bower et al., 2015).
To a first approximation, the Earth is a sphere: the Earth's polar diameter is about 43 km shorter than its equatorial diameter, a negligible difference of about 0.3 %.As a consequence, modelling physical processes which take place in the planet require the discretisation of a sphere.Furthermore, because core dynamics occur on vastly difference timescales than mantle dynamics, mantle modelling usually leaves the core out, thereby requiring simulations to be run on a hollow sphere mesh (with the noticeable exception of Gerya and Yuen, 2007).
Although so-called latitude-longitude grids would seem appealing, they suffer from the convergence of meridians at the poles (resulting in oversampling at poles) and the juxtaposition of triangles near the poles and quadrilaterals elsewhere.As a consequence, more regular, but more complex, grids have been designed over the years which tessellate the surface of the sphere into triangles or quadrilaterals (sometimes overlapping).There is the "cubed sphere" (Ronchi et al., 1996;Choblet et al., 2007), the yin-yang grid (Kageyama and Sato, 2004;Yoshida and Kageyama, 2004;Kameyama et al., 2008;Tackley, 2008;Crameri andTackley, 2014, 2016), the spiral grid (Hüttig and Stemmer, 2008), an icosahedron-based grid (Baumgardner and Frederickson, 1985;Tabata and Suzuki, 2002) or a grid composed of 12 blocks further subdivided into quadrilaterals (Zhong et al., 2000) as used in the CitcomS code.Note that Oldham et al. (2012) have also presented a method for generating a numerical grid on a spherical surface which allows the grid to be based on several different regular polyhedrons (including the octahedron, cube, icosahedron and rhombic dodecahedron).
How such meshes are built is often not discussed in the literature.It is a tedious exercise of three-dimensional geometry and it can be time-consuming, especially the connectivity array generation.In this paper, I present an open-source mesh generator for three hollow sphere meshes: the "cubed sphere" mesh, the CitcomS mesh and the icosahedral mesh.
I first present the basic workflow which has been implemented to arrive at such meshes; then, I showcase its efficiency and how accurate surfaces and volumes are represented.Finally, I provide a simple example of gravity and  gravity potential calculations on such meshes and compare the obtained values with the analytical solution derived in the Appendix.

Building the hollow sphere meshes
The open-source code library GHOST allows three different types of hollow sphere meshes to be built, i.e. meshes bounded by two concentric spheres: -The cubed sphere (HS06) is composed of six blocks which are themselves subdivided into N b × N b quadrilateral-shaped cells (Sadourny, 1972;Ronchi et al., 1996;Hernlund and Tackley, 2003;Burstedde et al., 2013).Four types of cubed sphere meshes have been proposed: the conformal, elliptic, gnomonic and spring types (Putman and Lin, 2007).However, only gnomonic meshes are considered here: these are obtained by inscribing a cube within a sphere and expanding to the surface of the sphere.The cubed sphere has recently been used in large-scale mantle convection simulation in conjunction with adaptive mesh refinement (Alisic et al., 2012;Burstedde et al., 2013).
Given the regularity and symmetry of these meshes, determining the location of the mesh nodes in space is a relatively straightforward task.Building the mesh connectivity in an efficient manner is where the difficulty lies.
The approach to building all three meshes is identical: 1.A reference square or triangle is populated with cells, as shown in Fig. 1, parameterised by a level l: the square is subdivided into l × l quadrilaterals, while the triangle is subdivided into l 2 triangles (implemented in the block_node_layout subroutine).
2. This reference square or triangle is then replicated nblock times (6, 12 or 20) and mapped onto a portion of a unit sphere.The blocks are such that their union covers a full sphere, but they cannot overlap except at the edges; see Fig. 2.This takes place in the map_blocks subroutine.
3. All block meshes are then merged together to generate a shell mesh.This task is rather complex as duplicate nodes must be removed and all connectivity arrays of the blocks must then be mended accordingly.This task is carried out in the merge_blocks subroutine.
4. Shell meshes are replicated nlayer+1 times outwards with increasing radii.The nlayer shells are then merged together to form a hollow sphere mesh, as shown in Fig. 3.This is carried out in build_hollow_sphere.
More information on these steps is available in the manual of the code.In Table 1, the number of nodes and cells for a variety of resolutions for all three mesh types is reported.Looking at the CitcomS literature of the past 20 years, we find that the mesh data presented in this table cover the various resolutions used, e.g. 12 × 48 3 (McNamara and Zhong, 2004;Arrial et al., 2014), 12×64 3 (Bull et al., 2014) 12×96 3 (Bull et al., 2010) and 12 × 128 3 (Becker, 2006;Weller and Lenardic, 2016;Weller et al., 2016).Note that, in the case of the HS06 and HS12 meshes, the mesh nodes are mapped out to the 6 or 12 blocks following either an equidistant or equiangular approach as shown in Fig. 6 (see Putman and Lin, 2007 for details on both approaches).
It is also worth mentioning that the aspect ratio of the cells usually influences the accuracy of the computed solution.As such, the user should strive to generate a mesh where cells' aspect ratios are (on average) close to unity by making nlayer a function of l.

Mesh generation performance
The total time to generate the coordinates and the connectivity of the final hollow sphere mesh was timed for all three mesh types and is shown in Fig. 4. The reported times scale ideally with the number of nodes, i.e. linearly or O(N ), up to 100 million mesh nodes and a mesh containing a million or so nodes can be built in less than a second on a laptop.

Areas and volume measurements
The area of the shell of unit radius can be calculated by summing the areas of each cell.In the case of triangles, Heron's formula is used, which states that the area of a triangle whose sides have lengths a, b and c is where s is the semi-perimeter of the triangle: (2) In the case of quadrilaterals, the four points composing each of them are not necessarily coplanar and the definition of the surface is ill-posed.Each quadrilateral is therefore decomposed into four triangles sharing a vertex in the middle given by the barycentre of the four points.The area of each quadrilateral is then approximated by the sum of the areas of all four inscribed triangles.
The relative error on the shell outer area, e A = (A meas − A th )/V th , where A meas is the total measured area and A th = 4π R 3 2 , is shown in Fig. 5, and the error is found to decrease linearly with the number of points for all three types of meshes.
Figure 6 shows shells which approximately count the same number of nodes.We see that the equiangular projection of the nodes for the HS06 and HS12 meshes yields cells whose area is homogenous in value compared to when the equidistant projection is used.
Although the volume of the hexahedra could have been computed directly with the formula of Grandy (1997), a Gaussian quadrature rule is used since it is also needed in the next section.The total volume of the spherical mesh is given by where stands for the volume inside the spherical shell; c is the volume of a cell.The sum runs over all the cells c = 1, . ..N el , and a 2 × 2 × 2 quadrature rule is used in each cell.
The measured hollow sphere mesh volume and its relative error e V = (V meas − V th )/V th , where V meas is the total measured volume and V th = 4π(R 3 2 − R 3 1 )/3, are shown in Fig. 7a, b, and the error is found to decrease with the number of points for all three types of meshes with a −2/3 exponent.

Gravity field and potential measurements
The gravity potential U can be computed by means of the Poisson equation (∇ 2 U = 4π Gρ), where ρ is the density and G is the gravitational constant (see Turcotte and Schubert, 2012, Chap. 5).The analytical solution of this equation in the case of a constant density spherical shell is given in Appendix A.
The gravity vector and potential at location r can also be computed by means of volume integrals: and both expressions can then be computed using the numerical quadrature described in the previous section.
In what follows, I set the inner radius to R 1 = 1 and the outer radius to R 2 = 2; the density of the material is ρ 0 = 10 6 and G = 6.673848 × 10 −11 m 3 kg −1 s −2 .Figures 8 and 9 show the gravity g = |g| and potential U measured at 256 points along a line between r = 0 and r = 4R 2 for all three mesh types.Data points align along the analytical curves, and the error is found to decrease as a function of the mesh resolution.4 Application: visualisation of a tomography dataset The S40RTS model is one of the most widely used tomographic models of the mantle (Ritsema et al., 2011).
It is based on a 20 million Rayleigh wave dispersion, 500 000 shear-wave travel time and 1100 normal-mode splitting function measurements.The data are widely available, for instance, on the website of the main author (http: //jritsema.earth.lsa.umich.edu//Research.html,last access: 7 October 2018) or as part of the SPECFEM 3Dglobe code (https://github.com/geodynamics/specfem3d_globe,last access: 7 October 2018).I have adapted the Fortran interface provided with the dataset and for each node of the grid the shear-velocity variation δ ln V s is computed, as well as the relative density variation δ ln ρ = ξ δ ln V s , where ξ = 0.25 is assumed to be constant with depth for simplicity (see Fig. 6 of Steinberger and Calderwood, 2006).The absolute density variation with regards to the PREM model (Dziewonski and Anderson, 1981) is then obtained as follows: δρ = ρ PREM × δ ln ρ. Results are shown in Fig. 10.

Conclusions
The three types of hollow sphere meshes presented in this work are currently in use in the ELEFANT code (http:// cedricthieulot.net/elefant.html,last access: 7 October 2018).Furthermore, the HS12 mesh was recently used in Thieulot (2017) in which a family of analytical solutions for viscous incompressible Stokes flow in a spherical shell is presented.
Following the example of CitcomS, each block of the final mesh could actually be built and used by a different Message Passing Interface (MPI) thread in the context of parallel calculations (Burstedde et al., 2013).Each block could then subsequently be divided to allow for more threads to be used than the original number of blocks.
Other tomography models than S40RTS (Ritsema et al., 2011) could have been chosen, such as UUP07 (van der Meer et al., 2017;Hall and Spakman, 2015) and other geophysical databases could have been coupled with it, such as the crust and lithospheric model, Litho1.0 (Pasyanos et al., 2014), to arrive at a more complete high resolution of the Earth.Gravity (anomaly) and geoid measurements could then be carried out.
Finally, this library is aimed at students and researchers alike.It provides the essential building block for many geophysical applications: a non-overlapping tessellation of the mantle.It also provides the starting point for any finiteelement or finite-volume method-based geodynamical code.Which of the three mesh types is best may be problemspecific, and potential users are encouraged to explore various combinations of resolutions and mesh types.However, under the assumption that the quality of the numerical solution correlates with the uniformity of the mesh cells' volume, it appears that an equiangular projection should always be preferred to an equidistant one.

Figure 1 .
Figure 1.Reference square and triangle meshes at level 5 .

Figure 2 .
Figure 2. From left to right: HS06, HS12 and HS20 shells coloured by block number.

Figure 4 .Figure 5 .
Figure 4. Measured times to build and assemble the mesh as a function of its number of nodes (N ).

Figure 7 .Figure 8 .
Figure 7. (a) Measured total volume for all three shells as a function of the number of nodes (N).(b) Corresponding total volume relative error as a function of N.

Figure 9 .
Figure 9. Gravity potential U as a function of the radial coordinate r.The gray area symbolises the spherical shell.

Figure 10 .
Figure 10.HS06 grid with 64 layers and level = 64.(a) δ ln V s at the surface of the model (R = 6346 km) and at depth down to the core mantle boundary (R = 3480 km); (b) absolute density variation δρ (kg m −3 ).