Weighted essentially nonoscillatory scheme on unstructured quadrilateral and triangular meshes for hyperbolic conservation laws
Abstract
In this paper, a thirdorder weighted essentially nonoscillatory (WENO) scheme is developed for hyperbolic conservation laws on unstructured quadrilateral and triangular meshes. As a starting point, a general stencil is selected for the cell with any local topology, and a unified linear scheme can be constructed. However, in the traditional WENO scheme on unstructured meshes, the very large and negative weights may appear for the mesh with lower quality, and the very large weights make the WENO scheme unstable even for the smooth tests. In the current scheme, an optimization approach is given to deal with the very large linear weights, and the splitting technique is considered to deal with the negative weights obtained by the optimization approach. The nonlinear weight with a new smooth indicator is proposed as well, in which the local mesh quality and discontinuities of solutions are taken into account simultaneously. Numerical tests are presented to validate the current scheme. The expected convergence rate of accuracy is obtained, and the absolute value of error is not affected by mesh quality. The numerical tests with strong discontinuities validate the robustness of current WENO scheme.
keywords:
Thirdorder WENO scheme, Unstructured meshes, Hyperbolic conservation laws, Finite volume method.1 Introduction
In recent decades, there have been continuous interests and efforts on the development of highorder schemes. For engineering applications, the construction of highorder numerical schemes on unstructured meshes becomes extremely demanding due to the complex computational domains. There are a gigantic amount of publications about the introduction and survey of highorder schemes, including discontinuous Galerkin (DG), spectral volume (SV), spectral difference (SD), correction procedure using reconstruction (CPR), essential nonoscillatory (ENO), weighted essential nonoscillatory (WENO), Hermite WENO (HWENO) and exact method, etc.
The DG scheme was first proposed DG1 () to solve the neutron transport equation. The major development of DG method was carried out for the compressible Euler equations DG2 (); DG3 (); DG5 (). Highorder accuracy is achieved by means of highorder polynomial approximation within each element rather than by means of wide stencils. Because only neighboring elements interaction is included, it becomes efficient in the application of complex geometry. Recently, a correction procedure via reconstruction framework (CPR) was developed CPR1 (); CPR2 (); CPR3 (); CPR4 (). The CPR formulation is based on a nodal differential form, with an elementwise continuous polynomial solution space. By choosing certain correction functions, the CPR framework can unify several well known methods, such as the DG, SV SV (), and SD SD () methods, and lead to simplified versions of these methods, at least for linear equations. The ENO scheme was proposed in Harten (); ENOShu () and successfully applied to solve the hyperbolic conservation laws and other convection dominated problems. The “smoothest” stencil is selected among several candidates in reconstruction to achieve high order accuracy and keep essentially nonoscillatory near discontinuities. Based on unstructured triangular meshes, the ENO scheme was developed as well Hartenunstructured (); Abgrall (); Sonar (). Following the ENO scheme, the WENO scheme WENOLiu (); WENOJS (); WENOM (); WENOZ (); Zhao () was developed. With the nonlinear convex combination of candidate polynomials, WENO scheme achieves higher order of accuracy and keeps essentially nonoscillatory property. WENO scheme improves upon ENO scheme in robustness, smoothness of fluxes, steadystate convergence, provable convergence properties, and more efficiency. The classical development of WENO scheme includes WENOJS scheme WENOJS (), WENOM scheme WENOM (), WENOZ scheme WENOZ (), et al. Based on the Hermite interpolation, HWENO schemes were developed HWENO1 (); HWENO2 (); HWENO3 () for the hyperbolic conservation laws. On the unstructured meshes, the WENO scheme with the same order of accuracy as the fundamental ENO scheme was constructed by Friedrichs Friedrich (). Compared with the fundamental ENO scheme, the WENO scheme is very smooth and stable in smooth regions, while the optimal order of accuracy is not obtained WENOJS (). Hu and Shu HuShu () presented the thirdorder and fourthorder WENO schemes on unstructured triangular meshes. The highorder of accuracy is obtained by the combination of lower order polynomials, which is similar with the onedimensional WENO scheme. However, its successful application was limited by the appearance of the negative linear weights, and the problem appears commonly on the unstructured meshes. For a mesh that is close to the regular meshes, such WENO scheme works well by a regrouping approach to avoid negative weights. However, for a mesh with lower quality, the large linear weights appears and WENO schemes become unstable even for the smooth flows. A splitting technique dealing with the negative linear weights was proposed in splittingweights (). This technique makes the WENO schemes HuShu () generally work well on better quality meshes. However, the linear weights can be very large with the lower quality mesh, and the WENO schemes become unstable even if the splitting technique is adopted. In order to avoid the negative linear weights and very large linear weights, a class of WENO schemes were proposed in Dumbser (). Due to the very large stencils, the boundary treatment becomes complex and consuming computer memory, and the computational efficiency is reduced because of the inefficient combination in reconstruction. Based on the two classes of WENO schemes in HuShu (); Dumbser (), a hybrid scheme is presented in hybridWENO () to deal with the problem with very large linear weights, while the scheme contains artificial parameter and the reconstruction is more complex. Another type of higherorder finite volume scheme is based on the socalled exact reconstruction kexact1 (); kexact2 (); kexact3 (). To guarantee the nonsingularity of the reconstruction, the exact reconstruction relies on the leastsquare approach.
In this paper, a thirdorder WENO scheme is developed on unstructured quadrilateral and triangular meshes for the hyperbolic conservation laws. As a starting point of WENO reconstruction, a general stencil is selected for a cell with any local topology. With the selected stencil, a unified linear scheme was constructed. With the linear candidate polynomials from the candidate substencils, the linear wights can be obtained. However, the appearance of very large linear wights on unstructured meshes causes the WENO schemes to be unstable HuShu (); hybridWENO (). For both quadrilateral and triangular meshes, an optimization approach is introduced to deal with the very large linear weights. The splitting technique is considered to deal with the negative weights obtained by the optimization approach. The nonlinear weight with a new smooth indicator is defined for the solutions with discontinuities. With the optimization approach for very large weights and the splitting technique for negative weights, the scheme becomes more robust compared with the previous WENO scheme HuShu () with the low quality meshes. A variety of numerical tests from the accuracy test to the solutions with strong discontinuities are presented to validate the accuracy and robustness of the current scheme.
This paper is organized as follows. In Section 2, the finite volume scheme on unstructured meshes are introduced. The thirdorder WENO schemes on unstructured quadrilateral and triangular meshes are presented in section 3. Section 4 includes numerical tests to validate the current algorithm. The section 5 is the conclusion.
2 Finite volume method
In this paper, the twodimensional Euler equations are considered
(1) 
where , , , and is the specific heat ratio. For a polygon cell , the boundary can be expressed as
where is the number of cell interface for cell . Integrating Eq.(1) over the cell , the semidiscretized form of finite volume scheme can be written as
(2) 
where is the cell averaged value over cell , is the area of , and is the outer normal direction of .
In this paper, a thirdorder spatial reconstruction will be introduced, and the line integral over is discretized according to Gaussian quadrature as follows
(3) 
where is the length of the cell interface , are the Gaussian quadrature weights, and the Gaussian quadrature points for are defined as
where are the endpoints of .
For the twodimensional Euler equations, the rotational invariance property RiemannSolver () can be expressed as
(4) 
where is the rotation matrix
The numerical fluxes in Eq.(3) at the Gaussian quadrature points in the global coordinate can be given with the fluxes in the local coordinate according to Eq.(4), and the flux in the local coordinate can be approximated by Riemann solvers
where and are the reconstructed variables at both sides of Gaussian quadrature points of cell interface. In this paper, the HLLC Riemann solver RiemannSolver () is adopted for numerical tests.
3 Thirdorder WENO scheme on unstructured meshes
The reconstruction procedure is the main theme of this paper, and it will be presented on both unstructured quadrilateral and triangular meshes in this section.
3.1 Linear reconstruction on quadrilateral meshes
As a starting point of WENO reconstruction, a linear reconstruction will be presented. For a piecewise smooth function over cell , a polynomial with degree can be constructed to approximate as follows
In order to achieve the thirdorder accuracy and satisfy conservative property, the following quadratic polynomial on the cell is constructed
(5) 
where is the cell average value of over cell and are basis functions, which are given as follows
(6) 
In order to construct the polynomial uniquely, the stencil for reconstruction is given in Figure.1. For arbitrary unstructured quadrilateral meshes, the cells must exist, and the cells may not exist since it is common that only three or four cells share a common node. Thus, preconditioning is needed for a general quadrilateral mesh, and the procedure of selecting stencil for cell is given as follows

initialize the parameter , which is corresponding to the cell ;

read the neighboring cell of cell from the mesh file;

read the neighboring cell of cell , if each of them is not exist, record and is replaced by ,
read the neighboring cell of cell , if each of them is not exist, record and is replaced by ; 
read the neighboring cell of cell , if each of them is not exist, record and is replaced by ,
read the neighboring cell of cell , if each of them is not exist, record and is replaced by .
To determine the polynomial , the following constrains need to be satisfied for the cells
(7) 
where is the cell averaged value over . According to the definition of , Eq.(7) can be written into the following linear system
where , , and is given as
Traditionally, the following least square method is used to obtain the coefficient
(8) 
However, the matrix may not be invertible. In order to obtain the coefficients in the quadratic polynomial uniquely, a more robust least square method is adopted
(9) 
where I is the unit matrix, is a small number, and takes in the computation. The matrix is invertible, and the quadratic polynomial can be determined uniquely. Meanwhile, if the matrix is invertible, the difference between the solution of Eq.(9) and Eq.(8) is . With the coefficient , the point value of at the Gaussian quadrature point can be written as a linear combination of the cell averaged values
(10) 
where if .
Remark 1
Denote , and the inverse of can be expressed as
(11) 
Multiplying on both side of Eq.(11), we have
According to the HamiltonCaylay theorem, can be expressed as
(12) 
where are given by
Substituting Eq.(12) into Eq.(11), the linear system for can be written as
If the matrix is invertible, and the solution of the linear system above can be given as
where corresponds to the solution of the system with . Substituting c into Eq.(11), and supposing the solution of Eq.(8) is , the solution of Eq.(9) can be given as
Thus, the difference between solutions of Eq.(9) and Eq.(8) is the order of .
Similar with the standard WENO reconstruction WENOLiu (); HuShu (), twelve substencils are selected from the large stencil given in Figure.1
Twelve candidate linear polynomials corresponding to the candidate substencils can be constructed as follows
where is the th cell in the substencils , and is the cell averaged value over cell . With the robust least square method introduced above, twelve linear polynomials can be fully determined. Similarly, the point value of at the Gaussian quadrature point can be written as a linear combination of the cell averaged values . However, in order to calculate the linear weights, it must be taken into account that the cell may not exist for arbitrary unstructured quadrilateral meshes. According to the local topology, the candidate polynomials is rewritten as follows

For , the linear polynomial can be expressed as

For , the linear polynomial can be expressed as
For the linear scheme, the linear combination of
satisfies
(13) 
where is the Gaussian quadrature point, and is the linear weights only depending on the local geometry of mesh. Substituting into at , we have
(14) 
Comparing the coefficient of the cell average value in Eq.(10) and Eq.(14), the following thirteen linear equations can be obtained
and its matrix form is
(15) 
where , , and . Similar with the analysis in HuShu (), the linear system for Eq.(15) contains thirteen linear equations, and generally the system is underdetermined with . and reproduce the linear function exactly and Eq.(13) is valid for , and under the following identical constraint
Two equations can be eliminated from the linear system with thirteen equations. The traditional least square becomes underdetermined, and the following least square method are used to obtain the linear weights
3.2 Linear reconstruction on triangular meshes
For the reconstruction on unstructured triangular meshes, the procedure is similar with the reconstruction on quadrilateral meshes, except the choice for stencil and substencils. The stencil is given in Figure.2, and the quadratic polynomial Eq.(5) can be constructed on cell as well. For the cell , the neighboring cell must exist and their neighboring cells may not exist, and the similar preconditioning is needed as well, and the procedure is given as follows

initialize the parameter , which is corresponding to the cell ;

read the neighboring cell of cell from the mesh file;

read the neighboring cell of cell , if each of them is not exist, record and is replaced by ;

read the neighboring cell of cell , if each of them is not exist, record and is replaced by .
In order to obtain the quadratic polynomial, the following constrains need to be satisfied
and the matrix form can be written as
where , , and is given as
Similar with the quadrilateral meshes, the polynomial can be determined. For the triangular meshes, nine candidate substencils are selected
and nine linear polynomials can be fully determined. Similarly, the point value of at the Gaussian quadrature point can be written as a linear combination of the cell averaged values . However, the cell may not exist for arbitrary unstructured meshes. In order to calculate linear weights later, the candidate polynomials is given according to the local topology as follows

For , the linear polynomial can be expressed as

For , the linear polynomial can be expressed as
For the linear scheme, the linear combination
satisfies
(16) 
where is the Gaussian quadrature point, and is the linear weights. The procedure of determining the linear weights for quadrilateral meshes can be extended to triangular meshes directly.
3.3 Optimization approach for linear weights
With the procedure introduced above, the linear weights for WENO reconstruction can be obtained. To deal with the flow with discontinuities, the nonlinear weights need to be constructed, which is one of the most important step in the WENO schemes. However, for the unstructured meshes, the very large linear weights and negative linear weights appear commonly, which makes the WENO schemes unstable even for the smooth flow problems. In order to construct a robust WENO schemes on unstructured meshes, the very large linear weights and negative linear weights need to be treated carefully. To illustrate the very large weights, the quadrilateral and triangular meshes in Figure.4 and Figure.4 are considered, where the meshes with cell size are given. The meshes in the left column are closer to the uniform quadrilateral and triangular meshes, and they denoted as the regular meshes. Meanwhile, the meshes in the right column are denoted as irregular meshes. The maximum of the liner weights for these meshes are given in Table.2, where the maximum becomes larger with the mesh refinement. Such very large weights will trigger the instability of WENO scheme.
The meshes with lower quality will result in the very large linear weights, and an optimization approach is introduced to deal with them. In the current optimization approach, the weighting parameters is introduced for each cell, such that the ill cell contributes little to the quadratic polynomial in Eq.(10). For quadrilateral meshes, the weighting parameters are defined
For triangular meshes, the weighting parameters are defined
The optimized coefficients of the quadratic polynomial in Eq.(10) are given by the following weighted linear system
where for the quadrilateral meshes, and for the triangular meshes. The maximum of optimized liner weights for the meshes in Figure.4 and Figure.4 are given in Table.2. With the mesh refinement, the maximum of liner weights is reduced greatly. This linear system is solved with the least square method proposed above. With the procedure above, the maximum becomes the order .
3.4 Nonlinear weights
Before implementing the WENO reconstruction, the splitting technique splittingweights () is considered for the negative linear weights which is obtained by the optimized approach as follows
where is taken in numerical tests, and the scaled nonnegative weights are given by
With the nonnegative weights, the splitting polynomial is written as
which satisfies
Based on the nonnegative weights obtained by splitting approach, the nonlinear weights are obtained as
where is a small positive number, is a new smooth indicator defined on unstructured meshes
and is defined as
where is a multiindex and is the derivative operator. In the new smooth indicator, the accuracy keeps the original order in smooth regions with . In addition, the linear weights are considered in constructing the new smooth indicator. Suppose that is a very large linear weight and the corresponding contains the discontinuity. The ratio of nonlinear weights of and becomes
The contribution of the candidate polynomial becomes very little, and the essentially nonoscillatory property retains for the very large linear weight. The final reconstructed value at the Gaussian quadrature points can be written as
4 Numerical tests
In this section, the numerical tests will be presented to validate the current WENO scheme. For the temporal discretization, the classical thirdorder TVD RungeKutta method TVDRK () is used for Eq.(2), which can be written as follows
where is the operator corresponding to the spatial discretization.
In order to eliminate the spurious oscillation and improve the stability, the reconstruction for Euler system can be performed on the characteristic variables. The Jacobian matrix at one quadrature point is , and the Roe’s mean matrix is used, which can be diagnoalized by the right eigenmatrix . The variables for reconstruction are defined as . With the reconstructed values, the conservative variables can be obtained by the inverse projection.
4.1 Accuracy tests
The first case is the advection of density perturbation, and the initial condition is given as follows
The computational domain is , and the periodic boundary conditions are applied in both directions. To keep the thirdorder accuracy, is used for the thirdorder scheme where is the scale of meshes. In order to validate the accuracy of the current scheme, both quadrilateral and triangular meshes are tested. The meshes with mesh size in Figure.4 and Figure.4 are shown as example, and the density distribution with cell size are given in Figure.6 and Figure.6. The and errors and convergence orders are presented in Table.4 and Table.4 for quadrilateral and triangular meshes at . The expected order of accuracy are obtained, and the order of error is not affected by the quality of meshes.
The second accuracy test is the isentropic vortex propagation problem. The mean flow is . An isentropic vortex is added to the mean flow with the following perturbations
where , , and the vortex strength . The computational domain is , and periodic boundary condition is applied to all boundaries. The regular and irregular quadrilateral meshes, regular and irregular triangular meshes are tested respectively. The and errors and convergence orders for current scheme are presented in Table.6 and Table.6 for quadrilateral and triangular meshes at . The expected order of accuracy are obtained, and the order of error is not affected by the quality of meshes as well. With the fixed mesh size , the longtime evolution of the vortex is computed. The density distribution at center line with and time periods on quadrilateral and triangular irregular meshes are given in Figure.9 with mesh size . The dissipation is well controlled by the current scheme even after long time evolution. Due to more number of cell with the same mesh scale, the scheme performs better on triangular meshes than quadrilateral meshes.
4.2 One dimensional Riemann problem
In this case, four onedimensional Riemann problems are tested by the thirdorder WENO scheme on the quadrilateral and triangular meshes. The first one is the Sod problem, and the initial condition is given as follows
The second one is the Lax problem, and the initial condition is given as follows
For these two cases, the computati