3.1. Motivation
Consider the one-dimensional linear advection problem
which represents one of the most straightforward hyperbolic problems. Given the initial condition
, its exact solution is simply
, which makes it a suitable starting point for the assessment of numerical schemes. It must be remarked that the Euler equations reduce to linear advection of density in the absence of velocity and pressure gradients.
Averaging Equation (20) over the
ith cell results in the ordinary differential equation (ODE) for the cell average
Notice that the right-hand side of Equation (21) involves the left-biased approximations of
u(
x) at the cell interfaces, which can be obtained through the WENO schemes described in the previous sections. In the present study, this ODE was solved in time using the explicit third-order TVD Runge–Kutta time-marching scheme [
13].
With the numerical schemes defined, consider the initial condition of two constant states separated by discontinuities at
x = 0 and
x = ±0.5
This problem was solved on the periodic domain
discretized into
N = 100 uniform cells using fifth- and seventh-order WENO schemes until
t = 100 s (100 cycles) using CFL = 0.1. Note that the mapped WENO methods using
gM(
ω) and
are denoted by WENO-M and WENO-IM(2,0.1), respectively. The final profiles are shown in
Figure 3.
It can be observed from
Figure 3a that at fifth order, the mapped WENO schemes WENO-M and WENO-IM(2,0.1) are able to capture the discontinuities with lesser numerical diffusion compared to the WENO-JS scheme. Recall that mapped WENO methods were developed to achieve optimal accuracy near smooth critical points. Yet, they seem to be able to capture sharper discontinuities than the original unmapped WENO-JS scheme. Between the mapped WENO schemes, WENO-IM(2,0.1) performs better as it captures the flat regions on both sides of the discontinuities more accurately. In contrast, WENO-M suffers a loss of accuracy to the right of the discontinuities and performs similar to WENO-JS scheme. On the other hand, the results become different at seventh order, as seen from
Figure 3b. Both the mapped WENO schemes deliver a similar, and at some regions even worse, performance compared to WENO-JS scheme. This can be determinantal when simulating the propagation of discontinuities and sharp gradients over long distances.
WENO-M and WENO-IM(2,0.1) both suffer a particularly significant loss in accuracy at the flat regions slightly to the left of the discontinuities. To better understand this loss in accuracy, the mapped weights were analysed in detail. The mapped weights
of the leftmost and rightmost stencils, i.e.,
j = 0 and
j = 2, at the end of the first cycle (
t = 1 s) are plotted for fifth-order WENO-M and WENO-IM(2,0.1) schemes in
Figure 4. Comparing
Figure 4a,b, a common trend can be discerned—the leftmost stencil weight
is larger than the rightmost stencil weight
in the region to the left of the discontinuity and vice versa to the right. Although the regions on either side of the discontinuity appear smooth and uniform, they are in fact highly non-smooth in the sense that the “flatness” of the profile changes vastly between one stencil to the next. This results in the smoothness indicators
ISj of adjacent stencils being an order of magnitude apart or more. Since a single WENO reconstruction involves multiple stencils (three for fifth order, four for seventh order), the resulting minimum and maximum values of
ISj are several orders of magnitudes apart. This explains the behaviour of the weights mentioned earlier. Interestingly, across the discontinuity itself, the weights remain close to the optimal weights
dj since the discontinuity has been sufficiently “smoothened” by numerical diffusion.
A close examination of the weights near discontinuities suggests that the fifth-order WENO-IM(2,0.1) scheme allocates
smaller weights to the non-smooth stencils (
j = 2 and
j = 0 to the left and right of the discontinuity, respectively) than the fifth-order WENO-M scheme at the end of first cycle, as shown in
Figure 4a,b. This is quite unexpected since WENO-IM(2,0.1) scheme generally allocates
greater weights to non-smooth stencils than WENO-M scheme, as evident from the shapes of their mapping functions shown in
Figure 2. The only conclusion is that WENO-M scheme produces a smoother, more diffused profile close to and across the discontinuity compared to WENO-IM(2,0.1) scheme at the end of the first cycle. This poses the question as to why the fifth-order WENO-IM(2,0.1) scheme captures a much sharper discontinuity compared to the fifth-order WENO-M scheme. It could be related to the behaviour of the mapping process across the discontinuity itself. Closer examination of the mapped weights across the discontinuity revealed that the mapped weights from the fifth-order WENO-IM(2,0.1) remained very close to their optimal values while those from the fifth-order WENO-M scheme deviated slightly. This is due to
having a wider flat portion in the neighbourhood of
than
(see
Figure 2). As a result,
attempts to map a wider range of weights to the optimal weights. Therefore, it is postulated that the fifth-order WENO-IM(2,0.1) scheme captures a sharper discontinuity because it is quicker to revert to the upstream central scheme at the smeared discontinuity and because it amplifies small weights
to a greater extent in general.
Moving on to seventh-order WENO schemes, the mapped WENO weights of the leftmost (
j = 0) and the rightmost (
j = 3) stencils WENO-M and WENO-IM(2,0.1) schemes at the end of the first cycle (
t = 1 s) are similarly shown in
Figure 5. The behaviour of the seventh-order mapped weights is similar to the fifth-order ones across the discontinuity and to its right. To the left of the discontinuity, however, the mapped weights
are relatively larger; for WENO-IM(2,0.1), they are even comparable to
at some locations. This is a crucial difference between the fifth- and seventh-order mapped WENO schemes. It is believed that this difference is due to the small value of
. Referring to the close-up view in
Figure 2c, it can be seen that both mapping functions
and
result in particularly severe amplifications to the small weights compared to the identity map (WENO-JS). In fact, the slope of
at
ω = 0
is strongly proportional to the inverse of
d. Feng et al. [
7] have pointed this out and recommended the use of the WENO-M scheme instead of the WENO-IM(2,0.1) scheme for seventh order since
. However, it is clear from
Figure 3b and
Figure 5 that the WENO-M scheme, while somewhat better than the WENO-IM(2,0.1) scheme, is still inadequate in overcoming the amplification problem for seventh order.
3.2. Proposed Solution
Based on the discussion in the previous section, it appears that limiting the amplification of small weights is important for seventh-order mapped WENO schemes. This idea was pursued in a couple of earlier studies [
14,
15] by ensuring that the slope of the mapping function vanishes at
ω = 0. This actively suppresses near-zero weights. In the present paper, a different, rather simpler approach is suggested to overcome the amplification problem. Consider a general family of rational mapping functions
given as
where
k is a positive even integer,
m is a positive integer and
s is a positive scaling factor. The corresponding mapped WENO scheme will be referred to as WENO-RM(
k,
m,
s).
gRM includes
gIM as a special case for
m = 1 and
. Since
is a positive real value less than or equal to 1,
ω(1 −
ω) is a positive quantity with an upper bound of
. Therefore, increasing the value of
allows the term [
ω(1 −
ω)]
m in the denominator to diminish and, therefore, causes the mapping function to move closer to the identity map
g(
ω)
= ω. In fact, for
m > 1,
, i.e.,
gRM functions are able to closely follow the identity map near
ω = 0 as shown in the inset plots in
Figure 6. This means that, in contrast to previous studies in which small weights were actively suppressed, small weights are preserved with little or no amplification upon mapping with
gRM functions. Therefore, for
m > 1, the proposed mapping function
gRM is more suitable for mapped WENO schemes at seventh order and above.
Increasing the values of
k and
s have the opposite effect of increasing
m. Increasing the value of
k causes more derivatives of
gRM to vanish at
ω =
d. Similarly, increasing the value of
s would increase the denominator of the rational expression causing the constant term
d to dominate. In both cases, the flat portion of the mapping function in the neighbourhood of
ω =
d is widened, thereby extending the range over which weights get mapped to optimal weights. Feng et al. [
7] did not investigate
mapping functions beyond
k = 2 because increasing
k severely amplifies the small weights (see Equation (23)). However, the introduction of parameter
m overcomes this restriction and, therefore, a greater range of
k values can be explored for mapping functions, which couple better with higher-order WENO methods with
m > 1.
Before assessing the performance of any mapping function, it is crucial to verify that it satisfies all the three properties (a)–(c) mentioned earlier. It can be easily verified that
gRM satisfies properties (a) and (c). However, it does not satisfy property (b) unconditionally.
gRM can be non-monotone and determining the exact conditions for monotonicity is not very straightforward. Nevertheless, it can be shown that
k ≥
m − 1 is a sufficient condition for monotonicity regardless of the value of
s and
d (refer to
Appendix A). Therefore, the present discussion will be restricted to
k ≥
m − 1. It is useful to mention that such restriction is considered reasonable because there are valid values of
k that satisfy the monotonicity condition for every value of
m.
3.3. Parametric Study
A parametric study was conducted to determine good combinations of parameters (
k,
m,
s) using a simple one-dimensional linear advection problem [Equation (20)] with the initial condition given as
where
,
and
δ = 0.005. From left to right, the initial profile consists of a Gaussian, a square wave, a triangle wave and half-ellipse. This initial condition was first used by Jiang and Shu [
5] in their classic paper on WENO schemes. Henceforth, it has become a popular choice for benchmarking numerical methods for hyperbolic problems.
The problem was solved on the domain
at two resolutions,
N = 400 and 800 cells. The simulation was run until
t = 200 s (100 cycles) using CFL = 0.1. The parameter combinations tested are listed in
Table 1. The corresponding
L1 error norms are plotted in
Figure 7 and
Figure 8 for
N = 400 and 800 cells, respectively.
A consistent trend can be observed in all cases. For any given combination of
k and
m, the errors tend to be large for small values of
s. Then, the errors reduce with increasing
s before increasing slightly at very large values of
s. The explanation for this trend is quite intuitive—increasing the value of
s extends the range of weights which are mapped to their optimal values. Too small a value of
s would not map the weights to optimal weights ‘rapidly’ enough resulting in greater numerical diffusion (see discussion in
Section 3.1). On the other hand, too large a value would cause the weights to be mapped to optimal weights over-aggressively producing spurious oscillations. This effect is illustrated in
Figure 9 for
k = 4 and
m = 2 at three different values of
s. When
s = 2, all four profiles suffer distortions. Upon increasing
s to 200, the distortions are nearly completely eliminated, and all four profiles are captured well. A small spurious oscillation occurs at the foot of the half-ellipse at
x ≈ 0.6. When increased further to
s = 2 × 10
4, prominent spurious oscillations form at the square wave and at the foot of the triangle wave.
Another notable trend is that increasing the value of m almost always results in larger errors at lower values of s. Increasing m causes the mapping function to follow the identity map more closely and decreases the range over which weights are mapped to optimal weights. As such, a large value of s is required to counteract its effect and widen the range before the errors start to reduce. Finally, increasing the value of k tends to have a favourable effect, suggesting that a wide flat region about ωj = dj is crucial with the only exception taking place when k is increased from 4 to 6 with m = 2.
Among the various combinations of parameters
k,
m and
s tested,
k = 6 and
m = 3 consistently results in very low errors over the entire range of
s values at both resolutions. However, a suitable value for
s appears to be dependent on the problem and resolution. For instance, among the values of
s tested,
s = 2 × 10
3 results in the smallest
L1 error norm when using
N = 400 cells. It reduces to
s = 2 when using
N = 800 cells. However, the choice of the value of
s may not be too critical since it is evident from
Figure 7 and
Figure 8 that the error norms remain small over a wide range of values of
s when
k = 6 and
m = 3.
To conclude this section, let us return to linear advection of sharp discontinuities, which motivated the development of this new family of mapping functions.
Figure 10a,b show the mapped weights for this problem at
t = 1 s (1 cycle) for two different values of
s with
k = 6 and
m = 3. Notice that in both cases, the mapped weights behave in the desired manner without any severe amplification in the flat regions on either side of the smeared discontinuity. It can also be seen that the mapped weights for the non-smooth stencils (
j = 0 to the left of the discontinuity and
j = 3 to its right)
reduce when
s is increased. This indicates that the discontinuity is captured with less dissipation when using
s = 2 × 10
3. The final profiles at
t = 100 s (100 cycles) using both values of
s are shown in
Figure 10c. In contrast to the results shown in
Figure 3b, the discontinuities are captured well without any distortions for both values of
s even at relatively long output times. The close-up view shows that using
s = 2 × 10
3 indeed results in a slightly sharper discontinuity compared to
s = 2. The value of
s was chosen to be 2 × 10
3 for all the test cases reported in the following section with the exception of the double-Mach reflection Euler problem for which a value of
s = 2 was chosen to ensure a stable run.