Manual
User Manual:
Open the PDF directly: View PDF 
.
Page Count: 54
| Download | |
| Open PDF In Browser | View PDF |
The Finite Element Method in Geodynamics
C. Thieulot
December 17, 2018
Contents
1 Introduction
1.1 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Essential literature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.3 Installation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
3
3
3
2 The physical equations of Fluid Dynamics
2.1 the Boussinesq approximation: an Incompressible flow . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
5
3 The
3.1
3.2
3.3
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
6
6
6
6
7
7
4 Additional techniques
4.1 The method of manufactured solutions . . . . .
4.2 Sparse storage . . . . . . . . . . . . . . . . . . .
4.3 Mesh generation . . . . . . . . . . . . . . . . .
4.4 The value of the timestep . . . . . . . . . . . .
4.5 Tracking materials . . . . . . . . . . . . . . . .
4.6 Visco-Plasticity . . . . . . . . . . . . . . . . . .
4.7 Picard and Newton . . . . . . . . . . . . . . . .
4.8 The choice of solvers . . . . . . . . . . . . . . .
4.9 The SUPG formulation for the energy equation
4.10 Tracking materials and/or interfaces . . . . . .
4.11 Dealing with a free surface . . . . . . . . . . . .
4.12 Pressure normalisation . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
11
11
11
11
11
11
11
11
11
11
11
11
11
3.4
Finite Element Method
Numerical integration . . . .
The mesh . . . . . . . . . . .
Elements and basis functions
3.3.1 The Q1 space . . . . .
The penalty approach . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5 fieldstone: simple analytical solution
13
6 fieldstone: Stokes sphere
16
7 fieldstone: Convection in a 2D box
17
8 fieldstone: solcx benchmark
19
9 fieldstone: solkz benchmark
21
10 fieldstone: solvi benchmark
22
11 fieldstone: the indentor benchmark
24
12 fieldstone: the annulus benchmark
26
13 fieldstone: stokes sphere (3D) - penalty
28
14 fieldstone: stokes sphere (3D) - mixed formulation
29
15 fieldstone: consistent pressure recovery
30
1
16 fieldstone: the Particle in Cell technique (1) - the effect of averaging
31
17 fieldstone: solving the full saddle point problem
34
18 fieldstone: solving the full saddle point problem in 3D
36
19 fieldstone: The non-conforming Q1 × P0 element
39
20 fieldstone: The stabilised Q1 × Q1 element
40
21 fieldstone: compressible flow (1)
41
22 fieldstone: compressible flow (2)
43
A The main codes in computational
A.1 ADELI . . . . . . . . . . . . . . .
A.2 ASPECT . . . . . . . . . . . . .
A.3 CITCOMS and CITCOMCU . .
A.4 DOUAR . . . . . . . . . . . . . .
A.5 GAIA . . . . . . . . . . . . . . .
A.6 GALE . . . . . . . . . . . . . . .
A.7 GTECTON . . . . . . . . . . . .
A.8 ELVIS . . . . . . . . . . . . . . .
A.9 ELEFANT . . . . . . . . . . . . .
A.10 ELLIPSIS . . . . . . . . . . . . .
A.11 FANTOM . . . . . . . . . . . . .
A.12 FLUIDITY . . . . . . . . . . . .
A.13 LAMEM . . . . . . . . . . . . . .
A.14 MILAMIN . . . . . . . . . . . . .
A.15 PARAVOZ/FLAMAR . . . . . .
A.16 PTATIN . . . . . . . . . . . . . .
A.17 RHEA . . . . . . . . . . . . . . .
A.18 SEPRAN . . . . . . . . . . . . .
A.19 SOPALE . . . . . . . . . . . . .
A.20 STAGYY . . . . . . . . . . . . .
A.21 SULEC . . . . . . . . . . . . . .
A.22 TERRA . . . . . . . . . . . . . .
A.23 UNDERWORLD 1&2 . . . . . .
geodynamics
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
45
45
45
45
45
45
45
45
45
45
45
45
45
45
45
45
45
45
45
45
45
45
45
45
B fieldstone.py
46
C fieldstone stokes sphere.py
47
D fieldstone convection box.py
47
E fieldstone solcx.py
48
F fieldstone indentor.py
49
G fieldstone saddlepoint.py
50
2
1
Introduction
practical hands-on approach
as little as possible jargon
no mathematical proof
no optimised codes (readability over efficiency). avoiding as much as possible to have to look elsewhere. very
sequential, so unavoidable repetitions (jacobian, shape functions)
FE is one of several methods.
All the python scripts and this document are freely available at https://github.com/cedrict/fieldstone.
1.1
Acknowledgments
Jean Braun, Philippe Fullsack, Arie van den Berg. Lukas van de Wiel. Robert Myhill. Menno, Anne Too many BSc
and MSc students to name indivisually, although Job Mos did produce the very first version of fieldstone as part of
his MSc thesis. The ASPECT team in general and Wolfgang Bangerth and Timo Hesiter in particular.
1.2
Essential literature
a)
1.3
b)
c)
d)
Installation
python3.6 -m pip install --user numpy scipy matplotlib
3
e)
2
The physical equations of Fluid Dynamics
Symbol
t
x, y, z
v
ρ
η
λ
T
∇
∇·
p
ε̇(v)
α
k
Cp
H
βT
meaning
Time
Cartesian coordinates
velocity vector
mass density
dynamic viscosity
penalty parameter
temperature
gradient operator
divergence operator
pressure
strain rate tensor
thermal expansion coefficient
thermal conductivity
Heat capacity
intrinsic specific heat production
isothermal compressibility
unit
s
m
m· s−1
kg/m3
Pa· s
Pa· s
K
m−1
m−1
Pa
s−1
K−1
W/(m · K)
J/K
W/kg
Pa−1
Let us start from the heat transport equation as shown in Schubert, Turcotte and Olson [37]:
ρCp
DT
Dp
− αT
= ∇ · k∇T + Φ + ρH
Dt
Dt
with
∂T
Dp
∂p
DT
=
+ v · ∇T
=
+ v · ∇p
Dt
∂t
Dt
∂t
In order to arrive at the set of equations that ASPECT solves, we need to neglect the ∂p/∂t. WHY? Also, their
definition of the shear heating term Φ is:
Φ = kB (∇ · v)2 + 2η ε̇d : ε̇d
For many fluids the bulk viscosity kB is very small and is often taken to be zero, an assumption known as the Stokes
assumption: kB = λ + 2η/3 = 0. Note that η is the dynamic viscosity and λ the second viscosity. Also,
τ = 2η ε̇ + λ(∇ · v)1
but since kB = λ + 2η/3 = 0, then λ = −2η/3 so
2
τ = 2η ε̇ − η(∇ · v)1 = 2η ε̇d
3
[from aspect manual] We focus on the system of equations in a d = 2- or d = 3-dimensional domain Ω that describes
the motion of a highly viscous fluid driven by differences in the gravitational force due to a density that depends on
the temperature. In the following, we largely follow the exposition of this material in Schubert, Turcotte and Olson
[37].
Specifically, we consider the following set of equations for velocity u, pressure p and temperature T :
1
in Ω,
(1)
−∇ · 2η ε̇(v) − (∇ · v)1 + ∇p = ρg
3
∇ · (ρv) = 0
ρCp
∂T
+ v · ∇T
∂t
in Ω,
(2)
− ∇ · k∇T = ρH
1
+ 2η ε̇(v) − (∇ · v)1
3
1
: ε̇(v) − (∇ · v)1
3
+ αT (v · ∇p)
∂X
+ ρT ∆S
+ v · ∇X
∂t
where ε̇(u) = 21 (∇u + ∇uT ) is the symmetric gradient of the velocity (often called the strain rate).
4
(3)
in Ω,
In this set of equations, (67) and (68) represent the compressible Stokes equations in which v = v(x, t) is the
velocity field and p = p(x, t) the pressure field. Both fields depend on space x and time t. Fluid flow is driven by the
gravity force that acts on the fluid and that is proportional to both the density of the fluid and the strength of the
gravitational pull.
Coupled to this Stokes system is equation (69) for the temperature field T = T (x, t) that contains heat conduction
terms as well as advection with the flow velocity v. The right hand side terms of this equation correspond to
• internal heat production for example due to radioactive decay;
• friction heating;
• adiabatic compression of material;
• phase change.
The last term of the temperature equation corresponds to the latent heat generated or consumed in the process of
phase change of material. In what follows we will not assume that no phase change takes place so that we disregard
this term altogether.
2.1
the Boussinesq approximation: an Incompressible flow
[from aspect manual] The Boussinesq approximation assumes that the density can be considered constant in all
occurrences in the equations with the exception of the buoyancy term on the right hand side of (67). The primary
result of this assumption is that the continuity equation (68) will now read
∇·v =0
This implies that the strain rate tensor is deviatoric. Under the Boussinesq approximation, the equations are much
simplified:
−∇ · [2η ε̇(v)] + ∇p = ρg
ρ0 Cp
∂T
+ v · ∇T
∂t
in Ω,
(4)
∇ · (ρv) = 0
in Ω,
(5)
− ∇ · k∇T = ρH
in Ω
(6)
Note that all terms on the rhs of the temperature equations have disappeared, with the exception of the source term.
5
3
The Finite Element Method
3.1
Numerical integration
3.2
The mesh
3.3
Elements and basis functions
Let us for a moment consider a single quadrilateral element in the xy-plane, as shown on the following figure:
Let us assume that we know the values of a given field u at the vertices. For a given point M inside the element in
the plane, what is the value of the field u at this point? It makes sense to postulate that uM = u(xM , yM ) will be
given by
uM = φ(u1 , u2 , u3 , u4 , xM , yM )
where φ is a function to be determined. Although φ is not unique, we can decide to express the value uM as a weighed
sum of the values at the vertices ui . One option could be to assign all four vertices the same weight, say 1/4 so that
uM = (u1 + u2 + u3 + u4 )/4, i.e. uM is simply given by the arithmetic mean of the vertices values. This approach
suffers from a major drawback as it does not use the location of point M inside the element. For instance, when
(xM , yM ) → (x2 , y2 ) we expect uM → u2 .
In light of this, we could now assume that the weights would depend on the position of M in a continuous fashion:
u(xM , yM ) =
4
X
Ni (xM , yM ) ui
i=1
where the Ni are continous (”well behaved”) functions which have the property:
Ni (xj , yj ) = δij
or, in other words:
N3 (x1 , y1 )
=
0
(7)
N3 (x2 , y2 )
=
0
(8)
N3 (x3 , y3 )
=
1
(9)
N3 (x4 , y4 )
=
0
(10)
The functions Ni are commonly called basis functions.
Omitting the M subscripts for any point inside the element, the velocity components u and v are given by:
u(x, y) =
4
X
Ni (x, y) ui
i=1
v(x, y) =
4
X
Ni (x, y) vi
i=1
Rather interestingly, one can now easily compute velocity gradients (and therefore the strain rate tensor) since we
6
have assumed the basis functions to be ”well behaved” (in this case differentiable):
4
˙xx (x, y)
=
∂u X ∂Ni
=
ui
∂x
∂x
i=1
˙yy (x, y)
=
∂v X ∂Ni
=
vi
∂y
∂y
i=1
˙xy (x, y)
=
1 ∂u 1 ∂v
1 X ∂Ni
1 X ∂Ni
+
=
ui +
vi
2 ∂y
2 ∂x
2 i=1 ∂y
2 i=1 ∂x
(11)
4
(12)
4
4
(13)
How we actually obtain the exact form of the basis functions is explained in the coming section.
3.3.1
The Q1 space
N1 (r, s)
=
0.25(1 − r)(1 − s)
N2 (r, s)
=
0.25(1 + r)(1 − s)
N3 (r, s)
=
0.25(1 + r)(1 + s)
N4 (r, s)
=
0.25(1 − r)(1 + s)
∂N1
(r, s)
∂r
∂N2
(r, s)
∂r
∂N3
(r, s)
∂r
∂N4
(r, s)
∂r
∂N1
(r, s)
∂s
∂N2
(r, s)
∂s
∂N3
(r, s)
∂s
∂N4
(r, s)
∂s
3.4
=
−0.25(1 − s)
=
+0.25(1 − s)
=
+0.25(1 + s)
=
−0.25(1 + s)
=
−0.25(1 − r)
=
−0.25(1 + r)
=
+0.25(1 + r)
=
+0.25(1 − r)
The penalty approach
In order to impose the incompressibility constraint, two widely used procedures are available, namely the Lagrange
multiplier method and the penalty method [1, 20]. The latter is implemented in elefant, which allows for the
elimination of the pressure variable from the momentum equation (resulting in a reduction of the matrix size).
Mathematical details on the origin and validity of the penalty approach applied to the Stokes problem can for
instance be found in [8], [33] or [19].
The penalty formulation of the mass conservation equation is based on a relaxation of the incompressibility constraint and writes
p
∇·v+ =0
(14)
λ
7
where λ is the penalty parameter, that can be interpreted (and has the same dimension) as a bulk viscosity. It is
equivalent to say that the material is weakly compressible. It can be shown that if one chooses λ to be a sufficiently
large number, the continuity equation ∇ · v = 0 will be approximately satisfied in the finite element solution. The
value of λ is often recommended to be 6 to 7 orders of magnitude larger than the shear viscosity [14, 21].
Equation (14) can be used to eliminate the pressure in Eq. (??) so that the mass and momentum conservation
equations fuse to become :
∇ · (2η ε̇(v)) + λ∇(∇ · v) = ρg = 0
(15)
[27] have established the equivalence for incompressible problems between the reduced integration of the penalty
term and a mixed Finite Element approach if the pressure nodes coincide with the integration points of the reduced
rule.
In the end, the elimination of the pressure unknown in the Stokes equations replaces the original saddle-point
Stokes problem [2] by an elliptical problem, which leads to a symmetric positive definite (SPD) FEM matrix. This is
the major benefit of the penalized approach over the full indefinite solver with the velocity-pressure variables. Indeed,
the SPD character of the matrix lends itself to efficient solving stragegies and is less memory-demanding since it is
sufficient to store only the upper half of the matrix including the diagonal [18] . ToDo: list codes which use this
approach.
8
The stress tensor σ is symmetric (i.e. σij = σji ). For simplicity I will now focus on a Stokes flow in two dimensions.
Since the penalty formulation is only valid for incompressible flows, then ˙ = ˙ d so that the d superscript is
ommitted in what follows. The stress tensor can also be cast in vector format:
σxx
−p
˙xx
σyy = −p + 2η ˙yy
σxy
0
˙xy
˙xx + ˙yy
˙xx
= λ ˙xx + ˙yy + 2η ˙yy
0
˙xy
∂u
∂x
2 0 0
1 1 0
∂v
1 1 0 +η 0 2 0 ·
= λ
∂y
0 0 1
0 0 0
{z
}
{z
}
|
|
∂u
+ ∂v
K
Remember that
∂y
C
4
∂x
4
∂u X ∂Ni
=
ui
∂x
∂x
i=1
∂v X ∂Ni
=
vi
∂y
∂y
i=1
4
4
X
X
∂u ∂v
∂Ni
∂Ni
+
=
ui +
vi
∂y
∂x
∂y
∂x
i=1
i=1
so that
∂u
∂x
∂v
∂y
∂u
∂y
+
∂v
∂x
=
∂N1
∂x
0
∂N2
∂x
0
∂N3
∂x
0
∂N4
∂x
0
0
∂N1
∂y
0
∂N2
∂y
0
∂N3
∂y
0
∂N4
∂y
∂N1
∂y
∂N1
∂x
∂N2
∂y
∂N2
∂x
∂N3
∂y
∂N3
∂x
∂N3
∂y
∂N4
∂x
|
{z
B
·
}
|
u1
v1
u2
v2
u3
v3
u4
v4
{z
V
}
Finally,
σxx
~σ = σyy = (λK + ηC) · B · V
σxy
We will now establish the weak form of the momentum conservation equation. We start again from
∇·σ+b=0
For the Ni ’s ’regular enough’, we can write:
Z
Z
Ni ∇ · σdΩ +
Ωe
Ni bdΩ = 0
Ωe
We can integrate by parts and drop the surface term1 :
Z
Z
∇Ni · σdΩ =
Ωe
Ni bdΩ
Ωe
or,
Z
∂Ni
∂x
0
∂Ni
∂y
0
∂Ni
∂y
∂Ni
∂x
Ωe
1 We
Z
σxx
· σyy dΩ =
Ni bdΩ
Ωe
σxy
will come back to this at a later stage
9
Let i = 1, 2, 3, 4 and stack the resulting four equations on top of one another.
Z
∂N1
∂x
0
∂N1
∂y
∂N2
∂x
0
Ωe
Z
0
∂N2
∂y
∂N3
∂x
0
Ωe
Z
∂N1
∂x
∂N2
∂y
∂N2
∂x
∂N3
∂y
0
∂N3
∂y
∂N4
∂x
0
∂N3
∂x
∂N4
∂y
0
∂N4
∂y
∂N4
∂x
Ωe
Z
∂N1
∂y
0
Ωe
σxx
· σyy dΩ
σxy
σxx
· σyy dΩ
σxy
σxx
· σyy dΩ
σxy
σxx
· σyy dΩ
σxy
Z
=
N1
Ωe
Z
=
Ni
Ωe
Z
=
N3
Ωe
Z
=
N4
Ωe
bx
by
bx
by
bx
by
bx
by
dΩ
(16)
dΩ
(17)
dΩ
(18)
dΩ
(19)
We easily recognize B T inside the integrals! Let us define
NbT = (N1 bx , N1 by , ...N4 bx , N4 by )
then we can write
Z
σxx
B T · σyy dΩ =
Nb dΩ
Ωe
Ωe
σxy
Z
and finally:
Z
Z
T
B · [λK + ηC] · B · V dΩ =
Ωe
Nb dΩ
Ωe
Since V contains the velocities at the corners, it does not depend on the x or y coordinates so it can be taking outside
of the integral:
Z
Z
B T · [λK + ηC] · BdΩ · |{z}
V =
Nb dΩ
Ωe
Ωe
{z
} (8x1) | {z }
|
Ael (8×8)
Bel (8×1)
or,
Z
Z
Z
T
T
λB · K · BdΩ +
V =
Nb dΩ
ηB · C · BdΩ · |{z}
Ω
Ω
| Ωe
{z
} | e
{z
} (8x1) | e {z }
Aλ
el (8×8)
Aη
el (8×8)
INTEGRATION - MAPPING
1. partition domain Ω into elements Ωe , e = 1, ...nel .
2. loop over elements and for each element compute Ael , Bel
3. a node belongs to several elements
→ need to assemble Ael and Bel in A, B
4. apply boundary conditions
5. solve system: x = A−1 · B
6. visualise/analyse x
10
Bel (8×1)
4
Additional techniques
4.1
The method of manufactured solutions
4.2
Sparse storage
4.3
Mesh generation
4.4
The value of the timestep
4.5
Tracking materials
4.6
Visco-Plasticity
4.7
Picard and Newton
4.8
The choice of solvers
4.9
The SUPG formulation for the energy equation
4.10
Tracking materials and/or interfaces
4.11
Dealing with a free surface
4.12
Pressure normalisation
11
so much to do ...
write about impose bc on el matrix
Q1Q1-stab
Q2Q1
Q3Q2
full compressible
total energy calculations
constraints
compositions, marker chain
van keken initial value with deformed mesh
free-slip bc on annulus and sphere . See for example p540 Gresho and Sani book.
non-linear rheologies (punch, two layer brick spmw16, tosn15)
Picard vs Newton
markers
Schur complement approach
periodic boundary conditions
open boundary conditions
export to vtu
free surface
SUPG
produce fastest version possible for convection box
zaleski disk advection
all kinds of interesting benchmarks
Busse convection pb, compare with aspect
cvi !!!
pure elastic
including phase changes (w. R. Myhill)
derivatives on nodes
Nusselt
discontinuous galerkin
formatting of code style
navier-stokes ? (LUKAS)
pressure smoothing
compute strainrate in middle of element or at quad point for punch?
GEO1442 code
GEO1442 indenter setup in plane ?
in/out flow on sides for lith modelling
Problems to solve:
colorscale
better yet simple matrix storage ?
12
5
fieldstone: simple analytical solution
From [12]. In order to illustrate the behavior of selected mixed finite elements in the solution of stationary Stokes
flow, we consider a two-dimensional problem in the square domain Ω = [0, 1] × [0, 1], which possesses a closed-form
analytical solution. The problem consists of determining the velocity field v = (u, v) and the pressure p such that
−ν∆v + ∇p = b
∇·v =0
v=0
in Ω
in Ω
on Γ
where the fluid viscosity is taken as ν = 1. The components of the body force b are prescribed as
bx
(12 − 24y)x4 + (−24 + 48y)x3 + (−48y + 72y 2 − 48y 3 + 12)x2
=
+(−2 + 24y − 72y 2 + 48y 3 )x + 1 − 4y + 12y 2 − 8y 3
by
(8 − 48y + 48y 2 )x3 + (−12 + 72y − 72y 2 )x2
=
+(4 − 24y + 48y 2 − 48y 3 + 24y 4 )x − 12y 2 + 24y 3 − 12y 4
With this prescribed body force, the exact solution is
Note that the pressure obeys
features
R
Ω
u(x, y)
= x2 (1 − x)2 (2y − 6y 2 + 4y 3 )
v(x, y)
= −y 2 (1 − y)2 (2x − 6x2 + 4x3 )
p(x, y)
= x(1 − x) − 1/6
p dΩ = 0
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (no-slip)
• direct solver
• isothermal
• isoviscous
• analytical solution
13
vx
1.0
vy
1.0
0.010
0.8
0.8
0.005
0.4
0.4
0.2
0.6
0.8
0.01
x
0.6
0.8
1.0
yy
0.8
0.6
0.0000075
0.6
0.8
1.0
0.2
0.4
x
0.6
0.8
1.0
p pth
1.0
0.00010
0.0000050 0.8
0.00005
0.0000025
0.00000
0.6
0.00005
y
y
y
x
0.6
0.04
0.4
0.0000025
0.00010
0.2
0.00000500.2
0.00020
0.0
0.0
0.0000075
0.0
0.0
0.000010
0.4
0.04
0.4
0.000005
0.2
xy
0.0000000
0.000000
0.4
0.15
0.00
0.0
0.0
vy tyth
0.000005
1.0
0.02
1.0
0.8
0.8
0.2
0.03
0.6
0.6
0.4
0.02
x
x
0.6
0.01
0.4
0.4
0.02
0.01
0.2
0.2
0.8
0.02
0.000010
0.2
0.10
1.0
0.03
1.0
0.8
0.0
0.0
0.2
0.0
0.0
1.0
0.00
0.0
0.0
vx txth
1.0
0.8
0.2
0.03
0.4
0.6
0.4
0.02
0.2
x
0.6
0.01
0.2
0.4
y
y
0.00
0.4
0.2
0.8
0.02
0.6
0.010
1.0
0.03
0.8
0.0
0.0
0.0
0.0
1.0
xx
1.0
0.05
y
0.010
x
0.4
0.005
0.005
0.2
0.4
0.00
0.6
0.000
y
y
0.000
y
0.6
0.2
0.05
0.8
0.005
0.6
0.0
0.0
p
1.0
0.010
0.2
0.4
x
0.6
0.8
1.0
0.100000
0.00015
0.00025
0.2
0.4
x
0.6
0.8
1.0
velocity
pressure
x2
x1
0.010000
error
0.001000
0.000100
0.000010
0.000001
0.01
0.1
h
Quadratic convergence for velocity error, linear convergence for pressure error, as expected.
ToDo:
pressure normalisation?
different cmat, a la schmalholz
To go further:
14
1. make your own analytical solution
15
6
fieldstone: Stokes sphere
Viscosity and density directly computed at the quadrature points.
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (free-slip)
• direct solver
• isothermal
• non-isoviscous
• buoyancy-driven flow
vx
1.0
0.8
0.8
0.0005
0.6
0.4
0.0005
0.0010
0.2
0.4
x
0.6
0.8
1.0
xx
1.0
0.8
0.2
0.6
0.2
0.0
0.0
0.2
0.2
0.002
0.2
0.4
x
0.6
0.8
0.4
0.0
0.0
1.0
yy
0.2
0.010
0.0
0.0
1.0
y
0.0050
0.2
0.010
0.2
0.4
x
0.6
0.8
density
0.4
x
0.6
0.8
0.005
0.010
0.015
0.2
0.4
x
0.6
0.8
1.0
viscosity
2.0
1.0
0.8
1.8
0.8
80
0.6
1.6
0.6
60
0.4
1.4
0.4
40
0.2
1.2
0.2
20
1.0
0.0
0.0
0.0025
0.2
0.000
0.0
0.0
1.0
y
0.0125
0.0075
0.005
100
y
0.0150
0.4
0.010
0.2
1.0
0.0175
0.0100
0.015
0.005
0.2
0.6
xy
0.4
0.005
0.8
1.0
y
0.000
0.4
x
0.8
0.6
y
y
0.4
0.8
x
0.6
0.8
0.6
0.6
0.4
0.005
0.6
0.4
0.2
1.0
0.010
0.005
0.2
0.0
0.4
0.8
0.000
p
0.4
0.001
1.0
0.010
0.8
0.0
0.0
0.001
0.000
0.4
0.2
1.0
y
y
0.0000
0.002
y
0.6
0.0
0.0
vy
1.0
0.0010
0.0
0.0
0.2
0.4
x
0.6
0.8
16
1.0
0.2
0.4
x
0.6
0.8
1.0
7
fieldstone: Convection in a 2D box
This benchmark deals with the 2-D thermal convection of a fluid of infinite Prandtl number in a rectangular closed
cell. In what follows, I carry out the case 1a, 1b, and 1c experiments as shown in [3]: steady convection with constant
viscosity in a square box.
The temperature is fixed to zero on top and to ∆T at the bottom, with reflecting symmetry at the sidewalls (i.e.
∂x T = 0) and there are no internal heat sources. Free-slip conditions are implemented on all boundaries.
The Rayleigh number is given by
αgy ∆T h3 ρ2 cp
αgy ∆T h3
=
(20)
Ra =
κν
kµ
In what follows, I use the following parameter values: Lx = Ly = 1,ρ0 = cP = k = µ = 1, T0 = 0, α = 10−2 ,
g = 102 Ra and I run the model with Ra = 104 , 105 and 106 .
The initial temperature field is given by
T (x, y) = (1 − y) − 0.01 cos(πx) sin(πz)
(21)
The perturbation in the initial temperature fields leads to a perturbation of the density field and sets the fluid in
motion.
Depending on the initial Rayleigh number, the system ultimately reaches a steady state after some time.
The Nusselt number (i.e. the mean surface temperature gradient over mean bottom temperature) is computed as
follows [3]:
R ∂T
∂y (y = Ly )dx
(22)
N u = Ly R
T (y = 0)dx
Note that in our case the denominator is equal to 1 since Lx = 1 and the temperature at the bottom is prescribed to
be 1.
Finally, the steady state root mean square velocity and Nusselt number measurements are indicated in Table ??
alongside those of [3] and [?]. (Note that this benchmark was also carried out and published in other publications
[?, ?, ?, ?, ?] but since they did not provide a complete set of measurement values, they are not included in the table.)
4
Ra = 10
Ra = 105
Ra = 106
Vrms
Nu
Vrms
Nu
Vrms
Nu
Blankenbach et al
42.864947 ± 0.000020
4.884409 ± 0.000010
193.21454 ± 0.00010
10.534095 ± 0.000010
833.98977 ± 0.00020
21.972465 ± 0.000020
Tackley [?]
42.775
4.878
193.11
10.531
833.55
21.998
Steady state Nusselt number N u and Vrms measurements as reported in the literature.
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (free-slip)
• Boussinesq approximation
• direct solver
• non-isothermal
• buoyancy-driven flow
• isoviscous
• CFL-condition
17
0.4
20
0.2
0.2
0.4
x
0.6
0.8
40 0.0
0.0
1.0
0.4
100
0.8
x
0.6
0.8
1.0
yy
1.0
0
100
0.4
0.2
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
0.2
200 0.00.0
0.25
0.8
100
0.2
0.4
x
0.6
0.8
1.0
qy
1.0
1.00 0.4
3
0.2
1.25 0.2
2
y
0.4
x
0.6
0.8
1.0
1.50 0.0
0.0
1.75
1
0.2
0.4
x
0.6
0.8
1.0
0
0.6
0.4
0.4
y
0.6
0.0
400000 0.0
1.0
40
x
0.6
0.8
1.0
0.0
rho
0.998
0.996
0.6
0.4
0
0.994
10 0.2
0.2
0.4
x
0.6
0.8
1.0
20 0.0
0.0
30
80
7
70
6
0.2
0.4
x
0.6
0.8
1.0
0.992
5
50
40
30
4
3
20
2
10
0
0.4
0.8
20
10
0.2
0.2
1.0
30
60
5
4
0.4
0.8
xy
0.0
0.0
6
0.8
0.75 0.6
0.2
x
0.6
0.8
7
0.6
0.0
0.0
0.4
0.2
y
0.50
0.2
0.4
0
0.00
qx
1.0
0.0
0.0
0.6
0.4
100
0.8
200000
0.2
0.2
y
0.6
y
y
60
0.4
200 1.0
0.8
0.6
0
y
40
0.2
60
xx
1.0
0.6
20
0.2
T
2000000.8
Nu
0.4
0
y
y
0
0.8
20
0.6
400000
1.0
p
1.0
40
0.8
20
0.6
0.0
0.0
1.0
40
0.8
vrms
1.0
1.0
60
vy
y
60
vx
1
0.00 0.05 0.10 0.15 0.20 0.25 0.30
time
ToDo:
implement steady state criterion
reach steady state
do Ra=1e4, 1e5, 1e6
plot against blankenbach paper and aspect
look at critical Ra number
18
0.00 0.05 0.10 0.15 0.20 0.25 0.30
time
8
fieldstone: solcx benchmark
The SolCx benchmark is intended to test the accuracy of the solution to a problem that has a large jump in the
viscosity along a line through the domain. Such situations are common in geophysics: for example, the viscosity in a
cold, subducting slab is much larger than in the surrounding, relatively hot mantle material.
The SolCx benchmark computes the Stokes flow field of a fluid driven by spatial density variations, subject to a
spatially variable viscosity. Specifically, the domain is Ω = [0, 1]2 , gravity is g = (0, −1)T and the density is given by
ρ(x, y) = sin(πy) cos(πx)
(23)
Boundary conditions are free slip on all of the sides of the domain and the temperature plays no role in this benchmark.
The viscosity is prescribed as follows:
1
f or x < 0.5
µ(x, y) =
(24)
106 f or x > 0.5
Note the strongly discontinuous viscosity field yields a stagnant flow in the right half of the domain and thereby yields
a pressure discontinuity along the interface.
The SolCx benchmark was previously used in [15] (references to earlier uses of the benchmark are available there)
and its analytic solution is given in [43]. It has been carried out in [25] and [17]. Note that the source code which
evaluates the velocity and pressure fields for both SolCx and SolKz is distributed as part of the open source package
Underworld ([29], http://underworldproject.org).
In this particular example, the viscosity is computed analytically at the quadrature points (i.e. tracers are not
used to attribute a viscosity to the element). If the number of elements is even in any direction, all elements (and
their associated quadrature points) have a constant viscosity(1 or 106 ). If it is odd, then the elements situated at the
viscosity jump have half their integration points with µ = 1 and half with µ = 106 (which is a pathological case since
the used quadrature rule inside elements cannot represent accurately such a jump).
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (free-slip)
• direct solver
• isothermal
• non-isoviscous
• analytical solution
19
1.0
0.0010
0.8
0.002
0.8
0.6
0.001
0.6
0.6
y
0.0000
0.4
y
0.0005
0.0005
0.0010
0.2
0.0
0.0
0.0015
0.0
0.0
0.4
x
0.6
0.8
1.0
xx
1.0
0.0075
0.8
0.0025
y
0.0000
0.4
0.0
0.0
0.0050
0.0075
0.2
0.4
x
0.6
0.8
1.0
vx txth
1.0
0.4
x
0.6
0.8
0.0100
yy
0.0050
0.0075
0.4
x
0.6
0.8
1.0
vy tyth
0.6
0.8
1.0
0.005
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
p pth
0.010
0.020
0.4
0.00000200.4
0.005
0.0000025
0.2
0.0000030
0.0
0.0
0.010
0.010
0.000
y
0.0
0.0
0.6
0.000
0.2
0.005
0.000003 0.0
0.0
x
0.0100
0.005
0.4
0.6
0.0000015
0.000002 0.2
0.4
0.010
0.6
0.2
0.2
0.020
0.015
y
y
0.000001
0.025
0.0000005
0.8
0.0000010
0.000000
0.4
1.0
xy
0.0000000 1.0
0.000002 0.8
0.000001
0.8
0.6
0.0025
0.2
x
0.6
0.015
0.0025
0.2
0.4
0.8
0.0050
0.4
0.2
1.0
0.0100
0.0000
0.0
0.0
0.2
0.0
0.0
0.0075
0.000003 1.0
0.8
0.1
1.0
0.6
0.0025
0.2
0.2
y
0.6
0.0
0.4
0.002
0.8
0.0050
0.1
0.2
1.0
0.0100
0.2
0.001
0.2
0.2
0.000
0.4
p
1.0
0.003
y
0.8
vy
0.0015
y
vx
1.0
0.2
0.4
x
0.6
0.8
What we learn from this
20
1.0
0.015
0.2
0.4
x
0.6
0.8
1.0
0.020
9
fieldstone: solkz benchmark
The SolKz benchmark [34] is similar to the SolCx benchmark. but the viscosity is now a function of the space
coordinates:
µ(y) = exp(By) with B = 13.8155
(25)
It is however not a discontinuous function but grows exponentially with the vertical coordinate so that its overall
variation is again 106 . The forcing is again chosen by imposing a spatially variable density variation as follows:
ρ(x, y) = sin(2y) cos(3πx)
(26)
Free slip boundary conditions are imposed on all sides of the domain. This benchmark is presented in [43] as well and
is studied in [15] and [17].
0.2
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
y
0.6
0.4
0.2
0.2
0.4
x
0.6
0.8
vx txth
1.0
0.4
0.6
0.8
0.0
0.0
0.2
101
x
0.6
0.8
1.0
0.8
0.4
0.2
x
error
y
0.2
0.4
x
0.6
0.8
1.0
0.0008
0.0006
0.0004
0.0002
0.0000
0.0002
0.0004
0.0006
0.0008
0.6
0.8
10
4
10
5
0.4
10
6
0.2
10
7
0.0
0.0
1.0
0.00010000
0.00001000
0.00000100
0.00000010
0.01
p pth
0.6
0.00100000
0.00000001
1.0
0.8
velocity
pressure
x2
x1
0.01000000
0.8
0.0004
0.0003
0.0002
0.0001
0.0000
0.0001
0.0002
0.0003
0.0004
1.0
0.6
0.10000000
0.0
0.0
0.75
0.50
0.25
0.00
0.25
0.50
0.75
0.4
0.6
0.2
density
0.2
x
0.4
1.0
0.8
0.0
0.0
0.4
0.6
1.0
0.6
0.2
0.8
0.5
x
1.0
xy
1.0
0.5
0.4
0.8
0.2
1e 7
0.2
0.6
0.4
0.0
0.0
0.0
x
0.6
1.0
0.2
0.4
0.8
1.0
y
y
102
0.4
0.8
0.4
0.0
0.0
103
0.2
0.6
0.6
104
0.0
0.0
x
vy tyth
1.0
0.2
0.4
0.2
1.0
0.00075
0.00050
0.00025
0.00000
0.00025
0.00050
0.00075
0.8
105
0.4
1.0
yy
1.0
1.0
0.6
0.8
0.2
1.0
0.8
0.6
0.4
2
4
6
0.2
x
0.0
0.0
x
y
y
0.6
0.4
0.00006
0.6
1e 7
0.8
0.2
0.2
0.0
0.0
0.8
1.0
6
4
2
0
0.0
0.0
0.00004
0.4
0.075
0.050
0.025
0.000
0.025
0.050
0.075
0.4
0.00002
0.2
p
0.6
0.2
y
0.8
0.0
0.0
0.00000
0.4
1.0
0.00075
0.00050
0.00025
0.00000
0.00025
0.00050
0.00075
0.8
0.00002
0.6
xx
1.0
1.0
0.00004
y
0.4
0.8
0.00006
y
y
0.6
vy
1.0
y
0.8
0.000100
0.000075
0.000050
0.000025
0.000000
0.000025
0.000050
0.000075
0.000100
y
vx
1.0
0.1
h
21
0.2
0.4
x
0.6
0.8
1.0
10
fieldstone: solvi benchmark
Following SolCx and SolKz, the SolVi inclusion benchmark solves a problem with a discontinuous viscosity field, but
in this case the viscosity field is chosen in such a way that the discontinuity is along a circle. Given the regular nature
of the grid used by a majority of codes and the present one, this ensures that the discontinuity in the viscosity never
aligns to cell boundaries. This in turns leads to almost discontinuous pressures along the interface which are difficult to
represent accurately. [36] derived a simple analytic solution for the pressure and velocity fields for a circular inclusion
under simple shear and it was used in [10], [38], [15], [25] and [17].
Because of the symmetry of the problem, we only have to solve over the top right quarter of the domain (see Fig.
??a).
The analytical solution requires a strain rate boundary condition (e.g., pure shear) to be applied far away from the
inclusion. In order to avoid using very large domains and/or dealing with this type of boundary condition altogether,
the analytical solution is evaluated and imposed on the boundaries of the domain. By doing so, the truncation error
introduced while discretizing the strain rate boundary condition is removed.
A characteristic of the analytic solution is that the pressure is zero inside the inclusion, while outside it follows the
relation
µm (µi − µm ) ri2
cos(2θ)
(27)
pm = 4˙
µi + µm r2
where µi = 103 is the viscosity of the inclusion and µm = 1 is the viscosity of the background media, θ = tan−1 (y/x),
and ˙ = 1 is the applied strain rate.
[10] thoroughly investigated this problem with various numerical methods (FEM, FDM), with and without tracers,
and conclusively showed how various averagings lead to different results. [15] obtained a first order convergence for
both pressure and velocity, while [25] and [17] showed that the use of adaptive mesh refinement in respectively the
FEM and FDM yields convergence rates which depend on refinement strategies.
22
vx
1.0
vy
1.0
0.0
0.2
0.8
0.6
0.6
0.6
0.4
0.6
0.4
0.4
0.4
0.6
0.4
0.2
0.2
0.2
0.8
0.2
0.0
0.0
0.0
0.0
0.0
0.8
1.0
xx
1.0
0.8
1.75
1.0
1.50
0.8
1.25
0.6
0.25
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
0.00
vx txth
1.0
0.000
0.4
x
0.6
1.50
0.2
0.0
0.0
1.75
0.0
0.0
0.8
0.8
1.0
vy tyth
103
0.6
0.8
0.8
102
y
0.6
y
0.6
0.4
0.4
101
0.2
0.2
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
100
0.0
0.0
10.00000000
0.2
0.4
x
0.6
0.8
1.0
velocity
pressure
x2
x1
1.00000000
0.10000000
error
0.01000000
0.00100000
0.00010000
0.00001000
0.00000100
0.00000010
0.00000001
0.01
0.0
0.0
0.1
h
23
1.100
1.075
1.050
1.025
1.000
0.975
0.950
0.925
0.900
0.4
x
0.6
0.8
1.0
p pth
1.00
0.75
0.50
0.25
0.00
0.25
0.50
0.75
1.00
15
10
5
0
5
0.2
1.0
density
1.0
0.8
x
0.2
15
1.0
0.4
0.000
0.4
0.8
0.6
0.005
0.2
0.6
0.8
0.015
0.010
x
xy
1.0
0.020
0.0
0.0
1.0
1.0
0.6
y
y
0.2
0.2
x
0.4
0.4
1.25
0.4
0.2
0.6
0.75
0.2
10
0.8
0.50
0.2
0.020
0.0
0.0
0.25
1.00
5
1.0
0.00
0.4
0.015
0.2
0.0
0.0
1.0
0.6
0.010
0.4
0.8
0.8
0.005
0.6
0.6
yy
1.0
0.8
x
0.4
0.50
0.2
0.4
y
y
0.75
0.4
0.2
0.6
1.00
0
y
0.6
5
y
x
10
10
0.2
0.4
x
0.6
0.8
1.0
1.0
15
100
0.8
0.6
10
1
10
2
10
3
y
0.4
15
y
0.8
y
0.8
y
0.8
0.2
p
1.0
0.4
0.2
0.0
0.0
0.2
0.4
x
0.6
0.8
1.0
11
fieldstone: the indentor benchmark
0.5
0.4
0.3
0.2
0.1
0.0
0.0
0.4
x
0.6
0.8
1.0
1.0
10.0
7.5
5.0
2.5
0.0
2.5
5.0
1.0
20.0
17.5
15.0
12.5
10.0
7.5
5.0
2.5
xx
0.2
0.4
0.2
0.4
x
x
0.6
0.6
0.8
0.8
0.5
0.4
0.3
0.2
0.1
0.0
0.0
0.2
vy
0.0
0.2
0.4
y
y
0.5
0.4
0.3
0.2
0.1
0.0
0.0
0.6
0.2
0.4
x
0.6
0.8
1.0
0.4
x
0.5
0.4
0.3
0.2
0.1
0.0
0.0
0.6
0.8
1.0
5.0
2.5
0.0
2.5
5.0
7.5
10.0
0.5
0.4
0.3
0.2
0.1
0.0
0.0
6
p
5
4
3
2
0.2
0.4
1.0
yy
0.2
0.8
y
0.2
0.3
0.2
0.1
0.0
0.1
0.2
0.3
x
0.6
0.8
1.0
xy
0.2
0.4
x
0.6
0.8
1.0
103
0.5
0.4
0.3
0.2
0.1
0.0
0.0
102
101
100
0.2
0.4
x
0.6
0.8
1.0
y
0.5
0.4
0.3
0.2
0.1
0.0
0.0
vx
y
0.5
0.4
0.3
0.2
0.1
0.0
0.0
y
y
y
y
The punch benchmark is one of the few boundary value problems involving plastic solids for which there exists an
exact solution. Such solutions are usually either for highly simplified geometries (spherical or axial symmetry, for
instance) or simplified material models (such as rigid plastic solids) [24].
In this experiment, a rigid punch indents a rigid plastic half space; the slip line field theory gives exact solutions
as shown in Fig. ??a. The plane strain formulation of the equations and the detailed solution to the problem were
derived in the Appendix of [41] and are also presented in [16].
The two dimensional punch problem has been extensively studied numerically for the past 40 years [46, 45, 7, 6, 22,
42, 4, 32] and has been used to draw a parallel with the tectonics of eastern China in the context of the India-Eurasia
collision [39, 28]. It is also worth noting that it has been carried out in one form or another in series of analogue
modelling articles concerning the same region, with a rigid indenter colliding with a rheologically stratified lithosphere
[30, 9, 23].
Numerically, the one-time step punch experiment is performed on a two-dimensional domain of purely plastic
von Mises material. Given that the von Mises rheology yield criterion does not depend on pressure, the density of
the material and/or the gravity vector is set to zero. Sides are set to free slip boundary conditions, the bottom to
no slip, while a vertical velocity (0, −vp ) is prescribed at the top boundary for nodes whose x coordinate is within
[Lx /2 − δ/2, Lx /2 + δ/2].
The following parameters are used: Lx = 1, Ly = 0.5, µmin = 10−3 , µmax = 103 , vp = 1, δ = 0.123456789 and the
yield value of the material is set to k = 1.
The analytical solution predicts that the angle of the shear bands stemming from the sides of the punch is π/4,
that√the pressure right under the punch is 1 + π, and that the velocity of the rigid blocks on each side of the punch is
vp / 2 (this is simply explained by invoking conservation of mass).
10
1
10
2
10
3
0.5
0.4
0.3
0.2
0.1
0.0
0.0
1
20
15
10
5
0
5
10
15
20
5
p (nodal)
4
3
2
0.2
0.4
x
0.6
0.8
1.0
1
6
0.0
5
0.2
4
0.00
p
0.2
0.05
v
u
0.4
0.10
0.4
3
0.6
0.05
2
0.8
0.10
1
1.0
0.0
0.2
0.4
x
0.6
0.8
1.0
0.0
0.2
0.4
x
0.6
24
0.8
1.0
0.0
0.2
0.4
x
0.6
0.8
1.0
ToDo: smooth punch
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (no-slip)
• isothermal
• non-isoviscous
• nonlinear rheology
25
12
fieldstone: the annulus benchmark
This benchmark is based on Thieulot & Puckett [Subm.] in which an analytical solution to the isoviscous incompressible
Stokes equations is derived in an annulus geometry. The velocity and pressure fields are as follows:
vr (r, θ)
=
g(r)k sin(kθ),
(28)
vθ (r, θ)
=
f (r) cos(kθ),
(29)
p(r, θ)
=
kh(r) sin(kθ),
(30)
ρ(r, θ)
=
ℵ(r)k sin(kθ),
(31)
with
f (r)
g(r)
h(r)
ℵ(r)
A
B
= Ar + B/r,
B
C
A
r + ln r + ,
=
2
r
r
2g(r) − f (r)
=
,
r
g0
f
f0
g
= g 00 − − 2 (k 2 − 1) + 2 + ,
r
r
r
r
2(ln R1 − ln R2 )
,
= −C 2
R2 ln R1 − R12 ln R2
R22 − R12
= −C 2
.
R2 ln R1 − R12 ln R2
(32)
(33)
(34)
(35)
(36)
(37)
The parameters A and B are chosen so that vr (R1 ) = vr (R2 ) = 0, i.e. the velocity is tangential to both inner
and outer surfaces. The gravity vector is radial and of unit length. In the present case, we set R1 = 1, R2 = 2 and
C = −1.
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions
• direct solver
• isothermal
• isoviscous
• analytical solution
• annulus geometry
• elemental boundary conditions
26
10
velocity
pressure
x2
x1
1
error
0.1
0.01
0.001
0.0001
0.01
h
27
13
fieldstone: stokes sphere (3D) - penalty
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (free-slip)
• direct solver
• isothermal
• non-isoviscous
• 3D
• elemental b.c.
• buoyancy driven
resolution is 24x24x24
28
14
fieldstone: stokes sphere (3D) - mixed formulation
This is the same setup as Section 14.
features
• Q1 × P0 element
• incompressible flow
• mixed formulation
• Dirichlet boundary conditions (free-slip)
• direct solver
• isothermal
• non-isoviscous
• 3D
• elemental b.c.
• buoyancy driven
29
15
fieldstone: consistent pressure recovery
p = −λ∇ · v
q1 is smoothed pressure obtained with the center-to-node approach.
q2 is recovered pressure obtained with [44]. R
All three fulfill the zero average condition: pdΩ = 0.
y
0.000
0.4
0.2
0.4
x
0.6
0.8
vx txth
1.0
0.000002
x
0.6
0.8
0.8
0.4
0.2
0.10
0.0
0.0
0.15
x
0.6
0.8
1.0
0.6
y
0.6
0.4
0.2
0.2
0.4
x
0.6
0.8
1.0
0.10
0.2
0.15
0.2
0.4
x
0.6
0.8
1.0
q2 pth
1.0
0.00002
0.00000
0.00002
0.00004
0.00006
0.2
0.4
0.004
0.8
0.003
1.0
0.6
0.4
0.2
0.0
0.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.6
0.8
1.0
1.0
0.8
0.002
0.6
0.4
0.001
0.4
0.2
0.000
0.2
0.2
0.4
x
0.6
0.8
1.0
0.001
0.0
0.0
velocity
pressure (el)
pressure (q1)
pressure (q2)
x2
x1
0.010000
error
0.001000
0.000100
0.000010
0.01
0.8
0.8
0.6
0.0
0.0
0.100000
0.000001
x
0.6
y
0.8
0.007
0.006
0.005
0.004
0.003
0.002
0.001
0.000
0.05
0.4
q1 pth
1.0
0.00
0.0
0.0
0.15
1.0
0.05
y
y
0.05
0.4
q2
1.0
p pth
0.0
0.0
1.0
0.8
0.00
0.6
x
0.8
0.2
0.000002
0.8
0.6
0.4
0.000001
0.6
x
0.6
0.000000
0.4
0.4
0.8
0.000001
0.2
0.2
1.0
0.000002
1.0
0.05
0.2
vy tyth
0.0
0.0
1.0
0.0
0.0
1.0
0.2
q1
1.0
0.8
0.4
0.000002
0.4
0.6
0.6
0.000001
0.2
x
y
y
0.000000
0.4
0.4
0.8
0.000001
0.6
0.2
0.2
0.10
0.2
0.010
1.0
0.8
0.0
0.0
0.0
0.0
1.0
0.05
0.4
0.005
0.2
0.010
0.00
0.6
0.000
0.4
0.005
0.2
0.0
0.0
0.6
0.05
0.8
0.005
y
0.6
0.0
0.0
0.8
0.005
p
1.0
0.010
y
0.8
vy
1.0
0.010
y
vx
1.0
0.1
h
In terms of pressure error, q2 is better than q1 which is better than elemental.
QUESTION: why are the averages exactly zero ?!
TODO:
• add randomness to internal node positions.
• look at elefant algorithms
30
16
fieldstone: the Particle in Cell technique (1) - the effect of averaging
features
• Q1 × P0 element
• incompressible flow
• penalty formulation
• Dirichlet boundary conditions (no-slip)
• isothermal
• non-isoviscous
• particle-in-cell
After the initial setup of the grid, markers can then be generated and placed in the domain. One could simply
randomly generate the marker positions in the whole domain but unless a very large number of markers is used, the
chance that an element does not contain any marker exists and this will prove problematic. In order to get a better
control over the markers spatial distribution, one usually generates the marker per element, so that the total number
of markers in the domain is the product of the number of elements times the user-chosen initial number of markers
per element.
Our next concern is how to actually place the markers inside an element. Two methods come to mind: on a regular
grid, or in a random manner, as shown on the following figure:
In both cases we make use of the basis shape functions: we generate the positions of the markers (random or
regular) in the reference element first (rim , sim ), and then map those out to the real element as follows:
xim =
m
X
Ni (rim , sim )xi
yim =
i
m
X
Ni (rim , sim )yi
i
where xi , yi are the coordinates of the vertices of the element.
When using active markers, one is faced with the problem of transferring the properties they carry to the mesh on
which the PDEs are to be solved. As we have seen, building the FE matrix involves a loop over all elements, so one
simple approach consists of assigning each element a single property computed as the average of the values carried by
the markers in that element. Often in colloquial language ”average” refers to the arithmetic mean:
n
hφiam =
1X
φi
n
k
where < φ >am is the arithmetic average of the n numbers φi . However, in mathematics other means are commonly
used, such as the geometric mean:
!
n
Y
hφigm =
φi
i
and the harmonic mean:
n
hφihm =
1X 1
n i φi
!−1
Furthermore, there is a well known inequality for any set of positive numbers,
hφiam
≥
hφigm
31
≥
hφihm
which will prove to be important later on.
Let us now turn to a simple concrete example: the 2D Stokes sphere. There are two materials in the domain, so
that markers carry the label ”mat=1” or ”mat=2”. For each element an average density and viscosity need to be
computed. The majority of elements contains markers with a single material label so that the choice of averaging does
not matter (it is trivial to verify that if φi = φ0 then hφiam = hφigm = hφihm = φ0 . Remain the elements crossed by
the interface between the two materials: they contain markers of both materials and the average density and viscosity
inside those depends on 1) the total number of markers inside the element, 2) the ratio of markers 1 to markers 2, 3)
the type of averaging.
This averaging problem has been studied and documented in the literature [35, 10, 40, 31]
Nodal projection. Left: all markers inside elements to which the green node belongs to are taken into account.
Right: only the markers closest to the green node count.
The setup is identical as the Stokes sphere experiment. The bash script runall runs the code for many resolutions,
both initial marker distribution and all three averaging types. The viscosity of the sphere has been set to 103 while
the viscosity of the surrounding fluid is 1. The average density is always computed with an arithmetic mean. Root
mean square velocity results are shown hereunder:
random markers, elemental proj
regular markers, elemental proj
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
0.009
0.008
0.007
vrms
0.006
0.005
0.004
0.003
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
0.01
0.009
0.008
0.007
0.006
vrms
0.01
0.005
0.004
0.003
0.002
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
0.002
0.001
0.001
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0
0.01
0.02
0.03
0.04
h
0.05
0.06
0.07
0.08
0.09
0.1
h
random markers, nodal proj
regular markers, nodal proj
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
0.009
0.008
0.007
vrms
0.006
0.005
0.004
0.003
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
0.01
0.009
0.008
0.007
0.006
vrms
0.01
0.005
0.004
0.003
0.002
0.002
0.001
0.001
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0
0.01
0.02
0.03
0.04
h
0.05
0.06
0.07
0.08
0.09
0.1
h
random markers, nodal proj
regular markers, nodal proj
am,
am,
am,
am,
am,
am,
gm,
gm,
gm,
gm,
gm,
gm,
hm,
hm,
hm,
hm,
hm,
hm,
0.009
0.008
0.007
vrms
0.006
0.005
0.004
0.003
09
16
25
36
48
64
09
16
25
36
48
64
09
16
25
36
48
64
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
m
0.01
0.009
0.008
0.007
0.006
vrms
0.01
0.005
0.004
0.003
0.002
0.002
0.001
0.001
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
0
h
0.01
0.02
0.03
0.04
0.05
h
32
0.06
0.07
0.08
0.09
0.1
Conclusions:
• With increasing resolution (h → 0) vrms values seem to converge towards a single value, irrespective of the
number of markers.
• At low resolution, say 32x32 (i.e. h=0.03125), vrms values for the three averagings differ by about 10%. At
higher resolution, say 128x128, vrms values are still not converged.
• The number of markers per element plays a role at low resolution, but less and less with increasing resolution.
• Results for random and regular marker distributions are not identical but follow a similar trend and seem to
converge to the same value.
• At low resolutions, elemental values yield better results.
• harmonic mean yields overal the best results
64 markers per element
0.01
0.009
am, reg, el
gm, reg, el
hm, reg, el
am, rand, el
gm, rand, el
hm, rand, el
am, reg, nod1
gm, reg, nod1
hm, reg, nod1
am, rand, nod1
gm, rand, nod1
hm, rand, nod1
am, reg, nod2
gm, reg, nod2
hm, reg, nod2
am, rand, nod2
gm, rand, nod2
hm, rand, nod2
0.008
0.007
vrms
0.006
0.005
0.004
0.003
0.002
0.001
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
h
33
0.08
0.09
0.1
17
fieldstone: solving the full saddle point problem
The details of the numerical setup are presented in Section 5.
The main difference is that we no longer use the penalty formulation and therefore keep both velocity and pressure
as unknowns. Therefore we end up having to solve the following system:
K G
V
f
·
=
or,
A · X = rhs
P
h
GT 0
Each block K, G and vector f , h are built separately in the code and assembled into the matrix A and vector rhs
afterwards. A and rhs are then passed to the solver. We will see later that there are alternatives to solve this approach
which do not require to build the full Stokes matrix A.
Each element has m = 4 vertices so in total ndof V × m = 8 velocity dofs and a single pressure dof, commonly
situated in the center of the element. The total number of velocity dofs is therefore N f emV = nnp × ndof V while
the total number of pressure dofs is N f emP = nel. The total number of dofs is then N f em = N f emV + N f emP .
As a consequence, matrix K has size N f emV, N f emV and matrix G has size N f emV, N f emP . Vector f is of size
N f emV and vector h is of size N f emP .
features
• Q1 × P0 element
• incompressible flow
• mixed formulation
• Dirichlet boundary conditions (no-slip)
• direct solver (?)
• isothermal
• isoviscous
• analytical solution
• pressure smoothing
34
0.010
vx
1.0
0.010
vy
1.0
0.005 0.8
0.005 0.8
0.6
0.6
0.6
0.4
0.0
0.0
0.4
0.0050.2
0.2
0.2
0.4
x
0.6
0.8
1.0
0.00
0.05
0.10
y
y
y
0.000
0.05
0.0
0.010 0.0
0.4
0.15
0.0050.2
0.2
0.4
x
0.6
0.8
1.0
0.00
q
1.0
0.05
0.8
0.10
0.6
y
0.8
0.000
p
1.0
0.0
0.010 0.0
0.20
0.2
0.4
x
0.6
0.8
1.0
0.15
0.4
0.20
0.2
0.25 0.00.0
0.2
0.4
x
0.6
0.8
1.0
0.06
0.8
0.01
0.6
0.02
0.8
0.01
0.6
0.00
y
y
0.00
1.0
0.4
0.01
0.2
0.0
0.0
0.02
0.2
0.4
x
0.6
0.8
1.0
0.4
0.01
0.2
0.0
0.0
0.02
0.2
0.0000010
1.0
vx txth
1.0
0.03
1.0
0.4
x
0.6
0.8
1.0
0.03
0.0
0.0
0.0000005
0.2
0.4
x
0.6
0.8
1.0
0.0
0.0
0.0000010
0.2
0.4
x
0.6
0.8
1.0
0.4
x
0.6
0.8
1.0
0.04
0.6
0.03
0.0
0.0
0.0000010
0.4
0.02
0.2
0.04 0.0
0.0
0.2
0.00
p pth
0.02
0.04
0.06
0.08
0.0000005
0.2
0.2
0.05
0.8
0.4
0.06
y
0.6
0.0000000
0.4
y
0.6
0.0000000
0.4
y
0.6
0.2
0.02
0.2
0.8
0.0000005
0.0
0.0
0.00
0.4
0.8
0.0000005
0.2
0.02
0.6
0.8
0.4
0.04 1.0
0.8
0.0000010
1.0
vy tyth
0.06
y
0.02
xy
x
0.6
0.8
1.0
q pth
1.0
0.10
0.2
0.4
x
0.6
0.8
1.0
0.12
0.02
0.8
0.04
0.6
0.06
0.4
0.08
0.2
0.0
0.0
0.01
0.00
y
1.0
0.03
yy
y
0.03
xx
0.25
0.10
0.2
0.4
x
0.6
0.8
1.0
0.12
Unlike the results obtained with the penalty formualtion (see Section 5), the pressure showcases a very strong
checkerboard pattern, similar to the one in [13].
Left: pressure solution as shown in [13]; Right: pressure solution obtained with fieldstone.
Rather interestingly, the nodal pressure (obtained with a simple center-to-node algorithm) fails to recover a correct
pressure at the four corners.
35
18
fieldstone: solving the full saddle point problem in 3D
this does not work because of the Dirichlet b.c. on all 6 sides and all three components It does work well with Q2Q1.
This benchmark begins by postulating a polynomial solution to the 3D Stokes equation [11]:
x + x2 + xy + x3 y
y + xy + y 2 + x2 y 2
v=
(38)
−2z − 3xz − 3yz − 5x2 yz
and
p = xyz + x3 y 3 z − 5/32
(39)
While it is then trivial to verify that this velocity field is divergence-free, the corresponding body force of the Stokes
equation can be computed by inserting this solution into the momentum equation with a given viscosity µ (constant
or position/velocity/strain rate dependent). The domain is a unit cube and velocity boundary conditions simply use
Eq. (38). Following [5], the viscosity is given by the smoothly varying function
µ = exp(1 − β(x(1 − x) + y(1 − y) + z(1 − z)))
(40)
Choosing β = 0 yields a constant velocity µ = e1 (and greatly simplifies the right-hand side). One can easily show
that the ratio of viscosities µ? in the system follows µ? = exp(−3β/4) so that choosing β = 10 yields µ? ' 1808 and
β = 20 yields µ? ' 3.269 × 106 .
We start from the momentum conservation equation:
˙ =f
−∇p + ∇ · (2µ)
The x-component of this equation writes
fx
∂
∂
∂
∂p
+
(2µ˙xx ) +
(2µ˙xy ) +
(2µ˙xz )
∂x ∂x
∂y
∂z
∂p
∂
∂
∂
∂µ
∂µ
∂µ
= −
+ 2µ ˙xx + 2µ ˙xy + 2µ ˙xz + 2 ˙xx + 2 ˙xy + 2 ˙xz
∂x
∂x
∂y
∂z
∂x
∂y
∂z
= −
Let us compute all the block separately:
˙xx
=
1 + 2x + y + 3x2 y
˙yy
=
1 + x + 2y + 2x2 y
˙zz
= −2 − 3x − 3y − 5x2 y
2˙xy
=
(x + x3 ) + (y + 2xy 2 ) = x + y + 2xy 2 + x3
2˙xz
=
(0) + (−3z − 10xyz) = −3z − 10xyz
2˙yz
=
(0) + (−3z − 5x2 z) = −3z − 5x2 z
In passing, one can verify that ˙xx + ˙yy + ˙zz = 0. We further have
∂
2˙xx
∂x
∂
2˙xy
∂y
∂
2˙xz
∂z
∂
2˙xy
∂x
∂
2˙yy
∂y
∂
2˙yz
∂z
∂
2˙xz
∂x
∂
2˙yz
∂y
∂
2˙zz
∂z
=
2(2 + 6xy)
=
1 + 4xy
= −3 − 10xy
=
1 + 2y 2 + 3x2
=
2(2 + 2x2 )
= −3 − 5x2
= −10yz
=
0
=
2(0)
36
(41)
(42)
∂p
∂x
∂p
∂y
∂p
∂z
= yz + 3x2 y 3 z
(43)
= xz + 3x3 y 2 z
(44)
= xy + x3 y 3
(45)
The viscosity is chosen to be
µ(x, y, z) = exp [1 − β(x(1 − x) + y(1 − y) + z(1 − z))]
Constant viscosity Setting β = 0 yields
µ(x, y, z) = e ' 2.718
so that
∂
µ(x, y, z) = 0
∂x
∂
µ(x, y, z) = 0
∂y
∂
µ(x, y, z) = 0
∂z
(46)
(47)
(48)
and
fx
=
=
=
fy
=
=
=
fz
=
=
∂
∂
∂
∂p
+ 2µ ˙xx + 2µ ˙xy + 2µ ˙xz
∂x
∂x
∂y
∂z
−(yz + 3x2 y 3 z) + 2(2 + 6xy) + (1 + 4xy) + (−3 − 10xy)
−
−(yz + 3x2 y 3 z) + µ(2 + 6xy)
∂
∂
∂
∂p
+ 2µ ˙xy + 2µ ˙yy + 2µ ˙yz
−
∂y
∂x
∂y
∂z
−(xz + 3x3 y 2 z) + µ(1 + 2y 2 + 3x2 ) + µ2(2 + 2x2 ) + µ(−3 − 5x2 )
−(xz + 3x3 y 2 z) + µ(2 + 2x2 + 2y 2 )
∂p
∂
∂
∂
−
+ 2µ ˙xz + 2µ ˙yz + 2µ ˙zz
∂z
∂x
∂y
∂z
−(xy + x3 y 3 ) + µ(−10yz) + 0 + 0
3 3
= −(xy + x y ) + µ(−10yz)
(51)
(52)
(53)
(54)
finally
yz + 3x2 y 3 z
2 + 6xy
f = − xz + 3x3 y 2 z + µ 2 + 2x2 + 2y 2
xy + x3 y 3
−10yz
⇒ sign problem with Burstedde et al Eq 26 !!
Variable viscosity
∂
µ(x, y, z)
∂x
∂
µ(x, y, z)
∂y
∂
µ(x, y, z)
∂z
(49)
(50)
=
−(1 − 2x)βµ(x, y, z)
=
−(1 − 2y)βµ(x, y, z)
=
−(1 − 2z)βµ(x, y, z)
37
yz + 3x2 y 3 z
2 + 6xy
= − xz + 3x3 y 2 z + µ 2 + 2x2 + 2y 2
(55)
xy + x3 y 3
−10yz
2˙xx
2˙xy
2˙xz
−(1 − 2x)βµ(x, y, z) 2˙xy − (1 − 2y)βµ(x, y, z) 2˙yy − (1 − 2z)βµ(x, y, z) 2˙yz
2˙xz
2˙yz
2˙zz
2 3
yz + 3x y z
2 + 6xy
= − xz + 3x3 y 2 z + µ 2 + 2x2 + 2y 2
xy + x3 y 3
−10yz
2
2 + 4x + 2y + 6x y
x + y + 2xy 2 + x3
−3z − 10xyz
−3z − 5x2 z
− (1 − 2x)βµ x + y + 2xy 2 + x3 − (1 − 2y)βµ 2 + 2x + 4y + 4x2 y − (1 − 2z)βµ
2
−3z − 10xyz
−3z − 5x z
−4 − 6x − 6y − 10x2 y
f
at (x, y, z) = (0, 0, 0), µ = exp(1),
at (x, y, z) = (0.5, 0.5, 0.5), µ = exp(1 − 3β/4)
Viscosity ratio is given by
exp(1 − 3β/4)
= exp(−3β/4)
µ? =
exp(1)
By varying β between 1 and 22 we can get up to 7 orders of magnitude viscosity difference !
features
• Q1 × P0 element
• incompressible flow
• saddle point system
• Dirichlet boundary conditions (free-slip)
• direct solver
• isothermal
• non-isoviscous
• 3D
• elemental b.c.
• analytical solution
38
fieldstone: The non-conforming Q1 × P0 element
19
features
• Non-conforming Q1 × P0 element
• incompressible flow
• mixed formulation
• isothermal
• non-isoviscous
• analytical solution
• pressure smoothing
try Q1 mapping instead of isoparametric.
0.005 0.8
0.6
0.4
0.6
0.000
0.0
0.0
0.4
0.0050.2
0.2
0.2
0.4
x
0.6
0.8
1.0
0.05
0.00
0.005 0.8
y
y
0.000
p
1.0
0.05
0.6
0.10
y
0.8
0.010
vy
1.0
0.0
0.010 0.0
0.4
0.15
0.0050.2
0.2
0.4
x
0.6
0.8
1.0
0.00
q
1.0
0.05
0.8
0.10
0.6
y
0.010
vx
1.0
0.0
0.010 0.0
0.20
0.2
0.4
x
0.6
0.8
1.0
0.15
0.4
0.20
0.2
0.25 0.00.0
0.2
0.4
x
0.6
0.8
1.0
0.06
0.01
0.6
0.02
0.8
0.01
0.6
0.00
y
y
0.00
1.0
0.4
0.01
0.2
0.0
0.0
0.02
0.2
0.4
x
0.6
0.8
1.0
0.4
0.01
0.2
0.0
0.0
0.02
0.2
0.0000010
1.0
vx txth
1.0
0.03
1.0
0.4
x
0.6
0.8
1.0
0.03
0.0
0.0
0.0000005
0.2
0.4
x
0.6
0.8
1.0
0.0
0.0
0.0000010
0.2
0.4
x
0.6
0.8
1.0
0.4
x
0.6
0.8
1.0
0.04
0.6
0.03
0.0
0.0
0.0000010
39
0.4
0.02
0.2
0.04 0.0
0.0
0.2
0.00
p pth
0.02
0.04
0.06
0.08
0.0000005
0.2
0.2
0.05
0.8
0.4
0.06
y
0.6
0.0000000
0.4
y
0.6
0.0000000
0.4
y
0.6
0.2
0.02
0.2
0.8
0.0000005
0.0
0.0
0.00
0.4
0.8
0.0000005
0.2
0.02
0.6
0.8
0.4
0.04 1.0
0.8
0.0000010
1.0
vy tyth
0.06
y
0.02
0.8
xy
x
0.6
0.8
1.0
q pth
1.0
0.10
0.2
0.4
x
0.6
0.8
1.0
0.12
0.02
0.8
0.04
0.6
0.06
0.4
0.08
0.2
0.0
0.0
0.01
0.00
y
1.0
0.03
yy
y
0.03
xx
0.25
0.10
0.2
0.4
x
0.6
0.8
1.0
0.12
20
fieldstone: The stabilised Q1 × Q1 element
The details of the numerical setup are presented in Section 5.
40
21
fieldstone: compressible flow (1)
We first start with an isothermal Stokes flow, so that we disregard the heat transport equation and the equations we
wish to solve are simply:
1
−∇ · 2η ε̇(v) − (∇ · v)1 + ∇p = ρg
3
∇ · (ρv) = 0
in Ω,
(56)
in Ω
(57)
The second equation can be rewritten ∇ · (ρv) = ρ∇ · v + v · ∇ρ = 0 or,
1
∇ · v + v · ∇ρ = 0
ρ
Note that this presupposes that the density is not zero anywhere in the domain.
We use a mixed formulation and therefore keep both velocity and pressure as unknowns. We end up having to
solve the following system:
K
G
V
f
·
=
or,
A · X = rhs
P
h
GT + Z 0
Where K is the stiffness matrix, G is the discrete gradient operator, GT is the discrete divergence operator, V the
velocity vector, P the pressure vector. Note that the term ZV derives from term v · ∇ρ in the continuity equation.
Each block K, G , Z and vectors f and h are built separately in the code and assembled into the matrix A and
vector rhs afterwards. A and rhs are then passed to the solver. We will see later that there are alternatives to solve
this approach which do not require to build the full Stokes matrix A.
Remark: the term ZV is often put in the rhs (i.e. added to h) so that the matrix A retains the same structure as
in the incompressible case. This is indeed how it is implemented in ASPECT. This however requires more work since
the rhs depends on the solution and some form of iterations is needed.
In the case of a compressible flow the strain rate tensor and the deviatoric strain rate tensor are no more equal
(since ∇ · v 6= 0). The deviatoric strainrate tensor is given by2
1
1
˙ = (v)
˙
˙
− (∇ · v)1
˙ d (v) = (v)
− T r()1
3
3
In that case:
˙dxx
=
˙dyy
=
2˙dxy
=
∂u 1 ∂u ∂v
2 ∂u 1 ∂v
−
+
−
=
∂x 3 ∂x ∂y
3 ∂x 3 ∂y
1 ∂u ∂v
1 ∂u 2 ∂v
∂v
−
+
=−
+
∂y 3 ∂x ∂y
3 ∂x 3 ∂y
∂u ∂v
+
∂y
∂x
and then
˙ d (v) =
2 ∂u
3 ∂x
−
1 ∂v
3 ∂y
1 ∂u
2 ∂y
+
1 ∂v
2 ∂x
1 ∂u
2 ∂y
+
1 ∂v
2 ∂x
− 13 ∂u
∂x +
2 ∂v
3 ∂y
(58)
(59)
(60)
From ~τ = 2η~d we arrive at:
d
˙xx
τxx
2/3
τyy = 2η ˙dyy = 2η −1/3
τxy
0
˙dxy
−1/3
2/3
0
0
0 ·
1/2
∂u
∂x
∂v
∂y
∂u
∂v
∂y + ∂x
4/3
= η −2/3
0
or,
~τ = Cη BV
2 See
the ASPECT manual for a justification of the 3 value in the denominator in 2D and 3D.
41
−2/3
4/3
0
0
0 ·
1
∂u
∂x
∂v
∂y
∂u
∂v
∂y + ∂x
In order to test our implementation we have created a few manufactured solutions:
• benchmark #1 (ibench=1)): Starting from a density profile of:
ρ(x, y) = xy
(61)
We derive a velocity given by:
vx (x, y) =
With gx (x, y) =
1
x
Cx
Cy
, vy (x, y) =
x
y
(62)
and gy (x, y) = y1 , this leads us to a pressure profile:
4Cx
4Cy
+
+ xy + C0
p = −η
3x2
3y 2
(63)
This gives us a strain rate:
−Cx
−Cy
˙yy = 2
˙xy = 0
x2
y
In what follows, we choose η = 1 and Cx = Cy = 1 and for a unit square domain [1 : 2] × [1 : 2] we compute C0
so that the pressure is normalised to zero over the whole domain and obtain C0 = −1.
˙xx =
• benchmark #2 (ibench=2): Starting from a density profile of:
ρ = cos(x) cos(y)
(64)
We derive a velocity given by:
vx =
With gx =
1
cos(y)
and gy =
Cx
Cy
, vy =
cos(x)
cos(y)
(65)
1
cos(x) ,
this leads us to a pressure profile:
!
4Cx sin(x) 4Cy sin(y)
p=η
+
+ (sin(x) + sin(y)) + C0
3 cos2 (x)
3 cos2 (y)
˙xx = Cx
sin(x)
cos2 (x)
˙yy = Cy
sin(y)
cos2 (y)
(66)
˙xy = 0
We choose η = 1 and Cx = Cy = 1. The domain is the unit square [0 : 1] × [0 : 1] and we obtain C0 as before
and obtain
1
C0 = 2 − 2 cos(1) + 8/3(
− 1) ' 3.18823730
cos(1)
(thank you WolframAlpha)
• benchmark #3 (ibench=3)
• benchmark #4 (ibench=4)
• benchmark #5 (ibench=5)
features
• Q1 × P0 element
• incompressible flow
• mixed formulation
• Dirichlet boundary conditions (no-slip)
• isothermal
• isoviscous
• analytical solution
• pressure smoothing
PROBLEMS:
- odd vs even number of elements
- q is ’fine’ everywhere except in the corners - revisit pressure smoothing paper?
42
22
fieldstone: compressible flow (2)
Let us start with some thermodynamics. Every material has an equation of state. The equilibrium thermodynamic
state of any material can be constrained if any two state variables are specified. Examples of state variables include
the pressure p and specific volume ν = 1/ρ, as well as the temperature T .
After linearisation, the density depends on temperature and pressure as follows:
ρ(T, p) = ρ0 ((1 − α(T − T0 ) + βT P )
where α is the coefficient of thermal expansion, also called thermal expansivity:
1 ∂ρ
α=−
ρ ∂T p
α is the percentage increase in volume of a material per degree of temperature increase; the subscript p means that
the pressure is held fixed.
βT is the isothermal compressibility of the fluid, which is given by
1
1 ∂ρ
βT =
=
K
ρ ∂P T
with K the bulk modulus. Values of βT = 10−12 −10−11 Pa−1 are reasonable for Earth’s mantle, with values decreasing
by about a factor of 5 between the shallow lithosphere and core-mantle boundary. This is the percentage increase
in density per unit change in pressure at constant temperature. Both the coefficient of thermal expansion and the
isothermal compressibility can be obtained from the equation of state.
The full set of equations we wish to solve is given by
−∇ · 2η ˙ d (v) + ∇p
1
∇ · v + v · ∇ρ
ρ
∂T
ρCp
+ v · ∇T − ∇ · k∇T
∂t
= ρg
=
0
in Ω
(67)
in Ω
= ρH + 2η ˙ d : ˙ d + αT (v · ∇p)
(68)
in Ω,
(69)
Note that this presupposes that the density is not zero anywhere in the domain.
We use a mixed formulation and therefore keep both velocity and pressure as unknowns. We end up having to
solve the following system:
K
G+W
V
f
·
=
or,
A · X = rhs
P
h
GT + Z
0
Where K is the stiffness matrix, G is the discrete gradient operator, GT is the discrete divergence operator, V the
velocity vector, P the pressure vector. Note that the term ZV derives from term v · ∇ρ in the continuity equation.
Each block K, G , Z and vectors f and h are built separately in the code and assembled into the matrix A and
vector rhs afterwards. A and rhs are then passed to the solver. We will see later that there are alternatives to solve
this approach which do not require to build the full Stokes matrix A.
Remark 1: the terms ZV and WP are often put in the rhs (i.e. added to h) so that the matrix A retains the same
structure as in the incompressible case. This is indeed how it is implemented in ASPECT, see also appendix A of [26].
This however requires more work since the rhs depends on the solution and some form of iterations is needed.
Remark 2: Very often the adiabatic heating term αT (v · ∇p) is simplified as follows: If you assume the vertical
component of the gradient of the dynamic pressure to be small compared to the gradient of the total pressure (in
other words, the gradient is dominated by the gradient of the hydrostatic pressure), then −ρg ' ∇p and then
αT (v · ∇p) ' −αρT v · g. We will however not be using this approximation in what follows.
We have already established that
~τ = Cη BV
43
features
• Q1 × P0 element
• compressible flow
• mixed formulation
• Dirichlet boundary conditions (no-slip)
• isoviscous
• analytical solution
• pressure smoothing
44
A
The main codes in computational geodynamics
In what follows I make a quick inventory of the main codes of computational geodynamics, for crust, lithosphere
and/or mantle modelling.
A.1
ADELI
A.2
ASPECT
A.3
CITCOMS and CITCOMCU
A.4
DOUAR
A.5
GAIA
A.6
GALE
A.7
GTECTON
A.8
ELVIS
A.9
ELEFANT
A.10
ELLIPSIS
A.11
FANTOM
A.12
FLUIDITY
A.13
LAMEM
A.14
MILAMIN
A.15
PARAVOZ/FLAMAR
A.16
PTATIN
A.17
RHEA
A.18
SEPRAN
A.19
SOPALE
A.20
STAGYY
A.21
SULEC
SULEC is a finite element code that solves the incompressible Navier-Stokes equations for slow creeping flows. The
code is developed by Susan Ellis (GNS Sciences, NZ) and Susanne Buiter (NGU).
A.22
TERRA
A.23
UNDERWORLD 1&2
45
B
fieldstone.py
46
C
fieldstone stokes sphere.py
D
fieldstone convection box.py
47
E
fieldstone solcx.py
48
F
fieldstone indentor.py
49
G
fieldstone saddlepoint.py
50
References
[1] M. Albers. A local mesh refinement multigrid method for 3D convection problems with strongly variable viscosity.
J. Comp. Phys., 160:126–150, 2000.
[2] K.-J. Bathe. Finite Element Procedures in Engineering Analysis. Prentice-Hall, 1982.
[3] M. Benzi, G.H. Golub, and J. Liesen. Numerical solution of saddle point problems. Acta Numerica, 14:1–137,
2005.
[4] B. Blankenbach, F. Busse, U. Christensen, L. Cserepes, D. Gunkel, U. Hansen, H. Harder, G. Jarvis, M. Koch,
G. Marquart, D. Moore, P. Olson, H. Schmeling, and T. Schnaubelt. A benchmark comparison for mantle
convection codes. Geophys. J. Int., 98:23–38, 1989.
[5] H.H. Bui, R. Fukugawa, K. Sako, and S. Ohno. Lagrangian meshfree particles method (SPH) for large deformation
and failure flows of geomaterial using elasticplastic soil constitutive model. Int. J. Numer. Anal. Geomech.,
32(12):1537–1570, 2008.
[6] C. Burstedde, G. Stadler, L. Alisic, L.C. Wilcox, E. Tan, M. Gurnis, and O. Ghattas. Large-scale adaptive mantle
convection simulation. Geophy. J. Int., 192:889–906, 2013.
[7] Edmund Christiansen and Knud D. Andersen. Computation of collapse states with von mises type yield condition.
International Journal for Numerical Methods in Engineering, 46:1185–1202, 1999.
[8] Edmund Christiansen and Ole S. Pedersen. Automatic mesh refinement in limit analysis. International Journal
for Numerical Methods in Engineering, 50:1331–1346, 2001.
[9] C. Cuvelier, A. Segal, and A.A. van Steenhoven. Finite Element Methods and Navier-Stokes Equations. D. Reidel
Publishing Company, 1986.
[10] D.R. Davies, C.R. Wilson, and S.C. Kramer. Fluidity: A fully unstructured anisotropic adaptive mesh computational modeling framework for geodynamics. Geochem. Geophys. Geosyst., 12(6), 2011.
[11] P. Davy and P. Cobbold. Indentation tectonics in nature and experiment. 1. experiments scaled for gravity.
Bulletin of the Geological Institutions of Uppsala, 14:129–141, 1988.
[12] Y. Deubelbeiss and B.J.P. Kaus. Comparison of Eulerian and Lagrangian numerical techniques for the Stokes
equations in the presence of strongly varying viscosity. Phys. Earth Planet. Interiors, 171:92–111, 2008.
[13] C.R. Dohrmann and P.B. Bochev. A stabilized finite element method for the Stokes problem based on polynomial
pressure projections. Int. J. Num. Meth. Fluids, 46:183–201, 2004.
[14] J. Donea and A. Huerta. Finite Element Methods for Flow Problems. 2003.
[15] J. Donea and A. Huerta. Finite Element Methods for Flow Problems. John Wiley & Sons, 2003.
[16] Jean Donea and Antonio Huerta. Finite Element Methods for Flow Problems. John Wiley & Sons, 2003.
[17] T. Duretz, D.A. May, T.V. Gerya, and P.J. Tackley. Discretization errors and free surface stabilisation in the finite
difference and marker-in-cell method for applied geodynamics: A numerical study. Geochem. Geophys. Geosyst.,
12(Q07004), 2011.
[18] M. Gerbault, A.N.B. Poliakov, and M. Daignieres. Prediction of faulting from the theories of elasticity and
plasticity: what are the limits? Journal of Structural Geology, 20:301–320, 1998.
[19] Taras Gerya. Numerical Geodynamic Modelling. Cambridge University Press, 2010.
[20] T.V. Gerya, D.A. May, and T. Duretz. An adaptive staggered grid finite difference method for modeling geodynamic Stokes flows with strongly variable viscosity. Geochem. Geophys. Geosyst., 14(4), 2013.
[21] G.H. Golub and C.F. van Loan. Matrix Computations, 4th edition. John Hopkins University Press, 2013.
[22] M. Gunzburger. Finite Element Methods for Viscous Incompressible Flows: A Guide to Theory, Practice and
Algorithms. Academic, Boston, 1989.
[23] T.J.R. Hughes. The Finite Element Method. Linear Static and Dynamic Finite Element Analysis. Dover Publications, Inc., 2000.
51
[24] T.J.R. Hughes, W.K. Liu, and A. Brooks. Finite element analysis of Incompressible viscous flows by the penalty
function formulation. J. Comp. Phys., 30:1–60, 1979.
[25] Hoon Huh, Choong Ho Lee, and Wei H. Yang. A general algorithm for plastic flow simulation by finite element
limit analysis. International Journal of Solids and Structures, 36:1193–1207, 1999.
[26] L. Jolivet, P. Davy, and P. Cobbold. Right-lateral shear along the Northwest Pacific margin and the India-Eurasia
collision. Tectonics, 9(6):1409–1419, 1990.
[27] L.M. Kachanov. Fundamentals of the Theory of Plasticity. Dover Publications, Inc., 2004.
[28] M. Kronbichler, T. Heister, and W. Bangerth. High accuracy mantle convection simulation through modern
numerical methods . Geophy. J. Int., 191:12–29, 2012.
[29] W. Leng and S. Zhong. Viscous heating, adiabatic heating and energetic consistency in compressible mantle
convection. Geophy. J. Int., 173:693–702, 2008.
[30] W. Leng and S. Zhong. Implementation and application of adaptive mesh refinement for thermochemical mantle
convection studies. Geochem. Geophys. Geosyst., 12(4), 2011.
[31] D.S. Malkus and T.J.R. Hughes. Mixed finite element methods - reduced and selective integration techniques: a
unification of concepts. Comput. Meth. Appl. Mech. Eng., 15:63–81, 1978.
[32] P. Molnar and P. Tapponnier. Relation of the tectonics of eastern China to the India-Eurasia collision: Application
of the slip-line field theory to large-scale continental tectonics. Geology, 5:212–216, 1977.
[33] L. Moresi, S. Quenette, V. Lemiale, C. Mériaux, B. Appelbe, and H.-B. Mühlhaus. Computational approaches to
studying non-linear dynamics of the crust and mantle. Phys. Earth. Planet. Inter., 163:69–82, 2007.
[34] G. Peltzer and P. Tapponnier. Formation and evolution of strike-slip faults, rifts, and basins during the india-asia
collision: an experimental approach. J. Geophys. Res., 93(B12):15085–15177, 1988.
[35] A.E. Pusok, B.J.P. Kaus, and A.A. Popov. On the Quality of Velocity Interpolation Schemes for Marker-in-Cell
Method and Staggered Grids. Pure and Applied Geophysics, pages doi:10.1007/s00024–016–1431–8, 2016.
[36] T. Rabczuk, P.M.A. Areias, and T. Belytschko. A simplified mesh-free method for shear bands with cohesive
surfaces . Int. J. Num. Meth. Eng., 69:993–1021, 2007.
[37] J.N. Reddy. On penalty function methods in the finite element analysis of flow problems. Int. J. Num. Meth. Fluids, 2:151–171, 1982.
[38] J. Revenaugh and B. Parsons. Dynamic topography and gravity anomalies for fluid layers whose viscosity varies
exponentially with depth. Geophysical Journal of the Royal Astronomical Society, 90(2):349–368, 1987.
[39] H. Schmeling, A.Y. Babeyko, A. Enns, C. Faccenna, F. Funiciello, T. Gerya, G.J. Golabek, S. Grigull, B.J.P.
Kaus, G. Morra, S.M. Schmalholz, and J. van Hunen. A benchmark comparison of spontaneous subduction
models - Towards a free surface. Phys. Earth. Planet. Inter., 171:198–223, 2008.
[40] D.W. Schmid and Y.Y. Podlachikov. Analytical solutions for deformable elliptical inclusions in general shear.
Geophy. J. Int., 155:269–288, 2003.
[41] G. Schubert, D.L. Turcotte, and P. Olson. Mantle Convection in the Earth and Planets. Cambridge University
Press, 2001.
[42] J. Suckale, J.-C. Nave, and B.H. Hager. It takes three to tango: 1. Simulating buoyancy-driven flow in the
presence of large viscosity contrasts. J. Geophys. Res., 115(B07409), 2010.
[43] P. Tackley. Three-dimensional models of mantle convection: Influence of phase transitions and temperaturedependent viscosity. PhD thesis, California Institute of Technology, 1994.
[44] Paul Tapponnier and Peter Molnar. Slip-line field theory and large-scale continental tectonics. Nature, 264:319–
324, November 1976.
[45] M. Thielmann, D.A. May, and B.J.P. Kaus. Discretization errors in the Hybrid Finite Element Particle-In-Cell
Method. Pure and Applied Geophysics, 2014.
[46] C. Thieulot, P. Fullsack, and J. Braun. Adaptive octree-based finite element analysis of two- and three-dimensional
indentation problems. J. Geophys. Res., 113:B12207, 2008.
52
[47] R.A. Trompert and U. Hansen. On the Rayleigh number dependence of convection with a strongly temperaturedependent viscosity. Physics of Fluids, 10(2):351–360, 1998.
[48] X. Yu and F. Tin-Loi. A simple mixed finite element for static limit analysis. Computers and Structures, 84:1906–
1917, 2006.
[49] S. Zhong. Analytic solutions for Stokes flow with lateral variations in viscosity. Geophys. J. Int., 124:18–28, 1996.
[50] O. Zienkiewicz and S. Nakazawa. The penalty function method and its application to the numerical solution of
boundary value problems. The American Society of Mechanical Engineers, 51, 1982.
[51] O.C. Zienkiewicz, M. Huang, and M. Pastor. Localization problems in plasticity using finite elements with
adaptive remeshing. International Journal for Numerical and Analytical Methods in Geomechanics, 19:127–148,
1995.
[52] O.C. Zienkiewicz, C. Humpheson, and R.W. Lewis. Associated and non-associated visco-plasticity and plasticity
in soil mechanics . Géotechnique, 25(4):671–689, 1975.
53
Index
Q1 × P0 , 13, 19, 25, 31, 34, 42, 44
analytical solution, 13, 19, 34, 39, 42, 44
arithmetic mean, 31
basis functions, 6
Boussinesq, 5
bulk modulus, 43
bulk viscosity, 4
compressibility, 43
compressible flow, 44
dynamic viscosity, 4
geometric mean, 31
harmonic mean, 31
incompressible flow, 19, 25, 31, 34, 39, 42
isothermal, 13, 19, 25, 31, 34, 39, 42
isoviscous, 13, 34, 42, 44
mixed formulation, 34, 39, 42, 44
non-isoviscous, 19, 25, 31, 39
nonconforming Q1 × P0 , 39
nonlinear rheology, 25
particle-in-cell, 31
penalty formulation, 7, 13, 19, 25, 31
pressure smoothing, 34, 39, 42, 44
second viscosity, 4
thermal expansion, 43
weak form, 9
54
Source Exif Data:
File Type : PDF File Type Extension : pdf MIME Type : application/pdf PDF Version : 1.5 Linearized : No Page Count : 54 Page Mode : UseOutlines Author : Title : Subject : Creator : LaTeX with hyperref package Producer : pdfTeX-1.40.18 Create Date : 2018:12:17 20:07:00+01:00 Modify Date : 2018:12:17 20:07:00+01:00 Trapped : False PTEX Fullbanner : This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) kpathsea version 6.2.3EXIF Metadata provided by EXIF.tools