FINITE VOLUME METHODS

LONG CHEN

The finite volume method (FVM) is a discretization technique for partial differential
equations, especially those that arise from physical conservation laws. FVM uses a volume
integral formulation of the problem with a finite partitioning set of volumes to discretize
the equations. FVM is in common use for discretizing computational fluid dynamics equations.
1. G ENERAL FORM OF FINITE VOLUME METHODS
We consider vertex-centered finite volume methods for solving diffusion type elliptic
equation
(Ku) = f in ,

(1)

with suitable Dirichlet or Neumann boundary conditions. Here Rd is a polyhedral

domain (d 2), the diffusion coefficient K(x) is a d d symmetric matrix function that
is uniformly positive definite on with components in L (), and f L2 (). We have
discussed finite element methods based on the discretization of the weak formulation and
finite difference methods based on the classic formulation. We shall now present finite
volume methods based on the following balance equation
Z
Z
(2)

(Ku) n ds = f dx, b ,
b

where n denotes the unit outwards normal vector of b.

Finite volume methods are discretizations of the balance equation (2). The discretization consists of three approximations:
(1) approximate the function u by uh in a N -dimensional space V;
(2) approximate arbitrary domain b by a finite subset B = {bi , i = 1 : M };
(3) approximate boundary flux (Ku) n on bi by a discrete one (Kh uh ) n.
We then end with a method: to find uh V such that:
Z
Z
(3)

(Kh uh ) n dS =
f dx, bi , i = 1 : M.
bi

bi

We call any method in the form (3) finite volume methods (FVMs).
Since finite volume methods discretize the balance equation (2) directly, an obvious
virtue of finite volume methods is the conservation property comparing with finite element
methods based on the weak formulation. This property can be fundamental for the simulation of many physical models, e.g., in oil recovery simulations and in computational fluid
dynamics in general.
On the other hand, the function space and the control volume can be constructed based
on general unstructured triangulations for complex geometry domains. The boundary condition can be easily built into the function space or the variational form. Thus FVM is
more flexible than standard finite difference methods which mainly defined on the structured grids of simple domains.
1

LONG CHEN

2. C ELL - CENTERED FINITE VOLUME METHOD

Let T be a triangular or Cartesian grid of . We choose the finite dimensional space
V = {v L2 () : v| is constant for all T }. Then dim V = N T , the number of
elements of T . We also choose B = T . See
Figure 2(a).
To complete the discretization,
FINITE VOLUME
METHODS:
7
FINITE VOLUME METHODS: FOUNDATION
AND
ANALYSISFOUNDATION AND
7 ANALYSIS
we need to assign the boundary flux of each element.
This can be done in a finite difference
fashion.
For example,
for
an interior side e (an
2. Finite
volume (FV)
methodslaws
for nonlinear conservation laws
2. Finite volume (FV) methods
for nonlinear
conservation
edgeIninthe2-D
and
a
face
in
3-D)
shared
by
two
elements

and

,
we
can define
1Rd , is first
2
Inthe
thecomputational
finite volume method,
domain,

aRd , is first tessellated into a

finite volume method,
domain, the
computational
tessellated
into
of non overlapping
control
volumes
that completely
cover the domain. Notationally,
collection of non overlappingcollection
control volumes
thatucompletely
cover
the
domain.
Notationally,
|

u
|
h 2
h 1
denote
a with
tessellation
the domain
with that
control
volumes
(4) let T denote a tessellation of let

theT
domain
volumes
T , T such
T T
T = .T T such that T T T = .
h uh n
e :=controlof
cscale
c1 T ,with
2
hT denote
a length
associated
each
volume
Let hT denote a length scaleLet
associated
with
each control
volume
e.g. h
diam(T).
For T , e.g. hT diam(T). For
T control
distinct
control
volumes
Ti and
Tvector
T , oriented
the
intersection
is either
an oriented edge (2-D)
two the
distinct
control
volumes
T
and
T
in
T
,
the
intersection
is
either
edge
(2-D)
j in an
i
j
where
normal
vector
ntwo
is
the
outward
unit
normal
of
e
in

,
i.e.
pointing
e
1
ornormal
face (3-D)
eijelse
with
oriented
normalatmost
a set
of measure
or face (3-D) eij with oriented
ij or
a set
of measure
d 2.
In each
controlat most d 2. In each control
ij or else
fromvolume,
1 to an
2 integral
and cconservation
i ,volume,
i =law
1, an
2statement
are
points
in each law
element
suchis that
the line segment
i
integral
statement
then imposed.
isconservation
then imposed.
connecting c2 and c1 is orthogonal to the side e. By the symmetry, for rectangles or
cubes c is the mass center
of . For simplex,
c should belaw)
the circumcenters
whichlaw asserts that the
Definition 2.1
(Integral
conservation
integral
Definition 2.1 (Integral conservation
law)
An integral
conservation lawAn
asserts
thatconservation
the
imposes
restriction
on
the
triangulation.
The
error
analysis
can
be
carried
out
in
the
finitecontrol volume T is
rate
of
change
of
the
total
amount
of
a
substance
with
density
u
in
a
fixed
rate of change of the total amount of a substance with density u in a fixed control volume T is
to the
the
total
flux
the
substance
through the
boundary
T system.
equal to fashion
the total flux
of the equal
substance
through
the of
boundary
T stability
difference
by considering
truncation
error
and
of the
resulting
Theory and computation along this approach is summarized in the book [9].
Another approach to discretize the boundary flux is through
mixed
finite element meth!
!
!
!
2 d
1
ods. The gradient operator is dunderstood
as

:
L

H
.
Optimal
can
f (u) error
d = 0estimate
.
(15)
u dx +
f (u) d = 0 . u dx +
(15)
dt T
dt of
T mixed T
be easily obtained by using that
finite element
methods T
[14].
Since the control volume is the element (also called cell) of the mesh and the unknown
This
integral conservation
law statement
is readily
obtained
upon spatial integration of the
This integral
law
statement
is readily
obtained
upon spatial
integration
ofmethods
the
is associated
to conservation
each element/cell,
it is often
called
cell-centered
finite
volume
and
(1a) in theofregion
T and application
the divergence theorem. The
divergence equation (1a) indivergence
the region equation
T and application
the divergence
theorem. of
The
the choice
difference
scheme
(4)
is
also
known
as
cell
centered
finite
difference
methods.
choice ofiscontrol
tessellation
flexible in
theexample,
finite volume
of control volume tessellation
flexiblevolume
in the finite
volumeis method.
For
Fig. method. For example, Fig.

storage location

storage location

control volume

control volume

(a) Mesh and dual mesh of cell-centered FVM

(b) Mesh and b.
dual
mesh of vertex-centered FVM
Cell-centered
Vertex-centered
a. Cell-centered
b.a.Vertex-centered
Figure
1. in
Control
volume
used in the finite volume method:
1. Control
used
the
finite
volume
method:
F IGUREFigure
1. Mesh
andvolume
dual variants
mesh
of
two
FVMs.
Thevariants
unknowns
are asso(a) cell-centered
(b) vertex-centered
control volume tessellation.
(a) cell-centered and (b) vertex-centered
controland
volume
tessellation.
ciated to black nodes.

1 depicts
2-Dtypical
triangle
complex
and tessellations
two typical control
1 depicts a 2-D triangle complex
and atwo
control
volume
(amongvolume
many tessellations (among many
others)
used inInthe
volume
method.
In themethod
cell-centered
others) used in the3.finite
volume
thefinite
cell-centered
finite volume
shownfinite
in volume method shown in
V ERTEX
-method.
CENTERED
FINITE
VOLUME
METHOD
Fig. serve
1a, theastriangles
themselves
as control
volumes
with of
solution unknowns (degrees of
Fig. 1a, the triangles themselves
control volumes
withserve
solution
unknowns
(degrees
freedom)
stored
onvertex-centered
a per
In the vertex-centered
freedom)
on another
a per triangle
basis.
In the
finite
volume
method
shown finite
in wevolume
We
now stored
discuss
popular
choice
of
V triangle
and B. basis.
To simplify
the notation,
con- method shown in
1b, as
control
volumesdual
are to
formed
as a geometric
the triangle complex and solution
Fig. 1b, control volumes are Fig.
formed
a geometric
the triangle
complexdual
and to
solution
sider
two dimensional triangular
grids
and homogenous
Dirichlet
condition. We
unknowns
stored
a per triangulation
vertexboundary
basis.
unknowns stored on a per triangulation
vertexonbasis.

refer to [16] for the general treatment in high dimensional simplicial grid and [15] for
Encyclopedia
of Computational
Mechanics.
Edited
by Thomas
Erwin Stein,
Rene de Borst and Thomas J.R. Hughes.
Encyclopedia
of Computational Mechanics.
Edited
by Erwin Stein,
Rene de Borst
and
J.R. Hughes.
rectangle
grids.
c 2004 John Wiley & Sons, Ltd.
c 2004 John Wiley
!
& Sons, Ltd.!
2
Let R be a polygon and let T be a triangular grid of . Denoted by VT be the
linear finite element spaces of H01 () based on T :
VT = {v H01 () : v| P1 ( ), T },
where P1 ( ) is the linear polynomial space on . We shall choose V = VT . The dimension
N is the number of interior vertices of T .

FINITE VOLUME METHODS

The control volume will be given by another mesh B = {bi , i = 1, , M } satisfying

= M

i=1 bi , and bi bj = , i 6= j,
and to reflect to the Dirichlet boundary condition, we set

bi }.
B = {bi B,
The element bi of B is not necessary to be polygons. But for practical reasons, bi are
chosen as polygons such that the boundary integral is easy to evaluate.
Given a triangulation T , one construction of B is given as follows: for each triangle
T , select a point c . The point c can coincides with middle points of edges,
but not the vertices of triangles (to avoid the degeneracy of the control volume). In each
triangle, we connect c to three middle points on the boundary edges. This will divide
each triangle in T into three regions. For each vertex xi of T , we collect all regions
containing this vertex and define it as bi . In Figure 2 we only draw the control volume
for interior vertices since the Dirichlet boundary condition is build into the space VT and
the unknown is only associated to interior vertices. Obviously for Neumann boundary

condition, we should use B.

There are three common choices of c :
Type A: c is the barycenter of .
Type B: c is the middle point of the longest edge.
Type C: c is the circumcenter of .
Type A is preferable for triangulations composed by equilateral triangles. In this case
will be divided into three parts with equal area. This symmetric property is important to
get optimal rate of convergence in L2 norm. Type B is better for right triangles, and can
be easily obtained by the longest edge bisection method. Type C is suitable for Delaunay
triangulations. The edges of the control volumes will be orthogonal to the intersected edges
of triangles, and if the grid T is a Delaunay triangulation, B will be a Voronoi diagram.

(a) Type A

(b) Type B

(c) Type C

F IGURE 2. Three types of grids and dual grids. The gray areas are the
control volumes of interior nodes. Type A: The point c is the barycenter
of . Type B: The points c is the middle point of the longest edge. Type
C: The point c is the circumcenter of .
Since we associate control volumes and unknowns to vertices, it s called vertex-centered
finite volume methods. It is also known as box method [1, 10] (since the control volume
is called box in these work), finite volume element methods [3, 2, 11] (to emphasis the
approximation of u is from finite element space), and generalized finite difference methods
[13, 12]. High order finite volume methods can be found in [5, 12].
1

LONG CHEN

4. P ETROV-G ALERKIN FORMULATION

We shall follow Bank and Rose [1] to formulate the vertex-centered linear finite volume
method as a Petrov-Galerkin method by choosing different trial space and test space for an
appropriate bilinear form. Other approaches can be found at [10, 3, 2, 15, 12, 6, 7, 4].
We first introduce a function space for the dual mesh. Let B be the dual mesh of a
triangulation T constructed in the previous subsection. We define a piecewise constant
function space on B by:
(5)

VB = {v L2 () : v|bi = constant , bi B}.

The set of interior sides of the mesh B is denoted by E(B). For each e E(B), we shall
fix a unit normal direction ne of e. That is ne is independent of the element containing
e. Suppose e is shared by two control volumes bi and bj . Without loss of generality,
we suppose the outward normal direction of e in bi coincides with ne . For any function
v VB , the jump of v across e is denoted by [v] = v|bi v|bj .
We define a bilinear form on VT and VB
X Z
(6)
a
(u, v) =
(Ku) ne [v]dS, u VT , v VB ,
eE(B)

and formulate the linear finite volume method as: find u

VT such that
(7)

a
(
u, v) = (f, v) for all v VB .

Remark 4.1. For Neumann boundary condition, we shall choose

VT = {v H 1 () : v| P1 ( ), T }, and

VB = {v L2 () : v|b = constant, bi B}.

i

For e bi , the flux (Ku) ne will be given by the boundary condition and can
be moved to the right hand side. All algorithms and analysis in this notes can be applied
to Neumann boundary condition in a straightforward way.
Since the trial space VT and the test space VB are different, the weak formulation (7) is
known as Petrov-Galerkin method.
5. R ELATION TO F INITE E LEMENT M ETHODS
We shall recall the finite element method and show the close relation between them. Let
a(u, v) be the bilinear form
Z
(8)
a(u, v) = (Ku) v dx.

The linear finite element method is: find uL VT such that

(9)

a(uL , v) = (f, v) for all v VT .

For FEM, the trial space and the test space are the same, which is known as the Galerkin
method.
To see the close relation, we now formulate the corresponding matrix equations for (7)
and (9). Let N (T ) be the set of interior nodes of T and N = #N (T ). Then dim VB =
dim VT = N . A basis of VB can be chosen as the characteristic function of each bi , i =
1, , N :
(
1 x bi ,
i = bi (x) =
0 otherwise .

FINITE VOLUME METHODS

The nodal basis of linear finite element space VT is the standard hat function:
i VT , i (xj ) = ij , xj N (T ), i = 1, , N.
Let u
=
equation

PN

j=1

j j . Choosing v = i , i = 1, , N in (7), we obtain a linear algebraic

U
= F ,
AU

(10)
with
Aij =

(Kj ) n, Fi =

f dx.

bi

bi

PN
Let uL = j=1 Uj j . Choosing v = i , i = 1, , N in (9), we obtain another linear
algebraic equation
AU = F,

(11)
with

Z
(Kj ) i , Fi =

Aij =

f i dx.

We shall prove that when K(x) is piecewise constant on each triangle, then A = A;
See [1, 10, 16]. The solution vectors are point values for uL and u
at vertices.
R The only
difference is the different way to compute the right hand side. For FEM, Fi = i f i dx,
R
is a weighted average over the star of a vertex. For FVM, Fi = bi f dx is the average over
the control volume bi . When we choose type A control volume, i.e. choosing c to be the
barycenter of , Fi can be thought as an approximation of Fi using mass lumping. This
modification enables the solution of linear FVM to satisfy the conservation property. It is
interesting to note that on uniform grids, three methods (FDM, FEM and FVM) result the
same matrix (up to a scaling). The right hand side are chosen from very different perspective. And amazingly for all three choices of right-hand side, the resulting approximation
converges to the same solution with the same order.
Lemma 5.1. Given a polyhedron with k-sides in Rn , let |Fi | denote the (n1)-measure
of the face Fi and ni the outward normal of the i-th side. Then
k
X

(12)

|Fi |ni = 0.

i=1

Proof. We pick up any interior point x and define i (x) as the simplex spanned
by Fi and x. Note that |i (x)| is a linear function of x and |i (x)| = 0 for x Fi .
Thus the direction of the vector i is the normal direction of the plane i (x) = 0. The
magnitude can be computed by a simple difference formula of directional derivatives to
yield |i (x)| = |Fi |ni .
Pk
We then differentiate i=1 |i (x)| = || to get
k
X
i=1

|Fi |ni =

k
X

|i (x)| = 0.

i=1


The following theorem is critical which says FVM and FEM have the stiffness matrix.
Let us introduce an isomorphism between linear spaces G : VB VT by mapping i
PN
PN
i , 1 i N . Then for any u = i=1 ui i VB , Gu = i=1 ui i VT . Note that

LONG CHEN

u and Gu share the same vector representation U = (u1 , , uN )T . We also use a simple
notation u
to denote Gu.
Theorem 5.2. Assume K(x) is piecewise constant on each T and bi consists
of middle points of edges, then
a(u, v) = a
(u, v),

u, v VT .

Proof. Since we assume K(x) is piecewise constant, we need only show the local stiffness
for Poisson equation on one triangle coincides. That reduces to prove
a
(i , Gj ) = a(i , j ).
To be specific, let us take 1 , 2 as an example. Since 1 is a constant over the triangle
, we can pull it out the integral and apply Lemma 12 twice to get
Z
1 (n1 + n2 )dS = 1 (|e1 |ne1 + |e2 |ne2 )
e1 e2
Z
1
1
= 1 (|l1 |n1 + |l3 |n3 ) = 1 |l2 |n2 = 1 2 .
2
2


6. E RROR ANALYSIS
The error analysis of vertex centered linear FVM (7) relies on the close relation between
linear finite element method and the linear finite volume method. From previous section,
linear FVM approximation u
can be thought as a perturbation of the FEM approximation
uL . First order optimal convergence rate in the energy norm can be obtained using this
relation.
Note that the right hand sides may be quite different for type B dual mesh. For example,
let f = 1 and consider the control volume in Figure 2(b). Then Fi = |i |/3 while Fi =
|i |/4. Nerveless optimal first order convergence in H 1 norm can be still derived by
comparing them in H 1 norm [10].
Furthermore if we choose type A dual mesh, then optimal second order convergence in
L2 -norm can be also derived [10, 4, 7, 8]. Note that for general choice of control volumes,
finite volume approximation may not lead to optimal L2 -norm convergent rate. See [11]
for such an example on type B dual mesh. Optimal L norm estimate can be also obtained
by treating it as a perturbation of finite element methods; See [7, 8]. We are not going to
discuss L2 or L error estimate.
Theorem 6.1. Assume K(x) is piecewise constant on each T , the solution u of the
diffusion equation is in the space H01 () H 2 () and the mesh is quasi-uniform with
mesh size h, then the finite volume approximation uh has optimal approximation order
(13)

ku uh k1, . h(kuk2, + kf k).

Proof. For any f L2 (), we define h f V0T as

hh f, vh i = (f, Gvh ),

for all vh VT .

and Qh f V0T as
hQh f, vh i = (f, vh ),

for all vh VT .

FINITE VOLUME METHODS

Following the notation of Hackbusch [10], we denoted by uG

h as the standard Galerkin
(finite element) approximation and uB
is
the
box
(finite
volume)
approximation. The
h
equivalence of the stiffness matrices means
Lh uG
h = Qh f,

Lh uB
h = h f.

Therefore by the stability of L1

h , we have
B
|uG
h uh |1 = sup

vh VT

hQh f h f, vh i
.
|vh |1

By the definition
hQh f h f, vh i = (f, vh Gvh ).
Denote the support the hat basis function at xi as i . Note that bi i and the operator
I G preserve constant function in the patch i and thus
(f, vh Gvh )bi kf kbi kvh Gvh ki Chkf kbi |vh |1,i .
Summing up and using Cauchy Schwarz inequality, we get the first order convergence
B
|uG
h uh |1 Chkf k.

The estimate (13) the follows from the triangle inequality and the estimate of the finite
element method.

R EFERENCES
[1] R. E. Bank and D. J. Rose. Some error estimates for the box method. SIAM J. Numer. Anal., 24:777787,
1987.
[2] Z. Cai. On the finite volume element method. Numer. Math., 58(1):713735, 1990.
[3] Z. Cai, J. Mandel, and S. F. McCormick. The finite volume element method for diffusion equations on
general triangulations. SIAM J. Numer. Anal., 28:392402, 1991.
[4] P. Chatzipantelidis. Finite volume methods for elliptic PDEs: a new approach. M2AN Math. Model. Numer.
Anal., 36(2):307324, 2002.
[5] L. Chen. A new class of high order finite volume methods for second order elliptic equations. SIAM J.
Numer. Anal., 47(6):40214043, 2010.
[6] S.-H. Chou, D. Y. Kwak, and Q. Li. Lp error estimates and superconvergence for covolume or finite volume
element methods. Numer. Methods Partial Differential Equations, 19(4):463486, 2003.
[7] S.-H. 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., 69(229):103120, 2000.
[8] R. E. Ewing, T. Lin, and Y. Lin. On the accuracy of the finite volume element method based on piecewise
linear polynomials. SIAM J. Numer. Anal., 39(6):18651888, 2002.
[9] R. Eymard, T. Gallouet, and R. Herbin. Finite volume methods. In Handbook of numerical analysis, Vol.
VII, Handb. Numer. Anal., VII, pages 7131020. North-Holland, Amsterdam, 2000.
[10] W. Hackbusch. On first and second order box schemes. Computing, 41(4):277296, 1989.
[11] J. Huang and S. Xi. On the finite volume element method for general self-adjoint elliptic problems. SIAM J.
Numer. Anal., 35(5):17621774, 1998.
[12] R. Li, Z. Chen, and W. Wu. Generalized difference methods for differential equations, volume 226 of Monographs and Textbooks in Pure and Applied Mathematics. Marcel Dekker Inc., New York, 2000. Numerical
analysis of finite volume methods.
[13] R. H. Li. On the generalized difference method for elliptic and parabolic differential equations. In Proc. of
the Symposium on the Finite Element Method between China and France, Beijing, China, 1982.
[14] T. Russell and M. Wheeler. Finite element and finite difference method for continuous flows in porous
media. Frontiers in Applied Mathematics, 1:35, 1983.
[15] E. Suli. Convergence of finite volume schemes for Poissons equation on nonuniform meshes. SIAM J.
Numer. Anal., 28(5):14191430, 1991.
[16] J. Xu and Q. Zou. Analysis of linear and quadratic simplicial finite volume methods for elliptic equations.
Numer. Math., 111(3):469492, 2009.