## 1 Introduction  
Most of the computational fluid dynamics (CFD) solvers currently used in aerospace engineering are based on schemes using finite volume (FV) methods, which are more suitable than schemes using finite difference (FD) methods for unstructured meshes. However, FV schemes usually have only second-order spatial accuracy, due to the difficulty of handling irregular local stencils, whose sizes grows rapidly as the order of accuracy increases. In high-order FV schemes, the local stencil of a cell is made up of the cell and the cell's neighbors, and the neighbors' neighbors, and so on. To ease the development of high-order schemes for unstructured meshes, it is better to use finite element (FE) methods whose local stencils are much more compact than those of FV schemes. One family of such compact high-order schemes suitable for unstructured meshes are called the discontinuous Galerkin (DG) methods, which assume continuous approximation in each cell (like classic FE methods), but allow discontinuities to exist on cell boundaries (like classic FV methods). Such an assumption gives scheme designers more freedom on choosing the basis for each cell, such as orthogonal functions which lead to modal DG methods, and Lagrange polynomials which lead to nodal DG methods [1]. The original DG method was introduced by [PERSON] and [PERSON] [2] for solving the neutron transport equation, which is a linear problem contains only first-order spatial derivatives. For nonlinear conservation laws, [PERSON] and [PERSON] [3] made the first attempt of using a DG method to discretize in space and using the explicit Euler method for time discretization. This scheme is second-order accurate in space but only first-order accurate in time. To improve the accuracy of time discretization, [PERSON] and [PERSON] [4] replace the first-order Euler method with a special second-order Runge-Kutta (RK) method [5, 6], which is explicit and total variation diminishing (TVD). This is the first successful Runge-Kutta discontinuous Galerkin (RKDG) method, which is high-order accurate both in space and time, for scalar conservation laws. This RKDGmethod was soon extended to one-dimensional system case [7], multi-dimensional scalar case [8] and multi-dimensional system case [9], which includes the Euler system that models inviscid compressible flows. Nearly the same time (late 1990s), the DG methods were also extended to the Navier-Stokes system that governs viscous flows. [PERSON] and [PERSON] [10] made the first attempt to apply the DG idea to both the unknowns and their gradients for solving the Navier-Stokes equations. This method was then generalized by [PERSON] and [PERSON] [11] to the family of local discontinuous Galerkin (LDG) methods, which extends the RKDG methods from conversation laws to convection-dominated problems. One may refer to [12], which is a comprehensive review article on this topic. In the rest of this article, we will only study the RKDG methods for conservation laws that are purely hyperbolic. In Sections 2.1 and 2.3, we will give a matrix formed formulation of the RKDG method for solving three-dimensional conservation law systems. When solving hyperbolic problems, limiting procedures, or limiters for short, are necessary for suppressing numerical oscillations that might occur near discontinuities (known as the Gibbs phenomenon [13]). For FD and FV schemes, there is the monotonic upstream-centered scheme for conservation laws (MUSCL) [14; 15; 16; 17], which could achieve third-order accuracy when used with caution. However, it is not directly applicable to FE schemes. Besides the RKDG formulations, there is another important contribution made by [8; 9], which is proposing the generalized slope limiters for multi-dimensional problems. These minmod type limiters, together with TVD RK methods, makes the high-order solutions free from non-physical oscillation, but tend to reduce the order of accuracy in smooth regions.