- Email: [email protected]
An improved mapped weighted essentially non-oscillatory scheme
Applied Mathematics and Computation 232 (2014) 453–468
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
An improved mapped weighted essentially non-oscillatory scheme Hui Feng a, Cong Huang a, Rong Wang b,⇑ a b
School of Mathematics and Statistics, Wuhan University, Wuhan, Hubei, PR China Department of General Education, South University of Science and Technology of China, Shenzhen, Guangdong, PR China
a r t i c l e
i n f o
Keywords: WENO methods Hyperbolic conservation laws Mapping function
a b s t r a c t In this paper we develop an improved version of the mapped weighted essentially non-oscillatory (WENO) method in Henrick et al. (2005) [10] for hyperbolic partial differential equations. By rewriting and making a change to the original mapping function, a new type of mapping functions is obtained. There are two parameters, namely A and k, in the new mapping functions (see Eq. (13)). By choosing k ¼ 2 and A ¼ 1, it leads to the mapping function in Henrick et al., i.e.; the mapped WENO method by Henrick et al. actually belongs to the family of our improved mapped WENO schemes. Furthermore, we show that, when the new mapping function is applied to any ð2r 1Þth order WENO scheme for proper choice of k, it can achieve the optimal order of accuracy near critical points. Note that, if only one mapping is used, the mapped WENO method by Henrick et al., whose order is higher than five, can not achieve the optimal order of accuracy in some cases. Through extensive numerical tests, we draw the conclusion that, the mapping function proposed by Henrick et al. is not the best choice for the parameter A. A new mapping function is then selected and provides an improved mapped WENO method with less dissipation and higher resolution. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction The essential non-oscillatory (ENO) method was studied by Harten et al. in [6,8,9,7]. It owes its success to the ENO property which is considered to be very important in the numerical simulation for hyperbolic conservation laws because high order numerical methods often produce spurious oscillations, especially near shocks or other discontinuities. The finite-volume ENO spatial discretization was studied by Harten et al. in [7], later the finite-difference ENO method was developed by Shu et al. in [16,17]. Later weighted ENO (WENO) methods have been introduced in [12,11] to address potential numerical instabilities due to the process of choosing ENO stencils [13]. Unlike ENO methods, a convex combination of all the ENO candidate stencils are used in WENO methods such that the numerical flux is approximated with the higher order of accuracy in smooth regions, while WENO methods still retain the ENO property in regions near discontinuities. In this paper we will use WENO-JS to refer the WENO scheme in [11,1]. See [14,15] for an overview and other references. Recently Henrick et al. [10] noted that the fifth-order WENO method proposed by Jiang and Shu [11] is only third-order accurate near critical points of the smooth regions in general. By using a simple mapping function to the original weights in ⇑ Corresponding author. E-mail addresses: [email protected] (H. Feng), [email protected] (C. Huang), [email protected] (R. Wang). http://dx.doi.org/10.1016/j.amc.2014.01.061 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
454
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
[11], Henrick et al. developed a mapped WENO method to achieve the optimal order of accuracy; i.e., the fifth order near critical points. However, when this mapping function is applied to the WENO scheme of the order P 7, it can not achieve the optimal order of accuracy near critical points. In this paper we use WENO-HAP to refer the family of mapped WENO schemes in [11]. Lately Feng et al. [4] studied the mapped WENO method in [10] and found that, when it is used for solving the problems with discontinuities, the mapping function in [10] may amplify the effect from the non-smooth stencils and thus cause a potential loss of accuracy near discontinuities. This effect may be difficult to be observed for the fifth-order WENO method unless a long time simulation is desired. In order to overcome this problem, a new piecewise polynomial mapping function was proposed in [4] by adding two requirements, that is, g 0 ð0Þ ¼ g 0 ð1Þ ¼ 0 to the original criteria in [10]. Clearly the requirement g 0 ð0Þ ¼ 0 tends to decrease the effect from the non-smooth stencils and g 0 ð1Þ ¼ 0 tends to increase the effect from the smoothest stencil. In this paper we use WENO-PM to refer the family of mapped WENO schemes in [4]. Through a further study of the mapped WENO method, we find that, the mapping function in [10] can be written in a simple and more meaningful form. Make a change to the new form, then a new family of mapping functions, which includes the one in WENO-HAP, is obtained. This kind of new mapping functions still satisfies the original criteria in [10]. Furthermore, when the new mapping function is applied to any ð2r 1Þth order WENO scheme for proper choice of k (see Theorem 2 for more details), it can achieve the optimal order of accuracy near critical points. Note that, if only one mapping is used, the mapped WENO method by Henrick et al., which order is higher than five, can not achieve the optimal order of accuracy in some cases. In order to overcome this problem, two or more mappings need to be used. However, in [4] Feng et al. have pointed out, the more mappings it uses, the less accurate numerical solution it provides due to potential numerical instability. Extensive numerical tests show that, the mapping function in WENO-HAP is not the best mapped WENO method in this new family. Then a new mapping function is selected. It provides an improved mapped WENO method with less dissipation and higher resolution. The remainder of this paper proceeds as follows. In Section 2, we give a brief description of WENO-JS [11] and WENO-HAP [10]. In Section 3, we give the improved version of the mapped WENO-HAP method. In Section 4, we compare the performance of different WENO methods. 2. The fifth-order WENO method 2.1. WENO spatial discretization Consider the one-dimensional scalar hyperbolic conservation law
ut ¼ fx ðuÞ;
a < x < b; t > 0:
ð1Þ
Assume a uniform spatial mesh; i.e., xj ¼ a þ jDx; j ¼ 0; 1; . . . ; N, where Dx ¼ ba . The spatial derivative fx ðuÞ at the point N x ¼ xj is approximated using a conservative finite difference scheme
LðuÞ ¼
1 ^ F jþ1 F^j1 ; 2 2 Dx
ð2Þ
where F^j1 is the numerical flux at xj1 . The numerical flux must be consistent with f; i.e., F^ðu; . . . ; uÞ ¼ f ðuÞ. The numerical 2 2 flux determines the particular numerical method and its properties. The fifth-order WENO scheme [11] has a numerical flux which is a convex combination of three third-order fluxes associated with three substencils; i.e.,
F^jþ1 ¼ x0 F^0jþ1 þ x1 F^1jþ1 þ x2 F^2jþ1 ; 2
2
2
ð3Þ
2
where
F^0jþ1 2
F^1jþ1 2
F^2jþ1 2
2f j2 7f j1 þ 11f j ; ¼ 16 fj1 þ 5f j þ 2f jþ1 ; ¼ 16 2f j þ 5f jþ1 fjþ2 : ¼
1 6
ð4Þ
In smooth regions, the WENO approximation is designed such that it approximates the solution with the high order of accuracy. On the other hand, near discontinuities, it keeps the essentially non-oscillatory property of the ENO scheme. The nonlinear weights are computed as
a xi ¼ P2 i
J¼0 aJ
;
ai ¼
di ð þ ISi Þ2
;
i ¼ 0; 1; 2;
ð5Þ
where is a small positive number that is introduced to prevent the denominator becoming zero, and the optimal weights, di ; i ¼ 0; 1; 2, are given by
d0 ¼ 1=10;
d1 ¼ 6=10;
and d2 ¼ 3=10:
ð6Þ
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
455
ISi are the smoothness indicators, which are defined as follows
2 1 2 13 fj2 2f j1 þ fj þ fj2 4f j1 þ 3f j ; 12 4 2 1 2 13 IS1 ¼ fj1 2f j þ fjþ1 þ fj1 fjþ1 ; 12 4 2 1 2 13 IS2 ¼ fj 2f jþ1 þ fjþ2 þ 3f j 4f jþ1 þ fjþ2 : 12 4 IS0 ¼
ð7Þ
Henrick et al. [10] then introduced a mapping function of the form
ðg HAP Þi ðxÞ ¼
x di þ ðdi Þ2 3di x þ x2
;
2
ðdi Þ þ ð1 2di Þx
ð8Þ
where di are the optimal weights given in (6), x 2 ½0; 1 and i ¼ 0; 1; 2. This function is monotonically increasing in ½0; 1 with 0 00 finite slopes, ðg HAP Þi ð0Þ ¼ 0; ðg HAP Þi ð1Þ ¼ 1, ðg HAP Þi ðdi Þ ¼ di , and ðg HAP Þi ðdi Þ ¼ ðg HAP Þi ðdi Þ ¼ 0. The mapped weights are given by:
aHAP xHAP ¼ P2 i ; aHAP ¼ ðg HAP Þi ðxi Þ; i ¼ 0; 1; 2; i i HAP J¼0 aJ
ð9Þ
where xi are computed by (5) and (7). In smooth regions, this mapped WENO scheme; i.e., the WENO-HAP scheme, gives the fifth order of accuracy even near critical points. We refer to [10] for more details. 2.2. Time discretization After the spatial derivative is discretized with one of the WENO schemes, an Ordinary Differential Equation (ODE) system is obtained as follow
du ¼ LðuÞ; dt
ð10Þ
where the spatial operator LðuÞ is represented in (2). In all the numerical examples in this paper, the ODE system (10) is solved using the explicit, third-order, TVD, Runge–Kutta method [16,5] as follows
uð1Þ ¼ un þ DtLðun Þ; 3 1 1 uð2Þ ¼ un þ uð1Þ þ DtLðuð1Þ Þ; 4 4 4 1 n 2 ð2Þ 2 nþ1 u ¼ u þ u þ DtLðuð2Þ Þ: 3 3 3
ð11Þ
3. The new mapped WENO methods 3.1. A new mapping function The mapping function (8) can be written in a simple and more meaningful form:
ðg HAP Þi ðxÞ ¼ di þ
ðx di Þ
3
2
ðx di Þ þ xð1 xÞ
:
ð12Þ
We now make a change to it, and obtain a new family of mapping functions,
g i ðxÞ ¼ di þ
ðx d i Þ k
kþ1
A
ðx di Þ A þ xð1 xÞ
;
ð13Þ
where we require A > 0; k is a positive even integer. Note that, if one chooses A ¼ 1; k ¼ 2 in (13), we obtain the mapping function in [10] right away. That is, the mapping function by Henrick et al. belongs to the family of the new mapping functions (13). Theorem 1. The new mapping functions (13) have the following properties: (1) g 0i ðxÞ P 0; x 2 ½0; 1, (2) g i ð0Þ ¼ 0; g i ðdi Þ ¼ di ; g i ð1Þ ¼ 1, ðkÞ ðkþ1Þ (3) g 0i ðdi Þ ¼ ¼ g i ðdi Þ ¼ 0; g i ðdi Þ – 0.
456
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Proof. (1) The denominator in g i ðxÞ is greater than 0 for any x 2 ½0; 1 because A > 0; k is a positive even integer. Now we come to verify that g i ðxÞ is monotonically increasing in [0,1] with finite slopes. Firstly g 0i ðxÞ can be shown to have the following form: 2k
g 0i ðxÞ ¼
k
pnum ðxÞ ðx di Þ A2 þ ðx di Þ AqðxÞ ; ¼ 2 k pden ðxÞ ððx d Þ A þ xð1 xÞÞ
ð14Þ
i
where
qðxÞ ¼ ð1 kÞx2 þ ðk 2di Þx þ di : Note the qðxÞ is a quadratic polynomial with a negative coefficient, ð1 kÞ, for x2 . Use the fact that the minimum value of a quadratic polynomial function, which is defined in x 2 ½a; b, and has a negative leading coefficient, can only be found at x ¼ a or b. We conclude that
min qðxÞ ¼ min ðqð0Þ; qð1ÞÞ ¼ min ðdi ; 1 di Þ > 0:
x2½0;1
Immediatly we have pnum ðxÞ ¼ 0, when x ¼ di , and > 0, otherwise, because k is even. Therefore, we conclude that g 0i > 0 except g 0i ¼ 0 for x ¼ di . Thus g i ðxÞ is monotonically increasing in ½0; 1. (2) g i ð0Þ ¼ 0; g i ðdi Þ ¼ di and g i ð1Þ ¼ 1 are easily verified. k kþ1 ðkÞ ðkþ1Þ (3) Note ðx di Þ jpnum , but ðx di Þ -pnum , because qðxÞ > 0 in ½0; 1, thus g 0i ðdi Þ ¼ ¼ g i ðdi Þ ¼ 0; g i ðdi Þ – 0 are verified. h 3.2. The rate of convergence We denote the mapped WENO method with the new mapping function (13) as WENO-IM (k, A), so the WENO-HAP method in [10] is WENO-IM (2, 1). In the smooth regions, away from critical points, the WENO-JS schemes achieve the optimal order of accuracy, but near critical points, the schemes drop the order of accuracy, detailed analysis can be found in [10]. By using the mapping function, the mapped WENO schemes can recover optimal convergence rates. Let ncp denote the order of the critical point, that is, f 0 ¼ ¼ f ðncp Þ ¼ 0; f ðncp þ1Þ – 0. Now we will show that, the WENO-IM (k,A) scheme can recover the optimal convergence rates by choosing the parameter k for different values of ncp . Before giving Theorem 2, two lemmas are necessary as follows. Lemma 1. The sufficient condition for the ð2r 1Þ-order WENO scheme to achieve the optimal order of accuracy is
xi di ¼ OððMxÞr Þ; r ¼ 2; 3; . . . ; 6; i ¼ 0; 1; . . . ; r 1:
ð15Þ
Proof. Evaluation of Eq. (1) at x ¼ xj , we have
F^jþ1=2 F^j1=2 hjþ1=2 hj1=2 duj @f ¼ ¼ ; @x x¼xj dt Mx Mx
ð16Þ
where hðxÞ is defined according to
f ðxÞ ¼
1 Mx
Z
xþMx=2
hðnÞdn;
ð17Þ
xMx=2
and F^j1=2 is the numerical flux of ð2r 1Þ-order WENO scheme which is a convex combination of r-order fluxes associated with r substencils, that is,
F^j1=2 ¼
r1 X
xi F^ij1=2 :
ð18Þ
i¼0
Adding and subtracting
F^j1=2 ¼
Pr1 ^i i¼0 di F j1=2 from Eq. (18), we have
r1 r1 X X di F^ij1=2 þ ðxi di ÞF^ij1=2 : i¼0
i¼0
The first term satisfies r1 X di F^ij1=2 ¼ hj1=2 þ OðMx2r1 Þ; i¼0
ð19Þ
457
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
so the second term must be at least an OðMx2r Þ quantity in order for the derivative to be approximated at ð2r 1Þ order. Expanding the second term gives r1 r1 X X ðxi di ÞF^ij1=2 ¼ ðxi di Þðhj1=2 þ Ai Mxr þ OðMxrþ1 ÞÞ i¼0
i¼0
¼ hj1=2
r1 r1 r1 X X X ðxi di Þ þ Mxr Ai ðxi di Þ þ ðxi di ÞOðMxrþ1 Þ: i¼0
Because
Pr1
i¼0 ð
i¼0
ð20Þ
i¼0
xi di Þ ¼ 0, it is sufficient to require
xi di ¼ OððMxÞr Þ; i ¼ 0; 1; . . . ; r 1
ð21Þ
for the overall scheme to achieve the optimal order of accuracy. h Remark 1. The above proof is a natural extension to the one in [10]. The following lemma can be found in the statement of page 565 in [10]. We will now state it explicitly. Lemma 2. For different values of ncp , the weights xi in the ð2r 1Þ-order WENO-JS scheme satisfy
xi di ¼ OððMxÞr1ncp Þ; r ¼ 2; 3; . . . ; 6; ncp ¼ 0; 1; . . . ; r 1;
ð22Þ
then the rate of convergence is
rc ¼
2r 1;
if
2r 2 ncp ; if
ncp ¼ 0;
ð23Þ
ncp ¼ 1; 2; . . . ; r 1:
We now give Theorem 2 for the different order WENO-IM (k, A) method, the proof is almost identical with the one in [10]. Theorem 2. When ncp ¼ 1, the WENO-IM (k, A) schemes can achieve the optimal ð2r 1Þ-order of accuracy if the new mapping function is applied to the original weights in the WENO-JS scheme with k P 2 (except for the case r ¼ 2). Remark 2. Use Lemma 1 and 2, we can obtain the detailed convergence rates for the WENO-IM (k, A) scheme with r ¼ 2; . . . ; 6 for different values of ncp . Table 1 is a comparison between WENO-JS, WENO-HAP and WENO-IM (k, A). 3.3. The choices of parameters k and A The new mapping function in Eq. (13) and the effects of the parameters k and A are shown in Figs. 1 and 2. Decreasing A or increasing k means the mapped weights are more close to the ideal weights, thus the mapped WENO scheme approaches the Table 1 Rates of convergence, r c , for WENO-JS, WENO-HAP and WENO-IM (k, A) schemes. ncp
r c -WENO-JS
r c -WENO-HAP
r c -WENO-IM (k, A)
r=2
0 1
3 1
3 1
3ðk P 2Þ 1
r=3
0 1 2
5 3 2
5 5 2
5ðk P 2Þ 5ðk P 2Þ 2
r=4
0 1 2 3
7 5 4 3
7 7 6 3
7ðk P 2Þ 7ðk P 2Þ 7ðk P 4Þ 3
r=5
0 1 2 3 4
9 7 6 5 4
9 9 9 7 4
9ðk P 2Þ 9ðk P 2Þ 9ðk P 2Þ 9ðk P 4Þ 4
r=6
0 1 2 3 4 5
11 9 8 7 6 5
11 11 11 11 8 5
11ðk P 2Þ 11ðk P 2Þ 11ðk P 2Þ 11ðk P 2Þ 11ðk P 6Þ 5
458
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Fig. 1. New mapping functions for WENO-IM (2, A) with d2 ¼ 0:3, and A ¼ 0:1; 1 and 10 respectively.
Fig. 2. New mapping functions for WENO-IM (k, 1) with d2 ¼ 0:3, and k ¼ 2; 4 and 6 respectively.
upstream central scheme in a larger region. That is, in smooth regions, the mapped WENO scheme can get better solution by decreasing A or increasing k. On the other hand, the mapping function may amplify the effect from the non-smooth stencils and thus cause a potential loss of accuracy near discontinuities [4]. So we have to find a balance between them. The balance means that the mapped scheme can obtain more accurate solutions while still retaining the non-oscillatory property of WENO method. 3.3.1. The choice of parameter k In Table 1, we prove that, in theory, WENO-IM (k, A) can achieve more accurate numerical solutions when a larger k is chosen for higher order WENO method. However, extensive numerical tests show that, WENO-IM (k > 2, A) methods do not usually provide better solutions than WENO-IM (k ¼ 2, A) no matter which order of WENO method is used. This is an interesting observation. We will now give an explanation. From (14), because k is positive and even, qð0Þ ¼ di , a simple calculation demonstrates that
g 0i ð0Þ ¼ 1 þ
1 k1
Adi
:
ð24Þ
When increasing k; g 0i ð0Þ becomes larger, that means the mapping function amplifies the effect from the non-smooth stencils much faster. This effect will be much larger while the mapping function is applied to higher order WENO schemes, because
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
459
the smallest optimal weight which they possess is much smaller. For example, when the fifth-order WENO-IM (k, A) method is used, the new weight for the non-smooth stencil can be ð1 þ 10k1 =AÞ times as large as the one by the fifth-order WENO-JS. Furthermore, when the seventh-order WENO-IM (k, A) method is used, the new weight for the non-smooth stencil can be ð1 þ 35k1 =AÞ times as large as the one by the seventh-order WENO-JS. This is probably why the use of larger values in k dose not necessarily generate better numerical solutions. Therefore we only consider the case when k ¼ 2 in the rest of the paper. 3.3.2. The choice of parameter A While the new mapping function is applied to the seventh-order WENO schemes, WENO-IM (2, 1) which is the method in Henrick et al. [10], seems to be a good choice. That is, WENO-IM (2, A > 1) does not provide more accurate solutions than WENO-IM (2, 1) as A is increasing, and WENO-IM (2, A < 1) always generates the oscillatory by decreasing A, because the effects from the non-smooth stencils become much larger than the one in [10]. Once again, it verifies the additional requirement g 0i ð0Þ ¼ 0 in [4] is necessary, especially for higher order WENO schemes. While the new mapping function is applied to the fifth-order WENO schemes, extensive numerical experiments demonstrate that WENO-IM (2, A > 1) always provide less accurate solutions than WENO-IM (2, 1) as A is increasing, and WENO-IM (2, A < 1) always provide more accuracy solutions than WENO-IM (2, A ¼ 1) as A is decreasing. However while A is too small, oscillatory starts to be generated. Extensive numerical tests show that, no matter in short time or long time simulation, A ¼ 0:1 seems to be a very good choice. In other words, it provides an improved mapped WENO method with less dissipation
Fig. 3. New mapping functions for the fifth order WENO-IM (2, 0.1).
Fig. 4. Comparison for the mapping functions for WENO-IM (2, 0.1) and WENO-IM (2, 1) for d0 ¼ 0:1.
460
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
and higher resolution. In our numerical tests, we only consider the fifth-order WENO-IM (2, 0.1) method with the mapping function: 3
g i ð xÞ ¼ d i þ
0:1ðx di Þ 2
0:1ðx di Þ þ xð1 xÞ
ð25Þ
:
Fig. 3 gives the graph for the mapping functions for WENO-IM (2, 0.1). Fig. 4 gives the comparison for the mapping functions of WENO-IM (2, 0.1) and WENO-IM (2, 1) for d0 ¼ 0:1. Clearly WENO-IM (2, 0.1) tends to use the optimal weights in a larger range. This may be the reason why it can provide better numerical solutions when the mapping function is applied to the fifth order WENO scheme. However, when the mapping function is applied to higher order (order P 7) WENO methods, g 0i ð0Þ becomes very large. This starts to contaminate the error for the numerical solution. The mapping function (25) is unstable to provide the better solution than WENO-IM (2, 1), let alone WENO-PM6 [4]. 4. Numerical tests In this section, we apply WENO-JS ( ¼ 106 ), WENO-Z ( ¼ 1040 ), WENO-PM6 ( ¼ 1060 ) and WENO-IM (k, A) ( ¼ 1040 ) to solve the conservation laws in one dimension, where WENO-Z is used to refer to the methods in [2,3]. The Lax-Friedrichs flux splitting is used in the numerical tests. The errors in the L1 norm is computed by comparing the numerical solution u with the exact solution uexact , according to
kerrork1 ¼
N 1 X jui ðuexact Þi j; N þ 1 i¼0
where N is the number of intervals in space. We have mentioned in previous section, that for the seventh-order WENO-IM (k, A) methods, WENO-IM (2, 1) seems to be the better choice. However, it still can not compete with WENO-PM6 as shown in [4]. So we only consider the fifth-order WENO-IM (k, A) methods in this section. We now focus on the performance of fifth-order WENO-IM (2, 0.1) in this section. 4.1. Scalar conservation laws in one dimension Example 1. We solve the linear advection equation:
ut þ ux ¼ 0;
0 6 x 6 2;
ð26Þ
uðx; 0Þ ¼ sinðpxÞ;
together with the periodic boundary conditions. Table 2 shows the convergence rates and errors for numerical solutions of Example 1 in L1 norm at output time t ¼ 1. The first row is the error of L1 norm, and the second row is the convergence rates. The CFL number is chosen to be ðMxÞ2=3 in order that the error for the overall scheme is a measure of the spatial convergence only. We observe that, all the schemes reach the desired fifth order of accuracy. Example 2. We solve the linear advection equation:
ut þ ux ¼ 0;
0 6 x 6 2;
ð27Þ
4
uðx; 0Þ ¼ sin ðpxÞ; together with the periodic boundary conditions.
Table 2 Convergence study for various schemes as applied to Example 1 at output time t ¼ 1, all schemes are used with CFL¼ ðMxÞ2=3 . N
40
80
160
320
WENO-JS
2.2864e-005
WENO-Z
3.8994e-006
WENO-IM (2, 1)
3.8991e-006
WENO-PM6
3.8916e-006
WENO-IM (2, 0.1)
3.8924e-006
7.1493e-007 4.9991 1.2328e-007 4.9832 1.2328e-007 4.9831 1.2326e-007 4.9806 1.2326e-007 4.9809
2.2438e-008 4.9938 3.8770e-009 4.9909 3.8770e-009 4.9909 3.8769e-009 4.9907 3.8769e-009 4.9907
7.0199e-010 4.9983 1.2154e-010 4.9954 1.2154e-010 4.9954 1.2154e-010 4.9954 1.2154e-010 4.9954
461
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Table 3 shows the convergence rates and errors for numerical solutions of Example 2 in L1 norm at output time t ¼ 1. The first row is the error of L1 norm, and the second row is the convergence rates. The CFL number is chosen to be ðMxÞ2=3 in order that the error for the overall scheme is a measure of the spatial convergence only. In the aspect of convergence rates, the WENO-JS scheme achieves the desired fifth order of accuracy when N P 160, and all the other schemes do not reach it. In the aspects of accuracy, The WENO-IM (2,0.1) provides the most accurate numerical solutions and is more accurate than other schemes almost by a factor of two at least when N P 80. Example 3. We solve the linear advection equation:
ut þ ux ¼ 0;
1 6 x 6 1;
ð28Þ
uðx; 0Þ ¼ u0 ðxÞ; where
81 ðGðx; b; z dÞ þ Gðx; b; z þ dÞ þ 4Gðx; b; zÞÞ; > 6 > > > > 1; > < u0 ¼ 1 j10ðx 0:1Þj; > > 1 > > > 6 ðFðx; a; a dÞ þ Fðx; a; a þ dÞ þ 4Fðx; a; aÞÞ; > : 0;
0:8 6 x 6 0:6; 0:4 6 x 6 0:2; ð29Þ
0 6 x 6 0:2; 0:4 6 0:6; otherwise;
and 2
Gðx; b; zÞ ¼ ebðxzÞ ;
Fðx; a; aÞ ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi maxð1 a2 ðx aÞ2 ; 0Þ:
ð30Þ
together with the periodic boundary condition. The constants are taken as a ¼ 0:5; z ¼ 0:7; d ¼ 0:0005; a ¼ 10 and log2 b ¼ 36d 2 . The solution contains a smooth but narrow combination of Gaussians, a square wave, a sharp triangle wave and
a half ellipse. Table 4 and 5 show the convergence rates and errors for numerical solutions of Example 3 in L1 norm at output time t ¼ 20 and t ¼ 2000 respectively. The first row is the error of L1 norm, and the second row is the convergence rates. The CFL number is chosen to be 0:1. When t ¼ 20, all the schemes converge with first order accuracy, and WENO-IM (2, 0.1) gives more accurate numerical solutions than other schemes. When t ¼ 2000, WENO-JS, WENO-Z and WENO-IM (2, 1) reduce their convergence rates seriously, especially, WENO-Z gives the negative convergence rate on the mesh of N ¼ 400, however, WENO-PM6 and WENO-IM (2, 0.1) still converge with first order accuracy, and on the mesh of N ¼ 400; 800 cells, WENOPM6 gives the best numerical solutions, WENO-IM (2, 0.1) gives slightly less accurate solutions than WENO-PM6, but more accurate solutions than WENO-JS, WENO-Z and WENO-IM (2, 1). Figs. 5–8 show the performance of various schemes at output time t ¼ 2000 on the mesh of N ¼ 400; 800 cells. On the mesh of N ¼ 400 cells, WENO-PM6 gives the highest resolution of Gaussians, however generates a small oscillation near x ¼ 0:5, and WENO-IM (2, 0.1) gives the highest resolution of square wave. On the mesh of N ¼ 800 cells, WENO-PM6 still gives the highest resolution of Gaussians, however generates a small oscillation near x ¼ 0:7, at this stage, WENO-PM6 and WENO-IM (2, 0.1) give the same resolution of square wave, and do not smear this feature. Table 6 shows the CPU time (in seconds) of various schemes at output time t ¼ 10 on the mesh of N ¼ 200; 400; 800 and 1600 cells respectively. The Matlab Language is used to devise the program for all the schemes, all the results presented here are computed on a PC with an Intel (R) Core (TM) i7–3770 K desktop processor and with 16.0 GB of memory. In order to ensure fairness in the comparison of CPU time, all the schemes share the same subroutine calls and are compiled with the same compilation options including optimization ones. The only differences between the implementation of all the schemes are on the subroutines for computing the different related WENO weights. We observe that WENO-Z is the most efficient scheme. However, it provides less accurate numerical solution than WENO-IM (2, 0.1). WENO-IM (2, 0.1) is slightly
Table 3 Convergence study for various schemes as applied to Example 2 at output time t ¼ 1, all schemes are used with CFL¼ ðMxÞ2=3 . N
40
80
160
320
WENO-JS
3.6393e-003
WENO-Z
1.4177e-003
WENO-IM (2, 1)
1.6499e-003
WENO-PM6
2.2364e-003
WENO-IM (2, 0.1)
1.3977e-003
4.7793e-004 2.9288 1.0794e-004 3.7153 1.4942e-004 3.4649 1.5091e-004 3.8894 6.7117e-005 4.3802
1.4829e-005 5.0103 5.5845e-006 4.2726 7.7330e-006 4.2722 8.3183e-006 4.1813 2.8832e-006 4.5409
3.8676e-007 5.2608 2.6142e-007 4.4170 3.5976e-007 4.4259 3.8574e-007 4.4306 1.0827e-007 4.7349
462
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Table 4 Convergence study for various schemes as applied to Example 3 at output time t ¼ 20, all schemes are used with CFL¼ 0:1. N
100
200
400
800
WENO-JS
1.2152e-001
WENO-Z
9.4926e-002
WENO-IM (2, 1)
9.6111e-002
WENO-PM6
9.7839e-002
WENO-IM (2, 0.1)
9.4494e-002
5.7747e-002 1.0733 3.8757e-002 1.2923 4.0298e-002 1.2540 4.0478e-002 1.2733 3.7664e-002 1.3270
2.5255e-002 1.1932 1.7682e-002 1.1322 1.8845e-002 1.0965 1.7411e-002 1.2171 1.7079e-002 1.1410
1.4352e-002 0.8153 8.3518e-003 1.0821 9.5250e-003 0.9844 8.1661e-003 1.0923 8.0227e-003 1.0901
Table 5 Convergence study for various schemes as applied to Example 3 at output time t ¼ 2000, all schemes are used with CFL¼ 0:1. N
100
200
400
800
WENO-JS
3.1104e-001
WENO-Z
2.5132e-001
WENO-IM (2, 1)
2.5496e-001
WENO-PM6
2.4304e-001
WENO-IM (2, 0.1)
2.3976e-001
2.9819e-001 0.0608 1.5717e-001 0.6772 1.9836e-001 0.3621 1.1025e-001 1.1404 1.1019e-001 1.1215
2.8040e-001 0.0888 1.9006e-001 0.2742 1.6024e-001 0.3079 5.2225e-002 1.0780 5.6666e-002 0.9595
2.4749e-001 0.1801 1.3242e-001 0.5214 1.5803e-001 0.0201 2.2554e-002 1.2113 2.6075e-002 1.1198
Fig. 5. Performance of the fifth-order WENO-JS, WENO-Z and WENO-IM (2, 0.1) for Example 3 when t ¼ 2000. The mesh of N ¼ 400 cells and CFL¼ 0:1 are used for all schemes.
slower than WENO-IM (2, 1) because of the additional multiplication with A, but is much faster than WENO-PM6 since WENO-PM6 is not very suitable for the vector operation due to the implementation of logic operations. 4.2. Euler system in one dimension We solve one-dimensional Euler equations of gas dynamics, but with different initial and boundary conditions. For the one-dimensional Euler equations, we have the conservation of mass, momentum and energy:
0
1
0
1
q qu B C B C @ qu A þ @ qu2 þ p A ¼ 0; E
t
uðE þ pÞ
x
ð31Þ
463
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Fig. 6. Performance of the fifth-order WENO-IM (2, 1), WENO-PM6 and WENO-IM (2, 0.1) for Example 3 when t ¼ 2000. The mesh of N ¼ 400 cells and CFL¼ 0:1 are used for all schemes.
Fig. 7. Performance of the fifth-order WENO-JS, WENO-Z and WENO-IM (2, 0.1) for Example 3 when t ¼ 2000. The mesh of N ¼ 800 cells and CFL¼ 0:1 are used for all schemes.
where q; u; E; p are the density, velocity, total energy and pressure, respectively. The system of Euler equations is closed by the equation of state for an ideal polytropic gas
E¼
p
c1
1 þ qu2 ; 2
ð32Þ
where c is the ratio of specific heat. For the examples in this paper, we use c ¼ 1:4. The numerical flux F^jþ1 at xjþ1 in (2) is calculated as follow [14]: 2
2
(1) Compute an average state ujþ1 using a Roe average; 2 (2) Compute the right eigenvectors, the left eigenvectors and the eigenvalues of the Jacobian f 0 ðujþ1 Þ, and denote them by 2
R ¼ Rðujþ1 Þ; 2
R
1
1
¼ R ðujþ1 Þ; 2
K ¼ Kðujþ12 Þ;
(3) Transform the flux and the solution which are in the potential stencil of the WENO reconstructions for obtaining the flux F^jþ1 , to the local characteristic fields by 2
vj ¼ R
1
uj ;
g j ¼ R1 f ðuj Þ;
j in a neighborhood of i;
(4) Use the Lax-Friedrichs flux splitting and the WENO reconstruction procedure for each component of the characteristic variables, to obtain the corresponding component of the flux b ; g jþ1 2
464
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Fig. 8. Performance of the fifth-order WENO-IM (2, 1), WENO-PM6 and WENO-IM (2, 0.1) for Example 3 when t ¼ 2000. The mesh of N ¼ 800 cells and CFL¼ 0:1 are used for all schemes.
Table 6 CPU time (in seconds) study for various schemes as applied to Example 3 at output time t ¼ 10, all schemes are used with CFL¼ 0:1. N
200
400
800
1600
WENO-JS WENO-Z WENO-IM (2, 1) WENO-PM6 WENO-IM (2, 0.1)
10.6 10.1 11.6 14.3 12.3
29.4 27.7 32.3 41.1 33.8
95.9 88.2 105.0 136.7 110.3
379.5 350.4 412.5 530.1 429.4
(5) Transform back into physical space by using
b F jþ1 ¼ R b g jþ1 ; 2
2
(6) Form the flux by taking
b F jþ1 ¼ b F þjþ1 þ b F jþ1 2
2
2
Example 4 (Shock tube problem with a sonic point). We solve the Euler Equations (31) on ½5; 5 with following initial condition [20]
ðq; u; pÞ ¼
ð1:0; 0:75; 1:0Þ; if x < 2; ð0:125; 0:0; 0:1Þ; if x P 2:
ð33Þ
This is a modified version of Sod’s problem [18], its solution has a right shock wave, a right traveling contact wave and a left sonic rarefaction wave. We compute the solution at the output time t ¼ 2 on the mesh of N ¼ 100; 200; 400 and 400 cells. The CFL number is chosen to be 0:4. Table 7 shows the errors of density for various WENO schemes in L1 norm, and WENO-IM (2, 0.1) provides the most accurate numerical solutions on different meshes. Fig. 9 is the plot for the density on the mesh of N ¼ 100. We observe that, all the schemes perform well at the sonic point.
Table 7 Performance of different fifth-order WENO schemes on shock tube problem with a sonic point, CFL¼ 0:4. N
100
200
400
800
WENO-JS WENO-Z WENO-IM (2, 1) WENO-PM6 WENO-IM (2, 0.1)
8.6385e-003 6.9880e-003 7.6275e-003 7.5799e-003 6.8783e-003
4.5876e-003 3.6937e-003 4.0102e-003 3.9629e-003 3.6330e-003
2.4777e-003 2.0126e-003 2.1624e-003 2.1260e-003 1.9772e-003
1.3338e-003 1.0828e-003 1.1562e-003 1.1345e-003 1.0636e-003
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
465
Fig. 9. Density plot for the shock tube problem with a sonic point when t ¼ 2. The mesh of N ¼ 100 cells and CFL¼ 0:4 are used for all schemes.
Fig. 10. Performance of the fifth-order WENO-IM (2, 0.1) for modified shock/turbulence interaction when t ¼ 5. The mesh of N ¼ 1500 cells and CFL¼ 0:4 are used.
Fig. 11. Comparison of the fifth-order WENO-JS and WENO-IM (2, 0.1) for modified shock/turbulence interaction when t ¼ 5. The mesh of N ¼ 1500 cells and CFL¼ 0:4 are used for all schemes.
466
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Fig. 12. Comparison of the fifth-order WENO-Z and WENO-IM (2, 0.1) for modified shock/turbulence interaction when t ¼ 5. The mesh of N ¼ 1500 cells and CFL¼ 0:4 are used for all schemes.
Fig. 13. Comparison of the fifth-order WENO-IM (2, 1) and WENO-IM (2, 0.1) for modified shock/turbulence interaction when t ¼ 5. The mesh of N ¼ 1500 cells and CFL¼ 0:4 are used for all schemes.
Fig. 14. Comparison of the fifth-order WENO-PM6 and WENO-IM (2, 0.1) for modified shock/turbulence interaction when t ¼ 5. The mesh of N ¼ 1500 cells and CFL¼ 0:4 are used for all schemes.
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
467
Example 5 (Modified shock/turbulence interaction). We solve the Euler equations (31) on ½5; 5 with following initial condition [19,21]
ðq; u; pÞ ¼
ð1:515695; 0:523346; 1:80500Þ; if x < 4:5; ð1 þ 0:1 sinð20pxÞ; 0:0; 1Þ;
if x P 4:5:
ð34Þ
The initial condition is set by a Mach 1:1 shock interacting with a high-frequency density perturbation. As time evolves, the shock moves into this perturbation, which spreads upstream. We compute the solution at the output time t ¼ 5 on the mesh of N ¼ 1500 cells. The CFL number is chosen to be 0:4. The density computed by WENO-IM (2, 0.1) is given in Fig. 10, the reference solution is an approximation to the exact solution using WENO-JS with N ¼ 10000. Fig. 11–14 show the comparisons of WENO-IM (2, 0.1) with other schemes on ½2; 1, we observe that WENO-IM (2, 0.1) gives the best resolution. 4.3. Summary of performance of fifth order WENO-IM (2, 0.1) We summarize the performance of our method as follow. (1) All the improved WENO methods perform better than the original WENO method. (2) Our method performs better than other improved WENO methods on smooth problems with large gradients. As presented in Example 2, the solution of our method is more accurate than other methods, in fact by a factor of two at least. (3) The most important advantage of our method presents on problems with discontinuities. The improvement is obvious if a long simulation is desired or if the solutions have many high-frequency waves or dense discontinuities. Example 3 demonstrates the improvement of our method if a long simulation is desired (see Table 5). And Example 5, the modified shock/turbulence interaction, demonstrates the improvement of our method by the mass high-frequency density perturbations. Note that the high-frequency perturbations are actually numerical discontinuities if a coarse mesh is used. Our method performs much better than WENO-Z and WENO-IM (2, 1) in these two examples. (4) Our method gives comparable solution with WENO-PM6 in Example 3, and gives more accurate solution than WENOPM6 in Example 5. In terms of CPU time, our method runs much faster then WENO-PM6. (5) The mapping function of our method is smoother than the one in WENO-PM6, because our mapping functions use rational functions but the mapping functions in WENO-PM6 use piecewise polynomial. (6) If the major parts of the numerical solutions are smooth while with a few discontinuities, and if a short time simulation is desired; e.g. Table 4 for Example 3 and Table 7 for Example 4, the Sod-type problem, the improvement of our method may be only minor comparing to other modified WENO methods. 5. Conclusions In this paper, a new mapping function is proposed, then an improved version of the mapped WENO method in Henrick et al. [10] is developed, we call it WENO-IM (k, A). The improved mapped WENO method can achieve the optimal order of accuracy near critical points by proper choice of parameter k. When the new mapping function is applied to the seventhorder WENO method, WENO-IM (2, 1), which is the same method proposed by Henrick et al. [10], seems to be a good choice. However, when the new mapping function is applied to the fifth-order method, WENO-IM (2, 1) is not the best choice any more. Numerical tests show that, no matter in short time or long time simulation, WENO-IM (2, 0.1) seems to be a better choice for the fifth-order WENO method. It provides better numerical solutions with less dissipation and higher resolution than WENO-JS, WENO-Z and WENO-IM (2, 1), special in long time simulation, and provides comparable numerical solutions with WENO-PM6, but it is more efficient than WENO-PM6 in terms of CPU time. Acknowledgments The authors would like to thank the anonymous referees for valuable suggestions and comments. The work of H. Feng was supported by National Natural Science Foundation of China (Nos. 91130022, 10971159 and 11161130003) and NCET of China. References [1] D.S. Balsara, C.W. Shu, Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy, J. Comput. Phys. 160 (2000) 405–452. [2] R. Borges, M. Carmona, B. Costa, W.S. Don, An improved weighted essentially non-oscillation scheme for hypebolic conservation laws, J. Comput. Phys. 227 (2008) 3191–3211. [3] M. Castro, B. Costa, W.S. Don, High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws, J. Comput. Phys. 230 (2011) 1766–1792. [4] H. Feng, F.X. Hu, R. Wang, A new mapped weighted essentially non-oscillatory scheme, J. Sci. Comput. 51 (2012) 449–473. [5] S. Gottlieb, C.W. Shu, E. Tadmor, Strong stability-preserving time discretization methods, SIAM Rev. 43 (2001) 89–112. [6] A. Harten, ENO schemes with subcell resolution, J. Comput. Phys. 83 (1987) 148–184. [7] A. Harten, B. Engquist, S. Osher, S. Chakravarthy, Uniformly high order essentially non-oscillatory schemes III, J. Comput. Phys. 71 (1987) 231–303.
468
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
[8] A. Harten, S. Osher, Uniformly high order essentially non-oscillatory schemes I, SIAM J. Numer. Anal. 24 (1987) 279–309. [9] A. Harten, S. Osher, B. Engquist, S. Chakravarthy, Some results on uniformly high order accurate essentially non-oscillatory schemes, Appl. Numer. Math. 2 (1986) 347–377. [10] A.K. Henrick, T.D. Aslam, J.M. Powers, Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points, J. Comput. Phys. 207 (2005) 542–567. [11] G.S. Jiang, C.W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. [12] X.D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994) 200–212. [13] C.W. Shu, Numerical experiments on the accuracy of ENO and modified ENO schemes, J. Sci. Comput. 5 (1990) 127–149. [14] C.W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, in: Advanced numerical approximation of nonlinear hyperbolic equations, Lecture Notes in Mathematics, vol. 1697, Springer-Verlag, Berlin, 1998, pp. 325–432. [15] C.W. Shu, High order weighted essentially non-oscillatory schemes for convection dominated problems, SIAM Rev. 51 (2009) 82–126. [16] C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77 (1988) 439–471. [17] C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes II, J. Comput. Phys. 83 (1989) 32–78. [18] G.A. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws, J. Comput. Phys. 107 (1978) 1–31. [19] V.A. Titarev, E.F. Toro, Weno schemes based on upwind and centred tvd fluxes, Comput. Fluids 34 (2005) 705–720. [20] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, third ed., Springer-Verlag, Berlin Heidelberg, 2009. [21] E.F. Toro, V.A. Titarev, Tvd fluxes for the high-order ader schemes, J. Sci. Comput. 24 (2005) 285–309.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
An improved mapped weighted essentially non-oscillatory scheme Hui Feng a, Cong Huang a, Rong Wang b,⇑ a b
School of Mathematics and Statistics, Wuhan University, Wuhan, Hubei, PR China Department of General Education, South University of Science and Technology of China, Shenzhen, Guangdong, PR China
a r t i c l e
i n f o
Keywords: WENO methods Hyperbolic conservation laws Mapping function
a b s t r a c t In this paper we develop an improved version of the mapped weighted essentially non-oscillatory (WENO) method in Henrick et al. (2005) [10] for hyperbolic partial differential equations. By rewriting and making a change to the original mapping function, a new type of mapping functions is obtained. There are two parameters, namely A and k, in the new mapping functions (see Eq. (13)). By choosing k ¼ 2 and A ¼ 1, it leads to the mapping function in Henrick et al., i.e.; the mapped WENO method by Henrick et al. actually belongs to the family of our improved mapped WENO schemes. Furthermore, we show that, when the new mapping function is applied to any ð2r 1Þth order WENO scheme for proper choice of k, it can achieve the optimal order of accuracy near critical points. Note that, if only one mapping is used, the mapped WENO method by Henrick et al., whose order is higher than five, can not achieve the optimal order of accuracy in some cases. Through extensive numerical tests, we draw the conclusion that, the mapping function proposed by Henrick et al. is not the best choice for the parameter A. A new mapping function is then selected and provides an improved mapped WENO method with less dissipation and higher resolution. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction The essential non-oscillatory (ENO) method was studied by Harten et al. in [6,8,9,7]. It owes its success to the ENO property which is considered to be very important in the numerical simulation for hyperbolic conservation laws because high order numerical methods often produce spurious oscillations, especially near shocks or other discontinuities. The finite-volume ENO spatial discretization was studied by Harten et al. in [7], later the finite-difference ENO method was developed by Shu et al. in [16,17]. Later weighted ENO (WENO) methods have been introduced in [12,11] to address potential numerical instabilities due to the process of choosing ENO stencils [13]. Unlike ENO methods, a convex combination of all the ENO candidate stencils are used in WENO methods such that the numerical flux is approximated with the higher order of accuracy in smooth regions, while WENO methods still retain the ENO property in regions near discontinuities. In this paper we will use WENO-JS to refer the WENO scheme in [11,1]. See [14,15] for an overview and other references. Recently Henrick et al. [10] noted that the fifth-order WENO method proposed by Jiang and Shu [11] is only third-order accurate near critical points of the smooth regions in general. By using a simple mapping function to the original weights in ⇑ Corresponding author. E-mail addresses: [email protected] (H. Feng), [email protected] (C. Huang), [email protected] (R. Wang). http://dx.doi.org/10.1016/j.amc.2014.01.061 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
454
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
[11], Henrick et al. developed a mapped WENO method to achieve the optimal order of accuracy; i.e., the fifth order near critical points. However, when this mapping function is applied to the WENO scheme of the order P 7, it can not achieve the optimal order of accuracy near critical points. In this paper we use WENO-HAP to refer the family of mapped WENO schemes in [11]. Lately Feng et al. [4] studied the mapped WENO method in [10] and found that, when it is used for solving the problems with discontinuities, the mapping function in [10] may amplify the effect from the non-smooth stencils and thus cause a potential loss of accuracy near discontinuities. This effect may be difficult to be observed for the fifth-order WENO method unless a long time simulation is desired. In order to overcome this problem, a new piecewise polynomial mapping function was proposed in [4] by adding two requirements, that is, g 0 ð0Þ ¼ g 0 ð1Þ ¼ 0 to the original criteria in [10]. Clearly the requirement g 0 ð0Þ ¼ 0 tends to decrease the effect from the non-smooth stencils and g 0 ð1Þ ¼ 0 tends to increase the effect from the smoothest stencil. In this paper we use WENO-PM to refer the family of mapped WENO schemes in [4]. Through a further study of the mapped WENO method, we find that, the mapping function in [10] can be written in a simple and more meaningful form. Make a change to the new form, then a new family of mapping functions, which includes the one in WENO-HAP, is obtained. This kind of new mapping functions still satisfies the original criteria in [10]. Furthermore, when the new mapping function is applied to any ð2r 1Þth order WENO scheme for proper choice of k (see Theorem 2 for more details), it can achieve the optimal order of accuracy near critical points. Note that, if only one mapping is used, the mapped WENO method by Henrick et al., which order is higher than five, can not achieve the optimal order of accuracy in some cases. In order to overcome this problem, two or more mappings need to be used. However, in [4] Feng et al. have pointed out, the more mappings it uses, the less accurate numerical solution it provides due to potential numerical instability. Extensive numerical tests show that, the mapping function in WENO-HAP is not the best mapped WENO method in this new family. Then a new mapping function is selected. It provides an improved mapped WENO method with less dissipation and higher resolution. The remainder of this paper proceeds as follows. In Section 2, we give a brief description of WENO-JS [11] and WENO-HAP [10]. In Section 3, we give the improved version of the mapped WENO-HAP method. In Section 4, we compare the performance of different WENO methods. 2. The fifth-order WENO method 2.1. WENO spatial discretization Consider the one-dimensional scalar hyperbolic conservation law
ut ¼ fx ðuÞ;
a < x < b; t > 0:
ð1Þ
Assume a uniform spatial mesh; i.e., xj ¼ a þ jDx; j ¼ 0; 1; . . . ; N, where Dx ¼ ba . The spatial derivative fx ðuÞ at the point N x ¼ xj is approximated using a conservative finite difference scheme
LðuÞ ¼
1 ^ F jþ1 F^j1 ; 2 2 Dx
ð2Þ
where F^j1 is the numerical flux at xj1 . The numerical flux must be consistent with f; i.e., F^ðu; . . . ; uÞ ¼ f ðuÞ. The numerical 2 2 flux determines the particular numerical method and its properties. The fifth-order WENO scheme [11] has a numerical flux which is a convex combination of three third-order fluxes associated with three substencils; i.e.,
F^jþ1 ¼ x0 F^0jþ1 þ x1 F^1jþ1 þ x2 F^2jþ1 ; 2
2
2
ð3Þ
2
where
F^0jþ1 2
F^1jþ1 2
F^2jþ1 2
2f j2 7f j1 þ 11f j ; ¼ 16 fj1 þ 5f j þ 2f jþ1 ; ¼ 16 2f j þ 5f jþ1 fjþ2 : ¼
1 6
ð4Þ
In smooth regions, the WENO approximation is designed such that it approximates the solution with the high order of accuracy. On the other hand, near discontinuities, it keeps the essentially non-oscillatory property of the ENO scheme. The nonlinear weights are computed as
a xi ¼ P2 i
J¼0 aJ
;
ai ¼
di ð þ ISi Þ2
;
i ¼ 0; 1; 2;
ð5Þ
where is a small positive number that is introduced to prevent the denominator becoming zero, and the optimal weights, di ; i ¼ 0; 1; 2, are given by
d0 ¼ 1=10;
d1 ¼ 6=10;
and d2 ¼ 3=10:
ð6Þ
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
455
ISi are the smoothness indicators, which are defined as follows
2 1 2 13 fj2 2f j1 þ fj þ fj2 4f j1 þ 3f j ; 12 4 2 1 2 13 IS1 ¼ fj1 2f j þ fjþ1 þ fj1 fjþ1 ; 12 4 2 1 2 13 IS2 ¼ fj 2f jþ1 þ fjþ2 þ 3f j 4f jþ1 þ fjþ2 : 12 4 IS0 ¼
ð7Þ
Henrick et al. [10] then introduced a mapping function of the form
ðg HAP Þi ðxÞ ¼
x di þ ðdi Þ2 3di x þ x2
;
2
ðdi Þ þ ð1 2di Þx
ð8Þ
where di are the optimal weights given in (6), x 2 ½0; 1 and i ¼ 0; 1; 2. This function is monotonically increasing in ½0; 1 with 0 00 finite slopes, ðg HAP Þi ð0Þ ¼ 0; ðg HAP Þi ð1Þ ¼ 1, ðg HAP Þi ðdi Þ ¼ di , and ðg HAP Þi ðdi Þ ¼ ðg HAP Þi ðdi Þ ¼ 0. The mapped weights are given by:
aHAP xHAP ¼ P2 i ; aHAP ¼ ðg HAP Þi ðxi Þ; i ¼ 0; 1; 2; i i HAP J¼0 aJ
ð9Þ
where xi are computed by (5) and (7). In smooth regions, this mapped WENO scheme; i.e., the WENO-HAP scheme, gives the fifth order of accuracy even near critical points. We refer to [10] for more details. 2.2. Time discretization After the spatial derivative is discretized with one of the WENO schemes, an Ordinary Differential Equation (ODE) system is obtained as follow
du ¼ LðuÞ; dt
ð10Þ
where the spatial operator LðuÞ is represented in (2). In all the numerical examples in this paper, the ODE system (10) is solved using the explicit, third-order, TVD, Runge–Kutta method [16,5] as follows
uð1Þ ¼ un þ DtLðun Þ; 3 1 1 uð2Þ ¼ un þ uð1Þ þ DtLðuð1Þ Þ; 4 4 4 1 n 2 ð2Þ 2 nþ1 u ¼ u þ u þ DtLðuð2Þ Þ: 3 3 3
ð11Þ
3. The new mapped WENO methods 3.1. A new mapping function The mapping function (8) can be written in a simple and more meaningful form:
ðg HAP Þi ðxÞ ¼ di þ
ðx di Þ
3
2
ðx di Þ þ xð1 xÞ
:
ð12Þ
We now make a change to it, and obtain a new family of mapping functions,
g i ðxÞ ¼ di þ
ðx d i Þ k
kþ1
A
ðx di Þ A þ xð1 xÞ
;
ð13Þ
where we require A > 0; k is a positive even integer. Note that, if one chooses A ¼ 1; k ¼ 2 in (13), we obtain the mapping function in [10] right away. That is, the mapping function by Henrick et al. belongs to the family of the new mapping functions (13). Theorem 1. The new mapping functions (13) have the following properties: (1) g 0i ðxÞ P 0; x 2 ½0; 1, (2) g i ð0Þ ¼ 0; g i ðdi Þ ¼ di ; g i ð1Þ ¼ 1, ðkÞ ðkþ1Þ (3) g 0i ðdi Þ ¼ ¼ g i ðdi Þ ¼ 0; g i ðdi Þ – 0.
456
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Proof. (1) The denominator in g i ðxÞ is greater than 0 for any x 2 ½0; 1 because A > 0; k is a positive even integer. Now we come to verify that g i ðxÞ is monotonically increasing in [0,1] with finite slopes. Firstly g 0i ðxÞ can be shown to have the following form: 2k
g 0i ðxÞ ¼
k
pnum ðxÞ ðx di Þ A2 þ ðx di Þ AqðxÞ ; ¼ 2 k pden ðxÞ ððx d Þ A þ xð1 xÞÞ
ð14Þ
i
where
qðxÞ ¼ ð1 kÞx2 þ ðk 2di Þx þ di : Note the qðxÞ is a quadratic polynomial with a negative coefficient, ð1 kÞ, for x2 . Use the fact that the minimum value of a quadratic polynomial function, which is defined in x 2 ½a; b, and has a negative leading coefficient, can only be found at x ¼ a or b. We conclude that
min qðxÞ ¼ min ðqð0Þ; qð1ÞÞ ¼ min ðdi ; 1 di Þ > 0:
x2½0;1
Immediatly we have pnum ðxÞ ¼ 0, when x ¼ di , and > 0, otherwise, because k is even. Therefore, we conclude that g 0i > 0 except g 0i ¼ 0 for x ¼ di . Thus g i ðxÞ is monotonically increasing in ½0; 1. (2) g i ð0Þ ¼ 0; g i ðdi Þ ¼ di and g i ð1Þ ¼ 1 are easily verified. k kþ1 ðkÞ ðkþ1Þ (3) Note ðx di Þ jpnum , but ðx di Þ -pnum , because qðxÞ > 0 in ½0; 1, thus g 0i ðdi Þ ¼ ¼ g i ðdi Þ ¼ 0; g i ðdi Þ – 0 are verified. h 3.2. The rate of convergence We denote the mapped WENO method with the new mapping function (13) as WENO-IM (k, A), so the WENO-HAP method in [10] is WENO-IM (2, 1). In the smooth regions, away from critical points, the WENO-JS schemes achieve the optimal order of accuracy, but near critical points, the schemes drop the order of accuracy, detailed analysis can be found in [10]. By using the mapping function, the mapped WENO schemes can recover optimal convergence rates. Let ncp denote the order of the critical point, that is, f 0 ¼ ¼ f ðncp Þ ¼ 0; f ðncp þ1Þ – 0. Now we will show that, the WENO-IM (k,A) scheme can recover the optimal convergence rates by choosing the parameter k for different values of ncp . Before giving Theorem 2, two lemmas are necessary as follows. Lemma 1. The sufficient condition for the ð2r 1Þ-order WENO scheme to achieve the optimal order of accuracy is
xi di ¼ OððMxÞr Þ; r ¼ 2; 3; . . . ; 6; i ¼ 0; 1; . . . ; r 1:
ð15Þ
Proof. Evaluation of Eq. (1) at x ¼ xj , we have
F^jþ1=2 F^j1=2 hjþ1=2 hj1=2 duj @f ¼ ¼ ; @x x¼xj dt Mx Mx
ð16Þ
where hðxÞ is defined according to
f ðxÞ ¼
1 Mx
Z
xþMx=2
hðnÞdn;
ð17Þ
xMx=2
and F^j1=2 is the numerical flux of ð2r 1Þ-order WENO scheme which is a convex combination of r-order fluxes associated with r substencils, that is,
F^j1=2 ¼
r1 X
xi F^ij1=2 :
ð18Þ
i¼0
Adding and subtracting
F^j1=2 ¼
Pr1 ^i i¼0 di F j1=2 from Eq. (18), we have
r1 r1 X X di F^ij1=2 þ ðxi di ÞF^ij1=2 : i¼0
i¼0
The first term satisfies r1 X di F^ij1=2 ¼ hj1=2 þ OðMx2r1 Þ; i¼0
ð19Þ
457
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
so the second term must be at least an OðMx2r Þ quantity in order for the derivative to be approximated at ð2r 1Þ order. Expanding the second term gives r1 r1 X X ðxi di ÞF^ij1=2 ¼ ðxi di Þðhj1=2 þ Ai Mxr þ OðMxrþ1 ÞÞ i¼0
i¼0
¼ hj1=2
r1 r1 r1 X X X ðxi di Þ þ Mxr Ai ðxi di Þ þ ðxi di ÞOðMxrþ1 Þ: i¼0
Because
Pr1
i¼0 ð
i¼0
ð20Þ
i¼0
xi di Þ ¼ 0, it is sufficient to require
xi di ¼ OððMxÞr Þ; i ¼ 0; 1; . . . ; r 1
ð21Þ
for the overall scheme to achieve the optimal order of accuracy. h Remark 1. The above proof is a natural extension to the one in [10]. The following lemma can be found in the statement of page 565 in [10]. We will now state it explicitly. Lemma 2. For different values of ncp , the weights xi in the ð2r 1Þ-order WENO-JS scheme satisfy
xi di ¼ OððMxÞr1ncp Þ; r ¼ 2; 3; . . . ; 6; ncp ¼ 0; 1; . . . ; r 1;
ð22Þ
then the rate of convergence is
rc ¼
2r 1;
if
2r 2 ncp ; if
ncp ¼ 0;
ð23Þ
ncp ¼ 1; 2; . . . ; r 1:
We now give Theorem 2 for the different order WENO-IM (k, A) method, the proof is almost identical with the one in [10]. Theorem 2. When ncp ¼ 1, the WENO-IM (k, A) schemes can achieve the optimal ð2r 1Þ-order of accuracy if the new mapping function is applied to the original weights in the WENO-JS scheme with k P 2 (except for the case r ¼ 2). Remark 2. Use Lemma 1 and 2, we can obtain the detailed convergence rates for the WENO-IM (k, A) scheme with r ¼ 2; . . . ; 6 for different values of ncp . Table 1 is a comparison between WENO-JS, WENO-HAP and WENO-IM (k, A). 3.3. The choices of parameters k and A The new mapping function in Eq. (13) and the effects of the parameters k and A are shown in Figs. 1 and 2. Decreasing A or increasing k means the mapped weights are more close to the ideal weights, thus the mapped WENO scheme approaches the Table 1 Rates of convergence, r c , for WENO-JS, WENO-HAP and WENO-IM (k, A) schemes. ncp
r c -WENO-JS
r c -WENO-HAP
r c -WENO-IM (k, A)
r=2
0 1
3 1
3 1
3ðk P 2Þ 1
r=3
0 1 2
5 3 2
5 5 2
5ðk P 2Þ 5ðk P 2Þ 2
r=4
0 1 2 3
7 5 4 3
7 7 6 3
7ðk P 2Þ 7ðk P 2Þ 7ðk P 4Þ 3
r=5
0 1 2 3 4
9 7 6 5 4
9 9 9 7 4
9ðk P 2Þ 9ðk P 2Þ 9ðk P 2Þ 9ðk P 4Þ 4
r=6
0 1 2 3 4 5
11 9 8 7 6 5
11 11 11 11 8 5
11ðk P 2Þ 11ðk P 2Þ 11ðk P 2Þ 11ðk P 2Þ 11ðk P 6Þ 5
458
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Fig. 1. New mapping functions for WENO-IM (2, A) with d2 ¼ 0:3, and A ¼ 0:1; 1 and 10 respectively.
Fig. 2. New mapping functions for WENO-IM (k, 1) with d2 ¼ 0:3, and k ¼ 2; 4 and 6 respectively.
upstream central scheme in a larger region. That is, in smooth regions, the mapped WENO scheme can get better solution by decreasing A or increasing k. On the other hand, the mapping function may amplify the effect from the non-smooth stencils and thus cause a potential loss of accuracy near discontinuities [4]. So we have to find a balance between them. The balance means that the mapped scheme can obtain more accurate solutions while still retaining the non-oscillatory property of WENO method. 3.3.1. The choice of parameter k In Table 1, we prove that, in theory, WENO-IM (k, A) can achieve more accurate numerical solutions when a larger k is chosen for higher order WENO method. However, extensive numerical tests show that, WENO-IM (k > 2, A) methods do not usually provide better solutions than WENO-IM (k ¼ 2, A) no matter which order of WENO method is used. This is an interesting observation. We will now give an explanation. From (14), because k is positive and even, qð0Þ ¼ di , a simple calculation demonstrates that
g 0i ð0Þ ¼ 1 þ
1 k1
Adi
:
ð24Þ
When increasing k; g 0i ð0Þ becomes larger, that means the mapping function amplifies the effect from the non-smooth stencils much faster. This effect will be much larger while the mapping function is applied to higher order WENO schemes, because
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
459
the smallest optimal weight which they possess is much smaller. For example, when the fifth-order WENO-IM (k, A) method is used, the new weight for the non-smooth stencil can be ð1 þ 10k1 =AÞ times as large as the one by the fifth-order WENO-JS. Furthermore, when the seventh-order WENO-IM (k, A) method is used, the new weight for the non-smooth stencil can be ð1 þ 35k1 =AÞ times as large as the one by the seventh-order WENO-JS. This is probably why the use of larger values in k dose not necessarily generate better numerical solutions. Therefore we only consider the case when k ¼ 2 in the rest of the paper. 3.3.2. The choice of parameter A While the new mapping function is applied to the seventh-order WENO schemes, WENO-IM (2, 1) which is the method in Henrick et al. [10], seems to be a good choice. That is, WENO-IM (2, A > 1) does not provide more accurate solutions than WENO-IM (2, 1) as A is increasing, and WENO-IM (2, A < 1) always generates the oscillatory by decreasing A, because the effects from the non-smooth stencils become much larger than the one in [10]. Once again, it verifies the additional requirement g 0i ð0Þ ¼ 0 in [4] is necessary, especially for higher order WENO schemes. While the new mapping function is applied to the fifth-order WENO schemes, extensive numerical experiments demonstrate that WENO-IM (2, A > 1) always provide less accurate solutions than WENO-IM (2, 1) as A is increasing, and WENO-IM (2, A < 1) always provide more accuracy solutions than WENO-IM (2, A ¼ 1) as A is decreasing. However while A is too small, oscillatory starts to be generated. Extensive numerical tests show that, no matter in short time or long time simulation, A ¼ 0:1 seems to be a very good choice. In other words, it provides an improved mapped WENO method with less dissipation
Fig. 3. New mapping functions for the fifth order WENO-IM (2, 0.1).
Fig. 4. Comparison for the mapping functions for WENO-IM (2, 0.1) and WENO-IM (2, 1) for d0 ¼ 0:1.
460
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
and higher resolution. In our numerical tests, we only consider the fifth-order WENO-IM (2, 0.1) method with the mapping function: 3
g i ð xÞ ¼ d i þ
0:1ðx di Þ 2
0:1ðx di Þ þ xð1 xÞ
ð25Þ
:
Fig. 3 gives the graph for the mapping functions for WENO-IM (2, 0.1). Fig. 4 gives the comparison for the mapping functions of WENO-IM (2, 0.1) and WENO-IM (2, 1) for d0 ¼ 0:1. Clearly WENO-IM (2, 0.1) tends to use the optimal weights in a larger range. This may be the reason why it can provide better numerical solutions when the mapping function is applied to the fifth order WENO scheme. However, when the mapping function is applied to higher order (order P 7) WENO methods, g 0i ð0Þ becomes very large. This starts to contaminate the error for the numerical solution. The mapping function (25) is unstable to provide the better solution than WENO-IM (2, 1), let alone WENO-PM6 [4]. 4. Numerical tests In this section, we apply WENO-JS ( ¼ 106 ), WENO-Z ( ¼ 1040 ), WENO-PM6 ( ¼ 1060 ) and WENO-IM (k, A) ( ¼ 1040 ) to solve the conservation laws in one dimension, where WENO-Z is used to refer to the methods in [2,3]. The Lax-Friedrichs flux splitting is used in the numerical tests. The errors in the L1 norm is computed by comparing the numerical solution u with the exact solution uexact , according to
kerrork1 ¼
N 1 X jui ðuexact Þi j; N þ 1 i¼0
where N is the number of intervals in space. We have mentioned in previous section, that for the seventh-order WENO-IM (k, A) methods, WENO-IM (2, 1) seems to be the better choice. However, it still can not compete with WENO-PM6 as shown in [4]. So we only consider the fifth-order WENO-IM (k, A) methods in this section. We now focus on the performance of fifth-order WENO-IM (2, 0.1) in this section. 4.1. Scalar conservation laws in one dimension Example 1. We solve the linear advection equation:
ut þ ux ¼ 0;
0 6 x 6 2;
ð26Þ
uðx; 0Þ ¼ sinðpxÞ;
together with the periodic boundary conditions. Table 2 shows the convergence rates and errors for numerical solutions of Example 1 in L1 norm at output time t ¼ 1. The first row is the error of L1 norm, and the second row is the convergence rates. The CFL number is chosen to be ðMxÞ2=3 in order that the error for the overall scheme is a measure of the spatial convergence only. We observe that, all the schemes reach the desired fifth order of accuracy. Example 2. We solve the linear advection equation:
ut þ ux ¼ 0;
0 6 x 6 2;
ð27Þ
4
uðx; 0Þ ¼ sin ðpxÞ; together with the periodic boundary conditions.
Table 2 Convergence study for various schemes as applied to Example 1 at output time t ¼ 1, all schemes are used with CFL¼ ðMxÞ2=3 . N
40
80
160
320
WENO-JS
2.2864e-005
WENO-Z
3.8994e-006
WENO-IM (2, 1)
3.8991e-006
WENO-PM6
3.8916e-006
WENO-IM (2, 0.1)
3.8924e-006
7.1493e-007 4.9991 1.2328e-007 4.9832 1.2328e-007 4.9831 1.2326e-007 4.9806 1.2326e-007 4.9809
2.2438e-008 4.9938 3.8770e-009 4.9909 3.8770e-009 4.9909 3.8769e-009 4.9907 3.8769e-009 4.9907
7.0199e-010 4.9983 1.2154e-010 4.9954 1.2154e-010 4.9954 1.2154e-010 4.9954 1.2154e-010 4.9954
461
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Table 3 shows the convergence rates and errors for numerical solutions of Example 2 in L1 norm at output time t ¼ 1. The first row is the error of L1 norm, and the second row is the convergence rates. The CFL number is chosen to be ðMxÞ2=3 in order that the error for the overall scheme is a measure of the spatial convergence only. In the aspect of convergence rates, the WENO-JS scheme achieves the desired fifth order of accuracy when N P 160, and all the other schemes do not reach it. In the aspects of accuracy, The WENO-IM (2,0.1) provides the most accurate numerical solutions and is more accurate than other schemes almost by a factor of two at least when N P 80. Example 3. We solve the linear advection equation:
ut þ ux ¼ 0;
1 6 x 6 1;
ð28Þ
uðx; 0Þ ¼ u0 ðxÞ; where
81 ðGðx; b; z dÞ þ Gðx; b; z þ dÞ þ 4Gðx; b; zÞÞ; > 6 > > > > 1; > < u0 ¼ 1 j10ðx 0:1Þj; > > 1 > > > 6 ðFðx; a; a dÞ þ Fðx; a; a þ dÞ þ 4Fðx; a; aÞÞ; > : 0;
0:8 6 x 6 0:6; 0:4 6 x 6 0:2; ð29Þ
0 6 x 6 0:2; 0:4 6 0:6; otherwise;
and 2
Gðx; b; zÞ ¼ ebðxzÞ ;
Fðx; a; aÞ ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi maxð1 a2 ðx aÞ2 ; 0Þ:
ð30Þ
together with the periodic boundary condition. The constants are taken as a ¼ 0:5; z ¼ 0:7; d ¼ 0:0005; a ¼ 10 and log2 b ¼ 36d 2 . The solution contains a smooth but narrow combination of Gaussians, a square wave, a sharp triangle wave and
a half ellipse. Table 4 and 5 show the convergence rates and errors for numerical solutions of Example 3 in L1 norm at output time t ¼ 20 and t ¼ 2000 respectively. The first row is the error of L1 norm, and the second row is the convergence rates. The CFL number is chosen to be 0:1. When t ¼ 20, all the schemes converge with first order accuracy, and WENO-IM (2, 0.1) gives more accurate numerical solutions than other schemes. When t ¼ 2000, WENO-JS, WENO-Z and WENO-IM (2, 1) reduce their convergence rates seriously, especially, WENO-Z gives the negative convergence rate on the mesh of N ¼ 400, however, WENO-PM6 and WENO-IM (2, 0.1) still converge with first order accuracy, and on the mesh of N ¼ 400; 800 cells, WENOPM6 gives the best numerical solutions, WENO-IM (2, 0.1) gives slightly less accurate solutions than WENO-PM6, but more accurate solutions than WENO-JS, WENO-Z and WENO-IM (2, 1). Figs. 5–8 show the performance of various schemes at output time t ¼ 2000 on the mesh of N ¼ 400; 800 cells. On the mesh of N ¼ 400 cells, WENO-PM6 gives the highest resolution of Gaussians, however generates a small oscillation near x ¼ 0:5, and WENO-IM (2, 0.1) gives the highest resolution of square wave. On the mesh of N ¼ 800 cells, WENO-PM6 still gives the highest resolution of Gaussians, however generates a small oscillation near x ¼ 0:7, at this stage, WENO-PM6 and WENO-IM (2, 0.1) give the same resolution of square wave, and do not smear this feature. Table 6 shows the CPU time (in seconds) of various schemes at output time t ¼ 10 on the mesh of N ¼ 200; 400; 800 and 1600 cells respectively. The Matlab Language is used to devise the program for all the schemes, all the results presented here are computed on a PC with an Intel (R) Core (TM) i7–3770 K desktop processor and with 16.0 GB of memory. In order to ensure fairness in the comparison of CPU time, all the schemes share the same subroutine calls and are compiled with the same compilation options including optimization ones. The only differences between the implementation of all the schemes are on the subroutines for computing the different related WENO weights. We observe that WENO-Z is the most efficient scheme. However, it provides less accurate numerical solution than WENO-IM (2, 0.1). WENO-IM (2, 0.1) is slightly
Table 3 Convergence study for various schemes as applied to Example 2 at output time t ¼ 1, all schemes are used with CFL¼ ðMxÞ2=3 . N
40
80
160
320
WENO-JS
3.6393e-003
WENO-Z
1.4177e-003
WENO-IM (2, 1)
1.6499e-003
WENO-PM6
2.2364e-003
WENO-IM (2, 0.1)
1.3977e-003
4.7793e-004 2.9288 1.0794e-004 3.7153 1.4942e-004 3.4649 1.5091e-004 3.8894 6.7117e-005 4.3802
1.4829e-005 5.0103 5.5845e-006 4.2726 7.7330e-006 4.2722 8.3183e-006 4.1813 2.8832e-006 4.5409
3.8676e-007 5.2608 2.6142e-007 4.4170 3.5976e-007 4.4259 3.8574e-007 4.4306 1.0827e-007 4.7349
462
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Table 4 Convergence study for various schemes as applied to Example 3 at output time t ¼ 20, all schemes are used with CFL¼ 0:1. N
100
200
400
800
WENO-JS
1.2152e-001
WENO-Z
9.4926e-002
WENO-IM (2, 1)
9.6111e-002
WENO-PM6
9.7839e-002
WENO-IM (2, 0.1)
9.4494e-002
5.7747e-002 1.0733 3.8757e-002 1.2923 4.0298e-002 1.2540 4.0478e-002 1.2733 3.7664e-002 1.3270
2.5255e-002 1.1932 1.7682e-002 1.1322 1.8845e-002 1.0965 1.7411e-002 1.2171 1.7079e-002 1.1410
1.4352e-002 0.8153 8.3518e-003 1.0821 9.5250e-003 0.9844 8.1661e-003 1.0923 8.0227e-003 1.0901
Table 5 Convergence study for various schemes as applied to Example 3 at output time t ¼ 2000, all schemes are used with CFL¼ 0:1. N
100
200
400
800
WENO-JS
3.1104e-001
WENO-Z
2.5132e-001
WENO-IM (2, 1)
2.5496e-001
WENO-PM6
2.4304e-001
WENO-IM (2, 0.1)
2.3976e-001
2.9819e-001 0.0608 1.5717e-001 0.6772 1.9836e-001 0.3621 1.1025e-001 1.1404 1.1019e-001 1.1215
2.8040e-001 0.0888 1.9006e-001 0.2742 1.6024e-001 0.3079 5.2225e-002 1.0780 5.6666e-002 0.9595
2.4749e-001 0.1801 1.3242e-001 0.5214 1.5803e-001 0.0201 2.2554e-002 1.2113 2.6075e-002 1.1198
Fig. 5. Performance of the fifth-order WENO-JS, WENO-Z and WENO-IM (2, 0.1) for Example 3 when t ¼ 2000. The mesh of N ¼ 400 cells and CFL¼ 0:1 are used for all schemes.
slower than WENO-IM (2, 1) because of the additional multiplication with A, but is much faster than WENO-PM6 since WENO-PM6 is not very suitable for the vector operation due to the implementation of logic operations. 4.2. Euler system in one dimension We solve one-dimensional Euler equations of gas dynamics, but with different initial and boundary conditions. For the one-dimensional Euler equations, we have the conservation of mass, momentum and energy:
0
1
0
1
q qu B C B C @ qu A þ @ qu2 þ p A ¼ 0; E
t
uðE þ pÞ
x
ð31Þ
463
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Fig. 6. Performance of the fifth-order WENO-IM (2, 1), WENO-PM6 and WENO-IM (2, 0.1) for Example 3 when t ¼ 2000. The mesh of N ¼ 400 cells and CFL¼ 0:1 are used for all schemes.
Fig. 7. Performance of the fifth-order WENO-JS, WENO-Z and WENO-IM (2, 0.1) for Example 3 when t ¼ 2000. The mesh of N ¼ 800 cells and CFL¼ 0:1 are used for all schemes.
where q; u; E; p are the density, velocity, total energy and pressure, respectively. The system of Euler equations is closed by the equation of state for an ideal polytropic gas
E¼
p
c1
1 þ qu2 ; 2
ð32Þ
where c is the ratio of specific heat. For the examples in this paper, we use c ¼ 1:4. The numerical flux F^jþ1 at xjþ1 in (2) is calculated as follow [14]: 2
2
(1) Compute an average state ujþ1 using a Roe average; 2 (2) Compute the right eigenvectors, the left eigenvectors and the eigenvalues of the Jacobian f 0 ðujþ1 Þ, and denote them by 2
R ¼ Rðujþ1 Þ; 2
R
1
1
¼ R ðujþ1 Þ; 2
K ¼ Kðujþ12 Þ;
(3) Transform the flux and the solution which are in the potential stencil of the WENO reconstructions for obtaining the flux F^jþ1 , to the local characteristic fields by 2
vj ¼ R
1
uj ;
g j ¼ R1 f ðuj Þ;
j in a neighborhood of i;
(4) Use the Lax-Friedrichs flux splitting and the WENO reconstruction procedure for each component of the characteristic variables, to obtain the corresponding component of the flux b ; g jþ1 2
464
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Fig. 8. Performance of the fifth-order WENO-IM (2, 1), WENO-PM6 and WENO-IM (2, 0.1) for Example 3 when t ¼ 2000. The mesh of N ¼ 800 cells and CFL¼ 0:1 are used for all schemes.
Table 6 CPU time (in seconds) study for various schemes as applied to Example 3 at output time t ¼ 10, all schemes are used with CFL¼ 0:1. N
200
400
800
1600
WENO-JS WENO-Z WENO-IM (2, 1) WENO-PM6 WENO-IM (2, 0.1)
10.6 10.1 11.6 14.3 12.3
29.4 27.7 32.3 41.1 33.8
95.9 88.2 105.0 136.7 110.3
379.5 350.4 412.5 530.1 429.4
(5) Transform back into physical space by using
b F jþ1 ¼ R b g jþ1 ; 2
2
(6) Form the flux by taking
b F jþ1 ¼ b F þjþ1 þ b F jþ1 2
2
2
Example 4 (Shock tube problem with a sonic point). We solve the Euler Equations (31) on ½5; 5 with following initial condition [20]
ðq; u; pÞ ¼
ð1:0; 0:75; 1:0Þ; if x < 2; ð0:125; 0:0; 0:1Þ; if x P 2:
ð33Þ
This is a modified version of Sod’s problem [18], its solution has a right shock wave, a right traveling contact wave and a left sonic rarefaction wave. We compute the solution at the output time t ¼ 2 on the mesh of N ¼ 100; 200; 400 and 400 cells. The CFL number is chosen to be 0:4. Table 7 shows the errors of density for various WENO schemes in L1 norm, and WENO-IM (2, 0.1) provides the most accurate numerical solutions on different meshes. Fig. 9 is the plot for the density on the mesh of N ¼ 100. We observe that, all the schemes perform well at the sonic point.
Table 7 Performance of different fifth-order WENO schemes on shock tube problem with a sonic point, CFL¼ 0:4. N
100
200
400
800
WENO-JS WENO-Z WENO-IM (2, 1) WENO-PM6 WENO-IM (2, 0.1)
8.6385e-003 6.9880e-003 7.6275e-003 7.5799e-003 6.8783e-003
4.5876e-003 3.6937e-003 4.0102e-003 3.9629e-003 3.6330e-003
2.4777e-003 2.0126e-003 2.1624e-003 2.1260e-003 1.9772e-003
1.3338e-003 1.0828e-003 1.1562e-003 1.1345e-003 1.0636e-003
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
465
Fig. 9. Density plot for the shock tube problem with a sonic point when t ¼ 2. The mesh of N ¼ 100 cells and CFL¼ 0:4 are used for all schemes.
Fig. 10. Performance of the fifth-order WENO-IM (2, 0.1) for modified shock/turbulence interaction when t ¼ 5. The mesh of N ¼ 1500 cells and CFL¼ 0:4 are used.
Fig. 11. Comparison of the fifth-order WENO-JS and WENO-IM (2, 0.1) for modified shock/turbulence interaction when t ¼ 5. The mesh of N ¼ 1500 cells and CFL¼ 0:4 are used for all schemes.
466
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
Fig. 12. Comparison of the fifth-order WENO-Z and WENO-IM (2, 0.1) for modified shock/turbulence interaction when t ¼ 5. The mesh of N ¼ 1500 cells and CFL¼ 0:4 are used for all schemes.
Fig. 13. Comparison of the fifth-order WENO-IM (2, 1) and WENO-IM (2, 0.1) for modified shock/turbulence interaction when t ¼ 5. The mesh of N ¼ 1500 cells and CFL¼ 0:4 are used for all schemes.
Fig. 14. Comparison of the fifth-order WENO-PM6 and WENO-IM (2, 0.1) for modified shock/turbulence interaction when t ¼ 5. The mesh of N ¼ 1500 cells and CFL¼ 0:4 are used for all schemes.
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
467
Example 5 (Modified shock/turbulence interaction). We solve the Euler equations (31) on ½5; 5 with following initial condition [19,21]
ðq; u; pÞ ¼
ð1:515695; 0:523346; 1:80500Þ; if x < 4:5; ð1 þ 0:1 sinð20pxÞ; 0:0; 1Þ;
if x P 4:5:
ð34Þ
The initial condition is set by a Mach 1:1 shock interacting with a high-frequency density perturbation. As time evolves, the shock moves into this perturbation, which spreads upstream. We compute the solution at the output time t ¼ 5 on the mesh of N ¼ 1500 cells. The CFL number is chosen to be 0:4. The density computed by WENO-IM (2, 0.1) is given in Fig. 10, the reference solution is an approximation to the exact solution using WENO-JS with N ¼ 10000. Fig. 11–14 show the comparisons of WENO-IM (2, 0.1) with other schemes on ½2; 1, we observe that WENO-IM (2, 0.1) gives the best resolution. 4.3. Summary of performance of fifth order WENO-IM (2, 0.1) We summarize the performance of our method as follow. (1) All the improved WENO methods perform better than the original WENO method. (2) Our method performs better than other improved WENO methods on smooth problems with large gradients. As presented in Example 2, the solution of our method is more accurate than other methods, in fact by a factor of two at least. (3) The most important advantage of our method presents on problems with discontinuities. The improvement is obvious if a long simulation is desired or if the solutions have many high-frequency waves or dense discontinuities. Example 3 demonstrates the improvement of our method if a long simulation is desired (see Table 5). And Example 5, the modified shock/turbulence interaction, demonstrates the improvement of our method by the mass high-frequency density perturbations. Note that the high-frequency perturbations are actually numerical discontinuities if a coarse mesh is used. Our method performs much better than WENO-Z and WENO-IM (2, 1) in these two examples. (4) Our method gives comparable solution with WENO-PM6 in Example 3, and gives more accurate solution than WENOPM6 in Example 5. In terms of CPU time, our method runs much faster then WENO-PM6. (5) The mapping function of our method is smoother than the one in WENO-PM6, because our mapping functions use rational functions but the mapping functions in WENO-PM6 use piecewise polynomial. (6) If the major parts of the numerical solutions are smooth while with a few discontinuities, and if a short time simulation is desired; e.g. Table 4 for Example 3 and Table 7 for Example 4, the Sod-type problem, the improvement of our method may be only minor comparing to other modified WENO methods. 5. Conclusions In this paper, a new mapping function is proposed, then an improved version of the mapped WENO method in Henrick et al. [10] is developed, we call it WENO-IM (k, A). The improved mapped WENO method can achieve the optimal order of accuracy near critical points by proper choice of parameter k. When the new mapping function is applied to the seventhorder WENO method, WENO-IM (2, 1), which is the same method proposed by Henrick et al. [10], seems to be a good choice. However, when the new mapping function is applied to the fifth-order method, WENO-IM (2, 1) is not the best choice any more. Numerical tests show that, no matter in short time or long time simulation, WENO-IM (2, 0.1) seems to be a better choice for the fifth-order WENO method. It provides better numerical solutions with less dissipation and higher resolution than WENO-JS, WENO-Z and WENO-IM (2, 1), special in long time simulation, and provides comparable numerical solutions with WENO-PM6, but it is more efficient than WENO-PM6 in terms of CPU time. Acknowledgments The authors would like to thank the anonymous referees for valuable suggestions and comments. The work of H. Feng was supported by National Natural Science Foundation of China (Nos. 91130022, 10971159 and 11161130003) and NCET of China. References [1] D.S. Balsara, C.W. Shu, Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy, J. Comput. Phys. 160 (2000) 405–452. [2] R. Borges, M. Carmona, B. Costa, W.S. Don, An improved weighted essentially non-oscillation scheme for hypebolic conservation laws, J. Comput. Phys. 227 (2008) 3191–3211. [3] M. Castro, B. Costa, W.S. Don, High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws, J. Comput. Phys. 230 (2011) 1766–1792. [4] H. Feng, F.X. Hu, R. Wang, A new mapped weighted essentially non-oscillatory scheme, J. Sci. Comput. 51 (2012) 449–473. [5] S. Gottlieb, C.W. Shu, E. Tadmor, Strong stability-preserving time discretization methods, SIAM Rev. 43 (2001) 89–112. [6] A. Harten, ENO schemes with subcell resolution, J. Comput. Phys. 83 (1987) 148–184. [7] A. Harten, B. Engquist, S. Osher, S. Chakravarthy, Uniformly high order essentially non-oscillatory schemes III, J. Comput. Phys. 71 (1987) 231–303.
468
H. Feng et al. / Applied Mathematics and Computation 232 (2014) 453–468
[8] A. Harten, S. Osher, Uniformly high order essentially non-oscillatory schemes I, SIAM J. Numer. Anal. 24 (1987) 279–309. [9] A. Harten, S. Osher, B. Engquist, S. Chakravarthy, Some results on uniformly high order accurate essentially non-oscillatory schemes, Appl. Numer. Math. 2 (1986) 347–377. [10] A.K. Henrick, T.D. Aslam, J.M. Powers, Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points, J. Comput. Phys. 207 (2005) 542–567. [11] G.S. Jiang, C.W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126 (1996) 202–228. [12] X.D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994) 200–212. [13] C.W. Shu, Numerical experiments on the accuracy of ENO and modified ENO schemes, J. Sci. Comput. 5 (1990) 127–149. [14] C.W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, in: Advanced numerical approximation of nonlinear hyperbolic equations, Lecture Notes in Mathematics, vol. 1697, Springer-Verlag, Berlin, 1998, pp. 325–432. [15] C.W. Shu, High order weighted essentially non-oscillatory schemes for convection dominated problems, SIAM Rev. 51 (2009) 82–126. [16] C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, J. Comput. Phys. 77 (1988) 439–471. [17] C.W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes II, J. Comput. Phys. 83 (1989) 32–78. [18] G.A. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws, J. Comput. Phys. 107 (1978) 1–31. [19] V.A. Titarev, E.F. Toro, Weno schemes based on upwind and centred tvd fluxes, Comput. Fluids 34 (2005) 705–720. [20] E.F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, third ed., Springer-Verlag, Berlin Heidelberg, 2009. [21] E.F. Toro, V.A. Titarev, Tvd fluxes for the high-order ader schemes, J. Sci. Comput. 24 (2005) 285–309.