Next Article in Journal
Ensemble Model of the Financial Distress Prediction in Visegrad Group Countries
Next Article in Special Issue
Application of Said Ball Curve for Solving Fractional Differential-Algebraic Equations
Previous Article in Journal
Inertial Extragradient Methods for Solving Split Equilibrium Problems
Previous Article in Special Issue
Local Dynamics of Logistic Equation with Delay and Diffusion
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Alternating Polynomial Reconstruction Method for Hyperbolic Conservation Laws

1
College of Meteorology and Oceanography, National University of Defense Technology, Changsha 410000, China
2
Department of Mathematics, National University of Defense Technology, Changsha 410000, China
*
Author to whom correspondence should be addressed.
Mathematics 2021, 9(16), 1885; https://doi.org/10.3390/math9161885
Submission received: 9 July 2021 / Revised: 4 August 2021 / Accepted: 6 August 2021 / Published: 8 August 2021
(This article belongs to the Special Issue Recent Advances in Differential Equations and Applications)

Abstract

:
We propose a new multi-moment numerical solver for hyperbolic conservation laws by using the alternating polynomial reconstruction approach. Unlike existing multi-moment schemes, our approach updates model variables by implementing two polynomial reconstructions alternately. First, Hermite interpolation reconstructs the solution within the cell by matching the point-based variables containing both physical values and their spatial derivatives. Then the reconstructed solution is updated by the Euler method. Second, we solve a constrained least-squares problem to correct the updated solution to preserve the conservation laws. Our method enjoys the advantages of a compact numerical stencil and high-order accuracy. Fourier analysis also indicates that our method allows a larger CFL number compared with many other high-order schemes. By adding a proper amount of artificial viscosity, shock waves and other discontinuities can also be computed accurately and sharply without solving an approximated Riemann problem.

1. Introduction

In this paper, we aim to solve the hyperbolic conservation law
u t + f ( u ) x = 0 ,
where u and f ( u ) can be either scalars or vectors.
The classical higher-order numerical schemes such as the finite difference (FD) and the finite volume (FV) frameworks usually involve only one degree of freedom (DOF) per cell. That means to achieve the kth order discretization of the first-order spatial derivative, a stencil of at least k + 1 cells is required. The compactness of a scheme can help simplify the polynomial interpolation in unstructured meshes and reduce communication costs between computational nodes on modern supercomputers. In this paper, the “compact” is used specifically for numerical schemes with a minimal stencil, which excludes so-called Padê schemes [1,2] since its approximation of the derivative involves solving a diagonal system of the full difference stencil.
A natural thought for the design of a compact high-order scheme is involving more than one DOF per cell so as to construct higher-order polynomials under the same stencil width. A series of multi-moment schemes is developed based on this idea. The constrained interpolation profile (CIP) scheme [3,4], which uses the semi-Lagrangian method to advance in time, individually updates the physical variable and its derivatives on each cell by the governing equations in different forms. The interpolated differential operator (IDO) method [5] is similar to CIP but uses an explicit time-advancing approach. For the IDO method, a stencil containing three adjacent cells suffices to approximate the first-order spatial derivative with fifth-order accuracy in one dimension, which makes it possible to obtain comparably accurate results with spectral methods in direct numerical simulation [6]. For IDO, updating the derivative-value variables induces the problem of computing the time derivative of u x , which involves f u u . However, for the system of conservation laws, f u u has a fairly complicated form and can be difficult to compute. Goodrich et al. proposed a similar Hermite method [7] on staggered grids.
Neither CIP nor IDO are conservative schemes, which means they are not ideal for shock-capturing problems and long-term simulations. To address the non-conservativeness, so-called CIP-conservative semi-Lagrangian (CIP-CSL) schemes [8,9] and the IDO scheme in conservative form (IDO-CF) [10,11] add the cell-averaged value as an extra independent model variable to evolve. The cell-averaged variable is advanced by the flux-form conservation law and is hence exactly conserved, while other kinds of variables are updated using the differential form of conservation laws. However, these conservative CIP algorithms need to find polynomials that match cell-averaged variables across cells, which can be time-consuming when encountered with non-uniform meshes. A framework termed the multi-moment finite volume method (MM-FVM) has been proposed on the basis of CIP-CSL. Compared with CIP-CSL, MM-FVM can construct the high-order interpolation on a local base and is more flexible when treating unstructured grids, as shown in [12]. Nonetheless, MM-FVM uses the semi-Lagrangian method in the time integration of point-based variables, which is difficult to implement in the multi-dimensional case.
Another class of high-order conservative schemes based on local flux reconstructions also uses a compact stencil and has become increasingly attractive in recent decades. Examples include the well-known discontinuous Galerkin (DG) [13,14,15], spectral volume [16,17], and spectral difference [18] methods and, more recently, the flux reconstruction/correction via procedure (FR/CPR) method [19,20]. These methods evolve k DOFs in each cell for a kth order spatial accuracy (in one dimension) and compute the numerical flux on cell boundaries to represent the interaction between adjacent cells. Thus, polynomial interpolation across cells is not needed, which makes these methods compact and flexible in handling unstructured meshes. The multi-moment constrained finite volume (MCV) method proposed by Ii and Xiao [21] can be viewed as a combination of CIP and the flux reconstruction approach. The MCV method introduces high-order spatial derivatives of numerical flux and reconstructs the local flux function more directly. MCV is generalized as multi-moment constrained flux reconstruction (MMC-FR) [22]. However, the nonlinear aliasing phenomenon [23,24] in treating nonlinear flux functions, which causes instability and accuracy loss, is one typical demerit for these methods based on local flux reconstruction.
This paper presents a novel conservative multi-moment approach termed the AltPoly method, which uses the alternating polynomial reconstruction as the time integration method. AltPoly divides the computational domain into cells and adopts point-based and cell-averaged values as model variables. The point-based variables including the physical variable and its spatial derivatives are updated together, which is different from the conventional multi-moment methods that update different types of model variables according to governing equations in different forms. AltPoly then solves a constrained least-squares problem to reconstruct the solution within the cell and update the model variables, which preserves conservativeness exactly. Solving Riemann problems on cell boundaries is not needed. Furthermore, AltPoly can also handle non-uniform meshes in one dimension efficiently with the aid of the coordinate transform.
The remainder of this paper is organized as follows: Section 2 describes the formulation of AltPoly and constructs artificial viscosity for discontinuous problems. Section 3 studies the accuracy and stability of AltPoly by Fourier analysis. Section 4 validates the numerical convergence rate and capability of solving some widely used benchmark problems involving both the scalar equation and Euler systems in 1D. Finally, a short summary in Section 5 ends this paper.

2. Formulation of AltPoly Method

This section gives the formulation of the AltPoly method for a one-dimensional conservation law. For concreteness and ease of comprehension, we only present the detailed formulation for AltPoly with 2 point-value variables and one extra cell-averaged variable per cell, which is termed AltPoly-3. As for the general AltPoly-r with r model variables within a cell, the formulation is similar to AltPoly-3 and is hence briefly illustrated in Appendix A.
To begin with, we first introduce some notation. The bold font lower-case letter denotes a column vector, e.g., a , u . The transpose a is denoted by a T . e k denotes the kth standard unit vector in a space in proper dimensions. a k k = 1 m denotes a column vector in the form of a 1 , , a m T . Bold font capital letters such as A denote matrices. The submatrix of A constructed from rows { r 1 , r 2 , , r k } and columns { c 1 , c 2 , , c l } is denoted as A [ r 1 , r 2 , , r k ; c 1 , c 2 , , c l ] . a denotes the 2-norm of the vector a , while A for a matrix A is defined by A = sup x 0 A x x .

2.1. Spatial Discretization and Model Variables

Suppose the spatial domain [ x min , x max ] is divided into N non-overlapping cells, among which the j-th cell is I j : = [ x j 1 / 2 , x j + 1 / 2 ] , j = 1 , , N . It should be noted that the cells do not necessarily have a uniform length. In the j-th element, we consider computing the time integration of three model variables: u , u x at the middle point of the element, and additionally the element-integrated value u ¯ j = 1 Δ x j x j 1 / 2 x j + 1 / 2 u d x .

2.2. Updating Model Variables

Denote the model variables at time level t = t 0 as u j t 0 , ( u x ) j t 0 , u ¯ j t 0 . This part illustrates the procedure to obtain u j t 1 , ( u x ) j t 1 , u ¯ j t 1 with t 1 = t 0 + Δ t , where Δ t is the time step length. In a nutshell, this procedure consists of 3 steps:
  • Hermite interpolation across cells using point-based variables;
  • Computing solution values on selected solution points based on reconstructed polynomials and updating solution values and cell-averaged variables via the Euler forward method;
  • Reconstructing the solution polynomial that matches updated solution values under the constraint of cell-averaged variables to preserve conservativeness.
Details are presented in the following and algorithm illustrations are presented in Figure 1 for convenience of illustration.

2.2.1. Step 1: Hermite Interpolation across Cells

In the first step, we find the interpolation polynomial h j 1 2 ( x ) on the interval [ x j 1 , x j ] that matches u j 1 t 0 , ( u x ) j 1 t 0 , u j t 0 , ( u x ) j t 0 as shown in Figure 1a. This step does not involve cell-averaged variables u ¯ j t 0 . This interpolation polynomial is of the third order since there are 4 conditions in total. Before constructing h j 1 2 ( x ) , we introduce a local coordinate
s ( x ) : = x 1 2 ( x j + x j + 1 ) x j + 1 x j [ 1 , 1 ]
for x [ x j , x j + 1 ] as is usually done in the finite element method [25]. The local coordinate can unify the interpolation templates so as to reduce the computational cost. The polynomial h j 1 2 ( x ) on the local coordinate s is denoted by h ˜ j 1 2 ( s ) : = h j 1 2 ( x ( s ) ) . Notice that h ˜ s and h x can be connected by the following rule:
h ˜ ( s ) s = h ( x ( s ) ) s = h x x s = ( x j x j 1 ) h x .
Then the coefficients of h ˜ j 1 2 ( s ) = i = 0 3 a i s i can be determined by the following system:
h ˜ j 1 2 ( 1 ) = k = 0 3 ( 1 ) k a k = u j 1 t 0 , h ˜ j 1 2 ( 1 ) s = k = 0 3 ( 1 ) k 1 k a k = ( x j x j 1 ) · ( u x ) j 1 t 0 , h ˜ j 1 2 ( 1 ) = k = 0 3 a k = u j t 0 , h ˜ j 1 2 ( 1 ) s = k = 0 3 k a k = ( x j x j 1 ) · ( u x ) j t 0 .
The matrix form is
Ha = d ,
where
H = 1 1 1 1 0 1 2 3 1 1 1 1 0 1 2 3 , a = a 0 a 1 a 2 a 3 , and d = u j 1 t 0 ( x j x j 1 ) · ( u x ) j 1 t 0 u j t 0 ( x j x j 1 ) · ( u x ) j t 0 .
In practice, the inverse of H is computed and stored in advance so as to reduce the computation cost in solving (4).

2.2.2. Step 2: Updating Solution Values on Solution Points

Based on Step 1, we can now approximate and update the solution in cell I j using polynomials h j 1 2 ( x ) and h j + 1 2 ( x ) . For AltPoly-3, four symmetric solution points x j , k = x j + ξ k Δ x j / 2 , k = 1 , , 4 on I j are needed in this step with ξ k [ 1 , 1 ] . In this paper, we choose ξ 1 = 1 , ξ 2 = 1 + σ , ξ 3 = 1 σ and ξ 4 = 1 as solution points with σ = 0.03 . To update the solution values on these solution points, we need to compute
u j , k = h j 1 2 ( x j , k ) , if x j , k x j ; h j + 1 2 ( x j , k ) , if x j , k > x j , and ( u x ) j , k = x h j 1 2 ( x j , k ) , if x j , k x j ; x h j + 1 2 ( x j , k ) , if x j , k > x j ,
to approximate u and u x . For convenience, we transfer the solution points ξ k into local coordinate s ( ξ k ) for k = 1 , 4 . Combining the conservative law (1) we can then compute u j , k ( t 1 ) by the Euler forward method:
u j , k t 1 = u j , k f u ( u j , k ) · ( u x ) j , k Δ t ,
where the function f u stands for f ( u ) u .
The cell-averaged variable u ¯ j ( t 1 ) is updated by the flux form of the conservation law:
u ¯ j t 1 = u ¯ j ( t 0 ) f ( u j + 1 2 ) f ( u j 1 2 ) Δ x j Δ t .
where u j + 1 2 is the function value at the cell boundary x j + 1 / 2 computed by
u j + 1 / 2 = h j + 1 2 ( x j + 1 2 ) = h ˜ j + 1 2 ( s ( x j + 1 2 ) ) .
When the computing grid is uniform, (10) is simply u j + 1 / 2 = h ˜ j + 1 2 ( 0 ) .
Figure 1b describes this step.
Remark 1.
The Euler time integration used for updating solution values on solution points is simple and easy to implement. However, the semi-Lagrangian time integration scheme [26] might be robust with a larger time step.

2.2.3. Step 3: Updating Model Variables

The cell-averaged variable u ¯ j has been updated (9). As for the point-value model variables, the four updated solution values on I j suffice to reconstruct a Lagrange interpolation polynomial l j ( x ) , so as to obtain the updated model variables u j t 1 and ( u x ) j t 1 . However, such a direct routine may cause the mismatch of u j and u ¯ j .
Similar to Step 1, we introduce a local coordinate ξ [ 1 , 1 ] :
ξ ( x ) : = x x j Δ x j , for x I j .
Let l ˜ j ( ξ ) = i = 0 3 b i ξ i be l j ( x ) under the local coordinate ξ . Then all conditions for l ˜ j ( ξ ) are listed as
l ˜ j ( ξ k ) = i = 0 3 b i ξ i = u j , k ( t 1 ) , for k = 1 , , 4 ,
1 2 1 1 l ˜ j ( ξ ) d ξ = b 0 + 1 3 b 2 = u ¯ j ( t 1 ) .
This is an over-determined linear system since it contains 5 conditions and 4 unknowns. However, (13) must hold for ensuring the conservation law. In other words, we need to find the optimal coefficients b i i = 1 , , 4 to fit (12) under the constraint (13). This is a typical constrained least-squares problem. To solve this, we rewrite the constraint (13) as
b 2 = 3 ( u ¯ j ( t 1 ) b 0 ) ,
and plug it into (12). Then the constrained least-squares problem is cast into the standard least-squares problem as
A b = y ,
where
A = 1 ξ 1 2 ξ 1 ξ 1 3 1 ξ 4 2 ξ 4 ξ 4 3 , b = b 0 b 1 b 3 , y = u j , 1 ( t 1 ) u j , 4 ( t 1 ) 3 u ¯ j ( t 1 ) ξ 1 2 ξ 4 2 .
Then the solution of (15) is
b = A y = ( A T A ) 1 A T y ,
where A is the generalized inverse or the Moore–Penrose inverse [27] of non-square matrix A .
From the polynomial expression of l j ( x ) , we can now retrieve the point-value model variables for the next time level:
u j ( t 1 ) = l j ( x j ) = l ˜ j ( 0 ) = b 0 ,
( u x ) j ( t 1 ) = l j ( x j ) x = l ˜ j ( 0 ) ξ ξ x = b 1 Δ x j .
So far, (9), (17) and (18) together give the updating rule for all model variables.
This polynomial reconstruction of l j ( x ) based on solving constrained squares problem (13) can be efficiently implemented if we compute and store A in advance.

2.3. Runge–Kutta Time Integration

Up to this point, we have only presented AltPoly using the Euler forward method, which is merely of first-order temporal accuracy. Since the presented AltPoly is not in the semi-discretization form, the Runge–Kutta algorithm cannot be directly applied for AltPoly. Therefore, a slight modification of the Runge–Kutta algorithm is adopted here.
To proceed, we introduce the following notation. Let u be a compounded list of all model variables, and its j-th entry is a vector u j = u j , ( u x ) j , u ¯ j . Denote the result of AltPoly upon u after one step with length Δ t of the temporal forward Euler method by E ( u , Δ t ) and the corresponding increment part is denoted by
D ( u , Δ t ) = E ( u , Δ t ) u .
Then, one time step of the Runge–Kutta procedure for AltPoly is
R ( u t 0 ) = u ( t 0 ) + 1 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) ,
where
k 1 = D ( u ( t 0 ) , Δ t ) ,
k 2 = D ( u ( t 0 ) + 1 2 k 1 , 1 2 Δ t ) ,
k 3 = D ( u ( t 0 ) + 1 2 k 2 , 1 2 Δ t ) ,
k 4 = D ( u ( t 0 ) + k 3 , Δ t ) .

2.4. Artificial Viscosity

To tackle problems that involve discontinuities such as shocks and contact discontinuities, we add artificial viscosities to the proposed method. The critical step is depicting the smoothness of a cell and determining whether the discontinuity exists. Our shock sensor borrows the idea from the sub-cell shock capturing for DG [28].
Within each cell, we define the following truncated cell average based on the moments for AltPoly-r:
u ¯ j trun = k = 0 k ˜ 1 2 k + 1 2 k u x 2 k j ( Δ x j ) 2 k ,
where 2 k ˜ is the largest even number smaller than r. It can be seen that u ¯ j trun utilizes the even-order partial differential point-value moments and is an approximation of u ¯ j with the error of magnitude O Δ x j 2 k ˜ + 2 . Therefore, we expect that the solution is continuous and then u ¯ j trun is close enough to u ¯ j , and hence declaim that a discontinuity exists in I j if | u ¯ j u ¯ j trun | exceeds some threshold. Based on this, the following smoothness indicator is defined:
θ j = | u ¯ j u ¯ j trun | u range + ϵ ,
where u range = u max u min and the parameter ϵ = 10 10 is added to avoid dividing zero. This smoothness indicator has a simple form and can be calculated directly from the model variables. Then the amount of viscosity at x j + 1 2 is calculated by the following smooth function,
ν j = 0 , if ϑ j ( , ϑ ) ; ν ˜ 2 1 + sin π ( ϑ j ϑ ) κ ϑ if ϑ j [ ϑ , κ ] ; ν ˜ if ϑ j ( κ , + ) ,
where ϑ j = log 10 θ j , the parameters ϑ ( Δ x j ) 2 k ˜ , ν ˜ = Δ x j / 4 , and κ is large enough to maintain a sharp but smooth shock profile.
So far, we have obtained the viscosity at the middle of each cell. Next we use the linear interpolation to construct the viscosity function between the middle points of two adjacent cells. Then, to solve the problem with this viscosity function, one only needs to replace f ( u ) with f vis ( u , u x ) = f ( u ) ν ( x ) u x , and modify (8) and (9), respectively.

3. Fourier (von Neumann) Analysis

This part gives the Fourier analyses of stability and accuracy for AltPoly in solving a linear advection equation. We only give a detailed analysis for AltPoly-3 for concreteness. AltPoly-r with other choices of r can be analyzed following a similar routine.

3.1. Formulation of the Amplification Matrix

Consider the linear scalar equation
u t + u x = 0 .
Suppose the computational domain is [ 0 , L ] with a uniform mesh spacing Δ x . The initial condition is periodic: u init ( x ) = e i w x / Δ x where i denotes the imaginary unit, and w = 2 π k Δ x / L [ 0 , π ) is the scaled wavenumber. The initial values of model variables on the cell I j = [ x j 1 / 2 , x j + 1 / 2 ] and its neighbor cells are given by
u j = e i w x j / Δ x , u j ± 1 = e ± i w u j ,
( u x ) j = i w Δ x e i w x j / Δ x , ( u x ) j ± 1 = e ± i w ( u x ) j ,
u ¯ j = 1 i w ( e i w / 2 e i w / 2 ) e i w x j / Δ x , u ¯ j ± 1 = e ± i w u ¯ j .
The discrete Fourier transform of the series u j is
u ^ k | u ^ k = j u j e 2 π i k x j / L .
u ^ k is not 0 only when k = k under the given initial condition. Hence, we only need to consider the Fourier coefficient u ^ k and denote it as u ^ in the following text for simplicity. u x ^ and u ¯ ^ are defined similarly for series ( u x ) j and u ¯ j .
The key procedure in our analysis is rewriting the AltPoly-3 formulation in Section 2 into the matrix form:
u ^ ( t 1 ) = S u ^ ( t 0 ) ,
where u ^ t 0 is defined as u ^ t 0 = [ u ^ t 0 , u x ^ t 0 Δ x / 2 , u ¯ ^ j t 0 ] , and S denotes the amplification matrix connecting u ^ t 1 at the next time level t = t 1 = t 0 + Δ t .
The coefficient vector a j ± 1 / 2 of the reconstructed polynomial h ˜ ( s ) in Step 1 of AltPoly-3 is
a j ± 1 / 2 = H 1 d j ± 1 / 2 = H 1 P u ^ t 0 e i w x j / Δ x
where d j 1 / 2 = u j 1 t 0 , ( u x ) j 1 t 0 Δ x / 2 , u j t 0 , ( u x ) j t 0 Δ x / 2 T ,
P = e i w 0 0 0 e i w 0 1 0 0 0 1 0 .
According to (26), we have d j + 1 / 2 = e i w d j 1 / 2 and a j + 1 / 2 = e i w a j 1 / 2 according to the periodicity of the spatial profile. According to (10), the function values at the cell boundaries x j ± 1 / 2 are
u j 1 / 2 ( 0 ) = e 1 T a j 1 / 2 , u j + 1 / 2 ( 0 ) = e i w u j 1 / 2 ,
where e k denotes the k-th standard unit vector.
Now we turn to Step 2 of AltPoly-3. As the solution points x j , 1 and x j , 2 are in the domain of h j 1 / 2 ( x ) , we transfer them into the local coordinate in [ x j 1 , x j ] . Let s 1 = s ( x j , 1 ) = 0 and s 2 = s ( x j , 2 ) = α . Then we can derive
u j , 1 u j , 2 = 1 s 1 s 1 2 s 1 3 1 s 2 s 2 2 s 2 3 a j 1 2 , ( u x ) j , 1 ( u x ) j , 2 = 2 Δ x 0 1 2 s 1 3 s 1 2 0 1 2 s 2 3 s 2 2 a j 1 2 .
Similarly, the variable value and first spatial derivatives on solution points x j , 3 , x j , 4 can be computed by
u j , 3 u j , 4 = 1 s 3 s 3 2 s 3 3 1 s 4 s 4 2 s 4 3 a j + 1 2 , ( u x ) j , 3 ( u x ) j , 4 = 2 Δ x 0 1 2 s 3 3 s 3 2 0 1 2 s 4 3 s 4 2 a j + 1 2 ,
where s 3 = α and s 4 = 0 .
Combining a j + 1 / 2 = e i w a j 1 / 2 , the updating of solution values on solution points can be expressed in the following matrix form:
[ u j , m ] 1 m 4 = W B a j 1 / 2 , [ ( u x ) j , m ] 1 m 4 = 2 Δ x W C a j 1 / 2
where the items of B , C are, respectively:
B m , n = s m n 1 , C m , n = ( n 1 ) s m n 2 for 1 m , n 4 ,
with 0 0 defined as 1, and W is a diagonal matrix with 1 , 1 , e i w , e i w on the diagonal line.
Since f u = 1 , the updated solution values here are computed by
[ u j , m t 1 ] m = 1 4 = WB 2 c WC a j 1 2 ,
where the CFL number c = | u f ( u ) | Δ t / Δ x = Δ t / Δ x for this linear problem.
As for Step 3, the updated cell average over the jth cell is
u ¯ j t 1 = u ¯ j c ( u j + 1 / 2 u j 1 / 2 ) = u ¯ j c ( e i w 1 ) e 1 T a j 1 / 2 .
Solving a constrained least-squares problem we have
b = A [ u j , m t 1 ] 1 m 4 3 [ ξ m 2 ] 1 m 4 u ¯ j Δ t .
From b we can easily obtain point-value variables of the next time level:
u j ( t 1 ) = e 1 T b , ( u x ) j ( t 1 ) = 2 Δ x e 2 T b .
Gathering (28), (30), (35), (37) and (38), we finally obtain the amplification matrix S in (27):
S = e 1 T e 2 T A ( WB 2 c WC + 3 c ξ ( e i w 1 ) e 1 T ) H 1 P 3 ξ e 3 T ) e 3 T c ( e i w 1 ) e 1 T H 1 P ,
which is the basis for the accuracy and the stability analysis as follows.

3.2. Accuracy Analysis

At the beginning, we assumed AltPoly-3 only had third-order spatial accuracy since all reconstruction polynomials are of the fourth-order, which has third-order spatial accuracy for first-order derivatives. Instead, however, we observed fourth-order spatial accuracy in numerical tests. Analyzing only truncation error is not enough to explain this phenomenon. We believe that the cell-averaged moment u ¯ j and the constrained polynomial reconstruction help improve the spatial accuracy.
The amplification matrix S can be divided into two parts,
S = S 0 + c S 1 ,
where S 0 and S 1 are independent from c. Let u ^ = [ 1 , i w / 2 , 1 i w ( e i w / 2 e i w / 2 ) ] . The principal eigenvalue of S 0 is 1 (see Appendix B for details). We denote the corresponding eigenvector as ψ , namely
S 0 ψ = ψ ,
and we make the last items of ψ and u ^ equal, namely ψ 3 = 1 i w ( e i w / 2 e i w / 2 ) .
With the assistance of Wolfram Mathematica ® 12, we decompose the following terms into Taylor series
ψ = u ^ + ε 1 w 4 e 1 + O ( w 5 )
and
S 1 ψ = i w ψ + ε 2 w 4 e 2 + O ( w 5 ) ,
where the scalars ε 1 , ε 2 , and the vectors v 1 , v 2 are independent with w.
Suppose we choose a proper c that ensures the stability, namely S i M for i N and some positive constant M. Divide the time domain 0 , T uniformly with step length Δ t and denote n = T / Δ t . It is obvious that the exact solution at t = T for (53) with the initial condition u init ( x ) = e i w x / Δ x is u exact = e i w ( x T ) / Δ x = e i n c w u init ( x ) . Hence the error of the numerical solution by AltPoly-3 is
S n u ^ e i n c w u ^ = S n ( u ^ ψ + ψ ) e i n c w ( u ^ ψ + ψ ) S n ψ e i n c w ψ + S n ( u ^ ψ ) + e i n c w ( u ^ ψ ) = S n ψ e i n c w ψ + O ( w 4 ) .
However,
S n ψ = S n 1 ( S 0 + c S 1 ) ψ = S n 1 ( S 0 ψ + c S 1 ψ ) = S n 1 ψ i c w ψ + ε 2 w 4 c e 2 + c O ( w 5 ) = S n 1 ( e i c w ψ + O ( ( c w ) 2 ) + ε 2 w 4 c e 2 + c O ( w 5 ) ) = e i c w S n 1 ψ + S n 1 O ( c 2 w 2 ) + c S n 1 ( ε 2 w 4 e 2 + O ( w 5 ) ) = = e i n c w ψ + k = 1 n e i ( n k ) c w S k 1 O c 2 w 2 + c k = 1 n e i ( n k ) c w S k 1 ε 2 e 2 w 4 + O ( w 5 )
According to the bound result (A16) in Appendix B, S k e 2 μ 1 k + ( k 1 ) c μ 2 k + O ( w ) with some positive constants μ 1 , μ 2 ( 0 , 1 ) . Then one can deduce that
S n ψ e i n c w ψ k = 1 n S k 1 O ( c 2 w 2 ) + c k = 1 n ( μ 1 k + ( k 1 ) c μ 2 k + O ( w ) ) w 4 + O ( w 5 ) n O ( c 2 w 2 ) + c μ 1 1 μ 1 + c μ 2 2 ( 1 μ 2 ) 2 O w 4 + n c O ( w 5 ) .
Remember that w Δ x , c = Δ t / Δ x , and n = T / Δ t , then we have
n O ( c 2 w 2 ) = T Δ t O ( Δ t 2 ) = O ( Δ t ) ,
n c O ( w 5 ) = T Δ t Δ t Δ x O ( Δ x 5 ) = O ( Δ x 4 ) ,
c O ( w 4 ) = Δ t Δ x O ( Δ x 4 ) = O ( Δ t Δ x 3 ) .
Combining (44)–(47), we can see that the error order of AltPoly-3 is
S n u ^ e i n c w u = O ( Δ t + Δ t Δ x 3 + Δ x 4 )
To ensure stability, the order of magnitude of Δ t is usually not higher than that of Δ x , namely c O ( 1 ) . Therefore, the truncated error (48) can be simply reduced to O ( Δ t + Δ x 4 ) . That means AltPoly-3 using the Euler time scheme has fourth-order spatial accuracy and first-order temporal accuracy.
When using RK4 as time integration strategy (20), the amplification matrix S RR 4 for this linear equation can be expressed as
S RK 4 = I + 1 6 K 1 + 2 K 2 + 2 K 3 + K 4 ,
where
K 1 = ( S I ) ,
K 2 = ( S I ) ( I + K 1 / 2 ) ,
K 3 = ( S I ) ( I + K 2 / 2 ) ,
K 4 = ( S I ) ( I + K 3 ) .
Expanding S RK 4 , we have
S RK 4 = 1 24 I + 8 S + 6 S 2 + S 4 .
One can easily show that S RK 4 has fifth-order truncated temporal accuracy and then the total error becomes
S RK 4 n u e i n c w u = O ( Δ t 4 + Δ t Δ x 3 + Δ x 4 ) .

3.3. Stability Analysis

To obtain stable time integration, the CFL number c should be carefully chosen such that the magnitude of the principal eigenvalue of the amplification matrix is not greater than 1 for all wavenumbers k. The principal eigenvalue of the amplification matrix S RK 4 for AltPoly with various moments can be seen in Figure 2. It is remarkable that all schemes have a fairly large range of c except for AltPoly-5. In particular, AltPoly-3, which has fourth-order spatial accuracy, is stable even with c > 1 , which is better than many other existing explicit high-order schemes using a compact stencil. As more moments are used, the stability condition becomes more stringent. For AltPoly-6, the tolerated CFL number is approximately c 0.2 , 0.4 . However, the maximum principal eigenvalue is around 1 + 10 2 when c 0.2 for AltPoly-6, which means short-term time integration can still maintain stability.

4. Numerical Results

In this part, we solve several types of one-dimensional hyperbolic law to illustrate the performance of the proposed method. We choose the four-stage Runge–Kutta method (20) for the time discretization.
For all numerical experiments, AltPoly-r with the moment number r varying from 3 to 6 are tested on uniform meshes. In addition, we check the accuracy of the proposed method on non-uniform meshes. Non-uniform grids are generated by adding a 40% random perturbation upon the element boundary of the uniform mesh. Specifically, suppose x j + 1 / 2 and Δ x are separately the cell boundary and cell length for the uniform grid, then the element boundary generated for the non-uniform grid is
x ˜ j + 1 2 = x j + 1 2 + β ζ j Δ x , ζ j Unif ( 1 , 1 ) ,
where ζ j are independently identically distributed random variables, Unif ( 1 , 1 ) denotes the uniform distribution over [ 1 , 1 ] , and β is the amount of perturbation, which is set as 0.2 in our experiments.
Special attention might be paid to specify the initial values of model variables. This is trivial if the initial profile is analytically given. Otherwise, these unknowns can be specified by implementing high-order interpolations and numerical integration.
Since the temporal accuracy is only fourth-order, we take Δ t ( Δ x ) ( 2 r 2 ) / 4 for AltPoly-r in the accuracy tests.

4.1. Linear Scalar Equation

In this subsection, we test the performance of AltPoly schemes on the following scalar equation [29]:
u t + u x = 0 .
We first check the accuracy of the proposed method through grid refinement tests. The initial condition is given as u ( x , 0 ) = sin ( π x ) with a periodical boundary condition on x [ 0 , 2 π ) . The exact solution is u ( x , t ) = sin ( x t ) . We refine the grid and record numerical errors after one period ( t = 2 ).
Two typical norms for errors are adopted here, i.e., E 1 = 1 N j = 1 N | u j u j true | , and E = max 1 j N | u j u j true | . The errors and the convergence rates are presented in Table 1. We can see that the spatial accuracy order is 2 r 2 for the AltPoly-r scheme, which matches our theoretical analysis in Section 3. Additionally, the order of accuracy is also maintained well for non-uniform meshes.
Next, we use a complicated initial condition from [30] to evaluate the ability of the proposed method in handling both continuous and discontinuous profiles. This initial condition contains four piece-wise functions:
u ( x , 0 ) = 1 6 ( G ( x , β , z δ ) + G ( x , β , z + δ ) + 4 G ( x , β , z ) ) for x [ 0.8 , 0.6 ] , 1 for x [ 0.4 , 0.2 ] , 1 | 10 ( x 0.1 ) | , for x [ 0 , 0.2 ] , 1 6 ( F ( x , α , α δ ) + F ( x , α , α + δ ) + 4 F ( x , α , a ) ) for x [ 0.4 , 0.6 ] , 0 , otherwise ,
where the involved functions are
G ( x , β , z ) = exp ( β ( x z ) 2 ) )
and
F ( x , α , a ) = max ( 1 α 2 ( x a ) 2 , 0 ) ,
with the constants set as α = 10 , a = 0.5 , z = 0.7 , δ = 0.005 , and β = log 2 / ( 36 δ 2 ) . As for the initial profile, we suggest using the routine used in the discontinuous Galerkin method.
A smaller time step can help reduce the oscillations around discontinuities. We take Δ t = 0.3 Δ x , 0.1 Δ x , 0.075 Δ x , and 0.01 Δ x , respectively, for AltPoly-3, AltPoly-4, AltPoly-5, and AltPoly-6.
The computational domain is divided into 200 elements. We first compute this problem without adding artificial viscosity, for which higher-order methods can produce better results in both continuous and discontinuous regions. The results at t = 2 are illustrated in Figure 3 and Figure 4 separately. Specifically, AltPoly-3 and AltPoly-4 generate significant oscillations around x = 0.4 and x = 0.2 . These oscillations are reasonable as our method is based on the smoothness assumption and the polynomial does have its limitations on presenting discontinuities. However, it is interesting that spurious oscillations are significantly reduced and the transitions around strong discontinuities are quite sharp for AltPoly-5 and AltPoly-6. After adding a proper amount of artificial viscosity, the proposed method is able to smooth the artificial oscillations for AltPoly method of various orders while not losing high accuracy in the smooth regions.

4.2. Nonlinear Burgers Equation

This subsection considers the 1D inviscid Burgers equation
u t + u 2 2 x = 0 ,
The initial condition is u ( x , 0 ) = 0.5 + sin ( x ) for x [ 0 , 2 π ) with the periodical boundary condition.
Although the initial profile is smooth, shock waves exist when t is large enough since the flux function is nonlinear. To evaluate the order of accuracy, we calculate the numerical solution before the shock wave occurs and compute the errors in respect of grid refinement. Table 2 lists the errors at t = 0.5 when the solution is still smooth. The numerical convergence rate results are consistent with our theoretical analysis.
The numerical results for the Burgers equation at t = 1.5 are shown in Figure 5 when a shock has already emerged. To handle the shock discontinuity, we add a proper amount of artificial viscosity according to Section 2.4. We can see that our schemes suppress the oscillations well around the discontinuity while still maintaining a sharp and accurate profile.

4.3. Euler Equations

We then use our scheme to solve the one-dimensional Euler system:
t ρ m e + x m ρ u 2 + p e u + p u = 0 ,
where ρ refers to the density, u the velocity, m = ρ u the momentum, e the total energy, and p the pressure:
p = ( γ 1 ) e 1 2 ρ u 2
with γ = 1.4 for a perfect dry gas. Just like in the scalar case, we introduce point-valued and cell-averaged moments for each variable. The cell-averaged moments are updated by the flux-form Equation (56). In Step 2 of AltPoly, solution values on auxiliary points are updated using the following non-conservative form equations:
ρ t = m x ,
m t = ( m x u + m u x + p x ) ,
e t = ( e x u + e u x + p x u + p u x ) ,
where
u x = x m ρ = m x ρ m ρ x ρ 2 , p x = ( γ 1 ) e x 1 2 ( m x u + m u x ) .
Neither characteristic decomposition nor Riemann solver are needed here.
We first check the order of accuracy for this nonlinear system by simulating advection of the density perturbation. The initial conditions are ρ ( x 0 ) = 1 + 0.2 sin x , v ( x , 0 ) = 1 , and p ( x , 0 ) = 1 for x [ 0 , 2 π ) with the periodic boundary condition. The exact solution is a traveling wave: ρ ( x , t ) = 1 + 0.2 sin ( x t ) , v ( x , t ) = 1 and p ( x , t ) = 1 . The error results of the density are listed in Table 3 and a satisfactory convergence rate is obtained.
We then consider the classical shock tube problem, which is described by Euler equations with the following Riemann-type initial distribution:
ρ , v , p = ρ L , v L , p L , 0 x 0.5 ; ρ R , v R , p R , 0.5 < x 1 , x [ 0 , 1 ] .
Two initial conditions are tested here. The first one is known as Sod’s problem and the initial condition is defined as
ρ L , v L , p L = ( 1 , 0 , 1 ) , ρ R , v R , p R = ( 0.125 , 0 , 0.1 ) .
The second one is Lax’s problem, which starts from the following initial condition:
ρ L , v L , p L = ( 0.445 , 0.698 , 3.258 ) , ρ R , v R , p R = ( 0.5 , 0 , 0.571 ) .
We compute these shock tube problems with 100 cells using the AltPoly schemes equipped with the artificial viscosity. Generally, a certain degree of spurious oscillation still exists but is in a tolerant range. As shown in Figure 6 and Figure 7, all AltPoly schemes can depict shocks sharply using only one or two cells. The higher-order schemes (AltPoly-5 and AltPoly-6) capture the contact discontinuities with only one or two cells, while more cells are needed for the lower-order AltPoly schemes.
To show the accuracy of our smoothness indicator (23), Figure 8 illustrates the cells where the artificial viscosity takes effect on the x-t plane. Generally, the number of cells recognized as discontinuous for lower-order schemes is significantly smaller than that for higher-order schemes. This might be because lower-order schemes have a larger dissipation and the smoothness indicator is less accurate than that of high-order schemes. It is remarkable that the artificial viscosity vanishes in most areas and is in effect only around the shock wave. The contact discontinuity does not trigger the smoothness indicator as it can be smoothed by the artificial viscosity in the first few time steps.

5. Conclusions

We have presented a high-order conservative solver termed AltPoly for 1D hyperbolic conservation laws. This method can be categorized as a multi-moment scheme since it adopts both point-values at the middle of the cell and cell-averaged values as the model variables. AltPoly updates the model variables via alternately implementing two polynomial reconstructions. The point-based variables are used to reconstruct a high-order accurate approximation of the solution via Hermite interpolation. The cell-averaged variables updated by the flux-form equations serve as the constraint to guarantee numerical conservation. The solution within a cell is corrected by polynomial reconstruction via solving a constrained least-squares problem. In other words, our method directly reconstructs the solution to maintain the conservation, which is different from existing schemes based on the local flux reconstruction. Fourier analysis shows the accuracy and the admissible range of the CFL number for the linear scalar equation. Further, our scheme can be extended to other kinds of PDEs such as parabolic equations, advection–diffusion equations, etc., which will be our future work.

Author Contributions

S.L., conceptualization, methodology, software, validation, writing—original draft preparation; Q.L., methodology, formal analysis, validation, writing—review and editing; H.L., supervision, methodology, writing—review and editing; J.S., project administration, resources. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by National Natural Science Foundation of China No. 41605070.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Outline of AltPoly-r

Steps 1 and 2 of AltPoly-r are almost the same as those of AltPoly-3 and are hence omitted here. In this part, we focus on introducing the constrained polynomial reconstruction in Step 3 of AltPoly-r for r 3 .
Following the notations in Section 2.2.3, the core of Step 3 is reconstructing the polynomial l ˜ i ( ξ ) = i = 0 2 r 3 b i ξ i from the following conditions:
l ˜ j ( ξ k ) = i = 0 2 r 3 b i ξ k i = u j , k ( t 1 ) , for k = 1 , , 2 r 2 ,
1 2 1 1 l ˜ j ( ξ ) d ξ = b 0 + 1 3 b 2 + + 1 2 r 3 b 2 r 4 = u ¯ j ( t 1 ) ,
where u j , k t 1 and u ¯ j ( t 1 ) are obtained in Step 2. Since the hyperbolic conservation law must be obeyed, we rewrite the constraint condition (A2) into the expression of b 2 r 4 ,
b 2 r 4 = ( 2 r 3 ) u ¯ j ( t 1 ) b 0 1 3 b 2 1 2 r 5 b 2 r 6 ,
and plug it into (A1) to arrive at
b 0 + b 1 ξ k + + b 2 r 5 ξ k 2 r 5 + ( 2 r 3 ) ( u ¯ j t 1 b 0 1 3 b 2 1 2 r 5 b 2 r 6 ) ξ k 2 r 4 + b 2 r 3 ξ k 2 r 3 = u j , k ( t 1 ) .
We rewrite it in matrix form:
Ab = y ,
where
A = A 1 A 2 ,
A 1 = 1 ξ 1 ξ 1 2 ξ 1 2 r 4 ξ 1 2 r 3 1 ξ 2 ξ 2 2 ξ 2 2 r 4 ξ 2 2 r 3 1 ξ 2 r 2 ξ 2 r 2 2 ξ 2 r 2 2 r 4 ξ 2 r 2 2 r 3 ,
A 2 = ( 2 r 3 ) ξ 1 2 r 4 0 1 3 ξ 1 2 r 4 0 1 5 ξ 1 2 r 4 0 1 2 r 5 ξ 1 2 r 4 0 0 ξ 2 2 r 4 0 1 3 ξ 2 2 r 4 0 1 5 ξ 2 2 r 4 0 1 2 r 5 ξ 2 2 r 4 0 0 ξ 2 r 2 2 r 4 0 1 3 ξ 2 r 2 2 r 4 0 1 5 ξ 2 r 2 2 r 4 0 1 2 r 5 ξ 2 r 2 2 r 4 0 0 ,
b = b 0 b 1 b 2 r 5 b 2 r 3 , y = u j , 1 t 1 u j , 2 t 1 u j , 2 r 2 t 1 ( 2 r 3 ) ξ 1 2 r 4 ξ 2 2 r 4 ξ 2 r 2 2 r 4 .
This linear system is overdetermined since it has 2 r 3 unknowns and 2 r 2 equations. The least-squares solution is
b = A y = ( A T A ) 1 A T y ,
where A denotes the generalized inverse or Moore–Penrose inverse [27] of A . From b we can obtain the updated point-based model variables.
In practice, we find that A has a large condition number when r 5 and we calculate it via symbolic computation using mathematical software to reduce the error. Throughout this paper, we used the symmetrical solution points 1 , 1 + σ 1 , 1 + σ 2 , , 1 + σ r 2 , 1 σ r 2 , 1 σ r 1 , , 1 σ 1 , 1 , for AltPoly-r with the choices of σ i shown in Table A1. Through numerical experiments, we found that smaller σ i can generally induce a more stable numerical scheme, which allows for a larger CFL number. However, when σ i approaches 0, the condition number of A tends to infinity, which will cause a loss of accuracy. Based on numerical trials, the choices of σ i are determined as a compromise between stability and accuracy.
Table A1. Choice of σ i for solution points.
Table A1. Choice of σ i for solution points.
rSolution Points
3 σ 1 = 0.03
4 σ 1 = 0.05 , σ 2 = 0.1
5 σ 1 = 1 34 , σ 2 = 2 11 , σ 3 = 1 4
6 σ 1 = 0.02 , σ 2 = 0.05 , σ 3 = 0.1 , σ 4 = 0.12

Appendix B. Some Technical Results Used in Section 3

Via symbolic computation by Mathematica, we have the following Taylor expansions in decimal form:
S 0 e 1 = ( 0.267 + 0.256 cos w ) e 1 + 0.379 sin w e 2 = 0.523 e 1 + w v 1 + O ( w 2 )
S 0 e 2 = 0.129 sin w e 1 + ( 0.508 0.246 cos w ) e 2 + = 0.129 w e 1 + 0.262 e 2 + O ( w 2 )
S 1 e 1 = 0.739 i sin w e 1 + 4.523 sin 2 ( w 2 ) e 1 i sin w e 3 = w v 2 + O ( w 2 )
S 1 e 2 = 0.521 e 1 + 1.261 i w e 2 + O ( w 2 )
where the vectors v 3 and v 4 are independent with w.
We can see that the magnitudes of all entries in S 0 [ 1 : 2 ; 1 : 2 ] are smaller than 1 when w is small enough. In addition, from (39) we know that
S 0 e 3 = e 3 .
Therefore the principal eigenvalue of S 0 is 1 as long as w is small enough.
Now we try to bound S n e 1 and S n e 2 . We first assume that S i M for any i 1 , which naturally holds as long as the scheme is stable. Then we have
S i e 1 = S i 1 ( ( S 0 + c S 1 ) e 1 ) = S i 1 ( 0.523 e 1 + w v + O ( w 2 ) ) = = ( 0.523 ) n e 1 + k = 1 i ( 0.523 ) n i S k 1 ( w v 1 + O ( w 2 ) ) 0 . 523 i + M ( w v 1 + O ( w 2 ) ) = 0 . 523 i + O ( w ) .
Hence,
S i e 2 = S i 1 0.262 e 2 + O ( w ) + c ( 0.521 ) e 1 + O ( w ) = S i 2 0 . 262 2 e 2 + ( 0.262 + S ) ( O ( w ) + c ( 0.521 e 1 + O ( w ) ) ) = = 0 . 262 i e 2 + k = 1 n ( 0 . 262 k 1 S i k ) O ( w ) + c O ( w ) A 1 0.521 c k = 1 i 0 . 262 k 1 S i k e 1 A 2
Since the magnitude of c = Δ t / Δ x is not greater than O ( 1 ) and S i M , then A 1 can be estimated by
A 1 k = 0 i 1 0 . 262 k 1 M O ( w ) M O ( w ) 1 0.262 = O ( w ) .
As for A 2 , the following bound can be deduced based on (A12),
A 2 k = 1 i 1 0 . 262 k 1 S i k e 1 k = 1 i 1 0 . 262 k 1 ( 0 . 523 i k + O ( w ) ) ( i 1 ) 0 . 523 i 1 + O ( w ) .
Combining (A14) and (A15), S i e 2 can be bounded by
S i e 2 0 . 262 i + ( i 1 ) 0 . 523 i c + O ( w ) .

References

  1. Lele, S.K. Compact finite difference schemes with spectral-like resolution. J. Comput. Phys. 1992, 103, 16–42. [Google Scholar] [CrossRef]
  2. Mohebbi, A.; Dehghan, M. High order compact solution of the one-space-dimensional linear hyperbolic equation. Numer. Methods Part. Differ. Equ. Int. J. 2008, 24, 1222–1235. [Google Scholar] [CrossRef]
  3. Yabe, T.; Ishikawa, T.; Wang, P.; Aoki, T.; Kadota, Y.k.; Ikeda, F. A universal solver for hyperbolic equations by cubic-polynomial interpolation II. Two-and three-dimensional solvers. Comput. Phys. Commun. 1991, 66, 233–242. [Google Scholar] [CrossRef]
  4. Yabe, T.; Xiao, F.; Utsumi, T. The constrained interpolation profile method for multiphase analysis. J. Comput. Phys. 2001, 169, 556–593. [Google Scholar] [CrossRef]
  5. Aoki, T. Interpolated differential operator (IDO) scheme for solving partial differential equations. Comput. Phys. Commun. 1997, 102, 132–146. [Google Scholar] [CrossRef]
  6. Imai, Y.; Aoki, T. Stable coupling between vector and scalar variables for the IDO scheme on collocated grids. J. Comput. Phys. 2006, 215, 81–97. [Google Scholar] [CrossRef]
  7. Goodrich, J.; Hagstrom, T.; Lorenz, J. Hermite methods for hyperbolic initial-boundary value problems. Math. Comput. 2006, 75, 595–630. [Google Scholar] [CrossRef] [Green Version]
  8. Tanaka, R.; Nakamura, T.; Yabe, T. Constructing exactly conservative scheme in a non-conservative form. Comput. Phys. Commun. 2000, 126, 232–243. [Google Scholar] [CrossRef]
  9. Yabe, T.; Tanaka, R.; Nakamura, T.; Xiao, F. An exactly conservative semi-Lagrangian scheme (CIP–CSL) in one dimension. Mon. Weather Rev. 2001, 129, 332–344. [Google Scholar] [CrossRef]
  10. Imai, Y.; Aoki, T.; Takizawa, K. Conservative form of interpolated differential operator scheme for compressible and incompressible fluid dynamics. J. Comput. Phys. 2008, 227, 2263–2285. [Google Scholar] [CrossRef]
  11. Onodera, N.; Aoki, T.; Kobayashi, H. Large-eddy simulation of turbulent channel flows with conservative IDO scheme. J. Comput. Phys. 2011, 230, 5787–5805. [Google Scholar] [CrossRef]
  12. Ii, S.; Shimuta, M.; Xiao, F. A 4th-order and single-cell-based advection scheme on unstructured grids using multi-moments. Comput. Phys. Commun. 2005, 173, 17–33. [Google Scholar] [CrossRef]
  13. Cockburn, B.; Shu, C.W. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework. Math. Comput. 1989, 52, 411–435. [Google Scholar]
  14. Cockburn, B.; Lin, S.Y.; Shu, C.W. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One-dimensional systems. J. Comput. Phys. 1989, 84, 90–113. [Google Scholar] [CrossRef] [Green Version]
  15. Qiu, J.; Shu, C.W. Runge–Kutta discontinuous Galerkin method using WENO limiters. SIAM J. Sci. Comput. 2005, 26, 907–929. [Google Scholar] [CrossRef]
  16. Wang, Z.J. Spectral (finite) volume method for conservation laws on unstructured grids. basic formulation: Basic formulation. J. Comput. Phys. 2002, 178, 210–251. [Google Scholar] [CrossRef] [Green Version]
  17. Wang, Z.J.; Zhang, L.; Liu, Y. Spectral (finite) volume method for conservation laws on unstructured grids IV: Extension to two-dimensional systems. J. Comput. Phys. 2004, 194, 716–741. [Google Scholar] [CrossRef]
  18. Liu, Y.; Vinokur, M.; Wang, Z.J. Spectral difference method for unstructured grids I: Basic formulation. J. Comput. Phys. 2006, 216, 780–801. [Google Scholar] [CrossRef]
  19. Huynh, H.T. A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods. In Proceedings of the 18th AIAA Computational Fluid Dynamics Conference, Miami, FL, USA, 25–28 June 2007; p. 4079. [Google Scholar]
  20. Wang, Z.J.; Gao, H. A unifying lifting collocation penalty formulation including the discontinuous Galerkin, spectral volume/difference methods for conservation laws on mixed grids. J. Comput. Phys. 2009, 228, 8161–8186. [Google Scholar] [CrossRef]
  21. Ii, S.; Xiao, F. High order multi-moment constrained finite volume method. Part I: Basic formulation. J. Comput. Phys. 2009, 228, 3669–3707. [Google Scholar] [CrossRef] [Green Version]
  22. Xiao, F.; Ii, S.; Chen, C.; Li, X. A note on the general multi-moment constrained flux reconstruction formulation for high order schemes. Appl. Math. Model. 2013, 37, 5092–5108. [Google Scholar] [CrossRef]
  23. Kirby, R.M.; Karniadakis, G.E. De-aliasing on non-uniform grids: Algorithms and applications. J. Comput. Phys. 2003, 191, 249–264. [Google Scholar] [CrossRef] [Green Version]
  24. Kanevsky, A.; Carpenter, M.H.; Hesthaven, J.S. Idempotent filtering in spectral and spectral element methods. J. Comput. Phys. 2006, 220, 41–58. [Google Scholar] [CrossRef]
  25. Zienkiewicz, O.C.; Taylor, R.L.; Zhu, J.Z. The Finite Element Method: Its Basis and Fundamentals; Elsevier: Amsterdam, The Netherlands, 2005. [Google Scholar]
  26. Staniforth, A.; Côté, J. Semi-Lagrangian integration schemes for atmospheric models—A review. Mon. Weather Rev. 1991, 119, 2206–2223. [Google Scholar] [CrossRef] [Green Version]
  27. Ben-Israel, A.; Greville, T.N. Generalized Inverses: Theory and Applications; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2003; Volume 15. [Google Scholar]
  28. Persson, P.O.; Peraire, J. Sub-cell shock capturing for discontinuous Galerkin methods. In Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, USA, 9–12 January 2006; p. 112. [Google Scholar]
  29. Dehghan, M.; Shokri, A. A meshless method for numerical solution of a linear hyperbolic equation with variable coefficients in two space dimensions. Numer. Methods Part. Differ. Equ. Int. J. 2009, 25, 494–506. [Google Scholar] [CrossRef]
  30. Jiang, G.S.; Shu, C.W. Efficient implementation of weighted ENO schemes. J. Comput. Phys. 1996, 126, 202–228. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Algorithm illustration.
Figure 1. Algorithm illustration.
Mathematics 09 01885 g001
Figure 2. Contour plots of the principal eigenvalues of the amplification matrices of AltPoly schemes using fourth-order RK time integration for the linear advection equation.
Figure 2. Contour plots of the principal eigenvalues of the amplification matrices of AltPoly schemes using fourth-order RK time integration for the linear advection equation.
Mathematics 09 01885 g002aMathematics 09 01885 g002b
Figure 3. Numerical results of the advection of a complicated profile at t = 2 via AltPoly without artificial viscosity using N = 200 cells.
Figure 3. Numerical results of the advection of a complicated profile at t = 2 via AltPoly without artificial viscosity using N = 200 cells.
Mathematics 09 01885 g003
Figure 4. Numerical results of the advection of a complicated profile at t = 2 via AltPoly with artificial viscosity using N = 200 cells.
Figure 4. Numerical results of the advection of a complicated profile at t = 2 via AltPoly with artificial viscosity using N = 200 cells.
Mathematics 09 01885 g004
Figure 5. Numerical results of the Burgers equation at t = 1.5 with N = 80 cells and Δ t = 0.15 Δ x .
Figure 5. Numerical results of the Burgers equation at t = 1.5 with N = 80 cells and Δ t = 0.15 Δ x .
Mathematics 09 01885 g005aMathematics 09 01885 g005b
Figure 6. Numerical results of Sod’s problem at t = 0.13 with N = 100 cells. Δ t = 0.2 Δ x , 0.15 Δ x , 0.1 Δ x , and 0.05 Δ x , respectively.
Figure 6. Numerical results of Sod’s problem at t = 0.13 with N = 100 cells. Δ t = 0.2 Δ x , 0.15 Δ x , 0.1 Δ x , and 0.05 Δ x , respectively.
Mathematics 09 01885 g006
Figure 7. Numerical results of Lax’s problem at t = 0.13 with N = 100 cells. Δ t = 0.1 Δ x , 0.1 Δ x , 0.08 Δ x , and 0.05 Δ x , respectively.
Figure 7. Numerical results of Lax’s problem at t = 0.13 with N = 100 cells. Δ t = 0.1 Δ x , 0.1 Δ x , 0.08 Δ x , and 0.05 Δ x , respectively.
Mathematics 09 01885 g007aMathematics 09 01885 g007b
Figure 8. The cells where artificial viscosity takes effect in Lax’s problem.
Figure 8. The cells where artificial viscosity takes effect in Lax’s problem.
Mathematics 09 01885 g008
Table 1. Accuracy test for linear hyperbolic law with initial condition u ( x , 0 ) = sin ( x ) at t = 2 π .
Table 1. Accuracy test for linear hyperbolic law with initial condition u ( x , 0 ) = sin ( x ) at t = 2 π .
rN L 1 ErrorOrder L ErrorOrder L 1 ErrorOrder L ErrorOrder
Uniform GridNon-Uniform Grid
3101.79 × 10 3 2.76 × 10 3 2.10 × 10 3 3.33 × 10 3
201.15 × 10 4 3.951.82 × 10 4 3.931.32 × 10 4 4.002.09 × 10 4 3.99
302.31 × 10 5 3.963.63 × 10 5 3.972.68 × 10 5 3.934.25 × 10 5 3.93
407.35 × 10 6 3.981.16 × 10 5 3.978.46 × 10 6 4.001.35 × 10 5 3.99
601.46 × 10 6 3.992.29 × 10 6 3.991.68 × 10 6 3.982.67 × 10 6 3.99
4101.32 × 10 4 2.05 × 10 4 1.33 × 10 4 2.06 × 10 4
202.06 × 10 6 6.013.22 × 10 6 5.992.06 × 10 6 6.013.23 × 10 6 5.99
301.81 × 10 7 5.992.84 × 10 7 5.981.82 × 10 7 5.992.86 × 10 7 5.98
403.23 × 10 8 6.005.06 × 10 8 6.003.24 × 10 8 6.005.08 × 10 8 6.00
602.84 × 10 9 6.004.45 × 10 9 5.992.85 × 10 9 5.994.48 × 10 9 5.99
5101.92 × 10 4 2.97 × 10 4 1.92 × 10 4 2.98 × 10 4
207.60 × 10 7 7.981.18 × 10 6 7.977.60 × 10 7 7.981.18 × 10 6 7.98
302.95 × 10 8 8.014.63 × 10 8 8.002.95 × 10 8 8.014.63 × 10 8 7.99
402.97 × 10 9 7.994.64 × 10 9 7.992.96 × 10 9 7.994.65 × 10 9 7.99
601.16 × 10 10 8.001.81 × 10 10 8.001.16 × 10 10 8.001.82 × 10 10 8.00
6201.57 × 10 6 2.46 × 10 6 1.57 × 10 6 2.46 × 10 6
302.73 × 10 8 10.004.28 × 10 8 9.992.73 × 10 8 9.994.28 × 10 8 9.99
401.54 × 10 9 9.992.41 × 10 9 9.991.54 × 10 9 9.992.41 × 10 9 9.99
602.49 × 10 11 10.173.96 × 10 11 10.142.28 × 10 11 10.393.82 × 10 11 10.22
Table 2. Accuracy test for Burgers equation with u ( x , 0 ) = 0.5 + sin ( x ) at t = 0.5 .
Table 2. Accuracy test for Burgers equation with u ( x , 0 ) = 0.5 + sin ( x ) at t = 0.5 .
rN L 1 ErrorOrder L ErrorOrder L 1 ErrorOrder L ErrorOrder
Uniform GridNon-Uniform Grid
3208.41 × 10 6 1.32 × 10 5 2.01 × 10 5 5.46 × 10 5
301.59 × 10 6 4.102.50 × 10 6 4.112.99 × 10 6 4.709.14 × 10 6 4.41
404.93 × 10 7 4.087.71 × 10 7 4.088.52 × 10 7 4.372.54 × 10 6 4.45
609.52 × 10 8 4.051.50 × 10 7 4.041.90 × 10 7 3.705.71 × 10 7 3.68
802.98 × 10 8 4.034.69 × 10 8 4.035.39 × 10 8 4.381.74 × 10 7 4.12
4201.98 × 10 5 1.21 × 10 4 2.11 × 10 5 1.29 × 10 4
302.18 × 10 6 5.431.83 × 10 5 4.662.13 × 10 6 5.651.73 × 10 5 4.96
404.24 × 10 7 5.703.60 × 10 6 5.654.16 × 10 7 5.683.53 × 10 6 5.52
603.58 × 10 8 6.092.94 × 10 7 6.173.65 × 10 8 6.002.99 × 10 7 6.09
806.38 × 10 9 5.995.41 × 10 8 5.896.54 × 10 9 5.975.28 × 10 8 6.03
5208.48 × 10 5 5.80 × 10 4 9.04 × 10 5 5.55 × 10 4
304.07 × 10 6 7.493.04 × 10 5 7.283.90 × 10 6 7.753.01 × 10 5 7.19
404.41 × 10 7 7.733.40 × 10 6 7.614.39 × 10 7 7.603.38 × 10 6 7.60
601.83 × 10 8 7.841.42 × 10 7 7.831.84 × 10 8 7.821.43 × 10 7 7.80
801.81 × 10 9 8.041.49 × 10 8 7.851.81 × 10 9 8.061.48 × 10 8 7.87
6204.54 × 10 4 2.74 × 10 3 4.68 × 10 4 2.73 × 10 3
301.30 × 10 5 8.771.00 × 10 4 8.151.29 × 10 5 8.861.00 × 10 4 8.15
407.54 × 10 7 9.895.87 × 10 6 9.877.62 × 10 7 9.835.85 × 10 6 9.87
601.32 × 10 8 9.981.02 × 10 7 9.991.32 × 10 8 10.001.03 × 10 7 9.97
807.38 × 10 10 10.036.09 × 10 9 9.817.67 × 10 10 9.906.09 × 10 9 9.82
Table 3. Accuracy test for the one-dimensional Euler equation with u ( x , 0 ) = 0.5 + sin ( x ) at t = 2 π .
Table 3. Accuracy test for the one-dimensional Euler equation with u ( x , 0 ) = 0.5 + sin ( x ) at t = 2 π .
rN L 1 ErrorOrder L ErrorOrder L 1 ErrorOrder L ErrorOrder
Uniform GridNon-Uniform Grid
3103.26 × 10 4 5.03 × 10 4 3.88 × 10 4 6.22 × 10 4
156.09 × 10 5 4.139.55 × 10 5 4.107.22 × 10 5 4.151.16 × 10 4 4.15
201.84 × 10 5 4.162.85 × 10 5 4.212.21 × 10 5 4.123.50 × 10 5 4.16
303.45 × 10 6 4.135.41 × 10 6 4.094.32 × 10 6 4.036.93 × 10 6 3.99
401.07 × 10 6 4.091.67 × 10 6 4.091.30 × 10 6 4.182.09 × 10 6 4.17
4101.12 × 10 5 1.72 × 10 5 1.20 × 10 5 2.40 × 10 5
159.81 × 10 7 6.001.54 × 10 6 5.969.91 × 10 7 6.161.55 × 10 6 6.75
201.76 × 10 7 5.972.76 × 10 7 5.981.78 × 10 7 5.972.79 × 10 7 5.97
301.55 × 10 8 5.992.43 × 10 8 5.991.57 × 10 8 5.992.47 × 10 8 5.98
402.77 × 10 9 5.994.35 × 10 9 5.982.81 × 10 9 5.984.42 × 10 9 5.98
5102.55 × 10 7 3.94 × 10 7 2.55 × 10 7 3.96 × 10 7
151.02 × 10 8 7.951.59 × 10 8 7.911.02 × 10 8 7.951.59 × 10 8 7.92
201.02 × 10 9 7.991.58 × 10 9 8.031.02 × 10 9 7.981.58 × 10 9 8.03
303.96 × 10 11 8.016.22 × 10 11 7.983.97 × 10 11 8.016.24 × 10 11 7.98
403.98 × 10 12 7.986.49 × 10 12 7.853.99 × 10 12 7.996.47 × 10 12 7.88
6105.08 × 10 7 7.85 × 10 7 5.08 × 10 7 7.85 × 10 7
128.35 × 10 8 9.901.27 × 10 7 9.988.35 × 10 8 9.901.27 × 10 7 9.98
158.90 × 10 9 10.031.40 × 10 8 9.918.91 × 10 9 10.031.40 × 10 8 9.90
205.03 × 10 10 9.997.90 × 10 10 9.981.44 × 10 9 10.012.25 × 10 9 10.02
248.34 × 10 11 9.861.33 × 10 10 9.775.02 × 10 10 9.977.92 × 10 10 9.90
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lin, S.; Luo, Q.; Leng, H.; Song, J. Alternating Polynomial Reconstruction Method for Hyperbolic Conservation Laws. Mathematics 2021, 9, 1885. https://doi.org/10.3390/math9161885

AMA Style

Lin S, Luo Q, Leng H, Song J. Alternating Polynomial Reconstruction Method for Hyperbolic Conservation Laws. Mathematics. 2021; 9(16):1885. https://doi.org/10.3390/math9161885

Chicago/Turabian Style

Lin, Shijian, Qi Luo, Hongze Leng, and Junqiang Song. 2021. "Alternating Polynomial Reconstruction Method for Hyperbolic Conservation Laws" Mathematics 9, no. 16: 1885. https://doi.org/10.3390/math9161885

APA Style

Lin, S., Luo, Q., Leng, H., & Song, J. (2021). Alternating Polynomial Reconstruction Method for Hyperbolic Conservation Laws. Mathematics, 9(16), 1885. https://doi.org/10.3390/math9161885

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop