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 a 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 (2r − 1) 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 (r + 1) 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 (2r − 1) 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 ε 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 ε on the convergence property of WENO scheme in detail, and their results showed that ε = 10 −2 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. ε = Ω(Δx m ), 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 lowdispersion 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 symmetrypreserving 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 symmetrypreserving 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 is used, and the positions of cell center and cell boundary are represented as so that 1 Δx ðh iþ1=2 − h i − 1=2 Þ is exactly equal to the first-order derivative of f at cell center i, i.e. f 0 i . Without loss of generality, the fifth-order WENO (WENO5) scheme is taken as an example to illustrate the WENO reconstruction procedure. In linear scheme, a fifthorder approximation of h(x) only needs a five-point stencil S 5 , while in WENO5, S 5 is divided into three sub-stencils {S 0 , S 1 , S 2 }, and each sub-stencil contains three adjacent grid cells, shown in Fig. 1. In the smooth region, WENO5 combines each sub-stencil with corresponding ideal weights to obtain a fifth-order approximation of h(x). 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 h(x), and the three sub-stencils of WENO5 givê where the superscript is used to distinguish different sub-stencils. h k i − 1=2 can be obtained via shifting each index in h k iþ1=2 by −1.Then theĥ iAE1=2 of stencil S 5 is obtained by a weighted average of h k iAE1=2 , In the smooth region, the weights approach their ideal values, which are

WENO5-JS scheme
In the classical fifth-order WENO (WENO5-JS), proposed by Jiang and Shu [5], the nonlinear weights are defined as ε is a small positive parameter to avoid denominator being zero, and the power parameter p ≥ 1 controls the amount of dissipation, usually taken as 2. β k is the smoothness indicator of k-th sub-stencil, defined as In the smooth region, in order to achieve optimal order of accuracy, i.e., satisfying the relation the weights of sub-stencils need to meet certain constraints. A strong sufficient condition is used in [6], given as and a weaker sufficient condition is derived in [14], Taylor series expansions of Eq.(6) at cell center x i give If f ' ≠ 0, one can have If f ' = 0 and f '' ≠ 0 (critical point), one can have If β k in Eq.(6) takes the form of where D is a non-zero constant independent of sub-stencils, and m is a positive parameter. Then Eqs. (5) and (13) give and which lead to 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, ω k obtained by WENO5-JS satisfies the weaker sufficient condition for achieving optimal order, i.e. Eq. (9). At a critical point, ω k obtained by WENO5-JS only satisfies ω k = d k + Ο(Δx), resulting in the accuracy loss at this 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 where τ 5 = | β 0 − β 2 | is the higher order smoothness indicator. τ 5 has two important properties, which are If S 5 is smooth, then τ 5 = Ο(Δx 5 ). If S 5 is discontinuous, then τ 5 = Ο(1) ≤ max(β 0 , β 1 , β 2 ).
Substitute Eq.(17) into Eq.(5), and the new weights are obtained by which lead to the fifth-order WENO-Z (WENO5-Z). At a critical point with f ' = 0 and f '' ≠ 0, β k = Ο(Δx 4 ) in Eq.(10), leading to τ 5 β k þε ¼ ΟðΔxÞ and thus ω k = d k + Ο(Δx). Therefore, Eq.(18) does not solve the problem of accuracy loss at a critical point. However, due to the larger assignment of weights to the less-smooth sub-stencils, the numerical results of WENO5-Z are substantially improved compared to WENO5-JS, and slightly better compared to WENO5-M, details shown in [9].
Castro et al. [10] then fixed the problem of accuracy loss and extended WENO5-Z to higher order. For (2r − 1) th-order WENO-Z, the weights are obtained by The power parameter p ≥ 1 in Eq. (19) plays two roles. On one hand, p controls the amount of numerical dissipation, and the scheme is more dissipative with a larger p. On the other hand, the accuracy at a critical point can be recovered by increasing p. The former is intuitive and easy to understand. As for the latter, at a critical point, ð τ 5 β k þε Þ p ¼ Ο ðΔx p Þ leads to ω k = d k + Ο(Δx p ), and thus the weaker sufficient condition, Eq.(9), for achieving optimal order is satisfied if p ≥ 2. Changing the value of p to alter convergence rate at critical points is a unique feature of WENO-Z, while in WENO-JS and WENO-M, p only controls the amount of numerical dissipation. In order to balance accuracy and stability, for (2r − 1) th-order scheme, p is generally suggested to be r − 1 [10].
3 The symmetry-preserving mapping method

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: Once the weight of WENO5-JS in Eq.(5) is obtained, the new unnormalized weight is given by Taylor series expansion of Eq. (21) gives 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, with periodic boundary conditions, is taken as an example for illustration. Consider a single square wave, and the initial condition is The third-order TVD Runge-Kutta method is used to do the discretization in time, expressed as Time is integrated to t = 20, 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 (b). Consequently, as is shown in Fig. 2 (a), the improvement of WENO5-M on the right side of the discontinuity is obviously less than that on the left side. Due to the asymmetry near the discontinuity, the whole error curve of WENO5-M is no longer symmetric about the center of domain x = 0, shown 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.
According to Eq.(6) and due to the symmetry, smoothness indicators of left and right stencils have the following relation, 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 θ, σ approach 0, it is approximately obtained that 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, where Fig. 3 Schematic of the analysis model and ; else: 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), besides the weight ω k , there is another quantity related to each sub-stencil, i.e. the smoothness indicator β k . In the smooth region, different from the weight ω k , the smoothness indicator β k is equal for each sub-stencil, and thus its ideal value is unique. In order to keep proportional to ω k , consider the inverse of smoothness indicator with power parameter p, i.e. 1 ðβ k þεÞ p . But for mapping, it is necessary to normalize them to a limited interval, The normalized variable, λ k , is the new mapping object, and its unique ideal value is In general, the ideal value of λ k for (2r − 1) th-order WENO scheme is 1 r . Substitute Eq.(12) into Eq.(38), and it is obtained at a critical point that It is straightforward to apply the mapping function, originally designed for ω k , to do the mapping for λ k . Considering the function by Henrick, the mapping function for λ k is unique and written as Similar to Eq.(22), after mapping, it is obtained that Then, new weight is evaluated and satisfies Eqs. (43) and (24) indicate that at a critical point, mapping λ k can achieve the same effect as mapping ω k .
According to Eqs. (28) and (38), λ k for the left and right stencils in Fig. 3 satisfies After mapping with function Eq.(41), the mapped value still satisfies According to Eq.(43), it is obtained that ; From Eqs. (28) and (29), one has According to Eqs.(33), (47), (48) and (49), it is obtained that 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 (b), the error curve of WENO5-FM is symmetric about the discontinuity, and almost symmetric about the center x = 0 on the whole domain. Consequently, compared to WENO5-JS, the improvement of WENO5-FM on both sides of a discontinuity is in a symmetric fashion. To present an intuitive comparison, the linear advection equation initialized by Eq. (26), is still numerically tested here. Results of the first way, where Eq.(34) is used to map ω k , are also presented, n = 6 as suggested in [7]. The corresponding scheme is abbreviated as WENO5-MP. As is shown in Fig. 4, with longer output time t = 100, the asymmetric improvement is more evident compared the case of t = 20. 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.
It is observed from Fig. 5 that at t = 250, both sides of the discontinuity are influenced in the result of WENO5-M, and similar to the right side of the discontinuity, the result on the left side also becomes distorted. On the right side of the discontinuity, the results of WENO5-MP and WENO5-FM are still close to each other, similar to the case of t = 100. However, on the left side, there is obvious overshoot in the result of WENO5-MP, shown with a blue line in Fig. 5. Taking the result at t = 400 in Fig. 6 into consideration together, it can be inferred that the overshoot of WENO5-MP is caused by numerical instability. From the results at t = 100, 250, 400, WENO5-FM is always able to obtain non-oscillatory results and the improvement is symmetric about the discontinuity.
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  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 Figs. 4, 5 and 6, it can be concluded that even with larger magnification of weight near strong discontinuity, the numerical stability of WENO5-FM is still ensured. By contrast, there is no oscillation observed in the result of WENO5-MP at t = 100, while evident oscillation occurs at t = 250 when the initial strong discontinuity is smoothed out and becomes 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 n = 4 and n = 2 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. Figure 8 shows the results of WENO5-MP with n = 6, 4, 2. In the results of n = 4, 2, there is no numerical oscillation observed, while n = 2 leads to a slightly more dissipative result. This further supports the inference 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. In terms of function design, using λ k to map has two main advantages over using ω k . On one hand, for (2r − 1) th-order scheme, there are r different ideal weights and thus r different functions need to be designed. While there is only one function to be designed for using λ k to map due to the unique ideal value, which greatly simplifies the function design. On the other hand, the minimum ideal weight is closer to 0 with higher order, e.g., d 0 ¼ 1 10 ; 1 35 ; 1 126 for fifth-, seventh-and ninth-order schemes, respectively. Ideal value close to 0 brings difficulty in balancing accuracy recovery and numerical stability. By contrast, λ ¼ 1 3 ; 1 4 ; 1 5 for fifth-, seventh-and ninth-order schemes respectively, which significantly alleviates this dilemma.  By introducing a global smoothness indicator, the relevance of less-smooth substencil 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 t = 100 for WENO5-Z, shown with a blue line in Fig. 10, while the distortion at t = 250 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 g 0 k ð0Þ ¼ 0 , 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 symmetricpreserving 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 The ideal value of λ Z k is the same as λ, i.e., Therefore, the mapping function is also the same as Eq.(41). Through mapping, it is obtained that Then, the new weight is obtained by 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 t = 100, 250, 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, q ; z = − 0.7, δ = 0.005, β ¼ log2 36δ 2 ; a = 0.5, α = 10. The initial condition consists of a Gaussian function, a square function, a piecewise linear triangle function and an ellipse function. It is shown in Fig. 12 that at a short time t = 4, there is no evident distortion in the results of WENO5-JS and WENO5-Z, and WENO5-Z has slight improvement near a discontinuity compared to WENO5-JS. 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 Figs. 13 and 14, at t = 100, 250, there is evident distortion in both results of WENO5-JS and WENO5-Z, and it is more severe with WENO5-JS. As expected, all the results of WENO5-ZM free from the distortion problem, and significantly improved results near discontinuities are observed. For comparison, the result of WENO5-FM is also presented in this test case. As is shown in these figures, WENO5-FM produces results very similar to that of WENO5-ZM, Specifically, WENO5-FM shows slightly more dissipation compared to WENO5-ZM, see the zoomed-in solutions.
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 wavenumber range, roughly from ξ = 1 to ξ = 1.5. However, as is shown in Figs. 12, 13 and 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 prediscrete 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 Step 2. The discrete sequence of λ k can be expressed as λ k, j = 0 + jΔλ, j = 0, …, N, and corresponding discrete sequence of mapped value can be obtained by 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 g(λ k, index ) at λ gives Equations (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 ε is set to be as small as 10 −20 , 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 10 4 for all cases.

One-dimensional Euler equations
The strong conservative form of one-dimensional Euler equations of gas dynamics is where ρ, u, p, E represent the density, velocity, pressure and total energy, respectively. The equation of state, is supplemented to close the system of equations. γ is the specific heat ratio and taken as 1.4 if not specified.

The Riemann problem of Sod
Zero-gradient boundary conditions are imposed at x = ± 0.5, and the numerical solutions are shown in Fig. 16. The solution consists of a left-going rarefaction wave, a central contact discontinuity and a right-going shock wave. All results are non-oscillatory, and due to less dissipation, WENO5-ZM shows better resolution near discontinuities than WENO5-Z, clearly observed in the zoomed-in solutions.

Mach 3 shock-density wave interaction
The initial condition is a Mach 3 shock interacting with a perturbed density field, and given as At x = ± 5, the zero-gradient boundary conditions are applied, and the numerical solutions are shown in Fig. 17. It is observed that in regions with discontinuities and high-frequency wave, WENO5-ZM shows the best resolution. While due to the excessive dissipation, the high-frequency wave is nearly smoothed out in the result of WENO5-JS.

Interacting blast waves
The test consists of two interacting blast waves. The strong shocks in the solution are computationally hard to solve, and schemes with unstable tendencies often fail to converge in this test. The initial condition on domain [0, 1] is ρ; u; p ð Þ¼ 1; 0; 1000 ð Þ ; x < 0:1; 1; 0; 100 ð Þ ; x > 0:9; 1; 0; 0:01 ð Þ ; else: Reflective boundary conditions are used at x = 0 and x = 1, and the numerical solutions are shown in Fig. 18. Obviously shown in the zoomed-in solutions, the result of WENO5-ZM shows the best resolution near discontinuities without any oscillation.

Two-dimensional Euler equations
For two-dimensional problem, the strong conservative form of Euler equations is where y, v represent the other direction and the velocity component in this direction, respectively. The two-dimensional equation of state is

Rayleigh-Taylor instability problem
The Rayleigh-Taylor instability (RTI) occurs at the interface between two fluids with different density when the heavier fluid accelerates to the lighter fluid. The computational domain is [0, 0.25] × [0, 1], and the initial condition is , is the speed of sound. The interface of two fluids with different density is at y = 0.5. The specific heat ratio in Eq.(68) is specified as γ = 5/3 in this case. At x = 0 and x = 0.25, reflective boundary conditions are used. ρ = 2, u = v = 0, p = 1 are set for y = 0, while ρ = 1, u = v = 0, p = 2.5 are set for y = 1. To exert the effect of gravity, the source terms ρ, ρv are explicitly added to the right hand side of third and fourth equations of Eq.(67).
As is shown in Fig. 19, all schemes can obtain symmetric results. Due to the Kelvin-Helmholtz instability, vortices are rolled up around the interface indicated by clustered contour lines. The least dissipation of WENO5-ZM leads to a result with the most details of vortices, while small-scale vortices are nearly smoothed out in the result of WENO5-JS due to its over-dissipation.

Double Mach reflection problem
Double Mach reflection (DMR) describes a Mach 10 moving shock with an angle of 60 ∘ with respect to the wall. The computational domain is [0, 4] × [0, 1], and the initial condition is Fig. 19 Contour lines of density of the RTI problem at t = 1.95. This figure is drawn with 15 density contours between 0.9 and 2.3. A uniform grid with 120 × 480 cells is used, Δx ¼ Δy ¼ where x 0 = 1/6 is the leading position of the wall, and θ = π/6 is the angle between initial oblique shock and y-axis. Supersonic inflow and outflow boundary conditions are used at x = 0 and x = 4 respectively. At the lower boundary y = 0, post-shock values are used at (0, x 0 ), while a reflective boundary condition is used at (x 0 , 4). As for the upper boundary y = 1, the exact solution of the Mach 10 moving oblique shock is applied. It is observed from Fig. 20 that all schemes can give results without any oscillation, and the main difference is reflected near the slip line. Vortices rolled up around the slip line can be observed in the result of WENO5-Z, while in the result of WENO5-ZM, they are more obvious, indicating the improvement of the latter.

2D Riemann problem
The 2D Riemann problem [33] Figure 21 shows the density contour obtained with WENO5-JS, WENO5-Z and WENO5-ZM. It is generally thought that capturing more small-scale structures indicates better resolution in this problem. Therefore, from the comparison of the three results, WENO5-ZM obviously has the best resolution.

Isentropic vortex propagation problem
The initial condition on a square domain [0, 10] × [0, 10] is given by where κ = 5 is the vortex strength, and (x c , y c ) = (5, 5) is the initial location of the vortex center. Periodic boundary conditions are used for all boundaries. Referring to [34], this test is tested under 10 periods and 100 periods with t = 100 and t = 1000, respectively. After 10 periods propagation with t = 100, 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 t = 1000. 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%

Conclusions
In the original mapping method proposed by Henrick, the weight of each substencil 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.