A compact subcell WENO limiting strategy using immediate neighbors for Runge-Kutta Discontinuous Galerkin Methods

TL;DR

This paper introduces a compact subcell WENO limiter for Discontinuous Galerkin Methods that uses only immediate neighbors, improving accuracy and efficiency in solving hyperbolic conservation laws.

Contribution

It presents a novel CSWENO limiter utilizing only neighboring cells, enhancing performance over traditional WENO limiters in high-order DG methods.

Findings

Slightly better performance than parent WENO limiter at higher orders

Effective in 1D and 2D Burgers and Euler equations

Maintains accuracy with a compact stencil

Abstract

A compact subcell WENO (CSWENO) limiter is proposed for the solution of hyperbolic conservation laws with Discontinuous Galerkin Method which uses only the immediate neighbors of a given cell. These neighbors are divided into the required stencil for WENO reconstruction and an existing WENO limiting strategy is used. Accuracy tests and results for one-dimensional and two-dimensional Burgers equation and one-dimensional and two-dimensional Euler equations for Cartesian meshes are presented using this limiter. Comparisons with the parent WENO limiter are provided wherever appropriate and the performance of the current limiter is found to be slightly better than the parent WENO limiter for higher orders.

Figures

Tables

Table 1. Table 1: 1D Burgers equation with the initial condition u ​ ( x , 0 ) = 0.5 + sin ⁡ x 𝑢 𝑥 0 0.5 𝑥 u(x,0)=0.5+\sin x , with periodic boundary conditions, t = 0.5 𝑡 0.5 t=0.5 , Nonuniform mesh with K 𝐾 K elements, L 1 subscript 𝐿 1 L_{1} and L ∞ subscript 𝐿 L_{\infty} errors for 𝐏 1 superscript 𝐏 1 \mathbf{P}^{1} , 𝐏 2 superscript 𝐏 2 \mathbf{P}^{2} and 𝐏 3 superscript 𝐏 3 \mathbf{P}^{3} based DG

		DG without limiter				DG with limiter
	$K$	$L_1$ error	Order	$L_{\infty }$ error	Order	$L_1$ error	Order	$L_{\infty }$ error	Order
${𝐏}^1$	20	5.32E-03		1.24E-02		9.16E-03		2.75E-02
	40	1.31E-03	2.02	3.22E-03	1.95	2.53E-03	1.92	6.08E-03	2.18
	80	3.18E-04	2.04	3.22E-03	1.95	6.52E-04	1.96	9.59E-04	2.66
	160	7.82E-05	2.02	2.13E-04	1.94	1.67E-04	1.96	2.62E-04	1.87
	320	1.91E-05	2.03	5.68E-05	1.91	4.45E-05	1.91	6.42E-05	2.03
${𝐏}^2$	20	3.12E-04		4.14E-03		4.89E-04		4.5E-03
	40	4.61E-05	2.76	5.25E-04	2.98	6.12E-05	2.99	4.56E-04	2.94
	80	6.42E-06	2.84	6.89E-05	2.93	8.14E-06	2.91	6.32E-05	2.85
	160	9.24E-07	2.79	9.87E-06	2.8	1.15E-06	2.82	1.02E-05	2.63
	320	1.31E-07	2.82	1.68E-06	2.55	2.05E-07	2.49	2.47E-06	2.05
${𝐏}^3$	20	2.02E-05		4.13E-04		2.05E-05		4.16E-04
	40	1.21E-06	4.06	3.42E-05	3.59	1.23E-06	4.06	3.43E-05	3.60
	80	7.62E-08	3.99	2.83E-06	3.60	7.67E-08	4.00	2.85E-06	3.59
	160	4.87E-09	3.97	2.81E-07	3.33	4.92E-09	3.96	2.81E-07	3.34
	320	2.98E-10	4.03	1.72E-08	4.03	3.06E-10	4.01	1.73E-08	4.02

Table 2. Table 2: 2D Burgers equation with the initial condition u ​ ( x , 0 ) = 0.5 + sin ⁡ ( x + y ) 𝑢 𝑥 0 0.5 𝑥 𝑦 u(x,0)=0.5+\sin(x+y) , with periodic boundary conditions, t = 0.25 𝑡 0.25 t=0.25 , Nonuniform mesh with K × K 𝐾 𝐾 K\times K elements, L 1 subscript 𝐿 1 L_{1} and L ∞ subscript 𝐿 L_{\infty} errors for 𝐏 1 superscript 𝐏 1 \mathbf{P}^{1} , 𝐏 2 superscript 𝐏 2 \mathbf{P}^{2} and 𝐏 3 superscript 𝐏 3 \mathbf{P}^{3} based DG

		DG without limiter				DG with limiter
	$K\times K$	$L_1$ error	Order	$L_{\infty }$ error	Order	$L_1$ error	Order	$L_{\infty }$ error	Order
${𝐏}^1$	20$\times$20	8.72E-03		1.56E-01		1.25E-02		2.9E-01
	40$\times$40	2.04E-03	2.10	4.93E-02	1.66	3.01E-03	2.05	7.41E-02	1.97
	80$\times$80	4.76E-04	2.10	1.37E-02	1.85	7.14E-04	2.08	1.96E-02	1.92
	160$\times$160	1.02E-04	2.22	3.89E-03	1.82	1.62E-04	2.15	5.34E-03	1.88
${𝐏}^2$	20$\times$20	9.1E-04		5.97E-02		9.2E-04		6.13E-02
	40$\times$40	1.13E-04	3.01	8.05E-03	2.89	1.24E-04	2.89	8.21E-03	2.9
	80$\times$80	1.76E-05	2.68	1.28E-03	2.65	1.96E-05	2.66	2.02E-03	2.02
	160$\times$160	2.5E-06	2.82	1.92E-04	2.74	2.68E-06	2.87	2.88E-04	2.81
${𝐏}^3$	20$\times$20	2.14E-04		9.71E-03		2.18E-04		9.92E-03
	40$\times$40	1.34E-05	4.00	6.17E-04	3.98	1.37E-05	3.99	6.31E-04	3.97
	80$\times$80	8.35E-07	4.00	4.05E-05	3.93	8.36E-07	4.03	4.11E-05	3.94
	160$\times$160	5.19E-08	4.01	2.13E-06	4.25	5.19E-08	4.01	2.13E-06	4.27

Table 3. Table 3: 2D Euler equations with the initial condition ρ ​ ( x , y , 0 ) = 1 + 0.2 ​ sin ⁡ ( x + y ) 𝜌 𝑥 𝑦 0 1 0.2 𝑥 𝑦 \rho(x,y,0)=1+0.2\sin(x+y) , u ​ ( x , y , 0 ) = 0.7 𝑢 𝑥 𝑦 0 0.7 u(x,y,0)=0.7 , v ​ ( x , y , 0 ) = 0.3 𝑣 𝑥 𝑦 0 0.3 v(x,y,0)=0.3 and p ​ ( x , y , 0 ) = 1.0 𝑝 𝑥 𝑦 0 1.0 p(x,y,0)=1.0 , with periodic boundary conditions, t = 2 ​ π 𝑡 2 𝜋 t=2\pi , Nonuniform mesh with K × K 𝐾 𝐾 K\times K elements, L 1 subscript 𝐿 1 L_{1} and L ∞ subscript 𝐿 L_{\infty} errors for density with 𝐏 1 superscript 𝐏 1 \mathbf{P}^{1} , 𝐏 2 superscript 𝐏 2 \mathbf{P}^{2} and 𝐏 3 superscript 𝐏 3 \mathbf{P}^{3} based DG

		DG without limiter				DG with limiter
	$K\times K$	$L_1$ error	Order	$L_{\infty }$ error	Order	$L_1$ error	Order	$L_{\infty }$ error	Order
${𝐏}^1$	20$\times$20	2.89E-03		7.42E-03		8.31E-03		3.24E-02
	40$\times$40	4.87E-04	2.57	2.11E-03	1.81	9.2E-04	3.18	7.52E-03	2.11
	80$\times$80	7.11E-05	2.78	6.89E-04	1.62	1.13E-04	3.03	9.13E-04	3.04
	160$\times$160	1.82E-05	1.97	2.16E-04	1.67	3.42E-05	1.72	2.45E-04	1.90
${𝐏}^2$	20$\times$20	9.88E-05		8.4E-04		1.25E-04		8.87E-04
	40$\times$40	1.25E-05	2.98	1.35E-04	2.64	1.64E-05	2.93	1.32E-04	2.75
	80$\times$80	1.54E-06	3.02	1.62E-05	3.06	2.02E-06	3.02	1.81E-05	2.87
	160$\times$160	1.84E-07	3.07	2.09E-06	2.95	2.31E-07	3.13	2.34E-06	2.95
${𝐏}^3$	20$\times$20	2.65E-06		6.87E-05		6.25E-06		7.01E-05
	40$\times$40	1.67E-07	3.99	4.42E-06	3.96	4.15E-07	3.91	6.24E-06	3.49
	80$\times$80	1.05E-08	3.99	2.85E-07	3.96	2.63E-08	3.98	4.32E-07	3.85
	160$\times$160	6.54E-10	4.01	1.65E-08	4.11	1.51E-09	4.12	3.12E-08	3.79

Table 4. Table 4: 2D Euler equations for the Isentropic Vortex problem with periodic boundary conditions, t = 2 𝑡 2 t=2 , Nonuniform mesh with K × K 𝐾 𝐾 K\times K elements, L 1 subscript 𝐿 1 L_{1} and L ∞ subscript 𝐿 L_{\infty} errors for density with 𝐏 1 superscript 𝐏 1 \mathbf{P}^{1} , 𝐏 2 superscript 𝐏 2 \mathbf{P}^{2} and 𝐏 3 superscript 𝐏 3 \mathbf{P}^{3} based DG

		DG without limiter				DG with limiter
	$K\times K$	$L_1$ error	Order	$L_{\infty }$ error	Order	$L_1$ error	Order	$L_{\infty }$ error	Order
${𝐏}^1$	20$\times$20	2.45E-03		1.89E-01		4.37E-03		2.81E-01
	40$\times$40	5.69E-04	2.11	4.71E-02	2.01	1.02E-03	2.10	7.18E-02	1.97
	80$\times$80	1.14E-04	2.32	1.29E-02	1.87	2.49E-04	2.03	1.85E-02	1.96
	160$\times$160	2.21E-05	2.37	3.34E-03	1.95	5.69E-05	2.13	4.23E-03	2.13
${𝐏}^2$	20$\times$20	6.43E-04		9.5E-02		2.82E-03		1.55E-01
	40$\times$40	7.87E-05	3.03	1.17E-02	3.02	2.98E-04	3.24	1.83E-02	3.08
	80$\times$80	8.92E-06	3.14	1.32E-03	3.15	3.76E-05	2.99	2.57E-03	2.83
	160$\times$160	1.17E-06	3.10	1.72E-04	2.94	4.18E-06	3.17	3.49E-04	2.88
${𝐏}^3$	20$\times$20	4.79E-06		8.35E-05		7.23E-06		8.37E-05
	40$\times$40	4.63E-07	3.96	5.98E-06	3.80	4.63E-07	3.97	6.33E-06	3.73
	80$\times$80	3.12E-08	3.99	4.16E-07	3.85	3.12E-08	3.89	4.24E-07	3.90
	160$\times$160	1.81E-09	4.04	2.72E-08	3.93	1.81E-09	4.11	3.25E-08	3.71

Equations

$$\frac{\partial u}{\partial t}+\frac{\partial f(u)}{\partial x}=0,x\in [x_L,x_R]=\mathbf{D}$$

$$u(x,0)=u_0(x),$$

$$u_h^k(x,t)=\sum _{n=0}^N{\hat{u}}_n^k(t){\psi }_n^k(x)\forall x\in {\mathbf{I}}^k$$

$$f_h^k(u_h^k)=\sum _{n=0}^N{\hat{f}}_n^k(t){\psi }_n^k(x)\forall x\in {\mathbf{I}}^k$$

$$x(r)=x_l^k+\frac{1+r}{2}{h}^k,h^k=x_r^k-x_l^k\forall r\in \mathbf{I}=[-1,1]$$

$$\frac{d}{dt}{\hat{u}}_h^k=({\mathbf{M}}^k{)}^{-1}({\mathbf{S}}^k{)}^T{\hat{f}}_h^k(u_h^k)-({\mathbf{M}}^k{)}^{-1}(f^{\ast }{|}_{r_{N_p}}{e}_{N_p}-f^{\ast }{|}_{r_1}{e}_1)$$

$${\mathbf{M}}^k=\left[M_{ij}^k\right]=\left[{\int }_{x_l^k}^{x_r^k}{\psi }_i^k(x){\psi }_j^k(x)\text{dx}\right]$$

$${\mathbf{S}}^k=\left[S_{ij}^k\right]=\left[{\int }_{x_l^k}^{x_r^k}{\psi }_i^k(x)\frac{d{\psi }_j^k(x)}{dx}\text{dx}\right]$$

$$Q(x_G)=\sum _{i=0}^N{\gamma }_i{p}_i(x_G)$$

$$\left[L_i^2\right]=\left[T_{ij}^3\right]\left[C_j^2\right]$$

$$\left[L_i^2\right]={\left[\begin{array}{cccc}{p}_0(x_G)& p_1(x_G)& p_1(x_G)& Q(x_G)\end{array}\right]}^T$$

$$\left[C_j^2\right]={\left[\begin{array}{ccccc}{\hat{u}}_0^{j-2}& {\hat{u}}_0^{j-1}& {\hat{u}}_0^j& {\hat{u}}_0^{j+1}& {\hat{u}}_0^{j+2}\end{array}\right]}^T$$

$$\left[T_{ij}^3\right]=\left[\begin{array}{ccccc}-\frac{1}{4}& \frac{13}{12}& \frac{1}{6}& 0& 0\\ 0& \frac{1}{2}& \frac{2}{3}& -\frac{1}{6}& 0\\ 0& 0& \frac{13}{6}& -\frac{23}{12}& -\frac{3}{4}\\ -\frac{1}{10}& \frac{21}{30}& \frac{17}{30}& -\frac{13}{60}& \frac{1}{20}\end{array}\right]$$

$${\gamma }_0=\frac{2}{5},{\gamma }_1=\frac{24}{45},{\gamma }_2=\frac{1}{15}$$

$$\left[L_i^2\right]=\left[T_{ij}^4\right]\left[C_j^2\right]$$

$$\left[T_{ij}^4\right]=\left[\begin{array}{ccccc}-\frac{3+6\sqrt{5}}{60}& \frac{5+18\sqrt{5}}{60}& \frac{58-12\sqrt{5}}{60}& 0& 0\\ 0& \frac{2\sqrt{5}-1}{30}& \frac{16}{15}& -\frac{1+2\sqrt{5}}{30}& 0\\ 0& 0& \frac{58+12\sqrt{5}}{60}& \frac{5-18\sqrt{5}}{60}& \frac{6\sqrt{5}-3}{60}\\ \frac{15-69\sqrt{5}}{3000}& \frac{63\sqrt{5}-29}{600}& \frac{163}{150}& -\frac{63\sqrt{5}+29}{600}& \frac{15+69\sqrt{5}}{3000}\end{array}\right]$$

$${\gamma }_0=\frac{235-33\sqrt{5}}{950},{\gamma }_1=\frac{48}{95},{\gamma }_2=\frac{235+33\sqrt{5}}{950}$$

$$\left[L_i^2\right]=\left[T_{ij}^5\right]\left[C_j^2\right]$$

$$\left[T_{ij}^5\right]=\left[\begin{array}{ccccc}-\frac{3}{4}& -\frac{23}{12}& \frac{13}{6}& 0& 0\\ 0& -\frac{1}{6}& \frac{2}{3}& \frac{1}{2}& 0\\ 0& 0& \frac{1}{6}& \frac{13}{12}& -\frac{1}{4}\\ \frac{1}{20}& -\frac{13}{60}& \frac{17}{30}& \frac{21}{30}& -\frac{1}{10}\end{array}\right]$$

$${\gamma }_0=\frac{1}{15},{\gamma }_1=\frac{24}{45},{\gamma }_2=\frac{2}{5}$$

$$\left[L_i^2\right]=\left[T_{ij}^6\right]\left[C_j^2\right]$$

$$\left[T_{ij}^6\right]=\left[\begin{array}{ccccc}\frac{6\sqrt{5}-3}{60}& \frac{5-18\sqrt{5}}{60}& \frac{58+12\sqrt{5}}{60}& 0& 0\\ 0& -\frac{1+2\sqrt{5}}{30}& \frac{16}{15}& \frac{2\sqrt{5}-1}{30}& 0\\ 0& 0& \frac{58-12\sqrt{5}}{60}& \frac{5+18\sqrt{5}}{60}& -\frac{3+6\sqrt{5}}{60}\\ \frac{15+69\sqrt{5}}{3000}& -\frac{63\sqrt{5}+29}{600}& \frac{163}{150}& \frac{63\sqrt{5}-29}{600}& \frac{15-69\sqrt{5}}{3000}\end{array}\right]$$

$${\gamma }_0=\frac{235+33\sqrt{5}}{950},{\gamma }_1=\frac{48}{95},{\gamma }_2=\frac{235-33\sqrt{5}}{950}$$

$${\beta }_i=\sum _{l=1}^N{\int }_{I_j}\Delta x_j^{2l-1}{\left(\frac{{\partial }^l}{\partial x^l}{p}_i(x)\right)}^{2}dx$$

$${\beta }_0=\frac{100({\hat{u}}_0^j{)}^2+\left(204{\hat{u}}_0^{j-2}-404{\hat{u}}_0^{j-1}\right){\hat{u}}_0^j+433({\hat{u}}_0^{j-1}{)}^2-462{\hat{u}}_0^{j-1}{\hat{u}}_0^{j-2}+129({\hat{u}}_0^{j-2}{)}^2}{12}$$

$${\beta }_1=\frac{43({\hat{u}}_0^{j+1}{)}^2+\left(70{\hat{u}}_0^{j-1}-156{\hat{u}}_0^j\right){\hat{u}}_0^{j+1}+156({\hat{u}}_0^j{)}^2-156{\hat{u}}_0^{j-1}{\hat{u}}_0^j+43({\hat{u}}_0^{j-1}{)}^2}{9}$$

$${\beta }_2=\frac{129({\hat{u}}_0^{j+2}{)}^2+\left(204{\hat{u}}_0^j-462{\hat{u}}_0^{j+1}\right){\hat{u}}_0^{j+2}+433({\hat{u}}_0^{j+1}{)}^2-404{\hat{u}}_0^j{\hat{u}}_0^{j+1}+100({\hat{u}}_0^j{)}^2}{12}$$

$${\omega }_i=\frac{{\bar{\omega }}_i}{{\sum }_i{\bar{\omega }}_i},\overline{{\omega }_i}=\frac{{\gamma }_i}{{\sum }_i(ϵ+{\beta }_i{)}^2}$$

$$u_G=\sum _{j=0}^N{\omega }_i{p}_i(x_G)$$

$${\hat{u}}_i^j=\Delta x_j\sum _G{w}_{G}u(x_G){\psi }_i^j(x_G)i=1,\dots ,N$$

Full text

Keywords: discontinuous Galerkin method, troubled cell indicator, limiting, WENO reconstruction, Quadrature Points

1 Introduction

In this paper, we look at the solution of hyperbolic conservation laws with the Runge Kutta Discontinuous Galerkin (RKDG) method [1]. The main advantage of this method is that, in a given cell, to advance the solution in time, we need the information from only the immediate neighbors. For higher orders, to control spurious oscillations near discontinuities, a limiter is used and weighted essentially non oscillatory (WENO) limiters are preferred as they maintain the order of the scheme. However, standard WENO limiters need the information from neighbors of the neighboring cells and the advantage of RKDG method is lost.

Zhong and Shu [2] addressed this issue and used the whole DG polynomial in a given cell for WENO reconstruction using only the immediate neighbors. Dumbser et al [3] used a different strategy where the target cell is divided into subcells and an a posteriori limiting strategy is used based on the Multi-dimensional Optimal Order Detection (MOOD) approach. This gives very good results but it is quite expensive computationally.

We propose a different strategy, where the immediate neighbors are divided into subcells based on the order of the scheme to get the required stencil, appropriate values are assigned to the new cells and the framework given in Qiu and Shu [4] is used for limiting. We call the limiting strategy as compact subcell WENO limiter or CSWENO limiter in short.

The paper is organized as follows. We describe the formulation of the discontinuous Galerkin method used for all our results in section 2, the proposed limiting procedure is described in section 3 and finally the testing of the limiter and the results are described in section 4 and we conclude the paper in section 5.

2 Formulation of Discontinuous Galerkin Method

Consider the nonlinear scalar conservation law as given below

[TABLE]

with the initial condition

[TABLE]

We look at solving (1) using the Discontinuous Galerkin method. We approximate the domain $\mathbf{D}$ by $K$ non overlapping elements whose domain is given by $\mathbf{I}^{k}=[x_{l}^{k},x_{r}^{k}]$. We will approximate the local solution as a polynomial of order $N=N_{p}-1$, where $N_{p}$ is the number of degrees of freedom of the approximation. This is termed to be $\mathbf{P}^{N}$ based Discontinuous Galerkin method. The approximation is given as:

[TABLE]

Here, $u_{h}^{k}(x,t)$ is the approximate local polynomial solution, $\psi_{n}^{k}(x)$ is the local polynomial basis of approximation and $\hat{u}^{k}_{n}(t)$ are the degrees of freedom.

Similarly, we will also approximate the flux $f(u)$ in the domain $\mathbf{D}$ as given below:

[TABLE]

We have used the orthonormalized Legendre polynomials as the local polynomial basis as suggested by Hesthaven et al [5]. The following affine mapping is employed.

[TABLE]

The degrees of freedom $\hat{u}^{k}_{n}$ can be advanced in time by the following scheme obtained from the weak form of the governing equation:

[TABLE]

Here, $\hat{u}_{h}^{k}=[\hat{u}_{0}^{k}\ldots\hat{u}_{N}^{k}]^{T}$, $e_{i}$ is a vector of dimension $N_{p}$ which has zero entries everywhere except at the $i$th location, and $\mathbf{M}^{k}$ is the local mass matrix which is given as:

[TABLE]

and $\mathbf{S}^{k}$ is the local stiffness matrix which is given by:

[TABLE]

Also, $f^{*}$ is the monotone numerical flux at the interface which is calculated using an exact or approximate Riemann solver. We have used the Lax-Friedrichs flux for all the test cases given below.

3 Proposed limiting Procedure

The common method for limiting in Discontinuous Galerkin method is:

1) Identify the cells which need to be limited. They are known as troubled cells.

2) Replace the solution polynomial in the troubled cell with a new polynomial that is less oscillatory but with the same cell average and order of accuracy.

For the first step, we have used the KXRCF troubled cell indicator [6] for all the calculations done in this paper as it is rated highly by Qiu and Shu in [7] on the basis of it’s performance in detecting the discontinuities in various test problems. The second step is where we do the limiting process. We will follow the method given in [4], where we use the finite volume WENO reconstruction procedure in the troubled cell, but with a new approach.

After identifying the troubled cells, we would like to reconstruct the moments, i.e., the values of $\hat{u}_{n}^{j}$ as given in (2) for the troubled-cell $\mathbf{I}_{j}$ for $n=1,\ldots,N$. That is, we retain the cell average $\hat{u}_{0}^{j}$ and reconstruct all the other degrees of freedom. This is done as given in [4] except for a modification which we suggest is that the stencil used for the reconstruction consists of only the immediate neighbors. For the $\mathbf{P}^{1}$ based DGM, we will need to use WENO3 reconstruction and we use the two-point Gauss quadrature rule to reconstruct the solution polynomial. Here, the procedure for limiting remains the same as in [4] which is a standard WENO reconstruction using the cell averages.

For the $\mathbf{P}^{2}$ based DGM, we have to use WENO5 reconstruction method. We use the four-point Gauss-Lobatto quadrature rule to reconstruct the solution polynomial, but instead of a stencil containing five cells as given in [4], we propose the following compact stencil. We take just three cells and divide the neighbors of the troubled-cell on both sides in half. This way, we get the new five point stencil required for the WENO5 reconstruction procedure as shown in Figure 1. Here, the neighbors $\mathbf{I}_{j-1}$ and $\mathbf{I}_{j+1}$ are divided in half to obtain the new cells $\mathbf{I}_{j1}$, $\mathbf{I}_{j2}$, $\mathbf{I}_{j3}$ and $\mathbf{I}_{j4}$ respectively. Now, we use the DG polynomial for $\mathbf{I}_{j-1}$ and take its average in the new cells $\mathbf{I}_{j1}$ and $\mathbf{I}_{j2}$ to obtain the cell averages for the new cells. Similarly, we can calculate the cell average for the cells $\mathbf{I}_{j3}$ and $\mathbf{I}_{j4}$. We can use the same procedure to divide the neighboring cells for $\mathbf{P}^{N}$ based DGM for any $N$ based on the required number of quadrature points.

Note that the grid size is halved for the newly formed cells $\mathbf{I}_{j1}$, $\mathbf{I}_{j2}$, $\mathbf{I}_{j3}$ and $\mathbf{I}_{j4}$. Now, we will proceed as in [4]. The step wise details of the WENO reconstruction are given below:

Step 1: For the $\mathbf{P}^{N}$ based DGM, we need a Gauss or Gauss-Lobatto quadrature rule which is accurate to at least $O(h^{2N+2})$. This means that the order of accuracy for the WENO reconstruction has to be at least $2N+1$. Therefore, for the $\mathbf{P}^{1}$ based DGM, we need the two point Gauss quadrature points given by $r=-1/\sqrt{3}$ and $r=1/\sqrt{3}$. Similarly, for $\mathbf{P}^{2}$ based DGM, we need the four point Gauss-Lobatto quadrature points given by $r=-1$, $r=-1/\sqrt{5}$, $r=1/\sqrt{5}$ and $r=1$ and for $\mathbf{P}^{3}$ based DGM, we use the four point Gauss quadrature points given by $r=-\sqrt{525+70\sqrt{30}}/35$, $r=-\sqrt{525-70\sqrt{30}}/35$, $r=\sqrt{525-70\sqrt{30}}/35$ and $r=\sqrt{525+70\sqrt{30}}/35$.

Step 2: The WENO reconstruction as given in [8] is described now which is to be performed in the given troubled cell $\mathbf{I}_{j}$. We identify $N+1$ small stencils $S_{i}$, $i=0,\ldots,N$ such that $\mathbf{I}_{j}$ belongs to each of them. We set $S_{i}=\bigcup_{l=0}^{N}\mathbf{I}_{j+i-l}$ with the understanding that for the $\mathbf{P}^{2}$ based DGM, $\mathbf{I}_{j-2}$, $\mathbf{I}_{j-1}$, $\mathbf{I}_{j+1}$ and $\mathbf{I}_{j+2}$ are replaced by $\mathbf{I}_{j1}$, $\mathbf{I}_{j2}$, $\mathbf{I}_{j3}$ and $\mathbf{I}_{j4}$ respectively as explained above. We also have the larger stencil $\mathbf{T}=\bigcup_{i=0}^{N}S_{i}$ which contains all the cells from the smaller stencils $S_{i}$.

Now, we have a polynomial of degree $N$, $p_{i}(x)$ corresponding to the stencil $S_{i}$ such that it’s cell average in each of the cells of the stencil $S_{i}$ agrees with the given cell average of $u$. We also have a polynomial of degree $2N$ reconstruction denoted by $Q(x)$ associated with the larger stencil $\mathbf{T}$, such that the cell average of $Q(x)$ in each of the cells of the stencil $\mathbf{T}$ agrees with the cell average of $u$ for that cell. The details of the construction of $p_{i}(x)$ and $Q(x)$ are given in [8].

Step 3: Next, we find the linear weights denoted by $\gamma_{0},\ldots,\gamma_{N}$, which satisfy

[TABLE]

where $x_{G}$ is a Gauss or Gauss-Lobatto quadrature point. A set of linear weights for each of the quadrature points is obtained. The value of the functions $Q(x)$ and $p_{i}(x)$ for each $i$ can be written as a function of the cell average of each cell in the stencil. This is used in WENO reconstruction. For example, for the $\mathbf{P}^{2}$ based DGM, for the Gauss-Lobatto quadrature point $r=-1$, we have:

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

The linear weights are given by

[TABLE]

For the Gauss-Lobatto quadrature point $r=-1/\sqrt{5}$, we have:

[TABLE]

where

[TABLE]

The linear weights are given by

[TABLE]

Similarly, for the Gauss-Lobatto quadrature point $r=1$:

[TABLE]

where

[TABLE]

The linear weights are given by

[TABLE]

Finally, for the Gauss-Lobatto quadrature point $r=1/\sqrt{5}$,

[TABLE]

where

[TABLE]

The linear weights are given by

[TABLE]

Step 4: As given by [8], we compute the smoothness indicator for each stencil $S_{i}$:

[TABLE]

For example, for the $\mathbf{P}^{2}$ based DGM, the smoothness indicators are given as:

[TABLE]

[TABLE]

[TABLE]

Step 5: We compute the nonlinear weights as given below:

[TABLE]

Here $\epsilon$ is a small number which is usually taken to be $10^{-6}$. The final WENO approximation is given by

[TABLE]

Step 6: We obtain the reconstructed degrees of freedom based on the reconstructed point values $u(x_{G})$ at the Gauss or Gauss-Lobatto quadrature points $x_{G}$ and a numerical integration as

[TABLE]

where $w_{G}$’s are the Gaussian quadrature weights for the points $x_{G}$. For a non-orthonormal basis, we define

[TABLE]

, $\mathbf{B}^{j}=\left[D_{1}^{j}-\hat{u}_{0}^{j}\left[M_{01}^{j}\right]\quad\ldots\quad D_{N}^{j}-\hat{u}_{0}^{j}\left[M_{0N}^{j}\right]\right]^{T}$, $\mathbf{X}^{j}=\left[\hat{u}_{0}^{j}\ldots\hat{u}_{N}^{j}\right]^{T}$ and

[TABLE]

Here, the terms $\left[M_{mn}^{j}\right]$ are given by (6). Then the reconstructed degrees of freedom are given by $\mathbf{X}^{j}=(\mathbf{A}^{j})^{-1}\mathbf{B}^{j}$. This will work for any polynomial basis. Now, we can get the reconstructed polynomial solution in $\mathbf{I}_{j}$ by (2). This completes the WENO limiting procedure. We call this limiting procedure the compact subcell WENO limiting or CSWENO limiting in short.

When you have a system of equations of the form (1) to solve, in order to achieve better results, the limiter is used with a local characteristic field decomposition as explained in [8].

For the two dimensional case, we reconstruct the values of the required function $u$ in the troubled cells at the tensor-product Gauss or Gauss-Lobatto quadrature points for the rectangular structured grids. For this, we get two polynomials, one each in $x$ and $y$ directions. Now, we apply the procedure given in [2] for time integration. For solving a system of equations, we use this with a local characteristic field decomposition.

Now, the semi-discrete scheme given in (5) along with the limiter is discretized in time by using the TVD Runge-Kutta time discretization introduced in [9]. We have used a third order TVD Runge-Kutta time discretization for all our calculations.

4 Results

In this section, we look at some of the results obtained to demonstrate the performance of the limiter (called the compact subcell WENO limiter or CSWENO limiter) described above. All the results are obtained using RKDG method and the CSWENO limiter with a third order TVD Runge-Kutta scheme for time integration unless otherwise specified.

4.1 Accuracy Tests

We test the accuracy of the schemes with the CSWENO limiter for scalar and system problems for both one-dimensional and two-dimensional test cases. We present the results of the accuracy tests using one and two-dimensional Burgers equations and one and two-dimensional nonlinear Euler equations. We used both uniform and non-uniform meshes for all the test cases. The non uniform meshes are obtained using a $10\%$ random perturbation of each node of the uniform mesh. We show only the results with nonuniform meshes as representative test cases.

Example 1: We solve the one dimensional nonlinear scalar inviscid Burgers equation:

[TABLE]

with the initial condition $u(x,0)=0.5+\sin x$, with periodic boundary conditions. The solution is smooth till $t=1.0$. The exact solution can be obtained using the Newton-Raphson method as given in [10]. The errors and numerical orders of accuracy are calculated at $t=0.5$ and are presented in Table 1. We can see that the CSWENO limiter maintains the order and magnitude of accuracy of the original DG method.

Example 2: We solve the two dimensional nonlinear scalar inviscid Burgers equation:

[TABLE]

with the initial condition $u(x,y,0)=0.5+\sin(x+y)$, with periodic boundary conditions. The solution is smooth till $t=0.5$. The exact solution is one-dimensional based on $x+y$ and can be calculated in a similar way to the one-dimensional problem. The errors and numerical orders of accuracy are calculated at $t=0.25$ and are presented in Table 2. Again, we can see that the CSWENO limiter maintains the order and magnitude of accuracy of the original DG method.

Example 3: We solve the two dimensional Euler equations:

[TABLE]

where $\textbf{U}=(\rho,\rho u,\rho v,E)^{T}$, $\textbf{f(U)}=u\textbf{U}+(0,p,0,pu)^{T}$ and $\textbf{g(U)}=v\textbf{U}+(0,0,p,pv)^{T}$ with $p=(\gamma-1)(E-\frac{1}{2}\rho(u^{2}+v^{2}))$ and $\gamma=1.4$. Here, $\rho$ is the density, $(u,v)$ is the velocity, $E$ is the total energy and $p$ is the pressure. The initial conditions are given by $\rho(x,y,0)=1+0.2\sin(x+y)$, $u(x,y,0)=0.7$, $v(x,y,0)=0.3$ and $p(x,y,0)=1.0$ and we use periodic boundary conditions. The exact solution is given by $\rho(x,y,0)=1+0.2\sin(x+y-t)$, $u(x,y,0)=0.7$, $v(x,y,0)=0.3$ and $p(x,y,0)=1.0$. The errors in density and numerical orders of accuracy are calculated at $t=2\pi$ and are presented in Table 3. Again, we can see that the CSWENO limiter maintains the order and magnitude of accuracy of the original DG method.

Example 4: We again solve the two dimensional Euler equations given by (21) for the Isentropic Euler Vortex problem suggested by Shu [8] as a test case.The exact solution is given by:

$\rho=\left(1-\left(\frac{\gamma-1}{16\gamma\pi^{2}}\right)\beta^{2}e^{2(1-r^{2})}\right)^{\frac{1}{\gamma-1}}$, $u=1-\beta e^{(1-r^{2})}\frac{y-y_{0}}{2\pi}$, $v=\beta e^{(1-r^{2})}\frac{x-x_{0}-t}{2\pi}$, and $p=\rho^{\gamma}$, where $r=\sqrt{(x-x_{0}-t)^{2}+(y-y_{0})^{2}}$, $x_{0}=5$, $y_{0}=0$, $\beta=5$ and $\gamma=1.4$. We initialize with the exact solution at $t=0$ and use periodic boundary conditions at the edges of the domain. The errors in density and numerical orders of accuracy are calculated at $t=2$ and are presented in Table 4. Again, we can see that the CSWENO limiter maintains the order and magnitude of accuracy of the original DG method.

4.2 Test Cases With Shocks

We will now test the CSWENO (compact subcell WENO) limiter for problems with solutions having shocks. Here also, we have used both uniform and nonuniform meshes and obtained similar results. We will only show the results with uniform meshes. In all the results, we have compared the solution obtained with the CSWENO limiter against that obtained with the WENO limiter given in [4] (also called the parent WENO limiter) and the simple WENO (also called SWENO) limiter described in [2]. As the described limiter is identical to WENO limiter for the $\mathbf{P}^{1}$ based DGM, we have only showed the comparisons for $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM. The WENO limiter uses a standard WENO reconstruction while the simple WENO limiter reconstructs the whole DG polynomial. In order to not clutter the comparison of the results, we have used fixed intervals between all the marked solution points and represented them only with lines otherwise.

Example 5: We solve the same nonlinear Burgers equation given in (19) as in Example 1 with the same initial condition $u(x,0)=0.5+\sin x$, with periodic boundary conditions. We now plot the results at $t=1.5$ when a shock has already appeared in the solution. The computed solution obtained at $t=1.5$ using 80 elements while using the CSWENO limiter for $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM is compared and plotted against the solution obtained using the WENO limiter, simple WENO limiter and the exact solution in Figures 2a and 2b. From the Figures, we can see that the performance of SWENO limiter is better than the CSWENO limiter and WENO limiter and the performance of CSWENO limiter is much better than the parent WENO limiter for both $\mathbf{P}^{2}$ based DGM and $\mathbf{P}^{3}$ based DGM.

Example 6: We solve the nonlinear nonconvex scalar Buckley-Leverett problem

[TABLE]

with the initial condition $u=1$ for $-1/2\leq x\leq 0$ and $u=0$ everywhere else. The exact solution is a shock-rarefaction-contact discontinuity mixture. The computed solution obtained at $t=0.4$ using 80 elements while using the CSWENO limiter for $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM is compared and plotted against the solution obtained using the parent WENO limiter, simple WENO limiter and the exact solution in Figures 3a and 3b. From the Figures, we can see that the performance of all three limiters is quite similar for $\mathbf{P}^{2}$ based DGM. For $\mathbf{P}^{3}$ based DGM, performance of the SWENO limiter is slightly better than the CSWENO limiter and the parent WENO limiter and the performance of the CSWENO limiter is much better than the parent WENO limiter.

Example 7: We now consider a one-dimensional system particularly the Euler equations for an ideal gas given by

[TABLE]

where $\textbf{U}=(\rho,\rho u,E)^{T}$ and $\textbf{f(U)}=u\textbf{U}+(0,p,pu)^{T}$ with $p=(\gamma-1)(E-\frac{1}{2}\rho u^{2})$ and $\gamma=1.4$. We will solve the Riemann problem in the domain $0\leq x\leq 1$ with the initial conditions given by Sod [11] as $(\rho_{L},u_{L},p_{L})=(1,0,1)$ for $x<0.5$ and $(\rho_{R},u_{R},p_{R})=(0.125,0,0.1)$ for $x\geq 0.5$. The computed solution for density obtained at $t=0.2$ using 200 grid points while using the CSWENO limiter for $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM is compared and plotted against the solution obtained using the parent WENO limiter, simple WENO limiter and the exact solution in Figures 4a and 4b. From the Figures, we can see that for $\mathbf{P}^{2}$ based DGM, performance of the SWENO limiter is slightly better than the CSWENO limiter and the parent WENO limiter and the performance of the CSWENO limiter is very similar to the parent WENO limiter. Also, the performance of all three limiters is quite similar for $\mathbf{P}^{3}$ based DGM.

Example 8: We again solve the Euler equations for an ideal gas as given by (23) for the Riemann problem in the domain $0\leq x\leq 1$ with the initial conditions given by Lax [12] as $(\rho_{L},u_{L},p_{L})=(0.445,0.698,3.528)$ for $x<0.5$ and $(\rho_{R},u_{R},p_{R})=(0.5,0,0.571)$ for $x\geq 0.5$. The computed solution for density obtained at $t=0.1$ using 200 elements while using the CSWENO limiter for $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM is compared and plotted against the solution obtained using the parent WENO limiter, simple WENO limiter and the exact solution in Figures 5a and 5b. Here, we can see that the simple WENO limiter performs slightly better than the parent WENO and the CSWENO limiter for both $\mathbf{P}^{2}$ based DGM and $\mathbf{P}^{3}$ based DGM and the CSWENO limiter performs quite better than the parent WENO limiter.

Example 9: We solve the problem of shock interaction with entropy waves as proposed in [13]. We solve the Euler equations (23) with a moving shock interacting with sine waves in density in the domain $0\leq x\leq 1$ with the initial conditions given as $(\rho,u,p)=(3.857143,2.629369,10.333333)$ for $x<0.125$ and $(\rho,u,p)=(1.0+0.2\sin(16\pi x),0,1)$ otherwise. The computed solution for density obtained at $t=0.178$s using 200 elements while using the CSWENO limiter for $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM is compared and plotted against the solution obtained using the parent WENO limiter, simple WENO limiter and the exact solution in Figures 6a and 6b. From the Figures, we can see that for both $\mathbf{P}^{2}$ based DGM and $\mathbf{P}^{3}$ based DGM, the simple WENO limiter performs slightly better than the parent WENO and the CSWENO limiter. Also, the CSWENO limiter performs much better than the parent WENO limiter.

Example 10: We solve the interaction of blast waves of Euler equation as proposed in [12]. We solve the Euler equations (23) in the domain $0\leq x\leq 1$ with the initial conditions given as $(\rho,u,p)=(1,0,1000)$ for $0\leq x<0.1$, $(\rho,u,p)=(1,0,0.01)$ for $0.1\leq x<0.9$ and $(\rho,u,p)=(1,0,100)$ otherwise with reflecting boundary conditions on both sides. The computed solution for density obtained at $t=0.038$s using 200 elements while using the CSWENO limiter for $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM is compared and plotted against the solution obtained using the parent WENO limiter, simple WENO limiter and the exact solution in Figures 7a and 7b. From the figures, for $\mathbf{P}^{2}$ based DGM, we can see that the performance of the CSWENO limiter and the simple WENO limiter is quite similar and both perform better than the WENO limiter. For $\mathbf{P}^{3}$ based DGM, the performance of all three limiters seems to be similar.

Example 11: As a test problem for the two-dimensional case, we solve the double Mach reflection problem which is given in [12]. We solve the two-dimensional Euler equations (21) in the computational domain $[0,4]\times[0,1]$. Initially, right moving Mach 10 shock is positioned at $x=1/6,y=0$ and it makes an angle $60^{0}$ with the $x$-axis. For the bottom boundary, we impose the exact post shock conditions from $x=0$ to $x=1/6$ and for the rest of the $x$-axis, we use reflective boundary conditions. For the top boundary, we set conditions to describe the exact motion of a Mach 10 shock. We compute the solution upto time $t=0.2$ for two different uniform meshes with $960\times 240$ and $1920\times 480$ cells in each mesh. The full solution using the CSWENO limiter for the most refined mesh (containing $1920\times 480$ cells) for $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM has been shown in Figure 8. A zoom-in view of the density contours near the double Mach stem has been shown in Figures 9, 10 and 11 using the CSWENO limiter, the parent WENO limiter and the SWENO limiter respectively for the two different mesh sizes. We have also added a positivity preserving component to the described limiter as given in [14].

We can see that the solution for all three limiters is quite well comparable to the solution obtained in [12].

5 Conclusions:

We have developed a different WENO limiting strategy for the solution of hyperbolic conservation laws using Discontinuous Galerkin method based on the limiter developed by Qiu and Shu [4]. Here, we identify the troubled cells and use only the immediate neighbors by dividing them into subcells. This is different from the subcell limiting strategy of Dumbser et al [3] which is much more accurate but quite complicated. These new cells are used for the reconstruction of the WENO polynomial. We termed this limiting procedure as the compact subcell WENO limiter (CSWENO limiter). We have tested the accuracy of this limiter using various standard test cases containing smooth solutions and calculating the numerical order of accuracy. We have also provided numerical results with shocks and compared the results obtained using this limiter with that obtained from the parent WENO limiter [4] and the simple WENO limiter proposed by Zhong and Shu [2]. We can conclude from the results that the performance of the simple WENO limiter is better in most cases than the CSWENO limiter. Also, the CSWENO limiter performs better than the parent WENO limiter for most of the examples discussed. Also, the CSWENO limiter uses only the cell averages for the WENO reconstruction and uses a very compact stencil like the simple WENO limiter which is highly beneficial near the boundaries. Implementation of this limiter for unstructured meshes is ongoing.