Compact Higher-order Gas-kinetic Schemes with Spectral-like Resolution for Compressible Flow Simulations

In this paper, a class of compact higher-order gas-kinetic schemes (GKS) with spectral resolution will be presented. Based on the high-order gas evolution model in GKS, both the interface flux function and conservative flow variables can be evaluated explicitly from the time-accurate gas distribution function. As a result, inside each control volume both the cell-averaged flow variables and their cell-averaged gradients can be updated within each time step. The flow variable update and slope update are coming from the same physical solution at the cell interface. Different from many other approaches, such as HWENO, there are no additional governing equations in GKS for the slopes or equivalent degrees of freedom independently inside each cell. For nonlinear gas dynamic evolution, the above compact linear reconstruction from the symmetric stencil can be divided into sub-stencils and apply a biased nonlinear WENO-Z reconstruction. In GKS, the time evolution solution of the gas distribution function at a cell interface covers a physical process from an initial non-equilibrium state to a final equilibrium one. The GKS evolution models unifies the nonlinear and linear reconstructions in gas evolution process for the determination of a time-dependent gas distribution function. This dynamically adaptive model helps to solve a long lasting problem in the development of high-order schemes about the choices of the linear and nonlinear reconstructions. Compared with discontinuous Galerkin (DG) scheme, the current compact GKS uses the same local and compact stencil, achieves the 6th-order and 8th-order accuracy, uses a much larger time step with CFL number $\geq 0.3$, and gets accurate solutions in both shock and smooth regions without introducing trouble cell and additional limiting process. At the same time, the current scheme solves the Navier-Stokes equations.


Introduction
High-order methods with spectral resolutions and shock capturing capability are needed in many engineering applications, such as turbulent flows, aeroacoustics, and various complex flows with shock and boundary layer interactions. With the great potentials of the highorder methods in solution accuracy and computational efficiency, extensive effort has been paid to the development of high-order schemes in the past decades. However, the spectral accuracy and shock capturing seem to be contradictory in the current CFD algorithms, where the linear schemes for the smooth flow and nonlinear schemes for discontinuous flow play different roles under different flow conditions [1]. It seems hard to possess both properties in a single high-order method.
For smooth flows, the compact schemes [26,34] and discontinuous Galerkin methods (DG) [43,10,11] are very attractive. The compact finite difference scheme constructs implicit relation between derivatives and flow variables on compact stencils. Dispersion and dissipation properties of the scheme have been fully analyzed [26], and the scheme influences greatly on the development of linear schemes. The reason for designing compact scheme is its flexibility in complex geometries and effective parallelization. The DG has a secondorder scheme stencil, but achieves high accuracy by using high-order piecewise polynomials within elements and evolving the multiple degrees of freedom (DOFs). Each element in DG only interacts with its neighboring elements and the scheme becomes very efficient for parallel computation [46]. But, both compact and DG methods seem very successful for the smooth flow, and problems are still remained for flow simulation with shock discontinuity. Theoretically, compact scheme is physically sounded because under the CFL condition, the domain of dependence is indeed only related to the neighboring cells. All non-compact schemes must be associated with dynamical weakness because they seek for help from the cells far away which may not have any physical connection in terms of the wave propagating speed.
Nonlinear schemes have been designed for the flows with discontinuous shocks. The successful nonlinear schemes include total variation diminishing (TVD) [16], essentially nonoscillation (ENO) [17], and weighted ENO (WENO) [29,23]. In past twenty years, the WENO-based methods have received the most attention among nonlinear schemes. The central ingredients in WENO scheme is to construct several low order polynomials and to design smoothness indicators to adaptively assemble them to get a higher one. Most current effort is about the selection of optimal stencil and the design of weighting functions. WENO scheme can achieve very high-order accuracy in the smooth region and maintain nonoscillatory property across shock waves [45]. But, their numerical dissipation is still very higher than the linear schemes [19,48], and the modified WENO schemes, such as WENO-M and WENO-Z [18,7], have been proposed recently. The hybrid linear and nonlinear schemes have been investigated as well [44,56,8].
In the existing DG and WENO methods, the first-order Riemann flux plays a key role for the flow evolution of compressible flow. Recently, beyond the 1st-order Riemann solver, many schemes based on the time-dependent flux functions have been developed, such as the Lax-Wendroff method [25], the generalized Riemann problem (GRP) [3,5,4], and ADER framework [50,13]. An outstanding method is the two-stage fourth-order scheme for the Euler equations [27], where both the flux and its time derivative are used in the construction of higher-order scheme. The two-stage fourth-order algorithms have been developed under the multi-stage multi-derivative (MSMD) framework [15,42,9]. Similar to the GRP method, the gas-kinetic scheme (GKS) is also based on a time accurate flux function at a cell interface [54,53,55]. The flux evaluation in GKS is based on a time evolution solution of kinetic model equation, which provides a physical process for the gas evolution from the initial nonequilibrium state to an equilibrium one. The comparison between GRP and GKS has been presented in [28]. Under the same MSMD framework, based on the WENO reconstruction a two-stage fourth-order (S2O4) GKS has been successfully developed for the Euler and Navier-Stokes equations [39,35]. The robustness of the S2O4 GKS is as good as the secondorder shock capturing scheme [36]. Due to the use of only two reconstructions in the 4thorder scheme, the GKS becomes very efficient in comparison with schemes based on the Runge-Kutta time stepping technique. By combining the second-order or third-order GKS fluxes and MSMD technique again, a family of higher-order gas-kinetic methods has been constructed [22,21].
The time dependent gas-distribution function in GKS at the cell interface provides not only the flux and its time derivative, but also time accurate flow variables at the cell in-terface. The design of compact GKS based on the cell averaged and cell interface values has been conducted [52,37,38,21]. In this paper, we are going to develop 6th-order and 8th-order compact GKS in rectangular mesh for compressible flow simulations. The construction of the compact high-order schemes is based on the following ingredients. Firstly, the high-order gas evolution model at cell interface is used for the flux and flow variable evaluation. Consequently, the cell averaged flow variables and cell averaged slopes can be updated. Secondly, based on the time-dependent flux and its time derivative, the two-stage fourth-order MSMD technique is used for the 4th-order temporal accuracy. Thirdly, based on the flow variables and their slopes within each cell, compact nonlinear and linear 6th-order and 8th-order spatial reconstructions will be designed for the determination of the initial non-equilibrium state and the evolved equilibrium one. Based on the GKS evolution model, in the smooth flow region the linear reconstruction will be achieved quickly in the determination of interface gas distribution function, and the spectral-like resolution can be obtained. At the same time, in the discontinuous region the nonlinear reconstruction will persist in the determination of gas evolution due to long existence of the non-equilibrium state, and the scheme becomes robust in capturing shock and discontinuous solutions. Since the gas evolution model in GKS is based on the physical process from the initial non-equilibrium state to the final equilibrium one, the nonlinear and linear reconstructions are naturally unified in a time evolution process with automatically identified non-equilibrium and equilibrium regions. As a result, both discontinuous shocks and smooth aeroacoustic waves can be accurately captured by the current compact 6th-order and 8th-order GKS. The compact GKS provides the Navier-Stokes solutions directly. Even for the 6th-order and 8th-order accuracy, the current compact scheme can use a time step determined by the CFL condition with a large CFL number (≥ 0.5) in all test cases in this paper. The current 6th-order and 8th-order compact GKS take advantages of both nonlinear WENO-Z reconstruction for the introduction of numerical dissipation in discontinuous region and linear reconstruction with much reduced dissipation in resolving small amplitude waves. The current compact GKS are ideal schemes for the study of multiscale and complex flow interactions, such as the flow transition and turbulence.
The idea of using flow variables and their gradients has been investigated in finite volume or finite difference methods. The finite volume Hermite WENO (HWENO) schemes update both flow variables and their first derivative [40,41,59]. The HWENO approach is also used in the hybrid schemes [2], and a monotonicity preserving strategy for detecting troubled zones is proposed. The HWENO schemes on unstructured grids have been extensively investigated [30,31,32]. The HWENO approach needs additional limiters as both DG and reconstructed DG methods. The major advantage of these HWENO approaches is the compactness of stencils in the reconstructions. The current compact GKS is different from the HWENO approach is mainly on the update of gradients. In GKS, there is no explicit evolution equation for the gradients and the gradients are obtained rigorously from the evolution solutions at cell interfaces. Therefore, the gradients in GKS are not introduced as additional degrees of freedom with their own evolution equations. The cell-averaged gradients in GKS like cell averaged conservative variables, which are obtained from the same gas distribution functions at cell interfaces. The robustness and improved resolution in GKS is due to the physically reliable gradients from the evolution solution of the original governing equation. In other words, the updated gradients in GKS are not coming from the "equations" by taking additional derivatives to the original ones, such as the HWENO and multi-layer compact (MLC) schemes [1]. In the HWENO and MCL, the gradients have their own evolution equations which are only reliable for the smooth flow, where new governing equations can be continuously created by taking spatial derivatives to the original equation.
This paper is organized as follows. The GKS and MSMD method will be introduced in Section 2. In Section 3, the compact linear and nonlinear 6th-order and 8th-order reconstructions will be presented. The dissipation and dispersion will be analyzed as well. In Section 4, the compact GKS will be tested in wide range of flow problems from the strong shock interaction to the linear acoustic wave propagation. The last section is the conclusion.

Gas-kinetic scheme and multi-stage multi-derivative method
In the past years, the gas-kinetic scheme (GKS) has been developed systematically [54,53]. The gas evolution model used in this paper for the flux and interface flow variable evaluation is almost identical to that of the fourth-order compact scheme [21]. Therefore, only a brief introduction about GKS will be presented in this section.

Gas-kinetic scheme
The gas kinetic evolution model in GKS is based on the BGK equation [6], where f is the gas distribution function, g is the corresponding equilibrium state that f approaches, and τ is defined as the collision time. In two-dimensional case, the equilibrium state is the Maxwellian distribution where λ = m/2kT , and m, k, T represent the molecular mass, the Boltzmann constant, and temperature, K is the number of internal degrees of freedom, i.e. K = (4 − 2γ)/(γ − 1) for two-dimensional flow, and γ is the specific heat ratio. ξ is the internal variable with ξ 2 = ξ 2 1 + ξ 2 2 + ... + ξ 2 K . The collision term satisfies the following compatibility condition where ψ = (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) T = (1, u, v, The macroscopic mass ρ, momentum (ρU, ρV ), and energy ρE at a cell interface can be 5 evaluated from the gas distribution function, and the flux in the x direction can be evaluated as well, Therefore, the central point in GKS is to evaluate the time-dependent gas distribution function f at a cell interface. By direct modeling on the mesh size scale [54], the conservations of mass, momentum and energy in a control volume get to written as where W ij is the cell averaged conservative variables, F i±1/2,j (t) and G i,j±1/2 (t) are the time dependent fluxes at cell interfaces in x and y directions. Like the traditional high-order finite volume scheme, the Gaussian quadrature points are used in order to get a high accuracy. In this paper, three-points Gaussian quadrature will be used on each cell interface in the 2D case.
To model the interface gas distribution function, the integral solution of BGK equation (1) is used where (x i+1/2 , y j ) = (0, 0) is the location of cell interface, and x i+1/2 = x + u(t − t ) and y j = y + v(t − t ) are the trajectory of particles. Here f 0 is the initial gas distribution function and g is the equilibrium state in space and time. The integral solution basically states a physical process from the particle free transport in f 0 in the kinetic scale to the hydrodynamic flow evolution in the integral term of g. The contributions from f 0 and g in the determination of f at the cell interface depend on the ratio of time step to the local particle collision time, i.e., exp(−∆t/τ ), and the modeling of f 0 and g in space and time is needed to evaluate the solution f from Eq. (5). For the continuum flow simulation, such as the NS solutions, the determination of f 0 and g depend only on the macroscopic flow variables and their initial reconstructions. In this paper, the WENO-Z method will be used as a nonlinear reconstruction in the determination of f 0 , and the linear reconstruction is adopted in the determination of g. Therefore, the above integral solution not only incorporates a physical evolution process from non-equilibrium to an equilibrium state, but alo from a nonlinear reconstruction to a linear one. This fact is critically important for the current scheme to capture both nonlinear shock and linear acoustic wave accurately in a single computation. Based on the integral solution, a simplified third-order gas distribution function can be obtained [58], where the terms related to g 0 are from the integral of the equilibrium state and the terms related to g l and g r are from the initial term f 0 in the Eq.(5). All other terms, such as A, a tt , a, ..., are coming from the derivatives of a Maxwellian distribution. All these coefficient in the above equation can be determined from the initially reconstructed macroscopic flow variables.

Multi-stage multi-derivative time stepping method
With the above time accurate gas distribution function at a cell interface, in order to improve the temporal accuracy of the scheme, the multi-stages multi-derivative (MSMD) technique will be used for the evaluation of the flux function and the interface values [21]. The implementation of MSMD needs the flux function and its time derivative. For conservation laws, the semi-discrete finite volume scheme is written as where L ij (W ) is the numerical operator for spatial derivative of flux, F and G are obtained fluxes at the Gaussian quadrature points. A fourth-order temporal accurate solution for W (t) at t = t n + ∆t can be obtained by where L and ∂ ∂t L are related to the fluxes and the time derivatives of the fluxes evaluated from the time-dependent gas distribution function f (t) at cell interfaces. And the middle state W * is obtained at time t * = t n + ∆t/2. Again, with the time accurate gas distribution function f (t), along the same line with MSMD technique the gas distribution function f at a cell interface at t n+1 becomes where f * is for the middle state at time t * = t n + ∆t/2, Therefore, based on the cell interface f n+1 j+1/2 , the flow variables W n+1 j+1/2 can be explicitly obtained, i.e., W n+1 j+1/2 = ψf n+1 j+1/2 dΞ, from which the slope inside each cell, such as in the 1D case, can be updated as More detailed formulation for the above GKS part can be found in [21].

Compact WENO reconstruction
In the classical high-order finite volume method [17,29], the pointwise value is reconstructed based on cell averages, and the Riemann solver is used for the flux evaluation at the interface, and the Runge-Kutta time stepping method is used for the time accuracy. However, for the gas kinetic scheme, as presented in the last section the cell-averaged conservative flow variables and cell-averaged slopes inside each cell can be updated. Based on these updated data, a 6th-order and a 8th-order compact linear and nonlinear reconstruction will be designed consistently in this section, which can be subsequently used in the determination of equilibrium state and initial non-equilibrium one in the gas kinetic scheme.

Compact linear reconstruction
In the section 2, the value W i+1/2 is provided by taking moments of the gas distribution function at the cell interface. The averaged slope in the cell I i can be obtained by the Gaussian formula where W i is the cell averaged slope in the cell I i . By taking advantage of the cell averaged variables and their slopes, the compact stencil Fig.1 will be used for the reconstruction of the value at the cell interface x i+1/2 . On each cell I i , the cell averages and their slopes are denoted as W i and W i . Based on the large stencil S i+1/2 , a series of high-order polynomials can be constructed, where n is the order of the polynomial and n = 5, 7 correspond to the 6th-order and 8th-order schemes, respectively. The polynomial P 7 (x) can be uniquely determined with the condition In order to determine the P 5 (x), the above condition has to be modified. For the cells I i and I i+1 , the following equations are strictly enforced and for the cells I i−1 and I i+2 the following conditions are satisfied in the sense of the least square Thus the linear reconstruction P n (x i+1/2 ) at the cell interface x = x i+1/2 can be written as where P n (x i+1/2 ) has the truncation error of O(∆x n+1 ). The first-order and second-order derivatives of the linear reconstruction at the cell interface can be obtained from P n (x), In order to eliminate the spurious oscillation and improve the stability, the above reconstruction is based on the characteristic variables. The dimension-by-dimension strategy is applied for the reconstruction in the two-dimensional case. The above compact reconstruction is mainly used to determine the equilibrium state in the normal direction in the gas kinetic scheme. Due to multi-dimensionality in the GKS evolution model, in the tangential direction a linear 5th-order reconstruction is used in both 6th-order and 8th-order schemes. The above linear reconstruction is used to determine the equilibrium state in GKS. In the reconstruction of the initial non-equilibrium state, the WENO-based compact nonlinear reconstruction will be designed in later subsection.

Resolution analysis
In this section, the resolution analysis of the above compact spatial reconstruction will be presented. The linear advection equation is used for the resolution analysis, The integral form on the cell I i becomes where W i is the cell averaged W (x) in the cell I i .
Remark 1. In order to make Fourier expansion for the cell averages W i , a continuous function W (x) is introduced [17] as follows.
Thus the cell average W i becomes the value of the continuous function W (x) at cell center In order to analyse the resolution of spatial discretization in the current compact schemes, we can compare the Fourier coefficients of the derivative given by the spatial discretization with the Fourier coefficients of exact derivative. Using the reconstructed value at the cell interface, the spatial discretization of (W i ) x in Eq. (19) becomes For the purpose of Fourier analysis, the function W (x) is assumed to be periodic over the domain[0, L], and W (x) can be decomposed into its Fourier coefficients It is convenient to introduce new variables ω = 2πk∆x/L = 2πk/N and s = x/∆x [26]. Then the Fourier expansion can be written as Fourier expansion, we get the modified wavenumber In addition to the modified wavenumber, the error in the phase speed can be alternatively used to evaluate the dispersive error. For the linear advection equation, the phase speed for a wavenumber ω is given by the current spatial discretization (exact time advancement is assumed) as  Figure 2: Plots of modified wavenumber and phase speed vs wavenumber for different schemes. cs-6th is the compact 6th-order tridiagonal scheme studied by Lele [26], and the scheme has the best resolution in the series of 6th-order scheme; cs-8th is the compact 8th-order pentadiagonal scheme [26], and the cs-8th uses the same formula in spatial discretization as the compact 8th-order GKS.
The comparison of the dispersion characteristics of the current compact GKS with the compact schemes by Lele [26] is shown in Fig.2. The traditional compact schemes use only the node values [26] or cell averages [12] as independent values in spatial discretization, and the curves of modified wavenumber and phase speed in the range ω ∈ (π/2, π) deviate greatly from the curves corresponding to the analytical solution. Therefore, these schemes have poor resolution for the wave with ω ∈ (π/2, π), because the wavenumber range (π > ω > π/2) corresponds to the wavelength (2∆x, 4∆x), and such wave cannot be determined by the cell averages (or point values) alone with less than four points in a wavelength. However, the compact GKS can retain a good resolution because of the use of slopes in each cell, and these slopes are evaluated from the gas evolution solutions at the cell interfaces which are independent from the cell averaged values. In the traditional compact schemes by Lele, the node values and derivatives are not fully independent. The coupling increases the error since the node values cannot resolve high wavenumber solution at all. As an example, we test the schemes for the initial condition (ρ, U, p) = (1+0.2 sin(2πx), 1, 1) in a computational domain [0, 2]. The results of linear density wave at t = 1 with meshes ∆x = 1/2 and ∆x = 1/4 are shown in Fig. 3. These results demonstrate that the current compact schemes have favorable resolution in comparison with the non-compact 7th-order linear upwind scheme.

WENO reconstruction
In gas dynamics, shock and contact discontinuities can appear. To deal with discontinuities, WENO reconstruction [29,23] can be used. In section 3.1, the linear reconstruction gives a unique reconstructed value and its slope at the cell interface. However, in order to capture the possible discontinuity at a cell interface, the values at the left and right sides of the cell interface have to be valuated separately. Therefore, for the nonlinear reconstruction the sub-stencils have to be defined first. For simplicity, the WENO reconstruction procedure is given in detail for the construction of the left side value of the cell interface x = x i+1/2 . The procedure for the right side value can be obtained similarly according to the symmetric property, which will be omitted here. The left side value by WENO reconstruction is given by the candidate polynomials as follows where W n i+1/2 is the nth − order WENO reconstruction for the left value of the cell interface x i+1/2 , l n is the number of candidate polynomials, δ n k is the WENO weight, and w n k,i+1/2 is the point value of the candidate polynomial w n k (x) at x i+1/2 . In the later presentation, the sets of candidate polynomials w n k (x) are the same for different nth−order scheme, such as n = 6 and 8, so the superscript is omitted for simplicity. In the current compact WENO reconstruction, three principles are considered. Firstly, if a discontinuity is located at one of the cell interfaces, theoretically the point-wise value at the interface cannot be reliable to determine the cell averaged slope on the left and right side cells. Secondly, because the smooth sub-stencils can play a dominant role in dealing with discontinuities appearing in the large stencil, the averaged slope used in the sub-stencil should be possibly away from the interface for the data reconstruction in order to make the sub-stencil sufficiently smooth. Thirdly, the order of the candidate polynomial can be higher without increasing the spatial size of sub-stencils in the compact schemes, and such higherorder candidate polynomials can achieve better resolution to solve discontinuities without spurious oscillation. Based on the above principles, the following sub-stencils are designed, where the 6th-order and 8th-order reconstructions are based on the same sub-stencils, w 2 (x) is a cubic polynomial, w 5 (x) and w 6 (x) are fourth-order polynomials, and others are quadratic polynomials. According to the similar reconstruction conditions in Eq.(11), w k,i+1/2 can be obtained as In the smooth region, the convex combination with δ n k = d n k recovers the reconstruction in Eq. (14) and Eq. (15), which can be the condition to get the linear weights [23]. The linear weight d 6 k of the 6th-order reconstruction can be obtained as The WENO-Z nonlinear weights are used in the current compact schemes and they are defined as [7] δ n k = where Z n ref is the local higher-order reference value, which is related to the accuracy of the scheme and is given in detail later. Here β k is the smooth indicator and defined as [23] where q k is the order of w k (x). In smooth region, the first two candidate polynomial w 0 (x) and w 1 (x) can be combined into a cubic polynomial which is symmetric counterpart of w 2 (x). Then, current sub-stencils can become basically symmetric for interface x i+1/2 . In order to maintain the symmetry for the nonlinear schemes, the smooth indicator β 1 of w 1 (x) corresponding to S 1 is replaced by the indicator of the cubic polynomial w 1 (x) on S 1 = {W i−1 , W i , W i+1 , W i−1 }. Even without showing in this paper, some tests demonstrate that the current choice β 1 can present a better resolution in the numerical results with excellent robustness, The detailed formulae for all β k , k = 0, ..., l n are given in Appendix.

Accuracy of the compact reconstruction
In this subsection, the accuracy of the compact reconstruction is analysed. The nonlinear reconstruction in Eq.(26) can be rewritten as The reconstruction can be split into two terms [18], i.e., the linear term and nonlinear term. For the linear part W n,opt i+1/2 ≡ ln k=0 d n k w k,i+1/2 , the error can be given as where W (x i+1/2 ) is the exact solution at x = x i+1/2 . In the second term in Eq.(30), the point value w k,i+1/2 approximates W (x i+1/2 ) as follows Substituting w k,i+1/2 into Eq.(30) and taking ln k=0 δ n k = ln k=0 d n k = 1 into account, we have To achieve the formal order of accuracy for the 6th-order and 8th-order reconstructions, the following sufficient condition is proposed.
For the current nonlinear weights with WENO-Z weighting functions, with the following formulation of the local higher-order reference value Z n ref , the sufficient condition Eq.(31) can be satisfied when In order to prove that the sufficient condition can be met, according to the Taylor series of β k , k = 0, 3, 4,

Numerical tests
In this section, we are going to test the 6th-order and 8th-order compact gas kinetic schemes. The cases include linear acoustic waves, blast wave, shock-shock interactions, shock acoustic wave interaction, as well as viscous flow computations. The mesh used in this paper is rectangular one and the time step is determined by the CFL condition with a CFL number (≥ 0.3) in all test cases. In all tests, the same linear reconstruction is used for the equilibrium state and the nonlinear reconstruction for the non-equilibrium state. There is no any additional "trouble cell" detection or any limiter designed for specific test. The gas kinetic evolution model basically presents a dynamical process from the nonlinear to linear one, and the convergence rate exp(−∆t/τ ) depends on the flow condition. The collision time τ for inviscid flow at a cell interface is defined by where ε = 0.01, C = 1, and p l and p r are the pressures at the left and right sides of a cell interface. For the viscous flow, the collision time is related to the viscosity coefficient, where µ is the dynamic viscosity coefficient and p is the pressure at the cell interface. In smooth flow regions, it will reduce to τ = µ/p. The reason for including pressure jump term in the particle collision time is to add artificial dissipation in the discontinuous region to keep the non-equilibrium dynamics in the shock layer.

Accuracy tests
The one-dimensional advection of density perturbation is tested first, and the initial conditions are given as follows The periodic boundary condition is adopted, and the analytic solutions are ρ(x, t) = 1 + 0.2 sin(π(x − t)), U (x, t) = 1, p(x, t) = 1.
With the rth-order spatial reconstruction and 2-stage 4th-order temporal discretization, the leading term of the truncation error is O(∆x r + ∆t 4 ). To keep the rth-order accuracy in the test, ∆t = C∆x r/4 needs to be used for the rth-order scheme. In the computation, the uniform meshes with N mesh points are used. The L 1 , L 2 and L ∞ errors and convergence orders at t = 2 for the 6th-order and 8th-order compact GKS are presented in Table 1 to Table 4, respectively. For the 8th-order scheme, the errors reduce quickly, and the order of accuracy does not attain 8 with N = 80 due to the limited round-off error.

One-dimensional Riemann problems
To test the performance of the schemes for capturing high frequency waves, the Shu-Osher problem for density-wave shock interaction is tested [47]. The initial condition is given by (ρ, U, p) = (3.857134, 2.629369, 10.33333), The computational domain is [0, 10] and 200 uniform mesh points are used. The nonreflecting boundary condition is used at both ends. The computed density profile and local enlargement at t = 1.8 from different schemes are shown in Fig.5. Current compact 6thorder and 8th-order schemes perform well with such a coarse mesh. Due to the improvement of order of accuracy, the 8th-order scheme performs better then the 6th-order one.
As an extension of the Shu-Osher problem, the Titarev-Toro problem is tested as well  The computational domain is [−5, 5]. The non-reflecting boundary condition is imposed on left end, and the fixed wave profile is given on the right end. Both compact 6th-order and 8th-order schemes are tested for this case. The computed density profile with 1000 mesh points at t = 5, local enlargement, and the exact solution for the Titarev-Toro problem are shown in Fig.6. In order to show the importance and accuracy of the 6th-order and 8th-order compact reconstructions, the results from the same GKS, but with the non-compact standard 7th-order WENO-JS reconstruction using the cell-averaged flow variables only, are included as well. As shown in Fig.6, even with 7th-order WENO reconstruction the dissipation and dispersion errors are much larger than those from the 6th-order compact GKS. Based on this observation, it clearly indicates that the use of high-order evolution model for the update of slope is favorable in the design of high-order schemes. The compactness of the scheme is also physically sounded because the CFL condition constrains the signal propagation locally to the neighboring cells only within a time step. To design a reliable compact scheme depends on the high-order gas evolution model at the cell interface. Any scheme based on the first order Riemann solver cannot achieve such a goal to have scheme with the properties of compactness, large CFL number, robustness, and high efficiency. The third test is the Woodward-Colella blast wave problem [51], and the initial conditions

One-dimensional acoustic wave
The one-dimensional acoustic wave propagation in x-direction was proposed by Bai et al. in [1]. The test demonstrates the high order and high resolution of the compact GKS to compute acoustic wave propagating through a long-distance. The initial conditions are given as follows where = 10 −5 is the magnitude of initial perturbation, and ω = 6π is the wavenumber of initial perturbations in velocity. The specific heat ratio is γ = 1.4. The acoustic wave given above is approximately linear because of the very small , and an analytical solution [1] is given from the approximate acoustic wave equation, The computational domain is [0, 1/3]. Periodic boundary conditions on both sides are adopted.
In the computation, we first obtain the cell average of the primitive variables, i.e.ρ i = i+1/2 i−1/2 ρ(x)dx exactly, then convert them to the conservative variables. Similarly, we could obtain the cell average of the derivatives, i.e.ρ x,i = i+1/2 i−1/2 ρ x (x)dx exactly, from which the derivatives of conservative variables can be obtained by the chain rules. Theoretically, this has a 2nd-order accuracy in space. However, it is enough for the comparison since the analytical solution in Eq. (33) is also an approximate solution. Fig. 8 and Fig. 9 show the density distributions at t = 0.01 (left) and t = 1.0 (right) with 20 and 40 uniform mesh points. The computational time is greatly larger than λ ρ 0 /a ∞ = 4.8 × 10 −4 , where λ ρ 0 is the initial wavelength of density perturbation. Thus the distance of acoustic wave propagation is greatly larger than the initial wavelength of density perturbation. The results from GKS with linear and nonlinear reconstructions are in good agreement. The results with 40 mesh points agree well with the analytical solution, where the small deviation near the extreme point is caused by the approximate analytical solution [1]. This case also shows that the 6th-order GKS cannot give the accurate solution with 20 mesh points after a long time wave propagation (100 periods) due to the relative large numerical error in comparison with the 8th-order scheme.

One-dimensional shock crossing velocity perturbation
The interaction between strong shock wave and turbulence is an important flow problem in gas dynamics, where the shock can cross the local wave packets which propagate faster or slower than the mean velocity. To validate the current schemes for this kind of problems, we propose a test of shock with a Mach number M a = 10 crossing a velocity perturbation. The initial condition is given as follows (ρ, U, p) = (8, 8.25, 116.5), x ≤ 1, (1.4, 0.2 sin(2π(x − 1)), 1), 1 < x ≤ 10. The computational domain is [0,10]. In order to compare the resolution for capturing the multi-scale structure from different schemes, a relative coarse mesh with 400 points is used in the computation. The non-reflecting boundary condition is used at both ends. The computed density profile and the local enlargement from different schemes at t = 0.8 are shown in Fig.10, and the distribution of velocity and its local enlargement are shown in Fig.11. The current compact 6th-order and 8th-order schemes perform well for resolving the multi-scale flow structures, where the non-compact WENO-7JS GKS has a relative large error. The compactness seems a preferred choice for a scheme with a better wave resolving power.

Double Mach reflection problem
In this subsection, the double Mach reflection problem is tested. The test was extensively studied by Woodward and Colella [51] for the inviscid flow. The computational domain is The reflecting boundary condition is used at the wall, and the exact post-shock condition is imposed for the rest of bottom boundary. At the top boundary, the flow variables are set to describe the exact motion of the Mach 10 shock. In this case, the compact 6th-order and 8th-order compact GKS are tested. The density distributions with 960 × 240 and 1920 × 480 uniform mesh points at t = 0.2 are shown in Fig. 12 and Fig. 13, and the corresponding local density enlargements are shown in Fig. 14. The 6th-order and 8th-order schemes resolve the flow structure under the triple Mach stem very well.

Lid-driven cavity flow
Here we test the performance of the compact GKS for the capturing of viscous flow solutions. The lid-driven cavity problem is one of the most important benchmarks for validating incompressible Navier-Stokes flow solvers. The fluid is bounded by a unit square and is driven by a uniform translation of the top boundary. In this case, the flow is simulated with  Fig. 18. The benchmark data [14] for Re = 1000 and 3200 are also presented, and the simulation results match well with these benchmark data. The results from 8th-order compact GKS are almost identical to the reference solutions even for Re = 3200 case. The cavity case fully validates the higher-order accuracy of compact GKS for viscous flow.

Two-dimensional acoustic wave
The case is an interaction of a shock wave with a single vortex as studied by Inoue and Hattori [20], where the fluid is viscous. The computational domain is [−20, 8] × [− 12,12]. The initial counterclockwise vortex (the case C in [20]) is set as follows where U θ and U r are the tangential and radial velocity respectively. Mach number M v of the vortex is M v = 0.25. The Mach number of shock wave is M s = 1.2. The Reynolds number is Re = 800 defined by Re = ρ ∞ a ∞ /µ ∞ , where the subscript ∞ denotes the quantity downstream of the shock wave. The initial location of vortex is (x v , y v ) = (1, 0), and the stationary shock is at x = 0. In the computation, the supersonic inflow boundary conditions   : Two-dimensional acoustic problem: radial distribution of the sound pressure p obtained by 6th-order and 8th-order compact GKS at t = 6 with 700 × 600 mesh points. The vortex center at t = 6 is approximately located at (−3.9, 0.0). at x = 8 as well as the periodic boundary conditions at y = ±12 are imposed. The nonreflecting boundary conditions are adopted at x = −20.
In the case, the compact 6th-order and 8th-order GKS are tested on the uniform mesh with 700 × 600 cells. The sound pressure at t = 6 are shown in Fig. 19, where the sound pressure is defined as p = (p − p ∞ )/p ∞ . The pressure distribution in Fig. 19 clearly shows that the incident shock wave and two reflected shock waves are connected at the triple point, and the reflected shock waves emanate from the compression region of the incident shock wave. In addition, both the first and second sound waves have the clear quadruple nature with opposite sign, which is similar with the result in [20]. Both the 6th-order and 8th-order GKS perform well to resolve the shock wave and density wave, and the solution is sufficiently smooth to resolve small amplitude waves. The radial distribution of the sound pressure p obtained by 6th-order and 8th-order compact GKS at t = 6 with 700 × 600 mesh points are shown in Fig. 20. The radial distribution of the sound pressure is subtracted from the vortex center with the angle θ = 135 degree with the positive x direction. The vortex center at t = 6 is approximately located at (−3.9, 0.0). The results are concordant well with the results in [20], where the reference solution in [20] is obtained with 1044 × 1170 mesh points and local mesh refinement with approximate 1/10 cell size reduction at the shock region.

Conclusion
In this paper, a class of compact high-order gas-kinetic scheme with WENO reconstruction has been proposed for the compressible Euler and Navier-Stokes solutions. Based on the gas-kinetic theory, the current scheme depends solely on the time-accurate solution of a gas distribution function at a cell interface for the flux and conservative flow variable evaluation. With the update of pointwise values at cell interfaces, besides the cell-averaged conservative flow variables, their cell-averaged gradients can be updated as well. Therefore, the compact linear and nonlinear reconstructions can be obtained from the local cell averaged flow variables and their slopes. With symmetrical stencils, the 6th-order and 8th-order compact linear reconstructions are given and used in the 1-D Fourier analysis. The analysis elucidates spectral resolution of the current spatial discretization, even for very large wave number. The GKS flux makes the scheme stable even for the initial reconstruction from symmetrical stencils. In order to capture shock and other discontinuities, the nonlinear compact WENO reconstruction has been constructed as well. Compared with other compact schemes, the current compact GKS is fully local and has explicit high-order discretization, without using any formulation for the globally connected flow variable and their slopes. Equipped with both linear and nonlinear compact reconstruction, the GKS gas evolution model can make a smooth transition in dynamics from the initial nonlinear reconstruction for non-equilibrium state to the final linear reconstruction for equilibrium state. The transition rate depends on the local flow conditions. As a result, the current compact high-order scheme can naturally capture both shock discontinuities and small amplitude acoustic waves in a single computation without using any trouble cell detection and additional limiters. Due to the existence of the time-derivative of the flux function, the multi-stage and multi-derivative time stepping technique has been applied here for the achieving high-order temporal accuracy of the scheme. More specifically, the two stage fourth-order time accuracy has been achieved. Both GKS and DG methods have the same compact stencils. But, the GKS uses a much large CFL number (≥ 0.3) for the 6th-and 8th-order schemes, and it is much more robust in capturing shock discontinuity and is much more efficient in simulating NS equations than the DG method. The extensive tests in this paper clearly demonstrate that the current 6th-order and 8th-order compact schemes advance the development of high-order CFD methods to a new level of maturity. The success of the current compact schemes depends solely on the high gas evolution model at a cell interface, which keeps and traces flow dynamics as close as possible to a real physical time evolution process. Any scheme based on the first order Riemann solver, which gives a time-independent gas evolution at a cell interface from two constant states whatever the high-order initial reconstruction is, may have difficulty to possess all these good properties of a high-order scheme, such as the compactness, large CFL number, robustness, high efficiency, and the dynamic unification of linear and nonlinear reconstructions.