We have implemented a Newton solver for ASPECT (Kronbichler et al., 2012) 2

N ^EWTON S ^OLVER S TABILIZATION FOR S ^TOKES S ^{OLVERS IN} G ^EODYNAMIC P^ROBLEMS

Fraters M. R. T. (1), Bangerth W. (2), Thieulot C. (1), Glerum A. (1,3), Spakman W. (1,4)

(1) DEPARTMENT OF E^ARTH S^CIENCES, U^TRECHT U^NIVERSITY, N^ETHERLANDS, (2) DEPARTMENT OF M^ATHEMATICS, C^OLORADO S^TATE U^NIVERSITY, U.S.A, (3) GFZ G^ERMAN R^ESEARCH C^ENTER

FOR G^EOSCIENCES, P^OTSDAM, G^ERMANY, (4) C^{ENTRE OF} E^ARTH EVOLUTION AND D^YNAMICS (CEED), UNIVERSITY OF O^SLO, N^ORWAY

A^BSTRACT

1. We have implemented a Newton solver for ASPECT (Kronbichler et al., 2012)

2. We only use solver libraries to solve linear systems (no Trilinos NOX or PETSC SNES) 3. We have implemented line search and

oversolving-prevention ourselves

4. The Jacobian is not always Symmetric Pos- itive Definite (SPD)

5. We force the Jacobian to be SPD in a cheap and optimal way

6. This allows for more complex rheologies to be used with a Newton solver

P^ROBLEM S^TATEMENT

We are interested solving the Stokes equations:

−∇ ·

 2η



˙

(u) − 1

3 (∇ · u)1



+ ∇p = ρg in Ω,

∇ · (ρu) = 0 in Ω,

The weak form of the Newton linearisation:

J_k^uu J_k^up J_k^pu 0

 δU_k δ ¯P_k



= F_k^u F_k^p



where for the incompressible case the Jacobian elements are (ignoring the pressure scaling):

J_k^uu

ij = (A_k)_ij +



ε(ϕ^u_i ), 2 ∂η(ε(u_k), p_k)

∂ε : ε(ϕ^u_j )

ε(u_k)

 , J_k^up

ij = B_ij^T +



ε(ϕ^u_i ), 2 ∂η(ε(u_k), p_k)

∂p ϕ^p_j 

ε(u_k)

 , J_k^pu

ij = B_ij.

In general J^uu neither symmetric, nor positive definite, which can be very bad for solvers.

R^ESTORING S^YMMETRY

Approx. J^uu with the eq. below where ^{Esym = 1}₂ ^{(Emnpq + Epqmn)} and ^{E(ε(uk))mnpq = h2ε(u)mn ∂η(ε(u),p)}_∂εpq ^{i :}

J_k^uu

ij ≈ (A_k)_ij +



ε(ϕ^u_i ),  ∂η(ε(u_k), p_k)

∂ε : ε(ϕ^u_j )

ε(u_k)

 +



ε(ϕ^u_j ),  ∂η(ε(u_k), p_k)

∂ε : ε(ϕ^u_i )

ε(u_k)



= (A_k)_ij + 

ε(ϕ^u_i ), E^sym(ε(u_k))ε(ϕ^u_j )

F^ORCING SYMMETRY AND P^OSITIVE DEFINITENESS

The positive definiteness of J^uu is determined by tensor H:

J^uu =



ε(ϕ^u_i ), 2η(ε(u))ε(ϕ^u_j )

 +



ε(ϕ^u_i ), E^sym(ε(u))ε(ϕ^u_j )



= ε(ϕ^u_i ), 2η(ε(u))I ⊗ I + E^sym(ε(u))

| {z }

=:H

ε(ϕ^u_j )

We can force tensor H to be SPD by scaling E^sym with a factor α:

H^spd = 2η(ε(u))I ⊗ I + αE^sym(ε(u))

with 0 < α ≤ 1. If ^{α = 0} J^uu will always be SPD, but to have optimal convergence we want ^α to be as large as possible (max one). It can be proven (Fraters et al., in prep) that the optimal value for ^α is (where ^{a = ε(u)}and ^{b = ∂η(ε(u),p)}_∂ε ):

α =







1 if h

1 − _kakkbk^b:a i2

kakkbk < 2η(ε(u))

2η(ε(u))



1− b:a kakkbk

2

kakkbk

otherwise.

R^EALISTIC 3^D C^ASE

P^RELIMINARY C^ONCLUSIONS

The Newton solver without stabilisation is:

• fast;

• prone to numerical breakdowns;

• very sensitive to the precise tweaking of parameters such as mini- mum linear tolerance, line search iterations, etc.

The Newton solver with stabilisation is:

• faster than Picard, slower than without stabilisation;

• almost immune to numerical breakdowns;

• very insensitive to tweaking of these parameters.

BENCHMARKING

We show here results of a setup based on Spiegelman et al. (2016), on a mesh that has been refined either 4 (64x16 cells) or 8 times

(1024x256 cells) uniformly.