SOLUTION OF NONLINEAR STOKES EQUATIONS
DISCRETIZED BY HIGH-ORDER FINITE ELEMENTS
ON NONCONFORMING AND ANISOTROPIC MESHES,
WITH APPLICATION TO ICE SHEET DYNAMICS∗

arXiv:1406.6573v3 [cs.NA] 9 Jul 2015

TOBIN ISAAC†, GEORG STADLER,† ‡ AND OMAR GHATTAS† §
Abstract. Motivated by the need for efficient and accurate simulation of the dynamics of
the polar ice sheets, we design high-order finite element discretizations and scalable solvers for the
solution of nonlinear incompressible Stokes equations. In particular, we focus on power-law, shear
thinning rheologies commonly used in modeling ice dynamics and other geophysical flows. We use
nonconforming hexahedral meshes and the conforming inf-sup stable finite element velocity-pressure
disc
pairings Qk × Qdisc
k−2 or Qk × Pk−1 , where k ≥ 2 is the polynomial order of the velocity space. To
solve the nonlinear equations, we propose a Newton-Krylov method with a block upper triangular
preconditioner for the linearized Stokes systems. The diagonal blocks of this preconditioner are sparse
approximations of the (1,1)-block and of its Schur complement. The (1,1)-block is approximated using
linear finite elements based on the nodes of the high-order discretization, and the application of its
inverse is approximated using algebraic multigrid with an incomplete factorization smoother. This
preconditioner is designed to be efficient on anisotropic meshes, which are necessary to match the
high aspect ratio domains typical for ice sheets. As part of this work, we develop and make available
extensions to two libraries—a hybrid meshing scheme for the p4est parallel adaptive mesh refinement
library, and a modified smoothed aggregation scheme for PETSc—to improve their support for solving
PDEs in high aspect ratio domains. In a comprehensive numerical study, we find that our solver yields
fast convergence that is independent of the element aspect ratio, the occurrence of nonconforming
interfaces, and of the mesh refinement, and that depends only weakly on the polynomial finite element
order. We simulate the ice flow in a realistic description of the Antarctic ice sheet derived from field
data, and study the parallel scalability of our solver for problems with up to 383 million unknowns.
Key words. Viscous incompressible flow, nonlinear Stokes equations, shear-thinning, high-order
finite elements, preconditioning, multigrid, Newton-Krylov method, ice sheet modeling, Antarctic ice
sheet.

1. Introduction. We design high-order finite element discretizations and scalable solvers for incompressible nonlinear Stokes equations describing creeping flows
of power-law rheology fluids. Applications include ice sheet dynamics [31], mantle
convection [53], magma dynamics [44] and other problems involving non-Newtonian
fluids [26]. Among the main challenges for the solution of these problems are the presence of local features that emerge from the nonlinear constitutive relation, the strongly
varying and anisotropic coefficients arising upon linearization, the incompressibility
condition leading to indefinite matrix problems, complex geometry and boundary
conditions, a wide range of length scales that may require highly-adapted meshes
with high aspect ratios, and large problem sizes that necessitate parallel solution on
large supercomputers. Our approach to cope with these challenges uses adaptively
† Institute for Computational Engineering & Sciences, The University of Texas at Austin, Austin,
TX, USA
‡ Courant Institute of Mathematical Sciences, New York University, New York, NY, USA
§ Department of Geological Sciences and Department of Mechanical Engineering, The University
of Texas at Austin, Austin, TX, USA
∗ Support for this work was provided by the U.S. Department of Energy Office of Science (DOESC), Advanced Scientific Computing Research (ASCR), Scientific Discovery through Advanced Computing (SciDAC) program, under award numbers DE-FG02-09ER25914, DE-11018096, and DEFC02-13ER2612, and the U.S. National Science Foundation (NSF) Cyber-enabled Discovery and
Innovation (CDI) program under awards CMS-1028889 and OPP-0941678. Allocations of computing
time on TACC’s Stampede under XSEDE, TG-DPP130002, and on ORNL’s Titan, which is supported by the Office of Science of the DOE under DE-AC05-00OR22725 are gratefully acknowledged.

1

refined nonconforming meshes, high-order inf-sup stable finite elements, and iterative
Newton-Krylov solvers combined with multilevel preconditioning techniques. We focus in particular on the construction of efficient solvers and preconditioners for discrete
systems resulting from high-order discretizations.
High-order finite element methods for partial differential equations (PDEs) are
attractive because, in many situations, the discrete solution converges rapidly to the
true solution as the approximation order k is increased or the characteristic mesh
size h is decreased. However, this increased accuracy per degree of freedom compared to low-order methods does not automatically translate into increased accuracy
per unit of computational work. This is due to the fact that matrices arising from
high-order discretizations are denser and, thus, more expensive to apply and to solve
systems with. The cost of applying a matrix arising from a high-order discretization
can be reduced drastically if the work is shifted from memory operations to floating
point operations. This can be achieved using matrix-free implementations and tensorproduct approximation spaces and element operations on hexahedral finite element
meshes. To precondition matrices arising from high-order discretizations, low-order
preconditioners based on the nodes of the high-order discretization have proven efficient [6, 11, 15, 29, 48]. These preconditioners allow fast construction and the use of
methods established for low-order discretizations.
Our approach to solving the nonlinear Stokes equations is an inexact NewtonKrylov method, with a block preconditioning strategy for the linearized equations,
built from preconditioners for the (1,1)-block and for its Schur complement. We consider a power-law rheology that involves the second invariant of the strain rate tensor,
for which the Newton linearization results in a fourth-order anisotropic tensor viscosity. We pay particular attention to the interplay between discretization and solver:
seemingly minor differences in either the discretization or the low-order preconditioner
can vastly impact the performance of conventional solution methods for both diagonal
blocks of the preconditioner.
Our driving application is the simulation of the dynamics of continental-scale
ice flows, which is a critical component of coupled climate modeling. Predicting the
contribution of ice sheets to sea-level rise is difficult because of the complexity of
accurately modeling ice sheet dynamics for the full polar ice sheets and the large
uncertainties in unobservable parameters governing these dynamics [45, Chapter 10,
Appendix 6]. To address these uncertainties, significant effort has been focused on
the development of inverse methods to infer ice sheet model parameters from observations [46, 50]. These inverse methods require the repeated solution of ice flow
equations for numerous parameter fields that may vary over wide ranges, and many
also require the repeated solution of related adjoint ice flow equations. Hence, inverse methods particularly stress the efficiency and robustness of solvers for nonlinear
Stokes equations.
A particular difficulty in ice sheet simulations is the high aspect ratio of the
computational domains, which is inherited by the discretization, leading to anisotropic
meshes. Discretizations with high-aspect ratio elements (and problems with highly
anisotropic material properties, which have many of the same properties) are known
to be challenging for implicit solvers and preconditioners. The development of robust
solvers for high aspect ratio domains is also important in other Earth science and
climate modeling problems. In ocean flow models, for instance, three-dimensional
implicit PDE models are now being used [36], whereas in the past they were often
replaced by two-dimensional approximations.
2

Related work. Several recent articles develop scalable solvers for Stokes problems
with varying viscosity [8,10,20,22,25,27,43]. These contributions use low-order stable
or stabilized finite elements for the discretization of the Stokes equations and address
nonlinearity mostly via a Picard fixed point approach. Scalable solvers for high-order
discretizations of nonlinear scalar problems and extensions to linear incompressible
flow problems are studied in [6]. Various asymptotics-based approximations of the
Stokes equations are used for ice sheet and glacier modeling, which reduce the indefinite Stokes equations to positive definite elliptic systems. These simplifications
are justified by the large differences between horizontal and vertical components of
the velocity; we refer to [30] for a comparison and discussion of the validity of these
different models. Ice sheet simulations using the full nonlinear Stokes equations can
be found, e.g., in [23, 41, 49].
Contributions. One of the main contributions of this paper is the design of discretizations, solvers and preconditioners that allow the fast and scalable iterative
solution of nonlinear Stokes problems. In particular, we obtain convergence that, for
a large class of realistic problems, is independent of the mesh size, the presence of
nonconforming interfaces in the mesh, and the element aspect ratio, and depends only
weakly on the polynomial order. Another contribution is the extension of low-order
preconditioners for high-order discretized problems to meshes with nonconforming
interfaces and high aspect ratio elements.
In addition to analyzing our solver techniques on workstation-sized model problems, we also demonstrate their performance and scalability on a series of larger
problems requiring a distributed memory parallel implementation, including a simulation of the dynamics of the Antarctic ice sheet. The simulation uses a geometry
and temperature field derived from field data and constitutes what we believe to be
the first highly resolved nonlinear Stokes-based continental scale simulations of the
Antarctic ice sheet.
We have also developed publicly available tools for discretizing and solving PDEs
on high aspect ratio domains, such as (but not limited to) those occurring for ice
sheets. One, presented in section 3.1, is an extension to the p4est parallel adaptive
mesh refinement (AMR) library [9] that allows it to construct favorable meshes for
these domains. The other one, discussed in section 5.1, is a plugin for PETSc’s [4]
generic algebraic multigrid (GAMG) preconditioner to improve its effectiveness for
these types of discretizations.
Limitations. We use algebraic multigrid (AMG) for preconditioning. An alternative would be to use geometric multigrid (GMG), which builds a hierarchy through
geometric coarsening of the mesh. GMG can be tailored to only coarsen in certain
directions while leaving others unchanged (semicoarsening), which can be useful for
anisotropic geometries [56]. We use the easier-to-use AMG, in which we incorporate
a minimal amount of geometric information.
Our simulations of ice sheet dynamics use a fixed temperature field and geometry. Simulations of evolving ice sheets would require coupling of the nonlinear Stokes
equations with a time-dependent advection-diffusion equation for the evolving temperature, and with a kinematic equation for the evolution of the ice sheet surface.
However, the solvers presented in this paper carry over as important components in a
time-stepping procedure for the simulation of time-evolving nonlinear viscous flows.
Overview. This paper is organized as follows. In section 2, we discuss the form
of the nonlinear Stokes equations and boundary conditions that are the focus of this
work, their variational formulation, and their linearization. In section 3, we present
3

stable mixed-space finite element discretizations and a discussion of adaptive mesh
refinement. In section 4, we give an overview of our approach to solving the resulting
discrete system of nonlinear equations and the linearized counterparts. The preconditioner for the (1,1)-block of the linearized Stokes systems that arise at each Newton
iteration, which is a critical component of the linear solver, is presented in detail in
section 5, followed by a discussion on the preconditioner for the Schur complement of
the (1,1)-block in section 6. We test our solver on a model geometry in section 7, and
then on a discretization of the Antarctic ice sheet in section 8, where we also study
the scalability of our method. We conclude with a discussion in section 9.
2. Nonlinear incompressible Stokes equations. After specifying the Stokes
equations with strain rate thinning power-law rheology in section 2.1, we present the
corresponding variational form and argue existence of a unique solution in section 2.2.
The linearization of the nonlinear equations is presented in section 2.3.
2.1. Problem statement. On an open, bounded domain Ω ⊂ R3 we consider
the incompressible Stokes equations
−∇ · σ = b,

x ∈ Ω,

(2.1a)

∇ · u = 0,

x ∈ Ω,

(2.1b)

where u is the flow velocity and b is a body force. The Cauchy stress tensor σ depends
on the strain rate tensor D(u) = 21 (∇u + ∇uT ), its second invariant DII (u) :=
1
2 D(u) : D(u) and, possibly, other physical quantities such as a temperature field.
Here, “:” denotes the Frobenius product
between second-order tensors A = (Aij ) and
P
B = (Bij ) defined by A : B = i,j Aij Bij . In the ice sheet problem, which is our
driving application, the stress tensor is given by Glen’s flow law
σ = −pI + B(T )(DII (u) + ε)

1−n
2n

Du,

(2.2)

where p is the pressure, B(T ) a positive-valued function of temperature T , n ≥ 1 is
the strain rate exponent, and ε > 0 a small regularization parameter that prevents
infinite effective viscosity for n > 1. For n = 1, eq. (2.2) reduces to a linear rheology,
and it describes a strain-rate weakening non-Newtonian fluid for n > 1. A common
value used for modeling the flow of glacial ice is n = 3. To complete the definition of
the boundary value problem eq. (2.1), it remains to specify the boundary conditions.
In ice sheet simulations, different parts of the boundary require different combinations of Dirichlet, Neumann, and Robin-type boundary conditions. This makes ice
flow a good problem for developing methods for other creeping flow problems with
complicated boundary conditions. At the ice-air interface, the homogeneous Neumann condition σn = 0 holds. At the ice-water interface, the normal stress matches
the hydrostatic water pressure, i.e., σn · n = −pw and the tangential components
of the boundary traction vanish. At the base of the ice sheet, complex interactions
occur between ice, water, rock, and till. In cold regions, the ice sticks to the bedrock,
while in temperate regions, water accumulates at the base and the ice can slide but is
subject to some amount of friction. A general way to describe these phenomena is to
use a Dirichlet condition in normal direction to describe melting and freezing at the
base of the ice sheet, combined with a Robin-type sliding law relating the tangential
component of velocity Tk u = (I − n ⊗ n)u to the tangential component of the stress
through a function β(·, ·, . . .), i.e.,
Tk σn = −β(|Tk u|, T, . . . ).
4

(2.3)

Physically realistic descriptions of sliding must include the dependence of β on |Tk u|
[61], which makes eq. (2.3) a nonlinear boundary condition for the flow; in this work,
however, we consider the linear case β(Tk u) = β(x)Tk u, and no basal freezing or
melting. To summarize, we use the following boundary conditions for the base ΓR :
Tk σn + β(x)Tk u = 0

x ∈ ΓR ,

(2.4a)

u·n=0

x ∈ ΓR .

(2.4b)

Polar ice sheets have a characteristic depth of less than 5 kilometers, while they
extend horizontally for thousands of kilometers. Because of this difference between
length scales, modelers often simplify eq. (2.1) using asymptotic expansions that require assumptions about the magnitude of the velocities and stresses in the ice sheet,
for instance, the shallow ice approximation [31] and the hydrostatic approximation [5].
The assumptions justifying these simplifications do not hold for the entire ice sheet,
which has led to approaches that combine simplified models in the interior with Stokes
equations at outlet glaciers [54]. To avoid these complications, we do not use simplified
models and focus on the efficient solution of the Stokes equations (2.1) instead.
2.2. Variational formulation. Here, we define a variational form of eqs. (2.1)
and (2.4) that defines the fields (u, p) as the unique solution in a vector space V × M
to
Z
Z
[µ(u)D(v) : D(u) − p∇ · v − q∇ · u] dx +
βTk v · Tk u ds = f (v),
(2.5)
Ω

ΓR
1−n

for all (v, q) ∈ V × M, where µ(u) = µ(u, T ) = B(T )(DII (u) + ε) 2n , and f is the
sum of the effects of body and boundary forces. We assume that B(T ) ∈ L∞ (Ω) is
uniformly bounded from below, that β ∈ L∞ (ΓR ) is nonnegative and that ∂Ω is Lipschitz continuous. For simplicity, we assume that the Dirichlet boundary conditions
are homogeneous and are incorporated into the space V.
For a similar problem, Jouvet and Rappaz [35] show that a unique solution (u, p)
0
exists in the Dirichlet-conforming subspace of [W1,r (Ω)]3 × Lr (Ω), where r = 1 + 1/n
and r0 = 1 + n. In appendix A we define a pair of spaces V and M, in which eq. (2.5)
is well-posed. This pair is only slightly modified from the pair above to account for
the linear Robin boundary condition.
2.3. Newton linearization. The Newton linearization of eq. (2.5) about a
velocity-pressure pair (u, p) are equations whose solution (ũ, p̃) ∈ V × M satisfies
Z
Z
0
[D(v) : (µ (u)D(ũ)) − p̃∇ · v − q∇ · ũ] dx +
βTk v · Tk ũ ds = −r(u, p, v, q)
Ω

ΓR

(2.6)
for all (v, q) ∈ V × M. Here, µ0 (u) is an anisotropic 4th-order tensor given by


n − 1 D(u) ⊗ D(u)
∂µ(u)
µ0 (u) = µ(u)I +
D(u) ⊗ D(u) = µ(u) I −
, (2.7)
∂DII (u)
2n
DII (u) + ε
and r(·) is the residual of eq. (2.5). Here, “⊗” denotes the outer product between two
second-order tensors. Compared to the Newton linearization eq. (2.6), the commonly
used Picard linearization of eq. (2.5) neglects the anisotropic part of the fourth-order
tensor µ0 (u). Using a finite element discretization of eqs. (2.6) and (2.7) is only
marginally more complex than the Picard linearization, as the action of µ0 (u) on
5

D(ũ), which is all that is required in a finite element method, can be computed using
Frobenius products with trial and test functions. The operator µ0 (u) is also found in
the adjoint equations corresponding to eq. (2.5), which are used in inverse methods
to infer uncertain parameters from observations [46, 50].
3. Discretization. Our goal for discretizing eqs. (2.5) and (2.6) is to obtain
discrete solutions that converge to the continuous solutions rapidly in terms of the
number of unknowns and time-to-solution. Our approach combines locally refined
meshes with high-order finite element approximation spaces. The adaptive meshes we
use are discussed in section 3.1, the finite element approximation spaces for the Stokes
equations are described in section 3.2, and computational aspects of the discretization
are addressed in section 3.3.
3.1. Meshing. We use a hierarchical approach to mesh refinement, which starts
with a coarse mesh of conforming hexahedra. This coarse mesh is expected to roughly
describe the geometry of the domain Ω. The fine mesh used for the finite element
discretization is obtained by hierarchical refinement of this coarse mesh. Refinement
can be used to improve the resolution of the geometry of Ω or the resulting small
hexahedra can simply inherit the geometry. For practical as well as numerical reasons,
we require our refined meshes to obey a 2:1 condition, where neighboring hexahedra
can differ by only one level of refinement, as illustrated in fig. 3.2. We use the p4est
library to manage refinement, to enforce this 2:1 condition and to partition the mesh
between processes in parallel computations [9, 32]. We use an extension of the p4est
library, developed by the first author, for meshing three-dimensional problems in
nearly two-dimensional domains. This extension uses a forest-of-quadtrees to manage
independently-refinable columns of hexahedra. Meshes created with this extension
have three key properties that forest-of-octrees meshes lack.
(1) Elements—and by extension degrees of freedom—are organized into columns
that span the height of the mesh. This organization can be exploited by
solvers: see section 5.1.
(2) Each column of elements is assigned to a single process. Split columns, which
can appear in meshes created by the forest-of-octrees approach, can have
negative consequences for solvers: see section 5.1.
(3) Hexahedra within a column may be locally refined in the vertical direction,
allowing for more flexible refinement than the purely isotropic refinement of
the forest-of-octrees approach. We exploit this fine control over the element
aspect ratio in our meshes of the Antarctic ice sheet: see section 8.
This type of mesh refinement is illustrated in fig. 3.1. This extension to the p4est
library appears as collection of “p6est” data types and functions, so named because
they build on elements of both the “p4est” interface for forests-of-quadtrees and
the “p8est” interface for forests-of-octrees. Documentation of the p6est interface is
available online.1
Meshes for the simulation of ice flow must address the different length scales
inherent in the problem. To accurately capture the vertical variations of the state
variables, a minimum vertical mesh resolution is necessary. Most ice sheet models use
∼10 nodes in each vertical column and have a horizontal resolution of 5 km. Since
the average thickness of the polar ice sheets is ∼2 km, the width-to-height aspect
1 For documentation of p6est, see http://p4est.github.io/api/group__p6est.html, as well as
the example/p6est/test/test all.c example distributed with the library, which demonstrates the
major I/O, mesh refinement, and visualization functions.

6

Fig. 3.1: Illustration of the type of mesh refinement provided by the p6est extension of
the p4est parallel AMR library. The hexahedral elements are organized into columns, which
are always assigned to a single process when the mesh is partitioned. Note that horizontal
refinement is decoupled from vertical refinement.

ratio φ of these discretizations is ∼25.2 The Antarctic ice shelves, moreover, are
typically ∼500 m thick, so the same horizontal resolution results in φ ≈ 100. We seek
discretizations and methods for the Stokes equations that support these aspect ratios.
3.2. Finite element discretization. In this section, we describe the finite element spaces used to discretize the velocity and pressure spaces V and M.
3.2.1. Discrete velocity spaces. Given a mesh T of hexahedra {Ki } with possibly nonconforming interfaces, we define a finite-dimensional subspace of C0 (Ω) using
isoparametric Lagrange finite elements. The nk := (k +1)3 nodes Ξk = {ξ rst }0≤r,s,t≤k
that define our Lagrange finite elements are the tensor-product Gauss-Lobatto nodes
of polynomial degree k on the reference domain K̂ = [−1, 1]3 . These basis functions
span Qk (K̂), the space of functions on K̂ that are the univariate polynomials of degree at most k in each of the coordinate directions. We map K̂ to an element Ki
by ϕi ∈ [Qk (K̂)]3 and use the tensor-product Gauss nodes of order k for numerical
quadrature. We define the finite-dimensional space
VT,k = {v ∈ C0 (Ω) : ∀Ki ∈ T, v ◦ ϕi ∈ Qk (K̂)},
and the velocity space WT,k = [VT,k ]3 ∩ V. For a conforming mesh, the set of element
nodes Ξk naturally defines a set of global nodes XT,k for VT,k by the images of element
nodes. For a nonconforming mesh T, however, not all element nodes correspond to
global nodes, as shown in fig. 3.2. To construct a set of global nodes at a nonconforming interface, we thus ignore some element nodes: if an element is more refined than
its neighbor, the element nodes on the interface with that neighbor do not contribute
to the set of points that define the global nodes. Instead, these nodes are known as
hanging or dependent nodes.
In general, function values at the nodes of element Ki must be interpolated from
the global vector of nodal values by a restriction matrix Ri . If Ki has no hanging
nodes, then Ri is simply a one-to-one association of Ξk to a subset of XT,k ; if Ki has
hanging nodes, then Ri interpolates values as described in fig. 3.2. We use identical
trial and test spaces, so
P a global nodal matrix A is assembled from element nodal
matrices {Ai } by A = i RiT Ai Ri .
3.2.2. Discrete pressure spaces. We use inf-sup stable mixed finite element
spaces to avoid the artificial compressibility that can be introduced by stabilized discretizations of incompressible flow. Additionally, to satisfy element-wise incompress2 In glaciology, one often uses the thickness-to-width aspect ratio  = φ−1 as the relevant limit
in asymptotic expansions is  → 0. In this work, however, we prefer using φ because we consider
the thickness of an ice sheet to be its characteristic length, in which we also measure its horizontal
extent.

7

(a)

(b)

Fig. 3.2: (a) A two-dimensional Q3 mesh with a 2:1 nonconforming interface. (b)
The nodal values along a nonconforming interface. Shown in green are the Gauss-Lobatto
nodes of the smaller element, which do not align with those of the larger element: function
values at these nodes depend on the values at the nodes of the larger element, so they are
not included in the global nodal basis. The matrix Ri that interpolates a function to the
nodes of the smaller element must interpolate cubic polynomials to the hanging nodes. This
polynomial interpolation is dense: the value at each of the hanging nodes is dependent on all
of the independent nodes. For two-dimensional nonconforming interfaces, Ri is defined by
tensor-product polynomial interpolation.

ibility, we favor piecewise discontinuous pressure spaces MT,k . This mass conservation
is particularly important for ice sheet simulations, where the change of the mass of
the ice sheet is an important quantity of interest in climate projections.
The two most common choices for approximation on the reference cube are
Pk−1 (K̂), which is the space of polynomials on K̂ of degree at most (k − 1), and
Qk−2 (K̂). We study two possibilities for inf-sup stable velocity-pressure finite element spaces. The pairing Qk (K̂) × Pk−1 (K̂) has an optimal order of convergence, and
has an inf-sup constant that is independent of k and of the type of hierarchical local
mesh refinement we use [28]. Its inf-sup stability, however, degrades with increasing
aspect-ratio φ [3]. We find (see section 6) that this degradation can be significant for
the element aspect ratios in our meshes. An alternative pairing is Qk (K̂) × Qk−2 (K̂),
which has a suboptimal order of convergence, but its inf-sup stability is uniform with
respect to boundary layer refinement making it appropriate for large values of φ [57].
For this pairing, the inf-sup constant decreases as O(k −1 ); however, for the moderate
values of k used in this work, this dependence is not problematic.
3.3. Computational aspects. For a mesh with Nel elements, the number of
degrees of freedom Ndof in a k-order finite element discretization is O(Nel k 3 ), and the
number of nonzero entries in the matrix for a system of equations defined on that space
is O(Nel k 6 ) = O(Ndof k 3 ). This increasing density means that sparse matrix-vector
products based on globally assembled element matrices are not efficient (in terms of
memory operations) for large values of k. We therefore compute nonlinear residuals and apply linear operators using a matrix-free approach to finite elements, where
only the coefficients and fields that define an operation are stored in memory and the
operation’s application to a specific vector is assembled from all element contributions at the time of application. This approach requires O(Ndof ) memory operations
per matrix-vector product or residual calculation. For high-order elements, it is thus
better suited to modern computer architectures, where the bandwidth for memory
operations is much narrower than the bandwidth for floating-point operations. This
reduction in memory operations comes at the expense of more floating point operations, but the tensor structure of Qk (K̂) allows for all such element computations to
be reduced to repeated applications of one-dimensional compute kernels, which can
8

be highly optimized [6, 16]. The same one-dimensional kernels are used to apply the
restriction operators Ri for hanging nodes.
4. Newton-Krylov method for nonlinear Stokes equations. Our goal is
to design a robust and scalable solver for the nonlinear Stokes equations eqs. (2.1)
and (2.2) with boundary conditions that include eq. (2.4). Ideally, the convergence
should be stable with respect to: (1) the element size and the mesh refinement pattern, (2) parameters in the rheology µ(u), (3) the Robin coefficient field β, (4) the
polynomial order k, and (5) the element aspect ratio φ. Here, we propose an inexact Newton-Krylov solver, outlined in sections 4.1 and 4.2 below. In section 4.3, we
present a test problem that allows us to study the effects of the factors listed above.
Variations of this problem will be used in sections 5 and 6 to analyze and optimize
the convergence of our linear solver, and in section 7 for a nonlinear Stokes problem.
4.1. Newton’s method for nonlinear Stokes equations. Given a velocity
and pressure pair (u, p), we (approximately) solve eq. (2.6) for a search direction
(ũ, p̃). We then conduct a line search in the direction (ũ, p̃) using the weak Wolf conditions [47] to guarantee that the nonlinear residual decreases. Each Newton update
is computed inexactly via a Krylov-space iterative method, but with a tolerance that
decreases in subsequent steps so as to guarantee quadratic convergence of the Newton
iterations [18] close to the solution. We ensure that discretization and differentiation
commute, so that the Jacobian obtained from discretizing eq. (2.6) is equivalent to differentiating the discretization of eq. (2.5). Discretization with one of the stable finite
element pairs discussed in section 3.2 results in a linear system with the symmetric
saddle-point system matrix

A(u) =


F (u) B T
,
B
0

(4.1)

where the (1,1)-block F (u) is the discretization of the sum of the terms involving β
and µ0 (u) in eqs. (2.6) and (2.7), and B is the discretized divergence operator.
4.2. Preconditioned Krylov method for linearized Stokes equations.
We solve systems involving A(u) using preconditioned Krylov space methods, typically restarted GMRES, or, if the preconditioner is not constant, its flexible variant
FGMRES [51]. As is well known, the performance of Krylov methods critically depends on the availability of an efficient preconditioner Ã for A(u). In the following, we
use the notation A = A(u) and F = F (u), i.e., in our notation we neglect the dependence of F and A on u. Due to the elliptic nature of F , a purely local preconditioner
for A cannot provide h-independent convergence and a multilevel preconditioner is
required. There are two main approaches for multilevel preconditioners for incompressible flow problems, namely monolithic and block preconditioning approaches.
The former approximates the saddle point system on each level of a hierarchy and
employs smoothers that are based on approximate local saddle point solutions (i.e.,
Vanka-type smoothers) [14, 34]. This approach typically requires a geometric mesh
hierarchy or involves stabilized discretizations. In contrast, block preconditioners are
built from preconditioners for F and for the Schur complement with respect to the
(1,1)-block, S := −BF −1 B T . They allow to build on existing solvers for elliptic systems and do not impose restrictions on the discretization underlying A. Due to this
flexibility, we follow this latter approach and use an upper-triangular block precondi9

tioner Ã, such that the preconditioned system becomes


−1 
F̃ B T
I + εF̃
F BT
−1
AÃ =
=
B 0
0 S̃
B F̃ −1


εF̃ B T S̃ −1
,
I + εS̃

(4.2)

where εF̃ = I − F F̃ −1 and εS̃ = I − (−B F̃ −1 B T )S̃ −1 . Here, the matrix F̃ is an
approximation of the (1,1)-block F , and S̃ is an approximation to the Schur complement S with respect to the (1,1)-block. If εF̃ and εS̃ vanished, AÃ−1 would be
a lower triangular matrix with all eigenvalues clustered at 1, and a preconditioned
Krylov method would converge in two iterations. Thus, our target is to devise F̃ −1
and S̃ −1 such that εF̃ and εS̃ are small and Ã-preconditioned Krylov methods converge quickly. Additionally, it is important that the setup time for the preconditioner
Ã is small, because Ã is recomputed as the Jacobian F = F (u) changes. In sections 5
and 6, we develop block preconditioners F̃ and S̃ and study their properties.
In sections 7 and 8, when we apply our preconditioner to problems with dimensional coefficients, we combine the preconditioner Ã that we describe in this section with a diagonal rescaling of the linearized Stokes system of the type described
in [43, Section 2.6]. This rescaling tends to reduce the round-off errors that result
from matrix entries that differ by several orders of magnitude due to different units.
4.3. Test problem setup. Our test problems are based on problem C from the
Ice Sheet Model Intercomparison Project [49]: Ω is a cutout of an “infinite slab”,
i.e., a sheet that is periodic in x- and y-directions. We use homogeneous Neumann
boundary conditions on the top surface of Ω; on its base, we employ homogeneous
Dirichlet conditions in the normal direction and Robin-type conditions in the tangential directions, as in eq. (2.4). We use the same boundary conditions when testing
systems in just the (1,1)-block F . In linear test problems, we use µ = 1 and β = 1
for the material and the sliding coefficients, unless specified otherwise.
To study the behavior of our solvers in the presence of nonconforming faces, we
use meshes for Ω that have nonconforming faces throughout, as illustrated in fig. 4.1.
We use two base meshes: Txy , which has nonconforming interfaces in only the x- and
y-directions, and Txyz , which has nonconforming interfaces in all directions. Note
that the latter mesh cannot be decomposed into columns, so it represents a mode of
mesh refinement that appears when using octree-based meshes, but not when using
the hybrid AMR approach discussed in section 3.1. In all test problems, Ω has unit
thickness, but we vary the length and width L of the domain. The finite elements
stretch with the domain, so that the number and pattern of elements is the same, but
the element aspect ratio φ varies between 1 and 100; see fig. 4.1b. These meshes have
576 elements, and the number of degrees of freedom in the Stokes systems discretized
on these meshes is roughly 2300(k − 1)k 2 .
For linear test problems, we test the effectiveness of our methods when the residual
contains multiple length scales. To achieve this, we compute right hand sides from
manufactured solutions, i.e., b = A(u∗ )(u∗ , p∗ )> when testing the linear Stokes solver
and b = F (u∗ )u∗ when testing the (1,1)-block solver. For that purpose, we create
scalar-valued, spatially variable fields s as the sum of a Fourier series with random
coefficients and pointwise random component:


cos ωjx cos ωky
N
X
 cos ωjx sin ωky 
−γ

s(x, y, z) =
(aj,k , bj,k , cj,k , dj,k )T 
 sin ωjx cos ωky  |(j, k)| + e(x, y, z), (4.3)
j,k=0
sin ωjx sin ωky
(j,k)6=(0,0)
10

(a)

(b)

(c)

(d)

Fig. 4.1: Meshes and manufactured solution used for the test problems. Elements are
colored according to their level of refinement. (a) Mesh Txy , which contains nonconforming
faces in x- and y-directions, for φ = 1. (b) Mesh Txy for φ = 10. (c) Mesh Txyz , which
contains nonconforming faces in all directions. (d) A manufactured solution u∗ , constructed
to contain variations with several length scales.

where N = 10, ω = 2π/L, γ = 3/2, and i = 1, 2, 3. The coefficients aj,k , bj,k , cj,k and
dj,k are randomly chosen from [−1, 1], but the Fourier coefficients decay because of
the |(j, k)|−γ := (j 2 + k 2 )−γ/2 term, making γ a control of the smoothness of s. The
extra term e is a random value from [−1/4, 1/4], added at each node of the discrete
vector. To generate the vector field u∗ that is the manufactured solution of our test
problems for the (1,1)-block solver, we generate a field s for each component of u∗ .
The magnitude of such a velocity field is shown in fig. 4.1d. In all our tests, the norm
we use to report the convergence of iterative solvers is the `2 -norm of the residual.
5. Preconditioning the (1,1)-block. The (1,1)-block F occurring in the Stokes
system is similar to the operator arising in linear elasticity. If we neglect boundary
conditions, its nullspace is given by the rigid-body modes. To approximately apply
the inverse of F , we use algebraic multigrid (AMG) and, in particular, we use the
smoothed aggregation multigrid (SA-AMG), which has theoretically proven convergence bounds [59]. SA-AMG uses Galerkin projections to create coarse approximation
spaces, where the coarse space is embedded in the fine space by a prolongation matrix
P . P is constructed by first creating a projector P̃ that projects coarse “aggregate”
nodes onto disjoint sets of fine nodes, followed by creating P from P̃ by applying a local smoothing operation based on F , while ensuring that the near nullspace of F (the
nullspace in the absence of boundary conditions) is well-approximated in the coarse
space [60]. To construct the prolongation and coarse matrices of the hierarchy, we use
SA-AMG as implemented by PETSc’s GAMG preconditioner. To build smoothers
for each level, we use PETSc’s KSP and PC objects for defining Krylov methods and
preconditioners. In this framework, we have three main design parameters, namely:
(1) The matrix F̃ used to construct the multigrid hierarchy. As multigrid is only
used as preconditioner and not as solver, F̃ can be based on a lower-order
11

(and thus sparser) approximations of the high-order element discretization of
the (1,1)-block used in the residual computation.
(2) How the hierarchy is coarsened. If fine nodes are only grouped into aggregates
when they are strongly connected to each other, then there will be more
aggregates and a less aggressively coarsened hierarchy. In standard SA-AMG,
connections between fine nodes are ignored when forming aggregates if they
are weaker than some threshold θ.
(3) The choice of smoothing operation on each level of the hierarchy.
These three choices cannot be made independently. The aggregation strategy and
the smoothers are related because the error components not damped by the smoother
must be corrected on the coarser levels. Moreover, the choice of F̃ affects the hierarchy
and the effectiveness of different smoothers. In this section, we study these parameters
on linear systems involving the (1,1)-block F rather than the full Stokes system A.
We note that two additional parameters, the size of the coarsest grid and the solver
used for the corresponding system, can also impact the scalability and performance
of multigrid preconditioners. In this work, however, we have used the default coarse
grid size of GAMG (which is less than 50 unknowns), and we use a direct solver for
the corresponding system in all of our numerical experiments. Note that in parallel
simulations, we repartition the coarsened multigrid levels to subsets of processors,
and, in particular, the system for the coarsest grid only uses a single processor.
5.1. SA-AMG aggregation and smoothers. The convergence of multigrid
using pointwise smoothers such as Jacobi and Gauss-Seidel is known to degrade for
anisotropic problems when isotropic coarsening is used [58, Chapter 4]. This slowdown
can also be seen in fig. 5.1a, where, for different element aspect ratios φ, we show
the convergence based on a symmetric Gauss-Seidel (e.g., symmetric successive overrelaxation, SSOR) smoother for hierarchies created with aggregation threshold θ =
0. This behavior occurs because for large φ, horizontal variations in the error have
less energy in the operator norm than vertical variations, and are thus not reduced
sufficiently. There are two primary approaches to addressing this issue. One can use
pointwise smoothers and use smaller aggregates to construct an operator hierarchy
such that the undamped error components are well represented on the coarser meshes
and can be corrected there (semicoarsening), or one can continue to use full-sized
aggregates and use a smoother that more effectively dampens all error components
(non-pointwise smoothing).
Although some form of semicoarsening could be effective for our problems, we
have not had success with it. In standard SA-AMG, semicoarsening is accomplished
by only aggregating strongly connected nodes, for which the corresponding matrix
entries satisfy
|Fij |2 ≥ θ|Fii ||Fjj |.

(5.1)

While this heuristic is generally good for scalar elliptic problems, for the vector system
F it is problematic. It is known that for a fixed threshold θ > 0, a sufficiently large
aspect ratio φ will cause the standard SA-AMG algorithm to break down for operators
whose nullspace contains the rigid-body modes [24, Section 5.3]. Even before this
breakdown, we observed poor convergence rates with θ > 0 and pointwise smoothers.
A better heuristic for anisotropic problems than eq. (5.1), which is reported to result
in good convergence for anisotropic elasticity problems, can be found in [37, Section
4], but even in this case, semicoarsening produces larger coarse operators and deeper
hierarchies, which increase the cost of storing and applying the preconditioner. We
12

(a) krj k/kr0 k: k = 3,
SSOR smoothing
100

(b) krj k/kr0 k: k = 3
IC(0) smoothing
100

φ=1
φ = 10
φ = 100

10−8

10−8

10−16

φ=1
φ = 10
φ = 100

0

20
40
Krylov iteration j

10−16

60

0

20
40
Krylov iteration j

60

Fig. 5.1: Convergence histories for the (1,1)-block test problem described in section 4.3,
discretized with k = 3 on mesh Txy (fig. 4.1a) at different aspect ratios. The solver is GMRES preconditioned by one SA-AMG V-cycle. The solver in (a) uses a Chebyshev(1)/blockJacobi/SSOR smoother, i.e., one optimally-damped application of block-Jacobi/SSOR; the
solver in (b) is the same, except that it uses an incomplete Cholesky factorization (IC(0))
instead of SSOR.

therefore choose not to use thresholding when constructing our hierarchies, i.e., we
do not consider the magnitudes of matrix entries when constructing our aggregates,
only the nonzero pattern of the matrix. This choice results in aggressively coarsened
hierarchies with smaller operator complexities, and has the additional advantage that
it allows us to reuse the projection matrices for multiple Newton steps, reducing the
cost of the subsequent preconditioner setups.
Since we must select a smoother that is compatible with our aggressive coarsening,
we use approaches based on incomplete factorizations of the matrix F̃ . In parallel,
we only compute the incomplete factorization for each process’s diagonal block of
the distributed F̃ matrix, which amounts to a block-Jacobi smoother. To build a
stationary smoother from this block-Jacobi/IC(0) smoother, an estimate of the largest
eigenvalue of the smoothed operator is required, either to calculate a single damping
parameter or the coefficients of a Chebyshev polynomial. We numerically estimate
the largest eigenvalue using an iterative method, which adds to the setup cost of
the multigrid preconditioner. To avoid this setup cost, one can alternatively use a
non-stationary smoother that does not require damping, such as a few iterations of a
Krylov method. This is a good choice when only a small reduction in the residual is
required, as for instance in the early iterations of an inexact Newton-Krylov method.
We note that in certain parameter regimes—large values of β or no-slip Dirichlet
conditions, and a fixed horizontal resolution—incomplete block-Jacobi factorizations
have been observed to work well as stand-alone preconditioners, even with fine vertical
resolution [7, Section 3.1]. However, we are not aware of a systematic study of the
parameter regimes where this good performance can be observed.
Using local Fourier analysis, one can show that geometric multigrid smoothed
by incomplete factorization can give φ-independent convergence for scalar elliptic
problems [58, Chapter 7]. A key component of this analysis is the ordering of degrees
of freedom in the incomplete factorization: tightly-coupled degrees of freedom within a
column must be sequential. We can order the degrees of freedom in our original system
this way, but the standard techniques for generating aggregates, such as GAMG’s
13

`
0
1
2

(a) standard
β = 10 1−2
0.14
0.14
0.17
0.27
0.20
0.51

MIS
1−4
0.57
0.75
0.82

1−8
0.63
0.78
0.83

`
0
1
2

(b) column-preserving MIS
β = 10 1−2 1−4 1−8
0.14
0.14 0.30 0.35
0.17
0.20 0.39 0.47
0.19
0.25 0.48 0.54

Fig. 5.2: The convergence factor for SA-AMG applied to the (1,1)-block test problem
for increasing mesh refinement and decreasing boundary condition strength. The aspect ratio
of the elements is φ = 100; the polynomial elements order is k = 3; the meshes are defined by
` levels of isotropic refinement of the mesh Txy . We show the convergence factor (estimated
over the number of iterations needed to reduce the residual by a factor of 10−14 or 100
iterations) for decreasing magnitudes of the basal Robin coefficient β. Table (a) shows the
results for a hierarchy constructed using standard randomized MIS aggregate construction,
and table (b) shows the same for a hierarchy constructed using the column-preserving MIS
aggregate construction described in the text.

default method of randomized maximal independent set (MIS) selection [2], result in
coarse discretizations having no special spatial structure.
Because of this, the multigrid hierarchies created by randomized MIS aggregation
are efficient for solving the (1,1)-block equations at a high aspect ratio only if the
boundary conditions are strong enough, which we demonstrate in fig. 5.2a. If the
basal boundary has full Dirichlet conditions, or a strong Robin coefficient β (β ≥ 1
in our non-dimensional model problem), we see that GMRES preconditioned by such
a V-cycle converges with near h-independence. If, however, we have β ≈ 0, then the
V-cycles become less efficient and lose h-independence. The culprits for this loss of
efficiency are error modes which are highly oscillatory in the x- and y-directions but
nearly constant in the z-direction. At high aspect ratios, the stress energy of these
modes is low in the absence of a strong boundary condition, so the smoother does
little to dampen them. If the coarse mesh correction is poor, then it will introduce
error components in these modes.
Our solution to this problem of weak boundary conditions is to use meshes whose
degrees of freedom form columns such as the hybrid AMR meshes discussed in section 3.1, and to modify the aggregation technique used by GAMG in order to preserve the column-structure of these degrees of freedom. We first construct aggregates
of columns using a standard aggregation technique such as randomized MIS, and
then partition each aggregated column into the final node aggregates. When subdividing a column, we try to ensure that each node aggregate is at least three nodes
tall, because three is usually the minimum diameter of aggregates used in standard
SA-AMG. The multigrid hierarchies created by this method are much closer to achieving h-independent convergence for weak boundary conditions, as we demonstrate in
fig. 5.2b. We note that this method appears to be effective even when the number
of degrees of freedom in each column varies as they do in our test meshes. We also
note that this technique does not depend on our particular choice of hexahedral elements, but could also be applied, e.g., to meshes with triangular prism elements,
which are commonly used in modeling ice sheets and other high aspect ratio domains.
We implement this aggregation technique in the DofColumns plugin3 for GAMG.

3 The

DofColumns plugin is available from https://bitbucket.org/tisaac/dofcolumns/.
14

5.2. Construction of low-order approximation F̃ . AMG requires assembled
matrices for the construction of the coarse grid hierarchy. For matrices arising from
high-order discretizations, this assembly requires significant memory and computation
compared to low-order discretizations [6,29]. In three dimensions, the cost to construct
an element matrix Fi is O(k 7 ) per element, or in terms of the number of nodes per
7/3
element O(nk ). Instead of the true element matrix Fi for element Ki , we construct
an approximation F̃i that treats the variables associated with nodes of the Qk (K̂)
finite element as variables for a k × k × k grid of Q1 (K̂) finite elements with corners
at the high-order node locations. This matrix F̃i can be constructed in O(nk ) steps.
Aspects of this lower-order preconditioning technique have been studied for simple
problems: spectral equivalence between high-order and lower-order discretizations of
the Laplacian has been proven in two dimensions [39] and is demonstrated numerically
in three dimensions. The (1,1)-block F differs from the operators used in previous
studies in that it involves variable coefficients, high-aspect ratio elements, and nonconforming interfaces between elements. As we will show below, the interaction of
these factors affects the stability and effectiveness of low-order preconditioning.
5.2.1. Influence of quadrature on low-order preconditioning. The effectiveness of the AMG preconditioner for F also depends on the choice of the quadrature
used to construct F̃ . Let us denote by F̃G and F̃GL the low-order matrices computed
with Gauss-Legendre and Gauss-Lobatto quadrature, respectively. We find that F̃GL
leads to a better and more robust preconditioner, particularly when combined with
an incomplete factorization smoother. The greater stability of Gauss-Lobatto quadrature comes from its diagonal-lumping behavior [17]. In tensor-product elements, the
directional derivative operators that are used to construct the element matrices can
be decomposed as Kronecker products of the form B ×B ×D, where B maps 1D nodal
values to 1D values at quadrature points and D maps 1D nodal values to derivatives
at quadrature points. When we refer to Gauss-Lobatto quadrature as “diagonallumping”, we mean that for Gauss-Legendre quadrature, B has off-diagonal entries,
whereas for Gauss-Lobatto quadrature B is the identity. This does not make F̃GL diagonal, but does increase sparsity: F̃GL contains ∼30% fewer nonzero entries than F̃G .
The factors of an incomplete factorization can become ill-conditioned when the original matrix is far from being diagonally dominant, particularly if there is no pivoting
performed during the factorization [52, Chapter 10.5]. The diagonal-lumping inherent
in Gauss-Lobatto quadrature increases the magnitude of diagonal entries relative to
off-diagonal entries. An additional advantage of the sparsity of F̃GL relative to F̃G is
the reduced cost of the incomplete factorization and the hierarchy construction.
5.2.2. Influence of nonconforming meshes on low-order preconditioning. As discussed in section 3.2, C 0 -conforming discretizations on nonconforming
meshes require element restriction matrices Ri that contain dense blocks for nonconforming interfaces. For an element with hanging nodes, the product RiT F̃i Ri can
5/3
at best be computed in O(k 5 ) time, or in terms of nk , O(nk ). One possibility for
recovering O(nk ) element assembly is to replace Ri with a matrix R̃i , in which the
values of hanging nodes only depend on the nearest independent nodes, as illustrated
in fig. 5.3. The convergence plots in figs. 5.1 and 5.2 were generated on the mesh Txy ,
which is the default mesh for our test problems, using these sparse-approximation restriction matrices {R̃i }: in this case, their use does not affect the convergence. If we
use the mesh Txyz (fig. 4.1c), which is generated by standard isotropic octree-based
refinement, the errors incurred by these sparse restriction operators increase as φ in15

(a)

(b)

Fig. 5.3: Illustration of different constructions of the low-order matrix F̃ at nonconforming interfaces. F̃ can be thought of as the matrix for a finite element space whose functions
are piecewise linear between the nodes of the original high-order finite element space. Because these nodes do not align at nonconforming interfaces (see fig. 3.2), the functions in this
low-order space have discontinuities. The nature of these discontinuities depends on how the
nodal values of the smaller element, shown in green, are interpolated from the global nodes,
which coincide with the nodes of the larger element, shown in red. (a) If the same restriction
operator Ri is used as for the original high-order finite element space, dependent hanging
nodes are interpolated as though the independent nodes represent high-order polynomials.
Each hanging node depends on every independent node along the nonconforming edge/face.
(b) If the approximate restriction operator R̃i is used, dependent hanging nodes are interpolated linearly (or bilinearly for two-dimensional faces) from the closest independent nodes.
This results in sparser matrices and faster global matrix construction.

krj k/kr0 k: k = 3, Txyz
100

Ri ,
Ri ,
Ri ,
R̃i ,
R̃i ,
R̃i ,

10−8

10−16

0

20
40
Krylov iteration j

φ=1
φ = 10
φ = 100
φ=1
φ = 10
φ = 100

60

Fig. 5.4: Relative residual krj k/kr0 k versus number of Krylov iterations for the test
problems described in section 4.3 for the mesh Txyz (fig. 4.1c), which has nonconforming
faces in all directions.We compare two possibilities of handling nonconforming interfaces in
the construction of the low-order preconditioner F̃ at different element aspect ratios. The
first possibility is to use {Ri }, the same hanging node restriction matrices as in the high-order
discretization. An alternative is to use {R̃i }, a piecewise linear approximation of {Ri }.

creases. Using the true {Ri } for Txyz improves the convergence, but we still observe
slight φ-dependence (see fig. 5.4).
To understand this behavior, note that in essence we are comparing two approximations to the high-order element matrix RiT Fi Ri : one, RiT F̃i Ri , where the element
matrix is replaced by a low-order approximation, and another, R̃iT F̃i R̃i , where additionally the high-order interpolation of the nodal values at nonconforming faces is
replaced by a low-order interpolation. The latter approximation only affects nonconforming faces; the fact that nonconforming interfaces are much larger in Txyz than
16

(a) F uj = λj F̃ uj , Txy
10

(b) F uj = λj F̃ uj , Txyz

2

Ri
R̃i

λj

101
100

10−1
10−2

0

1,000

2,000

0

j

1,000

2,000

j

Fig. 5.5: Comparison of the generalized eigenvalues λ that satisfy F u = λF̃ u for two
meshes, Txy and Txyz , and for F̃ constructed with either the true restriction operator Ri or
its sparse approximation R̃i . (a) The generalized eigenvalues for the mesh Txy are almost
the same whether Ri or R̃i is used. (b) The generalized eigenvalues for the mesh Txyz are
sensitive to the use of R̃i .

in Txy for large φ seems to explain the different convergence behavior for problems
discretized on the two meshes.
To investigate the influence of nonconforming faces on the low-order preconditioner numerically, we consider the generalized eigenvalue equation F u = λF̃ u for
φ = 100 on the meshes Txy and Txyz ; the results are shown in fig. 5.5. For eigenvectors u with eigenvalues far from 1, F̃ is a poor approximation of F . For the mesh Txy ,
the errors incurred by the sparse restriction matrices {R̃i } due to the nonconforming interfaces are small. For Txyz , errors incurred by {R̃i } are more significant, and
inspection of the extremal eigenvalues shows that they are associated with vectors u
for which F u is large only on nonconforming interfaces that are normal to the z-axis.
5.2.3. Residual computation for smoother on the finest mesh. When the
matrix F̃ used to generate a multigrid hierarchy differs from the true matrix F , one
has two possibilities for the residual computation in the smoothing step on the finest
mesh. Namely, one can use the high-order discretized operator, i.e., r = b − F x, or its
low-order approximation, i.e., r = b − F̃ x. The convergence in fig. 5.6a, which used F̃
to define the AMG residuals, should be compared to fig. 5.1b, which uses F to define
residuals and converges faster. Thus, in the following, we use the high-order operator
for the residual computation in the smoother on the finest mesh.
5.3. Convergence for different orders k. In fig. 5.6b we show the convergence
of our solver (GMRES, preconditioned by column-preserving SA-AMG with a damped
block-Jacobi/IC(0) smoother) for polynomial orders k=3, 6, and 9 on a mesh with
φ = 100. As can be seen, the iteration number is independent of the polynomial
order. To illustrate that the preconditioner is also k-independent with respect to
computational work, note that the operator complexities of the AMG hierarchies
with respect to F̃ (i.e., the sum of nonzero matrix entries in all operators in the AMG
hierarchy divided by the number of nonzero entries in F̃ ) reported by PETSc are
1.22, 1.26, and 1.27 for k = 3, 6 and 9, respectively, and that the average numbers
of nonzeros per row in F̃ are 86.83, 82.33, and 82.36. This demonstrates that the
cost to construct F̃ and the coarse hierarchy are proportional to the problem size,
17

(a) krj k/kr0 k: k = 3,
F̃ -based AMG residual
101

(b) krj k/kr0 k: φ = 100
100

φ=1
φ = 10
φ = 100

10−7

10−8

10−16

10−15
0

20
40
Krylov iteration j

k=3
k=6
k=9

60

0

20
40
Krylov iteration j

60

Fig. 5.6: Relative residual krj k/kr0 k versus number of Krylov iterations for the
test problems described in section 4.3 on mesh Txy using SA-AMG smoothed by blockJacobi/IC(0), (a) when the residual in the smoother on the finest level of the AMG hierarchy
is computed with F̃ instead of F , and (b) for varying orders of approximation.

but independent of k. The same can be said of the cost of computing the incomplete
factorization smoother on each level as we do not allow fill-in. Thus, the overall
computational complexity for preconditioning F is independent of k.
6. Preconditioning the Schur complement of the (1,1)-block. In this
section, we analyze how different approximations of the Schur complement S of the
(1,1)-block and different choices of the basis for the pressure space M affect the convergence rate for Stokes linear problems preconditioned as in eq. (4.2). For all our
tests, the preconditioner for the (1,1)-block F is a single multigrid V-cycle with the
parameter choices described above in section 5. We test the effectiveness of our preconditioner on Stokes problems whose setup is discussed in section 4.3. As discussed
in section 3.2, the discrete inf-sup constant for the Qk × Pk−1 mixed element is φdependent. In fig. 6.1a, we demonstrate this φ-dependence numerically. Because of
this instability we prefer the Qk × Qk−2 mixed element for anisotropic problems.
For constant viscosity µ > 0, S is known to be spectrally equivalent to the scaled
pressure mass matrix −µ−1 M . Because of this equivalence, a common choice is to
approximate S by −M̃ (µ−1 ), which is a diagonally-lumped approximation to the
µ−1 -weighted mass matrix. This Schur complement approximation does not take
into account the anisotropic part of the 4th-order tensor µ0 (u) defined in (2.7). As an
alternative to the weighted mass matrix, the least-squares commutator, also known as
the BFBt preconditioner, has proven to be a good Schur complement approximation,
in particular in the presence of strongly varying coefficients [19, 43].
We have implemented both of these Schur complement preconditioners, and have
found −M̃ (|µ|−1 ) to be the most efficient for the problems targeted in this paper. One
reason for this good performance is that we choose a basis for the discontinuous finite
element pressure space that nearly diagonalizes the mass matrix M for the pressure
space. As a consequence, the effect of replacing −M (|µ|−1 ) with its mass-lumped
counterpart −M̃ (|µ|−1 ) is minimal. Such a basis for Qk−2 (K̂) is given by a Lagrange
basis for the tensor-product Gauss nodes. For a mapped element Ki , the mass matrix
remains nearly diagonal provided the mapping is moderately nonaffine. To illustrate
the effect of the choice of the basis, in fig. 6.1b we compare the convergence when
18

(a) krj k/kr0 k: k = 3

(b) krj k/kr0 k: k = 3

100

100

10−7

10−7

10−14

0

20
40
Krylov iteration j

Qk−2 , φ = 10
Pk−1 , φ = 10

10−14

60

Gauss
Gauss-Lobatto

0

20
40
Krylov iteration j

60

Qk−2 , φ = 100
Pk−1 , φ = 100

Fig. 6.1: Relative residual versus number of Krylov iterations for a Stokes test problem,
posed on the same domain as the test problems in section 4.3, but with a uniformly refined
mesh. The convergence obtained with different discretizations or bases for the pressure space
are compared. (a) Comparison between the pressure spaces Qk−2 and Pk−1 . (b) Comparison
between Gauss points and Gauss-Lobatto points for nodal bases of the pressure space Qk−2 .

(a) krj k/kr0 k: k = 3
100

(b) krj k/kr0 k: φ = 100
100

φ=1
φ = 10
φ = 100

10−8

10−16

k=3
k=6
k=9

10−8

0

20
40
Krylov iteration j

10−16

60

0

20
40
Krylov iteration j

60

Fig. 6.2: Relative residual versus number of Krylov iterations for Stokes test problems.
(a) The convergence for different aspect ratios φ is compared. (b) The convergence for different polynomial orders k is compared.

using this Lagrange basis for the Gauss nodes to the convergence when using the
Lagrange basis for Gauss-Lobatto nodes, which is more commonly used as a basis for
tensor-polynomial finite elements.
In fig. 6.2a, we observe that the preconditioner for the Stokes system A results
in convergence independent of the aspect ratio φ; however, as shown in fig. 6.2b, we
find a dependence of the convergence on k. Note that this differs from our findings
obtained for the (1,1)-block shown in fig. 5.6b, which shows independence of k for
solving systems with F . The O(k −1 ) decay of the lower bound for the discrete inf-sup
constant for Qk × Qk−2 suggests that this k-dependence cannot be avoided.
7. Nonlinear ice stream problems with smooth and rough beds. To
test the nonlinear solver for eq. (2.5), we adapt a model problem from [13]. As in
19

(a)

(b)

Fig. 7.1: Ice stream model problem on a 400×400 km periodic domain of 1 km thick
ice with a 0.5◦ slope to the right. Shown in (a) is the Robin coefficient field β beneath the
ice and in (b) the magnitude of the velocity u at the top surface. A fast flowing ice stream
develops due to the shear thinning rheology.

section 4.3, the domain is a cutout of an infinite slab that is periodic in the horizontal
dimensions, but the pitch of the domain relative to the direction of gravity is θ = 0.5◦ ,
so that a flow is induced. The Robin coefficient field β is shown in fig. 7.1a. Although
β varies smoothly, the nonlinearity of the rheology causes the velocity u to develop a
narrow region of fast flow similar to an ice stream, as shown in fig. 7.1b. The constants
in the constitutive relationship eq. (2.2) are n = 3 and B(T ) ≡ 2.15 × 105 Pa a1/3 ,
which equals A−1/3 , where A = 10−16 Pa−3 a−1 , which is taken from [49]. Because
the top surface of the periodic domain is flat, we can convert the body force due to
gravity into a constant tangential surface force with magnitude ρg sin θ: we use an
ice density of ρ = 910kgm−3 , and the acceleration due to gravity is g = 9.81m sec−2 .
Note that with this change in the forcing, the pressure p is now interpreted as the
variation from hydrostatic pressure. We use ε = 1 × 10−6 a−2 , which has a negligible
effect for stresses of 105 Pa or greater, following the recommendation in [31, Chapter
2]. The periodic domain is 400 km × 400 km × 1 km, and we again use Txy as our
mesh so that the elements are stretched to φ = 100.
In fig. 7.2, we show the convergence behavior of the inexact Newton method for
k = 3, 4, and 5, and compare with the convergence of an inexact Picard method for k =
3. As can be seen, the inexact Newton method converges faster than Picard’s method,
both in terms of the number of nonlinear iterations and in terms of the total number
of Krylov iterations. As each nonlinear iteration requires a preconditioner setup (or
update), the superiority of Newton’s method compared to the Picard method is even
more pronounced if we consider time-to-solution rather than the number of Krylov
iterations. Although not shown in fig. 7.2, we did test Picard’s method with tighter
tolerances on the linear solves, but found that this did not improve the convergence
rate in terms of total Krylov iterations.
We next test our method on the same problem, but with the Robin coefficient
field reduced to 1% of the previous field, and with rough bed geometries instead of
the flat slab used in the previous test. These modifications make the problem more
realistic, but more challenging to solve. We generate bed topographies using random
coefficients in a truncated Fourier series, as in eq. (4.3). By changing the exponent
γ, which controls the decay of the Fourier coefficients, we are able to control how
rough the generated topography is. In figs. 7.3a and 7.3b, we show two topographies,
generated with ten Fourier modes and γ = 1.5 and γ = 1.0, respectively. In fig. 7.3c we
20

krj k/kr0 k: φ = 100, ISMIP-C, stream variant,
CG(2)/block-Jacobi/IC(0)
k
k
k
k
k
k
k
k

100

10−4

10−8

= 3,
= 3,
= 4,
= 4,
= 5,
= 5,
= 3,
= 3,

Linear
Newton [9]
Linear
Newton [10]
Linear
Newton [10]
Linear
Picard [50]

10−12
0

10

20
30
40
50
60
cumulative Krylov iteration j

70

80

Fig. 7.2: Convergence of the inexact Newton solver for the flat bed topography ice stream
model problem, for different polynomial orders k. Diamonds correspond to the nonlinear
residual at the start of a Newton step. The solid lines show the linear residuals of the Krylov
iterations. The dots show where the nonlinear residuals are evaluated and the preconditioners
are recomputed. The Krylov method is FGMRES(30). The smoother for the (1,1)-block AMG
preconditioner is two iterations of CG, preconditioned by IC(0) on each processor block. As
can be observed, the linear and the nonlinear convergence rates coincide close to the solution,
as desired for an inexact Newton-Krylov method. Far from the solution, where the linear
and nonlinear residuals are very different, the inexactness in the method avoids over-solving
of the linear systems, and is thus efficient in terms of Krylov iterations [18]. The bracketed
numbers in the legend indicate the total number of nonlinear iterations.

show the convergence behavior of our method on domains with these bed topographies.
In these nonlinear problems, the incomplete factorization of the (1,1)-block approximation F̃ sometimes encounters zero or negative pivots on the diagonal, which
can lead to poor convergence or can cause the solver to fail. Zero or negative pivots occur more likely when there are regions of rough topography and sharp solution
gradients. We have considered several strategies to make the preconditioner robust
in these cases. One remedy to avoid these bad pivots is to increase the sophistication
of the incomplete factorization, for instance by increasing the level of fill-in (e.g., use
IC(1) or IC(2) instead of IC(0)), or by using a drop-tolerance-based approach. However, these approaches incur additional setup and storage costs, and it is difficult to
anticipate a priori what amount of fill (or what drop tolerance) is needed to avoid
bad pivots.
PETSc implements strategies for modification of the matrix being factorized to
avoid bad pivots. The default is to use Manteuffel’s shifting strategy [42], in which
F̃ is replaced by F̃ + αI. Our experience is that when the magnitudes of the entries
on the diagonal of F̃ are highly variable, which can be due to variable coefficients or
variable element sizes, this shifting strategy is detrimental to the effectiveness of the
incomplete factorization as a smoother. In particular, it can shift the eigenvalues of the
high-frequency modes away from the region that is optimally damped by a Chebyshev
polynomial smoother, leading to stagnation in the convergence. We find that shifting
21

(a)

(c) krj k/kr0 k: φ = 100, k = 3, ISMIP-C, stream variant,
CG(2)/block-Jacobi/IC(0)[with shifting “in blocks”], two bed topographies
bed (a)
bed (b)
10−3

(b)
10−8

10−13
0

10

20
30
40
50
60
cumulative Krylov iteration j

70

80

Fig. 7.3: (a), (b) Bed topographies with different roughness for the ice stream model
problem. The vertical scale is exaggerated by a factor of 100. The maximum variation in
the bed topography in (b) is about half of the ice thickness. (c) Convergence of the inexact
Newton solver for the ice stream model problem as described in fig. 7.1, but with 1% of the
basal friction β and the rough bed topographies.

“in blocks”—increasing pivots to a lower bound when they become too small—results
more consistently in factorizations that are appropriate for use as smoothers. The
corresponding convergence behavior for problems with polynomial order k = 3 is
shown in fig. 7.3. Note that while compared to fig. 7.2, the efficiency in terms of
overall Krylov is similar, the inexact Newton method requires more iterations before
it reaches the asymptotic convergence regime. This is likely due to a combination of
a greater degree of nonlinearity in the problem caused by the more complex geometry
and the larger variations in the Robin coefficient β.
8. Antarctic ice sheet problem. We now demonstrate the performance of
our inexact Newton-Krylov solver for the simulation of the dynamics of the Antarctic ice sheet. Below, we detail how satellite and radar data are used to define the
computational domain and describe the mesh generation. In section 8.2, we study
the performance of our solver for this problem. Finally, in section 8.3, we compare to
results presented in [33], which are obtained with a similar solver that uses a SSOR
smoother in the preconditioner instead of an incomplete factorization smoother.
8.1. Problem Description. For the following simulations, we define the ice
density to be ρ = 917 kg/m3 , the pre-exponential in Glen’s power law (2.2) to be
B(T ) = 4.1×105 Pa a1/3 , and the regularizing constant that prevents infinite effective
viscosity to be  = 9.95×10−6 a−2 . The assumed sliding coefficient field β is computed
by first taking the ratio of the driving stress due to gravity and the observed surface
velocity, and then imposing maximum and minimum values. This results in β being
∼ 0.3 MPa a km−1 over most of the ice sheet and becomes almost zero (with a
minimum value of ∼ 10−14 MPa a km−1 ) in ice streams.
The body force due to gravity is −ρgz, where z is a unit vector that is normal to
22

the reference ellipsoid. Due to the complexity of the ice sheet geometry, this forcing
cannot be exactly represented by the discrete function spaces we use. Consequently,
we find that when we solve eq. (2.5) with the right-hand side computed by directly
discretizing the body force due to gravity, the velocity in our solution exhibits gridscale oscillations. To ameliorate this, we set the body force to be b = −ρgz − grad p̂,
where p̂ is a hydrostatic-like pressure solving
−∆p̂ = 0

x ∈ Ω,
x ∈ ΓN ,

p̂ = 0
∇p̂ · n = −ρgz · n

x ∈ ΓD ,

∇p̂ · n = −ρgz · n

x ∈ ΓR .

Note that p̂ satisfies homogeneous Dirichlet boundary conditions where the velocity
field has stress-free boundary conditions, and p̂ has Neumann boundary conditions
where the normal velocity component satisfies a Dirichlet condition. The pressure p
that is solved for is then a variation from p̂. Note that the computation of p̂ is an
upfront computation and is not included in any of the performance results that follow.
The geometric description of the ice sheet is constructed from the ALBMAP
dataset [40]. Elevation values in the ALBMAP dataset are given relative to the
EIGEN-GL04C geoid [21]: we convert these values to elevations relative to the WGS84
ellipsoid [1] using the software library GeographicLib [38], and then map the resulting
(latitude, longitude, elevation) geodetic coordinates into Cartesian coordinates. While
this results in a usable geometry for the ice sheet, the resulting geometry for the ice
shelves, the extension of the ice sheet onto the surface of the ocean, is far from
hydrostatic equilibrium, resulting in flow velocities several orders of magnitude too
fast. We were unable to correct this behavior using local geometry adjustments. While
we believe that a good geometry can be obtained, this was not further pursued for
these numerical studies and we have limited the domain to the grounded ice sheet.
From the ice thickness data, given on a latitude-longitude grid, we obtain a polygon describing the lateral boundaries of the ice sheet. We create a quadrilateral mesh
from this polygon by first using the triangular mesh generator Triangle [55], and
then splitting the triangles into quadrilaterals, to which we apply mesh smoothing
to improve the element quality. This resulting coarse quadrilateral mesh (fig. 8.1a)
contains ∼27,000 elements. Using the quadtree-based refinement for the horizontal
directions within our p4est extension for anisotropic domains, we refine this mesh to
construct the footprints for the columns of our final hexahedral mesh. We use several
refinement criteria: we require that the elements of our mesh have a footprint smaller
than (2.5 km)2 at this grounding line; we refine any column whose thickness varies by
more than a factor of 1.5; in keeping with other ice sheet models, we keep the aspect
ratio φ of the elements below 25. Once we have constructed a mesh that satisfies
these constraints (fig. 8.1b), we use uniform refinement (i.e., replace each hexahedron
by eight children) from this mesh when we perform scaling studies.
The Antarctic ice sheet contains some very thin regions which, because we constrain the aspect ratio of our elements, would require a large number of elements if we
modeled the true thickness of the ice sheet. To control the mesh size, we employ an
artificial minimum thickness of 200 m, enforced by modifying the bedrock topography.
8.2. Solver performance and scalability.
Nonlinear convergence for different orders k. We first test our solver for different
polynomial orders on a mesh obtained by one level of uniform mesh refinement. For
23

(a)

(b)

Fig. 8.1: Antarctic ice sheet problem: (a) The coarse quadrilateral mesh description of
the Antarctic ice sheet. (b) An oblique view of the refined hexahedral mesh, obtained using
the “p6est” extension of the p4est library for anisotropic AMR. The refinement near the
center is the result of mesh refinement to reduce the aspect ratio of elements in the mesh:
this type of refinement is not possible using purely octree-based refinement.

krj k/kr0 k: Antarctica, GMRES(2)/block-Jacobi/IC(0)
k
k
k
k
k
k

100

10−4

= 3,
= 3,
= 4,
= 4,
= 5,
= 5,

Linear
Newton [10]
Linear
Newton [9]
Linear
Newton [15]

10−8

10−12
0

50
100
150
200
cumulative Krylov iteration j

250

Fig. 8.2: Antarctic ice sheet problem: convergence of the incomplete Newton-Krylov
solver for different polynomial element orders k. We enforce a minimal ice thickness of 200
m, and the mesh is distributed over 1024 processes. The bracketed numbers are the total
number of nonlinear iterations.

smoothing in the (1,1)-block multigrid preconditioner, we use 2 GMRES iterations of
our block-Jacobi/IC(0) smoother. We discretize eq. (2.5) using Qk × Qk−2 elements
for k = 3, 4, and 5, resulting in problems with 51M, 121M, and 238M degrees of
freedom. The experiments were conducted on TACC’s Stampede supercomputer,
with each MPI process assigned to one Sandy Bridge Xeon core. Each discretization
is distributed across 1024 MPI processes. The convergence of our inexact NewtonKrylov solver is shown in fig. 8.2. Note that the convergence is similar to the results for
the rough bed topography model problem shown in fig. 7.3c. The increased number
of overall Krylov iterations (by about a factor of 2) is likely due to the more complex
geometry and boundary conditions.
Parallel scalability. We next test the parallel scalability of our solver on the
Q3 × Q1 discretization. We use three different levels of uniform mesh refinement,
` = 0, 1, and 2. In conducting our study, we found that our GMRES(2)/block24

Table 8.1: Antarctic ice sheet problem: scaling of the full nonlinear solver for k = 3 for
the geometry with minimal ice thickness of 200 m, using GMRES(10)/block-Jacobi/IC(0) as
a smoother in the (1,1)-block multigrid preconditioner. For different numbers of processes P ,
we report the total time (in seconds) to solve the nonlinear problem to a relative tolerance
of 10−12 in the `2 -norm from a zero initial guess, the parallel efficiency (eff.) as well as
the number of Newton iterations (#N) and overall Krylov iterations (#K) performed during
the solution. We also report the average time and efficiency for a single Stokes operator
application (Op), a preconditioner application (PC), and the total preconditioner setup time
(setup). Factorization of the coarsest operator (which has always less than 50 unknowns) is
included in the setup times, and applying the inverse of the factored coarse operator is included
in preconditioner application times. We report efficiency in reference to the smallest problem
size on the smallest number of processors.

P

Solve

128
256
512

67.7
36.9
20.3

1,024
2,048
4,096

78.0
44.2
30.2

8,192
16,384

108.0
74.7

eff.

#N #K
Op
eff.
PC
` = 0 uniform refinement, Ndof =7M
1.00
7
66
0.0275 1.00 0.981
0.91
7
67
0.0146 0.94 0.577
0.83
7
65
0.0080 0.86 0.299
` = 1 uniform refinement, Ndof =51M
0.86 11
75
0.0276 0.99 0.987
0.76 10
75
0.0152 0.90 0.561
0.56 10
75
0.0087 0.79 0.383
` = 2 uniform refinement, Ndof =383M
0.62 13
91
0.0295 0.93
1.13
0.45 10
89
0.0168 0.82
0.80

eff.

Setup

eff.

1.00
0.93
0.82

7.361
4.255
2.722

1.00
0.86
0.67

0.99
0.87
0.64

8.657
5.810
5.406

0.85
0.63
0.34

0.87
0.61

15.58
17.22

0.47
0.21

Jacobi/IC(0) smoother for the (1,1)-block multigrid preconditioner was sometimes
insufficient: because of the non-smooth coefficients obtained for some Newton steps,
the prolongation of coarse solutions onto fine grids would introduce errors that could
not be sufficiently damped by just a couple of smoother iterations. For this reason,
we employ a stronger GMRES(10)/block-Jacobi/IC(0) smoother in this scaling study.
We report the results in table 8.1. The table includes timings for each of the main
components of the nonlinear solver: the matrix-vector product, the preconditioner
application, and the preconditioner reconstruction.
The time to apply the preconditioner is almost entirely spent in applying the
AMG V-cycle to the (1,1)-block of the Stokes operator. This step also requires more
communication than the matrix-vector product, both to project and restrict vectors
in the hierarchy and to apply the smoothers. Note that GMRES requires global
reductions to compute the required inner products, which would not be the case for
a stationary smoother (e.g. Chebyshev).
It is well known that the setup phase in parallel implementations of algebraic
multigrid requires significant communication to properly aggregate degrees of freedom
across processor boundaries, compute prolongation matrices and repartition coarse
matrices. As a consequence, the setup is often the least scalable component of the
solver [8, 12], which can also be seen in table 8.1. Here, the reported times include
both the prolongator constructions in the initial setup and the Galerkin projections
to reconstruct the coarse matrices during subsequent Newton steps, so the increase in
Newton iterations for the larger problems affects the reported efficiency. In all cases,
the contribution of the initial prolongator setup is roughly equivalent to the cost of
25

one subsequent preconditioner reconstruction.
Note that though the total number of Krylov iterations increases by ∼ 37% from
the smallest to the largest problem, the number of outer inexact Newton steps almost
doubles from 7 to 13. Recent work [7] on solving the hydrostatic approximation to
the Stokes equations for ice sheet dynamics has demonstrated the effectiveness of
grid continuation in obtaining an initial guess that is near the region of asymptotic
convergence of Newton’s method. The hierarchical mesh refinement we use lends itself
naturally to grid continuation, so this approach could improve the efficiency of our
nonlinear solver.
8.3. Comparison with performance of SSOR-based smoothing. In [33],
we have used a similar nonlinear and linear solver framework as the one developed in
this paper to infer the Robin coefficient field β in the Stokes boundary value problem
such that the velocity fields of the solution closely matches satellite observations of
the Antarctic ice sheet’s surface (see fig. 8.3). This inference uses methods of adjointbased PDE-constrained optimization to find an optimal β, which requires the solution
of the Stokes problem described in this work, as well as linearized adjoint problems,
whose operators are similar to the linearized Stokes operator in (2.6).
The difference between the solver used in [33] and the one presented here lies in
the preconditioning for the (1,1)-block. In [33], a Chebyshev(2)/block-Jacobi/SSOR
smoother combined with the standard SA-AMG scheme (not the column-preserving
SA-AMG presented in section 5.1) is used. In table 8.2, we reproduce some of the
algorithmic scalings from [33], and compare them to scalings obtained for the same
problem with the solvers developed here: using a Chebyshev(2)/block-Jacobi/IC(0)
smoother and column-preserving SA-AMG. We present only algorithmic scalability
because the results in [33] were computed on Oak Ridge National Laboratory’s Titan supercomputer, while the new scalings were computed on TACC’s Stampede
supercomputer. As a consequence, the timings are not directly comparable. We
note, however, that the cost of applying an IC(0) preconditioner, in terms of memory
movement and floating point operations, is almost identical to the cost of applying
an SSOR preconditioner, so the only significant additional computation in our new
scaling results is in the computation of the IC(0) factorization, which is on par with a
single application of the IC(0) preconditioner. Thus, the algorithmic speedup shown
in table 8.2 is representative of the runtime speedup we would find, had these scaling
results been produced on the same architecture.
We note that the total number of iterations required to solve these problems is
lower than for the problems solved earlier in fig. 8.2, even though both are posed on
the same domain. This is because the aspect ratio of the elements is smaller, with all
element aspect ratios φ less than 10, and because the basal friction is larger, reaching
∼ 3 MPa a km−1 in some areas (versus the ∼ 0.3 MPa a km−1 used in fig. 8.2), and
having a minimum value in the ice streams that is three orders of magnitude larger
than the minimum used in fig. 8.2.
9. Conclusions. Several issues related to high-order, adaptive mesh discretizations and solvers for the simulation of nonlinear Stokes flow in three-dimensional
anisotropic domains are addressed in this work. Our main target problem has been the
nonlinear Stokes boundary value problem arising in ice sheet dynamics. We demonstrate an extension to the p4est library for adaptive mesh refinement of anisotropic
domains that combines quadtree-based refinement in the horizontal directions with
columns of elements to achieve a flexible approach to mesh refinement, with local control over the aspect ratio of elements in the mesh. We demonstrate that high-order
26

Fig. 8.3: Antarctic ice sheet problem: magnitude of the surface velocity field optimized to
match satellite observations. The inversion procedure to infer the basal coefficient is described
in [33].
Table 8.2: Antarctic ice sheet problem: comparison of algorithmic scalability of the
solver used in [33] (which uses SSOR and standard SA-AMG within the (1,1)-block preconditioner) versus the one presented here (which uses IC(0) and column-preserving SA-AMG).
Shown are the number of overall degrees of freedom Ndof , the number of processes P , the
number of Newton iterations #N and the number of cumulative Krylov iterations #K. The
data on the left is reproduced from [33].

Ndof
38M
268M

P
1,024
8,192

SSOR, standard SA-AMG
#N
#K
8
147
9
243

IC(0), column-preserving SA-AMG
#N
#K
7
85
8
98

finite elements discretized on these meshes are well-approximated by low-order approximations, which can be used for preconditioning. We present an efficient solver for the
linearized Stokes equations, with particular emphasis on the design of algebraic multigrid solvers for high-order discretizations, anisotropic domains, and hanging nodes.
Using incomplete factorization-based smoothing for the (1,1)-block yields efficient and
fast convergence. When using this type of smoothing, we demonstrate that columnpreserving SA-AMG, as implemented by our DofColumns plugin for PETSc’s GAMG
preconditioner, significantly improves over standard SA-AMG in its effectiveness, especially in the presence of weak boundary conditions. The numerical experiments on
our discretization of the Antarctic ice sheet show that, up to a point, the incomplete
factorization process can be made robust against bad pivots that occur in problems
with variable coefficients and rough topography by using only local shifting and by
combining the incomplete factorization with a non-stationary Krylov method as a
smoother.
10. Acknowledgments. We would like to thank the three referees who reviewed
this work, whose reviews greatly improved the final version of this work. In particular, we would like to thank the two initial referees, Ray Tuminaro and one who
remains anonymous: the column-preserving SA-AMG that we have developed in the
DofColumns plugin was not included in the initial submission of this work, but was
developed in response to their insightful comments and suggestions.
Appendix A. Well-posedness of eq. (2.5). The main difference
between
R
eq. (2.5) and the variational form in [35] is the boundary integral ΓR βTk v · Tk u ds:
27

for u ∈ W1,r (Ω), this form is not meaningful because the trace of u is not necessarily in
L2 (ΓR ). We compensate for this by bringing the boundary integral into the definition
of V. Let [C∞ (Ω)]30 be the space of smooth vector-valued functions in Ω that satisfy
the homogeneous Dirichlet part of the boundary conditions and let r = 1 + n1 . Under
the above assumptions, the functional
I(u) =

Z

r

1/r

|∇u| dx

Z
+

Ω

1/2
β|Tk u| ds
2

(A.1)

ΓR

defines a seminorm on [C∞ (Ω)]30 . We assume a problem setup in which I(u) also
defines a norm, which amounts to requiring that ifR u ∈ [C∞ (Ω)]30 is a rigid-body
motion that satisfies the Dirichlet conditions, then ΓR β|Tk u|2 ds > 0. We define V
to be the closure of [C∞ (Ω)]30 in this norm.
Following the same steps as in [35], one can show that the minimization problem
Z
Z
1+n
2n
1
2n
B(T )(DII (u) + ε)
Tk u · Tk u ds − f (u)
(A.2)
dx + 2
min
u
1+n Ω
ΓR
is well-posed in Vdiv , the subspace of V containing only the divergence-free functions.
To prove that eq. (2.5) is well-posed, we additionally need to choose a space M for
which we can prove the inf-sup condition
R
q∇ · u dx
inf sup Ω
≥ γ > 0.
(A.3)
q∈M u∈V kqkM kukV
0

This inequality still holds for M = Lr (Ω). To see this, consider the subspace Ṽ =
{u ∈ V : u|ΓR = 0}: for u ∈ Ṽ, kukV = kuk[W1,r (Ω)]3 , and so
R
inf0

sup

q∈Lr (Ω) u∈V

R
q∇ · u dx
q∇ · u dx
Ω
≥ inf0
sup
.
kqkM kukV
q∈Lr (Ω) u∈Ṽ kqkLr0 (Ω) kuk[W1,r(Ω) ]3
Ω

(A.4)

As long as ΓN 6= ∅, the term on the right is bounded from below as a particular case
of the inf-sup condition in [35]. Thus, (2.5) is well posed and has a unique solution.
REFERENCES
[1] Department of Defense World Geodetic System 1984: Its definition and relationships with
local geodetic systems, Tech. Rep. NGA TR8350.2, U.S. National Geospatial-Intelligence
Agency, 1984. Revised in 2004.
[2] M. F. Adams, A parallel maximal independent set algorithm, in Proceedings 5th copper mountain conference on iterative methods, 1998.
[3] M. Ainsworth and P. Coggins, The stability of mixed hp-finite element methods for Stokes
flow on high aspect ratio elements, SIAM Journal on Numerical Analysis, 38 (2000),
pp. 1721–1761.
[4] S. Balay, J. Brown, K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G.
Knepley, L. C. McInnes, B. F. Smith, and H. Zhang, PETSc users manual, Tech. Rep.
ANL-95/11 - Revision 3.3, Argonne National Laboratory, 2012.
[5] H. Blatter, Velocity and stress fields in grounded glaciers: a simple algorithm for including
deviatoric stress gradients, Journal of Glaciology, 41 (1995), pp. 333–344.
[6] J. Brown, Efficient nonlinear solvers for nodal high-order finite elements in 3D, Journal of
Scientific Computing, 45 (2010), pp. 48–63.
[7] J. Brown, B. Smith, and A. Ahmadia, Achieving textbook multigrid efficiency for hydrostatic
ice sheet flow, SIAM Journal on Scientific Computing, 35 (2013), pp. B359–B375.
28

[8] C. Burstedde, O. Ghattas, G. Stadler, T. Tu, and L. C. Wilcox, Parallel scalable adjointbased adaptive solution for variable-viscosity Stokes flows, Computer Methods in Applied
Mechanics and Engineering, 198 (2009), pp. 1691–1700.
[9] C. Burstedde, L. C. Wilcox, and O. Ghattas, p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees, SIAM Journal on Scientific Computing, 33
(2011), pp. 1103–1133.
[10] M. Cai, A. Nonaka, J. B. Bell, B. E. Griffith, and A. Donev, Efficient variablecoefficient finite-volume Stokes solvers, Communications in Computational Physics, 16
(2014), pp. 1263–1297.
[11] C. Canuto, P. Gervasio, and A. Quarteroni, Finite-element preconditioning of G-NI spectral methods, SIAM Journal on Scientific Computing, 31 (2010), pp. 4422–4451.
[12] E. Chow, R. D. Falgout, J. J. Hu, R. S. Tuminaro, and U. M. Yang, A survey of parallelization techniques for multigrid solvers, in Parallel Processing for Scientific Computing,
M. A. Heroux, P. Raghavan, and H. D. Simon, eds., Society for Industrial and Applied
Mathematics, Philadephia, PA, 2006, pp. 179–201.
[13] S. L. Cornford, D. F. Martin, D. T. Graves, D. F. Ranken, A. M. Le Brocq, R. M.
Gladstone, A. J. Payne, E. G. Ng, and W. H. Lipscomb, Adaptive mesh, finite volume
modeling of marine ice sheets, Journal of Computational Physics, 232 (2012), pp. 529–549.
[14] H. Damanik, J. Hron, A. Ouazzi, and S. Turek, A monolithic fem-multigrid solver for nonisothermal incompressible flow on general meshes, Journal of Computational Physics, 228
(2009), pp. 3869–3881.
[15] M. Deville and E. Mund, Finite-element preconditioning for pseudospectral solutions of elliptic problems, SIAM Journal on Scientific and Statistical Computing, 11 (1990), pp. 311–342.
[16] M. O. Deville, P. F. Fischer, and E. H. Mund, High-Order Methods for Incompressible
Fluid Flow, vol. 9 of Cambridge Monographs on Applied and Computational Mathematics,
Cambridge University Press, Cambridge, UK, 2002.
[17] M. Durufle, P. Grob, and P. Joly, Influence of Gauss and Gauss-Lobatto quadrature rules
on the accuracy of a quadrilateral finite element method in the time domain, Numerical
methods for partial differential equations, 25 (2009), pp. 526–551.
[18] S. C. Eisenstat and H. F. Walker, Choosing the forcing terms in an inexact Newton method,
SIAM Journal on Scientific Computing, 17 (1996), pp. 16–32.
[19] H. Elman, V. Howle, J. Shadid, R. Shuttleworth, and R. Tuminaro, Block preconditioners based on approximate commutators, SIAM Journal on Scientific Computing, 27 (2006),
pp. 1651–1668.
[20] H. C. Elman, D. J. Silvester, and A. J. Wathen, Finite Elements and Fast Iterative Solvers
with applications in incompressible fluid dynamics, Oxford University Press, Oxford, 2005.
[21] C. Förste, R. Schmidt, R. Stubenvoll, F. Flechtner, U. Meyer, R. König, H. Neumayer, R. Biancale, J.-M. Lemoine, S. Bruinsma, et al., The Geoforschungszentrum
Potsdam/groupe de Recherche de Geodesie Spatiale satellite-only and combined gravity
field models: EIGEN-GL04S1 and EIGEN-GL04C, Journal of Geodesy, 82 (2008), pp. 331–
346.
[22] M. Furuichi, D. A. May, and P. J. Tackley, Development of a Stokes flow solver robust to
large viscosity jumps using a Schur complement approach with mixed precision arithmetic,
Journal of Computational Physics, 230 (2011), pp. 8835–8851.
[23] O. Gagliardini and T. Zwinger, The ISMIP-HOM benchmark experiments performed using
the finite-element code Elmer, The Crysophere, 2 (2008), pp. 67–76.
[24] M. W. Gee, J. J. Hu, and R. S. Tuminaro, A new smoothed aggregation multigrid method
for anisotropic problems, Numer. Linear Algebra Appl, 16 (2009), pp. 19–37.
[25] T. Geenen, M. ur Rehman, S. P. MacLachlan, G. Segal, C. Vuik, A. P. van den Berg,
and W. Spakman, Scalable robust solvers for unstructured FE geodynamic modeling applications: Solving the Stokes equation for models with large localized viscosity contrasts,
Geochemistry Geophysics Geosystems, 10 (2009), p. Q09002.
[26] R. Glowinski and J. Xu, Numerical Methods for Non-Newtonian Fluids: Special Volume,
vol. 16 of Handbook of numerical analysis, North-Holland, 2011.
[27] P. P. Grinevich and M. A. Olshanskii, An iterative method for the Stokes-type problem with
variable viscosity, SIAM Journal on Scientific Computing, 31 (2009), pp. 3959–3978.
[28] V. Heuveline and F. Schieweck, On the inf-sup condition for higher order mixed fem on
meshes with hanging nodes, ESAIM: Mathematical Modelling and Numerical Analysis, 41
(2007), pp. 1–20.
[29] J. Heys, T. Manteuffel, S. McCormick, and L. Olson, Algebraic multigrid for higher-order
finite elements, Journal of computational Physics, 204 (2005), pp. 520–532.
[30] R. C. A. Hindmarsh, A numerical comparison of approximations to the Stokes equations used
29

in ice sheet and glacier modeling, Journal of Geophysical Research, 109 (2004), p. F01012.
[31] K. Hutter, Theoretical Glaciology, Mathematical Approaches to Geophysics, D. Reidel Publishing Company, 1983.
[32] T. Isaac, C. Burstedde, and O. Ghattas, Low-cost parallel algorithms for 2:1 octree balance, in Proceedings of the 26th IEEE International Parallel & Distributed Processing
Symposium, IEEE, 2012. http://dx.doi.org/10.1109/IPDPS.2012.47.
[33] T. Isaac, N. Petra, G. Stadler, and O. Ghattas, Scalable and efficient algorithms for
the propagation of uncertainty from data through inference to prediction for large-scale
problems, with application to flow of the Antarctic ice sheet, Journal of Computational
Physics, 296 (2015), pp. 348–368.
[34] A. Janka, Smoothed aggregation multigrid for a Stokes problem, Computing and Visualization
in Science, 11 (2008), pp. 169–180.
[35] G. Jouvet and J. Rappaz, Analysis and finite element approximation of a nonlinear stationary Stokes problem arising in glaciology, Advances in Numerical Analysis, 2011 (2012).
ID 164581.
[36] Y. Kanarska, A. Shchepetkin, and J. McWilliams, Algorithm for non-hydrostatic dynamics
in the regional oceanic modeling system, Ocean Modelling, 18 (2007), pp. 143–174.
[37] E. Karer and J. Kraus, Algebraic multigrid for finite element elasticity equations: Determination of nodal dependence via edge-matrices and two-level convergence, International
Journal for Numerical Methods in Engineering, 83 (2010), pp. 642–670.
[38] C. Karney, Algorithms for geodesics, Journal of Geodesy, 87 (2013), pp. 43–55.
[39] S. Kim, Piecewise bilinear preconditioning of high-order finite element methods, Electronic
Transactions on Numerical Analysis, 26 (2007), pp. 228–242.
[40] A. M. Le Brocq, A. J. Payne, and A. Vieli, An improved Antarctic dataset for high resolution
numerical ice sheet models (ALBMAP v1), Earth System Science Data, 2 (2010), pp. 247–
260.
[41] W. Leng, L. Ju, M. Gunzburger, S. Price, and T. Ringler, A parallel high-order accurate
finite element nonlinear Stokes ice-sheet model and benchmark experiments, Journal of
Geophysical Research, 117 (2012).
[42] T. A. Manteuffel, An incomplete factorization technique for positive definite linear systems,
Mathematics of computation, 34 (1980), pp. 473–497.
[43] D. A. May and L. Moresi, Preconditioned iterative methods for Stokes flow problems arising
in computational geodynamics, Physics of the Earth and Planetary Interiors, 171 (2008),
pp. 33–47.
[44] D. McKenzie, The generation and compaction of partially molten rock, Journal of Petrology,
25 (1984), pp. 713–765.
[45] G. A. Meehl, T. F. Stocker, W. D. Collins, A. T. Friedlingstein, A. T. Gaye, J. M.
Gregory, A. Kitoh, R. Knutti, J. M. Murphy, A. Noda, S. C. B. Raper, I. G. Watterson, A. J. Weaver, and Z.-C. Zhao, Global climate projections, in Climate Change
2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University
Press, 2007, pp. 747–845.
[46] M. Morlighem, H. Seroussi, E. Larour, and E. Rignot, Inversion of basal friction in
antarctica using exact and incomplete adjoints of a higher-order model, Journal of Geophysical Research: Earth Surface, 118 (2013), pp. 1746–1753.
[47] J. Nocedal and S. J. Wright, Numerical Optimization, Springer Verlag, Berlin, Heidelberg,
New York, 1999.
[48] L. Olson, Algebraic multigrid preconditioning of high-order spectral elements for elliptic problems on a simplicial mesh, SIAM Journal on Scientific Computing, 29 (2007), pp. 2189–
2209.
[49] F. Pattyn, L. Perichon, A. Aschwanden, B. Breuer, B. de Smedt, O. Gagliardini,
G. H. Gudmundsson, R. C. A. Hindmarsh, A. Hubbard, J. V. Johnson, T. Kleiner,
Y. Konovalov, C. Martin, A. J. Payne, D. Pollard, S. Price, M. Ruckamp, F. Saito,
O. Soucek, S. Sugiyama, and T. Zwinger, Benchmark experiments for higher-order and
full-Stokes ice sheet models (ISMIP-HOM), The Cryosphere, 2 (2008), pp. 95–108.
[50] N. Petra, H. Zhu, G. Stadler, T. J. R. Hughes, and O. Ghattas, An inexact Gauss-Newton
method for inversion of basal sliding and rheology parameters in a nonlinear Stokes ice
sheet model, Journal of Glaciology, 58 (2012), pp. 889–903.
[51] Y. Saad, A flexible inner-outer preconditioned GMRES algorithm, SIAM Journal on Scientific
Computing, 14 (1993), pp. 461–469.
[52] Y. Saad, Iterative methods for sparse linear systems, Society for Industrial and Applied Mathematics, Philadelphia, PA, second ed., 2003.
30

[53] G. Schubert, D. Turcotte, and P. Olsen, Mantle Convection in the Earth and Planets,
Cambridge University Press, 2001.
[54] H. Seroussi, H. B. E. N. Dhia, M. Morlighem, E. Larour, E. Rignot, and D. Aubry,
Coupling ice flow models of varying orders of complexity with the Tiling method, Journal
of Glaciology, 58 (2012).
[55] J. R. Shewchuk, Triangle: Engineering a 2D quality mesh generator and delaunay triangulator, in Applied Computational Geometry: Towards Geometric Engineering, M. C. Lin and
D. Manocha, eds., vol. 1148 of Lecture Notes in Computer Science, Springer-Verlag, 1996,
pp. 203–222. From the First ACM Workshop on Applied Computational Geometry.
[56] I. K. Tezaur, M. Perego, A. G. Salinger, R. S. Tuminaro, and S. F. Price, Albany/felix:
a parallel, scalable and robust, finite element, first-order Stokes approximation ice sheet
solver built for advanced analysis, Geoscientific Model Development, 8 (2015), pp. 1197–
1220.
[57] A. Toselli and C. Schwab, Mixed hp-finite element approximations on geometric edge and
boundary layer meshes in three dimensions, Numerische Mathematik, 94 (2003), pp. 771–
801.
[58] U. Trottenberg, C. Oosterlee, and A. Schüller, Multigrid, Academic Press, London,
2001.
[59] P. Vaněk, M. Brezina, and J. Mandel, Convergence of algebraic multigrid based on smoothed
aggregation, Numerische Mathematik, 88 (2001), pp. 559–579.
[60] P. Vaněk, J. Mandel, and M. Brezina, Algebraic multigrid by smoothed aggregation for
second and fourth order elliptic problems, Computing, 56 (1996), pp. 179–196.
[61] A. Zmitrowicz, Glaciers and laws of friction and sliding, Acta mechanica, 166 (2003), pp. 185–
206.

31

