Journal of Computational Physics 230 (2011) 7670–7686

Contents lists available at ScienceDirect

Journal of Computational Physics

journal homepage: www.elsevier.com/locate/jcp

Mesh quality effects on the accuracy of CFD solutions on

unstructured meshes
Aaron Katz ⇑, Venkateswaran Sankaran
US Army Aeroflightdynamics Directorate (AMRDEC), Moffett Field, CA 94035, United States

a r t i c l e i n f o a b s t r a c t

Article history: The order of accuracy and error magnitude of node- and cell-centered schemes are exam-
Received 22 October 2010 ined on representative unstructured meshes and flowfield solutions for computational
Received in revised form 22 May 2011 fluid dynamics. Specifically, we investigate the properties of inviscid and viscous flux dis-
Accepted 22 June 2011
cretizations for isotropic and highly stretched meshes using the Method of Manufactured
Available online 14 July 2011
Solutions. Grid quality effects are studied by randomly perturbing the base meshes and cat-
aloguing the error convergence as a function of grid size. For isotropic grids, node-centered
Keywords:
approaches produce less error than cell-centered approaches. Moreover, a corrected node-
Unstructured grids
Mesh quality
centered scheme is shown to maintain third order accuracy for the inviscid terms on arbi-
Manufactured solutions trary triangular meshes. In contrast, for stretched meshes, cell-centered schemes are
Error convergence favored, with cell-centered prismatic approaches in particular showing the lowest levels
High-order methods of error. In three dimensions, simple flux integrations on non-planar control volume faces
lead to first-order solution errors, while second-order accuracy is recovered by triangula-
tion of the non-planar faces.
Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction

The need to resolve high Reynolds number viscous flows around complex geometry places stringent demands on com-
putational fluid dynamics (CFD) algorithms. In recent years, unstructured grid schemes have achieved mainstream use for
such problems mainly due to their ability to automatically discretize complex domains. Unstructured CFD codes work on
a variety of cell types including tetrahedra, hexahedra, pyramids, and prisms, and many algorithms exist to compute flow
on these cell types, including cell-centered schemes, node-centered dual schemes, and hybrid schemes [1]. The general
objectives of the current article are to critically assess the accuracy characteristics of different unstructured discretizations
for resolving inviscid and viscous phenomena on a variety of mesh types ranging from isotropic to highly stretched high-
aspect ratio meshes. The specific objectives are to identify optimal schemes that maintain desirable accuracy levels even
for highly perturbed and irregular meshes.
Accuracy assessment of CFD schemes is generally referred to as verification and validation. Verification addresses the de-
gree to which the discretization scheme faithfully represents the continuous equations, while validation assesses the ability
to correctly predict the physics of the flow. The focus of the current study is on the the verification of unstructured discret-
ization schemes. Essentially, two questions need to be answered by verification. First, how fast errors are decreased as the
mesh is refined, which reveals the order of accuracy of a scheme. Second, what the absolute level of error is for a given solu-
tion and mesh, which is useful for comparing the relative accuracy of different schemes. A key aspect of this work is to inves-
tigate the relationship between error convergence and mesh quality since practical CFD meshes often contain irregular mesh

⇑ Corresponding author. Address: M/S 215-1, Ames Research Center, Moffett Field, CA 94035, United States. Tel.: +1 650 604 3697; fax: +1 650 604 5173.
E-mail address: akatz@merlin.arc.nasa.gov (A. Katz).

0021-9991/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved.
doi:10.1016/j.jcp.2011.06.023
 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686 7671

topologies. In this regard, we consider triangles and quadrilaterals in 2D and prismatic meshes in 3D. Mesh quality issues are
addressed by systematically perturbing regular meshes and quantifying the error behavior as a function of these perturba-
tions. We analyze the properties of both scalar model equations as well as the coupled fluid dynamics system. However, we
focus exclusively on mesh effects given a particular solution and do not specifically address system stiffness issues, such as
low Mach number effects.
Key to the process of verification is the definition and measurement of the ‘‘error.’’ In order to compute the error precisely
it is necessary to have an exact solution to the continuous equations. Such exact solutions are difficult to obtain or may not
exist for certain equations. A useful alternative is the Method of Manufactured Solutions (MMS) [2]. With MMS, a ‘‘manu-
factured’’ solution is arbitrarily prescribed, after which an analytic source term is added to the original continuous equations
such that the manufactured solution becomes an exact solution of the modified equation. MMS has achieved widespread use
for code verification both as a means of detecting coding errors and to evaluate the order of accuracy [3,4]. In addition,
researchers have used the method to evaluate grid quality effects [5,6]. In this work we use MMS in two ways – first to eval-
uate the formal order of accuracy of schemes, and second, to compare the magnitude of error for various cell types and for
different forms of manufactured solutions.
Solution errors are obtained by numerically solving the manufactured equations on a fixed domain for a series of increas-
ingly refined meshes. Formally solution errors are defined as the difference between the numerical solution and the exact
manufactured solution, and they represent the so-called solution error of the numerical scheme. In contrast, the truncation
errors, which are obtained by substituting the exact solution into the discrete residual are unreliable indicators of accuracy
since they typically converge at a rate lower than solution errors for unstructured finite volume schemes [7,8]. It is also
worth noting that, in calculating the solution errors, we fix the domain size and refine the mesh within the domain in a reg-
ular way, as opposed to the method of downscaling [9], where the domain size is systematically reduced while the number of
cells remains constant. Since the error convergence behavior in the downscaling method does not always match that ob-
tained with the traditional approach, we use the traditional method of mesh refinement exclusively in this work.
The focus of the current study is the identification of schemes which maintain their design order of accuracy regardless of
mesh perturbations, and exhibit low levels of error for the types of meshes and solutions commonly encountered in CFD.
Specifically, we test the ability of node-centered and cell-centered schemes for both inviscid and viscous flux discretizations.
Our findings generally align well with previous studies [10–12,5,13,14]. For instance, we confirm that the median-dual
approximation used in node-centered schemes is only first order accurate on general quadrilateral grids [12]. In the present
work, we further extend the investigations to address several additional issues such as the proper discretization of the MMS
source terms in solution error analysis, the formulation of formally third-order accurate inviscid node-centered schemes, and
the treatment of non-planar faces in three dimensions. Based on these studies, we make recommendations for the best ap-
proach (or approaches) for practical CFD modeling and simulation.
The paper is organized as follows: first, we briefly describe the spatial discretization schemes tested here. Next we dem-
onstrate how MMS can be used to assess order of accuracy and error, including the proper treatment of the source term. We
then present results of MMS tests for isotropic and stretched grids of various types for inviscid and viscous terms in two
dimensions. Finally, we present key results in three dimensions and offer some conclusions.

2. Description of schemes

The accuracy tests in this work involve four schemes: a standard 2D node-centered scheme, a novel 2D corrected scheme,
a standard 2D cell-centered scheme, and finally an improved 3D prismatic scheme. These schemes all approximate a partial
differential equation of the form

@Q
þ r  F ¼ r  Fv ; ð1Þ
@t
where Q is the vector of conserved variables, F is the inviscid flux vector, and Fv is the viscous flux vector. While Euler and
Navier–Stokes definitions of Q, F, and Fv have been implemented, most of the results in this work pertain to the linear scalar
case, where F = aQ, and Fv = lrQ, with a and l being constants. We consider both inviscid and viscous discretizations for the
2D node-centered and cell-centered schemes, while restricting the corrected node-centered and 3D prismatic schemes to an
inviscid discretization presently. The schemes are described briefly with the aid of Fig. 1.

2.1. 2D node-centered scheme

The 2D node-centered scheme is derived from a Galerkin finite element method on linear triangles [15]. Referring to
Fig. 1(a), the discretization for node 0 can be expressed as

@Q 0 X 1 X1
V0 þ ðF 0 þ F i Þ  A0i ¼ l 1 rQ iþ12  Aiþ12 ; ð2Þ
@t i
2 i
2 iþ2

where A0i is the projected area representing the sum of the dual facets which touch the edge connecting nodes 0 and node i,
and Aiþ1 is the outward pointing projected area of the face connecting nodes i and i + 1. In the above formula, the scalar vis-
2
7672 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686

i+1

j+1
i i
0 0
j

i−1
(a) node-centered stencil (b) cell-centered stencil (c) prismatic cell
Fig. 1. Cell-centered and node-centered stencils in two dimensions. Three-dimensional prismatic cell.

cous flux has been used, where liþ1 is the average value of l in the triangle formed from nodes (0, i, i + 1), and rQ iþ1 is the
2 2
gradient of Q in the same triangle, obtained with a trapezoidal integration around the triangle perimeter. This viscous for-
mulation follows the standard Galerkin procedure, and is exact for linear functions. The Navier–Stokes fluxes are of the same
form, and no generality is lost by showing the discretization for the scalar case. In order to enforce upwinding, an artificial
diffusion term of the form
X  
D0 ¼  d0i Q R0i ; Q L0i ; A0i ð3Þ
i

is added to the left-hand side of Eq. 2. For the Navier–Stokes equations, we use the CUSP formulation of Jameson [16,17],
although other formulations for the artificial diffusion may equivalently be used. To obtain orders of accuracy higher than
unity, the approximations Q R0i and Q L0i must be reconstructed using solution data at the surrounding nodes.In this work a
weighted linear least squares method is used to compute the gradients of the reconstruction. Thus the scheme remains exact
for linear functions. The weights of the least squares procedure are inversely proportional to the distance between nodes.
While the linear node-centered scheme resembles a finite volume scheme, it may also be derived using a Galerkin-
weighted weak form of the governing equation on linear triangles. Details of the equivalence of this method to a Galerkin
method, including more details on the viscous discretization are described by Barth [15]. While we have illustrated the
scheme for triangles, the inviscid scheme has also been implemented for quadrilaterals. The dual face areas are formed in
the same way as for triangles, i.e., by connecting the centroids of the cells to the face centers.

2.2. 2D corrected node-centered scheme

Along with the node-centered scheme just described, we also consider a new scheme, referred to as the ‘‘corrected’’
scheme. The corrected scheme is so-named because it involves a correction to the linear node-centered scheme to construct
a formally third order accurate scheme for arbitrary triangular meshes. The condition for the correction and the artificial dif-
fusion term is that all gradients are computed using a weighted quadratic least squares method instead of a linear method.
We note that so far this correction has only been implemented for the inviscid terms, and the viscous terms remain formally
second order accurate. The correction may be interpreted as a flux reconstruction, where the inviscid terms in Eq. 2 are mod-
ified to become
X1 X 1 1 1
 X1
ðF L þ F R Þ  A0i ¼ F 0 þ Dr 0i  rF 0 þ F i  Dr0i  rF i  A0i ¼ ðF 0 þ F i Þ  A0i  C 0 ;
i
2 i
2 2 2 i
2

where the correction, C0, is

X1
C0 ¼ ðDr 0i  ðrF i  rF 0 ÞÞ  A0i : ð4Þ
i
4

Here, Dr0i is the vector from node 0 to node i, rF0 and rFi are the gradients of the flux at nodes 0 and i, and A0i is the median-
dual directed face area associated with the edge connecting nodes 0 and i. These parameters are shown in Fig. 2. The form of
the artificial diffusion remains unchanged for the corrected scheme. No higher order quadrature is used as for cell-centered
quadratic schemes [13]. The only additional expense over the linear scheme is the computation of the gradient weights using
a quadratic approximation, which may be done as a preprocessing step, and the formation of the correction term, which in-
volves additional gradients of the flux. A quadratic least squares approximation to the flux may be constructed by finding a
quadratic function, F, that minimizes the function (see Fig. 3)
X
CðrF 0 Þ ¼ wi ðFðri Þ  F i Þ;
i

where wi are inverse distance weights of the least squares procedure. In other words, we seek to minimize the difference
between some quadratic function and the discrete values of F. The form of F is
 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686 7673

Fig. 2. Illustration of correction technique to obtain third order accuracy on arbitrary triangular grids.

Fig. 3. Node-centered grid in one dimension.

1
F ¼ F 0 þ Dr T0i rF 0 þ Dr T0i r2 F 0 Dr 0i ;
2
where r2F0 is the Hessian matrix of second derivatives. While the result of this procedure is an estimate of both first and
second derivatives of F, we emphasize that second derivative terms are never used and need not be stored. Only the first
derivative terms enter into the correction terms in Eq. 4. The formal third-order accuracy of the corrected inviscid flux will
be demonstrated later.

2.3. 2D cell-centered scheme

Along with the nodal schemes just described, a traditional 2D cell-centered scheme is also tested. The formulation of the
cell centered scheme is similar to that of Eqs. 2 and 3, but with the convective fluxes using left and right states as in
X 1  L   
F Q 0i þ F Q R0i  A0i : ð5Þ
i
2

In the cell-centered formulation, the state variables are located at the cell-centers, and the directed area, A0i, is the area of the
face separating cells i and i + 1 in Fig. 1(c). The gradients of the reconstruction are formed by first obtaining nodal values of
the conserved variables with an inverse-distance weighted linear least squares procedure. These nodal values are used in a
Green-Gauss gradient integration around the perimeter of each triangle. This hybrid cell-nodal approach has been used often
in the literature for achieving second order accuracy on unstructured grids [10,18–20]. While not included in this paper, we
have also implemented other gradient estimation methods, including a direct cell least squares method. This method gave
nearly identical results as the hybrid cell-nodal approach and is therefore not discussed further in this work. The viscous
terms for the cell-centered scheme are discretized as
X  
F v Q jþ1 ; rQ jþ1  A0i ; ð6Þ
2 2
j

where Q jþ1 and rQ jþ1 are the solution and gradient midway between nodes j and j + 1 at the quadrature point. The solution
2 2
Q jþ1 , at the quadrature point is the average of the nodal values at j and j + 1 obtained with the inverse-distance weighted
2
linear least square method. The gradient at the quadrature point is found by enforcing the conditions
rQ jþ12  Dr 0i ¼ Q i  Q 0 ; rQ jþ12  Dr jjþ1 ¼ Q jþ1  Q j ; ð7Þ

where r0i is the vector joining cell centers 0 and i, rjj+1 is the vector joining nodes j and j + 1. The cell-centered scheme works
with both triangle and quadrilateral cell types.

2.4. 3D cell-centered prismatic scheme

The 3D prismatic solver is also based on a cell-centered formulation. The solution variables are stored at the cells and are
obtained at the nodes with an inverse-distance weighted linear least squares method. The nodal values are then used to ob-
tain the gradients at each cell using the Green-Gauss method. The viscous terms for the 3D prismatic code have not yet been
implemented. A subtle, but important finding in this work deals with the method of inviscid flux integration. Polyhedral vol-
umes in three dimensions, other than tetrahedra, are composed of non-planar faces in general. Usually, the flux integration
7674 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686

scheme uses a single quadrature point at the center of each non-planar face. As an alternative, we also implement a flux inte-
gration scheme in which each non-planar face is triangulated, with a quadrature point placed at the center of each resulting
triangular facet. For quadrilateral faces encountered in prismatic grids, this results in two quadrature points on each non-
planar face. Triangulation of non-planar faces is shown in Fig. 1(c). Both single point and triangulated quadrature schemes
are implemented in the 3D prismatic solver for purposes of comparison.

3. Method of Manufactured Solutions

In order to formally assess the accuracy of the discretization methods listed in the previous section, exact solutions are
required. Most exact solutions for compressible viscous flows are quite simple and do not exercise all the terms in the gov-
erning equations. To address this problem, Roache [21] proposed the Method of Manufactured Solutions, or MMS. Here, in-
stead of solving Eq. 1 directly, we solve the equation augmented with an analytic source term,
@Q
þ r  F  r  F v ¼ S: ð8Þ
@t
The MMS procedure involves selecting an arbitrary ‘‘manufactured’’ solution and then substituting the manufactured solu-
tion into the original continuous differential equation. The manufactured solution will not in general identically satisfy the
differential equation, and the remainder is set equal to the source term [4]. Thereby, the manufactured solution represents
the exact solution of the modified continuous differential equation, i.e., the original equation with the source term added. We
note that the source term is not a function of Q, but is only a function of the independent variables and parameters of the
PDE.
Certain subtleties that are not well-understood in the literature arise when using MMS in conjunction with grid refine-
ment as a means to measure error convergence. In particular, we discuss the difference between truncation and solution er-
ror as well as the notion of source term discretization. We also show grid convergence results using an exact solution of the
Euler equations, a manufactured solution of the Euler equations, and a manufactured scalar solution. The objective of these
results is to verify the MMS procedure so that we will have confidence in the method to evaluate grid quality effects in our
subsequent tests.

3.1. Truncation error versus solution error

One common misconception in error analysis for unstructured grid schemes stems from the difference between trunca-
tion error and solution error. The two types of error do not necessarily converge at the same rate, a fact which has sometimes
led to confusion and erroneous conclusions [5,14]. Truncation error may be defined as the residual present when an exact
solution is introduced into a discretization. If no residual were present, the discretization would be exact. For example, con-
sider a discretization of a linear form of Eq. 8 at steady state. If known Dirichlet conditions Qb are supplied at some boundary
values, a linear discretization may be expressed as

DQ h ¼ BQ b þ CS; ð9Þ
where D represents the discrete scheme operating on the discrete solution, Qh, B incorporates the Dirichlet conditions, and C
operates on the known source terms, S, obtained from the MMS procedure. A common treatment of the source term is to set
C = I, where I is the identity matrix. However, other source term treatments are possible and even necessary to obtain correct
results using manufactured solutions, as we will show. Additionally, B is a matrix which contains zero-valued entries for
interior degrees of freedom, with non-zero-valued entries for degrees of freedom touching Dirichlet boundaries. In this case
the entries are dictated by the particular discretization method. Note that the discrete solution exactly satisfies Eq. 9 while
the exact solution, Q, applied to the discrete equations yields an error,
DQ ¼ BQ b þ CS þ Et ; ð10Þ
where Et is the truncation error.
On the other hand, the solution error is defined as the difference between the exact and discrete solutions. Using the
above definitions, we can relate the solution error to the truncation error by

Es ¼ Q  Q h ¼ D1 Et ð11Þ
or
DEs ¼ Et : ð12Þ
It should be emphasized that this result holds only for linear equations. For general equations, an additional linearization
step would be required first. However, the result is included here to emphasize the difference between truncation and solu-
tion error. The solution error is assumed to be dominated by terms of leading order p such that
p
Es ¼ Q  Q h ¼ K Dh þ H:O:T:; ð13Þ
 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686 7675

where K consists of constants independent of the mesh size, h, and H.O.T. stands for high order terms.
A few observations can be made from from the preceeding analysis:

 The solution error is a solution to the PDE governing Q when driven by the local truncation error as a source term. It is
appropriate to think of the truncation error as driving the solution error. If the truncation error were zero at every node,
the solution to Eq. 12 would be Es = constant for all nodes, and since some of the nodes must necessarily be Dirichlet con-
ditions for which the constant is zero, the constant must be zero everywhere, resulting in zero solution error. However,
for the general case, the truncation error is non-zero, and would drive a non-zero solution error.
 The order of convergence of solution error need not be the same as the order convergence of truncation error for finite
volume schemes [22,7,8,23]. Therefore, truncation error used in this way does not yield a reliable estimate of the order
of accuracy. For example, the truncation error for the linear node-centered scheme is dominated by first order terms for
arbitrary meshes. However, the solution error improves to second order accuracy.
 We must actually solve the system of equations to determine the solution error, which is of primary interest. Simply com-
puting the residual from the substitution of the exact solution only gives truncation error.

3.2. Source term discretization

Despite the fact the MMS source is known exactly at every location in the domain, a source term discretization is still
necessary, as symbolized by the operator C in Eq. 9. Often C = I, a simple point source discretization, is used since the source
is known analytically at each location in the flow. However, this simple source term discretization can lead to incorrect error
convergence computations, especially for high-order methods. While we desire to understand the effect of the source term
discretization on solution error, it is insightful to observe its effects on truncation error, which itself acts as a source term for
the solution error, as shown in Eq. 12.
As an example, consider the linear advection equation in one dimension, which only has a non-trivial solution with an
MMS source term:

aQ x ¼ sðxÞ: ð14Þ
Discretizing using a node-centered central difference approach with artificial diffusion to introduce upwinding, the approx-
imation becomes
a
ðQ  Q i1 Þ  Di ¼ C i S; ð15Þ
2Dxi iþ1
where Di consists of artificial diffusion/upwinding terms, CiS is the source term discretization at node i, and
Dxi ¼ 12 ðDxi1=2 þ Dxiþ1=2 Þ, where D xi+1/2 = xi+1  xi. The truncation error can be determined by Taylor series analysis for var-
ious source term discretizations, Ci. If we use a simple point source discretization CiS = si, the truncation error becomes
! ! !
Dx2iþ1  Dx2i1 aQ  Dx3iþ1 þ Dx3i1 Dx4iþ1  Dx4i1 aQ 
2x 4x
Et ¼  2 2
 2 2
ðaQ 3x Þ  2 2
þ Di þ H:O:T:
Dxi 4 12Dxi D xi 48

If we use linearly reconstructed states for the artificial diffusion operator then Di  O(h3) for a regular mesh and Di  O(h1)
for a general irregular mesh. Therefore, the truncation error for the point source discretization becomes Et  O(h2) for a reg-
ular mesh and Et  O(h1) for a general irregular mesh.
Alternately, since the node-centered finite volume scheme actually corresponds to a linear Galerkin finite element
scheme, we could use the Galerkin framework to derive the source term discretization. The Galerkin-weighted source term
discretization at node i is

2 1  
CiS ¼ si þ si1 Dxi1 þ siþ1 Dxiþ1 : ð16Þ
3 6Dxi 2 2

The truncation error using the Galerkin-weighted source term discretization becomes
! !
Dx2iþ1  Dx2i1 aQ sx
 Dx4iþ1  Dx4i1 aQ s3x

2x 4x
Et ¼  2 2
  2 2
 þ Di þ H:O:T:;
Dxi 4 6 Dxi 48 36

using the fact that aQ3x  s2x = 0 from the derivative of Eq. 14. Again, using linearly reconstructed states for the artificial dif-
fusion term, the truncation error for the point source discretization becomes Et  O(h3) for a regular mesh and Et  O(h1) for a
general irregular mesh. By using the Galerkin-weighted source term discretization we obtain a truncation error which is one
order more accurate than that for a regular mesh. Our numerical experiments indicate that the solution error likewise im-
proves by one order of accuracy on a regular mesh by using the Galerkin source discretization. Importantly, we will see that
the corresponding MMS results display the same trends as the solution error obtained using exact solutions.
Table 1 summarizes the effect of the source term discretization on an equivalent 2D scalar advection problem for which
non-trivial exact solutions exist with no MMS source term. The mesh types used for this study are shown in Fig. 4. This
7676 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686

Table 1
Truncation and solution error convergence rates for different source term discretization methods for exact and manufactured solutions using regular and
perturbed grids composed of equilateral triangles.

Source discretization Truncation error Solution error

Exact MMS Exact MMS
– Point Galerkin – Point Galerkin
Regular grid 3 2 3 3 2 3
Perturbed grid 1 1 1 2 2 2

experiment highlights several items of interest. First, different convergence rates are obtained for truncation and solution
error for the same series of meshes and solutions. The difference only appears for irregular grids. Second, the Galerkin meth-
od of MMS source term discretization leads to error convergence rates that match those of the exact solution with no MMS
source terms. On the other hand, the point source method does not accurately reflect the true order of accuracy for regular
grids. Finally, the convergence rate of the solution error degrades by one order for perturbed grids. More detailed studies of
these phenomena are considered later.
In summary, MMS solution error calculations appear to give the same order of convergence as exact solutions (no MMS
source term) only when a source term discretization consistent with the overall discretization method for the PDE is used.
While no rigorous proof has been devised yet, it appears that a sufficient condition on the source term discretization is equal-
ity of the order of the truncation error for the MMS modified PDE in Eq. 8 and the original PDE in Eq. 1. This is a convenient
check on the correctness of any source term discretization, since truncation errors can easily be checked, either analytically
or computationally (by simply computing the residual) without resolving the discrete solution. A dependence of error con-
vergence on the source term discretization was found to be true by Pautz [3] in his experiments, but we are not aware of
studies of the design of optimal source discretizations. We further note that proper discretization of the source terms
may play an important role for using MMS with high-order schemes.

3.3. Verification of the MMS procedure

For a more detailed verification of the MMS procedure itself, we perform an error convergence study using an exact solu-
tion to the Euler equations, a manufactured solution to the Euler equations, and a manufactured solution to the scalar con-
vection equation in two dimensions. The error convergence characteristics for each of these solutions should be the same if
we are to trust MMS as a means of verification. The exact solution to the Euler equations we use is Ringleb flow [24], a por-
tion of which is shown in Fig. 5(a). The Euler manufactured and scalar manufactured solutions are shown in Fig. 5(b) and (c).
The Euler manufactured solution is chosen similar to the solution of Roy [4] to be
a  a   
qx px qy py aqxy pxy
qðx; yÞ ¼ q0 þ qx sin þ qy sinþ qxy cos ;
L L L2
a px a py  
ux uy auxy pxy
uðx; yÞ ¼ u0 þ ux cos þ uy cos þ uxy cos ;
L L L2
a px a py  
av xy pxy
v ðx; yÞ ¼ v 0 þ v x sin v x þ v y sin
vy
þ v xy cos ;
L L L2
a px a py  
Px Py aPxy pxy
Pðx; yÞ ¼ P 0 þ Px sin þ Py cos þ P xy sin :
L L L2
The scalar manufactured solution is of a similar trigonometric form:

Fig. 4. Regular and perturbed isotropic grids composed of equilateral triangles used for MMS verification study.
 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686 7677

Fig. 5. Exact and manufactured solutions used for error convergence studies.

a  
/x p x a
/y py a/xy pxy
/ðx; yÞ ¼ /0 þ /x sin þ /y sin þ /xy cos :
L L L2
Systematic grid convergence tests are performed using the node-centered and corrected schemes on grids composed of iso-
tropic equilateral triangles. Both regular and perturbed triangles are tested. In all instances in this work, perturbed meshes
are created from corresponding regular meshes. The regular locations of the nodes are moved a random distance between
zero and 25% of the minimum distance between adjacent nodes, and in a random direction between zero and 360 degrees
in the plane of the mesh. Each mesh on each refinement level is independently perturbed randomly from its regular parent
mesh in the same manner and is scaled by the same fraction of its mesh size (statistically), such that mesh quality remains
constant throughout the refinement process. As such, meshes are not given the benefit of being of the same ‘‘family’’ as is
typically done for structured meshes, in which every other mesh line is removed, for example. This unstructured grid refine-
ment procedure is similar to that used by Diskin [10,11] and reveals the true behavior of schemes. The results of the veri-
fication test for each scheme, solution type, and mesh are shown in Table 2. As the table shows, the exact Euler solution,
MMS Euler solution, and MMS scalar solutions all lead to identical estimates of the solution error convergence rate on per-
turbed and uniform grids.
Consistent with the discussion in Section 3.2, the results shown in Table 2 indicate third order accuracy for the node-cen-
tered linear scheme on unperturbed meshes. When the MMS source term does not receive the proper Galerkin weighting,
only a second order result is obtained. This is likely the reason for the second order results observed by Diskin et al [12].
But the scheme is truly third order for regular meshes, a fact that is confirmed by the exact Ringleb flow results. A further
significant result is that the corrected scheme maintains third order accuracy even for perturbed grids. Unconditional third
order accuracy was previously obtained using schemes with high-order quadrature and second derivative computation [13].
In the case of the corrected scheme presented here neither of these expensive procedures are required. These results will be
confirmed with further tests in a later section.

Table 2
Verification of the MMS procedure for the node-centered and corrected schemes using Ringleb, manufactured Euler, and manufactured scalar solutions. Table
values are solution error convergence rates for regular and perturbed grids.

Ringleb Euler MMS Scalar MMS

Regular Perturbed Regular Perturbed Regular Perturbed
Node-centered 3 2 3 2 3 2
Corrected 3 3 3 3 3 3
7678 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686

The results of these tests support the notion that MMS is an accurate and useful method for verification as long as the
proper source term discretization is used. Consequently, the remaining tests in this work make use of manufactured solu-
tions and scalar equations to represent error characteristics of the Euler or full Navier–Stokes equations. This allows us
the freedom to choose solution forms of practical interest, such as the boundary layer type solutions of Fig. 5(d) and (e). Sim-
ilar to the manufactured solution of Sun et al. [6], the solution of Fig. 5(d) has the form
pffiffiffiffiffiffiffiffiffiffi
0 ffi
ðyy Þ

Qðx; yÞ ¼ 1  e clðxx0Þ
ð17Þ
6
with c = .59 and l = 10 . The solution of Fig. 5(e) has a similar form, expressed in polar coordinates as
pffiffiffiffiffiffiffiffiffiffiffi ðrr0 Þ

Qðr; hÞ ¼ 1  e clðh0 hÞ : ð18Þ

Verification results of the different schemes for different grid and solution types are given in the following section.

4. Results

We now present detailed results for isotropic grids, stretched grids, and 3D prismatic grids. The goal the isotropic grid
tests are to verify solution error convergence rates. The section on stretched meshes aims to examine how these error char-
acteristics change as a function of the mesh stretching, which is controlled by the wall cell aspect ratio. Finally, the 3D pris-
matic grid studies examine the effects of non-planar faces on order of accuracy.

4.1. Isotropic grids

The first way in which we use MMS is to verify the order of accuracy of the node-centered linear, node-centered corrected,
and cell-centered schemes on a variety of two dimensional grids. For the order of accuracy tests, we use the isotropic grids
shown in Fig. 6 and the scalar manufactured solution shown in Fig. 5(c). Both uniform and randomly perturbed grids com-
posed of quadrilaterals, equilateral triangles, and right triangles are tested. A series of six meshes of each type are created by
successively reducing the cell size by a factor of two. The cell size for a given mesh and scheme is defined as
 1
V total d
ds ¼ ; ð19Þ
ndof
where Vtotal is the total volume of all cells in the domain, ndof is the number of degrees of freedom in the mesh, and d is the
number of spatial dimensions. For a given refinement test, Vtotal remains constant and ndof increases as the mesh is refined.
For node-centered schemes, ndof is set to the number of nodes in the mesh, while for cell-centered schemes ndof is set to the
number of cells. This distinction is important in comparing the relative accuracy per degree of freedom for node- versus cell-
centered schemes. This is because there are roughly twice the number of cells as nodes in 2D triangular meshes, and five to

Fig. 6. Isotropic grids used for order of accuracy study.

 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686 7679

six times the number of cells as nodes in 3D tetrahedral meshes. Quadrilateral and hexahedral meshes contain the same
number of cells as nodes, excluding boundary effects, while prismatic meshes may contain differing numbers depending
on the surface topology. The perturbed meshes are constructed in the manner described in Section 3.3.
The results of the order of accuracy studies on the isotropic grids of Fig. 6 for scalar manufactured solutions are shown in
Fig. 7. We test the inviscid and viscous terms in isolation to examine the order of accuracy of each discretization. In each plot
of Fig. 7, the L2 norm of the solution error defined as Qcomputed  Qexact is plotted on the y-axis, and the inverse of the cell size
given in Eq. 19 is plotted on the x-axis. The figures show the most important aspects of this study. A more comprehensive
summary of the reported orders of accuracy for all schemes tested are given in Tables 3 and 4 for inviscid and viscous
schemes, respectively. Solution error convergence rates are shown for regular and perturbed grids.
As Table 3 shows, the inviscid node-centered scheme on quadrilateral meshes reduces to first order accuracy compared to
the identical scheme on equilateral or right triangles, which is second order accurate. This is consistent with the finding of
Diskin et al. [12]. Diskin further shows that second order accuracy can be recovered on quadrilateral grids by using flux inte-
gration on the individual facets of the dual control volume, but with a dramatic increase in computational expense and
complexity.
Table 3 also shows that the inviscid node-centered linear scheme on both equilateral and right triangles is actually third
order accurate for regular meshes. This result is somewhat surprising since, in the literature, we observe that linear schemes
are usually reported to be second order accurate. Furthermore, we observe that the corrected scheme maintains third order
accuracy, even for randomly perturbed meshes. It accomplishes this without additional quadrature or the need to retain sec-
ond derivative information for quadratic reconstruction, as is the case for cell-centered quadratic schemes. As stated previ-

-1 -1
10 10

-2 -2
10 10
node-cent. Qaud, pert. cell-cent. Qaud, pert.

10-3 10-3
L2 error

L2 error

1
1 1
1
10-4 10-4 2
2
node-cent. Qaud, reg.
cell-cent. Eq. Tri, pert.
10-5 node-cent. Eq. Tri, pert. 10-5
3 cell-cent. Rt. Tri., pert.
corrected, Eq. Tri., pert. 1
-6 -6
10 10
500 1000 1500 2000 500 1000 1500 2000
-1 -1
ds ds

-1 -1
10 10

cell-cent. Quad., pert.

-2 -2
10 10
node-cent. Rt. Tri., pert. cell-cent. Rt. Tri., pert.

10-3 10-3
L2 error

L2 error

1 1
2 2
10-4 10-4

node-cent. Rt. Tri., pert. cell-cent. Eq. Tri., pert.

10-5 10-5

-6 -6
10 10
500 1000 1500 2000 500 1000 1500 2000
-1 -1
ds ds

Fig. 7. Solution error convergence for invscid and viscous MMS scalar schemes on regular and perturbed isotropic meshes composed of quadrilaterals,
equilateral triangles, and right triangles.
7680 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686

Table 3
Solution error convergence rates for all inviscid schemes tested as predicted by the isotropic scalar MMS test for regular and perturbed grids.

Quadrilateral Equilateral triangle Right triangle

Regular Perturbed Regular Perturbed Regular Perturbed
Node-centered 2 1 3 2 3 2
Corrected – – 3 3 3 3
Cell-centered 2 2 2 2 2 2

Table 4
Solution error convergence rates for all viscous schemes tested as predicted by the isotropic scalar MMS test for regular and perturbed grids.

Quadrilateral Equilateral triangle Right triangle

Regular Perturbed Regular Perturbed Regular Perturbed
Node-centered 2 2 2 2 2 2
Cell-centered 2 2 2 2 2 2

ously, the only additional expense for the corrected scheme over the linear node-centered scheme is the expense of comput-
ing the gradients of the flux. The gradients of the flux may be computed in a manner identical to the gradients of the solution
variables themselves, which are needed anyway in the artificial diffusion formulation. These results seem to indicate that for
inviscid schemes on isotropic meshes, node-centered schemes are preferred over cell-centered schemes.
A final result from the isotropic order of accuracy tests is that both the cell-centered and node-centered viscous discret-
izations are second order accurate for both uniform and perturbed meshes of all cell types, as shown in Table 3. Note that a
viscous scheme for node-centered quadrilaterals was not tested due to the first order results obtained for the inviscid equa-
tion. In other words, the node-centered linear scheme on quadrilaterals is no longer considered a competitive method since
it is easily outperformed by other methods.

4.2. Stretched grids

We next turn our attention to the performance of the schemes on stretched meshes of various types. A total of six mesh
types are used, shown in Fig. 8. Flat and curved surface boundary layer type meshes are tested, along with the corresponding
scalar manufactured solutions shown in Fig. 5(d) and (e) and described in Eqs. 17 and 18. Three levels of mesh stretching are
tested, including wall cell aspect ratios of 102, 104, and 106. These mesh shapes, cell types, and manufactured solutions are

Fig. 8. Stretched grids used to assess grid quality effects.

 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686 7681

representative of typical high Reynolds number Navier–Stokes flows. As in the isotropic tests, purely inviscid and purely vis-
cous discretizations are tested independently to isolate the effects of mesh stretching on each discretization type.
The solution error convergence for inviscid terms using the flat surface meshes and scalar boundary layer type solutions
are shown in Fig. 9. Fig. 9(a) shows error convergence results for the node-centered and corrected schemes for flat surface
meshes composed of right triangular meshes of wall cell aspect ratio 106. The corrected scheme maintains its third order
accuracy despite the stretching. A similar study is performed for the cell-centered inviscid scheme using the same scalar
boundary layer type solution and flat stretched mesh, showing uniform second order accuracy in Fig. 9(b). For the cell-cen-
tered inviscid scheme, quadrilateral grids produce slightly less error than triangular grids.
A comparison of the inviscid node-centered and corrected schemes for right triangle flat surface grids of various aspect
ratios is shown in Fig. 9(c). The corrected scheme again shows third order accuracy, while the linear node-centered scheme
shows second order accuracy for all wall cell aspect ratios. As the wall cell aspect ratio increases so does the error, likely due
in part to increased stretching since the domain remains fixed.
While the inviscid study on stretched meshes seems to indicate only a slight advantage of the cell-centered scheme over
the node-centered scheme the situation becomes more pronounced when the error in the viscous discretization is consid-
ered. Fig. 10 shows solution error convergence results for the viscous terms using flat surface stretched meshes and a scalar
boundary-layer type solution. Fig. 10(a) shows that for an aspect ratio of 106, right triangles and equilateral triangles perform
similarly for the node-centered viscous discretization. However, for the cell-centered viscous scheme, Fig. 10(b) shows vary-
ing levels of error for different cell types. Quadrilaterals produce significantly lower error, followed by right triangles, then

0 0
10 10
6
node-cent. Rt. Tri., A=10 cell-cent. Rt. Tri., A=10 6
-1 -1
10 10

10-2 10-2
1
L2 error

L2 error

2 1
10-3 10-3 2
6
cell-cent. Eq. Tri., A=10
3
-4 -4
10 10 6
1
6
cell-cent. Quad., A=10
-5
corrected Rt. Tri., A=10 -5
10 10

-6 -6
10 10
500 1000 1500 2000 500 1000 1500 2000
-1
ds ds-1

0
10

node-cent. Rt. Tri.

-1
10 1
2
-2
10
L2 error

3
10-3
A=10 2
1 }
-4
A=10 4
10 A=10 6

10
-5
A=10 2
A=10 4
A=10 6
}
corrected Rt. Tri.
-6
10
500 1000 1500 2000
-1
ds

Fig. 9. Effects of mesh stretching for flat surface meshes on node-centered, corrected, and cell-centered inviscid schemes using a scalar boundary layer type
MMS solution and various cell types. Cell types include right triangles, equilateral triangles, and quadrilaterals (cell-centered only).
7682 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686

0 0
10 10
6
node-cent. Rt. Tri., A=10 cell-cent. Eq. Tri., A=10 6

cell-cent. Rt. Tri., A=10 6

-1 -1
10 10
L2 error

L2 error
1
10
-2
10
-2 1
2
2

10-3 node-cent. Eq. Tri., A=10 6 10-3

cell-cent. Quad., A=10 6

500 1000 1500 2000 500 1000 1500 2000

-1 -1
ds ds

A=10 2 A=10 2
10
0 A=10 4 10
0 A=10 4
6 6
A=10 A=10
A=10 2 A=10 2
A=10 4 A=10 4
6 6
A=10 A=10
-1 -1
10 10
node-cent. Eq. Tri. node-cent. Rt. Tri.
L2 error

L2 error

1 1
-2 -2
10 10
2 2

10-3 }} 10-3
}}
cell-cent. Eq. Tri. cell-cent. Rt. Tri.
500 1000 1500 2000 500 1000 1500 2000
-1 -1
ds ds

Fig. 10. Effects of mesh stretching for flat surface meshes on node-centered, and cell-centered viscous schemes using a scalar boundary layer type MMS
solution and various cell types. Cell types include right triangles, equilateral triangles, and quadrilaterals (cell-centered only).

equilateral triangles. Even more pronounced are the differences between the node-centered and cell-centered schemes as
the wall cell aspect ratio increases, shown in Fig. 10(c) and (d) for equilateral and right triangles, respectively. Both figures
show that the cell-centered viscous formulation produces roughly four times less error than the node-centered formulation.
This evidence points to an advantage of the cell-centered scheme over the node-centered scheme for stretched meshes
where viscous terms dominate. This is in contrast to isotropic results, which indicate an advantage for the node-centered
schemes. We also note that there is a slight decrease in accuracy as the wall cell aspect ratio increases for both cell- and
node-centered schemes.
Finally, a similar study is performed using stretched meshes over a curved surface instead of a flat surface. The corre-
sponding boundary layer type scalar MMS solution following the curvature of the mesh in Eq. 18 is used. Node-centered
and cell-centered schemes are compared using right triangles and quadrilaterals typically used in practical CFD applications.
Unfortunately, the corrected scheme fails to converge in many of the curved surface tests, likely due to ill-conditioning of the
least squares procedure. Further investigation is needed to analyze the stability of the corrected scheme on curved grids with
high aspect ratio. Similar to the flat surface test, the meshes have wall cell aspect ratios of 102, 104, and 106. Fig. 11(a) shows
the solution error convergence in the inviscid terms for the node-centered and cell-centered schemes on right triangles of
varying aspect ratio. Very little difference can be seen between node- and cell-centered discretizations. This is also the case
for the inviscid error for quadrilaterals and right triangles of cell aspect ratio 106 using the cell-centered scheme, as shown in
Fig. 11(b). However, when considering the corresponding error in the viscous terms, shown in Fig. 11(c) and (d), the cell-cen-
tered scheme shows a marked advantage, similar to the previous flat surface case. Importantly, cell-centered quadrilaterals
 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686 7683

10-1 10-1

10-2 10-2
1 1
2 2
L2 error

L2 error
10-3 10-3
6
Node, RT, A=10 2 Cell, RT, A=10
Node, RT, A=10 4 Cell, Q, A=10 6
Node, RT, A=10 6
10-4 Cell, RT, A=10 2 10-4
Cell, RT, A=10 4
Cell, RT, A=10 6
-5 -5
10 10

500 1000 1500 500 1000 1500

-1 -1
ds ds

6
Node, RT, A=10 2 Cell, RT, A=10
0 Node, RT, A=10 4 Cell, Q, A=10 6
10 10
0
Node, RT, A=10 6
Cell, RT, A=10 2
Cell, RT, A=10 4
10-1 Cell, RT, A=10 6
10-1
L2 error

L2 error

-2
10 2 10
-2

2
1
1

10-3 10-3

-4
10 10
-4

500 1000 1500 500 1000 1500

ds-1 ds-1

Fig. 11. Effects of mesh stretching for curved surface meshes on node-centered, and cell-centered inviscid and viscous schemes using a scalar boundary
layer type MMS solution and various cell types. Cell types include right triangles and quadrilaterals (cell-centered only).

perform significantly better for the viscous discretization than do right triangles. This result would suggest that cell-centered
prismatic meshes are to be preferred for highly stretched meshes used to capture viscous boundary layer effects. Additional
studies have confirmed that second order accuracy is preserved in both theta and radial directions independently at the
viscous and inviscid limits.
Other studies [25,26] have reported excessive diffusion when using node-centered schemes on stretched triangular
meshes. It has been hypothesized that inaccuracies arise due to the misalignment of dual control volume faces with the
one dimensional approximate Riemann flux methods commonly used. These reports are consistent with the findings in this
work. The evidence points to the fact that node-centered schemes on stretched triangular meshes produce higher levels of
error than cell-centered schemes.

4.3. 3D prismatic grids

The results of the error convergence tests on stretched meshes indicate potential advantages of cell-centered formula-
tions over node-centered formulations for resolving boundary layers. Furthermore, stretched quadrilateral cells lead to lower
7684 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686

errors than stretched triangular cells, as shown in Fig. 10(b) and Fig. 11(d). In three dimensions, this would correspond to the
use of stretched prismatic cells near the body surface. Along with greater accuracy, prismatic cells fill space more efficiently
than tetrahedra. Since a triangular prism decomposes into three tetrahedra, and a quadrilateral prism (hexahedron) decom-
poses into five tetrahedra, triangular and quadrilateral prisms are three and five times more efficient at filling a given space,
respectively. This becomes even more significant considering that a large percentage of the total number of cells is usually
devoted to the boundary layer regions in CFD simulations.
Because of our interest in prismatic meshes in three dimensions, we present an order of accuracy study that concerns a
phenomenon particular to three dimensional grids – the occurence of non-planar control volume faces. We test the inviscid
scheme of Section 2.4 in three dimensions using a scalar MMS solution. Regular and perturbed prismatic grids in a unit cube
are used in a grid refinement study, shown in Fig. 12(a) and (b). Note that for the regular grids, control volume faces remain
planar, while for perturbed grids, the faces become non-planar.
Two methods of flux integration are used: single point quadrature, and a triangulated face quadrature. With single point
quadrature, the convective flux is evaluated at the face center of each quadrilateral face separating adjacent control volumes.
With the triangulated face quadrature, each quadrilateral face in the prismatic mesh is first triangulated, with one quadra-
ture point placed at the center of each of the two resulting triangular facets for a total of two quadrature points per quad-
rilateral face. For the artificial diffusion terms, the single point at the face center is used for both convective flux integration
methods, with no corresponding decrease in accuracy observed.
The results of this error convergence study using an isotropic scalar manufactured solution are shown in Fig. 12(c), and
summarized in Table 5. The results indicate that the simple single point quadrature is insufficient to retain second order
accuracy for three dimensional meshes with non-planar control volume surfaces. The triangulated face quadrature maintains
second order accuracy, even for perturbed grids with non-planar faces. This result is not widely appreciated in the literature,

-6
10

1
-7 1
10
L2 error

-8
10

-9
10
single pt. quad., reg 1
single pt. quad., pert.
tringulated face, reg. 2
triangulated face, pert.
10-10

20 40 60 80
-1
ds

Fig. 12. Isotropic grids used to assess the effect of non-planar faces in three dimensions, along with error convergence results for the inviscid terms only.
 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686 7685

Table 5
Order of solution error convergence in the inviscid terms using a scalar MMS solutio for planar and non-planar faces in three dimensions on cell-centered
prismatic grids.

Planar faces Non-planar faces

Single quadrature point 2 1
Triangulated face 2 2

although it has been pointed out by Delanaye and Liu [13] and Liu and Vinokur [27]. Further tests are needed to understand
the implications of non-planar faces for viscous flux integration.

5. Conclusions and future work

Results are presented which examine order of accuracy and error magnitude for node- and cell-centered schemes on
common mesh types and realistic solutions. The Method of Manufactured Solutions is found to be an effective method of
evaluating scheme performance as long as the MMS source term is discretized in a manner consistent with the base scheme.
Our studies indicate that, for isotropic grids, node-centered approaches produce slightly less error than cell-centered ap-
proaches for comparable cell size. For regular triangular meshes, the linear node-centered scheme is actually third order
accurate, while the cell-centered scheme is only second order accurate. Both schemes reduce to second order accuracy
for arbitrary triangles. Furthermore, a new corrected inviscid scheme is formulated which is third order accurate on trian-
gular meshes, even for arbitrary perturbations or with severe mesh stretching. One disadvantage of linear node-centered
scheme is the reduction to first order accuracy for quadrilateral grids. By extension, first order accuracy would be observed
whenever the scheme departs from the linear Galerkin approximation, such as for containment dual schemes on non-
Delaunay triangulations.
In contrast to the isotropic case, for stretched grids, we found that the cell-centered approximation produces less error
than the node-centered schemes, especially when considering the accuracy of the viscous terms on coarse meshes. The
cell-centered scheme produced the lowest error for quadrilateral meshes and the highest error for equilateral meshes.
Due to the good performance of the cell-centered scheme on quadrilateral meshes, we extended our studies to three-
dimensional prismatic meshes. In order to maintain second order accuracy for general prismatic meshes, non-planar control
volume faces were first triangulated. The triangular facets were then used to integrate the convective terms using one quad-
rature point per triangular facet. This resulted in two convective flux evaluations per quadrilateral face. However, it was suf-
ficient to only use one evaluation of the more expensive diffusive flux per quadrilateral face to maintain second order
accuracy.
Future work will focus on the effects of mesh stretching in the presence of grid curvature. Also, more three dimensional
results are needed to more fully understand the implications of non-planar faces for viscous flux integration.

Acknowledgements

Development was performed at the HPC Institute for Advanced Rotorcraft Modeling and Simulation (HIARMS) located at
the US Army Aeroflightdynamics Directorate at Moffett Field, CA, which is supported by the Department of Defense High
Performance Computing Modernization Office (HPCMO). Material presented in this paper is a product of the CREATE-AV Ele-
ment of the Computational Research and Engineering for Acquisition Tools and Environments (CREATE) Program sponsored
by the US Department of Defense HPC Modernization Program Office.

References

[1] D.J. Mavriplis, Unstructured mesh discretizations and solvers for computational aerodynamics, in: AIAA paper 2007-3955, AIAA 18th CFD Conference,
Miami, June 2007.
[2] P. Roache, Code verification by the method of manufactured solutions, Transactions of the ASME 124 (2002) 4–10.
[3] S. Pautz, Verification of Transport Codes by the Method of Manufactured Solutions: The ATTILA Experience, Technical Report, LA-UR-01-1487, Los
Alamos, 2001.
[4] C. Roy, Review of code and solution verification procedures for computational simulation, Journal of Computational Physics 205 (2005) 131–156.
[5] E. Luke, S. Hebert, D. Thompson, Theoretical and practical evaluation of solver-specific mesh quality, in: AIAA paper 2008-0934, AIAA 46th ASM, Reno,
January 2008.
[6] H. Sun, Darmofal, R. Haimes, On the impact of triangle shapes for boundary layer problems using high-order finite element discretization, in: AIAA
paper 2010-0542, AIAA 48th ASM, Orlando, January 2010.
[7] B. Despres, Lax theorem and finite volume schemes, Mathematics of Computation 73 (2003) 1203–1234.
[8] B. Wendroff, Supraconvergence in Two Dimensions, Technical Report, LA-UR-95-3068, Los Alamos, 1995.
[9] J. Thomas, B. Diskin, C. Rumsey, Towards verification of unstructured-grid solvers, in: AIAA paper 2008-666, AIAA 46th ASM, Reno, January 2008.
[10] B. Diskin, J. Thomas, Comparison of node-centered and cell-centered unstructured finite-volume discretizations: inviscid fluxes, in: AIAA paper 2010-
1079, AIAA 48th ASM, Orlando, January 2010.
[11] B. Diskin, J. Thomas, E. Nielsen, H. Nishikawa, Comparison of node-centered and cell-centered unstructured finite-volume discretizations. Part 1:
Viscous Fluxes, in: AIAA paper 2009-0597, AIAA 47th ASM, Orlando, January 2009.
[12] B. Diskin, J. Thomas, Accuracy analysis for mixed-element finite-volume discretization schemes, in: NIA Report 2007-08, National Institute of
Aerospace, 2007.
7686 A. Katz, V. Sankaran / Journal of Computational Physics 230 (2011) 7670–7686

[13] M. Delanaye, Y. Liu, Quadratic reconstruction finite volume schemes on 3D arbitrary unstructured polyhedral grids, in: AIAA paper 1995-3259, AIAA
14th CFD Conference, Norfolk, June 1999.
[14] M. Svard, J. Gong, J. Nordstrom, An accuracy evaluation of unstructured node-centered finite volume methods, in: NIA Report 2005-04, 2005.
[15] T.J. Barth, Numerical aspects of computing viscous high Reynolds number flows on unstructured meshes, in: AIAA paper 1991-0721, AIAA 29th ASM,
Reno, January 1991.
[16] A. Jameson, Analysis and design of numerical schemes for gas dynamics 1 artificial diffusion, upwind biasing, limiters and their effect on accuracy and
multigrid convergence, International Journal of Computational Fluid Dynamics 4 (1995) 171–218.
[17] A. Jameson, Analysis and design of numerical schemes for gas dynamics 2 artificial diffusion and discrete shock structure, International Journal of
Computational Fluid Dynamics 5 (1995) 1–38.
[18] A. Jameson, J. Vassberg, A Vertex-Centroid (V-C) Scheme for the Gas-Dynamics Equations, Technical Report, International Conference on CFD, Kyoto,
Japan, July 2000.
[19] N. Frink, Upwind scheme for solving the euler equations on unstructured tetrahedral meshes, AIAA Journal 30 (1) (1992) 70–77.
[20] J. Batina, Implicit flux-split euler schemes for unsteady aerodynamic analysis involving unstructured dynamic meshes, in: AIAA paper 1990-0936, AIAA
31st Structural Dynamics and Materials Conference, Long Beach, CA, April 1990.
[21] P. Roache, S. Steinberg, Symbolic manipulation and computational fluid dynamics, AIAA Journal 22 (1984) 1390–1394.
[22] M. Giles, Accuracy of node-based solutions on irregular meshes, Lecture Notes in Physics 323 (1989) 273–277.
[23] W. Oberkampf, C. Roy, Verification and Validation in Scientific Computing, Cambridge, 2010.
[24] A.H. Shapiro, The Dynamics and Thermodynamics of Compressible Fluid Flow, vol. 2, The Ronald Press Company, 1954.
[25] M. Aftosmis, D. Gaitonde, T. Tavares, Behavior of linear reconstruction techniques on unstructured meshes, AIAA Journal 33 (1995) 2038–2049.
[26] H. Luo, D. Sharov, J. Baum, R. Löhner, On the computation of compressible turbulent flows on unstructured grids, in: AIAA paper 2000-0926, AIAA 38th
ASM, Reno, January 2000.
[27] Y. Liu, M. Vinokur, Exact integrations of polynomials and symmetric quadrature formulas over arbitrary polyhedral grids, Journal of Computational
Physics 140 (1998) 122–147.