Journal of Applied Mathematics

Volume 2013 (2013), Article ID 428681, 12 pages

http://dx.doi.org/10.1155/2013/428681

Research Article

Optimized Weighted Essentially Nonoscillatory Third-Order Schemes for Hyperbolic Conservation Laws

A. R. Appadu¹ and A. A. I. Peer²

¹Department of Mathematics and Applied Mathematics, University of Pretoria, Pretoria 0002, South Africa

²Department of Applied Mathematical Sciences, University of Technology, Mauritius, La Tour Koenig, Pointe-aux-Sables, Mauritius

Received 12 April 2013; Revised 25 June 2013; Accepted 25 June 2013

Academic Editor: Mohamed Fathy El-Amin

Copyright © 2013 A. R. Appadu and A. A. I. Peer. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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.

1. Introduction

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 [1]. 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 [1]. 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 [2]. All linear numerical schemes are either dispersive or dissipative [3]. 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 [4]. 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 [5]. 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 [6]. However, they generate a great deal of computational noise in highly steep gradients. An iteration smoothing process was proposed by Forester [7] 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 [8] and van Leer [9]. The approach of flux-corrected transport is an early version of the concept of limiters [10]. Boris and Book [8] introduced a flux-corrected method which improves the accuracy of a variety of finite difference algorithms. The second version of limiters is geometric limiters [10]. One example is van Leer’s scheme which blends Fromm’s method with a geometric limiter [9]. The scheme naturally detects a discontinuity and changes its behaviour accordingly.

Essentially Nonoscillatory (ENO) and WENO schemes have been developed to capture shocks efficiently [11]. 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 [5]. 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 [12]. 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 [4].

Kasahara and Washington [13] 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 [14] have combined the Leap-Frog scheme with the Trapezoidal Rule. Liska and Wendroff [4] 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

The smoothness indicators (13) of the nonlinear weights (11) and (12) turn out to be

with linear coefficients and .

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.

Figure 1: Plot of modulus of amplification factor versus phase angle versus cfl number for the WENO3 scheme.

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 .

(a)

(b)

Figure 2: Plots of AFM versus phase angle at some cfl numbers for the WENO3 scheme.

(a)

(b)

Figure 3: Plots of RPE versus phase angle at some cfl numbers for the WENO3 scheme.

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 [4];(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 [15], 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 .

The eeldld denotes exponential error for low dispersion and low dissipation [15, 16].

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 [15], 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 [17], we use the technique to understand why not all composite methods can be effective to control dispersion and dissipation in regions of shocks. In [16], we consider the family of third-order methods proposed by Takacs [2] 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.

Figure 4: Plot of IEELDLD versus cfl number for the WENO3 scheme.

6. Optimisation of WENO3 Scheme Discretizing the 1D Burgers’ Equation

Ascher and McLachlan [18] have obtained the dispersion relation for the equation

They considered the linearized version of (43) which is

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.

7. Quantification of Errors from Numerical Results [2, 15, 16]

In this section, we describe how Takacs [2] 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.

Table 1: Errors for Problem I at time and .

error
order

40

—

80

160

Table 2: Errors for Problem I at time and .

error
order

40

—

80

160

Table 3: Errors for Problem I at time and .

error
order

40

—

80

160

Table 4: Errors for Problem I at time and .

error
order

40

—

80

160

Table 5: Errors for Problem I at time and .

error
order

40

—

80

160

Table 6: Errors for Problem I at time and .

error
order

40

—

80

160

Table 7: Errors for Problem I at time and .

error
order

40

—

80

2.96

160

2.70

8.2. Problem II

We solve with the initial condition being for and elsewhere [19]. 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: Errors for Problem II at time .

cfl
Dissipation error
Dispersion error
Total error
eeldld

0.05

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

1.1

1.2

Figure 5: Results of Problem II at cfl number 0.05 with number of spatial steps, .

Figure 6: Results of Problem II at cfl number 1.1 with number of spatial steps, .

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

The time . We perform the experiment using WENO3 over 200 cells at some cfl numbers. The errors, dispersion errors, and dissipation errors are shown in Table 9. Results are shown in Figures 7 and 8.

Table 9: Errors for Problem III at time .

cfl
error
Dissipation error
Dispersion error
Total error
eeldld

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

1.1

1.2

Figure 7: Solutions for Problem III with 200 cells at using WENO3 scheme at cfl number 0.50.

Figure 8: Solutions for Problem III with 200 cells at using WENO3 scheme at cfl number 1.3.

We observe that the errors are least at cfl 0.50 and greatest at cfl 1.2. Figure 9 shows the plots of error, and total error and eeldld, all versus the cfl number. These three errors are least at low cfl numbers.

(a)

(b)

(c)

Figure 9: Plots of errors versus cfl number using WENO3 scheme for Problem III.

9. Conclusion

In this work, we study the spectral analysis of the dissipation and dispersion errors of the WENO3 scheme at some different cfl numbers and verify numerically the order of convergence of the WENO3 scheme. It is observed that low cfl numbers are preferred to minimise the error, dispersion and dissipation errors, and eeldld. An extension of this work is to obtain the optimal cfl of some WENO schemes discretizing some 1D nonlinear equations such as Korteweg-de-Vries and also 2D equations such as scalar advection, convection-diffusion, and Korteweg-de-Vries equations.

Nomenclature

: Time step

: Spatial step

: Time level

: Advection velocity

: / number

: Phase angle in 1D

: Relative phase error per unit time step

: Amplification factor

.

Acknowledgments

The authors thank the anonymous referees whose comments helped to improve the work and the presentation of this paper. This work was initiated when Dr. Appadu was on a research visit at the University of Technology which was funded by the Research Development Programme (RDP) of the University of Pretoria. The period of support for the RDP funds is from January 2012 to January 2014.

References

D. A. Randall, “An introduction to atmospheric modeling,” Tech. Rep., 2007.

L. L. Takacs, “A two-step scheme for the advection equation with minimized dissipation and dispersion errors,” Monthly Weather Review, vol. 113, no. 6, pp. 1050–1065, 1985.

H. Trac and U.-L. Pen, “A primer on eulerian computational fluid dynamics for astrophysics,” Publications of the Astronomical Society of the Pacific, vol. 115, no. 805, pp. 303–321, 2003.

R. Liska and B. Wendroff, “Composite schemes for conservation laws,” SIAM Journal on Numerical Analysis, vol. 35, no. 6, pp. 2250–2271, 1998.

Z. Jiang, “On dispersion-controlled principles for non-oscillatory shock-capturing schemes,” Acta Mechanica Sinica, vol. 20, no. 1, pp. 1–15, 2004.

S. Okamoto, K. Hara, A. Takei, and F. Masuda, “A study on numerical methods for air quality simulation,” Tokyo University of Information Sciences, vol. 5, no. 1, pp. 65–72, 2001.

C. K. Forester, “Higher order monotonic convective difference schemes,” Journal of Computational Physics, vol. 23, no. 1, pp. 1–22, 1977.

J. P. Boris and D. L. Book, “Flux-corrected transport. I. SHASTA, a fluid transport algorithm that works,” Journal of Computational Physics, vol. 11, no. 1, pp. 38–69, 1973.

B. van Leer, “Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme,” Journal of Computational Physics, vol. 14, no. 4, pp. 361–370, 1974.

P. Colella and E. Puckett, “Modern methods for fluid flow,” Tech. Rep., 1994.

X.-D. Liu, S. Osher, and T. Chan, “Weighted essentially non-oscillatory schemes,” Journal of Computational Physics, vol. 115, no. 1, pp. 200–212, 1994.

J. E. Fromm, “A method for reducing dispersion in convective difference schemes,” Journal of Computational Physics, vol. 3, no. 2, pp. 176–189, 1968.

A. Kasahara and W. M. Washington, “Thermal and dynamical effects of orography on the general circulation of the atmosphere,” in Proceedings of the WMO/IUGG Symbosium on Numerical Weather Prediction, pp. 47–56, Tokyo, Japan, November-December 1968.

A. F. Shchepetkin and J. C. Mcwilliams, “Quasi-monotone advection schemes based on explicit locally adaptive dissipation,” Monthly Weather Review, vol. 126, no. 6, pp. 1541–1580, 1998.

A. R. Appadu and M. Z. Dauhoo, “The concept of minimized integrated exponential error for low dispersion and low dissipation schemes,” International Journal for Numerical Methods in Fluids, vol. 65, no. 5, pp. 578–601, 2011.

A. R. Appadu, “Some applications of the concept of minimized integrated exponential error for low dispersion and low dissipation,” International Journal for Numerical Methods in Fluids, vol. 68, no. 2, pp. 244–268, 2012.

A. R. Appadu, “Investigating the shock-capturing properties of some composite numerical schemes for the 1-D linear advection equation,” International Journal of Computer Applications in Technology, vol. 43, no. 2, pp. 79–92, 2012.

U. M. Ascher and R. I. McLachlan, “Multisymplectic box schemes and the Korteweg-de Vries equation,” Applied Numerical Mathematics, vol. 48, no. 3-4, pp. 255–269, 2004.

A. A. I. Peer, M. Z. Dauhoo, A. Gopaul, and M. Bhuruth, “A weighted ENO-flux limiter scheme for hyperbolic conservation laws,” International Journal of Computer Mathematics, vol. 87, no. 15, pp. 3467–3488, 2010.