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 EARTH SCIENCES, UTRECHT UNIVERSITY, NETHERLANDS, (2) DEPARTMENT OF MATHEMATICS, COLORADO STATE UNIVERSITY, U.S.A,
(3) GFZ GERMAN RESEARCH CENTER FOR GEOSCIENCES, POTSDAM, GERMANY, (4) CENTRE OF EARTH EVOLUTION AND DYNAMICS (CEED), UNIVERSITY OF OSLO, NORWAY
ABSTRACT
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. We found that the Jacobian is not always Symmetric Positive 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
7. Our Newton solver also works for com- pressible models
PROBLEM STATEMENT
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:
Jkuu Jkup Jkpu 0
δUk δ ¯Pk
= Fku Fkp
where for the incompressible case the Jacobian elements are (ignoring the pressure scaling):
Jkuu
ij = (Ak)ij +
ε(ϕui ), 2
∂η(ε(uk), pk)
∂ε : ε(ϕuj )
ε(uk)
, Jkup
ij = BijT +
ε(ϕui ), 2 ∂η(ε(uk), pk)
∂p ϕpj
ε(uk)
, Jkpu
ij = Bij.
This is in general neither Symmetric, nor Posi- tive Definite, which can be very bad for solvers.
BENCHMARKING THE NEWTON METHOD
RESTORING SYMMETRY FOR PRECONDITIONING
Approximate Jacobian with Jkup ≈ BT and
Jkuu
ij ≈ (Ak)ij +
ε(ϕui ), ∂η(ε(uk), pk)
∂ε : ε(ϕuj )
ε(uk)
+
ε(ϕuj ), ∂η(ε(uk), pk)
∂ε : ε(ϕui )
ε(uk)
= (Ak)ij +
ε(ϕui ), Esym(ε(uk))ε(ϕuj )
where Esym = 12 (Emnpq + Epqmn) and E(ε(uk))mnpq = h2ε(u)mn ∂η(ε(u),p)∂εpq i.
DEMONSTRATING THE WORLD GENERATOR
The next step is to apply the Newton solver to complex 3D geodynamic settings, such as constructed by our World Generator (Fraters et al., in prep), which generates inititial conditions with little input:
Available features: different plate types, faults and variations in curvature and thickness within slabs.
RESTORING POSITIVE DEFINITENESS
The positive definiteness of the Jacobian is determined by tensor H:
Juu =
ε(ϕui ), 2η(ε(u))ε(ϕuj )
+
ε(ϕui ), Esym(ε(u))ε(ϕuj )
= ε(ϕui ) 2η(ε(u))I ⊗ I + Esym(ε(u))
| {z }
=:H
We can force tensor H to be SPD by scaling Esym with a factor α:
Hspd = 2η(ε(u))I ⊗ I + αEsym(ε(u))
with 0 < α ≤ 1. If α = 0 the Jacobian 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 − kakkbkb:a i2
kakkbk < 2η(ε(u))
2η(ε(u))
1− b:a kakkbk
2
kakkbk
otherwise.
