Los Alamos National Laboratory

Cray Research, Inc. and Los Alamos National Laboratory

Abstract

An efficient multigrid solver is the heart of a 3-D spherical finite element
code named TERRA developed for studying the mantle dynamics of the terrestrial
planets. The multigrid solver exploits a natural hierarchy of nested, almost
uniform, triangular grids constructed from the projection of the regular
icosahedron onto the sphere. The 3-D finite elements are hexagonal prisms
centered on the nodes of the mesh. Basis functions are Cartesian products of
piecewise linear spherical barycentric functions for the tangential directions
and normal piecewise linear functions for the radial direction. The dynamical
equations are written in terms of the primitive variables of velocity,
temperature, and pressure. In the infinite Prandtl number limit, appropriate
for creeping flow of silicate rock in planetary mantles, the momentum equation
reduces to a balance among buoyancy forces, viscous forces, and the pressure
gradient. Given the temperature field to each time step, a conjugate gradient
scheme built around the multigrid solver is used to find the velocity and
pressure simultaneously under the anelastic constraint that
\nabla \cdot (\rho_r **u** = 0, where \rho_r is a radially symmetric
reference density and **u** is velocity. The time step is adjusted
dynamically such that only one or two V-cycles are required per step.
Eulerian advection is utilized in advancing the temperature equation.

Applied to the Earth, TERRA shows how subduction of ocean lithosphere along the margins of the supercontinent Pangea naturally leads to Pangea's breakup and to the migration of the fragments toward their present locations. Global simulations with resolution on the order of 50 km are now practical on parallel machines.