1 Introduction

In recent decades, high-order finite-difference methods have been widely used in different areas of science and engineering, e.g., large eddy simulation (LES)^[1-5], direct numerical simulation (DNS)^[6-8] of turbulence, computational aeroacoustics (CAA)^[9-11], computational electromagnetics (CEM)^[12], and magnetic hydrodynamics (MHD)^[13]. On a Cartesian mesh, it is not difficult to construct a high-order interior scheme. However, it may become challenging to derive proper boundary closures to keep the whole scheme simultaneously be time stable and of the same order as the interior^[14].

In this paper, we will study a seventh-order dissipative compact finite-difference scheme^[15], which possesses good dispersive and dissipative properties and thus is suitable for the simulation of complex flow^[16-17]. On a uniform mesh, the fourth-order boundary closures^[15] are usually used for this scheme to ensure the time stability. However, to achieve a global convergence rate of the whole scheme as the interior, the boundary closures must be at most one order lower than the interior for first-order hyperbolic problems^[18]. Therefore, fourth-order boundary closures usually result in a global fifth-order convergence rate, which is not the same as the interior. If the boundary closures are elevated to be the sixth-order, one will obtain a scheme that is time stable on a uniform mesh^[19]. However, the scheme may become unstable if a grid transform is used. To rectify this, a new type of boundary closures has been proposed in Ref. [19], and a scheme is obtained, which is not only globally conservative but also time stable even when grid transforms are applied.

Since there is a dissipative free parameter in the dissipative compact scheme, Fourier analysis is usually used to optimize it. However, in Ref. [19], the optimization is only done for the interior scheme on a uniform mesh, and the same optimized parameter is used for the boundary closures. Since there is a difference between the interior scheme and the boundary closures, the optimized parameter for the interior will not be the optimized one for the near boundary schemes. In other words, the free parameters appearing in the boundary closures should be optimized, and so are the different near boundary schemes. In this paper, we attempt to improve the stability property of the global seventh-order scheme^[19] by an optimization with the consideration of the near boundary schemes.

Intuitively, to improve the stability property of a scheme, the maximum real part of the eigenvalues of the coefficient matrix is expected to be as small as possible. In our case, several free parameters are involved when we take the boundary closures into account. Since a genetic algorithm^[20] is an appropriate candidate to optimize the boundary closures^[21], the optimization will be conducted for two cases. One is to only optimize the boundary closures while maintaining the interior dissipative parameter as it is used in Ref. [19]. The other is to optimize all the interpolations involving free parameters.

The rest of this paper is organized as follows. In Section 2, we recall the global seventh-order dissipative compact finite-difference scheme^[19]. In Section 3, we state the optimization problem solved by a generic algorithm. We analyze the optimized results in terms of the eigenvalue distribution, the ε-pseudospectra^[22-25], the short time behavior, and the Fourier analysis. In Section 4, examples are numerically carried out to confirm the effectiveness of the recommended scheme. Finally, we draw the conclusions in Section 5.

2 Global seventh-order scheme

In the following, we recall the global seventh-order dissipative compact finite-difference scheme^[19], and state the corresponding optimization problem.

Consider the scalar one-dimensional conservation law in the computational domain ξ ∈ [0, 1] as follows:

(1)

where u=u(ξ, t) is the conservative quantity, and f(u) the flux. For convenience of the algorithm description, let us assume that f′(u) > 0, i.e., the information propagates from the left to the right. In this case, to ensure a unique solution to Eq. (1), one has to specify an initial condition u(ξ, 0) and a left boundary condition u(0, t).

As shown in Fig. 1, grids staggered by the flux points ξ_j-1/2 and the solution points ξ_j are considered. Suppose that the domain is divided into N units with the length h=1/N. Then, the uniform flux points can be expressed as follows:

(2)

However, for the solution points, we have three nonuniform solution points at both the left and the right boundaries, which are determined by using a global conservation restriction for the difference scheme (more details can be found in Ref. [19]). To be precise, the solution points are

(3)

We call the mesh used here the semi-uniform mesh, and distinguish it with the traditional uniform one.

Let u_j and u_j-1/2 be the approximations to u(ξ_j, t) and u(ξ_j-1/2, t), respectively. Then, at each time level of the algorithm, we have the values u_1/2 and u_j according to the aforementioned initial and boundary conditions. The algorithm can be summarized briefly as the follows:

(ⅰ) Apply an interpolation scheme to get the values u_j-1/2 and then the corresponding fluxes f_j-1/2=f(u_j-1/2).

(ⅱ) Calculate the derivatives f(u_j)_ξ with a difference scheme denoted as f_j′.

(ⅲ) Use a time integration scheme to solve the resulting ordinary differential equations

Now, we are ready to present the interpolation and the difference schemes used in Ref. [19] for Steps (ⅰ) and (ⅱ).

2.1 Interpolation scheme

Since the values u_1/2 and u_j are known, the following seventh-order dissipative compact scheme can be used to obtain the required values u_j-1/2 at the flux points:

(4)

where α_j ∈ [0, 1] are the dissipative parameters, and the other unspecified notations are constant coefficients which can be determined directly by Lagrange interpolations^[19]. Technically, for 4 ≤ j ≤ N-2, one can obtain the coefficients b_{j, k}⁽¹⁾ and c_{j, k} first by requiring the scheme to be eighth-order when α_j=0, and then determine the coefficients b_{j, k}⁽²⁾ by the seventh-order requirement of the terms containing α_j.

Introduce the vectors

Then, the interpolation scheme (4) can be written as a compact form CU_h= , i.e., U_h=C^-1, where B and C are (N+1)× (N+1) matrices which can be determined from Eq. (4).

2.2 Difference scheme

To approximate f(u)_ξ, explicit eighth-order difference schemes are used for the interior solution points, and sixth-order schemes are used for the near boundary points, i.e.,

(5)

where a_{j, k} are coefficients which can be determined by Lagrange interpolations^[19] as the interpolation scheme.

Introduce the vectors

Then, one can write Eq. (5) as F′=DF_h/h, where D is the N×(N+1) difference matrix following Eq. (5).

2.3 Linear system

Now, consider a linear flux f(u)=u in Eq. (1), and apply the interpolation and difference schemes presented above for the spatial direction. Then, we obtain the following ordinary differential equations:

(6)

where U=(u₁, u₂, ..., u_N)^T, and is defined at the end of Subsection 2.1. For the purpose of stability analysis, we can confine the study to the case with homogeneous boundary conditions, i.e., u(0, t)=0. Then, if an N × N matrix A is introduced to denote the matrix attained by eliminating the first column of -DC^-1B, we can express Eq. (6) as follows:

(7)

Since the coefficient matrix A contains the free controlled parameters α_j appearing in the interpolation scheme (4), we have to determine them before we try to solve Eq. (7).

For the case 7≤ j ≤ N-5, the interpolation scheme (4) does not involve nonuniform solution points. Therefore, the coefficients are the same. We can just use the optimized parameter α_j=0.554 674 for this case as that used in Ref. [19]. However, for the cases 4≤ j≤ 6 or N-4≤ j≤ N-2, the interpolation scheme (4) involves nonuniform solution points, and hence results in different coefficients with the previous case. Even though the same optimized parameter as the interior can still be used^[19], it is not an optimized choice. Strictly speaking, we should apply additional optimizations to determine the free parameters for the boundary closures to get a more stable scheme, which is the task of the next section.

3 Optimization and analysis of the results

In this section, we assume that the free parameters α_j for the interior and boundary closures are different. However, we use a same dissipative parameter for the interior scheme (7≤ j ≤ N-5). For convenience of the notations, we denote

To have a more stable system (7), we require the maximum real part of the eigenvalues (denoted by λ_j) of A as small as possible for a fixed N. Therefore, we want to choose appropriate φ_k (1≤ k ≤ 7) to minimize the following objective function:

(8)

This optimization problem contains seven parameters, which can be solved by a genetic algorithm^[20].

3.1 Results of a genetic algorithm

To minimize the objection function (8), we resort to the genetic algorithm toolbox in MATLAB. To implement the simulation, we have to specify some inputs. Here, the population size is set to be 500, the selection size is set to be 50, the max generation is set to be 5 000, and the constraint tolerance is set to be 10^-8. For the purpose of comparisons, two situations are implemented, i.e., the cases with fixed φ₄ of 0.554 674 or free φ₄. The optimized results are shown in Table 1. As we can see, the parameters near the boundaries are close to one.

3.2 Eigenvalue distribution

To see the optimized results more clearly, we plot the distributions of the eigenvalues of the coefficient matrices for different cases. As shown in Fig. 2, the maximum real parts of the eigenvalues of the two optimized results are both less than those obtained in Ref. [19], indicating a better stability property. In addition, the result of the case with free φ₄ fixed seems to be the best.

To compare further, we consider the spectral radius of the coefficient matrix A for the three different cases. In general, a smaller spectral radius indicates that a larger time step can be used for a given time scheme. As shown in Table 2, the radii do not differ a lot for different cases even though the result of the last case with free φ₄ is the best.

It is noted that the coefficient matrix A in Eq. (7) is nonnormal since the boundary closures are used. To get deep insight into the nonnormal matrix A, we can resort to the so-called ε-pseudospectra^[23]. One definition of the ε-pseudospectra is stated in terms of the perturbations of A^[24] as follows:

where denotes the complex space, E is a random matrix, and the norm ‖·‖ is chosen to be the 2-norm through out this paper.

At first, we take a look at the ε-pseudospectra to see whether the eigenvalues of A are sensitive to perturbations or not. As shown in Fig. 3, the ε-pseudospectra of A are still distributed on the left half complex plane for the optimized schemes with different perturbations. Therefore, the schemes remain to be time stable when a small perturbation is involved.

3.3 Short time behavior

For a linear system with a nonnormal coefficient matrix, it has been pointed out that the simple eigenvalue analysis cannot reveal the short time tendency of the solution^[24]. For some nonnormal matrices, the linear system may blow up in a short time. Therefore, it makes sense to check the short time behavior of the solution in addition to its asymptotic tendency.

For the linear system (7), the solution can be expressed as u(ξ, t)=e^tA/hu(ξ, 0). Therefore, we have

(9)

where ‖·‖ is chosen to be 2-norm here. It is clear that the term ‖e^tA/h‖ converges to zero for t→ ∞ if the system (7) is time stable. For the consideration of a short time behavior, we plot ‖e^tA/h‖ as a function of t for different cases (see Fig. 4). We can see that the three considered cases display a similar short time behavior. Before decaying monotonically, the term ‖e^tA/h‖ increases a little bit for all of the cases. Fortunately, since the maximum height is less than two, we expect that the system will not blow up in a short time. Therefore, all of them are safe for short time simulations.

3.4 Fourier analysis

One of the advantages of high-order compact schemes is that they admit a very good spectral property. From this point of view, we analyze in this section the resolution of the schemes with the Fourier analysis. For the analysis of compact schemes, it is best to consider the system as a whole and use a global Fourier analysis method^[21].

At first, we assume

(10)

where ω is the wavenumber, ω_j^* is the modified wavenumber of the point ξ_j, and is the imaginary unit. Substitute this assumption into the spatial discretization (see Eq. (6)). Then, we have

(11)

Introduce ζ=(ξ₁, ξ₂, ..., ξ_N-1, ξ_N), and =(ξ_1/2, ξ₁, ξ₂, ..., ξ_N-1, ξ_N). Then, we can rewrite Eq. (11) as follows:

(12)

For the purpose of comparisons, we can plot the real and imaginary parts of a middle modified wave number ω_j^* as a function of ω for different cases (see Fig. 5). We can see that the results of the case with fixed φ₄ are similar to those presented in Ref. [19], while the results of the case with free φ₄ seem to have a dispersive plot closer to the exact one. However, if zooming in, we can see that there is a small overshooting interval for the case with free φ₄. Therefore, we recommend the reader to use the results with fixed φ₄ in simulations.

4 Numerical experiments

To show the effectiveness of the optimized scheme with fixed φ₄ of 0.554 674, we consider the examples governed by the one-dimensional (1D) Euler equations or the two-dimensional (2D) Euler equations. In these nonlinear cases, we choose a Roe flux splitting method^[26] to calculate the nonlinear fluxes for the scheme. In the first two cases, proper source terms are chosen to ensure that the equations admit stationary solutions. Therefore, we compute the solutions until the residuals in space reach the machine zero. The last test case is a vortex convection with periodic boundary conditions. For the test of numerical convergence rate, the following L₂ error is used for the 1D case:

(13)

where Δx_j = x_j+1/2-x_j-1/2. Similar definition is used for the 2D case. For different numbers of the grid points N₁ and N₂, denote the corresponding errors as e₁ and e₂, respectively. Then, the numerical order of accuracy (or convergence rate) is defined by

(14)

4.1 1D nozzle flow

Consider the nozzle flow governed by the following Euler equations with a source term:

(15)

where

(16)

Here, the total energy per unit mass is defined by E=p/(γ-1)+ρu²/2 with the specific ratio γ=1.4, the area distribution of the nozzle is chosen to be A(x)=1+2.2(x-1.5)². In the simulation, the computational domain is set to be x∈ [0, 3] with the same boundary condition (ρ, u, p)=(1.0, 0.08, 1.0/γ) at both sides such that Eq. (15) admits a stationary solution.

To solve Eq. (15), a third-order TVD Runge-Kutta method^[28] is used until the residuals reach the machine zero (see Fig. 6 for the numerical stationary solution for the grid number N=80). Since the problem is isentropic^[27], i.e., p/ρ^γ= constant, we can use the errors of the entropy to test the convergence rate of the scheme. As shown in Table 3, the L₂ errors of the entropy decay at a seventh-order convergence rate as expected.

4.2 2D Euler equations

We consider the 2D Euler equations with a source term as follows:

(17)

where

Here, E=p/(γ -1)+ρ (u²+v²)/2, and S(U) is the source term. Since our purpose is to test the order of accuracy in the spatial direction, we assume that Eq. (17) has the solution

(18)

Then, we choose the source term that can be computed analytically by Eq. (18) such that Eq. (18) is actually a stationary solution to Eq. (17).

Here, the computational domain is chosen to be x × y∈ [-0.2, 0.2] × [-0.5, 0.5]. In the spatial direction, the Roe splitting method is used to compute the nonlinear fluxes. To improve the efficiency, an implicit Euler scheme is implemented for the time direction until the residuals reach the machine zero, rather than an explicit one used in the previous example. For simplicity, we just take a look at the errors of the density ρ. As shown in Table 4, it is evident that the scheme achieves the designed seventh-order accuracy for the density ρ.

4.3 2D vortex convection

We consider a 2D vorticity wave convection governed by Eq. (17) with a vanishing source term, i.e., S(U)=0. The initial value condition is

(19)

with

(20)

where

The computational domain is chosen to be x× y ∈ [-0.5, 1]× [-0.75, 0.75] with the periodic boundary conditions in both the x- and y-directions. The exact solution to this problem can be obtained by simply moving the initial value (19) to the right periodically with the velocity u_∞^[11].

Even though we are considering a case with periodic boundary conditions, the boundary closures derived in this paper is still used for the spatial discretization, while the third-order TVD Runge-Kutta method^[28] used in Subsection 4.1 is used for the temporal integration. We compute the solution until it propagates for one period, and obtain the test of accuracy as shown in Table 5. From the table, we can see that the scheme achieves the designed seventh-order accuracy approximately.

5 Conclusions

We have investigated a global seventh-order dissipative compact finite-difference scheme that contains undetermined free dissipative parameters. With a generic algorithm, we optimize the scheme by assuming that the boundary closures use different dissipative parameters with the interior. We obtain the results for the case with a fixed or free interior parameter.

At first, from the comparison of the eigenvalue distributions, we find that the two optimized results improve the stability property obviously. The optimized scheme with a free interior parameter is slightly better than that with a fixed interior parameter in terms of the maximum real part of the eigenvalues. Furthermore, ε-pseudospectra are introduced to see whether the nonnormal coefficient matrices are sensitive to perturbations or not. We find that the optimized schemes remain to be time stable for small perturbations. Then, we consider the short time behavior of the linear systems with nonnormal coefficient matrices. We observe that the norms of these systems increase in a short time interval, but will not blow up since the increments are small. Hence, the schemes are safe in a short time interval as well.

In practice, we prefer the case with a fixed interior dissipative parameter in terms of the Fourier analysis. The reason is that the case with free interior dissipative parameter exhibits an overshooting dispersion behavior in a small wavenumber interval. To show the effectiveness of the optimized scheme, two examples governed by 1D and 2D Euler equations with source terms and a 2D vortex convection case are computed numerically. All the results show that the semi-dicretized nonlinear systems are time stable and the scheme achieves the designed seventh-order accuracy.