Mathematical and Numerical Modeling of Flow and Transport 2013View this Special Issue
Research Article | Open Access
A. R. Appadu, A. A. I. Peer, "Optimized Weighted Essentially Nonoscillatory Third-Order Schemes for Hyperbolic Conservation Laws", Journal of Applied Mathematics, vol. 2013, Article ID 428681, 12 pages, 2013. https://doi.org/10.1155/2013/428681
Optimized Weighted Essentially Nonoscillatory Third-Order Schemes for Hyperbolic Conservation Laws
We describe briefly how a third-order Weighted Essentially Nonoscillatory (WENO) scheme is derived by coupling a WENO spatial discretization scheme with a temporal integration scheme. The scheme is termed WENO3. We perform a spectral analysis of its dispersive and dissipative properties when used to approximate the 1D linear advection equation and use a technique of optimisation to find the optimal cfl number of the scheme. We carry out some numerical experiments dealing with wave propagation based on the 1D linear advection and 1D Burger’s equation at some different cfl numbers and show that the optimal cfl does indeed cause less dispersion, less dissipation, and lower errors. Lastly, we test numerically the order of convergence of the WENO3 scheme.
Spatial discretization methods in order to solve partial differential equations can be broadly classified as finite difference, finite volume, finite element, and spectral methods. All these methods combined with explicit or implicit time integration schemes can be effectively applied to solve partial differential equations.
Computational dispersion arises from space differencing and must not be confused with instability . A noisy solution is not necessarily unstable. Both dispersion and instability can lead to noise. In the case of dispersion, the waves do not grow in amplitude but become separated from one another, each at its own speed . Computational dispersion causes waves of different wavelengths to travel at different speeds.
One of the most common ways of measuring the relative merit of a numerical scheme for advection is to consider the scheme’s dispersion and dissipation . All linear numerical schemes are either dispersive or dissipative . We now briefly survey some techniques implemented to minimise dispersion and/or dissipation, in regions of shocks.
Antidiffusion has been proposed in order to improve the resolving power of dissipative schemes . Antidiffusion does well with shock-type discontinuities, but at the same time, distortions are produced in smooth regions. Artificial viscosity does work to suppress oscillations and capture shock waves . However, it ruins the resolution of numerical schemes in recognising shock waves. Filters are also used to control computational noise. The Galerkin finite element methods are useful to solve the advection equation . However, they generate a great deal of computational noise in highly steep gradients. An iteration smoothing process was proposed by Forester  in order to smooth the computational values from the results of the Galerkin finite element method.
The concept of limiters was introduced by Boris and Book  and van Leer . The approach of flux-corrected transport is an early version of the concept of limiters . Boris and Book  introduced a flux-corrected method which improves the accuracy of a variety of finite difference algorithms. The second version of limiters is geometric limiters . One example is van Leer’s scheme which blends Fromm’s method with a geometric limiter . The scheme naturally detects a discontinuity and changes its behaviour accordingly.
Essentially Nonoscillatory (ENO) and WENO schemes have been developed to capture shocks efficiently . These are high-order accurate finite difference schemes designed for problems with piecewise smooth solutions containing discontinuities. The fundamental notion resides at the approximation level, where a nonlinear adaptive procedure is used to choose the locally smoothest stencil, thus avoiding crossing discontinuities in the interpolation procedure as much as possible. WENO methods suffer from a tendency for the weights of the various stencils to revert quite slowly to the underlying high-order method as the numerical resolution of a smooth function improves.
The effects of dispersion terms on numerical solutions have been done since a long time . However, it is only since the last two decades that extensive studies on the dispersion-controlled dissipative scheme were reported. Work on scheme dispersion reduction was first reported by Fromm . Fromm’s scheme in 1D is made up of a linear combination of Lax-Wendroff and the Beam-Warming second-order upwind schemes as these two dispersive schemes have phase errors in opposite directions. We can also combine dispersive and dissipative schemes to obtain composite schemes. This idea was implemented in meteorological codes which compose the oscillatory second-order Leap-Frog scheme with the dissipative backward Euler scheme .
Kasahara and Washington  use a three-level Leap-Frog scheme for 50 time steps which is inherently unstable followed by one cycle of the Lax-Wendroff scheme. The Leap-Frog scheme when applied for 50 time steps is unstable. The Lax-Wendroff scheme is essentially dispersive but does possess some damping property. The damping properties of Lax-Wendroff render the combined procedure of the application of 50 time steps of the three-level Leap-Frog scheme and one cycle of Lax-Wendroff scheme stable. Shchepetkin and McWilliams  have combined the Leap-Frog scheme with the Trapezoidal Rule. Liska and Wendroff  have obtained a composite scheme which comprises Lax-Wendroff as the dispersive scheme and the two-step Lax-Friedrichs scheme as the dissipative scheme to obtain the Lax-Wendroff/Lax-Friedrichs scheme.
The paper is organised as follows. In Sections 2 and 3, we construct the WENO reconstruction and the time integration schemes. The originality in this work resides in the fact that expressions for the amplification factor and the relative phase error have been obtained for the WENO3 scheme discretizing the 1D linear advection equation. Very few papers have examined these properties for WENO schemes as it is cumbersome to obtain the expression for the amplification factor of WENO schemes in general. Spectral analysis of the WENO3 scheme is discussed in Section 4. In Section 5, we compute the optimal cfl number of the WENO3 scheme by using a technique of optimisation that minimizes both the dispersion and dissipation errors. Such an optimal cfl number is important as the cfl number is chosen quite arbitrarily for WENO schemes, and, hence, the efficient shock-capturing properties of WENO methods are not exploited. Section 7 gives a description of how to quantify errors from numerical results into dispersion and dissipation. Some numerical experiments are performed in Section 8 and errors are compared at different cfl numbers. Section 9 highlights the salient features of the paper.
2. WENO Spatial Reconstructions
WENO schemes are commonly used for the solution of hyperbolic conservation laws of the form where is the conserved variable and is the flux function.
Integrating (1) over the cell and denoting its cell average by , we get the equation
Upwind schemes sample the solution at the cell centres , and in this case we can write (2) in the semidiscrete form as follows: where the exact flux is . Since the flux is required on cell boundaries, the discontinuities which may be present prevent the use of a quadrature formula to integrate in time. A Riemann solver is needed to identify the direction of propagation of the waves, and then an ODE solver is applied to advance (3) in time.
The relation (3) requires a reconstruction to recover the approximation of point values . Assuming that the cell averages are known, the following piecewise polynomial reconstruction is done: where is the characteristic function of the cell . Denoting the order of accuracy by , is a polynomial of degree at most on the cell , such that and maintains conservation; that is,
The approximations at the cell boundaries for the semidiscrete upwind scheme (3) are given by and . The evolution step consists of approximating the point values from and , by identifying the waves according to their direction of propagation, and to find the numerical flux, . Here we consider the monotone Lax-Friedrichs flux which is given by
A WENO reconstruction consists of associating a weight to all of the candidate stencils of an Essentially Nonoscillatory (ENO) reconstruction. All the possible approximations at cell boundaries for are evaluated using , and then where the weights satisfy If is smooth over all the candidate stencils, then some linear weights are used, such that The linear weights for are found by carrying out a Taylor series expansion of (10).
In case there is a discontinuity present in one of the candidate stencils, its corresponding weight turns to be nearly 0, such that only a smooth approximation is considered. The weights are given by where The parameter in (12) is used to avoid divisions by zeros during computations, and as per the majority of literature in the field we will take as .
The parameters are called the smoothness indicators of the candidate stencils. They make the weights of smooth stencils dominant over those of discontinuous stencils. The smoothness indicators are measures of the -norms of the derivatives of ; that is, If a stencil is smooth then its corresponding smoothness indicator , but if the stencil contains a discontinuity then . This strategy reduces the contribution of nonsmooth stencils.
Similarly, the WENO principle is used for , where for , and the terms satisfy .
3. Time Integration Schemes
In this section, we describe some time integration schemes to advance the solution of (2) in time. The time step is taken as where the CFL number is necessary for stability.
Let be an approximation to the spatial derivative and consider the semidiscrete formulation (3) as a system of initial value problem of ODE
We mention a result that motivates the use of high-order Strong Stability-Preserving (SSP) Runge-Kutta (RK) methods in time. If the first-order forward Euler time stepping method is SSP, that is, for , then the same strong stability is preserved for higher order SSPRK methods under the restriction where is the CFL coefficient for the high-order time integration.
The explicit stage Runge-Kutta method has the general form For nonnegative coefficients and , the Runge-Kutta method (20) is a convex combination of Euler forward schemes. The scheme should satisfy .
Lemma 1. If the forward Euler method (17) is strongly stable under the CFL restriction , , then the Runge-Kutta method (20) with is SSP provided the following CFL restriction (19) is fulfilled:
For the special class of optimal stage th-order SSPRK schemes, the CFL coefficient is maximized according to Lemma 1. Optimal second-order and third-order SSPRK schemes are, respectively, given by
4. Spectral Analysis of WENO3
In this paper, we focus on third-order WENO3 reconstruction. Setting in (8), we get where
For the sake of finding amplification factors of numerical schemes, we consider the linear advection case into (1), where is a constant. This results in the numerical flux (7) simplifying to and, therefore, (16) is approximated by
The fully discretized method is given by where .
On substituting by where is the amplification factor, , and is the phase angle, we obtain the amplification factor as where .
The real and imaginary parts of are given by respectively.
The modulus of the amplification factor, AFM, is computed as and the region of stability is obtained as as depicted in Figure 1.
The relative phase error, , is obtained as
Plots of the AFM versus the phase angle at six different cfl numbers, namely, 0.25, 0.5, 0.75, 1.0, 1.25, and 1.47, are shown in Figure 2. The scheme is least dissipative at . At low cfl numbers, the scheme is less dissipative especially for phase angles, . Plots of the RPE versus phase angle at the six different cfl numbers are shown in Figure 3. There is no big change in the dispersive character at cfl 0.25, 0.5, 0.75, and 1.0, and at these cfl numbers, phase lag behaviour is exhibited at all phase angles, . At phase angles greater than 0.5, phase lead behaviour is exhibited at .
The modified equation of WENO3 when discretized by the 1D linear advection equation is given by indicating that the scheme is third-order accurate in space and essentially dissipative in nature.
5. Optimisation of WENO3 Scheme Discretizing the 1D Linear Advection Equation
The technique of Minimized Integrated Exponential Error for Low Dispersion and Low Dissipation (MIEELDLD) has been introduced in Appadu and Dauhoo (2011). It basically enables us to choose the optimal parameters from two conditions; namely,(i)small amounts of dissipation when added can help to curb numerical dispersion ;(ii)the dissipation and dispersion errors must both be small in a numerical scheme to yield efficient shock-capturing properties.We now describe the technique of minimized integrated exponential error for low dispersion and low dissipation (MIEELDLD). Suppose that the amplification factor of the numerical scheme when applied to the 1D linear advection equation, given by is The quantities and measure the dispersion and dissipation errors, respectively. For a numerical scheme to have low dispersion and low dissipation, we require Also when dissipation neutralises dispersion optimally, we have Thus on combining these two conditions, we get the following condition necessary for dissipation to neutralise dispersion and for low dispersion and low dissipation character to be satisfied: where eldld denotes error for low dispersion and low dissipation.
If we plot the quantity, eldld versus RPE versus AFM, as done in , we can see that when (no dispersion) and (no dissipation), and this makes sense. However, the eldld takes a constant value of 2 when independent of the value of the AFM, and this presents a drawback of the measure. Therefore, we have presented a modification to the quantity, eldld, which is and this quantity goes to zero when and .
5.1. Only One Parameter Involved
If the number is the only parameter, we compute the integrated error for low dispersion and low dissipation, IELDLD, for a range of , and this integral will be a function of . The optimal is the one at which the integral quantity is closest to zero.
5.2. Two Parameters Are Involved
Suppose that we now have two parameters, say and . In that case, we can compute and this integral will be a function of and we can obtain the optimal value of .
We can also compute and this integral will consist of and . From there, we can obtain the optimal values of both and .
Considerable and extensive work on the technique of minimised integrated exponential error for low dispersion and low dissipation has been carried out in [15–17]. In , we have obtained the optimal cfl of some explicit methods like Lax-Wendroff, Beam-Warming, and Upwind Leap-Frog when applied to the 1D advection equation. In , we use the technique to understand why not all composite methods can be effective to control dispersion and dissipation in regions of shocks. In , we consider the family of third-order methods proposed by Takacs  where we optimize two parameters, namely, the cfl number and another variable which also controls dispersion and dissipation.
In this work, we use the technique, MIEELDLD, to compute the optimal cfl number of the WENO3 scheme. The plots of the integrated error versus cfl number are shown in Figure 4, which indicates that the optimal cfl is close to zero.
6. Optimisation of WENO3 Scheme Discretizing the 1D Burgers’ Equation
By considering a plane wave solution of the form , where and are the wavenumber and dispersion relation, respectively, the exact dispersion relation is
The exact phase velocity can then be deduced. A numerical scheme can then be used to discretize (44), and, thus, the numerical phase velocity can be obtained. The relative phase error is then computed as the ratio of the numerical phase velocity to the exact phase velocity.
Since a linearised version of the 1D Burgers’ equation, is not possible, we cannot obtain the relative phase error, and thus it is not possible to obtain the variation of the quantity, IEELDLD, versus to deduce the optimal cfl number. To estimate the optimal cfl number, we can run our numerical experiment at some cfl numbers, compute the error, dispersion, dissipation errors, and eeldld, and then guess the optimal cfl number. This approach has been done in Section 8.3.
In this section, we describe how Takacs  quantifies errors from numerical results into dispersion and dissipation errors.
The Total Mean Square Error is calculated as where represents the analytical solution and is the numerical (discrete) solution at a given grid point, .
The Total Mean Square Error can be expressed as Next, The Total Mean Square Error can be further expressed as
The expression in (50) can be rewritten as where and denote the variance of and , respectively, and and , denote the mean values of and respectively.
Thus, the Total Mean Square Error is given by which on further simplification yields Thus, we have But the correlation coefficient, , is given by . Hence, the Total Mean Square Error can be written as which simplifies to On putting , we get . Thus, we define as the dispersion error as correlation coefficient in statistics is analogous with phase lag or phase lead in Computational Fluid Dynamics.
Consequently, measures the dissipation error. The total error is defined as the sum of the dissipation and dispersion errors.
Since our domain , the error is calculated as where and are the computed and exact values, respectively, and and are the step length and number of cells, respectively.
8. Numerical Experiments
8.1. Problem I
We solve with the initial conditions for . We display the errors and convergence rates at time at some different cfl numbers, namely, 0.05, 0.1, 0.25, 0.5, 1.0, 1.25, and in Tables 1, 2, 3, 4, 5, 6, and 7. The order of convergence obtained numerically at the seven values of cfl is approximately equal to 3. At low cfl numbers, the order is in general closer to 3. This shows that the choice of the cfl number affects to a certain extent the order of convergence.
8.2. Problem II
We solve with the initial condition being for and elsewhere . We use some different values of cfl and cells, . We use . The various types of errors, namely, dissipation error and dispersion error, are tabulated for some values of cfl numbers in Table 8. It can be seen clearly that the cfl number influences the shock-capturing property of the scheme to a great extent. Results are depicted in Figures 5 and 6.
Table 8 shows the dissipation, dispersion, total errors, and eeldld at some values of cfl. We observe that the total error and eeldld are least at low cfl numbers. These two quantities increase monotonically as the cfl number increases.
8.3. Problem III
We solve Burger’ equation with initial conditions being