Abstract: Iterative Solvers for Modeling Mantle Convection with Strongly Varying Viscosity
Iterative Solvers for Modeling Mantle Convection with Strongly Varying Viscosity
Christoph Köstler
Abstract
This dissertation presents significant improvements to the spherical finite-element discretization and to the iterative solver of the Stokes equations within the high-performance mantle convection code Terra.
For this purpose, a stabilized Q1-Q1 finite-element discretization of the Stokes equations in a two-dimensional square domain has been studied in terms of evaluating its spectral properties depending on grid spacing, viscosity model and viscosity contrast. It could be shown that the spectrum of the Schur complement S becomes independent of the grid spacing when the stabilization as proposed by Dohrmann and Bochev (2004) is applied. To get this spectrum also independent of the viscosity contrast, S has to be scaled by the diagonal of a viscosity-weighted pressure mass matrix which is also spectral equivalent to the stabilization matrix.
The above-mentioned finite-element discretization has been extended to a Stokes solver study tool (SSST). which has been used to compare three Krylov subspace methods: Pressure Correction, MINRES and a CG method using a block-triangular preconditioner, proposed by Bramble and Pasciak (1988). Except MINRES, all solvers have been transformed to a restarted version, using an inner-outer scheme with inner and outer stopping criteria derived from eigenvalue estimates. In the comparison, emphasis was on performance and on robustness with respect to viscosity and to iteration parameter choices. The study revealed that the difference between the Krylov solvers was less than a factor of two. However, the pressure correction algorithm showed slightly the best performance while being the simplest method to implement.
To improve the spherical finite-element discretization in Terra, the stabilization matrix C has been included as in SSST. It is ready to use on grids with at least 84.5 millions of nodes, on coarser grids it must be weighted as the projection of the pressure to a piecewise constant function leads to a higher maximum divergence error. An adaptive weighting of C has been implemented into the solver of Terra.
From the findings of the two-dimensional study, the pressure correction algorithm of Terra has been refined and prior to solving, also S is scaled as in SSST, which leads to iteration numbers considerably less dependent on the viscosity variation than before. By applying the variable-viscosity mass matrix scaling, the total number of multigrid iterations in the first ten time steps of a convection simulation could be reduced by a factor of four in the presence of strong lateral viscosity variations. This improvement could be even larger, up to a factor of 20, if the convergence of the multigrid solver, which is used for calculating velocities and velocity search directions, would not depend that much on the viscosity variations. Volume-weighted harmonic viscosity averaging has been introduced to Terra to apply cellwise constant viscosities. These improvements allow to model the convection of Earth's mantle more realistically.
Citation: Christoph Köstler. Iterative Solvers for Modeling Mantle Convection with Strongly Varying Viscosity, Ph.D. Thesis, Friedrich-Schiller-Univ. Jena, 2011.