1 Introduction

It is well known that high-order methods have superior resolution properties and can be computationally efficient when applied to simulate complex flow fields with broadband length scales, such as direct numerical simulation (DNS) and large eddy simulation (LES) of turbulence^[1] and aeroacoustics^[2]. Especially, since high-order compact finite difference methods admit spectral-like resolution properties^[3], they are usually more accurate than explicit ones with the same order when applied to solve problems with smooth solutions^[4-5].

In this paper, we focus on a seventh-order dissipative compact finite difference scheme^[6]. In the previous work, to stabilize the scheme for long time simulations, fourth-order boundary closures are usually used^[7-9], resulting in an only fifth-order convergence rate for first-order hyperbolic conservation laws according to the result of Ref. [10]. To let the scheme achieve an optimal seventh-order convergence rate, we need to elevate the boundary closures to be at least sixth-order. However, we may encounter a time instability issue in that case.

Time stability^[11] (also called strict stability^[12-13]) is a property that is very important for long time simulations. To ensure time stability of the seventh-order scheme with sixth-order boundary closures, we intruduce the simultaneous approximation terms (SATs), which are successfully used by the summation-by-parts community^[14-15]. We will discuss the implementation of SATs for linear advection equations, coupled linear system, and Euler equations.

The rest of this paper is arranged as follows. In Section 2, we revisit the seventh-order dissipative compact scheme with fourth-order boundary closures^[6]. Then, we elevate the boundary closures to be sixth-order in Section 3, such that seventh-order global accuracy is achieved. Moreover, time instability is also discussed. In Section 4, the SATs are introduced to stabilize the scheme. The implementation of the improved scheme for one-dimensional Euler equations is discussed in Section 5. Finally, we conclude this paper in Section 6.

2 Revisiting seventh-order compact dissipative scheme

To state the seventh-order dissipative compact scheme proposed in Ref. [6], let us consider the following one-dimensional scalar conservation law:

(1)

where u=u(ξ, t) denotes the conservative variable, and f(u) stands for the flux. For convenience of presentation, we assume that f'(u)> 0 for all ξ∈[0, 1]. Therefore, to ensure a unique solution to Eq. (1), we specify an initial condition u(x, 0)=u₀(x) and a left boundary condition u(0, t)=g(t).

The grid adopted here is staggered by solution points ξ_j=(j-1)h with 1≤j≤N, and flux points ξ_j+1/2=(j-1/2)h with 1≤j≤N-1, as shown in Fig. 1, where h=1/(N-1) is the spatial step. Let u_j and f'_j be approximations to u(ξ_j, t) and f(u(ξ_j, t))_ξ, respectively. Similar notations are also introduced for flux points. Then, we can solve Eq. (1) by using the following procedure. (ⅰ) Compute the solution values u_j+1/2 at flux points by using interpolation schemes from the known values u_j. (ⅱ) Derive difference schemes to approximate the flux derivatives f'_j from the obtained flux values f_j=f(u_j) and f_j+1/2=f(u_j+1/2). (ⅲ) Apply some time-marching schemes to solve the resulting ordinary differential system,

(2)

Here, u₁=g(t). in the following, we will present the schemes proposed in ref. [6] for the steps (ⅰ) and (ⅱ), respectively.

2.1 Interpolation scheme

For the interior flux points ξ_j+1/2(3≤j≤N-3), the following seventh-order compact dissipative interpolation scheme is used:

(3)

where α is a constant controlling dissipation. Here, α=-0.3 is used, which can be determined by using an optimization procedure through the Fourier analysis^[6].

At the points ξ_5/2 and ξ_N-3/2, the fifth-order compact dissipative interpolation scheme is applied,

(4)

where the dissipative parameter is still chosen to be α=-0.3 for simplicity.

At the pointsξ_3/2 and ξ_N-1/2, fifth-order explicit interpolation schemes are derived to be

(5)

(6)

It is observed that Eqs. (3)-(6) result in a linear system with a tridiagonal matrix. Here, the traditional Thomas algorithm is used to solve the system.

2.2 Difference scheme

For the interior points ξ_j with 4≤j≤N-3, the following eighth-order central difference scheme is used:

(7)

For the other solution points, we have to design special near boundary difference schemes. At the solution points ξ₃ and ξ_N-2, we use the sixth-order central difference scheme,

(8)

For the solution points ξ₂, ξ_N-1, and ξ_N, the following fourth-order difference schemes are used:

(9)

(10)

2.3 Accuracy test

Since the finite difference scheme presented above is only fourth-order at solution points near domain boundaries, it is expected that the scheme can achieve at most fifth-order global accuracy for the conservation law (1) according to the classical result presented in Ref. [10]. To verify this numerically, we consider a numerical test for the one-dimensional scalar linear advection problem,

(12)

which admits the exact solution u(x, t)=sin(π(x-t))+sin(2π(x-t)).

Let us denote the resulting ordinary differential system after discretizing in space for Eq. (12) as

(13)

where U=[u₂, u₃, …, u_N]^T. Then a traditional fourth-order Runge-Kutta scheme^[16] for Eq. (13) can be written as

(14)

where ∆t denotes the fixed time step. To measure the errors between numerical results and exact solutions, we introduce the L^∞ norm and L² norm as follows:

(15)

Suppose e₁ and e₂ are the corresponding errors for different numbers of solution points N₁ and N₂, respectively. Then, we define the numerical order (O) of accuracy as

(16)

For numerical test, we consider N=20, 40, 80, and 160, and take the time step ∆t=1/10 000, which is small enough to make the errors in the time direction negligible here. As shown in Table 1, the numerical test for the problem (12) demonstrates that the finite difference scheme presented in Ref. [6] can only achieve fifth-order global accuracy as expected, even though the interior scheme is designed to be seventh-order. In the next section, we aim to elevate the global convergence rate to seventh-order by designing sixth-order boundary closures.

3 Global seventh-order scheme

To achieve seventh-order global accuracy, it is required that boundary closures have to be at least sixth-order^[10]. Here, we will modify the scheme presented in Section 2 to satisfy this requirement by using uniformly seventh-order interpolations and sixth-order difference schemes near the boundaries.

3.1 Near boundary interpolation schemes

For the interior points ξ_j+1/2 with 3≤j≤N-3, we still use the seventh-order dissipative compact interpolation scheme (3). We only need to modify the interpolations for the rest flux points, ξ_3/2, ξ_5/2, ξ_N-1/2, and ξ_N-3/2.

For the points ξ_3/2 and ξ_5/2, we can use the following seventh-order explicit interpolation scheme:

(17)

where the coefficients can be determined directly by using the Lagrangian interpolations,

(18)

(19)

For the points ξ_N-3/2 and ξ_N-1/2, the following seventh-order explicit interpolation scheme is used:

(20)

Again, by the Lagrangian interpolations, we can determine the coefficients,

(21)

(22)

3.2 Near boundary difference schemes

For the interior solution points ξ_j with 4≤j≤N-3, we still use the eighth-order scheme (7), while for the points ξ₃ and ξ_N-2, the sixth-order scheme (8) is applied. For the point ξ₂, the following sixth-order explicit difference scheme is designed:

(23)

For the other points ξ_N-1 and ξ_N, we derive the following sixth-order difference schemes:

(24)

(25)

3.3 Accuracy test

To show that the scheme designed here can achieve seventh-order global accuracy, we consider the initial boundary value problem (12). As shown in Table 2, we see that the modified scheme achieves seventh-order global accuracy approximately. It is noted that the final time t=0.1 is considered here, rather than t=1 considered in Subsection 2.3. The reason is that the modified scheme is actually not suitable for long time simulations. If we perform the computation for larger time, the errors arising from time integration would dominate the errors in the spatial direction. In that case, we would not observe a correct convergence rate of spatial discretization. This phenomenon can be better observed in Fig. 2. It is shown that large point-wise errors appear near the left boundary at time t=2, and the L²-norm error grows exponentially with increasing time. We can understand this phenomenon more deeply in the following section by eigenvalue analysis.

3.4 Eigenvalue analysis

As we have seen above, the modified scheme encounters a time instability issue. Here, we attempt to implement eigenvalue analysis to study the stability property of the scheme. To do this, we need to write the corresponding semi-discretized system (2) for Eq. (12) as a vector form.

At first, introduce the vector Û=[u_3/2, u_5/2, …, u_N-1/2]^T. Then, the interpolations (3), (17), and (20) can be written as

(26)

where Ũ=[u₁, U^T]^T with U defined in Eq. (13), and K and C are coefficient matrices that can be written down straightforwardly by observing the interpolations. To be more specific, the size of K is (N-1)×(N-1) while that of C is (N-1)×N.

Let U'=[u'₂, u'₃, …, u'_N]^T. Then, according to Eqs. (7), (8), and (23)-(25), we can express the difference scheme as a compact form,

(27)

where we have used Eq. (26) to obtain the last equality. Here, the sizes of the coefficient matrices B₁ and B₂ are (N-1)×N and (N-1)×(N-1), respectively. Therefore, we obtain the following semi-discretized system corresponding to Eq. (12):

(28)

For stability analysis of this linear system, we can confine the study to the case with a homogeneous boundary condition, i.e., u₁=0. Therefore, we have

(29)

where D is an (N-1)×(N-1) matrix that can be obtained from Eq. (28) by letting u₁=0.

Time stability of the system (29) requires that there is no eigenvalue of D in the right half complex plane. However, in Fig. 3, there are two eigenvalues with positive real parts, indicating that the solution of the system (29) will increase exponentially with time. This fact explains what we observed in Fig. 2. However, in practical calculations, we usually need to simulate problems for a long time. Therefore, we need to do some work to stabilize the modified scheme. In the next section, the SATs will be introduced to accomplish the job.

4 Time stable global seventh-order scheme

To stabilize the modified scheme presented in the last section. We will introduce the SATs to impose the boundary condition u(0, t)=g(t) here.

4.1 Construction of SATs

To impose the boundary condition u(0, t) weakly, we assume that u₁ is unknown when we consider spatial discretization. Therefore, in addition to the near boundary difference schemes presented in Subsection 3.2, we derive a sixth-order difference scheme at ξ₁,

(30)

Here, we impose the boundary condition by adding penalty terms (also called SATs) to the difference schemes that involve the value of u₁ (see Eqs. (30), (23), (8), and (7)), rather than by setting u₁=g(t) directly. To be more precise, we modify these difference schemes to be

(31)

(32)

(33)

(34)

Here, τ_j(j=1, 2, 3, 4) are penalty coefficients that can be tuned to stabilize the scheme. It is noted that we have chosen the signs of τ_j as the opposite signs of the coefficients of u₁ appearing in the original schemes, such that we can set here to obtain a stabilized scheme, even though they may be optimized further.

Remark 1 We try other choices of the signs of τ_j, but in the end we find that the selection presented in Eqs. (31)-(34) give the best stability property for a fixed value of τ, which will be analyzed in the following.

4.2 Stability analysis

Similar to the description in Subsection 3.4, we can write the improved semi-discretized system with SATs (see Eqs. (31)-(34)) as

(35)

where Ũ is defined in Eq. (29), is a modified coefficient matrix corresponding to D in Eq. (35), and E=[1, -1, 1, -1, 0, …, 0]^T. To analyze the time stability property of the system (35), we can set g=0 and write the result as

(36)

where with E₀=[1, 0, …, 0]^T. Here, we still choose N=40 for analysis (similar results can be observed for other values of N). From Fig. 4, we can see that the parameter τ can be tuned to let all eigenvalues of lie in the left half complex plane. For examples, the distributions of the eigenvalues of for two values of τ are shown in Fig. 5. For convenience, we just choose τ=5 hereinafter.

Remark 2 We observe that the maximum real part of the eigenvalues of the matrix is positive for all τ≤0. Thus, we omit the part for τ in Fig. 4.

Since non-uniform grids are often used in physical domains for practical problems, the above stability analysis based on uniform grids may not be sufficient. It is necessary to have a look at the cases when coordinate transforms are performed. To do so, let us consider the linear advection equation in physical coordinate x,

(37)

As an example, we consider the following physical grid generated by

(38)

where the stretching function φ is chosen as in Ref. [17] as follows:

(39)

with β being a parameter that controls stretching strength. For β=0, the transform (39) becomes linear, while for β≠0 it is nonlinear.

In the computational coordinate ξ, Eq. (37) reads

(40)

Then, a semi-discretized system corresponding to Eq. (36) can be written as

(41)

where P=diag{1/(x_ξ)₁, 1/(x_ξ)₂, …, 1/(x_ξ)_N}. As we can see from Fig. 6, the maximum real part of the eigenvalues of the matrix remains to be negative for different values of β with τ=5. Thus, we expect that the improved scheme with the SATs is suitable for the cases with nonuniform grids used in the physical space.

4.3 Numerical examples 4.3.1 One-dimensional linear advection equation

Here, we still consider the initial boundary value problem (12). As we can see from Table 3, seventh-order global accuracy is observed for the improved scheme with SATs. This demonstrates that the SATs do not affect the convergence rate of the original scheme. For long time simulations, as shown in Fig. 7, no exponential increase is observed, indicating that the improved scheme is time stable for this case.

4.3.2 Two-dimensional linear advection equation

We consider an initial boundary value problem for a two-dimensional linear advection equation,

(42)

This problem admits the exact solution (x, y, t)=sin(π(x+y-2t))+sin(2π(x+y-2t)). The spatial steps for the two spatial directions are both chosen to be h=1/(N-1) for simplicity. As we can see from Table 4, seventh-order global accuracy is observed for the improved scheme with the SATs. In addition, no exponential increase can be seen for the long time simulation as shown in Fig. 8, indicating that the scheme is time stable for this two-dimensional case.

4.3.3 Coupled linear system

Here, in the domain 0≤x≤1, we consider the following initial boundary value problem^[11]:

(43)

which admits the exact solution u(x, t)=sin(2π(x-t)), v(x, t)=-sin(2π(x+t)). Even though u and v are decoupled in the governed equations, they are coupled by the boundary conditions. This problem is often used to test time stability property of schemes by the summation-by-parts community.

The semi-discretized system corresponding to Eq. (43) can be written as

(44)

(45)

where Ṽ=[v₁, v₂, …, v_N]^T, is a derivative matrix corresponding to but with upwind direction swaped, E is defined in Eq. (35), and Ẽ=[0, …0, 1, -1, 1, -1]^T. It is noted that the coupled boundary conditions has been imposed weakly by using the SATs. The results illustrated in Table 5 show that seventh-order global accuracy is also achieved in this case. Moreover, the time stability is confirmed by the result of the long time simulation as shown in Fig. 9.

5 Generalization to Euler equations

In this section, we intend to apply the improved scheme with SATs to solve the following one-dimensional Euler equations:

(46)

where U=[ρ, ρu, E]^T is the conservative variable, and F=[ρu, ρu²+p, (p+E)u]^T is the flux. Here, ρ denotes the density, u denotes the velocity, and p denotes the pressure. The equation of state is set to be E=p/(γ-1)+ρu²/2 with γ being 1.4 for air. In particular, we consider the computational domain 0≤ξ≤1 with boundary conditions denoted to be

(47)

For the Euler equations (46), we can use the uniformly seventh-order interpolation schemes (see Subsection 3.1) to obtain left and right values at ξ_j+1/2, denoted by U_j+1/2^L and U_j+1/2^R, respectively. Then, we can obtain a corresponding numerical flux F_j+1/2(U_j+1/2^L, U_j+1/2^R) by using an approximate Riemann solver. Here, the Roe's Riemann solver^[18] is used. By using these obtained fluxes, we can compute the derivatives of the flux F'_j at ξ_j by using the improved difference scheme (see Subsection 3.2), yielding an ordinary differential system

(48)

As the linear case, we shall not expect the system (48) to be time stable. The SATs are needed to stabilize it.

5.1 Boundary treatment

Similarly to the linear case, the SATs should be introduced for the Euler equations (46) when boundary closures involve boundary points ξ₁ or ξ_N. To be specific, we should impose the SATs for Eq. (48) when 1≤j≤4 and N-3≤j≤N.

To introduce appropriate SATs, it is beneficial to express Eq. (46) as a quasi-linear form^[19],

(49)

where A(U)=F'(U) is the Jacobian matrix. It is well known that A(U) can be diagonalized as

(50)

where the sound speed , and the total specific enthalpy H=u²/2+γp/((γ-1)ρ). We first split A into two terms: A=A⁺ + A^- with A^±=KΛ^±K^-1. Here, Λ⁺=diag[(u-a)⁺, u⁺, (u+a)⁺] and Λ^-=diag[(u-a)^-, u^-, (u+a)^-] with the notations λ^±=(λ±|λ|)/2. Then, we impose the SATs to the system (48) by using the following modified near boundary schemes:

(51)

(52)

where the penalty coefficient τ is still fixed to be 5, G^L and G^R are boundary conditions (47), and A_j^± denote the frozen values of A^± at ξ_j, evaluated by the solution values from the last time step when time-marching schemes are employed (see Ref. [20] for more explanations).

5.2 Numerical example

Consider the following Euler equations with a source term^[21]:

(53)

where the source term has the form

(54)

with A(x)=1+2.2(x-1.5)². The boundary conditions for the left and right sides and initial values are all set to be [ρ, u, p]=[1, 0.08, 1/γ], such that the considered problem (53) admits a stationary solution. Since this problem is isentropic, i.e., p/ρ^γ is a constant, we choose the error of the entropy to evaluate the convergence rate of numerical results. We compute the problem until the residues reach the machine zero. The numerical results are presented in Table 6, showing that the scheme achieves the designed seventh-order convergence rate.

6 Conclusions

To improve the convergence rate of a seventh-order dissipative compact finite difference scheme, we elevate the boundary closures to sixth-order, such that seventh-order global accuracy is ensured for the first-order hyperbolic conservation laws. In addition, the SATs are introduced to impose boundary conditions weakly, yielding a time stable scheme that is suitable for long time simulations. The time stability of the improved scheme is shown by numerical eigenvalue analysis. Numerical experiments for linear advection equations are also implemented to validate the improved scheme. Moreover, we also generalize the scheme to solve one-dimensional Euler equations. In the future, we intend to use the improved scheme to solve two-dimensional Euler equations and discuss the use in multi-block domains.