High-order compact gas-kinetic schemes for three-dimensional flow simulations on tetrahedral mesh

A general framework for the development of high-order compact schemes has been proposed recently. The core steps of the schemes are composed of the following. 1). Based on a kinetic model equation, from a generalized initial distribution of flow variables construct a time-accurate evolution solution of gas distribution function at a cell interface and obtain the corresponding flux function; 2). Introduce the WENO-type weighting functions into the high-order time-derivative of the cell interface flux function in the multistage multi-derivative (MSMD) time stepping scheme to cope with the possible impingement of a shock wave on a cell interface within a time step, and update the cell-averaged conservative flow variables inside each control volume; 3). Model the time evolution of the gas distribution function on both sides of a cell interface separately, take moments of the inner cell interface gas distribution function to get flow variables, and update the cell-averaged gradients of flow variables inside each control volume; 4). Based on the cell-averaged flow variables and their gradients, develop compact initial data reconstruction to get initial condition of flow distributions at the beginning of next time step. A compact gas-kinetic scheme (GKS) up to sixth-order accuracy in space and fourth-order in time has been constructed on 2D unstructured mesh. In this paper, the compact GKS up to fourth-order accuracy on three-dimensional tetrahedral mesh will be further constructed with the focus on the WENO-type initial compact data reconstruction. Nonlinear weights are designed to achieve high-order accuracy for the smooth Navier-Stokes solution and keep super robustness in 3D computation with strong shock interactions. The fourth-order compact GKS uses a large time step with a CFL number 0.6 in the simulations from subsonic to hypersonic flow. A series of test cases are used to validate the scheme. The high-order compact GKS can be used in 3D applications with complex geometry.

Sections 4 and 5 will present the linear and nonlinear compact reconstructions for the determination of piecewise high-order polynomials inside each control volume. In Section 6, the compact GKS will be tested in a wide range of cases in three-dimensional space on tetrahedral mesh. The last section is the conclusion.

Time-accurate gas-kinetic evolution model
For the construction of high-order compact GKS, the use of high-order gas evolution model beyond the 1st-order Riemann solution is necessary. In this section, we will briefly present the time-accurate evolution solution of a gas distribution at a cell interface.
The gas-kinetic evolution is based on the kinetic BGK equation [19], where u = (u, v, w) is the particle velocity, f is the gas distribution function, g is the corresponding equilibrium state that f approaches, and τ is the particle collision time. The equilibrium state g is a Maxwellian distribution, where = 1/2RT , and R and T are the gas constant and temperature, respectively. K is the number of internal DOFs, i.e. K = (5 − 3γ )/(γ − 1) for three-dimensional flow, and γ is the specific heat ratio. ξ is the internal variable with ξ 2 = ξ 2 1 + ξ 2 2 + ... + ξ 2 K . At a relatively low temperature without exciting the vibrational mode, a diatomic molecule in a three-dimensional flow has two rotational DOFs in ξ , such as K = 2 . U = (U, V , W ) is the macroscopic flow velocity which is the same velocity in the Navier-Stokes (NS) equations. Due to the conservation of mass, momentum and energy during particle collisions, f and g satisfy the compatibility condition, at any point in space and time, where ψ = (ψ 1 , ψ 2 , ψ 3 , ψ 4 , ψ 5 ) T = (1, u, v, w, 1 2 (u 2 + ξ 2 )) T , d� = dudvdwdξ 1 ...dξ K .
The macroscopic conservative flow variables W = (ρ, ρU , ρV , ρW , ρE) can be evaluated from the gas distribution function, The corresponding fluxes for mass, momentum and energy in the i-th direction are given by (1) (4) F i = u i f ψd , with (u 1 , u 2 , u 3 ) = u in the 3-D case. According to the Chapman−Enskog expansion [20][21][22], the zeroth-order truncation f = g corresponds to the invicid Euler equations, and the first-order truncation f = g − τ (u · ∇ u g + g t ) gives the NS equations.
In GKS, the evolution solution W(t) and F(t) at a cell interface are determined by the time-accurate gas distribution function f. Based on the integral solution of BGK equation and the modeling for the initial state and equilibrium state distribution in local space and time [21], the time-accurate distribution function at a cell interface becomes where x 0 is the numerical quadrature point at the cell interface, and is the particle trajectory. Here f 0 is the initial state of gas distribution function f at t = 0 . In order to explicitly obtain the solution f of Eq. (5), both f 0 and g in Eq. (5) need to be modeled. The second-order accurate solution for f is [21] 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 Eq. (5). All the coefficients in Eq. (6), such as a kl and a kr (k = 1, 2, 3) , can be determined from the initially reconstructed macroscopic flow variables at the left and right sides of the cell interface. The above time evolution solution is distinguishable from the generalized Riemann problem (GRP) solver [23] in the following aspects. (1). The above distribution function takes a physical process from the flux vector splitting transport to the NS solution at the cell interface x 0 ; (2). The flow evolution has multidimensional mechanism with the participation of ∂x i terms in the solution, while the GRP solver with multidimensional properties has been developed for the wave equations [24]; (3). In smooth flow region, the evolution solution in GKS gets back to the Lax-Wendroff type central difference method; (4). The cell interface evolution solution is obtained explicitly without iterations.

Solution updates
The discrete conservation law in a control volume j is, where W(x, t) is the conservative flow variable in a control volume j , and F(t) is the corresponding flux across the cell interface ∂� j . The above integral conservation law is valid in any flow regimes from the rarefied to the continuum one once the dynamics of F(t) can be properly modeled [25,26]. The accuracy of the updated solution depends critically on the time-dependent interface flux function F(t) , which depends on the initial condition W j (t n ) and the evolution model for the cell interface flux F(t).
As a compact scheme, with the time-accurate interface flow variable W(t) , the cellaveraged gradient of the flow variable can be updated as well by the Gauss's law, With the consideration of possible discontinuous flow distribution around the cell interface, the above flow variable W(t n+1 ) is the value at the inner side on the cell interface of the control volume. In other words, W(x 0 , t n+1 ) may have multiple values, such as W l (x 0 ) and W r (x 0 ) at both sides of a cell interface. The outstanding example is that a stationary shock is exactly standing on the cell interface. The evolution solution of F(t) and W l,r (x 0 ) will be presented next.

Update of cell-averaged flow variable
The conservation law in Eq. (7) can be written as where W j is the cell-averaged flow variable defined as The surface integral in L j (t) is discretized by a q-point Gaussian quadrature rule, where |Ŵ l | is the face area of the cell, l 0 is the number of cell faces, n l is the unit outer normal vector, and q and ω k are the number of quadrature points and weight of the Gaussian quadrature rule.
With the consideration of possible discontinuous flux function F across the cell interface, such as a moving shock with speed v s passing through the cell interface v s = [F n ]/[W n ] , the flux function F n (t) on a cell interface may be a discontinuous function of time. In order to capture such a dynamic evolution without introducing oscillation, same as nonlinear reconstruction polynomial of flow variables in space, the nonlinearly limited flux function in time has be developed as well. Using a fourth-order time-accurate flux function in which terms of quadratic and above are limited, Eq. (9) is discretized as where in the second stage from t n to t n+1 a nonlinear limiter for the higher-order time derivatives of the flux function is introduced [14]. The nonlinearly limited operator L j is given as where ω t l is a nonlinear weight for the lth interface of the cell and ω t l ∈ [0, 1] . The nonlinear weight ω t is defined as where α is a positive integer, and α takes α = 3 in this paper. τ t Z is the local high-order reference value to indicate smoothness of the large stencil in the spatial reconstruction. IS L,R s and IS L,R d are smooth indicators corresponding to a smooth and a possibly discontinuous candidate polynomial in the nonlinear compact spatial reconstruction in the cells on both sides of the cell interface. The determination of these parameters are given in Eq. (29) at the end of Section 5. The above fourth-order scheme in time can be reduced to a second-order one-step method for discretizing Eq. (9), which is equivalent to the Lax-Wendroff method.

Update of cell-averaged gradient
The RHS of Eq. (8) can be discretized by the same q-point Gaussian quadrature rule as that in Eq. (11). The update of cell-averaged gradient becomes where W n+1 (x k ) is the value at the inner surface of j , which may be different from the value at the other side of the surface. In order to obtain a high-order accurate flow variable at the quadrature point, the macroscopic flow variable is evolved by two stages To provide the flow variable at both sides of a cell interface [14], the update model for W(x, t) is given by The evolution solution W e is given by the moments of the time-accurate distribution function in Eq. (6), and W l 0 and W r 0 are obtained from Eq. (6) as well with the assumptions smooth initial condition on both sides of the cell interface separately. Eq. (18) is the modeling for the real flow physics. The solutions of flow variables on both sides of the interface may become discontinuous, since the appropriate approach in a shock-capturing scheme is to assume a continuous subcell flow distribution and contribute all possible discontinuity to the cell interface. The weighting function e −�t/τ 0 is constructed from a physical relaxation model, while in the smooth flow region under �t ≫ τ 0 , a single flow variable at the cell interface is recovered. The relaxation time τ 0 is defined as where p l and p r are the pressures at both sides of the cell interface, and ε diss is a constant coefficient with a uniform value ε diss = 5 in all test cases in this paper. The same modeling of relaxation time in case of numerical shock wave is firstly proposed in GKS [21].

Stencils and linear compact reconstruction
The compact linear reconstruction from second to fourth order of accuracy is given in this section. Although the second-order scheme in one-and two-dimensional cases is well developed, the construction of a second-order scheme on tetrahedral mesh is not trivial due to its characteristics of the geometry. The first-level neighbors of a tetrahedron may not fill up the space around it, as shown in Fig. 1, and the domain of influence to the central tetrahedron cannot be fully covered by the direct neighboring mesh cells. As a result, on tetrahedral mesh with the stability requirement a secondorder finite volume scheme should use a relatively large stencil which includes the first-level and second-level neighboring cells [27]. In order to present a complete picture for the reconstruction, a second-order compact reconstruction is also presented in this section.

Compact reconstruction
A schematic of reconstruction stencils of second-to fourth-order compact GKS are shown in Fig. 1. The cell with black edges is the reconstructed cell. The cells with blue edges are the first-level neighboring cells of the reconstructed cell, and the cells with green edges are the second-level neighbors of the reconstructed cell. For the sake of simplicity and clarity of illustration, only two cells among the second-level neighboring cells are shown. In general, a first-level neighboring cell is connected to three second-level neighboring cells.
A reasonable stencil consistent with the physical domain of dependence should consist of the reconstructed cell and all neighboring cells sharing common nodes with it. The definition of compact scheme in this study is different from that of conventional compact scheme. In this study, the compactness means that the reconstruction stencil is consistent with the physical domain of dependence. Considering the simplicity of the algorithm and a finite order of accuracy in reconstruction, many subset stencils with different sizes can be defined from the largest complete and compact stencil. In the current compact GKS, inside each control volume, one cell-averaged flow variable and 3 cell-averaged derivatives are available. The stencils and adopted data for the compact reconstruction are determined from the following consideration.
1. The compact stencils are the subsets of the largest compact set consisting of the reconstructed cell and its neighboring cells with common nodes. 2. For the rth-order compact reconstruction, the number of adopted data is about N p = 2N DOF , where N DOF is the number of DOFs of the reconstructed rth-order polynomial. 3. The data for reconstruction is selected in the sequence from the reconstructed cell to the first-level neighbors and then to the second-level neighbors. 4. For the same level neighboring cells, the cell-averaged variable has a higher priority than the cell-averaged gradient.
For the reconstructed cell 0, the four first-level neighboring cells are denoted as 1 to 4. In general, each first-level cell is connected to other three second-level neighboring cells through the surfaces, and the numbers of the second-level cells for each first-level cell j (j = 1, 2, 3, 4) are defined as 5 + 3(j − 1), 6 + 3(j − 1) and 7 + 3(j − 1) . For the third-and fourth-order compact reconstruction, the stencils and adopted data can be obtained from where n l is the unit vector along the direction l . Suppose the second-level neighboring cell j 2 is connected to the first-level neighboring cell j 1 , and j 1 and 0 are connected by face Ŵ 0−j 1 . For ∇ l Q j 2 , n l is taken as the unit outer normal vector of the face Ŵ 0−j 1 . At the boundary of the computational domain, the first-level ghost cells are constructed based on the boundary condition. But the next level neighboring cells of the first-level ghost cells are not constructed. Therefore, the reconstruction scheme based on the stencils on the right figure of Fig. 1 cannot be directly applied due to the absence of cells. In this paper, a third-order compact reconstruction will be developed for the boundary cells and the scheme is stable in numerical examples.
Besides the high-order reconstruction, to get back to the second-order nonlinear reconstruction is critical for the flow simulation at discontinuities. Due to the use of the cell-averaged flow variables and their gradients, the second-order reconstruction will be different from the traditional methods with limiters. Biased stencils will be used instead of a central stencil to get second-order compact reconstruction. A total of four biased compact stencils for second-order reconstruction can be obtained.
Since the DOFs (zeroth-order and first-order derivatives) of a linear polynomial can be determined by the cell average and its cell-averaged gradient, the above biased stencils can fully determine linear polynomials. Then, the linear combination of the four linear polynomials can give a second-order compact reconstruction.

Linear system of compact reconstruction
The polynomial function used in the compact reconstruction is written as where a k is the DOF of the reconstruction polynomial. ϕ k is the zero-mean basis defined by where 0 ≤ l, m, n ≤ r and max{l + m + n} = r . ϕ k can make p r automatically satisfy the conservation condition in the reconstruction. h x , h y and h z are the characteristic scales of 0 along the three directions of axes, which take the values h x = h y = h z = h = |� j | 1/3 for isotropic mesh. The division by h in the expansion is to make the condition number of the matrix in the linear system of a k small. p r (x) is constrained by the following conditions where l indicates the directional derivative along n l . Based on the above constraints for p r (x) , the linear system for a k is obtained by the least square (LS) method or constrained least square (CLS) method. In the CLS method, some constraints, such as the one for Q j ( j = 1, 2, 3, 4 ), are strictly satisfied, and others are satisfied in the sense of least square.
For p 1 and p 2 , the CLS method is adopted, where the cell-averaged values Q j (j = 1, 2, 3, 4) are strictly satisfied. While the LS method is adopted for p 3 , which makes the fourth-order compact GKS have better stability on irregular meshes. If the CLS method is adopted for p 3 and the cell-averaged values of the first-level neighbors are strictly satisfied, the linear system of a k is more sensitive to the constraints from the first-level neighboring cells, due to the use of a smaller effective numerical domain of dependence. The CLS problem in the determination of p 1 and p 2 can be solved by the Lagrangian factor method and the linear system for a k can be obtained. A general form of the linear system for a k has been given in [18].
The linear polynomial of the second-order reconstruction can be the linear combination of four p 1 k (x).

Nonlinear compact reconstruction
To deal with discontinuities, nonlinear compact reconstruction with WENO method is adopted. WENO reconstruction can adaptively achieve high-order accuracy in smooth region and essentially non-oscillatory property in discontinuity region. The simplified WENO method is developed for simple implementation and good robustness on unstructured mesh [17,18]. The basis of the adaptivity in the simplified WENO is the nonlinear combination of a high-order polynomial and several lower-order polynomials. The four compact reconstructions p 1 k (k = 1, 2, 3, 4) on the biased stencils can be used as the lower-order candidate polynomials in the WENO method.

Simplified WENO reconstruction
The extension of the one-dimensional WENO reconstruction to unstructured mesh is difficult [28], especially for high-order (order ≥ 4 ) one. The simplified WENO reconstruction has been developed and implemented on triangular mesh [18]. The simplified WENO reconstruction is given as The nonlinear weight w k is and the linear weight d k is where ǫ is a small positive number and ǫ = 1 × 10 −15 is taken for all numerical tests in this paper. n is the number of the sub-stencils. A large number α improves the robustness of the scheme through better identifying the less smooth polynomials from all candidate polynomials. α = 3 is taken in the current compact reconstruction. τ Z is the local high-order reference value to indicate smoothness of the large stencil and it is given by IS k , where IS k is obtained by the conventional definition in [29,30]. The parameters C and C k are required to satisfy n k=1 C k = 1, C > 0 . C = n and C k = 1/n are taken in this paper. The compact scheme based on the simplified WENO reconstruction is insensitive to the values of C and C k .
Based on the concept for reconstruction without crossing discontinuity of the solution, the biased lower-order polynomials p k in the WENO reconstruction take four biased linear polynomials p 1 k and one zeroth-order polynomial p 0 . The smooth indicator of p 0 is obtained based on the local smoothest slope W j,x , W j,y and W j,z ( j = 0, 1, 2, 3, 4 ) by the definition in [30]. Such a construction guarantees that the smoothness indicator of p 0 corresponds to an auxiliary smoother linear reconstruction than linear reconstruction p 1 k . The calculation of IS based on the definition in [30] will be complicated for the highorder polynomial (order ≥ 3 ) on the tetrahedral mesh. In this paper, numerical quadrature is adopted to calculate IS. For a third-order polynomial p 3 , its IS can be where x 0 is the centroid of j . |α| is a multi-index, for example, when |α| = 1 , there are cases (α 1 , α 2 , α 3 ) = (1, 0, 0), (0, 1, 0) and (0, 0, 1). Apparently, two simplified calculations are made in Eq. (24). Firstly, the numerical quadrature with second-order accuracy is adopted for the integral. Secondly, terms related to the third derivative by taking |α| = 3 in IS are ignored, and only terms related to the first and second derivatives are included. In smooth regions, the obtained IS approximates the analytical one as where A i are the parameters dependent on cell's geometry, and W i are three first derivatives. The degradation of IS accuracy will not affect the accuracy of the final nonlinear reconstruction.

WENO weight with improved adaptivity
Adaptive variation of the magnitude of τ Z in WENO weight is important for the accuracy and robustness of the WENO scheme. τ Z is required to adaptively satisfy τ Z /IS k = O(h r )(r > 0) for the candidate polynomials for smooth solutions, and it satisfies τ Z /IS k = O(1) for the candidates crossing discontinuities. Thus the definition of τ Z is directly related to the accuracy and robustness of the scheme. For the current high-order simplified WENO reconstruction (accuracy order is higher than 2), the orders of the candidate polynomials include zero, one and two or three. As a result, the separation of order of the candidate polynomials results in a large deviation in the values of smoothness indicators. For example, the fourth-order nonlinear reconstruction is obtained from p 3 , p 1 and p 0 . A small first derivative of p 1 or large second and third derivatives of p 3 will lead to large deviations in the values of their IS k even in smooth region, and a small first derivative of p 1 is common at critical point.
In this paper, τ Z is defined by the hierarchical differences of IS k . τ Z consists of two parts. The first part is the difference of two IS k which are more likely to have small deviation in values in smooth region. To further identify discontinuity, the second part is obtained by two IS k whose values are more likely to have large deviation. τ Z can be uniformly written as where β is an adaptive power. For the fourth-order reconstruction, τ Z can be given as IS 0 is from the candidate polynomial p 0 , IS * is the auxiliary one with smaller deviation from IS 0 in values, and IS * * is the auxiliary one with larger deviation with IS * in the vicinity of discontinuities. IS * is given by Eq. (24) based on p 2 which is obtained from the third-order compact reconstruction. IS * * is taken as IS * * = ( 4 k=1 IS k − max 1≤k≤4 IS k − min 1≤k≤4 IS k )/2 which corresponds to smooth candidate polynomials (assuming only one discontinuity exists at four faces of the reconstructed cell) but not the smoothest one whose IS k may be close to 0. The first term of τ Z in Eq. (27) is going to be a higher-order small value, even at critical points. While the second term of τ Z in Eq. (27) can be O(1) in the vicinity of discontinuities. For the thirdorder and second-order reconstructions, the same τ Z as in Eq. (27) is used, but different IS * is defined. In this case, IS * is taken as IS * = max 1≤k≤4 {IS k } . β is defined as where the free-parameters are uniformly taken as C 1 = 1 and C 0 = min{sgn(ρ thres − ρ 0 ) + 1, sgn(p thres − p 0 ) + 1, 1} . ρ 0 and p 0 are the minimum cellaveraged density and pressure in the reconstructed cell and its first-level neighboring cells. ρ thres and p thres are uniformly taken as ρ thres = p thres = 5.0 × 10 −2 in this paper.
The improved adaptivity of τ Z is achieved by the hierarchical difference of indicators and the adaptive power β . The high-order part τ HO Z of τ Z can be a high-order small value in smooth region, and the lower-order part τ LO Z of τ Z can be reduced to a high-order small value in smooth region by taking β = 2 . At the same time, discontinuities of variables on the stencil can be identified in the current reconstruction sensitively. The first part and the second part will take values of order O(1) in discontinuities, which make the nonlinear weights of smoother polynomials candidates have values O(1). In particular, C 0 in the definition of β improves the robustness of the scheme in the regions with a very small value of density or pressure.
In Section 3, the S2O4 temporal discretization based on nonlinearly limited flux function is introduced, where the higher-order time derivatives of the flux function are nonlinearly limited. The nonlinear weight in Eq. (14) depends on the smoothness indicators in the spatial WENO reconstruction, and these indicators are given by

Numerical examples
In this section, a few test cases on 3D tetrahedral meshes will be conducted. The time step is determined by the CFL condition. In all test cases, CFL number is taken as C fl ≥ 0.5 , except C fl = 0.35 in the Taylor-Green vortex flow. For viscous flows, the time step is also limited by the viscous term as �t = C fl h 2 /(3ν) , where h is the cell size and ν is the kinematic viscosity coefficient. The collision time τ for inviscid flow at a cell interface is defined by where ε = 0.05 , ε diss = 10 , 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 region, it will reduce to τ = µ/p . The reason for including pressure jump term in the particle collision time is to increase the shock wave thickness numerically to the order of cell size [21].

Accuracy test
The three-dimensional advection of density perturbation is used to verify the order of accuracy of compact GKS. The initial condition is given by  The computational domain is [0, 2] 3 . The periodic boundary conditions are applied on all domain boundaries. The tetrahedral mesh is used in the computation. The tetrahedral mesh is obtained by dividing each regular hexahedral cell in the Cartesian mesh into six tetrahedral cells. The coarsest mesh used in the computation is shown in Fig. 2. The L 1 and L ∞ errors and convergence orders obtained by the compact GKS with linear reconstruction at t = 2 are presented in Table 1. Due to the non-uniform mesh cell, the accuracy order of L ∞ cannot reflect the true convergence order of the scheme. From the numerical results listed in Table 1, it can be seen that the accuracy orders of L 1 are almost the same as the theoretical ones. The errors versus mesh DOFs is given in Fig. 3 for linear and nonlinear compact GKS, where the abscissa represents the mesh DOFs in one direction. Compared with the second-order scheme, the high-order scheme has a smaller error under the same mesh DOFs; at the same time, under the same error requirement, the mesh DOF used ρ(x, y, z) = 1 + 0.2 sin(π(x + y + z)), p(x, y, z) = 1, U(x, y, z) = V (x, y, z) = W (x, y, z) = 1. by the high-order scheme is far less than that of the second-order one. For example, when the error limit is 10 −3 , the required mesh DOFs by the linear fourth-order scheme are nearly two orders of magnitude less than that by the linear second-order scheme, and the required mesh DOFs by the nonlinear fourth-order scheme is 1/35 of that by the nonlinear second-order scheme. Figure 4 lists the relationship between error and CPU time. The computation is performed by an OpenMP parallel code using 48 threads on a 2.2 GHz Intel(R) Xeon(R) workstation. The curves of errors versus CPU time are basically similar to the curves of errors versus mesh DOFs in one direction. The results indicate that the increase in CPU time relevant to the increase in algorithm complexity of high-order compact scheme accounts for a small proportion in the overall CPU time in comparison with the secondorder scheme.

Taylor-Green vortex flow
Taylor-Green vortex flow is a popular test to assess high-order schemes [31,32], which is used to verify the accuracy of the schemes for viscous flow. Starting from a smooth initial flow distribution, small scale flow structure in the flow field will emerge and evolve continuously. The initial condition is set as where U 0 = 1 , ρ 0 = 1 . The Mach number is Ma = 0.1 and it is determined by Ma = U 0 / √ γ p 0 /ρ 0 . The Reynolds number is Re = 280 defined by Re = ρ 0 U 0 /µ , where µ is the dynamic viscosity coefficient. The computational domain is [−π , π ] 3 . This test case is used to verify the accuracy and the linear stability of the current compact GKS, and only the linear compact GKS is adopted in the computation. In order to compare the performance of second-order and fourth-order compact GKS, a coarse mesh with 20 3 × 6 cells and a fine mesh with 40 3 × 6 cells are used in the computation, with the same tetrahedral mesh used in the previous accuracy test. The averaged kinetic energy is defined as where is the computational domain. E k is calculated by numerical quadrature. For the fourth-order scheme, the five-point Gaussian quadrature formula is used, where the values of flow variables at the quadrature points are obtained by the same fourth-order compact reconstruction. For the second-order scheme, the mid-point quadrature formula is adopted. The dissipation rate of kinetic energy is given by ε(E k ) is calculated by the central difference method with second-order accuracy by considering the small difference in the time steps. The integrated enstrophy is defined as where ζ is calculated by numerical quadrature, and the velocity derivative values at quadrature points are obtained by the compact reconstruction. Figure 5 presents the iso-surface of the second invariant of the velocity gradient tensor Qv = 0.18 at t = 15 colored by pressure. The left and right figures in Fig. 5 are the results from the second-order and fourth-order compact schemes, respectively. The fourth-order compact GKS has a better resolution than that from the second-order scheme. The time history of E k , ε(E k ) and ζ are shown in Fig. 6. The reference solution is from [31]. The results of E k and ε(E k ) obtained by the second-order compact GKS on the fine mesh is slightly better than that obtained by the fourth-order compact GKS on the coarse mesh. The results of E k and ε(E k ) obtained by the fourth-order compact GKS on the fine mesh has good agreement with the reference solution. For the enstrophy ζ , the result of the fourth-order scheme on the fine mesh is slightly better than that of the second-order scheme, but they are obviously better than that of the fourth-order scheme on the coarse mesh.   Fig. 8. In order to verify that the nonlinear scheme in this paper has similar accuracy as the linear scheme, the fourth-order compact GKS with WENO reconstruction is also used in this test case. The quantitative results of the compact GKS and other computations are given in Table 2, where the drag and lift coefficients, the wake length L, and the total number of mesh cells are listed for comparison. T-mesh, H-mesh and M-mesh represent tetrahedral mesh, hexahedral mesh and hybrid mesh respectively. The linear fourth-order GKS gives the same value of C D as direct DG scheme [34] and 5th-order hybrid scheme of FR and DG [35], but uses fewer total DOFs, where (r + 1)(r + 2)(r + 3)/6 DOFs are used in each cell for (r + 1)th-order DG scheme. The nonlinear fourth-order GKS gives the C D with a 0.17% error from the linear scheme, which validates the similar high-order accuracy of nonlinear fourth-order compact GKS to the linear one.

Riemann problems and blast wave problem
To validate the robustness of the compact GKS, the scheme is applied to two 1-D Riemann problems and the blast wave problem for compressible inviscid flow. The computational domain is given by [0, 1] × [0, 0.025] 2 . The tetrahedral mesh is used in the computation, and the average side length of tetrahedral cells is about 1/200. The initial conditions and the output times of the two Riemann problems are listed in Table 3. For the Lax shock tube problem [36], the fixed inflow condition and free boundary condition are used for the boundaries with x = 0 and x = 1 respectively, and the inviscid wall condition is imposed on the other boundaries. For the large pressure ratio problem [37], the free boundary condition is used on both ends at x = 0 and x = 1 , and the inviscid wall condition is imposed on the other boundaries.
The initial condition of the blast wave problem [38] is given as The output time is t f = 0.038 . The inviscid wall boundary condition is imposed on all boundaries for the blast wave problem.
(1, 0, 100), 0.9 ≤ x ≤ 1. The local enlargement of the mesh and the 3-D contours of density are shown in Fig. 9. The wave structures are accurately obtained by the fourth-order compact GKS on the tetrahedral meshes. The computed density and velocity profiles of these cases are shown in Figs. 10, 11 and 12, respectively. No numerical oscillations are observed near the shock fronts. The compact GKS performs well on the large pressure jump and the blast wave problems. On the current coarse mesh, especially tetrahedral one, the result of large pressure jump problem is reasonable in comparison with the results in [37], even the exact Riemann solver with shock fitting may not be helpful in this case on the tetrahedral mesh [39]. The 3-D density distribution and density distribution along lines are shown in Fig. 13, where the reference solution is from [40]. Line 1, line 2 and line 3 are determined by connecting the origin to (0, 1.2, 1.2), (1.2, 0, 1.2) and (1.2, 1.2, 0). The strong robustness of the 4th-order compact GKS is demonstrated by its use of a large CFL number CFL = 0.6 without additional limiting technique. In addition, the high resolution of the compact GKS for strong shock waves is verified by the high post-shock density peak, and the numerical shock wave remains sharp and spans only two mesh cells. The noncompact high-order GKS gives a numerical shock wave that is wider and has a lower post-shock peak [41].

3-D explosion and implosion problems
The 3-D Noh problem is an implosion test to model the gas compression with constant radial velocity towards a spherical center, where a moving strong shock wave is formed [42]. The computational domain is [0, 0.3] 3 . The initial density and pressure are ρ = 1 and p = 1 × 10 −4 , and the velocity is (U , V , W ) = (−x, −y, −z)/ x 2 + y 2 + z 2 . The ratio of the specific heat is γ = 5/3 . The inviscid wall condition is adopted along the boundaries x = 0 , y = 0 and z = 0 . The supersonic inflow boundary condition is imposed on the other boundaries with the same pressure and velocity as the initial conditions and the analytical solution of density,  0.3, 0.3, 0). Again, in this test case, a large CFL number CFL = 0.6 without additional limiting technique is used. The accuracy of the compact GKS for strong shock waves is verified by the precise post-shock density solution which is comparable to the result of the non-compact 5th-order GKS on structured mesh [43].

Supersonic and hypersonic inviscid flow around a sphere
The inviscid flow at high Mach numbers around a sphere is used to verify the strong robustness of the fourth-order compact GKS. The incoming inviscid flows have Mach numbers Ma = 2 and Ma = 20 separately. The adiabatic reflective boundary condition is imposed on the wall of the sphere. The same computational mesh as shown in Fig. 7 is used in the current computation. The 3-D pressure and 2-D Mach number distributions are presented in Fig. 15 when the flow reaches a steady state. In this test ρ = 64, r < t/3, (1 + t/r) 2 , r > t/3. (down). The fourth-order compact GKS is used. The mesh used in the computation is the same as shown in Fig. 7 case, the exact same code and parameter settings are used as the nonlinear reconstruction-based GKS in the accuracy test of Fig. 4.

Supersonic flow over the YF-17 fighter model
The supersonic flow over the YF-17 model is simulated by the fourth-order GKS under complex geometry. The mesh of the YF-17 model is provided at https:// cgns. github. io/ CGNSF iles. html, which is shown in Fig. 16. The total number of tetrahedral cells is about 1.7 × 10 5 . The incoming Mach number is set as Ma ∞ = 2.0 , and the angle of attack is AoA = 0 . The inviscid wall boundary condition is adopted on the model surface and on the symmetry plane.  Figure 17 presents the pressure and Mach number distributions on the model at a steady state. A low pressure area occurs in the downwind area at the extreme end of the model fuselage. Shock waves are generated at the nose of the model fuselage and the front of the wing.

Conclusion
In this paper, high-order compact GKS from second-up to fourth-order have been constructed on three-dimensional tetrahedral mesh for compressible flow simulations. The compact scheme works very well from the subsonic smooth flow to the hypersonic compressible flow simulations. The scheme shows the accuracy/efficiency in the smooth viscous flow computation and excellent robustness for the complex flow with discontinuities. More importantly, a large CFL number, such as CFL = 0.6, can be used for the fourth-order compact GKS, even at a Mach number 20 flow computation and Noh problem on tetrahedral mesh. The success of the compact scheme is mainly coming from the following fact in the algorithm development. A high-order time-accurate gas evolution model at a cell interface is adopted and it provides the evolution solution of flow variables and fluxes at cell interfaces, which can be used to update cell-averaged flow variables and their gradients. More importantly, in the high-order evolution model a nonlinear limiter on the high-order time derivative of flux function through WENO formulation is introduced in its temporal evolution, which makes the status of "nonlinear limiters" on the equivalent footing in the spatial reconstruction and temporal evolution for capturing the propagation of discontinuous solution in space and time. As a result, even for the fourth-order compact scheme, a large CFL number can be used, such as CFL = 0.6 in the current scheme in comparison with the 0.14 in the DG method. At the same time, based on the cell-averaged flow variables and their gradients improved WENO weighting functions in the compact reconstruction have been developed. The extension of the current compact GKS to three-dimensional mixed unstructured mesh is straightforward. It is expected that the high-order compact GKS will play an important role in engineering applications.