The Cubed Sphere Grid		Tetrahedral steady solution		Time-dependent convection

Numerical study of thermal convection in planetary mantles is complicated by their spherical shape and highly variable properties. These two difficulties are best approached using a grid or mesh based numerical method such as finite differences or finite elements, so that variable properties can be dealt with locally. Non-local numerical approaches, such as spectral schemes, often have difficulties dealing with highly variable properties such as diffusivities or viscosities, and are of little use in the more realistic applications to planetary interiors. Spherical polar grids (like the latitude and longitude system) are commonly used to cover the sphere, however the grid is highly non-uniform and exhibits a singularity at the north and south poles. Both of these issues complicate numerical implementation (eg. multigrid acceleration is made more difficult), and give results having highly spatially dependent resolution and Courant limited time steps. In order to develop a local finite volume/finite difference code for spherical geometries, we are making use of the cubed sphere coordinate system.

The cubed-sphere is a coordinate system obtained by projecting the sides of a cube onto a sphere (hence the name). It has six sections (blocks) on which an equi-angular grid can be constructed, where the grid lines on the sphere form great circles. The blocks are made to overlap one another and they are connected by interpolating a neighboring block's interior points to the exterior "ghost points" for each block. Each block is geometrically identical, due to the cubic symmetry, which greatly simplifies calculations. I am using the cubed sphere as a composite coordinate system (along with a contravariant basis) in the calculations, not simply as a mesh generator. Equations are solved by transforming the usual vector differential operators into the cubed sphere coordinate system using tensor calculus and applying the appropriate numerical approximations.

The code currently uses multi-grid acceleration methods for solution of the elliptic Stokes equations, and an upwind conservative MPDATA scheme for advective transport with or without diffusion. The multi-grid algorithm is based on the Full Approximation Scheme (FAS) and includes pre-defined (a priori) grid refinements as a natural part of this multi-grid structure when combined with local high order interpolation. In the near future a fully adaptive dynamically determined grid refinement scheme will be implemented, and tau-extrapolation (a form of Richardson extrapolation) will be used to increase the order of approximation for a certain class of problems.

The multigrid solver will feature several high performance features that allow it to incorporate very large viscosity variations relative to competing codes. This is motivated by the fact that the commonly used inter-grid transfer operators (for restriction/prolongation) in most mantle convection codes are actually of first-order, because of the highly variable or discontinuous nature of the fine grid operator coefficients. The usual claim that a bilinear interpolation is of second-order in this scenario is highly inaccurate, and leads to convergence problems when viscosity varies greatly between dicretized points. Use of the Galerkin coarse-grid operator in combination with an operator-dependent inter-grid transfer routine also results in a proper upscaling of the effective viscosity on the coarse-grid. This combined routine will converge in a number of operations that is entirely independent of the type of viscosity variations employed on the fine grid, i.e. discontinuous variations in viscosity of many orders of magnitude will converge at exactly the same rate as isoviscous distributions.

This code will eventually become an improved three-dimensional, spherical analogue of STAG3D, which is Paul Tackley's well-known finite-difference staggered grid convection code for Cartesian geometry. We are calling this code "CSSTAG" for the time being, unless we can find a better acronym. This is an important extension of our current codes since many of the problems our group is interested in, such as self-consistent plate-like behavior or thermo-chemical stability of dense piles at the base of the lower mantle, depend very strongly on geometry and will exhibit different behavior in a sphere that cannot be captured in a box. I am currently in the process of refining and verifying the code for variable viscosity scenarios, and will be conducting benchmark calculations with numerous other codes in the near future.

Click here for a movie of my rotational advection test. A steep 3-D Gaussian bell centered on each block is advected 2 pi around the sphere. You'll notice that MPDATA does a wonderful job of handling numerical dissipation as well as distortion due to non-orthogonality of the coordinates. After one full rotation you can notice a very slight shrinking of the equatorial isosurfaces, reflecting the remaining numerical diffusion due to higher order terms. Note that the polar gaussian shape has undergone pure rotation, with no apparent numerical dissipation or distortion. This reflects the fact that numerical diffusion is roughly proportional to the number of grid cells that the advected structure passes through. Because a much greater degree of thermal diffusion also acts on scalar fields advected using MPDATA (i.e. a low Peclet number), this technique is more than sufficiently accurate for treatment of the temperature field.

References:

• Hernlund, J.; Tackley, P.J.; Three-Dimensional Spherical Shell Convection at Infinite Prandtl Number Using the 'Cubed Sphere' Method. Book of Conference Proceedings, Second M.I.T. Conference on Computational Fluid and Solid Mechanics, 2003 (in press). Click here for a reprint (PDF, 393K).