High order semi-implicit weighted compact nonlinear scheme for the all-Mach isentropic Euler system

The computation of compressible flows at all Mach numbers is a very challenging problem. An efficient numerical method for solving this problem needs to have shock-capturing capability in the high Mach number regime, while it can deal with stiffness and accuracy in the low Mach number regime. This paper designs a high order semi-implicit weighted compact nonlinear scheme (WCNS) for the all-Mach isentropic Euler system of compressible gas dynamics. To avoid severe Courant-Friedrichs-Levy (CFL) restrictions for low Mach flows, the nonlinear fluxes in the Euler equations are split into stiff and non-stiff components. A third-order implicit-explicit (IMEX) method is used for the time discretization of the split components and a fifth-order WCNS is used for the spatial discretization of flux derivatives. The high order IMEX method is asymptotic preserving and asymptotically accurate in the zero Mach number limit. One- and two-dimensional numerical examples in both compressible and incompressible regimes are given to demonstrate the advantages of the designed IMEX WCNS.


Introduction
Low-Mach number flows in fluid dynamics that are slow compared with the speed of sound can be described by the following scaled isentropic Euler equations.
where ρ(x, t) > 0, u(x, t) = (u(x, t), v(x, t)) and q = ρu(x, t) are the density, velocity and momentum of the fluid, respectively, and ε > 0 is the scaled Mach number that is a measure of compressibility of the fluid. For example, ε = 1 corresponds to the fully compressible regime. The operators ∇, ∇· and ⊗ are the gradient, divergence and tensor product operators, respectively. The pressure p in Eq. (1) is given by the equation of state where > 0 and γ ≥ 1 are constants. Equation (1) are hyperbolic with wave speeds in direction n λ 1,3 = u · n ± c, λ 2 = u · n, ( 3 ) and the sound speed c = p (ρ)/ε. Klainerman and Majda [1,2] proved that the solution of Eq. (1) converges to the incompressible counterpart when ε ∼ 0. The asymptotic nature of the solution to Eq. (1) can be studied by plugging-in the ansatz for variables ρ, p and u. After performing a scale analysis, we obtain an incompressible Euler system for the velocity based on well-prepared initial data, i.e., and appropriate boundary conditions [3], e.g., periodic, wall or homogeneous Neumann boundary conditions. Taking the divergence of the first equation in (5) and using the incompressibility, we obtain In low Mach number flows, the challenge faced by standard hyperbolic numerical methods for Eq. (1) is that the system becomes stiff. The stiffness is caused by the stringent Courant-Friedrichs-Levy (CFL) restriction [3][4][5][6], i.e. t ≤ x/max i |λ i | = O(ε x) with the spatial and temporal grid sizes x and t. The restriction results in an increasingly large computational time for ε → 0. Moreover, the strong stability requirement for ε → 0 results in the excessive numerical viscosity of the standard hyperbolic solvers that is proportional to max i |λ i | = O(ε −1 ). This makes the standard hyperbolic solvers suffer from severe loss of accuracy due to the creation of spurious waves. Therefore, the stability and accuracy of the standard hyperbolic numerical methods highly depend on ε [4].
In the literature, some all-speed schemes [3][4][5][6][7][8][9][10][11][12][13][14][15][16] have been derived for the compressible-incompressible limit problems based on the idea of the asymptoticpreserving (AP) methodology introduced by Jin [17,18] for relaxation systems. The AP all-speed schemes are uniformly consistent with the low-Mach number limit and are uniformly stable. Although a robust AP scheme should allow both explicit and implicit time-marching techniques, the explicit strategies suffer from the CFL stability conditions and thus affect the efficiency of the method for small Mach numbers. The implicit strategies avoid the CFL restriction, but they need to solve highly nonlinear equations and add an excessive numerical dissipation to the slow wave, degrading the accuracy of the method [3]. To overcome the drawbacks, the semi-implicit or implicit-explicit (IMEX) time discretization strategies are widely applied based on the acoustic/advection splitting of the Euler equations into the stiff/non-stiff parts (or fast and slow fluxes) [3][4][5]10]. These strategies treat with the stiff part implicitly and the non-stiff part explicitly to solve the isentropic Euler equations effectively and efficiently. The IMEX time discretization combined with the shock-capturing schemes, e.g., the first-order local Lax-Friedrichs (or Rusanov) scheme [4], the second-order finite volume method [3,19], the arbitrary high order accurate discontinuous Galerkin method [20] and the fifth-order finite difference weighted essentially non-oscillatory (WENO) scheme [10], has shown its good performances in compressible and incompressible regimes.
Recently a class of weighted compact nonlinear schemes (WCNS's) as high order accurate shock-capturing schemes have been widely adopted to investigate compressible flows with shock waves and discontinuities [21][22][23][24][25][26][27]. The WCNS's are first proposed by Deng and Zhang [21] based on the cell-centered compact scheme of Lele [28] and are developed in [22][23][24]29]. The expression of the classical cell-centered compact scheme is that a linear combination of the unknown approximations of the first derivative of a function at the grid points in a stencil is equal to another linear combination of the function itself at the grid points in the same stencil. Compared to noncompact schemes, e.g., the classical WENO scheme [30], the stencil can be more compact in the compact scheme for the same order of accuracy [22]. WCNS's that integrate WENO interpolation into the cell-centered compact scheme for the first derivative of a function are variants of the WENO schemes and have the good performance in terms of accuracy, robustness and efficiency for compressible flows. Compared to the classical WENO schemes, the WCNS's feature a better spectral resolution at the same order of accuracy [21,22]. They are more flexible in the choice of flux splitting methods, e.g., Roe's flux difference splitting method, van Leer's flux vector splitting and Liou's advection upstream splitting method, and flux derivatives can be computed by different explicit or implicit compact finite difference methods [23,29]. In this work, in view of the good performance of the IMEX time discretization method and the WCNS's, we design a high order IMEX WCNS for the all-Mach isentropic Euler system based on the acoustic/advection splitting strategy [10]. The stiff part governing fast acoustic and the non-stiff part modeling slow nonlinear advection effects are solved by the high order implicit and explicit time discretization methods, respectively. In the zero Mach number limit, the high order IMEX scheme has the AP property and is asymptotically stable under a convective CFL condition independent of the Mach number ε. Oneand two-dimensional numerical examples are given to demonstrate good performance of the designed IMEX WCNS in both compressible and incompressible regimes.

WCNS method
In this section, a fifth-order WCNS is considered for the spatial discretization of flux derivatives. The WCNS procedure consists of three components, i.e., midpoint-to-node differencing, flux evaluation at the midpoint and node-to-midpoint weighted averaging interpolation of fluxes. Here, we briefly describe the framework of WCNS discretization for the first derivative, f x , in the one-dimensional scalar case u t +f (u) x = 0. For simplicity, we take a uniform grid defined by x i = i x = ih, i = 0, ..., N, where h is the uniform grid spacing and N is the number of cells.
The classical cell-centered compact scheme for the first derivatives [28] has no dissipation, but the dispersive errors dominate the scheme, which may destroy the numerical solutions. A class of hybrid cell-node and cell-edge WCNS's proposed in [29] overcome the disadvantages and are written as which introduce dissipation by approximating the cell-edge fluxes f i±(m−1/2) based on nonlinear weighted averaging interpolation of fluxes. The coefficients k , η l and ζ m in Eq. (8) can be acquired by matching the Taylor series coefficients of various orders. If the coefficients k = 0, a class of explicit WCNS's can be derived. The explicit WCNS's are preferred because of their low cost and simplicity for parallelization and vectorization, compared to the implicit WCNS's with a certain coefficient k * = 0. Here, we consider one of explicit WCNS's and use a sixth-order explicit central-type scheme obtained from Eq. (8) to approximate f x , i.e., By the Taylor series expansion of (9), we obtain where the truncation term O(h r ) comes from the numerical approximations of the celledge fluxes f i±1/2 , denoted byf i±1/2 . The accuracy of the scheme (9) is min (6, r). Here, the numerical fluxesf i±1/2 approximate f i±1/2 to the fifth order and thus a fifth-order WCNS is designed for the spatial discretization.
Various flux splitting methods [23] can be used in WCNS to introduce correct upwinding. For simplicity, we use the Lax-Friedrichs flux splitting, i.e., the flux is split into positive and negative parts: where α = max u |f u | is the maximum characteristic wave speed. The numerical fluxeŝ f i±1/2 at the cell centers in (9) are calculated aŝ A fifth-order WENO interpolation technique is considered to obtain the positive numerical fluxf + i+1/2 . The interpolation process of the negative numerical fluxf − i+1/2 is symmetric with respect to x i+1/2 and will not be shown here. For simplicity, the "+" sign in the superscript is removed.
The numerical fluxf i+1/2 on a cell center can be interpolated locally from the values on a five-point upwind-biased stencil S = {x i±2 , x i±1 , x i }. The stencil S can be divided into three sub-stencils S k = {x i+k−2 , x i+k−1 , x i+k }, k = 0, 1, 2. On each sub-stencil a thirdorder linear interpolation of f at cell center x i+1/2 can be derived, i.e., wheref k i+1/2 , k = 0, 1, 2, are the approximated values at cell center x i+1/2 on S k , k = 0, 1, 2, respectively. On the stencil S a fifth-order upwind-biased linear interpolation of f at cell center x i+1/2 can be derived, i.e., which can also be constructed linearly from the three third-order interpolationsf k i+1/2 : with the linear weights The central-type scheme (9) with the upwind-biased linear interpolation (16) or (17) will generate oscillations around shock waves and discontinuities. To capture shock waves and discontinuities, the numerical fluxf i+1/2 should be computed by using the nonlinear combination of upwind-biased interpolations, i.e., where ω k , k = 0, 1, 2, are the upwind-biased nonlinear weights that are associated with the relative smoothness of f on each sub-stencil. In smooth areas the nonlinear weights ω k should be close to the linear weights d k to achieve the designed order and be close to zero near discontinuities to suppress spurious numerical oscillations near discontinuities.
There are several types of nonlinear weights ω k in the literature, e.g., the WENO-JS [30], WENO-M [31] and WENO-Z [32][33][34] weights. The WENO-Z weights place a larger weight to the less smooth sub-stencil, and thus they are less dissipative than the WENO-JS weights and more efficient than the WENO-M weights. Here, the WENO-Z weights ω k proposed by Borges et al. [32] are chosen.
where the second-order smoothness indicators β k , k = 0, 1, 2, derived by Hu [35] are adopted to measure the smoothness of f on the sub-stencils S k . with In (20), the global higher order smoothness indicator τ = |β 0 − β 2 | and the parameter δ is a small constant to avoid denominators becoming zero.

Remark 1
The value of δ in (20) will influence the convergence, accuracy and stability of the schemes. Here, a parameter-free algorithm proposed by Zheng et al. [36] is used to calculate the value of δ as follows.
where β ave = 1 To avoid denominators becoming zero in (20) and (22), β k + δ in (20) and β std in (22) are replaced with max(β k + δ, err) and max(β std , err), respectively. Here, err is the square root of the smallest positive number allowed for a machine.

IMEX method
In this section, we briefly describe the framework of IMEX time discretization designed for the numerical integration of the stiff ordinary differential equation (ODE) as follows.
wheref andḡ are the non-stiff and stiff parts of the ODE, respectively. The IMEX time discretization treats the non-stiff part explicitly and the stiff part implicitly in order to overcome the stringent time-step restriction for stability. There are many IMEX methods designed for various problems in the literature [37][38][39][40][41]. Among them, the IMEX Runge-Kutta (IMEX-RK) schemes as combinations of diagonally implicit A-or L-stable RK (DIRK) methods and explicit RK methods offer a precise and robust approach to define high order AP schemes. The computational cost per time step is proportional to the number of stages in the IMEX-RK schemes and the accuracy is proportional to the time-step size and the order of accuracy.
The IMEX-RKp(s, σ , p) scheme is characterized by the s × s matrices A = (a ij ),Ã = (ã ij ), the vectors b,b ∈ R s , and the vectors c,c ∈ R s . This scheme with s stages of the implicit scheme for the stiff term and σ (= s + 1) stages of the explicit scheme for the nonstiff term has p th -order temporal accuracy. Applying this scheme to solve Eq. (23) from t n to t n+1 yields

Definition 1
The IMEX-RKp(s, σ , p) scheme is of type ARS, i.e., the matrices A,Ã for the implicit and explicit schemes are described respectively as with the submatrices

IMEX-ARSp wCNS
In this section, we design the IMEX-ARSp WCNS scheme for the isentropic Euler system. Here, the fifth-order WCNS described in the previous section is used for the spatial discretization and the p th -order IMEX-ARSp method is used for the time dis-cretization. We first split the Euler system into two parts based on the acoustic/advection splitting strategy [10] and obtain where U = (ρ, q), F = −(0, ∇ · q⊗q ρ ) is the non-stiff nonlinear term and G = −(∇ · q, 1 ε 2 ∇p) is the stiff linear term.
The p th -order IMEX-ARSp scheme that consists of (24) and (25) is applied to each component of (28). The intermediate stages are described as where k = 1, 2, ..., s, and the numerical solutions at time t n+1 are described as The intermediate stages (29) and (30) consist of two implicit steps and can be simplified by eliminating q (k) between (29) and (30) [10], and obtain a nonlinear elliptic equation for ρ (k) : Based on (30) and (35), q (k) can be calculated as

Remark 3
To update the solutions from t n to t n+1 by the p th -order IMEX-ARSp WCNS, we first solve the elliptic equation (33) to obtain ρ (k) and then calculate explicitly q (k) using (36). Here, (33) is a nonlinear equation for ρ (k) , which is solved by inexact Newtongeneralized minimum residual (Newton-GMRES) algorithms. Refer to [42,43] for the detailed description of the algorithms. Finally, the updates ρ n+1 and q n+1 are calculated explicitly using (31) and (32) with the intermediate values ρ (k) and q (k) . For the GSA IMEX-ARSp scheme, ρ n+1 = ρ (s) and q n+1 = q (s) .

Remark 4
The spatial discretization of both convective terms in (33)- (35) and the gradient terms in (34) and (35) are solved by the fifth-order WCNS in the component-bycomponent and dimension-by-dimension forms. Here, in the interpolation process of the cell-edge pressures, the Lax-Friedrichs flux with zero numerical viscosity is considered, i.e., p = 1 2 p + 1 2 p, otherwise the numerical diffusion will become large when ε → 0 [10].
Remark 5 For the explicit time discretization, the CFL stability condition is described as which leads to unaffordable computational costs in the low Mach number regime. For the IMEX-ARSp time discretization, the convective CFL stability condition is obtained as due to semi-implicit nature of the acoustic/advection splitting strategy. In order to avoid an excessively large time-step sizes when the velocity of the fluid is too small, Boscarino et al. [3] recommended the following stability restriction The stability conditions (38) and (39) of the IMEX-ARSp method are independent of the Mach number ε and thus are less restrictive contrary to a standard explicit method. Therefore, the IMEX-ARSp method is more efficient than the explicit method for low-Mach number flows.

Remark 6
For one-dimensional case, the term p = p xx in (33) is approximated by a fourth-order central-type difference scheme: For two-dimensional case, this central-type difference scheme in a dimension-bydimension form is applied to p = p xx + p yy .

AP property
Definition 3 A scheme for the low Mach number limit of the Euler equations is asymptotic preserving, if it serves as a consistent and stable discretisation, independent of the Mach number ε. Proposition 1 Consider the system (1) subject to well-prepared initial data (i.e., ρ 0 0 = const, ∇ · u 0 0 (x) = 0) and periodic, wall or homogeneous Neumann boundary conditions. The spatial derivatives are assumed to be continuous. The GSA IMEX-ARSp scheme ((29)- (32)) (see Definitions 1 and 2) is used for the time discretization of the system (1) from t = t 0 to t 1 with the time-step size t = t 1 − t 0 given by the stability conditions (38) or (39) independent of ε, and we then have and Here, u e is the exact solution of Eq. (5) with initial data u e (x, 0) = u 0 0 (x).

Remark 8 The AP property of the GSA IMEX-ARSp scheme is proved based on the assumption that the Euler equations are subject to well-prepared initial data. For the SA IMEX-ARSp scheme, the AP property can also be proved using a similar procedure under the same assumption. The proof of the AP property of the GSA IMEX-ARSp scheme combined with the fifth-order WCNS is relatively complex and has not been performed.
In the next section, we will verify numerically the AP property of the IMEX-ARSp WCNS. In addition, the IMEX-ARSp scheme combined some other shock capturing schemes, e.g., the MUSCL (monotone upstream-centered scheme for conservation laws) type finite volume scheme proposed in [44], can also maintain numerically the AP property.

Numerical examples
The third-order GSA IMEX-ARS3 WCNS is applied to solve the following numerical examples. In some convergence tests, the third-order SA IMEX-ARS3 WCNS is also tested. See Appendix for the characteristics of the GSA and SA IMEX-ARS3 scheme. In addition, the third-order explicit total variation diminishing Runge-Kutta time-stepping method [30] coupled with the fifth-order WCNS (TVD-RK3 WCNS) and the thirdorder GAS IMEX-ARS3 scheme coupled with the third-order MUSCL type finite volume scheme (IMEX-ARS3 MUSCL) are used to solve some examples for comparisons.
The order of temporal accuracy is numerically calculated as where E t 1 obtained with t 1 = O(h) and E t 2 obtained with t 2 = t 1 /2 are the global errors in the L ∞ and L 1 norms, respectively. The space-step size h deceases with the timestep size t. In the convergence tests, to test the temporal order of accuracy the time-step size is taken as t = CFL · h (CFL = 0.1), while it is set as t = CFL · h 2 (CFL = 0.1) to test the spatial order of accuracy. In the other tests, the time-step sizes of the IMEX-ARS3 WCNS and the IMEX-ARS3 MUSCL are computed by the stability conditions (38) or (39), while the time-step size of the TVD-RK3 WCNS is calculated by the stability condition (37).
where γ = 2, c = 2 √ γ /ε, L = 5 and the final time is fixed to T f = 0.  [13]. Tables 5 and 6 report the corresponding results at t = 0.05, which demonstrate a fifth-order spatial convergence of the GSA IMEX-ARS3 WCNS at different Mach numbers. [4,6,10] This example that consists in several interacting Riemann problems is used to validate the asymptotic stability of the IMEX-ARS3 WCNS. The computational domain is [ 0, 1] with periodic boundary conditions. The initial data is described as

Example 2 One-dimensional Riemann problem
and the pressure is p(ρ) = ρ 2 . The time-step size is computed by (38) Fig. 1a and b). The utilization of WENO interpolation based on the characteristic variables instead of the flux functions can effectively suppress spurious oscillations near discontinuities. Compared to the TVD-RK3 WCNS, the IMEX-ARS3 WCNS and the IMEX-ARS3 MUSCL capture the limit incompressible solution faster when ε is small. When ε approaches 0, the TVD-RK3 WCNS can not capture the correct solution, due to the excessive numerical viscosity.  putational domain is [ 0, 1] with zero Neumann boundary conditions. The pressure is p(ρ) = ρ. The non-well prepared initial data including discontinuity is

Example 3 Sod shock tube problem [5] This problem with non-well prepared initial data is used to investigate the performance of the GSA IMEX-ARS3 WCNS in the intermediate and incompressible regimes. The com-
The time-step size is computed by (39) Fig. 2a and b, it is observed that when ε = 0.3, the numerical solution obtained with the IMEX-ARS3 WCNS has a good agreement with the reference solution. When ε = 0.03 in the incompressible regime, nonphysical oscillations arise near discontinuities based on the IMEX-ARS3 WCNS (see Fig. 2c and  d). Nonphysical oscillations are a common phenomenon when shocks develop in the low Mach number regime [4,5] and can be eliminated by decreasing t and taking a smaller CFL = 0.01. [19] This example is used to test the temporal order of accuracy of the GSA IMEX-ARS3 WCNS for two-dimensional case. The computational domain is [ 0, 1] ×[ 0, 1] with periodic boundary conditions. The pressure is p(ρ) = 1 2 ρ 2 . The final time is T f = 0.02. The initial data is given by the following exact solution  ρ(x, t) = 1 + ε 2 1.5 2 16π 2 δ(r)(k(r) − k(π)), u(x, t) = 0.6 + 1.5(1 + cos(r))δ(r)(0.5 − y), v(x, t) = 1.5(1 + cos(r))δ(r)(x − 0.6t − 0.5),  zero Neumann boundary conditions. The pressure is p(ρ) = ρ 1.4 . This initial data includes four shock waves and is described as

Example 4 Convergence test for two-dimensional case
The time-step size is computed by (38) with CFL = 0.4. The values of the Mach number are chosen as ε = 1 and 2 to test the robust shock-capturing capability of the IMEX-ARS3 WCNS in the compressible regime. Figure 3 plots the surfaces of the density profile obtained with the IMEX-ARS3 WCNS and the IMEX-ARS3 MUSCL on a mesh 50 × 50 at T f = 0.1. From this figure, we see that the IMEX-ARS3 WCNS has a higher resolution, compared to the IMEX-ARS3 MUSCL.

Example 6
Cylindrical explosion problem [6,11] The computational domain is [ −1, 1] ×[ −1, 1] with periodic boundary conditions and the pressure is given by p(ρ) = ρ 2 . The initial density of the fluid is given by where r = x 2 + y 2 is the distance to the center (0, 0) of the domain. The initial velocity of the fluid is given by    (38) with CFL = 0.4 for ε = 1 and by (39) with CFL = 0.03 for ε = 0.001, respectively. Figure 4 plots the surfaces of the density profile and the velocity fields for ε = 1 obtained with the GSA IMEX-ARS3 WCNS at different times. Figure 5 shows the density profile and the discrete divergence of the velocity ∇ · u at time t = 0.05 for a small value of ε = 0.001. It is observed that when ε = 0.001, the density converges to the constant value 1 and the divergence ∇ · u is close to 0.
Example 7 Shear flow problem [3,10,16] This example is used to verify numerically the AP property of the IMEX-ARS3 WCNS in the low Mach number limit. The computational domain is [ 0, 2π] ×[ 0, 2π] with periodic The time-step size is computed by (39) with CFL = 0.2. The Mach number is ε = 10 −4 and the final time is T f = 10. Figure 6 demonstrates the discrete vorticity ω = v x − u y on a mesh 50 × 50 at different times. From this figure, we see the vorticity quickly develops into roll-ups with smaller and smaller scales with the time evolution. Figure 7 shows the time history of the L 1 -norm error of the discrete divergence of the velocity ∇ · u. It is observed that the L 1 -norm error increases with time, which is in agreement with the phenomenon in [10] and is caused by the appearance of the smaller and smaller scale structures in the shear flow. The maximum error is smaller than that obtained with a finer mesh in [10].

Example 8 Kelvin-Helmholtz instability problem [10]
This problem is solved in the domain [ 0, 2π] ×[ 0, 4π] with periodic boundary conditions. The pressure is p(ρ) = ρ 2 and the initial data reads ρ(x, 0) = 1, u(x, 0) = cos(y), v(x, 0) = 0.03 sin(0.5x). (69) The time-step size is computed by (39) with CFL = 0.2. The Mach number is ε = 10 −3 and the final time is T f = 40. The discrete vorticity ω = v x − u y on a mesh 50 × 50 at different times is depicted in Fig. 8, which shows the emergence and subsequent evolution of small scale flow features. The time history of the L 1 -norm error of the discrete divergence of the velocity ∇ ·u is shown in Fig. 9, from which we also observe that the L 1 -norm error increases with time and the maximum error is smaller than that obtained with a finer mesh in [10]. In this work, we have designed a high order IMEX WCNS for the compressibleincompressible limit problems described as scaled isentropic Euler equations. Based on the acoustic/advection splitting strategy, the semi-implicit method combines the thirdorder IMEX-ARS3 scheme for the time integration with the fifth-order WCNS for the spatial discretization. The GSA IMEX-ARS3 scheme has been proven to be asymptotic preserving, i.e., it converges to a consistent scheme in the zero Mach number limit, and asymptotically stable under a convective CFL condition. The IMEX-ARS3 WCNS has been verified numerically that it can achieve a third-order temporal accuracy for different values of the Mach number, from compressible to incompressible regimes. Numerical results show that the IMEX-ARS3 WCNS with the AP property can capture shocks and contact discontinuities with high resolution in the compressible regime, while it can capture incompressible features in the imcompressible regime. The future work will focus on the implementation of the proof of the AP property of the fully discrete scheme and the extension to compressible Navier-Stokes flows at any Mach number.

Appendix
This appendix is used to characterize the GSA and SA IMEX-ARSp schemes adopted in this paper. The third-order GSA IMEX-RK3 scheme is characterized by