Matrix-free, robust linear solvers for realistic, large scale, high-contrast forward mantle convection simulations

In order to accurately capture the physics of mantle convection, nonlinear, compressible flow models and high resolution numerical simulations are required. The viscosity inside the Earth's mantle exhibits strong gradients and in general depends on the pressure, temperature and strain rate, leading to numerically challenging high-contrast problems. Since the parameters governing geophysical models for mantle convection are often fraught with uncertainties, a major goal is to solve inverse problems (usually facilitated via the adjoint method), which require computationally expensive feedback loops. These requirements put a significant emphasis on the development of robust, iterative solvers for the associated linear systems and a highly scalable, parallel implementation usable on distributed memory HPC systems. Additionally, the prohibitive memory requirements of matrix assembly necessitate a matrix-free approach in order to make high resolution models feasible.

In this talk, we employ the finite element framework HyTeG that has successfully been applied to problems with 10^12 unknowns and consider two compressible forward mantle convection models, respectively with the truncated anelastic liquid approximation (TALA) and projected density approximation (PDA). The models feature high-contrast temperature-dependent viscosities together with realistic plate motion reconstructions as boundary conditions at Earth's surface. For the PDA model, we also consider the take-up and release of the latent heat that occurs at mineral phase boundaries by means of effective thermal expansivity and effective specific heat capacity stored on high-resolution lookup tables precomputed with a thermodynamics framework.

For the TALA model, as presented in [Burkhart et al. 2026], we consider the quasi-stationary Stokes system coupled with a time dependent advection diffusion equation, solved in a Gauss-Seidel like alternating fashion as part of a semi implicit scheme. To handle the advection dominated energy conservation equation we use a novel operator splitting approach, combining the BDF2 method with a particle method, resulting in an overall second order time discretisation. Due to handling the compressibility term implicitly (in contrast to a frozen velocity approach) the Stokes system admits an asymmetrical generalised saddle point problem, which we solve by the application of an outer FMGRES solver loop combined with an Uzawa type block precondtioner. As part of the block preconditioner we use geometric multigrid approaches and a new kind of BFBT type Schur complement approximation.

Due to the high contrast in thermal expansivity and specific heat capacity along mineral phase transitions, we discuss adaptations of this solving strategy for the PDA model and aim for a comparison with the entropy formulation introduced in [Dannberg et al. 2022].