Numerical Modeling of Earth Systems PDF
Numerical Modeling of Earth Systems PDF
Numerical Modeling of Earth Systems PDF
Thorsten W. Becker
Department of Earth Sciences,
University of Southern California, Los Angeles CA, USA
and
Boris J. P. Kaus
University of Mainz, Germany
March 8, 2016
1 Preface 7
1.1 Purpose of these lecture notes . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.3 Abbreviations used . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.4 Typesetting conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5 Other resources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
I Introduction 12
2 Introduction to numerical geodynamics 13
2.1 Numerical methods in the Earth Sciences . . . . . . . . . . . . . . . . . . . . 13
2.1.1 Philosophy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.2 Goals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.1.3 Overview of applications of numerical methods for Earth sciences . 14
2.1.4 Classification of numerical problems & solution methods . . . . . . . 18
2.2 Examples of applications for numerical methods . . . . . . . . . . . . . . . . 19
2.2.1 Linear inverse problems . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.2 Ordinary differential equations . . . . . . . . . . . . . . . . . . . . . . 19
2.2.3 Partial differential equations . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.4 Numerical solution methods . . . . . . . . . . . . . . . . . . . . . . . 22
2.3 Computing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3.1 Hardware issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3.2 Software - Computer Languages . . . . . . . . . . . . . . . . . . . . . 26
2.3.3 Elements of a computer program . . . . . . . . . . . . . . . . . . . . . 28
2.3.4 Guiding philosophy in writing a computer program . . . . . . . . . 28
2.3.5 Guidelines for writing efficient code . . . . . . . . . . . . . . . . . . . 29
2.4 Scaling analysis and non-dimensional numbers . . . . . . . . . . . . . . . . . 33
2.4.1 Scaling analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.4.2 Non-dimensionalization . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.4.3 Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
Bibliography 213
Index 219
Preface
7
CHAPTER 1. PREFACE
1.2 Acknowledgments
Parts of the original course notes are inspired by, or partially follow very closely, the treat-
ment in the numerical analysis lecture notes by Spiegelman (2004), the notes by Schmeling
(1994), and the textbook by Press et al. (1993) on numerical analysis, and the textbook by
Hughes (2000) on finite element methods.
Some of the finite difference exercises and parts of the scaling homework assignment
are loosely based on an ETH Zürich course taught by Yuri Podladchikov. Daoyuan Sun
provided the wave propagation exercise, and most of the finite element exercises are built
directly on the MILAMIN matlab software which is openly distributed by Dabrowski et al.
(2008). The multigrid problem set is based on an exercise by Zhong (2008), and one of
the elasticity, finite element exercises is inspired by a Numerical Analysis class taught
by Harro Schmeling at Frankfurt University in 1997 (Schmeling, 1994). A web blog at
http://gwlab.blogspot.com/2012/08/inverse-transformation-of-parametric.html
was helpful in coming to grips with inverse transformations.
Students who took the class at USC Earth Sciences in the Fall of 2005, 2008, 2013, and
2016 provided valuable feedback.
USC PhD candidate Francois Cadieux helped greatly with the type setting of the origi-
nal lecture notes in 2010, and also provided some very useful, additional material on high
Reynolds number fluid dynamics algorithms.
Partial funding for course development was provided by the US National Science
Foundation under CAREER grant EAR-0643365.
FD Finite differences
FE Finite elements
IC Initial condition
BC Boundary conditions
Mathematical symbols are denoted by a for scalars, a for vectors, and A for ten-
sors/matrices.
This text is meant to be fairly self-contained, but there is some related material which
we refer to often. This includes the following online resources:
• The finite element method by Hughes (2000), for further detailes on the finite element
method.
For finite elements, for example, there are, of course, a number of good textbooks.
Those include Kwon and Bang (1996), which provides a clear, step-by-step, introduction
with many MATLAB program examples, and the classic by Bathe (2007) for a comprehen-
sive reference. For computational geodynamics, we started when no appropriate text-
books were out there. Since then, Ismail-Zadeh and Tackley (2010) provided a broader, but
less-detailed, overview of general methods, and Gerya (2009) a very detailed, MATLAB
based approach that is focused on finite differences only, however.
Background material Short discussions of basic math, continuum mechanics, and MAT-
LAB are found in secs. 6, 7, and 8, respectively. Readers not familiar with this material
may wish to review those chapters.
Availability and contact A PDF of the lecture notes and MATLAB exercises as used for
a graduate class at USC (GEOL557) can be found on course web site
http://geodynamics.usc.edu/~becker/teaching-557.html
Introduction
12
Chapter 2
2.1.1 Philosophy
• Avoid black boxes (e.g. commercial codes) in general. They may or may not do
what you like them to do; if they don’t, you’re out of luck because if you cannot
modify the source code. Exception for “good” black boxes are matrix solver and
linear algebra packages, generally speaking (but see sec. 5.4).
• There are no “Easy” or “Model my field data” buttons, but you don’t have to be a
math-whiz either!
2.1.2 Goals
• Provide you with a basic understanding of numerical modeling, using solid Earth
science problems as an example.
• Allow you to solve simple research problems using tools presented here, and, more
importantly, allow you to write new programs and solve different scientific ques-
tions yourself independently.
13
CHAPTER 2. INTRODUCTION TO NUMERICAL GEODYNAMICS
• be comprehensive
Differential equations of the partial (PDE) or ordinary (ODE) kind, which can be solved
with
Inverse problems where a structural or physical model of the Earth is inferred from (a
potentially very large) set of data. Methods may include
• linear matrix inversion, least squares (with challenges related to those in finite
element methods)
• Monte Carlo, simulated annealing
• Genetic Algorithms
Figure 2.1: Estimated slip distribution on southern California faults for shear stress (left) and
Coulomb stress free, simultaneous slip, boundary element method computations (modified from
Becker and Schott, 2002). This method (e.g. Crouch and Starfield, 1983) uses Greens functions (here
computed following Okada, 1992) and either summation for estimating τ given slip distributions,
u, on faults, or solution of a large, linear system of equations to infer u for specified stress on the
faults.
Figure 2.2: Spherical, visco-plastic mantle convection computation results showing plate-like sur-
face motions (left: velocity vectors on top of viscosity, right: cold temperature isosurface in the in-
terior showing one-sided subduction due to slab rollback), modified from Foley and Becker (2009).
We use the CitcomS (Zhong et al., 2000) finite element software to solve the Stokes and energy
equations for an infinite Prantl number, incompressible fluid (cf. Zhong et al., 2007).
Figure 2.3: Global mantle circulation computation using mantle tomography to infer density dis-
tributions in the mantle and solving Stokes equation (modified from Becker and Faccenna, 2011),
see (cf. Hager and O’Connell, 1981; Zhong et al., 2000; Becker, 2006). Visualization with Paraview
(Kitware, Inc., 2006).
3. solve
• dimensional analysis
• analytical solution
– check if this is a standard problem someone else has solved
– check if terms can be neglected to simplify
– check if equations can be linearized
– numerical solution
Inverse problem
1. Formulate a model (e.g. Earth’s mantle wave speed variations are smooth)
3. Collect data
∂y
= f (y) (2.1)
∂t
∂x
= v (2.2)
∂t
→ Runge-Kutta, Burlisch-Stoer integration methods, for example (this will be dealt
with next, see 3.2).
∂2 f ∂2 f ∂2 f ∂2 f
k1 + k 2 + k 3 + k 4 + (2.3)
∂x2 ∂y2 ∂z2 ∂t2
∂f ∂f ∂f ∂f
p1 + p2 + p3 + p4 + (2.4)
∂x ∂y ∂z ∂t
f ( x, y, z, t) = g( x, y, z, t). (2.5)
elliptic PDEs : B2 < 4AC, or all k i are non-zero and have the same sign;
hyperbolic PDEs : B2 > 4AC, or all k i are non-zero and one has a different sign from the
others;
parabolic PDEs : B2 = 4AC, or all k i but one are non-zero and have the same sign.
Distinguishing between these different types of PDEs is important because the types
determine which boundary conditions and which initial conditions are neded. For example,
elliptic problems require boundary conditions, but parabolic and hyperbolic problems als
require initial conditions.
As for the boundary types and conditions, the following rules apply
type boundary conditions
hyperbolic open Cauchy
parabolic open Dirichlet or Neumann
elliptic closed Dirichlet or Neumann
where the boundary condition type means
Dirichlet : specifying f itself;
∂f ∂f
Neumann : specifying ∂n where n is the normal derivative on the boundary, ∂n = ∇ f · n,
where n is the normal at each point;
∂f
Cauchy : where both f and ∂n need to be specified.
• gravity potential equation: Gauss’ law states that the gravitational acceleration, g,
due to density in the volume, ρ, is
∇ g = −4πGρ. (2.8)
The gravitational field is conservative, and can hence be expressed as a potential Φ,
g = −∇φ such that
∇2 Φ = 4πGρ (2.9)
results.
For an isotropic and homogeneous medium, these equations correspond to the Poisson
2 ∂2 f ∂2 f ∂2 f
∇ f = 2 + 2 + 2 = g( x, y, z), (2.10)
∂x ∂y ∂z
or the Laplace equation
∇2 f = 0. (2.11)
∂2 f ∂f
∇2 f = g( x, y, z, t) + µ + γ , (2.12)
∂t2 ∂t
where boundary conditions may also be a function of time t, and we need additional
initial conditions which provide f ( x, y, z, t0 ) = Ψ ( x, y, z).
For µ = 0, this is a parabolic PDE, with the classic example of the time-dependent heat
equation, an example of a diffusion equation
∂T
ρc p = k∇2 T + ρH (2.13)
∂t
∂T H
= κ ∇2 T + . (2.14)
∂t cP
with thermal diffusivity
k
κ= (2.15)
ρc p
based on conductivity, k, density, ρ, and heat capacity, c P , at constant pressure, and inter-
nal heat production per volume, H.
If γ = 0 and µ > 0, the equation is a hyperbolic PDE, in the case of g = 0 a wave
equation,
∂2 f
= c2 ∇2 f , (2.16)
∂t2
where c is a velocity.
If we are interested in periodic solutions of eq. (2.16) only, we can use normal modes
and turn the hyperbolic PDE in a eigenvalalue problem by assuming
f = F ( x, y, z)eiωt , (2.17)
which yields
− ω 2 F = c2 ∇2 F. (2.18)
The resulting frequencies ω are the eigenfrequencies, e.g. of the vibrational modes of the
whole Earth induced by earthquakes.
In case both µ 6= 0 and γ 6= 0, we have a damped wave equation. In this case, two
initial conditions are needed,
f ( x, y, z, t0 ) = Ψ ( x, y, z) (2.19)
∂f
( x, y, z, t0 ) = ξ ( x, y, z) (2.20)
∂t
Cons
∂2 u
Z b Z b
− dx w 2 − dx ws = 0, (2.24)
a ∂x a
where h = ∂u
∂x ( b ). This is the weak form of the PDE, and picking w as a low order polyno-
mial forms the basis of the FE method.
Pros
Cons
u = ∑ a j eijx (2.27)
s = ∑ bj eijx (2.28)
(2.29)
from which
bj
aj = − 2 . (2.31)
j
The approach then consists in solving by computing b j by FFT of s, then use eq. (2.31),
then find f from inverse FFT of a j .
Pros
Cons
• non-local
• poor performance for lateral variations in material properties (need to iterate)
2.3 Computing
Here, we very briefly review some hardware related issues, give a few general program-
ming tips, and then move on to get started programming using the MATLAB language
in the next section. You can also refer to the hardware notes of Press et al. (1993) for some
background on machine architecture.
Memory
1 MB (megabyte) corresponds 1024 × 1024 bytes; 1 GB = 1024 MB. As of 2008, your PC
will have likely have at least ∼ 2 GB of Random Access Memory or RAM (as opposed to
hard drive space) meaning you can store how many floats and doubles? To increase the
available memory, one can use formerly called “supercomputers”. Those consist these
days mainly of
Distributed memory machines e.g. 200 × 2×quadcore (8 Central Processing Units or CPUs)
×8 GB RAM machines which need specially designed software to make use of par-
allelism, e.g. Message Passing Interface or MPI.
Shared memory machines This is the more expensive, old school approach where sev-
eral CPUs can share a larger than normal (e.g. 256 GB) memory. Compilers can
sometimes help make your code make use of “parallelism”, i.e. having the compu-
tational time decrease by using more than one core or CPU. Right now, typical PCs
can be considered shared memory (multi-core, i.e. CPU) machines.
• convenient debugging
Cons:
Libraries
• NETLIB
• BLAS
• LAPACK
• PETSc
Cons:
• Need to compile
% This is the main program. Notice the ’%’ symbol - it means this line is
% a comment and will be ignored at run time.
i = 0; % assign integer variable for loop
n = 100; % some number of elements
x = zeros(n,1); % allocate and initialize a vector x[] with n elements
y = 1;
for i = 1:n % loop from i = 1, 2, ..., n
x(i) = y^2; % assign some value
y = y+2; % increment variable
end % close loop
% notice the statements inside the loop are indented.
i = 1;
while (i <= n) % different loop construct
x(i) = mysin(x(i)); % function call
i = i+1;
printf("%g\n", x(i)); % output statement
end
• Break the task down into small into small pieces that can be reused within the
same program or in another program
• Test each part well before using it in a larger project to make code more robust.
• ensure that each subroutine gives error messages, in case non sensible input
arguments are given.
• do not ignore compiler warnings
Don’t use special tricks/packages that might not be available on other platforms.
3. Comment
• Add explanatory notes for each major step, strive for a fraction of comments to
code
≥ 30%
This will help re-usability, should you or someone else want to modify the code
later.
6. Visualize you intermediate results often (But don’t print it all out in color!)
Bugs in the code can often be seen easily when output is analyzed graphically, and
may show up as, e.g.
Object oriented programming forces you to follow rules 1 & 4 (not so much 2). Editors
and advanced development environments (such as the MATLAB DE) help with 3 & 6.
2. Use nested loops that are sorted by the fastest/major index, because memory access
is faster that way. The storage depends on the computer larguage (C vrs. FOR-
TRAN). e.g. in MATLAB , you would write
Even better, in MATLAB (and languages such as FORTRAN90) you can vectorize, i.e.
write symbolically for a vector x
x = x .* y
for the example above. MATLAB internally takes take that the looping is taking care
of in the most efficient way. This can make a huge difference, vectorize whenever
you can in.
if optional and usually zero, comment it out using pre-processor directives. I.e. in
C, you would write the code like so
#ifdef DEBUG
% code here for debugging version of program
#else
% code here for the regular version of program
#end
gcc -DDEBUG
depending on if you want those statements to be executed when the program runs.
for i = 1:n
x(i) = x(i)/180*pi;
end
It is better to do
fac = pi/180;
for i = 1:n
x(i) = x(i)*fac;
end
because it entails one less division per step. In MATLAB , it’s better still to use the
vectorized version, x=x.*fac.
The more eyes, the less bugs, and the better the performance.
1 http://github.com/
Reading
• Turcotte and Schubert (2002), Google, and Wikipedia for reference and material pa-
rameters
τ = 2η ε̇ (2.32)
where ε̇ is the strain-rate (in units of s−1 ). Say, we wish to estimate the typical amplitudes
of shear stress in a part of the crust that we know is being sheared at some (e.g. plate-)
velocity u over a zone of width L. The strain-rate in 3-D is really a tensor, ε̇, with 3 × 3
components that depends on the spatial derivatives of the velocity like so
!
1 ∂ui ∂u j
ε̇ = ε̇ ij = + , (2.33)
2 ∂x j ∂xi
and has to be either constrained by kinematics or inferred from the full solution. How-
ever, for our problem, we only need a “characteristic” value, i.e. correct up to a factor of
ten or so. Strain-rate is physically the change in velocity over length, and the characteris-
tic strain-rate is then given by
v
ε̇ ∝ (2.34)
L
where ∝ means “proportional to”, or “scales as”, to indicate that eq. (2.34) is not exact.
Assuming we know the viscosity η, we can then estimate the typical stress in the shear
zone to be
v
τ ∝ 2η . (2.35)
L
If you think about the units of all quantities involved (“dimensional analysis”), then this
scaling could not have worked out any other way. Viscosity is Pa s (stress times time),
velocity m/s (length over time), so stress=Pa s m/(s m)=Pa as it should be. 2
Note I: Whenever you work out, or type up, a new equation, it is always a good idea
to check if the units on both sides make sense.
Note II: The scaling of velocities and stress for a buoyancy driven problem, such as the
Stokes sinker discussed below, is entirely different!
2.4.2 Non-dimensionalization
A complementary approach that also takes into account the order of magnitude of vari-
ables is to simplify the governing equations by defining “characteristic” quantities and
then dividing all properties by those to make them “non-dimensional”. This way, the
non-dimensional quantities that enter the equation on their own should all be of order
unity so that the resulting collection of parameters in some part of the equation measures
their relative importance.
A classic example for this is based on the Navier Stokes equation for an incompress-
ible, Newtonian fluid. When body forces driving flow are due to temperature T fluctua-
tions in (the Earth’s) gravitational field
Dv
ρ = −∇ pd + η ∇2 v + ρ0 αTg (2.36)
Dt
where D is the total, Lagrangian derivative operator
D ∂
= + v · ∇, (2.37)
Dt ∂t
v velocity, ∇ the Nabla derivative operator ∇ = {∂/∂x, ∂/∂y, ∂/∂z}, t time, pd the dy-
namic pressure (without the hydrostatic part), η the viscosity, ρ0 reference density, α, and
g gravitational acceleration. One can now choose (as mentioned before for the Lorenz
equations) typical quantities that can be derived from the given parameters such as a ∆T
temperature difference, a fluid box height d, and some choice for the timescale. All other
characteristic values for physical properties can then be derived from those choices (see,
e.g., discussion in Ricard, 2007).
2 We will always use SI units unless it’s inconvenient for Earth applications, where we might use multi-
ples of SI units such as cm/yr instead of m/s for velocities. Also note that one year has roughly π · 107 s
(accurate up to 0.5%), i.e. 1 cm/yr is ≈ 3.2 · 10−10 m/s, and that you should account for leap years for
geological time conversions, meaning that 365.25 × 24 × 60 × 60 is the right number of seconds per year.
Let us assume that we are dealing with a box heated from below and cooled from
above, i.e. held at constant temperature difference of ∆T = Ttop − Tbot , with no internal
heating (Rayleigh-Benard problem). A typical choice for a characteristic timescale is to
use the diffusion time that can be constructed from the thermal diffusivity, κ, in the energy
equation
DT
= κ ∇2 T (2.38)
Dt
(no heat sources) that couples with the momentum equation, eq. (2.36). Because κ has
units of length2 /time, any diffusion-related time scale td for a given length l has to work
out like
d2
td ∝ , (2.39)
κ
by dimensional analysis. This relationship is hugely important for all diffusional pro-
cesses.
Using the characteristic quantities f c which result from this scaling and using l = d,
for all variables in eq. (2.36) and eq. (2.38), all other properties can be derived, e.g.
d vc
vc = ε̇ c = τc = η ε̇ c Tc = ∆T (2.40)
tc d
we divide all variables (spatial and temporal derivatives are dealt with like space and
time variables) to make them unit-less, non-dimensional f 0 = f / f c , and eq. (2.36) can
then be written as
1 D 0 v0
0 0
= −∇0 p0 + (∇0 )2 v0 − RaT 0 ez (2.41)
Pr D t
where we have used g = − gez . I.e., all material parameters have been collected in two
unit-less numbers after non-dimensionalization, the Prandlt number,
η ν
Pr = = (2.42)
ρκ κ
and the Rayleigh number 3
ρ0 gα∆Td3
Ra = . (2.44)
κη
3 Inthe derivation above, we have assumed that the system is heated from below and viscosity is con-
stant. The Rayleigh number is therefore valid for this bottom-heated case only. In Earth’s mantle, internal
heating (due to decay of radioactive elements) is at least equally important (e.g. Jaupart et al., 2007; Lay et al.,
2008). For the case of pure internal hearing, the Rayleigh-number is given by
ρ0 gαH 0 d5
Ra H 0 = , (2.43)
kκη
where k is conductivity and H 0 is the rate of internal heat generation per volume (H 0 = ρ0 H where H is per
unit mass). We can identify ∆T in (2.44) with H 0 d2 /k which makes sense, since the total heat flux, Q, should
scale as H 0 d3 and Q ∝ k ∆T
d . Also, rock viscosity depends on a range of quantities, including temperature
and strain-rate, making it imperative to properly (log) average viscosity when computing effective Rayleigh
numbers (e.g. Christensen, 1984).
The second way of expressing Pr uses the kinematic viscosity, ν = η/ρ, which, like κ, has
units of m2 /s; this makes it clear that Pr measures diffusion of momentum vrs. diffusion
of heat. Ra measures the vigor of convection by balancing buoyancy forces associated
with advection against diffusion and viscous drag.
Exercise: Verify this recasting of the Navier-Stokes equation by plugging in the non-
dimensionalized quantities.
Often, we then just drop the primes and write the equation like so
1 Dv
= −∇ p + ∇2 v − RaTez (2.45)
Pr Dt
where it is implied that all quantities are used non-dimensionalized (also see sec. 3.3.2).
This equation may still be hard to solve, but at least we now have sorted all material
parameters into two numbers, Ra and Pr.
Note I: The non-dimensional versions of the equations are also the best choice if you
want to write a computer program for a physical problem. Using non-dimensionalized
equations, all terms should be roughly of order unity, and the computer will not have
to multiply terms that are very large in SI units (e.g. η) with those that are very small,
reducing round-off error (e.g. v, what is the order of magnitude of η and of |v| for mantle
convection?).
Note II: This also means that when some geophysicist’s convection code spits out, say,
velocities, you will have to check what units those have, and more often than not you’ll
have to multiply by the vc from above to get back m/s, which you’ll then convert to
cm/yr.
Note III: You’ll also note that a few geodynamics papers will not provide the scaled
quantities used so that you can go back to SI units; sometimes this is because the values
used for the parameters in the models stray significantly from typical Earth values.
Particularly the Rayleigh number is key for mantle convection, because we typically
use the infite Prandtl number approximation, (Pr → ∞, why?). In this case, eq. (2.36)
simplifies to the Stokes equation,
η ∇2 v = ∇ pd − ρ0 αTg. (2.46)
Both Pr and Ra are discussed below. Fluid dynamics is full of these non-dimensional
numbers which are usually named after some famous person because they are so power-
ful. Any fluid that has the same Ra and Pr number as another fluid will behave exactly
the same way in terms of the overall style of dynamics, such as the resulting average
temperature structure and up and downwelling morphology.
The actual time scales of convection, e.g., may, however, be very different for two
systems at the same Rayleigh number (because of vc being different). This behavior al-
lows, for example, to conduct analog, laboratory experiments of mantle convection (e.g.
Jacoby and Schmeling, 1981; Faccenna et al., 1999). When conducting such experiments,
care needs to be taken that all relevant non-dimensional numbers agree between the real
Earth problem and the laboratory experiment (e.g. Weijermars and Schmeling, 1986). Also,
when changing length scales and material, different physical effects such as surface ten-
sion may matter in the lab, while they are irrelevant for mantle convection in general (see,
e.g., sec. 6.7 of Ricard, 2007, for a discussion of Mahagoni convection).
From an analytical point of view, if the non-dimensional quantities are either very
large or very small, we can simplify the full equations to more tractable special cases. For
a nice and more comprehensive treatment of this section, you may want to refer to Ricard
(2007).
2.4.3 Problems
1. For all of the following non-dimensional numbers, discuss briefly (2-3 sentences)
the processes which these numbers measure, e.g. by contrasting system behavior for
Th = 0 and Th = ∞, where Th is some non-dimensional number.
For each number, give numerical estimates for the Earth, at the present day. Docu-
ment your choices (i.e. providing references) for individual parameters before com-
puting joint quantities, mention where you got the estimates from, and what the
implications for Earth in terms of the dynamics are. A neat way to organize this
might be to use a table for each dimensionless number with appropriate headings
(e.g. parameter, estimate, reference).
You might have to look up definitions and other reference material, e.g. in a geody-
namics text, or on Google (note: don’t trust everything on the web . . . ). There are no
unique answers for this part of the problem set, and you will often have to decide on
an example problem for which you’ll pick a characteristic quantity such as length.
Some answers are actively debated in the literature.
(e) Deborah number for the subducting oceanic lithosphere, and for a laboratory
experiment on rock deformation. The Deborah number can be defined as
tr
De = (2.49)
tp
where you can use a Maxwell time
η
tM = (2.50)
µ
for the relaxation time tr , and t p is the time scale of observation. The Maxwell
time measures the visco-elastic relaxation time of a body with viscosity η and
shear modulus µ (think post-glacial rebound).
• What are characteristic Maxwell times for the crust? The upper mantle?
2. (a) Consider a solid, sinking sphere of radius a in a fluid of viscosity η and gravi-
tational pull g, and a density Stokes velocity contrast between sphere and fluid
of ∆ρ. Solve for the approximate sinking velocity of this “Stokes” sinker by
equating the gravitational pull force FP = ∆Mg = V∆ρg with the shear force
acting on the sphere’s area A, FS ∝ τA ∝ Aη ε̇ c . Here, I’ve used ∆M for the
mass anomaly, and V for the volume of the sphere. All equations follow from
F = ma and stress = force / area and some geometry.
Note: The full equation for a Stokes sinker is only very weakly dependent on
the viscosity of the sinker, ηs , itself, but scales mainly with the ambient vis-
cosity η. For ηs /η → ∞ and ηs /η → 0, the prefactor changes from 2/9 to
1/3, respectively (see further discussion in Becker and Faccenna, 2009, for the
subduction context).
(b) For flow induced by a Stokes sinker, does the stress scale with η and/or ∆ρ?
How does that compare with the velocities?
Note: The scaling of v and τ with combinations of ∆ρ and η are among the
most fundamental relations of mantle dynamics (velocities v might be the plate
velocities, dynamic topography scales with τ, for example).
(c) Estimate the Stokes velocity by dimensional analysis as in (a), but now assum-
ing that the viscosity of the fluid obeys a power-law,
τ n ∝ η ε̇ (2.51)
(for rocks, n ∼ 3) instead of
τ ∝ η ε̇ (2.52)
for Newtonian creep as assumed above. (These equations are written sloppily
and don’t have the right units. For correct units, consider a relationship like
τ (τ/µ)n−1 = η ε̇, where η is a material parameter, but you may use eq. (2.51)
for the scaling analysis.)
H h ρs
g
ρm ηm
(d) Estimate the rise velocity of a plume head large enough to cause the Deccan
traps.
3. You are moving the top of a fluid layer of height d at constant speed v(z = d) = v0 ,
and the fluid is held fixed at the bottom at z = 0. In this case, the laminar solution
for the flow velocity is a linear decrease of velocity with depth to v(0) = 0 at the
bottom.
(a) What material parameters set the stress in the fluid? What determines the
strain-rate and how does it vary with depth?
(b) Now consider two fluid layers, with the top fluid viscosity larger than the bot-
tom one by a factor of two. Sketch the solution for the dependence of v(z).
4. Using dimensional analysis, such as used above for the Stokes sinker, estimate the
velocity of a volcanic eruption (see Figure below for parameters).
Hint: You might want to proceed by first using the equations for laminar, pressure-
driven (look up “Hagen-Poiseuille”) flow in a pipe of radius R, and then estimate
the pressure difference from Figure 2.4.
40
Chapter 3
3.1 Introduction
Reading: Press et al. (1993), Chap. 17; Spiegelman (2004), Chap. 4; Spencer and Ware (2008),
sec. 16.
ODE An equation that involves the derivative of the function we want to solve for, and
that has only one independent variable (else it’s a PDE).
For example:
dy
= f ( x ), which can be solved by integration (3.1)
dx Z
y= f ( x ) dx + C (3.2)
d2 y dy
2
+ q( x ) = r(x) (3.3)
dx dx
is “second order”. However, we can always reduce ODEs to sets of first order equations.
For eq. (3.3), define
dy
= z( x ), then (3.4)
dx
dz
= r ( x ) − q ( x ) z ( x ). (3.5)
dx
41
CHAPTER 3. SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS
Or, in general
dyi ( x ) dy
= f i ( x, y, . . . , y N ) or = f ( x, y) (3.6)
dx dx
is a system for N coupled ODEs, all dependent on the independent variable x, which is
typically time, t. y is the solution vector we want to solve for. The actual solution of ODEs
will depend on the types of boundary conditions on y and the initial conditions.
We can distinguish between initial value and two point boundary values problems.
Dc dc
= + v · ∇c = f (3.10)
Dt dt
We can break down the integral into step sizes h from ti to ti + h with n = (t f − t0 )/h
partial integrals such that we only need to solve
Z ti + h
I= f (t, y(t)) dt (3.13)
ti
becomes the rule to advance y from ti to ti + h. This is the Euler method, and a really bad
idea. Consider the graphical representation in Figure 3.1, which shows that (3.14) is just
a simple extrapolation of y based on the slope at ti , which is given by equation (3.11). If y
has some curvature to it, the Euler scheme will lead to large errors quickly!
We can Taylor expand (eq. 6.4) y around t0 to get
h2 d2 y
y ( t i + h ) ≈ y ( t i ) + h · f ( y ( t0 ), t0 ) + +... (3.16)
2 dt2
Notice that the error of the Euler scheme goes as O (“order of”)(h2 ), and the scheme itself
is therefore, by definition, only accurate to first order. This means that tiny time steps
would have to be taken for a good solution. There are several improvements to the Euler
method.
1. The midpoint method of Figure 3.2 evaluates the derivative of y w.r.t. to t first at
half the Euler step
dy dy h
t + h/2, y (ti ) + (t ) (3.17)
dt i dt i 2
k 1 = h f ( ti , yi ) (3.20)
h 1
k 2 = h f ti + , yi + k 1 (3.21)
2 2
y i +1 = y i + k 2 + O h 3 , (3.22)
where it can be shown that this method is second order accurate. Note that higher
accuracy has come at a cost, f now needs to be evaluated twice and once at a y value
different from yi , and there are overall more operations per time step. However,
since the error is now O h3 , we can take larger time steps.
There are several avenues to refine the midpoint method further, but in general the
k 1 = h f ( ti , yi )
h k1
k 2 = h f ti + , yi +
2 2
h k2
k 3 = h f ti + , yi +
2 2
k4 = h f (ti + h, yi + k3 )
k k2 k3 k4
y i +1 = y i + 1 + + + + O h5 . (3.23)
6 3 3 6
The next improvement to the Runge Kutta method is to adapt the step size h during
forward integration, according to an a comparison of an on-the-fly error estimate with
some desired accuracy (e.g. Press et al., 1993, sec. 16.2). If you have a simple problem and
want to implement your own method, as in the problem set below, one way to test the
behavior of the forward routine is to successively half h and keep track of the deviations
of |y| at x f to make sure things converge. This is of course assuming that the parameters
for f remain within a reasonable range as used for the test, and adaptive step size is
preferred.
Moreover, especially tricky functions f require special methods, and Press et al. (1993),
chap. 16, discusses these, for example. Some of the more fancy methods are implemented
in MATLAB , read the help material for how to use the built-in ODE solvers. MATLAB
implementation of ODE solvers are discussed in Spencer and Ware (2008), sec. 16.5. Typi-
cally, you want to try ode45, and if that fails, ode113, or ode155.
We previously discussed the 4th order Runge Kutta method as a simple method to
solve initial value problems where the task is to forward integrate a vector y(t) from an
initial condition y0 (t = t0 ) to some time t f while the time derivatives of y are given by
dy
(t) = f (t, y, C ) (3.24)
dt
we made the dependence of f on constant parameters explicit in the C.
Numerically, this is done by successively computing yn+1 for time t + h from the last
known solution for yn at time t with time step h
∇·v = 0 (3.26)
∂T
+ v · ∇ T = κ ∇2 T (3.27)
∂t
∂v 1 ρ
+ (v · ∇) v = ν∇2 v − ∇ P + g. (3.28)
∂t ρ0 ρ0
∂T 0
0
+ v 0 · ∇ T 0 = ∇2 T 0 (3.29)
∂t
∂T 0
1 ∂ω 0 2
+ v · ∇ ω = ∇ ω − Ra , (3.30)
Pr ∂t0 ∂x 0
where the primed quantities are now non-dimensionalized. If chosen right, such a trans-
formation makes the individual terms of the equations that depend on the primed vari-
ables of order unity, such that the material parameter dependent factors (here Ra and Pr)
measure the importance of the terms during solution (see sec. 2.4).
Eq. (3.30) is eq. (3.28) rewritten in terms of vorticity ω, which is defined as
ω = ∇×v (3.31)
∂w ∂u
ω= − . (3.32)
∂x ∂z
for convective cells with wavelength 2/a. This is an example of a spectral method where
spatial variations in properties such as T are dealt with in the frequency domain, here
with one harmonic. Such an analysis is also common when examining barely super-
critical convective instabilities.
3.3.3 Problems
The resulting equations for the time dependent parameters of the approximate Lorenz
convection equations are
dW
= Pr ( T1 − W ) (3.36)
dt
dT1
= −WT2 + rW − T1
dt
dT2
= WT1 − bT2
dt
where b = 4/(1 + a2 ), r = Ra/Rac with the critical Rayleigh number Rac .
Hint: You can code this any way you want, but consider the following (Figure 3.4):
• You will want to separate a “driver” that deals with initial conditions, control-
ling the total time steps, plotting, etc., and an actual routine that computes the
Runge Kutta step following the formula we discussed in class. Those should
be separate m-files, or at least separate functions.
50
45
40
35
30
25
2
T
20
15
10
0
30
20 20
10 15
10
0 5
0
−10 −5
−10
−20 −15
T1
W
Figure 3.3: Solution to one of the problem set questions visualizing the behavior of the Lorenz
equations (the Lorenz attractor).
• You will want to make the Runge Kutta stepper independent of the actual func-
tion that is needed to compute dy/dt so that you can reuse it for other prob-
lems later. This can be done in MATLAB by defining a function myfunc that
computes the derivatives, and then passing the function name myfunc as an ar-
gument to the Runge Kutta time stepper. Within the time stepper, the function
then then has to be referred to as
@func. Alternatively, the function that computes the derivatives can be made
into its own “.m” file, in the same directory as the other subroutines, making it
available to all subroutines in that folder.
• If you need some inspiration on how to do this, download the m-file fragments
we provide for this sections problem set, lorenz_dy.m, lorenz.m, and rkstep.
m. There is also a complete Python implementation, ode.py and lorenz.py, if
you are curious.
2. Use initial condition y0 = {0, 0.5, 0.5}, parameters b = 8/3, Pr = 10 and solve
for time evolution for all three variables from t = 0 to t = 50, using a time step
h = 0.005. Use r = 2 and plot T1 and T2 against time. Comment on the temporal
character of the solution, what does it correspond to physically?
3. Change r to 10, and then 24. Plot T1 and T2 against time, and also plot the “phase
space trajectory” of the system by plotting y in W, T1 , and T2 space using MATLAB
plot3. Comment on how these solutions differ.
4. Increase r to 25 and plot both time behavior of T and the phase space trajectory.
What happened? Compare the r = 25 solution with the r = 24 solution from the
last question. Do you think r = 24 will remain steady for all times? Why? Why not?
5. Use r = 30 and show on one plot how T1 evolves with time for two different
initial conditions, the y0 from before, {0, 0.5, 0.5}, and a second initial condition
{0, 0.5, 0.50001}. Comment.
6. Compare your solution with h = 0.005 for T1 and an initial condition of your choice
in the r = 30 regime with the MATLAB -internal ODE solver you deem most appro-
priate. Plot the absolute difference of the solutions against time. Comment.
For help with making simple plots with MATLAB , see Spencer and Ware (2008), for
example. It is very easy to get such plots while developing code and debugging, but
often hard to generate publication quality results from MATLAB.How to do this is dis-
cussed in numerous online resources, including some helpful routines at http://geoweb.
princeton.edu/people/simons/software.html, for example. While it is generally pre-
ferred to remain within one framework, you might want to consider plotting MATLAB
generated results with external graphics packages such as Gnuplot2 (mainly for x-y type
plots, but some great extra capability) or GMT3 (a very versatile plotting software that
does, however, often require scripting).
dx
= −mx + yz (3.37)
dt
dy
= −my + (z − a) x
dt
dz
= 1 − xy
dt
2 http://www.gnuplot.info/
3 Wessel
and Smith (1998), http://gmt.soest.hawaii.edu/, also see Becker and Braun (1998), http://
geodynamics.usc.edu/~becker/igmt
% initial values
y= [...];
time =0;tstop=50;
h = 0.005;% timestep
save_each = 1;
nstep=0;save_step=0;
while(time < tstop) % loop while time is smaller than tstop
if(mod(nstep,save_each)==0) % only save every save_each step
save_step=save_step+1;
ysave(save_step,:)=y;
...
end
% advance the y(1:3) solution by one 4th order Runge Kutta step
y = y + rkstep(....);
nstep=nstep+1;
time=time+h;
end
function dy = rkstep(.... )
%
%
% perform one 4th order Runge Kutta timestep and return
% the increment on y(t_n) by evaluating func(time,y,parameters)
%
% ... parts need to be filled in
%
%
% input values:
% h: time step
% t: time
% y: vector with variables at time = t which are to be advanced
% func: function which computes dy/dt
% parameters: structure with any parameters the func function might need
% save computations
h2=h/2;
k1 = h .* dydt(...);
k2 = h .* dydt(...);
....
% return the y_{n+1} timestep
dy = ....
Figure 3.4: Suggested program structure for the Lorenz equation ODE solver exercise. Available
online as lorenz.m, rkstep.m, and dydt.m.
3. Examples from our own research where we have used simple ODE solutions, in-
clude some work on parameterized convection (Loyd et al., 2007), a method that
goes back at least to Schubert et al. (1980), see Christensen (1985). In this case, the
box is the mantle, and the total heat content of the mantle, as parameterized by the
mean temperature, is the property one solves for. I.e. we are averaging the PDEs
governing convection spatially, to solve for the time-evolution of average mantle
temperature.
The idea is that a convective system with Rayleigh number (see eq. (2.44))
ραTgh3
Ra = (3.38)
κη
transports heat at Nusselt number Nu following a scaling of
Q
Nu = = aRa β (3.39)
cT
(with some debate about β, see, e.g. Korenaga, 2008, for a review). The energy balance
for the mantle is
dT
Cp = H (t) − Q (3.40)
dt
where C p is the total heat capacity, H the time-dependent heat production through
radiogenic elements, and Q the heat loss through the surface. If viscosity is a func-
tion of temperature,
H
η = η0 exp (3.41)
RT
then the equations couple such that
1+ β
η ( T0 ) β
dT T
Cp = H ( t ) − Q0 . (3.42)
dt T0 η (T )
Figure 3.5: Poincare sections in y for the period doubling sequence to chaos for the spring-slider
system, eq. (3.43), as a function of normalized spring stiffness, κ 0 . Bottom figure shows zoom into
the dashed rectangular region highlighted on top (modified from Becker, 2000).
54
Chapter 4
Finite differences
∂ f ( x0 )
f ( x0 + h ) = f ( x0 ) + h + (4.1)
∂x
∂2 f ( x0 ) h2 ∂3 f ( x0 ) h3 ∂ n f ( x0 ) h n
+ + . . . + +
∂x2 2! ∂x3 3! ∂x n n!
O(hn+1 ).
∂ f ( x0 ) f ( x0 + h ) − f ( x0 ) ∂2 f ( x0 ) h ∂3 f ( x0 ) h2
= − − ... (4.2)
∂x h ∂x2 2! ∂x3 3!
If we now only compute the first term of this equation as an approximation, we can write
a discretized version
∂ f ( xi ) f − fi
= i +1 + O(h2 ) (4.3)
∂x h
55
CHAPTER 4. FINITE DIFFERENCES
∂ f ( x0 ) ∂2 f ( x0 ) h2 ∂3 f ( x0 ) h3
f ( x0 − h ) = f ( x0 ) − h+ − + ... (4.4)
∂x ∂x2 2! ∂x3 3!
In this case, the first, backward difference can be obtained by
∂ f ( xi ) f − f i −1
= i + O(h2 ). (4.5)
∂x h
Proceeding in a similar fashion, we can derive higher order derivatives. Introducing
the abbreviations
∂f ∂2 f
f0 = , f 00 = 2 . . . (4.6)
∂x ∂x
we can find, for example,
f i0+1 − f i0
f i00+1 = + O(h2 ) (4.7)
h
f i +2 − f i +1
− f i+1h− f i
= h
+ O(h2 ) (4.8)
h
f i +2 − 2 f i +1 + f i
= + O(h2 ) (4.9)
h2
which is the first order accurate, forward difference approximation for second order deriva-
tives around f i+1 .
If we wish to improve on accuracy, we can proceed by taking higher order terms of
the Taylor expansion, and using first order accurate estimates for the derivatives. For
example,
f ( x + h) − f ( x ) h 00
f 0 (x) = − f (x) + . . . (4.10)
h 2
f ( x + h) − f ( x )
= −
h
h f ( x + 2h) − 2 f ( x + h) + f ( x )
+ O(h ) + O(h3 )
2
or
2 h2
− f i +2 + 4 f i +1 − 3 f i
f i0+1 = + O(h3 ). (4.11)
2h
Alternatively, we can form the average of the first order accurate forward and back-
ward schemes, i.e. adding eqs. (4.3) and (4.5) and dividing by two. The result is the central
difference approximation of the first derivative
f i +1 − f i −1
f i0 = + O(h3 ) (4.12)
2h
and also second order accurate.
Note that eq. (4.12) involves fewer function evaluations than eq. (4.11), which is why
eq. (4.12) is preferred for actual implementations. Also, both equations now require
knowledge of f over three lateral grid points (three point stencil), rather than two as
was needed for first order accuracy. Moreover, improving accuracy by taking higher and
higher order polynomial expansions into account only works if the function f is actually
smooth (i.e. differentiable) in that way, and not “weird” or jumpy. In this case, lower order
approximations may do just as well.
By adding eqs. (4.1) and (4.4) an approximation of the second derivative is obtained
f i +1 − 2 f i + f i −1
f i00 = + O(h3 ). (4.13)
h2
A different way to derive the second derivative of second order accuracy is by com-
puting, by central differences which we now know should yield second accuracy, the first
derivatives (theoretically!) as evaluated at the points in between, i + 1/2 and i − 1/2,
and then computing the second derivative at i from this by using the central difference of
those two first derivatives:
f i +1 − f i
f i0+1/2 =
h
f i − f i −1
f i0−1/2 =
h
f i +1 − f i f i − f i −1
f i+1/2 − f i0−1/2
0
− f i +1 − 2 f i + f i −1
f i00 = = h h
= (4.14)
h h h2
(This is important for derivatives with variable coefficients, cf. sec. 4.1.3.) If h is not con-
stant, it should be computed as xi+1 − xi . The second order derivative in this case is
identical to the equations above with xi+1 − xi = xi − xi−1 = h.
Similarly, we can derive higher order derivatives, and higher order accuracy (but only
if f is of polynomial form). Both require more input values, a larger stencil. A general
approach to forming interpolations of f and dn f /dx n can be found in Fornberg (1996).
Note that the highest order derivative that usually occurs in geodynamics is the 4th -order
derivative.
where k i−1/2 is evaluated between ui and ui−1 , to maintain the second order accuracy of
the central difference approach for second derivatives (eq. 4.14)
00f 0 ( x + ∆x 0 h
2 ) − f (x − 2 )
f = , (4.26)
h
i.e. the first derivatives are premultiplied with the cofficients in between the nodes. If k
is spatially varying, the following, common approximations are therefore inadequate to
maintain second order accuracy:
u −u u −u
k i+1 i+1h i − k i i h i−1
∂ ∂u
k =
∂x ∂x i h
∂T ∂2 T
=κ 2 (4.28)
∂t ∂x
where
k
κ= (4.29)
ρc p
is the thermal diffusivity (a common value for rocks is κ = 10−6 m2 s−1 ; also see discussion
in sec. 2.4).
We are interested in the temperature evolution versus time, T ( x, t), which satisfies
eq. (4.28), given an initial temperature distribution (Fig. 4.1A). An example would be the
intrusion of a basaltic dike in cooler country rocks. How long does it take to cool the
dike to a certain temperature? What is the maximum temperature that the country rock
experiences?
The first step in the finite differences method is to construct a grid with points on
which we are interested in solving the equation (this is called discretization, see Fig. 4.1B).
The next step is to replace the continuous derivatives of eq. (4.28) with their finite differ-
ence approximations. The derivative of temperature versus time ∂T ∂t can be approximated
with a forward finite difference approximation as
i,n+1
T(x,0) i,n-1
time
Dt
Dx
L
x
space
Figure 4.1: A) Setup of the thermal cooling model considered here. A hot basaltic dike intrudes
cooler country rocks. Only variations in x-direction are considered; properties in the other di-
rections are assumed to be constant. The initial temperature distribution T ( x, 0) has a step-like
perturbation, centered around the origin with [−W/2; W/2] B) Finite difference discretization of
the 1D heat equation. The finite difference method approximates the temperature at given grid
points, with spacing ∆x. The time-evolution is also computed at given times with time step ∆t.
The last step is to specify the initial and the boundary conditions. If for example the
country rock has a temperature of 300◦ C and the dike a total width W = 5 m, with a
magma temperature of 1200◦ C, we can write as initial conditions:
In addition we assume that the temperature far away from the dike center (at | L/2|) re-
mains at a constant temperature. The boundary conditions are thus
The MATLAB code in Figure 4.2, heat1Dexplicit.m, shows an example in which the
grid is initialized, and a time loop is performed. In the exercise, you will fill in the ques-
tion marks and obtain a working code that solves eq. (4.33).
4.2.1 Exercises
1. Open MATLAB and an editor and type the MATLAB script in an empty file; alter-
natively use the template provided on the web if you need inspiration. Save the file
under the name heat1Dexplicit.m. If starting from the template, fill in the question
marks and then run the file by typing heat1Dexplicit in the MATLAB command
window (make sure you’re in the correct directory). (Alternatively, type F5 to run
from within the editor.)
2. Study the time evolution of the spatial solution using a variable y-axis that adjusts
to the peak temperature, and a fixed axis with range axis([-L/2 L/2 0 Tmagma]).
Comment on the nature of the solution. What parameter determines the relationship
between two spatial solutions at different times?
Does the temperature of the country rock matter for the nature of the solution? What
about if there is a background gradient in temperature such that the country rock
temperature increases from 300◦ at x = − L/2 to 600◦ at x = L/2?
3. Vary the parameters (e.g. use more grid points, a larger or smaller time step). Com-
pare the results for small ∆x and ∆t with those for larger ∆x and ∆t. How are these
solutions different? Why? Notice also that if the time step is increased beyond a
certain value, the numerical method becomes unstable and does not converge – it
grows without bounds and exhibits non-physical features.
Investigate which parameters affect stability, and find out what ratio of these pa-
rameters delimits this scheme’s stability region. This is called the CFL condition,
see von Neumann stability analysis in (cf. chap 5 of Spiegelman, 2004).
%heat1Dexplicit.m
%
% Solves the 1D heat equation with an explicit finite difference scheme
clear
%Physical parameters
L = 100; % Length of modeled domain [m]
Tmagma = 1200; % Temperature of magma [C]
Trock = 300; % Temperature of country rock [C]
kappa = 1e-6; % Thermal diffusivity of rock [m2/s]
W = 5; % Width of dike [m]
day = 3600*24; % # seconds per day
dt = 1*day; % Timestep [s]
% Numerical parameters
nx = 201; % Number of gridpoints in x-direction
nt = 500; % Number of timesteps to compute
dx = L/(nx-1); % Spacing of grid
x = -L/2:dx:L/2;% Grid
time = 0;
for n=1:nt % Timestep loop
% Plot solution
figure(1), clf
plot(x,Tnew);
xlabel(’x [m]’)
ylabel(’Temperature [^oC]’)
title([’Temperature evolution after ’,num2str(time/day),’ days’])
drawnow
end
Figure 4.2: MATLAB script heat1Dexplicit.m to solve eq. (4.28) (once the blanks indicated by the
questions marks are filled in . . . ).
4. Record and plot the temperature evolution versus time at a distance of 5 m from
the dikecountry rock contact. What is the maximum temperature the country rock
experiences at this location and when is it reached? Assume that the country rock
was composed of shales, and that those shales were transformed to hornfels above
a temperature of 600◦ C. What is the width of the metamorphic aureole?
5. Think about how one would write a non-dimensionalized version of the tempera-
ture solver.
6. Add a test with an analytical solution for diffusion and plot error vrs. resolution. A
good reference for analytical solutions for heat conduction problems is Carslaw and
Jaeger (1959), or see sec. 4.7.4.
The spatial discretization should be second order for a second order scheme.
7. Derive a finite-difference approximation for variable k (and variable ∆x allowing for
uneven spacing between grid points should you so desire). Test the solution for the
case of k = 10 inside the dike, and k = 3 in the country rock.
Limited stability and numerical aliasing/dissipation are two major drawbacks of ex-
plicit finite difference codes such as the one presented in sec. 4.2. Next, we will discuss
methods that do not have these limitations.
∂T ∂2 T
=κ 2 (4.38)
∂t ∂x
if material parameters are homogeneous.
In explicit finite difference schemes, the temperature at time n + 1 depends only on
the already known temperature at time n. The explicit finite difference discretization of
eq. (4.38) is
Tin+1 − Tin T n − 2Tin + Tin−1
= κ i +1 , (4.39)
∆t (∆x )2
using central differences for the spatial derivatives (subscript i indicating the x location
in 1-D, superscripts indicating the time). Eq. (4.39) can be rearranged in the following
manner (with all quantities at time n + 1 on the left and quantities at time n on the right-
hand-side)
n +1 n
Tin+1 − 2Tin + Tin−1
Ti = Ti + κ∆t (4.40)
(∆x )2
Since we know Tin+1 , Tin and Tin−1 , we can compute Tin+1 . This is schematically shown on
Figure 4.3a, and an algorithm based on eq. (4.40), such as the one of last section’s problem
set, is called a forward time, centered space (FTCS) because of the way it is computed.
The major advantage of explicit finite difference methods is that they are relatively
simple, only one solution for T needs to be stored, and the method is computationally
Figure 4.3: A) Explicit finite difference discretization. B) Implicit finite difference discretization.
C) Crank-Nicolson finite difference discretization.
fast for each time step. However, the main drawback is that stable solutions are obtained
only when
2κ∆t (∆x )2
0< ≤ 1 or ∆t ≤ for given ∆x. (4.41)
(∆x )2 2κ
If this condition is not satisfied, the solution becomes unstable, starts to wildly oscillate,
or “blow up”.
This can be shown by von Neumann stability analysis (Press et al., 1993, chap. 19.2),
analyzing the growth of of eigenmodes of the finite difference equation. However, phys-
ically, the stability condition eq. (4.41) means that the maximum time step needs to be
smaller than the time it takes for an anomaly to diffuse across the grid (nodal) spacing ∆x
(cf. diffusion time in sec. 3.3). The explicit solution, eq. (4.40), is an example of a condition-
ally stable method that only leads to well behaved solutions if a criterion like eq. (4.41) is
satisfied.
Note that eq. (4.41) can only hold for κ∆t > 0; having a negative diffusivity, or using
a time-reversed ∆t < 0, will invariably lead to blow up since small features will get em-
phasized rather than smoothed out. This is an issue if one wishes to reconstruct diffusive-
advective processes (such as mantle convection), going from the present-day temperature
field back in time (cf. Ismail-Zadeh and Tackley, 2010).
We will revisit an FTCS scheme similar to eq. (4.39) for advection that involves sin-
gle derivatives in space later. Unlike eq. (4.39), the FTCS scheme for advection is always
unstable. Not a good idea.
Even if the FTCS for diffusion can be made stable, the stability condition leads to
numerical convenience issues. Given that we are typically interested in spatial features
with wavelength, λ, within the solution that are much larger than ∆x, λ ∆x, because
we want to resolve the solution features at least with a few nodes, the explicit scheme,
eq. (4.39), will require
2
λ
1 (4.42)
∆x
steps per relevant time scale for the evolution of λ features, which is usually prohibitive.
An alternative approach is an implicit finite difference scheme, where the spatial deriva-
tives ∂2 T/∂x2 are evaluated (at least partially) at the new time step. The simplest implicit
discretization of eq. (4.38) is
Therefore, the fully implicit scheme will always yield the right equilibrium solution but
may not capture small scale, transient features.
It turns out that the fully implicit method described by eq. (4.43) is second order ac-
curate in space but only first order accurate in time, i.e. the error goes as O((∆x )3 , ∆t2 ).
It is possible to write down a scheme which is second order accurate both in time and
in space (i.e. O((∆x )3 , (∆t)3 )). One such scheme is the Crank-Nicolson scheme (see ex-
ercises, Fig. 4.3C), which is unconditionally stable. Note the analogy with the previous
derivation of spatial derivatives: forward or backward differences were only first order
accurate, while the central difference approach achieved second order accuracy O((∆x )3 ).
The Crank-Nicolson method is the time analog of central spatial differences. However,
any partially implicit method is more tricky to compute as we need to infer the future
solution at time n + 1 by solution (inversion) of a system of linear equations based on the
known solution at time n. This is discussed next.
T1 = Tle f t (4.55)
Ax = b. (4.57)
T1n+1
T2n+1
T3n+1
x= , (4.59)
T4n+1
T5n+1
T6n+1
Note that matrix A will have a unity entry on the diagonal and zero else for each node
where Dirichlet (fixed temperature) boundary conditions apply; see derivation below and
eqs. (4.73) and (4.74) for how to implement Neumann boundary conditions.
Matrix A also has an overall peculiar form because most entries off the diagonal are
zero. This “sparseness” can be exploited by specialized linear algebra routines, both in
terms of storage and speed. By avoiding computations involving zero entries of the ma-
trix, much larger problems can be handled than would be possible if we were to store
the full matrix. In particular, the fully implicit FD scheme leads to a “tridiagonal” system
of linear equations that can be solved efficiently by LU decomposition using the Thomas
algorithm (e.g. Press et al., 1993, sec. 2.4).
A = sparse(nx,nx);
for i=2:nx-1
A(i,i-1) = -s;
A(i,i ) = (1+2*s);
A(i,i+1) = -s;
end
A(1 ,1 ) = 1;
A(nx,nx) = 1;
(Exercise: Improve on the loop formulation for A assembly by using MATLAB vector
functionality.)
Once the coefficient matrix has been constructed, its structure can be visualized with
the command
>>spy(A)
(Try it, for example by putting a “break-point” into the MATLAB code below after assem-
bly.)
The right-hand-side vector b can be constructed with
b = zeros(nx,1);
b(2:nx-1) = Told(2:nx-1);
b(1) = Tleft; b(nx) = Tright;
The only thing that remains to be done is to solve the system of equations and find x.
MATLAB does this with
x = A\b;
The vector x is now filled with new temperatures T n+1 , and we can go to the next time
step. Note that, for constant ∆t, κ, and ∆x, the matrix A does not change with time. There-
fore we have to form it only once in the program, which speeds up the code significantly.
Only the vectors b and x need to be recomputed. (Note: Having a constant matrix helps a
lot for large systems because operations such as x = A\b can then be optimized further
by storing A in a special form.)
4.4.4 Exercises
1. Save the script heat1Dexplicit.m from last section as heat1Dimplicit.m. Program
the implicit finite difference scheme explained above. Compare the results with
results from last section’s explicit code.
2. Time-dependent, analytical solutions for the heat equation exists. For example, if
the initial temperature distribution (initial condition, IC) is
x 2
T ( x, t = 0) = Tmax exp − (4.61)
σ
where Tmax is the maximum amplitude of the temperature perturbation at x = 0 and
σ its half-width of the perturbance (use σ < L, for example σ = W). The solution is
then
− x2
Tmax
T ( x, t) = √ exp (4.62)
1 + 4tκ/σ2 σ2 + 4tκ
(for T = 0 BCs at infinity). (See Carslaw and Jaeger, 1959, for useful analytical solu-
tions to heat conduction problems).
Program the analytical solution and compare the analytical solution with the nu-
merical solution with the same initial condition. Compare results of the implicit
and FTCS scheme used last section to the analytical solution near the instability
region of FTCS,
κ∆t 1
s= 2
< . (4.63)
(∆x ) 2
√
Note: Eq. (4.62) can be derived using a similarity variable, x̃ = x/xc with xc ∝ κt.
Looks familiar?
(b) Write down a finite difference discretization of ∂2 T/∂x2 = 0 and solve it. (See
the limit case consideration above.)
Employ both methods to compute steady-state temperatures for Tle f t = 100◦ and
Tright = 1000◦ . Derive the analytical solution and compare your numerical solu-
tions’ accuracies. Use the implicit method for part (a), and think about different
boundary conditions, and the case with heat production.
4. Apply no flux boundary conditions at | x | = L/2 and solve the dike intrusion prob-
lem in a fully implicit scheme. Eqs. (4.73) and (4.74) need to replace the first and last
columns of your A matrix.
5. Derive and program the Crank-Nicolson method (cf. Figure 4.3C). This “best of both
worlds” method is obtained by computing the average of the fully implicit and fully
explicit schemes:
n + 1 n + 1 n + 1 n − 2T n + T n
Tin +1 n
− Ti κ Ti +1 − 2Ti + Ti −1 + Ti +1 i i −1
= 2
. (4.64)
∆t 2 (∆x )
This scheme should generally yield the best performance for any diffusion problem,
it is second order time and space accurate, because the averaging of fully explicit
and fully implicit methods to obtain the time derivative corresponds to evaluating
the derivative centered on n + 1/2. Such centered evaluation also lead to second
order accuracy for the spatial derivatives.
Compare the accuracy of the Crank-Nicolson scheme with that of the FTCS and fully
implicit schemes for the cases explored in the two previous problems, and for ideal
values of ∆t and ∆x, and for large values of ∆t that are near the instability region of
FTCS.
Hint: Proceed by writing out eq. (4.64) and sorting terms into those that depend on
the solution at time step n + 1 and those at time step n, as for eq. (4.44).
6. Bonus question: Write a codefor the thermal equation with variable thermal con-
ductivity k: ρc p ∂T ∂ ∂T
∂t = ∂x k ∂x . Assume that the grid spacing ∆x is constant. For
extra bonus, allow for variable grid spacing and variable conductivity.
Figure 4.4: Discretization of the numerical domain with fictitious boundary points, that are em-
ployed to set flux boundary conditions.
i.e. we are simply extrapolating from T2 to T0 with the slope given by c1 . Substituting
into eq. (4.69) yields
2. To apply this general formulation in a fully implicit scheme, we can again rearrange
terms from eq. eq. (4.71) to bring known quantities on the right-hand-side:
These equations now only involve grid points that are part of the computational
grid, and eqs. (4.73) and (4.74) can be incorporated into the matrix A and right-
hand-side b (compare with eq. 4.44).
4.6.1 Example
We consider a case of fluid flow in a porous media (governed by the Darcy equation)
whose diffusivity κ is a function of the fluid pressure (high fluid pressure increases per-
meability, p ↑→ κ ↑). In a 1-D, vertical (z) column the governing equation shall be
∂P ∂ ∂P
= κ ( P) (4.75)
∂t ∂z ∂z
where P is the fluid pressure, and κ ( P) the hydraulic diffusivity. The equation is nonlinear
because the diffusivity can be written as a function of the fluid pressure P, which is related
to the effect of dilation and cracking under enhanced fluid pressure.
To solve eq. (4.75), we need a constitutive law, and we assume that the hydraulic dif-
fusivity is given by
κ ( P) = κ0 + cPm (4.76)
where κ0 is the background diffusivity, and c and m (semi-empirical) constants.
We will use a fully implicit scheme, so that the discretization is done in analogy (sec-
ond order accurate second spatial derivative) to the implicit 1-D thermal diffusion prob-
lem:
P n +1 − P n +1 P n +1 − P n +1
Pin+1 − Pin κin++1/2
1 i +1
∆x
i
− κin−+1/2
1 i
∆x
i −1
= (4.77)
∆t ∆x
where the material parameters are evaluated in between nodes, for example by comput-
ing an average
κin+1 + κin±+11
κin±+1/2
1
= . (4.78)
2
(If diffusivity were merely heterogeneous (such as in the previous explicit heat equa-
tion example), but not dependent on the solution itself, we could use a “staggered” grid
where κ would be specified at nodes located in between the locations where P is to be
computed. )
The implicit equations can again be solved in matrix form as
AP = b (4.79)
for Pn+1 . The problem, however, is that κ depends on Pn+1 . Therefore we have to perform
iterations for the true Pn+1 and recompute A at each time step before advancing time. The
general recipe is
4. Use this new pressure estimate to recalculate diffusivities and the coefficients in A.
5. Return to step 2 and continue until the pressure Pn+1 stops changing, which in-
dicates that the solution has converged. Use as an indication of convergence the
following error estimate:
If convergence is reached (e.g. relative error < 10−4 ), continue to the next time step.
Exercise
• Write a program that solves the equations described above on the domain z ∈ [0; 1]
from t = 0 to t = 0.2. Assume that we have zero flux boundary conditions (i.e.
gradient ∂P/∂z = 0 on top and bottom). Use non-dimensional parameter values
κ0 = 0.05, c = 1, m = 2 and time-step 0.005. The initial pressure is to be unity
within [0.4; 0.6] and zero else. At each time step, compare the nonlinear solution to
the linear one, obtained by setting c = 0, to visualize the effect of the non-linear
terms.
i+1,j
i-1,j
Dz
Dx
L
x
where, ρ is density, c p heat capacity, k x,z the thermal conductivities in x and z direction,
and Q radiogenic heat production.
If the thermal conductivity is isostropic (k x = k z ) and constant, we can rewrite
∂2 T ∂2 T
∂T Q
=κ + 2 + . (4.82)
∂t ∂x2 ∂z ρc p
where ∆x and ∆z indicates the node spacing in both spatial directions, and there are now
two indices for space, i and j for zi and x j , respectively (Figure 4.5). Rearranging gives
n ∆t
Qi,j
n +1 n n n n
Ti,j = Ti,j + sx Ti,j +1 − 2Ti,j + Ti,j −1 + sz Tin+1,j n
− 2Ti,j + Tin−1,j + , (4.84)
ρc p
where
κ∆t κ∆t
sx = and sz = . (4.85)
(∆x )2 (∆z)2
Boundary conditions can be set the usual way. A constant (Dirichlet) temperature on
the left-hand side of the domain (at j = 1), for example, is given by
A constant flux (Neumann BC) on the same boundary at {i, j = 1} is set through fictitious
boundary points
∂T
= c1 (4.87)
∂x
Ti,2 − Ti,0
= c1
2∆x
Ti,0 = Ti,2 − 2∆xc1
which can then be plugged into eq. (4.84) so that for j = 1, for example,
n +1 n n n
Ti,1 = Ti,1 + s x 2Ti,2 − 2( Ti,1 + ∆xc1 )
n ∆t
Qi,1
sz Tin+1,1 n
Tin−1,1
+ − 2Ti,1 + + (4.88)
ρc p
2κ∆t
≤ 1. (4.90)
min((∆x )2 , (∆z)2 )
nz 29 30 31 32 33 34 35
22 23 24 25 26 27 28
15 16 17 18 19 20 21
z
8 9 10 11 12 13 14
1 2 3 4 5 6 7
nx
If n x are the number of grid points in x-direction and nz the number of points in z-
direction, we can write eqs. (4.93) and (4.94) in a more general way as:
∂2 T 1
|i,j = T − 2T(i−1)nx + j + T(i−1)nx + j−1 (4.95)
∂x2 (∆x )2 (i−1)nx + j+1
∂2 T 1
| = T − 2T + T( i −2) n x + j . (4.96)
∂z2
i,j
(∆z)2 i·nx + j ( i −1) n x + j
T1n+1 = T1,1
T2n+1 = T1,2
.
.
.
T n +1
( i −1) n x + j
= Ti,j
x = n +1 , (4.98)
T(i−1)n + j+1 = Ti,j+1
x
..
.
n + 1
Tn n −1 = Tnz ,nx −1
x z
Tnnx+n1z = Tnz ,nx
and the load (right hand side) vector is given by (Q = 0 for simplicity)
Tbottom
Tbottom
..
.
T(ni−1)n + j
x
b= n . (4.99)
T
( i −1) n x + j +1
..
.
Ttop
Ttop
1. Fill in the question marks in the script heat2Dexplicit.m (Figure 4.7), by program-
ming the explicit finite difference scheme. Employ zero flux boundary conditions
∂T
∂x = 0 on the left and on the right-side of the domain (outside the top and bot-
tom edges), and constant temperature conditions on the top and bottom. Ignore the
effects of radioactive heat.
2. Finish the code heat2Dimplicit.m (Figure 4.8), by programming the implicit finite
difference approximation of the 2D temperature equation.
% Numerical parameters
nx = 101; % # gridpoints in x-direction
nz = 51; % # gridpoints in z-direction
nt = 500; % Number of timesteps to compute
dx = L/(nx-1); % Spacing of grid in x-direction
dz = H/(nz-1); % Spacing of grid in z-direction
[x2d,z2d] = meshgrid(-L/2:dx:L/2, -H:dz:0); % create grid
% Compute stable timestep
dt = min([dx,dz])^2/kappa/4;
% Setup initial linear temperature profile
T = abs(z2d./H)*Tbot;
% Imping plume beneath lithosphere
ind = find(abs(x2d(1,:)) <= Wplume/2);
T(1,ind) = Tplume;
time = 0;
for n=1:nt
Figure 4.7: MATLAB script heat2D_explicit.m to solve the 2D heat equation using the explicit
approach.
−( x2 + z2 )
T ( x, z, t = 0) = Tmax exp (4.100)
σ2
Figure 4.8: MATLAB script heat2D_explicit.m to solve the 2D heat equation using the implicit
approach.
USC GEOL557: Modeling Earth Systems 82
CHAPTER 4. FINITE DIFFERENCES
−( x2 + z2 )
Tmax
T ( x, z, t) = √ exp . (4.101)
1 + 4tκ/σ2 σ2 + 4tκ
(As for the 1D√example, note that this uses a characteristic, time-dependent length-
scale of lc ∝ κt, as expected for a diffusion problem (cf. sec. 2.4), and see Carslaw
and Jaeger (1959) for more analytical solutions.)
Program the analytical solution and compare it with the explicit and fully implicit
numerical solutions with the same initial conditions at each time step. Comment on
the accuracy of both methods for different values of dt.
4. Add the effects of radioactive heat to the explicit/implicit equations above. Use
Turcotte and Schubert (2002) or Google to find typical values of Q, ρ, c p for rocks.
5. Bonus: Write a code for the thermal equation with variable thermal conductivity k.
Assume that the grid spacing ∆x is constant. This type of code is not only relevant
for thermal problems, but also for problems like hydro-geological problems (Darcy
flow, e.g. how far did the chemical waste go into the aquifer?), fluid movements
through the crust and through fault zones (which is related to the creation of ore de-
posits), magma migration through the mantle, geochemistry and mineral reactions
at grain-boundary scale, and, aftershocks and fluids.
Figure 4.9: Snapshots of a bottom heated thermal convection model with a Rayleigh-number of
5 × 105 and constant viscosity (no internal heating). Temperature is advected through a fixed
(Eulerian) grid (circles) with a velocity (arrows) that is computed with a Stokes solver.
Heat sources would lead to additional terms on the right hand side. Since temperature
variations lead to buoyancy forces, the energy equation is coupled with the Stokes (con-
servation of momentum) equation from which velocities v can be computed to close the
system needed for a convection algorithm.
Mantle convection codes typically deal with advection of a temperature field assum-
ing that there is significant diffusion at the same time, κ > 0, and will at times produce
non-physical artifacts (e.g. temperature ringing, overshoot “waves” or “halos” around
advected strong temperature contrast) in cases that are advection-dominated.
One example would be if a chemical composition C is to be treated akin to T with a
typical field method,
∂C
+ v · ∇C = κc ∇2 C. (4.106)
∂t
Chemical diffusivities are for mantle purposes zero (i.e. compositional mixing happens
by stirring, not molecular diffusion), κc ≈ 0, and special tricks are required to use field
methods to solve
∂C
+ v · ∇C = 0 (4.107)
∂t
(e.g. Lenardic and Kaula, 1993; Kronbichler et al., 2012), as discussed below (cf. sec. 7.05.5.1
of Zhong et al., 2007).
f 2 −1
f1
ηeff = + , (4.109)
η1 η2
and at the high (“strong”) end by the artihmetic mean viscosity,
ηeff = f 1 η1 + f 2 η2 . (4.110)
Those two correspond to the horizontal simple and pure shear deformation of a horizon-
tally layered η1 and η2 sandwich, respectively. The 1-D rheological element models are
Figure 4.10: Exploration of different tracer averaging and processing schemes (from Tackley and
King, 2003). Figure shows composition, C, and temperature, T, of a simple, thermo-chemical
convection test akin to van Keken et al. (1997), for different methods of estimating composition.
(Colorbars for C and T go from 0 (blue) to 1 (red).) The absolute method uses Ci = ANi /Vi to
compute the composition of each unit volume, Vi , from the number of particles of the “special”
(e.g. dense) species, Ni , within the element. Here, A is a constant, and tracers are originally only
placed in the region where C = 1 (here, the dense lower layer). Truncated absolute indicates filtering
of unphysical values which may arise following Lenardic and Kaula (1993). The ratio method assigns
tracers everywhere and defines their original composition, c, as unity or zero depending on tracers
being with the material of interest or outside. Then, Ci = Ni (1)/ ( Ni (0) + Ni (1)), where the Ni ( x )
are the number of tracers within Vi of composition c = x. Number of tracers per element increases
with each row. Note that the ratio method is generally preferred and appears to allow for much
smaller number of tracers per element (Tackley and King, 2003).
dashpots in series (with constant stress and additive strain-rates, eq. 4.109) and dashpots
in parallel (with constant-strain rate and additive stress, eq. 4.110), respectively.
Since it is usually impractical to decide on the deformation state a priori, the interme-
diate case of the geometric mean viscosity, or, equivalenty the log-average viscosity,
f f
ηeff = η11 η22 (4.111)
log ηeff = f 1 log η1 + f 2 log η2 (4.112)
(in practice, use eq. 4.112) is preferred, as explored by Schmeling et al. (2008). Different
choices in averaging methods at an elemental level can have drastically different effects
for certain problems (Schmeling et al., 2008). The question of an appropriate average vis-
cosity also arises in the case of defining an average Rayleigh number, e.g. for temperature
and strain-rate dependent viscosity. In these cases, the log-average of eq. (4.112) is also
generally preferred (Christensen, 1984).
Tracer approaches have gotten more popular over the last decades since they combine
intuitively appealing aspects such as the transport of different material (e.g. crust vrs.
peridotite in thermo-chemically convecting mantle with fractionation) with natural im-
plementation of path dependence (“memory”, as required, e.g., for visco-elasticity) with
the benefits of allowing a finite difference or finite element mesh to remain fixed in an
Eulerian frame (“marker in cell” methods, e.g. Gerya and Yuen, 2003; Moresi et al., 2003).
However, the book keeping that is involved in efficiently conducting the operations
that are involved in the numerical implementation (e.g. “find all markers within this unit
cell”) and other necessary steps such as visualization are somewhat involved, particu-
larly on distributed, multi-processor approaches, which is why we do not discuss them at
length here. Existing implementations and libraries for standard computations should be
consulted. However, we will return to a hybrid field/tracer approach for advection (the
semi-Lagrangian scheme) below.
FTCS method
In 1-D, the simplest way to discretize eq. (4.113) is by employing a central difference
scheme in space, and go forward in time (another example of a forward-time, central
space, FTCS, scheme):
Tin+1 − Tin T n − Tin−1
= −v x,i i+1 , (4.115)
∆t 2∆x
where v x,i is the v x velocity at location i.
Exercise 1 We will consider a 1-D problem, with constant v x velocity in which an ex-
ponential pulse of temperature is getting advected along the x axis (see Figure 4.11 and
exercise_1_ftcs.m).
• Program the FTCS method in the code of Figure 4.11 and watch what happens.
• Change the sign of the velocity.
• Change the time step and grid spacing and compute the non-dimensional parameter
|v x |∆t/∆x.
• When do unstable results occur? Put differently, can you find a ∆t small enough to
avoid blow-up?
As you can see from the exercise, the FTCS method does not work so well . . . In fact, it is
a nice example of a scheme that looks logical on paper, but looks can be deceiving. The
FTCS method is unconditionally unstable, blows up for any ∆t, as can be shown by von
Neumann stability analysis (cf. chap 5 of Spiegelman, 2004). The instability is related to the
fact that this scheme produces negative diffusion, which is numerically unstable.
Lax method
The Lax approach consists of replacing the Tin in the time-derivative of eq. (4.115) with
( Tin+1 + Tin−1 )/2. The resulting equation is
Exercise 2
• Program the Lax method by modifying the script of the last exercise.
• Try different velocities and ∆t settings and compute the Courant number, α, which is
given by the following equation:
v x ∆t
α= (4.117)
∆x
%
% FTCS advection schem
%
clear all
nx = 201;
W = 40; % width of domain
Vel = -4; % velocity
sigma = 1;
Ampl = 2;
nt = 500; % number of timesteps
dt = 1e-2; % timestep
dx = W/(nx-1);
x = 0:dx:W;
% Initial Gaussian T-profile
xc = 20;
T = Ampl*exp(-(x-xc).^2/sigma^2);
% Velocity
Vx = ones(1,nx)*Vel;
abs(Vel)*dt/dx
cfac = dt/(2*dx);
% Central finite difference discretization
for itime=1:nt
% central fin. diff
for ix=2:nx-1
Tnew(ix) = ???
end
% BCs
Tnew(1) = ???
Tnew(nx) = ???
% Update Solution & time incremement
T = Tnew;
time = itime*dt;
% Analytical solution for this case
T_anal = ???
figure(1),clf, plot(x,T,x,T_anal), ...
legend(’Numerical’,’Analytical’)
xlabel(’x’)
ylabel(’temperature’)
drawnow
end
α2 ∂2 ∂ n α2 (∆x )2 ∂ n
n +1 n
Mx − ( T − Ti ) + α∆x T − T =0 (4.118)
(∆x )2 ∂x2 i ∂x i 2 ∂x i
This formulation gives us much better accuracy (O(∆t2 , (∆x )2 ) by using a higher or-
der discretization in both time and space. But what is its stability range in terms of
Figure 4.12: Illustration of the Courant criterion (from Press et al., 1993, chap 19.1).
Courant number? Notice the difference in terms of artificial diffusion, and oscilla-
tions with respect to the simple Lax method.
As you saw from exercise 2, the Lax method does not blow up, but does have a lot of
numerical diffusion for α 6= 1 (which is hard to attain for realistic problems, as v will vary
in space and time). In fact, the Lax criterion stabilized the discretized advection equation
by adding some artificial diffusion. So, it is an improvement but it is far from perfect, since
you may now lose the plumes of Figure 4.9 around mid-mantle purely due to numerical
diffusion. As for the case of the implicit versus explicit solution of the diffusion equation,
you see that there are trade-offs between stability and accuracy. There is no free lunch,
and numerical modeling is also a bit of an art.
The stability requirement
|V |∆t
α= ≤1 (4.120)
∆x
is called the Courant criterion (Figure 4.12).
Note that we have replaced central with forward or backward derivatives, depending on
the flow direction. The idea is that the flux into the local cell at xi will only depend on the
gradient of temperature in the direction “upstream”, i.e. where the inflowing velocity is
coming from.
Exercise 3
• Program the upwind scheme method.
• Try different velocity distributions (not just constant) and compute the Courant
numbers α.
• Is the numerical scheme stable for all Courant numbers?
The upwind scheme also suffers from numerical diffusion, and it is only first order
accurate in space. For some applications, particularly if there’s also diffusion, it might
just be good enough because the simple trick of doing FD forward or backward is closer
to the underlying physics of transport than, say, FTCS. There are some mantle convection
codes that use streamline upwind schemes.
So far, we employed explicit discretizations. You are probably wondering whether
implicit discretizations will save us again this time. Bad news: they are not well-suited
for this type of problem (try it and see). Implicit schemes behave like parabolic partial
differential equations (e.g. the diffusion equation) in that a perturbation at node ( j, n)
will affect the solution at all nodes at time level n + 1. With hyperbolic PDEs like the
advection equation or the wave equation, disturbances travel at a finite speed (the speed
of the material displacement) and will not affect all nodes at time level n + 1. So we have
to come up with something else.
Modified Crank-Nicolson
One approach to solving the advection equation is the previously introduced Crank-
Nicolson semi-implicit scheme. Here we modify it slightly by introducing a general mass
operator Mx = {δ, 1 − 2 δ, δ}.
(a) (b)
Figure 4.13: Illustration of the upwind (a) and leapfrog (b) schemes (from Press et al., 1993,
chap 19.1).
Setting the mass operator to δ = 0 gives us the previously seen Crank-Nicolson semi-
implicit finite difference discretization, while setting δ = 61 gives us the finite element
formulation. Below is eq. (4.122) written out with δ = 16 .
1 1 ∆t 1 1 1 ∆t
− v Tin−+11 − 1− Tin +1
+ + v Tin++11
6 4 ∆x 3 6 4 ∆x
1 1 ∆t n 1 n 1 1 ∆t
= + v Ti−1 + 1 − Ti + − v Tin+1
6 4 ∆x 3 6 4 ∆x
(4.123)
The finite element Crank-Nicolson advection scheme is stable for α ≤ 1 and provides
an improvement over previous schemes in that it is accurate to O(∆t, (∆x )3 ). This allows
us to reduce the number of grid points to reach the same accuracy as the other schemes
presented, as long as ∆t is kept small enough.
Staggered leapfrog
The explicit discretizations discussed so far were second order accurate in time, but only
first order in space. We can also come up with a scheme that is second order in time and
space
Tin+1 − Tin−1 T n − Tin−1
= −v x,i i+1 , (4.124)
2∆t 2∆x
called staggered leapfrog because of the way it is centered in shifted space-time (Figure 4.13b).
The computational inconvenience in this scheme is that two time steps have to be stored,
T n−1 and T n .
Exercise 4
• Program the staggered leapfrog method (assume that at the first time step T n−1 =
T n ).
• Try with different values of the Courant number α and compare the accuracy and
stability of the different methods.
• Bonus: Also program the two formulations of the Crank-Nicolson method with δ =
0 and δ = 61 .
The staggered leapfrog method works quite well regarding the amplitude and trans-
port phase as long as α is close to one. If, however, α 1 and the length scale of the
to-be-transported quantity is small compared to the number of grid points (e.g. we have a
thin plume), numerical oscillations again occur (those are due to the lack of communica-
tion between cells, which can be remedied by artificial diffusion). The conditions where
leapfrog does not work well are typically the case in mantle convection simulations (cf.
Figure 4.9). Onward ever, backward never.
Similarly, the Crank-Nicolson method works well for v ≤ 0.1 and α ≤ 0.1, and elim-
inates the staggered problem. But what happens for α ≥ 0.1? What about the finite
element formulation? What about computational time? Is Crank-Nicolson’s increased
accuracy worth the extra work? Is it well-suited for mantle convection problems?
MPDATA
This is a technique that is frequently applied in (older) mantle convection codes. The idea
is based on Smolarkiewicz (1983) and represents an attempt to improve on the upwind
scheme by adding some anti-diffusion, which requires iterative corrections. The results
are pretty good, but MPDATA is somewhat more complicated to implement. Moreover
we still have a restriction on the time step (given by the Courant criterion), for details see
Spiegelman (2004).
Basic idea
The basic idea of the semi-Lagrangian method is illustrated in Figure 4.14A for one di-
mension in space, x, and is given by the following, simplified scheme. Instead of allowing
the numerical scheme to transport noise in from unknown regions, the semi-Lagrangian
method uses transport by going back one (e.g. Euler) step.
For each point i at xi and time tn :
1. Assume that the future velocity v x (tn+1 , xi ) at xi is known. Under the assumption
that the velocity at the old time step is close to the future velocity
v x ( t n +1 , x i ) ≈ v x ( t n , x i ) (4.125)
and that velocity does not vary spatially
v x ( t n , x i −1 ) ≈ v x ( t n , x i ) ≈ v x ( t n , x i +1 ), (4.126)
we can compute the location X where the particle came from by
X = xi − ∆tv x (tn+1 , xi ). (4.127)
Note 2: Most of the MATLAB interpolation functions will by default not extrapolate outside
the
[min( xi ), max( xi )]
range and return NaN (not a number). If extrapolation is desired, ’extrap’ needs to be set
as an option when calling the interp1 function.
3. Assume that T (tn+1 , xi ) = T (tn , X ), i.e. temperature has been transported (along
“characteristics”) without any modification (e.g. due to diffusion).
This scheme assumes that no heat-sources were active during the advection of T from
T (tn , X ) to T (tn+1 , xi ). If heat sources are present and are spatially variable, some extra
care needs to be taken (Spiegelman, 2004, sec. 5.6.1).
Exercise 5
• Program the semi-Lagrangian advection scheme illustrated in Figure 4.14A. Is there
a Courant criterion for stability?
true path
A n+1 vx(n+1,i)
B n+1 vx(n+1,i)
time
time
n+1/2
Dt Dt x
n vx(n,i-1) vx(n,i) n vx(n,i-2) x vx(n,i-1) vx(n,i)
T(n,i-2)
x T(n,i-1) T(n,i) x
space space
Figure 4.14: Basics of the semi-Lagrangian method. See text for explanation.
Exercise 6
• Program the semi-Lagrangian advection scheme with the centered midpoint method
as illustrated in Figure 4.14B (cf. Spiegelman, 2004, p. 67).
Some care has to be taken if point X is outside of the computational domain, since
MATLAB will return NaN for the velocity (or temperature) of this point. If no extrap-
olation is desired, use the velocity v x (tn+1 , xi ) in this case. A pseudo-code is given
by
if isnan(Velocity)
Velocity = Vx(i)
end
2D advection example
The semi-Lagrangian method is likely a good, general advection algorithm (except in the
case of pseudo spectral methods), so this is the one we will implement in 2D.
Assume that velocity is given by (the rotational field)
v x ( x, z) = z (4.130)
vz ( x, z) = − x. (4.131)
Moreover, assume that the initial temperature distribution is Gaussian and given by
( x + 0.25)2 + z2
T ( x, z) = 2 exp (4.132)
0.12
Exercise 7
clear all
% Velocity
Vx = z;
Vz = -x;
for itime=1:nt
Vx_cen = Vx(iz,ix);
Vz_cen = Vz(iz,ix);
% for ??
% X =?
% Z = ?
if isnan(Vx_cen)
Vx_cen = Vx(iz,ix);
end
if isnan(Vz_cen)
Vz_cen = Vz(iz,ix);
end
% end
% X = ?;
% Z = ?;
% Interpolate temperature on X
% T_X = interp2(x,z,?,?,?, ’cubic’);
if isnan(T_X)
T_X = T(iz,ix);
end
Tnew(iz,ix) = T_X;
end
end
Tnew(1,:) = T(1,:);
Tnew(nx,:) = T(nx,:);
Tnew(:,1) = T(:,1);
Tnew(:,nx) = T(:,nx);
T = Tnew;
time = itime*dt;
figure(1),clf
pcolor(x,z,T), shading interp, hold on, colorbar
contour(x,z,T,[.1:.1:2],’k’),
hold on, quiver(x,z,Vx,Vz,’w’)
axis equal, axis tight
drawnow
pause
end
T̃ n+1 − T n ∂T
+ vx =0 (4.133)
∆t ∂x
for example by using a semi-Lagrangian advection scheme. Then solve the diffusion
equation
∂ T̃ ∂ ∂ T̃
ρc p = k + Q. (4.134)
∂t ∂x ∂x
For this, we assumed that Q is spatially constant; if not, one should consider to slightly
improve the advection scheme by introducing source terms. A good general method
would be to combine Crank-Nicolson for diffusion with a semi-Lagrangian solver for
advection (Spiegelman, 2004, sec. 7.2), but we will try something simpler first:
Exercise 8
Figure 4.16: Staggered grid definition. Properties such as viscosity and density inside a control
volume (gray) are assumed to be constant. Moreover, a constant grid spacing in x and z-direction
is assumed.
∂v x ∂vz
+ = 0 (4.135)
∂x ∂z
∂σxx ∂σxz
+ = 0 (4.136)
∂x ∂z
∂σxz ∂σzz
+ − ρg = 0 (4.137)
∂x ∂z
∂v x
σxx = − p + 2µ (4.138)
∂x
∂vz
σzz = − p + 2µ (4.139)
∂z
∂v x ∂vz
σxz = µ + , (4.140)
∂z ∂x
where ρ is density and g = {0, g} the gravitational acceleration. The density is where
these continuum and force balance equations (eqs. 4.135 to 4.137) couple to the the energy
equation, e.g. the diffusion and advection of temperature for mantle convection, discussed
in the previous sections.
It has been suggested that a particularly nice way to solve these equations is to use a
staggered grid (more about this later) and to keep as variables v x , vz and p (Gerya and Yuen,
2003; Gerya, 2009).1 Since there are three variables, we need three equations. Substituting
eqs. (4.138)-(4.140) into eq. (4.136) and eq. (4.137) leads to:
∂v x ∂vz p
+ = (4.141)
∂x ∂z γ
∂P ∂ ∂v x ∂ ∂v x ∂vz
− +2 µ + µ + = 0 (4.142)
∂x ∂x ∂x ∂z ∂z ∂x
∂P ∂ ∂vz ∂ ∂v x ∂vz
− +2 µ + µ + = ρg (4.143)
∂z ∂z ∂z ∂x ∂z ∂x
Note that we added the term γP to the incompressibility equations. This is a “trick” called
the penalty method, which ensures that the system of equations does not become ill-
posed. For this to work, γ should be sufficiently large (∼ 104 or so), so that the condition
of incompressibility (conservation of mass, eq. 4.135) is approximately satisfied.
4.9.3 Exercise
1. Discretize eqs. (4.141)-(4.143) on a staggered grid as shown on Figure 4.16.
2. A MATLAB subroutine is shown on Figure 4.17. The subroutine sets up the grid,
the node numbering and discretizes the incompressibility equations.
1 For a comparison of different finite difference approaches, see Deubelbeiss and Kaus (2008), for example.
% Solve the 2D Stokes equations on a staggered grid, using the Vx,Vz,P formulation.
clear
% Material properties
% phase #1 phase #2
mu_vec = [1 1 ];
rho_vec = [1 2 ];
% Input parameters
Nx = 20;
Nz = .9*Nx;
W = 1;
H = 1;
g = 1;
% Setup the interface
x_int = 0:.01:W;
z_int = cos(x_int*2*pi/W)*1e-2 - 0.5;
% Setup the grids----------------------------------------------------------
dz = H/(Nz-1);
dx = W/(Nx-1);
[X2d,Z2d] = meshgrid(0:dx:W,-H:dz:0);
XVx = [X2d(2:end,:) + X2d(1:end-1,:)]/2; % Vx
ZVx = [Z2d(2:end,:) + Z2d(1:end-1,:)]/2;
XVz = [X2d(:,2:end) + X2d(:,1:end-1)]/2; % Vz
ZVz = [Z2d(:,2:end) + Z2d(:,1:end-1)]/2;
XP = [X2d(2:end,2:end) + X2d(1:end-1,1:end-1)]/2; % p
ZP = [Z2d(2:end,2:end) + Z2d(1:end-1,1:end-1)]/2;
% Compute material properties from interface, properties are computed in the center of a control volume
Rho = ones(Nz-1,Nx-1)*rho_vec(2);
Mu = ones(Nz-1,Nx-1)*mu_vec(2);
z_int_intp = interp1(x_int,z_int,XP(1,:));
for ix = 1:length(z_int_intp)
ind = find(ZVz(:,1)<z_int_intp(ix));
Rho(ind(1:end-1),ix) = mu_vec(1);
Mu(ind(1:end-1),ix) = rho_vec(1);
fac = (z_int_intp(ix) - ZVz(ind(end),1))/dz;
Rho(ind(end),ix) = fac*rho_vec(1) + (1-fac)*rho_vec(2);
Mu(ind(end),ix) = fac*mu_vec( 2) + (1-fac)*mu_vec( 2);
end
% Setup numbering scheme----------------------------------------------------
Number_Phase = zeros(Nz + Nz-1, Nx + Nx-1); % Create the general numbering scheme
Number_ind = zeros(Nz + Nz-1, Nx + Nx-1); % Create the general numbering scheme
Number_Vx = zeros(Nz-1,Nx );Number_Vz = zeros(Nz ,Nx-1);
Number_P = zeros(Nz-1,Nx-1);
for ix=1:2:Nx+Nx-1, for iz=2:2:Nz+Nz-1, Number_Phase(iz,ix) = 1; end; end % Vx equations
for ix=2:2:Nx+Nx-1, for iz=1:2:Nz+Nz-1, Number_Phase(iz,ix) = 2; end; end % Vz equations
for ix=2:2:Nx+Nx-1, for iz=2:2:Nz+Nz-1, Number_Phase(iz,ix) = 3; end; end % P equations
num = 1;
for ix=1:size(Number_Phase,2)
for iz=1:size(Number_Phase,1)
if Number_Phase(iz,ix)~=0
Number_ind(iz,ix) = num;
num = num+1;
end
end
end
num_eqns = num-1;
ind_Vx = find(Number_Phase==1); Number_Vx(find(Number_Vx==0)) = Number_ind(ind_Vx);
ind_Vz = find(Number_Phase==2); Number_Vz(find(Number_Vz==0)) = Number_ind(ind_Vz);
ind_P = find(Number_Phase==3); Number_P (find(Number_P ==0)) = Number_ind(ind_P );
% Setup the stiffness matrix
A = sparse(num_eqns,num_eqns);
Rhs_vec = zeros(num_eqns,1);
% Setup the incompressibility equations------------------------------------
ind_list = [];ind_val = [];
[ind_list,ind_val] = Add_coeffs(ind_list,ind_val, Number_Vx(:,2:end ), ( 1/dx));%dVx/dx
[ind_list,ind_val] = Add_coeffs(ind_list,ind_val, Number_Vx(:,1:end-1), (-1/dx));
[ind_list,ind_val] = Add_coeffs(ind_list,ind_val, Number_Vz(2:end,: ), ( 1/dz));%dVz/dz
[ind_list,ind_val] = Add_coeffs(ind_list,ind_val, Number_Vz(1:end-1,:), (-1/dz));
% Add local equations to global matrix
for i=1:size(ind_list,2)
A = A + sparse([1:size(ind_list,1)].’,ind_list(:,i),ind_val(:,i),num_eqns,num_eqns);
end
num_incomp = length(ind_list);
% Perform testing of the system of equation, setup some given matrixes
mu = mu_vec(1);
Vx = cos(XVx).*sin(ZVx);
Vz = -sin(XVz).*cos(ZVz);
P = 2*mu*sin(XP ).*sin(ZP );
C = zeros(num_eqns,1);
C(Number_Vx(:)) = Vx(:);
C(Number_Vz(:)) = Vz(:);
C(Number_P(:)) = P(:);
Rhs = A*C;
% Check whether the compressibility equations are implemented correctly
max(abs(Rhs(1:num_incomp)))
Figure 4.17: MATLAB script Staggered_Stokes.m that sets up numbering, matrix A and that
solves the incompressibility equations.
USC GEOL557: Modeling Earth Systems 102
CHAPTER 4. FINITE DIFFERENCES
Figure 4.19: Staggered grid definition with the boundary points. Within the purple domain, the
finite difference scheme for center points can be applied. At the boundaries, we have to apply a
special finite difference scheme which employ fictious boundary nodes.
Add the discretization of the force balance equations (including the effects of grav-
ity) into the equation matrix A. Assume that the viscosity is constant and µ = 1 in a
first step, but density is variable.
An example is given in how to verify that the incompressibility equation is incorpo-
rated correctly. This is done by assuming a given (sinusoidal) function for, let’s say,
v x (e.g. v x = cos(ωx ) cos(ωz)). From the incompressibility equation (eq. 4.135) a so-
lution for vz than follows. By setting those solutions in the c vector, we can compute
Ac and verify that rhs for those equations is indeed zero.
3. Add free-slip boundary conditions on all sides (which means vz = 0, σxz = 0 on
the lower and upper boundaries and σxz = 0, v x = 0 on the side boundaries). Use
fictious boundary points to incorporate the σxz boundary conditions.
4. Assume a model domain x = [0; 1], z = [0; 1], and assume that the density below
z = 0.1 cos(2πx ) + 0.5 is 1, whereas the density above it is 2. Compute the velocity
and pressure, and plot the velocity vectors.
5. Write the code for the case of variable viscosity (which is relevant for the Earth since
rock properties are a strong function of temperature).
Figure 4.20: Discretization for the streamfunction approach. The boundary conditions are set
through fictious boundary points.
∂v x ∂vz
+ = 0 (4.144)
∂x ∂z
∂σxx ∂σxz
+ = 0 (4.145)
∂x ∂z
∂σxz ∂σzz
+ − ρg = 0 (4.146)
∂x ∂z
∂v x
σxx = − p + 2µ (4.147)
∂x
∂vz
σzz = − p + 2µ (4.148)
∂z
∂v x ∂vz
σxz = µ + (4.149)
∂z ∂x
By substituting eqs. (4.147)-(4.149) into eqs. (4.144)-(4.146), we obtain (compare sec. 4.9)
∂v x ∂vz
+ = 0 (4.150)
∂x ∂z
∂p ∂ ∂v x ∂ ∂v x ∂vz
− +2 µ + µ + = 0 (4.151)
∂x ∂x ∂x ∂z ∂z ∂x
∂p ∂ ∂vz ∂ ∂v x ∂vz
− +2 µ + µ + = ρg (4.152)
∂z ∂z ∂z ∂x ∂z ∂x
We can eliminate pressure from eqs. (4.151) and (4.152) by taking the derivative of eq. (4.151)
versus z and subtracting eq. (4.152) derived versus x. This results in:
∂2 ∂2
∂v x ∂vz
2 µ −2 µ +
∂x∂z ∂x ∂x∂z ∂z
∂2 ∂2
∂v x ∂vz ∂v x ∂vz ∂
2
µ + − 2 µ + = − ρg. (4.153)
∂z ∂z ∂x ∂x ∂z ∂x ∂x
We can also use the incompressibility constraint (4.150) to simplify things a little bit more:
∂2
∂vz
−4 µ +
∂x∂z ∂z
∂2 ∂2
∂v x ∂vz ∂v x ∂vz ∂
2
µ + − 2 µ + = − ρg (4.154)
∂z ∂z ∂x ∂x ∂z ∂x ∂x
Now we introduce a variable Ψ (the stream function) which is defined by its relation-
ship to the velocities as
∂Ψ
vx = (4.155)
∂z
∂Ψ
vz = − (4.156)
∂x
Note that Ψ satisfies incompressibility by plugging eqs. (4.155) and (4.156) into eq. (4.144).
By using Ψ, we can write eq. (4.154) as:
∂2
2
∂ Ψ
4 µ +
∂x∂z ∂x∂z
∂2
2
∂ Ψ ∂2 Ψ ∂2
2
∂ Ψ ∂2 Ψ
∂
2
µ 2
− 2 − 2 µ 2
− 2 = − ρg. (4.157)
∂z ∂z ∂x ∂x ∂z ∂x ∂x
Note that this equation now has 4th order derivatives for Ψ (easier to see for constant µ,
where we can pull the viscosity out of the derivatives.) The challenge is to solve eq. (4.157)
for Ψ given then density gradients.
4.10.3 Exercise
1. Discretize eq. (4.157) on a grid as shown on Figure 4.20.
2. A MATLAB subroutine is shown on Figure 4.21. The subroutine sets up the grid and
the node numbering. Finish the code by programming the discretized eq. (4.157).
To start simple, assume that viscosity is constant.
3. Add free-slip boundary conditions on all sides (which means vz = 0, σxz = 0 on the
lower and upper boundaries and σxz = 0, v x = 0 on the side boundaries; you’ll have
to write these equations in terms of Ψ and employ fictious boundary points).
4. Assume a model domain x = [0; 1], z = [0; 1], and assume that the density below
z = 0.1 cos(2πx ) + 0.5 is 1, whereas the density above it is 2. Compute the velocity,
and plot the velocity vectors.
5. Write the code for the case of variable viscosity (which is relevant for the Earth since
rock properties are a strong function of temperature).
Figure 4.21: Code Streamfunction_Stokes.m that initializes the grid and node numbering for the
2D streamfunction approach.
We briefly discuss two examples for solving wave propagation type problems with
finite differences, the acoustic and the seismic problem.
∂2 p
1
= K∇ · ∇p . (4.162)
∂t2 ρ
If we assume that density is constant, we can introduce a parameter, v B , which turns out
to be a velocity of propagation
∂2 p
= v2B ∇2 p, (4.163)
∂t2
where the bulk sound velocity is s
K
vB = . (4.164)
ρ
Simplifying to a 2D case, we have
∂2 p ∂2 p ∂2 p
= v2B + 2 . (4.165)
∂t2 ∂x2 ∂z
The equation for propogation of SH waves, the transverse components of S waves, in
seismology has a similar form as eq. (4.165):
∂2 u
2
∂ u ∂2 u
2
= vSH + 2 , (4.166)
∂t2 ∂x2 ∂z
where u is the displacement and VSH is the velocity of the SH component.
Likewise, a similar equation also applies for tsunami waves at long wavelengths, in
the “shallow water approximation”,
∂2 ξ
2
∂2 ξ
2 ∂ ξ
= vSW + 2 . (4.167)
∂t2 ∂x2 ∂z
Here, ξ is the height of the tsunami wave, vSW is the velocity controlled by the water
depth H as p
vSW = gH. (4.168)
To solve eqs. (4.165)-(4.167), with finite differences, we use the mesh shown in Fig. 4.23.
n = P (i∆h, j∆h, n∆t ) and v
Here, we have pi,j i,j = v (i∆h, j∆h ) (meaning bulk velocity,
v B ). Applying the 2nd -order, second derivative formula to the acoustic wave equation
eq. (4.165),
n −1 n + p n +1
" n n + pn n n + pn
#
pi,j − 2pi,j i,j pi−1,j − 2pi,j i + 1,j pi,j − 1 − 2pi,j i,j + 1
= v2i,j + . (4.169)
∆t2 ∆h2 ∆h2
where
∆h2
ai,j = v2i,j . (4.171)
∆t2
Then, the pressure or displacement at time step n + 1 can be derived explicitly from time
step n and n − 1 as in eq. (4.170), though two solutions have to be stored. Note that we
use 2nd -order second derivatives in eq. (4.170).
Two considerations are required for choosing suitable time step ∆t and spatial step ∆h:
grid dispersion and stability:
• To achieve an accurate solution, we need at least 12 points per wavelength for space
for a scheme with 2nd order accuracy. For a 4th order scheme, a minimum of 6.5
points per wavelength are required. For a fixed frequency, this minimum wave-
length is determined by the minimum velocity (vmin ), so the accuracy of the system
is governed by (vmin ). Following a stability analysis, we can derive the stability
requirement here as:
1 ∆h
∆t ≤ √ (4.172)
2 vmax
where vmax is the maximum velocity on the grid.
Exercise 1
• Program the 2D acoustic wave propagation in standard grid scheme as in Fig. 4.23
(wave_acoustic_2D.m). Study the wavefield and seismograms with different choices
of ∆t and ∆h and demonstrate how ∆t and ∆h affect the stability and grid dispersion
in the program.
• Introduce heterogeneities in the velocities, such as a thin layer with half velocity, and
describe the difference from the isotropic model, especially how this layer affects the
observed seismograms. Run the code with the velocity inside the thin layer being
zero and explain the result.
∂2 u x ∂τxx ∂τxz
ρ 2
= + , (4.173)
∂t ∂x ∂z
∂2 u z ∂τxz ∂τzz
ρ 2
= + , (4.174)
∂t ∂x ∂z
∂u x ∂uz
τxx = (λ + 2µ) +λ , (4.175)
∂x ∂z
∂uz ∂u x
τzz = (λ + 2µ) +λ , (4.176)
∂z ∂x
and
∂u x ∂uz
τxz = µ + . (4.177)
∂z ∂x
Here, (u x , uz ) are the particle displacements. For seismic waves, they are typically called
radial and vertical components, respectively, if they are recorded at surface. Further, τ is
the stress tensor, and λ and µ the elastic, Lamé coefficients (µ is shear modulus).
Typically, those equations are solved for particle velocities as U = ∂u x
∂t and V = ∂t .
∂uz
Then, the system is transformed into the frst-order hyperbolic system, introducing the
abbreviations Σ = τxx , T = τzz , Λ = τxz ,
∂U ∂Σ ∂Λ
=b + , (4.178)
∂t ∂x ∂z
∂V ∂Λ ∂T
=b + , (4.179)
∂t ∂x ∂z
∂Σ ∂U ∂V
= (λ + 2µ) +λ , (4.180)
∂t ∂x ∂z
∂T ∂V ∂U
= (λ + 2µ) +λ , (4.181)
∂t ∂z ∂x
∂Λ ∂U ∂V
=µ + , (4.182)
∂t ∂z ∂x
with the buoyancy, b = 1/ρ.
A typical seismic wave propagation problem needs to deal with medium with variable
Poisson’s ratio, ν, which can be defined as
λ
ν= . (4.183)
2( λ + µ )
For the special case of λ = µ, ν = 1/4, and many rocks have Poisson’s ratios not far
from 1/4. For liquids, ν → 0.5. For seismic wave propagation, this is particularly impor-
tant when ocean water or the outer core of the Earth are needed to be considered in the
problem, which is hard to be resolved with the traditional set up of grid as in Fig. 4.23.
To satisfy both requirements for stability and grid dispersion at those problems, a P-
SV staggered-grid scheme is applied. Note the structure of the elastic wave problem,
eq. (4.178)-eq. (4.182), they allow the stress and particle velocity to be spatially interlaced
on the grids as in Fig. 4.24. The staggered-grid scheme allows the spatial derivative to be
computed to a much higher accuracy (e.g. Levander, 1988). This computational aspect is
similar to the staggered grid finite difference approach to the Stokes problem, discussed
in sec. 4.9.
To add the complexity, the stress and velocity field can also staggered in time. We
follow the explicit scheme and first update the velocities from time half-steps k − 1/2 to
k + 1/2, i.e. centered on time, k∆t, using second order, finite difference equations for the
first derivatives in eqs. (4.178)-(4.182) (Figure 4.24). Introducing
∆t
S= , (4.184)
∆h
we find
Uik++1/2,j
1/2
= Uik+−1/2,j
1/2
+ bi+1/2,j S Σik+1,j k
− Σi,j
+ bi+1/2,j S Λki+1/2,j+1/2 − Λki+1/2,j−1/2 , (4.185)
Vi,jk++1/2
1/2 = Vi,jk−+1/2
1/2
k k
+ bi,j+1/2 S Λi+1/2,j+1/2 − Λi−1/2,j+1/2
k k
+ bi,j+1/2 S Ti,j +1 − Ti,j . (4.186)
k +1
Ti,j = k
Ti,j + (λ + 2µ)i,j S Vi,jk++1/2
1/2 − Vi,jk+−1/2
1/2
+ λi,j S Uik++1/2,j
1/2
− Vik−+1/2,j
1/2
, (4.188)
+1 k+1/2 k+1/2
Λki+ 1/2,j+1/2 = Λk
i +1/2,j+1/2 + µ i +1/2,j+1/2 S Vi +1,j+1/2 − Vi,j+1/2
k+1/2 k +1/2
+ µi+1/2,j+1/2 S Ui+1/2,j+1 − Ui+1/2,j . (4.189)
Exercise 2
• Program the 2D elastic wave propagation in staggered grid scheme as in Fig. 4.24
(wave_elastic_staggered_2D.m). Choose ∆t and ∆h and describe the wavefield
(both vertical and horizontal components) for the model with uniform velocities.
Identify the first P and SV arrivals on the recorded seismograms.
• Include a thin liquid layer (vS = 0) in the model and explain the result. Note for a
typical wave propagation problem, the input models are v P , vS , and density ρ, so
the conversions to λ and µ are required in the program.
(see sec. 4.8.2 and, e.g., Mozco et al., 2004, p. 33ff for a discussion of averaging schemes).
Finite elements
• Textbooks
– The recommended background reading for this, finite element, part of the class
is Hughes (2000). The Hughes (2000) approach is widely followed in mantle con-
vection modeling, and aspects of the classic codes ConMan (King et al., 1990) and
CitcomS (Moresi and Solomatov, 1995; Zhong et al., 2000) (both now maintained
and developed by CIG, geodynamics.org) follow the approach and notation of
Hughes (2000) (see also Zhong et al., 2007).
– For additional reference, you might want to consider the MATLAB -based finite
element course by Kwon and Bang (1996) and the comprehensive and high-level
treatment by Bathe (2007).
115
CHAPTER 5. FINITE ELEMENTS
u|∂Ω = g (5.2)
∂u
ni = n · ∇u = h (5.3)
∂xi
If the PDE is, for example, an elastic deformation problem, then u would be displace-
ments, and Dirichlet conditions of g = 0 correspond to “no-slip”, i.e. no deformation
at the boundaries. Likewise, Neumann conditions for h = 0 would correspond to zero
tractions (the derivative of displacement times modulus are stresses).
The FE analysis then proceeds by two steps:
1. Converting the governing PDE from the regular, “strong” form (which we used for
FD) to the weak integral form (see below).
We will provide a highly abbreviated treatment, lacking any mathematical rigor. More-
over, we will omit the detailed discussion of different element types, or shape functions,
as well as implementation issues such as order of integration. Those issues are very im-
portant in practice, as choices in shape functions and element type may strongly affect
solution robustness and accuracy. However, we regrettably do not have the time to delve
into this (see, e.g. Hughes, 2000, for a more comprehensive treatment). For a specific man-
tle convection code pertaining to the commonly used ConMan and Citcom, see also Zhong
et al. (2007).
To contrast finite elements with finite difference methods, Table 5.1 summarizes the
main differences and how they affect usage.
Table 5.1: Comparison of the characteristics of finite difference and finite element approaches
∂2 T H (x)
2
− = 0. (5.5)
∂x κ
Mathematically, we require s to be smooth for a solution for u to exist. Additionally,
we will require
and
∂u(1)
=h (Neumann) (5.7)
∂x
boundary conditions (BCs), which closes the system for a two point boundary value prob-
lem. This formulation of the PDE with all original derivatives in place is called the strong
form. Eqs. (5.4) and (5.6) can, of course, be solved analytically by integration.
∂2 u
+x =0 (5.8)
∂x2
then by integration
∂u 1 2
+ x + c1 = 0 (5.9)
∂x 2
and integrating again
1
u + x 3 + c1 x + c2 = 0 (5.10)
6
1
u ( x ) = − x 3 − c1 x − c2 (5.11)
6
and use BCs eq. (5.6) and eq. (5.7)
∂u(1) 1
⇒ u(0) = −c2 = g; = − − c1 = h (5.12)
∂x 2
1 3 1
⇒ u( x ) = − x + h + x+g (5.13)
6 2
The analytical solution eq. (5.13) will be used in the numerical problem set to test the
approximate, numerical solutions. For more complicated, realistic problems, typically
no analytical solutions can be found (which is why we do numerical analysis in the first
place).
From the strong form, we will now move to a variation, or “weak form”. We require
that
(a) the trial solutions of u, among all possible solutions satisfy the essential BC
u (0) = g (5.14)
∂2 u
Z Z
− w − ws = 0 (5.18)
∂x2
from the integration by parts rule
Z Z
0
a b = ab | − b0 a (5.19)
Z 1 1 Z
∂w ∂u ∂u
⇒ − w − w s = 0. (5.20)
0 ∂x ∂x ∂x 0
∂u(1)
With ∂x = h and w(0) = 0, this can be written as
Z 1 Z 1
∂w ∂u
= h w (1) + w s dx (5.21)
0 ∂x ∂x 0
This is the weak form of the PDE. Equations of this type in mechanics are called “vir-
tual work”, or virtual displacement formulations (the w are the virtual displacements).
It can be shown that the weak and the strong form are identical (Hughes (2000), sec. 1.4)
and the FE method proceeds from eq. (5.21) by assuming u and w can be taken from a
simplified functional space, typically based on low order polynomials.
It is useful to define a shorthand notation
Z 1 Z 1
∂w ∂u
a(w, u) = dx and (w, s) = w s dx (5.22)
0 ∂x ∂x 0
Then, we can write eq. (5.21) as
ũ = ṽ + g̃ (5.27)
where g̃ is a given function such that g̃(0) = g, which satisfies the BCs because
where we solve for the LHS and the RHS is determined by BCs.
This is an example of a weighted residual method, there are other approaches such as
the Petrov-Galerkin method. The Galerkin method used here is the simplest because it
assumes that ṽ and w̃ are from the same function space, i.e. the same shape functions (see
below) are used for the solution ũ and the weights w̃.
NA (0 ) = 0 ∀A (5.31)
such that w̃(0) = 0 can be fulfilled, and NA ∈ W . If we also introduce a special shape
function for the boundary
then
such that ũ(0) = g. If we substitute eqs. (5.30) and (5.34) into eq. (5.29), then
! ! " # !
a ∑ c A NA , ∑ db NB = ∑ c A NA , s + ∑ c A NA (1 ) h−a ∑ c A NA , g N̂1 .
A B A A A
(5.35)
Because of bilinearity, we can write ∑ A c A G A = 0 with
GA = ∑ a( NA , NB ) dB − ( NA , s) − NA (1) h + a( NA , N̂1 ) g. (5.36)
B
The Galerkin equation (5.29) is supposed to hold for all w, therefore all c A , which means
that G A = 0. So, for all A
where the element size h A may vary, and a general grid spacing may be defined as h =
max (h A ).
We can then choose interior shape functions for 2 ≤ A ≤ n as
1
h A−1 ( x − x A−1 ) for x A −1 ≤ x ≤ x A ,
NA ( x ) = 1 (5.45)
h A ( x A +1 − x ) for x A ≤ x ≤ x A+1 ,
0 otherwise.
1
N̂1 ( x ) = ( x2 − x ) for x1 ≤ x ≤ xn+1 and (5.46)
h1
1
N̂n+1 ( x ) = ( x − xn ) for xn ≤ x ≤ xn+1 . (5.47)
hn
An illustration of interior and boundary shape functions is shown in Figure 5.2; note that
NA = 1 at x = x A and zero for other nodes.
With this choice, the shape functions are zero outside the vicinity of A. They have local
support, which means that K is a sparse matrix because the
Z 1
∂NA ∂NB
a( NB , NA ) = dx (5.48)
0 ∂x ∂x
integral is zero for B > A + 1. The K matrix is a banded matrix and the bandwidth
depends on how the nodes are numbered (leading to an optimization problem during
mesh design) and what basis functions are used. Besides symmetry and bandedness, K
is also positive definite, which means that
cT K c ≥ 0 (5.49)
Figure 5.3: Example of 1-D, linear shape functions in global (left, N ( X )) and element-local (right,
N (ξ ), with ξ ∈ [−1; 1]) coordinate systems.
1 1
x (ξ ) = ( h A ξ + x A + x A +1 ) ⇔ ξ (x) = (2x − x A − x A+1 ) (5.50)
2 hA
u ( x ) = ∑ d A NA ( x ) ⇔ u(ξ ) = N1 (ξ ) d1 + N2 (ξ ) d2 . (5.51)
We can express the global shape function NA of eq. (5.45) in a local coordinate system as
1
Na (ξ ) = (1 + ξ a ξ ) for a = 1, 2 (5.52)
2
1 1
i.e. N1 (ξ ) = (1 − ξ ); N2 (ξ ) = (1 + ξ ) with ξ ∈ [−1; 1] . (5.53)
2 2
where x ea are the global nodes that belong to the element e. For the assembly of the stiff-
ness matrix, derivatives of Na and x e with respect to ξ are required. We note that
∂Na ξa (−1) a
= = , (5.55)
∂ξ 2 2
∂x he
= , (5.56)
∂ξ 2
and −1
∂x e
∂ξ 2
= = . (5.57)
∂x ∂ξ he
−1
∂x e
Note 1: For higher dimensional problems, the term ∂ξ will be a matrix inverse.
Note 2: The choice of shape function is determined by the type of element. E.g. with a
two node element, shape functions can only be linear. If the solution to be approxi-
mated is u, then u will vary linearly over the element, and derivatives, ∂ x u, will be
constant.
The integrals over the problem domain [0; 1] can be written as summations over ele-
ments, therefore
n
K= ∑ Ke with K e = [K eAB ] (5.62)
e =1
(5.63)
F= ∑F e
with e
F = { FAe } (5.64)
e
Z
∂NA ∂NB
K eAB = dΩe dx (5.65)
∂x ∂x
∂NA ∂ N̂1
Z Z
FAe = NA s dx + h δe,n δA,n+1 − g dx (5.66)
Ωe Ωe ∂x ∂x
where the element domain Ωe = [ x1e ; x2e ].
Since the NA only have local support K eAB = 0 if A 6= {e or e + 1} or B 6= {e or e + 1},
and FAe = 0 if A 6= {e or e + 1}, and we can obtain the global stiffness matrix and force
vector by summing up elemental contributions K e and f e
Z
e ∂Na ∂Nb
Kab = a( Na , Nb ) = dx (5.68)
Ωe ∂x ∂x
(5.69)
e g
Z Ka,1 for e = 1,
fa = Na s dx + 0 for e = 2, . . . , n − 1, (5.70)
Ωe
δa,n+1 h for e = h.
The assembly proceeds as symbolized in Figure 5.4, and placing the K e element-local
matrix into the global stiffness matrix requires use of an assignment operator or array.
This is discussed in the worked example of the problem set.
∂Na
Since ∂ξ is independent of element data, the computation only has to be performed once,
∂ξ
independent of the actual nodal values for the solution. The derivatives ∂x
∂ξ and ∂x do,
however, depend on the element shape and need to be computed for each geometric
(mesh) configuration.
For the source term, we use the approximation
2
s̃ = ∑ sa Na (5.80)
a =1
i.e. the source term is assumed to vary linearly across the element with Na . Then, we can
write
Z 1
∂x (ξ )
Z
Na ( x ) s̃( x ) dx = Na ( x (ξ )) s̃( x (ξ )) dξ (5.81)
Ωe −1 ∂ξ
he 2
Z 1
=
2 ∑ Na (ξ ) Nb (ξ ) dξ sb . (5.82)
b =1 −1
R1
Since −1 Na Nb = 13 (1 + δab ),
he
e 2 1 s1
s = + boundary terms (5.83)
6 1 2 s2
(5.84)
he
2 s1 + s2
= + boundary terms (5.85)
6 s1 + 2 s2
Here, NA are the shape functions in the interior, B is another global node, and N̂1 the
boundary shape function for the essential boundary condition g. This can be further
abbreviated by
Z 1
∂NA ∂NB
K AB = a( NA , NB ) = dx (5.89)
∂x ∂x0
FA = ( NA , s) + NA (1)h − a( NA , N̂1 ) g (5.90)
Z 1 Z 1
∂NA ∂ N̂1
= NA sdx + δA,n+1 h − dx g, (5.91)
0 0 ∂x ∂x
where we have used the definitions of the bi-linear forms a(·, ·) and (·, ·) from before, and
the Kronecker delta (
1 if i = j
δi,j = (5.92)
0 else
is used for the flux boundary condition (also see Hughes, 2000, chap 1). Note that it is
sometimes helpful to not think about nodes but about elements. The weak form of the
equations are satisfied on average per element, and by constructing an appropriate map-
ping/numbering we can easily go from a single element to a 1D or 2D domain.
The approximate solution of u after discretization of the weak form is given by
n +1
ũ( x ) = ∑ d A NA ( x ) + N̂1 g = ∑ d A NA ( x ), (5.93)
A =2
where the latter summation implies choosing the boundary shape function and BC if
needed. The vector d = {d A } values have to be obtained by solution of the matrix equa-
tion
Kd = F, (5.94)
with K = {K AB } and F = { FA }.
We discussed previously how the integration over the domain can be broken down
into summation over integrals over each element (see sec. 3.2). This integration is most
easily performed in a local coordinate system −1 ≤ ξ ≤ 1 between the two nodes of each
element, which has a mapping to the corresponding, global coordinate interval [ x A ; x A+1 ].
We can also express the shape functions as x (ξ ),
The global K matrix and the F vector are then assembled by looping over all elements
1 ≤ e ≤ n and adding each element’s contribution for shared nodes. By change of inte-
gration variables and the chain rule, those elemental contributions follow as
e 1 1 −1
k = , (5.96)
∆x −1 1
where ∆x is the element size, x A+1 − x A , and for the force term
( k g
a1 for e = 1
∆x 2s1 + s2
fe = + δa,n+1 for e = n (5.97)
6 s1 + 2s2
0 else
where we have assumed that the source function s varies linearly over each element, and
s1 and s2 are the contributions from each local node a within the element from s( x ). After
assembly, one needs to ensure that each row of the global K matrix that corresponds to a
fixed value (Dirichlet) boundary condition, will only have a diagonal entry and the other
columns for this row are zero.
5.3.2 Exercises
(a) Download heat1dfe.m, and all helper routines for this section. Read through the
implementation of what is summarized above, heat1dfe.m, and understand this
code.
(b) Fill in the blanks in heat1dfe.m and experiment with a solution of eq. (5.86) for
n = 3 elements.
(a) Print out the stiffness matrix (full(stiffness)) to appreciate its banded struc-
ture. Does this look familiar to you?
(b) Choose the MATLAB solver (solver=0) and plot the finite element solution at
the nodes, as interpolated within the elements, and compare with the analytical
solution.
(c) Reads the following section on solving linear systems of equations, if you have
time, and work on the examples (optional).
Notes:
• Symmetry means K = KT , where KT is the transpose, KijT = K ji .
• Positive definite means that c T Kc > 0 for any non-zero c. Graphically, this corre-
sponds for a 2 × 2 matrix to a well defined minimum (lowest) point in a curved
landscape, which is important for iterative methods (e.g. Shewchuk, 1994).
Ax = b (5.101)
in the sense that |Ax − b| = min, i.e. deviations from the true solution are minimized,
for a matrix A that may be under-determined, i.e. not simply invertible. It can be
shown that the general least squares solution x LS is given by
−1
T
x LS = A · A · AT · d, (5.102)
−1
where AT · A is the generalized inverse (which exists even if the inverse of A,
−
A , does not exist because A is singular). AT · A is symmetric and positive definite,
1
meaning that Cholesky is also the method of choice for direct approaches to find
x LS .
where the summation is over all l but for l = j. The Jacobi method following eq. (5.104)
is implemented in jacobi.m. It serves mainly illustrative purposes but is guaranteed
Gauss-Seidel method
An improvement over the Jacobi method is the Gauss-Seidel (GS) approach, where the
iterative rule is
(D + K)di+1 = F − Udi . (5.105)
The main benefit is that di+1 can be computed from di directly, without having to store a
full previous solution, following
!
1
dij+1 = Fj − ∑ K jl dil+1 − ∑ K jl dil . (5.106)
K jj l<j l>j
Note that this operation can be done “in place”, and is implemented in gauss seidel.m.
The GS method will converge (somewhat faster than the Jacobi method) for diagonally
dominant, or positive definite and symmetric K.
Setting w = 1 will reduce SOR to the GS method. The optimal value of w is dependent
upon the matrix K, but setting w = 0.5 is a good starting point. The method has been
rigorously shown to converge for symmetric, positive definite matrices K for 0 < w < 2.
Exercise 1
• Plot the Jacobi, GS, and SOR solutions for 32 elements and a tolerance of 10−4 , 10−5 ,
and 10−6 on one plot each; comment on the accuracy and number of iterations re-
quired. Can you improve the definition of tolerance for the Jacobi method?
• Choose a tolerance of 10−6 , and record the number of iterations required to solve the
1-D FE example problem using the Jacobi and GS methods for increasing number of
elements, e.g. for 8, 16, 32, 64, and 128 elements. (You might want to automate these
computations and not wait until convergence and record the results by hand.)
• Plot the number of iterations against number of elements for both methods.
Conjugate gradient
You have now seen that while the Gauss-Seidel (GS) method is an improvement on the
Jacobi approach, it still requires a large number of iterations to converge. This makes
both methods impractical in real applications and other approaches are commonly used.
One of those is the conjugate gradient (CG) method which works for positive-definite,
symmetric, square (n × n) matrices. The CG method is explained in a nice, geometric
fashion by Shewchuk (1994). We cannot explore the motivation behind CGs in detail,
but conjgrad.m provides a pretty straightforward MATLAB implementation which you
should check out.
The CG method provides an exact solution after n iterations, which is often a pro-
hibitively large number for real systems, and approximate solutions may sometimes be
reached for a significantly smaller number of iterations. There are numerous tweaks in-
volving modifications to the conjugate gradient method that pertain to “pre-conditioners”
where we solve
M−1 Kd = M−1 F, (5.108)
for some M which approximates K but is simpler to handle than K. The best choice of
these is, for some applications, an active area of research (e.g. May and Moresi, 2008).
For sparse least-squares problems, such as for the typically-mixed determined seismic
tomography problem, the LSQR approach of Paige and Saunders (1982) is a popular choice
that is used by many researchers for linear inversions. The robustness of the iterative
solution compared to direct solvers was explored by Boschi and Dziewoński (1999), for
example (it works!).
Exercise 2
• Switch the solver from the GS method to conjugate gradient and increase the max-
imum iteration number stepwise from a fraction of n to the full n (as determined
by the number of elements which you should choose large, e.g. 200, for this exer-
cise). Test different initial guesses for di (e.g. all zero, random numbers), record the
convergence and comment on the solution.
Multigrid method
An interesting philosophy to solving PDEs of the type we are considering for the 1-D fi-
nite element example is by using several layers of variable resolution grids (e.g. Press et al.,
1993, sec. 19.6). The insight is based on the observation that the Gauss-Seidel method is
very good at reducing short-wavelength residuals in the iterative solution for d (“smooth-
ing”), but it takes a long time to reduce the largest wavelength components of the residual.
(You should try to plot successive solutions of the GS method compared to the analytical
solution for different starting d0 to visualize this behavior.)
For the multigrid (MG) method, the idea is to solve the equations to within the de-
sired tolerance only at a very coarse spatial discretization, where only a few iterations
are required. Then, the solution is interpolated up to finer and finer levels where only
a few GS iterations are performed at each level to smooth the solution. One then cycles
back and forth until convergence is achieved at the finest, true solution level. There are
several different approaches that are all called “multigrid” methods and basically only
share the same philosophy. Differences are, for example found in terms of the way the
cycling between fine and coarse resolutions are conducted (e.g. Briggs et al., 2000), and
we will only discuss the “V cycle” method. Multigrid methods are now implemented in
most 3-D finite element methods (Zhong et al., 2007) because MG has, ideally, the perfect
scaling of O( N ) where N is the size of the problem. MG methods areas of active research
(e.g. algebraic multigrid, which is related to adaptive mesh refinement).
The multigrid method is based on expressing the PDE on L MG levels of resolution
where the number of nodes in each level, ni , is given by
n i = b × 2i − 1 + 1 for i = 1, 2, . . . , L, (5.109)
where b is the base, typically a small number such as 2 or 4. At each ith level, we need
to construct separate stiffness matrices, Ki , and the corresponding force vector where the
resolution for the i = L solution is the best approximation to Kd = F, and the forcing is
only needed to be specified at FL (see below).
An example implementation may proceed like so (see, e.g. Press et al., 1993, sec. 19.6 for
some alternatives): We start at the highest level, L, and perform only a few, fixed number
of GS iterations for an rough approximate d L from
K L d L = FL (5.110)
to remove the short wavelength misfit starting from an initial guess d L = 0. The residual
is then given by
R L = FL − Kd L . (5.111)
We then project, or restrict, the residual to a coarser grid at L − 1 by a projection operator
P
R L−1 = PL→ L−1 R L . (5.112)
P will be some stencil giving more weight to the fine resolution nodes that are closer to
the coarse resolution node to which we project. We next GS iterate
Ki δdi = Ri (5.113)
for i = L − 1 for another small number of iterations (initializing di again with 0), per-
forming another “smoothing” step, reducing short wavelength fluctuations. Note that
eq. (5.113) now operates on the residual and not the load vector F such that we are com-
puting corrections of d, δd. We then repeat the smoothing and projection steps down to
i = 1 where eq. (5.113) can be solved quickly and exactly. This completes the downward
leg of the V cycle where the longest wavelength residual has been addressed.
Next, we have to propagate the correction δd1 from i = 1 to i = 2 and higher res-
olutions by means of a “prolongation”, i.e. an interpolation to higher resolution by an
interpolation operator I
δdi+1 = Ii→i+1 δdi . (5.114)
I may be a linear interpolation, for example, which is easy to compute for the mesh struc-
ture eq. (5.109). This upward interpolated δdi+1 can then be smoothed by using it as a
starting guess for a fixed number of GS iterations for
with δRi+1 = −Ki+1 δdi+1 and weighting α = (δRi+1 · Ri+1 )/|δRi+1 |2 . We continue by
projecting δdi in this fashion up to i = L, where we update d L = d L + δd L , which com-
pletes the upward leg of the V. The whole V cycle is then repeated until the desired toler-
ance for d L is reached at which point d L = d. Details of the implementations of the MG
method, such as the smoothing, restriction, and prolongation operations, depend on the
problem and the boundary conditions (e.g. Press et al., 1993; Briggs et al., 2000).
Exercise 3
• Compare the number of iterations needed for the MG solver with that of the GS
method for 32, 64, 128, 256 numbers of elements.
• Plot the scaling of the number of iterations, or time spent in the multigrid subrou-
tine, with the number of elements.
Temperature
1 1
0.9 0.9
0.8 0.8
0.7 0.7
0.6 0.6
0.5 0.5
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0 0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Figure 5.5: Coarse finite element mesh with solution for temperature, allowing for an elliptical
inclusion and a boundary at mid-model height.
We will now consider the solution of 2D boundary value problems using finiteele-
x1 x
ments, which can be easily expanded to three dimensions. We write x = =
x2 z
for a location vector xi with i = 1, 2 and a normal vector n on the boundary Γ of the do-
main Ω. As an example problem we will now revisit the linear heat conduction problem.
where repeated indices imply summation, q is heat flux, κij = κ ji is the conductivity
matrix (we use k normally for diffusivity but we wish to distinguish it from the stiffness
matrix), and T is temperature. In vector notation,
q = −κ · ∇ T κ = [κij ] (5.120)
∇·q = H (5.122)
∂qi ∂κij ∂T
or = =H for Ω (Poisson Eq.)
∂xi ∂xi ∂x j
with BCs:
T=g on Γg and
− qi ni = h on Γh (5.123)
where Γg and Γh are the parts of the boundary where fixed temperature or fixed flux
conditions apply, respectively. For isotropy, we recover
∂ ∂T ∂ ∂T
δij κ =κ =H (5.124)
∂xi ∂x j ∂xi ∂xi
2
∂ T ∂2 T
conductivity → κ + 2 =H. (5.125)
∂x2 ∂z
The weak form representation (first substitute the test function and integrate over the
domain, second integrate by parts) of eq. (5.123) is given by
Z Z Z
∂w
− dΩ q = dΩ wH + dΓ wh (5.126)
Ω ∂xi i Ω Γh
where the LHS is the diffusion term, the first term on the RHS corresponds to volumetric
heating, and the second term on the RHS is the flux through the boundary. See Hughes
(2000), sec. 2.3 for the derivation. Eq. (5.126) can then be expressed as
with
Z
∂w ∂T
a(w, T ) = dΩ κij
Ω ∂xi ∂x j
Z
(w, H ) = dΩ wH (area integral over Ω)
ZΩ
(w, h) Γ = dΓ wh. (line integral over Γ)
Γ
such that
Z
a(w, T ) = dΩ (∇w) T κ ∇ T (5.130)
Ω
1 0
with κ = κ for isotropy in 2-D.
0 1
Using the Galerkin approach of choosing the trial and weighting functions from the same
function space, we again posit for the solution
T̃ = ṽ + g̃ (5.131)
where we have again distinguished between interior nodes and shape functionm A ∈
I , and those on the Dirichlet boundary, with A ∈ B . Arguing as for the 1-D case, the
following assembly rules result
Kd=F (5.134)
K = K PQ K PQ = a( NA , NB ), (5.135)
where 1 ≤ P, Q ≤ neq and the number of free equations is given by the total number of
nodes minus the number of nodes on the Dirichlet boundary B .
P and Q can be computed from a 1D array that maps a global node A into a global equa-
tion number
P for A ∈ I
ID ( A) = (5.136)
0 for A ∈ B
ṽ( x) = ∑ NA d A (5.137)
and
F = { Fp } (5.138)
where
FP = ( NA , H ) + ( NA , h) Γh − ∑ a( NA , NB ) q B (5.139)
B∈B
Z
e e
K PQ = a( NA , NB ) = (∇ NA )T κ (∇ NB ) dΩ . (5.141)
Ωe
The RHS in eq. (5.141) corresponds to integrating over each element’s area.
nel
F= ∑ Fe F e = { FPe } (5.142)
e =1
Z Z
FPe =
Ωe
dΩ NA H +
Γhe
dΓ NA h − ∑ a( NA , NB )e q B (5.143)
B∈B
where Γhe is the part of the Neumann (flux) boundary within element e, and P = ID ( A), Q =
ID ( B). Within each element we compute for new nodes per element with 1 ≤ a, b ≤ nen
Z
e e e e
K = {Kab } Kab = a( Na , Nb ) = dΩ (∇ Na ) T κ (∇ Nb ) (5.144)
Ωe
f e = [ f ae ] (5.145)
Z Z nen
fa =
Ωe
Na f a dΩ +
Γhe
Na h dΓ − ∑ Kabe gbe (5.146)
b =1
where gbe = g( xeb ) for prescribed g and zero otherwise. It is convenient to write
Z
e
K = dΩ B T D B (5.147)
Ωe
where D is nsd × nsd (rows × columns); nsd , number of spatial dimensions. In our case
D is 2 × 2 and D = κ. B is nsd × nen such that B = { B1 , B2 , . . . , Bnen } and
Ba = ∇ Na (5.148)
is nsd × 1.
The B and D matrices’ general meaning is that of a gradient operator and that of a ma-
terial parameter matrix at an element level, respectively. For example, if the temperatures
at an element level are given by
d1e
de
2
de = {dea } = .. (5.149)
.
denen
then
nen
q ( x ) = −D ( x ) B ( x ) d e = −D ( x ) ∑ Ba dea (5.150)
a =1
can be used to compute the heat flux within each element. We will revisit this for the
elastic problem where B converts the nodal displacement solution into the strain.
are the same shape functions that are used to represent the solutions
nen
ṽ(ξ ) = ∑ Na (ξ ) dea . (5.154)
a =1
If the mapping from the element local, ξ, to real coordinate space, x, is differentiable, the
determinant j of the Jacobian J
!
∂x ∂x
∂x ∂z ∂z ∂x
j = det J = det ∂ξ
∂z ∂z
∂η
= − (5.155)
∂ξ ∂η ∂ξ ∂η ∂ξ ∂η
is j(ξ ) > 0 for all ξ within the element. j(ξ ) may, in practice, become very small, which in-
dicates that the element is greatly deformed (two edges almost align, for example) which
is to be avoided.
The practical use of j arises from element-local integration. Recall from the 1-D case
(eq. 5.71)
Z Z 1 Z 1
∂x ∂x
f ( x ) dΩ = f ( x (ξ )) dξ = f ( x (ξ )) dξ (5.156)
Ωe −1 ∂ξ −1 ∂ξ
The above equation is a result of the change of variables and can be used to evaluate the
a(., .) type integrals.
In 1-D, the objective is to optimally approximate (i.e. replace the integral over the element
with a summation)
Z 1 nint
−1
dξ g(ξ ) ≈ ∑ g(ξ̃ i ) Wi (5.158)
i =1
for a small number of integration points nint . The Wi are the weights for the function
values at the integration points ξ i . For example, the
Trapezoidal rule corresponds to nint = 2; ξ̃ 1 = −1; η̃2 = 1; Wi = 1 and is second order
accurate.
Z b
f ( a) + f (b)
f ( x ) dx ≈ (b − a) (5.159)
a 2
Gaussian quadrature is the optimal (fewest integration points for maximum accuracy)
strategy. For nint , it is defined by
2
Wi = ∂P 1 ≤ i ≤ nint (5.161)
nint
(1 − ξ̃ i2 ) ∂ξ (ξ̃ i )2
−1
dξ
−1
dη g ≈ ∑ ∑ g(ξ̃ i , η̃j ) Wi Wj . (5.163)
i =1 j =1
Finally, we often need to convert the derivatives of the shape functions with respect to the
global coordinates to local coordinates: ξ ( x). By means of the chain rule we obtain
∂Na ∂Na ∂ξ ∂Na ∂η
= + (5.164)
∂x ∂ξ ∂x ∂η ∂x
∂Na ∂Na ∂ξ ∂Na ∂η
= + . (5.165)
∂z ∂ξ ∂z ∂η ∂z
from which we can compute the derivatives of the coordinates versus natural coordinate
as
!
∂x ∂x
∂ξ ∂η
J = ∂ξ x = ∂z ∂z (5.169)
∂ξ ∂η
∂x ∂Na ∂z ∂Na e
e.g.
∂ξ
= ∑ ∂ξ xea ; ∂η = ∑ ∂η za .
It turns out that the Jacobian J is the inverse of eq. (5.166) (see eq. (6.69)):
! !
∂ξ ∂ξ
− 1 1
∂z
∂η − ∂x
∂η
∂x ∂z
∂η ∂η =J = (5.170)
∂x ∂z
j − ∂ξ∂z ∂x
∂ξ
with
∂x ∂x ∂z ∂x ∂z
j = det = − . (5.171)
∂ξ ∂ξ ∂η ∂η ∂ξ
Therefore
∂Na ∂Na ∂Na ∂Na
, = , J −1 (5.172)
∂x ∂z ∂ξ ∂η
Na (ξ ) = 12 (1 + ξ a ξ ) a = 1, 2
with two nodes at ξ 1 = −1; ξ 2 = 1.
∂ ξa (−1) a
∂ξ Na ( ξ ) = 2 = 2
Figure 5.10: Left: Quad element nodes in local coordinates. Right: Linear shape function
1
Na (ξ ) = Na (ξ, η ) = (1 + ξ a ξ )(1 + ηa η ) (5.176)
4
Quadrature in 2-D
nint visual ξ̃ i wi for integration
R1 R1
1 [×] {0, 0} 4 −1 dξ −1 dη
× × −1 √
2 {√ ; −1 } 1 ”
× × 3 3
1 − 1
{√ ; √ }
3 3
1 ”
{ √13 ; √13 } 1 ”
−1 √1
{√ 3
;
3
1} ”
For higher order quads, see Hughes (2000), sec. 3.7 (Figure 5.11).
Hughes (2000) p. 191 discusses the required level of Gaussian quadrature for adequate
convergence for different element types.
Triangular elements
N1 =r (2r − 1) N4 = 4 r s (5.182)
N2 =s(2s − 1) N5 = 4 s t (5.183)
N3 =t(2t − 1) N6 = 4 r t (5.184)
x= ∑ Na (ξ )xa , (5.185)
where ξ = {ξ, η } or ξ = {r, s}. However, the inverse transformation, going from an
arbitrary global x to the element-local ξ, is a bit more tricky. This conversion involves two
steps, finding the element in which x lies, and then projecting into the local coordinate
system.
Finding an element can be very fast, for example when the mesh consists of regular
quads, in which case this step only involves finding which row i and column j the point
x = { x, z} lies in by means of division of the grid spacings ∆x and ∆z, i = fix( x/(∆x )) + 1
and j = fix(z/(∆z)) + 1, where the element number might be n = i + ( j − 1)n x . Here, n x
is the number of elements in each row, and fix the MATLAB real to integer conversion
that does not round.
For irregular meshes, finding the element containing x can be tricky and care should
be taken to use optimally fast, geometric methods to find if the point lies within the
polygons defined, e.g. by an irregular triangular mesh. MATLAB provides the function
pointLocation for native Delaunay meshes and inpolygon for general meshes, for exam-
ple.
Triangular elements
Once the element is found, the type of projection depends on the type of element. For
triangular, 2D elements, the operation is fairly straightforward. If xi and zi are the nodal
coordinates (i = 1, 2, 3) and x = { x, z}, then
x2 z3 − x3 z2 − x2 z + xz2 + x3 z − xz3
r = (5.186)
x1 z2 − x2 z1 − x1 z3 + x3 z1 + x2 z3 − x3 z2
x z3 − x3 z1 − x1 z + xz1 + x3 z − xz3
s = − 1 (5.187)
x1 z2 − x2 z1 − x1 z3 + x3 z1 + x2 z3 − x3 z2
for the shape functions defined in eqs. (5.179) and (5.180).
Quadratic elements
If the mesh is regular, one can easily convert x to ξ = {ξ, η } given the distance from the
center of the element in question. If the element is irregular, different cases have to be
distinguished and those and the corresponding equations are discussed in Hua (1990).
This FE exercise and most of the following ones are based on the MILAMIN package
by Dabrowski et al. (2008) which provides a set of efficient, 2-D MATLAB-based FE rou-
tines including a thermal and a Stokes fluid solver. Given that the code uses MATLAB,
MILAMIN is remarkably efficient and certainly a good choice for simple 2-D research
problems that lend themselves to FE modeling. You may want to consider working on
expanding the MILAMIN capabilities, e.g. by adding advection to the thermal solver and
combining it with the Stokes solver for a convection code.
Over the next sections, we will discuss all of the issues described in Dabrowski et al.
(2008). This paper will be a good additional reference, and the original MILAMIN MAT-
LAB codes can be downloaded from http://milamin.org/ (the latter will not be of help
with the exercises).
Following, e.g., Hughes (2000), we use the Galerkin approach for which the resulting stiff-
ness matrix components, on an element level, is
Z
e e ∂Na ∂Nb ∂Na ∂Nb
Kab = κ + dΩ. (5.190)
Ωe ∂x ∂x ∂z ∂z
Here, a and b are node numbers local to element e, and integration Ωe is over the element
area.
If we express the spatial coordinates x = { x, z} in a node-local coordinate system ξ =
{ξ, η } and use Gaussian quadrature with N I NT points and weights Wi for integration,
we need to evaluate terms of the kind
Z 1 Z 1
∂Na ∂Nb ∂Na ∂Nb
e
Kab = dξ dηκ + J −1 j (5.191)
−1 −1 ∂ξ ∂ξ ∂η ∂η
N I NT
∂Na ∂Nb ∂Na ∂Nb
Kab = ∑ Wi κi
e
+ J −1 j (5.192)
i
∂ξ ∂ξ ∂η ∂η
where J −1 is the inverse and j = det(J ) = |J | the determinant of the Jacobian matrix
!
∂x ∂x
∂ξ ∂η
J= ∂z ∂z , (5.193)
∂ξ ∂η
respectively.
The load vector F has to be assembled on an element basis as
Z
Fae = e
Na HdΩ − Kab T̂b , (5.194)
Ωe
where the terms on the right hand side are due to heat sources, H, and a correction due to
prescribed temperatures on the boundaries T̂ (zero flux BCs need no specific treatment,
see Hughes, 2000, p. 69 and Dabrowski et al. (2008)). The global K and F are assembled by
looping through all elements and adding up the K e and F e contributions, while eliminat-
ing those rows that belong to nodes where essential boundary conditions (T̂) are supplied.
The solution is then obtained from solving
KT = F. (5.195)
Meshing
Download generate_mesh.m. Start by reading through this MATLAB code, it is a modi-
fication of the MILAMIN wrapper for triangle. Triangle is a 2-D triangular mesh gen-
erator by Shewchuk (2002). That work is freely available as C source code and a flexible,
production quality “Delaunay” mesh generator. A Delaunay mesh is such that all nodes
are connected by elements in a way that any circle which is drawn through the three
nodes of an element has no other nodes within its circumference.
Aside: The “dual graph” (sort of the graphic opposite) of a Delaunay mesh are the
Voronoi cells around each node. Those can be constructed based on the triangulation by
connecting lines that are orthogonal to each of the triangles’ sides and centered half-way
between nodes. Those two properties are important for computational geometry, inverse
theory, and interpolation problems.
A Delaunay triangulation is the best possible mesh for a given number of nodes in the
sense that the triangles are closest to equilateral. For FE analysis, we always strive for
nicely shaped elements (i.e. not distorted from their ideal, local coordinate system form)
so that the J −1 does not go haywire, and j = det(J ) remains positive.
Typically, meshers like triangle will allow you to refine the mesh (i.e. add more nodes)
for a given boundary structure and overall domain by enforcing minimum area and/or
angle constraints. Those refinements may also be iteratively applied based on an initial
solution of the PDE, e.g. to refine in local regions of large variations (adaptive mesh re-
finement, AMR).
Exercise Download a test driver for the triangle wrapper, mesher_test.m. You will have
to fill in the blanks after reading through generate_mesh.m, and make sure the triangle
binary (program) is installed on your machine in the directory you are executing your
MATLAB commands in.
Note: For this exercise and those below, please first inspect graphs on the screen while playing
with the code, and then only print out a few geometries.
(a) Create a triangular grid using three node triangles for the domain 0 ≤ x ≤ 1, 0 ≤
z ≤ 1 using minimum area constraint 0.1 and minimum angle 20◦ . Create a print
out plot of this mesh highlighting nodes that are on the outer boundary.
(c) Use second order triangles and an area constraint of 0.005 and minimum angle of
30◦ .
(d) Using the same quality constraints, create and print out a mesh plot of an elliptical
inclusion of radius 0.2, ellipticity 0.8, and 50 nodes on its perimeter. Color the ele-
ments of the inclusion differently from those of the exterior. Denote nodes on the
boundary of the inclusion.
(e) Create and plot a mesh with a circular hole and a circular inclusion of radius 0.1.
Thermal solver
Download thermal2d_std.m; this is a simplified version of the MILAMIN thermal solver
(thermal2d.m) which should be easier to read than the version of Dabrowski et al. (2008); it
also allows for heat production. Read through this MATLAB code and identify the matrix
assembly and solution method we discussed in sec. 5.6.1. Please make sure you take this
step seriously.
Also download and read through shp_deriv_triangle.m and ip_triangle.m which
implement linear (three node) and quadratic (six node), triangular shape functions and
derivatives, and weights for Gauss quadratures, respectively.
Exercise
(a) Download a rudimentary driver for the mesher and thermal solver,
thermal2d_test.m. You will need to fill in the blanks.
(b) Generate a regular mesh with area constraint 0.003 and solve the heat equation with
linear shape functions, without heat sources, given no flux on the sides, unity tem-
perature at the bottom, and zero temperature at the top. Plot your results. Use
constant conductivity.
(c) Place an elliptical inclusion with radius 0.4, ellipticity 0.8, and ten times higher con-
ductivity than the ambient material in the medium, and plot the resulting tempera-
tures. Experiment with variable resolutions and second order triangles. Comment
on the how the solution changes (visually only is OK).
(d) Set the heat production of the inclusion to 10 and 100, and plot the solution. Com-
pare with boundary conditions where zero temperatures are prescribed on all bound-
ary conditions.
(e) Compute the temperature as well as the geothermal gradient at a specific location.
This exercise is so that you gain experience using shape functions and derivatives
of shape functions and requires you to identify the N and ∂N equivalents.
Figure 5.14: Stress solution for a sheared, elastic box with an inclusion of different strength (see
problem set for details).
This FE exercise is again based on the MILAMIN package by Dabrowski et al. (2008).
Their “mechanical” solver (incompressible Stokes fluid, to be discussed in the next sec-
tion) was rewritten for the elastic problem, and simplified to reduce the dependency on
packages external to MATLAB.
A highly optimized version of the code that, for example, uses matrix reordering for
K is available from us (this one is closer to the original Dabrowski et al. (2008) code). When
inspecting the source codes, you should find many similarities (same mesher, same vari-
able structure, etc.) with last section’s 2-D heat equation exercise.
∇ · σ + f = 0, (5.196)
where σ = σij is the stress tensor and f a body force (such as due to gravity). (Note that
this equation is a general force balance equation in the absence of inertia. You can use it
for static elastic deformation, as we do here, or the Stokes fluid flow problem, as we will
discuss subsequently. The difference arises in the constitutive law.)
Written in component form as PDEs for the finite element domain Ω for each of the
three spatial coordinates i this is
∂σij
+ fi = 0 on Ω (5.197)
∂x j
u i = gi on Γgi (5.198)
σij n j = hi on Γhi . (5.199)
Here, Γh and Γhi , and similar for g, denotes that different components of the traction vector
may be specified on different parts of the domain boundary Γ.
In the case of linear, elastic behavior, the constitutive law linking dynamic with kine-
matic properties is given by the generalized Hooke’s law
Note 1: Notice the definition of the (i,j) derivative short-hand, e.g. operating on u to get
the tensor ε, like u(i,j) .
Note 2: C is a 4th order tensor and somewhat cumbersome to deal with. Noticing that
there are only 6 independent components in σ and ε, we can write the 21 indepen-
dent components of C in the Voigt notation, as a 6 × 6 matrix, CV . However, this
matrix has different definitions (see, e.g. Browaeys and Chevrot, 2004, for a discus-
sion), and is not a tensor anymore. I.e. you can do math with it, such as multiplying
CV · ε6 to get the stress state, where ε6 has the six independent components of ε, but
you cannot rotate CV anymore. For this, the full 4th order C has to be restored.
For an isotropic material, the constitutive law between total stress and strain thank-
fully simplifies to
σij = λε kk δij + 2µε ij = λ∆δij + 2µε ij , (5.202)
where µ and λ are the shear modulus and Lamè parameter, respectively; the latter speci-
fies how incompressible a body is. This formulation uses the isotropic dilation,
3
∆ = ε ii = ∑ ε ii , (5.203)
i =1
2ν νE
λ=µ = with E = 2µ(1 + ν), (5.204)
1 − 2ν (1 + ν)(1 − 2ν)
with the Poisson ratio ν and Young’s modulus E. (There are only two independent
material parameters for isotropic elasticity.) If a block is fixed at the base and loaded
in z-directions without constraints, then ν measures the deformation in the horizontal
ν = −ε xx /ε zz . E measures the stress exerted for the same experiment if the material is not
allowed to give way sideways (free-slip in z direction) by E = σzz /ε zz .
The incompressibility, K, is defined as
with
Z
a(w, u) = dΩ w(i,j) Cijkl u(k,l ) (5.210)
Z
(w, f ) = dΩ wi f i (5.211)
!
3 Z
(w, h) Γh = ∑ Γh
dΓ wi hi . (5.212)
i =1 i
Note that unlike the thermal problem, the solution function we wish to obtain using the
finite element method is a vector, u, rather than a scalar.
Matrix assembly
In the finite element approximation, we then solve for the nodal displacements d which
approximate u within the elements with shape functions N from
Kd = F. (5.213)
where a, b are local node numbers. The elemental force vector at local node a is given by
Z Z
f ie = dΩ Na f i + dΓ Na hi − ∑ k ab gb . (5.215)
Ωe Γhe b
i
We can represent the strain tensor ε as a strain vector e that can be computed from
displacements u by a gradient operator L, like
e = Lu or e j = L jk uk . (5.217)
Within each finite element the displacements can be obtained by summation over the
shape functions for each node a, Na , times the nodal displacements,
u = uk = Na d a = Na dka (5.219)
where d a is the displacement at the local node a, and dka is the k-th spatial component of
this displacement. Then,
(with τxz = 2σxy in analogy to γxy ), then the (symmetric) elasticity matrix D can be used
to obtain stresses from displacements like
s = De = DBa d a . (5.222)
D I J = Cijkl , (5.223)
where e(w) indicates applying the gradient operator to the virtual displacements, as op-
posed to e(u) as in eq. (5.217).
In the isotropic, 2-D plane strain approximation, D takes the simple form
ν
1 0
λ + 2µ λ 0 1− ν
E (1 − ν ) ν
D= λ λ + 2µ 0 = 1− ν 1 0 , (5.226)
( 1 + ν)(1 − 2ν) 0 1−2ν
0 2(1− ν )
0 0 µ
where plane strain means that no deformation is allowed in the y direction, ε yy = 0. For
the case of plane stress, where deformation is allowed and σyy = 0,
λ̄ + 2µ λ̄ 0 1 ν 0
E
D= λ λ̄ + 2µ 0 = ν 1 0 , (5.227)
1 − ν2 1− ν
0 0 µ 0 0 2
with
2λµ
λ̄ = . (5.228)
λ + 2µ
From eq. (5.228), it is apparent that plane stress reduces the effective, volumetric stiffness
of a body, for ν = 1/4, λ̄ = 2/3λ, because out of plane deformation is permitted.
Viscous equivalence
The constitutive law for linear viscous flow with viscosity η, and deviatoric stress σ, σ =
2η ε̇, is analogous to the elastic case with σ = 2µε, assuming the material is incompressible.
The latter can, in theory, be achieved by letting ν → 1/2 for which K/µ → ∞ such that
the linear FE approach can be used to solve simple fluid problems. In practice, however,
special care needs to be taken to allow for the numerical solution of the incompressible
elastic, or the Stokes flow case, which we discuss in sec. 5.8.
Exercises
(a) Make sure you have the MATLAB subroutines ip_triangle.m, shp_deriv_triangle.
m, generate_mesh.m, and the triangle binary from last section in your working di-
rectory. Both shape functions and the mesher will be reused.
(b) Download elastic2d_std.m, a simple linear elasticity solver, and calc_el_D.m which
assembles D. Also download the driver routine elastic2d_test.m. You will have
to fill in the blanks.
(c) Inspect elastic2d_std.m, compare with the notes above for linear elasticity, and the
heat solver from sec. 5.6.
(d) Download and inspect det2D.m, inv2D.m, and eig2d.m (for computing the determi-
nant, inverse, and eigen system of 2 × 2 matrices, respectively). Writing out these
(e) Consider a square, homogeneous elastic body with shear modulus µ = 1, Poisson’s
ratio ν = 1/4 and size 1 × 1 in x and z directions.
Figure 5.15: Load case sketches for some of the exercises, along with common symbols for kine-
matic boundary conditions.
(a) Assume the body is fixed at the base (zero displacement u for all z = 0),
and sheared with a uniform u x displacement of u0 = 0.1 at the top (z = 1)
(Load case a of Figure 5.15a). Assume the plane strain approximation and zero
density (i.e. zero body forces). What kind of geologic deformation state does
this correspond to? What kinds of displacements would you expect, and how
should the major (σ1 ) extensional and the major compressional (σ2 ) stress axis
align?
(b) Compute the displacements and stresses using the 2-D FE programs provided.
Use linear, three node triangles and experiment with the integration order. Use
a coarse mesh with area constraint 0.01 and angle constraint 25◦ .
For this and each subsequent problem, hand in three plots: i), of the deformed
mesh, indicating the shape of the deformed body, possibly exaggerating the
displacements of each node; ii) a plot where the background field (colored) is
the amplitude of displacement, and the foreground has displacement vectors,
plotted with origin at each original node location; and, iii), a plot of mean (nor-
mal) stress (colored in the background), along with extensional and compres-
sional stress axes vector-crosses. The MATLAB routines provided can, with
some alterations, perform all of these tasks.
(c) Compare the predicted stress and displacements for plane strain and plane
stress approximations. Comment.
(d) Compare the distorted mesh shape for linear triangles with that for six node,
quadratic shape functions. Increase the number of elements and compare the
predicted stress fields. Does the displacement and stress field agree with your
expectations for this load case?
(e) Consider Figure 5.15b and prescribe u x displacements linearly tapered from
u x = u0 at z = 1 down to u x = 0 at z = 0. Compare the predicted displace-
ments and stresses with load case Figure 5.15a. Comment on the stress and
displacement fields.
(f) Relax the kinematic boundary conditions on the sides and top and include
body forces with density ρ = 1 at a fixed (no slip) bottom boundary condi-
tion (Figure 5.15c). Compute the displacements and stresses, plot those, and
comment.
(g) Compute the body force load case of Figure 5.15d with free-slip (no horizontal
displacements, u x = 0, and no “vertical” shear stresses, τxz = 0) conditions
on the left and right sides. Compare the stress field with the previous, uncon-
strained case and comment.
(f) Consider the square elastic medium in 2-D plane strain plus a centered, spherical
inclusion with radius 0.2, shear modulus 0.001. Increase the resolution (e.g. use 100
nodes on the outline of the inclusion, 0.001 minimum element area, and 30◦ triangle
edge angle). Load the system as in Figure 5.15b, compute and plot the stress field,
and comment.
(g) Bonus: Write a subroutine that computes the stresses at the global node locations, as
opposed to the integration points within each element as is currently implemented.
Use the nodal stresses and trisurf to generate a plot of triangles colored according
to their normal stress. Compare with the previous plot.
so that the strong form of the incompressible elastic and fluid problems become
∂σij ∂σij
+ fi = 0 + fi = 0 (5.234)
∂x j ∂x j
∂ui ∂vi
=0 =0 (5.235)
∂xi ∂xi
u i = gi v i = gi (5.236)
σij n j = hi σij n j = hi (5.237)
where eqs. (5.234) and (5.235) hold in the domain Ω, eqs. (5.236) and (5.237) are boundary
conditions and hold on Γgi and Γhi respectively. u and v are displacement and velocity,
respectively, and
!
1 ∂u j ∂ui
ε ij = u(i,j) = + (5.238)
2 ∂xi ∂x j
!
1 ∂v j ∂vi
ε̇ ij = v(i,j) = + . (5.239)
2 ∂xi ∂x j
and if there are only displacement/velocity BCs, and f = 0, then pressures are only
determined up to a constant.
Mixed formulation
This is valid both for compressible and incompressible behavior, such that
∂σij
+ fi =0 on Ω (5.244)
∂x j
∂ui p
+ =0 on Ω (5.245)
∂xi λ
ui = gi on Γgi (5.246)
σij n j = hi on Γhi (5.247)
Equations in weak form
nsd Z
p
Z Z Z
∂ui
dΩ w(i,j) σij − dΩ q + = dΩ wi f i + ∑ dΓhi wi hi (5.248)
∂xi λ i
where w, g are virtual displacements and pressures, respectively, and nsd is the number
of dimensions. Special care must be taken in the next step: the choice of shape functions
for pressure and velocities/displacements (see Hughes, 2000, sec. 4.3), but in general, the
pressure shape functions should be lower order (e.g. linear) than the displacements (e.g.
quadratic).
Matrix formulation
K̄ G d F̄
= (5.249)
GT M p H
where the LHS is symmetric but not positive definite (it has negative eigen values or ∼
zero eigen values for pressure modes), p are the pressures at nodes (e.g. center of element),
d are the displacements at nodes (e.g. edges of elements), F̄ are the body forces, and H is
the divergence source in the boundary conditions.
This is called the “segregated d/p form” of the equations and is valid for the general
(including finite compressibility) case.
K̄ ā(w, d) stiffness matrix (symm, pos. def.) (5.250)
G −(∇ w, p) gradient operator (5.251)
GT −(q, ∇ · v) divergence operator (5.252)
p 1
M − q, symm, neg. def. for ν → and M → (5.253)
λ 2
In general, we can distinguish 3 cases:
solve for p
p = M −1 ( H − G T d ) (5.256)
substitute eq. (5.257) into eq. (5.254)
K̄ d + G M −1 ( H − G T d) = F̄
(K̄ − G M −1 G T )d = F̄ − G M −1 H (5.257)
which reduces to solving the following system of equations
Kd=F
where K is symmetric and positive definite and K = K̄ − G M −1 G T and F =
F̄ − G M −1 H from eq. (5.257). As the condition number of K is larger then that
of K̄, we generally need to use direct solvers for this approach.
If p is discontinuous on the elements, eq. (5.257) can be solved locally, on the ele-
ment level.
→ determine p from eq. (5.242).
(G T K̄ −1 G ) p = G T K̄ −1 F − H (5.258)
K p = F (5.259)
Kd= F (5.260)
u( x) = ∑ Na (x) da (5.261)
a
p( x ) = ∑ Ñã (x) pã (5.262)
ã
K ← Ke from element levels (5.263)
F ←fe (5.264)
K e = K̄ e − g e (m e )−1 (g e ) T
f e = f¯ − g e (m e )−1 he
pe = −(m e )−1 (g e ) T d
D here only has deviatoric terms, for the plane strain case
2 0 0
D̄ = µ 0 2 0 ,
0 0 1
Pressure components
1
Z
mãe b̃ = dΩ − Ñã Ñb̃ (5.265)
Ωe λ
Mixed (5.266)
Z
ga ã = − dΩ ∇ · ( Na e) Ñã (5.267)
heã = − ∑ gep ã gep (5.268)
in each element.
p0 = 0
while max (∆pi ) > tolerance
end
Figure 5.16: Pressure and velocity solution for a sinking, fluid slab impinging on viscosity contrast
problem.
This FE exercise is again based on the MILAMIN package by Dabrowski et al. (2008).
As for the heat and elasticity problems, we simplified their “mechanical”, incompressible
Stokes fluid solver to reduce the dependency on packages external to MATLAB.Dabrowski
et al. (2008) have a highly optimized version, which you can obtain from us or the origi-
nal authors; it uses, e.g., reordering of node numbers to improve matrix solutions which
comes an important memory issue for larger problems. The notation here is close to
Dabrowski et al. (2008), for simplicity, but Hughes (2000) has a somewhat clearer exposi-
tion.
secs. 4.9 and 5.7). These approximations transform the general, Navier-Stokes equation
for fluids into the Stokes equation, which is easier to solve, on the one hand, because there
is no turbulence. On the other hand, it is more complicated numerically as Stokes requires
implicit solution methods, whereas turbulent equations can often be solved in an explicit
manner.
The static force-balance equations for body forces due to gravity are given by
∂σij
∇ · σ = f = ρg or = ρgi , (5.273)
∂x j
where σ is the stress tensor, ρ density, and g gravitational acceleration (gi = gδiz ).
We assume that the medium is incompressible and a linear (Newtonian) fluid consti-
tutive law holds,
σij = − pδij + 2η ε̇0ij , (5.274)
where η is the viscosity, p pressure, and ε̇0 the deviatoric strain-rate tensor,
!
0 1 ∂vk 1 ∂vi ∂v j 1 ∂vk 1
ε̇ ij = v(i,j) − δij = + − δij , or ε̇0 = ε̇ − tr (ε̇), (5.275)
3 ∂xk 2 ∂x j ∂xi 3 ∂xk 3
where v are the velocities, and ε̇0 is the total strain-rate reduced by the isotropic part.
We can define
1 1 1
ė = ∑ ε̇ ii = ε̇ ii = tr (ε̇) (5.276)
3 i 3 3
in analogy to the pressure
1
3∑
p=− σii (5.277)
i
such that deviatoric stress and strain-rate are defined from the isotropic quantities as
Using the constitutive law, and assuming 2-D (x-z space), the Stokes equation can be
written as (also see sec. 7)
∂ 4 ∂v x 2 ∂vz ∂ ∂v x ∂vz ∂p
η − + η + − = 0 (5.281)
∂x 3 ∂x 3 ∂z ∂z ∂z ∂x ∂x
∂ 4 ∂vz 2 ∂v x ∂ ∂v x ∂vz ∂p
η − + η + − = ρgz . (5.282)
∂z 3 ∂z 3 ∂x ∂x ∂z ∂x ∂z
Often, we write the constitutive law for deviatoric quantities only,
∂vi
∇ · v = 0 or = 0, (5.284)
∂xi
which allows solving eq. (5.273) for the additional unknown, pressure. For ∇ · v = 0,
but we made the distinction between deviatoric and total strain-rate because we numer-
ically only approximate the incompressible continuity equation, eq. (5.284), by requiring
the divergence to be smaller than some tolerance.
There are several approaches to do this (e.g. penalty methods (sec. 4.9) or Lagrange
methods) which typically involve iterations to progressively introduce additional “stiff-
ness” to the medium (sec. 5.8). We shall allow for a finite, large bulk viscosity, κ, such that
eq. (5.284) is approximated by
∂v x ∂vz p
+ =− , (5.286)
∂x ∂z κ
the right hand side would → 0 for κ → ∞. Eq. (5.286) is valid for the incompressible and
the compressible cases. However, for the compressible case, where the constitutive law,
eq. (5.274), is replaced by
∂v
σij = κ k δij + 2η ε̇ ij , (5.287)
∂xk
p cannot be interpreted as the actual pressure, P = −σii /3, rather it is a pressure parame-
ter because
∂v
P = −(κ + 2η/3) i (5.288)
∂xi
and
∂vi
p = −κ . (5.289)
∂xi
Note: The general, compressible case is identical to the elastic formulation where v → u
and the constitutive law is
∂v
σij = λ k δij + 2µε̇ ij . (5.290)
∂xk
with nsd the number of spatial dimensions. We again use the Galerkin approach, which
leads to the matrix equations.
Matrix assembly
In analogy to the elastic problem, we define a (total) strain-rate vector ė = {ε̇ xx , ε̇ zz , γ̇xz =
2ε̇ xz } such that strain-rates on an element level can be computed from
ė = Bv, (5.296)
where v are velocities given at the element-local nodes, and B holds the derivatives, as
before. When expressed for the local node a and shape functions Na ,
∂Na
0
∂x
Ba = 0 ∂N a .
∂z
(5.297)
∂Na ∂Na
∂z ∂x
Likewise, deviatoric stresses can be computed from t = D ė, where the property matrix D
shall be given by 4 2
3 −3 0
D = η − 23 4
3 0
, (5.298)
0 0 1
for a plane-strain approximation (compare the elastic case). This allows to express the
stress vector with pressure part as
s = − pm + D ė, (5.299)
where m = {1, 1, 0}. The deviatoric-only version of D is
2 0 0
D0 = η 0 2 0 . (5.300)
0 0 1
v( x) ≈ ∑ Na (x)va . (5.301)
a =1
Given the incompressibility constraint, special care has to be taken in the choice of shape
functions, and we will use the seven-node, Crouzeix and Raviart (1973) triangle with quadratic
shape functions Na (cf. Dabrowski et al., 2008).
As detailed in Hughes (2000), one can either choose “conforming” elements for the
problem at hand and get a nice solution for the velocities and pressure right away (which
is what we do here), or choose theoretically inappropriate shape functions and later cor-
rect the pressure (e.g. for so-called “checkerboard modes”). The latter, rough-and-ready
approach may seem less appealing, but works just as well if done properly.
A departure from the elastic problem is that the pressure is treated differently from v,
and we use linear (constant) shape functions for
∇ · v = ε̇ v = Bv ve , (5.303)
and pe = −κBv ve .
The global system of equations for velocity, V , and pressure, P, at the nodes is given
by
A QT
V F
= , (5.304)
Q M P H
where F are the load vectors, e.g. due to body forces, and H is due to the divergence that
may be imposed traction loads for the compressible case (H = 0 for incompressible case).
On an element-level, the stiffness matrix is given by
QT
A
Z
e
kab = dΩ (5.305)
Ωe Q − κ1 M
T
Ba DBb −BvT Ñ T
Z
= dΩ ,
Ωe − ÑBv − κ1 Ña ÑbT
i.e. Q = − ÑBv , M = Ñ Ñ T , and A corresponds to the total stiffness matrix k in the elastic
case. We have omitted the dependence on the local node number in eq. (5.305). Note
that all operations involving Q and M involve the pressure, and not the velocity, shape
functions.
We avoid having to actually solve for the global p by using the “static condensation”.
This means that we locally (element by element) invert M to obtain the pressure from
p ≈ Ña0 p a0 = κ Ñ T M −1 Qve = −κBv ve . (5.306)
A0 V = f , (5.307)
which is to be solved for the nodal velocities V . Here, f = { f e } = {ρe g e } and (the Schur
complement)
A0 = A + κQT M −1 Q. (5.308)
A0 is now symmetric and positive-definite, and the regular, efficient matrix solution meth-
ods can be applied. (Note that the A is only symmetric if the Dirichlet boundary condi-
tions are applied carefully. If implemented straightforwardly, A is not symmetric.)
However, the matrix becomes ill-conditioned (hard to invert) for the desired large val-
ues of κ, which is why iterations for the velocity solution are needed in order to achieve
the incompressibility constraint. Our example code applies “Powell and Hestenes” iter-
ations for the global velocity and pressure vectors V and P (cf. Dabrowski et al., 2008), as
in
P0 = 0, i = 0 (5.309)
while max(∆Pi ) > tolerance
V i = (A0 ) −1 ( f − Q T P i )
∆Pi = M −1 QV i
Pi+1 = Pi + ∆Pi
i = i+1
end while
If and when the algorithm converges, the pressure correction ∆Pi = M −1 QV , which de-
pends on the divergence, M −1 QV , goes to zero. Above, all matrices are meant to be the
global, not element-local representation.
5.9.3 Exercises
(a) Make sure you have the common FE MATLAB subroutines from the earlier exercises
(ip_-triangle.m, shp_deriv_triangle.m, genereate_mesh.m), and the triangle
binary in your working directory.
(c) Compute the sinking velocity of a dense sphere (i.e. disk in 2-D) with radius 0.1 that
is centered in the middle of the 1 × 1 box with free-slip boundary conditions (no
shear stress tangentially to the boundary, no motion perpendicular to the boundary)
on all sides.
Ensure that the sphere is well resolved by choosing ∼ 50 points on its circumference
and using a high quality mesh. Use the second order triangles (six nodes on the
edges plus one added in the center), and six integration points.
(a) Note how boundary conditions are implemented in the MATLAB code, and
comment on essential and natural types.
(b) Compute the solution for the dense sphere with the same viscosity as the back-
ground. Plot the velocities on top of the pressure within the fluid. You may
choose whichever absolute parameter values you like but will have to be con-
sistent subsequently.
(c) Change the number of integration points to three, and replot. Change the type
of element to linear, replot. Comment on the velocity and pressure solution.
(d) The solver applies a finite bulk viscosity (it should be ∞ for an incompress-
ible fluid). For increasing sphere/medium viscosity contrasts upward of 103 ,
experiment with increasing the pseudo-incompressibility and comment on the
stability of the solution. After this experiment, reset to the starting value.
(e) The solver applies iterations to enforce the incompressibility constrain. Change
the tolerance criterion and comment on the resulting velocity and pressure so-
lutions.
(f) Change back to seven node triangles with six integration points. Plot the verti-
cal velocity, vz , along a profile for x ∈ [0; 1] at z = 0.5.
(g) Vary the radius of the sphere and comment on how the vz profiles are affected
by the size of the sinker relative to the box size. How small does the sphere
have to be to not feel the effect of the boundaries?
(h) Change the boundary conditions to no-slip (v = 0 on all domain edges), replot
the vertical velocity profile for a sphere of radius 0.1. Comment. Change back
to free-slip subsequently.
(i) Compute the sinking velocity of a dense sphere with radius 0.1 that is 0.001,
1, and 1,000 times the background viscosity. Define the sinking velocity as the
maximum velocity at the sphere’s origin at x = {0.5, 0.5}.
(j) Provide an analytical estimate for the sinking velocities and compare with the
numerical estimates.
(d) Compute the sinking velocities of a highly elliptical (choose ellipticity 0.975, radius
0.25) body whose viscosity is 1,000 times the background viscosity. Investigate the
case where this “needle” is oriented horizontally (i.e. perpendicular to the sinking
velocity at its center) and when it is oriented vertically (i.e. aligned with the sinking
velocity at its center). Comment on the difference in the maximum sinking velocity
between the two elliptical and the spherical cases.
(e) Bonus (somewhat involved): Compute the sinking velocity for a non-Newtonian, power-
law fluid with ε̇0I I ∝ τInI where n = 3, and I I indicated the second, shear invariants.
Hints: You will have to convert the constitutive law to a viscosity, for which you can
assume constant strain-rates. Then, you will have to modify the code to compute
the strain-rate tensor to obtain the second invariant, ε̇ I I . (You might want to check
the elastic exercise for the use of D and B to obtain strain and stress.) This strain-rate
will then enter the viscosity, and you will have to use a second iteration loop, start-
ing with a Newtonian viscosity, then updating the viscosity from the first velocity
solution, and repeat until velocities do not change by more than some tolerance.
∂
qi = −κij T (heat flux)
∂x j
∂ ∂
ρc p T+ q =H (5.310)
∂t ∂xi i
∂T
(ρc p − k ∇2 T = H for isotropic conductivity) (5.311)
∂t
Boundary conditions
T = g on Γg (essential)
−qi ni = h on Γh (natural)
Initial conditions
T ( x, t = 0) = T0 ( x) (5.312)
Weak form
where we have assumed that the spatial derivatives are now approximated by FE as ex-
pressed by v but time is still left continuous,
Matrix assembly
v( x, t) = ∑ NA ( x ) d A ( t ) (5.319)
with
The main difference with the static sets of equation for the heat equation is the introduc-
tion of the M matrix and the need to solve eq. (5.320) as an ODE.
M d˙ + K d = F (5.325)
with IC d = d0
Note that M, K are symmetric, M is positive definite and K is positive semi-definite (not
pos. def. anymore). A general approach to solve eq. (5.325) is by the generalized trape-
zoidal method (see Hughes, 2000, p. 459).
M v n +1 + K d n +1 = F n +1
dn+1 = dn + ∆t vn+α (5.326)
n+α n n +1
v = (1 − α)v + αv (5.327)
where dn and vn are the approximations to d(t = tn ) and d˙(t = tn ), respectively, with
tn+1 = tn + ∆t (5.328)
as for the finite difference method. For the following α’s the methods in the table below
are recovered.
α
0 forward Euler, fully explicit
0.5 midpoint, Crank-Nicolson
1 backward Euler, fully implicit
v - form implementation
(a) Start at t = t0 with d = d0 given for n = 0.
M v0 = F0 − K d0 (5.329)
Note that for the fully explicit case with α = 0 and a “lumped” M matrix (5.332) (i.e.
diagonal) does not involve any equation solving for time-stepping. M is lumped for ρ
and c p constant.
d - form implementation
Instead of eq. (5.332), we use (for α 6= 0)
1 1
(M + α ∆t K )dn+1 = F n+1 + M d˜n+1 (5.333)
α ∆t α ∆t
to obtain dn+1 , and then update
dn+1 − d˜n+1
v n +1 = (5.334)
α ∆t
The right hand side of eq. (5.333) is fast to compute for diagonal M.
1
The generalized trapezoidal methods are α < 2 conditionally stable for
2
∆t . (5.335)
(1 − 2α) h2
where h is the smallest grid spacing in the mesh (h , “mesh parameter”). For the fully
explicit method (α = 0), we recover
2
∆t ≤ (5.336)
h2
as in the finite difference method.
→ If the complete matrix inversion required for implicit schemes is not feasible, the
element-by-element approach of Hughes (2000) p. 484 (preconditioned conjugate
gradient with Crout factorization) can be used.
Appendix
180
Chapter 6
This section provides a few brief notes on math notation and concepts needed for this text.
Not all concepts and formula are presented in a mathematically rigorous way, and you
should refer to something like a Math for Engineers text for a more complete treatment.
For most of this text, it will be assumed that the reader is familiar with the material treated
in this chapter.
6.1 Calculus
6.1.1 Full and partial derivatives
In calculus, we are interested in the change or dependence of some quantity, e.g. u, on small
changes in some variable x. If u has value u0 at x0 and changes to u0 + δu when x changes
to x0 + δx, the incremental change can be written as
δu
δu = ( x0 )δx. (6.1)
δx
The δ (or sometimes written as capital ∆) here means that this is a small, but finite quan-
tity. If we let δx get asymptotically smaller around x0 , we of course arrive at the partial
derivative, which we denote with ∂ like
δu ∂u
lim ( x0 ) = . (6.2)
δx →0 δx ∂x
The limit in eq. (6.2) will work as long as u doesn’t do any funny stuff as a function of x,
like jump around abruptly. When you think of u( x ) as a function (some line on a plot)
that depends on x, ∂u/∂x is the slope of this line that can be obtained by measuring the
change δu over some interval δx, and then making the interval progressively smaller.
We call ∂u 0
∂x (we also write in shorthand ∂ x u ( x ) or u ( x ); if the variable is time, t, we also
use u̇(t) for ∂u/∂t) the partial derivative, because u might also depend on other variables,
e.g. y and z. If this is the case, the total derivative du at some { x0 , y0 , z0 } (we will drop (i.e.
181
CHAPTER 6. BASIC CALCULUS AND ALGEBRA REVIEW
not write down) the explicit dependence on the variables from now on) is given by the
sum of the changes in all variables on which u depends:
∂u ∂u ∂u
du = dx + dy + dz. (6.3)
∂x ∂y ∂z
Here, dx and similar are placeholders for infinitesimal changes in the variables. This
means that eq. (6.3) works as long as dx is small enough that a linear relationship between
δu and δx still holds. In fact, we can perform a Taylor approximation on any u( x ) around
x0 by
∂u ∂2 u ( x − x0 )2 ∂3 u ( x − x0 )3
u ( x ) = u ( x0 ) + ( x0 )( x − x0 ) + 2 ( x0 ) + 3 ( x0 ) ... (6.4)
∂x ∂x 2! ∂x 3!
2
Here, ∂∂xu2 is the second derivative, the change of the change of u with x. n! denotes the
factorial, i.e.
n! = 1 × 2 × 3 × . . . n. (6.5)
So, as long as dx = x − x0 is small, the derivative will work (for well behaved u). For ex-
ample, if better approximations are needed, e.g.when the strain tensor is not infinitesimal
anymore, quadratic and higher terms like the one that goes with the second derivative
in the series eq. (6.4) and so on need to be taken into account. Finite difference methods
essentially use Taylor approximations to approximate derivatives, as we will see later.
How to compute derivatives Here are some of the most common derivatives of a few
functions:
Sum rule:
( a f ( x ) + bg( x ))0 = a f 0 ( x ) + bg0 ( x ) (6.8)
Product rule:
( f ( x ) g( x ))0 = f 0 ( x ) g( x ) + f ( x ) g0 ( x ) (6.9)
Quotient rule:
g( x )
If f (x) = (6.10)
h( x )
g0 ( x ) h( x ) − g( x ) h0 ( x )
f 0 (x) = (6.11)
h ( x )2
If you need higher order derivatives, those are obtained by successively computing
derivatives, e.g. the third derivative of f ( x ) is
∂3 f ( x )
∂ ∂ ∂
= f (x) . (6.12)
∂x3 ∂x ∂x ∂x
Say, f ( x ) = x3 , then
∂3 x 3
∂ ∂ ∂ 3 ∂ ∂ 2 ∂
3
= x = 3x = 6x = 6. (6.13)
∂x ∂x ∂x ∂x ∂x ∂x ∂x
When applied to scalar field (a distribution of values that depends on spatial location),
such as a temperature distribution T ( x, y, z) (meaning T is variable with coordinates x, y,
and z, assumed implicitly for all properties from now on), the gradient operation
∂T
∂x
∂T
grad T = ∇ T = (6.14)
∂y
∂T
∂z
generates a vector from the scalar field which points in the direction of the steepest in-
crease in T.
Consider what ∇ can do to a vector field (i.e. vectors that vary in space, x). If
is a velocity field, then the divergence (grad dot product) operation on a vector field
div v = ∇ · v (6.16)
is equivalent to finding the dilatancy (volumetric) strain-rate 4̇ from the strain-rate tensor
components because
˙
∆V ∂v1 ∂v2 ∂v3
4̇ =
V
= tr (ε̇) = ∑ ε̇ ii = ε̇ 11 + ε̇ 22 + ε̇ 33 = ∂x1 + ∂x2 + ∂x3 = ∇ · v. (6.17)
i
Eq. (6.17) illustrates that the divergence has to do with sinks and sources, or volu-
metric effects. The volume integral over the divergence of a velocity field is equal to the
surface integral of the flow normal to the surface. (An electro-magnetics example: For the
magnetic field: div B = 0 because there are no magnetic monopoles, but for the electric
field: div E = q, with electric charges q being the “source”.)
If we take the vector instead of the dot product with the grad operator, we have the
curl or rot operation
curl v = ∇ ∧ v. (6.20)
The curl is a rotation vector just like ω. Indeed, if the velocity field is that of a the rigid
body rotation, v = ω ∧ r, one can show that ∇ ∧ v = ∇ ∧ (ω ∧ r ) = 2ω.
Second derivatives enter into the Laplace operator which appears, e.g. in the diffusion
equation:
∂2 T ∂2 T ∂2 T
∇2 T = 2 + 2 + 2 (6.21)
∂x ∂y ∂z
Some rules for second derivatives:
6.1.3 Integrals
Taking an integral Z
F(x) = f ( x )dx, (6.24)
in a general (indefinite) sense, is the inverse of taking the derivative of a function f ,
∂ f (x)
F = f (x) + c (6.25)
∂x
∂ ∂ f (x) ∂
F = ( f ( x ) + c ) = f 0 ( x ). (6.26)
∂x ∂x ∂x
Any general integration of a derivative is thus only determined up to an integration con-
stant, here c, because the derivative, which is the reverse of the integral, of a constant is
zero.
Graphically, the definite (with bounds) integral over f ( x )
Z b
f ( x )dx = F (b) − F ( a) (6.27)
a
along x, adding up the value of f ( x ) over little chunks of dx, from the left x = a to the
right x = b corresponds to the area under the curve f ( x ). This area can be computed by
subtracting the analytical form of the integral at b from that at a, F (b) − F ( a). If f ( x ) = c
(c a constant), then
F(x) = cx + d (6.28)
F (b) = cb + d (6.29)
F ( a) = ca + d (6.30)
F (b) − F ( a) = c ( b − a ), (6.31)
the area of the box (b − a) × c.
Here are the integrals (anti derivatives) of a few common functions, all only deter-
mined up to an integration constant C
x p +1
xp p +1 + C special case: f ( x ) = c = cx0 → F ( x ) = cx + C
ex e x +C
1/x ln(| x |) + C
sin( x ) − cos( x ) + C
cos( x ) sin( x ) + C
There are also a few very helpful definite integrals without closed-form anti derivatives,
e.g.
Z ∞ √
− x2 π
e dx = (6.32)
0 2
Linearity:
Z b Z b Z b
(c f ( x ) + dg( x )) dx = c f ( x )dx + d g( x ) (6.35)
a a a
Reversal: Z b Z a
f ( x )dx = − f ( x )dx (6.36)
a b
Zero length: Z a
f ( x )dx = 0 (6.37)
a
Additivity:
Z c Z b Z c
f ( x )dx = f ( x )dx + f ( x )dx (6.38)
a a b
Product rules:
1
Z
f 0 ( x ) f ( x )dx = ( f ( x ))2 + C (6.39)
Z 2 Z
f 0 ( x ) g( x )dx = f ( x ) g( x ) − f ( x ) g0 ( x )dx (6.40)
Quotient rule:
f 0 (x)
Z
dx = ln | f ( x )| + C (6.41)
f (x)
Gauß theorem The integral over the area Ω of the divergence of a vector field f is equiv-
alent to the boundary integral, ∂Ω, over the local normal (to the boundary), n, dotted with
f: Z Z
dA ∇ · f = ds n · f . (6.42)
Ω ∂Ω
where a and b are vectors of dimension n (n-dimensional, geometrical objects with a di-
rection and length, like a velocity) and the outcome of this operation is a scalar (a regular
number), c. In eq. (6.43), ∑in=1 means “sum all that follows while increasing the index i
from the lower limit, i = 1, in steps of of unity, to the upper limit, i = n”. In the examples
below, we will assume a typical, spatial coordinate system with n = 3 so that
a · b = a1 b1 + a2 b2 + a3 b3 , (6.44)
where 1, 2, 3 refer to the vector components along x, y, and z axis, respectively (ADD
FIGURE HERE). In the “Einstein summation” convention, we would rewrite ∑in=1 ai bi
simply as ai bi , where summation over repeated indices is implied, i.e. the ∑ is not written.
When we write out the vector components, we put them on top of each other
a1 ax
a = a2 = a y (6.45)
a3 az
or in a list, maybe with curly brackets, like so: a = { a1 , a2 , a3 }. Here, we will use a bold
font a to denote vectors as opposed to scalar a, but another common form is ~a, and on a
blackboard, you might also see vectors written as a because that’s easier.
We can write the amplitude (or: length, L2 norm) of a vector as
s
n q q
| a| = ∑ ai = a1 + a2 + a3 = a2x + a2y + a2z .
2 2 2 2 (6.46)
i
For instance, all of the basis vectors defining the Cartesian coordinate system, e x , ey , and
ez have unity length by definition, |ei | = 1. Those ei vectors point along the respective
axes of the Cartesian coordinate system so that we can assemble a vector from its compo-
nents like
a = { a x , ay , az } = a x e x + ay ey + az ez . (6.47)
For a spherical system, the er , eθ , and eφ unity vectors can still be used to express vectors
but the actual Cartesian components of ei depend on the coordinates at which the vectors
are evaluated.
We can restate eq. (6.43) and give another definition of the dot product,
where θ is the angle between vectors a and b. The meaning of this is that if you want to
know what component of vector a is parallel to b, you just take the dot product. Say, you
have a velocity v and want the normal velocity vn along a vector n with |n| = 1 that is
oriented at a 90◦ angle (perpendicular) to some plate boundary, you can use vn = v · n.
Also, eq. (6.47) only works because the basis vectors ei of any coordinate system are,
by definition, orthogonal (at right angle, perpendicular, at θ = 90◦ ) to each other and
ei · e j = 0 for all i 6= j. Likewise, ei · ei = 1 for all i since a · a = | a|2 , and basis vectors
have unity length by definition. Using the Kronecker δ
ei · e j = δij . (6.50)
c = a∧b (6.51)
that is at a right angle to both a and b (hence the right-hand-rule, with thumb, index, and
middle finger along a, b, and c, respectively). vector c’s length is given by
that is, c is largest when a and b are orthogonal, and zero if they are parallel. Compare
this relationship to eq. (6.48).
In 3-D,
a y bz − a z by
c = a ∧ b = a z b x − a x bz (6.53)
a x by − a y b x
(note that there is no i component of a or b in the i component of c, this is the aforemen-
tioned orthogonality property).
An example for a cross product is the velocity v at a point with location r in a body
spinning with the rotation vector ω, v = ω ∧ r. The rotation vector ω is different from,
e.g., r in that ω has a spin (a sense of rotation) to it (the other right-hand-rule, where your
thumb points along the vector and your fingers indicate the counter-clockwise motion).
ADD FIGURE
C = AB (6.58)
cij = ∑ aik bkj , (6.59)
k
where k goes from 1 to the number of columns in A, which has to be equal to the number
of rows in B. Note that, in general, AB 6= BA!
Quadratic matrices Have n × n rows and columns. All simple physical tensors, such as
stress or strain, can be written as quadratic matrices in 3 × 3.
Identity matrix = I, iij = δij , i.e. this matrix is unity along the diagonal, and zero for
all other elements.
Vector cross product based on the determinant The cross product c = a ∧ b (eq. 6.53)
can also be written as the determinant of the matrix
e x ey ez
a x ay az (6.63)
b x by bz
I IA = a11 a22 + a11 a33 + a22 a33 − a212 − a213 − a223 . (6.66)
These expressions arise when finding the eigenvectors and values of a tensor, eq. (6.73).
Transpose of a matrix (AT )ij = aijT = a ji , i.e. the transpose has all elements flipped by
row and column.
If the inverse exists, then (A−1 )−1 = A, (AT )−1 = (A−1 ) T , and (AB )−1 = B −1 A−1 . The
inverse only exists if det(A) 6= 0.
For the special case of a 2 × 2 matrix
a b
A= , (6.68)
c d
AAT = AT A = (6.70)
holds.
Eigenvalues and eigen vectors: Any n × n symmetric matrix A has n eigen vectors vi
that correspond to real eigenvalues λi such that
Avi = λi vi (6.71)
An example is the stress matrix which can be written in the principal axes system, where
the eigen vectors of the Cartesian representation of the stress matrix are the principal axes.
Eigenvalues can be found using
det(A − λI ) = 0 (6.72)
and eigen vectors subsequently by using the first property, which leads to
If a symmetric matrix A is transformed into the principal axes system, A0 , there are no
off-diagonal elements
a xx a xy a xz a1 0 0
A = ayx ayy ayz → A0 = 0 a2 0 (6.74)
azx azy azz 0 0 a3
where the a1 , a2 , and a3 correspond to the three eigenvalues λi . (The coordinate system
reference of A0 is then contained in the orientation of the eigen vectors vi .) For a matrix in
the principal axis system, the invariants are very easily computed:
See also sec. 7.2 for definitions of invariatns using deviators, such as for the deviatoric
stress tensor.
6.2.4 Tensors
The stress σ and strain ε are examples of second order (rank r = 2) tensors which, for
n = 3, 3-D operations, have 3r components and can be written as n × n matrices. You will
see tensors printed like so E, and the blackboard version again double underlining like ε,
making no distinction between tensors and matrices.
Tensors in a Cartesian space are defined by their properties under coordinate transfor-
mation. If a quantity v remains intact under rotation to a new coordinate system v0 such
that
3
vi0 = Lij v j = ∑ Lij v j (6.78)
j =1
holds, then v, a vector, is a first order tensor. Lij may be, for example, a rotation ma-
trix. Likewise, a second order tensor T is defined by remaining intact after rotation into
another coordinate system where it is expressed as T 0 such that
We will assume some familiarity with continuum mechanics as discussed in the context
of an introductory geodynamics course; a good reference for such problems is Turcotte and
Schubert (2002). However, here is a short and extremely simplified review of basic contin-
uum mechanics as it pertains to the remainder of the class. You may wish to refer to our
math review if notation or concepts appear unfamiliar, and consult chap. 1 of Spiegelman
(2004) for some clean derivations.
193
CHAPTER 7. CONTINUUM MECHANICS REVIEW
• Repeated indices indicate summation over these components (also called Einstein
summation convention).
3
∂vi ∂v ∂v1 ∂v2 ∂v3
∂xi
with i = 1, 2, 3 implies ∑ ∂xii = + +
∂x1 ∂x2 ∂x3
(7.1)
i =1
• In a Eulerian frame one uses a reference system for computations that is fixed in
space, for example a computational box in which we solve for advection of tem-
perature T in a velocity field v. Local changes in, e.g., temperature are then given
by
DT ∂T ∂T ∂T
= + v∇ T = + vi , (7.2)
Dt ∂t ∂t ∂xi
where D/Dt is the total derivative that we would experience if we were to ride on
a fluid particle in the convection cell (Lagrangian reference frame). D/Dt takes into
account local changes in a property with time (e.g. due to radioactive heating for T)
as well as advection of temperature anomalies by means of v in and out of our local
observation point.
σ11 σ12
σij = (2D) (7.3)
σ21 σ22
σ11 σ12 σ13
σij = σ21 σ22 σ23 (3D). (7.4)
σ31 σ32 σ33
• Meaning of the elements: Each row are components of the traction vectors acting on
the coordinate plane normal to the respective coordinate axis, the diagonal elements
are normal stresses, and off-diagonal elements are shear stresses.
σij : force/area (traction) on the i plane (plane with normal aligned with the i-th
coordinate axis) along the j direction.
• Special properties: Symmetric, i.e. σij = σji . This means that only six components
of σ need to be stored during computations since the other three can be readily
computed. Note: There are different convention for the order of storing elements
of σ (e.g. diagonal elements first, then off-diagonal; alternatively, upper right hand
side ordering within, for example, a finite element program).
• Cauchy’s formula: if you multiply the stress tensor (treated as a matrix) by a unit
vector, n j , which is normal to a certain plane, you will get the traction vector on this
plane (see above):
3
T(n)i = σn = σij n j = ∑ σij n j (7.5)
j =1
• In a model, the stress tensor is usually computed by solving the equilibrium equa-
tions.
Note: The number of equilibrium equations is less than the number of unknown
stress tensor components.
• Stress deviator
The deviatoric stress tensor invariants of τ are typically denoted as J (as opposed to
!
1 ∂ui ∂u j
ε ij = + (7.14)
2 ∂x j ∂xi
!
1 ∂vi ∂v j
ε̇ ij = + (7.15)
2 ∂x j ∂xi
∂v1 1 ∂v1 ∂v2 1 ∂v1 ∂v3
∂x1 2 ∂x2 + ∂x1 2 ∂x3
+ ∂x1
1 ∂v2 ∂v1 ∂v2 1 ∂v2 ∂v3
= + ∂x3 + (7.16)
2 ∂x1 ∂x2 ∂x2 2 ∂x2
1 ∂v3 ∂v1 1 ∂v3 ∂v2 ∂v3
2 ∂x1 + ∂x3 2 ∂x2 + ∂x3 ∂x3
• Note: The number of velocity components is smaller than the number of strain rate
components.
Here, λ, µ are elastic moduli (for an isotropic medium, there are two (bulk and shear)
independent moduli which can be related to all other commonly used parameters
such as Poisson’s ratio). η is (dynamic, shear) viscosity, bulk viscosities are usually
assumed infinite. Sometimes, kinematic viscosity ν = η/ρ is used.
• To solve a problem starting from the equilibrium equations for force balance, one
can replace stress by strain (rate) via the constitutive law, and then replace strain
(rate) by displacement (velocities). This results in a “closed” system of equations
in “fundamental” variables, meaning that the number of equations is equal to the
number of unknowns, the basic displacements (velocities).
• Material parameters for solid Earth problems can ideally be obtained by measuring
rheology in the lab. Alternatively, indirect inferences from seismology or geody-
namic modeling augmented by constraints such as post-glacial rebound need to be
used.
– Reversible elastic rheology at small stresses and strains over short time scales.
– Irreversible fluid flow (creep) at large strains and over long time scales. Ex-
amples are Newtonian viscous (rate-independent) or power-law (rate/stress
dependent) rheology; usually thermally activated. Intermediate stress levels.
– Rate-independent (instantaneous), catastrophic yielding at large, limit stresses.
Pressure sensitive, often temperature independent. Also called plastic, or fric-
tional (brittle), behavior. Important for cold material over long time-scales.
Conservation of energy
!
∂E ∂E ∂qi
+ vj + = ρQ (7.21)
∂t xj ∂xi
where E is energy, qi the energy flux vector, and Q an energy source (heat production).
E = c p ρT (7.24)
where c p is heat capacity, and T is temperature. If all material parameters are constant
(homogeneous medium), we can then write conservation of energy as
DT ∂T
= + v · ∇ T = κ ∇2 T + H (7.25)
Dt ∂t
or
∂T ∂T ∂ ∂ ∂2 T
+ vj =κ T+H = κ∑ 2 +H (7.26)
∂t ∂x j ∂xi ∂xi i ∂xi
with H = Q/ρ and the thermal diffusivity
k
κ= . (7.27)
ρc p
Equation of state 2: relationships for the isotropic parts of the stress/strain tensors
ρ = f ( T, p) (7.28)
∂ρ ∂ρvi
+ = 0 (7.30)
∂t ∂xi
!
∂vi ∂vi ∂σij
ρ + vj = + ρgi (7.31)
∂t xj ∂x j
!
∂E ∂E ∂q
+ vj + i = ρQ (7.32)
∂t xj ∂xi
E = c p ρT (7.33)
ρ = f ( T, P) (7.34)
ε̇˜ ij = R(σ̇˜ ij , σ̃ij ) (7.35)
∂T
qi = − k (7.36)
∂xi
where ρ is density, vi velocity, gi gravitational acceleration vector, E energy, qi heat flux
vector, Q an energy source (heat production, e.g. by radioactive elements), c is heat capac-
ity, T temperature, p pressure and k thermal conductivity. R indicates a general constitu-
tive law.
∂vi
=0 (7.37)
∂xi
∂σij
+ ρ 0 gi = 0 (7.38)
∂x j
!
∂T ∂T ∂ ∂T
ρ0 c p + vj = k + ρ0 Q (7.39)
∂t xj ∂xi ∂xi
σ̃ij
ε̇˜ ij = (7.40)
2η
σij = − pδij + σ̃ij (7.41)
Major simplifications: No inertial (Dρ/Dt) terms (infinite Prandtl number, see non-
dimensional analysis), incompressible flow, linear viscosity.
∂v x ∂vz
+ = 0 (7.42)
∂x ∂z
∂σxx ∂σxz
+ = 0 (7.43)
∂x ∂z
∂σxz ∂σzz
+ − ρg = 0 (7.44)
∂x ∂z
∂v x
σxx = − p + 2η (7.45)
∂x
∂vz
σzz = − p + 2η (7.46)
∂z
∂v x ∂vz
σxz = η + (7.47)
∂z ∂x
2
∂ T ∂2 T
∂T ∂T ∂T
ρ0 c p + vx + vz = k + 2 + ρ0 Q (7.48)
∂t ∂x ∂z ∂x2 ∂z
Introduction to MATLAB
Reading
8.1 Introduction
MATLAB is commercial software that provides a computing environment that allows for
sophisticated ways of developing and debugging computer code, executing programs,
and visualizing the output. MATLAB is also a computer language (sort of a mix between
C and Fortran) and this exercise for you to work through is mainly concerned with some
of the language aspects that we will use extensively throughout the book. Please read
through the more comprehensive and verbose MATLAB Intro and familiarize yourself
with MATLAB.
We will assume that your Windows, Mac, or Linux machine has MATLAB installed.
After starting up the program with the graphical user interface enabled, you will be pre-
sented with a number of windows, including an interactive window (“shell”) where you
can type in commands as we indicate below. Please also familiarize yourself with the
other components of the development environment, such as the built-in editor for MAT-
LAB programs, which are called “m-files”, so that you can be more efficient in writing
and debugging codes. There are numerous MATLAB -provided help resources accessible
through the environment, including video tutorials, access to the help pages, along with
extensive documentation on the web.
Also note that there is a free clone of MATLAB called octave. Given that MATLAB of-
ten uses freely available computational routines underneath the hood, it was fairly easy
to reproduce the computational basics of MATLAB.However, the MATLAB people also
added a bunch of proprietary visualization tools which are not available in octave. An-
other alternative is to use the freely available Python language and its MatPlotLib pack-
age, but we will not have time to explore such intriguing options.
202
CHAPTER 8. INTRODUCTION TO MATLAB
Matrix-vector multiplication:
5
1 4 5 10 = 130
D T bT =
2 3 6 142
17
Matrix-matrix multiplication:
1 2
1 4 5 4 3 = 42 44
DT D =
2 3 6 44 49
5 6
If you don’t know what’s going on here, and what the rules for such multiplications
are, please consult sec. 6.
Ac = Rhs
where A is a n × m matrix and Rhs is a n × 1 vector whose coefficients are both known,
and c is a m × 1 vector with unknown coefficients. If we take A = D and Rhs = b T , c is
(check!):
1
c =
2
2. Type the following commands and note how MATLAB deals with vectors
>>x=3
>>x=3;
>>x
>>y=x^2
>>x = [2, 5.6]
>>y=2 * x;
>>y=x^2;
>>y=x.^2
>>y = [3, 4]
>>x * y
>>x * y’
>>x .* y
>>pi
>>a=x*pi
3. Type demo and explore some examples. Also note the introductory tutorial videos
you might want to watch later.
4. Type help. You see a list of all help functions. Type help log10 to get information
about the log10 command. Type help logTAB where logTAB means typing log and then
pressing the TAB key without adding a white space. Notice the command completion
selection within the MATLAB shell. Note also that you can use the Up and Down arrows
to retrieve previous commands and navigate through your command history, and pUP will
bring up the last command line that started with a p. MATLAB also offers a graphical user
interface (GUI) to explore all of its features: click help in the menu bar, then product help.
Moreover, the function browser offers you a graphical way to find the suitable function
for what you are trying to accomplish. The function browser can be found under help,
functions browser, or can be brought up using a keyboard shortcut: Shift+F1.
>>x*y.’
>>x’*y
>>ones(1,5), zeros(6,1)
>>length(x)
>>whos
>>[x2d,y2d] = meshgrid(0:.1:2*pi,1:.1:2*pi)
>>size(x2d)
>>z2d = sin(x2d.*y2d)
>>surf(x2d,y2d,z2d)
>>mesh(x2d,y2d,z2d)
>>contour(x2d,y2d,z2d), colorbar
>>contourf(x2d,y2d,z2d), colorbar
>>[x2d,y2d,z2d] = peaks(30);
>>surf(x2d,y2d,z2d); shading interp
>>light; lighting phong
Some cool stuff (2): perform the example given at the end of
>>help coneplot;
12. Use the built in editor (or another text editor e.g. Emacs) and create a file “my-
surf.m”.
13. Type the plotting commands from the last section in the text file. A good program-
ming convention is to start the script with clear, which clears the memory of MATLAB.
Another good programming practice is to put lots of comments inside a MATLAB script.
A comment can be placed after %, e.g. % this is my first MATLAB script.
14. Start the script from within MATLAB by going to the directory where the text file
is saved. Type mysurf from within MATLAB and you should see the plot pop up in a new
figure window. Alternatively, within the MATLAB editor, you can press F5 to run. Also
note that there are various debugging features in the editor that are very helpful, such as
real-time syntax checking and addition of breakpoints.
8.3.5 Loops
Create an array na=100; a=sin(5*[1:na]/na); plot(a).
15. Ask instructions on using ”for”:
>>help for
16. Compute the sum of an array:
>>mysum=0; for i=1:length(a), mysum = mysum + a(i); end; mysum
17. Compare the result with the MATLAB inbuilt function sum
>>sum(a)
18. Exercise. Create x-coordinate array: dx=0.01; y=cos([0:dx:10]). Compute the inte-
gral of y=cos(x) on the x-interval 0 < x < 10. Use sum(y) and write a MATLAB -script.
Compare it with sin(10), the analytical solution.
21. Compare the results with the built in MATLAB function cumsum:
>>bednumber=1:length(depth)
>>plot(bednumber,depth,bednumber,cumsum(thickness))
22. What causes the discrepancy? Try to remove it, ask help cumsum
8.3.7 IF command
23. Ask help if. Find maxima of the above array thickness, and compare it with the in
built function max(thickness)
25. Find the number of beds with a maximum thickness less than 0.5.
8.3.10 Functions
MATLAB allows you to declare functions that return a value and use m-files to store those
functions. If you save
function xs = mysqr(x)
xs = x.^2;
as a mysqr.m in your working directory, you can then use your function just like a regular
MATLAB command.
y=[2,3,4]
mysqr(y);
The benefit of this is that you can now, for example, pass “quake” to functions and the
function will locally know that quake actually has the components lon, lat, and depth
which can be addressed within the subroutine.
26. Exercise: Write and test function that has two inputs, x and a polynomial. The
polynomial structure should have two entries, the order of the polynomial expansion n
and a vector a with n entries that hold the coefficients such that the function returns
n
y= ∑ a i x n −1 (8.1)
i =1
As an example for how the textbook can be used for a one semester course, we provide
the syllabus as Numerical Geodynamics was taught at the University of Southern Cal-
ifornia as GEOL540 in the Fall of 2008. Each week, the class met for a three hour slot
which typically consisted of some formal instruction by means of lectures and joint mat-
lab problem-set exercises in a computer lab. Some weeks, all class time is spent working
on problem sets. The class culminates in a three week final project part where students are
to either write their own code or combine codes used in class (for example, combine ad-
vection and diffusion solvers, and further with a Stokes solver to arrive at a self-contained
convection code).
210
CHAPTER 9. EXAMPLE SYLLABUS AS USC GEOL540 – 2008
(k) Finite elements IV (sec. 5.7): Compressible elastic problems. Elastic moduli, plane
stress, plane strain. Gradient operator, elasticity matrix, engineering strain convec-
tion. Visualization of stress states, eigensystems.
Problem set: 2-D FE elastic
(l) Finite elements V (sec. 5.9 & 5.8): Compressible and incompressible elasticity and
Stokes flow. Mixed formulation with discontinuous pressure. Powell-Hestenes iter-
ations.
Notes: Incompressible elastic/fluid problem Problem set: 2-D FE incompressible
Stokes
Becker, T. W. (2000), Deterministic chaos in two state-variable friction sliders and the ef-
fect of elastic interactions, in GeoComplexity and the physics of earthquakes, Geophys. Mono-
graph, vol. 120, edited by J. B. Rundle, D. L. Turcotte, and W. Klein, pp. 5–26, American
Geophysical Union, Washington, DC.
Becker, T. W., and A. Braun (1998), New program maps geoscientific data sets interac-
tively, Eos Trans. AGU, 79, 505.
Becker, T. W., and C. Faccenna (2009), A review of the role of subduction dynamics for re-
gional and global plate motions, in Subduction Zone Geodynamics, edited by F. Funiciello
and S. Lallemand, Int. J. Earth Sci., pp. 3–34, Springer.
Becker, T. W., and C. Faccenna (2011), Mantle conveyor beneath the Tethyan collisional
belt, Earth Planet. Sci. Lett., 310, 453–461.
Becker, T. W., and B. Schott (2002), On boundary-element models of elastic fault interac-
tion (abstract), Eos Trans. AGU, 83(47), NG62A–0925.
Boschi, L., and A. M. Dziewoński (1999), ‘High’ and ‘low’ resolution images of the Earth’s
mantle – Implications of different approaches to tomographic modeling, J. Geophys. Res.,
104, 25,567–25,594.
Briggs, W. L., V. E. Henson, and S. F. McCormick (2000), A multigrid tutorial, 2 ed., The
Society for Industrial and Applied Mathematics.
Browaeys, J., and S. Chevrot (2004), Decomposition of the elastic tensor and geophysical
applications, Geophys. J. Int., 159, 667–678.
Carslaw, H. S., and J. C. Jaeger (1959), Conduction of Heat in Solids, 2nd ed., Oxford Uni-
versity Press, London, p. 243.
213
BIBLIOGRAPHY
Christensen, U. R. (1985), Thermal evolution models for the Earth, J. Geophys. Res., 90,
2995–3007.
Clayton, R., and H. Engquist (1977), Absorbing boundary conditions for acoustic and
elastic wave equations, Bull. Seismol. Soc. Am., 67, 1529–1540.
Collino, F., and C. Tsogka (2001), Application of the perfectly matched absorbing layer
model to the linear elastodynamic problem in anisotropic heterogeneous media, Geo-
physics, 66, 294–307.
Crouch, S. L., and A. M. Starfield (1983), Boundary Element Methods in Solid Mechanics. With
Applications in Rock Mechanics, Allen and Unwin, London.
Crouzeix, M., and P. A. Raviart (1973), Conforming and nonconforming finite elements
methods for solving the stationary Stokes equation, Rev. Franc. d’Automat. Informat.
Rech. Opér., 3, 33–76.
Dahlen, F. A., and J. Tromp (1998), Theoretical Global Seismology, Princeton University
Press, Princeton, New Jersey.
Deubelbeiss, Y., and B. J. P. Kaus (2008), Comparison of Eulerian and Lagrangian numer-
ical techniques for the Stokes equations in the presence of strongly varying viscosity,
Phys. Earth Planet. Inter., 171, 92–111.
Foley, B., and T. W. Becker (2009), Generation of plate tectonics and mantle hetero-
geneity from a spherical, visco-plastic convection model, Geochem., Geophys., Geosys.,
10(Q08001), doi:10.1029/2009GC002378.
Gerya, T. V., and D. Yuen (2003), Characteristics-based marker-in-cell method with conser-
vative finite-differences schemes for modeling geological flows with strongly variable
transport properties, Phys. Earth Planet. Inter., 140, 293–318.
Golub, G. H., and C. F. Van Loan (1996), Matrix computations, 3 ed., Johns Hopkins Uni-
versity Press.
Gu, J.-C., J. R. Rice, A. L. Ruina, and S. T. Tse (1984), Slip motion and stability of a single
degree of freedom elastic system with rate and state dependent friction, J. Mech. Phys.
Solids, 32, 167–196.
Hager, B. H., and R. J. O’Connell (1981), A simple global model of plate dynamics and
mantle convection, J. Geophys. Res., 86, 4843–4867.
Ismail-Zadeh, A., and P. Tackley (2010), Computational Methods for Geodynamics, Cam-
bridge University Press.
Jacoby, W. R., and H. Schmeling (1981), Convection experiments and driving mechanism,
Geol. Rundschau, 24, 217–284.
Jaupart, C., S. Labrosse, and J.-C. Marechal (2007), Temperatures, heat and energy in the
mantle of the Earth, in Treatise on Geophysics, edited by G. Schubert and D. Bercovici,
pp. 253–303, Elsevier.
King, S. D., D. A. Raefsky, and B. H. Hager (1990), ConMan: vectorizing a finite element
code for incompressible two-dimensional convection in the Earth’s mantle, Phys. Earth
Planet. Inter., 59, 195–207.
Korenaga, J. (2008), Urey ratio and the structure and evolution of Earth’s mantle, Rev.
Geophys., 46, doi:10.1029/2007RG000241.
Kronbichler, M., T. Heister, and W. Bangerth (2012), High accuracy mantle convection
simulation through modern numerical methods, Geophys. J. Int., 191, 12–29.
Kwon, Y. W., and H. Bang (1996), The Finite Element Method Using Matlab, CRC Press.
Lay, T., J. Hernlund, and B. Buffett (2008), Core-mantle boundary heat flow, Nature Geosc.,
1, 25–32.
Lenardic, A., and W. M. Kaula (1993), A numerical treatment of geodynamic viscous flow
problems involving the advection of material interfaces, J. Geophys. Res., 98, 8243–8260.
Marone, C. (1998), Laboratory-derived friction laws and their application to seismic fault-
ing, Annu. Rev. Earth Planet. Sci., 26, 643–696.
May, D. A., and L. Moresi (2008), Preconditioned iterative methods for Stokes flow prob-
lems arising in computational geodynamics, Phys. Earth Planet. Inter., 171, 33–47.
Moresi, L. N., F. Dufour, and H.-B. Mühlhaus (2003), A Lagrangian integration point finite
element method for large deformation modeling of viscoelastic geomaterials, J. Comp.
Phys., 184, 476–497.
Mozco, P., J. Kristek, and L. Halada (2004), The finite-difference method for seismologists. An
introduction., Comenius University, Bratislava, online at http://www.spice-rtn.org/
events/workshops/venice2004/downloads/spicefdmcourse.tgz.
Okada, Y. (1992), Internal deformation due to shear and tensile faults in a half-space, Bull.
Seismol. Soc. Am., 82, 1018–1040.
Paige, C. C., and M. A. Saunders (1982), LSQR: an algorithm for sparse linear equations
and sparse least-squares, Trans. Math. Software, 8, 43–71.
Samuel, H., and M. Evonuk (2010), Modeling advection in geophysical flows with particle
level sets, Geochem., Geophys., Geosys., in press, doi:10.1029/2010GC003081.
Schmeling, H. (1989), Compressible convection with constant and variable viscosity: The
effect on slab formation, geoid, and Topography, J. Geophys. Res., 94, 12,463–12,481.
Schmeling, H. (1994), Skriptum: Numerische Methoden in der Geophysik, Institut für Meteo-
rologie und Geophysik, Universiät Frankfurt am Main.
Schubert, G., D. Stevenson, and P. Cassen (1980), Whole planet cooling and the radiogenic
heat source contents of the Earth and Moon, J. Geophys. Res., 85, 2531–2538.
Schubert, G., D. L. Turcotte, and P. Olson (2001), Mantle Convection in the Earth and Planets,
Cambridge University Press.
Smolarkiewicz, P. K. (1983), A simple positive definite advection scheme with small im-
plicit diffusion, Mon. Weather Rev., 111, 479–486.
Spencer, R. L., and M. Ware (2008), Introduction to Matlab, Brigham Young University,
available online, accessed 07/2008.
Spiegelman, M. (2004), Myths and Methods in Modeling, Columbia University Course Lec-
ture Notes, available online at http://www.ldeo.columbia.edu/~mspieg/mmm/course.
pdf, accessed 06/2006.
Suckale, J., J.-C. Nave, and B. H. Hager (2010), It takes three to tango 1: Simulating
buoyancy-driven flow in the presence of large viscosity contrasts, J. Geophys. Res., in
press, doi:10.1029/2009JB006916.
Tackley, P. J., and S. D. King (2003), Testing the tracer ratio method for modeling active
compositional fields in mantle convection simulations, Geochem., Geophys., Geosys., 4,
doi:10.1029/2001GC000214.
Turcotte, D. L., and G. Schubert (2002), Geodynamics, 2 ed., Cambridge University Press,
Cambridge.
van Keken, P. E., S. King, H. Schmeling, U. Christensen, D. Neumeister, and M.-P. Doin
(1997), A comparison of methods for the modeling of thermochemical convection, J.
Geophys. Res., 102, 22,477–22,495.
Weijermars, R., and H. Schmeling (1986), Scaling of Newtonian and non-Newtonian fluid
dynamics without inertia for quantitative modelling of rock flow due to gravity (in-
cluding the concept of rheological similarity), Phys. Earth Planet. Inter., 43, 316–330.
Wessel, P., and W. H. F. Smith (1998), New, improved version of the Generic Mapping
Tools released, Eos Trans. AGU, 79, 579.
Zhong, S., and B. H. Hager (2003), Entrainment of a dense layer by thermal plumes, Geo-
phys. J. Int., 154, 666–676.
Zhong, S. J., D. A. Yuen, and L. N. Moresi (2007), Numerical methods in mantle convec-
tion, in Treatise in Geophysics, vol. 7, edited by G. Schubert and D. Bercovici, pp. 227–252,
Elsevier.
219
INDEX