Reconstructed Interpolating Differential Operator Method with Arbitrary Order of Accuracy for the Hyperbolic Equation

: We propose a family of multi-moment methods with arbitrary orders of accuracy for the hyperbolic equation via the reconstructed interpolating differential operator (RDO) approach. Reconstruction up to arbitrary order can be achieved on a single cell from properly allocated model variables including spatial derivatives of varying orders. Then we calculate the temporal derivatives of coefﬁcients of the reconstructed polynomial and transform them into the temporal derivatives of the model variables. Unlike the conventional multi-moment methods which evolve different types of moments by deriving different equations, RDO can update all derivatives uniformly via a simple linear transform more efﬁciently. Based on difference in introducing interaction from adjacent cells, the central RDO and the upwind RDO are proposed. Both schemes enjoy high-order accuracy which is veriﬁed by Fourier analysis and numerical experiments.


Introduction
In the last few decades, high-order numerical methods are becoming increasingly popular in the research community due to the advantages over the low-order methods in achieving higher resolution more efficiently with lower computational cost [1]. The conventional high-order approaches based on high-order finite difference (FD) and the finite volume (FV) methods only use one degree of freedom (DOF) on each node or cell. For these methods, a high-order approximation requires a wide stencil which introduces huge communication costs between computational nodes on supercomputers.
To narrow the stencil, a natural and feasible approach is to include two or more DOFs per cell so as to construct a higher-order polynomial. The constrained interpolation profile (CIP) scheme [2,3] evolves two moments, i.e., the physical variable and the spatial derivative simultaneously and independently according to the governing equations in different forms. Specifically, CIP describes the spatial profile using the third-order Hermite interpolation and updates the profile according to the local analytic solution in via the semi-Lagrangian approach. Following this pioneering work, the so-called CIP-conservative semi-Lagrangian (CIP-CSL) schemes [4,5] add the cell-averaged value as an extra moment to ensure the conservativeness for the scalar conservative advection transport. Successive studies have resulted in a more general framework, the so-called CIP/multi-moment finite volume method (CIP/MM FVM) which has been applied to various fluid dynamic simulations [6][7][8]. The interpolated differential operator (IDO) method [9] is similar to CIP but IDO updates the model variables using the Eulerian method instead of the semi-Lagragian method. A conservative variant of IDO scheme also predicts the cell-averaged moment separately [10]. Both CIP and IDO require an additional equation for the spatial derivative moment which can be difficult to obtain when treating nonlinear and highdimensional problems.
A compact CIP/MM-FVM formulation that uses high-order derivative moments can reach arbitrary accuracy order has been presented in [11]. This approach predicts the highorder derivative moments in the Eulerian representation via solving a linearized high-order derivative Riemann derivative problem [12]. As a more flexible variant of CIP/MM-FVM, the multi-moment constrained finite volume (MCV) [13] treats the cell-averaged moments as constraints instead of explicit model variables. Following this principle, the MCV method has been generalized to the unstructured meshes [14,15] recently.
High-order schemes based on local flux reconstruction are another representative class of high-order schemes which also use compact stencil for spatial discretization and are increasingly attractive. Examples includes the discontinuous Galerkin(DG) [16], the spectral volume [17] and the flux reconstruction [18] method.
This article presents a new efficient compact multi-moment method termed reconstructed interpolating differential operator (RDO) method which does not involve the cross-cell reconstruction for the hyperbolic equations. This work can be viewed as an generalization of the original IDO scheme [9] to arbitrary orders, non-uniform meshes and higher dimensions. Moreover, this work proposes a novel updating rules for model variables via a direct linear transform, which does not require introducing additional equations for different types of moments. The core idea is simple and direct. In each cell, the solution is approximated by a polynomial to interpolate the moments including the point values and high-order derivatives at the cell boundaries. The next step computes the temporal derivatives of the polynomial coefficients with the aid of solution points. Then the temporal derivatives are retrieved from the temporal derivatives of polynomial coefficient via a linear transform. To deal with the non-uniqueness of temporal derivatives computed form different cells at cell boundaries, the central and the upwind RDO are proposed. Fourier analysis and numerical tests indicate that our scheme can achieve arbitrary orders of accuracy.
The remaining part of this paper is organized as follows. Section 2 gives the basic formulation of RDO in one dimension. Section 3 describes the generalization of RDO to higher dimensions. Then Section 4 investigates the accuracy order and stability by classical Fourier analysis. Numerical tests in Section 5 validate the Fourier analysis results and compare the computational efficiency and numerical performance of RDO of various orders. Section 6 ends this paper with a few conclusions and discussions.

Basic Formulation for the One-Dimensional Problem
Consider the 1D scalar hyperbolic equation ∂u ∂t where x is the spatial coordinate, t denotes time, u is the conserved variable and f (u) is the flux.
To begin with, we introduce some notations here. The bold lower-case letter such as a, u stands for a column vector. The bold capital letter like A denotes a matrix. The column vector [a 1 , . . . , a m ] T is denoted by [a k ] m k=1 for simplicity. A T denotes the transpose of A. The p norm of the vector a is denoted by a p and a is a simplified form of a 2 throughout this paper. The matrix norm of H is defined by H = max z =1 Hz / z .

Spatial Discretization
We divide the computational domain [x min , x max ] into N non-overlapping cells or elements, among which the j-th cell is I j := [x j−1/2 , x j+1/2 ]. For convenience, we denote the middle point of I j by x j = (x j−1/2 + x j+1/2 )/2 and the length of I j by ∆x j = x j+1/2 − x j−1/2 .
Model variables including point-based value and derivatives are independently evolved in our scheme. For simplicity, the RDO scheme that have M model variables in each cell is denoted by RDO-M. The allocation of model variables has a slight difference between M = 2K and M = 2K + 1. For RDO-2K, 2K model variables in each cell are all located at the cell boundaries. These model variables are {u j±1/2 , (u x ) j±1/2 , . . . , (u For RDO-(2K + 1), we add the function value at the middle of the cell, i.e., u j ≈ u(x j ) as an extra model variable and use a 2K-thorder polynomial to represent the solution polynomial. Specifically, the (2K + 1) model variables used in RDO-(2K + 1) are {u j±1/2 , (u x ) j±1/2 , . . . , (u x K−1 ) j±1/2 } ∪ u j .

Updating Model Variables
This section only presents the detailed formulation of RDO-2K (K ≥ 2). The procedure for RDO-(2K + 1) is almost the same and hence omitted here for simplicity. The procedure to compute time derivatives of model variables can be divided into three stages as follows.

First Stage: Reconstructing Polynomial from Model Variables
This stage finds a solution polynomial h j (x) on I j = [x j−1/2 , x j+1/2 ] to interpolate model variables defined at cell boundaries. The interpolation polynomial is of (2K − 1)-th order since there are totally 2K model variables. We first introduce a local coordinate to map I j to the standard cell [−1, 1], The usage of s(x) unifies the interpolation templates and hence reduces the computational cost in polynomial reconstruction. The solution polynomial h j (x) under the local coordinate s is thenh j (s) = h j (x(s)).
Supposeh j (s) = ∑ 2K−1 i=0 a i s i . Then the polynomial coefficients {a i } 2K−1 i=0 can be determined by the following linear system: for k = 0, . . . , K. This is a linear system about 2K unknowns {a 0 , . . . , a 2K−1 } and we can write it into the matrix form where H is the coefficient matrix and d j denotes the rescaled model variables The coefficients of the solution polynomial is then Through the high-order Hermite interpolation (4), the solution polynomial is (K − 1)th smooth globally. Introducing the local coordinate makes the inverse of the coefficient matrix H the same for all cells so as to reduce the computational costs in solving the linear system in each cell. When the condition number of this linear system is too high as K increases, one can use the symbolic tool in Maple or Matlab to obtain a sufficiently accurate inverse of the matrix to avoid the loss of accuracy.

Second Stage: Computing Temporal Derivatives on Solution Points
This part computes the temporal derivative of the obtained solution polynomial h j (x). For this, we choose 2K solution points x j,k = x j + ξ k ∆x j , k = 1, . . . , 2K on I j with k=1 are the same for each element. Through numerical and theoretical analysis later, we find that a simple choice of uniformly distributed solution points is sufficient for good numerical performance, i.e., We can now easily update the solution polynomials via updating solution values located on selected solution points. Current solution values and approximated spatial derivatives are, respectively, where dh j ds can be obtained by According to (1), the temporal derivatives of solution values become Then we use du j,k dt to retrieval the temporal derivatives of the coefficients of solution polynomials h j (s). This can be done by solving the following linear systems Note that s(x j,k ) = ξ k . This is a 2K × 2K linear system about da i dt and the coefficient matrix is the Vandermonde matrix This coefficient matrix remains the same for all cells. Therefore, we also can store the inverse in advance to reduce computational costs.

Third Stage: Retrieving the Temporal Derivatives
Using da i dt , the temporal derivatives of model variables can be directly retrieved by the following linear transform, The subscript j in the left-hand side of (14) indicates that this temporal derivative is computed from the cell I j alone. At the cell boundary x = x j+1/2 , we obtain the temporal derivatives from the I j and I j−1 at the same time, which can take different values. This leaves us the problem about how to define the proper temporal derivatives of (u x k ) j+1/2 K k=0 properly.
A natural thought is choosing the one from the upwind direction so as to obey the physical property of the hyperbolic equation. We denote dφ j+1/2 dt computed in the cell I j by dφ j+1/2 dt j for the model variable φ. Then the upwind temporal derivative of φ is where j up denotes the index in the upwind direction which is determined by the sign of f u (u j+1/2 ). To be more specific, This scheme is termed upwind RDO scheme due to the upwind property. Another idea to determine d dt φ j−1/2 is simply computing the algebraic average of the temporal derivatives from two cells just like the classical central difference scheme. Then the central temporal derivative of some model variables φ is One may consider a sum weighted by cell length in (17) when encountering the nonuniform mesh like the spectral element method [19]. However, (17) works sufficiently well, as shown in the numerical tests on non-uniform meshes.

Time Integration
We have obtained the spatial discretization and temporal derivatives of model variables, which can be described in the semi-discrete form where Φ denotes the collection of all model variables. Suppose the value of Φ at n-th time step is Φ (n) , then Φ (n+1) at the next time step is computed by the fourth-order Runge-Kutta method [20] to ensure the numerical accuracy in time. We have where ∆t is the time step length and

Discussion
The involved polynomial interpolation procedures (4) and (14) are implemented within the cell. In other words, no cross-cell polynomial interpolation is needed in the proposed schemes, which means the scheme is compact and makes it possible to solve the problems on the non-uniform mesh more flexibly than the conventional IDO methods.
We can also construct RDO-2 and RDO-3 by assigning merely function values at cell boundaries without derivative moments. Then central and the upwind RDO-2 are, respectively, the classical central difference and the upwind difference method. Besides, the spectral element method [21] with three element base function in one dimension is also the special case of central RDO-3.

Generalization to Higher Dimensions
Now we extend the algorithm to rectangular meshes to solve the two dimensional conservative equations and briefly explain how to extend this problem into three dimension. The core issue is to define proper model variables on cell boundaries and find a benign polynomial space to interpolate them. Furthermore, the reconstructed piecewise polynomials should also be at least differentiable at the cell boundaries. Therefore we also use Hermite interpolation to represent the solution polynomials.

Spatial Discretization
Suppose the computational domain is a rectangle [a, b] × [c, d] and is divided into N 1 × N 2 non-overlapping cells. We denote the cell at the i-th row and j-th column by To explain the intuition of spatial discretization, we first consider the bicubic Hermite interpolation on the cell I i,j . In this case, 16 DOFs are needed as 4 DOFs are assigned on each vertex. These 4 DOFs are u, u x , u y and the second-order mixed derivative u xy . A bit different from the 1D situation, u xy are requited here to guarantee that the reconstructed piecewise Hermite polynomials possesses continuous normal derivatives along cell edges. A similar allocation of model variables can be seen in the IDO method for solving the two-dimensional Poission Equations [22].
For general RDO-2K in 2D, the needed model variables on each vertex are ∂ i+j u ∂x i y j 0≤i,j≤K−1 , which means that each element has totally 4K 2 model variables. However, since the adjacent cells share the same model variables on common vertices, the number of effective DOFs on each element is only K 2 . We use the following polynomial spaces for representing the solution polynomial is the following full polynomial spaces [23], Its space dimension is dim(Q 2K ) = (2K) 2 . Q 2K is the tensor product of the 1D Hermite polynomials.
For RDO-(2K + 1), the reconstructed polynomial space becomes Q 2K+1 which has (2K + 1) 2 dimensions. Hence, each element needs additionally (2K + 1) 2 − (2K) 2 = 4K + 1 model variables, which include u i,j that denotes the function value on the barycenter, and K model variables, u, u n , u n 2 , · · · , u n K−1 on four edge centers. Here u n k denotes the kth-order normal derivative on the edge center.

First Stage: Reconstructing Solution Polynomial
The first step interpolates the solution polynomials h i,j (x, y) on I i,j . To begin with, we also introduce a local coordinate to map the cell Similar to the one-dimensional case, h i,j (x, y) is also transformed intoh i,j (r, s) under the local coordinate. Correspondingly, 4K 2 model variables serves as constraints. The solution polynomialh The inverse of the coefficient matrix of the linear system (23) about a α,β can also be solved in advance so as to effectively save computational costs in solving (23).

Second Stage: Computing Temporal Derivatives on Solution Points
This stage represents the solution polynomial on solution points and update it. The chosen solution points in 2D are tensor products of 1D points. Specifically, these 4K 2 points in the local coordinate are in the form of (s k , r l ) = (ξ k , ξ l ), k, l ∈ {1, . . . , 2K}.
Then the temporal derivatives of solution values become Then the temporal derivatives of coefficients of the solution polynomial, da α,β dt can be obtained by solving the following linear system,

Third Stage: Retrieving the Model Variables
From da i,j dt , we can obtain the temporal derivatives of model variables by the following relationships for RDO-2K for l, m ∈ {0, · · · , 2K − 1}. The subscript (i, j) in RHS of (30) specifies that the temporal derivatives are computed from I i,j and do not involve any external information. To properly introduce the external interaction, both central and upwind schemes can be reconstructed similar to the one dimensional case.

Three Dimensional Case
In three dimensions, RDO-2K on the cuboid mesh can also be constructed similarly via the polynomial space and chooses ∂ i+j+k u ∂x i y j z k 0≤i,j,z≤K−1 on each vertex as model variables. There are totally (2K) 3 model variables in each cell. The remaining procedure is exactly the same as in lower dimensional space and we do not repeat the formulation here.

Fourier Analysis
This part analyses the accuracy and the stability of presenting scheme in solving the linear advection equation. For concreteness, we only give a detailed analysis for RDO-4 since the generalized RDO-M can be analyzed following the same routine.

Formulation of the Amplification Matrix
Consider the linear scalar advection equation Suppose the spatial domain is [0, L] and is discretized by a uniform mesh spacing ∆x. The initial condition is periodic, u ini = e iwx/∆x where i denotes √ −1 and the scaled . Then the initial model variables for the upwind RDO-4 on I j = x j−1/2 , x j+1/2 are The core in Fourier analysis is reformulating the algorithm into the matrix form using the periodicity of u ini and S denotes the amplification matrix. We do this step by step as follows. Stage 1 described in Section 2 reconstructs the Hermite polynomial coefficients a j from the rescaled model variables d j in the cell I j , where H is defined by (5) and d j is the composition of the model variables after the coordinate transform, In Stage 2, the solution values and derivatives on solution points are, respectively, and where V and D can be found in (13) and (10). According to (31), the temporal derivatives on the solution points are d dt Then the temporal derivatives on solution points are transformed back to the temporal derivatives of model variables, which is the inverse problem of (37): where the subscript j in the left-hand side indicates that the temporal derivatives are computed locally form j-th cell, andS = −2HV −1 DH −1 .
In the linear case, V −1 D is an invariant about the choice of solution points. Furthermore, the matrix V −1 D has a simple structure Therefore, it can be easily concluded that RDO schemes using different sets of solution points are identical in solving the linear advection equation as the following theorem tells.
Moreover, it is noticed that the in-cell temporal derivatives (39) Then for the upwind discretization, the periodicity of the initial profile (42) tells that Hence, we can obtain the amplification matrix S up for the upwind scheme as follows, with Similarly, the temporal derivatives of d j for the central RDO-4 are Combining (39) and (47) we have

Accuracy Analysis
We can see that both the central and the upwind RDO-4 have fourth spatial local discretization errors. However, this does not suffice to ensure the third order accuracy in numerical experiments. In fact, the order accuracy is determined by the principal eigenvalues σ(w) according to [18], with σ(w) being the principle eigenvalue of the amplification matrix S. Other eigenvalues of S are spurious ones. Specifically, if σ(w) can be expanded as a Taylor series in w, then the scheme is m-th order accurate. Since it is difficult to calculate an expression manually in practice, we use the symbolic computation tools in Matlab .
For the upwind RDO-4, we can obtain the principal eigenvalue of the amplification matrix S upwind is This indicates that the upwind RDO-4 is third order accurate in space. Besides, the term − 1 72 w 4 introduces a small amount of dissipation which tends to smear the solution, which is certificated in the numerical experiments.
The principal eigenvalue of the amplification matrix S central for central RDO-4 is and we can see that only second order spatial accuracy is maintained for the central scheme.
It is also worthy noting that − i 16 w 3 introduces slight dispersion, which means that phase speed varies as w changes. Therefore, during the long-term simulation of a complicated initial profile that contains waves of different frequency, the profile will be distorted finally. This is also verified in numerical experiments in Section 5.
The accuracy analysis for higher-order RDO schemes is exactly the same as RDO-4. However, a concise Taylor expansion of principal eigenvalue like (51) and (52) is very difficult (if not impossible) to obtain even using the mathematical software. Here we try to calculate the truncation order of σ(w) of a scheme by the following approximation. First, choose a proper w, say w = π/4. Then the error is Next, we halve the wave number w and obtain the error corresponding to w/2 Noticing that for the m-th order accurate scheme, one has δ(w/2) ≈ 1 2 m+1 δ(w). (56) Hence the order of accuracy can be approximated by a simple logarithmic operation The modulus operation is required here since the error result can be complex valued. Errors and estimated accuracy orders for the upwind and the central schemes can be seen in Tables 1 and 2. The errors for the upwind schemes agree with the local discretization error well. The upwind scheme using M model variables per cell achieves (M − 1)-th accuracy order strictly. However, the accuracy orders for the central schemes varies quite irregularly. Specifically, when M = 5 and 9, the accuracy orders for central schemes remain consistent with that of upwind schemes. When M = 6 and 7, central schemes achieve higher accuracy orders than upwind schemes. Nevertheless, there also some choices of M like M = 4 and 8 such that the accuracy order of the central scheme is lower than the corresponding upwind scheme. It is remarkable that when M = 7, the central scheme enjoys two orders higher accuracy than that of the local discretization error.   Tables 1 and 2 and also contain abundant information about the dissipation and dispersion of the schemes. It should be noted that a positive real part of δ(w) will cause the scheme eventually blow up if no additional dissipation is added in time stepping. In Table 2, the real parts of δ(w) are negative for all upwind schemes which introduces a certain amount of spurious dissipation so as to stabilize the scheme. From Table 1 we can see that δ(w) are pure imaginary numbers with all choices of M, which means that the dispersion is dominant in the numerical error for the central schemes. However, it can be easily observed that the dispersion and dissipation error decays rapidly as M increases for the proposed schemes.

RDO-M Coarse Mesh Error (w = π/4) Fine Mesh Error (w = π/8)
In (39), one can see that the numerical error in computing H −1 dominates the error of temporal derivatives, since H and V −1 D have explicit expression. However, the error of computing H −1 is inevitably increasing as M increases, which indicates that there is a limit of accuracy for the presented method. Denote the numerical inverse of H byĤ, which is obtained by Matlab. It is well known that the error H −1 −Ĥ hinges on the condition number of H. Moreover, a more accurate estimate H −1 −Ĥ can be bounded by the following inequality, Then E inv determines the smallest error that RDO-M can reach through the grid refinement. The results of cond(H) and E inv with M varying form 4 to 11 are shown in Table 3. We can see that cond(H) and E inv grow quickly as M increases. Hence, there is a tradeoff between the higher convergence rate and E inv . In this work, we suggest that M ≤ 9 is ideal for the practical numerical experiments.

Stability Analysis
Now we investigate the maximum tolerant time step length when using the Runge-Kutta time integration scheme. Let the time step length be ∆t = c∆x and c is the Courant number for (31). The amplification matrix for the proposed schemes using the fourth order Runge Kutta method for time integration (19) can be simplified as To ensure the stability, the Courant number c should be chosen such that the spectral radius of the amplification matrix R, i.e., the largest absolute value of its eigenvalues, is not greater than 1 for all wavenumbers w. The spectral radius of R for the upwind schemes and the central schemes are presented in Figures 2 and 3, respectively, which can also tell the Courant number limitation related to the stability condition.
Both types of schemes have a fairly large permissible range of c. As more moments are used, the stability condition becomes more strict for both types of RDO schemes. Generally, the Courant number limitation of the central scheme is larger than that of the upwind schemes. For example, the Courant number limit for upwind RDO-5 is about 0.25 and this number is approximately 0.39 for the central RDO-5, although both schemes are fourth-order accurate. The Courant number limits for RDO schemes are presented in Table 4.

Numerical Results
This section provides several numerical tests to verify the formal accuracy and other numerical properties like diffusion and dispersion of the proposed schemes.

One-Dimensional Linear Case
We first test the accuracy when solving the scalar hyperbolic Equation (31) analysed in Section 4. The initial condition is u(x, 0) = sin(x) on [0, 2π) with the periodical boundary condition. The exact solution is u(x, t) = sin(x − t).
Both the central and the upwind RDO-M with the number of moments M varying from 4 to 9 are tested on uniform meshes. To measure the accuracy, two typical norms for errors are used here, We study the convergence rate by recording the errors at t = 2π and gradually refining the grids. The time step length is set as a very small number such that the temporal error does not influence the accuracy results. The error results are presented in Table 5. In general, the orders of accuracy computed from the grid refinement agree with the Fourier analysis results in Tables 1 and 2  Besides the uniform meshes, we also consider the non-uniform meshes that are randomly generated by setting cell boundaries as where N denotes the number of cells, ∆x = L/N, ζ j N j=1 are independently random variables, η represents the magnitude of perturbation and Unif(−1, 1) denotes the uniform distribution over [−1, 1]. We set η = 0.1 in the present work. In this setting, the maximal and the minimal possible values of ∆x j are, respectively, 0.8∆x and 1.2∆x.
The accuracy results of RDO schemes over non-uniform meshes are listed in Table 6. It can be observed that the lack of uniformity in the mesh causes slightly larger errors and does not influence the convergence rate essentially. Therefore, the proposed schemes works consistently well on both uniform and non-conform meshes.

Two-Dimensional Linear Case
Consider the following two-dimensional scalar hyperbolic equation The initial condition is u(x, y; 0) = sin(x + y) on [0, 2π] × [0, 2π] with the periodical boundary condition. The exact solution is u(x, y; t) = sin(x + y − 2t). The errors and the orders of accuracy for the central and the upwind RDO-M with M changing from 4 to 8 are shown in Table 7.
The errors and accuracy orders of the two-dimensional schemes are similar to that of the one-dimensional schemes, which illustrates that both the upwind and the central RDO schemes work well in two dimensions. Specifically, central RDO schemes achieves significantly lower errors compared with upwind RDO schemes when M = 6, 7. However, for any other presented values of M, upwind RDO performs better.

Nonlinear Case
This case considers solving the one-dimensional Burger's equation The initial condition is u(x, 0) = 0.5 + sin x with the periodical boundary condition for x ∈ [0, 2π]. We compute the numerical solution at t = 0.2π when the profile is still smooth and record the errors under grid refinement. The results can be seen in Table 8. Both schemes lose accuracy due to nonlinearity compared with results of the linear cases shown in Table 5. The accuracy orders of central RDO schemes coincide the accuracy orders of local polynomial reconstruction better when M ≥ 5. However, upwind schemes generally enjoy smaller errors in this test.

Advection of the Gaussian Profile
To investigate the diffusion and the dispersion of the proposed schemes and verify the long-term stability, a Gaussian profile used in [18] u( is set as the initial condition for the advection Equation (31) Tables 1 and 2 well. Table 1 indicates that the error of central RDO-M schemes are dominated by dispersion for all M. By contrast, Table 2 reveals that the upwind RDO-M suffers from diffusion mainly for all M with the only exception being M = 5. Correspondingly, the numerical solution by central RDO-4 and RDO-5 are remarkably distorted by the dispersion, especially in the long-term simulation. On the other hand, the upwind RDO schemes maintains a better and meanwhile a more flat profile than the central schemes do for all choices of M due to diffusion.

Comparison of Computational Costs
From the stability analysis, we can know that the higher-order schemes enjoy better accuracy at the expense of decreasing the CFL number and hence the efficiency. Therefore, it is necessary to compare the practical computational efficiency for RDO-M to achieve the same accuracy.
Then we consider simulating the advection of the Gaussian profile in the last test. All schemes use the fourth-order Runge-Kutta method with the temporal step length ∆t = 0.95c max ∆x for time integration, where c max denotes the largest tolerable Courant number. Furthermore, we take ∆t = 0.7c max for the non-uniform mesh since the stability condition is more stringent in the presence of mesh non-uniformness. The target accuracy is set as E 1 = 10 −5 . The number of cells N is increased gradually until the norm-1 error reaches the target accuracy. The minimal required N and the corresponding CPU time for RDO schemes are shown in Table 9. It should be remarked that the total DOFs are not MN since the model variables at each cell boundary are shared by two cells.
For both uniform and non-uniform meshes, the significant reduction can be easily observed in the number of cells, the computational storage (i.e., total DOFs) and CPU time as M and the accuracy orders increases. RDO schemes on non-uniform meshes generally require slightly more CPU time to reach the target accuracy as shown in Table 9. Hence, high-order schemes are more attractive in terms of computational efficiency.

Conclusions
A family of compact multi-moment schemes that use high-order derivative moments have been proposed, analyzed and tested in one dimension and two dimensions. Using the model variables at cell boundaries, the proposed schemes reconstruct the solution polynomial within a single cell and then retrieve the temporal derivatives of model variables efficiently. The central scheme also achieves excellent performance in several tests, although this scheme does not possess clear physical meaning as the upwind scheme does. This indicates that one may extend the scheme to the more complicated scenarios without considering finding the upwind direction.
The future work will focus on extending this formulation to the more generalized unstructured triangular mesh using the Argyris polynomials [24] which also have continuous gradient globally. Furthermore, designing the conservative counterpart of the current version is also on our schedule.
Author Contributions: S.L., conceptualization, methodology, software, software, validation, writingoriginal draft preparation; Q.L., formal analysis, writing-review and editing; H.L., supervision; J.S., project administration, resources. All authors have read and agreed to the published version of the manuscript.