Weighted BFBT Preconditioner for Stokes Flow Problems with
Highly Heterogeneous Viscosity∗

arXiv:1607.03936v2 [math.NA] 30 Jan 2017

Johann Rudi†

Georg Stadler‡

Omar Ghattas†§

January 31, 2017

Abstract
We present a weighted BFBT approximation (w-BFBT) to the inverse Schur complement of a
Stokes system with highly heterogeneous viscosity. When used as part of a Schur complement-based
Stokes preconditioner, we observe robust fast convergence for Stokes problems with smooth but highly
varying (up to 10 orders of magnitude) viscosities, optimal algorithmic scalability with respect to mesh
refinement, and only a mild dependence on the polynomial order of high-order finite element discretizations
(Qk × Pdisc
k−1 , order k ≥ 2). For certain difficult problems, we demonstrate numerically that w-BFBT
significantly improves Stokes solver convergence over the widely used inverse viscosity-weighted pressure
mass matrix approximation of the Schur complement. In addition, we derive theoretical eigenvalue
bounds to prove spectral equivalence of w-BFBT. Using detailed numerical experiments, we discuss
modifications to w-BFBT at Dirichlet boundaries that decrease the number of iterations. The overall
algorithmic performance of the Stokes solver is governed by the efficacy of w-BFBT as a Schur complement
approximation and, in addition, by our parallel hybrid spectral-geometric-algebraic multigrid (HMG)
method, which we use to approximate the inverses of the viscous block and variable-coefficient pressure
Poisson operators within w-BFBT. Building on the scalability of HMG, our Stokes solver achieves a
parallel efficiency of 90% while weak scaling over a more than 600-fold increase from 48 to all 30,000 cores
of TACC’s Lonestar 5 supercomputer.

1
1.1

Introduction
Motivation and governing equations

Many problems in science and engineering involve creeping flows of non-Newtonian fluids [19, 35]. Important
examples can be found in geophysical fluid flows, where the incompressible Stokes equations with power-law
rheology have become a prototypical continuum mechanical description for creeping flows occurring in mantle
convection [37], magma dynamics [32], and ice flow [22]. The linearization of the nonlinear momentum
equations, for instance within a Newton method, leads to incompressible Stokes-like equations with highly
heterogeneous viscosity fields.
In particular, simulations of Earth’s mantle convection at global scale [39] exhibit extreme computational
challenges due to a highly heterogeneous viscosity stemming from its dependence on temperature and strain
rate as well as sharp viscosity gradients in narrow regions modeling tectonic plate boundaries (six orders of
magnitude drop in ∼5 km) [6, 36]. This leads to a wide range of spatial scales since small localized features at
plate boundaries of size O(1 km) influence plate motion at continental scales of O(1000 km). The complex
character of the flow presents severe computational challenges for iterative solvers due to poor conditioning
of the linear systems that arise.
∗ This research was partially supported by NSF grants CMMI-1028889 and ARC-0941678 and U.S. Department of Energy
grants DE-FC02-13ER26128 and the Scientific Discovery through Advanced Computing (SciDAC) program.
† Institute for Computational Engineering and Sciences, The University of Texas at Austin, USA (johann@ices.utexas.edu,
omar@ices.utexas.edu).
‡ Courant Institute of Mathematical Sciences, New York University, USA (stadler@cims.nyu.edu).
§ Jackson School of Geosciences and Department of Mechanical Engineering, The University of Texas at Austin, USA.

1

w-BFBT Preconditioner for Stokes Flow Problems

2

Since we focus here on preconditioning linearized Stokes-like problems that arise at each step of a nonlinear
solver, we can simplify our problem setup by taking the viscosity to be independent of the strain rate, but
otherwise exhibiting severe spatial heterogeneity. Given a bounded domain Ω ⊂ Rd , d ∈ {2, 3}, right-hand side
forcing f (x), and spatially-varying viscosity µ(x) ≥ µmin > 0 for all x ∈ Ω, we consider the incompressible
Stokes equations with homogeneous Dirichlet boundary conditions


−∇ · µ(x) (∇u + ∇uT ) + ∇p = f in Ω,
(1a)
−∇ · u = 0
u=0

in Ω,

(1b)

on ∂Ω,

(1c)

where u and p are the unknown velocity and pressure fields, respectively.1

1.2

Discretization and computational challenges

Discretizing (1) leads to a linear algebraic system of equations of the form

   
f
A BT u
=
,
p
0
B 0

(2)

where A, B, and BT are matrices corresponding to discretizations of the viscous stress, divergence, and gradient
operators, respectively. The discretization is carried out by high-order finite elements on (possibly aggressively
adaptively refined) hexahedral meshes with velocity–pressure pairings Qk × Pdisc
k−1 of polynomial order k ≥ 2
with a continuous, nodal velocity approximation Qk and a discontinuous, modal pressure approximation
Pdisc
k−1 . These pairings yield optimal asymptotic convergence rates of the finite element approximation to the
infinite-dimensional solution with decreasing mesh element size, are inf-sup stable on general, non-conforming
hexahedral meshes with “hanging nodes,” and have the advantage of preserving mass locally at the element
level due to the discontinuous pressure [16,21,40]. While these properties have been recognized to be important
for geophysics applications (e.g., see [29, 30]), the high-order discretization, adaptivity, and discontinuous
pressure approximation present significant additional difficulties for iterative solvers (relative to low order,
uniform grid, continuous discretizations). Finally, a number of frontier geophysical problems, such as global
mantle convection with plate boundary-resolving meshes and continental ice sheet models with grounding
line-resolving meshes, result in billions of degrees of freedom, demanding efficient execution and scalability on
leading edge supercomputers [6, 36, 39].
The applications we target (such as global mantle convection) exhibit all of the difficulties described
above (severe heterogeneity, very large scale, need for aggressively-adapted meshes, need for high order, massconserving discretization) and thus demand robust and effective preconditioners for (2), resulting in iterative
solvers with optimal (or nearly optimal) algorithmic and parallel scalability. This paper describes the design
of such a preconditioner and its analysis and performance evaluation for problems with highly heterogeneous
viscosity. The preconditioner—which we call weighted BFBT (w-BFBT)—is of Schur complement type, and
we study its robustness as well as its algorithmic and parallel scalability.

1.3

Iterative methods and Schur complement approximations

An effective approximation of the Schur2 complement S := BA−1 BT is an essential ingredient for attaining
fast convergence of Schur complement-based iterative solvers for (2). More precisely, a sufficiently good
approximation of the inverse Schur complement S̃−1 ≈ S−1 is sought, which, together with an approximation
1 Note that Newton linearization (as opposed to, e.g., Picard) results in a fourth-order anisotropic viscosity tensor, with the
anisotropic component proportional to the tensor product of the second-order strain rate tensor with itself. This anisotropic
viscosity tensor also appears in adjoint Stokes equations, which are important for inverse problems [23, 34, 43]. For power-law
rheologies that are typical of many geophysical flows, the anisotropic component of the viscosity tensor is dominated by
the isotropic component µ(x)I [24, 36]. Thus for many forward and inverse modeling applications, the performance of our
preconditioner for isotropic viscosities is indicative of its performance for operators with anisotropic viscosity tensors. Hence we
focus on the former in this paper.
2 Strictly speaking, our definition is the negative Schur complement. However, as in [16], we prefer to work with positive-definite
operators and thus define the Schur complement to be positive rather than negative definite.

w-BFBT Preconditioner for Stokes Flow Problems

3

of the inverse viscous block, Ã−1 ≈ A−1 , is used in an iterative scheme with right preconditioning based on
an upper triangular block matrix:

A
B

BT
0


Ã
0

BT
S̃

−1    
ũ
f
=
.
p̃
0

(3)

Note that the original solution to (2) is recovered by applying the preconditioner once to the solution of (3).
For the preconditioned Stokes system (3), we use GMRES as the Krylov subspace solver. This particular
combination of Krylov method and preconditioner is known to converge in just two iterations for exact choices
of Ã−1 and S̃−1 [4].
The most widely used approximation of the Schur complement (for variable viscosity Stokes systems) is
the inverse viscosity-weighted
mass matrix of the pressure space [8, 24, 27, 30], denoted by Mp (1/µ), with
R
entries [Mp (1/µ)]i,j = Ω qi (x) qj (x)/µ(x) dx, where qi , qj ∈ Pdisc
k−1 are global basis functions of the finite
disc
disc
dimensional space Pk−1 . Since the basis functions of Pk−1 are modal and not orthogonal to each other, the
mass matrix is not diagonal, and thus Mp (1/µ) is typically diagonalized to further simplify its inversion. One
common way to obtain a diagonalized version is mass lumping. For nodal discretizations, the corresponding
diagonal elements are computed by summation of the entries of each matrix row, i.e., Mp (1/µ)1, where 1 is
the vector with ones in all entries. For modal discretizations, we generalize the lumping procedure by using
the coefficient vector, 1{qi }i , representing the constant function having value 1 in the associated basis {qi }i ,
i.e.,
M̃p (1/µ) := diag(Mp (1/µ)1{qi }i ).
(4)
Provided that µ is sufficiently smooth, Mp (1/µ) can be an effective approximation of S in numerical
experiments [7] and spectral equivalence can be shown [20]. However, it has been observed in applications
with highly heterogeneous viscosities µ (e.g., mantle convection [31, 36]) that convergence slows down
significantly due to a poor Schur complement approximation by Mp (1/µ). Therefore, we propose a new
approximation, w-BFBT, that remains robust when Mp (1/µ) fails.
Preconditioners based on BFBT approximations for the Schur complement were initially proposed in [13]
for the Navier–Stokes equations. Over the years, these ideas were refined and extended [14, 15, 17, 26, 38]
to arrive at a class of closely related Schur complement approximations: Pressure Convection–Diffusion,
BFBT, and Least Squares Commutator. The underlying principle, now in a Stokes setting, is that one seeks
a commutator matrix X such that the following commutator nearly vanishes,
AD−1 BT − BT X ≈ 0,

(5)

for a given diagonal matrix D−1 . The Navier–Stokes case differs from Stokes in that the viscous stress matrix
A contains an additional convection term. The motivation for seeking a near-commutator X is that (5)
can be rearranged by multiplying (5) with BA−1 from the left and, provided the inverse exists, with X−1
from the right to obtain S ≈ BD−1 BT X−1 , where the closer the commutator is to zero, the more accurate
the approximation [16]. The goal of finding a vanishing commutator can be recast as solving the following
least-squares minimization problem:
Find matrix X minimizing

AD−1 BT ej − BT Xej

2
C−1

for all j,

(6)

where ej is the j-th Cartesian unit vector and the norm arises from a symmetric and positive definite matrix
C. The solution is given by X = (BC−1 BT )−1 (BC−1 AD−1 BT ). Then the BFBT approximation of the
inverse Schur complement is derived by algebraic rearrangement of the commutator (5):
−1 T
S̃−1
B
BFBT := BC

−1

BC−1 AD−1 BT



BD−1 BT

−1

.

(7)

In the literature cited above (which addresses preconditioning of Navier–Stokes equations with constant
viscosity), the diagonal weighting matrices are chosen as C = D = M̃u , i.e., a diagonalized version of
the velocity space mass matrix; hence we call this the Mu -BFBT approximation of the Schur complement.
Mu -BFBT can be used for Stokes problems with constant viscosities providing convergence similar to
Mp (1/µ). However, the computational cost of applying (7) is significantly higher than the (cheap) application
of a possibly diagonalized inverse of Mp (1/µ). Moreover, Mu -BFBT is not an option for heterogeneous

Number of GMRES iterations

w-BFBT Preconditioner for Stokes Flow Problems

4

700
Mp (1/µ)
w-BFBT

600
500
400
300
200
100
0

0

4

8

12

16

20

24

28

Problem difficulty (number of sinkers) →

Figure 1: Left image shows the improvement in convergence obtained with the proposed w-BFBT preconditioner over a preconditioner using the inverse viscosity-weighted pressure mass matrix as Schur complement
approximation. The number of randomly placed sinkers (high viscosity inclusions in low viscosity medium)
increases along the horizontal axis. The vertical axis depicts the number of GMRES iterations required
for 106 residual reduction for the most popular Mp (1/µ) and the proposed w-BFBT preconditioner. Fixed
problem parameters are the dynamic ratio DR(µ) = max(µ)/min(µ) = 108 , discretization order k = 2, and
the mesh refinement level ` = 7, resulting in 1283 finite elements. Right image shows an example viscosity
field with 16 sinkers (blue spheres depict highly viscous centers of Gaussian-like sinkers, low viscosity medium
in red color ) and the streamlines of the computed velocity field.
viscosities because convergence becomes extremely slow or stagnates, as observed in [31]. Instead, in [31]
for finite element and in [18] for staggered grid finite difference discretizations, a re-scaling of the discrete
Stokes system (1) was performed, which essentially alters the diagonal weighting matrices C, D. By choosing
entries from A for these weighting matrices, it was possible to demonstrate improved convergence with BFBT
compared to Mp (1/µ) for certain benchmark problems with strong viscosity variations. Building on ideas
from [31], [36] chose the weighting matrices such that C = D = diag(A), which led to superior performance
compared to Mp (1/µ) for highly heterogeneous mantle convection problems. Hence we refer to this approach
as diag(A)-BFBT.
However, even diag(A)-BFBT can fail to achieve fast convergence for some problems/discretizations,
as shown below (Section 2.2). Moreover, choosing the weighting matrices as diag(A) is problematic for
high-order discretizations, where diag(A) becomes a poor approximation of A. These drawbacks lead us to
propose the following w-BFBT approximation for the inverse Schur complement:



−1 T −1
−1 T
T −1
S̃−1
BC−1
BD−1
,
(8)
wl ADwr B
wr B
w-BFBT := BCwl B
where Cwl = M̃u (wl ) and Dwr = M̃u (wr ) are lumped velocity space
p mass matrices (lumping analogously to
(4)) that are weighted by the square root of the viscosity, wl (x) = µ(x) = wr (x), x ∈ Ω.

1.4

Outline and summary of key results

After defining a class of benchmark problems (Section 2.1), we compare the convergence obtained with different
Schur complement approximations to motivate preconditioning with w-BFBT (Section 2.2). Theoretical
estimates for spectral equivalence of w-BFBT are derived in Section 3. This is followed by a detailed
numerical study showing when w-BFBT is advantageous over Mp (1/µ) (Section 4), and a discussion of
boundary modifications for w-BFBT that accelerate convergence (Section 5). In Section 6 we describe
an algorithm for w-BFBT-based Stokes preconditioning, which uses hybrid spectral-geometric-algebraic
multigrid (HMG). Finally, in Section 7, we provide numerical evidence for near-optimal algorithmic and
parallel scalability. In particular, we demonstrate that the preconditioner’s parallel efficiency remains high
when weak scaling out to tens of thousands of threads and even millions of threads.3
3 See [36], which demonstrated parallel scalability for a BFBT-type method. w-BFBT deviates from [36] because of a different
choice of diagonal weighting matrices in (7). However, the work per application of both preconditioners is the same, resulting in
comparable parallel scalability.

w-BFBT Preconditioner for Stokes Flow Problems

5

To motivate our study of w-BFBT, we give an example for a possible improvement in convergence in
Figure 1. There, a comparison is drawn between the Mp (1/µ) and w-BFBT approximations for the Schur
complement. The Stokes problem that is being solved is the multi-sinker benchmark problem from Section 2.1.
The difficulty of the problem can be increased by adding more and more high-viscosity inclusions, called
sinkers, into a low-viscosity background medium, which, as a result, introduces more variation in the viscosity.
As can be seen in the figure, the number of GMRES iterations remain flat when preconditioning with w-BFBT,
whereas the number of GMRES iterations increases significantly for higher sinker counts when using Mp (1/µ),
rendering Mp (1/µ) inefficient for these types of difficult problems. Therefore we propose w-BFBT as an
alternative Schur complement approximation for Stokes flow problems with a highly varying viscosity.

2

Benchmark problem and comparison of Schur complement approximations

This section further motivates the need for more effective Schur complement preconditioners. We first present
a class of benchmark problems that range from relatively mild viscosity variations to highly heterogeneous.
Then a challenging problem is used to compare Stokes solver convergence with different Schur complement
approximations to demonstrate the limitations of established methods and motivate the development of
w-BFBT.

2.1

Multi-sinker benchmark problem

The design of suitable benchmark problems is critical to conduct studies that can give useful convergence
estimates for challenging applications. We seek complex geometrical structures in the viscosity that generate
irregular, nonlocal, multiscale flow fields. Additionally, the viscosity should exhibit sharp gradients and its
dynamic ratio DR(µ) := max(µ)/min(µ) (also commonly referred to as viscosity contrast) can be six orders of
magnitude or higher in demanding applications. As in [29], we use a multi-sinker test problem with randomly
positioned inclusions (e.g., as in Figure 1, right image) to study solver performance. We find that the arising
viscosity structure is a suitable test for challenging, highly heterogeneous coefficient Stokes problems, and that
the solver performance observed for such models can be indicative of the performance for other challenging
applications.
In the (open) unit cube domain Ω = (0, 1)3 , we define the viscosity coefficient µ(x) ∈ [µmin , µmax ], x ∈ Ω,
0 < µmin < µmax < ∞, with dynamic ratio DR(µ) = µmax /µmin by means of rescaling a C ∞ indicator
function χn (x) ∈ [0, 1] that accumulates n sinkers via a product of modified Gaussian functions:
µ(x) := (µmax − µmin )(1 − χn (x)) + µmin , x ∈ Ω,


n

Y
ω 2
χn (x) :=
1 − exp −δ max 0, |ci − x| −
, x ∈ Ω,
2
i=1
where ci ∈ Ω, i = 1, . . . , n, are the centers of the sinkers, δ > 0 controls the exponential decay of the Gaussian
smoothing, and ω ≥ 0 is the diameter of a sinker where µmax is attained. Since all sinkers are equal in size,
inserting more of them inside the domain will eventually result in overlapping with each other and possible
intersections with the domain’s boundary. Throughout the paper, we fix δ = 200, ω = 0.1, and use the
same set of precomputed random points ci in all numerical experiments. Two parameters are varied: (i)
the number of sinkers n at random positions (the label Sn-rand indicates a multi-sinker problem with n
randomly positioned sinkers) and (ii) the dynamic ratio DR(µ) which in turn determines µmin := DR(µ)−1/2
and µmax := DR(µ)1/2 . The right-hand side of (2), f (x) := (0, 0, β(χn (x) − 1)), β = 10 constant, is such
that it forces the high-viscosity sinkers downward, similarly to a gravity that pulls on high-density inclusions
within a medium of lower density.

2.2

Comparison of Schur complement approximations

We compare convergence of the Stokes solver using the Schur complement approximation Mp (1/µ) with
diag(A)-BFBT and with the proposed w-BFBT. The problem parameters are held fixed to S16-rand and

w-BFBT Preconditioner for Stokes Flow Problems

kresk / kinit resk
kresk / kinit resk

fixed level ` = 5

fixed order k = 2

Mp (1/µ)

6

diag(A)-BFBT

w-BFBT

100

`=5
`=6
`=7

10−2
10−4
10−6

0

25

50

75 100 0

25

50

75 100 0

25

50

75 100

100

k
k
k
k

10−2
10−4
10−6

0

25

50

75 100 0

GMRES iteration

25

50

75 100 0

GMRES iteration

25

50

=2
=3
=4
=5

75 100

GMRES iteration

Figure 2: Comparison of Stokes solver convergence with Mp (1/µ) (left column), diag(A)-BFBT (middle
column), and w-BFBT (right column) preconditioning. We fix the problem S16-rand, DR(µ) = 108 while
varying mesh refinement level ` (top row) and discretization order k (bottom row). This comparison shows
that w-BFBT combines robust convergence of diag(A)-BFBT with improved algorithmic scalability when k
is increased.
DR(µ) = 108 . The numerical experiments are carried out using different levels of mesh refinement ` = 5, . . . , 7
(for fixed order k = 2) and different discretization orders k = 2, . . . , 5 (for fixed level ` = 5). A level `
corresponds to a mesh of 23` elements due to uniform refinement. Note that for these tests, the applications of


T −1
T −1
A−1 , BC−1
, and BD−1
are approximated using a multigrid method (introduced in Section 6).
wl B
wr B
These approximations are sufficiently accurate, such that the comparison is indicative of the effectiveness of
the different Schur complement approximations. In particular, improving the approximation does not change
the results, which are presented in Figure 2. In the left two plots, the poor Schur complement approximation
by Mp (1/µ) for this problem setup can be observed clearly. Convergence stagnates (similar results are found
in [29, 36]).
Preconditioner diag(A)-BFBT (Figure 2, middle) is able to achieve fast convergence for discretization
order k = 2. A limitation of diag(A)-BFBT is a strong dependence on the order k. This can be explained
by the decreasing diagonal dominance in the viscous block A with increasing order k: for higher k the
approximation of A by diag(A) deteriorates. Note that numerical experiments with Mu -BFBT are not
presented, because it performs poorly in the presence of spatially-varying viscosities. This leads to the
conclusion that the choice of the weighting matrices C, D in S̃−1
BFBT crucially affects the quality of the Schur
complement approximation.
The w-BFBT approximation delivers convergence that is nearly as fast as in the diag(A)-BFBT, k = 2
case, but without the severe deterioration when k is increased (see Figure 2, right). Thus, w-BFBT exhibits the
robustness of diag(A)-BFBT and additionally shows superior algorithmic scalability with respect to k. Having
illustrated the efficacy of w-BFBT for certain problem parameters, we next establish spectral equivalence
of w-BFBT (Section 3) and then analyze in more detail how crucial parameters influence convergence in
Section 4.

3

Spectral equivalence of w-BFBT

Before we show spectral equivalence, we introduce notation and basic definitions.

3.1

Basic definitions

Let Ω ⊂ Rd be a bounded domain with Lipschitz boundary ∂Ω. We denote as L2 (Ω) the class of realvalued square-integrable functions, equipped with the usual L2 -inner product (u , v)L2 (Ω) and induced norm

w-BFBT Preconditioner for Stokes Flow Problems

7

kukL2 (Ω) , u, v ∈ L2 (Ω). We also consider corresponding spaces of d-dimensional vector-valued functions
and (d × d)-dimensional tensor-valued functions with component-wise multiplication, denoted by L2 (Ω)d
and L2 (Ω)d×d . The subspace of L2 (Ω) that does not contain constant functions is denoted by L2 (Ω)/R.
A bounded function, say α = α(x), belongs to the space L∞ (Ω) by satisfying the following finite norm:
kαkL∞ (Ω) := ess supx∈Ω |α(x)| < ∞. We generalize the L2 -norms to classes of weighted L2α -norms for
functions f ∈ L2 (Ω)n , n ∈ {1, d, d × d}, defined by
1

kf kL2α (Ω)n := α 2 f

for α ∈ L∞ (Ω), 0 < α(x) a.e. in Ω.

L2 (Ω)n

Next, we introduce H m (Ω), with m ≥ 0, which is the Sobolev space of m derivatives in L2 (Ω), and for
m = 1 we use the inner product (u , v)H 1 (Ω) := (u , v)L2 (Ω) + (∇u , ∇v)L2 (Ω) , inducing the norm kukH 1 (Ω) .
Functions in H m (Ω) with vanishing trace on the boundary ∂Ω belong to the space H0m (Ω). Finally, we
say that a function belongs to the class of C ∞ (Ω) if it has partial derivatives of any order in Ω, and these
derivatives are continuous.
We transition from abstract definitions to fluid mechanics. The differential operators acting on velocity
u ∈ (H 1 (Ω))d and pressure p ∈ L2 (Ω) within the Stokes equations are defined in the sense of distributions:
∇s u :=

1
(∇u + ∇uT ),
2

Bu := −∇ · u,

B ∗ p := ∇p.

Moreover, assume a sufficiently regular, bounded viscosity µ ∈ H 1 (Ω) ∩ L∞ (Ω) such that 0 < µmin ≤ µ(x)
a.e. in Ω and then define the viscous stress tensor τ := 2µ∇s u. We denote the function space for velocity by

V := u ∈ (H 1 (Ω))d | n · u = 0 on ∂Ω ,
(9)
where n ∈ Rd is the outward unit normal vector at the boundary ∂Ω, and the function space for pressure by
Q := L2 (Ω)/R, and we introduce the viscous stress operator with a heterogeneous viscosity
Aµ : V → V 0 ,

Aµ u := −∇ · (2µ∇s u) = −∇ · τ .

Given exterior forces acting on the fluid f ∈ V 0 , we consider the incompressible Stokes problem with free-slip
and no-normal flow boundary conditions:
Aµ u + B ∗ p = f

in Ω,

(10a)

Bu = 0

in Ω,

(10b)

T [τ − pI] n = 0

on ∂Ω,

(10c)

u·n=0

on ∂Ω,

(10d)

in which we seek the velocity u ∈ V and pressure p ∈ Q. On the boundary, we have outward unit normal
vectors n ∈ Rd and tangential projectors T := I − nnT . For the theoretical analysis sections, we altered the
boundary conditions compared to (1c) to achieve a cleaner presentation.
For the definition of the w-BFBT approximation of the Schur complement, we introduce a Poisson operator
for higher regularity pressure functions
∗
Kw
: H 2 (Ω) → L2 (Ω),

∗
Kw
p := BwB ∗ p,

(11)

with an appropriate coefficient w (see below) and augmented with homogeneous Neumann boundary conditions,
∗
n · B ∗ p = 0. The L2 -adjoint of Kw
is denoted by Kw . Finally, we define the w-BFBT approximation of the
−1 ∗
Schur complement S = BAµ B by:
∗
S̃w-BFBT := Kw
(Bwl Aµ wr B ∗ )−1 Kwl ,
r

(12)

with sufficiently regular, bounded weight functions wl , wr ∈ H 1 (Ω)∩L∞ (Ω) such that 0 < wmin ≤ wl (x), wr (x)
a.e. in Ω. Note that the definitions of the w-BFBT weight functions in the discrete case (8) are reciprocal to
definitions in (11) and (12), because in the discrete case the weight functions were embedded into inverses of
mass matrices.

w-BFBT Preconditioner for Stokes Flow Problems

3.2

8

Main theorem on spectral equivalence of w-BFBT

One measure for the efficacy of a preconditioner consists of the ratio of the maximal to minimal eigenvalues
of the preconditioned system S̃−1
w-BFBT S. This section establishes inequalities for spectral equivalence of
w-BFBT by providing bounds on that ratio. The derivations are carried out in an infinite-dimensional setting.
We begin by stating the main result of this section in Theorem 3.1 and continue with proving this result
using a sequence of lemmas.
Theorem 3.1 (Main result). Let Q̂ = L2 (Ω)/R ∩ H 1 (Ω). If the left and right w-BFBT weight functions are
equal to
1
wl (x) = µ(x)− 2 = wr (x) for a.a. x ∈ Ω,
then the exact Schur complement is equivalent to the w-BFBT approximation such that




S̃w-BFBT q , q ≤ (Sq , q) ≤ Cw-BFBT S̃w-BFBT q , q
for all q ∈ Q̂,
where




1
2
2
2
Cw-BFBT := 1 + k∇µkL∞ (Ω)d CP,µ
+ 1 CK,µ
4
and the constants CP,µ , CK,µ > 0 stem from weighted Poincaré–Friedrichs’ and Korn’s inequalities, respectively
(see Remark 3.7 for more information); the viscosity µ assumes the role of the weight function in the weighted
inequalities.
If the viscosity and the w-BFBT weight functions are constant,
µ ≡ 1,

wl ≡ 1 ≡ wr ,

then the exact Schur complement is equivalent to the w-BFBT approximation such that




2
S̃w-BFBT q , q ≤ (Sq , q) ≤ (CP2 + 1) CK
S̃w-BFBT q , q
for all q ∈ Q̂
with constants CP , CK > 0 stemming from (the classical) Poincaré–Friedrichs’ and Korn’s inequalities,
respectively.

3.3

Proofs

The proof of Theorem 3.1 is established in the remainder of this section. In what follows, suprema are
understood over spaces excluding operator kernels that would cause a supremum to blow up. The following
basic, but hereafter frequently used, result is shown for completeness of the discussion.
Lemma 3.2 (sup-form of inverse operator). Let V be a complete Hilbert space and W ⊆ V be a dense
subspace. Assume the linear operator T : V → V 0 to be bounded, invertible, symmetric, and positive definite.
Then for any f ∈ V 0 follows
2

(w , f )
T −1 f , f = sup
.
w∈W (T w , w)
Proof. Let w ∈ W , then with Hölder’s inequality follows

2
2
(w , f ) = T 1/2 w , T −1/2 f ≤ T 1/2 w

2

T −1/2 f

2


= (T w , w) T −1 f , f .

Additionally, let v = T −1 f and since W ⊆ V dense, there exists a sequence {wk }k ⊂ W such that
wk → v = T −1 f ∈ V , hence
2
2
2

T −1 f , f
(wk , f )
(v , f )
→
=
= T −1 f , f ,
−1
(T wk , wk )
(T v , v)
(f , T f )
which shows that equality is achieved in the limit.

w-BFBT Preconditioner for Stokes Flow Problems

9

The next lemma establishes Schur complement properties that are essential for deriving lower and upper
bounds in the spectral equivalence estimates.
Lemma 3.3 (sup-form of Schur complement). With the definitions from Section 3.1, the following two
equalities hold:
2


(B ∗ p , wr B ∗ q)
S̃w-BFBT q , q = sup
for all q ∈ Q̂,
(13)
∗
∗
p∈P̂ (wl Aµ wr B p , B p)
where P̂ := C ∞ (Ω) and Q̂ := L2 (Ω)/R ∩ H 1 (Ω), and
2

(Sq , q) = sup
v∈V

(v , wr B ∗ q)
(wl Aµ wr v , v)

for all q ∈ Q.

(14)

Proof. For q ∈ Q̂, we use integration by parts on the left hand side of (13) to obtain

 Z

S̃w-BFBT q , q =
wr B ∗ (Bwl Aµ wr B ∗ )−1 Kwl q (B ∗ q) dx + b1 (q)
ZΩ

=
(Bwl Aµ wr B ∗ )−1 Kwl q (Kwr q) dx + b1 (q) + b2 (q)
Ω

with boundary terms
Z
b1 (q) := −
Z
b2 (q) :=


n · wr B ∗ (Bwl Aµ wr B ∗ )−1 Kwl q q dx,

∂Ω


(Bwl Aµ wr B ∗ )−1 Kwl q (n · wr B ∗ q) dx.

∂Ω

Using that P̂ ⊂ H 2 (Ω) is dense, application of Lemma 3.2 and further integration by parts yields



S̃w-BFBT q , q = sup
p∈P̂

2

(p , Kwr q)
+ b1 (q) + b2 (q)
(Bwl Aµ wr B ∗ p , p)
2

= sup
p∈P̂

((B ∗ p , wr B ∗ q) + b3 (p, q))
+ b1 (q) + b2 (q)
(wl Aµ wr B ∗ p , B ∗ p) + b4 (p)

with boundary terms
Z
b3 (p, q) := −

p (n · wr B ∗ q) dx,

∂Ω

Z
b4 (p) := −

(n · wl Aµ wr B ∗ p) p dx.

∂Ω
∗
Kw
r

Because the operator
from (11) is augmented with homogeneous Neumann boundary conditions,
∗
n · B p = 0, the boundary terms b1 (q), b2 (q), and b3 (p, q) vanish. In addition, p ∈ P̂ is sufficiently regular
for the term b4 (p) to be well-defined and it equals to zero because the velocity u satisfies n · u = 0 on ∂Ω.
Hence, (13) follows.
To show (14), let q ∈ Q, then for the exact Schur complement we apply integration by parts with a
vanishing boundary term



∗
−1 ∗
∗
−1 −1 −1
∗
∗
(Sq , q) = BA−1
µ B q , q = Aµ B q , B q = wr Aµ wl wl B q , wr B q
and (14) follows from Lemma 3.2.
A direct consequence of Lemma 3.3 is the following lower bound.

w-BFBT Preconditioner for Stokes Flow Problems

10

Corollary 3.4 (Lower bound, S̃w-BFBT . S). The exact Schur complement is bounded by the w-BFBT
approximation from below, i.e.,


S̃w-BFBT q , q ≤ (Sq , q) for all q ∈ Q̂,
(15)
where Q̂ := L2 (Ω)/R ∩ H 1 (Ω).
Proof. Let P̂ = C ∞ (Ω) and q ∈ Q̂. Since B ∗ maps P̂ into V , we combine (13) and (14) to get



S̃w-BFBT q , q = sup
p∈P̂

2

2

(B ∗ p , wr B ∗ q)
(v , wr B ∗ q)
≤ sup
= (Sq , q) ,
∗
∗
(wl Aµ wr B p , B p) v∈V (wl Aµ wr v , v)

and obtain the result (15).
We begin the derivation of an upper bound for the case of constant viscosity µ ≡ 1. Note that S̃w-BFBT
is scaling invariant with respect to constants multiplied to the w-BFBT weight functions wl , wr . Hence,
it always assumes the correct scaling of S independent of the viscosity constant. The result for constant
viscosity presented below in Lemma 3.5 is generalized to variable viscosity in Lemma 3.6. While Lemma 3.5 is
a special case of Lemma 3.6, we first prove the result for constant viscosity as the arguments are less technical
and easier to follow. In the proof of the result for variable viscosity, we build on some of the arguments from
the constant viscosity case and thus avoid unnecessary duplication.
Lemma 3.5 (Upper bound, S . S̃w-BFBT , for constant µ). Assume a constant viscosity µ ≡ 1 and constant
w-BFBT weight functions wl ≡ 1 ≡ wr , and, as before in Lemma 3.3, let Q̂ = L2 (Ω)/R ∩ H 1 (Ω). Then the
exact Schur complement is bounded by the w-BFBT approximation from above by


2
(Sq , q) ≤ (CP2 + 1) CK
S̃w-BFBT q , q
for all q ∈ Q̂
(16)
with constants CP , CK > 0 stemming from Poincaré–Friedrichs’ and Korn’s inequalities, respectively.
Proof. Let P̂ = C ∞ (Ω) and q ∈ Q̂, then due to (13) we can write


2
∗
2

(B ∗ p , B ∗ q) kB pk(H 1 (Ω))d
S̃w-BFBT q , q = sup
.
2
∗
∗
∗
p∈P̂ kB pk(H 1 (Ω))d (A1 B p , B p)

(17)

To estimate the second factor on the right-hand side of (17), note that
2

(A1 B ∗ p , B ∗ p) = 2 (∇s B ∗ p , ∇s B ∗ p) = 2 (∇B ∗ p , ∇B ∗ p) = 2 k∇B ∗ pk(L2 (Ω))d×d ,
where we used that n · B ∗ p = 0 on the boundary and that ∇B ∗ p is symmetric. Thus,
2

(A1 B ∗ p , B ∗ p) ≤ 2 kB ∗ pk(H 1 (Ω))d .

(18)

For the first factor on the right-hand side of (17), observe that for any v ∈ V there exists a sequence
{pi }i ⊂ P̂ ⊂ H 2 (Ω) such that K1∗ pi = BB ∗ pi → Bv, since K1∗ is invertible, where convergence is with respect
to the L2 -norm. Thus,
sup
p∈P̂

(B ∗ p , B ∗ q)
(v , B ∗ q)
=
sup
= kB ∗ qk(H −1 (Ω))d .
kB ∗ pk(H 1 (Ω))d
v∈V kvk(H 1 (Ω))d

(19)

Combining (18) and (19) provides the following estimate for the w-BFBT Schur complement approximation (17):

 1
2
S̃w-BFBT q , q ≥ kB ∗ qk(H −1 (Ω))d .
(20)
2
The exact Schur complement, on the other hand, in the form (14) from Lemma 3.3, can be bounded by
2

(Sq , q) = sup
v∈V

2

2
kvk(H 1 (Ω))d kB ∗ qk(H −1 (Ω))d
(v , B ∗ q)
≤ sup
.
2
(A1 v , v)
v∈V
2 k∇s vk(L2 (Ω))d×d

(21)

w-BFBT Preconditioner for Stokes Flow Problems

11

With Poincaré–Friedrichs’ inequality,
kvkL2 (Ω)d ≤ CP k∇vk(L2 (Ω))d×d ,
where the constant CP > 0 depends on the domain Ω, and Korn’s inequality,
k∇vk(L2 (Ω))d×d ≤ CK k∇s vk(L2 (Ω))d×d ,
with constant CK > 0, we obtain
1
2
2
2
2
kvk(H 1 (Ω))d ≤ k∇vk(L2 (Ω))d×d ≤ CK
k∇s vk(L2 (Ω))d×d ,
(CP2 + 1)
and substituting this into (21) gives
(Sq , q) ≤

2
(CP2 + 1) CK
2
kB ∗ qk(H −1 (Ω))d .
2

(22)

Together with (20), this yields the desired result (16).
We complete the presentation of spectral equivalence by deriving an upper bound for problems with
variable viscosities.
Lemma 3.6 (Upper bound, S . S̃w-BFBT , for variable µ). As before in Lemma 3.3, let Q̂ = L2 (Ω)/R∩H 1 (Ω).
If the left and right w-BFBT weight functions are equal to
1

wl (x) = µ(x)− 2 = wr (x)

for a.a. x ∈ Ω,

then the exact Schur complement is bounded by the w-BFBT approximation from above by





1
2
2
2
+ 1 CK,µ
S̃w-BFBT q , q
for all q ∈ Q̂
(Sq , q) ≤ 1 + k∇µkL∞ (Ω)d CP,µ
4

(23)

(24)

with constants CP,µ , CK,µ > 0 stemming from weighted Poincaré–Friedrichs’ and Korn’s inequalities, respectively, where µ assumes the role of the weight function.
Proof. Let the weight functions be equal, wl ≡ w ≡ wr , but (for now) otherwise arbitrary subject to the
condition 0 < wmin ≤ w(x) for a.a. x ∈ Ω. At the end of this proof we will argue the special role of the
choice (23) for the weight functions. In addition, let P̂ = C ∞ (Ω) and q ∈ Q̂, then due to (13) we can write


2
2

kB ∗ pk(H 1 (Ω))d
(B ∗ p , wB ∗ q)
S̃w-BFBT q , q = sup
.
2
∗
∗
∗
p∈P̂ kB pk(H 1 (Ω))d (wAµ wB p , B p)

(25)

We begin by estimating the second factor on the right-hand side of (25). For an arbitrary v ∈ (H 1 (Ω))d ,
observe that ∇s (wv) = w∇s v + ∇w ⊗ v, where “⊗” denotes the outer product of two vectors in Rd , and thus
√
√
2
(wAµ wv , v) = 2 (µ∇s (wv) , ∇s (wv)) = 2 k µw∇s v + µ∇w ⊗ vk(L2 (Ω))d×d .
Applying the triangle inequality and then Hölder’s inequality to the resulting terms,
√
√
2
2
2
k µw∇s vk(L2 (Ω))d×d ≤ k µwkL∞ (Ω) k∇s vk(L2 (Ω))d×d ,
√
√
2
2
2
k µ∇w ⊗ vk(L2 (Ω))d×d ≤ k µ∇wkL∞ (Ω)d kvkL2 (Ω)d ,
and thus we obtain the estimate
2

(wAµ wv , v) ≤ 2CA,µ,w kvk(H 1 (Ω))d ,
where

√
√
2
2
CA,µ,w := k µwkL∞ (Ω) + k µ∇wkL∞ (Ω)d .

(26)

w-BFBT Preconditioner for Stokes Flow Problems

12

Similarly to (19) and (20), we obtain the following estimate for the w-BFBT Schur complement approximation (25):


1
2
S̃w-BFBT q , q ≥
kwB ∗ qk(H −1 (Ω))d .
(27)
2CA,µ,w
Proceeding with the exact Schur complement, we obtain from (14) in Lemma 3.3 that
2
2
2
w−1 v (H 1 (Ω))d kwB ∗ qk(H −1 (Ω))d
w−1 v , wB ∗ q
≤ sup
.
(28)
(Sq , q) = sup
√
2
(Aµ v , v)
v∈V
v∈V
2
µ∇s v 2
d×d
(L (Ω))

We require a weighted Poincaré–Friedrichs’ inequality (see Remark 3.7 for details),
kvkL2

w−2

(Ω)d

≤ CP,w−2 k∇vkL2

w−2

(Ω)d×d

,

(29)

and also a weighted Korn’s inequality (see Remark 3.7 for more information),
k∇vkL2 (Ω)d×d ≤ CK,µ k∇s vkL2 (Ω)d×d .
µ

µ

With (29) and (30), we are able to bound (28) from above:


!
2
2
k∇vkL2 (Ω)d×d
CP,w
−2 + 1 CK,µ
2
w−2
sup
kwB ∗ qk(H −1 (Ω))d .
(Sq , q) ≤
2
v∈V k∇vk(L2 (Ω))d×d

(30)

(31)

µ

The supremum term in (31) and the constant CA,µ,w in (26) motivate the choice for the weight w to be
1

w := µ− 2 .
Then the supremum in (31) vanishes and (26) simplifies to
1
2
k∇µkL∞ (Ω)d .
4
Substituting this into (27) together with inequality (31) yields the desired result (24).
CA,µ,w = 1 +

Remark 3.7. In the proof of Lemma 3.6 we utilized a weighted Poincaré–Friedrichs’ inequality, for which the
optimal constant is
kvkL2µ (Ω)d
,
CP,µ = sup
v∈V k∇vkL2 (Ω)d×d
µ

where the viscosity takes the role of the weight function. While weighted Poincaré and Friedrichs’ inequalities
have been investigated in the literature numerous times, usually they are proven by contradiction and scaling
arguments, which does not provide information about the constants. If explicit constants are found, they
depend, in general, on the weight such that the resulting estimates are too pessimistic, e.g., CP,µ = O(DR(µ)).
Knowledge of constants that are robust with respect to weight functions is limited. In the context of a
posteriori error estimates for finite elements, weight-independent constants could be found for convex domains
and weights that are a positive power of a non-negative concave function [10]. These results were refined for
star-shaped domains under certain assumptions for the weights [42]. For another class of weights, namely
quasi-monotone piecewise constant weight functions, robust constants were derived in [33].
In addition to weighted Poincaré–Friedrichs’, we utilized a weighted Korn’s inequality in the proof of
Lemma 3.6. The optimal constant for this inequality is
CK,µ = sup
v∈V

k∇vkL2µ (Ω)d×d
k∇s vkL2µ (Ω)d×d

.

As for CP,µ , straightforward estimation results in an overly pessimistic weight-dependent constant, namely
CK,µ = O(DR(µ)), [25]. Other work utilizing weighted Korn’s inequalities usually aims to derive inequalities
for special domain shapes, e.g., [1].
In summary, robust constants for weighted Poincaré–Friedrichs’ and Korn’s inequalities for general weight
functions are difficult to obtain and limitations exist in the form of assumptions on the weights. Further
work on this topic could improve the constants for the spectral equivalence of w-BFBT but is beyond the
scope of this paper.

w-BFBT Preconditioner for Stokes Flow Problems

13

Table 1: Robustness classification for Schur complement approximations a Mp (1/µ) and b w-BFBT in terms
of number of GMRES iterations (10−6 residual reduction, GMRES restart every 100 iterations). Number of
randomly placed sinkers (#sinkers) is increased across rows, while dynamic ratio (DR(µ)) is increased across
columns. Discretization is fixed at k = 2, ` = 7.
(a) Mp (1/µ)
#sinkers \ DR(µ)
S1-rand
S4-rand
S8-rand
S12-rand
S16-rand
S20-rand
S24-rand
S28-rand

4

(b) w-BFBT

104

106

108

1010

29
53
64
70
85
84
117
108

31
63
79
86
167
167
286
499

31
71
93
99
231
380
3279
2472

29
80
165
180
891
724
5983
>10000

#sinkers \ DR(µ)
S1-rand
S4-rand
S8-rand
S12-rand
S16-rand
S20-rand
S24-rand
S28-rand

104

106

108

1010

29
39
38
38
40
34
31
29

29
41
40
40
45
36
32
31

29
42
41
43
47
37
39
42

30
44
44
45
48
38
55
60

Robustness of w-BFBT

In this section, we analyze the robustness properties of the widely used Schur complement approximation
Mp (1/µ) and the new w-BFBT via numerical experiments. Furthermore, we calculate the spectra for both
approaches and thus support the discussion in Section 3 about theoretical eigenvalue bounds with numerical
results. The comparison of Mp (1/µ) and w-BFBT is of particular importance, because of the widespread
use of the inverse viscosity-weighted mass matrix. It is therefore of interest to determine when convergence
with Mp (1/µ) deteriorates and using w-BFBT becomes beneficial. A comparison with diag(A)-BFBT was
not performed because Section 2.2 already showed that diag(A)-BFBT performs similarly or worse than
w-BFBT, hence there are no advantages in using diag(A)-BFBT over w-BFBT.
For the numerical experiments in this section, we return to the definitions and setup from Section 2.1. To
apply the inverse of Mp (1/µ), we diagonalize the mass matrix of the discontinuous, modal pressure space by
forming its lumped version (4). Moreover, to apply the approximate inverse of the viscous block in (3) we
use the same multigrid method for each of the two Schur approximations; this multigrid method is also used
for the inverse operators of w-BFBT in (8). The details of the multigrid method are provided in Section 6.
To compare the robustness, we vary two problem parameters: (i) the number of randomly placed sinkers n
and (ii) the dynamic ratio DR(µ). The parameter n influences the geometric complexity of the viscosity µ
while DR(µ) controls the magnitude of viscosity gradients.
Tables 1a and 1b present the number of GMRES iterations for a 10−6 residual reduction in the Euclidean
norm. Observe that for the S1-rand problem, the iteration count is essentially the same for both Mp (1/µ)
and w-BFBT, and that it stays stable across all dynamic ratios DR(µ) = 104 , . . . , 1010 . Hence for this
simple problem, w-BFBT has no advantages and its additional computational cost makes it less efficient than
Mp (1/µ). However, the limitations of the Mp (1/µ) approach become apparent by increasing the number
of randomly positioned sinkers. Two observations for Mp (1/µ) can be made from Table 1a. First, the
number of GMRES iterations rises with increasing number of sinkers (factor ∼80 increase for n = 1, . . . , 28,
DR(µ) = 108 ). Second, in a multi-sinker setup the dependence on DR(µ) becomes more severe (factor ∼50
increase for n = 24, DR(µ) = 104 , . . . , 1010 ). This demonstrates that Mp (1/µ) is a poor approximation of
the Schur complement for certain classes of problems with highly heterogeneous viscosities.
The advantages in robustness of the w-BFBT preconditioner are demonstrated in Table 1b. Compared
to Mp (1/µ), the number of GMRES iterations is stable and the increase over the whole range of problem
parameters is just a factor of 2. Only 60 iterations are needed for the most extreme problem, namely S28-rand,
DR(µ) = 1010 , for which convergence with Mp (1/µ) essentially stagnated.
More insight concerning the different convergence behaviors can be gained from the eigenvalues in Figure 3.
The plots in that figure are for two-dimensional multi-sinker problems, which are analogous to the threedimensional benchmark problems from Section 2.1. We discretize the problems on triangular meshes utilizing
the FEniCS library [28]. We choose Pbubble
× Pdisc
finite elements [11,12] because they represent a close analog
2
1
disc
to the Q2 ×P1 elements, which are employed on (three-dimensional) hexahedral meshes. Each plot shows the
eigenvalues of the exact Schur complement and the eigenvalues of the preconditioned Schur complement for the

S4-rand, DR(µ) = 104

10−8

10−8

S1-rand, DR(µ) = 106

S̃−1
w-BFBT S

M−1
p (1/µ) S
S

10−6

10−4

10−2

100

102

S̃−1
w-BFBT S

M−1
p (1/µ) S
S

10−6

10−4

14

10−8
S4-rand, DR(µ) = 106

S1-rand, DR(µ) = 104

w-BFBT Preconditioner for Stokes Flow Problems

10−2

Eigenvalue

100

102

10−8

S̃−1
w-BFBT S

M−1
p (1/µ) S
S

10−6

10−4

10−2

100

102

10−2

100

102

S̃−1
w-BFBT S

M−1
p (1/µ) S
S

10−6

10−4

Eigenvalue

Figure 3: Spectra of the Schur complement (gray), Mp (1/µ)-preconditioned Schur complement (red), and
w-BFBT-preconditioned Schur complement (blue); zero eigenvalues corresponding to the null space of the
Schur complement matrix are omitted. Results for viscosities with one sinker (S1-rand) are shown in top row,
and with four sinkers (S4-rand) in the bottom row of plots; DR(µ) = 104 in the left column and DR(µ) = 106
in the right column. The two-dimensional Stokes equations are discretized with Pbubble
× Pdisc
finite elements
2
1
on a uniform triangular mesh consisting of 512 triangles using the FEniCS library. As the problem difficulty
increases from one to four sinkers, the spreading of small eigenvalues for Mp (1/µ) becomes more severe,
which is disadvantageous for solver convergence. For w-BFBT, the spectrum remains largely unaffected by
increased sinker counts, which contributes to convergence that is robust with respect to viscosity variations.
Mp (1/µ) and w-BFBT approximations, where all inverse matrices, e.g., the viscous block matrix A and the
pressure Poisson matrices within w-BFBT, are inverted with a direct solver. Effective preconditioners exhibit
a strong clustering of eigenvalues, whereas the convergence of Krylov methods deteriorates if the eigenvalues
are spread out. We recognize different characteristics in the spectra associated with Mp (1/µ) and w-BFBT
preconditioning. For Mp (1/µ), the dominant eigenvalues are clustered around one while smaller eigenvalues,
i.e., eigenvalues  1, are spread out. The behavior for w-BFBT is the opposite: the dominant eigenvalues
are spread out and the smaller eigenvalues are tightly clustered around one. Now, as the problem difficulty
is increased by introducing more viscosity anomalies in the domain, the spreading of smaller eigenvalues
associated with Mp (1/µ) becomes more severe (compare in Figure 3 the top row of plots and the bottom row
of plots). We postulate that this is the property that is responsible for the deteriorating convergence with
Mp (1/µ) that was observed in Table 1a. With w-BFBT on the other hand, the spectrum remains largely
unaffected by increased sinker counts. The clustering of smaller eigenvalues around one remains stable, which
is likely the reason for the robustness of w-BFBT. The lower bound on the eigenvalues that we observe here
numerically supports the theoretical estimates on spectral equivalence in Section 3. Therefore we find the
lower bound to be sharp and, moreover, to be essential for the robustness of the w-BFBT preconditioner.
Remark 4.1. In addition to the Pbubble
× Pdisc
discretization used for the results in Figure 3, we also calculated
2
1
the spectra using P2 × P1 Taylor-Hood finite elements. We obtained very similar results for this discretization,
which uses continuous elements to approximate the pressure. Therefore, both the efficacy of w-BFBT as a
preconditioner and the issues with Mp (1/µ) seem to be largely unaffected by the specific type of discretization,
at least for the two cases that we tested.
Remark 4.2. The convergence of the Stokes solver with the Mp (1/µ) preconditioner can be improved by

w-BFBT Preconditioner for Stokes Flow Problems

15

approximating the heterogeneous viscosity µ(x) with elementwise constants, computed by averaging µ over
each element [3]. The benefit of faster convergence comes at the cost of slower asymptotic convergence
of the discrete finite element solution and an altered constitutive relationship, which might be a less
accurate representation of the physics; this is, however, problem-dependent. Moreover, for a nonlinear (e.g.,
power-law) rheology, elementwise averaging of µ can introduce non-physical, artificial disturbances in the
effective viscosity during Newton or Picard-type nonlinear solves. We observed such a behavior in mantle
convection simulations, which are governed by a nonlinear power-law rheology. Here, viscosity averaging led
to non-physical checkerboard-like patterns upon convergence of the nonlinear Newton solver.
Remark 4.3. In practice, the convergence of w-BFBT can be improved for coarse meshes, where the viscosity
variations over elements are large. This is achieved by alternative choices for the diagonal weighting matrices
Cwl and Dwr from (8) with the weights

1
2 4
= wr (x) for all x ∈ Ω,
wl (x) = µ2 (x) + |∇µ(x)|
where |·| denotes the Euclidean norm in Rd . These viscosity gradient-based w-BFBT weights have the
advantage of performing at least as well as the pure viscosity-based weights proposed in Section 1.3, but they
exhibit superior robustness on coarser meshes. They are, however, challenging to analyze theoretically.

5

Modifications for Dirichlet boundary conditions

In Section 2.2, deteriorating approximation properties of w-BFBT for increasing discretization order and mesh
refinement level could be observed. The numerical experiments in Figure 2, right did show slightly slower
convergence when k and ` were increased. This can stem from w-BFBT representing a poor approximation
to the exact Schur complement at the boundary ∂Ω in the presence of Dirichlet boundary conditions for the
velocity. This section investigates modifications to w-BFBT near a Dirichlet boundary and aims at obtaining
mesh independence and only a mild dependence on discretization order in terms of Stokes solver convergence.
Consider the commutator that leads to the w-BFBT formulation in an infinite-dimensional form: AB ∗ −
∗
B X ≈ 0, where A represents the viscous stress operator, B ∗ the gradient operator and X the sought
commuting operator. In case of an unbounded domain Ω = Rd and constant viscosity µ ≡ 1, this commutator
is exactly satisfied since (∇ · ∇)∇ − ∇(∇ · ∇) = 0. For Dirichlet boundary conditions, the commutator does
not, in general, vanish at the boundary. Therefore a possible source for deteriorating Schur complement
approximation properties of w-BFBT is a commutator mismatch for mesh elements that are touching the
boundary ∂Ω. A similar observation was also made in [16, 17]. A possible remedy is to modify the norm in
the least-squares minimization problem (6), which is represented by the matrix C−1 , such that a damping
factor is applied to the matrix entries near the boundary. By damping the influence of the boundary in the
minimization objective, more emphasis is given to the domain interior, and the w-BFBT approximation is
improved.
−1
Damping near Dirichlet boundaries can be incorporated by modifying the matrices C−1
wl or Dwr of the
w-BFBT inverse Schur complement approximation (8). A similar idea for BFBT in a Navier–Stokes setting
is presented in [17], where a damping to the weighting matrix D−1 in (7) is introduced to achieve mesh
independence (C−1 is not changed). There, damping affects the normal components of the velocity space
inside mesh elements touching ∂Ω and simply a constant damping factor of 1/10 is set regardless of mesh
refinement `. Also, only the discretization order k = 2 was considered (in addition to Q2 × Q1 and MAC
discretizations).
Now, we attempt to enhance our understanding of how modifications at a Dirichlet boundary ∂Ω S
influence
convergence and therefore the efficacy of w-BFBT as a Schur complement approximation. Let ΩD := e∈D Ωe ,
D := {e | Ωe ∩ ∂Ω 6= ∅} be the set of all mesh elements Ωe touching the Dirichlet boundary. Given values
√
al , ar ≥ 1, extend the previous definition of the weights wl = wr = µ (see Section 1.3) to a version with
boundary modification:
( p
( p
al µ(x) x ∈ ΩD ,
ar µ(x) x ∈ ΩD ,
wl (x) := p
and wr (x) := p
(32)
µ(x)
x∈
/ ΩD ,
µ(x)
x∈
/ ΩD .

w-BFBT Preconditioner for Stokes Flow Problems

16

Table 2: Influence of boundary modification factors al , ar on the Stokes solver convergence with w-BFBT
for discretizations: k = 2, ` = 5, . . . , 7 (see a, c, e) and k = 2, . . . , 5, ` = 5 (see a, b, d, f). Reported are the
number of GMRES iterations for 10−6 residual reduction for the problem S16-rand, DR(µ) = 106 . Colors
highlight solves within ∼5% of iterations above the lowest iteration count. Increase of mesh refinement level
` or discretization order k demands larger boundary amplification ar or al to maintain fast convergence.
(b) k = 3, ` = 5

(a) k = 2, ` = 5
al \ ar
1
2
4
8
16
32

1

2

4

8

16

32

33
33
33
34
34
34

33
33
34
34
34
34

34
34
34
36
38
39

34
34
36
39
43
44

34
34
38
43
46
49

35
34
39
44
49
53

al \ ar

1
2
4
8
16
32

(c) k = 2, ` = 6
al \ ar

1
2
4
8
16
32

1
2
4
8
16
32

2

4

8

16

32

38
37
38
38
39
39

37
38
40
42
44
45

37
38
42
47
50
51

37
39
44
50
53
56

37
39
46
51
56
59

(d) k = 4, ` = 5

1

2

4

8

16

32

37
34
33
34
34
34

34
34
33
34
34
34

33
34
34
35
36
37

34
34
35
38
39
39

34
34
36
39
40
41

34
34
37
39
41
42

al \ ar

1
2
4
8
16
32

(e) k = 2, ` = 7
al \ ar

1
41
38
37
36
37
37

1

2

4

8

16

32

44
39
36
36
36
36

39
39
39
40
41
41

36
39
43
47
49
50

36
40
47
52
56
58

36
41
49
56
60
63

36
41
51
58
63
66

(f) k = 5, ` = 5

1

2

4

8

16

32

45
37
34
34
34
34

37
36
36
36
36
36

34
35
38
39
40
41

34
36
39
42
44
44

34
36
40
44
45
46

34
36
41
44
46
47

al \ ar

1
2
4
8
16
32

1

2

4

8

16

32

63
53
47
44
43
44

53
51
51
51
52
53

46
51
55
59
62
64

43
51
59
65
69
72

43
52
62
69
75
78

44
53
64
72
78
82

We obtain matrices Cwl = M̃u (wl ) and Dwr = M̃u (wr ) in (8) that may differ at boundary elements in
ΩD due to possibly different values for al and ar . Note that amplifying the weight functions wl , wr at the
−1
boundary is similar to damping at the boundary after taking the inverses C−1
wl , Dwr .
The Stokes solver convergence under the influence of boundary amplifications al , ar is summarized in
Table 2. The table shows that the boundary amplification is most effective when performed non-symmetrically,
i.e., either al > 1 or ar > 1 but not both. Further, we deduce that with higher mesh refinement level `, the
boundary amplification should increase roughly proportional to 2` (or proportional to the reciprocal element
size, here h−1 = 2` ). Similar observations can be made for the discretization order k, i.e., amplification needs
to increase for larger k to avoid higher iteration counts. These implications were made based on extensive
numerical experiments for which Table 2 serves as a representative summary.
Remark 5.1. The theoretical derivations of spectral equivalence from Section 3 and the necessity for damping
at Dirichlet boundaries appear inconsistent. However, spectral equivalence was shown in infinite dimensions
whereas boundary damping is applied to the discretized problem. Therefore, we believe that the necessity for
damping is introduced through the discretization. It still remains an open question what might be causing
the slowdown in convergence that is avoided by damping.

6

Parallel hybrid spectral-geometric-algebraic multigrid (HMG) for
w-BFBT

Two aspects of the Stokes preconditioner with w-BFBT have not been discussed yet. One is the approximation
−1
of the inverse viscous block Ã−1 required in (3) and the other is the approximation of inverses K̃−1
wl ≈ Kwl :=

w-BFBT Preconditioner for Stokes Flow Problems

17

−1
T −1
−1 T −1
(BC−1
and K̃−1
in (8). These approximations are crucial for overall Stokes
wr ≈ Kwr := (BDwr B )
wl B )
solver performance and scalability and are addressed in this section. For brevity, we limit our discussion to
−1
K̃−1
wr ; the results also hold for K̃wl .
The approximation of the inverse viscous block Ã−1 is well suited for multigrid V-cycles. To this end,
in [36], we developed a hybrid spectral-geometric-algebraic multigrid (HMG) method, which exhibits extreme
parallel scalability and retains nearly optimal algorithmic scalability (see Section 7 for scalability results).
While traversing the HMG hierarchy shown in Figure 4, HMG initially reduces the discretization order (spectral
multigrid); after arriving at order one, it continues by coarsening mesh elements (geometric multigrid); once
the degrees of freedom fall below a threshold, algebraic multigrid (AMG) carries out further coarsening until a
direct solve can be computed efficiently ( [41]). Parallel forest-of-octrees algorithms, implemented in the p4est
parallel adaptive mesh refinement library [5, 9], are used for efficient, scalable mesh refinement/coarsening,
mesh balancing, and repartitioning in the geometric HMG phase. During parallel geometric coarsening,
the number of compute cores and the size of the MPI communicator is reduced successively to minimize
communication. Re-discretization of the differential equations is performed on each coarser spectral and
geometric level. The viscosity values in each element are stored at the quadrature points of the velocity
discretization, and are thus local to each element. The viscosity coarsening is done level-by-level during the
setup phase. The coarsening operator is the adjoint of the refinement operator, which performs element-wise
interpolation. This adjoint is computed with respect to the L2 -inner products, and since viscosity values are
not shared amongst elements, this does not require (an approximation of) a global mass matrix solve. The
transition from geometric to algebraic multigrid is done at a sufficiently small core count and small MPI
communicator. AMG continues to further reduce problem size (via Galerkin coarse grid projection) and the
number of cores down to a single core for the direct solver.
The operator Kwr can be regarded as a discrete, variable-coefficient Poisson operator on the discontinuous
pressure space Pdisc
k−1 with Neumann boundary conditions. Therefore, multigrid V-cycles can also be employed
to approximate the inverse K̃−1
wr . However, it turned out to be problematic to apply multigrid coarsening
directly due to the discontinuous, modal discretization of the pressure. We took a novel approach in [36] by
considering the underlying infinite-dimensional, variable-coefficient Poisson operator, where the coefficient
−1
is derived from the diagonal weighting matrix (here, C−1
wl or Dwr ). Then we re-discretize with continuous,
nodal high-order finite elements in Qk . An alternative would be to use Qk−1 , but we prefer to use Qk since
the corresponding data structures are readily available from the discretization of the velocity. Hence, this
choice avoids the setup cost related to discretization-specific parameters and their storage. Additionally, the
HMG hierarchy of the preconditioner acting on the velocity can be partially reused, again saving setup time
and memory. This continuous, nodal discretization of the Poisson operator is then approximately inverted
with an HMG V-cycle that is similar to the one described above for the inverse viscous block approximation
Ã−1 . Additional smoothing is applied in the discontinuous pressure space (Figure 4, green level ) to account
for high frequency modes in residuals that are introduced through projections between Qk and Pdisc
k−1 .
Our hybrid multigrid method combines high-order L2 -restrictions/interpolations and employs Chebyshevaccelerated point-Jacobi smoothers. This results in optimal or nearly optimal algorithmic multigrid performance (see Section 7), i.e., iteration numbers are independent of mesh size and only mildly dependent
on discretization order, while maintaining robustness with respect to highly heterogeneous coefficients. In
addition, the efficacy of the HMG preconditioner does not deteriorate with increasing core counts, because
the spectral and geometric multigrid is by construction independent of the number of cores and AMG is
invoked for prescribed small problem sizes on essentially fixed small core counts.
For all numerical experiments presented here, three pre- and post-smoothing iterations with a Chebyshev
accelerated point-Jacobi smoother are performed. AMG is always invoked on just one MPI rank after
geometrically coarsening the uniform mesh to ` = 2 with 64 elements (using a direct solver for such small
problems is also reasonable). GMRES is restarted after every 100 iterations throughout all experiments.
PETSc’s [2] implementations of Chebyshev acceleration, direct solver, AMG (called GAMG), and GMRES
are used.

w-BFBT Preconditioner for Stokes Flow Problems

18
HMG V-cycle

HMG hierarchy

modal to
nodal proj.

pressure space

discont. modal

spectral
p-coarsening

cont. nodal
high-order F.E.

p-MG

geometric
h-coarsening

trilinear F.E.
decreasing #cores

h-MG
AMG

algebraic
#cores < 1000
coars. small MPI communicator
single core

direct

high-order
L2 -projection
linear
L2 -projection
linear
projection

Figure 4: Hybrid spectral-geometric-algebraic multigrid (HMG). Left: Illustration of multigrid hierarchy.
From top to bottom, first, the multigrid levels are obtained by spectral coarsening (dark blue). Next, the
mesh is geometrically coarsened and repartitioned on successively fewer cores to minimize communication
(light blue). Finally, AMG further reduces problem size and core count (light red). The multigrid hierarchy
for the pressure Poisson operator Kwr additionally involves smoothing in the discontinuous, modal pressure
space (green). The projection from the discontinuous, modal to a continuous finite element nodal basis uses a
lumped mass matrix in the nodal space to avoid the global mass matrix system solve. Right: The multigrid
V-cycle consists of smoothing at each level of the hierarchy (circles) and intergrid transfer operators (arrows
downward for restriction and arrows upward for interpolation). To enhance efficacy of the the V-cycle as a
preconditioner, different types of projection operators are employed for these operators depending on the
phase within the V-cycle.

7

Algorithmic and parallel scalability for HMG + w-BFBT Stokes
preconditioner

After establishing the robustness of the Stokes solver with w-BFBT preconditioning in theory (Section 3)
and numerically (Section 4), addressing issues associated with Dirichlet boundary conditions (Section 5),
and introducing an effective and scalable multigrid method (Section 6), in this section we finally study the
scalability of the Stokes solver building on HMG+w-BFBT. One aspect of scalability is algorithmic scalability,
i.e., the dependence of Krylov iterations on the mesh resolution and the discretization order. The second
aspect is parallel scalability of the implementation, i.e., runtime measured on increasing numbers of compute
cores. Studying both aspects is required to fully assess the performance of a solver at scale.
The algorithmic scalability in Table 3 shows results for the Stokes solver as well as its individual components
by reporting iteration numbers for solving the systems Au = f and Kwr p = g. Studying the individual
components allows us to observe HMG performance in isolation and to compare it to the algorithmic scalability
of the full Stokes system, which is indicative of the quality of the w-BFBT Schur complement approximation.
All systems are solved with preconditioned GMRES down to a relative tolerance of 10−6 . The preconditioners
for A and Kwr are HMG-V-cycles as described in Section 6. For the w-BFBT preconditioner, we set a
constant left boundary amplification al = 1 and vary the right boundary amplification ar according to results
from Section 5. The iteration counts in Table 3a show textbook mesh independence when increasing the
level of refinement `. This holds for each component, A and Kwr , and also the whole Stokes solver, and
hence we conclude that the Schur complement approximation by w-BFBT is mesh-independent. When the
discretization order k is increased, the iteration counts presented in Table 3b increase mildly. The convergence
of both components A and Kwr exhibits a moderate dependence on k. Since the increase in number of
iterations is sightly larger for the full Stokes solve than for A and Kwr , we suspect a mild deterioration of
w-BFBT as a Schur complement approximation.
The following parallel scalability results complement the results presented in [36], where scalability to
millions of threads was demonstrated on IBM’s BlueGene/Q architecture, which is specifically designed for
large-scale supercomputers and thus differs from conventional clusters. The parallel scalability results here
were obtained on the full Lonestar 5 peta-scale system, which represents a rather conventional cluster, housed

w-BFBT Preconditioner for Stokes Flow Problems

19

Table 3: Algorithmic scalability for Stokes solver with HMG+w-BFBT preconditioning while a varying mesh
refinement level ` and b varying discretization order k (problem S16-rand, DR(µ) = 106 ). Computational
cost is expressed in number of GMRES iterations (abbreviated by It.) for full Stokes solve (10−6 residual
reduction). Left boundary amplification for Cwl is fixed to al = 1; right boundary amplification ar for
Dwr varies. Additionally, the numbers of GMRES iterations for solving only the sub-systems Au = f and
Kwr p = g are given for demonstration of HMG efficacy (here, f is the right-hand side of the momentum
equation and g is the discreet representation of ∇ · f ; however, random right-hand sides would give similar
convergence results).
(b) Algorithmic scalability (fixed level ` = 5)

(a) Algorithmic scalability (fixed order k = 2)
`

ar

u-DOF
[×106 ]

It.
A

p-DOF
[×106 ]

It.
Kwr

4
5
6
7
8
9
10

1
2
4
8
16
32
64

0.11
0.82
6.44
50.92
405.02
3230.67
25807.57

18
18
18
18
18
18
18

0.02
0.13
1.05
8.39
67.11
536.87
4294.97

8
7
6
6
6
6
6

DOF
[×106 ]
0.12
0.95
7.49
59.31
472.12
3767.54
30102.53

It.
Stokes

k

ar

40
33
33
34
34
34
34

2
3
4
5
6
7
8

2
4
8
16
32
64
128

u-DOF
[×106 ]

It.
A

p-DOF
[×106 ]

0.82
2.74
6.44
12.52
21.56
34.17
50.92

18
20
20
23
23
22
22

0.13
0.32
0.66
1.15
1.84
2.75
3.93

It.
Kwr
7
8
7
12
12
10
10

DOF
[×106 ]
0.95
3.07
7.10
13.67
23.40
36.92
54.86

It.
Stokes
33
37
36
43
50
54
67

at the Texas Advanced Computing Center (TACC). The Lonestar 5 supercomputer entered production in
January 2016 and is a Cray XC40 system consisting of 1252 compute nodes. Each node is equipped with two
Intel Haswell 12-core processors (Xeon E5-2680v3) and 64 GBytes of memory. Inter-node communication is
based on an Aries Dragonfly topology network that provides dynamic routing and thus enables optimal use
of the system bandwidth.
Results for weak scalability (DOF/core fixed to ∼1 million) in Figure 5a show that the Stokes solver with
w-BFBT (blue curve) maintains 90% parallel efficiency over a 618-fold increase in degrees of freedom along
with cores. Even for the setup of the Stokes solver (green curve), which mainly involves generation of the
HMG hierarchy, we observe 71% parallel efficiency. These are excellent results for such a complex implicit
multilevel solver with optimal algorithmic performance (when the mesh is refined, or nearly algorithmically
optimal when the order is increased) and with convergence that is independent of the number of cores.
Finally, Figure 5b reports strong scalability results (overall DOF fixed to 59 million) and how the number
of OpenMP threads substituting MPI ranks influences speedup.4 Over the 78-fold increase from 48 to 3744
cores, efficiency reduces moderately, to a worst-case 68% for 24×OMP1. However, note that in the largest
run with 29,640 cores, the granularity is only ∼2000 DOF/core, which is extremely challenging for strong
scalability. In this case, due to the increased communication volume, overlapping with decreased amounts of
computation becomes impossible and communication dominates the runtime. This behavior is expected for
an implicit solver, especially for a multilevel method that does not sacrifice algorithmic optimality for parallel
scalability.

8

Conclusions

For Stokes flow problems with highly heterogeneous viscosity, the commonly used inverse viscosity-weighted
mass matrix approximation of the Schur complement can be insufficient. As a consequence, convergence
of Schur complement-based Krylov solvers can be extremely slow. While each iteration with the weighted
BFBT (w-BFBT) Schur complement approximation proposed in this paper is computationally more costly,
we observe that it results in robust and fast convergence even for complex viscosity structures and for up
to ten orders of magnitude viscosity contrasts—properties than occur, for instance, in problems involving
non-Newtonian geophysical fluids. Our numerical findings are supported by theoretical spectral equivalence
results. To the best of our knowledge, a similar analysis has not been shown for any BFBT-type method
before.
4 Even though the processors of Lonestar 5 support two threads per physical core (Intel Hyper-Threading Technology),
assigning more than one OpenMP thread per core did not improve performance.

Degrees of freedom per unit time

w-BFBT Preconditioner for Stokes Flow Problems

1010

Ideal weak scalability
Solve [DOF/(sec/iter)]
Setup [DOF/sec]

0.90
0.94

109
108

20

0.71

0.97
0.94

107
106
48

0.97
456

3744

29,640

Number of cores (∼ 106 DOF/core)

Speedup [baseline sec/sec]

(a) Weak scalability (k = 2, ` = 7, . . . , 10)

512

64

Ideal speedup
Solve 24×OMP1
Solve 6×OMP4
Solve 2×OMP12
Solve 1×OMP24

1.20
8

1
48
(1,235,675)

1.08
0.68

0.18

0.05

1.00

456
(130,071)

3744
(15,842)

29,640
(2009)

Number of cores (DOF/core)

(b) Strong scalability (k = 2, ` = 7)

Figure 5: Parallel scalability on Lonestar 5 for Stokes solver with HMG+w-BFBT preconditioning (problem
S16-rand, DR(µ) = 106 as in Table 3a). a Weak scalability of setup and solve phases (normalized w.r.t.
deviations from const. DOF/core). Numbers along the graph lines indicate weak parallel efficiency w.r.t. ideal
weak scalability (baseline is 48 cores result). DOF/core is ∼1 million; the largest problem size on 29,640
cores has 30 billion DOF. b Strong scalability of solve phase for different configurations of OpenMP threads
(OMP) substituting MPI ranks on each node consisting of 24 cores. Numbers along the graph lines indicate
strong efficiency w.r.t. ideal speedup (baseline is 48 cores result).

w-BFBT Preconditioner for Stokes Flow Problems

21

At Dirichlet boundaries, we use a modification of w-BFBT that is necessary for mesh-independent
convergence. In this modification, we dampen the influence of elements at these boundaries, which would
otherwise dominate in w-BFBT and degrade its efficacy as a preconditioner. Finding an alternative remedy
is a subject of current research.
Indefinite problems with highly heterogeneous coefficients and high-order discretizations present significant
challenges for efficient solvers, especially in parallel. Nevertheless, we have demonstrated that with careful
attention paid to algorithm design and scalable implementation, viscosity-robust and mesh-independent
solvers can be designed that exhibit nearly optimal algorithmic and parallel scalability over a wide range of
problems sizes and core counts.

Acknowledgments
We wish to offer our deep thanks to the Texas Advanced Computing Center (TACC) for granting us early
user access to Lonestar 5 and enabling runs on the full system. We appreciate many helpful discussions with
Tobin Isaac, who also allowed us to build on his previous work on Stokes solvers. We also thank Carsten
Burstedde for his dedicated work on the p4est library.

References
[1] Gabriel Acosta, Ricardo G. Durán, and Ariel L. Lombardi, Weighted Poincaré and Korn
inequalities for Hölder α domains, Mathematical Methods in the Applied Sciences, 29 (2006), pp. 387–400.
[2] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris
Buschelman, Lisandro Dalcin, Victor Eijkhout, William D. Gropp, Dinesh Kaushik,
Matthew G. Knepley, Lois Curfman McInnes, Karl Rupp, Barry F. Smith, Stefano
Zampini, and Hong Zhang, PETSc users manual, Tech. Report ANL-95/11 - Revision 3.6, Argonne
National Laboratory, 2015.
[3] Wolfgang Bangerth and Timo Heister, ASPECT: Advanced Solver for Problems in Earth’s
ConvecTion, Computational Infrastructure in Geodynamics, 2015.
[4] Michele Benzi, Gene H. Golub, and Jörg Liesen, Numerical solution of saddle point problems,
Acta Numerica, 14 (2005), pp. 1–137.
[5] Carsten Burstedde, Omar Ghattas, Michael Gurnis, Tobin Isaac, Georg Stadler, Tim
Warburton, and Lucas C. Wilcox, Extreme-scale AMR, in SC10: Proceedings of the International
Conference for High Performance Computing, Networking, Storage and Analysis, ACM/IEEE, 2010.
[6] Carsten Burstedde, Omar Ghattas, Michael Gurnis, Eh Tan, Tiankai Tu, Georg Stadler,
Lucas C. Wilcox, and Shijie Zhong, Scalable adaptive mantle convection simulation on petascale
supercomputers, in SC08: Proceedings of the International Conference for High Performance Computing,
Networking, Storage and Analysis, ACM/IEEE, 2008.
[7] Carsten Burstedde, Omar Ghattas, Georg Stadler, Tiankai Tu, and Lucas C. Wilcox,
Parallel scalable adjoint-based adaptive solution for variable-viscosity Stokes flows, Computer Methods in
Applied Mechanics and Engineering, 198 (2009), pp. 1691–1700.
[8] Carsten Burstedde, Georg Stadler, Laura Alisic, Lucas C. Wilcox, Eh Tan, Michael
Gurnis, and Omar Ghattas, Large-scale adaptive mantle convection simulation, Geophysical Journal
International, 192 (2013), pp. 889–906.
[9] Carsten Burstedde, Lucas C. Wilcox, and Omar Ghattas, p4est: Scalable algorithms for
parallel adaptive mesh refinement on forests of octrees, SIAM Journal on Scientific Computing, 33 (2011),
pp. 1103–1133.
[10] Seng-Kee Chua and Richard L. Wheeden, Estimates of best constants for weighted Poincaré
inequalities on convex domains, Proceedings of the London Mathematical Society, 93 (2006), pp. 197–226.

w-BFBT Preconditioner for Stokes Flow Problems

22

[11] Michel Crouzeix and Pierre-Arnaud Raviart, Conforming and nonconforming finite element
methods for solving the stationary Stokes equations. I, Rev. Française Automat. Informat. Recherche
Opérationnelle Sér. Rouge, 7 (1973), pp. 33–75.
[12] Jean Donea and Antonio Huerta, Finite Element Methods for Flow Problems, John Wiley & Sons,
2003.
[13] Howard C. Elman, Preconditioning for the steady-state Navier–Stokes equations with low viscosity,
SIAM Journal on Scientific Computing, 20 (1999), pp. 1299–1316.
[14] Howard C. Elman, Victoria E. Howle, John Shadid, Robert Shuttleworth, and Raymond S.
Tuminaro, Block preconditioners based on approximate commutators, SIAM Journal on Scientific
Computing, 27 (2006), pp. 1651–1668.
[15]

, A taxonomy and comparison of parallel block multi-level preconditioners for the incompressible
Navier–Stokes equations, Journal of Computational Physics, 227 (2008), pp. 1790–1808.

[16] Howard C. Elman, David J. Silvester, and Andrew J. Wathen, Finite elements and fast iterative
solvers: with applications in incompressible fluid dynamics, Oxford University Press, 2014.
[17] Howard C. Elman and Raymond S. Tuminaro, Boundary conditions in approximate commutator
preconditioners for the Navier–Stokes equations, Electronic Transactions on Numerical Analysis, 35
(2009), pp. 257–280.
[18] Mikito Furuichi, Dave A. May, and Paul 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.
[19] Roland Glowinski and Jinchao Xu, Numerical Methods for Non-Newtonian Fluids: Special Volume,
vol. 16 of Handbook of Numerical Analysis, North-Holland, 2011.
[20] Piotr P. Grinevich and Maxim A. Olshanskii, An iterative method for the Stokes-type problem
with variable viscosity, SIAM Journal on Scientific Computing, 31 (2009), pp. 3959–3978.
[21] Vincent Heuveline and Friedhelm 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.
[22] Kolumban Hutter, Theoretical Glaciology, Mathematical Approaches to Geophysics, D. Reidel
Publishing Company, 1983.
[23] Tobin Isaac, Noemi Petra, Georg Stadler, and Omar 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.
[24] Tobin Isaac, Georg Stadler, and Omar Ghattas, Solution of nonlinear Stokes equations discretized
by high-order finite elements on nonconforming and anisotropic meshes, with application to ice sheet
dynamics, SIAM Journal on Scientific Computing, 37 (2015), pp. B804–B833.
[25] Volker John, Kristine Kaiser, and Julia Novo, Finite element methods for the incompressible
Stokes equations with variable viscosity, ZAMM - Zeitschrift für Angewandte Mathematik und Mechanik,
96 (2016), pp. 205–216.
[26] David Kay, Daniel Loghin, and Andrew Wathen, A preconditioner for the steady-state NavierStokes equations, SIAM Journal on Scientific Computing, 24 (2002), pp. 237–256.
[27] Martin Kronbichler, Timo Heister, and Wolfgang Bangerth, High accuracy mantle convection
simulation through modern numerical methods, Geophysical Journal International, 191 (2012), pp. 12–29.

w-BFBT Preconditioner for Stokes Flow Problems

23

[28] Anders Logg, Kent-Andre Mardal, and Garth Wells, Automated Solution of Differential
Equations by the Finite Element Method: The FEniCS book, vol. 84, Springer Science & Business Media,
2012.
[29] Dave A. May, Jed Brown, and Laetitia Le Pourhiet, pTatin3D: High-performance methods
for long-term lithospheric dynamics, in SC14: Proceedings of the International Conference on High
Performance Computing, Networking, Storage and Analysis, IEEE Press, 2014, pp. 274–284.
[30] Dave A. May, Jed Brown, and Laetitia Le Pourhiet, A scalable, matrix-free multigrid preconditioner for finite element discretizations of heterogeneous Stokes flow, Computer Methods in Applied
Mechanics and Engineering, 290 (2015), pp. 496–523.
[31] Dave A. May and Louis Moresi, Preconditioned iterative methods for Stokes flow problems arising in
computational geodynamics, Physics of the Earth and Planetary Interiors, 171 (2008), pp. 33–47.
[32] Dan McKenzie, The generation and compaction of partially molten rock, Journal of Petrology, 25
(1984), pp. 713–765.
[33] Clemens Pechstein and Robert Scheichl, Weighted Poincaré inequalities, IMA Journal of Numerical
Analysis, 33 (2013), pp. 652–686.
[34] Noemi Petra, Hongyu Zhu, Georg Stadler, Thomas J. R. Hughes, and Omar 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.
[35] Kumbakonam R. Rajagopal, Mechanics of non-Newtonian fluids, in Recent Developments in Theoretical Fluid Mechanics, G.P. Galdi and J. Necas, eds., vol. 291, Longman’s Scientific and Technical,
1993, pp. 129–162.
[36] Johann Rudi, A. Crisiano I. Malossi, Tobin Isaac, Georg Stadler, Michael Gurnis, Yves
Ineichen, Costas Bekas, Alessandro Curioni, and Omar Ghattas, An extreme-scale implicit
solver for complex PDEs: Highly heterogeneous flow in earth’s mantle, in SC15: Proceedings of the
International Conference for High Performance Computing, Networking, Storage and Analysis, ACM,
2015, pp. 5:1–5:12.
[37] Gerald Schubert, Donald L. Turcotte, and Peter Olson, Mantle Convection in the Earth and
Planets, Cambridge University Press, 2001.
[38] David J. Silvester, Howard C. Elman, David Kay, and Andrew J. Wathen, Efficient preconditioning of the linearized Navier–Stokes equations for incompressible flow, Journal of Computational and
Applied Mathematics, 128 (2001), pp. 261–279. Numerical analysis 2000, Vol. VII, Partial differential
equations.
[39] Georg Stadler, Michael Gurnis, Carsten Burstedde, Lucas C. Wilcox, Laura Alisic, and
Omar Ghattas, The dynamics of plate tectonics and mantle flow: From local to global scales, Science,
329 (2010), pp. 1033–1038.
[40] Rolf Stenberg and Manil Suri, Mixed hp finite element methods for problems in elasticity and
Stokes flow, Numerische Mathematik, 72 (1996), pp. 367–389.
[41] Hari Sundar, George Biros, Carsten Burstedde, Johann Rudi, Omar Ghattas, and Georg
Stadler, Parallel geometric-algebraic multigrid on unstructured forests of octrees, in SC12: Proceedings
of the International Conference for High Performance Computing, Networking, Storage and Analysis,
Salt Lake City, UT, 2012, ACM/IEEE.
[42] Rüdiger Verfürth, A Posteriori Error Estimation Techniques for Finite Element Methods, Oxford
University Press, 2013.
[43] Jennifer Worthen, Georg Stadler, Noemi Petra, Michael Gurnis, and Omar Ghattas,
Towards adjoint-based inversion for rheological parameters in nonlinear viscous mante flow, Physics of
the Earth and Planetary Interiors, 234 (2014), pp. 23–34.

