Efficient and practical Newton solvers for non-linear Stokes systems in geodynamic problems

Many problems in geodynamic modelling result in a non-linear Stokes problem in which the viscosity depends on the strain rate and pressure (in addition to other variables). After discretization, the resulting non-linear system is most commonly solved using a Picard fixedpoint iteration. However, it is well understood that Newton's method - when augmented by globalization strategies to ensure convergence even from points far from the solution - can be substantially more efficient and accurate than a Picard solver. In this contribution, we evaluate how a straightforward Newton method must be modified to allow for the kinds of rheologies common in geodynamics. Specifically, we show that the Newton step is not actually well posed for strain rate-weakening models without modifications to the Newton matrix. We derive modifications that guarantee well posedness and that also allow for efficient solution strategies by ensuring that the top left block of the Newton matrix is symmetric and positive definite.We demonstrate the applicability and relevance of thesemodifications with a sequence of benchmarks and a test case of realistic complexity.