Finite Volume Methods: Foundation and Analysis

Timothy Barth1 , Raphaèle Herbin2 and Mario Ohlberger3

1 NASA Ames Research Center, Moffett Field, CA, USA

2 Aix-Marseille
Université, CNRS, Centrale Marseille, Marseille, France
3 Applied Mathematics Münster, CeNoS, and CMTC, University of Münster, Münster, Germany

ABSTRACT

Finite volume methods are a class of discretization schemes resulting from the decomposition of a
problem domain into nonoverlapping control volumes. Degrees of freedom are assigned to each control
volume that determine local approximation spaces and quadratures used in the calculation of control
volume surface fluxes and interior integrals. An imposition of conservation and balance law statements
in each and every control volume constrains surface fluxes and results in a coupled system of equations
for the unknown degrees of freedom that must be solved by a numerical method.

Finite volume methods have proved highly successful in approximating the solution to a wide
variety of conservation and balance laws. They are extensively used in fluid mechanics, meteorology,
electromagnetics, semiconductor device simulation, materials modeling, heat transfer, models of
biological processes, and many other engineering problems governed by conservation and balance
laws that may be written in integral control volume form.

This chapter reviews elements of the foundation and analysis of modern finite volume methods
for approximating hyperbolic, elliptic, and parabolic partial differential equations. These different
equations have markedly different continuous problem regularity and function spaces (e.g., L∞ , L2 ,
and H01 ) that must be adequately represented in finite-dimensional discretizations. Particular attention
is given to finite volume discretizations yielding numerical solutions that inherit properties of the
underlying continuous solutions such as maximum (minimum) principles, total variation control, L2
stability, global entropy decay, and local balance law conservation while also having favorable accuracy
and convergence properties on structured and unstructured meshes.

As a starting point, a review of scalar nonlinear hyperbolic conservation laws and the development
of high-order accurate schemes for discretizing them is presented. A key tool in the design and analysis
of finite volume schemes suitable for discontinuity capturing is discrete maximum principle analysis. A
number of mathematical and algorithmic developments used in the construction of numerical schemes
possessing local discrete maximum principles are reviewed in one and several space dimensions.
These developments include monotone fluxes, TVD discretization, positive coefficient discretization,
nonoscillatory reconstruction, slope limiters, strong stability preserving time integrators, and so on.
When available, theoretical results concerning a priori and a posteriori error estimates and convergence
to entropy weak solutions are given.

A review of the discretization of elliptic and parabolic problems is then presented. The tools needed

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
2 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

for the theoretical analysis of the two point flux approximation scheme for the convection diffusion
equation are described. Such schemes require an orthogonality condition on the mesh in order for the
numerical fluxes to be consistent. Under this condition, the scheme may be shown to be monotone. A
weak formulation of the scheme is derived, which facilitates obtaining stability, convergence, and error
estimate results. The discretization of anisotropic problems is then considered and a review is given
of some of the numerous schemes that have been designed in recent years, along with their properties.
Parabolic problems are then addressed, both in the linear and nonlinear cases.

A discussion of further advanced topics is then given including the extension of the finite volume
method to systems of hyperbolic conservation laws. Numerical flux functions based on an exact or
approximate solution of the Riemann problem of gas dynamics are discussed. This is followed by the
review of another class of numerical flux functions for symmetrizable systems of conservation laws
that yield finite volume solutions with provable global decay of the total mathematical entropy for a
closed entropy system, often referred to as entropy stability.

Finally, a detailed review of the discretization of the steady-state incompressible Navier–Stokes

equations using the Marker-And-Cell (MAC) finite volume method is then presented. The MAC
scheme uses a staggered mesh discretization for pressure and velocities on primal and dual control
volumes. After reformulating the MAC scheme in weak form, analysis results concerning stability,
weak consistency, and convergence are given.

key words: finite volume methods, conservation laws, elliptic and parabolic equations,
nonoscillatory approximation, discrete maximum principles, higher order schemes

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 3

Contents

1 Introduction 3

2 Scalar Nonlinear Hyperbolic Conservation Laws 6

2.1 The method of characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.2 Weak solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.3 Entropy weak solutions and vanishing viscosity . . . . . . . . . . . . . . . . . . 9

2.4 Measure-valued or entropy process solutions . . . . . . . . . . . . . . . . . . . . 10

3 Finite Volume Methods for Nonlinear Hyperbolic Conservation Laws 11

3.1 Finite volume discretization from cell averages via exact or approximate
Riemann problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3.2 Monotone fluxes and E-flux functions . . . . . . . . . . . . . . . . . . . . . . . . 14

3.3 Stability, convergence, and error estimates . . . . . . . . . . . . . . . . . . . . . 19

4 Higher Order Accurate Finite Volume Methods for Hyperbolic Problems 24

4.1 Higher order accurate finite volume methods for hyperbolic problems in one
dimension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

4.2 Higher order accurate finite volume methods for hyperbolic problems in multiple
dimensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.3 Higher order accurate finite volume methods for hyperbolic problems on
unstructured meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

4.4 Higher order accurate time integration schemes . . . . . . . . . . . . . . . . . . 45

5 Finite Volume Methods for Elliptic and Parabolic Problems 47

5.1 Convergence analysis for the steady state reaction convection diffusion equation 48

5.2 Discretization of anisotropic elliptic problems . . . . . . . . . . . . . . . . . . . 56

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
4 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

5.3 The parabolic case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

6 Advanced Topics 68

6.1 Extension to systems of nonlinear hyperbolic conservation laws . . . . . . . . . 68

6.2 The Marker-and-Cell scheme for fluid flows . . . . . . . . . . . . . . . . . . . . 76

7 Related Chapters 86

1. Introduction

Finite volume methods (FVMs) are a popular class of discretization schemes that are
well suited to approximating conservation and balance laws. These laws may yield partial
differential equations (PDEs) of different type (hyperbolic, elliptic, or parabolic) as well as
coupled systems of equations with individual equations of different type. Consequently, the
regularity of the solution to these equations may be quite different from one another and so
too the functional spaces in which the solutions of the continuous problems are sought, viz. L∞ ,
H01 , L2 , and so on. While the piecewise approximation spaces used in FVMs (e.g., piecewise
constant spaces for the simplest FVMs) are natural candidates for hyperbolic problems,
these approximation spaces are not natural candidates for elliptic problems in H01 . But as
is revealed by mathematical analysis, these piecewise approximation spaces, in particular
piecewise constant spaces, are still provably viable candidates for problems in H10 when the
mesh satisfies certain technical requirements.

A question that is often asked by a non-expert concerns the differences between the finite
volume method, finite element method (FEM), and the finite difference method (FDM).
The answer truly lies in the concepts of the methods, but in some cases, these methods
do yield similar schemes. This similarity may be seen using the simple example u00 = f
discretized by all three of the methods using a constant mesh spacing on the unit interval
[0, 1]. Roughly speaking, one could say that the FEM is based on a weak formulation coupled
with a convenient finite-dimensional approximation of the infinite-dimensional function spaces.
The FDM relies on an approximation of the differential operators using Taylor expansions. The
FVM is constructed from a balance equation, rather than the PDE itself, with a consistent
approximation of the fluxes defined on the boundary of the control volume on which the
balance equation is written.

Confusion between the FVM and the FDM arises from the fact that the FVM is sometimes
called an FDM when the fluxes on the boundary of each control volume are approximated by
finite differences. This is sometimes the case, for instance, in oil reservoir simulations utilizing
isotropic diffusion models that are discretized on Cartesian grids such that the diffusion flux can

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 5

be easily approximated by a simple difference quotient. Moreover, numerous schemes have been
designed for hyperbolic problems that are called FDMs, although they can also be interpreted
as FVMs with suitable approximation of the fluxes at the interfaces of the discretization
control volume (also sometimes referred to as a “cell”). Links between the FVM and the finite
element method (FEM) can also be found. In certain instances, the FVM can be interpreted as
an FEM using a particular integration rule. Conversely, there are instances where the FEM can
be interpreted as FVM. For example, the piecewise linear FEM discretization of the Laplace
operator on a triangular mesh satisfying the weak Delaunay triangulation condition yields a
matrix that is the same as that of the FVM on the dual Voronoı̈ mesh; see Eymard et al.,
2000 for details. As another example, the FVM is sometimes presented as a discontinuous
Galerkin method (DGM) of lowest order that uses a finite-dimensional approximation of the
continuous space that is nonconforming. This is mathematically insightful, but the tools used to
analyze DGMs of higher order accuracy do not seem to directly apply to FVMs of higher order
accuracy. Other families of FVMs have been developed such as vertex-centered schemes, box
or covolume schemes, and finite volume element (FVE) methods that facilitate the compact
discretization of various differential operators. Particular attention is given in this chapter to
cell centered schemes because these schemes are widely used in industrial codes and are well
suited to the discretization of conservation laws of the general form
∂t u + ∇ · f (u, ∇u) + s(u) = 0 (1)
where u is a function of space and time, f ∈ C 1 (R × Rd , Rd ) is the flux function, s ∈ C(R, R) is
a low-order source term, and d is the space dimension. This conservation law may be obtained
from the following balance equation written for a control volume K with exterior boundary
normal nK Z Z Z
∂t udx + f (u, ∇u) · nK ds + s(u)dx = 0 (2)
K ∂K K
by letting the size of K tend to zero. In the above integrals and in the sequel, dx represents
the integration symbol on a d-dimensional subset of Rd and ds on a d − 1-dimensional subset
of Rd . Note that conversely, the balance equation (2) may be obtained from the conservation
law (1) using integration over a control volume K and applying the Stokes formula. If the
control volume K is a polytope (a polygon in 2D or a polyhedron in 3D), then the boundary
is the union of faces (or edges in 2D), denoted here by σ, so that (2) may be written as
Z X Z Z
∂t udx + f (u, ∇u) · nK ds + s(u)dx = 0 (3)
K σ⊂∂K σ K

Replacing the continuous time derivative with an explicit Euler time discretization with
uniform time step δt yields
Z X Z Z
1
(un+1 − un )dx + f (un , ∇un ) · nK ds + s(un )dx = 0
δt K σ K
σ⊂∂K
n
where u denotes an approximation of u at time tn = nδt. For each time tn and control volume
K, the discrete unknown unK approximates u in the control volume K at time tn = nδt. To
obtain the approximate
R equations needed to solve for unK (which defines the numerical scheme),
the flux integrals σ f (un , ∇un ) · nK ds must be discretized. Let FK,σ (un ) denote a numerical
flux that approximates f . A nontrivial task in developing a new FVM scheme is to devise

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
6 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

a numerical flux so that properties such as discrete conservation, consistency, accuracy, and
convergence (discussed later) are obtained from the resulting discretization. To illustrate the
task of devising a numerical flux, consider the linear convection equation that is obtained from
(2) by setting f (u, ∇u) = vu, where v is a constant vector of Rd and s(u) = 0. The balance
equation then reduces to the following simple form
Z Z
∂t udx + uv · nK ds = 0 (4)
K ∂K

In order to approximate the flux uv · nK on the faces of each control volume, one needs to
approximate the value of u on these edges as a function of the discrete unknowns uK associated
to each control volume K. This may be done in several ways. A straightforward choice is to
approximate the value of u on the face σ = σKL separating the control volumes K and L by
the mean value 12 (uK + uL ). This yields the so-called “centered” numerical flux
(cv,c) 1
FK,σ = vK,σ (uK + uL )
2
R
where vK,σ = σ v · nK ds. This centered choice is known to lead to stability problems and is
therefore not used in practice. A popular choice is the so-called “upwind” numerical flux given
by
(cv,c) + −
FK,σ = vK,σ uK − vK,σ uL
where x+ = max(x, 0) and x− = − min(x, 0). Note that this formula is equivalent to

 vK,σ uK if vK,σ ≥ 0
(cv,c)
FK,σ =

vK,σ uL otherwise
This numerical flux results in schemes satisfying the desired properties mentioned above. A
linear convection diffusion reaction balance equation can be obtained from (2) by setting
f (u, ∇u) = ∇u + vu, v ∈ Rd and s(u) = bu, b ∈ R,
∂t u − ∆u + div(vu) + bu = 0 on Ω
The flux through a given edge is then given by
Z Z
f (u) · nK,σ ds = (−∇u · nK,σ + v · nK,σ u) ds
σ σ
R
so that the additional diffusion term σ −∇u · nK,σ ds, involving the normal derivative to the
boundary of a control volume, must now be discretized. On a Cartesian grid, a possible simple
discretization is obtained using the difference quotient between the value of u in K and an
adjacent control volume L, that is,
(d) |σ|
FK,σ = − (uL − uK ) (5)
dKL
where |σ| stands for the (d−1)dimensional Lebesgue measure of σ (area if d = 3, length if d = 2)
and dKL is the distance between some (well-chosen) points of K and L. The numerical flux is
known in the porous media community as the “two point” (TP) flux, and the resulting scheme
as the “two point flux approximation” (TPFA) scheme. If the points that are used to compute

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 7

the distance dKL are carefully chosen, then the resulting diffusion flux (5) is consistent in the
finite difference sense (note, however, that the resulting approximation of the second-order
diffusion operator may not necessarily be consistent in the finite difference sense). Results
presented in Section 5 reveal that TP fluxes for the discretization of the diffusion flux yield
accurate results if the mesh satisfies an orthogonality condition. This orthogonality condition
requires that there exists a family of points xK , such that for a given interface σKL between
the control volumes K and L, the line segment xK xL is orthogonal to the edge σKL . The
length dKL is then defined as the distance between xK and xL . Such a family of points exists,
for instance, in the case of triangles, rectangles, or Voronoı̈ meshes, but not for general meshes.
For general meshes and for anisotropic diffusion problems, a wide variety of schemes have been
introduced in the recent years that are reviewed in Section 5.

The convergence analysis of FVMs is rather recent. Since these methods are currently
employed in nonlinear problems where the regularity of the solution is not clear, one would like
to obtain theoretical results on the convergence of the scheme without (nonphysical) regularity
assumptions on the data or solution. The usual path to the proof of convergence of FVMs,
which was initiated in Eymard, Gallouët, and Herbin (1999), Champier et al., 1993, and now
been used in a wide number of FVM papers and also adapted to other schemes, is based on
establishing the following set of theoretical results:

1. Establish a priori estimates on the solution to the scheme in a mesh-dependent norm

and deduce the existence of a solution to the scheme.
2. Prove a compactness result.
3. Prove a realistic regularity property of any possible limit.
4. By a passage to the limit in the scheme, prove that any possible limit satisfies a weak
form of the original PDE.

Note that a by-product of this approach is an existence proof for the original PDE. Even
though the existence is sometimes known before attacking the discretization of the problem,
it can be and sometimes has been proved by a numerical approximation technique. This is the
case for the now historical result on the existence of the Laplace and biharmonic equations by
the convergence of a finite difference approximation in the famous paper Courant et al. (1928)
(see Courant, Friedrichs and Lewy (1967) for its English version). The next natural question
is the rate of convergence of the scheme. Theoretical results on error estimates are often linked
to a uniqueness result and are therefore not always accessible.

2. Scalar Nonlinear Hyperbolic Conservation Laws

Many problems arising in science and engineering lead to the study of nonlinear hyperbolic
conservation laws. Some examples include fluid mechanics, meteorology, electromagnetics,
semiconductor device simulation, and numerous models of biological processes. As a prototype
conservation law, consider a flux function f depending only on u, and the Cauchy initial value

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
8 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

problem
∂t u + ∇ · f (u) =0 in Rd × R+ (6a)
u(x, 0) = u0 (x) in Rd (6b)
Here u(x, t) : Rd × R+ → R denotes the dependent solution variable, f ∈ C 1 (R, Rd ) denotes
the flux function, and u0 (x) : Rd → R the initial data.

The function u is a classical solution of the scalar initial value problem if u ∈ C 1 (Rd × R+ )
satisfies (6) pointwise. An essential feature of nonlinear conservation laws is that, in general,
gradients of u blow up in finite time, even when the initial data u0 is arbitrarily smooth.
Beyond some critical time t0 classical solutions of (6) do not exist. This behavior will be
demonstrated shortly using the method of characteristics. By introducing the notion of weak
solutions of (6) together with an entropy condition, it then becomes possible to define a class
of solutions where existence and uniqueness are guaranteed for times greater than t0 . These
are precisely the solutions that are numerically sought in the finite volume method.

2.1. The method of characteristics

Let u be a classical solution of (6). Further, define the vector

a(u) = f 0 (u) = (f10 (u), . . . , fd0 (u))T
A characteristic Γξ is a curve (x(t), t) such that
x0 (t) = a(u(x(t), t)) for t > 0
x(0) = ξ
Since u is assumed to be a classical solution, it is readily verified that
d
u(x(t), t) = ∂t u + x0 (t) · ∇u
dt
= ∂t u + a(u) · ∇u = ∂t u + ∇ · f (u) = 0
Therefore, u is constant along a characteristic curve and Γξ is a straight line since
x0 (t) = a(u(x(t), t)) = a(u(x(0), 0))
= a(u(ξ, 0)) = a(u0 (ξ)) = constant
In particular, x(t) is given by
x(t) = ξ + ta(u0 (ξ)) (7)
This important property may be used to construct classical solutions. If x and t are fixed and
ξ determined as a solution of (7), then
u(x, t) = u0 (ξ)
This procedure is the basis of the so-called method of characteristics. On the other hand,
this construction shows that the intersection of any two straight characteristic lines leads to
a contradiction in the definition of u(x, t). Thus, classical solutions can only exist up to the
first time t0 at which any two characteristics intersect.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 9

2.2. Weak solutions

Since, in general, classical solutions only exist for a finite time t0 , it is necessary to introduce
the notion of weak solutions that are well defined for times t > t0 .
Definition 1. (Weak solution) Let u0 ∈ L∞ (Rd ). Then, u is a weak solution of (6) if
u ∈ L∞ (Rd × R+ ) and (6) holds in the distributional sense, that is,
Z Z Z
(u∂t φ + f (u) · ∇φ) dt dx + u0 φ(x, 0) dx = 0 for all φ ∈ C01 (Rd × R+ ) (8)
Rd R+ Rd

Note that classical solutions are weak solutions and weak solutions that lie in C 1 (Rd × R+ )
satisfy (6) in the classical sense.

It can be shown (Kruzkov, 1970; Oleinik, 1963) that there always exists at least one weak
solution to (6) if the flux function f is at least Lipschitz continuous. Nevertheless, the class of
weak solutions is too large to ensure uniqueness of solutions. An important class of solutions are
piecewise classical solutions with discontinuities separating the smooth regions. The following
lemma gives a necessary and sufficient condition imposed on these discontinuities such that
the solution is a weak solution; see, for example, Godlewski and Raviart (1991) and Kröner
(1997). Later a simple example is given where infinitely many weak solutions exist.
Lemma 1. (Rankine–Hugoniot jump condition) Assume that the space-time domain Rd × R+
is separated by a smooth hypersurface S into two parts Ql and Qr . Furthermore, assume u is
a C 1 -function on Ql and Qr , respectively. Then, u is a weak solution of (6) if and only if the
following two conditions hold:

1. u is a classical solution in Ql and Qr .

2. u satisfies the Rankine–Hugoniot jump condition, that is,
[u]s = [f (u)] · n on S (9)
T
Here, (n, −s) denotes a unit normal vector for the (space-time) hypersurface S and [u]
denotes the jump in u across the hypersurface S.

In one space dimension (i.e., f = f is a scalar function), it may be assumed that S is

parameterized by (σ(t), t) such that s = σ 0 (t) and n = 1. The Rankine–Hugoniot jump
condition then reduces to
[f (u)]
s= on S (10)
[u]
Example 1. (Non-uniqueness of weak solutions) Consider the one-dimensional Burgers’
equation, f (u) = u2 /2, with Riemann data: u0 (x) = ul for x < 0 and u0 (x) = ur for x ≥ 0.
Then, for any a ≥ max(ul , −ur ) a function u given by

 ul , x < s1 t

−a, s1 t < x < 0
u(x, t) = (11)
 a,
 0 < x < s2 t
ur , s2 t < x

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
10 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

is a weak solution if s1 = (ul − a)/2 and s2 = (a + ur )/2. This is easily checked since u is
piecewise constant and satisfies the Rankine–Hugoniot jump condition. This elucidates a one-
parameter family of weak solutions. In fact, there is also a classical solution whenever ul ≤ ur .
In this case, the characteristics do not intersect and the method of characteristics yields the
classical solution (u ,
l x<ut l
u(x, t) = x/t, ul t < x < ur t (12)
ur , ur t < x
This solution is the unique classical solution but not the unique weak solution. Consequently,
additional conditions must be introduced in order to single out one solution within the class
of weak solutions. These additional conditions give rise to the notion of a unique entropy weak
solution.

2.3. Entropy weak solutions and vanishing viscosity

In order to introduce the notion of entropy weak solutions, it is useful to first demonstrate
that there is a class of additional conservation laws for any classical solution of (6). Let u be
a classical solution and η : R → R a smooth function. Multiplying (6a) by η 0 (u), one obtains
0 = η 0 (u)∂t u + η 0 (u)∇ · f (u) = ∂t η(u) + ∇ · F (u) (13)
where F is any primitive of η 0 f 0 . This reveals that for a classical solution u, the quantity η(u),
henceforth called an entropy function, is a conserved quantity.
Definition 2. (Entropy–entropy flux pair) Let η : R → R be a smooth convex function and
F : R → Rd a smooth function such that
F 0 = η0 f 0 (14)
in (13). Then (η, F ) is called an entropy–entropy flux pair or more simply an entropy pair for
the equation (6a).
Note 1. (Kruzkov entropies) The family of smooth convex entropies η may be equivalently
replaced by the nonsmooth family of the so-called Kruzkov entropies, that is, ηκ (u) ≡ |u − κ|
for all κ ∈ R. The associated entropy flux is then F κ (u) = (F (u) − F (κ))sg(u − κ), where sg
denotes the sign function (see e.g., Kröner, 1997).

Unfortunately, the relation (13) cannot be fulfilled for weak solutions in general, as it would lead
to additional jump conditions that would contradict the Rankine–Hugoniot jump condition
lemma. Rather, a weak solution may satisfy the relation (13) in the distributional sense with
inequality. To see that this concept of entropy effectively selects a unique, physically relevant
solution among all weak solutions, consider the viscosity-perturbed equation
∂t u + ∇ · f (u ) = ∆u (15)
with  > 0. For this parabolic problem, it may be assumed that a unique smooth solution u
exists. Multiplying by η 0 and rearranging terms yields the additional equation
∂t η(u ) + ∇ · F (u ) = ∆η(u ) − η 00 (u )|∇u|2

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 11

Furthermore, since η is assumed convex (η 00 ≥ 0), the following inequality is obtained

∂t η(u ) + ∇ · F (u ) ≤ ∆η(u )
Taking the limit  → 0 establishes (Málek, Nečas, Rokyta and Røužička, 1996) that u
converges toward some u a.e. in Rd × R+ where u is a weak solution of (6) and satisfies
the entropy condition
∂t η(u) + ∇ · F (u) ≤ 0 (16)
in the sense of distributions on Rd × R+ .

By this procedure, a unique weak solution has been identified as the limit of the
approximating sequence u . The obtained solution u is called the vanishing viscosity weak
solution of (6). Motivated by the entropy inequality (16) of the vanishing viscosity solution, it
is now possible to introduce the notion of entropy weak solutions. This notion is weak enough
for the existence and strong enough for the uniqueness of solutions to (6).
Definition 3. (Entropy weak solution) Let u be a weak solution of (6). Then, u is called an
entropy weak solution if u satisfies for all entropy pairs (η, F )
Z Z Z
(η(u)∂t φ + F (u) · ∇φ) dt dx + η(u0 )φ(x, 0) dx ≥ 0 (17)
Rd R+ Rd

for all φ ∈ C01 (Rd + +

× R , R ).

From the vanishing viscosity method, it is known that entropy weak solutions exist. The
following L1 contraction principle guarantees that entropy solutions are uniquely defined; see
Kruzkov (1970).
Theorem 1. (L1 -contraction principle) Let u and v be two entropy weak solutions of (6) with
respect to initial data u0 and v0 . Then, the following L1 -contraction principle holds
ku(·, t) − v(·, t)kL1 (Rd ) ≤ ku0 − v0 kL1 (Rd ) (18)
for almost every t > 0.

This principle demonstrates a continuous dependence of the solution on the initial data and
consequently the uniqueness of entropy weak solutions. Finally, note that an analog of the
Rankine–Hugoniot condition exists (with inequality) in terms of the entropy pair for all entropy
weak solutions
[η(u)]s ≥ [F (u)] · n on S (19)

2.4. Measure-valued or entropy process solutions

The numerical analysis of conservation laws is facilitated by an even weaker formulation of

solutions to (6). For instance, the convergence analysis of finite volume schemes makes it
necessary to introduce the so-called measure-valued or entropy process solutions; see DiPerna
(1985) and Eymard et al. (2000).

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
12 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

Definition 4. (Entropy process solution) A function µ(x, t, α) ∈ L∞ (Rd ×R+ ×(0, 1)) is called
an entropy process solution of (6) if u satisfies for all entropy pairs (η, F )
Z Z Z 1 Z
η(µ)∂t (φ + F (µ) · ∇φ) dα dt dx + η(u0 )φ(x, 0) dx ≥ 0
Rd R+ 0 Rd

for all φ ∈ C01 (Rd × R+ , R+ ).

The most important property of such entropy process solutions is the following uniqueness and
regularity result (see Eymard et al., 2000, Theorem 6.3).
Theorem 2. (Uniqueness of entropy process solutions) Let u0 ∈ L∞ (Rd ) and f ∈ C 1 (R).
The entropy process solution µ of problem (6) is unique. Moreover, there exists a function
u ∈ L∞ (Rd × R+ ) such that u(x, t) = µ(x, t, α) a.e. for (x, t, α) ∈ Rd × R+ × (0, 1) and u is
the unique entropy weak solution of (6).

3. Finite Volume Methods for Nonlinear Hyperbolic Conservation Laws

In the FVM for hyperbolic conservation laws, the computational domain, Ω ⊂ Rd , is first
tessellated into a collection of nonoverlapping control volumes that completely cover the
domain. Notationally, let T denote a tessellation of the domain Ω with control volumes K ∈ T
such that ∪K∈T K = Ω. Let hK denote a length scale associated with each control volume K,
for example, hK ≡ diam(K). For two distinct control volumes K and L in T , the intersection
is either an oriented edge (2-D) or face (3-D) denoted by σK,L with oriented normal nK,L
or else a set of measure at most d − 2. For simplicity, it is assumed throughout that control
volumes K are time invariant (unchanging in time). As mentioned in the introduction (see
(2)), for each control volume, an integral conservation law statement holds.
Definition 5. (Integral conservation law) An integral conservation law asserts that the rate
of change of the total amount of a substance with density u in a time-invariant control volume
K is equal to the total flux of the substance through the boundary ∂K
Z Z
d
u dx + f (u) · n ds = 0 (20)
dt K ∂K

This integral conservation law statement is readily obtained upon spatial integration of the
divergence equation (6a) in the region K and application of the divergence theorem. The
choice of control volume tessellation is flexible in the FVM. For example, Figure 1 depicts a
2-D triangle complex and two typical control volume tessellations (among many others) used
in the FVM. In the cell-centered FVM shown in Figure 1(a), the triangles themselves serve as
control volumes with solution unknowns (degrees of freedom) stored on a per triangle basis.
In the vertex-centered FVM shown in Figure 1(b), control volumes are formed as a geometric
dual to the triangle complex and solution unknowns stored on a per triangulation vertex basis.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 13

Storage location

Control volume

(a) (b)

Figure 1. Control volume variants used in the finite volume method: (a) cell-centered and (b) vertex-
centered control volume tessellation.

3.1. Finite volume discretization from cell averages via exact or approximate Riemann
problems

An important class of FVMs for hyperbolic problems comes from the introduction of control
volume cell averages Z
1
uK (t) ≡ u(x, t) dx, ∀K ∈ T (21)
|K| K
as unknowns in the numerical solution procedure. Recall that for simplicity, the control volume
K is assumed to be unchanging in time unless otherwise stated. The FVM can then be
interpreted as producing evolution equations for approximations of cell average unknowns
Z
d
|K| uK + f (u) · n ds = 0, ∀K ∈ T (22)
dt ∂K

if the flux integral can somehow be evaluated. In the context of the gas dynamic equations,
Godunov (1959) pursued an interpretation of the cell averages as a piecewise constant
representation of the numerical solution. This interpretation readily applies to scalar hyperbolic
conservation laws as well. The piecewise constant representation renders the numerical solution
multivalued at control volume interfaces. It then becomes unclear how the flux integral
appearing in (22) should be interpreted so that discrete conservation is maintained. Discrete
conservation in the finite volume method demands that for two control volumes K and L that
share an interface σK,L , the amount of substance that fluxes out of K through σK,L must
exactly equal the amount of substance that fluxes into L through σL,K . Godunov solved the
problems of multivalued solution states and discrete conservation by finding a single solution
state u∗ (uK , uL ) with symmetry u∗ (uK , uL ) = u∗ (uL , uK ) so that discrete conservation is
automatically satisfied
Z Z
∗
f (u (uK , uL )) · n ds = − f (u∗ (uL , uK )) · n ds (23)
σK,L σL,K

In Godunov’s original work, this single solution state was obtained by solving the one-
dimensional Riemann problem of gas dynamics. For the scalar hyperbolic problem (6) this
amounts to finding the state u∗ (uK , uL ) equal to the Riemann state uR (uK , uL ; nK,L )
described in the following definition.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
14 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

Definition 6. (Scalar Riemann state) Given two solution states u and v and normal vector
n, let hn (u) ≡ f (u) · n denote the flux in (6) projected in the direction n. Calculate
the Riemann state uR (u, v; n) by first solving the one-dimensional Riemann problem for
w(ξ, τ ; u, v) : R × R+ 7→ R
∂ ∂
w+ hn (w) = 0 (24)
∂τ ∂ξ
subject to initial data 
u if ξ < 0
w(ξ, 0; u, v) = (25)
v if ξ > 0
Next, determine the desired Riemann state from the Riemann problem solution w(ξ, τ ; u, v) by
evaluating it at ξ = 0 for any positive time, that is,
uR (u, v; n) = w(0, τ > 0; u, v) (26)
Owing to the self-similarity of the Riemann problem in the single parameter ξ/τ , this state
value is independent of time for τ > 0.

More generally, given two solution states u and v that share an interface with normal n,
variants of the Godunov approach are obtained by supplanting the true flux at this interface
by a numerical flux function g(u, v; n) : R × R 7→ R, a Lipschitz continuous function of the two
interface states u and v for a given normal n. In higher space dimensions, the flux integral
appearing in (20) is then approximated by
Z X
f (u) · n ds ≈ g(uK , uL ; nK,L )|σK,L | (27)
∂K σK,L ⊂∂K

The numerical flux is assumed to satisfy the properties:

• Conservation: This property ensures that fluxes from adjacent control volumes sharing
a mutual interface exactly cancel when summed. This is achieved if the numerical flux
satisfies the identity
g(u, v; n) = −g(v, u; −n) (28a)
• Consistency: Consistency is obtained if the numerical flux with identical state arguments
reduces to the true flux of that same state projected in the direction n, that is,
g(u, u; n) = f (u) · n (28b)

Combining (22) and (27) yields the following semi-discrete FVM shown here for a time invariant
(fixed) mesh.
Definition 7. (Semi-discrete FVM for hyperbolic problems) The semi-discrete finite volume
approximation of the hyperbolic conservation law problem (6) using cell averaged solution
unknowns (21) with continuous in time derivatives (22) and numerical flux function quadrature
in space (27) for time-invariant control volumes K is given by
d 1 X
uK + g(uK , uL ; nK,L ) |σK,L | = 0 ∀K ∈ T (29)
dt |K|
σK,L ⊂∂K

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 15

This system of ordinary differential equations can be marched forward in time using a variety
of explicit and implicit time integration methods; see Section 4.4.

Let unK denote a numerical approximation of the cell average solution in the control volume
K at time tn ≡ n∆t. A particularly simple time integration method is the forward Euler
scheme
d un+1 − unK
uK ≈ K (30)
dt ∆t
which can be viewed as resulting from a piecewise constant representation of the solution in
time. Using the explicit time advancement formula together with the flux quadrature (27)
yields the following fully discrete finite volume formulation:
Definition 8. (Fully discrete FVM for hyperbolic problems) The fully discrete finite volume
approximation of the hyperbolic conservation law problem (6) using Euler explicit time
advancement (30) and numerical flux function quadrature in space (27) for time-invariant
control volumes K is given by
∆t X
un+1
K = u n
K − g(unK , unL ; nK,L ) |σK,L | ∀K ∈ T (31)
|K|
σK,L ⊂∂K

Unfortunately, the numerical flux conditions (28a) and (28b) are insufficient to guarantee
that stable numerical solutions will be produced that converge to entropy satisfying weak
solutions (17). Consequently, additional numerical flux restrictions are necessary. In the
following section, an early result concerning monotone schemes is presented that addresses the
question of convergence to entropy weak solutions. These results motivate the construction
of monotone flux and E-flux functions that guarantee a local maximum principle and global
maximum norm stability.

3.2. Monotone fluxes and E-flux functions

An early result by Harten, Hyman and Lax (1976) addresses the question of convergence of
the fully discrete one-dimensional finite volume scheme to weak entropy satisfying solutions.
For ease of notation in describing stencil operators, the shorthand notation unj ≡ unKj has been
adopted.
Theorem 3. (Monotone schemes and weak solutions) Consider a 1-D finite volume
discretization of (6) with 2k + 1 stencil on a uniformly spaced mesh in both time and space
with corresponding mesh spacing parameters ∆t and ∆x
un+1
j = Hj (unj+k , . . . , unj , . . . , unj−k )
∆t n
= unj − (g n
− gj−1/2 ) (32)
∆x j+1/2
and consistent numerical flux of the form
gj+1/2 = g(uj+k , . . . , uj+1 , uj , . . . , uj−k+1 ) (33)

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
16 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

that is monotone in the sense

∂Hj
≥ 0, ∀ |l| ≤ k (34)
∂uj+l
Assume that unj converges boundedly almost everywhere to some function u(x, t); then as ∆t
and ∆x tend to zero with ∆t/∆x = constant, this limit function u(x, t) is an entropy satisfying
weak solution of (6).

Note that this theorem assumes convergence in the limit, which was later proved to be the
case in multidimensions by Crandall and Majda (1980).

The monotone scheme conditions (34) may be achieved by the use of Lipschitz continuous
monotone fluxes satisfying the following conditions
∂gj+1/2
≥0 if l = j (35a)
∂ul
∂gj+1/2
≤0 if l 6= j (35b)
∂ul
together with a CFL (Courant–Friedrichs–Lewy) like condition
 
∆t ∂gj+1/2 ∂gj−1/2
1− − ≥0 (36)
∆x ∂uj ∂uj

These monotone flux functions are considered further in the following section for the semi-
discrete FVM (29) and the fully discrete FVM (31). Example monotone flux functions are
then presented.

3.2.1. Monotone flux functions. The numerical flux function (33) accommodates large mesh
stencils unlike the numerical flux functions considered in Section 3.1 that only utilize two
states. Given two states u and v that share an interface with normal n, the flux monotonicity
conditions (35a) and (35b) reduce to

∂g(u, v; n)
≥0 (37a)
∂u
∂g(u, v; n)
≤0 (37b)
∂v
Some examples of two state monotone fluxes for the hyperbolic problem (6a) with convex
flux, f 00 > 0, are

• (Riemann flux), see Section 3.1


 min f (w) · n if u < v ,
w∈[u,v]
g R (u, v; n) = (38)
 max f (w) · n if u > v .
w∈[u,v]

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 17

• (“local” Lax–Friedrichs flux), see Shu and Osher (1988)

1 1
g LF (u, v; n) = (f (u) + f (v)) · n − sup |f 0 (w) · n|(v − u) (39)
2 2 w∈[u,v]
• (Upwind flux with sonic point modification), see Godlewski and Raviart (1991)

 f (u) · n if f 0 (w) · n ≥ 0, ∀w ∈ [u, v] ,
g upwind
(u, v; n) = f (v) · n if f 0 (w) · n ≤ 0, ∀w ∈ [u, v] , (40)
 LF
g (u, v; n) otherwise .
Note that this flux is not monotone in the usual sense introduced above, since the function
g upwind is not continuous at the sonic point. However, the notion of monotone flux can
be extended to noncontinuous functions in the following way:
Definition 9 (Monotone numerical flux.) Eymard et al., 2016 Let f be a locally
Lipschitz continuous function on R, and let A, B ∈ R such that A ≤ B. A function
g : R2 → R is said to be a monotone numerical flux for f on [A, B] if it satisfies the
following assumptions:
– There exists Lg ∈ R+ such that

|g(a, b) − f (a)| ≤ Lg |a − b|
∀(a, b) ∈ (A, B)2 , (41)
|g(b, a) − f (a)| ≤ Lg |a − b|
– g(s, s) = f (s), for all s ∈ [A, B],
– the function g : (a, b) 7→ g(a, b), from [A, B]2 to R, is nondecreasing with respect to
a and nonincreasing with respect to b.
The flux g upwind (·, ·; n) defined by (40) can be shown to be monotone for the function
f = f ·n. It is then easy to show that proofs of convergence that are written for Lipschitz
continuous monotone fluxes are also valid for monotone fluxes in the sense of Definition
41.

Monotone flux functions have played an enormously important role in the development of
FVMs for hyperbolic problems but can sometimes be difficult to construct or evaluate. In the
following section, a slightly weaker condition called the E-flux condition is presented. Monotone
fluxes satisfy this E-flux condition.

3.2.2. E-flux functions. Another class of numerical fluxes arising frequently in analysis and
practical implementations was introduced by Osher (1984). These fluxes are called E-fluxes due
to the relationship to Oleinik’s well-known E-condition, which characterizes entropy satisfying
discontinuities. E-fluxes satisfy the inequality
g E (u, v; n) − f (w) · n
≤ 0, ∀w ∈ [u, v] (42)
v−u
E-fluxes can be characterized by their relationship to the Riemann flux. Specifically, E-fluxes
are those fluxes such that
g E (u, v; n) ≤ g R (u, v; n) if v < u (43a)
E R
g (u, v; n) ≥ g (u, v; n) if v > u (43b)

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
18 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

Viewed another way, note that any numerical flux can be written in the form
1 1
g(u, v; n) = (f (u) + f (v)) · n − Q(u, v; n)(v − u) (44)
2 2
where Q(·) denotes a viscosity for the scheme. When written in this form, E-fluxes are those
fluxes that contribute at least as much viscosity as the Riemann flux, that is,

QR (u, v; n) ≤ QE (u, v; n) (45)

The most prominent E-flux is the Enquist–Osher flux

Z v
1 1
g EO (u, v; n) = (f (u) + f (v)) · n − |f 0 (w) · n| dw (46)
2 2 u

which was generalized to systems of conservation laws in Osher and Solomon (1982), see Sect.
6.1.2.

Monotone flux and E-flux functions provide the needed properties for proving discrete
maximum principles in the FVM. These results are presented in the following section. Section
6.1.3 shows that the E-flux condition (42) has a natural extension to systems of symmetrizable
hyperbolic conservation laws and plays an important role in proving energy/entropy stability
of those systems using the FVM.

3.2.3. Discrete maximum principles and stability using monotone flux and E-flux functions.
A compelling motivation for the use of monotone flux and E-flux functions in the FVMs (29)
and (31) is to obtain discrete maximum principles in the resulting numerical solutions. A
standard analysis technique is to first construct local discrete maximum principles that can
then be applied successively to obtain global maximum principles and maximum norm stability
results.

The following two lemmas concern the boundedness of local extrema and a discrete maximum
principle for FVMs that can be written in nonnegative coefficient form. The first lemma
addresses the evolution of local extrema using the semi-discrete finite volume discretization
(29) when rewritten in nonnegative coefficient form.
Lemma 2. (Semi-discrete LED property) The semi-discrete scheme for each K ∈ T
duK 1 X
= CK,L (uh )(uL − uK ) (47)
dt |K|
σK,L ⊂∂K

with uh ≡ {uK1 , uK2 , . . .} is local extremum diminishing (LED), that is, local maxima are
nonincreasing and local minima are nondecreasing, if

CK,L (uh ) ≥ 0, ∀σK,L ⊂ ∂K (48)

The second lemma addresses a local maximum principle using the fully discrete finite volume
discretization (31) when rewritten in nonnegative coefficient form.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 19

Lemma 3. (Fully discrete local maximum principle) The fully discrete scheme for the time
slab increment [tn , tn+1 ] and each K ∈ T
∆t X
un+1
K = unK + CK,L (unh )(unL − unK ) (49)
|K|
σK,L ⊂∂K

with uh ≡ {uK1 , uK2 , . . .} exhibits a local discrete maximum principle for each n = 0, 1, 2, . . .
min (unL , unK ) ≤ un+1
K ≤ max (unL , unK ) (50)
σK,L ⊂∂K σK,L ⊂∂K

if
CK,L (unh ) ≥ 0, ∀σK,L ⊂ ∂K (51)
and ∆t is chosen such that the CFL-like condition is satisfied
∆t X
1− CK,L (unh ) ≥ 0 (52)
|K|
σK,L ⊂∂K

The results of lemmas 2 and 3 require showing that the semi-discrete finite volume scheme
(29) and fully discrete finite volume scheme (31) can be placed in nonnegative coefficient form.
This can be easily shown when monotone flux and E-flux functions are used by rewriting the
flux divergence terms in the following equivalent forms
X X
g(uK , uL ; nK,L ) |σK,L | = CK,L (uh ) (uL − uK ) (53a)
σK,L ⊂∂K σK,L ⊂∂K
X g(uK , uL ; nK,L ) − f (uK ) · nK,L
= |σK,L | (uL − uK ) (53b)
uL − uK
σK,L ⊂∂K
X ∂g(uK , ũK,L )
= |σK,L | (uL − uK ) (53c)
∂uK
σK,L ⊂∂K

for appropriately chosen ũK,L ∈ [uK , uL ]. Nonnegativity of the coefficients CK,L in the (53a)
right-hand-side summands is achieved in the (53b) right-hand-side summands whenever the
numerical flux satisfies the E-flux condition (42) and similarly in (53c) whenever the numerical
flux satisfies the monotonicity conditions (37). These results are then used to make a statement
concerning stability in a maximum norm.
Theorem 4. (L∞ -stability) The fully discrete finite volume scheme (31) utilizing either the
monotone flux of Section 3.2.1 or the E-flux functions of Section 3.2.2 subject to a local CFL-
like condition as given in lemma 3 for each time slab increment [tn , tn+1 ] is L∞ -stable in the
following sense
inf u0 (x) ≤ unK ≤ sup u0 (x) (54)
x∈Rd x∈Rd

for all K ∈ T and time step tn , n = 0, 1, . . ..

The LED and local maximum principles discussed in lemmas 2 and 3 preclude the
introduction of spurious extrema and O(1) Gibbs-like oscillations that occur near solution

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
20 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

discontinuities computed using many numerical methods (even in the presence of grid
refinement). For this reason, FVMs for hyperbolic conservation laws that possess these
LED and discrete maximum principle properties have proved highly successful in practical
calculations.

3.3. Stability, convergence, and error estimates

Several stability results have been presented in Section 3.2.3 that originate from discrete
maximum principle analysis and are straightforwardly stated in the multidimensional setting
and on general unstructured meshes. In presenting results concerning convergence and error
estimates, a notable difference arises between one and several space dimensions. This is due to
the lack of a BV bound on the approximate solution in the multidimensional setting, except
in the case of Cartesian meshes.

3.3.1. Convergence results. The L∞ -stability bound (54) is an essential ingredient in the
proof of convergence of the fully discrete finite volume scheme (31). This bound permits the
extraction of a subsequence that converges against some limit in the L∞ weak-starred sense.
The primary task that then remains is to identify this limit with the unique solution of the
problem. To this end stronger estimates are needed, both for convergence and for deriving
convergence rates.

Let BV denote the space of functions with bounded variation, that is,
n o
BV = g ∈ L1 (Rd ) |g|BV < ∞
with Z
|g|BV = sup g∇ · ϕ dx
ϕ∈Cc1 (Rd )d Rd
kϕk∞ ≤1

From the theory of scalar conservation laws, it is known that, provided the initial data
is in BV, the solution remains in BV for all times. Therefore, it is desirable to have an
analog of this property for the approximate solution as well. Unfortunately, it is known
that finite volume schemes for hyperbolic problems on unstructured meshes may be non-
total variation diminishing (TVD), and in fact, the counterexample of Desprès, 2004 shows
that the approximate solutions of an upwind finite volume scheme on triangles may blow up
in the BV norm. However, the approximate finite volume solutions can be shown to fulfill a
weaker estimate, called a weak BV estimate; see Champier and Gallouët (1992), Champier
et al. (1993), Vila (1994), Cockburn, Coquel and Lefloch (1994), and Eymard et al. (1998).
Theorem 5. (Weak BV estimate) Let T be a regular triangulation, and let J be a uniform
partition of [0, τ ], for example, ∆tn ≡ ∆t. Assume that there exists some α > 0 such that
αh2 ≤ |K|, α|∂K| ≤ h. Let g be a monotone flux function, that is, a Lipschitz-continuous
function satisfying (35). Assume the following CFL-like condition for a given ξ ∈ (0, 1)
(1 − ξ)α2 h
∆t ≤
Lg

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 21

where Lg is the Lipschitz constant of the numerical flux function g. Furthermore, let u0 ∈
L∞ (Rd ) ∩ BV (Rd ) ∩ L2 (Rd ). Then, the numerical solution of the fully discrete solution (31)
fulfills the following estimate
X X p √
∆t χKL h|unK − unL |QKL (unK , unL ) ≤ C T |BR+h (0)| h (55)
n σKL ∈E

where C only depends on α, Lg , ξ, and the initial function u0 . In this formula QKL is defined
as
2gKL (u, v) − gKL (u, u) − gKL (v, v)
QKL (u, v) ≡
u−v
and χKL denotes the discrete cutoff function on BR (0) ⊂ Rd , that is,
n
χKL = 1 if (K ∪ L) ∩ BR (0) 6= ∅ ,
0 else

Note that in the case of a strong BV estimate (obtained√ in the case of Cartesian meshes),
the right-hand side of (55) would be O(h) instead of O( h). In the case of linear advection
problems, a somewhat stronger estimate was proven in Desprès and Lagoutière, 2010, which
implies the above weak BV estimate.

Another important property of monotone finite volume schemes is that they preserve the
L1 -contraction property (Theorem 1).
Theorem 6. (L1 -contraction property and Lipschitz estimate in time) Let uh , vh ∈ Vh0 be the
approximate monotone finite volume solutions corresponding to initial data u0 , v0 assuming
that the CFL-like condition for stability has been fulfilled. Then the following discrete L1 -
contraction property holds

kuh (·, t + τ ) − vh (·, t + τ )kL1 (Rd ) ≤ kuh (·, t) − vh (·, t)kL1 (Rd )

Furthermore, a discrete Lipschitz estimate in time is obtained

X X X
|K||un+1
K − unK | ≤ Lg ∆tn |σKL ||u0K − u0L |
K∈T K∈T σKL ∈EK

The principal ingredients of the convergence theory for scalar nonlinear conservation laws
are compactness of the family of approximate solutions and the passage to the limit within
the entropy inequality (17). In dealing with nonlinear equations, strong compactness is
needed in order to pass to the limit in (17). In one space dimension or in the case of
multidimensional Cartesian meshes, due to the BV estimate and the selection principle of
Helly, strong compactness is ensured and the passage to the limit is summarized in the well-
known Lax–Wendroff theorem; see Lax and Wendroff (1960).
Theorem 7. (Lax–Wendroff theorem) Let (um )m∈N be a sequence of discrete solutions defined
by the finite volume scheme in one space dimension with respect to initial data u0 . Assume that
(um )m∈N is uniformly bounded with respect to m in L∞ and um converges almost everywhere
in R × R+ to some function u. Then u is the uniquely defined entropy weak solution of (6).

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
22 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

With the lack of a BV estimate for the approximate solution in multiple space dimensions,
one cannot expect a passage to the limit of the nonlinear terms in the entropy inequality in the
classical sense. Nevertheless, the weak compactness obtained by the L∞ -estimate is enough
to obtain a measure-valued or entropy process solution in the limit. The convergence of a
subsequence of approximate solutions to the entropy weak solution is then obtained thanks to
the uniqueness result of Eymard et al., 1995. Note that this approach has been adapted for a
constrained hyperbolic equation in Andreianov et al., 2011.

The key theorem for this convergence result is the following compactness theorem of Eymard
et al. (2000), which is a convenient rewriting of a fundamental theorem on Young measures
for PDEs; see Tartar (1979) and Ball (1989).
Theorem 8. Let (um )m∈N be a family of bounded functions in L∞ (Rd ). Then, there exists a
subsequence (um )m∈N , and a function u ∈ L∞ (Rd × (0, 1)) such that for all functions g ∈ C(R)
the weak-? limit of g(um ) exists and
Z Z 1Z
lim g(um (x))φ(x) dx = g(u(x, α))φ(x) dx dα, for all φ ∈ L1 (Rd ) (56)
m→∞ Rd 0 Rd

In order to prove the convergence of a FVM, it now remains to be shown that the residual of the
entropy inequality (17) for the approximate solution uh tends to zero if h and ∆t tend to zero.
Before presenting this estimate for the finite volume approximation, a general convergence
theorem is given, which can be viewed as a generalization of the classical Lax–Wendroff result;
see Eymard et al. (2000).
Theorem 9. (Sufficient condition for convergence) Let u0 ∈ L∞ (Rd ) and f ∈ C 1 (R).
Further, let (um )m∈N be any family of uniformly bounded functions in L∞ (Rd × R+ ) that
satisfies the following estimate for the residual of the entropy inequality using the class of
Kruzkov entropy pairs (ηκ , F κ ) (see Note 1).
Z Z Z
(ηκ (um )∂t φ + F κ (um ) · ∇φ) dt dx + ηκ (u0 )φ(x, 0) dx ≥ −R(κ, um , φ) (57)
Rd R+ Rd

for all κ ∈ R and φ ∈ C01 (Rd + +

× R , R ) where the residual R(κ, um , φ) tends to zero for
m → ∞ uniformly in κ. Then, um converges strongly to the unique entropy weak solution of
(6) in Lploc (Rd × R+ ) for all p ∈ [1, ∞).
Theorem 10. (Estimate on the residual of the entropy inequality) Let (um )m∈N be a sequence
of monotone finite volume approximations satisfying a local CFL-like condition as given in (52)
such that h, ∆t tend to zero for m → ∞. Then, there exist measures µm ∈ M(Rd × R+ ) and
νm ∈ M(Rd ) such that the residual R(κ, um , φ) of the entropy inequality is estimated by
Z Z Z
R(κ, um , φ) ≤ (|∂t φ(x, t)| + |∇φ(x, t)|) dµm (x, t) + φ(x, 0) dνm (x)
Rd R+ Rd

for all κ ∈ R and φ ∈ C01 (Rd + +

× R , R ). The measures µm and νm satisfy the following
properties:

1. For all compact subsets Ω ⊂ Rd × R+ , limm→∞ µm (Ω) = 0.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 23

R
2. For all g ∈ C0 (Rd ) the measure νm is given by hνm , gi = Rd
g(x)|u0 (x) − um (x, 0)| dx.

These theorems are sufficient for establishing convergence of monotone finite volume schemes.
Corollary 1. (Convergence theorem) Let (um )m∈N be a sequence of monotone finite volume
approximations satisfying the assumptions of Theorem 10. Then, um converges strongly to the
unique entropy weak solution of (6) in Lploc (Rd × R+ ) for all p ∈ [1, ∞).

Convergence of higher order finite volume schemes can also be proved within the given
framework as long as they are L∞ -stable and allow for an estimate on the entropy residual
in the sense of Theorem 10; for details see Kröner, Noelle and Rokyta (1995) and Chainais-
Hillairet (2000).

3.3.2. Error estimates and convergence rates. There are two primary approaches taken
to obtaining error estimates for approximations of scalar nonlinear conservation laws. One
approach is based on the ideas of Oleinik and is applicable only in one space dimension; see
Oleinik (1963) and Tadmor (1991). The second approach, which is widely used in the numerical
analysis of conservation laws, is based on the doubling of variables technique of Kruzkov; see
Kruzkov (1970) and Kuznetsov (1976). In essence, this technique enables one to estimate the
error between the exact and approximate solution of a conservation law in terms of the entropy
residual R(κ, um , Φ) introduced in (57). Thus, an a posteriori error estimate is obtained.
Using a priori estimates of the approximate solution (see Section 3.2.3, and Theorems 5,
6), a convergence rate or an a priori error estimate is then obtained. The following theorem
gives a fundamental error estimate for conservation laws independent of the particular finite
volume scheme; see Eymard et al. (1998), Eymard et al. (2000), Chainais-Hillairet (1999), and
Kröner and Ohlberger (2000).
Theorem 11. (Fundamental error estimate) Let u0 ∈ BV (Rd ) and let u be an entropy weak
solution of (6). Furthermore, let v ∈ L∞ (Rd × R+ ) be a solution of the following entropy
inequalities with residual term R:
Z Z Z
(ηκ (v)∂t φ + Fκ (v) · ∇φ) dx + ηκ (u0 )φ(·, 0)dx ≥ −R(κ, v, φ) (58)
Rd R+ Rd

for all κ ∈ R and φ ∈ C01 (Rd × R+ , R+ ). Suppose that there exist measures µv ∈ M(Rd × R+ )
and νv ∈ M(Rd ) such that R(κ, v, φ) can be estimated independently of κ by
R(κ, v, φ) ≤ h|∂t φ| + |∇φ|, µv i + h|φ(·, 0)|, νv i (59)
d +
Let G ⊂⊂ R × R , ω ≡ Lip(f ), and choose T , R and x0 such that T ∈]0, (R/ω)[ and G lies
within its cone of dependence D0 , that is, G ⊂ D0 where Dδ is given as
[
Dδ := BR−ωt+δ (x0 ) × {t} (60)
0≤t≤T

Then, there exists a δ ≥ 0 and positive constants C1 , C2 such that u, v satisfy the following
error estimate
 p 
ku − vkL1 (G) ≤ T νv (BR+δ (x0 )) + C1 µv (Dδ ) + C2 µv (Dδ ) (61)

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
24 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

This estimate can be used either as an a posteriori control of the error, as the right-hand side
of the estimate (61) only depends on v, or it can be used as an a priori error bound if one is
able to estimate further the measures µv and νv using some a priori bounds on v. Finally, note
that comparable estimates to (61) are obtainable in an L∞ (0, T ; L1 (Rd ))-norm; see Cockburn
and Gau (1995) and Bouchut and Perthame (1998).

3.3.3. A posteriori error estimate. Based on the fundamental error estimate in Theorem 11,
the following theorem states an a posteriori error estimate that can be used to design self-
adaptive variants of finite volume schemes; see Kröner and Ohlberger (2000) and the review
article Ohlberger (2009).
Theorem 12. (A posteriori error estimate) Assume the conditions and notations as in
Theorem 11. Let v = uh be a numerical approximation to (6) obtained from a monotone finite
volume scheme that satisfies a local CFL-like condition as given in (52). Then the following
error estimate holds
Z
√

|u − uh |dx ≤ T ku0 − uh (·, 0)kL1 (BR+h (x0 )) + C1 η + C2 η (62)
G

where
X X
n+1
η ≡ |uK − unK |∆tn hdK
n∈I0 K∈M (tn )
X X
+ 2 ∆tn (∆tn + hKL ) × QKL (unK , unL )|unK − unL | (63)
n∈I0 σKL ∈E(tn )

with
2gKL (u, v)−gKL (u, u) − gKL (v, v)
QKL (u, v) ≡ (64)
u−v
and the index sets I0 , M (t), E(t) are given by
  
R+δ
I0 ≡ n | 0 ≤ tn ≤ min ,T ,
ω
M (t) ≡ {K | there exists x ∈ K such that (x, t) ∈ DR+δ },
E(t) ≡ {σKL | there exists x ∈ K ∪ L such that (x, t) ∈ DR+δ }

Furthermore, the constants C1 , C2 only depend on T , ω, ku0 kBV and ku0 kL∞ .

Note that this a posteriori error estimate is local, since the error on a compact set K is
estimated by discrete quantities that are supported in the cone of dependence DR+δ . Similar
results have been obtained for conservation laws on bounded domains Ohlberger and Vovelle
(2006) and for higher order discontinuous Galerkin generalizations Dedner et al. (2007).

3.3.4. A priori error estimate. Using the weak BV estimate (Theorem 5) and the Lipschitz
estimate in time (Theorem 6), the right-hand side of the a posteriori error estimate (62) can

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 25

be further estimated. This yields an a priori error estimate as stated in the following theorem;
for details, see Cockburn and Gremaud (1996, 1997, 1998a), Eymard et al. (1998), Chainais-
Hillairet (1999), and Eymard et al. (2000).
Theorem 13. (A priori error estimate) Assume the conditions and notations as in
Theorem 11 and let v = uh be the approximation to (6) given by a monotone finite volume
scheme that satisfies a local CFL-like condition as given in (52). Then there exists C ≥ 0 such
that Z
|u − uh | dx ≤ Ch1/4 (65)
G

Moreover, in the one-dimensional case or Cartesian multidimensional case, the optimal

convergence rate of h1/2 is obtained.

4. Higher Order Accurate Finite Volume Methods for Hyperbolic Problems

The FVMs for hyperbolic problems described in Section 3 are a powerful methodology that
has been used extensively in practical calculations. The error estimates of Section 3.3 predict
a deterioration from first-order accuracy for general problems and this deterioration can be
observed in some practical calculations. Nevertheless, for smooth solutions on smoothly varying
meshes, first-order accuracy is routinely observed. For problems that require very high solution
accuracy, first-order accuracy may still require a prohibitively large number of finite volume
cells. This has motivated the development of higher order accurate FVMs that require far
fewer finite volume cells while still providing a nonoscillatory resolution of discontinuities and
steep solution gradients. Unfortunately, Godunov (1959) has shown that all linear schemes that
preserve solution monotonicity are at most first-order accurate. Thus, higher order accurate
methods must utilize essential nonlinearity so that non-oscillatory resolution of discontinuities
and high-order accuracy away from discontinuities are simultaneously attained. These methods
are described below using both structured and unstructured meshes.

4.1. Higher order accurate finite volume methods for hyperbolic problems in one dimension

A significant step forward in the generalization of FVMs to higher order accuracy for hyperbolic
problems is due to van Leer (1979). In this work, van Leer generalized Godunov’s method by
employing linear solution reconstruction in each cell from given cell averages (Figures 2b).
This reconstruction yields a piecewise linear representation of the solution while still retaining
cell averages equal to the given data. The concept extends to higher order polynomials
such as depicted in Figure 2(c) using piecewise quadratic reconstruction; see Colella and
Woodward (1984). A close examination of Figure 2(b) reveals that this particular piecewise
linear reconstruction contains more local extrema than the underlying data. These spurious
extrema (oscillations) can seriously degrade the accuracy of a numerical solution in both space
and time. One approach for removing these spurious extrema during the reconstruction process

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
26 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

Figure 2. Piecewise polynomial approximation used in the finite volume method (a) cell averaging
of analytic data, (b) piecewise linear reconstruction from cell averages, and (c) piecewise quadratic
reconstruction from cell averages.

is to alter the slope of the reconstructed data in some or all of the cells. Another approach
is based on the idea of finding the smoothest piecewise polynomial reconstructions among
many possible candidates. The underlying mathematical principles and techniques used in
the development of higher order accurate finite volume approximations that control spurious
extrema are important technical contributions that are discussed in the remainder of this
section.

4.1.1. Total variation diminishing (TVD) finite volume methods. In considering the scalar
nonlinear conservation law (6) in one space dimension and time, Lax (1973) made the following
basic observation:

the total increasing and decreasing variations of a differentiable solution between

any pair of characteristics are conserved.

Furthermore, in the presence of shock wave discontinuities, information is lost and the total
variation decreases. For the 1-D nonlinear conservation law with compactly supported or
periodic solution data u(x, t), integrating along the constant time spatial coordinate at times
t1 and t2 yields
Z ∞ Z ∞
| du(x, t2 )| ≤ | du(x, t1 )|, t2 ≥ t1 (66)
−∞ −∞

This total variation property for scalar conservation laws motivated Harten (1983a) to consider
the discrete total variation in the design of numerical methods. Using a one-dimensional mesh
(assumed periodic here) with solution data uh ≡ {u1 , u2 , . . . , uN } at cell centroids of the
intervals [xj−1/2 , xj+1/2 ], j = 1, . . . , N , Harten considered the discrete total variations
X
TV(uh ) ≡ |∆j+1/2 uh |, ∆j+1/2 uh ≡ uj+1 − uj
j

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 27

and imposed the discrete total variation nonincreasing (TVNI) bound counterpart to (66)
TV(un+1
h ) ≤ TV(unh ) (67)
in the design of numerical discretizations for nonlinear conservation laws. A number of
theoretical results relating TVNI schemes and monotone schemes follow from analysis.
Theorem 14. (TVNI and monotone scheme properties, Harten, 1983a) (i) Monotone
schemes are TVNI. (ii) TVNI schemes are monotonicity preserving, that is, the number of
solution extrema is preserved in time.

Using the notion of discrete total variation, Harten (1983a) then constructed sufficient algebraic
conditions for achieving the TVNI inequality (67) in a fully discrete numerical method.
Theorem 15. (Harten’s explicit TVD criteria) The fully discrete 1-D discretization


un+1
j = unj + ∆t Cj+1/2 (unh )∆j+1/2 unh + Dj+1/2 (unh )∆j−1/2 unh , j = 1, . . . , N (68)
is TVNI if for each j
Cj+1/2 ≥ 0 (69a)
Dj+1/2 ≤ 0 (69b)


1 − ∆t Cj−1/2 − Dj+1/2 ≥ 0 (69c)

Note that although the inequality constraints (69) in Theorem 15 insure that the total variation
is nonincreasing, these conditions are often referred to as total variation diminishing (TVD)
conditions. Also note that inequality (69c) implies a CFL-like time step restriction that may
be different from the time step required for stability of the numerical method. The TVD
conditions are easily generalized to wider support stencils written in incremental form; see, for
example, Jameson and Lax (1986) and their corrected result in Jameson and Lax (1987).
Theorem 16. (Generalized explicit TVD criteria) The fully discrete explicit 1-D scheme
k−1
X (l)
un+1
j = unj + ∆t Cj+1/2 (unh )∆j+l+1/2 unh , j = 1, . . . , N (70)
l=−k

with integer stencil width parameter k > 0 is TVNI if for each j

(k−1)
Cj+1/2 ≥ 0 (71a)
(−k)
Cj+1/2 ≤ 0 (71b)
(l−1) (l)
Cj+1/2 − Cj−1/2 ≥ 0, −k + 1 ≤ l ≤ k − 1, l 6= 0 (71c)
 
(0) (−1)
1 − ∆t Cj−1/2 − Cj+1/2 ≥ 0 (71d)

While this simple Euler explicit time integration scheme may seem inadequate for applications
requiring true high-order space–time accuracy, special attention and analysis is given to this

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
28 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

fully discrete form because it serves as a fundamental building block for an important class
of high-order accurate Runge–Kutta time integration techniques discussed in Section 4.4 that,
by construction, inherit TVD (and later maximum principle) properties of the fully discrete
scheme (70).

The extension to implicit methods follows immediately upon rewriting the implicit scheme
in terms of the solution spatial increments ∆j+l+1/2 uh and imposing sufficient algebraic
conditions such that the implicit matrix acting on spatial increments has a nonnegative inverse.
Theorem 17. (Generalized implicit TVD criteria) The fully discrete implicit 1-D scheme
k−1
X (l)
un+1
j − ∆t Cj+1/2 (un+1
h )∆j+l+1/2 un+1
h = unj , j = 1, . . . , N (72)
l=−k

with integer stencil width parameter k > 0 is TVNI if for each j

(k−1)
Cj+1/2 ≥ 0 (73a)
(−k)
Cj+1/2 ≤ 0 (73b)
(l−1) (l)
Cj+1/2 − Cj−1/2 ≥ 0, −k + 1 ≤ l ≤ k − 1, l 6= 0 (73c)

Theorems 16 and 17 provide sufficient conditions for nonincreasing total variation of explicit
(70) or implicit (72) numerical schemes written in incremental form. These incremental
forms do not imply discrete conservation unless additional constraints are imposed on the
discretizations. A sufficient condition for discrete conservation of the discretizations (70) and
(72) is that these discretizations can be written in a finite volume flux balance form
k−1
X (l)
gj+1/2 − gj−1/2 = Cj+1/2 (uh )∆j+l+1/2 uh
l=−k

where gj±1/2 are the usual numerical flux functions. Section 4.1.2 provides an example of how
the discrete TVD conditions and discrete conservation can be simultaneously achieved. A more
comprehensive overview of finite volume numerical methods based on TVD constructions can
be found the books by Godlewski and Raviart (1991) and LeVeque (2002).

4.1.2. MUSCL finite volume methods. A general family of TVD FVMs with five-cell stencil is
the monotone upstream-centered scheme for conservation laws (MUSCL) discretization of van
Leer (1979) and van Leer (1985). MUSCL methods utilize a κ-parameter family of interpolation
formulas with slope limiter function Ψ(R) : R 7→ R
 
− 1+κ 1−κ 1
uj+1/2 = uj + Ψ(Ri )∆j−1/2 uh + Ψ ∆j+1/2 uh
4 4 Rj
 
1+κ 1 1−κ
u+
j−1/2 = uj − Ψ ∆j+1/2 uh − Ψ(Rj )∆j−1/2 uh (74)
4 Rj 4

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 29

where Rj is a ratio of successive solution increments

∆j+1/2 uh
Rj ≡ (75)
∆j−1/2 uh
These formulas are self-consistent with a piecewise linear representation of a solution
reconstructed from cell average values uj . When this reconstructed solution is cell averaged,
the linear corrections vanish identically (regardless of the value of the slope limiter Ψ(·)) and
the cell average values uj are obtained. The technique of incorporating limiter functions to
obtain nonoscillatory resolution of discontinuities and steep gradients can be found in the
flux corrected transport method developed earlier by Boris and Book (1973). For convenience,
the interpolation formulas (74) have been written for a uniformly spaced mesh, although the
extension to irregular mesh spacing is straightforward. The unlimited form of this interpolation
is obtained by setting Ψ(R) = 1. In this unlimited case, the truncation error for the
conservation law divergence in (6a) is given by
(κ − 1/3) ∂3
Truncation Error = − (∆x)2 3 f (u)
4 ∂x
Using the MUSCL interpolation formulas given in (74), sufficient conditions to be imposed on
the limiter function Ψ(·) to obtain the discrete TVD property are readily obtained.
Theorem 18. (MUSCL TVD FVM) The fully discrete 1-D scheme
∆t n
un+1 = unj − (g n
− gj−1/2 ), j = 1, . . . , N
j
∆xj j+1/2
with monotone Lipschitz continuous numerical flux function
gj+1/2 = g(u− +
j+1/2 , uj+1/2 )

utilizing the κ-parameter family of MUSCL interpolation formulas (74) and (75) is TVNI if
there exists a Ψ(R) such that ∀R ∈ R
3−κ 1+κ
0 ≤ Ψ(R) ≤ − (1 + α) (76a)
1−κ 1−κ
and
Ψ(R)
0≤ ≤2+α (76b)
R
with α ∈ [−2, 2 (1 − κ)/(1 + κ)] under the time step restriction
max
∆t 2 − (2 + α)κ ∂g
1− ≥0
∆xj 1−κ ∂u j

where
∂g
max  
≡ sup ∂g
u, u+ − ∂g
(u− e
e)
∂u e (e
∂u j+1/2 ) e
e
∂u j−1/2 , u
j ũ∈[u−j−1/2
,u−
j+1/2
]
+
˜
ũ∈[uj−1/2
,u+
j+1/2
]

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
30 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

Table 1. Members of the MUSCL TVD family of schemes.

κ Unlimited scheme βmax Truncation error

1/3 Third-order 4 0
1 ∂3
−1 Fully upwind 2 (∆x)2 3 f (u)
3 ∂x
3
1 2 ∂
0 Fromm’s 3 (∆x) f (u)
12 ∂x3
3
1 ∂
1/2 Low truncation error 5 − (∆x)2 3 f (u)
24 ∂x

From accuracy considerations away from extrema, it is desirable that the unlimited form
of the discretization is obtained. Consequently, the constraint Ψ(1) = 1 is also imposed
upon the limiter function. This constraint, together with the algebraic conditions (76a) and
(76b), is readily achieved using the well-known MinMod limiter with compression parameter
β determined from the TVD analysis
 
MM (3 − κ)
Ψ (R) = max(0, min(R, β)), β ∈ 1, (77)
(1 − κ)

Table 1 summarizes the MUSCL TVD scheme and maximum compression parameter β for a
number of values of κ. Another limiter due to van Leer that meets the technical conditions of
Theorem 18 and also satisfies Ψ(1) = 1 is given by
R + |R|
ΨVL (R) = (78)
1 + |R|
This limiter exhibits differentiability away from R = 0, which improves the iterative
convergence to steady state for many algorithms. Numerous other limiter functions are
considered and analyzed in Sweby (1984).

Unfortunately, TVD schemes locally degenerate to piecewise constant approximations at

smooth extrema, which locally degrades the accuracy. This is an unavoidable consequence of
the strict TVD condition.
Theorem 19. (TVD critical point accuracy, Osher, 1984) The TVD discretizations (68), (70)
and (72) all reduce to at most first-order accuracy at nonsonic critical points, that is, points
u∗ at which f 0 (u∗ ) 6= 0 and u∗x = 0.

The next class of finite volume methods addresses this deterioration at smooth extrema.

4.1.3. ENO/WENO finite volume methods. To address the degradation in accuracy of

TVD methods at critical points, Harten et al., 1986 proposed a new class of finite volume

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 31

discretizations in 1-D based on a weaker form of total variation control, see also Harten et al.,
1987; Harten, 1989. These discretizations utilize polynomial reconstruction from cell averages.
Let Pp (K) denote polynomials of degree at most p in the control volume K and Vhp the broken
space of piecewise p-order polynomials for a tesselation T , that is,
Vhp = {v v|K ∈ Pp (K), ∀K ∈ T } (79)
Let Rp0 : Vh0 7→Vhpdenote a polynomial reconstruction operator that maps cell averages to the
broken space of p-order polynomials. The evolution of solution total variation in the FVM may
be understood by first constructing the following operator composition form of the FVM for
a time slab increment [tn , tn+1 ]
uhn+1 = A · E(τ ) · Rp0 unh (80)
In this equation, E(τ ) is the evolution operator for the PDE, and A is the cell averaging
operator. Since A is a nonnegative operator and E(τ ) represents exact evolution over short
time, the evolution of the solution total variation satisfies
TV(un+1
h ) = TV(A · E(τ ) · Rp0 unh ) ≤ TV(Rp0 unh )
This indicates that the control of solution total variation depends entirely on total variation
properties of the reconstruction operator. The requirements of non-oscillatory high-order
accuracy for smooth solutions and discrete conservation suggests the following additional
design objectives for the reconstruction operator
• Rp0 (x; uh ) = u(x) + e(x)∆xp+1 + O(∆xp+2 ) (81a)
to insure accuracy whenever the infinite-dimensional solution u is smooth,
• A|K Rp0 uh = uh |K = uK (81b)
to insure discrete conservation,
• TV(Rp0 unh ) ≤ TV(unh ) + O(∆xp+1 ) (81c)
to insure an essentially nonoscillatory (ENO) reconstruction.
The last design objective (81c) permits a small O(∆xp+1 ) increase in total variation, thus
yielding
TV(un+1
h ) ≤ TV(Rp0 unh ) ≤ TV(unh ) + O(∆xp+1 )
which motivates the ENO name given to the reconstruction. Harten points out that the error
function e(x) may not be Lipschitz continuous at certain points so that the cumulative error in
the scheme is O(∆xp ) in a maximum norm but remains O(∆xp+1 ) in an L1 -norm. To achieve
the ENO reconstruction design criterion (81), Harten and coworkers considered breaking the
task into two parts:

1. Polynomial reconstruction from a given stencil of cell averages

2. Construction of the “smoothest” polynomial approximation by an adaptive stencil
selection algorithm

The resulting ENO reconstruction is then used to calculate numerical flux function states
u± ±
j+1/2 and uj−1/2 in the finite volume discretization

d 1
uj = − (g(u− + − +
j+1/2 , uj+1/2 ) − g(uj−1/2 , uj−1/2 )), j = 1, . . . , N
dt ∆xj

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
32 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

The semi-discrete equations are then evolved forward in time using some form of high-order
accurate time stepping that also prevents the introduction of spurious oscillations (Section 4.4).

In the following section, a commonly used reconstruction technique from cell averages is
presented. This is followed by a description of the adaptive stencil algorithm for obtaining the
smoothest polynomial approximation.

4.1.4. ENO reconstruction via primitive function. Given cell averages uj of a piecewise
smooth function u(x), one can readily evaluate pointwise values using the notion of a primitive
function U (x) given by Z x
U (x) = u(ξ) dξ
x0

by exploiting the relationship

j
X
∆xi ui = U (xj+1/2 ) = Uj+1/2
i=1

Let Hp+1 (x; U ) denote a (p+1)-order polynomial interpolant of a function U defined piecewise
for intervals [xj−1/2 , xj+1/2 ], j = 1, . . . , N . Since
d
u(x) ≡ U (x)
dx
a p-order reconstruction operator is obtained from
d
Rp0 (x; uh ) = Hp+1 (x; U (uh ))
dx
and consequently for smooth data
dl dl
R p (x; uh ) = u(x) + O(∆xp+1−l )
dxl dxl
By virtue of the use of the primitive function U (x), it follows that
A|K Rp0 (x; uh ) = uK
and from the polynomial interpolation problem for smooth data
Rp0 (x; uh ) = u(x) + O(∆xp+1 )
as desired.

4.1.5. ENO smoothest polynomial approximation. The piecewise polynomial interpolants

Hp+1 (x; U ) described in Section 4.1.4 depend on the particular choice of stencil for pointwise
values of the primitive function U for each cell [xj−1/2 , xj+1/2 ]. Consequently, the resulting
reconstructions will not generally satisfy the oscillation requirement (81c). This motivated
Harten and coworkers to consider a new algorithm for smoothest polynomial interpolation

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 33

using adaptive stencils. When used to interpolate pointwise values of the primitive function,
the resulting reconstruction satisfies (81a–c). Specifically, a high-order accurate interpolant is
constructed
dk 0 dk
k
Rp (x; uh ) = u(x) + O(∆xp+1−k ) , 0 ≤ k ≤ p
dx dxk
which avoids having Gibbs oscillations at discontinuities in the sense
T V (Rp0 (x; uh )) ≤ T V (u) + O(∆xp+1 )
The strategy pursued by Harten and coworkers was to construct an ENO polynomial
ENO
Hp+1 (x; U (uh )) in each interval [xj−1/2 , xj+1/2 ], j = 1, . . . , N , which interpolates the primitive
function U (x) at the p + 2 successive points {xi−1/2 }, ip (j) ≤ i ≤ ip (j) + p + 1 that include
xj−1/2 and xj+1/2 . This describes p + 1 possible polynomials depending on the choice of ip (j)
for an interval [xj−1/2 , xj+1/2 ]. The ENO strategy selects the value ip (j) for each interval that
produces the “smoothest” polynomial interpolant for a given input data. Information about
smoothness of U (x) is extracted from a table of divided differences defined recursively
U [xi−1/2 ] = U (xi−1/2 )
U [xi+1/2 ] − U [xi−1/2 ]
U [xi−1/2 , xi+1/2 ] =
xi+1/2 − xi−1/2
..
.
U [xi+1/2 , . . . , xi+k−1/2 ] − U [xi−1/2 , . . . , xi+k−3/2 ]
U [xi−1/2 , . . . , xi+k−1/2 ] =
xi+k−1/2 − xi−1/2
The stencil producing the smoothest interpolant is then chosen hierarchically by setting
i1 (j) = j
and for k = 1, . . . , p

ik (j) − 1 if |U [xik (j)−3/2 , . . . , xik (j)+k−1/2 ]| < |U [xik (j)−1/2 , . . . , xik (j)+k+1/2 ]|
ik+1 (j) =
ik (j) otherwise
(82)
Harten et al. (1986) demonstrated that this interpolant is monotone in any cell interval
containing a discontinuity and the resulting reconstruction satisfies the design objectives (81a-
c).

4.1.6. WENO reconstruction. The solution adaptive nature of the ENO stencil selection
algorithm (82) yields nondifferentiable fluxes that may impede the performance of solution
algorithms Jiang and Shu (1996). Recall that the p-order ENO reconstruction considers
p + 1 possible polynomial stencils. The stencil selection algorithm chooses only one of these
possible stencils and other slightly less smooth stencils may give similar O(∆xp+1 ) accuracy.
(0) (2) (p)
However, using a linear combination of all p+1 possible polynomials {Pp , Pp , . . . , Pp } with
2p+1
optimized weights ωk , k = 0, . . . , p potentially yields a more accurate O(∆x ) interpolant
for smooth enough data
p
X p
X
Pp (x) = ωk Pp(k) (x) + O(∆x2p+1 ) ωk = 1
k=0 k=0

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
34 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

For example, optimized weights for p = 0, 1, 2 are readily computed

p = 0: ω0 = 1
2 1
p = 1: ω0 = , ω1 =
3 3
3 3 1
p = 2: ω0 = , ω1 = , ω 2 =
10 5 10
In the weighted essentially nonoscillatory (WENO) schemes of Jiang and Shu (1996) and Shu
(1999), approximate weights ω ek , k = 0, . . . , p are devised such that for smooth solutions
ek = ωk + O(∆xp )
ω
so that the O(∆x2p+1 ) accuracy is still retained using these approximations
p
X p
X
Pp (x) = ek Pp(k) (x) + O(∆x2p+1 ) ,
ω ek = 1
ω
k=0 k=0

The approximate weights are constructed using the ad hoc formulas for k = 0, . . . , p
ωk αk
αk = ek = p
, ω
( + βk )2 X
αl
l=0

where  is an approximation to the square root of the machine precision and βk is a smoothness
indicator !2
Xp Z xj+1/2 (k)
2l−1 dl Pp (x)
βk = ∆x dx
xj−1/2 ∂xl
l=1

For a sequence of smooth solutions with decreasing smoothness indicator βk , these formulas
approach the optimized weights, ω ek → ωk . These formulas also yield vanishing weights
ek → 0 for stencils with large values of the smoothness indicator such as those encountered
ω
at discontinuities. In this way, the WENO construction retains some of the attributes of the
original ENO formulation but with increased accuracy in smooth solution regions and improved
differentiability often yielding superior robustness for steady-state calculations.

4.2. Higher order accurate finite volume methods for hyperbolic problems in multiple
dimensions

Although the one-dimensional finite volume discretizations given in previous sections may be
readily applied in multidimensions on a dimension-by-dimension basis, a result of Goodman
and LeVeque (1985) shows that TVD schemes in two or more space dimensions are only first-
order accurate.
Theorem 20. (Accuracy of TVD finite volume discretization on multidimensional Cartesian
meshes) Any two-dimensional finite volume scheme of the form
∆t ∆t
un+1 n
i,j = ui,j −
n
(gi+1/2,j n
− gi−1/2,j )− (hn − hni,j−1/2 ), 1 ≤ i ≤ M, 1 ≤ j ≤ N
|Ki,j | |Ki,j | i,j+1/2

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 35

with Lipschitz continuous numerical fluxes for integers p, q, r, s

gi+1/2,j = g(ui−p,j−q , . . . , ui+r,j+s )
hi,j+1/2 = h(ui−p,j−q , . . . , ui+r,j+s )

that is TVNI in the sense

T V (un+1
h ) ≤ T V (unh )
where
Xh i
T V (u) ≡ ∆yi+1/2,j |ui+1,j − ui,j | + ∆xi,j+1/2 |ui,j+1 − ui,j |
i,j

is at most first order accurate.

Motivated by the negative results of Goodman and LeVeque, weaker conditions yielding
solution monotonicity preservation have been developed from discrete maximum principle
analysis. These alternative constructions have the positive attribute that they extend to
unstructured meshes as well.

4.2.1. Positive coefficient FVMs on structured meshes. Lemma 3 considers FVMs of the form
∆t X
un+1
K = unK + CK,L (unh )(unL − unK )
|K|
σK,L ⊂K

with uh = {uK1 , uK2 , . . .} and establishes a local discrete maximum principle

min (unK , unL ) ≤ un+1
K ≤ max (unK , unL )
σK,L ⊂K σK,L ⊂K

for each K ∈ T and n = 0, 1, 2, . . . under a CFL-like condition on the time step parameter if
all coefficients CK,L (uh ) are nonnegative. Discretizations of this type are often called positive
coefficient discretizations or more simply positive discretizations. To circumvent the negative
result of Theorem 20, Spekreijse (1987) developed a family of high-order accurate positive
coefficient discretizations on two-dimensional structured meshes. For purposes of positivity
analysis, these methods are written in incremental form on a M × N logically rectangular 2-D
mesh

un+1
i,j = uni,j + ∆t Ani+1,j (uni+1,j − uni,j ) + Bi,j+1
n
(uni,j+1 − uni,j )

n
+ Ci−1,j (uni−1,j − uni,j ) + Di,j−1
n
(uni,j−1 − uni,j ) , 1 ≤ i ≤ M, 1 ≤ j ≤ N (83)
where the coefficients are nonlinear functions of the solution
Ani+1,j = A(. . . , uni−1,j , uni,j , uni+1,j , . . .)
n
Bi,j+1 = B(. . . , uni,j−1 , uni,j , uni,j+1 , . . .)
n
Ci−1,j = C(. . . , uni−1,j , uni,j , uni+1,j , . . .)
n
Di,j−1 = D(. . . , uni,j−1 , uni,j , uni,j+1 , . . .)
Once written in incremental form, the following lemma follows from standard positive
coefficient maximum principle analysis.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
36 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

Lemma 4. (Positive coefficient FVMs in multidimensions) The discretization (83) is a

positive coefficient FVM if for each 1 ≤ i ≤ M, 1 ≤ j ≤ N and time slab increment
[tn , tn+1 ], n = 0, 1, 2, . . .

Ani+1,j ≥ 0, Bi,j+1
n n
≥ 0, Ci−1,j n
≥ 0, Di,j−1 ≥0 (84)

and
1 − ∆t (Ani+1,j + Bi,j+1
n n
+ Ci−1,j n
+ Di,j−1 )≥0 (85)
with discrete maximum principle

min(uni,j , uni−1,j , uni+1,j , uni,j−1 , uni,j+1 ) ≤ un+1 n n n n n

i,j ≤ max(ui,j , ui−1,j , ui+1,j , ui,j−1 , ui,j+1 )

Using a procedure similar to that used in the development of MUSCL TVD FVMs in 1-D,
Spekreijse (1987) developed a family of monotonicity preserving MUSCL interpolations from
the positivity conditions of Lemma 4.
Theorem 21. (MUSCL positive coefficient FVM) Assume a fully discrete 2-D FVM
∆t ∆t
un+1 n
i,j = ui,j − (g n n
− gi−1/2,j )− (hn − hni,j−1/2 ), 1 ≤ i ≤ M, 1 ≤ j ≤ N
|Ki,j | i+1/2,j |Ki,j | i,j+1/2
for n = 0, 1, 2, . . . utilizing monotone Lipschitz continuous numerical flux functions

gi+1/2,j = g(u− +
i+1/2,j , ui+1/2,j )

hi,j+1/2 = h(u− +
i,j+1/2 , ui,j+1/2 )

and MUSCL extrapolation formulas

1
u−
i+1/2,j = ui,j + Ψ(Ri,j )(ui,j − ui−1,j )
2  
1 1
u+
i−1/2,j = ui,j − Ψ (ui+1,j − ui,j )
2 Ri,j
1
u−
i,j+1/2 = ui,j + Ψ(Si,j )(ui,j − ui,j−1 )
2  
1 1
u+
i,j−1/2 = ui,j − Ψ (ui,j+1 − ui,j )
2 Si,j
with
ui+1,j − ui,j ui,j+1 − ui,j
Ri,j ≡ , Si,j ≡
ui,j − ui−1,j ui,j − ui,j−1
This finite volume discretization satisfies the local maximum principle properties of Lemma 4
and is second-order accurate if the limiter Ψ = Ψ(R) has the properties that there exist
constants β ∈ (0, ∞), α ∈ [−2, 0] such that ∀R ∈ R
Ψ(R)
α ≤ Ψ(R) ≤ β, −β ≤ ≤2+α (86)
R
Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 37

with the constraint Ψ(1) = 1 and the smoothness condition Ψ(R) ∈ C 2 near R = 1 together
with a time step restriction for stability
!
n,max n,max
∆t ∂g ∂h
1 − (1 + β) + ≥0
|Ki,j | ∂u i,j ∂u i,j

where
∂g
max  
≡ sup ∂g
u, u+ − ∂g
(u− e
e) ≥0
∂u e (e
∂u i+1/2,j ) e
e
∂u i−1/2,j , u
i,j ũ∈[u−
i−1/2,j
,u−
i+1/2,j
]
+
˜
ũ∈[ui−1/2,j
,u+
i+1/2,j
]

∂h
max  
≡ sup ∂h
u, u+ − ∂h
(u− b
b) ≥0
∂u b (b
∂u i,j+1/2 ) b
b
∂u i,j−1/2 , u
i,j û∈[u−
i,j−1/2
,u−
i,j+1/2
]
ˆ
û∈[u+
,u+ ]
i,j−1/2 i,j+1/2

Many limiter functions satisfy the technical conditions (86) of Theorem 21. Some examples
include

• the van Leer limiter (78)

R + |R|
ΨVL (R) =
1 + |R|
• Koren limiter Koren (1988)
R + 2R2
ΨK (R) =
2 − R + 2R2

For smooth solutions, the Koren limiter results in state extrapolations identical to the most
accurate κ = 1/3 MUSCL interpolations of Section 4.1.2.

4.3. Higher order accurate finite volume methods for hyperbolic problems on unstructured
meshes

FVMs for hyperbolic problems have been extended to unstructured meshes using data
reconstruction from cell averages. The polynomial reconstruction operator in the FVM maps
cell-averaged data into the broken space (79) consisting of p-order piecewise R polynomials in
each control volume, Rp0 : Vh0 7→ Vhp , while preserving the cell average value, K Rp0 (x; uh ) dx =
uK |K| for all K ∈ T . Using this reconstruction operator, higher order accurate finite volume
discretizations on unstructured meshes using polynomial data reconstruction are of the general
form for time-invariant control volumes K
duK 1 X
=− ωq g(u− +
σ (xq,σ ; uh ), uσ (xq,σ ; uh ); nσ ) |σ|, ∀K ∈ T (87)
dt |K|
σ⊂∂K
1≤q≤Q

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
38 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

using Q-point quadrature with positive weights ωq ∈ R+ and locations xq,σ ∈ σ for
q = 1, . . . , Q. In this formula, each σ interface is assumed planar with unique normal nσ
and g(u, v; n) denotes any of the numerical fluxes described in Section 3. The σ interface
states used in the numerical flux quadrature are calculated from the reconstruction Rp0 (x; uh )

u± 0
σ (x; uh ) ≡ lim Rp (x ± nσ ; uh ) (88)
↓0

Due to the cell-wise piecewise structure of the reconstruction operator, the two states u±
σ (x; uh )
are generally distinct.

4.3.1. General p-exact reconstruction operators on unstructured meshes. The reconstruction

operator discussed in Section 4.1.4 exploits properties of a primitive function in 1-D. This
approach is problematic to extend to general unstructured meshes. Consequently, more general
reconstruction design principles have been developed and used in implementations. Abstractly,
the reconstruction operator Rp0 appearing in (88) serves as a finite-dimensional pseudoinverse
of the cell averaging operator A. The development of a general polynomial reconstruction
operator, Rp0 , that reconstructs p-degree polynomials from cell-averaged data on unstructured
meshes follows from the application of a small number of design principles:

1. (Conservation of the mean) Given solution cell averages uh ∈ Vh0 , the reconstruction
Rp0 uh is required to have the correct cell average, that is,

if v = Rp0 uh then uh = Av
More concisely,
ARp0 = I
so that Rp0 is a right inverse of the averaging operator A.
2. (p-exactness) A reconstruction operator Rp0 is p-exact if Rp0 A reconstructs polynomials
of degree p or less exactly, that is,
if u ∈ Pp and v = Au then Rp0 v = u
This can be written succinctly as
Rp0 A|Pp = I
so that Rp0 is a left inverse of the averaging operator A restricted to the space of
polynomials of degree at most p.
3. (Compact support) The reconstruction in a control volume K should only depend on
cell averages in a relatively small
neighborhood surrounding K. Recall that a polynomial
of degree p in Rd contains p+d d degrees of freedom. The support set for K is required
to contain at least this number of neighbors. As the support set becomes even larger
for fixed p, not only does the computational cost increase, but eventually, the accuracy
decreases as less valid data from further away is brought into the calculation.

Practical implementations of general p-order polynomial reconstruction operators fall into two
classes:

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 39

• Fixed support stencil reconstructions. These methods choose a fixed support set as a
preprocessing step. Various limiting strategies are then employed to obtain nonoscillatory
approximation; see, for example, Barth and Frederickson (1990) and Delanaye (1996) for
further details.
• Adaptive support stencil reconstructions. These ENO-like methods dynamically choose
reconstruction stencils based on solution smoothness criteria; see, for example, Harten
and Chakravarthy (1991), Vankeirsblick (1993), Abgrall (1994), Sonar (1997), and Sonar
(1998) for further details.

4.3.2. Higher order finite volume methods for hyperbolic problems on unstructured meshes
using linear reconstruction. Considerable simplification is possible when only linear
reconstruction polynomials are sought, R10 (x; uh ). The maximum principle analysis presented
next not only gives sufficient conditions for a discrete maximum principle but also explicitly
exposes the dependence of cell shape with respect to maximum allowable time step for the
fully discrete approximation. This geometrical shape parameter for a control volume K can
be calculated from the formula
−1
Γgeom
K = sup αK (θ) (89)
0≤θ≤2π

where 0 < αK (θ) < 1 represents the smallest fractional perpendicular distance from the
centroid to one of two minimally separated parallel hyperplanes with orientation θ and
hyperplane location such that all quadrature points in the control volume lie between or on the
hyperplanes as shown in Figure 3. Table 2 lists Γgeom values for various control volume shapes

Figure 3. Minimally separated hyperplanes hL (θ) and hU (θ) and the fractional distance ratio α(θ) for
use in the calculation of Γgeom .

in R1 , R2 , R3 , and Rd . As might be expected, those geometries that have exact quadrature

point symmetry with respect to the control volume centroid have geometric shape parameters
Γgeom equal to 2 regardless of the number of space dimensions involved.

Using standard maximum principle analysis, a discrete maximum principle is obtained under
a CFL-like time step restriction if the solution reconstruction in each and every control volume

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
40 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

Table 2. Reconstruction geometry factors for various control volume shapes utilizing midpoint flux
quadrature rule.

Control volume shape Space dimension Γgeom

Segment 1 2
Triangle 2 3
Parallelogram 2 2
Tetrahedron 3 4
parallelepiped 3 2
Simplex d d+1
Hyper-parallelepiped d 2

K can be bounded from above and below respectively by the neighboring cell averages. This
restriction is stated more precisely in the following theorem:
Theorem 22. (Finite volume maximum principle on unstructured meshes, R10 ) Let umin K and
umax
K denote the minimum and maximum value of solution cell averages for a given cell K and
all adjacent cell neighbors, that is,

umin
K ≡ min (uK , uL ) and umax
K ≡ max (uK , uL ) (90)
σK,L ⊂∂K σK,L ⊂∂K

The fully discrete finite volume scheme

∆t X
un+1
K = unK − ωq g(u− +
σ (xq,σ ; uh ), uσ (xq,σ ; uh ); nσ ) |σ|, ∀K ∈ T (91)
|K|
σ⊂∂K
1≤q≤Q

with monotone Lipschitz continuous numerical flux function, Q-point quadrature with
nonnegative weights, and linear reconstructions

u± 0
σ (x; uh ) ≡ lim R1 (x ± nσ ; uh ) (92)
↓0

exhibits the local maximum principle for each K ∈ T and n = 0, 1, 2, . . .

umin,n
K ≤ un+1
K ≤ umax,n
K

under the time step restriction

∆t geom X ∂g e
1− Γ sup (e e; nσ )|σ| ≥ 0
u, u (93)
|K| K min,n
,umax,n
e
e
∂u
σ ⊂ ∂K ũ∈[uK K ]
1≤q≤Q ũ∈[u
˜ min,n
,umax,n ]
K K

if the linear reconstruction evaluated at interface quadrature points satisfies

max(umin,n
K , umin,n
L ) ≤ u−,n max,n
σ (xq,σ ) ≤ min(uK , umax,n
L ), ∀σK,L ⊂ ∂K (94)

for all xq,σ ∈ σ, q = 1, . . . , Q.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 41

Equation (94) can be interpreted as stating that the data reconstruction in a cell, when
evaluated at quadrature points, should be bounded from above and below by cell-averaged
values of adjacent neighbors (including itself). This is completely consistent with the
interpretation originally given in van Leer (1979).

Comparing the maximum allowable time step using piecewise constant approximation (52)
with piecewise linear approximation (93) indicates that the maximum time step using piecewise
linear approximation is reduced by a factor of approximately 1/Γgeom K . Using shape regular
quadrilaterals (2-D) or hexahedron (3-D), this time step reduction is approximately 1/2. Using
triangles (2-D) or tetrahedra (3-D), this time step reduction is approximately 1/3 and 1/4,
respectively.

4.3.3. Slope limiters for linear reconstruction. Given a linear reconstruction R10 (x; uh ) that
does not necessarily satisfy the requirements of Theorem 22, it is possible to modify the
reconstruction so that the modified reconstruction does satisfy the requirements of Theorem 22.
For each control volume K ∈ T , assume a modified local reconstruction operator R e0 (x; uh )|K
1
of the form
Re10 (x; uh )|K = uK + αK (R10 (x; uh )|K − uK )
for αK ∈ [0, 1]. By construction, this modified reconstruction correctly reproduces the control
volume cell average for all values of αK , that is,
Z
1 e0 (x; uh ) dx = uK
R (95)
|K| K 1
The most restrictive value of αK for each control volume K is then computed on the basis of
the Theorem 22 constraint (94)
 min(umax ,umax )−uK


K L
R0 (xq,σ ;uh )|K −uK
if R10 (xq,σ ; uh )|K > min(umax max
K , uL )

 1
MM max(umin min
αK = min K ,uL )−uK (96)
σK,L ⊂∂K  R10 (xq,σ ;uh )|K −uK if R10 (xq,σ ; uh )|K < max(umin min
K , uL )
1≤q≤Q  
1 otherwise
where umax and umin are defined in (90). When the resulting modified reconstruction operator
is used in the extrapolation formulas (92), the discrete maximum principle of Theorem 22 is
attained under a CFL-like time step restriction. Utilizing the inequalities
max(uK , uL ) ≤ min(umax max min min
K , uL ) and min(uK , uL ) ≥ max(uK , uL )

it is straightforward to construct a more restrictive limiter function

 max(uK ,uL )−uK 0
 R10 (xq,σ ;uh )|K −uK if R1 (xq,σ ; uh )|K > max(uK , uL )



LM min(uK ,uL )−uK
αK = min 0 (97)
 R10 (xq,σ ;uh )|K −uK if R1 (xq,σ ; uh )|K < min(uK , uL )
σK,L ⊂∂K 
1≤q≤Q  
1 otherwise
that yields modified reconstructions satisfying the technical conditions of Theorem 22. This
simplified limiter (97) introduces additional slope reduction when compared to (96). This can

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
42 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

be detrimental to the overall accuracy of the discretization. The limiter strategy (97) and other
variants for simplicial control volumes are discussed further in Liu (1993), Wierse (1994), and
Batten et al. (1996). Note that in practical implementations, both limiters (96) and (97) require
some modification to prevent division by zero for constant solution data.

4.3.4. Linear reconstruction on simplicial control volumes. Linear reconstruction operators

on general control volumes that satisfy the cell averaging requirement often exploit the fact
that the cell average is also a pointwise value of any valid linear reconstruction evaluated at
the centroid of the control volume. This reduces the reconstruction problem to that of gradient
estimation given pointwise samples at the centroids. In this case, it is convenient to express
the reconstruction in the form
R10 (x; uh )|K = uK + (∇u)K · (x − xK ) (98)
where xK denotes the centroid for the control volume K and (∇u)K is the gradient
to be determined. Figure 4 depicts a 2-D simplex KO and three adjacent neighboring
simplices. Also shown are the corresponding four pointwise solution values {uA , uB , uC , uO }
located at centroids of each simplex. By selecting any three of the four pointwise solution

2
1
O

A
B

Figure 4. Triangle control volume KO (shaded) with three adjacent cell neighbors.

values, a set of four possible gradients are uniquely determined, that is, {∇(uA , uB , uC ),
∇(uA , uB , uO ),∇(uB , uC , uO ),∇(uC , uA , uO )}. A number of slope limited reconstruction
techniques are possible for use in the finite volume scheme (91) that meet the technical
conditions of Theorem 22.

1. Choose (∇u)KO = ∇(uA , uB , uC ) and limit the resulting reconstruction using (96) or
(97).
2. Limit four reconstructions corresponding to (∇u)KO equal to ∇(uA , uB , uC ),
∇(uA , uB , uO ), ∇(uB , uC , uO ), and ∇(uC , uA , uO ) using (96) or (97) and choose the
limited reconstruction with largest gradient magnitude. This technique is a generalization
of that described in Batten et al. (1996) wherein limiter (97) is used.
3. Evaluate four reconstructions corresponding to (∇u)KO equal to ∇(uA , uB , uC ),
∇(uA , uB , uO ), ∇(uB , uC , uO ), and ∇(uC , uA , uO ) and choose the largest gradient
magnitude that satisfies the maximum principle reconstruction bound inequality (94). If
all reconstructions fail the bound inequality, the reconstruction gradient is set equal to
zero; see Liu (1993).

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 43

4.3.5. Linear reconstruction on general control volume shapes. It is again convenient to

express the linear reconstruction in the form

R10 (x; uh )|K = uK + (∇u)K · (x − xK ) (99)

but now the shape and number of adjacent neighboring cells is irregular. Two common
techniques for simplified linear reconstruction include a simplified least squares technique and
a Green-Gauss quadrature technique.

Least squares linear reconstruction Again exploiting the fact that the cell average value is also
a pointwise value of the linear reconstruction evaluated at the centroid of a general control
volume shape, the task of linear reconstruction reduces to the problem of gradient estimation
given pointwise values. In the simplified least squares reconstruction technique, a triangulation
(2-D) or tetrahedralization (3-D) of centroids is constructed as shown in Figure 5. Referring

Figure 5. Triangulation of centroid locations showing a typical control volume K0 associated with
cyclically indexed graph neighbors Kj , j = 1, . . . , N0 .

to this figure, for each of the N0 edges of the simplex mesh incident to K0 , an edge projected
gradient constraint equation is constructed subject to a specified nonzero scaling wK0 for
j = 1, . . . , N0
wj (∇u)K0 · (xj − x0 ) = wk (uKj − uK0 )

The number of edges incident to a simplicial mesh vertex in Rd is greater than d away from
mesh boundaries thereby producing the following nonsquare matrix of constraint equations
   
w1 ∆x1 w1 ∆y1 w1 (uK1 − uK0 )
 .. ..   .. 
 . .  (∇u)K0 =  . 
wN0 ∆xN0 wN0 ∆yN0 wN0 (uKN0 − uK0 )

or in abstract form h i
~1
L ~ 2 ∇u = f~
L

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
44 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

This abstract form can be solved in a least squares sense using a variety of techniques, for
example, Gram-Schmidt, modified Gram-Schmidt, Householder rotations, SVD, and so on.
When least squares matrix conditioning is not an issue, the following symbolic solution may
be used !
1 ~ 1 · f~) − l12 (L
l22 (L ~ 2 · f~)
∇u = 2 ~ 2 · f~) − l12 (L
~ 1 · f~) (100)
l11 l22 − l12 l11 (L
~i ·L
with lij = L ~ j . The form of this solution in terms of scalar dot products over incident edges
suggests that the least squares linear reconstruction can be efficiently computed via an edge
data structure without the need for storing a nonsquare matrix.

Green–Gauss linear reconstruction Again referring to Figure 5, reconstruction gradients may

be approximated from mean value approximation and application of the Green–Gauss identity
Z Z N0 Z
X N0
X 1
|K0 |(∇u)K0 ≈ ∇u dx = u n dx = u n dx ≈ (uK0 + uKj ) ν 0j (101)
K0 ∂K0 j=1 σK0 ,Kj j=1
2
R
In this formula, ν 0j ≡ σK0 ,Kj
n dx which is independent of the shape of σK0 ,Kj and only
depends on its boundary shape (end points in 2-d). A notable property of this formula is
that the gradient calculation is exact whenever the numerical solution varies linearly over the
support of the reconstruction.

The slope limiting procedures (96) or (97) may then be applied to the least squares or Gauss-
Gauss reconstruction so that when used in the finite volume discretization (91) the discrete
maximum principle of Theorem 22 is obtained. Gradient reconstruction together with a slope
limitation on general meshes can be found in Buffard and Clain, 2010; Calgaro et al., 2010.
In these works, the limitation procedure is presented as a limitation of the slope defined by
the cell and face values on the basis of its comparison with other slopes defined by the values
taken by the solution in the neighborhood. Then, under geometric assumptions for the mesh,
this limitation may be shown to imply some conditions for the approximation at the face,
which again ensure, at least for pure convection problems, a local maximum principle Clain
and Clauzon, 2010; Clain, 2013; Berthon et al., 2014.

Another novel approach found in Tran, 2008; Michel et al., 2010 and Piar et al., 2013;
Therme, 2015 is based on the observation that for a linear convection term, stability conditions
may be exploited in the calculation of an admissible interval for the state value at the face.
This suggests a crude limitation process, which does not use any slope computation and simply
consists of performing a (one-dimensional) projection of the tentative affine reconstructed face
value on this interval. In addition, stability conditions are purely algebraic (in the sense that
they do not require any geometric computation) and thus work with arbitrary meshes.

4.3.6. Positive coefficient finite volume methods on unstructured meshes. Several related
positive coefficient schemes have been proposed for multidimensional simplicial meshes
using one-dimensional interpolation. The simplest example is the upwind triangle scheme as
introduced by Billey et al. (1987), Desideri and Dervieux (1988) and Rostand and Stoufflet

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 45

(1988) with later improved variants given by Jameson (1993) and Cournède, Debiez and
Dervieux (1998). These schemes are not finite volume methods in the sense described in Section
4.3.2 owing to the fact that a single multidimensional gradient is not obtained in each control
volume. Referring to Figure 6, the starting point for these methods is the semidiscrete FVM

b′

k′

k
b e′ a′
j e

j′

Figure 6. Triangle complex used in the upwind triangle schemes showing the linear extension of ejk
into neighboring triangle for the determination of points xj 0 and xk0 .

(29) from Section 3.1

d 1 X
uj = − g(uj , uk ; nj,k ) |σj,k | (102)
dt |Kj |
σj,k ⊂∂Kj

for each Kj ∈ T with numerical flux function

1 1
g(u, v; n) = (f (u) + f (v)) · n − |a(u, v; n)|(v − u) (103)
2 2

utilizing the mean value speed satisfying

(f (v) − f (u)) · n = a(u, v; n)(v − u) (104)

With further modifications at sonic points, the modified numerical flux can be shown to
be an E-flux. As discussed earlier, this FVM is at most first order accurate for hyperbolic
conservation laws. The main idea in the upwind triangle scheme is to add antidiffusion terms
to the numerical flux function (103) such that the sum of added diffusion and antidiffusion
terms in the numerical flux function vanishes entirely whenever the numerical solution varies
linearly over the support of the flux function. The amount of added antidiffusion is determined
from maximum principle analysis. The resulting method and maximum principle results are
summarized in the following theorem:

Theorem 23. (LED and maximum principles for the upwind triangle scheme) Referring to
simplex configuration in Figure 6, let uj denote the nodal solution value at a simplex vertex vj
in one-to-one correspondence with control volumes Kj ∈ T . Let g(uj 0 , uj , uk , uk0 ; njk ) denote

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
46 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

the numerical flux function with limiter function Ψ(·) : R 7→ R

1
g(uj 0 , uj , uk , uk0 ) ≡ (f (uj ) + f (uk )) · njk
2   
1 + hjk ∆j 0 j u
− a (uj , uk ; njk ) 1 − Ψ (uk − uj )
2 hj 0 j ∆jk u
  
1 hjk ∆kk0 u
+ a− (uj , uk ; njk ) 1 − Ψ (uk − uj )
2 hkk0 ∆jk u
utilizing the mean value speed a(uj , uk ; njk ) satisfying (104) and variable spacing parameter
hjk = |xk − xj |. The semidiscrete FVM
d 1 X
uj = − g(uj 0 , uj , uk , uk0 ; nj,k ) |σj,k |, ∀Kj ∈ T
dt |Kj |
σj,k ⊂∂Kj

with linearly interpolated values uj 0 and uk0 as depicted in Figure 6 is local extremum
diminishing (LED) in the sense of Lemma 2 and exhibits the local spatial maximum principle
at steady state u∗
min u∗k ≤ u∗j ≤ max u∗k
σj,k ⊂∂Kj σj,k ⊂∂Kj

if the limiter Ψ(R) satisfies ∀R ∈ R

[Ψ(R)]
0≤ , 0 ≤ Ψ(R) ≤ 2
R

Some standard limiter functions that satisfy the requirements of Theorem 23 include

• the MinMod limiter with maximum compression parameter equal to 2

ΨMM (R) = max(0, min(R, 2))
• the van Leer limiter
R + |R|
ΨVL (R) =
1 + |R|

Other limiter formulations involving three successive one-dimensional slopes are given in
Jameson (1993) and Cournède, Debiez and Dervieux (1998).

4.4. Higher order accurate time integration schemes

The derivation of finite volume schemes in Section 3 began with a semidiscrete formulation
(29) that was later extended to a fully discrete formulation (31) by the introduction of first-
order accurate forward Euler time integration. These latter schemes were then subsequently
extended to higher order accuracy in space using a variety of techniques. For many computing
problems of interest, first-order accuracy in time is then no longer enough. To overcome this

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 47

low-order accuracy in time, a general class of higher order accurate time integration methods
was developed that preserve stability properties of the fully discrete scheme with forward Euler
time integration. Following Gottlieb, Shu and Tadmor (2001) and Shu (2002), these methods
will be referred to as strong stability preserving (SSP) time integration methods.

Explicit SSP Runge–Kutta methods were originally developed by Shu (1988), Shu and Osher
(1988) and Gottlieb and Shu (1998) and called TVD Runge–Kutta time discretizations. In a
slightly more general approach, total variation bounded (TVB) Runge–Kutta methods were
considered by Cockburn et al. (1989), Cockburn et al. (1989), Cockburn, Hou and Shu (1990)
and Cockburn and Shu (1998b) in combination with the discontinuous Galerkin discretization
in space. Küther (2000) later gave error estimates for second-order TVD Runge–Kutta finite
volume approximations of hyperbolic conservation laws.

To present the general framework of SSP Runge–Kutta methods, consider writing the
semidiscrete FVM in the following form
d
U (t) = L(U (t)) (105)
dt
where U = U (t) denotes the solution vector of the semidiscrete FVM. Using this notation
together with forward Euler time integration yields the fully discrete form

U n+1 = U n − ∆t L(U n ) (106)

where U n is now an approximation of U (tn ). As demonstrated in Section 3.3, the forward

Euler time discretization is stable with respect to the L∞ -norm, that is,

kU n+1 k∞ ≤ kU n k∞ (107)

subject to a CFL-like time step restriction

∆t ≤ ∆t0 (108)

With this assumption, a time integration method is said to be SSP (Gottlieb, Shu and Tadmor,
2001) if it preserves the stability property (107), albeit with perhaps a slightly different
restriction on the time step
∆t ≤ c ∆t0 (109)

where c is called the CFL coefficient of the SSP method. In this framework, a general objective
is to find SSP methods that are higher order accurate, have low computational cost and storage
requirements, and have preferably a large CFL coefficient. Note that the TVB Runge-Kutta
methods can be embedded into this class if the following relaxed notion of stability is assumed

kU n+1 k∞ ≤ (1 + O(∆t)) kU n k∞ (110)

4.4.1. Explicit SSP Runge–Kutta methods. Following Shu and Osher (1988) and the review
articles by Gottlieb, Shu and Tadmor (2001) and Shu (2002), a general m-stage Runge–Kutta

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
48 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

method for integrating (105) in time can be algorithmically represented as

Ũ 0 := U n
l−1
X

Ũ l := αlk Ũ k + βlk ∆tL(Ũ k ) , αlk ≥ 0, l = 1, . . . , m (111)
k=0

U n+1 := Ũ m
Pl−1
To ensure consistency, the additional constraint k=0 αlk = 1 is imposed. If, in addition, all
βlk are assumed to be nonnegative, it is straightforward to see that the method can be written
as a convex (positive weighted) combination of simple forward Euler steps with ∆t replaced
by (βlk /αlk )∆t. From this property, Shu and Osher (1988) concluded the following lemma:
Lemma 5. If the forward Euler method (106) is L∞ -stable subject to the CFL condition (108),
then the Runge–Kutta method (111) with βlk ≥ 0 is SSP, that is, the method is L∞ -stable under
the time step restriction (109) with CFL coefficient
βlk
c = min (112)
l,k αlk

In the case of negative βlk , a similar result can be proven, see (Shu and Osher, 1988).

4.4.2. Optimal second- and third-order nonlinear SSP Runge–Kutta methods. Gottlieb, Shu
and Tadmor (2001) (Proposition 3.1) show that the maximal CFL coefficient for any m-stage,
mth order accurate SSP Runge–Kutta methods is c = 1. Therefore, SSP Runge–Kutta methods
that achieve c = 1 are termed ‘optimal’. Note that this restriction is not true if the number of
stages is higher than the order of accuracy; see Shu (1988).

Optimal second- and third-order nonlinear SSP Runge–Kutta methods are given in Shu and
Osher (1988). The optimal second-order, two-stage nonlinear SSP Runge–Kutta method is
given by
Ũ 0 := U n
1
:= Ũ 0 + ∆t L(Ũ 0 )
Ũ
1 0 1 1 1
U n+1 := Ũ + Ũ + ∆t L(Ũ 1 )
2 2 2
This method corresponds to the well-known method of Heun. Similarly, the optimal third-
order, three-stage nonlinear SSP Runge–Kutta method is given by
Ũ 0 := U n
Ũ 1
:= Ũ 0 + ∆t L(Ũ 0 )
3 0 1 1 1
Ũ 2 := Ũ + Ũ + ∆t L(Ũ 1 )
4 4 4
1 2 2
U n+1 := Ũ 0 + Ũ 2 + ∆tL(Ũ 2 )
3 3 3
Further methods addressing even higher-order accuracy or lower storage requirements are given
in the review articles of Gottlieb, Shu and Tadmor (2001) and Shu (2002) where SSP multistep
methods are also discussed.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 49

5. Finite Volume Methods for Elliptic and Parabolic Problems

As discussed in Section 1, hyperbolic conservation laws are often approximations to physical

problems with small or nearly vanishing viscosity. For some problems, the quantitative solution
effects of these small viscosity or diffusion terms are actually sought. For other problems (e.g.,
porous media problems) the diffusive terms are actually large and the diffusion coefficients (or
permeabilities) are often highly spatially dependent. For a linear convection diffusion problem,
the flux function f in (1) now becomes f (u, ∇u) = −λ∇u + vu, v ∈ Rd where λ > 0 is the
diffusion coefficient, which may be small or large depending on the problem, may depend on
x as in the case of heterogeneous media, or may be replaced by a matrix as in the case of
anisotropic media. The flux through a given edge σ of a control volume K becomes
Z Z
f (u) · nK,σ ds = (−λ∇u · nK,σ + v · nK,σ u) ds (113)
σ σ
R
so that the additional term σ −λ∇u · nK,σ ds must now be discretized.

Emphasis herein is given to analysis of TPFA methods. On the one hand, these methods
have specific monotonicity properties that allow one to cope with noncoercive problems and
irregular data. On the other hand, they assume a discretization of the Laplace operator on
so-called ∆-admissible meshes. Such meshes include the Voronoı̈ tessellation obtained from a
set of vertices of a mesh as depicted in Figure 7. Other methods that apply on general grids

Voronoi dual tessellation

Figure 7. Control volume obtained by the Voronoı̈ dual mesh. Edges of the Voronoı̈ dual are
perpendicular to the edges of the triangulation.

for anisotropic and heterogeneous diffusion operators are outlined in Section 5.2.

5.1. Convergence analysis for the steady state reaction convection diffusion equation

5.1.1. The continuous and discrete problems. Let Ω be an open bounded polygonal subset of
Rd , d = 2 or 3, f ∈ L2 (Ω), v ∈ Rd and b ∈ R, and consider the following steady state linear
reaction convection diffusion equation:
−∆u + div(vu) + bu = f on Ω (114)

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
50 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

together with homogeneous boundary conditions on ∂Ω. A weak formulation of this problem
is given by

1
 Find
Z u ∈ H0 (Ω) such
Z that Z Z
(115)
 ∇u · ∇φ dx + div(vu)φ dx + buφ dx = f φ dx ∀φ ∈ H01 (Ω)
Ω Ω Ω Ω

In order to discretize the above reaction convection diffusion equation, a “∆-admissible” mesh
of Ω is used in discretizing the Laplace operator. A ∆-admissible mesh T is composed of
control volumes with faces (or edges) denoted by E and a set of points P = (xK )K∈T chosen
such that the following orthogonality condition holds: for an interface σKL separating the
cells K and L, the line segment xK xL is orthogonal to this interface; see Figure 8. This
orthogonality condition on the mesh is used to prove the consistency of the flux as detailed
in the following section. Such a family of orthogonal points exists, for instance, in the case
of triangles, rectangles, or Voronoı̈ tessellations; see Eymard et al. (2000) for more details.
Unfortunately, this orthogonality condition is not always satisfied for general meshes so that
certain theoretical results given below are not valid. To address this problem, some schemes
designed specially for general meshes are discussed in Section 5.2.

|σ|

xK

xL
dKL
K
L σKL dK,σ

Figure 8. Notations for a control volume.

The classical flux form of the finite volume scheme is

X
FK,σ (uT ) + b|K|uK = |K|fK ∀K ∈ T (116)
σ∈EK

where |K| denotes the d-dimensional Lebesgue measure of K, fK is the mean value of f over

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 51

K, and FK,σ (uT ) is the numerical flux given by

( |σ| + −
− dKL (uL − uK ) + vK,σ uK − vK,σ uL if σ = σKL
FK,σ (uT ) = |σ| + (117)
− dK,σ (−uK ) + vK,σ uK if σ is an edge of K on ∂Ω
+ −
with vK,σ = |σ|v·nK,σ , vK,σ = max(vK,σ , 0) and vK,σ = − min(vK,σ , 0). One important feature
of FVMs is numerical flux consistency,
R namely that FK,σ (uT ) is a consistent approximation
∗
of the exact flux, F̄K,σ (u) = σ ∇u · nK,σ . Assuming a regular exact solution u, let FK,σ (u)
denote the flux obtained by replacing the discrete unknowns uK by the values of the exact
solution u(xK ). Then, the numerical flux (117) is consistent in the sense
∗
|FK,σ (u) − F̄K,σ (u)| → 0 as hT → 0

with hT = maxK∈T diamK. Note that the discrete Laplace operator ∆T uT resulting from
these numerical diffusion fluxes (see the precise formulas in (120) below) is not consistent in a
finite difference sense. This is demonstrated using a simple Laplace operator example with a
nonuniform one-dimensional mesh in Eymard et al., 2000, Chapter 2, Example 2.1. Even so,
this discretization is consistent in a weak sense. Specifically, thanks to conservation and the
consistency of the flux, it can be proved that the consistency error of the Laplace operator
tends to zero in a L∞ weak-star topology; see Eymard et al., 2000, Chapter 2, Remark 2.9.
The convergence analysis relies on the fact that the discrete Laplace operator is consistent for
some discrete dual H01 norm as addressed in Lemma 8.

Note that the system (116) and (117) always leads to a linear system containing N equations
and N unknowns with N = card(T ). This linear system may be written as
X
AK,L uL = |K|fK for all K ∈ T (118)
L∈T

with
X |σ|
+
AK,K = ( + |σ|vK,σ ) + b|K|
dσ
σ∈EK
|σ| −
AK,L = − − |σ|vK,σ with σ = K|L
dσ
AK,L = 0 if K and L do not share an interface

The finite volume scheme may also be written in the following equivalent weak form:

 Find uT ∈ HT (Ω) such thatZ Z
(119)
 [uT , φ]T + cT (uT , φ) + buT φ dx = f φ dx, ∀φ ∈ HT (Ω)
Ω Ω

where HT (Ω) is the space of piecewise constant functions on the control volumes of T . The
inner product [·, ·]T is defined by
X |σKL | X |σ|
[u, φ]T = (uL − uK )(φL − φK ) + u K φK
dKL dK,σ
σKL ∈Eint σ∈Eext

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
52 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

where Eint (resp. Eext , EK ) denotes the set of edges (or faces) included in Ω (resp. ∂Ω, ∂K),
|σ| the (d − 1)–dimensional Lebesgue measure of σ, dKL the distance between xK and xL
(Figure 8), and dK,σ the distance between xK and σ. In the first summation, σKL denotes
the edge separating the control volumes K and L, and in the last summation, the volume K
is the unique volume in which σ is an edge. To complete the weak formulation, the bilinear
convective form is defined by
" #
X X X
+ − +
cT (uT , φ) = φK (vK,σKL uK − vK,σKL uL ) + vK,σ uK
K∈T σKL ∈EK σ∈EK ∩Eext

Taking φ = 1K in (119), it is easily seen that (119) implies (116). Conversely, let φ ∈ HT (Ω).
Multiplying (116) by φK , summing the resulting equations for all K ∈ T , and reordering the
summations yields (119).

One may also define a discrete Laplace operator in HT in the following way. For φ ∈ HT ,
let ∆T φ ∈ HT be defined by
1 X (d)
(∆T φ)K = − FK,σ (φ) (120a)
|K|
σ∈EK
( |σ|
(d) (φK − φL ) if σ = σKL
FK,σ (φ) = d|σ|KL
(120b)
dKL (φK ) if σ ⊂ ∂Ω

Then, using the conservation property of the flux (FK,σ = −FL,σ if σ = σKL ), it follows that
Z Z
[u, φ]T = − ∆T u φ dx = − u ∆T φ dx, ∀ u, φ ∈ HT (Ω) (121)
Ω Ω

Using the following Poincaré inequality that holds for u ∈ HT (see e.g. Eymard et al., 2000,
Lemma 9.1):
kukL2 (Ω) ≤ diam(Ω) ||u||1,T (122)
a mesh dependent “discrete H01 norm” may be defined using the inner product introduced
above, that is,
1/2
X |σ| X |σ|
1/2
kuk1,T = ([u, u]T ) = (uL − uK )2 + u2K (123)
dKL dK,σ
σKL ∈Eint σ∈Eext

5.1.2. Monotonicity of the scheme. Elliptic PDEs are known to satisfy certain maximum
principle properties. An important property is positivity of the solution for a convection
diffusion problem obtained from Problem (114) with b ≥ 0 and f ≥ 0 together with a
homogeneous Dirichlet boundary condition. Let u be the solution to this problem; then u ≥ 0.
Note that this has been proved even in the noncoercive case; see Droniou, 2002. Moreover, a
more classical result states that if u is the solution to Problem (114) with b = 0, f = 0, and
div v = 0 (in that case the problem is coercive) with the nonhomogeneous Dirichlet boundary
condition
u = ub on ∂Ω (124)

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 53

where ub is a smooth function from Ω to R, then inf ∂Ω ub ≤ u ≤ sup∂Ω ub , see e.g. Evans,
2010, Chapter 6 .

These properties are important in physical applications; for such applications, numerical
schemes are sought that preserve these properties. This is indeed the case for the associated
finite volume scheme, which is given by (116) with the following modified numerical fluxes
(because of the nonhomogeneous boundary condition):
( |σ| + −
− dKL (uL − uK ) + vK,σ uK − vK,σ uL if σ = σKL
FK,σ (uT ) = |σ| + (125)
− dK,σ (ub,σ − uK ) + vK,σ uK if σ is an edge of K on ∂Ω

where ub,σ = ub (xσ ).

Lemma 6 (Positivity) If b ≥ 0, then the matrix of the scheme (116)-(125) is an inverse-
positive matrix, that is, it is invertible and its inverse is nonnegative (all its coefficients are
nonnegative). As a consequence, if fK ≥ 0 for all K ∈ T , and ub,σ ≥ 0 for all σ ∈ Eext , then
the solution (uK )K∈T of (116) and (117) (with b ≥ 0) satisfies uK ≥ 0 for all K ∈ T .

The proof of this result can be found in Eymard, Gallouët and Herbin (2002) in the case
div v = 0. In the case where the sign of div v is not known (and the problem is therefore not
coercive), the result is still true, and relies on the fact that the matrix A of the scheme is
irreducible and diagonally dominant by column, see for example, Fettah and Gallouët, 2013.

The positivity property immediately yields the existence and uniqueness of the solution of the
numerical scheme (119), which can also be proved directly thanks to an L2 a priori estimate
on the approximate solution.
Lemma 7 (Discrete maximum principle) If fK = 0 for all K ∈ T , b = 0 and div v = 0,
then the solution (uK )K∈T of (116)–(125) satisfies
min ub,σ ≤ uK ≤ max ub,σ for all K ∈ T
σ∈Eext σ∈Eext

The maximum principle can be deduced from the fact that the matrix A of the system is
diagonal dominant by row (note that this is only true if div v = 0).

5.1.3. Convergence results. The mathematical analysis of any numerical scheme must address
the question of existence of a solution (which is rather straightforward here since the problem
is linear) and the question of convergence (i.e., “does the approximate solution converge to
the solution of the continuous problem as the mesh size tends to 0?”). A related question
concerns obtaining a rate of convergence through error estimates that are usually dependent
on regularity assumptions on the continuous solution. The proof of the convergence of the
finite volume scheme for a semilinear equation generalizing (114) was first proven in Eymard,
Gallouët, and Herbin, 1999 (see also Eymard et al. (2000)). The result will be stated here
for the linear case and the main steps of the proof explained, since the techniques extend to
nonlinear problems.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
54 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

Under the assumptions of Lemma 6, it is easily seen that the system (119) (resp. (116))
has a unique solution uT ∈ HT (resp. (uK )K∈T ). Let (Tn )n∈N be a sequence of finite volume
discretizations satisfying the orthogonality condition and let hTn be the size of the mesh Tn ,
that is the maximum of the diameters of the control volumes of Tn . In this case, assuming
that hTn → 0 as n → +∞, the corresponding sequence (uTn )n∈N can be shown to converge in
L2 (Ω) to the unique solution of (115). The proof of this result may be decomposed into four
steps:

1. A priori estimates on the approximate solution in the HT norm and the L2 norm are
obtained which yield existence (and uniqueness) of uT solution of the scheme. These
estimates entail the weak convergence of (uTn )n∈N in L2 (Ω), up to a subsequence, to
some ū ∈ L2 (Ω).
2. Strong convergence and regularity of the limit, that is ū ∈ H01 (Ω), are obtained through
a discrete Rellich theorem, described below.
3. The fact that the limit ū is a weak solution of the continuous problem is obtained by a
passage to the limit in the scheme as hT tends to zero.
4. A classical uniqueness argument is then used to show that the whole sequence converges.

Note that there is no need to assume the existence of the solution to the continuous problem
as it is obtained as a by-product of the convergence of the scheme. In the present linear case,
this is not necessary, since existence is well known even in the non-coercive case; see Droniou,
2002. For more complicated nonlinear problems, obtaining the existence of the solution via the
convergence of the numerical scheme can be useful (see e.g. Bouillard et al., 2007).

For the sake of simplicity, these four steps will be detailed in the following paragraphs for
the pure diffusion operator. A sketch of the proof is provided for order h convergence in L2 and
HT norms assuming regularity conditions on the solution, namely u ∈ H 2 (Ω). Note that the
upwind scheme for the convection flux does not lead to any additional difficulty; see Eymard,
Gallouët, and Herbin (1999); Gallouët, Herbin and Vignal (2000). Order 2 convergence in
the L2 norm may be proved for the pure diffusion operator on uniform grids. However, the
same result on triangular meshes, which is observed in numerical experiments, remains an
open problem. Recall that higher convergence rates in weaker norms (including this special
case) are known and proved for most Galerkin methods via duality arguments (the so-called
Aubin-Nitsche lemma, Ciarlet (1991)).

5.1.4. A priori estimate.

Definition 10 (Discrete H −1 norm) Let ψ ∈ HT (Ω), then its discrete H −1 norm is defined
by Z
ψv dx
Ω
kψk−1,T = sup (126)
v∈HT (Ω),v6=0 kvk1,T
Note that, by the discrete Poincaré inequality (122),
kψk−1,T ≤ diam(Ω) kψkL2 (Ω)

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 55

Using the formulation (119) with v = 0 and b = 0, the finite volume scheme can written as
Z
[uT , v]T = f v dx, ∀v ∈ HT (Ω)
Ω

Choosing v = uT , by definition (126), one obtains

||uT ||1,T ≤ kf k−1,T (127)
Taking f = 0 yields the uniqueness (and therefore the existence) of the discrete solution. This
estimate also yields weak convergence of a subsequence of approximate solutions in L2 (Ω).

5.1.5. Convergence theorem. In order to prove strong convergence of the approximate

solutions, some control on solution oscillations is needed. In the finite element framework, the
family of approximate solutions is bounded in H 1 (Ω) so that one may use the Rellich theorem
to obtain compactness in L2 (Ω). This is too restrictive and not used here. However, the Rellich
theorem derives from the Kolmogorov theorem, which gives a necessary and sufficient condition
for a bounded family of Lp (Ω), p < +∞, to be relatively compact. Thus, the Kolmogorov
theorem is an adequate tool in finite volume analysis. In order to use it, some estimates on the
translates of functions of HT (Ω) are needed. Indeed, one may show in a way that is similar
to that of the continuous case (replacing the derivatives by differences) that for any function
v ∈ HT (Ω)
kv(· + η) − vk2L2 (Ω) ≤ |η| (|η| + 4hT )) ||v||21,T ∀η ∈ Rd
This estimate is then used in obtaining the following result:
Theorem 24 (Discrete Rellich theorem) Let (Tn )n∈N be a sequence of finite volume
discretizations satisfying the orthogonality condition, such that hTn → 0. Let (un )n∈N ⊂ L2 (Ω)
such that un ∈ HTn and kun k1,Tn ≤ C where C ∈ R. Then, there exists a subsequence (un )n∈N
and ū ∈ H01 (Ω) such that un → ū in L2 (Ω) as n → +∞.

From the discrete H 1 estimate (127), the above theorem yields the strong convergence of a
subsequence of the approximate solutions in L2 (Ω) to some function ū ∈ H01 (Ω).

5.1.6. Passage to the limit in the scheme. One now needs to show that the limit ū is a
solution to the continuous variational problem. Let (Tn ) be a sequence of discretizations such
that hTn → 0. For each mesh Tn , the finite volume scheme is given by
Z
[uTn , v]Tn = f v dx, ∀v ∈ HTn (Ω) (128)
Ω

Lemma 8 (Consistency of the discrete Laplace operator) Let T be a finite volume

mesh satisfying the orthogonality condition. Denote by PT and ΠT the following interpolation
operators:
PT : C(Ω) → HT (Ω), PT ϕ(x) = ϕ(xK ), ∀x ∈ K, ∀K ∈ T (129)
Z
1
ΠT : L2 (Ω) → HT (Ω), ΠT ϕ(x) = ϕ dx, ∀x ∈ K, ∀K ∈ T (130)
|K| K

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
56 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

For ϕ ∈ Cc∞ (Ω), define the consistency error R∆,T (ϕ) ∈ HT (Ω) on the discrete Laplace
operator by
R∆,T (ϕ) = ∆T PT ϕ − ΠT (∆ϕ)
Then, there exists Cϕ depending only on ϕ such that
kR∆,T (ϕ)k−1,T ≤ Cϕ hT (131)
for any hT sufficiently small (i.e., to say, smaller than the distance between the support of ϕ
and the boundary ∂Ω).

Proof: For ϕ ∈ Cc∞ (Ω), one has

kR∆,T (ϕ)k−1,T = sup X(v)
v∈HT (Ω),kvk1,T =1

with X
X(v) = |K| [(∆T PT ϕ)K vK − (ΠT (∆ϕ))K vK ]
K∈T

For hT small enough, ϕ vanishes in all the control volumes having an edge on the boundary
of the domain so that, by the definition of ∆T , PT , and ΠT
" Z #
X X
X(v) = vK FK,σ (PT ϕ) − ∇ϕ · nK,σ ds
K∈T σ∈EK σ
X (132)
= |σ| RK,σ (ϕ) (vK − vL )
σKL ∈Eint

where RK,σ (ϕ) is the consistency error on the fluxes, defined by

Z
1

RK,σ (ϕ) = FK,σ (PT ϕ) − ∇ϕ · nK,σ ds
|σ| σ
The property of consistency of the fluxes states that for a regular function ϕ, there exists
cϕ ∈ R depending only on ϕ such that
|RK,σ (ϕ)| ≤ cϕ hT
This result, proved in Eymard et al. (2000), is a central argument of the proof. It relies on
the orthogonality condition for the mesh, and is obtained by Taylor expansions. From the
Cauchy–Schwarz inequality and (132), it follows that
X(v) ≤ Cϕ hT kvk1,T
which concludes the proof.

An immediate consequence is the following corollary:

Corollary 2. Let (Tn )n∈N be a family of meshes satisfying the orthogonality property such
that hTn → 0. Let (uTn )n∈N ⊂ L2 (Ω) and ū ∈ H 1 (Ω) such that kuTn k1,T ≤ C, where C ∈ R+
and uTn → ū in L2 (Ω) as n → +∞, then
Z Z
uTn ∆Tn (PTn ϕ) dx → ū ∆ϕ dx as n → +∞, ∀ ϕ ∈ Cc∞ (Ω)
Ω Ω

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 57

A sketch of the proof of convergence of the scheme is now given. Set v = PTn ϕ in (128). From
(121), it follows that Z Z
− uTn ∆Tn (PTn ϕ) dx = f PTn ϕ dx
Ω Ω
R
Using Corollary 2 and the fact that the right-hand side converges to Ω
ϕ dx, passing to the
limit as n → +∞ yields Z Z
− ū∆ϕ dx = f ϕ dx
Ω Ω

Using the discrete Rellich theorem, it was shown previously that ū ∈ H01 (Ω); thus it can be
concluded that ū is the solution to (115).

5.1.7. Error analysis. An error estimate for convection diffusion equations was first obtained
in Herbin, 1995 in the case of continuous data and triangular meshes. It was extended to L2
data, general admissible meshes, and general boundary conditions in Gallouët, Herbin and
Vignal (2000). The key argument for the error analysis is the fact that Lemma (8) still holds
under regularity assumptions for the mesh for φ in H 2 (Ω). From the variational form of the
scheme (128), it follows that
Z
[uTn − PTn u, v]Tn = f v dx − [PTn u, v]Tn ∀v ∈ HTn (Ω)
Ω

where u is the solution to the continuous problem. Integrating the continuous equation
−∆u = f over each control volume, in order to replace the first term of the right-hand side of
the above relation, yields
Z
[uTn − PTn u, v]Tn = R∆,Tn (u)v dx ∀v ∈ HTn (Ω)
Ω

A first-order convergence result in the HT norm then follows from the stability estimate
(127). First-order convergence is also obtained in the L2 norm, thanks to the discrete Poincaré
inequality. Numerical evidence shows that second-order convergence is obtained for several
types of mesh. The mathematical proof is straightforward in the case of rectangular meshes,
since in this case, the consistency error on the flux is of order two. Second-order convergence
for general meshes remains an open question and is the object of ongoing works; see Omnes,
2011 and Droniou and Nataraj, 2016.

5.2. Discretization of anisotropic elliptic problems

A wide variety of schemes has been developed in the last few years for the numerical simulation
of anisotropic diffusion equations on general meshes, see Herbin and Hubert, 2008; Eymard
et al., 2011; Droniou, 2014; Lipnikov et al., 2014; Droniou et al., 2016 and references therein.
The rigorous analysis of these methods is useful to avoid designing numerical schemes that
seem to be well defined and robust, but which in the end do not converge to solutions of the
proper model. Such an example of a nonconverging TP flux method may be found in Faille,
1992b, Chapter III, Section 3.2 in the context of porous media.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
58 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

For simplicity, the following standard elliptic problem is discretized over a convex polygonal
domain Ω ⊂ Rd , d = 2 or 3:
−∇ · Λ∇u = f in Ω (133a)
u(x) = 0 on ∂Ω (133b)
for Λ ∈ Rd×d , a symmetric positive definite (SPD) matrix (assumed constant for the sake of
simplicity) and f ∈ L2 (Ω). The flux balance form of problem (133a) is given by
Z Z
− Λ∇u · nK ds = f dx (134)
∂K K

The solution to problem (133) is understood here in a weak sense. Specifically, there exists a
unique weak solution u ∈ H01 (Ω) satisfying
Z Z
Λ∇u(x) · ∇v(x) = g(x)v(x)dx ∀v ∈ H01 (Ω) (135)
Ω Ω

For heterogeneous anisotropic problems such as (133), a variety of recently developed

schemes are based on the definition of a numerical flux which uses the cell centered unknowns
together with some other unknowns that are located at the edges of the mesh. Some of these
methods permit the explicit elimination of the edge unknowns while ensuring the following
properties:

(P1) The schemes must apply on any type of grid: conforming or nonconforming, 2D or 3D
(see for instance the frameworks of kinetic formulations or financial mathematics), and
made with control volumes that are only assumed to be polyhedral (the boundary of
each control volume is a finite union of subsets of hyperplanes).
(P2) The matrices of the generated linear systems are expected to be sparse, symmetric, and
positive definite.
(P3) One should able to prove the convergence of the discrete solution and an associated
gradient to the solution of the continuous problem and its gradient with no regularity
assumption on the solution of the continuous problem as well as show error estimates if
the continuous solution is regular enough.

The idea of some of the schemes presented below is to find an approximation of the solution of
(133) by setting up a system of discrete equations for a family of values ((uK )K∈T , (uσ )σ∈E )
in the control volumes and on the interfaces. The values uσ on the interfaces are introduced
so as to allow for a natural consistent approximation of the normal fluxes in the case of an
anisotropic operator and a general (possibly nonconforming) mesh. Counting both cell and
interface unknowns yields card(T ) + card(E) as the total number of unknowns. But in some
cases, the unknowns (uσ )σ∈E can be eliminated. In fact, the edge unknowns can always be
eliminated by a suitable procedure (discussed below), but the resulting scheme is not always
optimal in terms of flux accuracy or robustness.

Following the idea of the finite volume framework, consider the flux balance form (134). The
integral boundary ∂K is decomposed as the sum of integrals over the interfaces of E which

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 59

compose ∂K
X  Z  Z
− Λ∇u(x) · nK,σ ds = f dx
σ∈EK σ K
R
The flux − σ Λ∇u(x) · nK,σ ds is approximated by a function FK,σ (u) of the values
((uK )K∈T , (uσ )σ∈E ) at the centers of the control volumes and on interfaces.

5.2.1. TPFA scheme. Assume first that Λ = λ Id where λ is a given scalar function so that
the diffusion is heterogeneous and isotropic. This is a problem of great interest in petroleum
reservoir engineering where the permeability of the medium is generally constant on each cell of
the grid (but may vary a lot from one cell to another). The so-called TPFA is a straightforward
extension of the finite volume scheme defined by the numerical fluxes (120b) for the isotropic
homogeneous problem to the heterogeneous diffusion problem. Begin by setting
uσ − uK
FK,σ (u) = −|σ|λK ∀K ∈ T , ∀σ ∈ EK , (136)
dK,σ
where λK is the mean value of λ on K. The unknowns uσ may be eliminated by using
conservation of the numerical flux, that is, FL,σ = −FK,σ if σ = σKL ⊂ Ω. If xL ∈
6 σ;
this yields the following value of uσ
1 λK,σ λL,σ

uσ = λK,σ λL,σ
uK + uL
+ dK,σ dL,σ
dK,σ dL,σ

If xL ∈ σ, then uσ = uK . Plugging these expressions in the flux (136) yields the following
numerical flux involving the harmonic mean of λ:
FK,σ = −τσ (uL − uK ) if σ ∈ Eint , σ = σKL
where
λK,σ λL,σ
τσ = |σ| if xσ 6= xK and xσ 6= xL
λK,σ dL,σ + λL,σ dK,σ
and
λK,σ
τσ = |σ| if xσ 6= xK and xσ = xL
dK,σ
In practice using this harmonic mean λ in the expression of the flux gives very good results.

This scheme is ideal for the ∆-admissible meshes introduced in Section 5.1. For instance, in
the case of rectangles, the properties obtained for the Laplace operator remain valid, that is,
the scheme is robust and converges to the correct solution of the problem. However, geological
grids often consist of very distorted (and sometimes degenerate) parallelograms for which
the orthogonality condition fails to hold. In this case, the numerical diffusion flux becomes
inconsistent. If the number of interfaces for which the orthogonality condition fails is large,
the approximate solution may be very far from the solution of the problem. In this case, a
natural idea to get a consistent numerical flux at an interface is to compute it by using more
points than just the two points of the neighboring cells. This is the principle used in the
methods developed in the works of Eymard et al., 2000 (see Section 3.1.1) and Faille, 1992a

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
60 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

where the values at the points needed to define the normal gradient were interpolated from the
cell unknowns. Alternatively, in (Coudière et al., 1999) a gradient is reconstructed on diamond
cells from the cell unknowns. In both cases, the scheme can only be proved to be stable under
some geometrical assumptions on the mesh. Consequently, several other methods have been
designed in recent years. Three of them are described below and a brief overview of some
others is also given.

5.2.2. Multiple Point Flux Approximation schemes. The class of multiple point flux
approximation (MPFA) methods was initiated in the works of Edwards and Rogers, 1994;
Edwards and Rogers, 1998; Aavatsmark et al., 1998a; Aavatsmark et al., 1998b. The idea
is to compute consistent fluxes by introducing an approximate gradient on subcells around
each vertex using some additional edge unknowns, which are then eliminated by exploiting
flux conservation. Numerous variants of MPFA schemes have been derived. The MPFA “O-
scheme” for a rectangular mesh is particularly simple to describe. For K ∈ T and v a vertex

ev
xσe σ v

VK,v σv

xσ
xK
σ
K nK,σ

v′

Figure 9. Notations for MPFA-O schemes defined on a rectangular mesh.

of K, let VK,v be the parallelepiped polyhedron whose faces are parallel to the faces of K
and that has xK , xσ and v as vertices, where xσ denotes the centroid of σ ∈ E. The discrete
flux balance over a cell K is established, but instead of computing the numerical flux over the
whole interface σ, two numerical fluxes FK,σ,v (u) and FK,σ,w (u) over the half edges σv and σw
(see Figure 9) are introduced. The discrete balance equation thus reads
X X Z
FK,σ,v (u) = f (x)dx
σ∈EK v∈Vσ K

The numerical fluxes over these half edges are computed by assuming an approximate solution
u, which is piecewise affine on the subcells VK,v so that its gradient ∇VK,v u is piecewise constant.
Then, observing that ∇VK,v u · (xK − xσ ) = uK − uσ and ∇VK,v u · (xK − xσe ) = uK − uσe , the
constant gradient of u is computed as
1 X
∇VK,v u = |σv |(vσ − vK )nK,σ
|VK,v |
σ∈EK,v

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 61

where EK,v is the set of interfaces so that the numerical flux through the half edge σv can be
computed as
1
FK,σ,v (u) = |σv |Λ∇VK,v u · nK,σ = |σv |2 (vσ − vK )AnK,σ · nK,σ
|VK,v |
In order to recover the complete scheme written only in terms of the cell unknowns (uK )K∈T ,
it remains to eliminate the face unknowns (uσ )σ∈Eint . This is accomplished by requiring the
conservation of the numerical fluxes FK,σ,v (u) through the half-edges σv . The following linear
system then needs to be solved for each vertex of the mesh:

FK,σ,v (u) = −FL,σ,v (u), ∀σ ∈ Ev , for any σ = K|L

Note that this local system is always invertible in the case of rectangular meshes considered
here, but not so for general distorted meshes. The MPFA O-scheme has been introduced and
used in the context of oil reservoir simulation. However, for highly nonconforming meshes, the
local systems may be difficult to invert and the numerical fluxes may become unrealistic.

The MPFA O-scheme may be constructed for general meshes; see Aavatsmark et al., 1998b.
Other variants have been proposed such as the so-called L scheme in Aavatsmark et al., 2008,
the U scheme in Aavatsmark et al., 1998b, and the scheme given in Agelas et al., 2010. See
Droniou, 2014 for a review of these methods and their properties.

Remark 1.[Gradient schemes] The numerical fluxes of the O-scheme were computed by
constructing an approximate gradient on subcells on the grid. The idea of using a discrete
gradient can be pursued further by using it in an approximate weak formulation of the scheme.
These so-called gradient schemes have been designed to deal with anisotropic elliptic problems
on polytope meshes in Eymard et al., 2012. The idea is to discretize the weak form of equation
(145a), thanks to a gradient discretization, which consists of a set of discrete unknowns, a
reconstruction operator, and a discrete gradient, which are both constructed from the set of
discrete unknowns such that the L2 norm of the discrete gradient is a norm on the set of
discrete unknowns. In order for the gradient scheme to be convergent for the linear problem,
the gradient discretization should satisfy the following properties:

• Coercivity ensures uniform discrete Poincaré inequalities for the family of gradient
discretizations. This is essential in obtaining a priori estimates for the solutions to
gradient schemes.
• Consistency ensures that the family of gradient discretizations “covers” the whole energy
space of the model (e.g., H01 (Ω) for the linear equation (133)).
• Limit-conformity ensures that the family of gradient and function reconstructions
asymptotically satisfies the Stokes formula.

See Eymard et al., 2012; Droniou et al., 2010; Droniou et al., 2016; Droniou et al., 2014 for more
precise definitions. Even though gradient schemes are not, in general, finite volume schemes (in
fact, they include conforming and nonconforming finite element schemes, mimetic methods,
and other schemes), some particular finite volume schemes are gradient schemes. Furthermore,

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
62 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

several recently developed schemes are in fact gradient schemes. This is the case for the O-
scheme on rectangular and simplicial meshes. This is also the case for some forms of thediscrete
duality finite volume (DDFV) scheme discussed next.

5.2.3. DDFV schemes. The design principle of DDFV schemes Hermeline (2000), Hermeline,
2003; Domelevo and Omnes, 2005; Andreianov et al., 2007; Boyer and Hubert, 2008; Hermeline,
2007; Andreianov et al., 2008; Coudière and Hubert, 2011 is again to introduce some discrete
gradients to compute the numerical diffusion flux on a face (or edge) of the mesh. However,
now the additional unknowns that are used to compute the fluxes are located at the vertices
rather than at the faces (or edges of the mesh). The approximate gradients are chosen piecewise
constant on the so-called “diamond cells” (dotted area in Figure 10) and they are computed
using the directions given by the vertices of the diamond cell in Figure 10, which are linearly
independent. The discrete gradient on the diamond cell Dσ associated with the interface

Figure 10. 2D DDFV: Control volumes K and L and diamond cell Dσ for the gradient reconstruction.

σ = σKL , whose vertices are v and w, is then computed such that

∇σ u · (xK − xL ) = uK − uL and ∇σ u · (v − w) = uv − uw

The approximate numerical flux through an interface σ belonging to the boundary of K can
then be defined as
FK,σ (u) = −|σ|Λσ ∇σ u · nK,σ
where Λσ is the mean value of Λ over the diamond cell Dσ . The discrete unknowns are the
values (uK )K∈T and (uv )v∈V where V is the set of vertices of the mesh.

In order to have as many equations as unknowns, the discrete balance equations are written
on the primal cells K as in the TP scheme and also on dual cells Pv associated with the
e : xK xL is easily computed
vertices, such as shown in Figure 10. The flux through the edge σ

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 63

using the discrete gradient

FPv ,eσ (u) = −|e
σ |Λσ ∇σ u · nPv ,eσ
The linear system to be solved is given by
X Z X Z
FK,σ (u) = f dx ∀K ∈ T and FPv ,eσ (u) = f dx ∀v ∈ V (137)
σ∈T K v∈V Pv

and is of order NT + NV , where NT is the number of cells of the mesh and NV the number
of vertices. The name DDFV was chosen as a reminder to the fact that the system may also
be written using a discrete divergence operator that can be deduced from the above discrete
gradient using a discrete Stokes formula. This property is used in the convergence analysis that
was performed in two space dimensions for a number of problems in Domelevo and Omnes,
2005; Boyer and Hubert, 2008; Andreianov et al., 2007.

The extension of the DDFV method to the three-dimensional setting was developed more
recently. Two main approaches exist: (i) the CeVe-DDFV method, which uses cell and vertex
unknowns as described in Hermeline, 2009; Coudière et al., 2009; Andreianov et al., 2010,
and (ii) the CeVeFE-DDFV method, which uses cell, vertex, faces and edges unknowns as
described in Coudière and Hubert, 2011; Coudière et al., 2011. The coercivity properties of the
two methods differ: the CeVe-DDFV method does not seem to be unconditionally coercive on
generic meshes, whereas the CeVeFE-DDFV method is unconditionally coercive; see Droniou,
2014. In fact, this latter method is a gradient scheme as described in Remark 1. The original
construction of this scheme necessitates four meshes. However, it can also be described using
only one mesh; see for example Droniou et al., 2016 for such a description, which relies on the
definition of a discrete gradient on the octahedral cells VK,v depicted in Figure 11. As shown in

v
F

D B
VK,v

xK
A C
A C
F
B=O
E

D
E

Figure 11. (a) Octahedral cell K for the CeVeFE-DDFV scheme. (b) Illustration of VK,v . (c)
Construction of a degenerate octahedron from a non-conforming hexahedral mesh in a heterogeneous
medium (CDEF is the intersection of the boundaries of two nonmatching hexahedral cells).

Andreianov et al., 2007 (Section IX.B) for 2D DDFV methods, the other three meshes may be
reconstructed from the “diamond” mesh, but they are not really needed to define the scheme

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
64 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

(nor to implement it). Note that the octahedral cells may be degenerate. This is the case, for
instance, when working with nonconforming hexahedral meshes used when modeling flows in
heterogeneous porous media.

5.2.4. HHM schemes. In the past decades, several schemes have been developed for elliptic
equations so as to satisfy some form of calculus formula at the discrete level. These schemes
are called mimetic finite difference (MFD) or compatible discrete operator (CDO) schemes.
See Lipnikov et al., 2014 for a review of MFD methods and Bonelle and Ern, 2014; Bonelle
et al., 2015 and references therein for CDO methods. Contrary to DDFV methods, which
construct discrete operators and duality products to satisfy fully discrete calculus formulas,
MFD/CDO methods are based on discrete operators satisfying a Stokes formula, which involves
both continuous and discrete functions. Depending on the choice of the location of the main
geometrical entities attached to the degrees of freedom (faces or vertices), two different
MFD/CDO families exist.

A first MFD method, hereafter called hybrid mimetic finite difference (HMFD), is obtained
by using a mixed form of (133) with fluxes through the mesh faces as initial unknowns, that
is, by introducing q = Λ∇u so that −∇ · q = f , and then discretizing this set of two equations.
The resulting scheme was proved in Droniou et al., 2010 to be embedded in a slightly larger
family, which also contains the hybrid finite volume (HFV) method in Eymard et al., 2010a
and the mixed finite volume (MFV) method in Droniou and Eymard, 2006. The schemes of
this family are called hybrid mimetic mixed (HMM) schemes. Each scheme in this family can
be written in three different ways depending on the approach considered: HMFD, HFV, or
MFV. The HFV formulation of a HMM scheme is very close to the weak formulation (135) of
the elliptic PDE. It consists of a weak formulation with a discrete gradient and a stabilization
term (bilinear form on (u, v)), see Eymard et al., 2010a. It was proved in Droniou et al.,
2013 that the discrete gradient can be modified to include the stabilization terms. Thus, all
HMM methods (and therefore also all HMFD methods) are part of the gradient scheme family
described in Remark 1.

5.2.5. Other schemes and topics of interest. Another type of scheme for diffusion problems
is the FVE method, which is based on finite element spaces with vertex unknowns and
flux balances on dual meshes around vertices; see Cai and McCormick, 1990; Cai, 1991; Cai
et al., 1991; Süli, 1991; Lazarov, Michev and Vassilevsky, 1996; Chatzipantelidis and Lazarov,
2005; Chatzipantelidis et al., 2013 and references therein. The FVE method can be recast in
Petrov–Galerkin form using a piecewise constant test space together with a conforming trial
space. To formulate and analyze the Petrov–Galerkin representation, two tessellations of Ω
are considered: a triangulation T with simplicial elements K ∈ T and a dual tessellation T ∗
with control volumes T ∈ T ∗ . In the class of conforming trial space methods such as the FVE
method, a globally continuous, piecewise p-order polynomial trial space with zero trace value
on the physical domain boundary is constructed

Xh = {v ∈ C 0 (Ω) v|K ∈ Pp (K), ∀K ∈ T and v|∂Ω = 0}

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 65

Median dual tessellation

Figure 12. Control volume for the finite volume element method formed from median dual segments
in each triangle.

using nodal Lagrange elements on the simplicial mesh. A dual tessellation T ∗ of the Lagrange
element is then constructed. See Figure 12 which shows a linear Lagrange element with dual
tessellation. These dual tessellated regions form control volumes for the FVM. The tessellation
technique extends to higher-order Lagrange elements in a straightforward manner. A piecewise
constant test space is then constructed using T ∗

Yh = {v v|T ∈ χ(T ), ∀T ∈ T ∗ }

where χ(T ) is a characteristic function in the control volume T . The FVE discretization of
(133a) then yields the following Petrov–Galerkin formulation: Find uh ∈ Xh such that
X Z Z 
wh Λ∇uh · n ds + wh f dx = 0, ∀wh ∈ Yh (138)
∀T ∈T ∗ ∂T T

The analysis of (138) by Ewing, Lin and Lin (2002) using linear elements gives an a priori
estimate in an L2 norm that requires the least amount of solution regularity when compared
to previous methods of analysis.
Theorem 25. (FVE a priori error estimate, Ewing, Lin and Lin (2002)) Assume a 2-D
quasi-uniform triangulation T with dual tessellation T ∗ such that ∃C > 0 satisfying
∗
C −1 h2 ≤ |T | ≤ Ch2 , ∀T ∈ T

Assume that u and uh are solutions of (133a) and (138) respectively with u ∈ H 2 (Ω),
f ∈ H β , (0 ≤ β ≤ 1). Then ∃C 0 > 0 such that the a priori estimate holds


ku − uh kL2 (Ω) ≤ C 0 h2 kukH 2 (Ω) + h1+β kf kH β (Ω) (139)

Unlike the finite element method, the error estimate (139) reveals that optimal order
convergence is obtained only if f ∈ H β with β ≥ 1. Moreover, numerical results show that the
source term regularity cannot be reduced without deteriorating the measured convergence rate.
Optimal convergence rates are also shown for the nonconforming Crouzeix–Raviart element
based FVM analyzed by Chatzipantelidis (1999) for u ∈ H 2 (Ω) and f ∈ H 1 (Ω).

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
66 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

Other FVMs or related schemes for elliptic boundary value problems have been proposed
and analyzed under a variety of names: box methods in Bank and Rose (1987), Bank and
Rose (1987),Croisille and Greff, 2005, Greff, 2007, covolume methods in Chou and Li (2000),
diamond cell methods in Coudière et al., 1999; Kútik and Mikula, 2015, and finite volume based
on the Crouzeix–Raviart element in Chatzipantelidis (1999); Ewing, Lin and Lin (2002).

A posteriori error estimates have been obtained for a number of the above mentioned
schemes; see for example, Lazarov and Tomov, 2002; Omnes, 2008; Di Pietro et al., 2011;
Ern and Vohralı́k, 2011; Chen and Wang, 2013; Cancès et al., 2014; Chen and Gunzburger,
2014; Arbogast et al., 2014; Erath, 2015. See the review article by Di Pietro and Vohralı́k,
2014 and references therein for more on this subject.

A number of recent methods for anisotropic problems can be found in Herbin and Hubert,
2008; Eymard et al., 2011 which give a numerical comparison of the methods for several 2D and
3D application benchmarks using various types of meshes. See also the recent review articles
in Droniou, 2014; Droniou, 2014; Di Pietro and Vohralı́k, 2014.

5.2.6. The transmissivity structure. In the theoretical study of linear or nonlinear PDEs, it
is sometimes useful to use a nonlinear function of the unknowns as a test function in the weak
formulation. This may be the case, for instance, to establish some of the properties of the
solutions or in order to show the existence of solutions to the systems. Here are four examples
where this approach is used to show

1. positivity of the solution for a linear elliptic problem with nonnegative right-hand side;
2. L∞ estimates for the solution of an elliptic problem for a right-hand side in Lp , p > d2 ;
3. existence of solutions for elliptic problems in a nonvariational framework, for instance,
if the right-hand side is only integrable; see Boccardo and Gallouët, 1989;
4. existence of solutions for noncoercive convection diffusion problems (and positivity if the
right-hand-side is nonnegative); see Droniou, 2002.

See Gallouët, 2007 for more details and references on these four examples. In the theoretical
study of discretization schemes for the same equations, the same technique (i.e., using a
nonlinear function of the unknown in a weak formulation) is also often used. However, in
the current state of the art, this technique is successful only if the schemes are of the form
X
τK,L (uK − uL ) = rhs
L∈N(K)

with τK,L ≥ 0. This is the case for the TPFA scheme on admissible meshes (see e.g., Gallouët
and Herbin, 1999; Droniou and Gallouët, 2002; Droniou, 2003; Droniou et al., 2003; Chainais-
Hillairet and Droniou, 2011) for noncoercive problems and/or problems with nonregular right-
hand-side. For general schemes, this is an open problem.

Some results (e.g., positivity, maximum principle) have been obtained for finite volume
schemes using more general meshes but the price to pay is that the scheme is nonlinear, even

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 67

in the case of a linear problem; see for example Le Potier, 2005; Le Potier, 2008; Le Potier,
2009; Genty and Le Potier, 2011; Droniou and Le Potier, 2011; Cancès et al., 2013 for the
original schemes and Droniou, 2014 for a review.

5.3. The parabolic case

5.3.1. The continuous problem. Next, consider a time-dependent convection diffusion

equation. Let T > 0, u0 ∈ L2 (Ω) and v ∈ Rd be given. The time-dependent convection
diffusion problem, u : Ω × [0, T ] → R, is given by

 ∂t u + div(vu) − ∆u = 0 in Ω × (0, T )
u = 0 in ∂Ω × (0, T ) (140)

u(., 0) = u0 in Ω
For the sake of simplicity, consider the case v = 0. A weak formulation of this problem is

−1

 Find u ∈ L2 (0, T ; H01 (Ω)) 2
Z such that ∂t u ∈ L (0, T ; H (Ω)) and

< ∂t u, ϕ >H −1 ,H01 + ∇u(x, ·) · ∇ϕ(x, ·)dx = 0 ∀ ϕ ∈ H01 (Ω) a.e. in (0, T ) (141)

 Ω
 u(·, 0) = u
0

Here, the duality product < ·, · >H −1 ,H01 is defined by

Z T
< ∂t u, ϕ >H −1 ,H01 = − u(·, t)∂t (·, t)ϕdt ∈ H01 (Ω)
0

Hence, identifying the space L2 (Ω) with its dual space (L2 (Ω))0 , the statement ∂t u ∈
L2 (0, T ; H −1 (Ω)) is to be understood as requiring that there exists v ∈ L2 (0, T ; H −1 (Ω))
such that Z T Z T
− u(·, t)∂t (·, t)ϕdt = v(·, t)ϕ(·, t)
0 0
| {z } | {z }
∈H01 (Ω) ∈H −1 (Ω)

which makes sense since H01 (Ω) can then be identified to its dual space H −1 (Ω). Note that if
u satisfies (141), then u ∈ C([0, T ], L2 (Ω)).

As in the steady state case, some estimates on u are obtained in order to get compactness
properties despite the lack of regularity of the approximate finite volume solution. In
the continuous framework, the natural estimates are in L2 (0, T ; H01 (Ω)) for u and in
L2 (0, T ; H −1 (Ω)) for ∂t u. These estimates give compactness in L2 (0, T ; L2 (Ω)) for a sequence
of approximate solutions (using, for instance, the Faedo–Galerkin approximation), which
proves the existence of a solution to (141). It is reasonable to therefore look for the same kind of
estimates in the discrete framework, which will also yield the compactness in L2 (0, T ; L2 (Ω))
of a sequence of approximate finite volume solutions.

5.3.2. The finite volume scheme. Let T be an admissible mesh of Ω, in the sense introduced
in Section 5.1, and let δt = T /M be the (uniform) time step. The finite volume scheme with

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
68 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

implicit Euler discretization in time is given by

 X
 un+1 − unK

 |K|
K
+ FK,σ (un+1
T )=0 0≤n≤M −1
δt
Z σ∈EK (142)

 u0K = 1
 u0 (x)dx
|K| K
|σ| n+1 −
with FK,σ (un+1
T )=− (u n+1
− uK +
) + vK,σ un+1 − vK,σ un+1
L .
dKL L K

The existence and uniqueness of a solution (unK )n∈N to (142) is easily deduced from
the steady-state case. Denote by HD (Ω × (0, T )) the set of functions of L2 (Ω × (0, T )),
which are piecewise constant on the subsets K × [tn , tn+1 ). Define an approximate solution
uD ∈ HD (Ω × (0, T )) by uD (x, t) = unK , ∀x ∈ K, ∀t ∈ [tn , tn+1 ). Using a variational technique
similar to the way the estimate (127) is established in the steady-state case, the following a
priori estimates on uD may be obtained:
kuD kL∞ (0,T ; L2 (Ω)) ≤ C (143)
and
M
X
δt kuD (·, tn )k21,T ≤ C (144)
n=1

where k · k1,T is the discrete H01 norm defined by (123) and C depends only on the initial
condition. Using equation (142), a discrete H −1 estimate on the discrete time derivative is
derived. For n = 0, . . . , M − 1, let δtn u = un+1δt−un ∈ HT (Ω); then the following estimate is
obtained:
MX−1
δt kδtn uk−1,T ≤ C
n=1

where the discrete dual norm k · k−1,T is defined by (126), and C depends only on the initial
condition. The compactness of a sequence of approximations may then be deduced from a
generalization of the Aubin–Simon lemma for sequences of embedded subspaces that can be
stated in the general Lp setting and used in the context of less regular solutions or nonlinear
problems.
Lemma 9 (Generalized Aubin–Simon lemma) Let 1 ≤ p < +∞. Let B be a Banach
space, (X` )`∈N be a sequence of Banach spaces included in B, and (Y` )`∈N be a sequence of
Banach spaces. Assume that the sequence (X` , Y` )`∈N satisfies the following conditions:

1. The sequence (X` )`∈N is such that any sequence (u` )`∈N , satisfying u` ∈ X` (for all
` ∈ N) and (ku` kX` )`∈N bounded, is relatively compact in B.
2. X` ⊂ Y` (for all ` ∈ N) and if the sequence (u` )`∈N is such that u` ∈ X` (for all ` ∈ N),
(ku` kX` )`∈N bounded and ku` kY` → 0 (as ` → +∞), then any subsequence converging in
B converges (in B) to 0.

Let T > 0 and (u` )`∈N be a sequence of Lp ((0, T ), B) satisfying the following conditions:

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 69

PN
1. For all ` ∈ N, there exists N ∈ N? and k1 , . . . , kN in R?+ such that i=1 ki = T and
u` (t) = vi for t ∈ (ti−1 , ti ), i ∈ {1, . . . , N }, t0 = 0, ti = ti−1 + ki , vi ∈ X` . (Of course,
the values N , ki and vi depend on `.)
2. The sequence (u` )`∈N is bounded in Lp ((0, T ), B).
3. The sequence (ku` kL1 ((0,T ),X` ) )`∈N is bounded.
4. The sequence (kδt u` kLp ((0,T ),Y` ) )`∈N is bounded, where the function δt u` is defined a.e.
by
vi − vi−1
δt u` (t) = for t ∈ (ti−1 , ti )
ki

Then there exists u ∈ Lp ((0, T ), B) such that, up to a subsequence, u` → u in Lp ((0, T ), B).

See, for example, Chénier et al., 2015 for its proof. The convergence in L2 (0, T ; L2 (Ω)) of uD
to some function ū ∈ L2 (0, T ; H01 (Ω)) is then obtained. As in the elliptic case, a passage to
the limit in the scheme then yields that ū = u, weak solution of (141). This analysis may be
generalized to the case of nonhomogeneous Dirichlet boundary conditions; see Bouillard et al.,
2007.

Error estimates for finite volume schemes applied to linear parabolic equations may also be
obtained; see for example, Bradji, 2008; Bradji and Fuhrmann, 2010; Chatzipantelidis et al.,
2013. Note that in the case of parabolic equations with L1 data, only the TPFA scheme has
been analyzed in Gallouët et al., 2012. As in the elliptic case, the analysis is based upon some
test functions that are nonlinear functions of the unknowns, as in the continuous case. As
already mentioned in Section 5.2.6, the extension of these results to numerical schemes on
general grids is still an open problem. Nonlinear parabolic systems have also been analyzed,
with several areas of applications: electrochemistry in Bradji and Herbin, 2008, porous media
in Bouillard et al., 2007; Chainais-Hillairet and Droniou, 2007, image processing in Mikula and
Ramarosy, 2001; Frolkovič and Mikula, 2007; Drblı́ková and Mikula, 2007; Drblı́ková et al.,
2009, to cite only a few.

5.3.3. Nonlinear conservation laws including diffusion terms. As a final scalar PDE model,
consider the nonlinear first-order conservation law with an added second-order Laplacian
diffusion term
∂t u + ∇ · f (u) − ε4u = 0 in Rd × R+ (145a)
u(x, 0) = u0 in Rd (145b)
Here, u(x, t) : Rd × R+ → R denotes the dependent solution variable, f ∈ C 1 (R, Rd ) the
hyperbolic flux, and ε ≥ 0 a small diffusion coefficient. An application of the divergence and
Gauss theorems to (145a) integrated in a region K yields the following integral conservation
law form Z Z Z
∂
u dx + f (u) · nK ds − ε∇u · nK ds = 0 (146)
∂t K ∂K ∂K
Using the finite volume approximation of the diffusion operator of the previous section, the
fully discrete form (31) of Section 3 is extended to the integral conservation law (146) by the

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
70 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

introduction of a numerical diffusion flux function FK,σ (uh ) for a control volume K ∈ T such
that Z X
ε∇u · nK ds ≈ FK,σ (uh )
∂K σ∈EK

where the discrete diffusive flux is chosen similar to the linear parabolic case, as FK,σ (um h ) :=
|σ| m m
dKL (uL − uK ). When combined with the general finite volume formulation (31) for hyperbolic
conservation laws, the following fully discrete scheme is produced
∆tn X

un+1
K = unK − g(unK , unL ; σ) |σ| − FK,σ (um
h ) , ∀K ∈ T (147)
|K|
σ∈EK

In this equation, the index m may be chosen either as n or n + 1, corresponding to an explicit

or implicit discretization.

Stability analysis reveals a CFL-like stability condition for the explicit scheme [choice m = n
in (147) ]
α3 (hnmin )2
∆tn ≤
αLg hnmin + ε
where Lg denotes the Lipschitz constant of the hyperbolic numerical flux, α is a positive mesh-
dependent parameter, and  is the diffusion coefficient. In constructing this bound, a certain
form of shape regularity is assumed such that there exists an α > 0 such that for all j, k with
hK ≡ diam(K)
αh2K ≤ |K| α|∂K| ≤ hK αhK ≤ dKL (148)
Thus, ∆tn is of the order h2 for large ε and of the order h for ε ≤ h. In cases where the diffusion
coefficient is larger than the mesh size, it is advisable to use an implicit scheme (m = n + 1).
In this latter situation, no time step restriction has to be imposed; see Eymard et al. (2002)
and Ohlberger (2001b).

In order to demonstrate the main difficulties when analyzing convection-dominated

problems, consider the following result from Feistauer et al. (1999) for a homogeneous diffusive
boundary value problem. In this work, a MFV, finite element method sharing similarities with
the methods described above is used to obtain a numerical approximation uh of the exact
solution u. Using typical energy-based techniques, they prove the following error bound.
Theorem 26. For initial data u0 ∈ L∞ (R2 ) ∩ W 1,2 (R2 ) and τ > 0 there exist constants
c1 , c2 > 0 independent of ε such that
ku(·, τ ) − uh (·, τ )kL2 (Ω) ≤ c1 hec2 τ / (149)

This estimate is fundamentally different from estimates for the purely hyperbolic problems
of Sections 3 and 4. Specifically, this result shows how the estimate strongly depends on the
small parameter ε, ultimately becoming unbounded as ε tends to zero.

In the context of convection dominated or degenerate parabolic equations, Kruzkov

techniques have been used by Carrillo (1999) and Karlsen and Risebro (2000) in proving the

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 71

uniqueness and stability of solutions. Utilizing these techniques, convergence of finite volume
schemes (uniform with respect to ε → 0) was proved in Eymard et al. (2002) and error estimates
were obtained for viscous approximations in Evje and Karlsen (2002) and Eymard, Gallouët
and Herbin (2002). In Ohlberger (2001a, 2001b) uniform a posteriori error estimates suitable
for adaptive meshing are given, see also Ohlberger and Rohde, 2002 for the case of weakly
coupled systems. Using the theory of nonlinear semigroups, continuous dependence results
were also obtained in Cockburn and Gripenberg (1999) (see Cockburn (2003) for a review).
FVMs for nonlinear systems of parabolic equations and degenerate parabolic equations have
received a lot of attention in the recent years in a wide area of applications such as porous
media, corrosion models, population models, and so on as discussed in Chainais-Hillairet and
Peng, 2004; Chainais-Hillairet and Droniou, 2007; Andreianov et al., 2010; Andreianov et al.,
2011; Chainais-Hillairet et al., 2011; Andreianov et al., 2011; Bessemoulin and Filbet, 2012;
Andreianov et al., 2010; Chainais-Hillairet et al., 2013; Bessemoulin et al., 2014; Chainais-
Hillairet et al., 2015 and references therein.

6. Advanced Topics

6.1. Extension to systems of nonlinear hyperbolic conservation laws

The early widespread use of FVMs comes, in part, from the relative ease in algorithmically
extending the numerical discretization of Section 3 to systems of nonlinear hyperbolic
conservation laws of the form
∂t u + ∇ · f (u) = 0 in Rd × R+ (150a)
u(x, 0) = u0 (x) in Rd (150b)
In these equations, u(x, t) : Rd ×R+ → Rm denotes the vector of dependent solution variables,
f : Rm → Rm×d denotes the flux vector, and u0 : Rd → Rm denotes the initial data vector
function. It is also assumed that the conservation law system (150) possesses an additional
scalar entropy equation in divergence form with convex entropy–entropy flux pair {U, F },
U : Rm 7→ R convex and F : Rm → R1×d
(U (u))t + ∇ · (F (u)) ≤ 0 (151)
with the right-hand side equal to zero for classical (smooth) solutions. As discussed further in
Section 6.1.3, the existence of this entropy extension equation with convex entropy function U
is sufficient to deduce that (150a) can be recast in symmetric hyperbolic form under a change
of variables u 7→ v where v ∈ Rm are the so-called entropy variables for the conservation
law system. This also implies hyperbolicity in the sense that eigenvalues of the flux jacobian
∂f
matrix A(u; n) ≡ ∂u · n are real (but not necessarily distinct). This eigenstructure is used in
the construction of upwind numerical flux functions in Sections 6.1.1 and 6.1.2.

6.1.1. Numerical flux functions for linear hyperbolic systems. The main task in extending
FVMs to systems of nonlinear conservation laws is the construction of a suitable numerical flux

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
72 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

function. To gain insight into this task, consider the one-dimensional linear Cauchy problem
∂t u + ∂x (A u) = 0 in R × R+
u(x, 0) = u0 (x) in R (152)
where A ∈ Rm×m is a constant matrix. Assume the matrix A has m has real eigenvalues,
λ1 ≤ λ1 ≤ · · · ≤ λm , and a complete set of right and left eigenvectors denoted by r k ∈ Rm
and lk ∈ Rm , respectively for k = 1, . . . , m. Furthermore, let X ∈ Rm×m denote the
matrix of right eigenvectors, X = [r 1 , . . . , r m ], and Λ ∈ Rm×m the diagonal matrix of
eigenvalues, Λ = diag(λ1 , . . . , λm ) so that A = XΛX −1 . The one-dimensional system (152) is
readily decoupled into scalar equations via the transformation to the characteristic variables
α ≡ X −1 u, solved as scalar equations, and recomposed into system form yielding
m
X
u(x, t) = lk · u0 (x − λk t) r k
k=1

Using this solution form, it is straightforward to solve the associated Riemann problem for
w(ξ, τ ; u, v) ∈ Rm
∂τ w + ∂ξ (Aw) = 0 in R × R+
with initial data 
u if ξ < 0
w(ξ, 0; u, v) =
v if ξ > 0
thereby producing the following Riemann-like numerical flux function
g(u, v) = Aw(0, τ > 0; u, v)
1 1
= (Au + Av) − |A| (v − u) (153)
2 2
with |A| ≡ X|Λ|X −1 . When this numerical flux function is used in the one-dimensional fully
discrete finite volume discretization with numerical solution unj at xj = j∆x and tn = n∆t
∆t
un+1
j = unj − (g(unj+1 , unj ) − (g(unj , unj−1 ))
∆x
∆t

= unj − A+ (unj − unj−1 ) + A− (unj+1 − unj ) A± ≡ XΛ± X −1 (154)
∆x
the resulting discretization reproduces the Courant–Isaacson–Rees (CIR) upwind discretization
for linear systems of hyperbolic equations given in Courant et al., 1952.

6.1.2. Numerical flux functions for nonlinear systems of conservation laws. In Godunov,
1959, exact solutions of the one-dimensional nonlinear Riemann problem of gas dynamics were
used in the construction of a numerical flux function
g R (u, v) = f (w(0, t > 0; u, v)) (155)
where w(ξ, τ ; u, v) ∈ Rm is a solution of a nonlinear Riemann problem
∂τ w + ∂ξ f (w) = 0 in R × R+ (156)

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 73

with initial data 

u if ξ < 0
w(ξ, 0; u, v) =
v if ξ > 0
Recall that solutions of the Riemann problem for gas dynamic systems with ideal gas law are
a composition of shock, contact, and rarefaction wave family solutions. For the gas dynamic
equations considered by Godunov, a unique solution of the Riemann problem exists for general
states u and v except those states producing a vacuum. Even so, the solution of the Riemann
problem is both mathematically and computationally nontrivial. Consequently, a number of
alternative numerical fluxes have been proposed that are more computationally efficient. These
alternative numerical fluxes can be sometimes interpreted as approximate Riemann solvers.
A partial list of alternative numerical fluxes is given here. A more detailed treatment of this
subject is given in Godlewski and Raviart (1991), Kröner (1997), and LeVeque (2002). In
∂f
describing these fluxes, it is assumed that the flux jacobian matrix A(u; n) ≡ ∂u · n is
m×m
diagonalizable via the matrix of right eigenvectors denoted by X ∈ R and the diagonal
matrix of real eigenvalues Λ ≡ diag[λ1 , . . . , λm ], so that A ≡ XΛX −1 .

• Osher–Solomon flux (Osher and Solomon, 1982). This numerical flux is a system
generalization of the Enquist–Osher flux (46) of Section 3. All wave families are
approximated in state space as rarefaction or inverted rarefaction waves with Lipschitz
continuous partial derivatives. The Osher-Solomon numerical flux is of the form
Z
1 1 v
g OS (u, v; n) = (f (u)+f (v)) · n− |A(w; n)| dw (157)
2 2 u
where |A| denotes the matrix absolute value via eigenvalues and eigenvectors. Integrating
on m rarefaction wave integral subpaths that are each parallel to a right eigenvector, a
system decoupling occurs on each subpath integration. Furthermore, for the gas dynamic
equations with ideal gas law, it is straightforward to construct m − 1 Riemann invariants
on each subpath thereby eliminating the need for path integration altogether. This
reduces the numerical flux calculation to purely algebraic computations with special
care taken at sonic points; see Osher and Solomon (1982).
• Roe flux (Roe, 1981). Roe’s numerical flux can be interpreted as approximating all wave
families as discontinuities. The numerical flux is of the form
1 1
g Roe (u, v; n) = (f (u) + f (v)) · n − |A(u, v; n)|(v − u)
2 2
where A(u, v; n) is the “Roe matrix” satisfying the matrix mean value identity
(f (v) − f (u)) · n = A(u, v; n)(v − u) (158)
with A(u, u; n) = A(u; n). For the equations of gas dynamics with ideal gas law, the
Roe matrix takes a particularly simple form. Steady discrete mesh-aligned shock profiles
are resolved with at most one intermediate point. The Roe flux does not preclude the
formation of entropy violating expansion shocks unless additional steps are taken near
sonic points.
• VFRoe flux (Gallouët and Masella, 1996). It is not always straightforward to find
matrices that satisfy the Roe condition (158). Another approach, using a linearized

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
74 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

problem but not requiring the Roe condition, is the so-called “VFRoe” scheme. The
idea is to use the Godunov flux (155) but with w as the solution of a linearized Riemann
problem rather than the exact Riemann problem (156). The VFRoe flux is given by
g VFRoe (u, v) = f (w(0, t > 0; u, v))
where w(ξ, τ ; u, v) ∈ Rm is a solution of a linearized Riemann problem
∂τ w + A(u, v)∂ξ w = 0 in R × R+ (159)
with initial data 
u if ξ < 0
w(ξ, 0; u, v) =
v if ξ > 0
where A(u, v) is an approximate jacobian matrix. Several choices of such matrices exist
as discussed in Masella et al., 1999; Brun et al., 2000; Buffard et al., 2000; Gallouët
et al., 2002; Gallouët et al., 2003. As in the case of the Roe scheme, solutions violating
an entropy inequality may occur. Entropy correction techniques have been designed to
avoid this problem (Helluy et al., 2010).
• Steger–Warming flux vector splitting (Steger and Warming, 1981). Steger and Warming
considered a splitting of the flux vector for the gas dynamic equations with ideal gas
law that exploited the fact that the flux vector is homogeneous of degree one in the
conserved variables. Euler’s identity for homogeneous functions of degree one then yields
f (u) · n = A(u; n)u. Steger and Warming then considered the matrix splitting
A = A+ + A− , A± ≡ XΛ± X −1
where Λ± is computed componentwise. From this matrix splitting, the final upwind
numerical flux function was constructed as
g SW (u, v; n) = A+ (u; n) u + A− (v; n)v
Although not part of their explicit construction, for the gas dynamic equations with ideal
gas law, the jacobian matrix ∂g SW /∂u has eigenvalues that are all nonnegative and the
jacobian matrix ∂g SW /∂v has eigenvalues that are all nonpositive whenever the ratio of
specific heats γ lies in the interval [1, 5/3]. The matrix splitting leads to numerical fluxes
that do not vary smoothly near sonic and stagnation points. Use of the Steger–Warming
flux splitting in the schemes of Sections 3 and 4 results in rather poor resolution of
linearly degenerate contact waves and velocity slip surfaces due to the introduction of
excessive artificial diffusion for these wave families.
• van Leer flux vector splitting. van Leer (1982) provided an alternative flux splitting for
the gas dynamic equations that produces a numerical flux of the form
g VL (u, v; n) = f − (u; n) + f + (v; n)
using special Mach number polynomials to construct fluxes that remain smooth near
sonic and stagnation points. As part of the splitting construction, the jacobian matrix
∂g SW /∂u has eigenvalues that are all nonnegative and the matrix ∂g SW /∂v has
eigenvalues that are all nonpositive. The resulting expressions for the flux splitting
are somewhat simpler when compared to the Steger–Warming splitting. The van Leer
splitting also introduces excessive diffusion in the resolution of linearly degenerate contact
waves and velocity slip surfaces.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 75

• System local Lax–Friedrichs flux. This numerical flux is the system equation counterpart
of the scalar local Lax–Friedrichs flux (39). For systems of conservation laws, the Lax–
Friedrichs flux is given by
1 1
g LF (u, v; n) = (f (u) + f (v)) · n − αmax (v − u)
2 2
where αmax is calculated from
αmax = max sup |λk (w; n)|
1≤k≤m w∈[u,v]

The system Lax-Friedrichs flux is usually not applied on the boundary of domains since it
generally requires an overspecification of boundary data. The system Lax–Friedrichs flux
introduces a relatively large amount of artificial diffusion when used in the schemes of
Section 3. Consequently, this numerical flux is typically only used together with relatively
high-order reconstruction schemes where the detrimental effects of excessive artificial
diffusion are mitigated.
• Harten–Lax–van Leer flux (Harten, Lax and van Leer, 1983). The Harten–Lax–van Leer
numerical flux originates from a simplified two wave model of more general m wave
systems such that waves associated with the smallest and largest characteristic speeds
of the m wave system are always accurately represented in the two-wave model. The
following numerical flux results from this simplified two-wave model
1
g HLL (u, v; n) = (f (u) + f (v)) · n
2
1 αmax + αmin αmax αmin
− (f (v) − f (u)) · n + (v − u)
2 αmax − αmin αmax − αmin
where
αmax = max (0, sup λk (w; n)), αmin = min (0, inf λk (w; n))
1≤k≤m w∈[u,v] 1≤k≤m w∈[u,v]

Using this flux, full upwinding is obtained for supersonic flow. Modifications of this flux
are suggested in Einfeldt et al. (1998) to improve the resolution of intermediate waves
as well.
• Entropy stable fluxes. An important class of numerical fluxes are motivated from energy
analysis of the FVM for symmetrizable systems of conservation laws. A brief introduction
to this theory and resulting flux functions are described in Section 6.1.3. Use of these
fluxes in the FVM guarantees a form of global entropy stability as characterized in
Theorem 27.
• Relaxation schemes (Jin and Xin, 1995). Another way to avoid solving the Riemann
problem (either because it is time consuming or because it not solvable) is to transform
the hyperbolic system into a simpler one, for instance, into a linear system with relaxation
term such as in Jin and Xin, 1995. The relaxation scheme can be reinterpreted as defining
a particular approximate Riemann solver for the original system conservation laws, see
LeVeque and Pelanti, 2001. Relaxation has also been used to replace a real equation of
state (EOS) by the perfect gas EOS in the Euler system in Coquel and Perthame, 1998.
Relaxation schemes are used in a wide area of applications such as pressureless gases in
Berthon et al., 2006, two-phase flows in Pelanti, 2011, discontinuous fluxes in Karlsen

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
76 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

et al., 2003, multicomponent flows in Dellacherie, 2003, Gallouët et al., 2010, Dorogan
et al., 2012, magnetohydrodynamics in Bouchut et al., 2009, or kidney physiology in
Perthame, 2015, to cite only a few.
• Well-balanced schemes (Greenberg and LeRoux, 1996). The notion of well-balanced
schemes was introduced by Greenberg and LeRoux, 1996 and Gosse, 1996 to deal with
source terms in hyperbolic systems in such a way that the steady state solutions are
correctly approximated. Well-balanced schemes have been used in several applications;
see Bouchut and Zeitlin, 2008 and Bouchut, 2004 for a thorough presentation of the
schemes and Berthon, 2016; Amadori and Gosse, 2016; Desprès, 2016 for some recent
work on the subject.

Further examples of numerical fluxes include the kinetic flux vector splitting due to Deshpande
(1986), the convective upwind and split pressure (CUSP) flux of Jameson (1993) and Tatsumi,
Martinelli and Jameson (1994), the advection upstream spitting method(AUSM) flux of Liou
and Steffen (1993), and the related internal energy-based solvers on staggered grids of Herbin
et al., 2013; Herbin et al., 2016b. Note that these latter staggered schemes have also been
implemented in implicit and semi-implicit time discretizations using staggered meshes in
Gastaldo et al., 2011 or collocated meshes in Herbin et al., 2016a in coping with low-Mach
number flows.

6.1.3. Entropy stability of the finite volume method for systems of hyperbolic conservation laws.
In the symmetrization theory for first-order conservation laws (see Godunov, 1961; Mock, 1980;
Tadmor, 1987; Harten, 1983a) one seeks a locally invertible mapping u(v) : Rm 7→ Rm applied
to (150a) so that when transformed
d
∂u ∂v X ∂f i ∂v
+ =0
∂v ∂t
|{z} i=1
∂v
|{z} ∂xi
SPD symmetric
∂f i
the matrix ∂u∂v is SPD and the matrices ∂v are symmetric. Mock, 1980 has proved that
the existence of an entropy equation (151) with convex entropy function U is sufficient
to guarantee that the system can be symmetrized via a change of variable with the new
variables v (commonly referred to as entropy variables) calculated from v T = ∂U
∂u . The entropy
variables also serve as nonlinear weights such that the weighted combination of conservation
law equations equals the entropy extension equation, that is,
v · (ut + ∇ · f ) = Ut + ∇ · F
with the right-hand side equal to zero for smooth solutions but not necessarily so for nonsmooth
solutions. Specifically, the inequality in the entropy equation (151) repeated here
Ut + ∇ · F ≤ 0
implies a form of stability and decay in a closed entropy system:

• The total mathematical entropy is nonincreasing (macroscopic Boltzmann H-theorem)

Z
d
U (u(x, t)) dx ≤ 0 (160)
dt Ω

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 77

• A two-sided bound on the total mathematical entropy for t0 ≤ t

Z Z Z
∗
U (u ) dx ≤ U (u(x, t)) dx ≤ U (u(x, t0 )) dx (161)
Ω Ω Ω

where u∗ is the minimum mathematical entropy state (maximum physical entropy state)
Z
1
u∗ = u(x, t0 ) dx
meas(Ω) Ω
that is invariant in time for a closed entropy system.
• System stability for t0 ≤ t
ku(x, t) − u∗ kL2 (Ω) ≤ (c−1
0 C0 )
1/2
ku(x, t0 )) − u∗ kL2 (Ω) (162)
where C0 and c0 are L2 -norm induced positive upper and lower bounds of the
∂v
transformation jacobian ∂u .

The last system stability statement is readily obtained from convexity of the entropy function
after first rewriting the two-sided total entropy bound using a Taylor series with integral
remainder formula
Z Z 1
∂2U
(1 − θ)(u(x, t) − u∗ ) · (u(θ)) (u(x, t) − u∗ ) dθ dx
Ω 0 ∂u2
Z Z 1
∂2U
≤ (1 − θ)(u(x, t0 ) − u∗ ) · (u0 (θ)) (u(x, t0 ) − u∗ ) dθ dx
Ω 0 ∂u2
where u(θ) ≡ u∗ + θ (u(x, t) − u∗ ) and u0 (θ) ≡ u∗ + θ (u(x, t0 ) − u∗ ). This shows how one can
go from a statement of boundedness for the scalar convex entropy function U to a statement
of boundedness for the vector of conservation variables u.

These stability results have motivated the construction of numerical discretizations that
discretely inherit or mimic these forms of entropy stability. The following theorem considers
semidiscrete and fully discrete finite volume discretizations and provides sufficient conditions
to be imposed on the numerical flux function so that discrete counterparts of the total entropy
bounds (160) and (161) are obtained.
Theorem 27. (Discrete entropy stability) Consider the hyperbolic conservation law system
(150) with convex entropy extension equation (151) that is symmetrizable via the locally
invertible change of variables u 7→ v and cell-averages uK ≡ u(v K ) for a problem domain
representing a closed entropy system.

1. The semidiscrete finite volume discretization

d X
uK |K| + g(uK , uL ; nK,L ) |σK,L | = 0, ∀K ∈ T (163)
dt
σK,L ⊂∂K

exhibits nonincreasing discrete total entropy

d X
U (uK )|K| ≤ 0 (164)
dt
K⊂T

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
78 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

2. The fully discrete finite volume discretization

X
(un+1
K − unK ) |K| + g(un+1 n+1
K , uL ; nK,L ) |σK,L | = 0, ∀K ∈ T , n = 0, 1, . . .
σK,L ⊂∂K
(165)
satisfies the two-sided discrete total entropy bound
X X X
U (u∗ )|K| ≤ U (unK )|K| ≤ U (u0K )|K| (166)
K⊂T K⊂T K⊂T

with X
1
u∗ = u0K |K| (167)
meas(Ω)
K⊂T

if the numerical flux function satisfies any of the following related sufficient conditions:

• The numerical flux function is of the following path integration form; see Barth, 1998
Z
1 1 1
g(u− , u+ ; n) = (f (u− ) + f (u+ )) · n − |A(u(v(θ)))|uv [v]+
− dθ (168)
2 2 0
where |A|uv is the matrix absolute value with respect to the matrix uv (see remark) and
v(θ) ≡ v − + θ [v]+
−.
• The numerical flux function is of the following form for Q ∈ Rm×m ; see Tadmor, 2004;
Fjordholm et al., 2012
g(u− , u+ ; n) = g ∗ (u− , u+ ; n) − Q(u− , u+ ; n)[v]+
− (169)
where g ∗ (u− , u+ ; n) satisfies
∗ ∗
[v]+ +
− · g (u− , u+ ; n) = [F · n]− (170)
∗
with F ∗ a primitive function that satisfies f T = ∂F
∂v and Q is any dissipation matrix
such that (see remarks)
[v]+ +
− · Q[v]− ≥ 0 (171)
• The numerical flux function satisfies the system generalization of Osher’s E-flux
condition (42) for scalar conservation laws, see Barth, 2006
[v]+
− · (g(u− , u+ ; n) − f (u(v(θ))) · n) ≤ 0, ∀θ ∈ [0, 1] (172)
with v(θ) ≡ v − + θ [v]+
−.

In these equations, [v]+

− ≡ v + − v − and the shorthand notation u± ≡ u(v ± ) have been used.

Remarks 1.(Theorem 27)

• Observe that the formulations (163) and (165) contain nonlinearity, not only from the
use of numerical flux functions g(·, ·; ·) but also from the use of entropy variables in
uK ≡ u(v K ).

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 79

• The sufficient conditions (168) and (172) are those originally obtained for the
discontinuous Galerkin finite element method that have been reduced to the present
FVM by choosing finite-dimensional trial and test spaces equal to piecewise constant
polynomials.
• The matrix absolute value |A|uv with respect to the matrix uv appearing in the path
integrated flux (168) can be efficiently computed as

|A|uv = X|Λ|X T (173)

where X is the matrix of scaled right eigenvectors that diagonalizes A such that

A = XΛX −1 , uv = XX T

using the constructive proof given in Barth, 1998.

• These entropy stability results have been extended to finite volume methods with higher
order reconstruction together with specific forms of the (169) dissipation matrix Q
discussed in Fjordholm et al., 2012 including the obvious choice |A|uv

Q = X|Λ|X T (174)
∂F ∗
• It is well known that the assumption of the existence of F ∗ satisfying f Ti = ∂vi
used in (170) is not always generically valid. For example, in ideal compressible
∂F ∗
magnetohydrodynamics, f T = ∂vi − ∂v ∂φ
B i where B is the magnetic induction field
∂φ
and ∂v are involution multipliers, so that sufficient conditions resulting from entropy
stability analysis are significantly different; see Barth, 2007.

6.2. The Marker-and-Cell scheme for fluid flows

The Marker-and-Cell (MAC) scheme, introduced in the mid-60s by Harlow and Welch, 1965, is
one of the most popular engineering methods for approximating the Navier-Stokes equations.
This is primarily due to the simplicity, efficiency and remarkable mathematical properties of
the method; see for example Patankar, 1980; Wesseling, 2001. The first error analysis seems
to be that of Porsching, 1978 in the case of the time-dependent Stokes equations on uniform
square grids. The mathematical analysis of the scheme was performed for the steady state
Stokes equations and Navier–Stokes equations in Nicolaı̈des, 1992; Nicolaı̈des and Wu, 1996
for uniform rectangular meshes with H 2 regularity assumption on the pressure. In the 90s,
using the tools that were developed for the finite volume theory that can be found in Eymard
et al., 2000, an order 1 error estimate for nonuniform meshes was obtained in Blanc, 1999, with
order 2 convergence for uniform meshes, under the usual regularity assumptions (H 2 for the
velocities, H 1 for the pressure). The convergence of the MAC scheme for the Stokes equations
with a right-hand side in H −1 (Ω) was later proved in Blanc, 2005.

6.2.1. The MAC scheme for the steady-state Navier–Stokes equations. Consider the
incompressible steady-state Navier-Stokes equations with Dirichlet boundary conditions on

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
80 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

a bounded domain Ω of Rd , d = 2 or 3. Denoting by ū the velocity and by p̄ the pressure, the

governing equations are given by

div ū = 0 in Ω (175a)
− ∆ū + (ū · ∇)ū + ∇p̄ = f in Ω (175b)
ū = 0 on ∂Ω (175c)

In this section, the notation (ū, p̄) denotes a solution to the continuous problem, while
the notation (u, p) denotes a solution to the approximate problem obtained from the MAC
discretization. Writing the balance form of these equations for a control volume K for the
incompressibility condition (175a) and a control volume D for the momentum equation (175a)
(both included in Ω) yields
Z
ū · n = 0 (176a)
∂K
Z Z Z
− ∇ū · n + ū(ū · n) + ∇p̄ = f (176b)
∂D ∂D D

2 d
where f ∈ L (Ω) . The balance form (176a) of the incompressibility condition is also given by
X Z
ū · nK,σ = 0
σ⊂∂K σ

where nK,σ denotes the unit vector normal to σ and outward K.

In the MAC scheme, a Cartesian mesh T of the domain is used with the set of edges (or
faces in 3D) denoted by E. The discrete velocity unknowns uσ are located on the edges (or
faces) σ of the mesh. These unknowns are approximations of ū · ei on the faces (or edges)
that are orthogonal to ei , the ith vector of the canonical basis of Rd . This is illustrated in
Figure 13 which depicts a 2-D MAC grid where the horizontal components of the velocities are
located at the × symbols on the vertical edges and the vertical components of the velocities
are located at the  symbols on the horizontal edges. Taking for K a pressure cell (dotted cell
on Figure 13) with edges R(or faces) orthogonal to the vectors of the canonical basis, a natural
discretization of the flux σ u · nK,σ is |σ|uσ K,σ where K,σ = nK,σ · ei = ±1 depending on
the direction of nK,σ and uσ is the discrete unknown associated with the face σ. If σ belongs
to the set E (i) of edges that are orthogonal to the ith unit vector of the canonical basis, then
uσ is an approximation of the mean value of u · ei over σ. The balance form (176a) of the
incompressibility condition may then be discretized as
X
|σ|uσ K,σ = 0
σ∈EK

where EK is the set of edges (or faces) of K. Since the discrete unknowns are approximations
of the components of ū on the corresponding edges (or faces), it is natural to introduce the
(i)
spaces HE , i = 1, . . . , d for the discrete velocity unknowns as the set of piecewise constant
functions on the velocity control volumes that are centered on the edges (or faces) σ ∈ E (i) of
the mesh. The resulting velocity control volumes are the North-West (i = 1) and North-East

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 81

   
× • × • × • × • ×
   

× • × • × • × • ×

   

× • × • × • × • ×

   
× • × • × • × • ×
   

Figure 13. MAC unknowns: •, pressure; ×, horizontal velocity; , vertical velocity. Control volumes:
pressure cell, dotted area; horizontal and vertical velocities, striped areas.

(i = 2) dashed cells in Figure 13. In order to take (partially) into account the homogeneous
Dirichlet boundary conditions, the spaces
n o
HE,0 = u ∈ HE , u(x) = 0 ∀x ∈ Dσ , σ ∈ Eeext , i = 1, . . . , d
(i) (i) (i)

Qd (i)
are defined. By introducing HE,0 = i=1 HE,0 , a vector function u ∈ HE,0 can be defined as
X
u = (u1 , . . . , ud ) with ui = uσ 11Dσ , where 11Dσ is the characteristic function of Dσ , that
σ∈E (i)
is 11Dσ (x) = 1 if x ∈ Dσ and zero otherwise. A discrete divergence divK u of u ∈ HE,0 can
then be defined on the cell K
1 X
divK u = |σ|uσ K,σ (177)
|K|
σ∈EK

The discrete divergence of u = (u1 , . . . , ud ) ∈ HE,0 may also be written as

d
X
divT (u) = (δi ui )K 11K (178)
i=1

where 11K is the characteristic function of the control volume K and (δi ui )K is the discrete
derivative of ui on K defined by
|σ| −→
(δi ui )K = (uσ0 − uσ ) with K = [σσ 0 ] and σ, σ 0 ∈ E (i) (179)
|K|
−→
where K = [σσ 0 ] means that σ and σ 0 are faces of K that are parallel and oriented. The
discrete derivatives and divergence are consistent in the sense that if ϕ = (ϕ1 , . . . , ϕd ) is a

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
82 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

smooth vector function over Ω and ΠE,i are “reasonable” interpolators (for i = 1, . . . , d), then
(i)
δi (ΠE ϕi ) tends to ∂i ϕ and divT (ΠE ϕ) tends to divϕ as the mesh size tends to 0. This is
stated precisely in the following lemma.
Lemma 10 (Discrete derivative and divergence consistency) Let D = (T , E) be an
MAC grid, and let ΠE be an interpolator from Cc∞ (Ω)d to HE,0 such that, for any ϕ =
(ϕ1 , . . . , ϕd )t ∈ (Cc∞ (Ω))d , there exists Cϕ ≥ 0 depending only on ϕ such that
 
(1) (d)
ΠE ϕ = Π(1)E ϕ 1 , · · · , Π
(d)
E ϕ d ∈ HE,0 × · · · × HE,0 , where
(i)
(180)
|ΠE ϕi (x) − ϕi (xσ )| ≤ Cϕ h2T , ∀x ∈ Dσ , ∀σ ∈ E (i) , ∀i = 1, · · · , d
Let η denote the regularity of the mesh
 
|σ| (i) 0 (j)
ηT = max , σ ∈ E , ∀σ ∈ E , i, j ∈ 1, . . . , d, i 6= j (181)
|σ 0 |
Then there exists Cϕ,η ≥ 0 such that
(i)
|δi ΠE ϕi (x) − ∂i ϕi (x)| ≤ Cϕ,η hT
for a.e. x ∈ Ω. As a consequence, if (Dn )n∈N = (Tn , En )n∈N is a sequence of MAC grids such
that ηn ≤ η for all n and hTn → 0 as n → +∞, then divTn (ΠEn ϕ) → divϕ uniformly as n →
+∞.

Examples of interpolators satisfying (180) are, for instance, the mean value over an edge (or
face) or the value of a function at the centroid of the edge (or face).

For the momentum equation, there is a need to distinguish the various components of the
velocity and momentum as well as associate a control volume with each unknown. To illustrate
this, Figure 13 shows a cell associated with the horizontal velocity (North-West dashed lines, ×
in the center) centered on a vertical edge belonging to E1 and a cell associated with the vertical
velocities (North-East dashed lines,  in the center) centered on a vertical edge belonging to
E2 . The set E of faces that are orthogonal to the ith unit vector ei can be decomposed as
(i) (i) (i) (i)
E (i) = Eint ∪ Eext where Eint (resp. Eext ) are the edges of E (i) that lie in the interior (resp. on
the boundary) of the domain Ω. Next, take for the control volume D in (176b) a velocity cell
Dσ , σ ∈ E (i) , corresponding to the ith velocity component. Each velocity grid consisting of such
cells is an admissible mesh for the Laplace operator. Therefore, the finite volume discretization
given in Section (5) can be directly applied. The diffusion fluxes as in (120) for each velocity
grid i = 1, . . . , d (replacing the cell K by a cell Dσ and the edge σ by an edge  of the set
e σ ) of edges (or faces) of the velocity cell Dσ are defined by
E(D
( ||
(d,i) − d (uσ0 − uσ ) if  = Dσ |Dσ0
Fσ, (ui ) = (182)
− ||
d (−uσ ) if  ⊂ ∂Ω
(c,i)
The numerical nonlinear convection flux Fσ, (u)(ui ) of the ith velocity component through an
0 0 (i)
R  = σ|σ separating the velocity cells Dσ and Dσ0 , σ, σ ∈ Eint , which is an approximation
edge
of  ui u · nσ, , can be defined by
(c,i) uσ + uσ0
Fσ, (u)(ui ) = ||uσ,
2
Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 83

where uσ, is an approximation of the velocity flux u · n through the edge . This flux
must be chosen carefully to obtain the L2 stability of the scheme. More precisely, a discrete
counterpart of the free divergence of u must be satisfied also on the dual cells. Two cases can
be distinguished:

- First case – the vector ei is normal to  and  is included in a primal cell K with
E (i) (K) = {σ, σ 0 }. The mass flux through  = σ|σ 0 is then given by
1
||uσ, = (−|σ|uK,σ + |σ 0 |uK,σ0 ) (183)
2

- Second case – the vector ei is tangent to  and  is the union of the halves of two primal
faces τ and τ 0 such that σ = K|L with τ ∈ E(K) and τ 0 ∈ E(L). The mass flux through
 is then given by
1
||uσ, = (|τ |uK,τ + |τ 0 |uL,τ 0 ) (184)
2

Using this definition, the usual finite volume property of local conservativity of the flux through
−−→
an interface σ|σ 0 is obtained, that is,

||uσ, = −||uσ0 , (185)

together with the following discrete free divergence condition on the dual cells:
X 1 X 1 X
||uσ, = |σ|uK,σ + |σ|uL,σ = 0 (186)
2 2
∈E(Dσ ) σ∈E(K) σ∈E(L)

Note that uσ, = 0 if  ⊂ ∂Ω, which is consistent with the boundary conditions (175c). A
discretization of each component of (176b) may now be written as
X X (i)
(d,i)
Fσ, (ui ) + CE (u)ui + |Dσ |δi p = |Dσ |fi,σ , i = 1, . . . , d
e σ)
∈E(D e σ)
∈E(D

where δi p denotes the discrete derivative of p given by the following (natural) differential
quotient for i = 1, . . . , d
1
(δi p)σ = (pL − pK ) for σ = K|L
dK,L
with K, L chosen such that nK,σ · ei = 1. Let LT denote the space P of functions that are
piecewise constant on the pressure control volumes and let p = P K∈T pK 11K . The discrete
gradient of p may be defined as ∇E p = (δ1 p, . . . , δd p)t where δi p = σ∈E (i) (δi p)σ 11Dσ . Using
int
this notation, the following discrete duality property can be succinctly written:
Z Z
q divT v + ∇E q · v = 0, ∀q ∈ LT ∀v ∈ HE,0 (187)
Ω Ω

In order to finish writing the scheme, let LT ,0 denote the set of functions that have zero mean
value and are piecewise constant on the pressure control volumes (i.e., the control volumes of

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
84 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

the rectangular mesh T , pointed cell in Figure 13). The ith component of the discrete Laplace
operator is classically defined as
(i) (i) (i)
−∆E : HE,0 −→ HE,0
X 1 X
ui 7−→ −∆E ui = − (d,i)
Fσ, (ui ) 11Dσ (188)
|Dσ |
σ∈E (i) ∈Ẽ(Dσ )

where Ẽ(Dσ ) denotes the faces of Dσ . Then the discrete Laplace operator of the full velocity
vector is defined by
−∆E : HE,0 −→ HE,0
(1) (d) (189)
v 7→ −∆E v = (−∆E u1 , . . . , −∆E ud )t
(i)
Finally, define the ith component CE (u) of the non linear convection operator by
(i) (i) (i)
CE (u) : HE,0 −→ HE,0
(i)
X 1 X
vi 7−→ CE (u)vi = F (c,i) (u)σ, (vi )11Dσ (190)
|Dσ |
σ∈Eeint
(i)
∈Ẽ(Dσ )
=σ|σ 0

and the full discrete convection operator C E (u), HE,0 −→ HE,0 by

(1) (d)
C E (u)v = (CE (u)v1 , . . . , CE (u)vd )t
With these notations, the MAC scheme for the discretization of (176) on a rectangular (non
uniform) staggered grid T can be written in following compact form for u ∈ HE,0 , p ∈ LT ,0 :
− ∆E u + C E (u)u + ∇E p = f E (191a)
divT u = 0 (191b)
(i)
where f E = (f1,E1 , . . . , fd,Ed ) and fi,Ei is the L2 projection on HE defined by
X  1 Z 
fi,Ei = fi (x) 11Dσ
(i)
|Dσ | Dσ
σ∈Eint

When performing the convergence analysis of the scheme (191), it can be shown that sequences
of approximate solutions tend to a weak solution of (176) as the mesh size tends to zero. A
weak formulation of (176) is given by
Find (ū, p̄) ∈ H01 (Ω)d × L20 (Ω) such that ∀(v, q) ∈ H01 (Ω)d × L20 (Ω)
Z Z Z Z
∇ū : ∇v + ((ū · ∇)ū) · v − p̄ divv = f ·v (192a)
ZΩ Ω Ω Ω

q divū = 0 (192b)
Ω

where L20 (Ω) stands for the subspace of L2 (Ω) of zero mean-valued functions. Similarly, a
weak form of the scheme equivalent to (191) is convenient for use in convergence analysis. As

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 85

in Section 5, a discrete H01 inner product is defined; however, here it must be defined on the
velocity grid. This inner product satisfies ∀(u, v) ∈ HE,0 2
Z Xd
−∆E u · v = [u, v]1,E,0 = [ui , vi ]1,E (i) ,0 (193)
Ω i=1

with
X || X ||
[ui , vi ]1,E (i) ,0 = (uσ − uσ0 ) (vσ − vσ0 ) + uσ vσ
d d
∈Eeint ∈Eeext
(i) (i)
−−→ ⊂∂(Dσ )
=σ|σ 0
(i) (i)
HE,0 × HE,0 → R HE,0 × HE,0 → R
The bilinear forms and are inner products on
(u, v) 7→ [ui , vi ]1,E (i) ,0 (u, v) 7→ [u, v]1,E,0
(i)
HE,0 and HE,0 respectively, which induce the following scalar and vector discrete H01 norms:
X || X || 2
kui k21,E (i) ,0 = [ui , ui ]1,E (i) ,0 = (uσ − uσ0 )2 + u for i = 1, . . . , d (194a)
d d σ
∈Eeint ∈Eeext
(i) (i)
−−→ ⊂∂(Dσ )
=σ|σ 0

and (194b)
d
X
kuk21,E,0 = [u, u]1,E,0 = kui k21,E (i) ,0 (194c)
i=1

M N
σ′
Dǫ
ǫ = σ|σ ′
σ = K|L

K L

Dσ

Figure 14. Full grid for definition of the derivative of the velocity.

When utilizing Cartesian grids, this inner product may be formulated as the L2 inner product
of discrete gradients. Consider the following discrete gradient of each velocity component ui :
X
∇E (i) ui = (δ1 ui , . . . , δd ui ) with δj ui = (δj ui )D 11D (195)
∈Ee(i)
⊥ej

uσ0 − uσ −−→
where (δj ui )D = with  = σ|σ 0 , and D =  × xσ xσ0 ; see Figure 14. This definition
d
is compatible with the definition of the discrete derivative (δi ui )K given by (179), since, if

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
86 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

 ⊂ K then D = K. With this definition, it follows that

Z
(i)
∇E (i) u · ∇E (i) v = [u, v]1,E (i) ,0 , ∀u, v ∈ HE,0 , ∀i = 1, . . . , d (196)
Ω

where [u, v]1,E (i) ,0 is the discrete H01 inner product defined by (193). Defining the following
gradient for u:
∇E u = (∇E (1) u1 , . . . , ∇E (d) ud )
it follows from the component formulas that
Z
∇E u : ∇E v = [u, v]1,E,0
Ω
With this formulation, the MAC scheme for the linear Stokes problem can be interpreted as
a gradient scheme in the sense introduced in Eymard et al., 2012 (see Eymard et al., 2015 for
more details on the generalization of this formulation to other schemes). Thanks to this result,
(strong) convergence of this discrete gradient to the gradient of the exact velocity as well as
strong convergence of the pressure can be shown.

Defining a weak form bE of the nonlinear convection operator

d
X (i)
bE (u, v, w) = bE (u, vi , wi ), ∀(u, v, w) ∈ HE,0 3 (197)
i=1
where for i = 1, . . . , d Z
(i) (i)
bE (u, vi , wi ) = CE (u)vi wi
Ω
it is now possible to introduce a weak formulation
Find (u, p) ∈ HE,0 × LT ,0 such that ∀(v, q) ∈ HE,0 × LT
Z Z Z
∇E u : ∇E v + bE (u, u, v) − p divT (v) = fE · v (198a)
ZΩ Ω Ω

divT u q = 0 (198b)
Ω
which is equivalent to the MAC scheme (191).

6.2.2. Convergence analysis of the MAC scheme. The proof of convergence of the MAC
scheme using this latter weak form then closely follows the proof of existence of a solution to
the Navier-Stokes equations; see for example Boyer and Fabrie, 2013. This analysis requires
estimates on the trilinear form bE .
Lemma 11 (Estimates on bE ) Let D = (T , E) be a MAC grid and let bE be defined by (197).
For d = 3, there exists CηT > 0, depending only on the regularity ηT of the mesh defined by
(181), such that
|bE (u, v, w)| ≤ CηT kukL4 (Ω)d kvk1,E,0 kwkL4 (Ω)d , ∀(u, v, w) ∈ E E × HE,0 2 (199)
and
|bE (u, v, w)| ≤ CηT kuk1,E,0 kvk1,E,0 kwk1,E,0 , ∀(u, v, w) ∈ E E × HE,0 2 (200)

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 87

An important property needed to obtain some a priori estimates on the velocity is that the
nonlinear convection term vanishes when taking u as test function in (192). This is also the
case for its discrete counterpart, as stated in the next lemma.
Lemma 12 (bE is skew-symmetric) Let (u, v, w) ∈ E E × HE,0 × HE,0 , then

bE (u, v, w) = −bE (u, w, v) (201)

and therefore
bE (u, u, u) = 0 ∀u ∈ E E (202)

From this property, using a so-called Fortin operator to obtain a discrete divergence free test
function, following estimates and existence result are obtained:
Theorem 28 (Existence and estimates) There exists a solution to (198) and CηT > 0,
depending only on the regularity ηT of the mesh and Ω, such that any solution of (198) satisfies
the stability estimate
kuk1,E,0 + kpkL2 (Ω) ≤ CηT kf kL2 (Ω)d (203)

A challenging part of the convergence analysis is the study of the nonlinear convection term.
One approach is to reconstruct a full grid velocity for each component, which converges as the
component itself.
Lemma 13 (Weak consistency of the nonlinear convection term) Let (Dn )n∈N with
Dn = (Tn , En ) be a sequence of meshes such that hTn = maxK∈Tn diam(K) → 0 as n → +∞.
Further, assume that there exists η > 0 such that ηTn ≤ η for any n ∈ N (with ηTn defined by
(181)). Let (v n )n∈N and (wn )n∈N be two sequences of functions such that

- v n ∈ HEn ,0 and wn ∈ HEn ,0 ,

- the sequences (v n )n∈N and (wn )n∈N converge in L2 (Ω)d to v̄ ∈ L2 (Ω)d and w̄ ∈ L2 (Ω)d
respectively.

Let (ΠEn )n∈N be a family of interpolators satisfying (180) and let ϕ ∈ Cc∞ (Ω)d . Then
bE (v n , wn , ΠEn ϕ) → b(v̄, w̄, ϕ) as n → +∞.

The above consistency result together with the estimates on the velocity and the pressure yield
the following convergence result:
Theorem 29 (Convergence of the scheme) Let (Dn )n∈N with Dn = (Tn , En ) be a
sequence of meshes such that hTn = maxK∈Tn diam(K) → 0 as n → +∞ . Further, assume
that there exists η > 0 such that ηTn ≤ η for any n ∈ N [with ηTn defined by (181)]. Let
(un , pn ) be a solution to the MAC scheme (191) or its weak form (198), for D = Dn , then
there exists ū ∈ H01 (Ω)d and p̄ ∈ L2 (Ω) such that, up to a subsequence

• the sequence (un )n∈N converges to ū in L2 (Ω)d ,

• the sequence (∇n un )n∈N converges to ∇ū in L2 (Ω)d×d ,

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
88 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

• the sequence (pn )n∈N converges to p̄ in L2 (Ω),

• (ū, p̄) is a solution to the weak formulation (192).

The proof of this result is obtained by taking interpolates of smooth functions as test functions
in the weak form of the scheme (198) and passing to the limit in the weak form of the scheme
Herbin et al., 2014; Gallouët et al., 2015a.

6.2.3. Further developments in the incompressible and compressible case. The time
dependent Navier–Stokes equations can also be easily discretized with the MAC scheme, either
with a time implicit scheme or else with a pressure correction scheme. In the case of the implicit
scheme, the scheme is proven to be convergent in Gallouët et al., 2015a. The case of the pressure
correction scheme remains an open question.

The MAC scheme can also be easily written and implemented for compressible flows. Some
theoretical results have been proven in the case of a perfect gas and the implicit Euler scheme
for the compressible Stokes equations (Eymard et al., 2010b) and in the case of the semi-
stationary Navier–Stokes equations (Gallouët et al., 2015c). Convergence results and error
estimates were also obtained for the same type of schemes on simplicial and quadrilateral
staggered grids using the Crouxeix–Raviart and the Rannacher–Turek finite element spaces
for the discretization of the diffusion operator; see Gallouët et al., 2015b,Gallouët et al., 2015d.

In the case of the Euler and Navier–Stokes equations, the convergence remains an open
question because of the lack of estimates. However, stability results exist. In particular, the
conservation of the discrete kinetic equation can be obtained (Gallouët et al., 2010; Herbin
and Latché, 2010). The weak consistency of the scheme has been proved for a decoupled
scheme in the case of the Euler (Herbin et al., 2013; Herbin et al., 2016b). Weak consistency
in this context means that if some estimates on the approximate solutions are assumed, then a
sequence of approximate solutions can be shown to converge to a weak solution of the system
as the mesh and time steps tend to zero (under appropriate CFL conditions). The limit of
the scheme may also be shown to be an entropy weak solution of the Euler system (for the
perfect gas EOS). Note that similar type schemes on simplicial and quadrilateral staggered
grids have been developed and studied using the Crouxeix–Raviart and the Rannacher–Turek
finite element spaces for the discretization of the diffusion operator. Both an implicit scheme
and a pressure correction scheme have been studied in Gallouët et al., 2008; Babik et al., 2011;
Gastaldo et al., 2011.

7. Related Chapters

(See also Finite Element Methods, Discontinuous Galerkin Methods for Computational Fluid
Dynamics, Aerodynamics)

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 89

REFERENCES

I. Aavatsmark, T. Barkve, O. Boe, and T. Mannseth. Discretization on unstructured grids

for inhomogeneous, anisotropic media. II. Discussion and numerical results. SIAM J. Sci.
Comput. 1998a; 19(5):1717–1736.
I. Aavatsmark, T. Barkve, O. Boe, and T. Mannseth. Discretization on unstructured grids for
inhomogeneous, anisotropic media. Part I: Derivation of the methods. SIAM Journal on
Sci. Comp. 1998b; 19:1700–1716.
I. Aavatsmark, E. Reiso, and E. Teigland. Control volume discretization for quadrilateral grids
with faults and local refinements. Comput. Geosci. 2001; 5:1–23.
I. Aavatsmark, G. T. Eigestad, B. T. Mallison, and J. M. Nordbotten. A compact
multipoint flux approximation method with improved robustness. Numer. Methods Partial
Differential Equations 2008; 24(5):1329–1360.
R. Abgrall. On essentially non-oscillatory schemes on unstructured meshes: analysis and
implementation. J. Comput. Phys. 1994; 114:45–58.
L. Agelas, D. A. Di Pietro, and J. Droniou. The G method for heterogeneous anisotropic
diffusion on general meshes. M2AN Math. Model. Numer. Anal. 2010; 44(4):597–625.
D. Amadori and L. Gosse. Trror Estimates for well-balanced and time-split schemes on a
locally damped wave equation. Math. Comp. 2016; 85(298):601–633.
B. Andreianov, M. Bendahmane, and K. Karlsen. A gradient reconstruction formula for finite-
volume schemes and discrete duality. In Finite volumes for complex applications V, ISTE,
London 2008, pp. 161–168.
B. Andreianov, F. Boyer, and F. Hubert. Discrete duality finite volume schemes for Leray-
Lions-type elliptic problems on general 2D meshes. Numer. Methods Partial Differential
Equations 2007; 23(1):145–195.
B. Andreianov, M. Bendahmane, and K. H. Karlsen. Discrete duality finite volume schemes
for doubly nonlinear degenerate hyperbolic-parabolic equations. J. Hyperbolic Differ. Eq.
2010; 7(1):1–67.
B. Andreianov, R. Eymard, M. Ghilani, and. N Marhraoui. Finite volume approximation of
degenerate two-phase flow model with unlimited air mobility. Numer. Methods Partial
Differential Equations 2013; 2(29):441–474.
B. Andreianov, M. Bendahmane, and M. Saad. Finite volume methods for degenerate
chemotaxis model. J. Comput. Appl. Math. 2011; 235(14):4015–4031.
B. Andreianov, M. Bendahmane, and R. Ruiz-Baier. Analysis of a finite volume method for a
cross-diffusion model in population dynamics. J. Comput. Appl. Math. 2011; 235(2):307–
344.
B. Andreianov, P. Goatin, and N. Seguin. Finite volume schemes for locally constrained
conservation laws. Numer. Math. 2010; 115(4):609–645.
T. Arbogast, D. Estep, B. Sheehan, and S. Tavener. A posteriori error estimates for mixed
finite element and finite volume methods for problems coupled through a boundary with
nonmatching grids. IMA J. Numer. Anal. 2014; 34(4):1625–1653.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
90 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

F. Archambeau, J.-M. Hérard, and J. Laviéville. Comparative study of pressure-correction

and Godunov-type schemes on unsteady compressible cases. Comput. & Fluids 2009;
38(8):1495–1509.
F. Babik, R. Herbin, W. Kheriji, and J.-C. Latché. Discretization of the viscous dissipation
term with the MAC scheme. In Finite volumes for complex applications. VI. Problems &
perspectives. Volume 1, 2, vol. 4 of Springer Proc. Math., Springer, Heidelberg 2011, pp.
571–579.
J. Ball. A version of the fundamental theorem for Young measures, PDEs and continuum
models of phase transitions (Nice, 1988). Lecture Notes in Phys. 1989; 344:207-215,
Springer, Berlin, 1989.
R. Bank and D. J. Rose. Some error estimates for the box method. SIAM J. Numer. Anal.
1987; 24:777–787.
T. J. Barth and P. Frederickson 1990. Higher order solution of the Euler equations on
unstructured grids using quadratic reconstruction. Report 90-0013, American Institute
for Aeronautics and Astronautics, 1990.
T. J. Barth and D. C. Jespersen 1989. The design and application of upwind schemes
on unstructured meshes. Report 89-0366, American Institute for Aeronautics and
Astronautics, 1989.
T. J. Barth. Numerical methods for gasdynamic systems on ynstructured meshes. In Kröner,
Ohlberger, & Rohde (eds.), An Introduction to Recent Developments in Theory and
Numerics for Conservation Laws, vol. 5 of Lecture Notes in Computational Science and
Engineering, pp. 195–285. Springer-Verlag, Heidelberg 1998.
T. J. Barth. On discontinuous Galerkin approximations of Boltzmann moment systems with
Levermore closure. Comput. Meth. Appl. Mech. Engrg. 2006; 195:3311–3330.
T. J. Barth. On the role of involutions in the discontinuous Galerkin discretization of Maxwell
and magnetohydrodynamics systems. IMA Vol. Math. Appl. 2007; Springer, 142:69–88.
P. Batten, C. Lambert, and D. M. Causon. Positively conservative high-resolution convection
schemes for unstructured elements. Int. J. Numer. Meth. Engrg. 1996; 39:1821–1838.
A. Berman and R. J. Plemmons. Nonnegative matrices in the mathematical sciences, volume 9
of Classics in Applied Mathematics. Society for Industrial and Applied Mathematics
(SIAM), Philadelphia, PA, 1994. Revised reprint of the 1979 original.
C. Berthon, Y. Coudière, and V. Desveaux. Second-order MUSCL schemes based on Dual Mesh
Gradient Reconstruction (DMGR). Mathematical Modelling and Numerical Analysis 2014;
48:583–602.
C. Berthon, M. Breuß, and M.-O. Titeux. A relaxation scheme for the approximation of
the pressureless Euler equations. Numer. Methods Partial Differential Equations 2006;
22(2):484–505.
C. Berthon. A fully well-balanced, positive and entropy-satisfying Godunov-type method for
the shallow-water equations. Math. Comp., 85(299):1281–1307, 2016.
M. Bessemoulin-Chatard, C. Chainais-Hillairet, and F. Filbet. On discrete functional
inequalities for some finite volume schemes. IMA J. Numer. Anal. 2015; 35(3):1125–
1149.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 91

M. Bessemoulin-Chatard, C. Chainais-Hillairet, and M.H. Vignal. Study of a finite volume

scheme for the drift-diffusion system. Asymptotic behavior in the quasi-neutral limit.
SIAM J. Numer. Anal. 2014; 52(4):1666–1691.
M. Bessemoulin-Chatard, C. Chainais-Hillairet, and F. Filbet. A finite volume scheme for
nonlinear degenerate parabolic equations. SIAM J. Sci. Comput. 2012; 34(5):B559–B583,
2012.
V. Billey, J. Périaux, P. Perrier, and B. Stoufflet. 2-D and 3-D Euler computations with finite
element methods in aerodynamics. vol. 1270 of Lecture Notes in Mathematics. Springer-
Verlag, Berlin, 1987.
P. Blanc. Error estimate for a finite volume scheme on a MAC mesh for the Stokes problem. In
Finite volumes for complex applications II, Hermes Sci. Publ., Paris, 1999, pp. 117–124.
P. Blanc. Convergence of a finite volume scheme on a MAC mesh for the Stokes problem with
right-hand side in H −1 . In Finite volumes for complex applications IV, ISTE, London,
2005, pp. 133–142.
L. Boccardo and T. Gallouët. Nonlinear elliptic and parabolic equations involving measure
data. J. Funct. Anal. 1989; 87(1):149–169.
J. Bonelle and A. Ern. Analysis of compatible discrete operator schemes for elliptic problems
on polyhedral meshes. ESAIM Math. Model. Numer. Anal. 2014; 48(2):553–581.
J. Bonelle, D. A. Di Pietro, and A. Ern. Low-order reconstruction operators on polyhedral
meshes: application to compatible discrete operator schemes. Comput. Aided Geom.
Design 2015; 35/36:27–41.
J. Boris and D. Book. Flux corrected transport: SHASTA, a fluid transport algorithm that
works. J. Comp. Phys. 1973; 11:38–69.
F. Bouchut & B. Perthame. Kruzkov’s estimates for scalar conservation laws revisited. Trans.
Am. Math. Soc. 1998; 350(7):2847–2870.
F. Bouchut. Nonlinear stability of finite volume methods for hyperbolic conservation laws and
well-balanced schemes for sources. Frontiers in Mathematics. Birkhäuser Verlag, Basel,
2004.
F. Bouchut, Christian Klingenberg, and Knut Waagan. An approximate Riemann solver for
ideal MHD based on relaxation. In Hyperbolic problems: theory, numerics and applications,
volume 67 of Proc. Sympos. Appl. Math., pages 439–443. Amer. Math. Soc., Providence,
RI, 2009.
F. Bouchut and Vladimir Zeitlin. A robust well-balanced scheme for multi-layer shallow water
equations. Discrete Contin. Dyn. Syst. Ser. B 2010; 13(4):739–758.
F. Bouchut and Vladimir Zeitlin. Finite volume schemes for the approximation via
characteristics of linear convection equations with irregular data. J. Evol. Equ. 2011;
113(4):687–724.
N. Bouillard, R. Eymard, R. Herbin, and Ph. Montarnal. Diffusion with dissolution and
precipitation in a porous medium: mathematical analysis and numerical approximation
of a simplified model. M2AN Math. Model. Numer. Anal. 2007; 41(6):975–1000.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
92 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

F. Boyer and F. Hubert. Finite volume method for 2D linear and nonlinear elliptic problems
with discontinuities. SIAM Journal on Numerical Analysis 2008; 46(6):3032–3070.
F. Boyer and P. Fabrie. Mathematical tools for the study of the incompressible Navier-Stokes
equations and related models, vol. 183 of Applied Mathematical Sciences. Springer, New
York, 2013.
A. Bradji. Some simple error estimates of finite volume approximate solution for parabolic
equations. Comptes Rendus Mathematique 2008; 346(9):571–574.
A. Bradji and R. Herbin. Discretization of coupled heat and electrical diffusion problems
by finite-element and finite-volume methods. IMA journal of numerical analysis 2008;
28(3):469–495.
A. Bradji and J. Fuhrmann. Error estimates of the discretization of linear parabolic
equations on general nonconforming spatial grids. Comptes Rendus Mathematique 2010;
348(19):1119–1122.
G. Brun, J.-M. Hérard, D. Jeandel, and M. Uhlmann. An approximate Roe-type Riemann
solver for a class of realizable second order closures. Int. J. Comput. Fluid Dyn. 2000;
13(3):233–249.
T. Buffard and S. Clain. Monoslope and multislope MUSCL methods for unstructured meshes.
Journal of Computational Physics 2010; 229:3745–3776.
T. Buffard, T. Gallouët, and J.-M. Hérard. A sequel to a rough Godunov scheme: application
to real gases. Comput. & Fluids 2000; 29(7):813–847.
Z. Q. Cai. On the finite volume element method. Numer. Math. 1991; 58(7):713–735.
Z. Q. Cai and S. McCormick. On the accuracy of the finite volume element method for diffusion
equations on composite grids. SIAM J. Numer. Anal. 1990; 27(3):636–655.
Z. Q. Cai, J. Mandel, and S. McCormick. The finite volume element method for diffusion
equations on general triangulations. SIAM J. Numer. Anal. 1991; 28(2):392–402.
C. Calgaro, E. Chane-Kane, E. Creusé, and T. Goudon. L∞ -stability of vertex-based MUSCL
finite volume schemes on unstructured grids: Simulation of incompressible flows with high
density ratios. Journal of Computational Physics 2010; 229:6027–6046.
C. Cancès, M. Cathala, and C. Le Potier. Monotone corrections for generic cell-centered
finite volume approximations of anisotropic diffusion equations. Numer. Math. 2013;
125(3):387–417.
C. Cancès, I. S. Pop, and M. Vohralı́k. An a posteriori error estimate for vertex-centered finite
volume discretizations of immiscible incompressible two-phase flow. Math. Comp. 2014;
83(285):153–188.
J. Carrillo. Entropy solutions for nonlinear degenerate problems. Arch. Ration. Mech. Anal.
1999; 147:269–361.
C. Chalons, and F. Coquel. Modified Suliciu relaxation system and exact resolution of isolated
shock waves. Math. Models Methods Appl. Sci. 2014; 24(5):937–971.
C. Chainais-Hillairet. Finite volume schemes for a nonlinear hyperbolic equation: convergence
towards the entropy solution and error estimates. M2AN Math. Model. Numer. Anal.
1999; 33:129–156.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 93

C. Chainais-Hillairet. Second-order finite-volume schemes for a non-linear hyperbolic equation:

error estimates. Math. Methods Appl. Sci. 2000; 23(5):467–490.
C. Chainais-Hillairet and Y.-J. Peng. Finite volume approximation for degenerate drift-
diffusion system in several space dimensions. SMath. Methods Appl. Sci. 2004; 14(3):461–
481.
C. Chainais-Hillairet and J. Droniou. Convergence analysis of a mixed finite volume scheme
for an elliptic-parabolic system modeling miscible fluid flows in porous media. SIAM J.
Numer. Anal. 2007; 45(5):2228–2258.
C. Chainais-Hillairet and F. Filbet. Asymptotic behaviour of a finite-volume scheme for the
transient drift-diffusion model. IMA J. Numer. Anal. 2007; 27(45):689–716.
C. Chainais-Hillairet and J. Droniou. Finite-volume schemes for noncoercive elliptic problems
with Neumann boundary conditions. IMA J. Numer. Anal. 2011; 31(1):61–85.
C. Chainais-Hillairet and M. Gisclon and A. Jüngel. A finite-volume scheme for the
multidimensional quantum drift-diffusion model for semiconductors. Numer. Meth.
Partial Diff. Eq. 2011; 27(6):1483–1510.
C. Chainais-Hillairet, S.. Krell and A. Mouton. Study of discrete duality finite volume schemes
for the Peaceman model. SIAM J. Sci. Comput. 2013; 35(6):A2928–A2952.
C. Chainais-Hillairet, S.. Krell and A. Mouton. Convergence analysis of a DDFV scheme
for a system describing miscible fluid flows in porous media. Numer. Methods Partial
Differential Equations 2015; 31(3):723–760.
S Champier and T Gallouët. Convergence d’un schéma décentré amont sur un maillage
triangulaire pour un problème hyperbolique linéaire. Modélisation mathématique et
analyse numérique 1992; 26(7):835–853.
S. Champier, T. Gallouët, and R. Herbin. Convergence of an upstream finite volume scheme for
a nonlinear hyperbolic equation on a triangular mesh. Numer. Math. 1993; 66(2):139–157.
P. Chatzipantelidis. A Finite volume method based on the Crouzeix-Raviart element for elliptic
problems. Numer. Math. 1999; 82:409–432.
P. Chatzipantelidis and R. D. Lazarov. Error estimates for a finite volume element method for
elliptic PDEs in nonconvex polygonal domains. SIAM J. Numer. Anal. 2005; 42(5):1932–
1958.
P. Chatzipantelidis, R. Lazarov, and V. Thomée. Some error estimates for the finite volume
element method for a parabolic problem. Comput. Methods Appl. Math. 2013; 13(3):251–
279.
L. Chen and M. Wang. Cell conservative flux recovery and a posteriori error estimate of
vertex-centered finite volume methods. Adv. Appl. Math. Mech. 2013; 5(5):705–727.
Q. Chen and M. Gunzburger. Goal-oriented a posteriori error estimation for finite volume
methods. J. Comput. Appl. Math. 2014; 265:69–82.
E. Chénier, R. Eymard, T. Gallouët and R. Herbin. An extension of the MAC scheme to
locally refined meshes : convergence analysis for the full tensor time-dependent Navier-
Stokes equations. Calcolo 2015; 52:69-107.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
94 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

S. Chou and Q. Li. Error estimates in L2 , H 1 and L∞ in covolume methods for elliptic and
parabolic problems: a unified approach. Math. Comp. 2000; 69:103–120.
P. G. Ciarlet. Basic error estimates for elliptic problems. In P. G. Ciarlet and J. L. Lions,
editors, Handbook of Numerical Analysis, Volume II, volume 2, pages 17–351. North
Holland, 1991.
S. Clain and V. Clauzon. L∞ stability of the MUSCL methods. Numerische Mathematik 2010;
116:31–64.
S. Clain. Finite volume maximum principle for hyperbolic scalar problems. SIAM J. Numer.
Anal. 2013; 51(1):467–490.
B. Cockburn, S. Y. Lin, and C. W. Shu. TVB Runge-Kutta local projection discontinuous
Galerkin finite element method for conservation laws. III. One-dimensional systems. J.
Comp. Phys. 1989; 84:90–113.
B. Cockburn and C. W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite
element method for conservation laws. II. General framework. J. Comp. Phys. 1989;
52:411–435.
B. Cockburn, S. Hou, and C. W. Shu. The Runge-Kutta local projection discontinuous Galerkin
finite element method for conservation laws. IV. The multidimensional case. Math. Comp.
1990; 54(190):545–581.
B. Cockburn, F. Coquel, and P. G. Lefloch. An error estimate for finite volume methods for
multidimensional conservation laws. Math. Comp. 1994; 63:77-103.
F. Coquel, J.-M. Hérard, K. Saleh, and N. Seguin. A robust entropy-satisfying finite volume
scheme for the isentropic Baer-Nunziato model. ESAIM Math. Model. Numer. Anal. 2014;
48(1):165–206.
F. Coquel, K. Saleh, and N. Seguin. A robust and entropy-satisfying numerical scheme for fluid
flows in discontinuous nozzles. Math. Models Methods Appl. Sci. 2014; 24(10):2043–2083.
B. Cockburn and H. Gau. A posteriori error estimates for general numerical methods for scalar
conservation laws. Comput. Appl. Math. 1995; 14:37–47.
B. Cockburn and P.-A. Gremaud. A priori error estimates for numerical methods for scalar
conservation laws. Part 1: The general approach. Math. Comput. 1996a; 65:533–573.
B. Cockburn and P.-A. Gremaud. A Priori error estimates for numerical methods for scalar
conservation laws. Part 2: Flux splitting montone schemes on irregular cartesian grids.
Math. Comput. 1996b; 66:547–572.
B. Cockburn and P.-A. Gremaud. A Priori error estimates for numerical methods for scalar
conservation laws. Part 3: Multidimensional flux-splitting monotone schemes on non-
cartesian grids. SIAM J. Numer. Anal. 1998a; 35:1775–1803.
B. Cockburn and C. W. Shu. The Runge-Kutta discontinuous Galerkin method for
conservation laws. V. Multidimensional systems. J. Comput. Phys. 1998b; 141(2):199–
224.
B. Cockburn and G. Gripenberg. Continuous dependence on the nonlinearities of solutions of
degenerate parabolic equations. J. Diff. Equations 1999; 151(2):231–251.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 95

B. Cockburn. Continuous dependence and error estimates for viscosity methods. Acta Numer.
2003; 12:127–180.
P. Colella and P. Woodward. The piecewise parabolic methods for gas-dynamical simulations.
J. Comp. Phys. 1984; 54:174–201.
F. Coquel and B. Perthame. Relaxation of energy and approximate Riemann solvers for general
pressure laws in fluid dynamics. SIAM J. Numer. Anal. 1998; 35(6):2223–2249.
Y. Coudière and F. Hubert. A 3D discrete duality finite volume method for nonlinear elliptic
equations. SIAM Journal on Scientific Computing 2011; 33(4):1739–1764.
Y. Coudière, F. Hubert, and G. Manzini. A CeVeFE DDFV scheme for discontinuous
anisotropic permeability tensors. In Finite volumes for complex applications VI, vol. 4 of
Springer Proc. Math., pp. 283–291. Springer, Heidelberg 2011.
Y. Coudière, C. Pierre, O. Rousseau, and R. Turpault. A 2D/3D discrete duality finite volume
scheme. Application to ECG simulation. Int. J. Finite Vol. 2009; 6(1).
Y. Coudière, J.-P. Vila, and P. Villedieu. Convergence rate of a finite volume scheme for a
two-dimensional convection-diffusion problem. M2AN Math. Model. Numer. Anal. 1999;
33(3):493–516.
R. Courant, K. Friedrichs, and H. Lewy. Über die partiellen Differenzengleichungen der
mathematischen Physik. Math. Ann. 1928; 100(1):32–74.
R. Courant, E. Isaacson, and M. Rees. On the solution of nonlinear hyperbolic differential
equations by finite differences. Comm. Pure Appl. Math. 1952; 5:243–255.
R. Courant, K. Friedrichs, and H. Lewy. On the partial difference equations of mathematical
physics. IBM J. Res. Develop. 1967; 11:215–234.
P.-H. Cournède, C. Debiez, and A. Dervieux. A positive MUSCL scheme for triangulations.
Report 3465, Institut National De Recherche En Informatique Et En Automatique
(INRIA), 1998.
M. Crandall and A. Majda. Monotone difference approximations of scalar conservation laws.
Math. Comp. 1980; 34:1–21.
J.-P. Croisille and I. Greff. An efficient box-scheme for convection-diffusion equations with
sharp contrast in the diffusion coefficients. Comput. & Fluids 2005; 34(4-5):461–489.
A. Dedner, C. Makridakis, and M. Ohlberger. Error control for a class of Runge-Kutta
discontinuous Galerkin methods for nonlinear conservation laws. SIAM J. Numer. Anal.
2007; 45(2):514–538.
M. Delanaye. Polynomial reconstruction finite volume schemes for the compressible Euler and
Navier-Stokes equations on unstructured adaptive grids. Ph.D. thesis, University of Liége,
Belgium.
S. Dellacherie. Relaxation schemes for the multicomponent Euler system. M2AN Math. Model.
Numer. Anal. 2003; 37(6):909–936.
S. M. Deshpande. On the Maxwellian distribution, symmetric form, and entropy conservation
for the Euler equations. Report TP-2613, NASA Langley, 1986.
J. A. Desideri and A. Dervieux. Compressible flow solvers using unstructured grids. von
Karman Institute Lecture Series 1988-05, von Karman Institute, Brussels, 1988.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
96 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

B. Desprès. An explicit a priori estimate for a finite volume approximation of linear advection
on non Cartesian grids. SIAM J. Numer. Anal. 2004; 42(2):484–504.
B. Desprès and C. Buet. The structure of well-balanced schemes for Friedrichs systems with
linear relaxation. Appl. Math. Comput. 2016; 272(2):440–459.
B. Desprès and F. Lagoutière. Generalized Harten formalism and longitudinal variation
diminishing schemes for linear advection on arbitrary grids. ESAIM: Mathematical
Modelling and Numerical Analysis 2010; 35(62):1149–1183.
D. A. Di Pietro, M. Vohralı́k, and C. Widmer. An a posteriori error estimator for a finite
volume discretization of the two-phase flow. In Finite volumes for complex applications.
VI. Problems & perspectives. Volume 1, 2, vol. 4 of Springer Proc. Math., pp. 341–349.
Springer, Heidelberg 2011.
D. A. Di Pietro and A. Ern. Mathematical aspects of discontinuous Galerkin methods, 69,
Springer, 2012.
D. A. Di Pietro and M. Vohralı́k. A review of recent advances in discretization methods, a
posteriori error analysis, and adaptive algorithms for numerical modeling in geosciences.
Oil and Gas Science and Technology 2014; 69(4):701–730.
R. J. DiPerna. Measure-valued solutions to conservation laws. Arch. Rational Mech. Anal.
1985; 88(3):223–270.
K. Domelevo and P. Omnes. A finite volume method for the Laplace equation on almost
arbitrary two-dimensional grids. M2AN Math. Model. Numer. Anal. 2005; 39(6):1203–
1249.
K. Dorogan and J.-Marc. Hérard. A two-dimensional relaxation scheme for the hybrid
modelling of gas-particle two-phase flows. Int. J. Finite Vol., 8:30, 2012.
O. Drblı́ková, A. Handlovičová, and K. Mikula. Error estimates of the finite volume scheme for
the nonlinear tensor-driven anisotropic diffusion. Appl. Numer. Math. 2009; 59(10):2548–
2570.
O. Drblı́ková and K. Mikula. Convergence analysis of finite volume scheme for nonlinear tensor
anisotropic diffusion in image processing. SIAM J. on Numer. Anal. 2007; 46(1):37–60.
J. Droniou. Non-coercive linear elliptic problems. Potential Analysis 2002; 17(2):181–203.
J. Droniou. Error estimates for the convergence of a finite volume discretization of convection-
diffusion equations. J. Numer. Math. 2003; 11(1):1–32.
J. Droniou. Finite volume schemes for diffusion equations: introduction to and review of
modern methods. Math. Models Methods Appl. Sci. 2014; 24(8):1575–1619.
J. Droniou and R. Eymard. A mixed finite volume scheme for anisotropic diffusion problems
on any grid. Numer. Math. 2006; 105(1):35–71.
J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. Gradient Schemes for
elliptic and parabolic problems. in preparation.
J. Droniou, R. Eymard, T. Gallouët, and R. Herbin. A unified approach to mimetic finite
difference, hybrid finite volume and mixed finite volume methods. Math. Models Methods
Appl. Sci. 2010; 20(2):265–295.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 97

J. Droniou, R. Eymard, T. Gallouët, and R. Herbin. Gradient schemes: a generic framework

for the discretisation of linear, nonlinear and nonlocal elliptic and parabolic equations.
Math. Models Methods Appl. Sci. 2013; 23(13):2395–2432.
J. Droniou, R. Eymard, and R. Herbin. Gradient schemes: generic tools for the numerical
analysis of diffusion equations. M2AN Math. Model. Numer. Anal., to appear.
J. Droniou and T. Gallouët. Finite volume methods for convection-diffusion equations with
right-hand side in H −1 . M2AN Math. Model. Numer. Anal. 2002; 36(4):705–724.
J. Droniou, T. Gallouët, and R. Herbin. A finite volume scheme for a noncoercive elliptic
equation with measure data. SIAM J. Numer. Anal. 2003; 41(6):1997–2031.
J. Droniou and C. Le Potier. Construction and convergence study of schemes preserving the
elliptic local maximum principle. SIAM J. Numer. Anal. 2011; 49(2):459–490.
J. Droniou and N. Nataraj in preparation. Improved L2 estimate for gradient schemes, and
super-convergence of finite volume methods. In preparation.
M. G. Edwards and C. F. Rogers 1994. A flux continuous scheme for the full tensor pressure
equation. In 4th European Conference on the Mathematics of Oil Recovery 1994.
M. G. Edwards and C. F. Rogers. Finite volume discretization with imposed flux continuity
for the general tensor pressure equation. Computational Geosciences 1998; 2(4):259–290.
B. Einfeldt, C. Munz, P. L. Roe, and B. Sjögreen. On Godunov-type methods near low
densities. J. Comp. Phys. 1992; 92:273–295.
C. Erath. A nonconforming a posteriori estimator for the coupling of cell-centered finite volume
and boundary element methods. Numer. Math. 2015; 131(3):425–451.
A. Ern and M. Vohralı́k. A unified framework for a posteriori error estimation in elliptic
and parabolic problems with application to finite volumes. In Finite volumes for complex
applications. VI. Problems & perspectives. Volume 1, 2, vol. 4 of Springer Proc. Math.,
pp. 821–837. Springer, Heidelberg 2011.
L. C. Evans. Partial differential equations, vol. 19 of Graduate Studies in Mathematics.
American Mathematical Society, Providence, RI, second ed. 2010.
S. Evje and K. H. Karlsen. An error estimate for viscous approximate solutions of degenerate
parabolic equations. J. Nonlin. Math. Phys. 2002; 9(3):262–281.
R. E. Ewing, T. Lin, and Y. Lin. On the accuracy of the finite volume element method based
on piecwise linear polynomials. SIAM J. Numer. Anal. 2002; 39(6):1865–1888.
R. Eymard, T. Gallouët, J. Vovelle. Limit boundary conditions for finite volume
approximations of some physical problems. J. Comput. Appl. Math. 2003; 161(2):349–369.
R. Eymard, T. Gallouët, M. Ghilani, and R. Herbin. Error estimates for the approximate
solution of a nonlinear hyperbolic equation given by finite volume schemes. IMA J. of
Numer. Anal. 1998; 18:563–594.
R. Eymard, T. Gallouët, R. Herbin, and A. Michel. Convergence of a finite volume scheme for
nonlinear degenerate parabolic equations. Numer. Math. 2002; 92(1):41–82.
R. Eymard, T. Gallouët and R. Herbin. Finite volume methods. In P. G. Ciarlet & J.-L. Lions
(eds.), Techniques of Scientific Computing, Part III, Handbook of Numerical Analysis,
VII, North-Holland, Amsterdam 2000, pp. 713–1020.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
98 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

R. Eymard, T. Gallouët and R. Herbin. Existence and uniqueness of the entropy solution to
a nonlinear hyperbolic equation. In Chinese Ann. Math. Ser. B,16(1):1–14.
R. Eymard, T. Gallouët, and R. Herbin. Finite volume approximation of elliptic problems and
convergence of an approximate gradient. Appl. Numer. Math. 2001; 37(1–2):31–53.
R. Eymard, T. Gallouët, and R. Herbin. Error estimates for approximate solutions of a
nonlinear convection-diffusion problem. Adv. Diff. Equations 2002; 7(4):419–440.
R. Eymard, T. Gallouët, and R. Herbin. Convergence of finite volume schemes for semilinear
convection diffusion equations. Numer. Math., 82(1):91–116, 1999.
R. Eymard, T. Gallouët, and R. Herbin. Discretization of heterogeneous and anisotropic
diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilization
and hybrid interfaces. IMA J. Numer. Anal. 2010a; 30(4):1009–1043.
R. Eymard, T. Gallouët, R. Herbin, and J.-C. Latché. Convergence of the MAC scheme for
the compressible Stokes equations. SIAM J. Numer. Anal. 2010b; 48(6):2218–2246.
R. Eymard, G. Henry, R. Herbin, F. Hubert, R. Klofkorn, and G. Manzini 2011. 3D Benchmark
on Discretization Schemes for Anisotropic Diffusion Problems on General Grids. In
Proceedings of Finite Volumes for Complex Applications VI, Praha. Springer, Springer
2011, pp. 895–930.
R. Eymard, C. Guichard, and R. Herbin. Small-stencil 3D schemes for diffusive flows in porous
media. M2AN Math. Model. Numer. Anal. 2012; 46:265–290.
R. Eymard, P. Féron, and C. Guichard 2015. Gradient schemes for the incompressible steady
Navier-Stokes problem. In 6th International conference on Approximation Methods and
Numerical Modelling in Environment and Natural Resources. Université de Pau 2015.
R. Eymard, T. Gallouët, and R. Herbin 2016. Finite volume methods: schemes and analysis.
in preparation.
I. Faille. A control volume method to solve an elliptic equation on a two-dimensional irregular
mesh. Comput. Methods Appl. Mech. Engrg. 1992a; 100(2):275–290.
I. Faille 1992b. Modélisation bidimensionnelle de la genèse et de la migration des hydrocarbures
dans un bassin sédimentaire. Ph.D. thesis, University Joseph Fourier, Grenoble I, 1992b.
M. Feistauer, J. Felcman, M. Lukáčová-Medvid’ová, and G. Warnecke. Error estimates for a
combined finite volume-finite element method for nonlinear convection-diffusion problems.
SIAM J. Numer. Anal. 1999; 36(5):1528–1548.
A. Fettah and T. Gallouët. Numerical approximation of the general compressible Stokes
problem. IMA J. Numer. Anal. 2013; 33(3):922–951.
F. Filbet. Convergence of a finite volume scheme for the Vlasov-Poisson system. SIAM J.
Numer. Anal. 2001; 39(4):1146–1169.
F. Filbet. A finite volume scheme for the Patlak-Keller-Segel chemotaxis model. Numer. Math.
2006; 104(4):457–488.
U. Fjordholm, S. Mishra, and E. Tadmor. Arbitrarily high-order accurate entropy stable
essentially nonoscillatory schemes for systems of conservation laws. SIAM J. Num. Anal.
2012; 50:544–573.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 99

P. Frolkovič and K. Mikula. Flux-based level set method: a finite volume method for evolving
interfaces. Applied numerical mathematics 2007; 57(4):436–454
T. Gallouët 2007. Nonlinear methods for linear equations. In T. Aliziane, K. Lemrabet,
A. Mokrane, & D. Teniou (eds.), Actes du 3eme colloque sur les Tendandes des
Applications Mathématiques en Tunise, Algérie, Maroc (14-18 avril 2007), AMNEDP-
USTHB 2007, pp. 17–22.
T. Gallouët, Ph. Helluy, J.-M. Hérard, and Julien Nussbaum. Hyperbolic relaxation models
for granular flows. M2AN Math. Model. Numer. Anal. 2010; 44(2):371–400.
T. Gallouët, L. Gastaldo, R. Herbin, and J.-C. Latché. An unconditionally stable pressure
correction scheme for the compressible barotropic Navier-Stokes equations. M2AN Math.
Model. Numer. Anal. 2008; 42(2):303–331.
T. Gallouët, J.-M. Hérard, and N. Seguin. Some approximate Godunov schemes to compute
shallow-water equations with topography. Comput. & Fluids, 32(4):479–513, 2003.
T. Gallouët, J.-M. Hérard, and N. Seguin. Some recent finite volume schemes to compute Euler
equations using real gas EOS. Internat. J. Numer. Methods Fluids 2002; 39(12):1073–
1138.
T. Gallouët and R. Herbin. Finite volume approximation of elliptic problems with irregular
data. In Finite volumes for complex applications II, Hermes Sci. Publ., Paris 1999, pp.
155–162.
T. Gallouët, R. Herbin, J.-C. Latché, and K. Mallem. Convergence of the Marker-and-Cell
scheme for the incompressible Navier-Stokes equations. submitted.
T. Gallouët, R. Herbin, J.-C. Latché, and D. Maltese. Convergence of the Marker-and-Cell
scheme for the compressible stationary Navier-Stokes equations. submitted.
T. Gallouët, R. Herbin, and J.-C. Latché. Kinetic energy control in explicit finite volume
discretizations of the incompressible and compressible Navier-Stokes equations. Int. J.
Finite Vol. 2010; 7(2):6.
T. Gallouët, R. Herbin, and M. H. Vignal. Error estimates on the approximate finite volume
solution of convection diffusion equations with general boundary conditions. SIAM J.
Numer. Anal. 2000; 37(6):1935–1972.
T. Gallouët, R. Herbin, D. Maltese, and A. Novotny. Convergence of the Marker-and-Cell
scheme for a semi-stationary compressible Navier-Stokes problem. under revision.
T. Gallouët, R. Herbin, D. Maltese, and A. Novotny. Error estimates for a numerical
approximation to the compressible barotropic NavierStokes equations. IMA Journal of
Numerical Analysis, 2015.
T. Gallouët, A. Larcher, and J. Latché. Convergence of a finite volume scheme for
the convection-diffusion equation with L1 data. Mathematics of Computation 2012;
81(279):1429–1454.
T. Gallouët and J.-M. Masella. Un schéma de Godunov approché. C. R. Acad. Sci. Paris Sér.
I Math. 1996; 323(1):77–84.
L. Gastaldo, R. Herbin, W. Kheriji, C. Lapuerta, and J.-C. Latché. Staggered discretizations,
pressure correction schemes and all speed barotropic flows. In Finite volumes for complex

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
100 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

applications. VI. Problems & perspectives. Volume 1, 2, vol. 4 of Springer Proc. Math.,
Springer, Heidelberg 2011, pp. 839–855.
A. Genty and C. Le Potier. Maximum and minimum principles for radionuclide transport
calculations in geological radioactive waste repository: comparison between a mixed hybrid
finite element method and finite volume element discretizations. Transp. Porous Media
2011; 88(1):65–85.
E. Godlewski and P.-A. Raviart. Hyperbolic Systems of Conservation Laws. Mathematiques
& Applications, Ellipses, Paris, 1991.
S. K. Godunov. A finite difference method for the numerical computation of discontinuous
solutions of the equations of fluid dynamics. Mat. Sb. 1959; 47:271–290.
S. K. Godunov. An interesting class of quasilinear systems. Dokl. Akad. Nauk. SSSR 1961;
139:521–523.
J. D. Goodman and R. J. L. Veque. On the accuracy of stable schemes for 2D conservation
laws. Math. Comp. 1985; 45(171):15–21.
S. Gottlieb and C.-W. Shu. Total variation diminishing Runge-Kutta schemes. Math. Comput.
1998; 67(221):73–85.
L. Gosse. A well-balanced scheme able to cope with hydrodynamic limits for linear kinetic
models. Appl. Math. Lett. 2015; 42(3):15–21.
L. Gosse. Un schéma-équilibre adapté aux lois de conservation scalaires non-homogènes. C.
R. Acad. Sci. Paris Sér. I Math. 1996; 323(5):543–546.
S. Gottlieb, C. W. Shu, and E. Tadmor. Strong stability-preserving high-order time
discretization methods. SIAM Rev. 2001; 43(1):89–112.
I. Greff. Nonconforming box-schemes for elliptic problems on rectangular grids. SIAM J.
Numer. Anal. 2007; 45(3):946–968.
J. M. Greenberg and A. Y. Leroux. A well-balanced scheme for the numerical processing of
source terms in hyperbolic equations. SIAM J. Numer. Anal. 1996; 33(1):1–16.
F. Harlow and J. Welch. Numerical calculation of time-dependent viscous incompressible flow
of fluid with a free surface. Physics of Fluids 1965; 8:2182–2189.
A. Harten. High resolution schemes for hyperbolic conservation laws. J. Comp. Phys. 1983a;
49:151–164.
A. Harten. On the symmetric form of systems of conservation laws with entropy. J. Comp.
Phys. 1983b; 49:357–393.
A. Harten. ENO schemes with subcell resolution. J. Comp. Phys. 1989; 83:148–184.
A. Harten, S. Osher, B. Engquist, and S. Chakravarthy. Some results on uniformly high order
accurate essentially non-oscillatory schemes. Appl. Num. Math. 1986; 2:347–377.
A. Harten, S. Osher, B. Engquist, and S. Chakravarthy. Uniformly high-order accurate
essentially nonoscillatory schemes III. J. Comput. Phys. 1987; 71(2):231–303.
A. Harten and S. Chakravarthy. Multi-dimensional ENO schemes for general geometries.
Report ICASE-91-76, Institute for Computer Applications in Science and Engineering
(ICASE), 1991.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 101

A. Harten, J. M. Hyman, and P. D. Lax. On finite-difference approximations and entropy

conditions for shocks. Comm. Pure and Appl. Math. 1976; 29:297–322.
A. Harten, P. D. Lax, and B. van Leer. On upstream differencing and Godunov-type schemes
for hyperbolic conservation laws. SIAM Rev. 1983; 25:35–61.
P. Helluy, J.-M. Hérard, H. Mathis, and S. Müller. A simple parameter-free entropy correction
for approximate Riemann solvers. C. R. Méc. Acad. Sci. Paris 2010; 338(9):493–498.
R. Herbin. An error estimate for a finite volume scheme for a diffusion-convection problem on
a triangular mesh. Numer. Methods Partial Differential Equations 1995; 11(2):165–173.
R. Herbin and M. Ohlberger. A posteriori error estimate for finite volume approximations
of convection diffusion problems. proceedings: Finite volumes for complex applications–
problems and perspectives, Porquerolles, Hermes Science Publications: Paris, 2002;
753-760.
R. Herbin and F. Hubert 2008. Benchmark on Discretization Schemes for Anisotropic
Diffusion Problems on General Grids for anisotropic heterogeneous diffusion problems.
In R. Eymard & J.-M. Hérard (eds.), Finite Volumes for Complex Applications V, Wiley,
2008, pp. 659–692.
R. Herbin and J.-C. Latché. Kinetic energy control in the MAC discretization of the
compressible Navier-Stokes equations. Int. J. Finite Vol. 2010; 7(2):6.
R. Herbin, J.-C. Latché, and T. T. Nguyen. Explicit staggered schemes for the compressible
Euler equations. In Applied mathematics in Savoie—AMIS 2012: Multiphase flow in
industrial and environmental engineering, vol. 40 of ESAIM Proc., EDP Sci., Les Ulis
2013, pp. 83–102.
R. Herbin, J.-C. Latché, and C. Zaza. A cell-centered pressure-correction scheme for the
compressible Euler equations. M3AS, under revision.
R. Herbin, J.-C. Latché and T. Nguyen. On some consistent explicit staggered schemes for
the shallow water and Euler equations. ESAIM: Mathematical Modelling and Numerical
Analysis, under revision.
R. Herbin, J.-C. Latché, and K. Mallem. Convergence of the MAC Scheme for the steady-
state incompressible Navier-Stokes equations on non-uniform grids. In J. Fuhrmann,
M. Ohlberger, & C. Rohde (eds.), Finite Volumes for Complex Applications VII-Methods
and Theoretical Aspects, vol. 77 of Springer Proceedings in Mathematics & Statistics,
Springer International Publishing 2014, pp. 343–351.
F. Hermeline. A finite volume method for the approximation of diffusion operators on distorted
meshes. J. Comput. Phys. 2000; 160(2):481–499.
F. Hermeline. Approximation of diffusion operators with discontinuous tensor coefficients
on distorted meshes. Computer methods in applied mechanics and engineering 2003;
192(16):1939–1959.
F. Hermeline. Approximation of 2-D and 3-D diffusion operators with variable full tensor
coefficients on arbitrary meshes. Comput. Meth. Appl. Mech. Engin. 2007; 196(21):2497–
2526.
F. Hermeline. A finite volume method for approximating 3D diffusion operators on general
meshes. J. Comp. Phys. 2009; 228(16):5763–5786.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
102 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

A. Jameson and P. D. Lax. Conditions for the construction of multipoint variation diminishing
difference schemes. Appl. Numer. Math. 1986; 2(3-5):235–345.
A. Jameson and P. D. Lax. Corrigendum: conditions for the construction of multipoint
variation diminishing difference schemes. Appl. Numer. Math. 1987; 3(3):289.
A. Jameson. Artificial diffusion, upwind biasing, limiters and their effect on accuracy and
convergence in transonic and hypersonic flows. Report AIAA-93-3359, American Institute
for Aeronautics and Astronautics, 1993.
G. Jiang and C. W. Shu. Efficient implementation of weighted ENO schemes. J. Comp. Phys.
1996; 126:202–228.
Shi Jin and Zhou Ping Xin. The relaxation schemes for systems of conservation laws in
arbitrary space dimensions. Comm. Pure Appl. Math. 1995; 48(3):235–276.
M. Küther. Error estimates for second order finite volume schemes using a TVD-Runge-
Kutta time discretization for a nonlinear scalar hyperbolic conservation law. East-West
J. Numer. Math. 2000; 8(4):299–322.
K. H. Karlsen and N. H. Risebro. On the uniqueness and stability of entropy solutions
of nonlinear degenerate parabolic equations with rough coefficients. Preprint 143,
Department of Mathematics, University of Bergen, 2000.
K. H. Karlsen, C. Klingenberg, and N. H. Risebro. Relaxation schemes for conservation laws
with discontinuous coefficients. In Hyperbolic problems: theory, numerics, applications,
Springer, Berlin, 2003, pp. 611–620.
B. Koren. Upwind schemes for the Navier-Stokes equations. In Proceedings of the Second
International Conference on Hyperbolic Problems, Vieweg, Braunschweig, 1988.
D. Kröner. Numerical Schemes for Conservation Laws. Wiley–Teubner, Stuttgart, 1997.
D. Kröner and M. Ohlberger. A posteriori error estimates for upwind finite volume schemes
for nonlinear conservation laws in multidimensions. Math. Comput. 2000; 69:25–39.
D. Kröner, S. Noelle, and M. Rokyta. Convergence of higher order upwind finite volume
schemes on unstructured grids for conservation laws in several space dimensions. Numer.
Math. 1995; 71:527–560.
S. N. Kruzkov. First order quasilinear equations in several independent variables. Math. USSR
Sbornik 1970; 10:217–243.
P. Kútik and K. Mikula. Diamond-cell finite volume scheme for the Heston model. Discrete
Contin. Dyn. Syst. Ser. S 2015; 8(5):913–931.
N. N. Kuznetsov. Accuracy of some approximate methods for computing the weak solutions of a
first-order quasi-linear equation. USSR, Comput. Math. and Math. Phys. 1976; 16(6):159–
193.
P. D. Lax and B. Wendroff. Systems of conservation laws. Comm. Pure Appl. Math. 1960;
13:217–237.
P. D. Lax. Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock
Waves. SIAM, Philadelphia, Penn., 1973.
R.D. Lazarov, I.D. Michev, and P.S. Vassilevsky. Finite volume methods for convection-
diffusion problems. SIAM J. Numer. Anal. 1996; 33:31–35.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 103

R. Lazarov and S. Tomov. A posteriori error estimates for finite volume element
approximations of convection-diffusion-reaction equations. Comput. Geosci. 2002; 6(3-
4):483–503.
C. Le Potier. Schéma volumes finis monotone pour des opérateurs de diffusion fortement
anisotropes sur des maillages de triangles non structurés. Comptes Rendus Mathématique
2005; 341(12):787–792.
C. Le Potier. Finite volume scheme satisfying maximum and minimum principles for
anisotropic diffusion operators. In Finite volumes for complex applications V, ISTE,
London 2008, pp. 103–118.
C. Le Potier. A nonlinear finite volume scheme satisfying maximum and minimum principles
for diffusion operators. Int. J. Finite Vol. 2009; 6(2).
R. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press
2002.
R. . J. LeVeque and M. Pelanti. A class of approximate Riemann solvers and their relation to
relaxation schemes. J. Comput. Phys. 2001; 172(2):572–591.
M. S. Liou and C. J. Steffen. A new flux-splitting scheme. J. Comput. Phys. 1993; 107:23–39.
K. Lipnikov, G. Manzini, and M. Shaskov. Mimetic finite difference method. J. Comput. Phys.
2014; 257(part B):1163–1227.
X.-D. Liu. A maximum principle satisfying modification of triangle based adaptive stencils
for the solution of scalar hyperbolic conservation laws. SIAM J. Numer. Anal. 1993;
30:701–716.
J. Málek, J. Nečas, M. Rokyta, and M. Røužička. Weak and measure-valued solutions
to evolutionary PDEs. Applied Mathematics and Mathematical Computation, vol.
13. Chapman and Hall: London, 1968; 44-177.
S. Martin and J. Vovelle. Convergence of implicit finite volume methods for scalar conservation
laws with discontinuous flux function. M2AN Math. Model. Numer. Anal. 2008; 42(5):699–
727
J.-M. Masella, I. Faille, and T. Gallouët. On an approximate Godunov scheme. Int. J. Comput.
Fluid Dyn. 1999; 12(2):133–149.
B. Merlet L∞ - and L2 -Error estimates for a finite volume approximation of linear advection
SIAM J. Numer. Anal. 2007; 46(1):124–150.
A. Michel, Q.H. Tran and G. Favennec. A genuinely one-dimensional upwind Scheme with
accuracy enhancement for multidimensional advection problems. ECMOR XII 12 th
European Conference on the Mathematics of Oil Recovery, 2010.
K. Mikula and N. Ramarosy. Semi-implicit finite volume scheme for solving nonlinear diffusion
equations in image processing. Numerische Mathematik 2001; 89(3):561–590.
M. S. Mock. Systems of conservation laws of mixed type. J. Diff. Eqns. 1980; 37:70–88.
R. Nicolaı̈des. Analysis and convergence of the MAC scheme I: The linear problem. SIAM J.
Numer. Anal. 1992; 29:1579–1591.
R. Nicolaı̈des and X. Wu. Analysis and convergence of the MAC scheme II, Navier-Stokes
equations. Math. Comp. 1996; 65:29–44.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
104 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

M. Ohlberger. A posteriori error estimates for finite volume approximations to singularly

perturbed nonlinear convection-diffusion equations. Numer. Math. 2001a; 87(4):737–761.
M. Ohlberger. A posteriori error estimates for vertex centered finite volume approximations
of convection-diffusion-reaction equations. M2AN Math. Modell. Numer. Anal. 2001b;
35(2):355–387.
M. Ohlberger and C. Rohde. Adaptive finite volume approximations for weakly coupled
convection dominated parabolic systems. IMA J. Numer. Anal. 2002; 22(2):253–280.
M. Ohlberger and J. Vovelle. Error estimate for the approximation of nonlinear conservation
laws on bounded domains by the finite volume method. Math. Comp. 2006; 75(253):113–
150 (electronic).
M. Ohlberger. A review of a posteriori error control and adaptivity for approximations of
non-linear conservation laws. Internat. J. Numer. Methods Fluids 2009; 59(3):333–354.
O. A. Oleinik. Discontinuous solutions of non-linear differential equations. Amer. Math. Soc.
Transl.(2) 1963; 26:95–172.
P. Omnes. An a posteriori error bound for the discrete duality finite volume discretization of
the Laplace equation and application to the adaptive simulation in a domain with a crack.
newblock In Finite volumes for complex applications V ISTE, London 2008, pp. 617–624.
P. Omnes. On the second-order convergence of a function reconstructed from finite volume
approximations of the Laplace equation on Delaunay-Voronoi meshes. ESAIM Math.
Model. Numer. Anal. 2011; 45(4):627–650.
S. Osher and F. Solomon. Upwind difference schemes for hyperbolic systems of conservation
laws. Math. Comput. 1982; 38(158):339–374.
S. Osher. Riemann solvers, the entropy condition, and difference approximations. SIAM J.
Numer. Anal. 1984; 21(2):217–235.
S. Patankar. Numerical heat transfer and fluid flow. Series in Computational Methods
in Mechanics and Thermal Sciences, vol. XIII. Washington - New York - London:
Hemisphere Publishing Corporation; New York. McGraw-Hill Book Company 1980.
M. Pelanti, F. Bouchut, and A. Mangeney. A Riemann solver for single-phase and two-phase
shallow flow models based on relaxation. Relations with Roe and VFRoe solvers. J.
Comput. Phys., 230(3):515–550, 2011.
B. Perthame, N. Seguin, Nicolas and M. Tournus. A simple derivation of BV bounds for
inhomogeneous relaxation systems. Commun. Math. Sci., 13(2):577–586, 2015.
L. Piar, F. Babik, R. Herbin, and J.-C. Latché. A formally second order cell centered scheme
for convection-diffusion equations on general grids. International Journal for Numerical
Methods in Fluids 2013; 71:873–890.
T. A. Porsching. Error estimates for MAC-like approximations to the linear Navier-Stokes
equations. Numer. Math. 1978; 29(3):291–306.
P. L. Roe. Approximate Riemann solvers, parameter vectors, and difference schemes. J.
Comput. Phys. 1981; 43:357–372.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
 FINITE VOLUME METHODS: FOUNDATION AND ANALYSIS 105

P. Rostand and B. Stoufflet. TVD schemes to compute compressible viscous flows on

unstructured meshes. Proceedings of the Second International Conference on Hyperbolic
Problems. Vieweg: Braunschweig, 1988.
E. Süli. Convergence of finite volume schemes for Poisson’s equation on nonuniform meshes.
SIAM J. Numer. Anal. 1991; 28:1419–1430.
C. W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing
schemes. J. Comput. Phys. 1988; 77:439–471.
C. W. Shu. Total-variation-diminishing time discretizations. SIAM J. Sci. Statist. Comput.
1988; 9:1073–1084.
C. W. Shu. High Order ENO and WENO Schemes, Lecture Notes in Computational Science
and Engineering, vol. 9. Springer-Verlag: Heidelberg, 1999.
C. W. Shu. A survey of strong stability preserving high order time discretizations. Collected
lectures on the preservation of stability under discretization (Fort Collins, CO, 2001),
SIAM: Philadelphia, 2002; 51-65.
T. Sonar. On the construction of essentially non-oscillatory finite volume approximations to
hyperbolic conservation laws on general triangulations: polynomial recovery, accuracy, and
stencil selection. Comput. Methods Appl. Mech. Eng. 1997; 140:157–181.
T. Sonar. On families of pointwise optimal finite volume ENO approximations. SIAM J.
Numer. Anal. 1998; 35(6):2350–2379.
S. P. Spekreijse. Multigrid solution of monotone second-order discretizations of hyperbolic
conservation laws. Math. Comput. 1987; 49:135–155.
J. L. Steger and R. F. Warming. Flux vector splitting of the inviscid gasdynamic equations
with application to finite difference methods. J. Comput. Phys. 1981; 40:263–293.
P. K. Sweby. High resolution schemes using flux limiters for hyperbolic conservation laws.
SIAM J. Numer. Anal. 1984; 21(5):995–1011.
A. Szepessy. An existence result for scalar conservation laws using measure valued solutions.
Commun. Partial Diff. Equations 1989; 14:1329–1350.
E. Tadmor. Entropy functions for symmetric systems of conservation laws. J. Math. Anal.
Appl. 1987; 122:355–359.
E. Tadmor. Local error estimates for discontinuous solutions of nonlinear hyperbolic equations.
SIAM J. Numer. Anal. 1991; 28:891–906.
E. Tadmor. Entropy stability theory for difference approximations of nonlinear conservation
laws and related time-dependent problems. Acta Numerica 2004; 13:451–512.
L. Tartar. Compensated compactness and applications to partial differential equations,
Nonlinear analysis and mechanics: Heriot-Watt Symposium, Vol. IV, Res. Notes in Math..
Pitman: Boston, Mass.-London, 1979.
S. Tatsumi, L. Martinelli, and A. Jameson. Design, Implementation, and Validation of Flux
Limited Schemes for the Solution of the Compressible Navier-Stokes Equations. Report
AIAA-94-0647, American Institute for Aeronautics and Astronautics, 1994.
N. Therme 2015. Schémas numériques pour la simulation de l’explosion. Ph.D. thesis, Aix-
Marseille Université.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.
106 ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS

Q. Tran. Second-order slope limiters for the simultaneous linear advection of (not so)
independent variables. Commun. Math. Sci. 2008; 6(3):569–593.
B. van Leer. Towards the ultimate conservative difference schemes V. A second order sequel
to Godunov’s method. J. Comput. Phys. 1979; 32:101–136.
B. van Leer. Flux-Vector Splitting for the Euler Equations. Report ICASE-82-30, Institute for
Computer Applications in Science and Engineering, NASA Langley, 1982.
B. van Leer. Upwind-Difference Schemes for Aerodynamics Problems Governed by the
Euler Equations, Lectures in Applied Mathematics 22, AMS: Providence, Rhode Island,
1985.
P. Vankeirsblick. Algorithmic Developments for the Solution of Hyperbolic Conservation Laws
on Adaptive Unstructured Grids. PhD thesis, Katholieke Universiteit Leuven, Belgium,
1993.
J. P. Vila. Convergence and error estimates in finite volume schemes for general multi-
dimensional scalar conservation laws I: Explicit monotone schemes. RAIRO, Model. Math.
Anal. Numer. 1994; 28:267-295.
C. Viozat, C. Held, K. Mer, and A. Dervieux. On Vertex-Center Unstructured Finite-Volume
Methods for Stretched Anisotropic Triangulations. Report 3464, Institut National De Rec-
herche En Informatique Et En Automatique (INRIA),
1998.
J. Vovelle. Convergence of finite volume monotone schemes for scalar conservation laws on
bounded domains. Numer. Math. 2002; 90(3):563–596.
P. Wesseling. Principles of Computational Fluid Dynamics. Springer 2001.
M. Wierse. Higher Order Upwind Scheme on Unstructured Grids for the Compressible Euler
Equations in Time Dependent Geometries in 3D. PhD thesis, University of Freiburg,
Germany, 1994.

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes.
c 2016 John Wiley & Sons, Ltd.