A generalized framework to construct third order WENO weights using weight limiter functions

TL;DR

This paper introduces a novel, generalized framework for constructing third order WENO weights using weight limiter functions, ensuring non-oscillatory solutions and third order accuracy, including at smooth extrema.

Contribution

A new framework for WENO weights based on weight limiter functions that guarantees third order accuracy and non-oscillatory behavior, distinct from existing methods.

Findings

New weights guarantee third order accuracy in smooth regions.

Proposed weights maintain non-oscillatory shock capturing.

Numerical tests confirm effectiveness across various problems.

Abstract

The main aim of this work is not to improve any existing non-linear weight but to give a generalized framework for the construction of non-linear weights to get non-oscillatory third order WENO schemes. It is done by imposing necessary conditions on weights to get non-oscillatory WENO reconstruction which give further insight on the structure of weights to ensure non-occurrence of oscillations and characterize the solution region for third order accuracy. This framework for WENO weights is new and completely different from the prevailing existing approach. New non-linear weights are designed using a function of smoothness parameter termed as weight limiter functions. Many such weight limiter functions are given and analyzed. These new weights are simple and by construction guarantee for exact third order accuracy in smooth solution region including smooth extrema away from critical point. Numerical results for various test problems are given and compared. Results show that proposed weights give third order accuracy without loosing the non-oscillatory shock capturing ability of the resulting scheme.

Tables

Table 1. Table 1: Rate of convergence for IC ( 4.4 ) at final time t = 0.5 𝑡 0.5 t=0.5 with C ​ F ​ L = 0.25 . 𝐶 𝐹 𝐿 0.25 CFL=0.25.

N	WENO3-${\omega }_0^1$	Rate	WENO3-${\omega }_0^1$	Rate
	$L^{\infty }$ error		$L^1$ error
80	1.24987e-02	-Inf	4.44919e-03	-Inf
160	3.81984e-03	1.71	8.45946e-04	2.39
320	5.91893e-04	2.69	9.35838e-05	3.18
640	7.22438e-05	3.03	8.53235e-06	3.46
1280	8.88347e-06	3.02	6.51676e-07	3.71
2560	9.28233e-07	3.26	4.18069e-08	3.96

Table 2. Table 2: Rate of convergence for IC ( 4.4 ) at final time t = 0.5 𝑡 0.5 t=0.5 with C ​ F ​ L = 0.25 𝐶 𝐹 𝐿 0.25 CFL=0.25 .

N	WENO3-${\omega }_{0,5}^1$	Rate	WENO3-${\omega }_{0,5}^1$	Rate
	$L^{\infty }$ error		$L^1$ error
80	8.83849e-03	-Inf	2.68444e-03	-Inf
160	2.53953e-03	1.80	4.69982e-04	2.51
320	3.37745e-04	2.91	4.75547e-05	3.30
640	3.65303e-05	3.21	4.15503e-06	3.52
1280	4.45903e-06	3.03	3.16314e-07	3.72
2560	4.65153e-07	3.26	2.03236e-08	3.96

Table 3. Table 3: Rate of convergence for IC ( 4.5 ) at final time t = 0.5 𝑡 0.5 t=0.5 with C ​ F ​ L = 0.25 𝐶 𝐹 𝐿 0.25 CFL=0.25 .

N	WENO3-${\omega }_0^1$	Rate	WENO3-${\omega }_0^1$	Rate
	$L^{\infty }$ error		$L^1$ error
80	4.76670e-02	-Inf	1.85105e-02	-Inf
160	1.83789e-02	1.37	4.16305e-03	2.15
320	5.61142e-03	1.71	7.71866e-04	2.43
640	9.96515e-04	2.49	9.43049e-05	3.03
1280	1.14091e-04	3.13	8.67465e-06	3.44
2560	1.41383e-05	3.01	7.30714e-07	3.57

Table 4. Table 4: Rate of convergence for IC ( 4.5 ) at final time t = 0.5 𝑡 0.5 t=0.5 with C ​ F ​ L = 0.25 . 𝐶 𝐹 𝐿 0.25 CFL=0.25.

N	WENO3-${\omega }_{0,5}^1$	Rate	WENO3-${\omega }_{0,5}^1$	Rate
	$L^{\infty }$ error		$L^1$ error
80	3.46967e-02	-Inf	1.07368e-02	-Inf
160	1.29620e-02	1.42	2.39369e-03	2.17
320	3.65421e-03	1.83	4.12535e-04	2.54
640	5.80597e-04	2.65	4.72786e-05	3.13
1280	5.80477e-05	3.32	4.20458e-06	3.49
2560	7.12878e-06	3.03	3.54383e-07	3.57

Equations

$$u_t+f(u{)}_x$$

$$u$$

$$\frac{d{u}_i(t)}{dt}=-\frac{1}{\Delta x}\left({\hat{f}}_{i+\frac{1}{2}}-{\hat{f}}_{i-\frac{1}{2}}\right)=:L(u)$$

$$f(u(x))=\frac{1}{\Delta x}{\int }_{x-\frac{\Delta x}{2}}^{x+\frac{\Delta x}{2}}h(\xi )d\xi .$$

$$\frac{\partial f}{\partial x}{\mid }_{x=x_i}=\frac{1}{\Delta x}\left(h_{i+\frac{1}{2}}-h_{i-\frac{1}{2}}\right)$$

$${\hat{f}}_{i\pm \frac{1}{2}}=h_{i\pm \frac{1}{2}}+O(\Delta x^r)$$

$$\frac{d{u}_i(t)}{dt}=-\frac{1}{\Delta x}\left(h_{i+\frac{1}{2}}-h_{i-\frac{1}{2}}\right)\approx -\frac{1}{\Delta x}\left({\hat{f}}_{i+\frac{1}{2}}-{\hat{f}}_{i-\frac{1}{2}}\right)$$

$$f(u)=f^{+}(u)+f^{-}(u)$$

$$f^{\pm }(u)=\frac{1}{2}(f(u)\pm \alpha u)$$

$${\hat{f}}_{i+\frac{1}{2}}={\hat{f}}_{i+\frac{1}{2}}^{+}+{\hat{f}}_{i+\frac{1}{2}}^{-}$$

$$S_k(i):=\{x_{i+k-1},x_{i+k}\},k=0,1$$

$${\hat{f}}_{i+\frac{1}{2}}^k=\sum _{l=0}^1{c}_{k,l}{f}_{i+k-1+l}$$

$${\hat{f}}_{i+\frac{1}{2}}^0=\frac{3}{2}{f}_i-\frac{1}{2}{f}_{i-1}$$

$${\hat{f}}_{i+\frac{1}{2}}^1=\frac{1}{2}{f}_i+\frac{1}{2}{f}_{i+1}$$

$${\hat{f}}_{i+\frac{1}{2}}^0=h_{i+\frac{1}{2}}-\frac{\Delta x^2}{4}{f}^{(2)}(0)+O(\Delta x^3)$$

$${\hat{f}}_{i+\frac{1}{2}}^1=h_{i+\frac{1}{2}}+\frac{\Delta x^2}{4}{f}^{(2)}(0)+O(\Delta x^3)$$

$${\hat{f}}_{i+\frac{1}{2}}=\sum _{k=0}^1{\omega }_k{\hat{f}}_{i+\frac{1}{2}}^k$$

$$\sum _{k=0}^1{\omega }_k=1,{\omega }_k\ge 0.$$

$$h_{i+\frac{1}{2}}=\sum _{k=0}^1{d}_k{\hat{f}}_{i+\frac{1}{2}}^k+O(\Delta x^3)$$

$$d_0=\frac{1}{3},d_1=\frac{2}{3}$$

$${\omega }_k=\frac{{\alpha }_k}{{\sum }_{l=0}^1{\alpha }_l},{\alpha }_k=\frac{d_k}{(ϵ+{\beta }_k{)}^p}$$

$${\beta }_k=\Delta x{\int }_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{\left(\frac{d{\hat{f}}^k}{dx}\right)}^{2}dx,k=0,1.$$

$${\beta }_0=(f_i-f_{i-1}{)}^2$$

$${\beta }_1=(f_{i+1}-fi{)}^2$$

$${\alpha }_k=d_k(1+\frac{\tau }{(ϵ+{\beta }_k)})$$

$$\tau =\left|\frac{{\beta }_0+{\beta }_1}{2}-{\beta }_3\right|,$$

$${\alpha }_k=d_k\left(1+\frac{{\tau }_p}{(ϵ+{\beta }_k)}+\lambda \frac{{\beta }_k+ϵ}{{\tau }_p+ϵ}\right)$$

$$r_i=\frac{{\Delta }_{-}{f}_i}{{\Delta }_{+}{f}_i}.$$

$${\hat{u}}_{i+\frac{1}{2}}={\omega }_0{\hat{u}}_{i+\frac{1}{2}}^0+{\omega }_1{\hat{u}}_{i+\frac{1}{2}}^1$$

$${\hat{u}}_{i+\frac{1}{2}}^0=\frac{3}{2}{u}_i-\frac{1}{2}{u}_{i-1}.$$

$${\hat{u}}_{i+\frac{1}{2}}^1=\frac{1}{2}{u}_i+\frac{1}{2}{u}_{i+1}.$$

$${\omega }_0$$

Full text

Keywords: Hyperbolic conservation laws, Third order WENO reconstructions, Limiter functions, Non-linear weights, Data dependent stability.

AMS subject classifications. 65M06, 65M06, 35L65

1 Introduction

The initial value problem for system of conservation laws which models flow phenomena in Gas Dynamics, Aero-dynamics, Astrophysical modeling, meteorology and weather prediction etc. can be written in one space dimension with given $f(u):\mathbb{R}^{m}\to\mathbb{R}^{m}~{}\text{and}~{}u_{0}:\mathbb{R}\to\mathbb{R}^{m}$ as

[TABLE]

where $u:\mathbb{R}\times(0,\infty)\to\mathbb{R}^{m}$ is the unknown conserved variable. In general the closed form solution is not known for such complex systems which normally have strongly irregular solutions (large jumps, discontinuities) and such problems may involve complicated smooth solution region structures. The numerical treatment of such problems faces typical complexity, for example, oscillations near shock and contact discontinuity which makes a high order numerical scheme unstable. On the other hand in certain situations, the time of evolution of these complex structures are so long that it is impractical to use low order numerical methods to obtain an acceptable resolution. Therefore, development and analysis of high order accurate and non-oscillatory shock capturing numerical schemes for problem (1.1) has been an active area of research. Among existing non-oscillatory shock capturing schemes, total variation diminishing (TVD) and weighted essentially non-oscillatory schemes (WENO) have been of great interest because of their ability to capture the discontinuity with high resolution.

The concept of TVD schemes is introduced first in [1] and was further utilized to develop and analyze high resolution TVD schemes in [2, 3, 4]. The basic idea therein is to use flux limiters to get high order reconstructed value of the flux. More details on flux limited TVD methods and applications can be found in [5, 6]. Unfortunately, though TVD schemes ensure for removal of spurious oscillations from the numerical approximations however they are criticized due to their degenerate accuracy at extrema [7]. This degeneracy causes cornered approximation to smooth solution due to clipping error.

On the other hand, essentially non-oscillatory (ENO) scheme based on cell-average is proposed by Harten et al. in [8]. The idea behind ENO schemes is to choose the smoothest stencil among several candidate stencil for high order accurate approximation of the flux function at cell boundaries and also to avoids oscillations near shocks. However, the cell-average approach to reconstruct point values from given cell average values in multidimensional case is computationally costly compared to the flux version of efficient ENO scheme (ENO schemes based on point values) given in [9, 10]. Later many improved version of ENO schemes has been developed see [11]. The most popular of them is cell averaged weighted ENO scheme (WENO introduced by Liu, Osher and Chan [12]. It uses non-linear convex combination of all interpolating polynomial obtained from candidate stencils of ENO scheme, which results in to a higher order accurate non-oscillatory scheme compared to ENO scheme using the same stencil. The most significant contribution of this technique is the construction of non-linear weights and smoothness indicator based on undivided differences. Later Jiang and Shu [13] introduced finite difference flux version WENO scheme popularly named as WENO-JS scheme by modifying the smoothness measurement and extended the scheme up to $5th$ order accuracy. Later, Henrick et al. found that non-linear weights in WENO-JS failed to recover optimal (ideal) order of accuracy at critical points and they presented mapped WENO scheme (WENO-M) [14]. In mapped WENO the non-linear weights are constructed through the construction of mapping function so that weights remain as close as possible to optimal (Ideal) weights except at highly non-smooth regions. Another approach was developed by Borges et al. [15] by using global smoothness measurement of fifth order WENO scheme named as WENO-Z scheme, having the same accuracy of WENO-M scheme but with less computational cost. Further many modified versions of WENO schemes have been developed to improve the order of accuracy or to get computationally cheaper results by changing the smoothness indicator of the non-linear weights [16, 17, 18]. In [19] energy stable WENO3 schemes are proposed.

Recently active research work has been carried out specific to improvisation of non-linear weights of third order WENO schemes. The classical WENO-JS non-linear weights are modified to achieve optimal third order accuracy regardless of critical point in [20]. These schemes are known as WENO-Z3 schemes. Further, improvement to get faster WENO schemes, WENO-Z3 wights are modified in [21] termed as WENO-N3 which are again modified in [22] to get third order accuracy at critical points and termed and WENO-NP3. In [23], WENO-NP3 weights are modified by changing the smoothness indicator named as WENO-F3. Very recently, in [24], WENO-Z3 weights are significantly modified to get ENO solution near strong discontinuity and termed as WENO-P+3 schemes. It is also shown that WENO-P+3 scheme outperforms WENO3 schemes using existing improved non-linear weights.

The main aim of this work is not to improve any existing weight but to present a generic framework to devise non-linear WENO weights which guarantee by construction for third order accuracy in smooth solution region. This is achieved via analyzing the weights of third order WENO reconstruction in the light of data dependent stability notion given in [25, 26] along with the flux limiters framework. More precisely, simple non-linear weights are designed by using weight limiter as a function of a smoothness measurement which is again a function of consecutive gradient ratio. An important feature is that weight limiter function controls the weights such that they achieve ideal weight for third order accuracy in smooth region including most of extreme points and give non-oscillatory approximation for discontinuity. Various examples of such weight limiter functions are given.

The organization of the paper is as follows: Section 2 provides brief review of WENO scheme and in particular third order WENO weights in section 2.1. The main detailed contribution of work is given section 3. In section 4 demonstrates the computational testing and quantitative performance of the third order WENO scheme using proposed non-linear weights. Conclusion and remark on ongoing work given in section 5.

2 Review of WENO scheme

In this section we briefly describe the flux version third order WENO finite difference scheme for one dimensional scalar conservation law proposed in [13]. Let $I_{i}$ be a partition of a given domain with $i^{th}$ cell $I_{i}=\left[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}\right]$, center of $I_{i}$ is $\displaystyle x_{i}=\frac{x_{i-\frac{1}{2}}+x_{i+\frac{1}{2}}}{2}$ and function value $f$ at node $x_{i}$ is given by $f_{i}=f(x_{i})$. For simplicity $\{x_{i+\frac{1}{2}}\}_{i}$ is uniformly spaced, notation $\Delta x=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$ indicate the size of $I_{i}$. Moreover, notation $u_{i}^{n}$ used for approximation to $u$ at the grid point $(x_{i},t^{n})$ where $t^{n}=n\Delta t$ is the discrete value in the time direction.

The semi-discretized approximation of the one -dimensional hyperbolic conservation law (1.1) can be expressed as:

[TABLE]

where $u_{i}(t)$ is the numerical approximation to the point values $u(x_{i},t)$ and the numerical flux $\hat{f}$ is a function of $(r+s)$ variables i.e., $\hat{f}_{i+\frac{1}{2}}=\hat{f}(u_{i-r},...,u_{i+s})$. The numerical flux function should satisfy Lipschitz continuity in each of its arguments, be consistent with the physical flux $f$, that is, $\hat{f}(u,u,...,u)=f(u)$.

To compute numerical flux $\hat{f}_{i\pm\frac{1}{2}}$, a function $h$ is defined implicitly by the following equation (see Lemma $2.1$ of [10]).

[TABLE]

differentiation of the above equation (2.2) at the point $x=x_{i}$ yields,

[TABLE]

which indicates that the numerical flux $\hat{f}$ approximates $h$ at cell boundaries $x_{i\pm\frac{1}{2}}$ with high order accuracy, that is,

[TABLE]

with $r$ depending on the degree of interpolation.

Using equation (2.3) in equation (2.1), we have

[TABLE]

In order to avoid entropy violating solution and ensure numerical stability flux $f(u)$ is splitted into the two components $f^{+}$ and $f^{-}$ such that

[TABLE]

where $\frac{df^{+}(u)}{du}\geq 0$ and $\frac{df^{-}(u)}{du}\leq 0$. Among many flux splitting methods, the following global Lax-Friedrichs splitting is heavily used for its simplicity and capability to produce very smooth fluxes, which is defined as

[TABLE]

where $\alpha=\max_{u}|f^{\prime}(u)|$. Let $\hat{f}^{+}_{i+\frac{1}{2}}$ and $\hat{f}^{-}_{i+\frac{1}{2}}$ be the numerical fluxes obtained from the positive and negative parts of $f(u)$ respectively and from (2.5) we have

[TABLE]

In the following approximation of $\hat{f}^{+}_{i+\frac{1}{2}}$ is describe for third order accuracy. Note that negative part $\hat{f}^{-}_{i+\frac{1}{2}}$ of (2.7) can be approximated accordingly as it is symmetric to positive part $\hat{f}^{+}_{i+\frac{1}{2}}$ with respect to $x_{i+\frac{1}{2}}$. Also, for simplicity we will drop the ${}^{\prime}+^{\prime}$ sign in the superscript.

2.1 Third order WENO

The construction of $\hat{f}^{+}_{i+\frac{1}{2}}$ for the classical third-order WENO scheme uses a $3$-point stencil $S(i):=\{x_{i+k-1},x_{i+k},x_{i+k+1}\},(k=0,1)$ which is subdivided into two candidate sub-stencils. Let

[TABLE]

be the sub-stencil consisting of $2$-points starting at $x_{i+k-1}$ and let

[TABLE]

be the first degree polynomial approximation constructed on the stencil $S_{k}(i)$ to approximate the value $h(x_{i+\frac{1}{2}})$ where $c_{k,l},~{}l=0,1$ are the Lagrange interpolation coefficients depending on the left-shift parameter $k$. The expression of $\hat{f}^{k}_{i+\frac{1}{2}}$ can be written as

[TABLE]

To define $\hat{f}^{j}_{i-\frac{1}{2}}$ each index needs to be shifted by $-1$. Moreover, the Taylor expansion of the equations in (2.10) are given by

[TABLE]

These functions are combined to define a new WENO approximation to the value $h_{i+\frac{1}{2}}$, that is,

[TABLE]

where non-linear weights $\omega_{k}$ satisfy the following convexity property

[TABLE]

To construct the non-linear weights $\omega_{k}$ it is first considered the case that the function $h$ is is smooth on all the stencil $S_{k}(i)$ with $k=0,1$. After that the constants $d_{k}$ are found such that its linear combination with $\hat{f}^{k}_{i+\frac{1}{2}}$ retains the third order convergence order to $h_{i+\frac{1}{2}}$, that is,

[TABLE]

The coefficients $d_{k}$ are called the ideal or linear weights.The specific values of $d_{k}$ are as follows [27]

[TABLE]

Note that each $d_{k}\geq 0$ and $\sum_{k=0}^{1}d_{k}=1$. The non-linear weights $\omega_{k}$ in (2.12) are constructed such that final reconstruction become essentially non-oscillatory using following designing criteria

• i.

In smooth regions the non-linear weights should converge to the linear weights with required order of accuracy.

• ii.

Weight corresponding to non-smooth or discontinuous region should tend to zero so that contribution from the non-smooth regions in the approximation of $\hat{f}_{i+\frac{1}{2}}$ is negligible

2.1.1 WENO-JS3 Weights [13]

The non-linear WENO-JS3 weights are defined as

[TABLE]

where $\epsilon$ is a positive small number which is set to be $\epsilon=10^{-6}$ to avoid division by zero, $p=2$ is chosen to increase the difference of scales of distinct weights ar non-smooth parts of the solution. Note that $\alpha_{k}$ are the unnormalized weights and $\omega_{k}$ are the normalized weights. In case of third order accuracy the smoothness indicator $\beta_{k}$ can be defined as

[TABLE]

which reduces to

[TABLE]

2.1.2 WENO-Z3 Weights [28]

The WENO-JS3 weights (2.16) were further modified to WENO-Z3. These WENO-Z3 weights are obtained by modifying $\alpha_{k}$ in (2.16) as

[TABLE]

where $\tau=|\beta_{0}-\beta_{1}|$ and $\beta_{k}$ are given by (2.18).

2.1.3 WENO-N3 Weights [21]

A more high resolution and efficient WENO-N3 weights compared to WENO-Z3 are proposed, which are obtained by using

[TABLE]

where $\beta_{3}=\frac{13}{12}(f_{i-1}-2f_{i}+f_{i+1})^{2}+\frac{1}{4}(f_{i-1}-f_{i+1})^{2}$. WENO-N3 weights are further improved to WENO-NP3 and WENO-F3 weights to achieve third order accuracy at critical points in [22, 23]. Though WENO3 scheme using these weights give improved approximation to smooth extrema they exhibits oscillations near discontinuity (See Results of WENO-NP3 in figure 4b)).

2.1.4 WENP-P+3 Weights [24]

In 2018, to improve further resolution of WENO3 scheme near strong discontinuity WENO-N3 weights are modified by defining

[TABLE]

where $\lambda=\Delta_{x}^{\frac{1}{6}}$ and $\tau_{p}=\left|\frac{\beta_{0}+\beta_{1}}{2}-\frac{1}{4}(f_{i-1}-f_{i+1})^{2}\right|.$ Note that though WENO-P+3 scheme outperformed other WENO3 schemes in capturing the strong discontinuity but degenerate to first order accuracy at critical points [24]. Also the shock capturing non-oscillatory behavior of WENP-P+3 scheme heavily depends on choice of parameter $\lambda$ (see Figure 5 in [24]) and parameter $\epsilon$ in (2.21) see Figure 4(a).

3 Construction of new non-linear weights

In this section a new approach is given to construct non-linear weights such that the weighted reconstruction (2.12) achieve essentially non-oscillatory property. Let us define the parameter of the ratio of consecutive gradients as

[TABLE]

It can be observe from (2.10) that in case of linear flux i.e., $f=au,\;a>0$ the reconstructed value at cell interface is a convex combination of second order upwind and centered flux. More precisely in this case (2.12) reduces to

[TABLE]

where $\omega_{0}\geq 0,\omega_{1}\geq 0\;\text{and}\;\omega_{0}+\omega_{1}=1$ . The upwind and centered flux are

[TABLE]

[TABLE]

In order to characterize the weights $\omega$ such that the reconstruction (3.2) be non-oscillatory, we recall our own results on data dependent stability given in [25, 26]

Lemma 3.1**.**

For the linear transport problem (4.1) the Forward in time scheme using second order centered flux (3.4) is data dependent stable and non-oscillatory in the solution data region where $r_{i}\in\mathcal{R}_{c}=(-\infty,-1)\cup\left[\frac{a\lambda}{2-a\lambda},\infty\right)$ for every $i$ where $a\lambda=a\frac{\Delta t}{\Delta x}\leq 1.$

Lemma 3.2**.**

For the linear transport problem (4.1), the forward in time scheme using second order upwind flux (3.3) is data dependent stable and non-oscillatory in the solution data region where $r_{i-1}\in\mathcal{R}_{up}=\left[-\frac{2-3a\lambda}{a\lambda},3\right)$ for every $i$ where $a\lambda=a\frac{\Delta t}{\Delta x}\leq\frac{1}{2}.$

On dropping out the subscript $i$, following can be deduced.

• Characterization of centred weight $\omega_{1}$: Lemma 3.1 gives the characteristics to be satisfied by the weight $\omega_{1}$ associated with the centered flux (3.4) . More precisely, note that for $a\lambda\rightarrow 0$ the non-oscillatory region $\mathcal{R}_{c}$ of centred flux $(-\infty,-1)\cup\left[\frac{a\lambda}{2-a\lambda},\infty\right)\rightarrow(-\infty,-1)\cup(0,\infty)$. It clearly concludes that centered flux is oscillatory at $r=0$ irrespective of choice of $\lambda$. Thus in order to avoid contribution of centered flux to construct non oscillatory third order WENO scheme using (3.2), a necessary condition is $\omega_{1}\rightarrow 0\;\text{for}\;r\rightarrow 0$. However for $r$ away from zero, $\omega_{1}$ can be nonzero i.e, $\omega_{1}\rightarrow 1\;\text{for}\;\rightarrow r\pm\infty.$ In particular, under CFL number $a\lambda=\frac{1}{2}$111Linear stability condition for second order scheme using upwind flux (3.3), since $\mathcal{R}_{c}=(\infty,-1)\cup(\frac{1}{3},\infty)$, therefore centered weight must be defined such that $\omega_{1}\rightarrow 0$ as $r\rightarrow 0$ in the interval $[0,1/3)$ however $\omega_{1}$ can take any value in $[0,1]$ for $\displaystyle r\geq\frac{1}{3}$.
• Characterization of upwind weight $\omega_{0}$: On the other hand, Lemma 3.2 charaterize the weight $\omega_{0}$ associated with the upwind flux (3.3). Note that for $\lambda\rightarrow 0$ non-oscillatory stability region of upwind flux (3.3) i.e., $\mathcal{R}_{up}=\left[-\frac{2-3a\lambda}{a\lambda},3\right)\rightarrow(-\infty,3)$ which suggests that weight $\omega_{0}\rightarrow 0$ for $r\rightarrow\infty$. In particular, under CFL number $a\lambda=\frac{1}{2}$, since $\mathcal{R}_{up}=[-1,3)$ upwind weight $\omega_{0}$ can take any value in $[0,1]$ for $r\in[-1,3]$.

Thus, necessary conditions for non-oscillatory approximation are $\omega_{1}\rightarrow 0$ for $r\rightarrow 0$ and $\omega_{0}\rightarrow 0$ for $r$ away from $[-1,3]$. In particular, for smooth region of solution (which corresponds to $r\approx 1$), convexity property (2.13) along with Lemma 3.1 and 3.2 enables to choose ideal weights $\omega_{0}=\frac{1}{3},\;\omega_{1}=\frac{2}{3}$ for third order accuracy [27]. Moreover, it also follows that third order non-oscillatory approximation can be achieved for $1/3<r<3$ which is the common stability region of centred and upwind flux.

Based on above observations, the following class of non-linear weights is constructed

[TABLE]

where $\chi(r)$ is termed as weight limiter function which must satisfy the following characteristics to yield non-oscillatory third order weights.

• i.

Non-Oscillatory Conditions:

• (a)

$\chi(0)=0$ to achieve $\omega_{0}=1$ and $\omega_{1}=0$ (Upwind only flux).

• (b)

$\chi(\pm\infty)=\frac{3}{2}$ to get $\omega_{0}=0$ and $\omega_{1}=1$ (Centered only flux).

• ii.

Third order Accuracy Condition: $\chi$ should be differentiable at $r=1$ and $\chi(1)=1$ to achieve ideal weights $\omega_{0}=\frac{1}{3},\;\omega_{1}=\frac{2}{3}$.

The above characterization paves the way to design weight limiter functions such that it achieve third order accuracy even for extrema away from critical points. For example

[TABLE]

In Figure 1, weight limiter functions $\chi_{i}$ (3.6) and corresponding weights $\omega_{0}^{i}$ (3.5) are given. Another class of limiter with better discontinuity capturing property can be defined as follows

[TABLE]

where $k\geq 1$ is a constants which guarantee that resulting scheme maintain the third order accuracy in the solution region $1/k\leq|r|\leq k$.

Note that choice $k>3$ violets the non-oscillatory conditions viz $r\leq 3$ of upwind flux in Lemma 3.2 and $r>\frac{1}{3}$ of centered flux in Lemma 3.1. This suggests that $k$ must satisfying $1\leq k\leq 3$. Numerical results also support this restriction (see Figure 3(b) and 6(c)). Thus it can be concluded weight corresponding to $k=3$ gives maximum possible region for third order accuracy without oscillations. In Figure 2, weight limiter functions $\chi_{5}^{k}$ (3.7) and corresponding weights $\omega_{0,5}^{k}$ (3.5) are given for different values of parameter $k$.

Note that all above limiters in (3.6) and (3.7) are constructed such that they satisfy $\chi(\pm 1)=1$ condition thus guarantee for third order accuracy of scheme at smooth solution region $r\approx 1$ including smooth extrema where $r\approx-1$.

Remark 3.1**.**

The weight limiters defined in (3.6) and (3.7) are functions of smoothness parameter (3.1) similar to the flux limiters used in construction of high resolution shock capturing schemes in [2, 6, 5]. Note that, the conditions on flux limiters in [2] are derived by imposing the total variation diminishing (TVD) property which demanded that flux limiter functions must vanish for all point of extrema i.e., $r\leq 0$ which cause degenerate accuracy of the TVD schemes at extrema. Numerically this degeneracy leads to clipping error while approximating smooth solution. However, weight limiter function $\chi$ in (3.5) does not necessarily needed to vanish for all $r<0$. Thus weight in (3.5) may retains third order reconstruction even for solution extrema $r\approx-1$ and thus clipping error can be reduce.

Remark 3.2**.**

The proposed weights using limiter functions 3.6 or 3.7 are free of parameters and gives consistent approximations

4 Numerical results

In this section various standard test problems are considered to analyze accuracy and the non-oscillatory behavior of third order WENO3 scheme with proposed weight limiters. The following name convention is used through out this section.

• •

WENO3-$\omega_{0}^{i}$ denotes the result obtained by weights (3.5) using function $\chi_{i}$ from (3.6).

• •

WENO3-$\omega_{0,5}^{l}$ denotes the result obtained by weights (3.5) using function $\chi_{5}^{k}$ in (3.7) for $k=l$.

• •

WENO-$*$ represents the results obtained by various weights (*) in section 2 e.g., WENO-$JS3$ represents WENO-JS weights (2.16).

4.1 Linear transport equation

Consider the linear transport equation,

[TABLE]

with $a=1$ and along with following initial conditions.

4.1.1 Test for Non-oscillatory property

Example 1 Consider the discontinuous initial condition

[TABLE]

Example 2 Consider the smooth initial condition with sharp turn

[TABLE]

where computational domain is $0\leq x\leq 1,\omega=5\pi$, $x_{c}=0.5$ and $\sigma=0.2$.

In Figure 3, numerical results for example 1 are given using weight limiter (3.7) for $k=1,3$. It can be seen in Figure 3(a) that numerical approximation with both choices of $k$ is non-oscillatory however choice $k=1$ yields more diffusive solution compared to the choice $k=3$. Figure 3(b) shows oscillatory behavior of weight limiter (3.7) for the choice $k=4$. We mention that the numerical solution obtained by weight limiter functions (3.6a)-(3.6d) are not shown as they falls between these two approximations.

In Figure 4, numerical solution by WENO3 scheme using various existing and newly proposed weights are given and compared. It can be noted that proposed WENO-$\omega_{0,5}^{k}$ weights gives better approximation to maxima without oscillations for $k=3$. The results with choice $k=10$ gives improved approximation for smooth maxima but exhibits small undershoot in bottom smooth region between $0.3\leq x\leq 0.4$ and $0.6\leq x\leq 0.7$. This behavior is also true for WENO-$NP3$ and WENO-P+3 for the choice $\epsilon=1e-6$ in (2.21). The result obtained by WENO-P+3 for the choice $\epsilon=1e-40$ in (2.21) does not show any such oscillations though there is a drop in the peak of smooth maxima with flatness compared to the choice $\epsilon=1e-6$. These results show that in order to achieve improved approximation without clipping error for smooth maxima one has to compromise on the non-oscillatory property.

4.1.2 Accuracy test

Consider equation (4.1) with smooth initial conditions

** Example 3**

[TABLE]

Example 4

[TABLE]

In Figure 5, numerical approximation using proposed weights (3.6a) and (3.7) for smooth initial conditions is given and compared with weight WENO-$Z3$ [20] . It can be seen from results that proposed weight $\omega_{0,5}^{3}$ gives significantly better approximation for smooth extrema compared to the WENO-$Z$.

In Tables 1 and 2, the convergence rate in both $L^{1}$ and $L^{\infty}$ error is given for the WENO3 scheme using proposed weights $\omega_{0}^{1}$ and $\omega_{0,5}^{1}$ for the problem (4.1) with initial condition (4.4). Similarly Tables 3 and 4 shows the convergence rate of these weights for the problem (4.1) with initial condition (4.5). We remark that, computationally weights corresponding to all the limiter functions in (3.6) and (3.7) give more or less similar third order accuracy however weight limiters (3.6a), (3.6b) and (3.7) for $k=1$ give better consistent third order convergence rate due to their smooth nature.

4.2 Inviscid Burgers’ equation

Consider the Burger equation

[TABLE]

with discontinuous initial condition

[TABLE]

The exact solution of Burger equation corresponding to initial condition (4.7) consists a left rarefaction and a steady shock at point $x=\frac{1}{3}$. Figure 6-(a) and (b) shows WENO3 scheme using weight limiter $\chi_{5}^{k},(k=3)$ capture the head and tail of the rarefaction more accurately compared to using $\chi_{1}$. Also for both choices of weights, resolution of steady shock is similar and with no oscillations. Figure 6-(c) shows the oscillatory behavior of weight limiter (3.7) $\chi_{5}^{k}$ for choice $k=4$.

4.3 1D Euler equations

The numerical simulations are performed on the one-dimensional Euler equations which are given by

[TABLE]

where $\rho,~{}u,~{}E,~{}p$ are the density, velocity, total energy and pressure respectively. The system (4.8) represents the conservation of mass, momentum and energy. The total energy for an ideal polytropic gas is defined as

[TABLE]

and the eigen values of the Jacobian matrix $A(U)=\frac{\partial F}{\partial U}$ are

[TABLE]

where U=\left(\begin{array}[]{c}\rho\\ \rho u\\ E\end{array}\right),~{}~{}F=\left(\begin{array}[]{c}\rho u\\ \rho u^{2}+p\\ u(E+p)\end{array}\right) and $\gamma$ is the ratio of specific heats and its values is taken as (4.9).

4.3.1 Sod’s shock tube problem

Sod’s shock tube problem [29] is one-dimensional Euler system (4.8) with Riemann data

[TABLE]

in the computational domain $[-5,5]$ at final time $t=1.3$. The density profile of the solution consists of a rarefaction region as well as a shock and contact discontinuity region. In Figure 7(a) and 7(b) results are given using weights $\omega_{0}^{1}$ and $\omega_{0,5}^{3}$ respectively. It can be seen that WENO3-$\omega_{0}^{1}$ captures the shock and contact discontinuities correctly without oscillations whereas WENO3-$\omega_{0,5}^{3}$ shows a small local undershoot at the foot of the rarefaction. It is needed to be mentioned here that this undershoot does not grow on further refinement of the mesh and is similar to the results obtained by WENO3-$YC$ scheme.

4.3.2 Lax’s shock tube problem

Lax’s shock tube problem [30] with initial condition

[TABLE]

is an one-dimensional Euler system of equation (4.8) in the computational domain $[-5,5]$ and is run upto $t=1.3$ with zero gradient boundary conditions. The numerical results of density profiles along with reference solutions are displayed in Figure 8(a) and 8(b). In this test WENO3-$\omega_{0,5}^{3}$ resolves the discontinuities more crisply compared to WENO3-$\omega_{0}^{1}$ with spurious oscillations.

4.3.3 Shu-Osher test

The Shu-Osher test [31] is the one-dimensional Euler system of equation (4.8) with initial condition

[TABLE]

in the computational domain $[-5,5]$and is run upto time $t=1.8$. Solution obtained using weight $\omega_{0,5}^{3}$ is given in In Figure 9(a) and compared in 9(b). It can be seen from Figure 9(b) that WENO3-$\omega_{0,5}^{3}$ resolve the smooth waves significantly better than WENO3-$YC$ and capture the shock crisply with no oscillations.

4.4 2D Riemann gas dynamic problem

Consider the two dimensional Euler equations of motion for gas dynamics

[TABLE]

where $\mathbf{f}$ and $\mathbf{g}$ are fluxes in the $x$ and $y$ direction and $\mathbf{w}:\mathbb{R}\times\mathbb{R}\times(0,\infty)\to\mathbb{R}^{m}$ is the unknown with

[TABLE]

Total energy $E$ an the pressure $p$ are related by the follwing equation

[TABLE]

Here $\rho,u,v,p$ and $E$ are density, components of velocity in the $x$ and $y$ coordinate directions,pressure, and total energy respectively.$\mathbf{w}$ is the vector of conservative variables,$\mathbf{f(w)}$ and $\mathbf{g(w)}$ is $x$ and $y$ direction wise flux component respectively.

4.4.1 Example 1:

The two-dimensional Riemann problem of gas dynamic [32] is defined by initial configuration as

[TABLE]

The numerical solution is computed on the computational square domain $[0,1]\times[0,1]$ with Dirichlet boundary conditions. The sqaure is divided into four quadrants by lines $x=0.5$ and $y=0.5$.And initial data as constant states are defined on each of the four quadrants and evolve upto time $t=0.25$ with $CFL=0.45$ for a grid $400\times 400$. Numerical result of WENO3-$\omega_{0}^{1}$ and WENO3-$\omega_{0,5}^{3}$ are given in Figure 10.

4.4.2 Example 2:

This two-dimensional Riemann problem from [33] is defined by initial constant which is divided by the lines $x=0.8$ and $y=0.8$ as

[TABLE]

The numerical solution is computed on the computational square domain $[0,1]\times[0,1]$ with Dirichlet boundary conditions. The sqaure is divided into four quadrants by lines $x=0.8$ and $y=0.8$.And initial data as constant states are defined on each of the four quadrants and evolve upto time $t=0.8$ with $CFL=0.15$ for a grid $400\times 400$. Numerical results of WENO3-$\omega_{0}^{1}$ and WENO3-$\omega_{0,5}^{3}$ are given in Figure 11.

4.4.3 Explosion problem

The explosion problem proposed in [31] is a circularly symmetric $2D$ problem with initial circular region of higher density and higher pressure. In particular we set the center of the circle to the origin, its radius to $0.4$ and compute on a quadrant $(x,y)\in(0,1.5)\times(0,1.5)$. Inside the circle initial data are $p=1,~{}\rho=1,~{}u=0,~{}v=0$ and outside it is $p=0.1,~{}\rho=0.125,~{}u=0,~{}v=0$. i.e. the gas is initially at rest and its gas constant is $\gamma=\frac{1}{4}$. This problem (evolution of unstable contact at later times) is sensitive to perturbations of the interface and as noted in [31] for the cells which are crossed by the initial interface circle one needs to use area weighted initial density and pressure. In Figure 12, the numerical results for explosion problem are given.

4.4.4 Implosion problem

This shock problem has been presented in [34].In this problem the gas is placed in a square box.The gas has initialy smaller density and pressure inside a smaller square centered at the center of the box than the rest of the box.We use the box $(x,y)\in(-0.3,0.3)\times(-0.3,0.3)$ and the smaller square with corners at $(\pm 0.15,0),(0,\pm 0.15)$. we did the computation only in the upper right quadrant $(x,y)\in(0,0.3)\times(0,0.3)$ of the box with diamond corner box $|x|+|y|<0.15$. Initial data inside the diamond corner box are $p=0.14,~{}\rho=0.125,~{}u=0,~{}v=0$ and outside are $p=1.0,~{}\rho=1,~{}u=0,~{}v=0$. i.e. initial velocities are zero. The gas constant is $\gamma=\frac{1}{4}$. On all four boundaries Reflecting boundary conditions are used. In Figure 13, the numerical results for implosion problem are given.

It can be observed from above four 2D Reimann Problem examples that, both the schemes WENO3-$\omega_{0}^{1}$, and WENO3-$\omega_{0,5}^{3}$ resolves the flow feature however scheme WENO3-$\omega_{0,5}^{3}$ gives significantly better resolution to discontinuities compared to scheme WENO3-$\omega_{0}^{1}$.

5 Conclusion & Future Work

A new framework is given to construct non-linear weight by utilizing necessary conditions for non-oscillatory WENO3 reconstruction. A characterization non-linear weights for constructing third order WENO schemes is done and various nonlinear weights using weight limiting functions are given. Computational results are given and in some cases compared with WENO3 schemes using other established weights. These results show third order accuracy of the resulting schemes and ability to resolve smooth as well discontinuous region of the solution. The work to propose similar framework for fifth order WENO scheme is under progress and partial computational results are obtained. However, an extensive theoretical investigation is further required to get suitable necessary conditions. This work will be reported separately.

Acknowledgment: Authors acknowledge the Science and Engineering Board, New Delhi, India for providing necessary financial support through funded projects File No. EMR/2016/ 000394 and MTR/2017/000187.