An Improved WENO-Z Scheme with Symmetry-Preserving Mapping

Since the classical weighted essentially non-oscillatory (WENO) scheme is proposed, various improved versions have been developed, and typical one is the WENO-Z scheme. Although better resolution is achieved, it is shown in this article that, the result of WENO-Z scheme suffers evident distortion in the long-time simulation of the linear advection equation. In order to fix the problem of WENO-Z, a symmetry-preserving mapping method is proposed in this article. In the original mapping method, the weight of each sub-stencil is used to map, which is demonstrated to cause asymmetric improvement about a discontinuity. This asymmetric improvement will lead to a distorted solution, more severe with longer output time. In the symmetry-preserving mapping method, a new variable related to the smoothness indicator is selected to map, which has the same ideal value for each sub-stencil. Using the new mapping method can not only fix the distortion problem of WENO-Z, but also improve the numerical resolution. Several benchmark problems are conducted to show the improved performance of the resultant scheme.


Introduction
Weighted essentially non-oscillatory (WENO) finite difference/volume scheme has become one of the most popular shock-capturing schemes to solve hyperbolic conservation laws, which are developed from the essentially non-oscillatory (ENO) schemes [1][2][3]. Unlike linear schemes using a single fixed stencil to interpolate, an r th-order ENO scheme chooses r stencils as candidates. The smoothest stencil is selected from r candidates with the aid of smoothness indicator. It is reasonable to abandon discontinuous stencils when there is a discontinuity, however, in the smooth region, r stencils together can reach ( ) 21 r − th-order accurate interpolation. Therefore, the ENO schemes do not make full use of all stencils' information in the smooth region, leading to the problem of accuracy losing and computation wasting.
In order to utilize the information of less-smooth stencils in ENO schemes, Liu et al. [4] introduced the WENO schemes. Using the smoothness indicator proposed by Liu et al., WENO scheme adopts a convex combination of r stencils, and an r th-order ENO scheme can be converted to an ( ) 1 r + th-order WENO scheme. Jiang and Shu [5] redefined the smoothness indicator of stencils, so that an r th-order ENO scheme can be converted to an ( ) 21 r − th-order WENO scheme, which is known as the WENO-JS scheme. However, Henrick et al. [6] pointed out that the actual accuracy order of WENO-JS scheme is less than the optimal one in smooth region with critical points. To recover the accuracy at critical points, they introduced an extra mapping process where the weights of stencils are mapped to more accurate values, and the resultant scheme is known as WENO-M. Then, Feng et al. [7,8] further improved the WENO-M scheme by optimizing the mapping function involved in the mapping process. To increase the relevance of less-smooth stencils, Borges et al. [9] took both local and global smoothness indicators into consideration, and the resolution of scheme in both discontinuous and smooth regions is significantly improved. Castro et al. [10] then extended this idea to higher order schemes, which is known as the WENO-Z. From another perspective, Fedkiw et al. [11,12] observed that a small parameter e used in formula to avoid the denominator being 0, has evident influence on the accuracy of WENO scheme. Shen et al. [13] investigated the effect of e on the convergence property of WENO scheme in detail, and their results showed that 2 10  − = can significantly improve the numerical solution for subsonic and transonic flows. Don et al. [14,15] theoretically and numerically demonstrated that if  is taken as a power function of the uniform grid space, i.e.
( ) m x  =   , the WENO-JS and WENO-Z schemes can retain the optimal order in the smooth region, regardless of critical points. However, a larger  leads to a less stable scheme, and thus the effect of adjusting  to recover accuracy is problem-dependent. Acker et al. [16] pointed out that to improve the resolution on a coarse grid, increasing the relevance of less-smooth stencils is much more important than recovering accuracy at a critical point.
By adding a new term that increases the relevance of less-smooth stencils to the original formula of weight, they proposed an improved version of the fifth-order WENO-Z scheme. Besides WENO-M and WENO-Z, other variants of WENO scheme have emerged and gained attention. Using a detector as weighting function to couple WENO reconstruction and linear compact scheme, hybrid method was proposed in [17,18] to obtain low-dissipation and low-dispersion scheme for turbulence and aeroacoustics. Then, Sun et al. [19,20] developed a linear minimum dispersion and controllable dissipation scheme to hybridize with WENO scheme to retain good performance in flow with shocks. However, hybrid function is critical in hybrid methods and still a big challenge.
To reduce the dissipation, Martin et al. [21] proposed WENO-SYMBO scheme through adding a downwind stencil and optimizing the ideal coefficients for better spectral property. Hu et al. [22] proposed a new weighting strategy and new smoothness indicator for the downwind stencil, leading to the adaptive central-upwind sixth-order WENO scheme (WENO-CU6). However, since the background central scheme cannot suppress disturbances from accumulated dispersion errors, spurious waves may appear at marginally-resolved wave numbers. Based on Lagrangian interpolation polynomial, Pan et al. [23,24] proposed a new smoothness indicator with more succinct form compared to the classical one by Jiang and Shu [5], leading to a scheme named WENO- . Don et al. [25] found that due to the numerical unstable form of the local smoothness indicators, asymmetry occurs in an otherwise problem for seventh and higher-order WENO schemes. Then, they introduced an equivalent but compact and numerical stable form of the smoothness indicators. Recently, Fu et al. [26][27][28] proposed a family of targeted ENO (TENO) schemes. By using a set of low-order candidate stencils with incrementally increasing width, discontinuities and small-scale fluctuations are efficiently separated, and the numerical dissipation is significantly diminished by an ENO-like stencil selection.
Due to better resolution while with almost the same computational cost as WENO-JS, WENO-Z has gained extensive attention. In this article, from the simulation of linear advection equation, it is observed that the numerical solution of WENO-Z suffers distortion problem near a discontinuity. To overcome this deficiency, a symmetry-preserving mapping method is proposed to improve the WENO-Z scheme. This article is organized as follows. Section 2 describes the classical WENO-JS scheme and WENO-Z scheme. In Section 3, the original mapping method and the symmetry-preserving mapping method are introduced in detail. The mapped WENO-Z is proposed in Section 4. Numerical experiments with one-and two-dimensional benchmark problems are presented in Section 5. Conclusions are drawn in Section 6.

WENO Reconstruction
Here, a grid with uniform space x D is used, and the positions of cell center and cell boundary are represented ,, S S S , and each sub-stencil contains three adjacent grid cells. In the smooth region, WENO5 combines each sub-stencil with corresponding ideal weights to obtain a fifth-order approximation of ( ) hx. While in the region with discontinuities, to avoid numerical oscillation, sub-stencils containing the discontinuity are suppressed by assigning small weights to them. Letĥ denote the approximation of ( ) hx, and the three sub-stencils of WENO5 give where the superscript is used to distinguish different sub-stencils.
In the smooth region, the weights approach their ideal values, which are , . 10 10 10
In the context of this article, a critical point refers to the point where the first-order derivative is zero while the second-order derivative is non-zero. From Eqs. (11) and (12), two observations can be obtained in the smooth region as follows.
 At a point which is not a critical point,

WENO5-Z scheme
Borges et al. [9] found that the increase of weights of less-smooth sub-stencils can substantially improve the numerical solution near both discontinuities and small-scale structures. Based on this, through the novel use of higher information that is already presented in the framework of classical WENO5-JS, they devised new smoothness indicators given as   Castro et al. [10] then fixed the problem of accuracy loss and extended WENO5-Z to higher order. For

The mapping method by Henrick
Henrick et al. [6] first pointed out the fact that WENO5-JS suffers accuracy loss at a critical point, and then they introduced an extra mapping process to recover the accuracy, where the weights are mapped to more accurate values. To be specific, the mapping is achieved through a set of mapping functions having the form of ( ) The function is smooth and has the following features: ( ) Then, new weight is obtained by Substitute Eq. (22) into Eq. (23), and new weight satisfies Eq. (24) indicates that the strong sufficient condition, Eq.(8), is satisfied. Therefore, the problem of accuracy loss at a critical point is fixed. The resultant mapped WENO scheme is generally abbreviated as WENO-M.

The Asymmetric improvement of WENO5-M
Here, the linear advection equation, 0, [ 1,1], uu x tx  + =  −  (25) with periodic boundary conditions, is taken as an example for illustration. Consider a single square wave, and the initial condition is se.
x uxt The third-order TVD Runge-Kutta method is used to do the discretization in time, expressed as Time is integrated to 20 t = , and the numerical results are shown in Fig. 2. Due to the less dissipation of WENO5-M compared to WENO5-JS, an improvement in the result near the discontinuity is observed in Fig. 2 (a). From the absolute pointwise errors in Fig. 2 (b), however, the improvement on both sides of the discontinuity is not in a symmetric fashion. Specifically, the error curve of WENO5-JS is almost symmetric about the discontinuity, shown with a red line in Fig. 2 (b), while WENO5-M has obvious larger error on the right side of the discontinuity, shown with a blue line in Fig. 2   In order to clarify the cause of this asymmetric improvement of WENO5-M, a simple model shown in Fig. 3 is used. On a six-point domain with a discontinuity at the center, there are two sets of five-point stencils symmetric about the discontinuity. A quantity with superscript L belongs to the left stencil, while with superscript R belongs to the right stencil.
Obviously, for the left stencil, 0-th sub-stencil is the smoothest, and 2-th sub-stencil is the least smooth. While for the right stencil, 2-th sub-stencil is the smoothest, and 0-th sub-stencil is the least smooth.
According to Eq.(5), the weight of each sub-stencil on the same stencil satisfies the relation, Assuming that , are weights assigned to the least-smooth sub-stencils of left and right stencils, the weight of each sub-stencil can be written as WENO scheme behaves like an ENO scheme at a discontinuity, and thus , are small values close to 0. Then, new weights are obtained through the mapping process, i.e., As , From Eq. (32), it is observed that after mapping, the magnification of weight of the least-smooth sub-stencil is not identical for both sides of the discontinuity, i.e.
, obviously larger on the right side. It is the mismatched magnification that results in the asymmetric improvement of WENO5-M around a discontinuity.

A symmetry-preserving mapping method
It is the derivative of mapping function at 0 that determines the magnification of weight of less-smooth sub-stencils, i.e., ( ) Due to the difference in ideal weights for each sub-stencil, mapping the weight with the function by Henrick inevitably causes the problem of asymmetric improvement. There are two ways to overcome the problem of mismatched magnification. One way is to redesign a set of mapping functions with identical derivatives at 0, and the other is to choose a new mapping object that has the same ideal value for each sub-stencil. An example of the former is the work in [7], where the new mapping function is smooth and piecewise, The new function has identical derivative at 0, In [7], the researchers noticed the distortion problem of WENO5-M, but they improperly attributed this to the numerical instability, which was thought to be caused by over-amplification of weight of less-smooth sub-stencil. Therefore, they proposed the function with derivative of 0 at 0, in the form of Eq. (34), to retain the ENO property at a strong discontinuity. Although the reason is not properly understood, due to restricting all derivatives to be identical at 0, using mapping function in the form of Eq.(34) is able to avoid the problem of asymmetric improvement. However, extra constraint on derivative brings difficulty in designing concise mapping function, and complicated function will result in significantly increased computational cost. The other way is choosing a new mapping object that has identical ideal value for each sub-stencil, so that extra constraint on derivative at 0 is not necessary. Furthermore, when the ideal value is the same, the mapping function is unique, and the function design is much simplified. Therefore, the second way is more reasonable.
In the formula of weight Eq. (5) Eqs.(43) and (24) indicate that at a critical point, mapping k  can achieve the same effect as mapping k  .
From Eqs. (28) and (29), one has Eq.(50) indicates that using the new mapping object k  , the problem of mismatched magnification of weights of less-smooth sub-stencils is fixed.
In this article, WENO5-JS with the symmetry-preserving mapping method is abbreviated as WENO5-FM.
As is shown in Fig. 2  where Eq.(34) is used to map k  , are also presented, 6 n = as suggested in [7]. The corresponding scheme is abbreviated as WENO5-MP. As is shown in Fig. 4, with longer output time 100 t = , the asymmetric improvement is more evident compared the case of 20 t = . On the left side of the discontinuity, three schemes have similar results, while on the right side, the result of WENO5-M is much worse than the others and suffers obvious distortion. As is mentioned before, to enhance the numerical stability near a discontinuity, the function used in WENO5-MP constrains the derivative at 0 to be 0, but numerical oscillation is observed in Fig. 5 and Fig. 6. To figure out the numerical instability, it is necessary to have an observation at the distribution of mapping functions.
As is shown in Fig. 7  It is observed that the magnification of weight near 0 is smaller with blue line, which is beneficial to capture strong discontinuities. However, for a weak discontinuity in numerical sense, i.e., the weight is not so close to 0 or 1 but deviates significantly from the ideal value, the adjustment of weight is obviously stronger with blue line.
From the results in Fig. 4, Fig. 5 and Fig. 6, it can be concluded that even with larger magnification of weight near strong discontinuity, the numerically stability of WENO5-FM is still ensured. By contrast, there is no oscillation observed in the result of WENO5-MP at 100 t = , while evident oscillation occurs at 250 t = when the initial strong discontinuity is smoothed out and become much weaker due to numerical dissipation. It can be inferred that, the numerical instability of WENO5-MP comes from the over-adjustment of weight near a weak discontinuity. As is mentioned previously, there is a controllable parameter n in the function defined by Eq. (34). Functions with 4 n = and 2 n = are plotted in Fig. 7, and it is shown that with smaller n , the adjustment of weight is weaker, especially in the region with weak discontinuities. Fig. 8 shows the results of WENO5-MP with 6,4,2 n = . In the results of 4,2 n = , there is no numerical oscillation observed, while 2 n = leads to a slightly more dissipative result.
This further supports the infer that the numerical instability is due to the over-adjustment of weight near a weak discontinuity. The above results suggest that weak discontinuity is also important for numerical stability. To ensure a stable mapping, excessive adjustment should be avoided near both strong and weak discontinuities. The core of mapping method lies in the mapping function, and different requirements on a scheme, such as strong stability, low dissipation, can be met with specifically designed function.   By introducing a global smoothness indicator, the relevance of less-smooth sub-stencil is increased in WENO5-Z, and thus less dissipation is achieved compared to WENO5-JS. However, it is shown in Fig. 10 that the distortion still exists in the results of WENO5-Z. Different from WENO5-JS, the distortion is already obvious at 100 t = for WENO5-Z, shown with a blue line in Fig. 10, while the distortion at 250 t = is less severe compared to WENO5-JS.

The distortion of WENO-Z in long-time simulation
The distortion problem of WENO5-JS and WENO5-Z, in the long-time simulation of linear advection equation, is also observed by Feng et al. [7]. By designing a new set of mapping functions that satisfy ( ) the distortion problem of WENO5-JS was fixed with an extra mapping process. Inspired by this, the mapping method can also be potential to overcome the distortion problem of WENO5-Z. However, based on the results in Section 3, it is necessary to adopt the symmetric-preserving mapping method.

Mapped WENO5-Z scheme
It is straightforward to do the symmetric-preserving mapping for WENO5-Z. Specifically, the new variable for WENO5-Z to map is Z 5 1. In this article, WENO5-Z with the symmetric-preserving mapping method is abbreviated as WENO5-ZM.
To demonstrate the effect of the mapping on WENO5-Z, the same numerical simulation in Section 4.1 is conducted. As is shown in Fig. 11, there is no more distortion observed in the results of WENO5-ZM at 100,250 t = , shown with the green and orange lines. With distortion problem fixed, the improvement near the discontinuity is significant compared to WENO5-Z. To further demonstrate the performance of WENO5-ZM in problems with multiple discontinuities, consider the following initial condition, ,  By contrast, the improvement of WENO5-ZM is much significant, clearly shown in the zoomed-in solutions with the green line. As is shown in Fig. 13 and Fig. 14  Following the approximated dispersion relation (ADR) analysis [29], the spectral properties of various schemes are shown in Fig. 15. With the symmetry-preserving mapping, the bias of overall scheme towards the linear scheme is increased. Thus, the spectral properties of WENO5-FM and WENO5-ZM are improved compared to WENO5-JS and WENO5-Z, respectively. Since WENO5-Z has better resolution than WENO5-JS in the intermediate-to-high wave-number range, WENO5-ZM shows lower dispersion and dissipation than WENO5-FM. Notice that WENO5-Z shows better spectral properties than WENO5-FM in the intermediate wave-number range, roughly from 1 However, as is shown in Fig. 12, Fig. 13 and Fig. 14, WENO5-FM produces much better results than WENO5-Z near discontinuities.

Pre-discrete mapping method
The disadvantage of mapping method is obvious, i.e., due to involving extra mapping process, the computational cost will be increased. As was reported in [10], WENO5-M led to 20%~30% extra cost compared to WENO5-JS. More complex function, like Eq.(34), will undoubtedly worsen this disadvantage.
In order to reduce the cost of mapping process, Hong et al. [30] proposed a pre-discrete mapping method.
The method consists of three steps, which are described as follows.
Step 1. Discretize the interval (0,1) , range of k  , with a uniform space 1 N  = , where N is the number of subintervals and satisfies 1 .
Step 2. The discrete sequence of k  can be expressed as , 0 , 0,..., Step 3. The index of the point closest to k  is obtained by where int denotes a rounding operation. Finally, the mapped value of k  is given by In fact, , k index  is a first-order approximation to k  , and in the smooth region it can be written as Also, Taylor series expansion of ( ) Eqs. (61) and (42) indicate that the pre-discrete mapping method can do the same mapping as the original method.
In the pre-discrete mapping method, the calculation of mapping function, which is the most time-consuming in mapping process, is only done once in step 1 before the iteration. Instead of calculating the function directly, the mapped value is obtained by one multiplication, one rounding operation and one addressing operation, which leads to lower computational cost. Moreover, this low cost is independent of the specific form of mapping function, and thus all kinds of functions can be used at the same cost, which frees the function design from the limit of computational cost.

Numerical results
In this section, several one-and two-dimensional Euler problems are conducted to validate the performance of WENO5-ZM. The Roe scheme [31] is used to solve a local Riemann problem at cell boundary. An explicit third-order TVD Runge-Kutta method, defined as Eq. (27), is used to solve the resulting set of ordinary differential equations in time. Referring to [14], the parameter e is set to be as small as 20 10 -, and the power parameter p is taken as 2 in all schemes. All numerical experiments are computed with double precision (64 bits). In the pre-discrete mapping method, N is taken as 4 10 for all cases.

One-dimensional Euler equations
The strong conservative form of one-dimensional Euler equations of gas dynamics is is supplemented to close the system of equations. g is the specific heat ratio and taken as 1.4 if not specified.

The Riemann problem of Sod
Sod problem [32] is used to demonstrate the ability of shock-capturing. The initial condition on

Mach 3 shock-density wave interaction
The initial condition is a Mach 3 shock interacting with a perturbed density field, and given as

Two-dimensional Euler equations
For two-dimensional problem, the strong conservative form of Euler equations is

Isentropic vortex propagation problem
The initial condition on a square domain [0,10] [0,10]  is given by Referring to [34], this test is tested under 10 periods and 100 periods with 100 t = and 1000 t = , respectively. After 10 periods propagation with 100 t = , it is shown in Fig. 23 (a) that the three results are very close to each other and in good agreement with the reference result. However, as is shown with the contours in Fig. 22, the difference becomes evident after 100 periods propagation with 1000 t = .
A quantitative comparison is given in Fig. 23 (b). It is observed that WENO5-JS and WENO5-ZM still keep the vortex at the center of the computational domain, while in the result of WENO5-Z, the vortex travels a little faster and thus is off center. Moreover, the result of WENO5-ZM is the closest to that of the optimal linear scheme but with slightly more dissipation.

The computational cost
Due to involving an extra mapping process, the computational cost of WENO5-ZM must be higher than that of WENO5-Z. If the improvement of WENO5-ZM comes at the cost of significantly increased computational time, then the attractiveness of the improved scheme will be greatly reduced. In order to give an intuitive comparison, the computational cost is numerically investigated here.
Take the DMR problem in Section 5.2.2 as the object to test, and to eliminate occasionality, the experiment is repeated three times under the same conditions. The computation is performed with double precision on a 2.40GHz Intel Xeon CPU E5-2676, and the code is developed in Fortran. In average, the computational cost of WENO5-ZM is about 6% higher than that of WENO5-Z. Compared to 20%~30% extra cost of WENO5-M reported in [10], such a low extra cost indicates the efficiency of WENO5-ZM.

Conclusions
In the original mapping method proposed by Henrick, the weight of each sub-stencil is used to map for recovering the accuracy at critical points. The numerical experiment of linear advection equation in this article shows that, the result around a discontinuity is improved but in an asymmetric fashion. With longer output time, this asymmetric improvement results in a solution with evident distortion. From a qualitative analysis, it is found that there is mismatched magnification of weight of less-smooth sub-stencil for both sides of a discontinuity, due to the difference in ideal weights of sub-stencils. To overcome the problem of asymmetric improvement, a symmetry-preserving mapping method is proposed in this article. Instead of the weight, a new variable related to the smoothness indicator is chosen as the mapping object, which has the same ideal value for each sub-stencil.
Moreover, since the function is unique and the ideal value is not close to 0, function design is much simplified in the new mapping method. In the long-time simulation of linear advection equation, it is found that the numerical solution of WENO5-Z suffers evident distortion, which seriously degrades the quality of solution. When applying the new mapping method to WENO5-Z, the distortion no longer exists, and thus the result is significantly improved.
Numerical experiments with one-and two-dimensional Euler equations are conducted to validate the performance of mapped WENO5-Z. In all problems, the results of mapped WENO5-Z show better resolution at both discontinuities and small-scale structures compared to WENO5-Z, and there is no oscillation observed in these results. The extra computational cost, brought by the mapping process, is moderate with a value of 6%, which is numerically investigated by the DMR problem.