H.T. Huynh

Glenn Research Center, Cleveland, Ohio

High-Order Space-Time Methods for

Conservation Laws

NASA/TM—2013-218077

September 2013

AIAA–2013–2432
NASA STI Program . . . in Profi le

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientifi c and Technical Information (STI) 

program plays a key part in helping NASA maintain 

this important role.

The NASA STI Program operates under the auspices 

of the Agency Chief Information Offi cer. It collects, 

organizes, provides for archiving, and disseminates 

NASA’s STI. The NASA STI program provides access 

to the NASA Aeronautics and Space Database and 

its public interface, the NASA Technical Reports 

Server, thus providing one of the largest collections 

of aeronautical and space science STI in the world. 

Results are published in both non-NASA channels 

and by NASA in the NASA STI Report Series, which 

includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major signifi cant phase 

of research that present the results of NASA 

programs and include extensive data or theoretical 

analysis. Includes compilations of signifi cant 

scientifi c and technical data and information 

deemed to be of continuing reference value. 

NASA counterpart of peer-reviewed formal 

professional papers but has less stringent 

limitations on manuscript length and extent of 

graphic presentations.

• TECHNICAL MEMORANDUM. Scientifi c 

and technical fi ndings that are preliminary or 

of specialized interest, e.g., quick release 

reports, working papers, and bibliographies that 

contain minimal annotation. Does not contain 

extensive analysis.

• CONTRACTOR REPORT. Scientifi c and 

technical fi ndings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientifi c and technical 

conferences, symposia, seminars, or other 

meetings sponsored or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientifi c, 

technical, or historical information from 

NASA programs, projects, and missions, often 

concerned with subjects having substantial 

public interest.

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientifi c and 

technical material pertinent to NASA’s mission.

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results.

For more information about the NASA STI 

program, see the following:

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question to help@sti.nasa.gov

• Fax your question to the NASA STI 

Information Desk at 443–757–5803

• Phone the NASA STI Information Desk at

 443–757–5802

• Write to:

 STI Information Desk

 NASA Center for AeroSpace Information

 7115 Standard Drive

 Hanover, MD 21076–1320
H.T. Huynh

Glenn Research Center, Cleveland, Ohio

High-Order Space-Time Methods for

Conservation Laws

NASA/TM—2013-218077

September 2013

AIAA–2013–2432

National Aeronautics and

Space Administration

Glenn Research Center

Cleveland, Ohio 44135

Prepared for the

21st Computational Fluid Dynamics Conference

sponsored by the American Institute of Aeronautics and Astronautics

San Diego, California, June 24–27, 2013
Acknowledgments

This work was supported by the Aeronautical Sciences Project of NASA. The author wishes to thank Drs. James DeBonis and 

Seth Spiegel for their thorough reviews and numerous valuable suggestions.

Available from

NASA Center for Aerospace Information

7115 Standard Drive

Hanover, MD 21076–1320

National Technical Information Service

5301 Shawnee Road

Alexandria, VA 22312

Available electronically at http://www.sti.nasa.gov

Level of Review: This material has been technically reviewed by technical management. 
NASA/TM—2013-218077 1 

High-Order Space-Time Methods for Conservation Laws 

H.T. Huynh 

National Aeronautics and Space Administration 

Glenn Research Center 

Cleveland, Ohio 44135

Abstract. Current high-order methods such as discontinuous Galerkin and/or flux reconstruction can 

provide effective discretization for the spatial derivatives. Together with a time discretization, such 

methods result in either too small a time step size in the case of an explicit scheme or a very large system 

in the case of an implicit one. To tackle these problems, two new high-order space-time schemes for 

conservation laws are introduced: the first is explicit and the second, implicit. The explicit method here, 

also called the moment scheme, achieves a CFL (Courant-Friedrichs-Lewy) condition of 1 for the case of 

one-spatial dimension regardless of the degree of the polynomial approximation. (For standard explicit 

methods, if the spatial approximation is of degree ݌, then the time step sizes are typically proportional to 

1/݌ଶ.) Fourier analyses for the one and two-dimensional cases are carried out. The property of super 

accuracy (or super convergence) is discussed. The implicit method is a simplified but optimal version of 

the discontinuous Galerkin scheme applied to time. It reduces to a collocation implicit Runge-Kutta (RK) 

method for ordinary differential equations (ODE) called Radau IIA. The explicit and implicit schemes are 

closely related since they employ the same intermediate time levels, and the former can serve as a key 

building block in an iterative procedure for the latter. A limiting technique for the piecewise linear 

scheme is also discussed. The technique can suppress oscillations near a discontinuity while preserving 

accuracy near extrema. Preliminary numerical results are shown. 

1. Introduction 

The discontinuous Galerkin (DG) method is currently among the most widely used high-order 

numerical methods for solving the compressible Navier-Stokes equations on unstructured meshes. It was 

introduced for the neutron transport equation by Reed and Hill (1973), analyzed by LaSaint and Raviart 

(1974) and developed and made popular for fluid dynamics equations by Cockburn, Shu, Bassi, Rebay, 

and others (see e.g., Bassi and Rebay 1997a,b, Cockburn, Karniadakis, and Shu 2000, Cockburn and Shu 

2005, 2009, Shu 2012, and the references therein). Efficient DG schemes using the Lagrange polynomials 

for nodal points as basis functions known as nodal DG methods can be found in (Hesthaven and 

Warburton 2008). Alternative methods employing the differential form (as opposed to DG, which uses the 

integral form) have been proposed. Kopriva and Kolias (1996) pioneered this approach with the 

staggered-grid method on quadrilateral meshes. The method was extended to triangular meshes by Liu, 

Vinokur, and Wang (2006) and named spectral difference (SD). Another class of schemes called spectral 

volume (SV) was presented by Wang, Zhang, and Liu (2004). 

Recently, an approach to high-order accuracy with the advantage of simplicity and economy called 

flux reconstruction (FR) was introduced in (Huynh 2007, 2009a). The approach employs the differential 

form and results in several new schemes with favorable properties. In addition, by recovering DG, SD, 

and SV, it unifies existing methods. Extensions of the FR method to 2D unstructured meshes were 

presented in (Wang and Gao 2009) and (Gao and Wang 2009) and was named lifting collocation penalty 

or LCP schemes. The involved authors later combined FR and LCP and called them correction procedure 

via reconstruction or CPR methods. Extension to 3D was presented in (Haga, Gao, and Wang 2010, 2011) 

 
NASA/TM—2013-218077 2 

 and (Wang, Gao, and Haga 2011). A mathematical foundation for FR was provided via energy 

stability proofs: for an SD scheme identified via FR by Jameson (2010), for a one-parameter class of FR 

schemes in 1D by Vincent, Castonguay, and Jameson (2011), for 2D triangular grids by Castonguay, 

Vincent and Jameson (2012), for linear-advection-diffusion problem in 1D by Castonguay et al. (2013) 

and for advection-diffusion in 2D, by Williams et al. (2011). Extension for a particularly simple FR 

scheme called ݃ଶ to the Navier-Stokes equations was presented in (Liang, Miyaji, and Zhang 2013). 

Applications of FR/CPR methods to turbulent internal flows for turbomachinaries were elaborated in (Lu, 

Yuan, and Dawes 2012 and Lu, Liu, and Dawes 2013). Interface elements dealing with non-conforming 

polynomials together with p-adaptation for viscous flow simulations were presented in (Cagnone and 

Nadarajah 2012, Cagnone, Vermeire, and Nadarajah 2013), and an implicit large eddy simulation (ILES) 

solver was developed for FR/CPR schemes by Vermeire, Cagnone, and Nadarajah (2013). 

The high-order methods discussed can provide effective discretization for the spatial derivatives. 

Together with a time discretization, however, such methods result in either too small a time step size in 

the case of an explicit scheme or a very large system in the case of an implicit one. Efforts to improve 

time stepping have been proposed by several authors. 

Concerning explicit schemes, a space-time DG method with a favorable CFL (Courant-FriedrichsLewy) condition employing a staggered mesh was introduced in Lowrie et. al (1995). For the case of a 

semi-discrete formulation, Warburton and Hagstrom (2008) presented a staggered-mesh scheme where 

the solution is projected onto a mesh staggered from the original one and then projected back onto the 

original mesh resulting in a solution ‘smoother’ than the original. They obtain a CFL condition 

proportional to 1/݌ as opposed to 1/݌ଶ where ݌ is the degree of the approximating polynomials. 

However, multi-dimensional extensions of these staggered-mesh approaches are both complex and costly. 

Additional efforts can be found under the heading “Taming the CFL number” of (Hesthaven and 

Warburton 2008) and the references therein. In addition to these efforts, Dumbser and Munz (2005) 

developed ADER-DG using the ADER finite-volume scheme (Arbitrary order using derivatives). Spacetime expansion DG (or STE-DG) as a modification of ADER-DG was introduced in (Lörcher, Gassner, 

and Munz 2007 and Gassner, Lörcher, and Munz 2008), but these methods have highly restricted CFL 

limits. 

Researchers in the finite-volume community have devised several explicit methods with a CFL 

condition of 1 for the case of one spatial dimension. In fact, in 1977, Van Leer presented a series of five 

schemes for convection with this property. Among the five, scheme I is the least accurate but becomes the 

most popular and, in the literature, is generally the scheme implied when the MUSCL approach is 

referred. Concerning the more accurate schemes such as III and V, the problem is, as stated in (Van Leer 

and Nomura 2005), “When trying to extend these schemes beyond advection, viz., to a nonlinear 

hyperbolic system like the Euler equations, the first author ran into insuperable difficulties because the 

exact shift operator no longer applies, and he abandoned the idea”. 

For scheme III, which can be considered as a piecewise-linear space-time DG method for convection, 

this difficulty was overcome in (Huynh 2006) where the exact shift operator was replaced by (a) an 

integration by parts (Gauss theorem), (b) the observation that the space-time fluxes contain the 

information provided by the shift operator, and (c) a successive procedure of updating the moments. The 

resulting method is called the upwind moment scheme. The approach was further analyzed and applied to 

nonlinear hyperbolic equations by Suzuki and Van Leer (2007). As briefly discussed in (Huynh 2007), the 

moment scheme can be extended to arbitrary order. Such an extension was independently obtained by Lo 

and was studied in combination with Van Leer’s recovery scheme for diffusion (2011). The description of 

the extension here and that of Lo thus have considerable similarities (the differences will also be 

discussed). Extensions for scheme V can be found in (Eymann and Roe 2013). 

Concerning implicit schemes, it was shown by LaSaint and Raviart (1974) via a quadrature of left

Radau type (Hildebrand 1987) that discretizing an ODE via discontinuous Galerkin results in an implicit 
NASA/TM—2013-218077 3 

Runge-Kutta (RK) method. In addition, the scheme is accurate to order 2݌ ൅ 1 if the solution is 

approximated by a polynomial of degree ݌ in time. Implicit space-time DG schemes for conservation 

laws have been studied by several authors (e.g., Bar-Yoseph and Elata 1990, Van der Vegt and Van der 

Ven 2002, Pesch and Van der Vegt 2008). The common difficulty is that the implicit system is large. 

Recent effort includes the large eddy simulation solver by Vermeire, Cagnone, and Nadarajah (2013). 

Another aspect of time stepping methods concerns their performance for large scale calculations using 

massively parallel supercomputers or graphics processing units (GPU). Here, the DG or FR/CPR spatial 

discretization is highly parallelizable. It is well known that this ease of parallelization is retained with an 

explicit time stepping, e.g., RK, but not an implicit one. 

In this paper, to tackle the time stepping problem, two new high-order space-time schemes for 

conservation laws are introduced: the first is explicit and the second, implicit. The explicit method, called 

the moment scheme, achieves a CFL condition of 1 for the case of one-spatial dimension regardless of the 

degree of approximation. Compared with standard approaches, the key difference is the built-in physics of 

advection. Such a CFL condition is a significant gain compared with a typical time step size proportional 

to 1/݌ଶ. For the case of multiple-dimensions, however, there is a reduction in the CFL limits. The current 

extension differs from that of Lo in (a) the case of two-spatial dimensions and its Fourier analysis are 

presented and (b) additional details, observations, and the relation with the implicit time stepping are 

discussed. It will be shown via Fourier analysis that for both one and two-dimensional cases, the moment 

schemes are accurate to order 2݌ ൅ 1 if the solution is approximated by polynomials of degree ݌. 

The implicit method here is a simplified but optimal DG scheme applied to time. It employs the FR 

approach and simplification is made possible by the fact that time always moves forward (as opposed to 

space where waves can travel in all directions). The method reduces to an implicit RK scheme using the 

right Radau points as collocation points called the Radau IIA (Lambert 1991, Hairer, Norsett, and Wanner 

1991, Hairer and Wanner 1993). It is also accurate to order 2݌ ൅ 1 if polynomial of degree ݌ is employed 

to approximate the solution. A key advantage of using the right Radau points (as opposed to the left 

Radau points employed by LaSaint and Raviart 1974) is that the method is simpler and involves fewer 

interpolations in the solution procedure. More importantly, the explicit and implicit schemes here are 

closely related since they use the same intermediate time levels, and the former can serve as a key 

building block in an iterative procedure for the latter. 

The DG derivation here also differs from those in the literature: the standard DG formulation is highly 

algebraic, and the behavior of the solution is hidden in a system of equations; the current formulation is 

geometric, and the behavior of the solution can be observed. 

In addition to the time stepping methods, since the numerical tests have discontinuities, a limiting 

technique that can suppress oscillations near discontinuities while preserving accuracy near extrema is 

discussed. Finally, preliminary numerical results are shown. 

This paper is essentially self-contained and organized as follows. In Section 2, we review Van Leer’s 

scheme III except that our derivation is for arbitrary order. The moment schemes, which have a CFL 

condition of 1 and can be considered as explicit space-time DG methods, are presented in Section 3. 

Extension to the case of two spatial dimensions is discussed in Section 4. Fourier analyses are carried out 

in Section 5 for the one- as well as two-dimensional cases. DG time discretization in the FR/CPR 

framework is presented in Section 6. Section 7 deals with limiting. Preliminary numerical solutions are 

shown in Section 8. Finally, conclusions and discussions can be found in Section 9. 

2. Exact Time Evolution and Projection 

In the fourth installment of a celebrated series of five papers entitled “Towards the ultimate 

conservative difference scheme”, Van Leer (1977) introduced five schemes for advection; among them, 

the first three are piecewise linear and, the last two, piecewise parabolic. Concerning the piecewise linear 
NASA/TM—2013-218077 4 

methods, scheme I is the least accurate but, due to its ease of extension to systems, becomes the most 

popular and is typically known as the MUSCL scheme (monotone upstream-centered schemes for 

conservation laws). As shown in (Huynh 2003), scheme II formulated for a staggered mesh results in the 

CE/SE scheme with ߳ ൌ 1/2 (conservation element/solution element, Chang 1995). Scheme III is 

remarkable since it is third-order accurate in spite of being piecewise linear. The property of having 

higher than expected order of accuracy is called super accuracy or super convergence. The last two 

schemes are piecewise parabolic. Here, scheme IV leads to the standard PPM (piecewise parabolic 

method). Scheme V stores and updates both the cell averages and interface values. Schemes III and V 

have the same error, which is only 1/9 that of scheme IV. As mentioned above, when trying to extend the 

more accurate schemes to the Euler equations, Van Leer “ran into insuperable difficulties”. However, 

scheme III plays a critical role in serving as a guide for the moment scheme (Huynh 2006), which can be 

extended to systems with relative ease. Below, we review Van Leer’s approach for scheme III, which 

amounts to a projection after an exact shift to account for time evolution. Our presentation of the scheme 

is for arbitrary order, not just piecewise linear. 

Consider the scalar advection equation 

ݑ௧ ൅ ܽݑ௫ ൌ 0 (2.1)

with initial condition 

ݑሺݔ, 0ሻ ൌ ݑ଴ሺݔሻ (2.2)

where ݐ is time, ݔ space, and ܽ the advection speed, a positive constant. By assuming that ݑ ଴is periodic 

or of compact support, boundary conditions are trivial and are omitted. The exact solution at time ݐ is 

obtained by shifting the data curve to the right a distance ܽݐ, 

ݑୣ୶ୟୡ୲ሺݔ, ݐሻ ൌ ݑ଴ሺݔ െ ܽݐሻ. (2.3)

Legendre Polynomials. To approximate the solution, we need some preparations. Let ݇ be a 

nonnegative integer. Denote by ࡼ ௞the space of polynomials of degree ൑ ݇. As is standard in the finitevolume community, for the reference interval, we use ܫ ൌ ሾെ1/2, 1/2ሿ. Let the inner product of any two 

functions ݒଵ and ݒଶ on ܫ be defined by 

ሺݒଵ, ݒଶሻ ൌ න ݒଵሺߦሻ ݒଶሺߦሻ

ଵ/ଶ

ିଵ/ଶ

(2.4) ,ߦ݀

and the ܮଶ norm of a function ݒ by

‖ݒ ‖ൌ ቆන ൫ݒሺߦሻ൯

ଶ ଵ/ଶ

ିଵ/ଶ ݀ߦቇଵ/ଶ

. (2.5)

Let ܮ ௞be the Legendre polynomial of degree ݇ defined by applying the Gram-Schmidt process to 

orthogonalize the monomials 1, ߦ, ߦଶ, ߦଷ, ... together with the condition that ܮ௞ሺ1/2ሻ ൌ 1. The first few 

Legendre polynomials on ܫ ൌ ሾെ1/2, 1/2ሿ are 

ܮ଴ሺߦሻ ൌ 1, 				ܮଵሺߦሻ ൌ 2ߦ, ܮଶሺߦሻ ൌ 6ߦଶ െ 1/2,

ܮଷሺߦሻ ൌ 20ߦଷ െ 3ߦ, ܮସሺߦሻ ൌ 70ߦସ െ 15ߦଶ ൅ 3/8. (2.6)
NASA/TM—2013-218077 5 

Their plots are shown in Fig. 2.1. Note that compared with the standard Legendre polynomials ܲ௞ defined 

on ሾെ1, 1ሿ, for ߦ in ሾെ1/2, 1/2ሿ, ܮ௞ሺߦሻ ൌ ܲ௞ሺ2ߦሻ. 

Fig. 2.1 Legendre Polynomials

Discretization. Let the domain of calculation be divided into non-overlapping cells or elements 

denoted by ܧ ௝ൌ ൣݔ௝ିଵ/ଶ, ݔ௝ାଵ/ଶ൧ of centers ݔ ௝and widths ݄௝. For each ܧ௝, let the local coordinate ߦ on 

ܫ ൌ ሾെ1/2, 1/2ሿ be 

ߦ ൌ ሺݔ െ ݔ௝ሻ/݄௝. (2.7)

Next, let ݌൒0 be an integer. On each cell ܧ௝, we wish to approximate the solution by a polynomial of 

degree ݌ denoted by ݓ௝. To this end, we need to define the projection operator. Consider the ݌൅1

Legendre polynomials ܮ௞, ݇ ൌ 0, ... , ݌. With ݐ fixed, the solution ݑሺݔ, ݐሻ can be approximated on ܧ݆ by its 

projection onto ࡼ௣: 

ሻߦሺ௞ܮ ௞ ,௝ݑ෍ ൌ ሻߦሺ௝ݓ

௣

௞ୀ଴

(2.8)

where the coefficients are given by 

ሻߦሺ௞ܮ ൯ݐ ,௝݄ߦ ൅ ௝ݔ൫ݑ න 		 ൌ ௞ ,௝ݑ

ଵ/ଶ

ିଵ/ଶ

(2.9) . ൗ ଶ௞‖ܮ‖ ߦ݀

Thus, on ܧ௝, ݓ ௝is the least square fit of degree ݌ for ݑ, and since the projection does not increase the 

norm, we have ฮݓ௝ฮ ൑ ‖ݑ‖. The collection of all ݓ ௝as ݆ varies forms a piecewise polynomial function 

denoted by, in the physical coordinate, ݓሺݔ, ݐሻ. The function ݓሺݔ, ݐሻ is generally discontinuous at the 

cell interfaces (see Fig. 2.2(a)). 

Note that ݑ௝, ଴ is the cell average value on ܧ௝, and ݑ௝, ଵ represents half of the normalized slope so that 

for ݌ൌ1, the value at the right interface is ݑ௝, ൅ ଴ ݑ௝, ଵ. In fact, for all ݌, the value at the right interface is 

ݑ௝, ൅ ଴ ݑ௝, ଵ ൅⋯൅ݑ௝,௣ , i.e., each ݑ௝, ௞ results in a ݇-th degree approximation to the right interface value. 
NASA/TM—2013-218077 6 

For simplicity but not necessity, assume that the mesh is uniform: ݄௝ ൌ ݄. Denote the time step by Δݐ

and the CFL number ܽΔݐ/݄ by ߪ. Then, since ܽ൒0, the CFL condition takes the form 

(2.10) 1. ൑ ߪ ൑ 0

That is, in one time step, the wave advects a distance no more than one cell width. It is this property that 

we wish to preserve when extending the method to one-dimensional systems of equations. 

Let the projection of the initial data be carried out and denoted by, on ܧ݆,

௝ݓ

଴ሺߦ, 0ሻ ൌ ෍ݑ௝,௞ 

ሻߦሺ௞ܮ ଴

௣

௞ୀ଴

.

Exact Time Evolution and Projection. Assume that at time ݐ௡, the projection ݓ௝ሺߦ, ݐ௡ሻ is known; 

i.e., ݑ݆, ݇

௡ ൌ ݑ௝, ௞ are known for all ݆ and all ݇൑݌) the superscript ݊ is understood). We wish to calculate 

݇ ݆,ݑ all

௡ାଵ at time ݐ௡ାଵ (the superscript ݊൅1 is retained). 

Given the piecewise polynomial data ݓሺݔ, ݐ௡ሻ ൌ ݓሺݔሻ ൌ ሼݓ௝ሺݔሻሽ, the exact solution at time ݐ௡ାଵ for 

this data is obtained by shifting the curve 	ݓሺݔሻ to the right a distance ܽΔݐ as shown in Fig. 2.2(b): 

ݒሺݔሻ ൌ ݓሺݔ െ ܽΔݐሻ.

That is, on ܧ௝, with the local coordinate ߦ given by (2.7), 

ݒ௝ሺߦሻ ൌ	ቊݓ௝ିଵሺߦെߪ൅1ሻ, if ߦ ൏ െ1/2 ൅ ߪ

(2.11) . ߪ ൅ 1/2െ ൒ ߦ if ,ሻߪെߦሺ௝ݓ

To evaluate the solution, we split the integral into two parts (again see Fig. 2.2(b)); with ݇ ൌ 0, ... , ݌, 

݇ ݆,ݑ

௡ାଵ ൌ 	ቈන ݓ݆െ1ሺߦെߪ൅1ሻ ܮ݇ሺߦሻ

ߪ൅1/2െ

െ1/2

ሻߦሺ݇ܮ ሻߪെߦሺ݆ݓ න ൅ ߦ݀

1/2

ߪ൅1/2െ

(2.12) . ൗ 2‖݇ܮ‖ ቉ߦ݀

Such a solution can easily be evaluated by a mathematical programming package such as Mathematica. 

For the case ݌ൌ0 or piecewise constant data, scheme III reduces to the first-order upwind scheme. 

For the case ݌ൌ1 or piecewise linear, the solution is 

0 ݆,ݑ

௡ାଵ ൌ	ݑ௝, ൅ ௞ ߪሺݑ௝ିଵ, ଴ ൅ ݑ௝ିଵ, ଵ െ ݑ௝, ଴ െ ݑ௝, ଵሻ൅ߪ2ሺെݑ௝ିଵ, ଵ ൅ ݑ௝, ଵሻ, (2.13)

1 ݆,ݑ

௡ାଵ ൌ	ݑ௝, ଵ ൅ ߪ൫െ3ݑ௝ିଵ, ଴ െ 3ݑ௝ିଵ, ଵ െ 3ݑ௝, ଴ െ 3ݑ௝, ଵ൯ ൅

																															ߪଶ൫3ݑ݆െ1, 0 ൅ 6ݑ݆െ1, 1 െ 3ݑ݆, 0൯൅ߪଷ൫െ2ݑ݆െ1, 1 ൅ 2ݑ݆, 1൯

(2.14)

In matrix form, 
NASA/TM—2013-218077 7 

ቆ

0 ݆,ݑ

௡ାଵ

1 ݆,ݑ

௡ାଵቇ ൌ	 ൬ ߪ ߪሺ1െߪሻ

െ3ߪሺ1െߪሻ െߪሺ3 െ 6ߪ ൅ 2ߪ2ሻ

൰ ቀ

଴ ,ଵି௝ݑ

ݑ௝ିଵ, ଵ ቁ	൅		

														൬ 1 െ ߪ െߪሺ1െߪሻ

	3ߪሺ1െߪሻ ሺ1െߪሻሺ1 െ 2ߪ െ 2ߪ2ሻ

൰ ቀ

	଴ ,௝ݑ

ݑ௝, ଵ ቁ

(2.15)

(a) Data (b) Solution 

Fig. 2.2 Scheme III (a) Piecewise linear data; (b) Solution obtained by shifting the data a distance ܽݐ߂

and calculating the moments of the discontinuous function. 

Note the term ߪଷ above for the case ݌ൌ1. For a general ݌, the function ݓ௝ሺߦെߪሻ	ܮ௣ሺߦሻ is of degree 

2݌ in ߦ ;as a result, the term of highest degree in the solution is ߪଶ௣ାଵ. It will be shown by Fourier 

analysis that the method is accurate to order 2݌ ൅ 1 and is stable for 0൑ߪ൑1. 

Also note that the method is clearly energy-stable: the ܮଶ-norm of the solution at time ݐ௡ାଵ is bounded 

by that of the data at time ݐ ௡since the solution is obtained by an orthogonal projection. 

As the final remark of this section, we discuss the extension of this approach. For the advection case, 

the discontinuity at an interface evolves via the exact shift operator. In the case of systems such as the 

Euler equations, in one spatial dimension, these discontinuities result in some combination of a shock, a 

contact, or a fan. Tracking these waves accurately is extremely difficult. Resolving these waves for the 

multi-dimensional cases appears to be an impossible task. 

3. Moment Scheme for the Case of One Spatial Dimension 

Formulation. To avoid the exact shift operator—which makes extensions difficult or impossible—we 

take a different path by using a space-time domain and integration by parts together with the observation 

that the space-time fluxes contain the information provided by the shift operator. Consider the scalar 

conservation law 

ݑ௧ ൅ ݂௫ ൌ 0 (3.1)

0.5 1.0 1.5 2.0 2.5 3.0 3.5

x

0.5

1.0

1.5

2.0

2.5

3.0

u

0.5 1.0 1.5 2.0 2.5 3.0 3.5

x

0.5

1.0

1.5

2.0

2.5

3.0

Cell j u

Cell ܽΔݐ jെ1 
NASA/TM—2013-218077 8 

where ݐ is time, ݔ spatial coordinate, and ݂ ൌ ݂ሺݑሺݔ, ݐሻሻ the flux. With the wave speed ܽሺݑሻ ൌ ݂݀/݀ݑ, 

(3.1) implies 

ݑ௧ ൅ ܽݑ௫ ൌ 0. (3.2)

On the cell ܧൌܧ௝, let ߶ be a test function, i.e., a polynomial of degree ݌. The conservation law in 

weak form is 

൅ ݔ݀ ሻݔሺ ߶ሻݐ ,ݔ௧ሺݑ	න

ா

න ݂௫ ߶ ݀ݔ

ா

ൌ 0. (3.3)

We will set ߶ൌܮ௞, the Legrendre basis function, later. Thus, the weak form amounts to projecting onto 

ࡼ௣. As in the DG formulation, by applying integration by parts to the second term, 

∂

൅ ݔ݀ ሻݔሺ	߶ሻݐ ,ݔሺݑන ݐ∂

ா

ሾ݂ ߶ሿ௫ೕషభ/మ

௫ೕశభ/మ െ න ݂ ߶௫ ݀ݔ

ா

ൌ 0.		 (3.4)

The fluxes ݂ for the boundary term ሾ݂	߶ሿ௫ೕషభ/మ

௫ೕశభ/మ are chosen to be the upwind ones. At each interface, this 

flux is common for the two adjacent cells and allows the data in these cells to interact. Without this 

common flux, there is no interaction, and the resulting method is inconsistent with physics. 

Integrating in time from ݐ ௡to ݐ) ∗ݐ ∗will be set to the intermediate time levels as well as the final time 

ݐ௡ାଵ later), we obtain 

	ݔ݀	ሻݔሺ	߶ሻ∗ݐ ,ݔሺݑන

ா

	െ නݑሺݔ, ݐ௡ሻ ߶ሺݔሻ ݀ݔ

ா

൅ න ሾ݂୳୮୵ ߶ሿ௫ೕషభ/మ

௫ೕశభ/మ

௧∗

௧೙

		ݐ݀	

െ 		න න ݂ ߶௫ ݀ݔ ݀ݐ

ா

௧∗

௧೙

ൌ 0.

(3.5)

Denote the surface and volume integral by, respectively, 

SI ൌ െ න ሾ݂୳୮୵ ߶ሿ௫ೕషభ/మ

௫ೕశభ/మ

௧∗

௧೙

(3.6) ݐ݀

and 

VI ൌ න න ݂ ߶௫ ݀ݔ ݀ݐ

ா

௧∗

௧೙

. (3.7)

The problem reduces to the evaluation of SI and VI, 

ݔ݀ ሻݔሺ	߶ሻ∗ݐ ,ݔሺݑන

ா

ൌ නݑሺݔ, ݐ௡ሻ ߶ሺݔሻ ݀ݔ

ா

൅ SI ൅ VI. (3.8)

Here, for ߶ൌܮ ௞and ݐ ∗ൌ ݐ௡ାଵ, the left hand side above in the local coordinate yields the solution 

௞ ,௝ݑ	ଶ௞‖ܮ‖௝݄

݊൅1 whereas the first term on right hand side, ݄௝‖ܮ‖௞ଶ	ݑ௝,௞ 

 .௞ ,௝ݑ	ଶ௞‖ܮ‖௝݄ ൌ ݊

Consider now the case of advection. So far, the current approach is equivalent to the exact evolution 

approach of the previous section. Indeed, with the piecewise polynomial data ݓ௝ሺݔሻ at time ݐ௡, the exact 
NASA/TM—2013-218077 9 

space-time solution is ݓሺݔ, ݐሻ ൌ ݓ௝ሺݔ െ ܽݐሻ, which has discontinuities along the characteristic lines 

joining the points ሺݔ௝ିଵ/ଶ, ݐ௡ሻ and ሺݔ௝ିଵ/ଶ ൅ ܽΔݐ, ݐ௡ାଵሻ as shown in Fig. 3.1. The two formulations, one 

by projecting ݓሺݔ, ݐ௡ାଵሻ, and the other by integrating on the space-time cell ܧ ௝ൈ ሾݐ௡, ݐ௡ାଵሿ yield 

identical solutions, which are exact for the data. At first glance, it appears that we have made the problem 

more difficult since the volume integral takes the form 

VI ൌ න න ܽݑ߶ ௫ ݀ݔ ݀ݐ

ா

௧∗

௧೙

. (3.9)

In the previous section, we only deal with a discontinuity at one point in space; here, we have to deal with 

a discontinuity along a line in a space-time domain. However, as will be shown, the term ߶௫ comes to the 

rescue in a procedure which ignores the discontinuity. 

Fig. 3.1. Space-time cells ܧ ௝ൈ ሾݐ௡, ݐ௡ାଵሿ and the paths of discontinuities (dark line segments) for the 

case of advection 

Time Integration. We proceed by discussing time integration. Due to the term ߪଶ௣ାଵ in the solution 

as noted in the previous section, we need a quadrature with a degree of precision of 2݌ in time, i.e., it is 

exact for polynomials of degree 2݌ or less (note that we only need a degree of precision of 2݌ since the 

vertical edge of ܧ ௝ൈ ሾݐ௡, ݐ௡ାଵሿ accounts for the additional power in ߪଶ௣ାଵ). Here, we use to the right 

Radau quadrature with ݌൅1 evaluation points—consistent with the spatial degree of freedom. Another 

reason for using this quadrature is that the DG method applied to time can be reduced to the Radau IIA 

method (Huynh 2009b), which is a collocation method corresponding to this quadrature. In other words, 

the moment scheme formulated here has the same intermediate time levels and is closely related to the 

implicit Runge-Kutta (RK) scheme resulting from a space-time DG discretization (see §6). 

A brief review of the right Radau quadrature and the associated Butcher matrix is in order (we can also 

employ the left Radau quadrature as in LaSaint and Raviart (1974); the resulting method is, however, 

more complex). Set ݏൌ݌൅1 where ݏ is the number of stages. Using notations from ODE, on the 

interval ሾ0, 1ሿ, let the ݏ right Radau points be denoted by 	ܿ௟, ݈ ൌ 1, ... , ݏ where ܿ௦ ൌ 1. Recall that these 

points result in a quadrature formula on ሾ0,1ሿ that is exact for polynomials of degree 2ሺݏ െ 1ሻ or less. For 

each ݈, let the corresponding Lagrange polynomials on ሾ0, 1ሿ be denoted by La௟ and defined by the 

condition that La௟ሺܿ௠ሻ ൌ 0 for all ്݈݉ and La௟ሺܿ௟ሻ ൌ 1:
NASA/TM—2013-218077 10 

La௟ሺ߬ሻ ൌ ෑ ߬െܿ௠

ܿ௟ െ ܿ௠

௦

௠ୀଵ

௠ஷ௟

. (3.10)

The plots of these functions for the case ݏൌ2 and ݏൌ3 ( where ݏൌ݌൅1 (are shown in Fig. 3.2. 

(a) ݏൌ2 (correspoding	to		݌ ൌ 1ሻ (b) ݏൌ3

Fig. 3.2 Lagrange polynomials for the right Radau points: (a) Two points 	ܿଵ ൌ 1/3, ܿଶ ൌ 1, (b) Three 

points ܿଵ ൌ 	0.155051, ܿଶ ൌ 0.644949	, and ܿଷ ൌ 1. 

Next, suppose the values of a function ݒ are known at these Radau points: ݒ ௠ൌ ݒሺܿ௠ሻ, ݉ ൌ 1, ... , ݏ. 

We wish to obtain a quadrature formula for ׬ ݒሺ߬ሻ݀߬ ௖೗

଴ , ݈ ൌ 1, ... , ݏ. To this end, the polynomial of degree 

݌ൌݏെ1 interpolating ݒ௠, ݉ ൌ 1, ... , ݏ, can be expressed as 

ݒሺ߬ሻ ൌ෍ ݒ ௠La௠ሺ߬ሻ

௦

௠ୀଵ

. (3.11)

Integrating from 0 to ܿ௟, we obtain 

න ݒሺ߬ሻ݀߬

௖೗

଴

௠ݒ ෍ൌ

௦

௠ୀଵ

න La௠ሺ߬ሻ ݀߬

௖೗

଴

. (3.12)

Denote the Butcher matrix of dimension ݏൈݏ by ܣ. The entries ܽ௟௠ of ܣ are given by 

ܽ௟௠ ൌ න La௠ሺ߬ሻ ݀߬

௖೗

଴

. (3.13)

These ܽ௟௠ can easily be obtained by using a software package such as Mathematica. The above two 

equations imply, for each ݈, 

න ݒሺ߬ሻ݀߬

௖೗

଴

௠ݒ ௠௟ܽ ෍ൌ

௦

௠ୀଵ

. (3.14)

Note that for each row ݈, setting ݒൌ1, we obtain ܿ௟ ൌ ∑ ܽ௟௠ ௦

௠ୀଵ . Also note that for ݈ ൌ 1, ... , ݌, the 

above quadrature has degree of precision of only ݌. However, for ݈ൌݏൌ݌൅1, i.e., for the integral from 
NASA/TM—2013-218077 11 

0 to 1, due to the definition of the Radau points, the corresponding quadrature has a degree of precision of 

2݌. The column vector ܥ ൌ ሼܿ௟ሽ and the matrix ܣ for the case ݌ൌ1 and ݌ൌ2 are, respectively 

ቂ

1/3

1 ቃ ൤

5/12 െ1/12

3/4 1/4 ൨			and					 ൥0.155051

0.644949

1

൩ ൥0.196815 െ0.0655354 		0.023771

0.394424 0.292073 െ0.0415488

0.376403 0.512486 		0.111111

൩. (3.15)

Returning to the conservation law, for the space-time domain ܧ ௝ൈ ሾݐ௡, ݐ௡ାଵሿ, denote the time level for 

each stage ݈ of the ݏ stages of the above right Radau quadrature by 

(3.16) .ݐΔ௟ ܿ௡ ൅ ݐ ൌ ௟ ,௡ݐ

These are the stages of a Runge-Kutta type scheme where ݐ௡, ௦ ൌ ݐ௡ାଵ. See Fig. 3.3(a) and (b). 

Space-Time Solution with No Interaction. Next, we discuss an approximation for the solution on the 

space-time cell ܧ ௝ൈ ሾݐ௡, ݐ௡ାଵሿ with no interaction (NI) i.e., no upwinding. This approximation denoted 

௝ݓ by

ேூሺݔ, ݐሻ is a polynomial of degree ݌, which can be a Taylor series around ൫ݔ௝,	 ݐ௡൯. The case of 

݌ൌ1 is trivial: ݓ௝ሺݔሻ provides ݑ௫ and (3.2), ݑ௧. For a general ݌, in each cell ܧ௝, at time ݐ௡, the solution 

is known and is a polynomial of degree ݌ in ݔ, 

ሻݔሺ௞ܮ ௞ ,௝ݑ෍ ൌ ሻݔሺ௝ݓ

௣

௞ୀ଴

. (3.17)

From (3.2), we can calculate all the time as well as mixed space-time derivatives up to degree ݌, i.e., ݑ௧, 

ݑ௫௧, 	ݑ௧௧, ... at ൫ݔ௝,	 ݐ௡൯ via the Cauchy-Kovalevsky procedure (Harten, Osher, Engquist, and 

Chakravarthy 1987). The space-time Taylor series expansion around ൫ݔ௝,	 ݐ௡൯ up to degree ݌ then yields 

௝ݓ

ேூሺݔ, ݐሻ on ܧ ௝ൈ ሾݐ௡,	 ݐ௡ାଵሿ. 

As an alternative to the above Taylor series expansion, we can use the nodal Lobatto points together 

with a RK time stepping with no interaction. For example, consider the case ݌ൌ3. See Fig. 3.3(b). On 

ܧ௝, at time ݐ௡, since the cubic ݓ௝ሺݔሻ is known, we can obtain the values at the four Lobatto points via a 

ሺ4 ൈ 4ሻ matrix multiplication. Next, applying the standard 4-stage RK scheme to (3.2) with no 

interaction, we can evaluate the nodal values at time ݐ௡ାଵ. From the known nodal values for ݑ at ݐ ௡and 

ݐ௡ାଵ, we can calculate the corresponding ݑ௧ again with no interaction by (3.2). Then, using a Hermite 

cubic interpolation with the two ends at ݐ ௡and ݐ௡ାଵ, we can evaluate the nodal values accurate to degree 

݌ൌ3 for all intermediate time levels ݐ௡, ௟ ൌ ݐ ൅ ௡ܿ	௟Δݐ, ݈ ൌ 1, 2, 3. Note that this step is inexpensive and 

the equations for either primitive or conservative variables can be employed with the former more 

economical than the latter. 

An observation concerning the case of advection with ܽ൒0 is in order. We focus, for the moment, on 

only one cell, say, ܧ ௝with the polynomial data ݓ௝. For the space-time domain ܧ ௝ൈ ሾݐ௡, ݐ௡ାଵሿ, along the 

two interfaces ሺݔ௝ିଵ/ଶ, ݐሻ and ሺݔ௝ାଵ/ଶ, ݐሻ, the fluxes are identical to those by method of characteristics 

applied to ݓ ௝assuming the domain for ݓ௝ሺݔሻ is ሺെ∞, ∞ሻ. Indeed, both the method of characteristics and 

space-time Taylor series produce the exact solution for the polynomial data. However, when interaction 

comes into play, whereas along the right interface, the upwind fluxes reduces to the result by method of 

characteristics for cell ݆, along the left interface, the values by ݓ ௝at ݔ௝ିଵ/ଶ െ ܽ	ܿ௟	Δݐ, namely, 

௝ݓܽ

ேூ൫ݔ௝ିଵ/ଶ, ݐ௡,௟ ൯ൌݓ௝൫ݔ௝ିଵ/ଶ െ ܽ	ܿ௟	Δݐ൯, are from outside the cell ܧ ௝and will be discarded by the 
NASA/TM—2013-218077 12 

upwind step. In other words, along ሺݔ௝ିଵ/ଶ, ݐሻ, the upwind procedure sorts out and employs the values by 

 .ଵି௝ݓ

Thus, for general conservation laws, via one of the above two procedures (Cauchy-Kovalevsky or 

explicit RK), for all intermediate time levels ݐ௡, ௟ of the right Radau type, at each cell interface ݔ௝ାଵ/ଶ, we 

have two values ݑ ௅and ݑோ (no interaction). See Fig. 3.3. These values in turn yield the upwind flux 

݂୳୮୵ሺݔ௝ାଵ/ଶ, ݐ௡,௟ ሻ. 

(a) ݌ൌ1 (b) ݌ൌ3

Fig. 3.3 Space-time nodal points with space via Lobatto and time via right Radau: (a) ݌ൌ1, and (b) 

݌ൌ3. The number of spatial points, ݌൅1, is the same as the number of Legendre basis functions. Along 

each interface, at each Radau intermediate time level, we can obtain the values ݑ ௅and ݑோ, which in turn 

yield the upwind flux. 

Moment Solution. We can now discuss the moment solution. On the space-time domain ܧ ௝ൈ

ሾݐ௡, ݐ௡ାଵሿ, with ݐ ∗ൌ ݐ௡,௟ , switching to the local coordinate ߦ on ܫ ൌ ሾെ1/2, 1/2ሿ, since ݀ݔ ൌ ݄௝݀ߦ,) 3.8 (

yields 

௞ ,௝ݑ	ଶ௞‖ܮ‖௝݄	

௞ ,௝ݑ	ଶ௞‖ܮ‖௝݄	ൌ ݈ ݊,

																				െ			න ሾ݂୳୮୵ ܮ௞ሿିଵ/ଶ

ଵ/ଶ

௧೙, ೗

௧೙

	ݐ݀	ߦ݀ ሻక௞ܮሺ ݂ න න ൅ ݐ݀

ூ

௧೙, ೗

௧೙

.	

(3.18)

For the above last two terms, integrate in (normalized) time from 0 to ܿ௟ via the right Radau quadrature 

(3.14), we obtain the following surface and volume integrals: with Legendre index ݇ and time index ݈, 

SI௞

௟ ൌ െ	Δݐ ෍ܽ௠௟ሾ݂୳୮୵ሺݔ௝ାଵ/ଶ, ݐ௡,௠ሻ ܮ௞ሺݔ௝ାଵ/ଶሻ െ ݂୳୮୵ሺݔ௝ିଵ/ଶ, ݐ௡,௠ሻ	ܮ௞ሺݔ௝ିଵ/ଶሻሿ

௦

௠ୀଵ

(3.19)

and 

VI௞

௟ ൌ Δݐ ෍	ܽ ௠௟න ݂ሺߦ, ݐ௡,௠ሻሺܮ௞ሻకሺߦሻ݀ߦ

ூ

௦

௠ୀଵ

. (3.20)

That is, 

௞ ,௝ݑ	ଶ௞‖ܮ‖௝݄

௞ ,௝ݑ ଶ௞‖ܮ‖௝݄ ൌ ݈ ݊,

݊ െ SI௞

௟ ൅ VI௞௟ . (3.21)

t

x

t

x

ோݑ ௅ݑ

ோݑ ௅ݑ

ோݑ ௅ݑ

௡ݐ

ାଵ௡ݐ

௡ݐ

ାଵ௡ݐ
NASA/TM—2013-218077 13 

The solution we are seeking corresponds to time level ݈ൌݏ and Legendre indices ݇ ൌ 0, ... , ݌. 

The surface integral is straightforward. The volume integral is where the difficulty lies. This integral 

must be carried out in a manner that two contrasting requirements hold: (a) the procedure can be extended 

with relative ease to the case of systems but (b) the space-time discontinuity in Fig. 3.1 is somehow 

accounted for. 

What comes to the rescue here is the fact that the degree of ሺܮ௞ሻక is one lower than that of ܮ ;௞as a 

result, we can successively calculate 	ݑ௝,଴ 

ଵ ,௝ݑ	 ,௟ ,௡ଶ ,௝ݑ	 ,௟ ,௡௣ ,௝ݑ	 and , ..., ௟,௡

௡, ௟ as follows. 

First, for ݇ൌ0, the volume integral vanishes since ሺܮ଴ሻక ൌ 0. This fact and (3.21) yield 	ݑ௝,଴ 

௡, ௟ for all ݈. 

଴ ,௝ݑ	 With

௡, ௟ known, we can calculate 	ݑ௝, ଵ௡, ௟ where ݇ൌ1 as follows. Since (2.6) implies ሺܮଵሻక ൌ 2, the 

volume integral VI௞

௟ of (3.20) involves only the projection of ݂ onto 2ܮ଴. Thus, we need to approximate ݂

up to ܮ ଴but no higher, which is given by ݂ሺ	ݑ௝,଴ 

௡, ௟ሻ. The solution 	ݑ௝, ଵ௡, ௟ can therefore be obtained for all ݈. 

Next, with 	ݑ௝,଴ 

௡, ௟ and 	ݑ௝, ଵ௡, ௟ known, we can calculate 	ݑ௝, ଶ௡, ௟ where ݇ൌ2 as follows. Since (2.6) implies 

ሺܮଶሻక ൌ 12ߦ ൌ 6ܮଵ, the volume integral VI௞

௟ of (3.20) involves the projection of ݂ onto 6ܮଵ. Thus, we 

only need to approximate ݂ up to ܮଵ. With 	ݑ௝,଴ 

௡, ௟ and 	ݑ௝, ଵ௡, ௟ known, we can approximate ݑ by 

଴ ,௝ݑ	ଵ ,௝ݑ	 ൅ ሻߦሺ଴ܮ	௟ ,௡

௡, ௟	ܮଵሺߦሻ, and ݂ by 	݂௝, ଴௡, ௟	ܮ଴ሺߦሻ ൅ 	݂௝, ଵ௡, ௟	ܮଵሺߦሻ (by whatever procedure that provides ݂ to 

the same degree as ݑ(. The solution 	ݑ௝, ଶ

௡, ௟ can then be obtained for all ݈. 

Assuming that 	ݑ௝,଴ 

ଵି ௞,௝ݑ	 ,... ,௟ ,௡

௡, ௟ are known, we now calculate 	ݑ௝,௞ ௡, ௟. Since ሺܮ௞ሻక is of degree ݇െ1, it 

can be expressed as a linear combination of ܮ଴,... , ܮ௞ିଵ. In fact, in general, 

ሺܮଶ௠ାଵሻక ൌ 2෍ ܮଶ௟

௠

௟ୀ଴

and ሺܮଶ௠ାଶሻక ൌ 2෍ ܮଶ௟ାଵ

௠

௟ୀ଴

, (3.22)

i.e., 

ሺܮଶ௠ାଵሻక ൌ 2ܮଶ௠ ൅ 2ܮଶ௠ିଶ ൅ ⋯ and ሺܮଶ௠ାଶሻక ൌ 2ܮଶ௠ାଵ ൅ 2ܮଶ௠ିଵ ൅ ⋯

Thus, the volume integral VI௞

௟ by (3.20) involves projecting ݂ onto ࡼ௞ିଵ, and we only need to 

approximate ݂ up to ܮ௞ିଵ. Since the solution can be approximated by 	ݑ௝,଴ 

൅⋯൅ ሻߦሺ଴ܮ	௟ ,௡

ଵି ௞,௝ݑ	

௡, ௟ 	ܮ௞ିଵሺߦሻ, the flux ݂ can be approximated by 	݂௝, ଴௡, ௟	ܮ଴ሺߦሻ ൅ ⋯ ൅ 	݂௝, ௞ିଵ௡, ௟ 	ܮ௞ିଵሺߦሻ (again, by 

whatever procedure that provides ݂ to the same degree as ݑ (and the solution 	ݑ௝,௞ 

௡, ௟ follows for all ݈. This 

completes the algorithm. 

Note that in the calculation for 	ݑ௝,௣ 

௡, ௟, namely ݇ൌ݌, the last Legendre component ݇ for the above 

successive calculation, we only need to evaluate 	ݑ௝,௣ 

௡, ௦; the quantities 	ݑ௝,௣ ௡, ௟ where ݈ ൌ 1, ... , ݌ are not 

needed (we would need them if we use the left Radau or Gauss quadrature or if we use the current 

procedure as a building block in an iterative procedure for the implicit right Radau method). 

Also note that if ݑ is approximated by ݓ ௝ൌ 	ݑ௝,	଴ ܮ଴ሺߦሻ ൅ ⋯ ൅ 	ݑ௝,௞ ିଵ	ܮ௞ିଵሺߦሻ, there are several ways 

to approximate ݂ up to 	ܮ௞ିଵ. One way is by evaluating ݂௠ ൌ ݂ሺݓ௝ሺߦ௠ሻሻ where ߦ ௠are the ݇ Gauss 

points and then projecting the polynomial interpolating ሼ݂௠ሽ onto ࡼ௞ିଵ. 

For the case of advection, the current moment approach is equivalent to the exact evolution approach 

of the previous section. This fact has also been verified by the author via a Mathematica program. 
NASA/TM—2013-218077 14 

For the piecewise linear case, the moment scheme is simple: ݓ௝

ேூሺݔ, ݐሻ is a linear function where, with 

ߦ on ሾെ1/2, 1/2ሿ and ߬ on ሾ0, 1ሿ, 

௝ݓ

 ߬.ߪଵ ,௝ݑ2 െ ߦଵ ,௝ݑ2 ൅ ଴ ,௝ݑ ൌ ሻ߬ ,ߦேூሺ

The extension of the piecewise linear case to the Euler equations is also simple (Huynh 2006). 

Algorithm. The moment scheme is summarized below. Assume that the data ݑ݆, ݇ at time ݐ ௡are given 

for all ݆ and ݇ ൌ 0, ... , ݌. 

(1) In each space-time cell ܧ ௝ൈ ሾݐ௡, ݐ௡ାଵሿ, obtain ݓ௝

ேூ with no interaction by a Cauchy-Kovalevsky 

procedure or an explicit RK with no interaction. 

(2) At each interface ݔ௝ିଵ/ଶ, for each right Radau intermediate time level, calculate the left and right 

values using ݓ௝ିଵ

ேூ and ݓ௝ேூ. Obtain upwind fluxes from these two values. 

(3) Successively evaluate ݑ௝,଴ 

ଵ ,௝ݑ ,௟ ,௡ଶ ,௝ݑ ,௟ ,௡௣ ,௝ݑ and , ..., ௟,௡

௡, ௟ via (3.21) by calculating the surface integral 

via (3.19) and volume integral via (3.20). This completes the algorithm. 

We conclude this section by the following remarks. 

If the upwind fluxes at the interfaces are given (e.g., in an iterative procedure for the implicit right 

Radau time stepping), the above algorithm provides the solution for the space-time cell ܧ ௝ൈ ሾݐ௡, ݐ௡ାଵሿ. 

This observation shows that the moment scheme can serve as a key building block in an iterative 

procedure for the implicit time-stepping method. 

One way to apply the moment approach to the nodal type of schemes is to transform the nodal values 

into the modal (Legendre) components, carry out the moment calculations, and then transfer back. The 

various transferring, however, involve some costs. The question is whether it is possible to use an 

iteration procedure in the nodal framework with no transferring. 

4 Moment Scheme for the Case of Two Spatial Dimensions 

For the case of two spatial dimensions, the approach of exact time evolution and projection of Section 

2 can be extended in a manner that if the flow is along the diagonal direction, say, from southwest to 

northeast, then the solution recovers the exact data of the southwest cell for the CFL number ൫ߪ௫, ߪ௬൯ ൌ

ሺ1,1ሻ. Such an approach must allow for interaction among cells that share a common edge as well as cells 

that share a common corner. Its key drawback is that extension to systems is extremely difficult. 

Therefore, the extension of the moment scheme here allows interaction only among immediate neighbors, 

i.e., cells that share a common edge. The price to pay for this simplification is a reduction in CFL 

condition when the flow is not along the edge direction. 

To proceed, with ݑ ൌ ݑሺݔ, ݕ, ݐሻ, the 2D conservation law takes the form 

ݑ௧ ൅ ݂௫ ൅ ݃௬ ൌ 0.

Let initial condition be ݑሺݔ, ݕ, 0ሻ ൌ ݑ୧୬୧୲ሺݔ, ݕሻ	. For the Euler equations, ݑ, ݂, and ݃ are vectors of four 

entries. With no loss of generality, we consider only the scalar case. Assume that the solution ݑ is 

periodic or of compact support so that boundary conditions are trivial. 
NASA/TM—2013-218077 15 

Denote ࢌ ൌ ሺ݂,  ݃ሻ and સ ൌ ቀ డ

డ௫ , 	 డడ௬ቁ. Then 	સ ∙ ࢌ ൌ ݂௫ ൅ ݃௬ , and the above can be written as 

ݑ௧ ൅ સ∙ ࢌ ൌ 0. (4.1)

Let the domain of calculation be divided into a non-overlapping unstructured mesh of cells ܧ ௝which 

are quadrilaterals or triangles. On each cell ܧ௝, let the solution ݑ	be approximated by ݑ ௝of degree ݌. If the 

cell is a quadrilateral, we can use either a polynomial in ߦ and ߟ of degree ݌ or less, or we can use tensor 

products of degree ݌ൈ݌. For example, if ݌ൌ1, we can use either a linear polynomial ܽ ൅ ܾߦ ൅ ܿߟ, or 

for the case of tensor products, a bilinear polynomial ܽ ൅ ܾߦ ൅ ܿߟ ൅ ݀ߟߦ. 

Using ݑ௝, the fluxes ݂ and ݃ can be approximated by ݂௝ and ݃௝ of the same degree and type as ݑ௝. The 

piecewise polynomial functions ሼݑ௝ሽ, ሼ݂௝ሽ, and ሼ݃௝ሽ can be and usually are discontinuous across cell 

interfaces. 

From here on, unless otherwise stated, we focus only on the approximate polynomials. For simplicity 

of notation, the subscript ݆ is understood: ݑ௝, ࢌ ௝, and ܧ ௝are abbreviated to ݑ, ࢌ, and ܧ when there is no 

ambiguity. 

Let ߶ be a test function on ܧ, i.e., a polynomial of the same type as ݑ. We require the solution to 

satisfy, for the moment, formally, 

ܧ݀ ߶ ௧ݑ	න

ா

൅ න ሺસ ∙ ࢌሻ߶ ݀ܧ

ா

ൌ 0.

As in the 1D case, integrating by parts and applying the upwind flux for the boundary term to account for 

data interaction, the solution is required to satisfy 

ܧ݀	߶ ௧ݑ	න

ா

൅ න ሺࢌ ∙ ࢔ሻ୳୮୵ ߶ ݀ݏ

డா

െ න સ߶ ∙ ࢌ ݀ܧ

ா

ൌ 0. (4.2)

Integrating in time from ݐ ௡to ݐ ∗ൌ ݐ௡,௟, 

ܧ݀	߶	ሻ∗ݐ ,ݕ ,ݔሺݑන

ா

ܧ݀ ߶ ሻ௡ݐ ,ݕ ,ݔሺݑන െ

ா

െ SI െ VI ൌ 0	 (4.3)

where the surface and volume integral are, respectively, 

SI	 ൌ െ න න ሺࢌ ∙ ࢔ሻ୳୮୵ ߶ ݀ݏ

డா

௧∗

௧೙

(4.4) ݐ݀

and 

VI ൌ න න સ߶ ∙ ࢌ ݀ܧ

ா

௧∗

௧೙

(4.5) .ݐ݀

The problem reduces to the evaluation of SI and VI. 

To this end, as in the 1D case, we first obtain ݓ௝

ேூ with no interaction by a Cauchy-Kovalevsky 

procedure or an explicit RK with no interaction up to appropriate order. For the case of tensor products, 

we would need to derive time derivatives of degree up to 2݌ if the Cauchy-Kovalevsky procedure is 
NASA/TM—2013-218077 16 

employed. Next, along each edge, for each right Radau intermediate time level, calculate the left and right 

values using appropriate ݓ௝

ேூ. Obtain normal upwind fluxes from these two values. 

The surface integral again poses no difficulty. As for the volume integral, we make use the fact that 

both components of સ߶, namely, ߶௫ and ߶௬ are of at least one degree lower than ߶. As in the 1D case, 

we can evaluate the volume integral and successively calculate the Legendre components ݑ௝,଴ ,଴

଴,ଵ ,௝ݑ , ௟ ,௡

௡, ௟ , 

ଵ, ଴,௝ݑ଴,ଶ ,௝ݑ , ௟ ,௡

௡, ௟ , ݑ௝, ଵ,ଵ

ଶ, ଴,௝ݑ , ௟ ,௡

௡, ௟ , ... This completes the algorithm. 

5 Fourier Analyses 

One-Dimensional Case. For the 1D case, consider the advection equation with positive speed ܽ, 

ݑ௧ ൅ ܽݑ௫ ൌ 0.

The cells are ܧ ௝ൌ ሾ݆ െ 1/2, ݆ ൅ 1/2ሿ. Denote by ࢁ ௝the column vector of ܭൌ݌൅1 components, 

൯௣,௝ݑ,...,ଵ,௝ݑ ,଴,௝ݑ൫ ൌ௝ ࢁ

்

where ܶ represents the transpose. Let Δݐ be the time step and set ߪ ൌ ܽΔݐ. Assume that the data ࢁ௝

௝ࢁ ൌ௡ 

are known. The solution ࢁ௝

௡ାଵ can be expressed as 

௝ࢁ

(5.1) ଵି௝ࢁଵି࡯ ൅ ௝ࢁ଴࡯ ൌ ାଵ௡

where ࡯ ଴and ࡯ିଵ are ܭൈܭ matrices. As examples, for ݌ൌ1, i.e., the linear case, 

࡯ିଵ ൌ 	 ൬ ߪ ߪሺ1െߪሻ

െ3ߪሺ1െߪሻ െߪሺ3 െ 6ߪ ൅ 2ߪଶሻ

൰ , ࡯ ଴ൌ ൬ 1െߪ ߪሺ1െߪሻ

െ3ߪሺ1െߪሻ ሺ1 െ ߪሻሺെ1 ൅ 2ߪ ൅ 2ߪଶሻ

൰.

and, for ݌ൌ2 or the parabolic case, 

	࡯ିଵ ൌ 	ቌ

ߪ ߪሺ1െߪሻ ߪሺ1െߪሻሺ1 െ 2ߪሻ

െ3ߪሺ1െߪሻ െߪሺ3 െ 6ߪ ൅ 2ߪଶሻ െ3ߪሺ1െߪሻሺ1 െ 3ߪ ൅ ߪଶሻ

5ߪሺ1െߪሻሺ1 െ 2ߪሻ 5ߪሺ1െߪሻሺ1 െ 3ߪ ൅ ߪଶሻ ߪሺ5 െ 30ߪ ൅ 50ߪଶ െ 30ߪଷ ൅ 6ߪସሻ

ቍ,		

and 

࡯ ଴ൌ ቌ

1 െ ߪ െߪሺ1െߪሻ െߪሺ1െߪሻሺ1 െ 2ߪሻ

3ߪሺ1െߪሻ ሺ1െߪሻሺ1 െ 2ߪ െ 2ߪଶሻ െ3ߪሺ1െߪሻሺ1െߪെߪଶሻ

െ5ߪሺ1െߪሻሺ1 െ 2ߪሻ 5ߪሺ1െߪሻሺ1െߪെߪଶሻ ሺ1െߪሻሺ1 െ 4ߪ െ 4ߪଶ ൅ 6ߪଷ ൅ 6ߪସሻ

ቍ.

Note the highest degree term for ݌ൌ1 is ߪଷ and, that for ݌ൌ2, ߪହ. 

Stability. Let the wave number be denoted by ݓ,– ߨ൏ ݓ൑ ߨ, and the imaginary number, by I

instead of the standard notation ݅ to avoid confusion with the cell index ݅ later. Assuming the solution is a 

harmonic such that ࢁ௝ିଵ ൌ ݁ିࡵ௪ࢁ௝, we obtain the amplification matrix 

(5.2) ଵି࡯௪ࡵ݁ି ൅ ଴࡯ ൌ ࡭

which satisfies, by (5.1), 
NASA/TM—2013-218077 17 

௝ࢁ

௡ାଵ ൌ ࢁ࡭௝.) 5.3(

With each value of ߪ and ݓ, the corresponding matrix ࡭ has ݌൅1ൌܭ complex eigenvalues. For 

stability, these eigenvalues must have magnitude no larger than 1. Numerical calculations for these 

eigenvalues for all ݌ tested by this author result in stable schemes provided that 0൑ߪ൑1. 

Figures 5.1 show the magnitude of the amplification factors (or eigenvalues) for the case ݌ൌ1 and 

Figs 5.2, the case ݌ ൌ 2.

(a) Principal (b) Spurious 

Figs. 5.1 Magnitude of the two amplification factors (or eigenvalues) for the case ݌ൌ1. 

(a) Principal (b) Spurious (c) Spurious 

Figs. 5.2 Magnitude of the three amplification factors for the case ݌ൌ2. 

Accuracy. To calculate the order of accuracy, denote the principal eigenvalue of ࡭ by ݁௣ ൌ ݁௣ሺߪ, ݓሻ. 

This value approximates ݁ିࡵఙ௪. A scheme is accurate to order ݉ if, for a fixed ߪ, ݁ ௣approximates ݁ିࡵఙ௪

to ܱሺݓ௠ାଵሻ for small ݓ, 

݁௣ െ ݁ିࡵఙ௪ ൌ ܱሺݓ௠ାଵሻ. (5.4)
NASA/TM—2013-218077 18 

In practice, it is difficult if not impossible to derive a Taylor series expression for ݁௣ when ݌൐1. 

Therefore, we obtain the order of accuracy of a scheme by a numerical calculation. First, set ߪൌߪ ଴ൌ

0.8, say, and ݓൌݓ ௖ൌ ߨ/4. We can calculate the coarse mesh error 

er௖ ൌ ݁௣ሺߪ ଴, ݓ௖ሻ െ݁ିࡵఙబ௪೎. (5.5)

By halving the wave number, ݓ ௙ൌ ݓ௖/2 ൌ ߨ/8) equivalent to doubling the number of mesh points), the 

corresponding fine mesh error is 

er௙ ൌ ݁௣ሺߪ ଴, ݓ௙ሻ െ݁ିࡵఙబ௪೑. (5.6)

For a scheme to be m-th order accurate, after one time step, 

ቤ

er௖

er௙

ቤ ൎ 2௠ାଵ. (5.7)

That is, 

݉ ൎ

Log ൬ฬer௖

er௙

ฬ൰

Logሺ2ሻ െ 1. (5.8)

Thus, for a scheme to be of order ݉, when we march to a certain final time, doubling the mesh results in 

an error smaller by a factor of 2௠. 

It can be verified, as will be shown in the more involved 2D cases, that a moment scheme using 

polynomials of degree ݌ is accurate to order 2݌ ൅ 1. 

Two-Dimensional Case. We next briefly describe Fourier (Von Neumann) analysis for the 2D case. 

On the domain (,) (,), consider the advection equation 

ݑ௧ ൅ ܽݑ௫ ൅ ܾݑ௬ ൌ 0.

The cells are the squares ܧ௜, ௝ൌ ሾ݅ െ 1/2, ݅ ൅ 1/2ሿ ൈ ሾ݆ െ 1/2, ݆ ൅ 1/2ሿ centered at ሺ݅, ݆ሻ. For the case of 

polynomials of degree ݌, the number of Legendre basis functions is ܭ ௣ൌ ሺ݌ ൅ 1ሻሺ݌ ൅ 2ሻ/2. For the case 

of tensor products, the number of Legendre basis functions again denoted by ܭ௣, is ܭ ௣ൌ ሺ݌൅1ሻଶ. In 

addition, denote by ࢁ௜, ௝the column vector of ܭ ௣components, 

ࢁ௜, ௝ൌ ൫ݑ௜,௝,଴,଴, ݑ௜,௝,ଵ,଴, ݑ௜,௝,଴,ଵ,...൯் (5.9)

As an example, for ݌ൌ1, 

ଵ൯,଴,௝,௜ݑ ,଴,ଵ,௝,௜ݑ ,଴,଴,௝,௜ݑ൫ ൌ௝ ,௜ࢁ

்

. 

Or, for the case of tensor products, 

ࢁ௜, ௝ൌ ൫ݑ௜,௝,଴,଴, ݑ௜,௝,ଵ,଴, ݑ௜,௝,଴,ଵ, ݑ௜,௝,ଵ,ଵ൯

் . 
NASA/TM—2013-218077 19 

Set ߪ௫ ൌ ܽΔݐ and ߪ௬ ൌ ܾΔݐ. The 2D moment solution can be written as 

௝,௜ࢁ

(5.10) .ଵି௝,௜ࢁଵ,ି଴࡯ ൅ ௝,ଵି௜ࢁ଴,ଵି࡯ ൅ ௝,௜ࢁ଴,଴࡯ ൌ ାଵ௡

Here, ࡯଴,଴,	 ࡯ିଵ,଴, and ࡯଴,ିଵ are ܭ ௣ൈ ܭ ௣matrices whose entries involve ߪ௫ and ߪ௬. As an example, for 

the linear case (5.9), 

࡯ିଵ,଴ ൌ ቌ

ߪ௫ ߪ௫ሺ1െߪ௫ሻ െߪ௫ߪ௬

െ3ߪ௫ሺ1െߪ௫ሻ െߪ௫ሺ3 െ 6ߪ௫ ൅ 2ߪ௫

ଶሻ ߪ௫ߪ௬ሺ3 െ 2ߪ௫ሻ

3ߪ௫ߪ௬ ߪ௫ߪ௬ሺ3 െ 2ߪ௫ሻ ߪ௫ሺ1 െ 2ߪ௬

ଶሻ

ቍ	,

࡯଴,ିଵ ൌ ൮ ߪ௬ െߪ௫ߪ௬ ߪ௬൫1 െ ߪ௬൯

3ߪ௫ߪ௬ ߪ௬ሺ1 െ 2ߪ௫

ଶሻ ߪ௫ߪ௬ሺ3 െ 2ߪ௬ሻ

െ3ߪ௬൫1 െ ߪ௬൯ ߪ௫ߪ௬൫3 െ 2ߪ௬൯ െߪ௬൫3 െ 6ߪ௬ ൅ 2ߪ௬

ଶ൯

൲	,			

and 

࡯଴, ଴ൌ ቌ

1െߪ௫ െ ߪ௬ ߪ௫ሺെ1 ൅ ߪ௫ ൅ ߪ௬ሻ ߪ௬ሺെ1 ൅ ߪ௫ ൅ ߪ௬ሻ

െ3ߪ௫ሺെ1 ൅ ߪ௫ ൅ ߪ௬ሻ 1 െ 3ߪ௫ ൅ 2ߪ௫

ଷ െ ߪ௬ ൅ 2ߪ௫ଶߪ௬ 2ߪ௫ߪ௬ሺߪ௫ ൅ ߪ௬ሻ

െ3ߪ௬ሺെ1 ൅ ߪ௫ ൅ ߪ௬ሻ 2ߪ௫ߪ௬ሺߪ௫ ൅ ߪ௬ሻ 1െߪ௫ െ 3ߪ௬ ൅ 2ߪ௫ߪ௬

ଶ ൅ 2ߪ௬ଷ

ቍ.

Stability. Let the wave number in the ݔ-direction be denoted by ݓ௫ and in the ݕ-direction, ݓ௬

(involving no partial derivative), and let the imaginary number be denoted by I. Similar to the 1D case, 

replacing ࢁ௜ିଵ,௝ and ࢁ௜,௝ିଵ respectively by ݁ିࡵ௪ೣࢁ௜, ௝and ݁ିࡵ௪೤ࢁ௜,௝, we obtain the amplification matrix 

(5.11) .ଵ,ି଴࡯೤௪ࡵ݁ି ൅ ଴,ଵି࡯ೣ௪ࡵ݁ି ൅ ଴,଴࡯ ൌ ࡭

Consider a fixed value ሺߪ௫, ߪ௬ሻ. For each wave number ሺݓ௫, ݓ௬ሻ where –ߨ൏ ݓ௫ ൑ ߨ and –ߨ൏ ݓ௬ ൑ ߨ, 

the corresponding matrix ࡭ has ܭ ௣eigenvalues. For stability, these eigenvalues must have magnitude no 

larger than 1 for all wave numbers. The domain of all values ሺߪ௫, ߪ௬ሻ satisfying this stability condition is 

called the stability region of the scheme. 

Figure 5.3 shows the stability region for a square mesh and polynomials of degree ݌ ൌ 1, ... , 4. As is 

well-known, the CFL limits along the diagonal for the first-order upwind scheme (݌ൌ0) is √2/2 ൎ

0.7071. In the case of the moment schemes, for ݌ ൌ 	1, ... ,4, they are, respectively, 

0.4714, 0.2859, 0.2190, and 0.1654. 
NASA/TM—2013-218077 20 

(a) ݌ൌ1 (b) ݌ൌ2

(c) ݌ൌ3 (d) ݌ൌ4

Figure 5.3. Stability regions for schemes using polynomials of degree ݌. Note that for ݌ൌ3 and 4, 

the stability conditions along the ݔ and ݕ directions are 1, but these conditions reduce very quickly when 

the flow slightly deviates from the axes directions. 

Figure 5.4 shows the stability region for a square mesh and polynomials via tensor products of degree 

݌ൈ݌ where ݌ ൌ 1,	 2, 3, and 4. The CFL limits along the diagonal for ݌ ൌ 1, ... ,4 are respectively 

0.4714, 0.2834, 0.2526, and 0.1702. 
NASA/TM—2013-218077 21 

(a) ݌ൌ1 (b) ݌ൌ2

(c) ݌ൌ3 (d) ݌ൌ4

Figure 5.4. Stability regions for schemes using polynomials via tensor products of degree ݌ൈ݌. 

Again, for ݌ൌ3 and 4, the stability conditions along the ݔ and ݕ directions are 1, but these conditions 

reduce very quickly when the flow slightly deviates from the axes directions. 

Thus, the reduction in the CFL condition for flows along the diagonal direction is significant. The 

main reason is that information does not directly communicate among cells that share only one corner. 

That is, for a flow from, say, southwest to northeast, information flows to the cells to the east and the 

north before it can reach the cell at the northeast corner. This difficulty concerning information flow can 

be remedied by using the moment scheme as a procedure to update the solution in an iterative process for 

an implicit space-time DG scheme as will be discussed in the next section. 
NASA/TM—2013-218077 22 

Accuracy. To calculate the order of accuracy, using the notation of the 1D case, denote the principal 

eigenvalue of the amplification matrix ࡭ defined in (5.11) by ݁௣ ൌ ݁௣ሺߪ௫, ߪ௬, ݓ௫, ݓ௬ሻ. This principal 

eigenvalue approximates ݁ିࡵሺఙೣ௪ೣାఙ೤௪೤ሻ

. To find the order of accuracy of the scheme, first, fix ሺߪ௫, ߪ௬ሻ, 

say, ൫ߪ௫, ߪ௬൯ ൌ ሺ0.2, 0.1ሻ. For the coarse mesh, set ൫ݓ௫, ݓ௬൯ ൌ ሺߨ/4, ߨ/6ሻ. We can calculate the coarse 

mesh error 

er௖ ൌ ݁௣൫ߪ௫, ߪ௬, ݓ௫, ݓ௬൯െ݁ିࡵሺఙೣ௪ೣାఙ೤௪೤ሻ

. (5.12)

Next, halving the wave number, set ൫ݓ௫, ݓ௬൯ ൌ ሺߨ/8, ߨ/12ሻ. A formula similar to the above yields the 

corresponding fine mesh error er௙. For a scheme to be m-th order accurate, (after one time step) 

Log ൬ฬer௖

er௙

ฬ൰

Logሺ2ሻ െ 1 ൎ ݉. (5.8)

As can be seen in Tables 5.1 and 5.2, the 2D moment schemes using polynomials either of degree ݌ or 

via tensor products of degree ݌ൈ݌ are accurate to order 2݌ ൅ 1. 

Polynomial Degree Coarse Mesh Error Fine Mesh Error Order of Accuracy 

1 1.88 ൈ 10ିଷ 1.28 ൈ 10ିସ 2.87

2 2.58 ൈ 10ିହ 4.26 ൈ 10ି଻ 4.92

3 1.83 ൈ 10ି଻ 7.47 ൈ 10ିଵ଴ 6.94

4 7.22 ൈ 10ିଵ଴ 7.49 ൈ 10ିଵଷ 8.91

Table 5.1. Errors and order of accuracy of 2D moment schemes using polynomials of degree ݌ ;

൫ߪ௫, ߪ௬൯ ൌ ሺ0.2, 0.1ሻ, coarse mesh ൫ݓ௫, ݓ௬൯ ൌ ሺߨ/4, ߨ/6ሻ, and fine mesh ൫ݓ௫, ݓ௬൯ ൌ ሺߨ/8, ߨ/12ሻ. 

Polynomial Degree Coarse Mesh Error Fine Mesh Error Order of Accuracy 

1ൈ1 8.3 ൈ 10ିସ 5.42 ൈ 10ିହ 2.94

2ൈ2 3.09 ൈ 10ି଺ 4.83 ൈ 10ି଼ 5.

3ൈ3 9.47 ൈ 10ିଽ 3.81 ൈ 10ିଵଵ 6.96

4ൈ4 1.14 ൈ 10ିଵଵ 1.41 ൈ 10ିଵସ 8.66

Table 5.2. Errors and order of accuracy of 2D moment schemes using tensor products (degree ݌ൈ݌ ;(

൫ߪ௫, ߪ௬൯ ൌ ሺ0.2, 0.1ሻ, coarse mesh ൫ݓ௫, ݓ௬൯ ൌ ሺߨ/4, ߨ/6ሻ, and fine mesh ൫ݓ௫, ݓ௬൯ ൌ ሺߨ/8, ߨ/12ሻ. 

Due to the larger degree of freedom, errors for schemes via tensor products are smaller than the 

corresponding errors for schemes via polynomials of degree ݌. 
NASA/TM—2013-218077 23 

6 Implicit Space-Time DG Schemes 

The moment scheme discussed in Section 3 is closely related to the implicit space-time DG method of 

this section. When applying the DG approach to discretize time, as shown by LaSaint and Raviart (1974) 

who employed a quadrature of left Radau type, the result is an implicit Runge-Kutta scheme. It was later 

shown in (Huynh 2009b) that using the FR framework, the right Radau points are a natural choice for 

quadrature points. The resulting DG method is simpler and is identical to a collocation method called 

Radau IIA (Lambert 1991, Hairer, Norsett, and Wanner 1991, Hairer and Wanner 1993). We review the 

DG method and apply the FR approach to time marching to show the close relation between the resulting 

implicit space-time method and the moment scheme. 

Consider the ODE 

ᇱݑ

ሺݐሻ ൌ ݂ሺݑ, ݐሻ, (6.1a)

with initial condition 

(b6.1 (.଴ݑሻൌ଴ݐሺݑ

Here, some of the notations are from ODE; as such, the above ݂ is not related to the flux ݂ of Section 4. 

Next, let Δݐ be the time step size and set ݐ ௡ൌ ݊Δݐ, ݊ ൌ 0, 1, 2, ...; the assumption of constant step size is 

only for convenience; in practice, time step sizes can vary. When dealing with time, we use ܫ ൌ ሾ0, 1ሿ and 

ݐൌݐ ൅ ௡߬Δݐ where ߬ varies on ܫ ;thus, ݐ varies on ܫ ௡ൌ ሾݐ௡, ݐ௡ାଵሿ. 

Assume that the data ݑ ௡at time ݐ ௡is known, we wish to calculate the solution ݑ௡ାଵ at ݐ௡ାଵ. For 

݊ൌ0, ݑ ௡is the initial condition ݑ ଴in (6.1b). 

The DG method seeks to approximate the solution for (6.1) by a function denoted by ݑ ௛which can be 

and usually is discontinuous across each ݐ௡. On ܫ௡, this function is a polynomial of degree ݏെ1, which 

is defined by its values at ݏ points at locations corresponding to the stages. Note that when applying these 

time-stepping methods to conservation laws, if ݌ is the degree of the spatial polynomial, the number of 

stages in time is 

 1.൅݌ൌܭൌݏ

Due to the property of super accuracy in time discussed later, we could combine parabolic approximation 

in space or ݌ൌ2 with linear approximation in time or ݏൌ2. Such combinations remain to be explored. 

Assume that the calculation for ܫ௡ିଵ ൌ ሾݐ௡ିଵ, ݐ௡ሿ is finished, and we have the solution ݑ ௛on that 

interval. Thus, the value just to the left of ݐ ௡is known, namely ݑ ௛ቀݐ–	௡ቁ Recall that the solution can be 

and usually is discontinuous across each ݐ௡. We assume time marches forward, thus, the value at ݐ ௡is 

that from the left (small square of Fig. 6.1), 

(6.2) .ሻି௡ ݐሺ௛ݑ ൌ௡ ݑ

We wish to calculate ݑ ௛on ܫ ௡ൌ ሾݐ௡, ݐ௡ାଵሿ, i.e., to obtain 

ݑ௛ሺݐ௡ାଵ ିሻ ൌ ݑ௡ାଵ. (6.3)
NASA/TM—2013-218077 24 

Fig. 6.1. DG solution ݑ ௛can be and usually is discontinuous across each ݐ௡. The value at ݐ ௡is chosen 

to be that from the left (small squares): ݑ ௡ൌ ݑ௛ሺݐ	௡ିሻ.

DG formulation. On ܫ ௡ൌ ሾݐ௡, ݐ௡ାଵሿ, let ߶௞ be a set of basis functions for ࡼ௄ିଵ, ݇ ൌ 1, ... , ܭ. The 

solution is approximated by 

௞߶ ௞ݑ෍ൌ௛ ݑ

௄

௞ୀଵ

.

Let ߶ be a test function, which is in ࡼ௄ିଵ; typically, ߶ is one of the basis functions ߶௞, ݇ ൌ 1, ... , ܭ. 

Again on ܫ௡, the weak form for ݑᇱ ൌ ݂ is, formally (the correct equation is (6.6) or (6.7)), 

ሺݑᇱ

, ߶ሻ ൌ ሺ݂, ߶ሻ. 

Let ݂௛ be the projection of ݂ onto ࡼ௄ିଵ. Since ߶ varies on ࡼ௄ିଵ, we have ሺ݂,߶ሻ ൌ ሺ݂௛, ߶ሻ, i.e., we can 

project first, then apply the inner product. The above implies 

ሺݑᇱ

, ߶ሻ ൌ ሺ݂௛, ߶ሻ. 

We wish to find ݑ ௛in ࡼ௄ିଵ such that for all ߶ in ࡼ௄ିଵ, again formally, 

௛ݑሺ

ᇱ , ߶ሻ ൌ ሺ݂௛, ߶ሻ. (6.4)

For the DG method, ݑ௛ሺݐ	௡ାሻ is allowed to be and is generally different from ݑ௛ሺݐ	௡ିሻ ൌ ݑ௡. By 

applying integration by parts to the left hand side above, 

ሾݑ߶௛ሿ௧೙

௧೙శభ

െ ሺݑ௛,߶ ᇱ

ሻ ൌ ሺ݂௛, ߶ሻ.

Concerning the boundary term, there is no ambiguity with test function ߶, but for ݑ௛, since it is 

discontinuous across ݐ ௡and ݐ௡ାଵ, a more precise expression is 

ሾݑ߶௛ሿ௧೙

௧೙శభ

ൌ ݑ௛ሺݐ௡ାଵିሻ߶ሺݐ௡ାଵሻ െ ݑ௛ሺݐ௡ାሻ߶ሺݐ௡ሻ.

To involve the initial condition, however, we replace ݑ௛ሺݐ௡ାሻ in the boundary term by ݑ௛ሺݐ	௡ିሻ ൌ ݑ ௡in 

(6.2). Thus, the DG solution ݑ ௛is required to satisfy, for all ߶ in ࡼ௄ିଵ, 

ݑ௛ሺݐ௡ାଵሻ߶ሺݐ௡ାଵሻ െ ݑ௛ሺݐ ௡ିሻ߶ሺݐ௡ሻ െ ሺݑ௛,߶ ᇱ

ሻ ൌ ሺ݂௛, ߶ሻ. (6.6)

t

u

u uh h

t

n

ሻି௡ݐሺ௛ݑ ൌ௡ ݑ

ݑ௡ାଵ ൌ ݑ௛ሺݐ௡ାଵିሻ

ଵି௡ݐ ାଵ௡ݐ
NASA/TM—2013-218077 25 

Applying integration by parts again to ሺݑ௛,߶ ᇱ

ሻ, 

െሾݑ௛ሺݐ	௡ିሻ െ ݑ௛ሺݐ ௡ାሻሿ ߶ሺݐ௡ሻ ൅ ሺݑ௛

ᇱ , ߶ሻ ൌ ሺ݂௛, ߶ሻ. (6.7)

Note that this time we use ݑ௛ሺݐ	௡ାሻ for the boundary term ݑ௛ሺݐ௡ሻ. 

Also note that the above, the result of integration by parts twice, is the ‘strong form’ whereas (6.6) is 

the ‘weak form’. 

Flux Reconstruction formulation. Our goal is to eliminate the test function in (6.7) so that the 

integral formulation results in a differential one. To this end, we raise the following question: can we find 

a polynomial ݃୐୆ ൌ ݃ on ܫ ௡which possesses the property that, for any ߶ of degree ܭെ1 or less, 

െ߶ሺݐ௡ሻ ൌ ሺ݃ᇱ

, ߶ሻ. (6.8)

If such a polynomial exists, we can combine ሺݑ௛

ᇱ , ߶ሻ in (6.7) with ሺ݃ᇱ

, ߶ሻ and eliminate the test function. 

From (6.7), since ݂௛ is of degree ܭെ1, we require ݃ᇱ

 to also be of degree ܭെ1; as a result, ݃ is of 

degree ܭ. Applying integration by parts to the right hand side of (6.8), we have 

െ߶ሺݐ௡ሻ ൌ ݃ሺݐ௡ାଵሻ߶ሺݐ௡ାଵሻ െ ݃ሺݐ௡ሻ߶ሺݐ௡ሻ െ ሺ݃, ߶ᇱ

ሻ. (6.9)

The above holds if ݃ satisfies, 

݃ሺݐ௡ሻ ൌ 1, ݃ሺݐ௡ାଵሻ ൌ 0, (6.10)

and, for all ߶ in ࡼ௄ିଵ, 

ሺ݃, ߶ᇱ

ሻ ൌ 0. (6.11)

Since ߶ is of degree ܭെ1, ߶ᇱ

 is of degree ܭെ2; moreover, ߶ᇱ spans ࡼ௄ିଶ as ߶ spans ࡼ௄ିଵ. Eq. (6.11) 

implies that ݃ is orthogonal to ࡼ௄ିଶ, i.e., for any polynomial ߮ of degree ܭെ2, 

ሺ݃, ߮ሻ ൌ 0. (6.12)

Orthogonality to ܭࡼെ2 provides ܭെ1 conditions; (6.10) provides the other two. 

Let the right Radau polynomial of degree ܭ be denoted by ܴோ, ௄ and defined by the conditions that 

ܴோ, ௄ሺݐ௡ሻ ൌ 1 and ܴோ, ௄ vanishes at the ܭ right Radau points on ܫ௡. See Fig. 6.2(a). We wish to show that 

݃ൌܴோ, ௄. (6.13)

That is, (6.10) and (6.12) hold with ݃ replaced by ܴோ, ௄. Indeed, (6.10) is immediate by the definition of 

ܴோ, ௄. Next, since ܴோ, ௄ is of degree ܭ, for any ߮ of degree ܭെ2, ܴோ, ௄߮ is of degree 2ሺܭ െ 1ሻ, and 

ܴோ, ௄߮ vanishes at the ܭ right Radau points. In addition, the right Radau quadrature is exact for ܴோ, ௄߮. 

Consequently, (6.12) holds, and (6.13) follows. 

The answer to the question posed for (6.8) is thus positive. Note that ݃ൌ݃୐୆, which corrects for the 

jump at the left boundary, is the right Radau polynomial (vanishing at ݐ௡ାଵ). Loosely put, the condition 

݃ሺݐ௡ሻ ൌ 1 deals with the jump at the left boundary whereas (6.12) together with ݃ሺݐ௡ାଵሻ ൌ 0 implies 

that away from the left boundary, ݃ approximates the zero function. 
NASA/TM—2013-218077 26 

We now return to our task of eliminating ߶. With ݃ defined above, ሺ݃′, ߶ሻ ൌ െ߶ሺݐ

݊ሻ. Consequently, 

by (6.7), 

ሾݑ௛ሺݐ	௡ିሻ െ ݑ௛ሺݐ ௡ାሻሿ ሺ݃′, ߶ሻ ൅ ሺݑ௛

ᇱ , ߶ሻ ൌ ሺ݂௛, ߶ሻ. (6.14)

What is crucial here is that ߶ can be factored out and canceled and, since ݑ ௡ൌ ݑ௛ሺݐ	௡ିሻ, 

ሺሾݑ ௡െ ݑ௛ሺݐ ௡ାሻሿ݃൅ݑ௛ሻ′ ൌ ݂௛. (6.15)

Therefore, ሾݑ ௡െ ݑ௛ሺݐ	௡ାሻሿ݃൅ݑ ௛is a polynomial of degree ܭ whose value at ݐ ௡ൌ ݐ௡ା is, by (6.10), 

ሾݑ௛ሺݐ	௡ିሻ െ ݑ௛ሺݐ ௡ାሻሿ݃ሺݐ௡ሻ ൅ ݑ௛ሺݐ ௡ାሻ ൌ ݑ௛ሺݐ ௡ିሻ ൌ ݑ௡.) 6.16(

Thus, to deal with the jump ሾݑ௛ሺݐ	௡ିሻ െ ݑ௛ሺݐ	௡ାሻሿ ൌ ݑ ௡െ ݑ௛ሺݐ	௡ାሻ at ݐ௡, we set 

ܷ ൌ ሾݑ௛ሺݐ ௡ିሻ െ ݑ௛ሺݐ ௡ାሻሿ݃൅ݑ௛.) 6.17(

See Fig. 6.2(b). Differentiating ܷ and employing (6.15), we obtain 

ܷᇱ

ሺݐሻ ൌ ݂௛ሺݐ, ݑ௛ሺݐሻሻ. (6.18)

The function ܷ is the key idea of FR. On ሾݐ௡, ݐ௡ାଵሿ, assume that ݑ ௛which does not match the 

boundary condition ݑ ௡is given. We can add the correction term ሾݑ ௡െ ݑ௛ሺݐ	௡ାሻሿ݃ to ݑ ௛as in (6.17) so 

that the resulting ܷ takes on the left boundary value ݑ௡, does not alter the value ݑ௛ሺݐ௡ାଵିሻ ൌ ݑ௡ାଵ, and 

retains the property of ‘best possible’ approximation to ݑ௛. The function ܷ is continuous across cell 

boundaries (as opposed to ݑ ௛which is discontinuous). To put it differently, ݑ ௛is of degree ݌ ;ܷ is of 

degree ݏൌ݌൅1 determined by ݌൅2 conditions; two conditions are known, namely, ܷሺݐ௡ሻ ൌ ݑ ௡and 

ܷሺݐ௡ାଵሻ ൌ ݑ௛ሺݐ௡ାଵ	ିሻ; the remaining ݌ൌݏെ1 conditions are given by the requirement that the 

projection of ܷ and ݑ ௛onto ࡼ௣ିଵ are the same. This last requirement is the result of the fact that ݃ is 

orthogonal to ݌ࡼെ1 by (6.12).

Denoting the ܭ right Radau points on ܫ ௡by ݐ௡,௞, ݇ ൌ 1, ... , ܭ where ݐ௡, ௄ൌ ݐ௡,௦ ൌ ݐ௡ାଵ. See also Fig. 

3.3. Since ݃ vanishes at these points, (6.17) implies 

ܷ൫ݐ௡,௞൯ൌݑ௛ሺݐ௡,௞ሻ. (6.19)

By employing these Radau points as collocation points (or quadrature points), (6.18) leads to 

ܷᇱ

ሺݐ௡,௞ሻ ൌ ݂ሺݐ௡,௞, ܷ௡,௞ሻ. (6.20)

Here, we approximate ݂௛, which is of degree ݌ൌܭെ1, by using the ܭ values ݂ሺݐ௡,௞,	 ܷ௡,௞ሻ, ݇ ൌ

 .ܭ , ... 1,

In summary, to calculate the DG solution, we first obtain the polynomial ܷ by solving the implicit 

equations (6.20); here, ܷ is of degree ܭൌݏ interpolating the ܭ൅1 values at ݐ ௡and at the ܭ right Radau 

points: ܷሺݐ௡ሻ ൌ ݑ ௡and ܷ൫ݐ௡,௞൯ൌܷ௡,௞, ݇ ൌ 1, ... , ܭ. The DG solution ݑ ௛of degree ܭെ1 is defined by 

the values ܷ at these ܭ right Radau points. 

The scheme defined by (6.20) is identical to an implicit collocation Runge-Kutta scheme using the 

right Radau points and accurate to order 2݌ ൅ 1 ൌ 2ܭ െ 1 called Radau IIA (Lambert 1991, Hairer, 

Norsett, and Wanner 1991). 
NASA/TM—2013-218077 27 

(a) Radau polynomial ܴோ, ଷ (b) ܷ and ݑ௛

Figure 6.2. (a) The correction function ݃, i.e., the right Radau polynomial of degree ݏൌ3 on 

ܫ ൌ ሾ0,1ሿ; (b) ݑ ௛is of degree ݏെ1ൌ2 and discontinuous across ݐ ;௡ܷ is of degree ݏൌ3, continuous 

across ݐ௡, and defined by: ܷሺݐ௡ሻ ൌ ݑ௡, ܷሺݐ௡ାଵሻ ൌ ݑ௛ሺݐ௡ାଵ	ିሻ, and the projections of ܷ and ݑ ௛onto 

ࡼ௦ିଶ (providing ݏെ1 conditions) are the same. 

We now discuss the relation between the current implicit time stepping and the moment scheme 

(explicit time stepping) of Section 3. For the conservation law ݑ௧ ൅ ݂௫ ൌ 0, with the same spatial 

discretization using the Legendre polynomial, the key difference between the two time stepping methods 

is the following. For the moment scheme, the Cauchy-Kovalevsky procedure, i.e., the space-time Taylor 

series expansion with no interaction denoted by ݓ௝

ேூ applies. For the implicit scheme, this expansion must 

be replaced by the solution ݑ௝,௞ 

௡, ௟ itself; recall that the Legendre components correspond to ݇ ൌ 0, ... , ݌, 

and the time stages, ݈ ൌ 1, ... , ݏ. As such, the moment method of updating the solution can be employed in 

an iterative manner to calculate the implicit solution. This subject, however, is beyond the scope of the 

current paper. 

7. A Limiting Technique 

As is well-known, high-order methods generate oscillations near a shock or a discontinuity. One way 

to avoid oscillations is to limit the data so that it is not ‘too steep’ before updating. The problem with 

limiting is that it can cause a loss of accuracy near extrema. We explore ideas that that can suppress 

oscillations near shocks while preserving accuracy near extrema here. The limiting procedure below is 

designed for the piecewise linear case with one spatial dimension. The procedure employs the techniques 

introduced in (Huynh 1995) and (Suresh and Huynh 1997) but has a stencil of only three cells and is 

simplified by requiring that the magnitude of the slope lies in a certain range. 

To proceed, assume that ݑ௝, ଴ and ݑ௝, ଵ are given for all ݆. We abbreviate the cell average values: 

 .௝ݑ ൌ଴ ,௝ݑ

As discussed after (2.9), the first moment ݑ௝, ଵ represents half of the normalized slope so that the value at 

the right interface of cell ݆ is ݑ௝, ൅ ଴ ݑ௝, ଵ. Next, set 

ݏ ௅ൌ	ݑ ௝െ ݑ௝ିଵ, ݏோ ൌ ݑ௝ାଵ െ ݑ௝,) 7.1(

and 

௡ݐ

ܷ

௛ݑ

ାଵ௡ݐ

௡ݑ
NASA/TM—2013-218077 28 

	ݏ୐୧୫ଵ ൌ 2 minሺ|ݏ|௅,| ݏோ|ሻ. (7.2)

The limit ݏ୐୧୫ଵ for the normalized slope is similar to Van Leer's limit (1977) except that we use the 

magnitudes. It can suppress oscillations near a discontinuity, but has the drawback of causing a loss of 

accuracy near extrema. With 2หݑ௝, ଵห representing the normalized slope, if 

2หݑ௝, ଵห൑ݏ୐୧୫ଵ, (7.3)

limiting does not alter the slope, and we move on to the next ݆ index. At most of the smooth regions of the 

flow, the above condition is satisfied; thus, limiting is typically needed only for a small number of cells. 

If (7.3) does not hold, the cell is near a discontinuity or an extremum, and the following calculations 

are carried out. Our goal is to enlarge the limit in a manner that near a shock, it reduces to (7.1), but near 

an extrema, it is larger and provides enough ‘room’ so that the original slope is not altered. To put it 

differently, at smooth regions, the results by the moment or DG schemes are highly accurate and should 

not be altered. 

As will be shown by numerical examples, the moment scheme produces, for each cell, a slope that is 

downwind biased if ߪ is near 0, and upwind biased is ߪ is near 1. Therefore, we take extra care to account 

for this slope variation. 

For any real numbers ݔ and ݕ, let minmodሺݔ, ݕሻ be the median of ݔ, ݕ, and 0, 

minmodሺݔ, ݕሻ ൌ 1

2

ሾsignሺݔሻ ൅ signሺݕሻሿ minሺ|ݔ|,| ݕ|ሻ.

At interface ݆ െ 1/2, set 

ݒଵଵ ൌ ൫ݑ ൅ ௝ݑ௝ିଵ൯/2,							ݒଵଶ ൌ ݑ௝ିଵ ൅ ݑ௝ିଵ, ଵ െ ݒଵଵ, ݒଵଷ ൌ ݑ ௝െ ݑ௝, ଵ െ ݒଵଵ. (7.4)

Loosely put, we use ݒଵଵ as a pivot; ݒଵଶ is from the left, and ݒଵଷ, the right; typically, near a smooth 

extrema, ݒଵଶ and ݒଵଷ are of the same sign and relatively close to each other. Set 

ݒଵସ ൌ minmodሺݒଵଶ, ݒଵଷሻ. (7.5)

At interface ݆ ൅ 1/2, set 

ݒଶଵ ൌ ൫ݑ ൅ ௝ݑ௝ାଵ൯/2,							ݒଶଶ ൌ ݑ௝ାଵ െ ݑ௝ାଵ, ଵ െ ݒଶଵ, ݒଶଷ ൌ ݑ ൅ ௝ݑ௝, ଵ െ ݒଶଵ, 

and 

ݒଶସ ൌ minmodሺݒଶଶ, ݒଶଷሻ. (7.6)

Near a smooth extremum, ݒଵସ and ݒଶସ are typically of the same sign and relatively close to each other. 

For the cell ݆, set 

ݒ୫୫ ൌ |minmodሺݒଵସ, ݒଶସሻ|, ݒ୫୶ ൌ maxሺ|ݒଵସ|, |ݒଶସ|ሻ. (7.7)

Thus, ݒ୫୫ and ݒ୫୶ are also relatively close to each other near a smooth extremum. Near a discontinuity, 

however, they are far apart. With ‘ra’ short for ratio and ‘ap’ for accuracy preserving, set 
NASA/TM—2013-218077 29 

ݒ୰ୟ ൌ ݒ୫୫

ݒ୫୶ ൅ 10ିଶ଴ , ݒୟ୮ ൌ 6ݒ୫୫ݒ୰ୟ . (7.7)

Note that ݒ୰ୟ is a non-dimensional quantity. Near a discontinuity, it is close to 0; near an extremum, it is 

relatively close to 1. Therefore, ݒୟ୮ provides ‘room’ for the interface value near a smooth extremum. 

Next, set 

	ݏ୫୧୬ ൌ minሺ|ݏ|௅,| ݏோ|ሻ , 				ݏ୐୧୫ଶ ൌ ݏ୫୧୬ ൅ 2ݒୟ୮, ݏ୐୧୫ଷ ൌ maxሺݏ୐୧୫ଵ, 	ݏ୐୧୫ଶሻ.	 (7.8)

Note the factor 2 accounts for transferring a quantity of type ݒ to type ݏ. Finally, limit ݑ௝, ଵ by 

ݑ௝, ଵ 	 ← 	 signሺݑ௝, ଵሻ min൫|ݑ௝, ଵ|, ݏ୐୧୫ଷ/2൯. (7.9)

Note that if (7.3) holds and we still apply the above limiting procedure, the result recovers the original 

ݑ௝, ଵ. Thus, (7.3) only serves the purpose of saving computing time; whether it is applied or not, the final 

solution is the same. 

8. Numerical Results 

The following preliminary numerical tests show the results by the piecewise linear moment (or upwind 

moment) scheme. Recall that for 1D advection, the scheme is identical to Van Leer’s scheme 3. 

The first test is for an advection with constant speed 1. The initial condition is a sine wave on the 

domain ሾ0, 1ሿ with 6 cells and periodic boundary conditions. The continuous (red) curve represents the 

exact solution, the dots the cell average numerical solutions, and the line segments, the piecewise linear 

solutions. Figs. 8.1 and 8.2 show the following behavior of the scheme: for each CFL number, the scheme 

has a ‘preferred mode’ or eigenvector. This eigenvector corresponds to a piecewise linear function with 

slopes that are downwind biased if ߪ is close to 0, and upwind biased if ߪ is close to 1. 

Figs. 8.1 shows the results with ߪ ൌ 0.1. For small ߪ, the solution ‘likes’ downwind biased slopes in 

the sense that no matter what kind of slopes we start with, the solution will turn them into downwind 

biased slopes after a few time steps. In (a), for each cell ݆, we start with slopes biased in the ‘wrong’ 

direction, i.e., upwind biased slopes defined by the exact interface value ݑ௝ିଵ/ଶ and exact cell average 

value ݑ௝. Fig. 8.1(c) shows the solution after 10 time steps; the wave travels a distance of one cell. At the 

maximum, the slope is now slanted downward, opposite to the initial slope in (a). In fact, in (b), after just 

4 time steps, the solution already exhibits this ‘preferred mode’ for the slopes: downwind biased for small 

ߪ. Fig. 8.1(d) shows the solution after 60 time steps; the wave advects 1 period. Comparing (a) and (d), 

again note the switch in bias of the slopes. 

Figs. 8.2 shows the results with a CFL number of 0.9677, which is the exact opposite compared with 

Figs. 8.1. For CFL number close to 1, the solution ‘likes’ upwind biased slopes. We again start with 

slopes biased in the ‘wrong’ direction, i.e., downwind biased slopes defined by the exact cell average 

value ݑ ௝and exact interface value ݑ௝ାଵ/ଶ as shown in Fig. 8.2(a). After 31 time steps, the wave travels a 

distance of 30 cells, i.e., the profile advects 5 periods. The solution is shown in Fig. 8.2(b), where the 

slopes are upwind biased. 

The above property of biasing toward the upwind or downwind sides holds true for higher-order as 

well (݌൐1(. In that case, the ݌൅1 right Radau and left Radau points come into play in a manner similar 

to the linear case above. (The upwind biased linear functions in Fig.8.1(a) interpolate the values at the two 

left Radau points for each cell, and the downwind biased ones in Fig.8.2(a), the values at the two right 

Radau points.) 
NASA/TM—2013-218077 30 

(a) Data (0 time steps) (b) Solution after 4 time steps 

(c) Solution after 10 time steps (d) Solution after 60 time steps 

Fig. 8.1. CFL number ߪ ൌ 0.1. The continuous (red) curve represents the exact solution. (a) Piecewise 

linear data with slopes biased in the ‘wrong’ direction, i.e., upwind biased (slanted upward near the 

maximum); (b) after just 4 time steps, solution exhibits a ‘preferred mode’ for the slopes, which is 

downwind biased since ߪ is close to 0; (c) solution after 10 time steps; the wave travels a distance of one 

cell; and (d) solution after 60 time steps; the wave advects 1 period. 

(a) Data (b) Solution 

Fig. 8.2. CFL number ߪ ൌ 0.9677, number of time steps 31, distance of advection 5 periods. The 

continuous (red) curve represents the exact solution. (a) Piecewise linear data with slopes biased in the 

‘wrong’ direction, i.e., downwind biased, or slanted downward near the maximum; (b) the solution 

exhibits a ‘preferred mode’ for the slopes, which is upwind biased for ߪ close to 1. 
NASA/TM—2013-218077 31 

Note that the solution in Fig. 8.1(d) is damped more than that in Fig. 8.2(b) in spite of the fact that the 

wave travels much further in the latter. The reason is that for Fig. 8.2(b), the error associated with the 

CFL number of 0.9677 is small (ߪൌ1 has no error), and the number of time steps is only half that of 

Fig. 8.1(d). In other words, the errors accumulate more in Fig. 8.1(d) due to the larger number of time 

steps. Since ߪൌ1 has no error, for the moment scheme (as opposed to the case of RK time stepping), to 

minimize error, we take time steps as large as possible. 

Note that if we use piecewise linear DG with RK time stepping, the solution behaves in a manner 

similar to that of Fig. 8.1, i.e., it ‘likes’ downwind biased slopes. The moment scheme, however, has this 

peculiar property of shifting the slopes it ‘likes’ to upwind biased when ߪ is close to 1. It appears difficult 

to devise a semi-discrete formulation with this shifting ability. 

The second test is again advection. The initial condition consists of a Gaussian, a rectangular, and a 

semi-ellipse wave on the domain ሾ0, 1ሿ with 100 cells and periodic boundary conditions. The CFL 

number is .8, and the final time is ݐൌ1, i.e., the profile advects 1 period after 125 time steps. The 

solutions by the standard upwind scheme using the Van Albada limiter and the upwind moment scheme 

using the limiting procedure in Section 7 are shown in Fig. 8.3. The dots represent the cell average values. 

Note that the standard upwind solution in (a) has a loss of accuracy (only first-order accurate) near 

extrema, and discontinuities are smeared. The upwind moment solution has no loss of accuracy near 

extrema, and discontinuities are resolved with four points. 

(a) Standard upwind scheme (b) Upwind moment scheme 

Fig. 8.3. Advection problem; 125 time steps with ߪ ൌ 0.8; distance travelled: 1 period; (a) solution by 

standard upwind scheme with Van Albada limiter; (b) solution by the moment scheme with the new 

limiter. 

The next test has both shocks and local extrema (Shu and Osher 1989). On ሾെ1, 1ሿ, a moving Mach 3 

shock interacts with sine waves in density described by 

ሺߩ, ݑ, ݌ሻ ൌ ൝ ሺ3.857, 2.629, 10.333ሻ, െ1 ൑ ݔ ൑ െ0.8

ሺ1 ൅ 0.2 sinሺ5ߨ xሻ , 0, 1ሻ, െ0.8 ൑ ݔ ൑ 1

.	 (7.9)

The final time is ݐ ൌ 0.36. Fig. 8.4 shows the results with 300 cells by standard upwind (with Van 

Albada limiter) and the upwind moment scheme. Here, for upwind moment, limiting is applied to the 

characteristic variables. The solid lines represent the solution with 3600 cells by the standard upwind 

scheme. Again, the upwind moment solution is highly accurate, and accuracy is preserved near extrema. 
NASA/TM—2013-218077 32 

Fig. 8.4. Shu and Osher’s problem. The upwind moment solution is highly accurate: the shock is well 

captured and there is no loss of accuracy near extrema. 

The last test is the (2D) oblique shock problem (Yee, Warming, and Harten 1983). The domain is 

ሾ0, 4ሿ ൈ ሾ0, 1ሿ. The boundary conditions are: left, inflow with ሺߩ, ݑ, ݒ, ݌ሻ ൌ ሺ1. , 2.9, 0. , 1/1.4ሻ; top, 

ሺ1.7, 2. 6193, െ.5063, 1.5282ሻ; bottom, solid wall; right, outflow. The rectangular uniform mesh is 

80 ൈ 20. Here, the limiter is applied in the ݔ and ݕ directions when needed (to the characteristic 

variables). Figure 8.5 shows (a) pressure distribution result by the upwind moment scheme and (b) 

pressures along ݕ ൌ 0.475 by the upwind moment and the standard upwind schemes. The solid line 

represents the exact solution. Again, the upwind moment solution is highly accurate. It appears that the 

slope update via moment calculation can ‘turn’ faster due to the compact stencil compared to the slope 

update via standard interpolation using neighboring data. For this case, compared with the standard 

(simplified) second-order upwind scheme with Van Albada’s limiter, computing cost increases by 

roughly a factor of five. 

Fig.8.5 (a) Oblique shock; pressure distribution by upwind moment scheme. 
NASA/TM—2013-218077 33 

Fig.8.5 (b) Oblique shock solutions along ݕ ൌ 0.475

4. Conclusions and Discussion 

In conclusion, two new high-order space-time schemes for conservation laws were introduced: the first 

is explicit and the second, implicit. The explicit method called the moment scheme achieves a CFL 

condition of 1 for the case of one-spatial dimension regardless of the degree of approximation. Fourier 

analyses for the one and two-dimensional cases were carried out. The moment scheme is accurate to order 

2݌ ൅ 1 if the solution is approximated by piecewise polynomials of degree ݌. The implicit method is a 

simplified but optimal version of the discontinuous Galerkin scheme applied to time discretization. It 

reduces to a collocation implicit Runge-Kutta method using the right Radau points called Radau IIA. The 

explicit and implicit schemes are closely related since they employ the same intermediate time levels, and 

the former can serve as a key building block in an iterative procedure for the latter. A limiting technique 

for the piecewise linear scheme that can suppress oscillations near discontinuities while preserving 

accuracy near extrema was presented. Preliminary numerical results were shown. 

Whereas the current findings are encouraging, further research on both time-stepping and limiting is 

needed. 

References 

P. Bar-Yoseph and D. Elata (1990), ‘An efficient ܮଶ Galerkin finite element method for multidimensional non-linear hyperbolic systems’, Inter. J. Numer. Methods in Engineering, Vol. 29, pp 1229-

1245. 

F. Bassi and S. Rebay (1997a), ‘A high-order accurate discontinuous finite element method for the 

numerical solution of the compressible Navier-Stokes equations’, J. Comput. Phys., 131 (1997), pp. 267-

279. 

F. Bassi and S. Rebay (1997b), ‘High-order accurate discontinuous finite element solution for the 2D 

Euler equations’, J. Comput. Phys., 138, pp. 251-285. 
NASA/TM—2013-218077 34 

J. S. Cagnone and S. K. Nadarajah (2012), ‘A stable interface element scheme for the p-adaptive 

lifting collocation penalty formulation’, J. Comput. Phys., Vol. 231, pp. 1615–1634. 

J.-S. Cagnone, B. C. Vermeire, and S. Nadarajah (2013), ‘A Polynomial Adaptive LCP Scheme for 

Viscous Compressible Flows’, AIAA-2013-1000. 

P. Castonguay, P. E. Vincent, A. Jameson (2012), ‘A New Class of High-Order Energy Stable Flux 

Reconstruction Schemes for Triangular Elements’, Journal of Scientific Computing, Volume 51, Number 

1, Pages 224-256. 

Castonguay, D. M. Williams, P. E. Vincent, A. Jameson (2013), ‘Energy Stable Flux Reconstruction 

Schemes for Advection-Diffusion Problems’, Computer Methods in Applied Mechanics and Engineering. 

Under Review. 

S.-C. Chang (1995), ‘The method of space-time conservation element and solution element. A new 

approach for solving the Navier-stokes and Euler equations’, J. Comp. Phys., 119, pp. 295-324. 

B. Cockburn, G. Karniadakis, and C.-W. Shu (2000), ‘The development of discontinuous Galerkin 

methods’, Springer, pp 3–50. 

B. Cockburn and C.-W. Shu (2005), ‘Foreword for the special issue on discontinuous Galerkin 

method’, Journal of Scientific Computing, 22-23, pp. 1–3. 

B. Cockburn and C.-W. Shu (2009), ‘Foreword for the special issue on discontinuous Galerkin 

method’, Journal of Scientific Computing, 40, pp. 1–3. 

M. Dumbser and C-D. Munz (2005), ‘ADER discontinuous Galerkin schemes for aeroacoustics’, C.R. 

Mecanique 333, pages 683-687. 

T. Eymann (2012), ‘Active Flux Schemes’, Ph.D. Dissertation, University of Michigan. 

H. Gao and Z.J. Wang (2009), ‘A High-Order Lifting Collocation Penalty Formulation for the NavierStokes Equations on 2-D Mixed Grids’, AIAA-2009-3784. 

G. Gassner, F. Lorcher, and C.D. Munz (2008), ‘A discontinuous Galkerin scheme based on a spacetime expansion’. II. Viscous flow equations in multi dimensions. Journal of Scientific Computing, 

34:260-286. 

T. Haga, H. Gao, Z.J. Wang (2010), ‘A High-Order Unifying Discontinuous Formulation for 3-D 

Mixed Grids’, AIAA Paper 2010-540. 

T. Haga, H. Gao and Z. J. Wang (2011), ‘A High-Order Unifying Discontinuous Formulation for the 

Navier-Stokes Equations on 3D Mixed Grids,’ Math. Model. Nat. Phenom., Vol. 6 (3), 28-56. 

E. Hairer, S. P. Norsett, and G. Wanner (1991), ‘Solving ordinary differential equations I’, 2nd 

edition, Springer-Verlag. 

E. Hairer and G. Wanner (1993), ‘Solving ordinary differential equations II’, Springer-Verlag. 

Harten, B. Engquist, S. Osher, and S. R. Chakravarthy (1987), ‘Uniformly high-order accurate 

essentially nonoscillatory schemes, III’, J. Comput. Phys., 71, pp. 231-303. 

J.S. Hesthaven and Tim Warburton (2008), ‘Nodal Discontinuous Galerkin Methods’ (Springer). 

F. B. Hildebrand (1987), ‘Introduction to Numerical Analysis’, Second Edition, Dover Books on 

Advanced Mathematics. 

H. T. Huynh (1995), ‘Accurate upwind methods for the Euler equations’, SIAM J. Numer. Anal. 32, 

1565. 
NASA/TM—2013-218077 35 

H.T. Huynh (2003), ‘Analysis and improvement of upwind and centered schemes on quadrilateral and 

triangular meshes’, AIAA Paper 2003-3541. 

H.T. Huynh (2006), ‘An upwind moment scheme for conservation laws’, Proceedings of the Third 

International Conference on Computational Fluid Dynamics, ICCFD3, Toronto, Springer, pages 761-766. 

H.T. Huynh (2007), ‘A flux reconstruction approach to high-order schemes including discontinuous 

Galerkin methods’, AIAA Paper 2007-4079. 

H.T. Huynh (2009a), ‘A Reconstruction Approach to High-Order Schemes Including Discontinuous 

Galerkin for Diffusion’, AIAA Paper 2009-403. 

H.T. Huynh (2009b), ‘Collocation and Galerkin Time-Stepping Methods’, AIAA Paper 2009-4323. 

A. Jameson (2010), ‘A proof of the stability of the spectral difference method for all orders of 

accuracy’, J. Sci. Comput. 45(1–3), 348–358. 

D.A. Kopriva and J.H. Kolias (1996), ‘A conservative staggered-grid Chebyshev multidomain method 

for compressible flows’, J. Comput. Phys. 125, 244. 

J.D. Lambert (1991), ‘Numerical methods for ordinary differential systems’, John Wiley and Sons. 

P. LaSaint and P.A. Raviart (1974), ‘On the finite element method for solving the neutron transport 

equation’, in C. de Boor, editor, ‘Mathematical aspects of finite elements in partial differential equations’, 

Academic Press, pp. 89-145. 

C. Liang, K. Miyaji and B. Zhang (2013), ‘An Efficient Correction Procedure via Reconstruction for 

Simulation of Viscous Flow on Moving and Deforming Domain’, AIAA Paper 2013-0998. 

Y. Liu, M. Vinokur, and Z.J. Wang (2006), ‘Discontinuous Spectral Difference Method for 

Conservation Laws on Unstructured Grids’, J. Comput. Phys., 216, 780–801. 

M. Lo (2011), ‘A Space-Time Discontinuous Galerkin Method for Navier-Stokes with Recovery’ 

Ph.D. Dissertation, University of Michigan. 

F. Lorcher, G. Gassner, and C.D. Munz (2007), ‘A discontinuous Galerkin scheme based on a spacetime expansion. I. Inviscid compressible flow in one space dimension’, Journal of Scientific Computing, 

32:175-199. 

R. B. Lowrie, P. L. Roe, B. van Leer (1995), ‘A space-time discontinuous Galerkin method for the 

time accurate numerical solution of hyperbolic conservation laws’, AIAA 95-1658. 

Y. Lu, X. Yuan, and W.N. Dawes (2012), ‘Investigation of 3D internal flow using new fluxreconstruction high order method’, Proceedings of ASME Turbo Expo 2012, GT2012-69270. 

Y. Lu, K. Liu, and W.N. Dawes (2013), ‘Large Eddy Simulations for 3D turbine blades using a high 

order flux reconstruction method’, ASME Paper GT2013-94707 

L. Pesch and J.J.W. van der Vegt (2008), ‘A discontinuous Galerkin finite element discretization of 

the Euler equations for compressible and incompressible fluids’, Journal of Computational Physics, 227, 

pp 5426-5446 

W.H. Reed and T.R. Hill (1973), ‘Triangular mesh methods for the neutron transport equation’, Los 

Alamos Scientific Laboratory Report, LA-UR-73-479. 

P. L. Roe (1986), ‘Characteristic-based schemes’, Ann. Rev. Fluid Mech., 18, pp 337--365. 

C-W. Shu (2012), ‘Discontinuous Galerkin method for time dependent problems: survey and recent 

developments’. Preprint. 
NASA/TM—2013-218077 36 

C-W. Shu and S. Osher (1989), ‘Efficient implementation of essentially non-oscillatory shockcapturing schemes’, J. Comp. Phys. 83, pp 32-78. 

A. Suresh and H. T. Huynh (1997), ‘Accurate Monotonicity-Preserving Schemes with Runge–Kutta 

Time Stepping’, J. Comp. Phys. 136, pp 83-99. 

Y. Suzuki and B. van Leer (2007), ‘An analysis of the upwind moment scheme and its extension to 

systems of nonlinear hyperbolic-relaxation equations’, AIAA Paper 2007-4468. 

G. D. van Albada, B. van Leer, and W. W. Roberts., Jr. (1982), ‘A comparative study of 

computational methods in cosmic gas dynamics’, Astronom. and Astrophys., 108, pp. 76-84. 

J.J.W. van der Vegt and H. van der Ven (2002), ‘Space-time discontinuous Galerkin finite element 

method with dynamic grid motion for inviscid compressible flows’, Part I. General formulation’, J. 

Comput. Phys., 182 546-585. 

B. Van Leer (1974), ‘Towards the ultimate conservative difference scheme: II. Monotonicity and 

conservation combined in a second-order scheme’, J. Comput. Phys., 14, pp. 361-.370. 

B. Van Leer (1977), ‘Towards the ultimate conservative difference scheme: IV. A new approach to 

nzmerical convection’, J. Comput. Phys., 23, pp. 276-298. 

B. Van Leer (1979), ‘Towards the ultimate conservative difference scheme: V. A second-order sequel 

to Godunov’s method’, J. Comput. Phys., 32, pp. 101-136. 

B. van Leer and S. Nomura (2005), ‘Discontinuous Galerkin for diffusion’, 17th AIAA-CFD 

Conference, , AIAA paper 2005-5108. 

B. C. Vermeire, J.-S. Cagnone, and S. Nadarajah (2013), ‘ILES Using the Correction Procedure via 

Reconstruction Scheme’, AIAA-2013-1001. 

P. E. Vincent, P. Castonguay, A. Jameson (2011), ‘A New Class of High-Order Energy Stable Flux 

Reconstruction Schemes’, J. Sci. Comput., Volume 47, Number 1, Pages 50-72. 

Z.J. Wang and H. Gao (2009). ‘A Unifying Lifting Collocation Penalty formulation including the 

discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids’, J. 

Comput. Phys., 228, No. 2, pp. 8161–8186. 

Z.J. Wang, H. Gao and T. Haga (2011). ‘A Unifying Discontinuous Formulation for Hybrid Meshes’, 

in ‘Adaptive High-Order Methods in Computational Fluid Dynamics’, World Scientific, edited by Z.J. 

Wang. 

Z.J. Wang,Y. Liu, G. May, A. Jameson (2007), ‘Spectral difference method for unstructured grids II: 

extension to the Euler equations’, J. Sci. Comput. 32, 45 

Z.J. Wang, L. Zhang and Y. Liu (2004), ‘Spectral (finite) volume method for conservation laws on 

unstructured grids IV: extension to two-dimensional Euler equations’, J. Comput. Phys., 194, No. 2, pp. 

716–741. 

T. Warburton and T. Hagstrom (2008), ‘Taming the CFL number for discontinuous Galerkin methods 

on structured meshes’, SIAM J. Numer. Anal., Vol 46, pp 3151-3180. 

D. Williams, P. Castonguay, P. E. Vincent, A. Jameson (2011), ‘An Extension of Energy Stable Flux 

Reconstruction to Unsteady, Non-linear, Viscous Problems on Mixed Grids’, AIAA Paper 2011-3405. 

H. C. Yee, R. F. Warming, and A. Harten (1983), ‘Implicit Total Variation Diminishing (TVD) 

Schemes for Steady-State Calculations’, AIAA 83-1902. 

REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 

control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY)

01-09-2013 

2. REPORT TYPE

Technical Memorandum 

3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE

High-Order Space-Time Methods for Conservation Laws 

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Huynh, H., T. 

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

WBS 794072.02.03.02.03 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

John H. Glenn Research Center at Lewis Field 

Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION

 REPORT NUMBER

E-18759 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration 

Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S

 ACRONYM(S)

NASA 

11. SPONSORING/MONITORING

 REPORT NUMBER

NASA/TM-2013-218077 

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited 

Subject Categories: 01 and 64 

Available electronically at http://www.sti.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802

13. SUPPLEMENTARY NOTES

14. ABSTRACT

Current high-order methods such as discontinuous Galerkin and/or flux reconstruction can provide effective discretization for the spatial 

derivatives. Together with a time discretization, such methods result in either too small a time step size in the case of an explicit scheme or a 

very large system in the case of an implicit one. To tackle these problems, two new high-order space-time schemes for conservation laws are 

introduced: the first is explicit and the second, implicit. The explicit method here, also called the moment scheme, achieves a CourantFriedrichs-Lewy (CFL) condition of 1 for the case of one-spatial dimension regardless of the degree of the polynomial approximation. (For 

standard explicit methods, if the spatial approximation is of degree p, then the time step sizes are typically proportional to 1/p2). Fourier 

analyses for the one and two-dimensional cases are carried out. The property of super accuracy (or super convergence) is discussed. The 

implicit method is a simplified but optimal version of the discontinuous Galerkin scheme applied to time. It reduces to a collocation implicit 

Runge-Kutta (RK) method for ordinary differential equations (ODE) called Radau IIA. The explicit and implicit schemes are closely related 

since they employ the same intermediate time levels, and the former can serve as a key building block in an iterative procedure for the latter. 

A limiting technique for the piecewise linear scheme is also discussed. The technique can suppress oscillations near a discontinuity while 

preserving accuracy near extrema. Preliminary numerical results are shown.

15. SUBJECT TERMS

Computational fluid dynamics (CFD); Discontinuous Galerkin; High-order methods; Conservation laws; Flux reconstruction; 

Explicit/implicit time stepping methods; Limiting technique 

16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF

 ABSTRACT

UU 

18. NUMBER

 OF 

 PAGES

44 

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email:help@sti.nasa.gov) 

a. REPORT

U 

b. ABSTRACT

U 

c. THIS 

PAGE

U 

19b. TELEPHONE NUMBER (include area code)

443-757-5802 

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18