A hybrid ENO reconstruction with limiters for systems of hyperbolic conservation laws

We consider a new class of essentially non-oscillatory (ENO) piecewise polynomial reconstructions together with interpolants based on Monotone Upwind Schemes for Conservation Laws. We improve the second-order ENO polynomial reconstruction by choosing an additional point inside the stencil in order to obtain the highest accuracy when combined with various non-linear limiters. The resulting algorithms are based on only one stencil selection, and we show that they can be efficiently implemented with largest allowable CFL numbers using optimal strong stability-preserving Runge-Kutta time evolution methods. The numerical results indicate that in some cases the schemes yield errors smaller in magnitude as compared to the fourth-order ENO scheme.


Introduction
In this paper, we review and modify the essentially nonoscillatory (ENO) reconstruction of [1] for the approximate solution of hyperbolic conservation laws with initial data u(x, 0) = u 0 (x). For the sake of simplicity, we first consider scalar conservation laws, that is, d = 1 in (1). The discretisation of (1) is done on a uniform spatial grid where the cell I j =[ x j−1/2 , x j+1/2 ] has a width h. Also, let x j = 1 2 (x j−1/2 + x j+1/2 ) be the mid-cell grid point of I j . Integrating (1) over the cell I j leads to the semi-discrete conservation equation that samples solutions at cell centres and can be formulated as where the cell average of u on I j is given byū j ≡ 1 h I j u(x, t) dx and the exact flux isf j+ 1 2 ≡ f (u(x j+ 1 2 , t)). The computation of the flux at cell boundaries requires a polynomial reconstruction that has been the subject of considerable work [1][2][3][4].
*Correspondence: apeer@umail.utm.ac.mu 1 Department of Applied Mathematical Sciences, University of Technology, Mauritius, La Tour Koenig, Pointe-aux-Sables, Mauritius Full list of author information is available at the end of the article ENO reconstructions aim to achieve high accuracy in smooth regions in addition to resolving discontinuities with correct positions by adaptively selecting the smoothest stencil among several candidates without introducing spurious oscillations. The authors in [5,6] showed that the weighted ENO (WENO) schemes [3,4,7], both upwind or central based, fail to satisfy the sign property. Furthermore, they showed that ENO reconstructions and interpolation procedures are stable.
However, low-order ENO schemes have been observed to smear discontinuities [8] and smoothing up of corners [9]. Higher-order ENO reconstructions decrease damping but generate significant oscillations when solving hyperbolic systems unless costly characteristic decompositions are used [10]. Also, Xu and Shu [11] noted that it is not possible to maintain the non-oscillatory property and at the same time avoid the smearing of discontinuities by high-order ENO schemes.
The works of [12][13][14] demonstrated that hybridized schemes perform better than the existing schemes by taking advantages of each of their components. For example, the authors in [15][16][17] proposed the use of limiters in conjunction with ENO reconstructions. The numerical experiments described in those works show that the modified scheme reduced smearing near discontinuities and gave good resolutions of corners and local extremas. Furthermore, different ENO-type methods have been extended to multi-dimensional problems, see for example http://www.iaumath.com/content/7/1/15 [15,[18][19][20][21], and Serna and Marquina [15] noted that some ENO methods combined with limiters performed better than conventional ENO methods of similar order for multi-dimensional problems where fine structures appear to be important. This paper proposes a new hybrid scheme that reduces the damping near discontinuities. Our technique avoids the costly characteristic decomposition and is based on modifying the second-order ENO polynomial reconstruction to obtain a higher-order polynomial combined with non-linear limiters to avoid spurious oscillations. Instead of using a reconstruction over two cells for finding the fluxes of (2), we use an additional point within the stencil by means of a MUSCL-type interpolant [22] with numerical derivatives. Those derivatives depend on nonlinear limiters that vary according to the smoothness of its vicinity and improve the behaviour of the scheme near discontinuities.
We give an outline of this paper. We first describe the new reconstructions with various interpolants and limiters. We then give the results of numerical experiments for our schemes, and we extend the new class of hybrid schemes to solve hyperbolic systems of conservation laws. The conclusions on this work are presented in the final section.

A hybrid ENO-flux limiter scheme
We use a finite volume approach for one-dimensional equations, and we approximate u(x j , t) by u j . Assuming that the cell averagesū j are known at time t, Harten et al. [1] proposed a piecewise polynomial reconstruction u(x, t) = j p j (x) χ j (x) to solve (2), where χ j (x) is the characteristic function of the cell I j . In the secondorder ENO2 reconstruction, p j (x) is a linear polynomial on the cell I j , and it maintains conservation, that is, The two-cell stencils of ENO2 are obtained by either choosing r = 0 or r = 1 cell to the right of I j . Let Then interpolating V (x) at the cell boundaries x j−r+i−1/2 for i = 0, 1, 2 by using Newton's interpolation formula and differentiating the new polynomial, we get where V is a divided difference. The 0th-degree divided The ENO polynomials approximate the cell boundaries of (2) by taking u − j+1/2 = p j (x j+1/2 ) and u + j+1/2 = p j+1 (x j+1/2 ). An evolution step consists of collecting the point values u(x j+1/2 , t) from u − j+1/2 and u + j+1/2 , andf j+ 1 , where h(., .) is a monotone flux. Some possible choices for these class of schemes are the Godunov, Lax-Friedrichs and Roe with entropy fix fluxes. The advantages and disadvantages associated with these fluxes are discussed by [8]. In this paper, we will use the Lax-Friedrichs flux. We start the new reconstruction by selecting a two-cell stencil as for the ENO scheme. Let y i = x j+i−r−1/2 for i = 0, 1 and 2 denote the cell boundaries of the stencil. The next step consists of adding another point y 3 = x j ∈ I j . The interpolating polynomial is then given by In order to compute the divided differences of (4), we need a polynomial that retains information within the cells I j . Nessyahu and Tadmor [23] used a second-order MUSCL interpolant to compensate for the numerical viscosity introduced by the Lax-Friedrichs scheme [24]. Since 1 h I j L j (x)dx =ū j , polynomial (5) retains conservation. The divided difference V [ y 2 , y 3 ] given by is computed using the MUSCL interpolant (5), such that is the characteristic function of I j . The divided differences are found to be From (4) and (6), we deduce the following approximations at the cell boundaries x j+ 1 u + (1) u + (0) For the numerical derivative u j , we first consider the non-linear limiter of [23]: where the differences are given by ū j+1/2 =ū j+1 −ū j and 2ū j =ū j+1 − 2ū j +ū j−1 . The MinMod limiter MM is set as The classical MinMod limiter with θ = 1 (MM1) is second-order accurate and gives the numerical derivative which is non-oscillatory in the sense that Since MM1 may oversmear some discontinuities [23], we can instead choose θ = 2 (MM2) to have a steeper slope, unless oscillations are present, in which case we let u j = 0. This mechanism is aimed to reduce spurious oscillations allowed by other reconstructions like ENO [1] and WENO [3]. A Taylor expansion of the cell averages about the cell boundaries shows that the reconstructions (7), (8), (9) and (10) combined with the different values of the Min-Mod limiter vary from first-order near discontinuities up to third order in smooth regions.
In the numerical examples, we also combine the newly developed reconstructions with the following two nonlinear limiters: (a) UNO limiter. The accuracy of (12) drops at the non-sonic critical grid values u j , where [25,26]. The UNO limiter [27] adds second-order differences to the MinMod limiter (12) to achieve high accuracy including at critical points This feature retains information about the slopes of the solution. However, it uses a wider stencil with respect to existing ENO reconstructions of comparable accuracy. This stencil allows the limiter to avoid discontinuities, and in case extremas cannot be avoided, the accuracy of the UNO limiter decreases until a non-oscillatory approximation is obtained. Contrary to the power limiter of [15] which influences only second-order derivatives, the UNO limiter fully adapts the accuracy of our proposed polynomial reconstruction to the approximated profile. A Taylor expansion of the cell averages of the approximations (7) to (10) combined with the UNO limiter gives up to third-order accuracy. (b) Harmod limiter [15] is the harmonic mean of two differences of the same sign

Higher-order interpolant
A smoother limiter can be obtained for the proposed scheme by using a higher-order piecewise polynomial, which in turn can be obtained using the third-order nonoscillatory reconstruction due to [28]. This reconstruction seeks quadratic polynomials of the form: such that the piecewise parabolic reconstruction satisfies the properties of conservation, 1 h I j q j (x)dx =ū j . It is also third-order accurate, that is, In addition, the quadratic polynomial interpolates the two neighbouring cell averagesū j±1 . These three constraints uniquely determine the coefficients as a j = ū j − 1 24 2ū j , b j = 0ūj and c j = 1 2 2ū j . Finally, in order to avoid spurious extremas at the cell interfaces, a convex modification of the form The limiter θ j is constructed in terms of the following cell quantities: , , http://www.iaumath.com/content/7/1/15 and is given in [29] as otherwise.
The quadratic interpolant can then be written as follows: where u j = θ j 0ūj , u j = θ j 2ū j and u j =ū j − 1 24 u j . Considering the additional point y 3 = x j ∈ I j to compute the divided differences V [ y 2 , y 3 ] and using (15), we get (6). Thus, as before, we get the approximations (7), (8), (9) and (10) at the cell boundaries, and where now the numerical derivatives dependent on (14). We recover third-order accurate ENO spatial reconstructions when the parameter θ j = 1 in smooth regions.
We remark that the hybrid reconstructions can be extended to ENO schemes of any order by interpolating the primitive values at cell boundaries along with the additional point x j . The first-order divided difference containing x j of the interpolating polynomial is evaluated using a conservative piecewise polynomial with numerical derivatives of similar accuracy as the ENO scheme. Those numerical derivatives are chosen to be non-oscillatory in the sense of (W)ENO schemes [1,3] or to damp spurious oscillations as in [23,30].

Scalar test problems
We describe the results of numerical experiments using some scalar test problems. The scalar test problems have periodic boundary conditions and are solved on the interval [ −1, 1]. The pth-order ENO scheme is denoted by ENOp, and the hybrid schemes by the limiter used. The higher-order hybrid reconstruction with the quadratic polynomial (15) is denoted as Quadratic. All the reconstructions are combined with the Lax-Friedrichs flux [8] and with strong stability-preserving Runge-Kutta (SSPRK) methods [31][32][33]. The experiments are carried out in MATLAB environment on a 1.91-GHz processor with 992 MB of RAM.

Problem 2.
We consider the linear advection equation u t +u x = 0 with initial condition given by the square wave u(x, 0) = 1 for |x| < 1/3 and 0 elsewhere. Figure 1 shows the results obtained at time T = 4 using 50 equally spaced cells. We observe that the hybrid schemes approximate the exact solution with high accuracy for c = 1. ENO3 and ENO4 become unstable and diverge for c = 1. Figure 1b,c gives the ENO3 and ENO4 solutions for c = 0.5, where the solutions are damped near the discontinuities in both cases.
Problem 3. Our next experiment consists of the inviscid Burgers' equation u t + (0.5u 2 ) x = 0 with the discontinuous initial condition u(x, 0) = 1 for |x| < 1/3 and 0 elsewhere. In Table 2, we give the L 1 error and CPU time of the different solutions for T = 0.64. The time evolution process for the different schemes is done with the third-order SSPRK scheme and c = 0.8. We see that as expected, ENO1 produces the largest error. ENO2 and the hybrid scheme with the limiter MM1 produce numerical solutions of comparable accuracy. The remaining hybrid schemes give results of comparable accuracy to ENO3 and ENO4. The hybrid schemes with the MM1 and Harmod limiters have comparable CPU times to ENO2, whereas using the MM2 and UNO limiters give timings of the same order as ENO3. The Quadratic limiter is again the slowest of the different hybrid schemes and diverges on refined grids for large CFL numbers.

Systems of conservation laws
We extend our scheme to solve one-dimensional hyperbolic systems of conservation laws for (1) of the type The Jacobian A(u) of the flux F(u) has distinct real eigenvalues. At present, we pay more attention to solve the Euler equations of gas dynamics for a polytropic gas: Here ρ, q, p and E are respectively the density, velocity, pressure and total energy of the conserved fluid, and the ratio of the specific heats γ = 1.4.
There are two methods to extend the numerical schemes considered, namely by doing a componentwise extension and using characteristic decomposition. Liu and Osher [10] noted that high-order ENO/WENO schemes generate significant oscillations when solving hyperbolic systems unless costly characteristic decompositions are used. In the present work, we show that our new scheme is still efficient while adopting the less expensive componentwise extension for the stencil selection for each variable of U. The eigenvalues of the Jacobian matrix A are where a = √ γ p/ρ is the sound speed. We then use the Lax-Friedrichs flux for systems to find the values at the cell boundaries by the methods described previously We test the numerical schemes for the Euler equations before the perturbations in the solutions reach the boundary of the computational domain on the interval [ −5, 5] with c = 0.8 using the third-order SSPRK scheme.

Sod problem
We solve the equations of gas dynamics (17) with the initial conditions given by [34] U(

Lax problem
Next, we solve (17) using the initial condition of [35] U( For this more severe shock tube problem, we show the different approximations on 200 cells at T = 1.5 in Figures 8, 9, 10, 11, 12 and 13. Similar to Sod's problem, the MM1 limiter with the hybrid scheme gives comparable results to ENO2. The hybrid scheme with the MM2 limiter produces some overshooting in the approximation of the density. For componentwise extension, the UNO and Harmod limiters give an overall good resolution with little smearing of the approximations of shocks for the velocity and pressure profiles, whereas ENO3 allows more oscillations in the numerical results of the velocity.
In Table 3, we give the CPU times for the different numerical solutions of the shock tube problems. We note that speed relationship of the schemes with componentwise extension is similar to scalar case.

Multi-dimension extension
In the 1D experiments, we demonstrated that the hybrid scheme with UNO limiter is similar to ENO3 on most problems for similar CPU times. Now, we extend the hybrid scheme with UNO limiter to solve 2D equations which can be generalized to multi-dimensional problems. Following the recommendations of [8], we use a finite difference approach for such problems. This technique is computationally less expensive than the finite volume technique for which quadrature rules are necessary. A more detailed explanation of finite difference schemes can be found in [3,8].

Conclusion
The implementation of ENO schemes requires many selection statements. In this paper, we have proposed ENO-flux limiter schemes which reduce the number of selection steps. These schemes are shown to use the large CFL numbers allowed by the SSPRK methods for the time evolution step. The new methods were based on only one stencil selection, and we recovered thirdorder ENO reconstructions on smooth regions, which otherwise were only obtained after more selection steps. For some discontinuous problems, our numerical results indicated that the hybrid ENO-flux limiter schemes performed better and were computationally quicker to run as compared to some higher-order ENO schemes. When the limiter MM1 was used with the hybrid scheme, we got comparable results to ENO2. Applying the MM2 limiter produced sharper results on shocks, but it generated overshooting in the profiles of hyperbolic systems because of the steep approximations of the slopes. The hybrid schemes with UNO and Harmod limiters achieved the best resolved approximations, produced sharp resolutions of shocks and reduced the oscillations which may be present in ENO3. For hyperbolic systems, componentwise extensions were adopted instead of characteristic decompositions. We would like to mention that according to existing literature, for example [10], characteristic decompositions may reduce oscillations but at the expenses of more computations.