Efficient WENO Schemes for Nonuniform Grids

1.1 Scope

Non-linear interpolation methods as essentially non-oscillatory (ENO) and weighted ENO (WENO) schemes were developed to approximate numerical solutions of partial differential equations (PDEs), see e.g. [13, 16, 20,21,22, 26], although they have also been used successfully in non-PDE applications as image processing (see e.g. [2, 4]), computer vision (see e.g. [23] mentioned by [20]) or finance (see e.g. [14]). Recently, a vast quantity of new WENO variants have been proposed to enhance the accuracy of the traditional WENO scheme proposed in [21, 22] as, for instance, Mapped WENO (WENO-M) by Henrick et al. [11], WENO-Z [8] or a WENO scheme with progressive order of accuracy close to discontinuities [1].

The idea behind the design of the WENO algorithm is quite powerful. It is based on the construction of a piecewise interpolator trying to avoid the points where the function is discontinuous. For this, a convex combination of interpolations of lower degree is proposed, using non-linear weights for the contribution of each interpolator. The final result is an interpolator of high-order accuracy if the data used are free of discontinuities and a lower order of accuracy interpolator if the utilized stencil crosses a singularity. The main ingredients of WENO algorithms are the optimal weights, which are the values needed to obtain the highest order of accuracy, and the smoothness indicators which indicate if a subset of points contains a discontinuity. There are certain situations where the optimal weights are negative and non-linear WENO weights do not satisfy the required conditions (see [20]). Some solutions for this issue are summarized in [20]. Levy et al. proposed in [15] a central scheme based on a centered version of WENO, named CWENO, which overcomes the mentioned matter.

On the one hand, in [5], the authors introduced a central WENO scheme based on a global average weight independent of the optimal weights. This method has two key ingredients: first, the global average weight, which is defined using the smoothness indicators through elementary operations, and second, the evaluation of the interpolator polynomial using the whole stencil. This method was developed for structured grids. However, in some applications (see [10, 18], mentioned in [20]), WENO reconstructions on non-uniform meshes are needed.

On the other hand, in [9], the authors presented an alternative family of smoothness indicators that are simpler and computationally cheaper than those originally proposed by Jiang and Shu in [13]. They proved that the computational time needed to obtain the proposed smoothness indicators increased linearly with respect to the accuracy order of the WENO scheme considered while for Jiang and Shu’s smoothness indicators the growth is quadratic.

In this paper, we combine and extend the ideas proposed in [5, 9] for nonuniform grids and prove theoretically that the new proposed scheme satisfies the desired accuracy. In addition, we show that this scheme efficiently solves, in some cases, the loss of order of accuracy near smooth extrema associated with the original WENO and CWENO schemes. We will present some mimetic examples to check these theoretical results. Another key feature of this method is that the cost of the computation of the smoothness indicators for a nonuniform stencil does not have a significant impact on the whole computational cost of the scheme, compared with the cost that the computation of the smoothness indicators has on the original WENO and CWENO schemes. We will apply the new WENO reconstruction method in finite volume schemes to obtain numerical solutions in shock problems from hyperbolic conservation laws to illustrate its performance.

1.2 Outline of This Paper

The remainder of the paper is divided as follows: in Sect. 2 we briefly review the classical Lagrange interpolation and the notation that will be used throughout this paper. We present linear reconstructions from point values and cell averages on nonuniform grids. The new WENO method will be introduced in Sect. 3 and some theoretical results will be provided. In particular, we will prove that the new WENO scheme satisfies the necessary properties to achieve maximum order of accuracy when the data do not cross any singularity and the resulting interpolator is essentially non-oscillatory. Section 4 is devoted to showing the computational aspects of the algorithm. In Sect. 5, we present some numerical experiments to validate the theoretical results and to test the performance of the new scheme in shock problems from hyperbolic conservation laws. Finally, some conclusions and remarks will be drawn in Sect. 6.

We start introducing basic definitions which will be used throughout this work. We also review the definition of linear pointwise reconstructions on nonuniform grids from given point values and cell averages.

Definition 2.1

In the remainder of the paper, we denote by \(\Pi _k\), \(k \in \mathbb {N}_0\), the space of polynomials of maximal degree k, i.e.

and by \(\bar{\Pi }_k\) the space of polynomials of exact degree k, \(k \in \mathbb {N}_0\), i.e.

Definition 2.2

For \(\alpha \in \mathbb {Z}\):

Since, for positive f, g,

it follows that

2.1 Reconstructions on Nonuniform Grids

Since we want to design WENO schemes working on nonuniform grids, the first step is to describe the procedure to obtain a reconstruction from given data at the desired order. We will focus both on the pointwise reconstructions from point values (the traditional Lagrange interpolation/extrapolation) and the reconstructions from cell averages.

2.1.1 Reconstructions from point values

Let \(S=\{x_0,\ldots ,x_{R-1}\}\) be a stencil, \(x_i<x_j\), \(i<j\), and \(f_i:=f(x_i)\), \(0\le i\le R-1\) for some function f. It is well-known that the polynomial \(p\in \Pi _{R-1}\) satisfying \(p(x_i)=f_i\) is given by the Lagrange form:

2.1.2 Reconstructions from cell averages

Now, let us consider a stencil consisting of cell interfaces given by \(S=\{x_0,\ldots ,x_R\}\), \(x_i<x_j\), \(i<j\), being the associated cells \(C_i=[x_i,x_{i+1}]\), \(0\le i\le R-1\) and \(\bar{f}_i\) the known values given by

Then, let us assume that we want to find \(p\in \Pi _{R-1}\) such that p(x) interpolates f(x). To do so, we impose

To find p, we use the strategy proposed in [3] for structured grids, which we next generalize to the nonuniform case. Let us consider

with \(p_i\in \Pi _{R-1}\) to be determined, \(0\le i\le R-1\). Then, imposing (2.1) on (2.2) and using linearity, we have for a fixed j,

Therefore, (2.3) is satisfied if we impose

with \(\delta _{i,j}\) the Kronecker delta.

Let \(P_i\in \Pi _R\) be a primitive of \(p_i\), namely, such that \(P_i'=p_i\). Then, for each \(0\le i\le R-1\), (2.4) reads

Hence, it suffices to impose

to write \(P_i\) in Lagrange form as

Thus, taking derivatives at both sides of (2.5) we get

and (2.2) then reads

The classical Lagrange methods introduced in Sect. 2.1 are useful to interpolate smooth data but some artifacts appear when they are used in applications with non-continuous inputs. To solve this problem, we present in this section a generalization of the WENO method for nonuniform grids.

3.1 Construction of the Method and Theoretical Results

Let \(h>0\), \(R\ge 3\) and let \(S_{R,h}=\{x_{0,h},\ldots ,x_{R-1,h}\}\) be a stencil, with \(x_{i,h}=z+c_ih\), for some fixed \(z\in \mathbb {R}\) and \(c_i\in \mathbb {R}\), \(0\le i\le R-1\), \(c_i<c_j\), \(i<j\). Let f be a function and define \(f_{i,h}=f(x_{i,h})\).

We denote by \(x_h^*\in \mathbb {R}\) the point in which we desire to reconstruct the data from the stencil \(S_{R,h}\). This point is located between \(x_{r-1,h}\) and \(x_{r,h}\) if R is even, with \(r=R/2\), or rather between \(x_{r-1,h}\) and \(x_{r,h}\) (left-biased) or between \(x_{r,h}\) and \(x_{r+1,h}\) (right-biased) if R is odd, with \(r=(R-1)/2\). Note that this is actually a generalization of the traditional centered or left/right-biased WENO/CWENO schemes on cartesian grids to nonuniform grids, respectively, (see e.g. [15, 16]).

Let us first consider \(p_{R,i,h}\) the reconstruction polynomials satisfying \(p_{R,i,h}(x_{j,h})=\mathcal {L}[f](x_{j,h})\), where the operator \(\mathcal {L}\) is defined depending on the framework chosen (point values or cell average reconstructions), i.e.,

with \(i\le j\le i+r\), for \(0\le i\le r'\), with \(r=\lfloor (R-1)/2\rfloor \) and \(r'=\lceil (R-1)/2\rceil \), and evaluate them at \(x_h^*\), where \(\lfloor \cdot \rfloor \) denotes the floor function and \(\lceil \cdot \rceil \) the ceil function.

Following the same idea based on the Yamaleev-Carpenter weight construction proposed in [25], we define \(d_{R,h}\) as

where \(f[x_{0,h},\ldots ,x_{R-1,h}]\) is the undivided difference defined as

The simplified smoothness indicators, akin to the idea proposed in [9], which in turn is the key to the simplicity of this algorithm, are defined by

Since in this case, in the context of an arbitrary grid, there is no point in computing ideal weights in order to build the reconstruction with optimal accuracy from the lower order reconstruction, we simply consider a uniform convex combination based on the \(r'+1\) values \(\frac{1}{r'+1}\) and define

with s an exponent to be determined (will be done so when the accuracy analysis is performed) and \(\varepsilon >0\) a small positive quantity in order to avoid divisions by zero. With this weight design, we can define the following weighted essentially non-oscillatory reconstruction:

Since this reconstruction does not attain the desired optimal order, R, we use the reconstruction associated with the whole stencil \(S_{R,h}\), \(p_{R,h}\) satisfying \(p_{R,h}(x_{i,h})=\mathcal {L}[f](x_{i,h})\), which we evaluate at \(x_h^*\), and the following global average weight, given by

Then, the final reconstruction with essentially non-oscillatory properties is given by

Now that we have introduced and defined all the components of the WENO method, we are in a position to prove the principal properties of the scheme in terms of the order of the weights and accuracy of the reconstructions obtained.

Proposition 3.1

For the standard case, and keeping the same notation, if f has a critical point of order k at z, \(0\le k\le R-2\), then it holds \(0\le \omega _{R,i,h}\le 1\), \(0\le \omega _{R,h}\le 1\) and

with \(m=\max \limits _{0\le i<R-1}m_i\), where

Proof

If \(f\in \mathcal {C}^R\), then \(d_{R,h}=\mathcal {O}(h^{2R-2})\). On the other hand, by all the considerations derived from [6, Section 2] and [9, Section 2] we conclude that, since \(r=1\) if \(3\le R\le 4\) and \(r>1\) if \(R\ge 5\), it holds \(I_{R,i,h}=\overline{\mathcal {O}}(h^{2m_i})\), with

Therefore, we have

Hence,

with \(m=\max \limits _{0\le i<R-1}m_i\).

Now, let us assume that a discontinuity crosses \(S_{R,h}\). Then, we define as A the set of indexes j such that the discontinuity does not cross \(S_{R,{j},h}\) and as B the set of indexes j such that the discontinuity crosses \(S_{R,{j},h}\). Note that, by construction of \(r'\), the stencils do not overlap, and thus \(A\ne \emptyset \). Then, on the one hand, if \(j\in A\):

On the other hand, if \(j\in B\):

Consequently, if \(i\in A\), we have:

with \(m':=\max \nolimits _{j\in A}m_j\), which satisfies \(m'\ge m_i\), since \(i\in A\) by assumption. However, if \(i\in B\), then:

Now we have to analyze \(\omega _{R,h}\). To do so, we first study the accuracy of \(J_{R,h,s}\):

Therefore, if \(f\in \mathcal {C}^R\), then \(d_{R,h}=\mathcal {O}(h^{2R-2})\) and

where the last equality holds since \(m'\le m\). Finally, if a discontinuity crosses \(S_{R,h}\), then \(d_{R,h}=\overline{\mathcal {O}}(1)\) and

\(\square \)

Theorem 3.1

Keeping the same notation, if \(s\ge (r+1)/2\), z is not a critical point of f for \(3\le R\le 4\) or z is not a critical point of f of order \(R-2\) for \(R\ge 5\), then the optimal accuracy is attained both in the case of smoothness and the case of discontinuities, namely:

Proof

Let us first assume that \(f\in \mathcal {C}^R\). Then, by Proposition 3.1:

with \(m=\max _{0\le i<R-1}m_i\), where

Now, taking into account our hypothesis involving the critical points, we have that \(m_i=k+1\), \(\forall \, 0\le i\le R-1\), \(\forall \, R\ge 3\). Moreover, \(k=0\) for \(3\le R\le 4\) and \(k\ne R-2\) for \(R\ge 5\), and thus we can assume \(k<R-2\) for any \(R\ge 3\) (the case \(k>R-2\) being trivial, as we will point out at the end of the proof). Therefore, taking into account that

and using that \(\sum \nolimits _{i=0}^{r'}\omega _i=1\), we have

Using (3.2) and (3.3):

where in the last equality we have used that \(R\ge r+1\). Therefore, if one wants the optimal R-th order accuracy, one must impose

On the other hand, since

it suffices to impose

Let us consider now the case where a discontinuity crosses \(S_{R,h}\), then, again by Proposition 3.1, we have

and now it holds

Therefore,

Combining in this case (3.4) and (3.5), we can conclude that

Then, in order to attain the optimal accuracy, one must impose

and therefore, it suffices with the choice

to attain optimal accuracy both in the smooth and the discontinuous case. \(\square \)

3.2 Remarks Regarding the Computational Cost

As stated in the introduction, one of the strongest points of the proposed WENO method is the ease of computation of the smoothness indicators, even taking into consideration the fact that it is defined in the context of a nonuniform grid.

Unlike the case of the traditional Jiang–Shu smoothness indicators, which in turn are used in the method proposed in [5], the indicators used in this work, based on those proposed in [9], are considerably cheaper to compute, especially when the grid is non-uniform.

To be more specific, given \(R\in \mathbb {N}\), a stencil of the form \(S_{R,h}=\{x_{0,h},\ldots ,x_{R-1,h}\}\), with \(x_{i,h}=z+c_ih\), \(c_i<c_j\), the set of nodal values \(\{f_{0,h},\ldots ,f_{R-1,h}\}\), with \(f_{i,h}=f(x_{i,h})\), for f a function to be approximated, and the reconstruction point \(x_h^*=z+c^*h\), the simplified expression for the classical Jiang–Shu smoothness indicators has the structure

with \(\alpha _{i,j,k}\) to be determined for each \(0\le i\le j\le R-1\), \(0\le k\le R-1\).

From (3.6), it can be clearly seen that the number of operations that are required to compute the Jiang–Shu smoothness indicators in terms of the accuracy order is \(\mathcal {O}(R^3)\). Therefore, any WENO method using these smoothness indicators has, at least, a cubic cost with respect to the aforementioned order of the method. This includes the techniques proposed by Levy, Puppo and Russo in [15] or those proposed in [5].

In contrast, from (3.1) the smoothness indicators proposed in this paper have linear cost with respect to the accuracy order, namely, \(\mathcal {O}(R)\). It is also important to mention that, depending on the coefficients \(c_i\) and the location of the reconstruction point \(x_h^*\), the ideal weights in the method proposed by Levy, Puppo and Russo can be negative. This issue was fixed in [5] and the method proposed here neither suffers from this phenomenon.

Since computing an interpolating polynomial under a non-uniform grid requires evaluating an expression of the form

with \(\beta _{i,k}\) to be determined, \(0\le i,k\le R-1\), the cost from that part of the algorithm requires a quadratic number of operations with respect to the accuracy order, namely, \(\mathcal {O}(R^2)\). Thus, taking into account that \(\lfloor \frac{R}{2}\rfloor \) polynomials and smoothness indicators have to be computed in the process, the computational cost of the WENO methods using the Jiang–Shu smoothness indicators have \(\mathcal {O}(R^4)\) cost when they are applied on nonuniform grids, while the WENO method proposed in this paper has \(\mathcal {O}(R^3)\) cost with respect to the accuracy order in a context of nonuniform grids.

In the particular case in which a static nonuniform grid is used during the whole time evolution, the terms \(\alpha _{i,j,k}\) and \(\beta _{i,k}\) can be precomputed for each local stencil and stored to be used on each time iteration. This implies that the computational cost with respect to the order of the expression in (3.6) is reduced to quadratic, while the cost in (3.7) is reduced to linear. Therefore, the global computational cost of the CWENO methods [5, 15] can be reduced to \(\mathcal {O}(R^3)\), while the cost associated with the method presented in this paper is reduced to \(\mathcal {O}(R^2)\).

In this Section, we summarize the algorithm of the proposed WENO method, including the considerations involving the items to attain unconditionally the optimal accuracy. We divide this section into two parts: First, we introduce the description for the reconstruction from point values and after we present the necessary changes to adapt it for reconstructions from cell average values.

4.1 Algorithm for Reconstructions from Point Values

Input data:

• Grid of nodes: \(S_{R,h}=\{x_{0,h},\ldots ,x_{R-1,h}\}\), with \(x_{i,h}=z+c_ih\), \(c_i<c_j\), \(i<j\), \(0\le i,j\le R-1\). It suffices to provide the \(c_i\) values, \(0\le i\le R-1\).

• Nodal values: \(\{f_{0,h},\ldots ,f_{R-1,h}\}\), with \(f_{i,h}=f(x_{i,h})\).

• Reconstruction point: \(x_h^*\), satisfying \(x_{R/2-1,h}\le x_h^*\le x_{R/2,h}\) if R is even or \(x_{(R-1)/2-1,h}\le x_h^*\le x_{(R-1)/2+1,h}\) if R is odd. It suffices to provide \(c^*\), with \(x_h^*=z+c^*h\).

• Small positive quantity: \(\varepsilon >0\).

Then, the next steps are to be followed:

1. (1)

   Compute

   $$\begin{aligned} r&=\left\lfloor \frac{R-1}{2}\right\rfloor ,\\ r'&=\Bigg \lceil \frac{R-1}{2}\Bigg \rceil ,\\ s&=\Bigg \lceil \frac{r+1}{2}\Bigg \rceil . \end{aligned}$$
2. (2)

   Compute the reconstruction polynomials at \(c_h^*\) of the corresponding substencils \(S_{R,i,h}\), \(0\le i\le r'\), \(p_{R,i,h}\), and the one associated to the global stencil, \(S_{R,h}\), \(p_{R,h}\), which are:

   $$\begin{aligned} p_{R,i,h}(c_h^*)&=\sum _{j=i}^{i+r}\left[ \prod _{l=i,l\ne j}^{i+r}\frac{(c_h^*-c_l)}{(c_j-c_l)}\right] f_j,\\ p_{R,h}(c_h^*)&=\sum _{j=0}^{R-1}\left[ \prod _{l=0,l\ne j}^{R-1}\frac{(c_h^*-c_l)}{(c_j-c_l)}\right] f_j. \end{aligned}$$
3. (3)

   Compute the smoothness indicators:

   $$\begin{aligned} I_{R,i,h}:=\sum _{j=i}^{r+i-1}\left( \frac{f_{j+1,h}-f_{j,h}}{c_{j+1}-c_j}\right) ^2,\quad 0\le i\le r'. \end{aligned}$$

   Optionally, the denominators can be omitted in both expressions without any impact on the accuracy of the scheme.

4. (4)

   Obtain \(d_{R,h}\):

   $$\begin{aligned} d_{R,h}=\left( (R-1)!\sum _{i=0}^{R-1}\frac{1}{\prod _{j=0,j\ne i}^{R-1}(c_i-c_j)}f_i\right) ^2. \end{aligned}$$
5. (5)

   Compute \(J_{R,h,s}\):

   $$\begin{aligned} J_{R,h,s}=\sum _{i=0}^{r'}\frac{1}{I_{R,i,h}^s+\varepsilon }. \end{aligned}$$
6. (6)

   Generate the weights:

   $$\begin{aligned} \omega _{R,i,h}&=\frac{\alpha _{R,i,h}}{\sum \nolimits _{j=0}^{r'}\alpha _{R,j,h}},\text { with }\alpha _{R,i,h}=\frac{1}{r'+1}\left( 1+\frac{d_{R,h}^s}{I_{R,i,h}^s+\varepsilon }\right) ,\\ \omega _{R,h}&=\frac{1}{1+d_{R,h}^sJ_{R,h,s}}. \end{aligned}$$
7. (7)

   Obtain the half-optimal-order essentially non-oscillatory reconstruction:

   $$\begin{aligned} {\tilde{q}}_{R,h}(x_h^*)=\sum _{i=0}^{r'}\omega _ip_{R,i,h}(x_h^*). \end{aligned}$$
8. (8)

   Compute the optimal order essentially non-oscillatory reconstruction:

   $$\begin{aligned} q_{R,h}(x_h^*)=\omega _{R,h}p_{R,h}(x_h^*)+(1-\omega _{R,h}){\tilde{q}}_{R,h}(x_h^*). \end{aligned}$$

4.2 Algorithm for Reconstructions from Cell Averages Values

Input data:

• Grid of cells: \(S_{R,h}=\{C_{0,h},\ldots ,C_{R-1,h}\}\), with \(C_i=[x_{i,h},x_{i+1,h}]\), \(0\le i\le R-1\), \(x_{i,h}=z+c_ih\), \(c_i<c_j\), \(i<j\), \(0\le i,j\le R\). It suffices to provide the \(c_i\) values, \(0\le i\le R\).

• Cell averages: \(\{\bar{f}_{0,h},\ldots ,\bar{f}_{R-1,h}\}\), with

  $$\begin{aligned} \bar{f}_{i,h}=\frac{1}{x_{i+1,h}-x_{i,h}}\int _{x_{i,h}}^{x_{i+1,h}}f(x)dx. \end{aligned}$$

• Reconstruction point: \(x_h^*\), satisfying \(x_{(R+1)/2-1,h}\le x_h^*\le x_{(R+1)/2,h}\) if R is odd or \(x_{R/2-1,h}\le x_h^*\le x_{R/2+1,h}\) if R is even. It suffices to provide \(c^*\), with \(x_h^*=z+c^*h\).

• Small positive quantity: \(\varepsilon >0\).

In step (2) reconstruction polynomials should be defined by:

In step (3) the smoothness indicators are:

and in step (4):

Remark 4.1

Note that, if one works with fixed grids, then all the terms involving computations related to the grid spacings, namely, \(c_i\), as well as the relative position of the reconstructions points \(c_h^*\), can be computed and stored only once, and then used to perform the computations, thus reducing to linear cost with respect to the order the computation of terms such as \(p_{R,i,h}\), \(p_{R,h}\) and \(d_{R,h}\), yielding a comparable performance with respect to the classical WENO schemes with structured grids.

In this section, we perform some numerical examples to check the accuracy properties of our schemes. We divide the experiments into two sections: Firstly, we apply our algorithm in mimetic algebraic numerical tests to assess the order of accuracy. Secondly, we numerically solve several PDEs. In both cases, we implement the method to non-continuous data.

In all the experiments the C++ wrapper [12] of the GNU MPFR library [24] is used with a precision of 332 bits (\(\approx 100\) digits) and taking \(\varepsilon =10^{-10^5}\).

5.1 Algebraic Numerical Experiments

The goal of this subsection is to check the behaviour of the proposed generalized WENO scheme in different scenarios. We present two different tests: one with a smooth function, and one for non-smooth functions computed both for point values and for cell average values for which the estimated order of accuracy is calculated.

In all the numerical experiments, we consider a mesh of the form \(x_{i,h}=z+c_ih_n\), \(0\le i<N\), with \(h_{n+1}=\frac{h_n}{2}\) and provide the table of the resulting experiments based on taking \(h_0=0.2\) and \(0\le n\le 20\). The order values are thus computed as

with \(x_{h_n}^*\) the reconstruction node and \(u_{h_n}^*\) the corresponding resulting WENO reconstruction.

Representation of the different nonuniform meshes considered in Tests 1 and 2 for \(h=1\), \(z=0\) and the values of \(c*\) indicated in each case. Note that to simplify the notation we have used \(x_i\) to denote \(x_{i,h}\). a Mesh for Test 1 and smooth initial data, b mesh for Test 1 and discontinuous initial data, c mesh for Test 2 and smooth initial data and d mesh for Test 2 and discontinuous initial data

Full size image

5.1.1 Test 1: Reconstruction of a Smooth Function

Let us consider the function \(f(x):=xe^x\) and the grid given in Fig. 1a, defined by the values \(c_i\), showed in Table 1 (column 1), with \(N=12\) and \(z=0\). We consider the reconstruction point \(x_h^*=z+c^*h\), with \(c^*=0\), hence in this case \(x_h^*\) is always located between the nodes \(x_{5,h}\) and \(x_{6,h}\).

First, we interpret the given data as point values, namely, we consider \(f_{i,h}=f(x_{i,h})\), \(0\le i\le 11\). In this case, the expected accuracy order should be the optimal twelfth-order accuracy that would attain the reconstruction procedure without WENO weights. By applying it to our WENO scheme, Table 2 shows the numerical errors and accuracy order obtained, which clearly shows the convergence to the aforementioned optimal order, as we have proved in Theorem 3.1.

Now, we interpret the data as cell averages, namely, we consider

using the same function f and the grid given in Fig. 1b defined by the values \(c_i\), in Table 1 (column 2). We consider \(x_h^*=z+ch\), with \(c^*=0\). As can be seen in Table 2, the behaviour of the numerical solutions is the same as before: the expected order of accuracy is reached.

5.1.2 Test 2: Reconstruction of a Discontinuous Function

Let us now consider now the function

and the grid given by the data shown in column 3 of Table 1 with \(c^*=2.3251\). Note that the reconstruction point is now located between \(x_{4,h}\) and \(x_{5,h}\). In this case, the discontinuity divides the stencil into two halves, where the first three nodes contain the information corresponding to the left side of the discontinuity, namely, \(xe^x\), and the last eight nodes contain the information related to the right side of the discontinuity, that is, \(2xe^x+1\). The expected behaviour is therefore an accuracy drop to sixth order (since \(\lfloor \frac{11+1}{2}\rfloor =6\)), which is confirmed in Table 3.

We repeat the experiment for cell averages with nodes given by the values \(c_i\), \(0\le i\le 11\) in the fourth column of Table 1 and \(c^*=0.5041\). In this case, the reconstruction point is located in the cell with bounds \(x_{5,h}\) and \(x_{6,h}\) as it can be seen in Fig. 1d. As before, the discontinuity separates the cells into two parts. The first one is given by the first five cells, and the second one contains six cells with the information related to the right side of the discontinuity. As we expected, in Table 3 we observe that the accuracy order reaches sixth order, as we have studied in Theorem 3.1.

5.2 Numerical Experiments Involving Conservation Laws

To illustrate the applicability of our proposed WENO interpolator in the context of numerical solutions of PDEs we next consider some experiments using one-dimensional hyperbolic conservation laws \(u_t+f(u)_x=0.\) First, we will test our method using some academic examples such as the linear advection equation, with \(f(u)=u\), and Burgers equation, which corresponds with \(f(u)=u^2/2\). Later, we will perform two more sophisticated experiments to compare the results obtained with the non-uniform approach and the uniform one in more challenging scenarios.

For this purpose, a finite volume scheme will be used (see e.g. [19]), which is suitable to attain a method with high order accuracy. Given a domain \(\Omega =[a,b]\) and a partition of \(\Omega \), \(\Big \{x_{i+\frac{1}{2}}\Big \}_{i=0}^N\), the semidiscrete form of the conservation law reads

where

To approximate (5.2) for approximations \(\overline{u}_{i}(t)\approx \overline{u}_t(x_i,t)\), we use the conservative difference formula:

where the numerical fluxes, \(\overline{f}_{i+\frac{1}{2}}\approx f(u(x_{i+\frac{1}{2}},t))\), are approximated using:

For the autonomous linear advection equation:

Left-biased upwind reconstructions are applied:

$$\begin{aligned} \overline{f}_{i+\frac{1}{2}}=q^-_{5,h}\left( x_{i+\frac{1}{2}}\right) , \end{aligned}$$

where \(q^-_{5,h}\) is the WENO reconstruction using five cells left-biased, i.e., \(\overline{u}_{i+l}\), \(-3\le l\le 2\).

For the rest of the numerical experiments in this work:

Local Lax-Friedrichs flux splitting is utilized:

$$\begin{aligned} \overline{f}_{i+\frac{1}{2}}=\frac{1}{2}\left( f\left( q^+_{5,h}\left( x_{i+\frac{1}{2}}\right) \right) +f\left( q^-_{5,h}\left( x_{i+\frac{1}{2}}\right) \right) -\alpha \left( q^+_{5,h}\left( x_{i+\frac{1}{2}}\right) -q^-_{5,h}\left( x_{i+\frac{1}{2}}\right) \right) \right) , \end{aligned}$$

where \(q^{\pm }_{5,h}\) are the WENO reconstructions using five cells left and right-biased and \(\alpha =\max _{u} |f'(u)|\).

The time-stepping solver that will be used along all the experiments is the second- or third-order Runge–Kutta TVD scheme [21]. The spatial discretization will always have order 5.

5.2.1 Construction of the Nonuniform Grid

For the sake of reproducibility, we provide all the details about the non-uniform grid construction in each experiment. Given \(n\in \mathbb {N}\), we define

with \(R(0,\xi )=R(n,\xi )=0\) and

for \(1\le j\le n-1\). The parameter \(\xi \) controls the fluctuation of the cell interfaces with respect to a uniform grid, and the function \(r:\mathbb {N}^3\rightarrow [0,1]\) is the Wichmann-Hill pseudo-random number generator given by

with seed update

In all the experiments, the selected starting seeds will be chosen as

In the experiments analyzing the order under smooth solutions, namely, those involving several simulations in which the number of cells is duplicated, the seeds are not restarted upon the next doubled resolution; instead, the last updated seed from the previous experiment with halved resolution is used as starting seed for the next one.

5.2.2 Test 3: Linear Advection Equation

First of all, we compute the approximations of the solution of the linear advection equation:

with initial and boundary conditions:

which exact solution is given by:

As the method used to solve the ODE is third-order accurate, for test purposes we choose \(\Delta t=(\widetilde{\Delta x})^{5/3}\) with \(\widetilde{\Delta x}=\min _{i}\{\Delta x_i\}\). With this selection, we guarantee that \(\Delta t^3=O(\Delta x^{5})\), with \({\Delta x}=\max _{i}\{\Delta x_i\}\) and \(\Delta t/\Delta x \le 1\). In this case, the fluctuation parameter is \(\xi =0.1\). We calculate the error for \(T=1\) and the numerical order in a similar way as we computed in Equation (5.1) for values of \(n=20\cdot 2^j\) with . The results are displayed in Table 4. It is clear that the expected order of accuracy is reached.

We consider now the discontinuous initial data given by:

\(\Delta t=0.9\widetilde{\Delta x}\) and the fluctuation parameter \(\xi =0.25\). The simulation is now run until \(T=1.5\) and with \(n=40\) and \(n=100\). We can see in Fig. 2 that the behaviour of the numerical solution is oscillation-free and does not present Gibbs phenomenon close to the discontinuities present in the solution.

Test 3: Numerical solution of the linear advection equation with discontinuous initial data for a \(n=40\), b \(n=100\)

Full size image

5.2.3 Test 4: Burgers Equation

If we repeat the same experiments as in Test 3 but for the Burgers equation with \(T=0.3\) for smooth initial data and \(T=1\) for non-smooth initial data, we obtain the results shown in Table 5 and in Fig. 3. Again, the expected approximation order, 5, is achieved for both norm 1 and norm infinity. The simulations in Fig. 3 show that for small values of n some numerical oscillations appear near the boundaries. However, this problem is overcome when we refine the mesh.

Test 4: Numerical solution of the Burgers equation with discontinuous initial data: a \(n=20\), b for \(n=100\)

Full size image

5.2.4 Test 5: Scalar Equation with Solution Approaching a Dirac Delta Function

It is the goal of this test to show a numerical experiment where the results obtained using the nonuniform approach significantly outperform those of the uniform one. We consider the equation (see [17])

with initial condition \(u(x,0)=1\), \(x\in [0,2\pi ]\), and periodic boundary conditions. By the method of characteristics, the exact solution of (5.3) is given by

For a fixed value of \(x\in [0,2\pi ]\), there holds

and for fixed \(t\ge 0\)

therefore, the solution of this equation behaves progressively like a Dirac delta function as time increases. This means that the solution will have a strong gradient near \(x=\pi \), while it would remain close to zero in the rest of the domain. The structure of the solution suggests that the use of a nonuniform grid with decreasing distance between nodes as we approach the point \(x=\pi \) could be convenient.

We set \(T=8\) and consider uniform grids with \(n=100 \cdot 2^{j}-1, j=1,\dots ,10\) cells. We define non-uniform grids of \(2\,m\) cells with cell boundaries defined by

and, symmetrically, \(x_{m+j+\frac{1}{2}}=2\pi -x_{m-j+\frac{1}{2}}\), \(j=1,\dots ,m\), with \(\kappa >1\).

The errors are computed as

where \(u_{j}^{N}\) is the result of N time steps of length \(\Delta t=T/N\), subject to the following CFL restriction

The ODE solver used in this test is TVD-RK2.

In Fig. 4 we show the numerical solutions obtained using both uniform and non-uniform meshes. Figure  4c and d show that the accuracy of the approximated solutions obtained using non-uniform grids is higher than the ones obtained using uniform grids, even when a significantly lower resolution is considered.

Table 6 displays the errors and orders obtained when using uniform grids. For comparison, the error obtained using a non-uniform grid with \(2m=198\) cells and \(q=1.1\) is \(2.87e-2\), which is comparable to the error obtained with a uniform grid of 25,599 cells. If we consider \(2m=398\) cells and \(q=1.04\) the computed error is \(2.67e-3\), which in this case is comparable to the error obtained with a uniform grid with 51,199 cells.

It is worth mentioning that although the number of cells on the non-uniform grids considered is much smaller than the number of cells on the uniform ones, the smallest cell size for a non-uniform grid with \(2m=198\) cells and \(q=1.1\) is \(2.51e-5\) while for the grid with \(2\,m=398\) and \(q=1.04\) is \(5.13e-5\). Those sizes are smaller than the cell size for the uniform grids with 25,599 and 51,199 cells, given by \(2.45e-4\) and \(1.23e-4\), respectively.

Test 5: Numerical solution of the scalar equation converging to a Dirac delta function in the interval [3.138, 3.145], obtained using: a uniform meshes with \(n=25{,}599, 51{,}199, 102{,}399\), b non-uniform mesh with \(n=198\) and \(q=1.1\). c and d are enlarged views of (a) and (b) respectively

Full size image

5.2.5 Test 6: Shu–Osher Problem

In this test, we compare the performance of our method against the CWENO-type methods. Given that the method proposed by Levy, Puppo and Russo [15] can have negative weights depending on the distribution of the grid and the relative position of the reconstruction point, we will use the alternative proposed in [5], which was designed to avoid this issue as one of its purposes. The main goal of this test is to compare the impact of the computation of the smoothness indicators in the global cost of the method, which has been discussed theoretically in Sect. 3.2.

We start by describing the model of the 1D Euler equations for gas dynamics. If we denote the density as \(\rho \), the velocity as v and the specific energy of the system by E, then the equations are given by \(\varvec{u} = ( \rho , \rho v, E)^{\textrm{T}}\) and \(\varvec{f} (\varvec{u}) = (\rho v, p+\rho v^2, v(E+p))^{\textrm{T}}\). The variable p represents the pressure and is defined by the equation of state

where \(\gamma \) is the adiabatic constant taken as \(\gamma =1.4\). We consider as spatial domain \(\Omega =(-5,5)\), and as initial condition

This initial condition specifies the interaction of a Mach 3 shock with a sine wave and is accompanied by left inflow and right outflow boundary conditions.

We simulate until \(T=1.8\) with a resolution of \(n=256\) cells, \(\text {CFL}=0.5\) and the same parameters as those stated in Sect. 5.2.1 to generate the non-uniform grid. Then, we compare the results against a reference solution computed with a resolution of \(N=65536\) cells. The results of the experiment can be seen in Fig. 5.

Test 6: Numerical solution of the Shu–Osher problem with CWENO-type schemes with uniform (blue) and non-uniform (red) grids with \(n=256\) cells

Full size image

In this simulation, our method is 1.13 times faster than CWENO for the fifth-order scheme and 8.24 times faster for the seventh-order scheme, mostly due to the considerably more complex computations required to obtain the classical Jiang–Shu smoothness indicators on a nonuniform grid, required by the chosen CWENO method. Furthermore, our method produces a better approximation than the CWENO method, as Fig. 5 reflects.

In this paper, we have presented a novel WENO approach capable of performing essentially non-oscillatory reconstructions in a general context of nonuniform grids. Both the theoretical results and the numerical experiments show that the scheme has the desired accuracy.

The efficiency of the algorithm relies strongly on the simplification of the smoothness indicators that can be performed in weight designs akin to the one proposed by Yamaleev–Carpenter [25], as shown in [9]. Another important factor is that the coefficients involving the grid spacing can be computed and stored previously if the numerical simulation uses a fixed grid, yielding in that case a similar efficiency to the uniform grid case. As for the cases in which the grid is not fixed, then Newton polynomials, together with the associated computation of the divided differences properly scaled, are recommended instead.

Our next goal is to tackle the challenge of redesigning the weights of the method so that it can be unconditionally optimal, regardless of the order of the critical point to which the stencil converges, for stencils with four or more nodes, similarly as done in [6] and [7] for the case of uniform stencils.

Another of our future research projects involves a generalization of the current approach into an interpolator capable of preserving the essentially non-oscillatory properties regardless of the relative position of the point to be interpolated or extrapolated with respect to the stencil. This approach can be especially interesting to handle boundary conditions in finite-difference schemes, among other possible applications.