Newton solver stabilisation for Stokes solvers in geodynamic problems

The most commonly used method by the geodynamical community for solving non-linear equations is the Picard
fixed-point iteration. However, the Newton method has recently gained interest within this community because
it formally leads to quadratic convergence close to the solution as compared to the global linear convergence of
the Picard iteration. In mantle dynamics, a blend of pressure and strain-rate dependent visco-plastic rheologies is
often used. While for power-law rheologies the Jacobian is guaranteed to be Symmetric Positive Definite (SPD),
for more complex (compressible) rheologies, the Jacobian may become non-SPD. Here we present a new method
for efficiently enforce the Jacobian to be SPD, necessary for our current highly efficient Stokes solvers, with
a minimum loss in convergence rate. Furthermore, we show results for both incompressible and compressible
models.