Performance evaluation of high-order compact and second-order gas-kinetic schemes in compressible flow simulations

TL;DR

This study compares high-order compact and second-order gas-kinetic schemes in compressible flow simulations, demonstrating that high-order schemes offer superior accuracy and efficiency, especially in turbulent flows with shocks and vortices.

Contribution

It provides the first verification of high-order compact gas-kinetic schemes' advantages in viscous flow simulations with discontinuities.

Findings

CGKS-5th achieves similar resolution to GKS-2nd at lower computational cost

CGKS-5th delivers higher accuracy in turbulent flows with shocks

Multi-GPU parallelization enables large-scale applications

Abstract

The trade-off among accuracy, robustness, and computational cost remains a key challenge in simulating complex flows. Second-order schemes are computationally efficient but lack the accuracy required for resolving intricate flow structures, particularly in turbulence. High-order schemes, especially compact high-order schemes, offer superior accuracy and resolution at a relatively modest computational cost. To clarify the practical performance of high-order schemes in scale-resolving simulations, this study evaluates two representative gas-kinetic schemes: the newly developed fifth-order compact gas-kinetic scheme (CGKS-5th) and the conventional second-order gas-kinetic scheme (GKS-2nd). Test cases ranging from subsonic to supersonic flows are used to quantitatively assess their accuracy and efficiency. The results demonstrate that CGKS-5th achieves comparable resolution to GKS-2nd at roughly an order of magnitude lower computational cost. Under equivalent computational resources, CGKS-5th delivers significantly higher accuracy and resolution, particularly in turbulent flows involving shocks and small-scale vortices. This study provides the first clear verification of the advantages of high-order compact gas-kinetic schemes in simulating viscous flows with discontinuities. Additionally, multi-GPU parallelization using CUDA and MPI is implemented to enable large-scale applications.

Tables

Table 1. Table 1: Subsonic Taylor-Green vortex flow: the comparison of computational efficiency between linear CGKS-5th and linear GKS-2nd.

Number of GPUs	MESH	Scheme	Computational time	Ratio
2	${96}^3$	linear CGKS-5th	876	1
2	${296}^3$	linear GKS-2nd	7746	8.8

Table 2. Table 2: High-speed Taylor-Green vortex flow: the comparison of computational efficiency between CGKS-5th and GKS-2nd.

Number of GPUs	MESH	Scheme	Computational time	Ratio
4	${116}^3$	nonlinear CGKS-5th	407	1
4	${384}^3$	nonlinear GKS-2nd	2842	7.0

Table 3. Table 3: Turbulent Flow Past a Cylinder: the evaluation of the computational efficiency of CGKS-5th and GKS-2nd for t ∈ [ 0 , 10 ] t\in[0,10] with a computational cost difference within 5 % \% .

Number of GPUs	MESH	Scheme	Time	Ratio
8	$352\times 176\times 176$	CGKS-5th	7502.3	1
8	$528\times 264\times 264$	GKS-2nd	7115.9	0.95

Table 4. Table 4: Transonic Jet Simulation: the comparison of L C L_{C} calculated by CGKS-5th and GKS-2nd with reference data.

Method	$L_C$
Exp. Arakeri et al. [44]	7.7D
Exp. Power et al. [43]	7.0D
LES Fosso Pouangué et al. [42]	6.2D
GKS-2nd	5.6D
CGKS-5th	6.6D

Equations

$$f_t+\mathbf{u}\cdot \nabla f=\frac{g-f}{\tau },$$

$$g=\rho {\left(\frac{\lambda }{\pi }\right)}^{(N+3)/2}{e}^{-\lambda [(\mathbf{u}-\mathbf{U}{)}^2+{\mathbf{\xi }}^2]},$$

$$\int \frac{g-f}{\tau }\mathbf{\psi }\text{d}\Xi =0,$$

$$\mathbf{Q}=\int f\mathbf{\psi }\text{d}\Xi ,$$

$$\mathbf{F}=\int \mathbf{u}f\mathbf{\psi }\text{d}\Xi ,$$

$$f({\mathbf{x}}_G,t,\mathbf{u},\mathbf{\xi })=\frac{1}{\tau }{\int }_0^{t}g({\mathbf{x}}^{\prime },t^{\prime },\mathbf{u},\mathbf{\xi })e^{-(t-t^{\prime })/\tau }\text{d}{t}^{\prime }+e^{-t/\tau }{f}_0(-\mathbf{u}t,\mathbf{\xi }).$$

$$f({\mathbf{x}}_G,t,\mathbf{u},\mathbf{\xi })=$$

$$+$$

$$\hat{f}(t)=f^n+t{f}_t^n.$$

$$f^n$$

$$f_t^n$$

$${\mathbf{F}}^n=\int \mathbf{u}{f}^n\mathbf{\psi }\text{d}\Xi ,{\mathbf{F}}_t^n=\int \mathbf{u}{f}_t^n\mathbf{\psi }\text{d}\Xi .$$

$${\mathbf{Q}}^n=\int f^n\mathbf{\psi }\text{d}\Xi ,{\mathbf{Q}}_t^n=\int f_t^n\mathbf{\psi }\text{d}\Xi .$$

$$\frac{\text{d}{\mathbf{Q}}_i}{\text{d}t}={\mathcal{L}}_i(\mathbf{Q}),$$

$${\mathcal{L}}_i(\mathbf{Q})=-\frac{1}{|{\Omega }_i|}{\int }_{\partial {\Omega }_i}\mathbf{F}(t)\mathbf{n}\text{d}S,$$

$${\int }_{{\sigma }_j}\mathbf{F}(t)\mathbf{n}\text{d}S=\sum _{k=1}^K{\omega }_k\mathbf{F}({\mathbf{x}}_k,t)S_j,$$

$$\nabla \mathbf{Q}(t)=\frac{1}{|\Omega |}{\int }_{\Omega }\nabla \mathbf{Q}(t)\text{d}V=\frac{1}{|\Omega |}{\int }_{\partial \Omega }\mathbf{Q}(t)\mathbf{n}\text{d}S,$$

$$({\partial }_l\mathbf{Q})(t)$$

$$=\frac{1}{\left|l_G\right|}\left(\mathbf{Q}({\mathbf{x}}_2,t)-\mathbf{Q}({\mathbf{x}}_1,t)\right),$$

$$S^{cell-5th}=\{{\Omega }_0,{\Omega }_{1,m},{\Omega }_{2,n}|m=1,\dots ,6,n=1,\dots ,12\}.$$

$$S^{data-5th}=S^0\cup S^{face}\cup S^{edge},$$

$$S^0$$

$$S^{face}$$

$$S^{edge}$$

$$P^4(\mathbf{x})={\mathbf{Q}}_0+\sum _{|\mathbf{d}|=1}^4{a}_{\mathbf{d}}{p}_{\mathbf{d}}(\mathbf{x}),$$

$$p_{\mathbf{d}}(\mathbf{x})=\frac{1}{d_1!d_2!d_3!}{\delta x}^{d_1}{\delta y}^{d_2}{\delta z}^{d_3}-\frac{1}{\left|{\Omega }_0\right|}{\int }_{{\Omega }_0}\frac{1}{d_1!d_2!d_3!}{\delta x}^{d_1}{\delta y}^{d_2}{\delta z}^{d_3}\text{d}V,$$

$$\delta x=\frac{1}{h_x}(x-x_0),\delta y=\frac{1}{h_y}(y-y_0),\delta z=\frac{1}{h_z}(z-z_0),$$

$$\begin{array}{rl}\frac{1}{\left|{\Omega }_k\right|}{\int }_{{\Omega }_k}{P}^4(\mathbf{x})\text{d}V& ={\mathbf{Q}}_k,{\Omega }_k\in S^{cell-5th},\\ \frac{h_1}{\left|{\Omega }_{1,m}\right|}{\int }_{{\Omega }_{1,m}}\nabla P^4(\mathbf{x})\text{d}V& =h_1\cdot (\nabla \mathbf{Q}{)}_{1,m},m=1,\dots ,6,\\ \frac{h_2}{\left|{\Omega }_{2,n}\right|}{\int }_{{\Omega }_{2,n}}\frac{\partial }{{\partial }_{\mathbf{\tau }}}{P}^4(\mathbf{x})\text{d}V& =h_2\cdot ({\partial }_{\mathbf{\tau }}\mathbf{Q}{)}_{2,n},n=1,\dots ,12,\\ \frac{h_3}{\left|l_G\right|}{\int }_{l_G}\frac{\partial }{{\partial }_l}{P}^4(\mathbf{x})\text{d}l& =h_3\cdot ({\partial }_l\mathbf{Q}{)}_{0,G},\end{array}$$

$$S_m^{cell-2}=\{{\Omega }_0,{\Omega }_{1,m}\}.$$

$$S_m^{data-2}=\{{\mathbf{Q}}_0,({\partial }_l\mathbf{Q}{)}_{0,G}^m,{\mathbf{Q}}_{1,m}\},m=1,\dots ,6,$$

$$P_m(\mathbf{x})={\mathbf{Q}}_0+\sum _{|\mathbf{d}|=1}{b}_{\mathbf{d}}^m{p}_{\mathbf{d}}(\mathbf{x}),$$

$$\begin{array}{rl}\frac{1}{\left|{\Omega }_{1,m}\right|}{\int }_{{\Omega }_{1,m}}{P}_m^1(\mathbf{x})\text{d}V& ={\mathbf{Q}}_{1,m},\\ \frac{1}{\left|l_G\right|}{\int }_{l_G}\frac{\partial }{{\partial }_l}{P}_m^1(\mathbf{x})\text{d}l& =({\partial }_l\mathbf{Q}{)}_{0,G}^m,\end{array}$$

Full text

1 Introduction

To address complex flow simulations in engineering applications, it is essential to balance resolution, robustness, and computational cost. Most commercial and industrial software currently relies on second-order numerical methods, such as MUSCL [1], owing to their robustness, low computational cost, and broad applicability. However, second-order schemes are inherently constrained by significant numerical dissipation and dispersion, which limit their accuracy in high-resolution scenarios. For example, accurately predicting the position and intensity of shock waves in high-speed flows, capturing wave propagation and reflection in acoustic problems, and resolving multi-scale turbulence structures in direct numerical simulations (DNS) all expose the limitations of second-order methods.

To overcome these challenges, high-order schemes have been developed to provide improved accuracy and resolution for complex flow problems. Over the past few decades, various high-order methods have been proposed, including finite difference (FD) methods [2, 3], essential non-oscillatory (ENO) schemes [4, 5], weighted essential non-oscillatory (WENO) schemes [6, 7], Hermite WENO (HWENO) schemes [8], discontinuous Galerkin (DG) methods [9, 10], flux reconstruction (FR) methods [11], correction procedure using reconstruction (CPR) methods [12], and compact gas-kinetic scheme (CGKS) [13, 14, 15, 16, 17]. Finite difference and finite volume methods typically achieve high-order accuracy by using expanded stencils that incorporate additional neighboring information for reconstruction. In contrast, DG-type methods improve accuracy by increasing the degree of interpolation polynomials within each cell, thereby enhancing the internal degrees of freedom and naturally exhibiting compactness. However, the widespread adoption of high-order schemes in industrial engineering is limited by lower robustness, higher computational costs, and greater implementation complexity. Achieving an optimal balance among resolution, robustness, and computational efficiency is essential to fully exploit the potential of high-order schemes in practical applications.

In recent decades, the gas-kinetic scheme (GKS) has been systematically developed for simulating flows across a wide range of regimes, from subsonic to supersonic [18, 19, 20]. GKS performs a one-step, second-order-accurate temporal evolution from kinetic to hydrodynamic scales and offers several intrinsic advantages: it computes both inviscid and viscous fluxes within a single unified formulation and adaptively transitions between equilibrium-state fluxes in smooth regions and non-equilibrium transport fluxes suitable for discontinuities. In conjunction with linear reconstruction and slope limiters, the conventional second-order GKS can be obtained. Building on this foundation, the high-order compact GKS (CGKS) further exploits the gas distribution function to achieve time-accurate evolution of flow variables at cell interfaces. As a result, both the cell-averaged conservative variables and their gradients within each control volume can be updated simultaneously at the next time level. Coupled with the two-stage temporal discretization method [21, 22, 23], high-order CGKS have been developed for high-fidelity simulation of compressible flows [13, 14, 15, 16, 17]. For structured meshes, a high-order CGKS was proposed in [13], demonstrating spectral-like resolution in one dimension using a dimension-by-dimension reconstruction strategy. However, this approach substantially increases algorithmic complexity in three-dimensional applications. High-order GKS in curvilinear coordinates have been developed using dimension-by-dimension reconstruction for both laminar [24] and turbulent flows [25]. Nevertheless, when extended to compact schemes in three dimensions, the dimension-by-dimension strategy still introduces significant implementation challenges. To achieve high resolution and efficiency while maintaining implementation simplicity, a multidimensional fifth-order CGKS has been newly developed for non-orthogonal structured meshes [17]. This method employs line-averaged derivatives within each control volume to provide additional degrees of freedom. Its resolution, robustness, and stability have been validated through a series of complex flow simulations.

In this paper, we conduct a systematic performance evaluation of the newly developed fifth-order compact gas-kinetic scheme (CGKS-5th) against the conventional second-order gas-kinetic scheme (GKS-2nd) on structured meshes. Two comparative strategies are employed: (i) assessing the computational cost required to achieve a target resolution, and (ii) evaluating accuracy and resolution under comparable computational resources. While the advantages of high-order schemes over second-order schemes have been extensively validated in complex smooth flows [26, 27, 28], such verification remains limited for flows involving discontinuities. Theoretically, the use of nonlinear reconstructions or limiters in discontinuous flows may impact the performance of high-order schemes. This study seeks to address this gap. To achieve this, a series of numerical tests, ranging from subsonic to supersonic turbulent flows, are conducted to compare the resolution and efficiency of the two schemes under these criteria. Furthermore, both schemes are implemented on multiple graphics processing units (GPUs) using the Compute Unified Device Architecture (CUDA) and the Message Passing Interface (MPI). Numerical experiments are conducted on NVIDIA GeForce RTX 4090 GPUs for achieving significant computational acceleration.

This paper is organized as follows: Section 2 introduces the gas-kinetic scheme and finite volume framework. Section 3 provides a brief review of the new fifth-order compact gas-kinetic scheme and the second-order gas-kinetic scheme. Numerical examples are presented in Section 4, and conclusions are drawn in the final section.

2 Time-accurate solution of GKS and finite volume method

In the gas-kinetic method, time-dependent numerical fluxes are computed from a physically modeled, space-time coupled solution of the BGK equation for the gas distribution function [18, 19]. In the gas-kinetic method, a physically modeled, space-time coupled solution of the BGK equation for the gas distribution function is employed to compute time-dependent numerical fluxes [18, 19]. The BGK equation is a simplification of Boltzmann equation, and the three-dimensional BGK equation [29, 30] can be written as

[TABLE]

where $\bm{u}=(u,v,w)$ is the particle velocity, $\tau$ is the collision time and $f$ is the gas distribution function. The equilibrium state $g$ is given by Maxwellian distribution

[TABLE]

where $\rho$ is the density, $\bm{U}=(U,V,W)$ is the macroscopic fluid velocity, and $\lambda=\displaystyle\frac{\rho}{2p}$, $p$ is the pressure. In the BGK equation, the collision operator involves a simple relaxation from $f$ to the local equilibrium state $g$. The variable $\bm{\xi}$ accounts for the internal degree of freedom of molecular motion, $\bm{\xi}^{2}=\xi_{1}^{2}+\dots+\xi_{N}^{2}$, where $N=(5-3\gamma)/(\gamma-1)$ is the number of internal degree of freedom and $\gamma$ is the specific heat ratio. The collision term satisfies the compatibility condition

[TABLE]

where $\displaystyle\bm{\psi}=\left(1,u,v,w,\frac{1}{2}(u^{2}+v^{2}+w^{2}+\bm{\xi}^{2})\right)^{T}$ and $\text{d}\Xi=\text{d}u\text{d}v\text{d}w\text{d}\xi_{1}\dots\text{d}\xi_{N}$. The macroscopic conservative variables $\bm{Q}=(\rho,\rho U,\rho V,\rho W,\rho E)$ can be calculated through the gas distribution function $f$

[TABLE]

and the corresponding fluxes can be given by taking moments of the gas distribution function

[TABLE]

where $\bm{F}$ represents either the Euler flux or the Navier-Stokes (NS) flux, depending on the order of approximation of $f$ to $g$. This study focuses solely on viscous flows described by the NS equations [18, 19]. In GKS, $\bm{F}$ at cell interface is determined by the gas distribution function $f$. Based on the integral solution of BGK equation Eq.(1), the gas distribution function $f(\bm{x}_{G},t,\bm{u},\bm{\xi})$ can be given by

[TABLE]

By modeling and approximating the unknown terms in the integral solution, a second-order accurate explicit solution for the distribution function is obtained as

[TABLE]

where $g_{0}$, $g_{l}$ and $g_{r}$ are determined by the conservative variables reconstructed at the interface, and the coefficients are determined by the macroscopic variables and their derivatives. The specific formulas can be referred to in [18, 19].

For the second-order evolution solution, the time accurate distribution function $f(t)$ at cell interface can be approximated through a linearization in time [31]

[TABLE]

The two coefficients $f^{n}$ and $f_{t}^{n}$ are calculated as follows

[TABLE]

where $\bar{f}(\Delta t)$ and $\bar{f}(\Delta t/2)$ are the time integrations of $f(t)$ over the interval $[t_{n},t_{n}+\Delta t]$ and $[t_{n},t_{n}+\Delta t/2]$, respectively. The numerical fluxes and their time derivatives can be obtained by taking moments of $\hat{f}(t)$ and $\hat{f}_{t}(t)$ at $t=t_{n}$

[TABLE]

Simultaneously, the flow variables and their time derivatives can be given by

[TABLE]

The obtained conservative variables at interfaces are used to achieve compact reconstructions in high-order compact GKS.

Taking moments of Eq.(1) and integrating with respect to space, the semi-discretized finite volume scheme can be obtained as

[TABLE]

where $\bm{Q}_{i}$ is the cell-averaged conservative variables over cell $\Omega_{i}$. The operator $\mathcal{L}_{i}$ is defined as

[TABLE]

where $|\Omega_{i}|$ is the volume of $\Omega_{i}$, $\partial\Omega_{i}$ represents the cell interfaces of $\Omega_{i}$, and $\bm{n}$ is the unit direction of the cell interface. To achieve the expected order of accuracy, the Gaussian quadrature is used for the flux integration

[TABLE]

where $S_{j}$ is the area of $\sigma_{j}$, $\sigma_{j}$ is one of the cell interfaces, $\omega_{k}$ is the Gaussian quadrature weight, $\bm{x}_{k}$ is the Gaussian quadrature point, $K$ is the total number of Gaussian quadrature points. The numerical flux $\bm{F}(\bm{x}_{k},t)$ at Gaussian quadrature point can be given by Eq.(3). In this study, structured meshes are utilized, where each face is a quadrilateral with $K=4$ Gaussian points for high-order schemes. For the second-order scheme, this integral reduces to a form that only involves the midpoints of the cell interfaces.

3 High-order compact and second-order gas-kinetic schemes

Second-order schemes have been widely used to solve a broad range of engineering problems owing to their robustness and computational efficiency. However, their accuracy is often inadequate for accurately capturing complex flow features in high-fidelity simulations. In contrast, high-order schemes, due to their superior accuracy, can faithfully resolve intricate details of the flow field, although this comes at the cost of increased algorithmic complexity. To enable a clear and fair assessment of the differences between high- and low-order methods in simulations, this study conducts a performance evaluation within the framework of the finite-volume method based on GKS. Specifically, we consider a recently proposed CGKS-5th [17] and a conventional GKS-2nd, as representatives of high-order and second-order approaches, respectively. A brief introduction to these two schemes is provided in this section.

3.1 Flow variables update

For both CGKS-5th and GKS-2nd, the cell-averaged conservative variable $\bm{Q}(t)$ within a control volume $\Omega$ is updated according to the conservation laws, as given in Eq.(2). In CGKS-5th, additionally, the cell-averaged and line-averaged derivatives of $\bm{Q}(t)$ are also updated.

Using Gauss’s theorem, the cell-averaged gradient of $\bm{Q}(t)$ over $\Omega$ can be expressed as

[TABLE]

where $\bm{n}$ is the unit normal vector on the cell interface. The line-averaged partial derivative along a segment defined by points $\bm{x}_{1}$ and $\bm{x}_{2}$ is given by

[TABLE]

where $\left|l_{G}\right|=\left|\bm{x}{2}-\bm{x}{1}\right|$, and $\bm{Q}(\bm{x}_{1},t)$ and $\bm{Q}(\bm{x}_{2},t)$ are the values of $\bm{Q}$ at $\bm{x}_{1}$ and $\bm{x}_{2}$ inside the cell, respectively. At each time step, $\bm{Q}(\bm{x}_{1},t)$ and $\bm{Q}(\bm{x}_{2},t)$ are obtained from the time-accurate solution of the gas distribution function in the GKS. This solution is also used to compute the numerical flux at the interfaces, and it is specified using interface values provided by spatial reconstruction.

3.2 Fifth-order compact gas-kinetic scheme

Figure 1(a) and Figure 1(b) illustrate the schematic of the compact stencils used for spatial reconstruction in CGKS-5th and GKS-2nd, respectively. In these figures, the red cubes denote face-neighboring cells, while the blue cubes indicate edge-neighboring cells. For the target cell $\Omega_{0}$, its interfaces are denoted by $F_{m}$ ($m=1,\dots,6$), and its edges are denoted by $E_{n}$ ($n=1,\dots,12$). The cell adjacent to $\Omega_{0}$ that shares the face $F_{m}$ is represented by $\Omega_{1,m}$, while the cell adjacent to $\Omega_{0}$ that shares only the edge $E_{n}$ is represented by $\Omega_{2,n}$.

For CGKS-5th, to achieve fifth-order accuracy, a compact stencil containing all the face-neighboring and edge-neighboring cells is defined as follows

[TABLE]

The following data are utilized for fifth-order reconstruction

[TABLE]

where

[TABLE]

$\bm{Q}_{0}$, $\bm{Q}_{1,m}$ and $\bm{Q}_{2,n}$ represent the cell-averaged conservative values over cells $\Omega_{0}$, $\Omega_{1,m}$ and $\Omega_{2,n}$, respectively. $(\partial_{l}\bm{Q})_{0,G}$ represents the line-averaged partial derivatives within the target cell $\Omega_{0}$. $(\nabla\bm{Q})_{1,m}$ denotes the cell-averaged gradient over cell $\Omega_{1,m}$, where $\nabla=(\partial_{x},\partial_{y},\partial_{z})$. $(\partial_{\bm{\tau}}\bm{Q})_{2,n}$ denotes the cell-averaged directional derivative over the cell $\Omega_{2,n}$, with $\partial_{\bm{\tau}}=(\partial_{\tau_{1}},\partial_{\tau_{2}})$, where $\partial_{\tau_{1}}$ and $\partial_{\tau_{2}}$ are linear combinations of $\partial_{x}$, $\partial_{y}$, and $\partial_{z}$. Further details on the computation of these quantities can be found in [17].

For the target cell $\Omega_{0}$, a quartic polynomial $P^{4}(\bm{x})$ can be reconstructed based on $S^{data-5th}$, defined by

[TABLE]

where $\bm{Q}_{0}$ is the cell-averaged variable over cell $\Omega_{0}$, the multi-index $\bm{d}=(d_{1},d_{2},d_{3})$, and $|\bm{d}|=d_{1}+d_{2}+d_{3}$. The zero-mean basis $p_{\bm{d}}(\bm{x})$ is defined as

[TABLE]

where

[TABLE]

$h_{x}$, $h_{y}$ and $h_{z}$ are the characteristic scales of $\Omega_{0}$ along the three directions of axes, and $(x_{0},y_{0},z_{0})$ is the centroid of $\Omega_{0}$. To determine the polynomial $P^{4}(\bm{x})$, the following system is solved using the constrained least-squares method

[TABLE]

where the parameters $h_{1}$, $h_{2}$, and $h_{3}$ are used to reduce the condition number of the reconstruction matrix, the choice of which can be found in [17]. In the constrained least-squares system, the first-row equation in Eq.(9) is satisfied exactly, while the remaining equations are satisfied in the least-squares sense.

To handle discontinuities, nonlinear reconstruction is necessary. For CGKS-5th, we adopt the generalized ENO (GENO) nonlinear reconstruction method [32] to obtain nonlinear reconstruction. At discontinuities, the GENO reconstruction is dominated by a second-order reconstruction with ENO property. The second-order reconstruction is determined based on sub-stencils. The sub-stencils $S_{m}^{cell-2}$ ($m=1,\dots,6$) of the CGKS-5th scheme are selected as

[TABLE]

The construction of several linear polynomials on these sub-stencils is based on the following two types of data

[TABLE]

where $\bm{Q}_{1,m}$ is the cell-averaged conservative value over cell $\Omega_{1,m}$, $(\partial_{l}\bm{Q})^{m}_{0,G}$ is the selected line-averaged partial derivative within the target cell $\Omega_{0}$. Detailed descriptions of this selection can be found in [17]. Using $S_{m}^{data-2}$ ($m=1,\dots,6$), linear polynomials $P_{m}^{1}(\bm{x})$ can be constructed of the form

[TABLE]

where the basis functions $p_{\bm{d}}(\bm{x})$ are defined as in Eq.(8). To determine these linear polynomials, the following constrains need to be satisfied

[TABLE]

where $\bm{Q}_{1,m}$ is the cell-averaged variable over cell $\Omega_{1,m}$. The resulting linear system is solved using a constrained least-squares method.

To accommodate non-uniform meshes while keeping memory usage low and the implementation simple, the CGKS-5th scheme adopts a simple transformation strategy [17]. Specifically, a coordinate mapping sends each target cell to a single reference cell, allowing the reconstruction to be represented in a unified polynomial form. Consequently, only one reconstructed polynomial on the reference cell needs to be stored, reconstructed values are then mapped back from the reference coordinates to the physical coordinates.

Using reconstructions on the compact stencil and sub-stencils, GENO is applied to achieve nonlinear reconstruction. It adaptively transitions from high-order linear reconstruction in smooth regions to second-order reconstruction near discontinuities, effectively balancing accuracy and robustness. With the reconstructed polynomial $P^{4}(\bm{x})$ and $P_{m}^{1}(\bm{x}),m=1,...,6$, the point-value $\bm{Q}(\bm{x}_{G})$ and the spatial derivatives $\partial_{x,y,z}\bm{Q}(\bm{x}_{G})$ at Gaussian quadrature point for the nonlinear reconstruction are computed as

[TABLE]

where

[TABLE]

with linear weights $d_{1}=\dots=d_{6}=\displaystyle\frac{1}{6}$. The small positive constant is set to $\epsilon=10^{-15}$. Details of the computation of $\chi$ can be found in [32, 17].

3.3 Second-order gas-kinetic scheme

For consistency, we retain the notation used in Section 3.2. For GKS-2nd, the reconstruction stencil for the target cell $\Omega_{0}$ is defined as

[TABLE]

which involves only face-neighboring cells, as shown in Figure 1(b). The second-order reconstruction can be achieved based on

[TABLE]

A linear polynomial $P^{1}(\bm{x})$ is constructed using $S^{data-2nd}$ and is defined as in Eq.(10). To determine this second-order polynomial, the following system is solved using the least-square method

[TABLE]

To handle discontinuities, the linear polynomial $P^{1}(\bm{x})$ is coupled with the Venkatakrishnan limiter [33] to enable a nonlinear reconstruction as follows

[TABLE]

The limiter value $\phi_{0}$ for the target cell $\Omega_{0}$ is taken as

[TABLE]

where $\phi_{F_{m}}$ denotes the Venkatakrishnan limiter [33] evaluated on face $F_{m}$.

3.4 Temporal discretization

According to Eq.(7), the discretization of the conservation laws over cell $\Omega_{i}$ can be represented as

[TABLE]

where $\bm{Q}_{i}^{n+1}$ represents the cell-averaged conservative variable over cell $\Omega_{i}$ at $t_{n+1}=t_{n}+\Delta t$. For GKS-2nd, a second-order temporal discretization is used

[TABLE]

For the CGKS-5th scheme, high-order temporal accuracy is achieved with a two-stage, fourth-order method [21, 22, 34]

[TABLE]

The cell-averaged gradients and line-averaged derivatives are evaluated simultaneously using conservative variables at cell interfaces. The evolution of the conservative variable on one side of an interface follows a two-stage update

[TABLE]

To provide conservative variables on both sides of an interface, which may differ near discontinuities, the update model for $\bm{Q}(t)$ is

[TABLE]

Further details are provided in [31, 15].

4 Numerical verification

In this section, the performance of CGKS-5th and GKS-2nd is evaluated using two approaches. The first approach quantifies the computational resources required to achieve a prescribed target resolution. The second approach evaluates the accuracy and resolution of the schemes under equivalent computational resources. Using these criteria, numerical tests are conducted ranging from subsonic to supersonic turbulent flows, thereby assessing performance across smooth and discontinuous solutions. For clarity, CGKS-5th and GKS-2nd with linear reconstructions are specifically referred to as linear CGKS-5th and linear GKS-2nd, respectively.

For evolution of flow fields, the time step $\Delta t$ is given by the CFL condition. For viscous flows, the time step is additionally constrained by the viscous term as $\Delta t=\text{CFL}\cdot h^{2}/(3\nu)$, where $h$ is the cell size and $\nu$ is the kinematic viscosity coefficient. For the inviscid flows, the collision time $\tau$ takes

[TABLE]

where $p_{l}$ and $p_{r}$ denote the pressure on the left and right sides of the cell interface, $\epsilon=0.05$ and $C=10.0$. For the viscous flows, the collision time is related to the viscosity coefficient,

[TABLE]

where $\mu$ is the dynamic viscous coefficient and $p$ is the pressure at the cell interface. In smooth flow regions, the collision time reduces to

[TABLE]

In the numerical cases presented in this section, CGKS-5th and GKS-2nd are implemented on multiple GPUs using the Compute Unified Device Architecture (CUDA) for parallel computations and the Message Passing Interface (MPI) for inter-process communication. This configuration enables large-scale simulations and is executed on NVIDIA GeForce RTX 4090 GPUs with double-precision.

4.1 Subsonic Taylor-Green vortex flow

The Taylor-Green vortex is a well-established benchmark in fluid mechanics, widely used to assess numerical methods in terms of turbulence generation, energy dissipation, and the evolution of vortex structures [35, 36, 37]. In this subsection, the subsonic Taylor-Green vortex is selected as a representative benchmark case, enabling a detailed performance evaluation of the linear CGKS-5th and linear GKS-2nd. In this case, the comparison focuses on analyzing the computational cost required by each scheme to achieve a prescribed target resolution.

The initial velocity field is given by

[TABLE]

and the initial pressure field is set as

[TABLE]

In the computation, $U_{0}=1$, $\rho_{0}=1$. The Reynolds number is $Re=1600$ defined by $Re=\rho_{0}U_{0}L/\mu_{0}$ , where $\mu_{0}$ is the dynamic viscosity coefficient. The Mach number is defined as $Ma_{\infty}=U_{0}/\sqrt{\gamma p_{0}/\rho_{0}}$, and this case is conducted at $Ma_{\infty}=0.1$, which lies within the incompressible limit. The Taylor-Green vortex flow is defined in a periodic domain defined as $[-\pi L,\pi L]^{3}$, $L=1$. A uniform mesh with $96^{3}$ cells is utilized for the linear CGKS-5th, while a uniform mesh with $296^{3}$ cells is employed for the linear GKS-2nd. It has been tested that the above grid configurations ensure comparable resolution for the two schemes.

Figure 2 presents the iso-surfaces of the Q-criterion ($Q=1.6$) at $t=10$ and $t=20$, computed using the linear CGKS-5th on a $96^{3}$ mesh, with the iso-surfaces colored by the Mach number. It can be observed that CGKS-5th is capable of capturing flow structures on a coarse mesh. In order to quantitatively evaluate the performance of CGKS-5th and GKS-2nd, several time-averaged quantities are calculated. The volume-averaged kinetic energy is given by

[TABLE]

where $\Omega$ is the computational domain and $E_{k}$ is calculated by numerical quadrature. The integrated enstrophy is defined as

[TABLE]

where the dynamic viscosity $\mu$ is computed by the Sutherland’s law, $\epsilon^{S}$ is calculated by numerical quadrature, and the velocity derivative values at quadrature points are obtained by the compact reconstruction. The comparison of the time histories of kinetic energy $E_{k}$ and enstrophy $\epsilon^{S}$ using the linear CGKS-5th on a $96^{3}$ mesh and the linear GKS-2nd scheme on a $296^{3}$ mesh is presented in Figure 3. The reference line in Figure 3 corresponds to a spectral solution computed on a $512^{3}$ mesh [38]. As observed in Figure 3, the two numerical results achieve equivalent resolution.

Table 1 summarizes the computation time of the linear CGKS-5th on a $96^{3}$ mesh and the linear GKS-2nd on a $296^{3}$ mesh, both using multi-GPU parallel computing. As shown in Table 1, CGKS-5th achieves nearly an order of magnitude improvement in computational efficiency compared to GKS-2nd for equivalent resolution. This demonstrates that, compared to the second-order GKS, the fifth-order CGKS can achieve higher resolution with significantly lower computational costs.

4.2 High-speed Taylor-Green vortex flow

In previous studies, the subsonic and transonic Taylor-Green vortex flows have been extensively analyzed [37, 39, 40]. However, the high-speed Taylor-Green vortex flow (Mach number $Ma_{\infty}\geq 2$) has not been widely adopted as a benchmark problem. Under such supersonic conditions, the original fluctuation form presented in Eq.(13) cannot be directly extended to the high-speed regime. To address this limitation, we propose a new initial condition, inspired by the classical form, that is tailored to high-speed flows. In this subsection, the supersonic Taylor-Green vortex flow with $Ma_{\infty}=2$ is considered, which exhibits pronounced shock-vortex interactions. CGKS-5th and GKS-2nd employ nonlinear combination and Venkatakrishnan limiter to handle discontinuities, respectively, whereas the nonlinear part of CGKS-5th being algorithmically more complicated. The objective is to evaluate whether the nonlinear CGKS-5th can still deliver superior computational efficiency under these more demanding conditions.

The supersonic Taylor-Green vortex flow is defined in a periodic domain defined as $[-\pi L,\pi L]^{3},L=1$. The initial condition is set as

[TABLE]

and

[TABLE]

In the computation, $\rho_{0}=1,\,p_{0}=1,\,U_{0}=\sqrt{\gamma}Ma_{\infty},\,Ma_{\infty}=2.0,\,T_{0}=1$, the Reynolds number $Re=1600$, and the Prandtl number is set as 0.7. The dynamic viscosity $\mu$ is computed through the Sutherland law

[TABLE]

where

[TABLE]

The kinetic energy integrated over the domain is given by

[TABLE]

where $\Omega$ is the computational domain and $E_{k}$ is calculated by numerical quadrature. The total viscous dissipation rate is defined as

[TABLE]

where $\epsilon^{S}$ is the solenoidal (enstrophy) dissipation, and $\epsilon^{D}$ is the dilatational dissipation, $\epsilon_{T},\epsilon^{S},\epsilon_{D}$ are calculated by numerical quadrature, and the velocity derivative values at quadrature points are obtained by the compact reconstruction.

A uniform mesh with $116^{3}$ cells is utilized for CGKS-5th, while a uniform mesh with $384^{3}$ cells is employed for GKS-2nd. Figure 5 presents the iso-surfaces of the Q-criterion ($Q=4$) at $t=5$ and $t=10$, obtained using CGKS-5th on the $116^{3}$ mesh. The iso-surfaces are colored by the Mach number, illustrating the flow structures. Figure 5 shows the density gradient magnitude on the centerline (x-y) plane at $z=0$ for $t=5$ and $t=10$, computed using CGKS-5th on a $116^{3}$ mesh. These results demonstrate the vortex decay over time and highlight the capability of CGKS-5th to accurately resolve flow structures with high resolution. The comparison of the time histories of kinetic energy $E_{k}$ and enstrophy dissipation $\epsilon^{S}$ obtained using CGKS-5th on a $116^{3}$ mesh and GKS-2nd on a $384^{3}$ mesh is presented in Figure 6. Additionally, the results from CGKS-5th on a $256^{3}$ mesh are used as the reference solution. As shown in Figure 6, the numerical results computed with CGKS-5th on a $116^{3}$ mesh achieve the same resolution as those computed with GKS-2nd on a $384^{3}$ mesh. Table 2 compares the computational time required for these simulations using multi-GPU parallelization. Notably, despite the higher algorithmic complexity of CGKS-5th caused by its nonlinear reconstruction, it still demonstrates a significant efficiency advantage over the limiter-based GKS-2nd, achieving an improvement close to an order of magnitude.

4.3 Transonic and Supersonic Jet Simulation

Jet simulation plays a critical role in advancing the understanding of jet flows and optimizing them for a wide range of engineering applications. In transonic and supersonic jet simulations, computation is particularly challenging due to the involvement of complex physical phenomena, such as the interaction of shear layers and vortex structures, as well as the transition from laminar to turbulent flow. To accurately and efficiently capture these intricate dynamics, a combination of high-order numerical schemes and large-scale computations is typically required to perform long-duration, scale-resolving simulations. As discussed in the previous subsection, achieving comparable resolution with a second-order scheme may require tens of times more cells than a fifth-order scheme, making the simulation prohibitively expensive at that scale. In this subsection, the accuracy and resolution of CGKS-5th and GKS-2nd are compared under the constraint of equal computational resources. Specifically, the performance of the schemes is evaluated at matched computational cost by employing different mesh sizes, with mesh distributions carefully designed to ensure a fair comparison.

For this case, the isothermal jet with Mach numbers $Ma_{\infty}=0.9$ and $Ma_{\infty}=1.5$ are considered, with the Reynolds number set as $Re=4\times 10^{5}$ based on a jet diameter of $D=2r_{0}=1$. None of the simulations include any nozzle geometry. The inflow axial velocity of the jet is initialized through a hyperbolic-tangent profile similarly to previous studies [41, 42]

[TABLE]

where $U_{c}=Ma_{\infty}\cdot c_{\infty}$ is the inflow centerline velocity, $c_{\infty}=1$, and the initial momentum thickness $\delta_{\theta}(0)$ is set as $0.03D$. The inlet temperature profile is calculated using the Crocco-Busemann relation, assuming constant pressure

[TABLE]

where $T_{c}$ and $T_{\infty}$ are the jet centerline and ambient temperatures, respectively, in the present isothermal case, $T_{c}=T_{\infty}$ and $\sigma={Pr}^{1/3}$ is the recovery factor, the Prandtl number $Pr$ is set to 0.7. Outside the inlet, the initial condition is given as follows,

[TABLE]

No disturbance is added to the field. The non-reflecting boundary condition and sponge zone are implemented in the outflow directions to filter out possible reflected waves.

In this case, a structured mesh consisting of $352\times 176\times 176$ cells is designed for CGKS-5th, as illustrated in Figure 7. In the flow direction, the mesh size is constant at $\Delta x=0.04D$ up to $x=9.6D$ to resolve the domain near the nozzle exit. Between $9.6D<x<15D$, the mesh size increases exponentially from $0.04D$ to $0.16D$. Beyond this, a stretched mesh forms a sponge zone with $\Delta x_{max}=1.6D$ to damp vortices and perturbations. Radially, the mesh size is $0.025D$ at the jet center, decreasing exponentially to $0.015D$ at $r=r_{0}$ to resolve small-scale structures near the jet boundary. Outside the jet, the mesh size increases back to $0.025D$ and gradually expands to $0.2D$ within $r<4.5D$. Beyond this, a sponge region is created with a maximum mesh size of $0.4D$. To achieve computational cost equivalence with the fifth-order scheme using the above mesh, while maintaining a consistent stretching design, a structured mesh consisting of $528\times 264\times 264$ cells is designed for GKS-2nd, where the size of each cell is approximately $1/1.5$ times the cell size of the original mesh at the corresponding locations. Table 3 presents the computational time required for CGKS-5th on a mesh consisting of $352\times 176\times 176$ cells and GKS-2nd on a mesh consisting of $528\times 264\times 264$ cells. It can be observed that the computational costs of the two schemes are approximately equivalent, with a cost difference within 5$\%$.

For the transonic jet simulation with $Ma_{\infty}=0.9$, Figures 9 and 9 provide instantaneous visualizations of the iso-surfaces of the Q-criterion, which are colored by the Mach number, alongside grayscale pressure maps illustrating pressure fluctuations relative to $P_{0}$. The results are obtained with CGKS-5th and GKS-2nd, respectively. The higher-order CGKS-5th clearly resolves more intricate turbulent structures than the second-order GKS-2nd. Figure 11 and Figure 11 present the mean axial velocity and root mean square (RMS) axial velocity along the centerline ($\bm{r}=0$). These results are computed using CGKS-5th on a mesh with $352\times 176\times 176$ cells and GKS-2nd on a mesh with $528\times 264\times 264$ cells. Theses figures compare the numerical results with experimental results from Arakeri et al.[44] and Power et al.[43], and with results obtained using the LES method combined with a sixth-order compact scheme on a multiblock structured grid containing 24.5 million cells (Fosso Pouangué et al.) [42]. Furthermore, in Figure 11, the position is normalized by the corresponding potential core length, $L_{C}$, defined as the distance from the jet inlet to the point where the centerline velocity drops to 95$\%$ of the inflow centerline velocity. The potential core length values obtained from CGKS-5th, GKS-2nd, and the reference data are presented in Table 4. As shown, $L_{C}$ computed by CGKS-5th agrees more closely with the reference data than that from GKS-2nd. As shown in Figure 11 and Figure 11, the mean axial velocity computed using CGKS-5th demonstrates excellent agreement with both the experimental data [44, 43] and the numerical results obtained from the sixth-order LES on a refined mesh [42]. In contrast, GKS-2nd exhibits significant discrepancies. Similarly, the RMS axial velocity predicted by CGKS-5th also aligns well with the LES results [42] and the experimental data from Arakeri et al. [44], whereas the second-order scheme displays clear deviations. Notably, in the region $x\in[4,10]$, where the turbulent structures are most complex, the results computed by GKS-2nd show pronounced discrepancies in both Figure 11 and Figure 11. These results demonstrate that, at matched computational cost, CGKS-5th delivers higher accuracy and resolution, yielding more reliable turbulence predictions.

For the supersonic jet simulation with $Ma_{\infty}=1.5$, Figure 13 and Figure 13 show instantaneous visualizations of the iso-surfaces of the Q-criterion, colored by the Mach number, along with grayscale pressure visualizations representing pressure fluctuations around $P_{0}$, obtained using CGKS-5th and GKS-2nd, respectively. The higher-order CGKS-5th again demonstrates superior capability in resolving flow-field details. Figure 14 presents the instantaneous shadowgraph contours obtained using CGKS-5th and GKS-2nd. Both results show the propagation of sound waves within the computational domain. Notably, CGKS-5th results reveal clearer small-scale structures and sharper internal shock features within the jet, evidencing the higher effective resolution of the high-order method. Figure 16 and Figure 16 compare the mean axial velocity and root-mean-square (RMS) axial velocity obtained using CGKS-5th and GKS-2nd. The oscillations in the mean axial velocity computed by CGKS-5th over $x\in[0,5]$ align with the locations of internal shock waves identified in Figure 14, whereas the second-order GKS-2nd fails to capture this signature. These results demonstrate that, at matched computational cost, CGKS-5th achieves substantially higher resolution for complex turbulent flows featuring shocks and small-scale vortical structures, underscoring its practical value for turbulence simulations.

5 Conclusion

This study presents a performance comparison between CGKS-5th and GKS-2nd on structured meshes for complex viscous flows, with a focus on computational efficiency and simulation accuracy. Two comparative approaches are employed: (i) the computational cost required to achieve comparable resolution, and (ii) the accuracy and resolution attainable under matched computational resources. Across both smooth and discontinuous flows, CGKS-5th achieves the same resolution at approximately a 7-9 times lower computational cost than GKS-2nd. Under equal resource constraints, CGKS-5th further delivers substantially higher accuracy and resolution, with particularly pronounced benefits for turbulent flows involving shocks and small-scale vortical structures. The advantages of high-order schemes for smooth compressible flows have been well-recognized in previous studies. However, their effectiveness in complex compressible flow simulations, where nonlinear reconstructions are introduced to handle discontinuities and may reduce the benefits of high-order schemes, has lacked thorough validation. For the first time, this study demonstrates, through simulations of turbulence benchmarks and jet flow problems, that high-order schemes, particularly compact schemes with high-resolution properties, offer significant advantages in both accuracy and efficiency over second-order schemes, underscoring their potential value for engineering applications.

Acknowledgements

The current research is supported by National Key R$\&$D Program of China (Grant Nos.2022YFA1004500), National Science Foundation of China (92371107, 12172316), and Hong Kong research grant council (16301222, 16208324).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.