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.