Hybrid BSQI-WENO Based Numerical Scheme for Hyperbolic Conservation Laws
In this paper, we intend to use a B-spline quasi-interpolation (BSQI) technique to develop higher order hybrid schemes for conservation laws. As a first step, we develop cubic and quintic B-spline quasi-interpolation based numerical methods for hyperbolic conservation laws in 1 space dimension, and show that they achieve the rate of convergence 4 and 6, respectively. Although the BSQI schemes that we develop are shown to be stable, they produce spurious oscillations in the vicinity of shocks, as expected. In order to suppress the oscillations, we conjugate the BSQI schemes with the fifth order weighted essentially non-oscillatory (WENO5) scheme. We use a weak local truncation based estimate to detect the high gradient regions of the numerical solution. We use this information to capture shocks using WENO scheme, whereas the BSQI based scheme is used in the smooth regions. For the time discretization, we consider a strong stability preserving (SSP) Runge-Kutta method of order three. At the end, we demonstrate the accuracy and the efficiency of the proposed schemes over the WENO5 scheme through numerical experiments.
Keywords:Weighted Essentially Non-Oscillatory Scheme B-spline Quasi-interpolation
In this work, we intend to use B-Spline Quasi-Interpolation (BSQI) to develop numerical schemes for the hyperbolic conservation laws
where and is the flux function, along with the initial condition
It is well-known that the classical solution of (1.1) may ceases to exist in finite time, even the initial data is sufficiently smooth. The higher order finite difference, spectral, and finite element methods are well suitable in approximating the solution in the smooth regions, but fail to provide a non-oscillatory solution in the presence of jump discontinuities. Thus, it is challenging to devise an higher order numerical scheme, which preserves the higher order accuracy in smooth regions and produce non-oscillatory solution in a vicinity of shocks.
The BSQI approximations are known (see Sablonnière ) to be very accurate in approximating smooth functions and their derivatives. This motivated many researchers to use this technique to develop numerical methods for some partial differential equations. For instance, Li et al.  developed a mesh-free method for the viscous Burgers’ equation, and Zhu and Kang  developed a finite difference method for the Burgers-Fisher equation. Recently, we (Kumar and Baskar ) have developed higher order schemes for some Sobolev type equations with special reference to equal width and BBMB equations. But these equations involve dissipative and dispersive effects that lead to smeared shocks and dispersive shocks, respectively. Zhu and Kang  made an initial effort in devising cubic B-spline quasi-interpolation (CBSQI) based numerical scheme for hyperbolic conservation laws. But their study is restricted only to scalar conservation laws and also they used quasi-linear form of the equation to develop the numerical scheme and presented numerical results for Burgers’ equation up to the critical time at which the shock forms. To our knowledge, there is no work on developing a BSQI based conservative numerical method for hyperbolic conservation laws.
In the present work, our interest is to develop higher order methods for (1.1)-(1.2) using cubic and quintic BSQI (denoted by CBSQI and QnBSQI, respectively) for approximating space derivative and the 3 order Runge-Kutta method for time derivative. Although the resulting schemes are conditionally stable, as they are higher order methods, they develop spurious oscillations near the shocks. It is known that such a problem is common for any higher order methods. There are at least two approaches to minimize the non-physical oscillations in the numerical solution, namely, using limiters and using essentially non-oscillatory (ENO) reconstructions. Out of these, the ENO approach can minimize the oscillations without much compromise in the order of accuracy.
An improved version of the ENO approach is the well known Weighted ENO (WENO) reconstruction introduced by Liu et al.  (also see Jiang and Shu ). However, this reconstruction algorithm is expensive. In order to minimize the expense, there are many hybrid schemes introduced in the literature for the hyperbolic conservation laws. For instance, a non-conservative hybrid method is introduced by Harabetian and Pego , the multi-domain methods in one and two space dimensions are developed using Fourier continuation and WENO scheme by Shahbazi et al. [32, 33]. Costa and Don [4, 5] obtained a high order hybrid method by combining the spectral and the central schemes with WENO. Cheng and Liu [3, 2] introduced multi-domain RKDG methods with proper boundary treatments. The characteristics hybrid compact-WENO scheme is proposed by Ren et al. .
In developing efficient hybrid schemes, choosing smooth indicator is an important component of the algorithm. Karni et al. , used a Weak Local Truncation Error (WLTE) based smooth indicator in developing an adaptive algorithm for hyperbolic systems. Further Kurganov and Liu  developed adaptive artificial viscosity methods for the system of Euler equations of motion using Weak Local Residual (WLR) (similar to WLTE). The WLR is used in many other problems, for instance, Mungkasi et al.  used WLR for the system of shallow water equations and Chen et al.  developed the adaptive artificial viscosity method for the Saint-Venant system using this indicator. To reduce the computational cost in the case of system of Euler equations, Dewar et al.  introduced WLR pressure based smooth indicator. A comparison study of many switch indicators for hyperbolic conservation laws is performed by Li and Qiu [21, 22] for uniform and curvilinear meshes.
In this paper, we propose higher order hybrid BSQI-WENO schemes for hyperbolic conservation laws. In a BSQI-WENO scheme, we conjugate the BSQI based numerical scheme with a WENO scheme. To detect the high gradient region of the numerical solution, we use a WLR based estimate of Kurganov and Liu . We use this information to capture shocks using WENO scheme, whereas the BSQI based scheme is used in the smooth region. For the time discretization, we consider a SSP Runge-Kutta method of order three. To validate the numerical scheme, we perform the numerical experiments using Burgers and Buckley-Leverett equations in the case of scalar conservation laws and Euler equations in the case of system of equations. It is also interesting and important to study the efficiency and the accuracy of the hybrid methods constructed using the BSQI based numerical schemes and WENO schemes.
The outline of the paper is as follows. In section 2, we discuss the BSQI and its application in approximating the derivative of a function. In section 3, we develop the semi-discrete form of hyperbolic conservation laws using CBSQI and QnBSQI approximations to the space derivative. We further perform the von Neumann stability analysis and provide the numerical evidence for the order of accuracy. We develop two hybrid numerical schemes, namely, the Hybrid4 and Hybrid6 by combining the CBSQI with WENO3 and QnBSQI with WENO5, respectively, in section 4. We perform some numerical experiments in section 5, to validate the schemes for scalar as well as system of equations. Finally, we test the accuracy and the efficiency of the proposed algorithms over the pure WENO in section 6.
2 B-Spline Quasi-Interpolation
In this section, we briefly discuss the construction of BSQI of degree upto 5 and give the numerical rate of convergence in each cases. For more details on the basic concepts of spline spaces and B-splines we refer to Schumaker , De boor , Nürnberger , and Kvasov .
Let the space interval be divide into subintervals , of equal length . The boundary points of the subintervals are called the knots and let
denotes the set of knots. The spline space of degree is defined by (see Schumaker )
where is the polynomial space of degree . The set forms a vector space of dimension and the set of B-splines of degree forms a basis for this space. B-splines can be constructed recursively by starting with the B-spline of degree 0 given by
The recursive formula for the higher degree B-splines takes the form
The general form of the B-spline quasi-interpolant (BSQI) of degree for a function is given by (see Sablonnière )
where the coefficients can be obtained by imposing the condition that the quasi-interpolant is exact if is a polynomial of degree less than or equal to . The coefficients are the linear combination of function (approximate) values at the nodes, where
if is odd, then the nodes are simply the knots, i.e. , and the set of nodes involved in calculating is
if is even, then the nodes are taken as and the set of nodes involved in calculating is
For more details, we refer to Sablonnière .
In case when support of B-spline
contains distinct knots, the coefficients for quadratic (QBSQI), cubic (CBSQI), quadric (QdBSQI), and quintic (QnBSQI) are given as follows:
We have the following theorem for general BSQI in approximating a function:
There exists a constant such that for all and for any set of knots on with step size , we have
From the construction we can see that a B-splines of degree is times continuously differentiable. Thus, we can also approximate the derivatives of a function with the corresponding derivatives of BSQI as long as it is possible, i.e.
The approximation of function and its first derivative at point are as follows:
Case I. When is even, we have
Case II. when is odd, we have
In the following example, we demonstrate the accuracy and the rate of convergence in approximating the first derivative of a smooth function.
Consider the function
whose first derivative is given by In order to make the comparison among the BSQI’s in terms of their accuracy, we consider the -error norm given by
where the superscript denotes the derivative and denotes the error associated with the derivative. In Figure 1(a), we depict the log-log plot of the -error and the number of mesh points, in the case of the first derivative approximation of the model problem (2.16). The QBSQI converges to the first derivative with rate two, whereas CBSQI and QdBSQI achieves the fourth order of convergence, which is clear from Figure 1(b). The QnBSQI has a higher order accuracy among the considered BSQI, which has the rate of convergence six in the case of the first derivative approximation.
3 BSQI Scheme
In this section, we discuss the semi-discrete form of (1.1) obtained using the BSQI based spatial approximation. For the sake of notation simplicity, we restrict to , and consider the scalar conservation laws of the form
The BSQI formulation for a general 1D system (1.1) can be done componentwise in a similar way.
The above methods can work for any degree of BSQI. But in this work, we intend to use only the CBSQI and QnBSQI formulations, which are given, respectively, by
for , where we have used the notation . The above semi-discrete formulations can be written in the conservation form as
where denotes the numerical flux given by
for the CBSQI scheme (3.21) and
for the QnBSQI scheme (3.22).
where and with being the wave-number, we get
where is given by
Note that the spectrum of the Runge-Kutta method of order three is known to contain some part of the imaginary axis. Thus, this method along with the CBSQI and QnBSQI space discretizations can be stable under the condition
where and (see Hirsch  for more details). This implies that, we can find the condition on such that the scaled spectrum of CBSQI and QnBSQI schemes are contained in the stability spectrum of the Runge-Kutta method of order three.
In the following example, we demonstrate the order of accuracy of the CBSQI and QnBSQI schemes in the case of linear advection equation.
Consider the linear advection equation
with the initial data
The rate of convergence of a numerical scheme, denoted by , for the error , is given by
where and denote the error and the number of mesh points, respectively. From Example 1, we expect the rate of convergence, with respect to the space variable, of the CBSQI and the QnBSQI schemes to be four and six, respectively. In order to achieve the rate numerically, we take the time step . The numerical solution is computed over the domain at time . Table 2 and Table 3, clearly shows that, as the number of mesh points increase, the numerical solution obtained with CBSQI and QnBSQI schemes converge to the exact solution with rate four and sixth, respectively.
We now show numerically that the CBSQI and QnBSQI schemes achieve the expected order of accuracy in the case of Burgers equation as long as the solution remains smooth.
Consider the Burgers’ equation
In order to show accuracy of the BSQI schemes, we consider (3.30) with the following initial data
In Table 4 and Table 5, we compare the -, -, and -error and their rate of convergence obtained using the CBSQI and the QnBSQI schemes, respectively. The numerical solution is computed over the at time , at which the solution remains smooth. In order to achieve the expected rate of convergence, we choose the time step significantly small to make time discretization error small. In both cases, we choose the time step . From the tables, it is clear that the CBSQI scheme converges to the exact solution with a rate more then four, whereas the QnBSQI scheme initial has order nearly four which increased further to order more then six. ∎
4 Hybrid Scheme
The CBSQI and the QnBSQI schemes developed in the previous section are stable with 3 order Runge-Kutta time discretization and also they achieve the expected order of accuracy when the solution remains smooth, as demonstrated in Example 2 and 3. But the main disadvantage of these schemes is that the numerical solution develops spurious oscillations in the case of non-smooth solution. In order to minimize the oscillations without much compromise in the order of accuracy, we devise an hybrid scheme, where we use the BSQI scheme in the region where the solution is smooth, whereas in the non-smooth region we propose to use WENO scheme. Thus, we propose the hybrid semi-discrete scheme for (3.18) to be
with being a smooth indicator, which is a constant over the cell i.e . Using this in (4.33), the above scheme can be written as
Here, denotes either CBSQI or QnBSQI flux given by (3.24) and (3.25), respectively. Whereas denotes the WENO flux of appropriate order. More precisely, we conjugate WENO flux of order 3 with CBSQI and form a hybrid scheme with numerical flux
and we combine QnBSQI with WENO flux of order 5 to get
To be self-content, we quickly recall the WENO fluxes and for more details, we refer to Jiang and Shu .
The flux can be decomposed into positive and negative parts as
such that and . For instance, these can be taken as
where . We recall that the stencil of WENO3 is and , and the corresponding numerical flux is written as
where the positive part is given by
and and . The negative part can be written in a similar way.
The stencil of WENO5 is
and the positive part of the numerical flux is written as
Here the smooth indicators are given by
and , and .
Thus, the hybrid scheme switches between WENO and BSQI scheme, depending upon the smooth indicator. Once the shock region is detected, the hybrid scheme switches to the expensive shock capturing WENO scheme, otherwise the solution is computed using the BSQI scheme.
For the smooth indicator, we use the weak local truncation error (WLTE) introduced by Karni and Kurganov . The WLTE is derived on the fact that a weak solution of (3.18) satisfies the integral form given by (see Godlewski and Raviart  )
for all . The variation in the value of from smooth to discontinuous regions and vice versa for computed solution can be taken as a measure of smooth indicator. The is referred as weak truncation error for with respect to . In practice, the calculation of weak truncation error is seem to be a difficult task, since is a general test function. To overcome this difficulty Kurganov and his co-workers (see [17, 14]) used test function with B-splines as
Karni and Kurganov  established an estimate for the WLTE (4.43) and remarked that the WLTE bound can be converted into error bound in the smooth region. Although their result is valid for one dimension scalar problem, they showed numerical that the WLTE behavior is similar in the case of 1D system of equations. With the aid of these estimates and further numerical experiments, Kurganov and Liu  proposed the following estimate and used it as a smooth indicator:
where , and is the accuracy of the numerical scheme. In our hybrid scheme (4.34), we propose to use the smooth indicator defined by
where is some empirically chosen nonnegative real number. In general, we take in our numerical experiments. To make smooth transition between smooth and rough regions, we extend the rough region to the neighboring node points of . That is, once the discontinuity is detected in a vicinity of , which is the case if , then we take
where is again an empirically chosen positive integer. In our numerical computation, we take . The smooth transition can also be made using a smaller value of , but that lead to over estimate of discontinuous regions.
5 Numerical Experiments
This section comprises the numerical implementation of hybrid BSQI-WENO schemes for linear and non-linear problems. The BSQI schemes are oscillatory in nature at shock positions, whereas the hybrid BSQI-WENO schemes developed in the previous section are non-oscillatory because of the usage of the WENO scheme in the non-smooth regions and the high order accuracy is maintained in the smooth regions using BSQI approximation. The developed algorithm works for any degree of BSQI, but to demonstrate the idea we use CBSQI and QnBSQI for space approximation in conjugation with WENO3 and WENO5, respectively. We study the accuracy of the CBSQI-WENO3 and the QnBSQI-WENO5 schemes in terms of rate of convergence. Along with the accuracy, we also demonstrate the efficiency of the proposed hybrid BSQI-WENO.
We start the numerical experiment with smooth initial data.
Example 4 (Smooth Initial Data)
In order to measure the accuracy and convergence rate of the hybrid algorithm, we consider the Example 3.1 and Example 3.2. In both cases, we take in smooth indicator. We fixed the time step and increase the mesh points in space. In Figure 2(a), we compare the accuracy of the Hybrid4 and the Hybrid6 schemes with the WENO5 scheme in case of linear advection equation at time . Both the Hybrid4 and the Hybrid6 schemes achieve their expected convergence rate. In terms of accuracy, we observe that the Hybrid6 scheme is better than both the WENO5 and the Hybrid4 schemes. Similar behavior is observed from Figure 2(b) in the case of Burgers’ equation at time .
As a next step, we perform the numerical experiment in the case of the linear advection equation with discontinuous initial data.
Example 5 (Linear advection equation)
Consider the linear advection equation (3.27) with initial data
which has jump discontinuities at points and moves towards right as time progress. The CBSQI and QnBSQI schemes capture the correct position of the discontinuities but with oscillations in the profile (not shown here). We now compute the solution using the Hybrid6 scheme whose numerical flux is given by (4.36). In Figure 3(a), we compare the numerical solution with the exact solution at time along with the initial data. The number of mesh points is taken to be and the CFL number is . In the figure, the symbol ‘’ represents the region where the WENO5 is used whereas the symbol ‘’ represents the region where QnCBSQI is used. We observe that the discontinuities are captured without oscillations. In Figure 3(b), the plot of smooth indicator over the plane is shown. It is clear from the figure that the smooth indicator moves along the characteristic direction and detect the discontinuity positions correctly. Also, we observe that the solution is computed with WENO5 only at discontinuity positions, whereas in the smooth regions QnBSQI is used to compute the solution.∎
We now test the Hybrid6 scheme in the case of some nonlinear scalar conservation laws. We first start with the inviscid Burgers equation, which is a common test case for convex flux.
Example 6 (Nonlinear equation with convex flux - Burgers equation)
To show the shock capturing ability of the Hybrid6 scheme (4.36), we consider the Burgers’ equation (3.30) with discontinuous initial data (5.46). The numerical solution is computed at time over the domain with mesh points. The CFL number is taken to be . The solution contains a rarefaction wave and a shock moving towards right. As expected the QnBSQI scheme is oscillatory at the shock position, which spreads over the smooth region of the solution as the time progresses, which can be observed from Figure 4(a).
Figure 4(b) depicts that the numerical solution at time obtained using the Hybrid6 scheme (4.36). Here, we observe that the scheme captures shock and rarefaction accurately without any oscillation. The smooth indicator over the -plane is shown in Figure 4(c), where we observe that the indicator moves along the discontinuity. The left state of the initial discontinuity at the point is connected with the right state by the rarefaction fan, where the singularity at the lower part of the rarefaction wave is smoothen after time due to numerical dissipation, which is well illustrated by the smooth indicator.
In Figure 4(d), we show the plot of the percentage of the computational domain in which WENO5 is used verses the number of grid points. We kept the time fixed at and increased the number of mesh points and observed that the percentage of the region where WENO5 scheme is used has been reduced drastically from to when the number of mesh points is increased from to . The low percentage of usage of the WENO scheme for large number of mesh points, makes the hybrid BSQI-WENO schemes efficient.∎
We now turn our attention to (3.18) with non-convex fluxs. Since vanishes at one or more points, (3.18) fails to be genuinely nonlinear which makes the solution more complicated as it may involve the formation of the composite waves. The composite waves are generally the combination of joined rarefaction and shock waves, for more detail see LeVeque . Many popular schemes fail to approximate the entropy solution of (3.18) with non-convex flux, for instance semi-discrete central-upwind schemes (see Kurganov et al. ). Here we test the Hybrid6 scheme (4.36) in case of the Buckley-Leverett problem.
Example 7 (Non-convex flux - Buckley-Leverett equation)
Consider the conservation law (3.18) with non-convex flux, called the Buckley-Leverett equation,
along with the initial data (5.46). In Figure 5(a), the numerical solution obtained using QnBSQI scheme at time is compared with the reference solution over the domain with and CFL is taken to be . In this figure, we observe that the QnBSQI scheme fails to produce non-oscillatory result. The rarefaction fan is missing and the numerical solution is oscillatory at the rarefaction position, which makes the solution to moves with different speed. As a result, the QnBSQI scheme fails to capture the correct positions of the shocks. In comparison to QnBSQI, the Hybrid6 scheme (4.36) produces non-oscillatory solution with correct positions of shocks as is apparent from Figure 5(b) and thus converges to the entropy solution. In Figure 5(c), we show the smooth indicator on -plane and the percentage of the usage of WENO5 scheme is depicted in Figure 5(d). As observed in the case of Burgers equation, here also we see that the usage of WENO5 is approximately when the number of mesh points approaches . This shows that, the Hybrid6 scheme captures the entropy solution with less computational complexity. ∎
We continue with the scalar conservation law (3.18) with another non-convex flux.
Example 8 (Another non-convex flux)
Another test case for non-convex flux is (see for instance Kurganov et al. )
which belongs to , whose convexity changes at . We consider the following initial data
The exact solution is given by