It C
Resume:
J Sci Comput (****) **: *** ***
A Study of Viscous Flux Formulations for a p-Multigrid
Spectral Volume Navier Stokes Solver
R. Kannan Z.J. Wang
Received: 19 July 2008 / Accepted: 7 January 2009 / Published online: 22 January 2009
Springer Science+Business Media, LLC 2009
Abstract In this paper, we improve the Navier Stokes ow solver developed by Sun et al.
based on the spectral volume method (SV) in the following two aspects: the development
of a more ef cient implicit/p-multigrid solution approach, and the use of a new viscous ux
formula. An implicit preconditioned LU-SGS p-multigrid method developed for the spectral
difference (SD) Euler solver by Liang is adopted here. In the original SV solver, the viscous
ux was computed with a local discontinuous Galerkin (LDG) type approach. In this study,
an interior penalty approach is developed and tested for both the Laplace and Navier Stokes
equations. In addition, the second method of Bassi and Rebay (also known as BR2 approach)
is also implemented in the SV context, and also tested. Their convergence properties are
studied with the implicit BLU-SGS approach. Fourier analysis revealed some interesting
advantages for the penalty method over the LDG method. A convergence speedup of up
to 2-3 orders is obtained with the implicit method. The convergence was further enhanced
by employing a p-multigrid algorithm. Numerical simulations were performed using all the
three viscous ux formulations and were compared with existing high order simulations (or
in some cases, analytical solutions). The penalty and the BR2 approaches displayed higher
accuracy than the LDG approach. In general, the numerical results are very promising and
indicate that the approach has a great potential for 3D ow problems.
Keywords Spectral volume LDG Penalty BR2 Implicit LU-SGS High-order
p-multigrid
1 Introduction
The spectral volume (SV) method was originally developed by Wang, Liu et al. for hyper-
bolic conservation laws on unstructured grids [26, 50 54]. The spectral volume method is a
subset of the Godunov type nite volume method, which is the starting block for the devel-
opment of a plethora of methods such as the k -exact nite volume [4, 16], MUSCL [48, 49],
R. Kannan Z.J. Wang
Department of Aerospace Engineering, Iowa State University, 2271 Howe Hall, Ames, IA 50011, USA
e-mail: abpwon@r.postjobfree.com
166 J Sci Comput (2009) 41: 165 199
and the essentially non-oscillatory (ENO) [1, 21] methods. The spectral volume method can
be viewed as an extension of the Godunov method to higher order by adding more degrees-
of-freedom (DOFs) in the form of subcells in each cell (simplex). The simplex is referred to
as a spectral volume (SV) and the subcells are referred to as control volumes (CV). Every
simplex in the SV method utilizes a structured pattern to construct the subcells (CVs).
As in the nite volume method, the unknowns (or DOFs) are the subcell-averaged solu-
tions. A nite volume procedure is employed to update the DOFs. The spectral volume
method shares many similar properties with the discontinuous Galerkin (DG) [14, 15] and
the spectral difference (SD) [24, 43] methods, such as discontinuous solution space and
compactness. They mainly differ on how the DOFs are updated. The DG being a derivative
of the nite element method has unknown the elemental nodal values as DOF. The spectral
volume being a derivative of the nite volume has subcell averages as its DOF while the
spectral difference has point wise values as DOF. In terms of complexity, DG requires both
volume and surface integrations. In contrast, SV requires only surface integrations and the
SD requires differentiations.
The SV method was successfully implemented for 2D Euler [53] and 3D Maxwell equa-
tions [26]. Recently, Chen [12, 13] developed many high-order SV partitions for simplexes
in 2D and 3D with relatively small Lebesgue constants. Comparisons between the SV and
DG methods were made by Sun and Wang [44] and by Zhang and Shu [56]. The SV method
was also applied to solve the 3D Euler and Navier Stokes equations by Haga et al. [20] on
Japan s Earth Simulator. More recently, Van den Abeele et al. [46, 47] performed Fourier
analysis for both 1D and 2D SV methods and identi ed a weak instability in several SV
partitions. New partitions are derived which showed improved stability properties.
The main goal of this research study is to extend the SV method for solving the Navier
Stokes equations. The proper discretization of viscous uxes has been studied extensively in
the DG methods and shown to be very important for accuracy and stability. In the classical
second-order nite volume method, the solution gradients at an interface are usually com-
puted by averaging the gradients of the neighbour cells sharing the face in consideration.
However, for higher order elements, special care has to be taken in computing the solu-
tion gradients. In the late 1970s and early 1980s, Arnold [3] and Wheeler [53] introduced
the discontinuous nite element methods known as penalty methods for elliptic operators.
More recently, many researchers [5, 6, 9, 10, 14, 15] have applied DG methods to diffusive
operators. One procedure was the local discontinuous Galerkin (LDG) method, developed
by Cockburn and Shu [14, 15]. This method dealt with rewriting a second-order equation
as a rst-order system and then discretize the rst-order system using the DG formulation.
Their simplicity and effectiveness have made them the main choice for discretizing the vis-
cous uxes. Baumann and Oden [9] came up with different DG methods for the second
order viscous terms. In 1997, Bassi and Rebay [5] came up with a symmetrical scheme to
discretize the viscous uxes. In 2000, the above researchers came up with a second con-
cept (also referred to as BR2 [7]). More details on the above methods will be discussed in
this paper. Recently, Zhang and Shu [54] conducted some Fourier analysis for the different
formulations.
Recently Sun et al. [45] implemented the SV method for the Navier Stokes equations
using the LDG [14] approach to discretize the viscous uxes. The LDG approach alternates
between the right and the left directions for obtaining the CV averaged gradient and the
residual respectively. More recent numerical tests indicated that the computational results
using LDG are somewhat dependent on how the faces are oriented especially for unstruc-
tured and non uniform grids.
In this paper, we study and implement two more viscous ux formulations in the SV
method setting. The rst is the penalty method and the second is the 2nd method of Bassi
J Sci Comput (2009) 41: 165 199 167
and Rebay (i.e. BR2). Fourier analysis was performed on all methods and yielded some
interesting results on accuracy and time step limit (for explicit simulations) and speed of
convergence.
The convergence properties of high order methods like DG, SV, SD is generally poor
with explicit time marching approaches. In this paper, we adopt a preconditioned implicit
LU-SGS approach. The original LU-SGS approach was developed by Yoon and Jameson
on structured grids [21], and later improved upon by Jameson and Caughley [23], Chen
and Wang [11] using an element Jacobi preconditioner. It was also used by Liang et al.
[24] to solve the Euler equations with a p-multigrid algorithm for the SD method. In addi-
tion, we also use a p-multigrid algorithm to enhance the convergence rate. The p-multigrid
method employs smoothing operators, which hasten the convergence by using the coarse
levels constructed by reducing the order of the interpolation polynomial. The p-multigrid
was initially proposed by Ronquist and Patera [35], and extended by Maday and Munoz
[28]. The acceleration of the convergence by the p-multigrid algorithm on Euler equations
was demonstrated by Bassi and Rebay [8], Nastase and Mavriplis [31] and Luo et al. [27].
Helenbrook [22] examined the performance of the Laplace equation and the convection
equations in two dimensions. Fidkowski et al. [17] developed a p-multigrid approach using
an effective line-implicit smoother for their high-order DG Navier Stokes solver. Luo et al.
[27] demonstrated a low storage p-multigrid approach for a 3D DG Euler solver, in which
an explicit Runge Kutta scheme is used for the highest order discretization, while implicit
SGS smoothers are employed for the lower-order operators. The above concept was put to
use for the SD method by Liang et al. [24]. Recently, Nastase and Mavriplis [31] devel-
oped an hp-multigrid approach for their high-order inviscid DG solver, and demonstrated
h-independent and p-independent convergence rates.
The paper is organized as follows. In the next section, we review the basics of the SV
method. The new penalty and the BR2 formulations are presented in Sect. 3. A detailed
linear analysis is performed for the three viscous uxes in Sect. 4. Time integration schemes
and the p-multigrid method are discussed in Sect. 5. Section 6 presents with the different
test cases conducted in this study. Finally conclusions from this study are summarized in
Sect. 7.
2 Basics of the Spectral Volume Method
Consider the general conservation equation
Q (fi (Q) fv (Q)) (gi (Q) gv (Q))
+ + = 0, (2.1)
t x y
in domain with appropriate initial and boundary conditions. In (2.1), x and y are the
Cartesian coordinates and (x, y), t [0, T ] denotes time, Q is the vector of conserved
variables, and fi and gi are the inviscid uxes in the x and y directions, respectively. fv and
gv are the viscous uxes in the x and y directions, respectively. Domain is discretized
into I non-overlapping triangular (2D) cells. In the SV method, the simplex grid cells are
called SVs, denoted Si, which are further partitioned into CVs, denoted Cij, which depend
on the degree of the polynomial reconstruction. Figure 1 shows linear, quadratic and cubic
partitions in 2D.
We need N unknown control volume solution averages (or DOFs) to construct a degree
k polynomial. N is calculated using the below formula (in 2D)
(k + 1)(k + 2)
N= (2.2),
2
168 J Sci Comput (2009) 41: 165 199
Fig. 1 Partitions of a triangular SV. Linear, quadratic and cubic reconstructions are shown in (a), (b) and (c)
respectively
where k is the degrees of the polynomial, constructed using the CV solution averages. The
CV averaged conserved variable for Cij is de ned as
1
Qi,j = j = 1, . . ., N, i = 1, . . ., I, (2.3)
QdV,
Vi,j Ci,j
where Vi,j is the volume of Cij . Given the CV averaged conserved variables, a degree k
polynomial can be constructed such that it is (k + 1)th order approximation to Q. In other
words, we can write the polynomial as
N
pi (x, y) = (2.4)
Lj (x, y)Qi,j,
j =1
where the shape functions Lj (x, y) satisfy
1
Ln (x, y)dV = j,n . (2.5)
Vi,j Ci,j
Equation (2.1) is integrated over the Cij . This results in the following equation
K
1
Q
+ (F .n)dA = 0, (2.6)
t Vi,j Ar
r =1
where F = (fi fv, gi gv ) where Ar represents the r th face of Cij, n is the outward
unit normal vector of Ar and K is the number of faces in Cij . The surface integration on
each face is done using a (k + 1)th order accurate Gauss quadrature formula. The uxes are
discontinuous across the SV interfaces. The inviscid uxes are handled using a numerical
Riemann ux such as the Rusanov ux [36], the Roe ux [34] or AUSM ux [25]. The
handling of the viscous uxes is discussed in the next section.
3 Spectral Volume Formulation for the Diffusion Equation
The following 2D diffusion equation is considered rst in domain with appropriate initial
and boundary conditions
u
( u) = 0, (3.1)
t
J Sci Comput (2009) 41: 165 199 169
where is a positive diffusion coef cient. We de ne an auxiliary variable
q = u. (3.2)
Equation (3.1) then becomes
u
( q ) = 0. (3.3)
t
Using the Gauss-divergence theorem, we obtain
K
q ij Vij = u ndA, (3.4)
Ar
r =1
K
d uij
Vij q ndA = 0, (3.5)
dt Ar
r =1
where q ij and uij are the CV averaged gradient and solution in Cij . As the solution u is
cell-wise continuous, u and q at SV boundaries are replaced by numerical uxes q and u.
The above equations thus become
K
q ij Vij = u ndA, (3.6)
Ar
r =1
K
d uij
Vij q ndA = 0. (3.7)
dt r =1
3.1 LDG Approach
The commonly used approach for obtaining the numerical uxes is the LDG approach. In
this approach, the numerical uxes are de ned by alternating the direction in the following
manner [45]
u = uL, (3.8)
q = qR, (3.9)
where uL and uR are the left and right state solutions of the CV face in consideration and qL
and qR are the left and right state solution gradients of the face (of the CV) in consideration.
Thus if the CV face lies on the SV boundary, uL = uR and qL = qR .
3.2 Penalty Approach
It can be seen (from (3.8) and (3.9)) that LDG is inherently unsymmetric. A symmetric
approach was given by Bassi and Rebay [5], in which the numerical uxes are de ned by
u = 0.5 (uR + uL ), (3.10)
q = 0.5 (qR + qL ). (3.11)
170 J Sci Comput (2009) 41: 165 199
Fig. 2 Stencils in 1D. Case (a): compact stencil; Case (b): non compact stencil
Analysis by Brezzi et al. [10] showed that the approach might be unstable in some situations.
In this paper, we suggest the following the penalty approach to obtain the numerical uxes:
u = 0.5 (uR + uL ), (3.12)
Ar
q = 0.5 (qR + q L ) + (uR uL )n (3.13),
Vij
where qL and qR are the left and right state solution gradients of the face (of the CV) in
consideration, Ar is the area of the face (of the CV) in consideration, Vij is the CV volume.
One can see a similarity between the above equation and an approximate Riemann (like
Roe, Rusanov or AUSM) ux. The approximate Riemann ux is obtained by averaging the
left and right state uxes and then adding a dissipation term. This dissipation term is
1. Proportional to the jump in the solution between the right and left states.
2. Proportional to the Jacobian term/matrix or its eigen values (the Jacobian term in 1D
f
is Q ). For instance, in Rusanov ux, it is the maximum eigen value of the Jacobian
matrix.
Equation (3.13) is obtained by averaging the left and right states and then penalizing it
with the penalty term. This is similar to the structure of the approximate Riemann ux. The
Ar
Jacobian term in this case has a dimension of 1/length. So we picked Vij as an approximation
to the eigen value. The penalty term has a sign that is opposite to the dissipation term. This
is because the dissipation terms come on the RHS.
It can be seen that the current formulation is still non compact. By de nition, a non
compact system is one wherein the residual of a cell is dependent on the solution of cells
which are neither the current cell nor its adjacent face-neighbors. Compactness is generally
a desired attribute as it has serious desirable properties:
1. Very desirable when the cells are clustered. This is because the reconstructions are done
from information obtained from a restricted region and hence is more accurate.
2. Ef cient for parallel applications: lesser information is needed w.r.t. storage and message
transfers
Figure 2 shows compact and non compact stencils in 1D. The red and blue colors indicate
the domain of dependence.
3.3 Second Approach of Bassi and Rebay (BR2)
The penalty approach is de nitely more symmetrical than LDG. However, the implementa-
tion results in a non-compact stencil. We now extend the compact BR2 formulation to the
SV method. The BR2 formulation is slightly different from the above two methods. Two
different gradients are evaluated for the residual computation:
I. The rst gradient is referred to as the self gradient. It is computed using the solution from
the cell in consideration. In other words, there is no contribution from the neighboring
cells. This gradient is utilized for computing the viscous uxes through the CV faces,
which lie on the SV boundary.
J Sci Comput (2009) 41: 165 199 171
Fig. 3 Two SVs having a
common SV face f
II. The second gradient is referred to as the averaged gradient. It is computed by using the
averaged solution obtained from the left and the right sides of the CV face. This gradient
is used for computing the viscous uxes through the CV faces which do not lie on the
SV boundary.
BR2 uses the concept of lifting functions associated with each of the three faces of a
SV. These lifting functions, r + and r (which are actually polynomials) are then used to
correct the gradients at the face. Every SV face has the above lifting functions, r + and
r associated with it. CV averaged lifting functions are rst computed for the SV face in
consideration. The actual procedure can be explained using an example. Consider the two
2nd order SVs, shown in Fig. 3. The SV boundary face in consideration is f and is colored
in red. This comprises of two CV faces, f 1 and f 2. The CVs of the SV in consideration
+
are also marked in the gure. The CV averaged lifting functions r1 and r1 for CV1 can be
obtained using (3.14) and (3.15)
(u u+ ) +
+
V + r1 = (3.14)
n dA,
2
f1
(u+ u )
V r1 = (3.15)
n dA,
2
f1
where + and refer to the current CV and the neighboring CV respectively, V + and V
are the volumes of the current and neighboring CVs respectively, n+ and n are the outward
normal w.r.t. SVs (V + and V respectively) for the face in consideration. Note: n+ = n .
+
Similarly, one could obtain the CV averaged lifting functions r2 and r2 for CV2. The
CV averaged lifting functions for the CV3 is zero because it does not have a sub-face on
++ +
face f . Given the CV averaged lifting functions, r1, r2 and r3 (which is zero!), a linear
polynomial can be constructed. This is the actual lifting functions r + corresponding to the
faces f 1 and f 2.
The above procedure is repeated for the other SV faces so as to obtain the lifting functions
corresponding to the CV faces located in the other two SV faces. The other lifting function
r can be obtained in a similar way. The gradients used for the CV faces located on the SV
boundary are given by
1
q = (qR + r + + qL + r ). (3.16)
2
172 J Sci Comput (2009) 41: 165 199
It must be noted that BR2 requires nearly twice the number of CV averaged gradient
computations and their reconstructions. In addition, the computations involved in evaluating
the lifting function have more oating point operations (FLOPs) than the penalty or the LDG
scheme. As expected, additional storage is also required for the BR2 scheme.
4 Time Integration Algorithms
With the high-order SV method, it is clearly possible to achieve lower error levels than the
traditional rst order or second order methods. However, most of the developments of the SV
methods employed explicit TVD Runge Kutta schemes [50 53] for the temporal advance-
ment and those numerical experiments are all performed on single grids. These explicit
temporal discretizations can be implemented with ease. In addition, they require limited
memory and can be vectorized and parallelized. However it is well known that the explicit
scheme is limited by the Courant Friedrichs Lewy (CFL) condition. When the polynomial
order of the SD method is high or the grid is stretched due to complex geometries, the con-
vergence rate of the explicit scheme slows down rapidly. Solution strategies to remedy this
problem include the implicit method [24] and the multigrid method [24, 27].
In general, the high-order spatial operators are much stiffer than lower-order ones. In
addition, for time accurate problems, the maximum allowable CFL number decreases with
increasing the order of accuracy for explicit schemes. This problem becomes even more
severe for viscous ow problems with highly clustered meshes: explicit high-order methods
are adversely limited by time step size. Hence the development of ef cient time integration
schemes for high-orders is essential. Some of the existing methods are reviewed in the next
two paragraphs.
Explicit Runge Kutta schemes have been one of the most widely used methods of time
integration. One particular class of the above is the strong-stability-preserving (SSP) Runge
Kutta schemes. This class of time integration method was originally developed by Shu
[39, 40], and Shu and Osher [41] and named TVD Runge Kutta schemes. It was also the
topic of research for many other researchers e.g. [19, 42]. This class of Runge Kutta schemes
preserves the stability properties of forward Euler in any norm or semi-norm. They have
been popular for high-order spatial operators because of its TVD or SSP property. The co-
ef cients of these schemes are not unique. However, optimal versions with maximum CFL
numbers have been derived in [39, 40] for the second and third-order schemes, and by Spi-
teri and Ruuth [42] for the fourth-order counterpart. More details can be found in a recent
review article by Gottlieb [18].
The limitations of the explicit methods were mentioned in the early part of this section.
It is obvious that implicit algorithms are necessary to overcome the time step limit suffered
by explicit algorithms especially for highly clustered viscous meshes. Many types of im-
plicit algorithms have been successfully developed for unstructured grid-based solvers in
the last two decades, e.g., the element Jacobi, Gauss Seidel, precondition GMRES [5, 37],
matrix-free Krylov [33], lower upper symmetry Gauss Seidel [2, 38], and line-implicit al-
gorithms [29]. In addition, these implicit algorithms can serve as smoothers for geometric
or p-multigrid approaches. Many of these algorithms have been successfully applied to high-
order spatial discretizations (e.g. [27, 33] for DG and [24] for SD). In almost all implicit
approaches, the non-linear system of equations is linearized and then solved with an itera-
tive solution algorithm. Even though these implicit methods offer extremely high speedup,
they need huge memory to store the associated matrices. One intelligent way to mitigate the
J Sci Comput (2009) 41: 165 199 173
above problem is to use the traditional multi stage Runge Kutta method as the higher level
smoother and the implicit scheme for the lower levels [27].
The following unsteady equation is considered in the current paper and the time integra-
tion schemes (either explicit or implicit schemes) can be used as a smoother.
Q
= Ric (Q), (4.1)
t
where Ric (Q) denotes the residual vector of the current cell and Q is the solution vector.
The remainder of the chapter gives the formulations for the three stage SSP Runge Kutta
schemes and the implicit LU-SGS method.
4.1 The Explicit Method
We mainly concentrate on the three stage SSP Runge Kutta scheme. The three-stage explicit
SSP Runge Kutta [39, 40] can be written as follows:
u(1) = un tRi (un );
i i
3n 1 (1)
u(2) =
i u+
[u tRi (u(1) )];
(4.2)
4i 4i
1 2 (2)
un+1
i = un +
[u tRi (u(2) )].
3i 3i
The main advantage of the Runge Kutta scheme is that it requires little memory for
storage. In addition, this method is inherently simple and so is easy to program. These are
the main reasons for it being one of the most preferred methods of time integration.
The main bottleneck associated with the Runge Kutta scheme is the limit imposed on the
time step. Euler (and Navier Stokes) equations for realistic geometries entail a rather strict
limit on the time step. Even though the above can be circumvented by using a very higher
order (several RK steps) scheme, it is seldom used as it required lots of storage and thus
adversely affects its simplistic nature. This prompts us to switch over to implicit schemes.
4.2 The Implicit Method
It is well-known that the explicit scheme is limited by the Courant Friedrichs Lewy (CFL)
condition. When the polynomial order of the SD method is high for or the grid is stretched
due to complex geometries, the convergence rate of the explicit scheme slows down rapidly.
The implicit methods are normally formulated by the linearization of a given set of equa-
tions. At each current cell c, the backward Euler difference can be written as
c c
Qn+1 Qn
[Rc (Qn+1 ) Rc (Qn )] = Rc (Qn ), (4.3)
t
where n refers to the current time level and n + 1 the next time level.
c c
Let Qc = Qn+1 Qn and linearizing the residual, we obtain
Rc Rc
Rc (Qn+1 ) Rc (Qn ) Qc + (4.4)
Qnb,
c Qnb
Q nb=c
174 J Sci Comput (2009) 41: 165 199
where nb indicates all the neighboring cells contributing to the residual of cell c. Therefore,
the fully linearized equations for (4.3) can be written as
I Rc Rc
Qc Qnb = Rc (Qn ). (4.5)
t Qnb
Qc nb=c
However, it costs too much memory to store the entire LHS implicit Jacobian matrices. For
instance, a 3rd order Euler system of equations in 2D would result in a [96 96] matrix per
cell. In addition, the cost of performing a Gaussian Elimination of that matrix is O(963 ) i.e.
O(900000) units. Therefore, we employ a LU-SGS scheme to solve (4.5), i.e., we use the
most recent solution for the neighboring cells,
I Rc Rc
c nb
Q(k+1) = Rc (Qn ) + Q .
(4.6)
t Qnb
Qc nb=c
The matrix
I Rc
D= (4.7),
t Qc
is the element (or cell) matrix, which is not quite small for 2nd to 3rd order SV schemes. For
instance, D is [24 24] for the 3rd order SV Euler/Navier Stokes system of equations. The
cost of performing a Gaussian Elimination for a matrix of this size is O(243 ) i.e. O(14000)
units. These matrices also consume much lesser storage space.
Rc
Since we do not want to store the matrices Qnb, (4.6) is further manipulated as follows
Rc Rc
nb c nb nb
Q = Rc (Qn, {Qn }) + Q
Rc (Qn ) +
Qnb Qnb
nb=c nb=c
Rc
c nb c nb
Rc (Qn, {Q }) Rc (Q, {Q }) Q
c
Qc
Rc
= Rc (Q ) Q . (4.8)
c
Qc
Q(k+1) = Q(k+1) Q . We can combine (4.6) and (4.8) together to obtain
2
Let c c c
Q
I Rc c
Q(k+1) = Rc (Q )
c
2
(4.9)
.
t t
Qc
Equation (4.9) is then solved with a direct LU decomposition solver and sweeping sym-
metrically forward and backward through the computational grid. Note that once (4.9) is
solved to machine zero, the unsteady residual is zero at each time step. For steady ow
Q
problem, we found that the term t c in the right-hand-side of (4.9) can be neglected and
leading to quicker convergence. Note that Q is the difference between the current solu-
c
tion Q and the solution at the previous time level Qn . In reality, the entire system is swept
c c
several times in order to proceed to the next time level. As a result, Q is in uenced by
c
the solution occurred several sweeps ago. This introduces an under-relaxation effect. Hence
Q
neglecting the t c term accelerates the convergence. We de ne the solver obtained using
Q c
(4.9) as implicit normal approach. If term is dropped, the iterative solver is de ned as
t
implicit simpli ed approach.
J Sci Comput (2009) 41: 165 199 175
5 Fourier Analysis for the Diffusion Equation
5.1 One Dimensional Analysis
In this analysis, we follow a technique described by Zhang and Shu [55] and focus on linear,
quadratic and cubic reconstructions. The SV is partitioned into two equal CVs for the sec-
ond order simulations. The CVs for the third and the fourth order are clustered toward the
SV boundaries. The locations of the CV faces (i.e. nodes in 1D) were based on the Gauss
quadrature points. For the sake of simplicity, let us rst consider a linear partition shown in
Fig. 4. In this case, all the formulations can be cast in the following form:
uj 1,1 uj +1,1 uj +2,1
uj,1 uj 2,1 uj,1
d
=A +B +C +D +E, (5.1)
uj 1,2 uj +1,2 uj +2,2
dt uj,2 uj 2,2 uj,2
where A, B, C, D and E are constant matrices. We now seek a general solution of the
following form
u(x, t) = uk (t)eikx,
(5.2)
where k is the index of modes (k = 1, 2, . . .) representing the wave number. Obviously, the
2
analytical solution for (3.1) is u(x, t) = eikx k t . The solution we are looking for can be
expressed as
uj,1 uk,1
= eikxj,3/2 . (5.3)
uj,2 uk,2
Substituting (5.3) into (5.1), we obtain the advancement equation:
uk,1 uk,1
= G(k, h) (5.4),
uk,2 uk,2
where the ampli cation matrix is given by
G = e 2ikh A + e ikh B + C + eikh D + e2ikh E. (5.5)
The above method can be easily extended to 3rd and 4th orders. In general, all but one
of the eigen values of G is made up of spurious modes and is damped rapidly. The error
associated with the scheme, convergence properties can be determined by analyzing the non
spurious mode.
Spatial Analysis Figure 5a shows the variation of the principal eigen value with respect
to the non dimensional frequency (= kh) for the 2nd order case. This is also referred to
as the Fourier footprint. Figure 5b shows the error at lower wave numbers. The error is the
difference between the principal (physical mode) eigen value and the analytical eigen value
(i.e. 2 ).
It is clear that the errors associated with the penalty formulation are lower than that of
LDG. This means that the system is damped faster in the penalty case and hence converges
faster. The situation is similar for the 3rd order as is seen from Figs. 6a and 6b. Thus in
Fig. 4 Linear spectral volume
in 1D
176 J Sci Comput (2009) 41: 165 199
Fig. 5 Second order Fourier analysis. Case (a): Fourier footprint; Case (b): error associated at low wavenum-
bers
Fig. 6 Third order Fourier analysis. Case (a): Fourier footprint; Case (b): error associated at low wavenum-
bers
general, we expect the penalty scheme to converge faster than the LDG for the 2nd and 3rd
orders. The principal eigen value lies on top of the y = x 2 parabola for the fourth order
(Fig. 7a). The fourth order errors are too negligible (Fig. 7b) to comment on convergence
phenomena. It must be noted that the matrices A and E are zero matrices for the LDG and
BR2 formulations.
Temporal Analysis In this section, we compute the time step requirements of the 3-stage
Runge Kutta scheme. We use the 3 stage Runge Kutta scheme to solve and study (4.4). We
employ the standard Taylor expansion to link [ukn ] (the CV averaged solution vector at the
current time step) to [ukn+1 ] (the CV averaged solution vector at the next time step):
[ukn ]
t 2 2 [ukn ]
t 3 3 [ukn ]
[ukn+1 ] = [ukn ] +
+ + (5.6)
t .
2 t 3
2 6
t t
J Sci Comput (2009) 41: 165 199 177
Fig. 7 Fourth order Fourier analysis. Case (a): Fourier footprint; Case (b): error associated at low wavenum-
bers
Table 1 Maximum
Method 2nd order 3rd order 4th order
non-dimensional time step
( = dt2 ) for obtaining stable
dx LDG 0.157 0.0266 0.0109
solutions for the LDG, penalty
Penalty 0.182 0.0322 0.0133
and BR2 methods in 1D
BR2 0.314 0.0784 0.0321
Combining (5.6) and (5.4), we obtain the following advancement equation:
[ukn+1 ] = [H ][ukn ],
(5.7)
where [H ] is
[G]2 2 [G]3 3
[H ] = [I ] + [G] t + t+ (5.8)
t,
2 6
with [I ] being the Identity matrix.
The eigen values of [H ] matrix were computed. The moduli of each of the above com-
puted eigen values need to be less than unity to ensure a stable system. Table 1 lists the
maximum non-dimensional time step ( = xt2 ) required for obtaining a stable solution. It
can be seen that the penalty and BR2 methods permit a higher time-step limit than the LDG
method. In fact, the time step permitted by the BR2 method is more than double of that of
the penalty method. Numerical experiments veri ed the above conclusion.
5.2 Two Dimensional Analysis
In this analysis, we follow a technique described by Van Den Abeele et al. [47] and focus
on the linear reconstruction. The SV is partitioned into 3 CVs. We had to use a basic unit,
for imposing periodicity. In a one-dimensional case, the basic unit was the spectral volume
itself. It comprises of 2 SVs in 2D. The basic unit (only SVs are shown) is shown in Fig. 8.
This solution process would result in solving a [6 6] set of equations. The rate of change
of the solution (in the basic unit) can be written as a linear combination of solutions from
178 J Sci Comput (2009) 41: 165 199
Fig. 8 Basic Building block
represented by bold lines; the
neighbors (represented by normal
lines) are used in analysis
the left (L), right (R ), top (T ), bottom (B ) and central (C ) units. The formulations can be
cast in the following form:
d
[ ui,j ] = L[ ui 1,j ] + R [ ui +1,j ] + C [ ui,j ] + B [ ui,j 1 ] + T [ ui,j +1 ], (5.9)
dt
where L, R, C, B and T are constant matrices and [ ui,j ] is the solution vector of the basic
building unit. We now seek general solution of the following form
u(x, t) = uk (t)eik(x cos +y sin,
(5.10)
where k is the index of modes (k = 1, 2, . . .) representing the wave number and is the di-
2
rection of wave propagation. Obviously, the analytical solution for (3.1) is u(x, t) = eikx k t .
We follow the method described in the 1D analysis section and obtain the following result:
[uk ] = G(k, h)[uk ],
(5.11)
where the ampli cation matrix is given by
G = e ikh cos L + eikh cos R + C + e ikh sin B + T eikh sin . (5.12)
Spatial Analysis There is an extra degree of freedom, . Due to symmetry of the system
and of the building unit, we need to focus on varying from 0 to /4 radians. In this paper,
we focus on these extreme cases. Figures 9a and 9b show the variation of the principal eigen
value with respect to the non dimensional frequency (= kh) for = 0 and /4 radians
respectively. Figures 10a and 10b show the errors at low frequencies. In general, the trends
are qualitatively similar to those obtained using 1D analysis. In addition, it can be seen that
the errors are much lower when = /4 radians. This is probably due to the way our SVs
are de ned (right triangles).
Temporal Analysis As in the 1D analysis, we stick to the 3-stage SSP Runge Kutta scheme.
Table 2 lists the maximum non dimensional time step required for obtaining a stable solu-
tion. The time steps mentioned in Table 2 are the same for both = 0 and /4 radians. Once
again, the penalty and BR2 methods permit a higher time-step limit than the LDG method.
The time step permitted by the BR2 method is more than that of the penalty method.
5.3 Discussion on Stability and Spurious Modes
In the previous sections, we analyzed the stability requirements of the 3 stage SSP Runge
Kutta smoother. It has been a common practice in CFD to use the maximum possible time
J Sci Comput (2009) 41: 165 199 179
Fig. 9 Second order Fourier footprint in 2D. Case (a): = 0 radians; Case (b): = /4 radians
Fig. 10 Error associated at low wavenumbers for the 2D 2nd order scheme. Case (a): = 0 radians; Case (b):
= /4 radians
Table 2 Maximum
LDG 0.0267
non-dimensional time step for
Penalty 0.0465
obtaining stable solutions for the
2nd order, 2D LDG, penalty and BR2 0.0635
BR2 methods
step to drive a system to convergence. This time step is limited by the stability requirements
(as determined in the earlier sections). A different trend was observed during the course
of this study. Consider the 2nd order LDG case: the ampli cation matrix G has the eigen
values 16/ h2 (spurious mode) and 2(1 cos / h2 (physical mode). Thus if
uj,+1
n n
uj,1
1
= H (, dt) (5.13)
.
uj,+1 n
n
uj,2
2
Obviously for stability, the eigen values of the growth matrix H (, dt ) need to lie be-
tween 1 and 1. Figure 11a shows the eigen values for this system for a non dimensional
time step of 0.157. The system satis es the requirement for stability. However, the spurious
180 J Sci Comput (2009) 41: 165 199
Fig. 11 Eigen values of a 2nd order LDG growth matrix for RK 3stage scheme. Case (a): non dimensional
time step = 0.157; Case (b): non dimensional time step = 0.138
Fig. 12 Damping occurring in a
2nd order LDG simulation using
two different time steps
mode was never damped. The numerical solution overshoots the analytical solution and one
can expect some oscillations. Therefore steady-state problems may never converge. The time
step needs to be lowered to reduce the magnitude of the eigen value of the spurious mode.
Figure 11b shows the eigen values for this system for a non dimensional time step of 0.138.
In other words, there exists a range, wherein the system converges to a wrong solution. In
general, this phenomenon can be seen in higher orders for LDG, Penalty and BR2 schemes.
The above conclusion was tested using a numerical simulation. Figure 12 shows the
variation of magnitude (L1 norm) of the solution as a function of time for a 2nd order LDG.
The simulation with a non dimensional time step = 0.157 is stable but has no damping. This
is the case wherein the spurious mode dominates over the physical mode. The simulation
with a non dimensional time step = 0.138 is both stable and consistent.
J Sci Comput (2009) 41: 165 199 181
Table 3 Non dimensional time step criteria for obtaining stable and consistent solutions for the LDG, penalty
and BR2 schemes in 1D
Method 2nd order 2nd order 3rd order 3rd order 4th order 4th order
(stability) (no overshoot) (stability) (no overshoot) (stability) (no overshoot)
LDG 0.157 0.138 0.0266 0.0266 0.0109 0.0106
Penalty 0.182 0.175 0.0322 0.0322 0.0133 0.0129
BR2 0.314 0.224 0.0784 0.0784 0.0321 0.0300
Fig. 13 Eigen values of a 2nd order LDG growth matrix for Crank-Nicolson scheme. Case (a): non dimen-
sional time step = 0.5; Case (b): non dimensional time step = 0.25
Table 3 lists the non dimensional time steps required for obtaining a stable solution as
well as a physically relevant solution (i.e. all the spurious modes are damped). Thus in an
explicit simulation, it is generally necessary to run the case with a time step lower than the
cut-off stability limit.
The above pattern was also observed while using the Crank-Nicolson scheme. There
were no stability issues as it is an implicit scheme. However, the spurious modes start to
dominate after a critical value. Figures 13a and 13b show the eigen values of the H matrix
using a non dimensional time step of 0.5 and 0.25 respectively. It is clear that the former has
inadequate damping in most parts of the spectrum, while the latter is suf ciently damped.
6 The p-Multigrid Method
The Gauss-Seidel or Jacobi iterations produce smooth errors when applied on the above
mentioned equations. The error vector has its high frequencies nearly removed in a few it-
erations using a higher order polynomial; but low frequencies are removed very slowly. The
key idea of the p-multigrid method is to solve the nonlinear equations using a lower order
polynomial such that smooth becomes rough and low frequencies act like high frequen-
cies. The p-multigrid method operates on a sequence of solution approximations of different
polynomial orders. Therefore it offers the exibility of switching between higher and lower
polynomial levels without changing the actual geometrical nodal grid points.
182 J Sci Comput (2009) 41: 165 199
p 1 p 2
To accomplish the communication between different levels, restriction (Ip, Ip 1 ) and
p 1
p
prolongation (Ip 1, Ip 2 ) operators are required in addition to the aforementioned relax-
ation scheme as a smoother. Restriction consists of moving solutions and their residuals
at the unknown points from a space of higher polynomial order to a lower polynomial or-
der. Prolongation refers to a reverse procedure in which lower order polynomial solution
correction is redistributed as corrections to the solutions of the unknown points at a higher
polynomial order.
Classical multigrid method begins with a two-level process. First, iterative relaxation is
applied using the higher order polynomial, thus basically reducing high-frequency errors.
Then, a coarse-grid correction is applied, in which the smooth error is determined at the
lower polynomial level. This error is interpolated to the higher polynomial level and used
to correct the existing higher order solutions. Applying this method recursively to solve the
lower polynomial level problems leads to multigrid.
De ning three polynomial levels from the highest to the lowest as p, p 1 and p 2,
we want to solve:
Rp (Qp ) = rp . (6.1)
The rhs rp is zero for the highest level polynomial.
Rp 1 (Qp 1 ) = rp 1, (6.2)
Rp 2 (Qp 2 ) = rp 2 . (6.3)
We employ the implicit LU-SGS schemes as the smoothers for all three levels. The fol-
lowing steps are used to update the solutions on the highest p level in one big V cycle.
Improve Qp by application of a few times the smoother similar as (6.1)
p n
I Rc
Qp = Rp (Q ) rp .
2
(6.4)
p c p
t Qc
Restrict the latest solution Qp to the coarser level for an approximate solution Qp 1
Q0 1 = Ip 1 (Qp ).
p
(6.5)
p
Compute the defect on the nest level
dp = rp Rp (Qp ) = Rp (Qp ). (6.6)
Compute the right hand side of (6.2) as
rp 1 = Rp 1 (Q0 1 ) + Ip 1 dp .
p
(6.7)
p
Improve Qp 1 by application of a few times the smoother
p 2 n
I Rc
Qp 2 = Rp 2 (Q 2 ) rp 2 .
2
(6.8)
p 2 c p
t Qc
p 1
Restrict the latest solution Qp 1 obtained by Qk+1 = Qk 1 + 2
to the coarser
Qc
p 1 p
level for an approximate solution Qs 2
p
Contact this candidate