The Development of a Volume-of-Fluid Interface

Tracking Method for Modeling Problems in

Mantle Convection

Jonathan Robey1

August 27, 2016

1
Funded in part by a GAANN Fellowship
1 Summary
In problems involving mantle convection, it is often important to track the move-
ment of a conserved fluid or material volume. Although in some cases this can
be done by modeling the advection equation for a compositional variable, often
the smearing/diffusion introduced by standard numerical methods for modeling
the advection equation obscures or eliminates essential information, such as the
material interface. In our work, we implement a Volume-of-Fluid (VoF) inter-
face tracking method in a large open source finite element code, ASPECT [4],
for modeling convection in the Earth’s mantle.
In section 2 we briefly discuss the applications of our method to geodynamic
problems. In section 3, the we describe the equations governing convection in
the Earth’s mantle and the numerical methodology we use to model them. In
section 4, we describe in detail the VoF methodology. In section 5, we present
numerical results including computations to validate the accuracy of the newly
implemented VoF algorithm in ASPECT and show applications of this new
algorithm to a problems in mantle convection. In section 6, we discuss future
directions for this work.

2 Geodynamics Background
Mantle convection is one of the primary areas of interest in the field of geo-
dynamics. Due to the nature of the problem, computational modeling is an
essential tool for examining wide variety of problems that are associated with
mantle convection. The underlying problem, thermally driven convection (i.e.
the Rayleigh-Bernard equations), is well established problem in fluid dynamics.
In a number of situations the ability to carefully track a fluid or material volume
is highly desirable, due to variations in fluid properties or interest in location of
the fluid interface. This suggests the use of an interface tracking method.

3 Computational Background
3.1 Introduction
Over the past forty years, a number of codes have been developed for the pur-
pose of modeling convection in the Earth’s mantle. ASPECT, a finite element
code designed for extensibility, is currently among the best codes in the field.
However, it does not have any interface tracking ability, other than with tracer
particles. Furthermore, to the best of our knowledge no other state-of-the-art
mantle convection codes have interface tracking capability. Thus, the addition
of a VoF interface tracking extension to ASPECT will be of significant value to
researchers in the mantle convection community.
The VoF method is an interface tracking algorithm designed to use the ad-
vection of an indicator or ‘color’ function to approximate the location of the
interface between two fluid volumes. The core of the method is similar to a

1
finite volume advection scheme, though specialized for to tracking the inter-
face of an indicator or color function. As such, VoF methods are inherently
conservative; i.e. their basic design will conserve the material being tracked.

3.2 An Overview of Current Interface Tracking Methods

In this section we discuss a general overview of several standard interface track-
ing methods.
The interface tracking problem is found in a wide variety of applications in
science and technology, and as such a number of approaches to modeling the
motion of an interface have been proposed. As with many numerical solutions
to PDEs, it is first necessary to choose between Eulerian and Lagrangian ap-
proaches. In a number of ways, the Lagrangian approach is the more intuitive
for the interface problem, but is subject to many difficulties especially in cases
where the interface distorts or fragments. While it is possible to find solutions
to these issues, the most robust solutions frequently require significant com-
putational resources in addition to presenting complex algorithmic issues. For
these reasons, an Eulerian approach is frequently the better choice.
The simplest Eulerian approach is to use a continuous approximation of
the indicator function. However, this has the interface smearing/diffusion issue
mentioned earlier. For this reason, it is often better to employ an interface track-
ing algorithm, which provide the user with sub-grid resolution of the material
interface.
One Eulerian interface tracking algorithm is the level set method[9], in which
the interface is given by the zero contour of smooth function φ. This generates a
sharp boundary. However for problems such as the ones we are interested in one
must solve a Hamilton-Jacobi equation to advance the location of the interface.
This requires a complex set of methodologies, originally developed for modeling
shock waves such as limiters and entropy conditions. This is still an active area
of research in the finite element community. Additionally, this approach does
not inherently conserve volume, since the level set function φ is not a conserved
quantity.
Another Eulerian approach is the VoF method. Like the level set method,
the VoF method has as a primary goal of maintaining a sharp interface between
the two fluids. In the case of the VoF method this is done using volume fraction
data fi , and then reconstructing the fluid interface on a sub-grid level in each cell
of the mesh where 0 < fi < 1. The interface reconstruction step only requires
access to the volume fraction in neighboring cells, and as such is inherently
parallelizable.

3.3 Equations
The equations solved by Aspect for the velocity field are the incompressible
Stokes equations, given below in equation 1, where f~(~x, t) is the body force,
in this case gravity. In this equation, P is pressure, µ is viscosity, ~u = (u, v)
is the velocity vector, and ρ is density. The temperature is handled using the

2
standard advection-diffusion equation as given in equation 2. Here k is the
thermal conductivity and cp is the specific heat at constant pressure. Any
continuous composition fields C are handled using the advection equation as
given in equation 3. The composition field that we wish to maintain a sharp
interface on using an interface tracking scheme is to be propagated using the
VoF method, which will be described in more detail below.

 2
∂ u ∂2u

∂P
− +µ + 2 =ρ~ex · f~(~x, t) (1a)
∂x ∂x2 ∂y
 2
∂2v

∂P ∂ v
− +µ + 2 =ρ~ey · f~(~x, t) (1b)
∂y ∂x2 ∂y
∇ · ~u =0 (1c)

∂T k
+ ~u · ∇T = ∆T (2)
∂t ρc
∂C
+ ~u · ∇C = 0 (3)
∂t
For the velocity solution, the Boussinesq approximation is used to close the
equations.

ρ(T ) = ρ0 (1 − αv (T − Tref )) (4)

This particular approximation combines the incompressible Stokes equations
with a buoyancy force based on a density variation which is linear in tempera-
ture. The relevant velocity field equations are then given by [from 10]

 2
∂ u ∂2u

∂P
− +µ + =0 (5a)
∂x ∂x2 ∂y 2
 2
∂2v

∂P ∂ v
− +µ + 2 =gρ0 αν (T − T0 ) (5b)
∂y ∂x2 ∂y
∇ · ~u =0 (5c)

Like standard advection equation schemes, the VoF method begins with the
integral form of the standard fluid advection equation (6) where Ωi denotes a
grid cell.
Z Z
d
Φ(~x, t)d~x + ~u · ∇Φ(~x, t)d~x = 0 (6)
dt
Ωi Ωi

At this point standard methods assume a continuous approximation f (~x, t)

of Φ allowing a reduction to a pointwise limit and obtain
R a differential equation
from that. By contrast, VoF methods track fi = Φ(~x, t)d~x, and continue
Ωi

3
from the integral form by making use of the additional constraint that Φ is an
indicator function in order to reconstruct an approximation to the true interface.
The reconstruction does require the assumption of some minimum length scale
for features of interest, and for pure VoF methods also depend on surrounding
cell volumes. The resultant approximation to the true interface Φ∗ then does not
include any of the numerical smearing that would occur for standard advection
methods, so the interface remains localized to a single cell.
To obtain the new state values {fin+1 } from n
R {fi }, the divergence theorem is
applied to the integral form 6. Letting fi = f (~x, tn )d~x, this is then integrated
n
Ωi
over time to obtain the equation
n+1 n+1
tZ Z tZ Z
fin+1 = fin + f ∇ · ~ud~xdt − f~u · ~nd~xdt (7)
tn Ωi tn ∂Ωi

This can then be discretized using average values computed from the finite
element solution for ∇·~u and ~u ·~n, while using a method of characteristics based
approach to find the donor region that is then used to compute the volume frac-
tions on the apertures (∂Ωi x[tn , tn+1 ]). For the incompressible Stokes equations
∇ · ~u = 0 and normally we can assume that this term may be dropped. How-
ever, for some volume fraction update algorithms, such as a dimensionally split
approach, the effective velocity field during the update step may not remain
divergence free and therefore the term is retained.
This discretization of equation 7 is equivalent to a Discontinuous Galerkin
method with zeroth order elements and a specialized flux term computation or
equivalently to a standard finite volume method. As such, the method can be
included in the Finite Element code ASPECT without requiring significant logic
to translate between assumptions.

4 Computational Approach
4.1 Implementation
The necessary finite element matrix assembler and solver for the Boussinesq ap-
proximation have already been implemented as part of ASPECT. It is therefore
sufficient to consider only the implementation of the volume of fluid method,
and the means by which any information necessary to the main algorithm may
be provided to the velocity computation.
To implement the VoF method, it is necessary to select both a a suitable
interface reconstruction method and an algorithm for advecting the volume
fractions based on the reconstructed interface and the computed velocity field.
For the interface reconstruction algorithm, we use the Efficient Least-Squares
VoF algorithm (ELVIRA) developed by Pilliod [5] and Pilliod and Puckett [6]
due to ease of implementation, low computational cost, and second-order accu-
racy. The algorithm is a modification of the full iterative least-squares interface

4
reconstruction method [7] that selects between a short list of candidates shown
to include the true interface for a linear interface, or a second order approxi-
mation for general interfaces given the eight surrounding cells and a maximum
interface curvature. To ensure that linear interfaces are reproduced exactly for
cells on the edge of the mesh, the normal direction of the interface in the cell
at the past timestep is appended to the list of candidates. In the cases where
this addition has any effect, it will result in the selected normal direction being
closer to that generated by the more accurate least-squares algorithm specified
by Puckett [7].
The piecewise linear interface algorithm ELVIRA constructs interfaces of the
form ~n · ~x = d for each cell, with d chosen such that the volume enclosed in the
cell is equivalent to the current volume fraction. For application to the finite
element mesh generated for ASPECT, this was modified to be ~n · x̂ = d where
x̂ is the coordinates for a standard unit cell.
This algorithm requires the computation of the volume fraction or inter-
face location given the other quantity and a particular normal vector rather
frequently. Due to the chosen method of computation, which makes use of the
same relation between the two for the two-dimensional case noted by Scardovelli
and Zaleski [8], the current implementation is limited to a parallelogram mesh.
Acceptance of this implmentation greatly reduces the required complexity of the
code for the moment and can also be expected to significantly reduce the com-
putational cost. Generalizing to permit more freedom in mappings is likely to
necessitate an iterative scheme for computing the appropriate d given a volume
fraction f , which can be expected to significantly increase the computational
cost of the algorithm. See section 6 ‘Future Work’.
In order to simplify the flux fraction calculation, one standard approach for
VoF methods is to use a dimensionally split algorithm with appropriate desig-
nations for the donor regions. This does require some additional consideration
of the method of handling the divergences in the split velocity field as discussed
in section 3.3. The particular method chosen in the case of this implementation
is to include the divergence term for the advection calculation with an implicit
discretization of the volume fraction for said term, ensuring that the volume
fraction remains bounded within [0, 1].
The final implementation problem is how to provide the VoF data to the
main solver. In order to best make use of the existing solver infrastructure,
the volume of fluid data is used to generate a finite element composition field
approximation to the color function that is then used in the standard manner
by the velocity equation solver. The problem is then to select the appropriate
means of making this approximation. Use of a discontinuous galerkin finite
element permits the approximation to be local to cell of interest. Considering
the relative structural scales and configurations, it was decided that polynomial
elements of degree 1 would be most reasonable for both ease of calculating the
approximation and not reaching beyond the level of information available. It
is possible to generate simple interpolations that satisfy any two of the obvious
criteria of bounded on [0, 1], volume fraction preserving, and that the fluid
interface be at C(~x) = 0.5, but attempting to satisfy all three of these is not

5
(a) Linear interface trans- (b) Circular interface (c) Circular interface ro-
lation problem translation problem tation problem, red dot is
center of rotation

Figure 1: Sketches of rate of convergence tests

possible with the given discretization. Considering the case for minimum L2
error for a linear mapping and a P−1 element, the requirements generated are
that both volume fraction and fluid moment match for the true interface and
the approximation. This suggests that the correct approach is to use the first
two criteria.

5 Results
5.1 Validation
In order to ensure that the implementation of the interface tracking method is
second-order accurate for smooth interfaces, we tested it on several problems
advecting a known interface under a prescribed velocity field.
These tests require a method by which to approximate the error of the VoF
method. The obvious approach is to consider the volume where the approximate
and exact interfaces differ. However, even for linear interfaces the calculation of
this quantity is not simple, and useful test cases cannot all have linear interfaces
everywhere.
As such, an approximation by quadrature is necessary, and the use of a
continuous approximation to the Heaviside function is likely to produce more
useful information. In order to select a reasonable approximation, we begin by
noting that the use of an iterated midpoint quadrature could equivalently be
describe as evaluating the differences in volume fractions on given subregions of
the cell with the volume fractions specified only by the type of fluid at the cell
center. When considered in this manner, an obvious improvement is to do the
volume fraction estimation by using linear approximations to the interfaces on
the given volumes, which results in a C 1 function (C 2 almost √everywhere) which
differs from the Heaviside function only in regions within h̃ 2 2 of the interface
where h̃ is the width of a subregion. Examination of the various cases shows that
on volumes with non-disjoint fluid volumes and linear non-intersecting interfaces

6
this is exact. This approximation will be reported in tables as the ‘L1 Interface
Error’.
Another approach to error estimation is to make use of the initialization
procedure to generate the ‘correct’ volume fractions for the discretization at
a given time, and then do a standard `1 norm of the difference between the
the correct and computed volume fractions. This approach does not take into
consideration the current interface reconstruction, since it only measures the
difference in the current volume fractions. This error measurement will be
referred to as the ‘L1 Field Error’.
Both error measurements are normalized by the correct interface length.
This value can then be intuitively considered as the average distance from the
correct interface.

k L1 Interface Error L1 Field Error

4 1.27552 · 10−16 1.32630 · 10−16
5 2.75932 · 10−16 2.74613 · 10−16
6 4.69302 · 10−16 4.69034 · 10−16
7 1.19555 · 10−15 1.19552 · 10−15
8 2.02016 · 10−15 2.01992 · 10−15

Table 1: Error after translation in [0, 1]2 of the linear interface y = 1 − x in the
constant velocity field ~u = (− 14 , − 41 ) for 1 time unit, sketch is in figure 1a

First, the implementation was confirmed to reproduce linear interfaces under

constant velocity fields to within machine epsilon (figure 1a). Due to the lack
of sufficient bordering cells near the boundaries, an exact reproduction did, as
expected, require the previously mentioned addition of an extrapolation of the
past interface normal to the reconstruction candidate normals. The particular
error values for the test are given in table 1, and are all on the order of ma-
chine epsilon. Though the reported error is growing, inspection suggests that
it is likely due to accumulation of machine precision subtraction errors over the
entirety of the mesh, which is likewise increasing in size.
The second, and more complex validation problem was translation of a cir-
cular interface under a constant mesh-aligned velocity field (figure 1b). In par-
ticular, using the domain [−2, 2]2 , a circular interface of radius 12 centered at
(−1, 0) was advected under the velocity field ~u = (1, 0) for two time units under
varying CFL numbers. The results of this test are given in table 2. For CFL
σ = 1, we observe the expected design convergence rate of 2 where the advection
error can be expected to be minimal.
Finally, the algorithm was tested on advecting a circular interface (figure 1c).
In particular, the test problem advects a circular interface of radius 18 centered
at ( 81 , 0) in the region [− 12 , 21 ]2 under the velocity field ~u = π(−y, x). The
results of this test are given in table 3. Again, the observed convergence rate is
second-order.

7
 CFL σ = 1
k L1 Interface Error Rate L1 Field Error Rate
5 1.89618 · 10−3 1.26228 · 10−3
6 3.82585 · 10−4 2.31 2.36274 · 10−4 2.42
7 9.07219 · 10−5 2.08 3.31707 · 10−5 2.83
8 2.07321 · 10−5 2.13 6.34625 · 10−6 2.39

Table 2: Error after mesh-aligned constant velocity field translation of a circular

interface for two diameters, sketch in figure 1b

k L1 Interface Error Rate L1 Field Error Rate

4 8.08431 · 10−3 8.08431 · 10−3
5 2.35671 · 10−3 1.78 2.33441 · 10−3 1.79
6 5.27489 · 10−4 2.16 5.23932 · 10−4 2.16
7 1.41767 · 10−4 1.90 1.40612 · 10−4 1.90
8 3.52731 · 10−5 2.01 3.50275 · 10−5 2.01

Table 3: Error after 1 full rotation of a circular interface offset from center of
rotation by 1 Radius, sketch in figure 1c

5.2 Application
In this section, we present three computations of problems that are relevant to
those that arise in mantle convection. The first two are a standard benchmarks
in the field, and the third is directly related to an active research project I am
engaged in with Professors Turcotte, Billen, and Kellogg in the departmen of
Earth and Planetary Sciences.

Figure 2: Contour and colormap for initial temperature as generated on a 192x64

mesh

8
 For the first problem the newly implemented VoF method is used to track the
uppermost 10% of the fluid volume under the following initial conditions. The
values used for the relevant constants are ρ0 = 1, cp = 1, Tref = 0.5, k = 10−5 ,
αv = 10−4 , and µ = 10−3 for a Rayleigh constant of Ra = 105 . The domain is
[−π, π]x[0, 1] so d = 1 and L = 2π with free-slip and no flow velocity boundary
conditions. The temperature boundary condition are zero Neumann for the
left and right boundaries, and the top and bottom use constant temperature
boundary conditions of T0 = 0 and T1 = 1 respectively. The temperature was
then initialized with perturbation amplitude A = 0.05 to the function given in
equation 8, which is also shown as a colormap and superimposed contour plot
in figure 2.

(
T0 + T1 (1 − yd ) + A sin( yπ xπ 9d
d ) cos( L ) y > 10 or y <
d
10
T (x, y, 0) = (8)
Tref otherwise

(a) Uniform 96x32 mesh

(b) AMR mesh with maximum resolution 384x128

Figure 3: Fluid state for Ra = 105 convection problem

The results for a uniform 192x64 mesh are shown in figure 3a. The general
structure of the fluid interfaces observed matches that described by Turcotte [10]

9
for simple convection rolls. The result of using a AMR refinement scheme to
obtain a slightly higher resolution is given in figure 3b. The contrast between
the two interfaces is likely due to a better resolution of the initial temperature
in the AMR version adjusting the velocity field and thus the length of the initial
timesteps.

Figure 4: Convection problem showing interface reconstruction algorithm re-

solving sub-mesh scale features as ‘bubbles’

One of the design benefits of the VoF method is the ability of the generated
interface to fragment without need for costly and complex special coding. As can
be seen in figure 4., the reconstruction algorithm has resolved a feature with an
effective scale too small to have a definite form as a series of disconnected fluid
regions rather than a connected stream due to beginning with the assumption of
a connected local region. This is a smaller scale loss of accuracy than would be
expected with any code which did not include interface tracking, though it may
be one possible case for comparison when considering reconstruction algorithms.
η = η0

(a) Uniform 128x128 mesh (b) AMR with maximum resolution of

512x512

Figure 5: Comparison of results for the falling block problem with constant
viscosity showing the mesh and the reconstructed interface (white contour)

The second problem is a standard benchmark in the field [2], the ‘Falling

10
 η = 10η0

(a) Uniform 128x128 mesh (b) AMR with maximum resolution of

512x512

Figure 6: Comparison of results for the falling block problem with constant
viscosity showing the mesh and the reconstructed interface (white contour)

Figure 7: Comparison between VoF (left), Bounded Discontinuous Galerkin [3]

(middle), and FEM Advection method in ASPECT (right)

11
Box’ problem. For this problem, an initially square fluid volume with density
ρ = 3.3 · 103 and viscosity η sinks through a background fluid of density ρ0 =
3.2 · 103 and viscosity η0 = 1021 . The results at the final timestep for η = η0 and
η = 10η0 for uniform meshes are shown in figures 5a and 6a. Visual comparison
shows similar final structures to the tracer particle methods in [2]. Note that
although in the figures on the left there is some rounding of sharp corners as
a result of the underresolution, the equivalent AMR tests with a maximum
resolution equivalent to a 512x512 mesh in figures 5b and 6b do not exhibit
this numerical artifact. Examination reveals a much better match to the tracer
approximations in the [2]. A comparsion to other methods for modeling this
problem available in ASPECT is shown in figure 7.
The final test problem is related to an active area of research. Heterogeneous
compositional mantle models are frequently invoked to explain some observa-
tions obtained from geochemistry and seismology. In particular, two regions in
the Earth’s lower mantle, under the Pacific and Africa, are often interpreted
as being piles of dense material, known as Large Low-Shear-Velocity Provinces
(LLSVPs). In collaboration with Professors Billen, Kellogg, and Turcotte we
are currently engaged in a research project to numerically model the thermo-
chemical convection of a two-layer, density stratified region in which each of the
materials obeys equations 1, 2 and 4. The particular configuration was initial-
ized with the division between the composition at z = d2 , a full depth Rayleigh
number of Ra ≈ 105 , and initial perturbations limited to boundary layers of size
d
10 , where d is the height of the domain. The tests were run for varying values
of the nondimensional Bouyancy parameter B = ρ0∆ρ α∆T , with varying initial
temperature perturbations. After a parameter survey, it was found that for a
particular parameter space, the change from stratified to full depth convection
did not take place until the interface between the two fluids showed some initial
perturbation, so the initial conditions were adjusted to permit the addition of
an initial perturbation to the fluid interface. Due to the scale of the behaviour
of interest and the number of test runs necessary, a relatively coarse mesh of size
96x32 was used. The results of the test run are given in figures 8a and 9a for
two interesting values of B. Figures 8 and 9 contain a comparison between the
experiments of Davaille [1] and our computations for both stratified convection
(figure 9) and full depth convection (figure 8) regimes. Note that despite the
long time and expected high flow rate near the composition boundary, there is
very little variation in the reconstructed interface for the case in the stratified
region. Similarly, as is the design goal, the interface remains quite sharp in the
case of the full depth convection mode. For the stratified region, a comparison
to data from Davaille [1] descibing the temperature as a function of depth is
shown in figure 10.

12
(a) Density colormap with temperature contours(white) and reconstructed inter-
face(black) superimposed for model with B = 0.4

(b) Full depth convection from experiment (B = Rρ = 0.2) by Dav Davaille [from 1,
Fig1b]

Figure 8: Comparison of full depth convection model to experiment by Davaille

[1]

13
(a) Density colormap with temperature contours(white) and reconstructed inter-
face(black) superimposed for model with B = 0.7

(b) Stratified convection from experiment (B = Rρ = 2.4) by Davaille [from 1, Fig1f]

Figure 9: Comparison of stratified convection model to experiment by Davaille

[1]

14
(a) Stratified convection from experiment (b) Stratified convection from VoF model
(B = Rρ = 2.4) by Davaille [1, Fig 1a] (B = 0.7)

Figure 10: Comparsion of upwelling temperature-depth plots to experimental

temperature-depth plots from Davaille [1]. Note: These are preliminary results
and no attempt has been made to match the experimental values of the bouyancy
d2
parameter B, viscosity ratio γ, and depth ratio d1 .

15
6 Future Work
• Implementation of a 3D interface reconstruction, and use of the 3D formu-
las for volume frction and interface location in Scardovelli [8] would allow
the examination of the behavior of 3D systems.
• Improve volume fraction calculation to permit use of non-parallelogram
meshes.

• Examination of other interface reconstruction approaches, such as the mo-

ment of fluid algorithm that allows the tracking of multiple materials.
• Examination of alternate advection algorithms; i.e. for flux fraction cal-
culation to the donor region approach commonly used. This is likely to
require significant additional information.

References
[1] Anne Davaille. “Simultaneous generation of hotspots and superswells by
convection in a heterogeneous planetary mantle”. In: Nature 402.6763
(Dec. 1999), pp. 756–760. doi: 10.1038/45461. url: http://dx.doi.
org/10.1038/45461.
[2] Taras V. Gerya and David A. Yuen. “Characteristics-based marker-in-cell
method with conservative finite-differences schemes for modeling geolog-
ical flows with strongly variable transport properties”. In: Physics of the
Earth and Planetary Interiors 140.4 (2003), pp. 293–318. issn: 0031-9201.
doi: http://dx.doi.org/10.1016/j.pepi.2003.09.006. url: http://
www.sciencedirect.com/science/article/pii/S0031920103001900.
[3] Ying He, Elbridge Gerry Puckett, and Magali I. Billen. “A Discontinuous
Galerkin Method with a Bound Preserving Limiter for the Advection of
non-Diffusive Fields in Solid Earth Geodynamics”. In: Physics of the Earth
and Planetary Interiors (submitted) (2016).
[4] M. Kronbichler, T. Heister, and W. Bangerth. “High Accuracy Mantle
Convection Simulation through Modern Numerical Methods”. In: Geo-
physics Journal International 191 (2012), pp. 12–29.
[5] James Edward Pilliod. “An Analysis of Piecewise Linear Interface Recon-
struction Algorithms for Volume-of-Fluid Methods”. MS Thesis. Gradu-
ate Group in Applied Mathematics, University of California, Davis, Sept.
1992.
[6] James Edward Pilliod and Elbridge Gerry Puckett. “Second-order accu-
rate volume-of-fluid algorithms for tracking material interfaces”. In: Jour-
nal of Computational Physics 199.2 (2004), pp. 465–502. issn: 0021-9991.
doi: http://dx.doi.org/10.1016/j.jcp.2003.12.023. url: http://
www.sciencedirect.com/science/article/pii/S0021999104000920.

16
 [7] Elbridge Gerry Puckett. “A volume-of-fluid interface tracking algorithm
with applications to computing shock wave refraction”. In: Proceedings of
the Fourth International Symposium on Computational Fluid Dynamics.
1991, pp. 933–938.
[8] Ruben Scardovelli and Stephane Zaleski. “Analytical Relations Connect-
ing Linear Interfaces and Volume Fractions in Rectangular Grids”. In:
Journal of Computational Physics 164.1 (2000), pp. 228–237. issn: 0021-
9991. doi: http://dx.doi.org/10.1006/jcph.2000.6567. url: http:
//www.sciencedirect.com/science/article/pii/S0021999100965677.
[9] J. A. Sethian. Level Set Methods and Fast Marching Methods. 2nd. Cam-
bridge University Press, 1999.
[10] Donald Turcotte and Gerald Schubert. Geodynamics. 3rd. Cambridge Uni-
versity Press, 2014.

17
Proposed Exam Syllabus
• Fluid Dynamics
– Derivation of the Incompressible Euler Equations
∗ Conservation of Mass
∗ Conservation of Momentum
– Rayleigh-Bernard Equations/Mantle convection
∗ Incompressible Stokes Equations
∗ Temperature Advection/Diffusion Equation
∗ Boussinesq approximation
– Nondimensionalization of the Stokes Equations
– References:
∗ Alexandre J. Chorin and Jerrold E. Marsden. A Mathematical
Introduction to Fluid Mechanics. Springer New York, 1993. doi:
10.1007/978-1-4612-0883-9. url: http://dx.doi.org/10.
1007/978-1-4612-0883-9
∗ Donald Turcotte and Gerald Schubert. Geodynamics. 3rd. Cam-
bridge University Press, 2014
• Finite Difference Methods
– Stability of Explicit Finite Difference Methods for the 1D Advection
Equation
– Consistency of Explicit Finite Difference Methods for the 1D Advec-
tion Equation
– Lax Equivalence Theorem
– Godunov’s Theorem and implications
– Slope Limiters (Eg Minmax, Van-Leer)
– Flux Corrected Transport
– References:
∗ Randall J. LeVeque. Numerical Methods for Conservation Laws.
Springer Science + Business Media, 1992. doi: 10.1007/978-
3-0348-8629-1. url: http://dx.doi.org/10.1007/978-3-
0348-8629-1
• Finite Element Methods
– Derivation of the Galerkin Method for Q1 on a 1D Uniform Mesh (eg
for Heat Equation)
– References:
∗ Fish and Belytschko, A First Course in Finite Elements

18
• Numerical Methods in Linear Algebra
– Jacobi Iteration
– Gauss-Seidel Iteration
– Multigrid
∗ V-cycle
∗ W-cycle
∗ Full multigrid
– Conjugate Gradient method
– References:
∗ William L. Briggs, Van Emden Henson, and Steve F. McCormick.
A Multigrid Tutorial, Second Edition. Society for Industrial &
Applied Mathematics (SIAM), Jan. 2000. doi: 10 . 1137 / 1 .
9780898719505. url: http : / / dx . doi . org / 10 . 1137 / 1 .
9780898719505
∗ Richard Barrett et al. Templates for the Solution of Linear Sys-
tems: Building Blocks for Iterative Methods. Society for Indus-
trial & Applied Mathematics (SIAM), Jan. 1994. doi: 10.1137/
1.9781611971538. url: http:/ /dx.doi.org /10.1137/1 .
9781611971538

19