On long enough time scales, the Earth mantle (the region between the rigid plates at the surface and the liquid metal outer core at depth) behaves like a non-Newtonian fluid with a strongly temperature dependent viscosity. While it moves only a few centimeters per year, the large length scales nevertheless lead to very large Rayleigh numbers and, consequently, very complex and expensive numerical simulations. At the same time, given the inaccessibility of the Earth mantle to direct experimental observation implies that numerical simulation is one of the few available tools to elucidate what exactly is going on the mantle, how it affects the long-term evolution of Earth's thermal and chemical structure, as well as what drives and sustains plate motion.

I will here review the approach we have taken in building the state-of-the-art open source solver ASPECT (see http://aspect.dealii.org) to simulate realistic conditions in the Earth and other celestial bodies. ASPECT is built using some of the most widely used and best software libraries for common tasks, such as deal.II for mesh handling and discretization, p4est for parallel partitioning and rebalancing, and Trilinos for linear algebra. In this talk, I will focus on the choices we have made regarding the numerical methods used in ASPECT, and in particular on the interplay between higher order discretizations on adaptive meshes, linear and nonlinear solvers, optimal preconditioners, and approaches to scale to thousands of processor cores. All of these are necessary for simulations that can answer geosphysical questions.