A critical analysis of some popular methods for the discretisation of
the gradient operator in finite volume methods
Alexandros Syrakos∗1 , Stylianos Varchanis1 , Yannis Dimakopoulos1 , Apostolos Goulas2 , and
John Tsamopoulos1
1

arXiv:1606.05556v6 [cs.NA] 29 Dec 2017

Laboratory of Fluid Mechanics and Rheology, Dept. of Chemical Engineering, University of Patras, 26500
Patras, Greece
2
Laboratory of Fluid Mechanics and Turbomachinery, Dept. of Mechanical Engineering, Aristotle University
of Thessaloniki, 54124 Thessaloniki, Greece

Abstract
Finite volume methods (FVMs) constitute a popular class of methods for the numerical simulation of fluid flows. Among the various components of these methods, the discretisation of the
gradient operator has received less attention despite its fundamental importance with regards to
the accuracy of the FVM. The most popular gradient schemes are the divergence theorem (DT) (or
Green-Gauss) scheme, and the least-squares (LS) scheme. Both are widely believed to be secondorder accurate, but the present study shows that in fact the common variant of the DT gradient
is second-order accurate only on structured meshes whereas it is zeroth-order accurate on general
unstructured meshes, and the LS gradient is second-order and first-order accurate, respectively.
This is explained through a theoretical analysis and is confirmed by numerical tests. The schemes
are then used within a FVM to solve a simple diffusion equation on unstructured grids generated
by several methods; the results reveal that the zeroth-order accuracy of the DT gradient is inherited by the FVM as a whole, and the discretisation error does not decrease with grid refinement.
On the other hand, use of the LS gradient leads to second-order accurate results, as does the use
of alternative, consistent, DT gradient schemes, including a new iterative scheme that makes the
common DT gradient consistent at almost no extra cost. The numerical tests are performed using
both an in-house code and the popular public domain PDE solver OpenFOAM.

This is the accepted version of the article published in:
Physics of Fluids 29, 127103 (2017); doi:10.1063/1.4997682

1

Introduction

Finite volume methods (FVMs) are used widely for the simulation of fluid flows; they are employed by
several popular general-purpose CFD (Computational Fluid Dynamics) solvers, both commercial (e.g.
ANSYS Fluent, STAR-CD) and open-source (e.g. OpenFOAM). One of the key components of FVMs
is the approximation of the gradient operator. Computing the gradient of the dependent variables is
needed for the FVM discretisation on non-Cartesian grids, where the fluxes across a face separating
two finite volumes cannot be expressed as a function of the values of the variables at the centres of
these two volumes alone, but additional terms involving the gradients must also be included [1]. The
gradient operator is even more significant when solving complex flow problems such as turbulent flows
modelled by the RANS and some LES methodologies [1] or non-Newtonian flows [2–5], where the
additional equations solved (turbulence closure, constitutive equation etc.) may directly involve the
velocity gradients. Apart from the main task of solving partial differential equations (PDEs), gradient
calculation may also be important in auxiliary activities such as post-processing [6].
∗

Corresponding author. E-mail address: alexandros.syrakos@gmail.com, syrakos@upatras.gr

1

The two most popular methods for calculating the gradient on grids of general geometry are
based on the use of (a) the divergence theorem (DT) and (b) least-squares (LS) minimisation. They
have both remained popular for nearly three decades of application of the FVM; for example, the
DT method (also known as “Green–Gauss gradient”) has been used in [1, 7–11], and the LS method
in [10–15]. Often, general-purpose commercial or open-source CFD codes present the user with the
option of choosing between these two methods. Their popularity stems from the fact that they are
not algorithmically restricted to a particular grid cell geometry but can be applied to cells with an
arbitrary number of faces. This property is important, because in recent years the use of unstructured
grids is becoming standard practice in simulations of complex engineering processes. The tessellation
of the complex geometries typically associated with such processes by structured grids is an arduous
and extremely time-consuming procedure for the modeller, whereas unstructured grid generation is
much more automated. Therefore, discretisation schemes are sought that are easily applicable on
unstructured grids of arbitrary geometry whereas early FVMs either used Cartesian grids or relied on
coordinate transformations which are applicable only on smooth structured grids.
In the literature, usually the gradient discretisation is only briefly discussed within an overall
presentation of a FVM, with only a relatively limited number of studies devoted specifically to it
(e.g. [6, 16–20]). This suggests that existing gradient schemes are deemed satisfactory, and in fact
there seems to be a widespread misconception that the DT and LS schemes are second-order accurate
on any type of grid. For the DT gradient, a couple of studies have shown that it is potentially zerothorder accurate, depending on the grid properties. Syrakos [21] noticed that it converged to wrong
values on composite Cartesian grids in the vicinity of the interfaces between patches of different
fineness. His analysis, given here in Sec. 5.4, concludes that this is due to grid skewness. Later,
Sozer et al. [20] tested a simpler variant of the DT gradient that uses arithmetic averaging instead
of linear interpolation and proved that in the one-dimensional case it converges to incorrect values
if the grid is not uniform. In numerical tests they also noticed that the scheme is inconsistent on
two-dimensional grids of arbitrary topology. In fact, as early as in [7] it was briefly mentioned that
the DT gradient fails to calculate exactly the gradients of linear functions, providing a hint to its
inconsistency. An important question is whether this inconsistency inhibits convergence of the FVM
to the correct solution. The results reported in [10] concerning the solution of a Poisson equation
suggest zeroth-order FVM accuracy when the DT gradient is employed on skewed meshes, although
the authors did not explicitly attribute this to inconsistency of the DT scheme.
The present paper presents mathematical analyses of the orders of accuracy of both methods on
various types of grid, backed by numerical experiments. The common DT method is proved to be
zeroth-order accurate on grids of arbitrary geometry and second-order accurate on smooth structured
meshes. This is shown to be true even if a finite number of “corrector” steps is performed (e.g. [10,22])
whereby an iterative procedure uses the gradient calculated at the previous iteration to improve the
accuracy. Theoretically, this procedure increases the order of accuracy only in the limit of infinite
iterations. Practically, the desired accuracy can be achieved with a finite number of iterations, but
this number is a priori unknown and increases with the grid fineness. Furthermore, corrector steps
make the DT gradient much more expensive than the LS. However, since the FVM solution typically
proceeds with a number of outer iterations, we exploit this fact to interweave the gradient iterations
with the FVM outer iterations to obtain a first-order accurate DT gradient at a cost that is nearly
the same as that of the uncorrected DT scheme. On the other hand, the LS method is shown to be
first-order accurate on arbitrary grids and second-order accurate on structured grids, while for the
default distance-based weighting scheme a particular non-integer exponent (−3/2) enlarges the set
of grid types for which the method is second-order accurate. The different order of accuracy of the
gradient schemes on structured versus unstructured grids is attributed to the fact that the former
become less skewed and more uniform as they are refined, whereas the latter usually do not.
First-order accurate gradients are, for most purposes, compatible with second-order accurate FVMs
because differentiation of a second-order accurate variable gives a first-order accurate derivative. In
Section 6 we proceed beyond the analysis of the gradient schemes themselves and test them as components of a FVM for the solution of a Poisson equation on various kinds of unstructured grids.
Experiments with both an in-house code and OpenFOAM show that the FVMs are zeroth-order ac-

2

N2

N1
N3

c01 m1
c1

P

c05

n1

c5

N4

Figure 1: Part of an unstructured grid, showing cell P and its neighbouring cells, each having a single
common face with P . The shaded area lies outside the grid. The faces and neighbours of cell P are numbered
in anticlockwise order, with face f separating P from its neighbour Nf . The geometric characteristics of its
face 1 which separates it from neighbouring cell N1 are displayed. The position vectors of the centroids of cells
P and Nf are denoted by the same characters but in boldface as P and Nf ; cf is the centroid of face f and
c0f is its closest point on the line segment connecting P and Nf ; mf is the midpoint between P and Nf ; nf
is the unit vector normal to face f , pointing outwards of P . The shown cell P also has a boundary face (face
5), with no neighbour on the other side; in this case c05 is the projection of P onto the boundary face.

curate when they employ the common DT gradient, whereas they are second-order accurate if they
employ instead the least-squares gradient, or the proposed iteratively corrected DT gradient, or an
alternative DT gradient scheme, achieving accuracy that is in fact only marginally worse than that
obtained on Cartesian grids.

2

Preliminary considerations

We will focus on two-dimensional problems, although the one- and three-dimensional cases are also
discussed where appropriate and the conclusions are roughly the same in all dimensions. Let x, y
denote the usual Cartesian coordinates, and let φ(x, y) be a function defined over a domain, whose
gradient ∇φ we wish to calculate. It is convenient to introduce a convention where subscripts beginning
with a dot (.) denote differentiation with respect to the ensuing variable(s), e.g. φ.x ≡ ∂φ/∂x, φ.xy ≡
∂ 2 φ/∂x∂y etc. If the variables x, y etc. are used as subscripts without a leading dot then they are used
simply as indices, without any differentiation implied. Therefore, we seek the gradient ∇φ = (φ.x , φ.y )
of the function φ. The domain over which the function is defined is discretised by a grid into a number
of non-overlapping finite volumes, or cells. A cell can be arbitrarily shaped, but its boundary must
consist of a number of straight faces, as in Fig. 1, each separating it from a single other cell or from the
exterior of the domain (the latter are called boundary faces). We assume a cell-centred finite volume
method, meaning that the values of φ are known only at the geometric centres of the cells and at the
centres of the boundary faces. The notation that is adopted in order to describe the geometry of the
grid is presented in Fig. 1. Also, we will denote vectors by boldface characters.
Our goal is to derive approximate algebraic gradient operators ∇a which return values ∇a φ(P ) ≈
∇φ(P ) at each cell centroid P , using information only from the immediate neighbouring cells and
a , φa ). The
boundary faces. The components of the approximate gradient are denoted as ∇a φ = (φ.x
.y
a
operators ∇ must be capable of approximating the derivative on grids of arbitrary geometry. To
study the effect of grid geometry on the accuracy, we must define some indicators of grid irregularity.
With the present grid arrangement, we find it useful to define three kinds of such grid irregularities,
which we shall call here “non-orthogonality”, “unevenness” and “skewness” (other possibilities exist,
see e.g. [23]). We will define these terms with the aid of Fig. 1. With the nomenclature defined in that
figure, we will say that face f of cell P exhibits non-orthogonality if Nf − P is not parallel to nf ; a
3

N2
c2 δ η

N1

δξ
c1

P
c3

c4

N3
N4

Figure 2: Part of a grid formed by equispaced parallel grid lines. See the text for details.

measure of non-orthogonality is the angle between the vectors (Nf −P ) and nf . Also, we will say that
face f exhibits unevenness if the midpoint of the line segment joining P and Nf , mf = (P + Nf )/2,
does not coincide with c0f (i.e. the cell centres are unequally spaced on either side of the face); a
measure of unevenness is kc0f − mf k/kNf − P k. Finally, we will say that face f exhibits skewness if
c0f does not coincide with cf (i.e. the line joining the cell centres does not pass through the centre of
the face); a measure of the skewness is kcf − c0f k/kNf − P k.
Before discussing the DT and LS gradient schemes, it is useful to examine a simpler scheme which
is applicable on very plain grids that are formed from two families of equispaced parallel straight
lines intersecting at a constant angle, as in Fig. 2. In this case all the cells of the grid are identical
parallelograms (Cartesian grids with constant spacing are a special case of this category). Figure
2 shows a cell P belonging to such a grid, and its four neighbouring cells. The vectors δ ξ and δ η
are parallel to the grid lines and span the size of the cells. Due to the grid properties it holds that
δ ξ = c1 − c3 = N1 − P = P − N3 and δ η = c2 − c4 = N2 − P = P − N4 . It can be assumed that two
variables, ξ and η, are distributed in the domain, such that in the direction of δ ξ the variable ξ varies
linearly while η is constant, and in the direction of δ η the variable η varies linearly while ξ is constant.
Then the grid can be considered to be constructed by drawing lines of constant ξ and of constant η,
equispaced (in the ξ, η space) by ∆ξ and ∆η, respectively. Let us assume also that the grid density
can be increased by adjusting the spacings ∆ξ and ∆η, but their ratio must be kept constant, e.g. if
∆ξ = h then ∆η = αh with α being a constant, independent of h. The variable h determines the grid
fineness. Therefore, the direction of the grid vectors δ ξ and δ η remains constant with grid refinement,
but their lengths are proportional to the grid parameter h.
This idealized grid exhibits no unevenness or skewness. It possibly exhibits non-orthogonality,
but this poses no problem as far as the gradient calculation is concerned. Since points N3 , P and
N1 are collinear and equidistant, the rate of change of any quantity φ in the direction of δ ξ can
be approximated at point P from the values at N3 and N1 using second-order accurate central
differencing. In the same manner, the rate of change in the direction of δ η at P can be approximated
from the values at N4 and N2 . So, let the grid vectors be written in Cartesian coordinates as
δ ξ = (δxξ , δyξ ) and δ η = (δxη , δyη ), respectively. Then, by expanding the function φ in a two-dimensional
Taylor series along the Cartesian directions, centred at point P , and using that to express the values
at the points N1 = P + δ ξ , N3 = P − δ ξ , N2 = P + δ η and N4 = P − δ η we get
φ(N1 ) − φ(N3 ) = ∇φ(P ) · (2δ ξ ) + O(h3 )
φ(N2 ) − φ(N4 ) = ∇φ(P ) · (2δ η ) + O(h3 )
This system can be solved for ∇φ(P ) ≡ (φ.x , φ.y ) to get
"
#
"
#"
#
φ.x (P )
δyη −δyξ
φ(N1 ) − φ(N3 )
1
=
·
2ΩP −δxη
φ.y (P )
φ(N2 ) − φ(N4 )
δxξ

+

1
2ΩP

"

δyη −δyξ
−δxη

#"
#
O(h3 )
·
O(h3 )
δxξ

where ΩP = δxξ δyη − δyξ δxη = kδ ξ × δ η k is the volume (area in 2D) of cell P . The last term in the above
equation, involving the unknown O(h3 ) terms, is of order O(h2 ) because ΩP = O(h2 ) and all the δ’s
4

(x2,1 , y2,1 )

ξ = ξ0 + ∆ξ
ξ = ξ0
(x1,1 , y1,1 )
η = η0 + ∆η

(x0,1 , y0,1 )
(x2,0 , y2,0 )
ξ = ξ0 + 2∆ξ
(x1,0 , y1,0 )

η = η0
(x0,0 , y0,0 )

Figure 3: Part of a grid (dashed straight lines) constructed from lines of constant ξ (red curves) and η
(blue curves), where ξ, η are variables distributed smoothly in the domain. The lines are equispaced with
constant spacings ∆ξ and ∆η, respectively. Grid node (i, j) is located at (ξi , ηj ) = (ξ0 + i ∆ξ, η0 + j ∆η) in the
computational space, where (ξ0 , η0 ) is a predefined point, and at (xi,j , yi,j ) in the physical space.

are O(h). So, carrying out the matrix multiplications we arrive at
" η
#
"
#
"
#
δy (φ(N1 )−φ(N3 )) − δyξ (φ(N2 )−φ(N4 ))
O(h2 )
φ.x (P )
1
+
=
∇φ(P ) ≡
2ΩP δxξ (φ(N2 )−φ(N4 )) − δxη (φ(N1 )−φ(N3 ))
O(h2 )
φ.y (P )
|
{z
}

(1)

∇s φ(P )

Finally, we drop the unknown O(h2 ) terms in Eq. (1) and are left with a second-order approximation
to the gradient, ∇s ; the dropped terms are called the truncation error of the operator ∇s . The fact
that this formula has been derived for a grid constructed from equidistant parallel lines may seem too
restrictive, but in fact the utility of Eq. (1) goes beyond this narrow context, as will now be explained.
Consider again structured grid generation based on a pair of variables ξ, η distributed smoothly in
the domain, where curves of constant ξ and η are drawn at equal intervals of ∆ξ and ∆η, respectively.
This time ξ and η are not required to vary linearly nor to be constant along straight lines. Therefore,
a curvilinear grid such as that shown in Fig. 3 may result, constructed by joining the points of
intersection of these two families of curves by straight line segments (the dashed lines in Fig. 3).
There is a one-to-one correspondence of coordinates (x, y) of the physical space to coordinates
(ξ, η) of the computational space. Therefore, not only are the computational coordinates functions of
the physical coordinates (ξ = ξ(x, y) and η = η(x, y)), but the physical coordinates (x, y) can also be
regarded as functions of the computational coordinates (x = x(ξ, η) and y = y(ξ, η)). Since the latter
vary smoothly in the domain, (x, y) can be expanded in Taylor series around a reference point (ξ0 , η0 ).
For example, for the x coordinate,
x(ξ0 + δξ, η0 + δη) = x(ξ0 , η0 ) + x.ξ δξ + x.η δη + O(h2 )

(2)

where the derivatives are evaluated at point (ξ0 , η0 ). In this way, we can express the coordinates of
all the grid vertices (xi,j , yi,j ) shown in Fig. 3 as functions of (x0,0 , y0,0 ) ≡ (x(ξ0 , η0 ), y(ξ0 , η0 )) and of
the derivatives x.ξ etc. there. Using these expansions, it is easy to show that




x1,1 − x1,0
y1,1 − y1,0
lim
=
lim
= 1
(3)
∆ξ,∆η→0 x0,1 − x0,0
∆ξ,∆η→0 y0,1 − y0,0
and


lim

∆ξ,∆η→0

x2,0 − x1,0
x1,0 − x0,0




=

lim

∆ξ,∆η→0

y2,0 − y1,0
y1,0 − y0,0


= 1

(4)

Equation (3) implies that, as the grid is refined by reducing the ∆ξ, ∆η spacings, neighbouring
grid lines of the same family become more and more parallel, so that grid skewness tends to zero.
5

Equation (4) implies that, as the grid is refined, neighbouring cells tend to become of equal size, so
that grid unevenness tends to zero. The conclusion is that on structured grids which are constructed
from smooth distributions of auxiliary variables (ξ, η), grid refinement1 causes the geometry of a cell
and its neighbours to locally approach that depicted in Fig. 2. This has an important consequence:
any numerical scheme for computing the gradient that reduces to Eq. (1) on parallelogram grids (such
as that of Fig. 2) is of second-order accuracy on smooth structured grids (such as that of Fig. 3)2 . It
so happens that both the DT scheme and the LS scheme belong to this category. Unfortunately, grid
refinement does not engender such a quality improvement when it comes to unstructured meshes.

3

Gradient calculation using the divergence theorem

For grids of more general geometry, like that of Fig. 1, a more general method is needed. Let ΩP
denote the volume of cell P and SP its bounding surface. The DT gradient scheme is based on a
derivative of the divergence theorem, which can be expressed as follows:
Z
Z
∇φ dΩ =
φ n ds
ΩP

SP

where n is the unit vector perpendicular to SP at each point, pointing outwards of the cell, while dΩ
and ds are infinitesimal elements of the volume and surface, respectively. The bounding surface SP
can be decomposed into F faces which are denoted by Sf , f = 1, . . . , F (F = 5 in Fig. 1). These faces
are assumed to be straight (planar, in 3 dimensions), as in Fig. 1, so that the normal unit vector n
has a constant value nf along each face f . Therefore, the above equation can be written as
!
Z
Z
F
X
∇φ dΩ =
nf
φ ds
(5)
ΩP

Sf

f =1

According to the midpoint integration rule [1,24], the mean value of a quantity over cell P (or face
f ) is equal to its value at the centroid P of the cell (or cf of the face), plus a second-order correction
term. Applying this to the mean values of ∇φ and φ over ΩP and Sf we get, respectively:
Z
Z
1
2
∇φ dΩ = ∇φ(P ) + O(h ) ⇒
∇φ dΩ = ∇φ(P ) ΩP + O(h4 )
(6)
ΩP ΩP
ΩP
1
Sf

Z

Z

2

φ ds = φ(cf ) + O(h ) ⇒
Sf

φ ds = φ(cf ) Sf + O(h3 )

(7)

Sf

where h is a characteristic grid spacing (we used that ΩP = O(h2 ) and Sf = O(h)). Substituting
these expressions into the divergence theorem, Eq. (5), we get
F
1 X
∇φ(P ) =
φ(cf ) Sf nf + O(h)
ΩP

(8)

f =1

The above formula is exact as long as the unknown O(h) term is retained, which consists mostly
of face contributions arising from Eq. (7), whereas the volume contribution from Eq. (6) is only
O(h2 ). Dropping this term would leave us with a first-order accurate approximation. However, this
approximation cannot be the final formula for the gradient because Eq. (8) contains φ(cf ), the φ
values at the face centres, whereas we need a formula that uses only the values at the cell centres.
The common practice is to approximate φ(cf ) by φ(c0f ), the exact values of φ at points c0f rather than
1

It is stressed that grid refinement must be performed in the computational space (ξ, η), by simultaneously reducing
the spacings ∆ξ and ∆η. Otherwise, if refinement is performed directly in the physical space (x, y), for example by
joining the centroid of each cell to the centroids of its faces thus splitting it into four child cells, then equations (3) – (4)
do not necessarily hold.
2
To be precise, this depends on the skewness decreasing fast enough with grid refinement – see Section 3. For
structured grids, using Taylor series such as Eq. (2), it can be shown that the skewness is O(h).

6

cf (see Fig. 1); these values also do not belong to the set of cell-centre values, but since c0f lies on the
line segment joining cell centres P and Nf , the value φ(c0f ) can in turn be approximated by linear
interpolation between φ(P ) and φ(Nf ), say φ̄(c0f ) (the overbar denotes linear interpolation):
φ̄(c0f ) ≡

kc0f − Nf k
kNf − P k

φ(P ) +

kc0f − P k
kNf − P k

φ(Nf ) = φ(c0f ) + O(h2 )

(9)

Linear interpolation is known to be second-order accurate, hence the O(h2 ) term in the above equation.
Thus, by using φ̄(c0f ) instead of φ(cf ) in the right hand side of Eq. (8) and dropping the unknown
O(h) term we obtain an approximation to the gradient which depends only on cell-centre values of φ:
F
1 X
φ̄(c0f ) Sf nf
∇ φ(P ) ≡
ΩP
d0

(10)

f =1

This is called the “divergence theorem” (DT) gradient. It applies in both two and three dimensions.
An important question is whether and how the replacement of φ(cf ) by φ̄(c0f ), that led from Eq.
(8) to the formula (10), has affected the accuracy. To answer this question we first need to deduce
how much φ̄(c0f ) differs from φ(cf ). This can be done by using the exact φ(c0f ) as an intermediate
value. We will consider structured grids and unstructured grids separately.
Structured grids
As discussed in the previous Section, the skewness of smooth structured grids diminishes with refinement, and in fact expressing the points involved in its definition as Taylor series of the form (2) it can
be shown that kcf − c0f k/kNf − P k = O(h), or cf − c0f = O(h2 ). Therefore, expanding φ(cf ) in a
Taylor series about point c0f gives φ(cf ) = φ(c0f ) + ∇φ(c0f ) · (cf − c0f ) + O(h2 ) = φ(c0f ) + O(h2 ). Next,
φ(c0f ) can be expressed in terms of φ̄(c0f ) according to Eq. (9). Putting everything together results in
the sought relationship:
φ(cf ) = φ(c0f ) + O(h2 ) = (φ̄(c0f ) + O(h2 )) + O(h2 ) = φ̄(c0f ) + O(h2 )

(11)

Substituting this into Eq. (8), and considering that ΩP = O(h2 ) and Sf = O(h), we arrive at
∇φ(P ) =

F
1 X
φ̄(c0f ) Sf nf + O(h)
ΩP

(12)

f =1

The first term on its right-hand side is just the approximate gradient ∇d0 , Eq. (10). Therefore, the
truncation error of ∇d0 is τ ≡ ∇d0 φ(P ) − ∇φ(P ) = O(h). Actually, the situation is even better
because the leading terms of the contributions of opposite faces to the O(h) term in Eq. (12) cancel
out leaving a net O(h2 ) truncation error, making the method second-order accurate. The proof is
tedious and involves writing analytic expressions for the truncation error contributions of each face;
it can be performed more easily for a grid such as that shown in Fig. 2 where the divergence theorem
procedure described in the present Section is an alternative path to arrive at the exact same formula,
Eq. (1), ∇d0 ≡ ∇s , with its O(h2 ) truncation error.
The error cancellation between opposite faces does not always occur; sometimes the worst-case
scenario of O(h) truncation error predicted by Eq. (12) holds. An example is at boundary cells. Fig. 4
shows a boundary cell P belonging to a Cartesian grid which exhibits neither skewness nor unevenness
(in fact, it does not even exhibit non-orthogonality). Yet cell P has no neighbour on the boundary
side, and so the centre of its boundary face, c3 , is used instead of a neighbouring cell centre. This
introduces “unevenness” in the x- (horizontal) direction because the distances kN1 − P k = h and
kc3 − P k = h/2 are not equal. The x- component of the divergence theorem gradient (10) reduces to
d0
φ.x
(P ) =

1
(φ(N1 ) + φ(P ) − 2φ(c3 ))
2h
7

(13)

N2

c3

N1

P

N4

h

h/2

Figure 4: A boundary cell P belonging to a Cartesian grid.

which is only first-order accurate, since expanding φ(N1 ) and φ(c3 ) in Taylor series about P gives
d0
φ.x
(P ) = φ.x (P ) +

1
φ.xx (P ) h + O(h2 )
8

This offers a nice demonstration of the effect of error cancellation between opposite faces. Equation
(10) is obtained from Eq. (8) by using interpolated values φ̄(c0f ) instead of the exact but unknown
values φ(cf ) at the face centres. Therefore, one might expect that since the cell of Fig. 4 has a
boundary face and the exact value at its centre, φ(c3 ), is used rather than an interpolated value, the
result would be more accurate; but the above analysis shows exactly the opposite: the error increases
from O(h2 ) to O(h). This is due to the fact that by dropping the interpolation error on the boundary
face the corresponding error on the opposite face 1 is no longer counterbalanced.
Another example where error cancellation does not occur and the formal O(h) accuracy predicted
by Eq. (12) holds is the common case of grid that consists of triangles that come from dividing each
cell of a smooth structured grid along the same diagonal. For example, application of this procedure
to the grid shown in Fig. 2 results in the grid of Fig. 5. The latter may be seen to exhibit neither
unevenness nor skewness as the face centres cf coincide with the midpoints of the line segments joining
cell centre P to its neighbours Nf . Thus Eq. (11) holds, leading to Eq. (12). This time, however,
a tedious but straightforward calculation where neighbouring φ(Nf ) values are expressed in Taylor
series about P and substituted in (12) shows that there is no cancellation and the truncation error
remains O(h). If the grid comes from triangulation of a curvilinear structured grid (Fig. 3) then the
skewness is not zero but it diminishes with refinement at an O(h) rate (Eq. (11)) and exactly the same
conclusions hold.

N1

N2

P

c1
c3

c2

N3

Figure 5: Grid of triangles constructed through bisection of the cells of the structured grid of Fig. 2 along
their same diagonal.

In fact, from the considerations leading to Eq. (11) it follows that if a grid generation algorithm
is such that skewness diminishes as O(hp ) then the order of accuracy of the DT gradient is at least
min{p, 1}.
8

Unstructured grids
Unstructured grids usually consist of triangles (or tetrahedra in 3D), or of polygonal (polyhedral in
3D) cells which are also formed by a triangulation process. They are typically constructed using
algorithms that are based on geometrical principles that do not depend on grid fineness, so that there
is some self-similarity between coarse and fine grids and refinement does not reduce the skewness, i.e.
kcf − c0f k/kNf − P k = O(1). This means that cf − c0f = O(h) and therefore instead of Eq. (11) we
now have
φ(cf ) = φ(c0f ) + O(h) = (φ̄(c0f ) + O(h2 )) + O(h) = φ̄(c0f ) + O(h)
(14)
Substituting this into the exact equation (8) we get
F

1 X
∇φ(P ) =
φ̄(c0f ) + O(h) Sf nf + O(h)
ΩP
f =1

=

F
F
1 X
1 X
φ̄(c0f ) Sf nf +
O(h) Sf nf + O(h)
ΩP
ΩP
f =1

f =1

which, considering that ΩP = O(h2 ) and Sf = O(h), means that
∇φ(P ) =

F
1 X
φ̄(c0f ) Sf nf + O(1) + O(h)
ΩP

(15)

f =1

Again, the first term on its right-hand side is the approximate gradient ∇d0 , Eq. (10). But this time,
unlike in the structured grid case, the truncation error of ∇d0 is τ ≡ ∇d0 φ(P ) − ∇φ(P ) = O(1).
Nor do the leading terms of the face contributions to the truncation error cancel out to increase its
order, because on unstructured grids the orientations and sizes of the faces are unrelated. This is
an unfortunate result, because it means that the approximation (10) is zeroth-order accurate, as the
error O(1) does not decrease with grid refinement, but instead ∇d0 φ(P ) converges to a value that is
not equal to ∇φ(P ) (see Section 5.4 for an example).
Acknowledging that the lack of accuracy is due to the bad representation of the φ(cf ) values in
Eq. (8) by φ̄(c0f ), application of the formula (10) is often followed by a “corrector step” where, instead
of the values φ̄(c0f ), the “improved” values φ̂(cf ) are used, defined as
φ̂(cf ) ≡ φ̄(c0f ) + ∇d0 φ(c0f )·(cf − c0f )

(16)

where ∇d0 φ(c0f ) is obtained by linear interpolation (Eq. (9)) between ∇d0 φ(P ) and ∇d0 φ(Nf ) at
point c0 (these were calculated in the previous step (10), the “predictor” step). This results in an
approximation which is hopefully more accurate than (10):
∇d1 φ(P ) ≡

F
1 X
φ̂(cf ) Sf nf
ΩP

(17)

f =1

The values φ̂(cf ) are expected to be better approximations to φ(cf ) than φ̄(c0f ) are, because Eq. (16)
tries to account for skewness by mimicking the Taylor series expansion
φ(cf ) = φ(c0f ) + ∇φ(c0f ) · (cf − c0f ) + O(h2 )

(18)

Unfortunately, since in Eq. (16) only a crude approximation of ∇φ(c0f ) is used, namely ∇d0 φ(c0f ) =
∇φ(c0f ) + O(1), what we get by subtracting Eq. (16) from Eq. (18) is φ̂(cf ) = φ(cf ) + O(h) which
may have greater accuracy but not greater order of accuracy than the previous estimate φ̄(c0f ) =
φ(cf ) + O(h), Eq. (14). Substituting this into Eq. (8) we arrive again at an equation similar to (15)
which shows that the error of the approximation (17) is also of order O(1).
9

Further correction steps may be applied in the same manner; a fixed finite number of such steps
may increase the accuracy, but the order of accuracy with respect to grid refinement will remain
zero. But if this procedure is repeated until convergence to an operator ∇d∞ , say, then ∇d∞ would
simultaneously satisfy both Eqs. (16) and (17), or combined in a single equation:
∇

d∞

F
F
1 X d∞
1 X
0
0
φ(P ) −
∇ φ(cf )·(cf − cf ) Sf nf =
φ̄(c0f ) Sf nf
ΩP
ΩP
f =1

(19)

f =1

where we have moved all terms involving the gradient to the left-hand side. An analogous equation
can be derived for the exact gradient, from Eqs. (8), (18), and (9):
∇φ(P ) −

F
F
1 X
1 X
∇φ(c0f )·(cf − c0f ) Sf nf =
φ̄(c0f ) Sf nf + O(h)
ΩP
ΩP
f =1

(20)

f =1

By subtracting Eq. (20) from Eq. (19) we get an expression for the truncation error τ ≡ ∇d∞ φ − ∇φ:
F
1 X
τ (c0f )·(cf − c0f ) Sf nf = O(h)
τ (P ) −
ΩP

(21)

f =1

The left-hand side is a linear combination of not only τ (P ) but also τ (Nf ) at all neighbour points.
Since ΩP = O(h2 ), Sf = O(h) and cf − c0f = O(h), the coefficients of this linear combination are O(1)
i.e. they do not depend on the grid fineness but only on the grid geometry (skewness, unevenness etc.).
Expression (21) suggests that τ = O(h), i.e. that ∇d∞ is first-order accurate, although to ascertain
this one would have to assemble Eqs. (21) for all grid cells in a large linear system Aτ = b ⇒ τ = A−1 b
(where τ stores the x- and y- (and z-, in 3D) components of τ at all cell centres and b = O(h)), select
a suitable matrix norm k · k, and ensure that kA−1 k remains bounded as h → 0.
So, iterating Eqs. (16) and (17) until convergence one would hope that a first-order accurate gradient ∇d∞ would be obtained. As reported in [19], convergence of this procedure is not guaranteed
and underrelaxation may be needed, leading to a large number of required iterations and great computational cost. In practice, the desired accuracy will be reached with a finite number of iterations,
but this number is not known a priori and must be determined by trials; furthermore, in order to
maintain the property that grid refinement improves the accuracy, the number of iterations should
increase with grid refinement in order to avoid accuracy stagnation. Considering that the main rival
of the DT gradient, the LS gradient, costs approximately as much as the DT gradient with a single
corrector step (see Sec. 6.1), it can be seen that this iterative procedure is impractical and may end
up consuming most of the computational time of a FVM solver. Alternatively, one could directly
solve all the equations (19) simultaneously in a large linear system; this is also proposed and tested
in [19], along with a second-order accurate method that additionally computes the Hessian matrix.
Obviously, this approach also involves a very large computational cost.
However, we would like to propose here a way to avoid the extra cost of the iterative procedure for
the gradient. The FVM generally employs outer iterations to solve a PDE (e.g. SIMPLE iterations,
as opposed to the inner iterations of the linear solvers employed). Outer iterations are necessary when
the equations solved are non-linear, but may be used also when solving linear problems if some terms
are treated with the deferred-correction approach. FVMs that use gradient schemes such as those
presently discussed are almost always iterative, with the gradients computed using the values of the
dependent variable from the previous outer iteration. The idea is then to exploit these outer iterations,
dividing the “gradient iterations” among them: at each outer iteration n a single “gradient iteration”
is performed according to
F

1 X  n−1 0
∇ φ(P ) =
φ
(cf ) + ∇d(n−1) φ(c0f )·(cf − c0f ) Sf nf
ΩP
dn

(22)

f =1

which means that one has to store not only the solution of the previous outer iteration φn−1 but also
the gradient of that iteration ∇d(n−1) φ, which many codes do already. Since only one calculation is
10

performed per outer iteration, the cost of this procedure is almost as low as that of the uncorrected
DT gradient (10), but it converges to the first-order accurate ∇d∞ gradient operator instead of the
zeroth-order accurate ∇d0 . For explicit time-dependent FVM methods, Eq. (22) could be applied
with n − 1 denoting the previous time step; in the very first time step it may be necessary to perform
several “gradient iterations” to obtain a consistent gradient to begin with.
Iterative schemes that are different from (22) can also be devised; a similar iterative technique
is used in [25] in order to determine the coefficients of reconstruction polynomials without solving
large linear systems. There, Jacobi, Gauss-Seidel, and SOR-type iterations are applied, with the
latter found to be the most efficient. The present scheme (22) is reminiscent of the Jacobi iterative
procedure (although not exactly equivalent, since ∇φ(P ) appears also at the right-hand side, within
the averaged gradient, evaluated from the previous iteration). A scheme reminiscent of the GaussSeidel iterative procedure would be one where in the right-hand side of Eq. (22) new values of the
gradient at neighbour cells (calculated at the current outer iteration) are used whenever available
instead of using the old values. However, in the present work only the scheme (22) is tested in Section
6.1.
At this point we would like to make a brief comment concerning the calculation of the gradient at
boundary cells. The above analysis has assumed that at boundary faces the values of φ at the face
centres are available. However, this is not always the case; in situations where these values are not
available, they must be approximated to at least second-order accuracy as otherwise the DT gradient
will be inconsistent, as the preceding analysis has shown. For example, in problems with Neumann
boundary conditions the directional derivative normal to the boundary, g say, is given rather than the
boundary values. If face 5 of Fig. 1 belongs to such a boundary, then expressing φ(P ) in a Taylor
series about point c05 gives
φ(c05 ) = φ(P ) + gkc05 − P k + O(h2 )

(23)

(where g is measured in the direction pointing out of the domain). Equation (23) provides a secondorder accurate approximation for φ(c05 ), but only a first-order accurate approximation for φ(c5 ). Therefore, if we just use the value φ̄(c05 ) ≡ φ(P ) + gkc05 − P k in the DT gradient formula then we will
get a zeroth-order accurate gradient. Note that this will hold even for structured grids if the grid
lines intersect the boundary at an angle (c05 6= c5 ). It is not difficult to show that a second order
approximation is φ(c5 ) ≈ φ̄(c05 ) + ∇φ(P ) · (c5 − c05 ). This introduces ∇φ(P ) also in the right-hand
side of the gradient expression (8), which poses no problem for the iterative procedure (22).
Finally, we note that alternative schemes can be derived from the Gauss divergence theorem that
are consistent on unstructured grids and are worth mentioning, although the present work focuses on
the standard method. For example, consider the application of the divergence theorem method not to
the actual cell P but to the auxiliary cell marked by dashed lines in Fig. 6. The endpoints of each face
of this cell are either cell centroids or boundary face centroids, and so the value of φ can be computed
at any point on the face to second-order accuracy, by linear interpolation between the two endpoints.
Therefore, in Eq. (8) the values φ(cf ) (cf being the face centroids of the auxiliary cell, marked by
empty circles in Fig. 6) are calculated to second-order accuracy rather than first-order, leading to
first-order accuracy of the computed gradient. Second-order accuracy is inhibited also by the fact that
the midpoint integration rule (Eq. (6)) requires that the gradient be computed at the centroid of the
auxiliary cell, which now does not, in general, coincide with point P (the situation changes if P and
the centroid of the auxiliary cell tend to coincide with grid refinement). This method, along with a
more complex variant, is mentioned in [7]; it is also tested in Section 6.1. Yet another method would
be to use cell P itself rather than an auxiliary cell, but, similarly to the previous method, to calculate
the values φ(cf ) in Eq. (8) by linear interpolation from the values at the face endpoints (vertices)
rather than from the values at the cell centroids straddling the face. An extra step must therefore
precede where φ is approximated at the cell vertices to second-order accuracy from its values at the cell
centroids. This adds to the computational cost and furthermore requires of the grid data structures
to contain lists relating each vertex to its surrounding cells. More information on this method and
further references can be found in [11, 26].

11

N2

N1
N3

P

c5
N4

Figure 6: The divergence theorem method can approximate ∇φ(P ) to first-order accuracy on arbitrarily
irregular grids if applied to the auxiliary cell bounded by dashed lines, rather than on the actual cell P . This
cell is formed by joining the centroids of the neighbouring cells and of the boundary faces of cell P by straight
line segments. The open circles denote the centroids of these segments.

4

Gradient calculation using least squares minimisation

The starting point of the “least-squares” (LS) method for calculating the gradient of a quantity φ at
the centroid of a cell P is the expression of the values of φ at the centroids of all neighbouring cells as
Taylor series expansions about the centroid of P . For convenience, we decompose the position vectors
in Cartesian coordinates as P = (x0 , y0 ), Nf = (xf , yf ) and denote ∆xf ≡ xf − x0 , ∆yf ≡ yf − y0 .
Then, for each neighbour f it follows from the Taylor expansion that
φ(Nf ) − φ(P ) = ∇φ(P ) · (Nf − P ) + ςf

(24)

where

1
1
φ.xx (∆xf )2 + φ.xy ∆xf ∆yf + φ.yy (∆yf )2 + O(h3 )
(25)
2
2
is an associated truncation error (the second derivatives are evaluated at P ). The basic idea of the
method is to drop the unknown ςf terms and solve the remaining linear system by least squares since,
in general, the number of equations (F = number of neighbouring cells) will be greater than the
number of unknowns (two, φ.x (P ) and φ.y (P )). Of course some of the “neighbours” may be boundary
face points such as the centroid c5 in Fig. 1, or, in the case of Neumann boundary conditions, the
projection point c05 . In the latter case φ(c05 ) can be approximated by Eq. (23). Furthermore, unlike for
the DT method, it is now easy to incorporate additional, more distant cells as neighbours, a practice
that may prove advantageous in some cases3 (see e.g. [26]), although in the present work we will
restrict ourselves to using only the cells that share a face with P .
It is advantageous to first multiply each equation f by a suitably selected weight wf (weighted
least squares). Then the system of all these equations can be written in matrix form as




 

 
w1
∆φ1
w1
∆x1 ∆y1 "
ς1
#











w
∆x
∆y
∆φ
w
2
2
2
2
2  φ.x (P )


 

 ς2 
· .  = 
· .
·
+  .  (26)




.
.
.
..
..
..  φ.y (P )

  .. 

  ..
 .. 
wF
∆φF
wF
∆xF ∆yF
ςF
|
|
| {z }
{z
} | {z }
{z
} |
{z
} | {z }
ςf =

0

0

0

0

W

b

W

A

z

ς

where ∆φf ≡ φ(Nf ) − φ(P ). The alternative matrix notation W b = W (Az + ς) is also introduced
above, to facilitate the discussions that follow. The above equations are exact, with the vector ς
3
It may even be necessary in some cases. For example, some of the triangular cells at the corners of the grid in Fig.
20(b) have only one neighbour cell; if the values of φ are not available at the boundary faces (e.g. if φ is the pressure in
incompressible flows) then at least one more distant cell must be used, because the LS calculation requires at least two
neighbour points.

12

ensuring that the system (26) has a unique solution, which is independent of the weights, despite
having more equations than unknowns. The redundant equations are just linear combinations of the
non-redundant ones. But once ς is dropped, the system will, in general, have no solution.
Systems with no solution can be solved in the “least squares” sense [27,28]. In matrix notation, if a
linear system Az = b with more equations than unknowns cannot be solved because the vector b does
not lie in the column space of the matrix A then the best that can be done is to find the vector z that
minimises the error b − Az. This error will be minimised when its projection onto the column space of
A is zero, i.e. when b − Az is perpendicular to each column of A, or AT (b − Az) = 0 ⇒ AT Az = AT b
(where T denotes the transpose). This latter system (called the normal equations) has a unique
solution provided that A has independent columns. The solution z is called the least squares solution
because it minimises the L2 norm kb − Azk, i.e. the sum of the squares of the individual components
of the error e ≡ b − Az
P(note that each individual error component ei is the error of the corresponding
equation i, ei = bi − j aij zj ).
In the weighted case W Az = W b the matrix A is replaced by the product W A and b by W b, so
that the corresponding solution to the normal equations is
z = (AT W T W A)−1 AT W T W b

(27)

Now the quantity minimised is the norm of the error of the weighted system, W (b − Az) = W e. Thus,
if equation i is assigned a larger weight than equation j then the method will prefer to make ei small
at the expense of ej .
For now, we proceed with the least squares methodology on Eq. (26) without discarding ς in order
to assess its impact on the error. We therefore first left-multiply the system (26) by (W A)T , and then
solve it to obtain
z = (AT W T W A)−1 AT W T W b + (AT W T W A)−1 AT W T W ς

(28)

The matrix AT W T W A is invertible provided that W A has independent columns, which requires that
A have independent columns because W is diagonal. The two columns of A will be independent if
there are also two independent rows, i.e. vectors Nf − P , so we need at least two neighbours Nf to
lie at different directions with respect to P , which is normally the case (in three dimensions A has
three columns and three independent neighbour directions are needed). Then, substituting for A, W ,
z, b and ς from Eq. (26) into Eq. (28) and performing the matrix multiplications and inversions, we
arrive at
z ≡ ∇φ(P ) = ∇ls φ(P ) + τ
(29)
where ∇ls φ(P ) is the first term on the right side of Eq. (28),
 F
  F

F
P
P
P
2
2
2
2
−
(∆yf ) wf
∆xf ∆yf wf  
∆xf ∆φf wf 
1 
f =1
 f =1
 f =1

ls
∇ φ(P ) ≡
 F
· F

F
P
P
D P



2
2
2
2
−
∆xf ∆yf wf
(∆xf ) wf
∆yf ∆φf wf
f =1

f =1

|

(30)

f =1

{z

} |

M

{z
βb

}

with (1/D)M = (AT W T W A)−1 and βb = AT W T W b, and τ is the second term:
1
τ =
M·
D

"

F
P
f =1

F
P

∆xf ςf wf2

f =1

#T
∆yf ςf wf2

=

1
M βς
D

(31)

with βς = AT W T W ς. In the 2D case, D = |M |, the determinant of M . Equation (29) gives the exact
gradient ∇φ(P ), but since the truncation error τ is unknown we drop it and use the expression (30)
alone as the approximate “least squares” gradient. In the 3D case the explicit formula is more involved
than Eq. (30), but again the least squares gradient is obtained by solving Eq. (27). If, just for the
purposes of discussing the 3D case, we denote the Cartesian components of the P
displacement vectors
1
2
3
T
T
as Nf − P = (∆xf , ∆xf , ∆xf ), then the (i, j) entry of the matrix A W W A is f ∆xif ∆xjf wf2 (it is
13

P
a symmetric matrix) and the i-th component of the vector AT W T W b is f ∆xif ∆φf wf2 . The normal
equations (27) are sometimes ill-conditioned, and in such cases it is better to solve the least-squares
system by QR factorisation of A [27].
The above derivation has produced the explicit expression (31) for the truncation error τ . To
analyse it further, we can assume that all the weights share the same dependency on the grid spacing,
namely wf = O(hq ) for some real number q (independent of f ), as is the usual practice. Then the
factors of Eq. (31) have the following magnitudes: Since ∆xf and ∆yf are O(h) the coefficients of
M have magnitude O(h2+2q ). Consequently, M being 2 × 2, its determinant will have magnitude
(O(h2+2q ))2 ⇒ 1/D = O(h−2(2+2q) ). Finally, considering that ςf = O(h2 ), the components of βς are
of O(h2q+3 ). Multiplying all these together, Eq. (31) shows that τ = O(h), independently of q. This
is not surprising, given that the approximation is based on Eq. (24) which assumes a linear variation
of φ in the vicinity of point P . So, ∇ls is at least first order accurate, even on grids of arbitrary
geometry.
The order of accuracy of ∇ls may be higher than one if some cancellation occurs between the
components of τ for certain grid configurations, similarly to the DT operator. In particular, a tedious
but straightforward calculation shows that when applied to a parallelogram grid such as that shown
in Fig. 2, ∇ls again reduces to the second-order accurate formula (1), provided only that the weights
of parallel faces are equal: w1 = w3 and w2 = w4 . This will hold due to symmetry if the weights are
dependent only on the grid geometry. Of course, as for the DT gradient, this has the consequence
that the LS gradient is second-order accurate also on smooth curvilinear structured grids.
The choice of weights
According to the preceding analysis, the least squares method is first-order accurate on grids of
arbitrary geometry and second-order accurate on smooth structured grids. This holds irrespective of
the choice of weights, and in fact it holds even for the unweighted method (wf = 1). The question
then arises of whether a suitable choice of weights can offer some advantage.
The weights commonly used are of the form wf = (∆rf )−q where ∆rf = kNf − P k is the distance
between the two cell centres and q is an integer, usually chosen as q = 1 [11, 12, 16] or q = 2 [6, 14].
The unweighted method amounts to q = 0. As noted, the least squares method finds the approximate
gradient that minimises kW (b − Az)k, which for the various choices of q amounts to minimising:
q=0:


2
X
∆φf
ls
∆rf
− ∇ φ(P ) · df
∆rf

(32)

f

q=1:

X  ∆φf

2

− ∇ φ(P ) · df

(33)

2
X  1  ∆φf
ls
− ∇ φ(P ) · df
∆rf ∆rf

(34)

f

q=2:

ls

∆rf

f

where df = (Nf − P )/∆rf is the unit vector in the direction from P to Nf , and so ∇ls φ · df is the
least squares directional derivative ∂φ/∂rf in that direction. In the q = 1 case, expression (33) shows
that the least squares procedure shows no preference in trying to set the directional derivative in each
neighbour direction f equal to the finite difference ∆φf /∆rf . On the other hand, in the unweighted
method (32) the discrepancies between the directional derivatives and the finite differences are weighted
by the distances ∆rf so that the method prefers to reduce ∆φ/∆rf − ∇ls φ · df along the directions of
the distant neighbours at the expense of the directions of the closer neighbours. The exact opposite
holds for the q = 2 case (34) where the discrepancies are weighted by 1/∆rf so that the result is
determined mostly by the close neighbours. Intuitively, this latter choice seems more reasonable as
the linearity of the variation of φ is lost as one moves away from P and thus at distant neighbours
the finite differences ∆φ/∆rf are less accurate approximations of the directional derivatives. Thus it
is not surprising that usually the q > 0 methods outperform the unweighted method.

14

These are the common weight choices, but it so happens that the particular non-integer exponent
q = 3/2 confers enhanced accuracy compared to the other q choices under special but not uncommon
circumstances. This fact does not appear to be widely known in the literature; we have seen it only
briefly mentioned in [18, 29]. So, consider the vector βς in the expression (31) for the error, and
substitute for ς from Eq. (25). Then the first component of βτ becomes
X

∆xf ςf wf2 =

f

X
φ.yy X
φ.xx X
(∆xf )3 wf2 + φ.xy
(∆xf )2 ∆yf wf2 +
∆xf (∆yf )2 wf2
2
2
f

+

X

f

f

O(h4 ) wf2

(35)

f

Now, if two neighbours, i and j say, lie at opposite directions to point P but at the same distance then
wi = wj , ∆xi = −∆xj and ∆yi = −∆yj so that their contributions in each of the first three sums in
the right hand side of Eq. (35) cancel out. If all neighbour points are arranged
in such pairs, like in
P
4
the grid of Fig. 2, then these three sums become zero leaving only the
O(h )wf2 term. The same
4−2q
holds for the second component of βς , so that βς = O(h
) overall because wf = O(h−q ). Then Eq.
(31) gives τ = (1/D)M βς = O(h2 ) i.e. the method is second-order accurate for any exponent q.
The particular choice wf = (∆rf )−3/2 amounts to dropping the wf ’s from the first three sums of
the right-hand side of Eq. (35) and replacing therein every instance of ∆xf by ∆xf /∆rf and every
instance of ∆yf by ∆yf /∆rf . These ratios are precisely cos θf and sin θf , respectively, where θf is
the angle that the direction vector Nf − P makes with the horizontal direction. If two neighbours,
i and j, lie at opposite directions then θi = θj + π so that cos θi = − cos θj and sin θi = − sin θj ,
and their contributions cancel out in the aforementioned three sums, irrespective of whether these two
neighbours lie at equal distances to P or not. The same holds for the second component of βς . Thus,
if all neighbour points are arranged in such pairs then again what remains of the right-hand side of
Eq. (35) is only the last term and so τ = O(h2 ). For example, the LS gradient with q = 3/2 is
second-order accurate at the boundary volume P of Fig. 4, whereas it is only first-order accurate with
any other choice of q. Furthermore, the same result will hold if the neighbours are not arranged in
pairs at opposite directions but tend to become so with grid refinement. This is the case with smooth
structured grids, as shown in Section 2, and therefore the LS gradient with q = 3/2 is second-order
accurate at boundary cells of all smooth structured grids.
Another property of the q = 3/2 LS gradient is that it is second order accurate if all neighbour
points are arranged at equal angles. Unfortunately, this property holds only with more than three
neighbour points, which limits its usefulness. A proof is provided in Appendix A.
Finally, a question that arises naturally is whether full 2nd-order accuracy can be achieved on
arbitrary grids by allowing non-diagonal entries in the weights matrix. It turns out that this is indeed
possible, but yields a method that is equivalent to the least squares solution of a system of Taylor
expansions (24) with terms higher than first-order included. The procedure is sketched in Appendix B.
We do not advocate it, because direct solution of the system of higher-order Taylor expansions would
have the added advantage of solving also for the second derivatives – see e.g. [30]. An alternative
second-order accurate method is described in [31]. The second-order accurate methods are much more
expensive than the present method.

5
5.1

Numerical tests on the accuracy of the gradient schemes
One-dimensional tests

The methods are first tested on a one-dimensional problem so as to examine the effect of unevenness,
isolated from skewness. The derivative of the single-variable function φ(x) = tanh x is calculated at
101 equispaced points spanning the x ∈ [0, 2] interval. The
against the exact
P results are compared
2
ls
solution φ.x = 1 − (tanh x) and the mean absolute error i |φ.x (x(i) ) − φ.x (x(i) )|/101 is recorded for
each method. In order to introduce unevenness, the neighbours of point xi are not chosen from this set
of equispaced points but are set at xi,f = xi + ∆xf /2r where the ∆xf belong to a predetermined set
15

(a) {∆x0f } = {0.05, 0.10}

(b) {∆x0f } = {−0.10, 0.05}

(c) {∆x0f } = {−0.10, 0.05, 0.15}

(d) {∆x0f } = {−0.20, −0.10, 0.05, 0.15}

Figure 7: Mean errors of the calculation of the derivative of φ(x) = tanh(x) with various methods – see the
text for details.

of displacements {∆xf }Ff=1 (common to all xi points) and the integer r is the level of grid refinement.
For example, if the chosen set of displacements is {−0.05, 0.1} (F = 2 neighbours), then the derivative
at point xi will be calculated using the values of φ at the three points {xi − 0.05/2r , xi , xi + 0.1/2r }.
The order of accuracy of each method is determined by incrementing the level of refinement r.
In order to test the methods thoroughly, several sets of initial displacements {∆x0f }Ff=1 were used.
In Fig. 7, each diagram corresponds to a different such set and the mean error is plotted as a function
of the number of displacement halvings, r. The slope of each curve reveals the order of accuracy of
the corresponding method. The methods tested are: (a) the DT method, denoted “d”, (b) the LS
methods with weight exponents q = 0, 1, 1.5, 2 and 3, indicated on each curve, (c) the non-diagonal
weights LS method of Appendix B denoted “ND”, and (d) a simpler variant of the DT method which
is sometimes used [20], denoted “da”, where the values at face centres are calculated by arithmetic
averaging, φ(c0f ) = (φ(P ) + φ(N ))/2, instead of linear interpolation (9). The DT methods are only
applicable in the case plotted in Fig. 7(b) because exactly two neighbours are required, one on each
side of xi . To apply the DT methods xi is regarded as the centroid of a cell of size equal to the
minimum distance between xi and any of its neighbours.
A result not shown is that when the displacements stencil is symmetric (e.g. {∆x0f } = {−0.1, 0.1})
all methods produce identical, second-order accurate results. This is consistent with the two-dimensional
16

case where on symmetric grids like that of Fig. 2 all methods reduce to the formula (1). When the
stencil is not symmetric, Fig. 7 shows that the LS methods with a diagonal weights matrix become
first-order accurate, except for the q = 3/2 method which retains second-order accuracy when there
are equal numbers of neighbours on either side (Figs. 7(b) and 7(d)). The unweighted method (q = 0)
is always the least accurate; the optimum accuracy is achieved with 1 ≤ q ≤ 2, while a further increase
of q is unprofitable (q = 3). The method of Appendix B is always second-order accurate. Concerning
the DT methods (Fig. 7(b)), the method of Sec. 3, indicated by “d” in the figure, gives identical results
with the q = 1 least squares method. On the other hand, the simplified method, indicated by “da”, is
zeroth-order accurate, which agrees with the findings reported in [20].

5.2

Uniform Cartesian grids

Next, the two-dimensional methods are used to calculate the gradient of the function φ(x, y) = tanh(x)·
tanh(y) on the unit square (x, y) ∈ [0, 1] × [0, 1] using uniform Cartesian grids of different fineness.
The exact gradient is φ.x = (1 − (tanh x)2 ) tanh y and φ.y = (1 − (tanh y)2 ) tanh x. All grid cells
are geometrically identical squares of side h = 0.25/2r where r is the level of refinement; however,
boundary cells are topologically different from interior cells because they have one or more boundary
faces where the function value at the face centre has to be used (Fig. 4). It so happens that the
function tanh has zero second derivative at the boundaries x = 0 and y = 0, which may artificially
increase the order of accuracy of the methods there. However, the general behaviour of the methods
at boundary cells can be observed at the x = 1 and y = 1 boundaries where no such special behaviour
of the tanh function applies.
On each grid the gradient of φ is calculated at all cell centres using the DT method ∇d0 (Eq.
(10)), and the LS methods ∇ls (Eq. (30)) with weight exponents q = 0, 1, 1.5 and 2. Since there is
no skewness (c0 = c in Eq. (16)), the application of corrector steps is meaningless. The methods are
evaluated by comparing the mean and maximum truncation errors, defined as
τmean ≡

Mr
1 X
k∇a φ(P ) − ∇φ(P )k
Mr

(36)

P =1

Mr

τmax ≡ max k∇a φ(P ) − ∇φ(P )k
P =1

(37)

where k · k denotes the L2 norm of a vector in a single cell, Mr is the number of cells of grid r, and ∇a
is any of the gradient schemes considered. These errors are plotted in Figs. 8(a) and 8(b), respectively.
The theory predicts that all methods should be second-order accurate at interior cells, i.e. τ =
O(h2 ), because they all reduce to formula (1) there. At boundary cells (Fig. 4) the favourable conditions that are responsible for second-order accuracy are lost, and the methods should revert to
first-order accuracy except for the LS method with q = 3/2 of which second-order accuracy is still
expected. Figure 8(b) confirms that the maximum error, which occurs at some boundary cell, is
τmax = O(h) for all methods, except the q = 3/2 method for which τmax = O(h2 ). Of the other
methods, the unweighted LS method (q = 0) is the least accurate and the q = 2 method is slightly
better than the q = 1 method. The DT method gives identical results with the q = 1 LS method, as
they both revert to Eq. (13) at boundary cells.
Figure 8(a) shows the mean errors. Since all methods give the same results at interior cells any
differences are due to the different errors at the boundary cells. The slope of the curves is the same,
and corresponds to τmean = O(h2 ). This does not contradict with the O(h) errors at the boundary
cells: the number of such cells along each side of the domain equals 1/h, so that in total there are
O(1/h) boundary cells, each contributing an O(h) error. Their total contribution to τmean in Eq. (36)
is therefore O(1/h)·O(h)/Mr = O(h2 ) since Mr = 1/h2 . The interior cell contribution is O(h2 ) also,
since τ = O(h2 ) at each individual interior cell.

5.3

Smooth curvilinear grids

Next, we try the methods on smooth curvilinear grids. The same function φ = tanh(x) · tanh(y) is
differentiated, but the domain boundaries now have the shapes of two horizontal and two vertical
17

(a) mean error, τmean

(b) maximum error, τmax

Figure 8: The mean (a) and maximum (b) errors (defined by Eqs. (36) and (37), respectively) of the various
methods for calculating the gradient of the function φ = tanh(x) tanh(y) on uniform Cartesian grids. The
abscissa r designates the grid; grid r has a uniform spacing of h = 0.25/2r . The exponent q used for the LS
weights is shown on each curve. The DT method produces identical results with the q = 1 LS method.

sinusoidal waves (Fig. 9), beginning and ending at the points (0, 0), (1, 0), (1, 1) and (0, 1). The grid
is generated using a very basic elliptic grid generation method [32]. In particular, smoothly varying
functions ξ(x, y) and η(x, y) are assumed in the domain, and the grid consists of lines of constant ξ
and of constant η. The left, right, bottom and top boundaries correspond to ξ = 0, ξ = 1, η = 0
and η = 1, respectively. In the interior of the domain ξ and η are assumed to vary according to the
following Laplace equations:
ξ.xx + ξ.yy = 0
η.xx + η.yy = 0
These equations guarantee that ξ and η vary smoothly in the domain, but their solution ξ =
ξ(x, y), η = η(x, y) is not much help in constructing the grid. Instead, we need the inverse functions
x = x(ξ, η), y = y(ξ, η) which explicitly set the locations of all grid nodes; node (i, j) is located at
(xi,j , yi,j ) ≡ (x(ξ = i ∆ξ, η = j ∆η), y(ξ = i ∆ξ, η = j ∆η)). With ξ, η ∈ [0, 1], the constant spacings
∆ξ and ∆η are adjusted according to the desired grid fineness. Therefore, using the chain rule of
partial differentiation it can be shown that the above equations can be expressed in inverse form as
g22 x.ξξ − 2g12 x.ξη + g11 x.ηη = 0
g22 y.ξξ − 2g12 y.ξη + g11 y.ηη = 0
where
g11 = x.ξ2 + y.ξ2
g22 = x.η2 + y.η2
g12 = x.ξ x.η + y.ξ y.η
In order to cluster the points near the boundaries in the physical domain, we accompany the above
equations with the following boundary conditions: at the bottom boundary we set x = 0.5+0.5 sin(π(ξ−
0.5)) and y = sin(2πx), and at the left boundary we set y = 0.5+0.5 sin(π(η−0.5)) and x = − sin(2πy).
At the top and right boundaries we set the same conditions, respectively, adding 1 to x at the right
boundary and 1 to y at the top boundary. Better results can be obtained by using a more elaborate
method such as described in [33], but this suffices for the present purposes.
18

(a) r = 0

(b) r = 1

(c) r = 2

(d) r = 3

Figure 9: The first four of the series of grids constructed on a domain with sinusoidal boundaries via elliptic
grid generation. See the text for more details.

Now (x, y) have become the dependent variables while (ξ, η) are the independent variables, which
acquire values in the unit square [0, 1] × [0, 1]. The grid equations were solved numerically with a finite
difference method on a 513 × 513 point uniform Cartesian grid. Note that the dependent variables
(x, y) are stored at the grid nodes, i.e. at the intersection points of the grid lines, instead of at the cell
centres. The derivatives are approximated by second-order accurate central differences; for example,
at point (i, j), x.ξ ≈ (xi+1,j −xi−1,j )/2h, x.η ≈ (xi,j+1 −xi,j−1 )/2h, x.ξξ ≈ (xi+1,j −2xi,j +xi−1,j )/h2 etc.
where h = 1/512 is the grid spacing. The resulting system of nonlinear algebraic equations was solved
using a Gauss-Seidel iterative method where in the equations of the (i, j) node all terms are treated
as known from their current values except for xi,j and yi,j which are solved for. The convergence of
the method was accelerated using a minimal polynomial extrapolation technique [34], and iterations
were carried out until machine precision was reached.
After obtaining x(ξ, η) and y(ξ, η), a series of successively refined grids in the physical domain were
constructed by drawing lines of constant ξ and lines of constant η at intervals ∆ξ = ∆η = 0.25/2r , for
r = 0, 1, . . . 7. The first four grids are shown in Fig. 9. The gradient calculation methods were then
applied on each of these grids. The errors of each method are depicted in Fig. 10.
These grids exhibit all kinds of grid irregularity, but the unevenness and the skewness diminish with
grid refinement, as explained in Section 2. In particular, Table 1 shows that the measures of both these
grid qualities have magnitudes of O(h). These measures were defined in Sec. 2: kc0f − mf k/kNf − P k
for unevenness and kcf − c0f k/kNf − P k for skewness (see Fig 1 for definitions). The values listed in
Table 1 are the average values among all faces of each grid, excluding boundary faces. Therefore, it is
19

(a) mean error, τmean

(b) maximum error, τmax

Figure 10: The mean (a) and maximum (b) errors (defined by Eqs. (36) and (37), respectively) of the gradient
schemes applied on the function φ = tanh(x) tanh(y) on smooth curvilinear grids (Fig. 9). The abscissa r
designates the grid; r = 0 is the coarsest grid (Fig. 9(a)), and grid r comes from subdividing every cell of grid
r − 1 into 4 child cells in the computational space (see text). The blue solid lines correspond to the LS methods
with q = 0, 1 and 2 as indicated on each curve; the blue dash-dot line corresponds to the LS method with
q = 3/2; and the red dashed lines correspond to the DT methods dc where c is the number of corrector steps.

expected that eventually the methods will behave as in the Cartesian case. Indeed, Fig. 10(b) shows
that as the grid is refined τmax tends to decrease at a first-order rate for all methods, except for the
least squares method with q = 3/2, for which it decreases at a second-order rate. Accordingly, Fig.
10(a) shows that τmean tends to decrease at a second-order rate for all methods, with the second-order
rate attained earlier by the q = 3/2 method. Therefore, like on the Cartesian grids, the methods are
second-order accurate at interior cells but revert to first-order accuracy at boundary cells, except for
the q = 3/2 method.
Of the LS methods the least accurate is the unweighted method, followed by the q = 2 method.
The undisputed champion is the q = 3/2 method because, as mentioned, it retains its second-order
accuracy even at boundary cells. Concerning the DT methods, the method with no corrections (Eq.
10) performs similarly to the unweighted LS method. Application of a corrector step (Eq. (17)) now
does make a difference, since skewness is present at any finite grid density, bringing the accuracy of
the method on a par with the best weighted LS methods, except of course the q = 3/2 method. We
also tried a second corrector step but it did not bring any noticeable improvement.
Next, on the same grids we also applied discrete gradient operators to calculate the gradient of a
linear and of a quadratic function. The schemes tested are the DT gradient operators ∇d0 , ∇d1 and
∇d∞ (the latter approximated with 100 correction steps), and the LS gradient operators with q = 0,
1 and 3/2, denoted as ∇ls0 , ∇ls1 and ∇ls3/2 , respectively, in Tables 1 and 2. “Exactness” is sometimes
used as an aid to either determine the order of accuracy of a method or to design a method to achieve
a desired order of accuracy (e.g. Appendix B): a first-order accurate gradient scheme would normally
be exact for linear functions, and a second-order accurate gradient scheme would normally be exact
for quadratic functions. However, the results listed in Tables 1 and 2 show that in the present case the
DT gradient is not exact even for linear functions while the LS gradient is exact for linear functions
but not for quadratic functions, despite both methods being second-order accurate.
In particular, Table 1 shows that the DT gradient without corrector steps (∇d0 ) is not exact for the
linear function (the errors are not zero) but converges to the exact gradient at a rate that approaches
second-order as the grid is refined. Performing a corrector step (∇d1 ) brings a significant improvement
in accuracy, with an observed convergence rate order of between 2 and 3, but still the operator is not
exact. This inexactness is anticipated since the grids are skewed and the DT scheme cannot cope with

20

skewness. However, grid refinement causes skewness to diminish and the DT accuracy to improve at
a second-order rate. In the limit of many corrector steps (∇d∞ ) the operator becomes exact, with the
errors at machine precision levels even at the coarsest grid. All of the LS schemes are also exact for
the linear function.
Concerning the quadratic function (Table 2), none of the schemes is exact but they are all secondorder accurate, with the q = 3/2 LS scheme being the most accurate, and the q = 0 LS and zerocorrection DT schemes being the least accurate. A single corrector step in the DT scheme (∇d1 ) brings
the maximum attainable improvement, since the error levels of ∇d1 and ∇d∞ are nearly identical.
The second-order convergence rates are due to the improvement of grid quality with refinement, as
explained in Sections 3 and 4.
Table 1: Mean errors (Eq. (36)) of various schemes for calculating the gradient ∇φ = (1, 2) of the linear
function φ(x, y) = x + 2y + 0.5 on the series of grids shown in Fig. 9. Also displayed are the measures of grid
skewness and unevenness (defined in Sec. 2), averaged over all faces of each grid excluding boundary faces.

Grid r

Skew.

Unev.

∇d0

0

1.33·10−1

1.10·10−1

4.16·10−1

1

−2

5.68·10

8.00·10

−2

−1

2

3.03·10−2

4.83·10−2

6.93·10−2

3

1.64·10−2

2.70·10−2

4

8.44·10−3

5

−3

4.29·10

7.41·10

6

2.16·10−3

3.78·10−3

6.64·10−4

7

1.09·10−3

1.91·10−3

1.86·10−4

1.86·10

∇d1

∇d∞

∇ls0

∇ls1

∇ls3/2

8.00·10−2

2.19·10−15

6.04·10−16

7.39·10−16

9.37·10−16

−2

−15

−15

−15

1.68·10−15

4.46·10

1.34·10

1.45·10

4.80·10−3

9.53·10−15

2.43·10−15

2.41·10−15

2.81·10−15

2.41·10−2

8.91·10−4

1.93·10−14

4.42·10−15

4.52·10−15

4.84·10−15

1.43·10−2

7.72·10−3

1.16·10−4

3.97·10−14

8.74·10−15

8.77·10−15

9.07·10−15

−3

−3

−5

−14

−14

−14

1.74·10−14

2.32·10

1.90·10

1.43·10

7.98·10

1.70·10

1.71·10

1.94·10−6

1.60·10−13

3.36·10−14

3.38·10−14

3.40·10−14

3.25·10−7

3.20·10−13

6.69·10−14

6.70·10−14

6.72·10−14

Table 2: Mean errors (Eq. (36)) of various schemes for calculating the gradient ∇φ = (2x + 2y, 2x − 2y) of the
quadratic function φ(x, y) = x2 + 2xy − y 2 on the series of grids shown in Fig. 9.

5.4

Grid r

∇d0

∇d1

∇d∞

∇ls0

∇ls1

∇ls3/2

0

4.95·10−1

2.28·10−1

2.29·10−1

3.54·10−1

1.94·10−1

1.47·10−1

1

2.20·10−1

1.27·10−1

1.29·10−1

1.72·10−1

8.16·10−2

5.81·10−2

2

−2

8.38·10

4.48·10

−2

−2

−2

8.09·10

3.13·10

−2

2.24·10−2

3

2.84·10−2

1.46·10−2

1.46·10−2

3.20·10−2

1.12·10−2

7.43·10−3

4

8.99·10−3

4.63·10−3

4.64·10−3

1.14·10−2

3.76·10−3

2.05·10−3

5

2.66·10−3

1.40·10−3

1.40·10−3

3.75·10−3

1.19·10−3

5.28·10−4

6

−4

7.45·10

4.00·10

−4

−4

−3

1.15·10

3.50·10

−4

1.35·10−4

7

2.07·10−4

1.10·10−4

3.29·10−4

9.74·10−5

3.41·10−5

4.47·10

4.01·10

1.10·10−4

Grids of localised high distortion

Structured grids that are constructed not by solving partial differential equations, as in Section 5.3,
but by algebraic methods may lack the property that unevenness and skewness diminish with grid
refinement. This is especially true if the domain boundaries include sharp corners at points other than
grid line endpoints. For example, the grid of Fig. 11 is structured, consisting of piecewise straight
lines. At the line joining the sharp corners, the intersecting grid lines change direction abruptly. This
causes significant skewness which is unaffected by grid refinement.

21

Figure 11: A structured grid where the grid lines belonging to one family change direction abruptly at the
dashed line joining the pair of sharp corners, where grid skewness and not diminishing with grid refinement.

(a) r = 0

(b) r = 1

(c) r = 2

Figure 12: A series of multi-level, or composite, grids. Each grid r comes from the previous grid r − 1 by
evenly subdividing each cell into four child cells.

A similar situation may occur when adaptive mesh refinement is used, depending on the treatment
of the interaction between levels. Figure 12 shows multi-level grids, which consist of regions of different
fineness. Such grids are often called composite grids [35]. One possible strategy is to treat the cells at
the level interfaces as topologically polygonal [8, 29, 36]. For example, cell P of Fig. 13(a) has 6 faces,
each separating it from a single other cell. Its face f1 separates it from cell N1 which belongs to the
finer level. Faces such as f1 , which lie on grid level interfaces, exhibit non-orthogonality, unevenness,
and skewness. If the grid density is increased throughout the domain, as in the series of grids shown
in Fig. 12, then these interface distortions remain insensitive to the grid fineness, like for the marked
line in Fig. 11. Alternative schemes exist which avoid changing the topology of the cells by inserting
a layer of transitional cells between the coarse and the fine part of the grid (e.g. [37, 38]) but they also
lead to high, non-diminishing grid distortions at the interface.
We computed the gradient of the same function φ(x, y) = tanh(x) tanh(y) on a series of composite
grids the first three of which are shown in Fig. 12. Figure 14 shows how τmean and τmax vary with grid
refinement. This time, τmean is defined a little differently than Eq. (36) to account for the different grid
levels: the error of each individual cell is weighted by the cell’s volume (i.e. the area, in the present
two-dimensional setting):
Mr
1 X
τmean ≡
ΩP k∇a φ(P ) − ∇φ(P )k
(38)
Ω
P =1

where Ω is the total volume of the domain and ΩP is the volume of cell P .
We can identify three classes of cells that are topologically different. Apart from the familiar
classes of interior and boundary cells, there is now also the class of cells that touch the level interfaces,
which shall be called interface cells (these belong to two sub-classes, coarse- and fine-level cells, as
in Figs. 13(a) and 13(b), respectively). Interface cells possess high skewness and unevenness that do
not diminish with grid refinement. The behaviour of the gradient-calculation methods at interior and

22

N4

N3

f4

f3
N2
f2

N5

N2
θ2 = ϑ

P

c2

f5
f1

N3

θ1 = −ϑ
N1

c3

P c
1
c3

c01

N1

f6
N4

N6
(a) coarse cell

(b) fine cell

Figure 13: Topological and geometrical characteristics of a coarse cell (a) and of a fine cell (b) adjacent to a
level interface, in a composite grid.

boundary cells has already been tested in Sections 5.2 and 5.3, so our interest now focuses on interface
cells. Skewness, the most detrimental grid distortion, is encountered only at those and therefore it is
there that the maximum errors (Fig. 14(b)) occur.
In Fig. 14(b) we observe that at the level interfaces all the LS methods converge to the correct
solution at a first-order rate, τmax being lowest for the q = 1 method, followed closely by the q = 3/2
and q = 2 methods. On the other hand, none of the DT methods converge to the correct solution
there, although nearly an order of magnitude accuracy improvement is obtained with each corrector
step. We can determine the operator to whom ∇d0 converges as follows. For the interface cell P of
Fig. 13(b), formula (10) amounts to the following series of approximations:
φ.x (P ) ≈
≈

φ(c1 ) − φ(c3 )
φ(c01 ) − φ(c3 )
≈
h
h
[α φ(N1 ) + (1 − α) φ(P )] − [0.5 φ(P ) + 0.5 φ(N3 )]
d0
≡ φ.x
(P )
h

(39)

where h is the length of the side of cell P . In the last step, the values of φ at points c01 and c3 were
approximated with linear interpolation between points N1 and P , and P and N3 , respectively; α is
an interpolation factor which equals α = 0.3 for the present geometry. We then substitute in (39)
φ(N1 ) and φ(N3 ) with their two-dimensional Taylor series about P , considering that if P = (x0 , y0 )
then N1 = (x0+3h/2, y0−h/2) and N3 = (x0−h, y0 ) – see Fig. 13(b). The following result is obtained:
d0
φ.x
(P ) =

3α + 1
α
φ.x (P ) −
φ.y (P ) + O(h)
2
2

d0 converges not to φ but to an operator that involves both φ and φ .
Therefore, as h → 0, φ.x
.x
.x
.y
Next we examine the mean error in Fig. 14(a). The plot can be interpreted by considering separately the error contributions of each class of cells. The contributions of interior and boundary cells
to τmean are both O(h2 ), as discussed in Section 5.2 (for the q = 3/2 method the boundary cell
contribution is O(h3 )).
At interface cells the LS methods behave similarly as on boundary cells, because they produce
O(h) errors there as well due to unevenness and skewness. The total length of the level interfaces is
constant, O(1). The number of interface cells is O(h−1 ) because it equals this constant length divided
by the cell size which is O(h). Their contribution to the mean error in Eq. (38) is (number of cells) ×
(volume of one cell) × (error at a cell) = O(h−1 ) × O(h2 ) × O(h) = O(h2 ). This is confirmed by Fig.
14(a), where all the LS methods converge to the exact solution at a second-order rate.

23

(a) mean error, τmean

(b) maximum error, τmax

Figure 14: The mean (a) and maximum (b) errors (defined by Eqs. (38) and (37), respectively) of the gradient
calculation methods applied to the function φ = tanh(x) tanh(y) on locally refined grids (Fig. 12). The abscissa
r designates the grid; r = 0 is the coarsest grid (Fig. 12(a)), and grid r comes from subdividing every cell of
grid r − 1 into 4 identical child cells. The blue solid lines correspond to the LS methods with weight exponents
q = 0, 1 and 2, which are indicated on each curve; the blue dash-dot line corresponds to the LS method with
q = 3/2; and the red dashed lines correspond to the DT methods dc where c is the number of corrector steps.

On the other hand, for the DT methods the contribution of interface cells to τmean is (number of
cells) × (volume of one cell) × (error at a cell) = O(h−1 ) × O(h2 ) × O(1) = O(h). Figure 14(a) shows
that for the d0 method (no corrector steps) this O(h) component is so large that it dominates τmean
even at coarse grids. For the d1 method (one corrector step), at coarse grids this O(h) component
is initially small compared to the bulk O(h2 ) component that comes from all the other cells, so that
τmean appears to decrease at a second-order rate up to a refinement level of r = 3; but eventually it
becomes dominant and beyond r = 5 the d1 curve is parallel to the d0 curve, with a first-order slope.
With two corrector steps, the O(h) component is so small that up to r = 7 it is completely masked
by the O(h2 ) component and the method appears to be second-order accurate. More grid refinements
are necessary to reveal its asymptotic first-order accuracy.
An observation that raises some concern in Fig. 14(a) is that the DT methods outperform the LS
methods at those grids where they have not yet degraded to first order. This suggests that there may
be some room for improvement in the latter. A potential source of the problem is suggested by Fig.
13(a). Along the horizontal direction, cell P has one cell on its left side (N5 ) and two cells on its right
side (N1 and N2 ). Since points N1 and N2 are quite close to each other there is some overlap in the
information they convey. Yet the weights of the LS method depend only on the distance of Nf from
P , while any clustering of the Nf points in some direction is not taken into account. Thus, points N1
and N2 , being closer to P than N5 may individually contribute equally (q = 1) or more (q > 1) to the
calculation of the gradient at P than point N5 does. Combined they contribute much more. So, the
horizontal component of the gradient is calculated using mostly information from the right of cell P ,
whereas information from its left is undervalued. The DT methods do not suffer from this deficiency
because they weigh the contribution of each point by the area of the respective face; faces f1 and f2
are half in size than f5 , and so points N1 and N2 together contribute to the gradient approximately
as much as N5 alone does.
In order to seek a remedy, we investigate the contributions of points N1 , N2 and N5 (Fig. 13(a))
to the vector βς of the error expression (31). As for Eq. (35), we substitute for ςf from Eq. (25) into
the expression for βς , but then proceed to substitute ∆xf = ∆rf cos θf , ∆yf = ∆rf sin θf . For the
particular points under consideration we have, with reference to Fig. 13(a), θ1 = −ϑ, θ2 = ϑ and
θ5 = π, while ∆r2 = ∆r1 . Finally, concerning the weights, since cells P and N5 both belong to the

24

(a) mean error, τmean

(b) maximum error, τmax

Figure 15: As for Fig. 14, but with the weights (41) applied at level interfaces in the case of the LS methods.

same grid level we choose not to tamper with the weight w5 and simply use the possibly second-order
accurate q = 3/2 scheme: w5 = (∆r5 )−3/2 . Due to symmetry we set w2 = w1 and our goal is to
suitably select w1 to achieve a small error. Putting everything together, the joint contribution of these
three points to βς is, neglecting higher order terms,




φ.xx (cos ϑ)3 (∆r1 )3 w12 − 12 + φ.yy (sin ϑ)2 cos ϑ (∆r1 )3 w12


(40)
2
3
2
2φ.xy cos ϑ (sin ϑ) (∆r1 ) w1
Since the values of the higher order derivatives can be arbitrary, it is obvious that the above contribution does not become zero for any choice of w1 ; thus, second-order
√ accuracy cannot be achieved.
The best that can be done is to choose w1 = (cos ϑ)−3/2 (∆r1 )−3/2 / 2 so that the term in parentheses
multiplied by φ.xx becomes zero. In fact, since this choice is not guaranteed to minimise the error
and since cos ϑ ≈ 1 for relatively small ϑ, we chose to drop the cos ϑ factor and applied the following
scheme to all least squares methods:
(
kNf − P k−q
face f does not touch a finer level
wf =
(41)
1
√ kNf − P k−q
face f touches a finer level
2
(In the three-dimensional
case, cell P of Fig. 13(a)
would have four fine-level neighbours on its right
√
√
and the 1/ 2 factor in (41) would become 1/ 4 = 1/2). Using this scheme, we obtained the results
shown in Fig. 15, which can be seen to be better than those of Fig. 14, especially for the q = 3/2 and
q = 2 methods which now rival the d2 method with respect to the mean error.
The modification (41) is easy and inexpensive. More elaborate methods could be devised, applicable to more general grid configurations, in order to properly deal with the issue of point clustering at certain angles. For example, a point Nf could be assigned to span a “sector” of angle
∆θf = |θf +1 − θf −1 |/2 (assuming that the points are numbered in either clockwise or anti-clockwise
order); the sum of all these sectors would then equal 2π, and they could be incorporated into the
weights such that points with smaller sectors would have less influence over the solution. Alternatively, the face areas could be incorporated into the weights, as in the DT method [18].

5.5

Grids with arbitrary distortion

As mentioned in Section 3, general-purpose unstructured grid generation methods result in unevenness
and skewness that are insensitive to grid fineness throughout the domain, not just in isolated regions
25

(a) r = 0

(b) r = 1

(c) r = 2

Figure 16: A series of excessively distorted grids. Grid r is constructed by random perturbation of the nodes
of Cartesian grid r + 1 of Section 5.2.

as for the grids of Section 5.4. Therefore, in the present Section the methods are tested under these
conditions; grids of non-diminishing distortion were generated by randomly perturbing the vertices
of a Cartesian grid. Using such a process we constructed a series of grids, the first three of which
are shown in Fig. 16. The perturbation procedure is applied as follows: Suppose a Cartesian grid
with grid spacing h. If node (i, j) has coordinates (xij , yij ) then the perturbation procedure moves
0 ) = (x + δ x , y + δ y ) where δ x and δ y are random numbers in the interval
it to a location (x0ij , yij
ij
ij
ij
ij
ij
ij
[−0.25h, 0.25h). Because all perturbations are smaller than h/4 in both x and y, it is ensured that all
grid cells remain simple convex quadrilaterals after all vertices have been perturbed. Grids based on
triangles as well as three-dimensional cases will be considered in Section 6.
The gradient of the same function φ = tanh(x) tanh(y) is calculated, and the errors are plotted in
Fig. 17. This time, all cells belong to a common category. Concerning the mean error, Fig. 17(a) is in
full agreement with the theory. The LS methods converge to the exact gradient at a first-order rate,
with the unweighted method being, as usual, the least accurate, while the differences between the
weighted methods are very slight. On the other hand, the DT methods, as expected, do not converge
to the exact gradient. Performing corrector steps improves things, but in every case the convergence
eventually stagnates at some grid fineness.
As far as the maximum errors are concerned, Fig. 17(b) shows that those of the LS methods
decrease at a first order rate, with the unweighted method being the least accurate. On the other
hand, the maximum errors of the DT methods actually increase with grid refinement. Presumably this
is due to the fact that as the number of grid nodes is increased the probability of encountering higher
degrees of skewness somewhere in the domain increases. Performing corrector steps reduces the error,
but it is interesting to note that grid refinement causes a somewhat larger error increase when more
corrector steps are performed. The deterioration of the ∇d0 method with grid refinement propagates
across the iterative correction procedure and eventually it is expected that the errors produced by
the ∇dc operator used as a predictor will become so large that it will provide less (rather than more)
accurate face centre values to the resulting “corrected” operator ∇d(c+1) , making the latter worse than
∇dc itself.
Similarly to Sec. 5.3 we also tested the gradient schemes on a linear function. The results are listed
in Table 3, together with grid quality metrics. Table 3 confirms that in the present case grid skewness
and unevenness are roughly independent of grid refinement, i.e. their measures are O(1), whereas they
were O(h) in the structured grid case (Table 1). As a result, the DT gradient schemes with a finite
number of corrector steps (∇d0 and ∇d1 in Table 3) now do not converge to the exact gradient. In
the limit of infinite corrector steps the DT gradient (∇d∞ ) becomes exact for linear functions, as are
all the LS gradient schemes (∇ls0 and ∇ls1 in Table 3).

26

(a) mean error, τmean

(b) maximum error, τmax

Figure 17: The mean (a) and maximum (b) errors (defined by Eqs. (36) and (37), respectively) of the gradient
calculation methods applied on the function φ = tanh(x) tanh(y) on the series of globally distorted grids (Fig.
16). The blue solid lines correspond to the least squares methods with weight exponents q = 0, 1 and 2, which
are indicated on each curve; the blue dash-dot line corresponds to the least squares method with q = 3/2; and
the red dashed lines correspond to the divergence theorem methods dc where c is the number of corrector steps.

6

Use of the gradient schemes within finite volume PDE solvers

So far we have examined the gradient schemes per se, examining their truncation error through mathematical tools and numerical experiments. Although there are examples of independent use of a
gradient scheme such as in post-processing, the application of main interest is within finite volume
methods (FVMs) for the solution of partial differential equations (PDEs). The gradient scheme is
but a single component of the FVM and how it affects the overall accuracy depends also on the PDE
solved as well as on the rest of the FVM discretisation. In the present section we provide some general
comments and some simple demonstrations. The focus is on the effect of the gradient scheme on
unstructured grids, where the DT scheme was shown to be inconsistent; on such grids the tests of Sec.
5.5 showed that the weights exponent of the LS method plays a minor role and so we only test the
q = 1 LS variant.
We begin with the observation that the approximation formula (10) is very similar to the formulae
used by FVMs for integrating convective terms of transport equations over grid cells. Therefore,
Table 3: Mean errors (Eq. (36)) of various schemes for calculating the gradient ∇φ = (1, 2) of the linear
function φ(x, y) = x + 2y + 0.5 on the series of grids shown in Fig. 16. Also displayed are the measures of grid
skewness and unevenness (defined in Sec. 2), averaged over all faces of each grid excluding boundary faces.

Grid r

Skew.

Unev.

∇d0

∇d1

∇d∞

∇ls0

∇ls1

0

4.79·10−2

4.94·10−2

1.99·10−1

1.55·10−2

3.83·10−15

9.84·10−16

1.05·10−15

1

5.46·10−2

5.32·10−2

2.43·10−1

1.96·10−2

7.29·10−15

1.55·10−15

1.59·10−15

2

4.97·10−2

5.14·10−2

2.17·10−1

1.61·10−2

1.48·10−14

3.07·10−15

3.15·10−15

3

5.16·10

−2

−2

−1

−2

−14

−15

6.18·10−15

4

5.22·10−2

5
6

5.16·10

2.29·10

1.81·10

5.17·10−2

2.32·10−1

1.89·10−2

5.87·10−14

1.22·10−14

1.23·10−14

5.18·10−2

5.17·10−2

2.32·10−1

1.86·10−2

1.18·10−13

2.45·10−14

2.45·10−14

5.18·10−2

5.14·10−2

2.32·10−1

1.87·10−2

2.34·10−13

4.86·10−14

4.86·10−14

27

2.94·10

6.11·10

according to the same reasoning as in Section 3, such formulae also imply truncation errors of order
O(1) on arbitrary grids; this is true even if an interpolation scheme other than (9) is used to calculate
φ̄(c0f ), or even if the exact values φ(c0f ) are known and used. The order of the truncation error can
be increased to O(h) by accounting for skewness through a correction such as (16), provided that the
gradient used is at least first-order accurate. These observations may raise concern about the overall
accuracy of the FVM; however, it is known that the order of reduction of the discretisation error with
grid refinement is often greater than that of the truncation error. Thus, O(1) truncation errors do
not necessarily imply O(1) discretisation errors; the latter can be of order O(h) or even O(h2 ). This
phenomenon has been observed by several authors (including the present ones [36]) and a literature
review can be found in [39]. The question then naturally arises of whether and how the accuracy of the
gradient scheme would affect the overall accuracy of the FVM that uses it. A general answer to this
question does not yet exist, and each combination of PDE / discretisation scheme / grid type must be
examined separately. For simple cases such as one-dimensional ones the problem can be tackled using
theoretical tools but on general unstructured grids this has not yet been achieved [39]. Thus we will
resort to some simple numerical experiments that amount to solving a Poisson equation.

6.1

Tests with an in-house solver

We solve the following Poisson equation:
∇ · (−k∇φ) = b(x, y)

on Ω = [0, 1]×[0, 1]

(42)

φ = c(x, y)

on SΩ

(43)

with k = 1, where SΩ is the boundary of the domain Ω (the unit square), and


b(x, y) = 2 tanh(x) tanh(y) 2 − (tanh(x))2 − (tanh(y))2
c(x, y) = tanh(x) tanh(y)

(44)
(45)

This is a heat conduction equation with a heat source term and Dirichlet boundary conditions. The
source term was chosen such that the exact solution of the above equation is precisely φ = c(x, y). The
domain was discretised by a series of progressively finer grids of 32 × 32, 64 × 64, . . . , 512 × 512 cells
denoted as grids 0, 1, . . . 4, respectively. The grids were generated by the same randomised distortion
procedure as those of Fig. 16, only that their boundaries are straight (Fig. 20(a)). Corresponding
undistorted Cartesian grids were also used for comparison. According to the FV methodology, we
integrate Eq. (42) over each cell, apply the divergence theorem and use the midpoint integration rule
to obtain for each cell P a discrete equation of the following form:
F
X

Df = b(P ) ΩP

(46)

f =1

where b(P ) is the value of the source term at the centre of cell P and
Z
Df =
−k∇φ · nf dS ≈ −Sf k ∇φ(cf ) · nf

(47)

Sf

is the diffusive flux through face f of the cell. We will test here two alternative discretisations of the
fluxes Df , both of which utilise some discrete gradient operator.
The first is the “over-relaxed” scheme [8, 40] which, according to [41], is very popular and is the
method of choice in commercial and public-domain codes such as FLUENT, STAR-CD, STAR-CCM+
and OpenFOAM. This scheme splits the normal unit vector nf in (47) into two components, one in
the direction Nf − P and one tangent to the face. The unit vectors along these two directions are
denoted as df and tf , respectively (Fig. 18(a)). Thus, if nf = αdf + βtf then by taking the dot
product of this expression with nf itself and noting that tf and nf are perpendicular (nf · tf = 0) we

28

t1
c01
c1
P

N1

d1

N1
P0

c01

n1

c1

P5

P

n1

N10

c5
n5

(a)

(b)

Figure 18: (a) Adopted notation for the “over-relaxed” diffusion scheme (Eq. (48)) and for the scheme for
boundary faces (Eq. (50)); (b) Notation for the custom scheme (49). Face 1 of cell P is an inner face, while its
face 5 is a boundary face.

obtain α = (df · nf )−1 . We can therefore split nf = d∗f + t∗f where d∗f = df /(df · nf ) and t∗f = βtf
calculated most easily as t∗f = nf − d∗f . Then the flux (47) is approximated as
d∗f

}|

{
1
· df )
− Sf k ∇φ(c0f ) · t∗f
(nf · df )
z

Df ≈

−Sf k (∇φ(c0f )

≈ −Sf k

φ(Nf ) − φ(P )
1
− Sf k ∇a φ(c0f ) · t∗f
kNf − P k (df · nf )

= −Sf k

φ(Nf ) − φ(P )
− Sf k ∇a φ(c0f ) · t∗f
(Nf − P ) · nf

(48)

where ∇a is a discretised gradient calculated at the cell centres (either a DT gradient ∇d or the LS
gradient ∇ls ) whose role we wish to investigate, and ∇a φ(c0f ) is its value calculated at point c0f by
linear interpolation (9). This scheme obviously deviates somewhat from the midpoint integration rule,
substituting ∇φ(c0f ) · nf instead of ∇φ(cf ) · nf in (47), i.e. it does not account for skewness.
We also try an alternative scheme which is the standard scheme in our code [42] and is a slight
modification of a scheme proposed in [1]:
φ(Nf0 ) − φ(P 0 )

Df ≈ −Sf k

kNf0 − P 0 k

(49)

where
φ(P 0 ) ≈ φ(P ) + ∇a φ(P ) · (P 0 − P )
φ(Nf0 ) ≈ φ(Nf ) + ∇a φ(Nf ) · (Nf0 − Nf )
and points P 0 and Nf0 (Fig. 18(b)) are such that the line segment joining these two points has length
kNf − P k, is perpendicular to face f , and its midpoint is cf . Thus this scheme tries to account for
skewness. We will refer to this as the “custom” scheme.
Finally, irrespective of whether scheme (48) or (49) is used for the inner face fluxes, the boundary
face fluxes are always calculated as
Df ≈ −Sf k

φ(cf ) − φ(Pf )
kcf − Pf k

29

(50)

where f is a boundary face, point Pf is the projection of P on the line which is perpendicular to face
f and passes through cf (see Fig. 18(a)), and
φ(Pf ) ≈ φ(P ) + ∇a φ(P ) · (Pf − P )
Grid non-orthogonality activates in all schemes, (48), (49) and (50), a component that involves the
approximate gradient ∇a (it is activated in (49) also by unevenness or skewness). It is not hard to show
that this component contributes O(1) to the truncation error if ∇a is first-order accurate, and O(1/h)
if it is zeroth-order accurate. On the undistorted Cartesian grids that were used for comparison the
∇a terms of the schemes are not activated; in fact, schemes (48) and (49) reduce to the same simple
formula there. The reason why we do not simply substitute ∇a φ(cf ) for ∇φ(cf ) in Eq. (47) instead
of using schemes such as (48) and (49) is that this would allow spurious high-frequency (cell-to-cell)
oscillations in the numerical solution φ. This occurs because the ∇a calculation is mostly based on
differences of φ values that are two cells apart, and so the gradient field ∇a φ can be smooth even
though φ oscillates from cell to cell. Schemes such as (48) and (49) express the diffusive flux mostly
as a function of the direct difference of φ across the face, φ(Nf ) − φ(P ), and thus do not allow such
oscillations to go undetected. The gradients ∇a are used in an auxiliary fashion, in order to make the
diffusive fluxes consistent. However, it has recently been shown [43–45] that schemes such as (48) and
(49) can also be interpreted as equivalent to substituting an interpolated value of ∇a φ for ∇φ(cf ) in
(47) and adding a damping term which is a function of the direct difference φ(Nf ) − φ(P ) and of the
gradients ∇a φ(P ) and ∇a φ(Nf ) and tends to zero with grid refinement.
The system of discrete equations (46) is linear, but for convenience it was solved with a fixed-point
iteration procedure where in each outer iteration a linear system is solved (by a few inner iterations
of a preconditioned conjugate gradient solver) whose equations involve only the unknowns at a cell
and at its immediate neighbours, thus avoiding extended stencils. The matrix of this linear system
is assembled only from the parts of Eq. (48) or (49) that directly involve φ(P ) and φ(Nf ), while the
terms involving the gradients are calculated using the estimate of φ from the previous outer iteration
and incorporated into the right-hand side of the linear system, as is customary. Outer iterations were
carried out until the magnitude of the residual per unit volume had fallen below 10−8 in every grid
cell, and the number of required such iterations for each diffusion flux scheme and gradient scheme are
listed in Table 4, where the operator ∇d∞ is obtained with the scheme (22), while the operator ∇dx is
the one obtained by applying the divergence theorem to the auxiliary cell of Fig. 6. Table 4 includes
the percentage of the total calculation time consumed in calculating the gradient. It may be seen that
the cost of corrector steps is quite significant, with the ∇d2 gradient requiring more than 22% of the
total computational effort. On the other hand, the scheme (22) is quite efficient, obtaining the ∇d∞
gradient at a cost that is almost as low as that of ∇d0 ; however, in the case of the over-relaxed scheme,
it also requires somewhat more outer iterations for convergence. The LS gradient costs about the same
as the ∇d1 gradient. The cost of the ∇dx gradient appears somewhat inflated due to a quick but not
very efficient implementation. Overall, it seems that the iterative convergence rate depends mostly on
the chosen flux discretisation scheme rather than on the choice of gradient scheme: in most cases the
over-relaxed scheme converges faster than the custom scheme, which is not surprising since the former
Table 4: Number of outer iterations for convergence (residual per unit volume below 10−8 ) of the over-relaxed
(48) and custom (49) schemes for solving the Poisson equation (42) – (45) on the 512 × 512 distorted grids,
employing various gradient schemes. Also shown is the percentage of CPU time spent on computing the gradient.

∇d0

∇d1

∇d2

∇d∞

∇dx

∇ls

iterations

901

617

601

786

573

585

CPU %

6.8

14.9

22.1

9.7

17.2

14.3

iterations

792

874

914

752

919

919

CPU %

7.4

16.5

23.6

10.5

17.9

15.2

Scheme
Over-relaxed

Custom

30

(a) mean error, mean

(b) maximum error, max

Figure 19: The mean (a) and maximum (b) discretisation errors of various FVM schemes to solve the diffusion
problem (42) – (45), on a series of highly distorted quadrilateral grids (Fig. 20(a)). Grids 0, 1, . . . , 4 have 32×32,
64 × 64, . . . , 512 × 512 volumes, respectively. DT0, DT1, DT2, DT, DTX and LS (solid lines) denote the “overrelaxed” scheme (48) with the ∇d0 , ∇d1 , ∇d2 , ∇d∞ , ∇dx and ∇ls gradient schemes, respectively. DT1c, DTc,
DTXc and LSc (dash-dot lines) denote the “custom” scheme (49) with the respective gradient schemes. C
denotes results on Cartesian grids, where all methods are identical.

is known for its robustness [41]. Nevertheless, the choice of gradient scheme may significantly affect
the iterative convergence rate in the FVM solution of other kinds of problems (depending of course
also on the choice of FVM discretisation and iterative solver). In [26] it is reported that for convection
problems unweighted LS gradients can result in much better iterative convergence than weighted ones,
and that including more distant neighbours in the LS gradient computation also improves convergence.
The discretisation errors with respect to grid refinement are plotted in Fig. 19 for various flux and
gradient scheme combinations. The discretisation error is defined as the difference between the exact
solution φ of the problem (42) – (43) and the numerical (FVM) solution; Fig. 19(a) plots the mean
absolute value of the discretisation error over all grid cells, and Fig. 19(b) plots the maximum absolute
value.
The first thing that one can notice from Figure 19 is that the common DT gradient (∇d0 , ∇d1 ,
d2
∇ ) leads to zeroth-order accuracy with respect to both the mean and maximum discretisation errors,
irrespective of whether the over-relaxed or custom flux scheme is used. Similar trends as those of Fig.
17 are observed: corrector steps reduce the error but are eventually unable to converge to the exact
solution. Although on coarse grids the DT1 and DT2 lines give the impression that they converge,
eventually grid refinement cannot reduce the error below a certain point. In fact, an increase of the
mean error can be observed on grid 4. Things are worse when the DT gradient is used in combination
with the custom scheme (line DT1c), presumably because the gradient ∇a plays a more significant
role in that scheme than in the over-relaxed scheme.
On the other hand, the ∇d∞ gradient obtained through the iterative procedure (22), the “auxiliary
cell” gradient ∇dx , and the LS gradient ∇ls , all lead to second-order convergence to the exact solution.
As long as the gradient scheme is consistent, the FVM accuracy seems to depend not on the gradient
scheme but on the flux scheme, i.e. lines DT, LS and DTX (over-relaxed scheme (48)) nearly coincide,
as do lines DTc, LSc and DTXc (custom scheme (49)). The mean discretisation error of the custom
scheme is about the same as that obtained on Cartesian grids (line C), or even marginally lower when
used in combination with the ∇d∞ gradient! The mean discretisation error of the over-relaxed scheme
is only about 40% higher, which is a very good performance given that it does not account for skewness
and that its iterative convergence is better (Table 4).

31

(a)

(b)

(c)

Figure 20: Samples from each grid category employed to solve the 2D Poisson problems using OpenFOAM:
(a) a grid of distorted quadrilaterals; (b) a grid constructed with the Netgen algorithm; (c) a grid constructed
with the Gmsh algorithm. All domains have unit length along each dimension.

6.2

Tests with OpenFOAM

In order to investigate further the severity of the problem associated with the use of the common DT
gradient we performed also a set of experiments using the popular public domain CFD solver OpenFOAM (https://openfoam.org), version 4.0, (see e.g. [46]) which we have been using recently in our
laboratory [47,48]. We solved again a Poisson problem (42) – (43) but instead of (45) and (44) we have
c(x, y) = sin(πx) sin(πy) and b(x, y) = 2π 2 c(x, y), respectively. The problem was solved on the same
set of grids composed of distorted quadrilaterals (Fig. 20(a)), but, as these are artificially distorted
grids, we also repeated the calculations using a couple of popular grid generation procedures which
are more typical of real-life engineering applications, namely the Netgen [49] (https://ngsolve.org)
and Gmsh Delaunay [50] (http://gmsh.info) algorithms implemented as plugins in the SALOME
preprocessor (www.salome-platform.org). Both generate grids of triangular cells, coarse samples of
which are shown in Fig. 20. The Netgen algorithm can be seen to be effective in producing smooth
grid, reminiscent of that shown in Fig. 5, in four distinct areas within the domain that meet at an
“X”-shaped interface where grid distortion is high. The Gmsh grids are less regular. The “LaplacianFoam” component of OpenFOAM was used to solve the equations with default options, which include
the DT gradient (10) as the chosen gradient scheme (option “gradSchemes” is set to “Gauss linear”).
We repeated the calculations with the gradient option set to LS (gradSchemes: leastSquares).
The discretisation errors are plotted in Fig. 21(a), against the mean cell size4 h. We note that
on the distorted quadrilateral grids and on the Gmsh grids the DT gradient (lines DTQ and DTG )
engenders zeroth-order accuracy, like the DT0 case of Fig. 19. On the other hand, the Netgen grids
are quite smooth (Fig. 20(b)), resembling grids that come from triangulation of structured grids in
most of the domain, and this has the consequence that even with the DT gradient (lines DTN ) the
mean error decreases at a second-order rate. With the LS gradient second-order accuracy is exhibited
on all grids, and furthermore the accuracy is nearly at the same level as on the Cartesian grids.
Three-dimensional simulations were also performed. The governing equations are again (42) – (43)
where we set c(x, y, z) = sin(πx) sin(πy) sin(πz) and b(x, y, z) = 3π 2 c(x, y, z). The domain is the unit
cube z, y, z ∈ [0, 1] and the grids were generated with the Netgen and Gmsh algorithms, where in the
latter we chose the Delaunay method to construct the grid at the boundaries and the frontal Delaunay
algorithm to construct the grid in the interior. The discretisation errors are plotted in Fig. 21(b).
This time it is observed that the DT gradient engenders zeroth-order accuracy even on the Netgen
grids. Presumably, in the 3D case although the Netgen algorithm produces relatively smooth surface
grids on the bounding surfaces of the cube, in the bulk of the domain it packs the tetrahedra in a
way that the skewness is large throughout. On the other hand, the LS gradient again always results
in second-order accuracy.
4

Unlike in previous plots we use the mean cell size instead of the grid level in Fig. 21 because in the case of the Netgen
and Gmsh grids each level does not have precisely 4 (in 2D) or 8 (in 3D) times as many cells as the previous level.

32

(a) mean error, 2D problem

(b) mean error, 3D problem

Figure 21: The mean errors of the OpenFOAM solutions of the 2D (a) and 3D (b) Poisson equations, for
various configurations. The abscissa h is the mean cell size, h = (Ω/M )1/d where Ω = 1 is the domain volume,
M is the total number of grid cells, and d = 2 (for 2D) or 3 (for 3D). The black solid line (C) corresponds to
the 2D problem solved on undistorted Cartesian grids. Red lines (DT) and blue lines (LS) correspond to results
obtained with the DT and LS gradients, respectively. Subscript Q (DTQ , LSQ – solid lines) corresponds to
grids of distorted quadrilaterals (Fig. 20(a)). Subscript N (DTN , LSN – dash-dot lines) corresponds to Netgen
grids (Fig. 20(b)). Subscript G (DTG , LSG – dashed lines) corresponds to Gmsh Delaunay grids (Fig. 20(c)).

7

Final remarks and conclusions

The previous sections showed that on arbitrary grids the DT and LS gradients are zeroth- and firstorder accurate, respectively, but higher orders of accuracy, up to second, are obtained on particular
kinds of grids – a summary of the present findings is given in Table 5.
Unfortunately, general-purpose unstructured grid generation algorithms, such as the popular Netgen and Gmsh algorithms tested in Sec. 6, retain high skewness at all levels of refinement (except
in the 2D Netgen case) resulting in zeroth-order accuracy for both the DT gradient and the FVM
solver that employs it. In our tests, OpenFOAM using the DT gradient was unable to reduce the discretisation errors by grid refinement even though the problems solved were very fundamental, namely
linear Poisson problems with sinusoidal solutions on the simplest possible domains, the unit square
and cube. Zeroth-order accuracy is very undesirable because it deprives the modeller of his/her main
tool for estimating the accuracy of the solution, i.e. the grid convergence study. Therefore the common
DT gradient should be avoided unless the grid refinement algorithm is known to reduce the skewness
towards zero. This is especially important nowadays that automatic unstructured grid generation
algorithms are preferred over structured grids. On the other hand, if the LS gradient is employed
instead then the FVM for the problems of Sec. 6 proved to be second-order accurate on all grids, irrespective of whether their skewness diminishes or not with refinement. The same holds for consistent
DT gradient schemes, such as the proposed iterative scheme (22), the scheme employing the auxiliary
cell of Fig. 6, or any other consistent scheme mentioned in Sec. 3.
On grids where both methods exhibit the same order of accuracy the optimal method is problemdependent. Usually schemes of the same order of accuracy will produce comparable errors on any
given grid but in some cases the errors can differ significantly. For example, in aerodynamics problems
concerning the boundary layer flow over a curved solid boundary, when structured grids consisting of
cells of very high aspect ratio are employed the DT method is known to significantly outperform the LS
method [16, 20] despite both being second-order accurate (or first-order, on triangulated grids). This
is because the DT method benefits from the alignment of the cell faces with the flow (see Appendix
C).
The present study has been restricted to the FVM solution of a Poisson problem. The effect of
33

Table 5: Order of the mean and maximum truncation errors of the gradient schemes on various grid types.
Common DT

Cartesian (§5.2)
Smooth Structured (§5.3)
Locally Distorted (§5.4)
Arbitrary (§5.5)

mean
2
2
1
0

max.
1
1
0
0

LS q 6= 3/2,
∇d∞ , ∇dx
mean max.
2
1
2
1
2
1
1
1

LS q = 3/2
mean
2
2
2
1

max.
2
2
1
1

the choice of gradient discretisation scheme for the FVM solution of flow problems, including nonNewtonian or turbulent ones where the role of the gradient is more significant, will be investigated in
a future study.

Acknowledgements
This research has been funded by the LIMMAT Foundation under the Project “MuSiComPS”.

Appendix A

Accuracy of the least squares gradient when the neighbouring points are arranged at equal angles

Consider the case that all angles between two successive neighbours are equal; this means that they
are equal to 2π/F (F being the total number of neighbours). Therefore,
θf = θ1 + (f − 1)

2π
F

(A.1)

As per Eq. (35), we will express βς in the error expression (31) substituting for ςf from Eq. (25). We
also substitute ∆xf = ∆rf cos θf and ∆yf = ∆rf sin θf . In the special case that all neighbours are
equidistant from P , i.e. all ∆rf ’s are equal, the weights wf are also equal and the products (∆rf )3 wf2
can be factored out of the sums:


1
 2 φ.xx


βτ = (∆rf )3 wf2 
1

2 φ.xx

F
P

(cos θf )3 + φ.xy

f =1
F
P

F
P

f =1

(cos θf

)2 sin θ

f

+ φ.xy

f =1

F
P


cos θf (sin θf )2 
f =1


F
F
P
P

cos θf (sin θf )2 + 21 φ.yy
(sin θf )3

(cos θf )2 sin θf +

f =1

1
2 φ.yy

(A.2)

f =1

where higher-order terms are omitted. If the neighbours are not equidistant but the weight scheme
wf = (∆rf )−3/2 is used, then all products (∆rf )3 wf2 = 1 are again equal and (A.2) still holds.
To proceed further, it will be convenient
to use complex arithmetic: cos θf = (eiθf + e−iθf )/2 and
√
iθ
−iθ
f
f
sin θf = (e − e
)/2i, where i ≡ −1. Then it is straightforward to show that the sums that
appear in the right hand side of Eq. (A.2) are equal to

34

F
X

F

(cos θf )3 =

f =1
F
X


1 X  i3θf
e
+ 3eiθf + 3e−iθf + e−i3θf
8

(A.3)

f =1
F

(cos θf )2 sin θf =

f =1

(A.4)

f =1

f =1
F
X


1 X  i3θf
+ eiθf − e−iθf − e−i3θf
e
8i
F


1 X  i3θf
−i3θf
−iθf
iθf
cos θf (sin θf ) = −
+ e
− e
− e
e
8
2

(A.5)

f =1

F
X

F


1 X  i3θf
(sin θf ) = −
e
− 3eiθf + 3e−iθf − e−i3θf
8i
3

f =1

(A.6)

f =1

P iθ P −iθ P i3θ
P −i3θ
f ; if these sums are zero
The above expressions involve the sums
e f,
e f,
e f and
e
then the sums on the left hand sides of (A.3) – (A.6) are also zero and it follows, in the same way as
for Eq. (35), that the LS gradient is second-order accurate. So, we consider the first of these sums,
substituting for θf from Eq. (A.1).
F
X

iθf

e

=

f =1

F
X

e

i(θ1 +(f −1) 2π
F )

f =1

= e

iθ1

F
X

2π

ei(f −1) F

(A.7)

f =1

We denote ζ = ei2π/F and note that it is the first F -th root of unity, i.e. ζ F = ei2πF/F = ei2π = 1.
Then the sum that appears at the end of Eq. (A.7) is
F
−1
X

ζ n = 1 + ζ + ζ 2 + · · · + ζ F −1 = 0

(A.8)

n=0

Equation (A.8) holds for any F -th root of unity except ζ = 1, as follows from the identity 1 − ζ F =
(1 − ζ)(1 + ζ + ζ 2 + · · · + ζ F −1 ): the left hand side is zero because ζ F = 1, so the right hand side must
be
as well, and since ζ 6= 1, it follows that (1 + ζ + ζ 2 + · · P
· + ζ F −1 ) = 0, i.e. Eq. (A.8) [28]. Thus
P zero
iθ
f
e = 0. In exactly the same manner it can
e−iθf = 0 as well.
Pbei3θshown that
We proceed in the same way for the sum
e f:
F
X
f =1

ei3θf =

F
X

2π
ei3(θ1 +(f −1) F ) = ei3θ1

f =1

F
X

2π

ei3(f −1) F

(A.9)

f =1

This time we denote ζ = ei3(2π)/F . It is also an F -th root of unity, since ζ F = ei3(2π)F/F = ei6π = 1.
So, it will also satisfy Eq. (A.8), unless ζ = 1. Now, ζ = ei2π(3/F ) will be 6= 1 precisely if 3/F is not
an
i.e. if F is neither 1 norP
3. Thus, for F 6= 3 Eq. (A.8) still holds, and Eq. (A.9)
that
P integer,
P says
i3θ
i3θ
−i3θ
f
f
f
e
is zero. But if F = 3 then
e
is not zero. The same holds for the last sum,
e
, as
can be shown in exactly the same way.
So, to summarise, if F > 3 then all the terms in (A.3) – (A.6) become zero and the leading term
of the truncation error vanishes, granting second-order accuracy to the LS gradient, provided that
the q = 3/2 weights scheme is used (or that the neighbouring points are equidistant from P ). But if
F = 3 then not all terms vanish and the leading part of the truncation error does not reduce to zero
– the gradient remains first-order accurate. This is unfortunate, as the F = 3 case is very common,
corresponding to triangular cells.

35

Appendix B

Using a LS weights matrix with off-diagonal elements

Consider the one-dimensional case. For a general matrix W the weighted LS method returns the value
ls (P ) which minimises the quantity
of φ.x

2
F


X X
ls

kW (b − Az)k2 =
wij ∆φj − ∆xj φ.x
(P ) 
i

(B.1)

j=1

It can be shown that if the method returns the exact derivative of any quadratic function φ then it is
second-order accurate. The Taylor expansion of a quadratic function φ is just
1
φ.xx (P ) (∆xf )2
2
We can use this equation to substitute for ∆φj into Eq. (B.1) to get
∆φf = φ.x (P ) ∆xf +

(B.2)


2
F
F
X
X 
X
1
ls
 φ.x (P ) − φ.x
kW (b − Az)k2 =
wij (∆xj )2 
(P )
wij ∆xj + φ.xx (P )
2
i

j=1

j=1

Apparently, if we can choose the weights so that j wij (∆xj )2 = 0 for each row i of W , then the
terms involving φ.xx (P ) will vanish and the remaining quantity will be minimised (actually, made
ls
zero)
P by the 2choice φ.x (P ) = φ.x (P ). Thus the procedure will produce the exact result. The equation
j wij (∆xj ) = 0 can be written in matrix notation as X wi = 0, where X is the 1 × F matrix with
elements X1j = (∆xj )2 and wi is the F × 1 column vector with elements wij . We can therefore set
each row wi of W equal to a basis vector of the null space of X. A simple choice is one where in
row f the first element is (∆x1 )−2 , the f +1 element is −(∆xf +1 )−2 , and the rest of the elements are
zero. W is no longer square but has size (F −1) × F , so the weighted system W Az = W b has one less
equation than the original system Az = b. Equation i of the system W Az = W b has the form


1
1
∆φ1
∆φi+1
ls
φ.x (P )
−
−
(B.3)
=
2
∆xi+1 ∆x1
(∆xi+1 )
(∆x1 )2
P

for i = 1, . . . , F − 1. A little manipulation shows that Eq. (B.3) is just Eq. (B.2) for f = i + 1, where
φ.xx (P ) has been eliminated by using again Eq. (B.2), but for f = 1, to express it as a function of ∆φ1
and ∆x1 . Therefore, this method is nearly equivalent to the unweighted LS solution of the system of
equations (B.2). Extension to the 2D or 3D cases follows the same lines.

Appendix C

The gradient schemes in aerodynamics problems

It is known in aerodynamics that the DT gradient is more accurate than the LS gradient for the
computation of boundary layer flow over a curved solid boundary when grid cells of very high aspect
ratio are employed (typically in excess of 1000) [16, 20]. This is due to a fundamental difference
between the methods: the LS gradient uses the directions in which the neighbours lie, represented by
the unit vectors df = (Nf − P )/kNf − P k, while the DT gradient uses the directions normal to the
faces, represented by the unit vectors nf . This allows the possibility of assisting the DT gradient by
aligning the cell faces with the direction of the exact gradient, when this direction is known a priori.
So, consider a structured grid of high aspect ratio cells over a curved boundary as in Fig. 22 and
assume a boundary layer flow where the contours of the dependent variable φ follow the shape of the
boundary so that φ(P ) = φ(N1 ) = φ(N3 ) ⇒ ∆φ1 = ∆φ3 = 0 and the exact ∇φ is directed normal to
the boundary with approximate magnitude (φ(N2 ) − φ(N4 ))/2∆y. The grid being structured all of
the examined methods are second order accurate, but the unweighted LS gradient is known to perform
very poorly: From expression (32) it follows that because ∆r1 and ∆r3 are much larger than ∆r2 and
∆r4 this scheme relies on information mostly in the directions d1 and d3 , and since the directional
derivatives ∆φ1 /∆r1 and ∆φ3 /∆r3 in these directions appear to be zero, the predicted gradient is
also near zero although the true gradient may be very large. This gross underestimation is known
36

to occur when the ratio of y-displacements γ/∆y of neighbours N1 and N3 to neighbours N2 and N4
(Fig. 22) is larger than one (i.e. when the terms (∂φ/∂y)ls γ are larger than the terms (∂φ/∂y)ls ∆y in
the minimisation expression), and in practical applications it is typically around 50 [16]. As this ratio
is proportional to (∆x)2 /(R ∆y) [16], reducing the size ∆x of the cells in the streamwise direction
rapidly reduces the error. Using a weighted LS method dramatically improves the accuracy, but this
remedy does not work if the grid is triangulated as in Fig. 5.
∆x
ey
ex

n2 , d2

P

d3
n3

γ
N3

N2

n4 , d4
N4

d1
n1

∆y
N1

Figure 22: Structured grid of high aspect ratio cells over a curved boundary (the aspect ratio ∆x/∆y is greatly
downplayed for clarity).

The inaccuracy of the LS gradients in this case is due to the invalidation of their fundamental
assumption that the variation of φ in the neighbourhood of the cell is linear, as the cell size is large
enough for the contours of φ to curve significantly across it. On the other hand, the DT gradient
benefits from the alignment of the normals n2 and n4 of the long faces (which dominate the calculation
due to weighting by face area) with the direction of the exact gradient. Thus, for this particular φ
distribution the given grid alignment assists this gradient scheme to achieve good accuracy. A similar
(somewhat less favourable) situation holds in the triangulated grid case.
Another issue is related to numerical stability. Certain discretisation schemes use the values of
the dependent variable at face centres, evaluated as φ(cf ) ≈ φ(P ) + ∇a φ(P ) · (c − P ). Consider the
same boundary layer flow as before, but for the moment assume a straight boundary so that the grid
looks like that of Fig. 2, on which all the examined gradient schemes reduce to ∇s , Eq. (1). Since
φ(N1 ) = φ(N3 ) = φ(P ), Eq. (1) returns a gradient in the direction (−δyξ , δxξ ) i.e. perpendicular to
the vector δ ξ . Thus φ(c1 ) ≈ φ(P ) + ∇s φ(P ) · (c1 − P ) = φ(P ) = φ(N1 ), which is reasonable. Now
consider a curved boundary case where cell P has exactly the same shape as before but its neighbours
N1 and N3 are tilted to follow the curvature as in Fig. 22. The DT gradient depends only on the face
normals nf which have not changed for cell P and thus it gives the same result for φ(c1 ) as before.
On the other hand, the vectors d1 and d3 have now changed and thus the LS gradient is affected by
the curvature and may no longer be perpendicular to the direction c1 − P ; this leads to a predicted
value for φ(c1 ) that is either larger or smaller than both φ(P ) and φ(N1 ), introducing an oscillatory
variation in the direction d1 that may be detrimental to numerical stability [18]. In other cases this
issue can occur with both the LS and the DT gradients, and so some solvers such as OpenFOAM offer
versions of these gradient schemes that use limiters in order to avoid this problem.

References
[1] J. H. Ferziger and M. Peric, Computational methods for fluid dynamics. Springer, 3rd ed., 2002.
[2] P. J. Oliveira, F. T. Pinho, and G. A. Pinto, “Numerical simulation of non-linear elastic flows with
a general collocated finite-volume method,” J. Non-Newtonian Fluid Mech., vol. 79, pp. 1–43,
1998.
[3] A. M. Afonso, M. S. N. Oliveira, P. J. Oliveira, M. A. Alves, and F. T. Pinho, “The finite volume
method in computational rheology,” in Finite-Volume Methods – Powerful Means of Engineering
Design, ch. 7, pp. 141–170, In-Tech Open Publishers, 2012.

37

[4] A. Syrakos, G. Georgiou, and A. Alexandrou, “Solution of the square lid-driven cavity flow of
a Bingham plastic using the finite volume method,” J. Non-Newtonian Fluid Mech., vol. 195,
pp. 19–31, 2013.
[5] A. Jalali, M. Sharbatdar, and C. Ollivier-Gooch, “An efficient implicit unstructured finite volume
solver for generalised Newtonian fluids,” International Journal of Computational Fluid Dynamics,
vol. 30, no. 3, pp. 201–217, 2016.
[6] C. D. Correa, R. Hero, and K.-L. Ma, “A comparison of gradient estimation methods for volume
rendering on unstructured meshes,” IEEE Trans. Visual Comput. Graphics, vol. 17, pp. 305–319,
2011.
[7] T. J. Barth and D. C. Jespersen, “The design and application of upwind schemes on unstructured
meshes,” in AIAA Paper 89-0366, 1989.
[8] H. Jasak, Error Analysis and Estimation for the Finite Volume Method with Application to Fluid
Flows. PhD thesis, Imperial College, London, 1996.
[9] Ž. Lilek, S. Muzaferija, M. Perić, and V. Seidl, “An implicit finite-volume method using nonmatching blocks of structured grid,” Numer. Heat Transfer, vol. 32, pp. 385–401, 1997.
[10] J. Wu and P. Traoré, “Similarity and comparison of three finite-volume methods for diffusive
fluxes computation on nonorthogonal meshes,” Numer. Heat Transfer B, vol. 64, pp. 118–146,
2014.
[11] F. Moukalled, L. Mangani, and M. Darwish, The Finite Volume Method in Computational Fluid
Dynamics. Springer, 2016.
[12] T. J. Barth, “A 3-D upwind Euler solver for unstructured meshes,” in AIAA Paper 91-1548-CP,
1991.
[13] S. Muzaferija and D. A. Gosman, “Finite-volume CFD procedure and adaptive error control
strategy for grids of arbitrary topology,” J. Comput. Phys., vol. 138, pp. 766–787, 1997.
[14] C. Ollivier-Gooch and M. Van Altena, “A high-order-accurate unstructured mesh finite-volume
scheme for the advection–diffusion equation,” J. Comput. Phys., vol. 181, pp. 729–752, 2002.
[15] F. Bramkamp, P. Lamby, and S. Müller, “An adaptive multiscale finite volume solver for unsteady
and steady state flow computations,” J. Comput. Phys., vol. 197, pp. 460–490, 2004.
[16] D. J. Mavriplis, “Revisiting the least-squares procedure for gradient reconstruction on unstructured meshes,” in AIAA Paper 2003-3986, 2003.
[17] B. Diskin and J. L. Thomas, “Accuracy of gradient reconstruction on grids with high aspect
ratio,” tech. rep., NIA Report No. 2008-12, 2008.
[18] E. Shima, K. Kitamura, and K. Fujimoto, “New gradient calculation method for MUSCL type
CFD schemes in arbitrary polyhedra,” in AIAA Paper 2010-1081, 2010.
[19] L. J. Betchen and A. G. Straatman, “An accurate gradient and Hessian reconstruction method for
cell-centered finite volume discretizations on general unstructured grids,” Int. J. Numer. Methods
Fluids, vol. 62, pp. 945–962, 2010.
[20] E. Sozer, C. Brehm, and C. C. Kiris, “Gradient calculation methods on arbitrary polyhedral
unstructured meshes for cell-centered CFD solvers,” in AIAA Paper 2014-1440, 2014.
[21] A. Syrakos, Analysis of a finite volume method for the incompressible Navier-Stokes equations.
PhD thesis, Aristotle University of Thessaloniki, 2006.

38

[22] M. S. Karimian and A. G. Straatman, “Discretization and parallel performance of an unstructured
finite volume Navier–Stokes solver,” Int. J. Numer. Methods Fluids, vol. 52, pp. 591–615, 2006.
[23] Y. Kallinderis and C. Kontzialis, “A priori mesh quality estimation via direct relation between
truncation error and mesh distortion,” J. Comput. Phys., vol. 228, pp. 881–902, 2009.
[24] R. L. Burden and J. D. Faires, Numerical Analysis. Brooks/Cole, 9th ed., 2011.
[25] Q. Wang, Y.-X. Ren, J. Pan, and W. Li, “Compact high order finite volume method on unstructured grids III: Variational reconstruction,” Journal of Computational Physics, vol. 337, pp. 1–26,
may 2017.
[26] B. Diskin and J. L. Thomas, “Comparison of node-centered and cell-centered unstructured finitevolume discretizations: Inviscid fluxes,” AIAA Journal, vol. 49, pp. 836–854, apr 2011.
[27] L. N. Trefethen and D. Bau III, Numerical Linear Algebra. SIAM, 1997.
[28] G. Strang, Linear Algebra and its Applications. Brooks/Cole, 4th ed., 2006.
[29] S. Muzaferija, Adaptive Finite Volume Method for Flow Prediction using Unstructured Meshes
and Multigrid Approach. PhD thesis, Imperial College, London, 1994.
[30] Y. Liu and W. Zhang, “Accuracy preserving limiter for the high-order finite volume method on
unstructured grids,” Computers & Fluids, vol. 149, pp. 88–99, 2017.
[31] J.-M. Vaassen, D. Vigneron, and J.-A. Essers, “An implicit high order finite volume scheme for
the solution of 3D Navier–Stokes equations with new discretization of diffusive terms,” J. Comput.
Appl. Math., vol. 215, pp. 595–601, 2008.
[32] J. F. Thompson, Z. U. Warsi, and C. W. Mastin, Numerical grid generation: foundations and
applications. North-holland Amsterdam, 1985.
[33] Y. Dimakopoulos and J. Tsamopoulos, “A quasi-elliptic transformation for moving boundary
problems with large anisotropic deformations,” J. Comput. Phys., vol. 192, pp. 494–522, 2003.
[34] A. Sidi, “Review of two vector extrapolation methods of polynomial type with applications to
large-scale problems,” Journal of Computational Science, vol. 3, pp. 92–101, 2012.
[35] U. Trottenberg, C. Oosterlee, and A. Schüller, Multigrid. Academic Press, 2001.
[36] A. Syrakos, G. Efthimiou, J. G. Bartzis, and A. Goulas, “Numerical experiments on the efficiency
of local grid refinement based on truncation error estimates,” J. Comput. Phys., vol. 231, pp. 6725–
6753, 2012.
[37] R. Schneiders, “Refining quadrilateral and hexahedral element meshes,” in Proceedings of the 5th
International Conference on Numerical Grid Generation in Computational Fluid Simulations,
pp. 679–688, Mississippi State University, 1996.
[38] N. Chatzidai, A. Giannousakis, Y. Dimakopoulos, and J. Tsamopoulos, “On the elliptic mesh
generation in domains containing multiple inclusions and undergoing large deformations,” J.
Comput. Phys., vol. 228, pp. 1980–2011, 2009.
[39] B. Diskin and J. L. Thomas, “Notes on accuracy of finite-volume discretization schemes on irregular grids,” Appl. Numer. Math., vol. 60, pp. 224–226, 2010.
[40] P. Traoré, Y. M. Ahipo, and C. Louste, “A robust and efficient finite volume scheme for the
discretization of diffusive flux on extremely skewed meshes in complex geometries,” J. Comput.
Phys., vol. 228, pp. 5148–5159, 2009.
[41] I. Demirdžić, “On the discretization of the diffusion term in finite-volume continuum mechanics,”
Numer. Heat Transfer B, vol. 68, pp. 1–10, 2015.
39

[42] A. Syrakos and A. Goulas, “Estimate of the truncation error of finite volume discretization of the
Navier-Stokes equations on colocated grids,” Int. J. Numer. Methods Fluids, vol. 50, pp. 103–130,
2006.
[43] H. Nishikawa, “Robust and accurate viscous discretization via upwind scheme – i: Basic principle,” Computers & Fluids, vol. 49, pp. 62–86, oct 2011.
[44] A. Jalali, M. Sharbatdar, and C. Ollivier-Gooch, “Accuracy analysis of unstructured finite volume
discretization schemes for diffusive fluxes,” Computers & Fluids, vol. 101, pp. 220–232, sep 2014.
[45] H. Nishikawa, Y. Nakashima, and N. Watanabe, “Effects of high-frequency damping on iterative
convergence of implicit viscous solver,” Journal of Computational Physics, vol. 348, pp. 66–81,
nov 2017.
[46] E. Robertson, V. Choudhury, S. Bhushan, and D. Walters, “Validation of OpenFOAM numerical methods and turbulence models for incompressible bluff body flows,” Computers & Fluids,
vol. 123, pp. 122–145, 2015.
[47] N. K. Lampropoulos, Y. Dimakopoulos, and J. Tsamopoulos, “Transient flow of gravity-driven
viscous films over substrates with rectangular topographical features,” Microfluid. Nanofluid.,
vol. 20, p. 51, 2016.
[48] G. Karapetsas, N. K. Lampropoulos, Y. Dimakopoulos, and J. Tsamopoulos, “Transient flow of
gravity-driven viscous films over 3D patterned substrates: conditions leading to Wenzel, Cassie
and intermediate states,” Microfluid. Nanofluid., vol. 21, p. 17, 2017.
[49] J. Schöberl, “NETGEN an advancing front 2D/3D-mesh generator based on abstract rules,”
Computing and Visualization in Science, vol. 1, pp. 41–52, 1997.
[50] C. Geuzaine and J.-F. Remacle, “Gmsh: A 3-D finite element mesh generator with built-in preand post-processing facilities,” Int. J. Numer. Methods Eng., vol. 79, pp. 1309–1331, 2009.

40

