On the numerical overshoots of shock-capturing schemes
OOn the numerical overshoots of shock-capturing schemes
Huaibao Zhang a , Fan Zhang b, ∗ a School of Physics, Sun Yat-sen University, Guangzhou 510006, China b Centre for mathematical Plasma-Astrophysics, Department of Mathematics, KU Leuven, Celestijnenlaan 200B, 3001 Leuven, Belgium
Abstract
This note introduces a simple metric for benchmarking shock-capturing schemes, and this metric is especially focus-ing on the shock-capturing overshoots, which may damage the robustness of numerical simulations, as well as thereliability of numerical results. The idea is to numerically solve the linear advection equation with square waves ofdifferent wavenumbers as initial solutions. Therefore, the exact overshoots after one step of temporal evolution canbe easily calculated, and shown as a function of CFL number and reduced wavenumber. The present metric is able toquantify the numerical overshoots of shock-capturing schemes, and some new features of several canonical schemesare shown.
Keywords: overshoot; shock-capturing; high-order; TVD; WENO
1. Introduction
Numerical simulations of compressible flows are a highly challenging topic because shock waves which causestrong density, velocity and pressure discontinuities, need to be tackled by using robust [1], accurate [2, 3] and efficient[4, 5] shock-capturing schemes, the development of which is still not trivial. This topic becomes especially challengingwhile high numerical resolution and computational robustness are expected to be provided by a unified framework ofnumerical schemes, since achieving high-resolution usually expects high-order polynomials for spatial approximation,and high-order polynomials tend to be oscillatory if they are crossing discontinuities (Gibbs phenomenon). Theamplitude of these kind of oscillatory solutions does not decrease even as the grid is refined. Therefore, in the pastseveral decades, numerous high-resolution shock-capturing schemes have been designed and investigated [1, 2, 3, 6, 7].In particular, total variation diminishing (TVD) schemes [1] enforces the condition of non-increase of the totalvariation in time, usually with using flux or slope limiters, thus achieving oscillation-free shock-capturing solutions.However, it is well known that TVD schemes have reduced order of accuracy at smooth extrema of the solutions. Uni-form high-order accuracy can be achieved by the essentially non-oscillatory (ENO) scheme which uses adaptive sten-cils to avoid or, at least, minimize reconstructions crossing discontinuities, thus reducing the oscillation or overshootessentially. An important improvement, the weighted essentially non-oscillatory (WENO) scheme [3], uses a convex-combination procedure of reconstructions carried out over given amount of candidate stencils, and the shock-capturingperformance is highly depending on the nonlinear weights which in essence represent the relative local smoothness of ∗ Corresponding author.
Email address: [email protected] (Fan Zhang)
Preprint submitted to X February 22, 2021 a r X i v : . [ phy s i c s . f l u - dyn ] F e b he solution on each candidate stencil. WENO schemes [3, 6, 7] maintain the ENO property and significantly improvethe overall resolution. However, ENO schemes or WENO schemes do not completely remove numerical oscillation orovershoot in shock-capturing computations [8]. Therefore, in order to guarantee the robustness of WENO schemes,great efforts have been made and the resulting schemes can pass severe test cases involving blast waves and/or nearvacuum [9, 10].However, in this note, there is not intention of improving or reviewing the robust high-order shock-capturingschemes. It should be noted that, to the Authors’ knowledge, the properties of numerical schemes mentioned abovewere still evaluated based on a case-by-case basis, and rather than a quantitative analysis, visualized investigation andcomparison of numerical results are the most common evaluation methodology. Only until recently, a quantitativemetric for shock-capturing error has been given by Zhao et al. [8], showing that the ENO and WENO schemes mayyield significant over-amplification of incoming waves at discontinuities, which may be caused by the activation oflinearly unstable stencils in the vicinity of the discontinuity. Especially, their analysis shows that the over-amplificationhappens within a certain range of wavenumber space, which means that the traditional evaluation with simple shocksimulations may not be able to trigger the over-amplification behaviour.Despite the success of the shock-capturing error metric of Zhao et al. [8], in this note we would like to further ex-plore more details of those popular shock-capturing schemes, by using a simple test. For example, in shock-capturingcomputations the CFL number is sometimes restricted to guarantee the robustness, and we would like to unveil thedetailed behaviour while the CFL number is changing. Moreover, we are interested in the behaviour of simulatingmultiple closely located shock waves, which was not given in Ref.[8]. Another fundamental property needs to be justi-fied is the relative magnitude of the overshoots produced by the shock-capturing schemes at discontinuities. Therefore,at least these issues need to be investigated but haven’t been introduced, to the Authors’ knowledge.The organization of the paper is as follows. The idea is introduced in section 2, and then in section 3 several typicalshock-capturing schemes are investigated by the present metric. Eventually, the concluding remarks are given in thelast section.
2. The numerical technique
In this note, we investigate the IVP of one-dimensional wave propagation in an unbounded domain, governed bythe scalar linear advection equation, i.e. ∂ u ∂ t + a ∂ u ∂ x = , (1)where the characteristic velocity a is assumed to be 1, without loss of generality. In practice, the spatial discretizationof Eq.(1) is given on an equally spaced grid, resulting an ordinary differential equation (ODE) system, i.e. du i dt = − ∂ u ∂ x | x = x i , ≤ x i ≤ L , i = , · · · , N , (2)and the periodic boundary condition is given in the left ( x =
0) and right ( x = L ) boundaries to model the unboundeddomain. To be more specific, in this work we have N =
500 and L = u ( x , T ) = . ( sin ( π x / λ n )) , (3)where T =
0, and λ n is the wave length, which is defined as λ n = L / n , n = , · · · , N / . (4)Therefore, the exact solutions of this linear advection problem can be given by simply shifting the initial solutions by adistance of aT , where T is the time at which the computations end. Moreover, the reduced wavenumber which relatesto the frequency of the initial square waves is defined as ϕ = ω ∆ x = ( π / λ ) ∆ x , (5)where ω is wavenumber. The reduced wavenumber is to be used to specify the performance. In this work, we focus on the shock-capturing schemes which are used for spatial solutions, i.e. the right hand sideof Eq.(2). The partial derivatives in x -direction are approximated by using the conservative finite difference formula,i.e. du i dt = − ∆ x ( h i + / − h i − / ) . (6)The flux function h i ± / at half points can be implicitly defined by f ( x ) = ∆ x (cid:90) x + ∆ x / x − ∆ x / h ( ξ ) d ξ , (7)and then the semi-discretized form can be written as du i dt ≈ − ∆ x ( ˆ f i + / − ˆ f i − / ) , (8)where ˆ f i + / = ˆ f ( f i − l , · · · , f i + r ) , (9)and ˆ f i − / can be easily achieved symmetrically. It should be noted that for shock-capturing schemes the function ˆ f is usually a nonlinear function, and l and r define the width of the stencil, which depends on the numerical schemesbeing used.In this work, the classical TVD schemes [1, 11] and the ENO/WENO schemes [1, 2, 3, 6, 7] are investigated, sincethey have several essential properties which are of interest. TVD schemes are oscillation-free in one-dimensional sim-ulations, and thus we would like to use the minmod (MM) limiter and superbee (SB) limiter, which are respectively themost and least restrictive second-order TVD schemes [11], as the baseline schemes to show the different behaviours ofhigh-order schemes. ENO/WENO schemes, however, can not guarantee providing oscillation-free shock-capturing re-sults, and Zhao et al. [8] proved that the amplitude of the oscillation relates to the frequency of the smooth disturbance3nteracting with the discontinuity. Zhao et al. [8] also indicated that the activation of linearly unstable stencils mightcontribute to the resulting exaggerated amplitude or overshoot of propagating smooth waves near the discontinuity.Here, we are especially interested in more detailed behaviours and try to quantify the overshoot caused by simplediscontinuities without interacting smooth waves, since this scenario may directly and specifically mimic high-orderreconstruction crossing shocks.Without any further discussion and modification on classical shock-capturing schemes, we directly use TVDschemes and ENO/WENO schemes as the flux function in Eq.(9) to solve the right hand side of the ODE system,i.e. Eq.(2), and then the temporal differential term in the left hand side is solved by using the third-order stronglystable Runge-Kutta method [12], of which more detailed formulas are also omitted here for simplicity. It is worthnoted that the temporal solutions, more or less, change the numerical overshoot of a specific spatial reconstructionscheme. However, in this work we specifically focus on the relative difference between those high-order schemes, andthus the commonly used temporal solution is applied without further examination. Theoretically, TVD schemes may provide oscillation-free shock-capturing results, and ENO/WENO schemes areless dissipative but may produce small oscillations. Therefore, compared with TVD schemes, ENO/WENO schemesmay produce sharp profile of discontinuities, and thus, although ENO/WENO schemes produce oscillation, the overallnumerical error of ENO/WENO schemes might be still smaller than that of TVD schemes. In this work, we are notinterested in the overall numerical error, especially we would like to focus on the overshoots near discontinuities whichprobably leads to negative pressure or density in gas dynamics simulations. Therefore, in the present framework, onlythe resulting exaggerated amplitude or overshoot will be investigated, and we ignore dissipative numerical behaviourssince they have been discussed in numerous literatures.In this work, we do not consider long time computations, and instead, only one step of temporal evolution will beconducted, with the initial field given in the subsection 2.1 and varying time steps T calculated based on CFL ∈ ( , . ] ,resulting ˆ u ( x , T ) = u ( x , T ) + u ε ( x , T )= . ( sin ( π ( x + aT ) / λ n )) + u ε ( x , T ) , (10)where a =
1, as mentioned above, and u ε ( x , T ) is the error function. It is straightforward to fetch the numericalovershoots from the error function as follows u overshoot ε ( x ) = | u ε ( x , T ) | , if | ˆ u ( x , T ) | > . , , otherwise , (11)which means that the numerical error caused by the dissipative effect is excluded from our discussion. For eachcombination of T and λ n , only the maximum overshoot, i.e. max ( u overshoot ε ( x )) , will be shown in the next section.It is well known that using large time steps will challenge the stability of numerical schemes, which usually doesn’t4uantitatively show the influence of involving spatial discontinuity. Therefore, the present simple metric is expectedto further explain numerical behaviours in shock-capturing computations.It should be noted again that in this work only one step of temporal evolution is taken into account, because forlong-time simulations numerical dissipation will inevitably smooth out all the discontinuities and smooth waves. Forinstance, Deng et al. [13] showed that WENO schemes may significantly smear smooth waves and discontinuities, andeven their low-dissipation high-order scheme cannot completely avoid this trend in a very long-time linear advectionsolution. At the meantime, even the most popular fifth-order WENO-JS scheme [6] may produce minor overshootin certain wavenumber range [8]. Therefore, by only considering one step of temporal evolution, we may be able toavoid mixing up all the numerical effects.
3. Results and discussion
As mentioned above, in this section we consider several representative shock-capturing schemes, including: (i)second-order TVD schemes [14, 11] with the minmod (MM) and superbee (SB) flux limiters; (ii) schemes of the ENOclass [15] with third-order and fifth-order of accuracy; (iii) the WENO-JS scheme with fifth-order [6] and seventh-order [9] of accuracy; (iv) the WENO-Z scheme with fifth-order [7] and seventh-order [16] of accuracy; and (v) themonotonicity preserving WENO (MPWENO) schemes with fifth-order and seventh-order of accuracy [9].As shown in Fig.1(a), at least for the present one-dimensional linear advection equation, the minmod limiterguarantees overshoot-free results even using a large CFL number. The superbee limiter produces minor overshoot ifthe wavenumber and CFL number are both large, as shown in Fig.1(b). It should be noted that the overshoot maychange a bit while using a different time stepping technique, but in this work we focus on the relative differencebetween numerical schemes, as mentioned above. In general, these two TVD limiters are able to eliminate oscillationnear discontinuities, and thus, the other TVD limiters, which are at least more restrictive than the superbee limiter,should be able to remove oscillation as well. (a) minmod (b) superbeeFigure 1: Numerical overshoots of TVD schemes.
It is also encouraging and surprising to find that the third-order ENO3 scheme shows overshoot-free solutions, as5hown in Fig.2(a). Whereas, the fifth-order ENO5 scheme produces significant error, as shown in Fig.2(b), which isexpectable since its five-point candidate stencils inevitably cross closely located discontinuities. The strong evidenceis that the overshoot is directly related to the reduced wavenumber. It can be found that if the reduced wavenumber isapproximately larger than 1 .
05, the overshoot error will suddenly increase regardless of the CFL number. However, ifthe CFL number is relative large, the numerical overshoot is still observable in low-wavenumber region ( ϕ < . (a) ENO3 (b) ENO5Figure 2: Numerical overshoots of ENO schemes. The result of the WENO-JS5 scheme, as shown in Fig.3(a), shows similar behaviours. In fact, as long as thediscontinuities are separated by a sufficiently large distance, in which the smoothness of the candidate stencils can beguaranteed, the overshoot error is rather small, even if a large CFL number is used. Moreover, the overshoot error isalso not monotonously depending on the CFL number. The WENO-JS7 scheme, in general, shows similar behaviours,as shown in Fig.3(b). However, it is worth noted that the WENO-JS7 scheme already shows noticeable overshoot ifthe reduced wavenumber is approximately larger than 0 .
8. We believe that the large stencil used by the WENO-JS7scheme is the reason causing this behaviour. It is also important that the amplitude of the overshoot produced by theWENO-JS7 scheme is significantly larger than that of the WENO-JS5 scheme.As shown in Fig.4, the results of the WENO-Z schemes have similar patterns as those of the WENO-JS schemes.However, the WENO-Z schemes, which are proved to be less dissipative compared with the WENO-JS schemes [16],produce more significant overshoots in general. In particular, the overshoots in the intermediate wavenumber regionare even more significant.The MPWENO schemes [9], which incorporate the monotonicity preserving bounds [17], have significantly re-duced the overshoot error, as shown in Fig.5. Especially, the MPWENO7 scheme, which is seventh-order accurate,still shows relative small overshoot even in the high-wavenumber region. Therefore, it is very effective to use themonotonicity preserving bounds for improving the shock-capturing capability of WENO schemes. Of course, the6 a) WENO-JS5 (b) WENO-JS7Figure 3: Numerical overshoots of WENO-JS schemes.(a) WENO-Z5 (b) WENO-Z7Figure 4: Numerical overshoots of WENO-Z schemes. (a) MPWENO5 (b) MPWENO7Figure 5: Numerical overshoots of MPWENO schemes.
More details are shown in Fig.6 and Fig.7. In these two figures, the previous results are projected into the CFLnumber space and the wavenumber space respectively. In Fig.6, the most interesting behaviour being shown is that sev-eral schemes show relative small overshoots if the CFL number is close to 0.6. In particular, the maximum overshootof the fifth-order WENO-Z5 scheme with CFL ≈ .
65, is approximately 30% smaller than the maximum overshootof the same scheme with CFL = .
4. It is worth noted that, although according to various canonical numerical testcases, we can expect that higher-order schemes are usually less robust compared with lower-order schemes, in thiswork we are able to quantitatively compare the overshoots produced by different shock-capturing schemes. Also, it isunsurprising that increasing the CFL number significantly affects the amplitude of overshoots, but by doing this testwe are able to quantitatively visualize the difference between these commonly used numerical schemes.In Fig.7, we can, again, easily find out the most robust high-order shock capturing schemes (MPWENO), asshown in the previous results. At the meantime, we can more clearly distinguish the influence of the wavenumber,or more specifically, the influence of closed located discontinuities. For instance, all the WENO-type schemes startto produce significant overshoots if the reduced wavenumber is approximately larger than 1.05, as also mentionedabove. This behaviour is especially clear in Fig.7(b). We don’t want to describe every details in these figures, butwe may realise that the robustness of these schemes encounters challenges if the full stencil crosses more than onediscontinuity. Therefore, in order to improve the robustness of a shock-capturing scheme, this situation must bespecifically investigated. On the other hand, for a single discontinuity all the schemes being investigated here arecapable to produce overshoot-free results. It is worth noted again that in Ref.[8], shock-capturing error involvinginteracting oscillating smooth waves were discussed. 8 a) L ∞ (b) L Figure 6: The overshoot variations in the CFL number space.(a) L ∞ (b) L Figure 7: The overshoot variations in the wavenumber space. . Conclusions In this note we have designed a simple numerical metric for assessing shock-capturing schemes. Despite that themethod introduced only involves discontinuities without smooth waves, one may still use it to quantitatively evaluatethe shock-capturing capability of numerical schemes. According to the present results, several simple conclusions canbe directly observed: (i) the amplitude of overshoots is significantly depending on the distance between discontinuities;(ii) The amplitude of overshoots is not monotonously varying while the CFL number is increasing; and (iii) Themaximum overshoots of different schemes may be produced in very different scenarios, with respect to the CFLnumber and/or the reduced wavenumber (representing the distance between discontinuities).These observations are quantitatively provided, and thus we expect that the present method can be a new errormetric for shock-capturing schemes which need to be sufficiently robust for critical situations. Moreover, besides theconventional case-by-case investigation, in which the resolution and robustness are usually assessed altogether, thepresent note further provides a tool with which one may be able to find out the weakness or advantage of a shock-capturing scheme clearly. Of course, there are more factors affecting shock-capturing computations, especially inmore practical numerical simulations with complex physics, and thus more investigation is necessary.
References [1] A. Harten, On a class of high resolution total-variation-stable finite-difference schemes, SIAM J. Numer. Anal.21 (1) (1984) 1–23. doi:10.1137/0721001 .[2] A. Harten, B. Engquist, S. Osher, S. R. Chakravarthy, Uniformly high order accurate essentially non-oscillatoryschemes, III, Journal of Computational Physics 71 (2) (1987) 231–303. doi:10.1016/0021-9991(87)90031-3 .[3] X. D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, Journal of Computational Physics115 (1) (1994) 200–212. doi:10.1006/jcph.1994.1187 .[4] H. Zhang, F. Zhang, J. Liu, J. McDonough, C. Xu, A simple extended compact nonlinear scheme with adaptivedissipation control, Communications in Nonlinear Science and Numerical Simulation 84 (2020) 105191. doi:10.1016/j.cnsns.2020.105191 .[5] C. Xu, F. Zhang, H. Dong, H. Jiang, Arbitrary high-order extended eno schemes for hyperbolic conservationlaws, International Journal for Numerical Methods in Fluids n/a (n/a). doi:10.1002/fld.4968 .[6] G.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, Journal of Computational Physics126 (1) (1996) 202–228. doi:10.1006/jcph.1996.0130 .[7] R. Borges, M. Carmona, B. Costa, W. S. Don, An improved weighted essentially non-oscillatory scheme forhyperbolic conservation laws, Journal of Computational Physics 227 (6) (2008) 3191–3211. doi:10.1016/j.jcp.2007.11.038 . 108] G. Zhao, M. Sun, A. Memmolo, S. Pirozzoli, A general framework for the evaluation of shock-capturing schemes,Journal of Computational Physics 376 (2019) 924–936. doi:10.1016/j.jcp.2018.10.013 .[9] D. S. Balsara, C.-W. Shu, Monotonicity preserving weighted essentially non-oscillatory schemes with increas-ingly high order of accuracy, Journal of Computational Physics 160 (2) (2000) 405 – 452. doi:10.1006/jcph.2000.6443 .[10] X. Y. Hu, N. A. Adams, C.-W. Shu, Positivity-preserving method for high-order conservative schemes solvingcompressible euler equations, Journal of Computational Physics 242 (2013) 169 – 180. doi:https://doi.org/10.1016/j.jcp.2013.01.024 .[11] P. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal.21 (5) (1984) 995–1011. doi:10.1137/0721062 .[12] S. Gottlieb, C. Shu, E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Review43 (1) (2001) 89–112. doi:10.1137/S003614450036757X .[13] X. Deng, Y. Shimizu, F. Xiao, A fifth-order shock capturing scheme with two-stage boundary variation diminish-ing algorithm, Journal of Computational Physics 386 (2019) 323–349. doi:10.1016/j.jcp.2019.02.024 .[14] S. Osher, S. Chakravarthy, High resolution schemes and the entropy condition, SIAM Journal on NumericalAnalysis 21 (5) (1984) 955–984. arXiv:https://doi.org/10.1137/0721060 , doi:10.1137/0721060 .[15] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes,ii, J. Com-put. Phys. 83 (1) (1989) 32–78. doi:10.1016/0021-9991(89)90222-2 .[16] W.-S. Don, R. Borges, Accuracy of the weighted essentially non-oscillatory conservative finite differenceschemes, Journal of Computational Physics 250 (2013) 347 – 372. doi:10.1016/j.jcp.2013.05.018 .[17] A. Suresh, H. Huynh, Accurate monotonicity-preserving schemes with runge–kutta time stepping, Journal ofComputational Physics 136 (1) (1997) 83 – 99. doi:10.1006/jcph.1997.5745doi:10.1006/jcph.1997.5745