# A compact subcell WENO limiting strategy using immediate neighbors for
Runge-Kutta discontinuous Galerkin methods for unstructured meshes

###### Abstract

In this paper, we generalize the compact subcell weighted essentially non
oscillatory (CSWENO) limiting strategy for Runge-Kutta discontinuous Galerkin
method developed recently in [1] for structured meshes to unstructured
triangular meshes. The main idea of the limiting strategy is to divide the
immediate neighbors of a given cell into the required stencil and to use a
WENO reconstruction for limiting. This strategy can be applied for any type of
WENO reconstruction. We used the WENO reconstruction proposed in [2] and
provided accuracy tests and results for two-dimensional Burgers’ equation and
two dimensional Euler equations to illustrate the performance of this limiting
strategy.

Keywords: discontinuous Galerkin method, WENO limiter, unstructured meshes

##  1 Introduction

In this paper, we look at the solution of hyperbolic conservation laws with
the Runge-Kutta discontinuous Galerkin (RKDG) method [3] on unstructured
meshes. To control spurious oscillations near discontinuities, a limiter is
used. For higher orders, weighted essentially non oscillatory (WENO) limiters
are preferred as they maintain the order of the scheme. We generalize the WENO
limiting strategy developed recently in [1] for structured meshes to
unstructured triangular meshes. This strategy works with any type of WENO
reconstruction in the target cell using only the immediate neighbors by
dividing them into subcells to obtain the required stencil. We used the WENO
reconstruction proposed in [2] as it is quite simple for implementation.

WENO limiting for RKDG method was first presented by Qiu and Shu for
structured meshes in [4]. Their limiter was extended to unstructured meshes by
Zhu et al. in [5]. This limiter requires neighbors of neighbors for limiting
in a given cell for higher orders. Various other limiters given in [6], [7],
[8] addressed this issue of using a wide stencil and used the polynomial in a
given cell along with the polynomials of the immediate neighbors to obtain a
limited polynomial using WENO reconstruction. On the other hand, Dumbser et al
[9] used a different strategy where the target cell is divided into subcells
and an a posteriori limiting strategy is used. This was further refined using
an Adaptive Mesh Refinement (AMR) technique by Giri and Qiu in [10].

We extend the compact limiting strategy used in [1] to unstructured meshes.
Here, we use the immediate neighbors for limiting in a given cell and these
neighbors are divided into subcells for higher order limiting. We assign
appropriate values to these subcells as explained in Section 3. This limiting
strategy can be used with any type of WENO reconstructions as given in [11]
(called type-I WENO reconstruction), or [5] (called type-II WENO
reconstruction) or [12] (mixed reconstruction) or the more recent methods
given in [13] and [2]. WENO schemes form a weighted combination of several
local reconstructions based on different stencils (called small stencils) and
use it as the final WENO reconstruction. Type-I reconstruction consists of
WENO schemes whose order of accuracy is not higher than that of the
reconstruction on each small stencil. These schemes require very wide stencils
for higher orders. On the other hand, type-II WENO schemes require smaller
stencils whose order of accuracy is higher than that of the reconstruction on
each small stencil. Each of these reconstructions have their pros and cons and
they are discussed in detail in [12]. We use the WENO reconstruction given by
Zhu and Shu in [2] as it is quite simple in implementation and extension to
higher orders is easy. These [2] schemes have linear weights that can be any
positive numbers on the condition that their sum is one and number of spatial
stencils is smaller than that of the same order accurate classical finite
volume WENO schemes [14] on triangular meshes. We emphasize that our procedure
works with any of the WENO reconstructions listed above by dividing the
neighbors of the cell in which limiting is needed as is required by the
reconstruction. We call this 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 the testing of the
limiter and the results are described in Section 4 and finally we conclude the
paper in Section 5.

##  2 Formulation of discontinuous Galerkin method on unstructured meshes

Consider a two-dimensional conservation law of the form for variable $u$

$\frac{\partial u}{\partial t}+\frac{\partial f(u)}{\partial x}+\frac{\partial
g(u)}{\partial y}=0,\qquad\qquad x,y\in\mathbf{D}$ (1)

with the initial condition at $t=0$ given by,

$u(x,y,0)=u_{0}(x,y),$

and the appropriate boundary conditions on the boundary $\partial\mathbf{D}$.
Here $f(u)$ and $g(u)$ are fluxes in $x$ and $y$ directions respectively.

We assume that the domain $\mathbf{D}$ can be triangulated by $K$ elements as

$\mathbf{D}=\bigcup_{k=1}^{K}\mathbf{I}_{k}$ (2)

where $\mathbf{I}_{k}$ is a straight sided triangle and the triangulation is
assumed to be geometrically conforming i.e., $\partial\mathbf{D}$ is
approximated by a piece wise linear polygon with each line a face of the
triangle.

We approximate the local solution as a polynomial of order $N$ as given by:

$u_{h}^{k}(x,y,t)=\sum_{n=0}^{N_{p}-1}\hat{u}^{k}_{n}(t)\psi_{n}^{k}(x,y)\qquad\forall
x,y\in\mathbf{I}_{k}$ (3)

where $\psi_{n}^{k}(x,y)$ is a two dimensional polynomial basis of order $N$,
$N_{p}$ is the number of degrees of freedom and $h$ is the size of the grid
(which is in general the average of the lengths of the sides of the triangle).
This is termed to be $\mathbf{P}^{N}$ based discontinuous Galerkin method
whose formal order of accuracy is $N+1$. The number of degrees of freedom
$N_{p}$ is given by:

$N_{p}=\frac{(N+1)(N+2)}{2}$

We use the orthonormal basis as given in [15] for a standard isoparametric
triangle $I$ such that

$I=\\{\vec{r}=(r,s)|(r,s)\geq-1;r+s\leq 0\\}$ (4)

The basis is given as:

$\psi_{m}(\vec{r})=\sqrt{2}\left[P_{i}(a)P_{j}^{(2i+1,0)}(b)\right](1-b)^{i}$
(5)

where

$m=j+(N+1)i+1-\frac{i}{2}(i-1),\qquad(i,j)\geq 0;\qquad i+j\leq N$

and

$a=2\frac{1+r}{1-s}-1;\qquad b=s$

and $P_{n}^{(\alpha,\beta)}(x)$ is the $n$th order Jacobi polynomial.

We approximate the fluxes $f(u(x,y,t),x,y,t)$ and $g(u(x,y,t),x,y,t)$ in the
domain $\mathbf{D}$ as

$f_{h}^{k}(u_{h}^{k})=\sum_{n=0}^{N_{p}-1}\hat{f}^{k}_{n}(t)\psi_{n}^{k}(x,y)\qquad\forall
x,y\in\mathbf{I}_{k}$ (6)

$g_{h}^{k}(u_{h}^{k})=\sum_{n=0}^{N_{p}-1}\hat{g}^{k}_{n}(t)\psi_{n}^{k}(x,y)\qquad\forall
x,y\in\mathbf{I}_{k}$ (7)

where $\hat{f}^{k}_{n}(t)=f(\hat{u}^{k}_{n}(t))$ and
$\hat{g}^{k}_{n}(t)=g(\hat{u}^{k}_{n}(t))$. Substituting equations (3), (6)
and (7) in (1) and integrating it by parts, we get the following scheme to
advance the degrees of freedom $\hat{u}^{k}_{n}(t)$ in time in an element
given by $\mathbf{I}_{k}=[x_{l}^{k},y_{l}^{k}]\times[x_{r}^{k},y_{r}^{k}]$:

$\displaystyle\frac{d}{dt}\hat{u}_{h}^{k}=(\mathbf{M}^{k})^{-1}(\mathbf{S_{x}}^{k})^{T}\hat{f}_{h}^{k}(\hat{u}_{h}^{k})+(\mathbf{M}^{k})^{-1}(\mathbf{S_{y}}^{k})^{T}\hat{g}_{h}^{k}(\hat{u}_{h}^{k})$
(8)
$\displaystyle-(\mathbf{M}^{k})^{-1}(f^{*}|_{r_{N_{p}}}e_{N_{p}}-f^{*}|_{r_{1}}e_{1})-(\mathbf{M}^{k})^{-1}(g^{*}|_{s_{N_{p}}}e_{N_{p}}-g^{*}|_{s_{1}}e_{1})=L(u_{h}^{k})$

Here, $\hat{u}_{h}^{k}=[\hat{u}_{0}^{k}\ldots\hat{u}_{N_{p}-1}^{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:

$\mathbf{M}^{k}=\left[M_{ij}^{k}\right]=\left[\int_{\mathbf{I}_{k}}\psi_{i}^{k}(x,y)\psi_{j}^{k}(x,y)~{}\text{dA}\right]$
(9)

and $\mathbf{S_{x}}^{k}$ and $\mathbf{S_{y}}^{k}$ are the local stiffness
matrices given by:

$\mathbf{S_{x}}^{k}=\left[S_{xij}^{k}\right]=\left[\int_{\mathbf{I}_{k}}\psi_{i}^{k}(x,y)\frac{\partial\psi_{j}^{k}(x,y)}{\partial
x}~{}\text{dA}\right]$ (10)

$\mathbf{S_{y}}^{k}=\left[S_{yij}^{k}\right]=\left[\int_{\mathbf{I}_{k}}\psi_{i}^{k}(x,y)\frac{\partial\psi_{j}^{k}(x,y)}{\partial
y}~{}\text{dA}\right]$ (11)

Also, $f^{*}$ and $g^{*}$ are the $x$ and $y$ components of 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 on unstructured meshes

In this section, we describe the details of the proposed limiting procedure
using WENO reconstruction and the division of cells into subcells for the
discontinuous Galerkin method. This step is a generalization of the procedure
in [1] for structured meshes. The common method for limiting in discontinuous
Galerkin method is:
1) Identify the cells which need to be limited, 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 [16] for
all the calculations done in this paper as it is rated highly by Qiu and Shu
in [17] on the basis of it’s performance in detecting the discontinuities in
various test problems. We give a brief description of the troubled cell
indicator in the following subsection.

###  3.1 KXRCF Troubled cell indicator

The troubled cell indicator used in this work was developed by Krivodonova et
al. [16]. This indicator is termed as KXRCF troubled cell indicator by using
the names of the authors of the paper [16].

For a given problem, partition the boundary $\partial\mathbf{I}_{j}$ of a
given cell $\mathbf{I}_{j}$ into portions $\partial\mathbf{I}_{j}^{-}$ and
$\partial\mathbf{I}_{j}^{+}$ where the flow is into ($\bar{v}.\bar{n}<0$) and
out of ($\bar{v}.\bar{n}>0$) $\mathbf{I}_{j}$, respectively. For the scalar
conservation laws, $\bar{v}$ is taken to be $(f^{\prime}(u),g^{\prime}(u))$
and for the Euler equations it is the velocity vector. The troubled cell
indicator is defined as

$R_{j}=\frac{|\int_{\partial\mathbf{I}_{j}^{-}}(U_{j}-U_{nbj})ds|}{h^{(N+1)/2}|\partial\mathbf{I}_{j}^{-}|\quad||U_{j}||}$
(12)

where $U_{j}$ is the discontinuous Galerkin value of $u$ on $\mathbf{I}_{j}$,
$U_{nbj}$ the value across the boundary, $||U_{j}||$ the standard $L^{2}$ norm
in the cell $\mathbf{I}_{j}$ and $N$ is the order of the polynomial basis. We
take $h$ to be the radius of the circumscribed circle in $\mathbf{I}_{j}$. Now
if $R_{j}>C_{K}$, $\mathbf{I}_{j}$ is identified to be a troubled cell. We
take $C_{K}$ to be 1.

###  3.2 Limiting procedure

After identifying the troubled-cells, we would like to reconstruct the values
of $\hat{u}_{n}^{j}$ for the troubled-cell $\mathbf{I}_{j}$ for
$n=1,\ldots,N_{p}-1$. That is, we retain the cell average $\hat{u}_{0}^{j}$
and reconstruct all the other degrees of freedom. To do that, we use the
quadrature points given in [18] (SCP-quadrature) and find the solution at
those points using WENO reconstruction. These quadrature points are termed as
$(x_{q},y_{q})$ for $q=0,\ldots,N_{p}-1$. We can use any of the WENO
reconstructions given in [11] (called type I WENO reconstruction), or [5]
(called type II WENO reconstruction) or [12] (mixed reconstruction) or the
more recent methods given in [13] and [2]. WENO schemes form a weighted
combination of several local reconstructions based on different stencils
(called small stencils) and use it as the final WENO reconstruction. Type-I
reconstruction consists of WENO schemes whose order of accuracy is not higher
than that of the reconstruction on each small stencil. These schemes require
very wide stencils for higher orders. On the other hand, type-II WENO schemes
require smaller stencils whose order of accuracy is higher than that of the
reconstruction on each small stencil. But, they sometimes require negative
linear weights and do not work properly for geometries of poor mesh quality.
Each of these reconstructions have their pros and cons and they are discussed
in detail in [12]. We use the WENO reconstruction given by Zhu and Shu in [2]
as it is quite simple in implementation and extension to higher orders is
easy. These [2] schemes have linear weights that can be any positive numbers
on the condition that their sum is one and number of spatial stencils is
smaller than that of the same order accurate classical finite volume WENO
schemes [14] on triangular meshes. We emphasize that our procedure works with
any of the WENO reconstructions listed above by dividing the neighbors of the
cell in which limiting is needed as is required by the reconstruction. For the
$\mathbf{P}^{1}$ based DGM, we now describe the procedure for the
reconstruction of the moments $\hat{u}_{1}^{k}$ and $\hat{u}_{2}^{k}$ in the
troubled cell $\mathbf{I}_{j}$ using the WENO reconstruction procedure for
triangles. We relabel the troubled cell as $\Delta_{0}$ and its immediate
neighbors as $\Delta_{1},\Delta_{2}$ and $\Delta_{3}$ as shown in Figure 1
where they are labeled according to their subscripts.

Figure 1: Troubled Cell and its immediate neighbors

For a third order WENO reconstruction, we also need the neighbors of the cells
$\Delta_{1},\Delta_{2}$ and $\Delta_{3}$ as suggested in [5]. Instead of using
their actual neighbors, we divide each of those triangles into three equal
parts by joining the centroid of a triangle to its vertices as shown in Figure
2a. Now, we get the new stencil required for the third order WENO
reconstruction. Now, we use the DG local polynomial for $\Delta_{1}$ and take
its average in the new cells $\Delta_{11}$, $\Delta_{12}$ and $\Delta_{13}$ to
obtain the cell averages for the new cells. Similarly, we can calculate the
cell average for the cells $\Delta_{21}$, $\Delta_{22}$, $\Delta_{23}$,
$\Delta_{31}$, $\Delta_{32}$ and $\Delta_{33}$ using the DG local polynomial
of cells $\Delta_{2}$ and $\Delta_{3}$. 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. As the order of the polynomial
increases, the number of triangles required to maintain the order of the
scheme increases. We split the neighbors in such a way that this will not
reduce the order of accuracy as explained in [14].

(a) Division of cell neighbors for $\mathbf{P}^{1}$ based DGM

(b) Division of cell neighbors for $\mathbf{P}^{2}$ based DGM

Figure 2: Division of cell neighbors for WENO reconstruction using
$\mathbf{P}^{1}$ and $\mathbf{P}^{2}$ based DGM Figure 3: Required stencils
for $\mathbf{P}^{1}$ based DGM

Now the step by step WENO reconstruction procedure to be used for
$\mathbf{P}^{1}$ based DGM is described below:

Step 1: We use stencils T1, T2 and T3 as shown in Figure 3 for WENO
reconstruction. For T1 = {$\Delta_{0}$}, we assign a zeroth order polynomial
$q_{1}(x,y)=\hat{u}_{0}^{j}$. Now on T2 =
{$\Delta_{0}$,$\Delta_{11}$,$\Delta_{21}$,$\Delta_{31}$}, we obtain a linear
polynomial $q_{2}(x,y)$ such that it has the same cell average as $u$ on the
troubled cell $\Delta_{0}$ and matches the cell averages of $u$ (obtained
earlier) on the other triangular cells in the set T2$\backslash${$\Delta_{0}$}
in a least square sense [14]. To obtain a quadratic polynomial $q_{3}(x,y)$,
we use the stencil T3 =
{$\Delta_{0}$,$\Delta_{11}$,$\Delta_{21}$,$\Delta_{31}$,$\Delta_{12}$,$\Delta_{13}$,$\Delta_{22}$,$\Delta_{23}$,$\Delta_{32}$,$\Delta_{33}$}
such that it has the same cell average as $u$ on the troubled cell
$\Delta_{0}$ and matches the cell averages of $u$ on the other triangular
cells in the set T3$\backslash${$\Delta_{0}$} in a least square sense [14].

Step 2: Now, we obtain equivalent expressions for these polynomials
$q_{1}(x,y)$, $q_{2}(x,y)$ and $q_{3}(x,y)$. We define:

$p_{1}(x,y)=q_{1}(x,y)$ (13)

$p_{2}(x,y)=\frac{1}{\gamma_{2,2}}q_{2}(x,y)-\frac{\gamma_{1,2}}{\gamma_{2,2}}p_{1}(x,y)\quad\text{with}\quad\sum_{l=1}^{2}\gamma_{l,2}=1\quad\text{and}\quad\gamma_{2,2}\neq
0$ (14)

and

$p_{3}(x,y)=\frac{1}{\gamma_{3,3}}q_{3}(x,y)-\sum_{l=1}^{2}\frac{\gamma_{l,3}}{\gamma_{3,3}}p_{l}(x,y)\quad\text{with}\quad\sum_{l=1}^{3}\gamma_{l,3}=1\quad\text{and}\quad\gamma_{3,3}\neq
0$ (15)

Here the $\gamma_{l,m}$ are linear weights. We choose them as in [2] as
$\gamma_{l,m}=\bar{\gamma}_{l,m}/\left({\sum_{l=1}^{m}\bar{\gamma}_{l,m}}\right)$
with $\bar{\gamma}_{l,m}=(1/h)^{l-1}$, $l=1,\ldots,t$ and $m=2,3$ ($t=m$)
where $h$ is the grid size (which is in general the average of the lengths of
the sides of the triangle). Here we changed the expression for
$\bar{\gamma}_{l,m}$ from $\bar{\gamma}_{l,m}=10^{l-1}$ (which is used in [2])
to $\bar{\gamma}_{l,m}=(1/h)^{l-1}$. This change was made so that
$\gamma_{l,m}$ = $O(h^{m-l})$. In this manner, we get larger linear weights
for high degree polynomials and smaller linear weights for low degree
polynomials. Also, using a Taylor series expansion of cell averages about the
barycenter of the triangle $\Delta_{0}$ as given in [19], we obtain
$u(x_{q},y_{q})-p_{l}(x_{q},y_{q})=O(h^{l})$, where $u(x_{q},y_{q})$ is the
exact solution at the quadrature point $q$.

Step 3: We now compute the smoothness indicator for each stencil denoted by
$\beta_{l}$ for $l=2,3$. As given by [20], we use:

$\beta_{l}=\sum_{|l|=1}^{N}|\Delta_{0}|^{|l|-1}\int_{\Delta_{0}}\left(\frac{\partial^{|l|}}{\partial
x^{l_{1}}\partial x^{l_{2}}}p_{i}(x,y)\right)^{2}dxdy,\quad l=2,3$ (16)

where $l=(l_{1},l_{2})$, $|l|=l1+l2$. For $\beta_{1}$, we magnify it from 0 to
a value defined in equation (18) as given by [2]. Using the stencils in Figure
3, we construct three polynomials
$p_{1,l}(x,y)\in\text{span}\\{\frac{x-x_{l}}{|\Delta_{l}|^{1/2}}\\}$,
$l=1,2,3$, satisfying $p_{1,1}(x_{12},y_{12})=\bar{u}_{12}-\bar{u}_{11}$,
$p_{1,1}(x_{13},y_{13})=\bar{u}_{13}-\bar{u}_{11}$,
$p_{1,2}(x_{22},y_{22})=\bar{u}_{22}-\bar{u}_{21}$,
$p_{1,2}(x_{23},y_{23})=\bar{u}_{23}-\bar{u}_{21}$,
$p_{1,3}(x_{32},y_{32})=\bar{u}_{32}-\bar{u}_{31}$,
$p_{1,3}(x_{33},y_{33})=\bar{u}_{33}-\bar{u}_{31}$, where $(x_{l},y_{l})$ are
centroids of $\Delta_{l}$, $l=11,12,13;21,22,23;31,32,33$, respectively. Now,
we use equation (16) to find $\beta_{1,l}$, $l=1,2,3$. We set
$\lambda_{1,1}=\lambda_{1,2}=\lambda_{1,3}=1/3$ and

$\sigma_{l}=\lambda_{1,l}\left(1+\frac{\left(\left[|\beta_{1,1}-\beta_{1,2}|^{2}+|\beta_{1,2}-\beta_{1,3}|^{2}+|\beta_{1,3}-\beta_{1,1}|^{2}\right]/3\right)^{2}}{\beta_{1,l}+\epsilon}\right),l=1,2,3$
(17)

with $\sigma=\sigma_{1}+\sigma_{2}+\sigma_{3}$, where $\epsilon=10^{-10}$.
Then we get

$\beta_{1}=\sum_{|l|=1}|\Delta_{0}|^{|l|}\left(\frac{\partial^{|l|}}{\partial
x^{l_{1}}\partial
x^{l_{2}}}\left(\sum_{l=1}^{3}\frac{\sigma_{l}}{\sigma}p_{1,l}(x,y)\right)\right)^{2}$
(18)

Again using a Taylor series expansion of cell averages about the barycenter of
the triangle $\Delta_{0}$, we obtain for $l=1,2,3$,
$\beta_{l}=d(x_{l},y_{l})h^{2}(1+O(h^{2}))$ for some function $d$ where
$(x_{l},y_{l})$ is the location of barycenter of $\mathbf{I}_{j}$. If a
particular stencil, say T2 is not smooth, then $\beta_{2}=O(1)$.

Step 4: We compute the nonlinear weights using the WENO-Z recipe [21]. We
first find $\tau$ using

$\tau=\frac{|\beta_{3}-\beta_{1}|^{2}+|\beta_{3}-\beta_{2}|^{2}}{4}$ (19)

Here, using the Taylor series expansions about the barycenter of the triangle
$\Delta_{0}$, we obtain $\tau=O(h^{6})$. Now, the nonlinear weights are given
by

$\omega_{l}=\frac{\bar{\omega}_{l}}{\sum_{l=1}^{3}\bar{\omega}_{l}},\quad\bar{\omega}_{l}=\gamma_{l,3}\left(1+\frac{\tau}{\epsilon+\beta_{l}}\right)\quad
l=1,2,3$ (20)

Again we take $\epsilon$ as $10^{-10}$. Now assuming $\epsilon<<\beta_{l}$, we
get $\omega_{1}-\gamma_{1,3}=O(h^{6})$, $\omega_{2}-\gamma_{2,3}=O(h^{5})$,
and $\omega_{3}-\gamma_{3,3}=O(h^{4})$. If a particular stencil, say T2 is not
smooth, then $\omega_{2}-\gamma_{l,2}=O(1)$. We note that the stencil T1 is
always smooth as it contains only one cell with a given cell average. Now, we
can write the final WENO approximation at the required quadrature point $q$ as

$u_{q}=\sum_{l=1}^{3}\omega_{l}p_{l}(x_{q},y_{q})\quad q=0,\ldots,N_{p}-1$
(21)

Now, we write

$u_{q}-u(x_{q},y_{q})=\omega_{0}p_{0}(x_{q},y_{q})+\omega_{1}p_{1}(x_{q},y_{q})+\omega_{2}p_{2}(x_{q},y_{q})-u(x_{q},y_{q})$
(22)

which can written as

$\displaystyle
u_{q}-u(x_{q},y_{q})=[(\omega_{3}-\gamma_{3,3})(p_{3}(x_{q},y_{q})-u(x_{q},y_{q}))]+[(\omega_{2}-\gamma_{2,3})(p_{2}(x_{q},y_{q})-u(x_{q},y_{q}))]$
$\displaystyle+[(\omega_{1}-\gamma_{1,3})(p_{1}(x_{q},y_{q})-u(x_{q},y_{q}))]+[\gamma_{3,3}(p_{3}(x_{q},y_{q})-u(x_{q},y_{q}))]$
$\displaystyle+[\gamma_{2,3}(p_{2}(x_{q},y_{q})-u(x_{q},y_{q}))]+[\gamma_{1,3}(p_{3}(x_{q},y_{q})-u(x_{q},y_{q}))]$
(23)

If all the stencils are smooth, from the right hand side of equation (3.2), we
can see that $u(x_{q},y_{q})-u_{q}=O(h^{3})$. If any one of the stencils (say
T2) is not smooth, then $\omega_{2}-\gamma_{2,3}=O(1)$ and the right hand side
of equation (3.2) will become $O(h^{2})$. This limit is the required order of
accuracy for the reconstruction of moments for $\mathbf{P^{1}}$ based DGM.

Step 5: Finally, we obtain the reconstructed degrees of freedom based on the
reconstructed point values $u_{q}$ at the quadrature points $(x_{q},y_{q})$
and a numerical integration as

$\hat{u}_{i}^{j}=|\Delta_{0}|\sum_{q}w_{q}u_{q}\psi_{i}^{j}(x_{q},y_{q})\quad
i=1,\ldots,N_{p}-1$ (24)

where $|\Delta_{0}|$ is the area of $\Delta_{0}$ and $w_{G}$’s are the
quadrature weights for the points $(x_{q},y_{q})$ as given by [18].

For the $\mathbf{P}^{2}$ based DGM, the procedure to construct the second
order moments $\hat{u}_{1}^{k}$, $\hat{u}_{2}^{k}$, $\hat{u}_{3}^{k}$,
$\hat{u}_{4}^{k}$ and $\hat{u}_{5}^{k}$ in the troubled cell $\Delta_{0}$ is
similar to the above procedure. Here, we require some extra neighbors for the
WENO reconstruction as we need a fifth order WENO reconstruction and they are
shown in Figure 2b. Here note that the new triangles $221,231$ etc., are
obtained by dropping a perpendicular to the opposite side as required. The
cell averages for the new cells are obtained using the same procedure
described above for $\mathbf{P}^{1}$ based DGM. The WENO reconstruction
procedure for this case is described below:

Figure 4: Required stencils for $\mathbf{P}^{2}$ based DGM

Step 1: We use stencils T1, T2, T3 and T4 as shown in Figure 4 for WENO
reconstruction. We obtain $q_{1}(x,y)$, $q_{2}(x,y)$ and $q_{3}(x,y)$ using
stencils T1, T2 and T3 as described in the procedure above. Now, using T4 =
{$\Delta_{0}$,$\Delta_{11}$,$\Delta_{21}$,$\Delta_{31}$,$\Delta_{12}$,$\Delta_{13}$,$\Delta_{121}$,$\Delta_{131}$,$\Delta_{22}$,$\Delta_{23}$,$\Delta_{221}$,$\Delta_{231}$,$\Delta_{32}$,$\Delta_{33}$,$\Delta_{321}$,$\Delta_{331}$},
we obtain a cubic polynomial $q_{4}(x,y)$, such that it has the same cell
average as $u$ on the troubled cell $\Delta_{0}$ and matches the cell averages
of $u$ on the other triangular cells in the set T4$\backslash${$\Delta_{0}$}
in a least square sense. Again a quartic polynomial $q_{5}(x,y)$ is obtained
on T5 such that it has the same cell average as $u$ on the troubled cell
$\Delta_{0}$ and matches the cell averages of $u$ on the other triangular
cells in the set T4$\backslash${$\Delta_{0}$} in a least square sense.

Step 2: We obtain $p_{1}(x,y)$, $p_{2}(x,y)$ and $p_{3}(x,y)$ using equations
(13), (14) and (15) respectively. We now define:

$p_{4}(x,y)=\frac{1}{\gamma_{3,4}}q_{4}(x,y)-\sum_{l=1}^{2}\frac{\gamma_{l,4}}{\gamma_{3,4}}p_{l}(x,y)\quad\text{with}\quad\sum_{l=1}^{3}\gamma_{l,4}=1\quad\text{and}\quad\gamma_{3,4}\neq
0$ (25)

$p_{5}(x,y)=\frac{1}{\gamma_{4,5}}q_{5}(x,y)-\sum_{l=1}^{3}\frac{\gamma_{l,5}}{\gamma_{4,5}}p_{l}(x,y)\quad\text{with}\quad\sum_{l=1}^{4}\gamma_{l,5}=1\quad\text{and}\quad\gamma_{4,5}\neq
0$ (26)

We choose $\gamma_{l,m}$ again as given above for $m=2,3$. For $m=4,5$, we set
$\gamma_{l,m}=\bar{\gamma}_{l,m}/\left(\sum_{l=1}^{m}\bar{\gamma}_{l,m}\right)$
with $\bar{\gamma}_{l,m}=(1/h)^{l-1}$, $l=1,\ldots,t$ and $m=4,5$ ($t=m-1$)
where $h$ is the grid size.

Step 3: Using the same procedure given in Step 3 for $\mathbf{P}^{1}$ based
DGM and equations (16), (18), we obtain $\beta_{l}$, for $l=1,2,3,4,5$.

Step 4: We find $\tau$ required for calculating nonlinear weights using

$\tau=\frac{|\beta_{4}-\beta_{1}|^{2}+|\beta_{4}-\beta_{2}|^{2}+|\beta_{4}-\beta_{3}|^{2}}{16}$
(27)

Now, we calculate the nonlinear weights using

$\omega_{l}=\frac{\bar{\omega}_{l}}{\sum_{l=1}^{4}\bar{\omega}_{l}},\quad\bar{\omega}_{l}=\gamma_{l,5}\left(1+\frac{\tau}{\epsilon+\beta_{l}}\right)\quad
l=1,2,3,4$ (28)

Again we take $\epsilon$ as $10^{-10}$. Now, we can write the final WENO
approximation at the required quadrature point $q$ as

$u_{q}=\sum_{l=1}^{4}\omega_{l}p_{l}(x_{q},y_{q})\quad q=0,\ldots,N_{p}-1$
(29)

Using equations (29) and (3.2), we can say that when the solution of smooth,
$u(x_{q},y_{q})-u_{q}=O(h^{5})$. If one of the stencils is not smooth, we have
$u(x_{q},y_{q})-u_{q}=O(h^{3})$ . This limit is the required order of accuracy
for the reconstruction of moments for $\mathbf{P^{2}}$ based DGM. We can
reconstruct the degrees of freedom in the same way by dividing the neighbors
and using WENO reconstruction for $\mathbf{P^{3}}$ based DGM.

Now, equation (24) gives us the degrees of freedom for $\mathbf{P}^{2}$ based
DGM. For a non-orthonormal basis, we define

$D_{i}^{j}=|\Delta_{0}|\sum_{q}w_{q}u_{q}\psi_{i}^{j}((x_{q},y_{q}))\quad
i=1,\ldots,N_{p}-1$ (30)

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

$\mathbf{A}^{j}=\left[A\right]=\begin{bmatrix}\left[M_{11}^{j}\right]&\ldots&\left[M_{1N_{p}-1}^{j}\right]\\\
\ldots&\ldots&\ldots\\\
\rule{0.0pt}{10.76385pt}\left[M_{N_{p}-11}^{j}\right]&\ldots&\rule{0.0pt}{10.76385pt}\left[M_{N_{p}-1N_{p}-1}^{j}\right]\end{bmatrix}$
(31)

Here, the terms $\left[M_{mn}^{j}\right]$ are given by equation (9). Then the
reconstructed degrees of freedom are given by
$\mathbf{X}^{j}=(\mathbf{A}^{j})^{-1}\mathbf{B}^{j}$. This formulation will
work for any polynomial basis. Now, we can get the reconstructed polynomial
solution in $\mathbf{I}_{j}$ by equation (3). We call this limiting procedure
the compact subcell WENO limiting or CSWENO limiting in short. For solving a
system of equations, we use this with a local characteristic field
decomposition with the corresponding Jacobians in the $x$ and $y$ directions
as explained in [6].

Now, the semi-discrete scheme given in equation (8) along with the limiter is
discretized in time by using the TVD Runge-Kutta time discretization
introduced in [22] which is described below briefly.

If $\\{t^{m}\\}_{m=0}^{M}$ is a partition of $[0,T]$ and $\Delta t^{m}$ =
$t^{m+1}-t^{m}$, $m=0,\ldots,M-1$, the time-marching algorithm reads as
follows:

1) Set $u_{h}^{0}$ = $u_{0h}$, the initial condition;

2) For $m=0,\ldots,M-1$ compute $u_{h}^{m+1}$ from $u_{h}^{m}$ as follows:
(i) set $u_{h}^{(0)}$ = $u_{h}^{m}$;
(ii) for $i=1,2,\ldots,d$ compute the intermediate functions:

$u_{h}^{(i)}=\left[\sum_{l=0}^{i-1}\alpha_{il}u_{h}^{(l)}+\beta_{il}\Delta
t^{m}L_{h}(u_{h}^{(l)})\right];$ (32)

where $d$ is the order of the time integration.

(iii) set $u_{h}^{m+1}$ = $u_{h}^{(d)}$;

Some Runge-Kutta time discretization parameters are given in Table 1.

Runge-Kutta discretization parameters
---
order |  $\alpha_{il}$ |  $\beta_{il}$ |  max$\\{\alpha_{il}/\beta_{il}\\}$
2 |  1 |  1 |  1
|  $\frac{1}{2}$ $\frac{1}{2}$ |  $0$ $\frac{1}{2}$ |
3 |  1 |  1 |  1
|  $\frac{3}{4}$ $\frac{1}{4}$ |  0 $\frac{1}{4}$ |
|  $\frac{1}{3}$ 0 $\frac{2}{3}$ |  0 0 $\frac{2}{3}$ |
Table 1: TVD Runge-Kutta discretization parameters for orders 2 and 3

We have used the 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 in Section 3. We used Gmsh 4.6.0 software [23] for the
generation of meshes for all our calculations. 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 the two-dimensional test cases. We present the results of
the accuracy tests using two-dimensional Burgers equations and the two-
dimensional Euler equations. For all the accuracy tests conducted, we have
marked all the cells as troubled cells. We used meshes which contain some
triangles of irregular size as shown in Figure 5 for all our calculations to
illustrate that the limiter retains the order of the scheme even for such
meshes.

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

$\frac{\partial u}{\partial t}+\frac{\partial(u^{2}/2)}{\partial
x}+\frac{\partial(u^{2}/2)}{\partial
y}=0,\qquad\qquad(x,y)\in[-2,2]\times[-2,2]$ (33)

with the initial condition $u(x,y,0)=0.5+\sin(\pi(x+y)/2)$, with periodic
boundary conditions in both directions. The exact solution is smooth till
$t=0.5/\pi$. A sample mesh used is shown in Figure 5. The errors and numerical
orders of accuracy are calculated at $t=0.5/\pi$ by marking all the cells as
troubled cells and are presented in Table 2. We can see that the CSWENO
limiter maintains the order and magnitude of accuracy of the original DG
method.

|  | DG without limiter | DG with limiter
---|---|---|---
$h$ |  $L_{1}$ error | Order |  $L_{\infty}$ error | Order |  $L_{1}$ error | Order |  $L_{\infty}$ error | Order
$\mathbf{P}^{1}$ | 4/20 | 2.31E-02 |  | 4.14E-01 |  | 6.64E-03 |  | 8.35E-02 |
4/40 | 4.89E-03 | 2.24 | 8.47E-02 | 2.29 | 1.70E-03 | 1.97 | 2.21E-02 | 1.92
4/80 | 1.01E-03 | 2.27 | 1.84E-02 | 2.20 | 4.37E-04 | 1.96 | 6.30E-03 | 1.81
4/160 | 2.34E-04 | 2.11 | 4.76E-03 | 1.95 | 1.12E-04 | 1.96 | 1.82E-03 | 1.79
$\mathbf{P}^{2}$ | 4/20 | 4.17E-04 |  | 9.39E-03 |  | 3.97E-04 |  | 9.66E-03 |
4/40 | 5.47E-05 | 2.93 | 1.34E-03 | 2.81 | 5.28E-05 | 2.91 | 1.82E-03 | 2.41
4/80 | 7.59E-06 | 2.85 | 2.08E-04 | 2.69 | 7.32E-06 | 2.85 | 2.84E-04 | 2.68
4/160 | 1.09E-06 | 2.80 | 4.14E-05 | 2.33 | 9.81E-07 | 2.90 | 4.02E-05 | 2.82
$\mathbf{P}^{3}$ | 4/20 | 3.69E-05 |  | 2.27E-03 |  | 3.87E-05 |  | 8.23E-04 |
4/40 | 2.51E-06 | 3.88 | 1.63E-04 | 3.80 | 2.50E-06 | 3.95 | 6.08E-05 | 3.76
4/80 | 1.73E-07 | 3.86 | 1.13E-05 | 3.85 | 1.63E-07 | 3.94 | 4.16E-06 | 3.87
4/160 | 1.21E-08 | 3.84 | 8.28E-07 | 3.77 | 1.05E-08 | 3.96 | 2.83E-07 | 3.88

Table 2: 2D Burgers equation with the initial condition
$u(x,y,0)=0.5+\sin(\pi(x+y)/2)$, with periodic boundary conditions in both
directions, $t=0.5/\pi$, Triangular mesh with size $h$, $L_{1}$ and
$L_{\infty}$ errors for $\mathbf{P}^{1}$, $\mathbf{P}^{2}$ and
$\mathbf{P}^{3}$ based DGM Figure 5: Burgers equation - Sample mesh where the
mesh points on the boundary are uniformly distributed with cell length $h$ =
4/20

Example 2: We solve the two dimensional Euler equations:

$\textbf{U}_{t}+\textbf{f(U)}_{x}+\textbf{g(U)}_{y}=0\qquad\qquad(x,y)\in[0,2]\times[0,2]$
(34)

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(\pi(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 in both directions. The
exact solution is given by $\rho(x,y,t)=1+0.2\sin(\pi(x+y-t))$,
$u(x,y,t)=0.7$, $v(x,y,t)=0.3$ and $p(x,y,t)=1.0$. A sample mesh used is shown
in Figure 6. The errors in density and numerical orders of accuracy are
calculated at $t=2$ by marking all the cells as troubled cells 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.

|  | DG without limiter | DG with limiter
---|---|---|---
$h$ |  $L_{1}$ error | Order |  $L_{\infty}$ error | Order |  $L_{1}$ error | Order |  $L_{\infty}$ error | Order
$\mathbf{P}^{1}$ | 2/20 | 6.69E-03 |  | 2.87E-02 |  | 2.72E-03 |  | 6.95E-03 |
2/40 | 1.51E-03 | 2.15 | 7.58E-03 | 1.92 | 6.90E-04 | 1.98 | 1.81E-03 | 1.94
2/80 | 3.13E-04 | 2.27 | 1.54E-03 | 2.30 | 1.74E-04 | 1.99 | 4.53E-04 | 2.00
2/160 | 8.05E-05 | 1.96 | 3.57E-04 | 2.11 | 4.44E-05 | 1.97 | 1.17E-04 | 1.95
$\mathbf{P}^{2}$ | 2/20 | 1.02E-04 |  | 2.27E-03 |  | 8.21E-05 |  | 3.31E-03 |
2/40 | 1.15E-05 | 3.15 | 3.11E-04 | 2.87 | 1.03E-05 | 2.99 | 5.97E-04 | 2.47
2/80 | 1.23E-06 | 3.22 | 4.66E-05 | 2.74 | 1.28E-06 | 3.01 | 9.19E-05 | 2.70
2/160 | 1.55E-07 | 2.99 | 6.29E-06 | 2.89 | 1.69E-07 | 2.92 | 1.24E-05 | 2.89
$\mathbf{P}^{3}$ | 2/20 | 2.31E-06 |  | 4.25E-05 |  | 8.32E-07 |  | 7.24E-06 |
2/40 | 1.50E-07 | 3.95 | 3.25E-06 | 3.71 | 4.92E-08 | 4.08 | 5.31E-07 | 3.77
2/80 | 9.91E-09 | 3.92 | 2.27E-07 | 3.84 | 2.99E-09 | 4.04 | 4.32E-08 | 3.62
2/160 | 7.07E-10 | 3.81 | 1.52E-08 | 3.90 | 1.88E-10 | 3.99 | 2.93E-09 | 3.88

Table 3: 2D Euler equations with the initial condition
$\rho(x,y,0)=1+0.2\sin(\pi(x+y))$, $u(x,y,0)=0.7$, $v(x,y,0)=0.3$ and
$p(x,y,0)=1.0$, with periodic boundary conditions in both directions, $t=2$,
Triangular mesh with size $h$, $L_{1}$ and $L_{\infty}$ errors for density
with $\mathbf{P}^{1}$, $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM Figure
6: 2D Euler equations with the initial condition
$\rho(x,y,0)=1+0.2\sin(\pi(x+y))$, $u(x,y,0)=0.7$, $v(x,y,0)=0.3$ and
$p(x,y,0)=1.0$ \- Sample mesh where the mesh points on the boundary are
uniformly distributed with cell length $h$ = 2/20

Example 3: We again solve the two dimensional Euler equations given by (34)
for the Isentropic Euler Vortex problem suggested by Shu [24] as a test case
in the domain $[0,10]\times[-5,5]$. The analytical 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 analytical solution at $t=0$ and use
periodic boundary conditions at the edges of the domain in both directions. A
sample mesh used is shown in Figure 7. The errors in density and numerical
orders of accuracy are calculated at $t=2$ by marking all the cells as
troubled cells 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.

|  | DG without limiter | DG with limiter
---|---|---|---
$h$ |  $L_{1}$ error | Order |  $L_{\infty}$ error | Order |  $L_{1}$ error | Order |  $L_{\infty}$ error | Order
$\mathbf{P}^{1}$ | 10/20 | 4.57E-02 |  | 1.45E-01 |  | 4.35E-02 |  | 1.35E-01 |
10/40 | 1.00E-02 | 2.19 | 3.58E-02 | 2.02 | 1.10E-02 | 1.98 | 3.75E-02 | 1.85
10/80 | 2.27E-03 | 2.14 | 1.11E-02 | 1.69 | 2.81E-03 | 1.97 | 1.01E-02 | 1.90
10/160 | 6.00E-04 | 1.92 | 3.23E-03 | 1.78 | 7.07E-04 | 1.99 | 2.67E-03 | 1.92
$\mathbf{P}^{2}$ | 10/20 | 6.67E-03 |  | 9.93E-02 |  | 4.89E-03 |  | 8.39E-02 |
10/40 | 7.67E-04 | 3.12 | 1.38E-02 | 2.85 | 6.20E-04 | 2.98 | 1.13E-02 | 2.89
10/80 | 9.26E-05 | 3.05 | 1.88E-03 | 2.88 | 7.91E-05 | 2.97 | 1.80E-03 | 2.65
10/160 | 1.19E-05 | 2.96 | 3.12E-04 | 2.59 | 1.02E-05 | 2.95 | 2.73E-04 | 2.72
$\mathbf{P}^{3}$ | 10/20 | 8.11E-04 |  | 8.49E-03 |  | 8.27E-04 |  | 7.82E-03 |
10/40 | 5.03E-05 | 4.01 | 6.49E-04 | 3.71 | 5.28E-05 | 3.97 | 5.85E-04 | 3.74
10/80 | 3.35E-06 | 3.91 | 5.28E-05 | 3.62 | 3.44E-06 | 3.94 | 4.53E-05 | 3.69
10/160 | 2.32E-07 | 3.85 | 4.09E-06 | 3.69 | 2.32E-07 | 3.89 | 3.34E-06 | 3.76

Table 4: 2D Euler equations for the Isentropic Vortex problem with periodic
boundary conditions in both directions, $t=2$, Triangular mesh with size $h$,
$L_{1}$ and $L_{\infty}$ errors for density with $\mathbf{P}^{1}$,
$\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM Figure 7: 2D Euler equations
for Isentropic Vortex Problem - Sample mesh where the mesh points on the
boundary are uniformly distributed with cell length $h$ = 10/20

###  4.2 Test Cases With Shocks

We now test the CSWENO (compact subcell WENO) limiter for problems with
solutions having shocks.

Example 4: We solve the problem of shock interaction with entropy waves as
proposed in [25] in the two-dimensional domain. We solve the 2D Euler
equations with a moving shock interacting with sine waves in density in the
domain $[0,1]\times[0,1]$ with the initial conditions given as
$(\rho,u,v,p)=(3.857143,2.629369,0.0,10.333333)$ for $x<0.125$ and
$(\rho,u,v,p)=(1.0+0.2\sin(16\pi x),0,0,1)$ otherwise. Non reflecting boundary
condition is applied at $x=0$ and $x=1$ and periodic boundary conditions are
applied at the other two boundaries. The computed solution for density
obtained at $t=0.178$s using $h=1/200$ at the $y=0.5$ line while using the
CSWENO limiter for $\mathbf{P}^{1}$, $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$
based DGM is compared and plotted against the exact solution in Figure 8a. A
surface color plot for density is shown in Figure 8b. We can see that good
resolution is obtained in the solution for this problem using CSWENO limiter.

(a) Comparison of density solution on $y=0.5$ line of Shock entropy wave
Problem for $\mathbf{P}^{1}$, $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM
with the exact solution

(b) Surface plot for density for $\mathbf{P}^{3}$ based DGM

Figure 8: Solution of Shock entropy wave Problem at $t=0.178$ with $h=1/200$

Example 5: We solve the double Mach reflection problem as given in [26]. We
solve the two-dimensional Euler equations 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 a mesh size of $h=1/200$
which contains 370,046 triangles. A sample mesh of size $h=1/20$ is shown in
Figure 9. The density contours for the solution obtained using the CSWENO
limiter for $\mathbf{P}^{1}$, $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM
is shown in Figure 10. A zoom-in view of the density contours near the double
Mach stem is shown in Figure 11 for $\mathbf{P}^{1}$, $\mathbf{P}^{2}$ and
$\mathbf{P}^{3}$ based DGM. We can see that the solution obtained using CSWENO
limiter is quite well comparable to the solution obtained in [26].

Figure 9: Double Mach reflection problem - Sample mesh where the mesh points
on the boundary are uniformly distributed with cell length $h$ = 1/20

(a) Density contours for the solution at t=0.2 with $\mathbf{P}^{1}$ based DGM

(b) Density contours for the solution at t=0.2 with $\mathbf{P}^{2}$ based DGM

(c) Density contours for the solution at t=0.2 with $\mathbf{P}^{3}$ based DGM

(d) Density Range

Figure 10: 50 equally spaced density contours for Double Mach reflection
solution using the CSWENO limiter with $h=1/200$ for $\mathbf{P}^{1}$(top),
$\mathbf{P}^{2}$(middle) and $\mathbf{P}^{3}$(bottom) based DGM

(a) $\mathbf{P}^{1}$ based DGM

(b) $\mathbf{P}^{2}$ based DGM

(c) $\mathbf{P}^{3}$ based DGM

(d) Density Range

Figure 11: Density variation for Double Mach reflection solution using the
CSWENO limiter for $\mathbf{P}^{1}$, $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$
based DGM in the region $[2,2.9]\times[0,0.6]$ using 50 equally spaced
contours

Example 6: As another test problem, we solve the flow over a forward facing
step which is given in [26]. We solve the two-dimensional Euler equations in a
flow set up which contains a right-going Mach 3 uniform flow in a wind tunnel
of width 1 unit and length 4 units. The step height is 0.2 units and is
located 1 units from the left hand end of the wind tunnel. The problem is
initialized by a uniform, right-going Mach 3 flow. Reflective boundary
conditions are applied along the walls of the tunnel and in-flow and out-flow
boundary conditions are applied at the entrance and the exit, respectively. We
compute the solution upto time $t=4.0$ for a mesh size of $h=1/100$ which
contains 78982 triangles. A sample mesh of size $h=1/20$ is shown in Figure
12. The density contours for the solution obtained using the CSWENO limiter
for $\mathbf{P}^{1}$, $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM is shown
in Figure 13. We can see that the solution obtained using CSWENO limiter is
quite well comparable to the solution obtained in [26].

Figure 12: Flow over a forward facing step - Sample mesh where the mesh points
on the boundary are uniformly distributed with cell length $h$ = 1/20

(a) Density contours for the solution at t=4.0 with $\mathbf{P}^{1}$ based DGM

(b) Density contours for the solution at t=4.0 with $\mathbf{P}^{2}$ based DGM

(c) Density contours for the solution at t=4.0 with $\mathbf{P}^{3}$ based DGM

(d) Density Range

Figure 13: 50 equally spaced density contours for solution of flow over a
forward facing step using the CSWENO limiter with $h=1/100$ for
$\mathbf{P}^{1}$(top), $\mathbf{P}^{2}$(middle) and $\mathbf{P}^{3}$(bottom)
based DGM

Example 7: As another test case, we look at the 2D Riemann problem of gas
dynamics which is one of the most extensively studied problem which also
contains a lot of intricate flow structures. We solve the two-dimensional
Euler equations in the domain $[0,1]\times[0,1]$ for 2D Riemann problem
configurations (3) and (12) as given by the nomenclature in [27]. The initial
conditions for configurations (3) and (12) are given respectively as

$\displaystyle(\rho,u,v,p)(x,y,0)=\begin{cases}(1.5,0,0,1.5)\quad\text{if
$x\geq 0.5$ and $y\geq 0.5$}\\\ (0.5323,1.206,0,0.3)\quad\text{if $x<0.5$ and
$y\geq 0.5$}\\\ (0.138,1.206,1.206,0.029)\quad\text{if $x<0.5$ and $y<0.5$}\\\
(0.5323,0,1.206,0.3)\quad\text{otherwise}\end{cases}$ (35)

$\displaystyle(\rho,u,v,p)(x,y,0)=\begin{cases}(0.5313,0,0,0.4)\quad\text{if
$x\geq 0.5$ and $y\geq 0.5$}\\\ (1,0.7276,0,1)\quad\text{if $x<0.5$ and $y\geq
0.5$}\\\ (0.8,0,0,1)\quad\text{if $x<0.5$ and $y<0.5$}\\\
(1,0,0.7276,1)\quad\text{otherwise}\end{cases}$ (36)

We compute the solution upto time $t=0.3$ for configuration (3) and till
$t=0.25$ for configuration (12). We use a mesh size of $h=1/200$ which
contains 92552 triangles. A sample mesh of size $h=1/20$ is shown in Figure
14. The density contours for the solution obtained using the CSWENO limiter
for $\mathbf{P}^{1}$, $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM are
shown in Figures 15 and 16 respectively for configurations (3) and (12). We
can see that the solution obtained using CSWENO limiter is quite well
comparable to the solution obtained in [27].

Figure 14: 2D Riemann problem of gas dynamics - Sample mesh where the mesh
points on the boundary are uniformly distributed with cell length $h$ = 1/20

(a) Density contours for the solution at t=0.3 with $\mathbf{P}^{1}$ based DGM

(b) Density contours for the solution at t=0.3 with $\mathbf{P}^{2}$ based DGM

(c) Density contours for the solution at t=0.3 with $\mathbf{P}^{3}$ based DGM

(d) Density Range

Figure 15: 50 equally spaced density contours for solution at $t=0.3$ for 2D
Riemann problem configuration 3 using the CSWENO limiter with $h=1/200$ for
$\mathbf{P}^{1}$, $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM

(a) Density contours for the solution at t=0.3 with $\mathbf{P}^{1}$ based DGM

(b) Density contours for the solution at t=0.3 with $\mathbf{P}^{2}$ based DGM

(c) Density contours for the solution at t=0.3 with $\mathbf{P}^{3}$ based DGM

(d) Density Range

Figure 16: 50 equally spaced density contours for solution at $t=0.25$ for 2D
Riemann problem configuration 12 using the CSWENO limiter with $h=1/200$ for
$\mathbf{P}^{1}$, $\mathbf{P}^{2}$ and $\mathbf{P}^{3}$ based DGM

##  5 Conclusions:

We have generalized the compact subcell WENO (CSWENO) limiting strategy for
the solution of hyperbolic conservation laws using discontinuous Galerkin
method proposed in [1] to unstructured triangular meshes. Using this strategy,
we identify the troubled cells and use only the immediate neighbors by
dividing them into subcells appropriately based on the order of accuracy of
the scheme. These new cells are used for the WENO reconstruction. This
formulation is different from the subcell limiting strategy of Dumbser et al
[9] and Giri et al [10] which is much more accurate but quite complicated as
they use subcells in an a posteriori limiting strategy. This limiting strategy
can be used with any of the WENO reconstructions given in [11] (called type I
WENO reconstruction), or [5] (called type II WENO reconstruction) or [12]
(mixed reconstruction) or the more recent methods given in [13] and [2] by
dividing the immediate neighbors appropriately. We use the WENO reconstruction
given in [2] as it is quite simple in implementation and extension to higher
orders is easy. We emphasize that our procedure works with any of the WENO
reconstructions listed above by dividing the neighbors of the troubled cell as
is required by the reconstruction. 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 for standard test cases which are solutions of two-
dimensional Euler equations to illustrate the performance of the limiter.

## References

  * [1] S. R. Siva Prasad Kochi and M. Ramakrishna, “A compact subcell WENO limiting strategy using immediate neighbours for Runge-Kutta discontinuous Galerkin methods.,” International Journal of Computer Mathematics, vol. 98, no. 3, pp. 608–626, 2021.
  * [2] J. Zhu and C.-W. Shu, “A new type of multi-resolution WENO schemes with increasingly higher order of accuracy on triangular meshes.,” Journal of Computational Physics, vol. 392, pp. 19–33, 2019.
  * [3] B. Cockburn and C.-W. Shu, “Runge-Kutta discontinuous Galerkin method for convection-dominated problems.,” J. Sci. Comput., vol. 16, pp. 173–261, 2001.
  * [4] J. Qiu and C.-W. Shu, “Runge-Kutta discontinuous Galerkin method using WENO limiters.,” SIAM J. Sci. Comput., vol. 26, pp. 907–929, 2005.
  * [5] J. Zhu, J. Qiu, C.-W. Shu, and M. Dumbser, “Runge-Kutta discontinuous Galerkin method using WENO limiters II: Unstructured meshes.,” Journal of Computational Physics, vol. 227, pp. 4330–4353, 2008.
  * [6] J. Zhu, X. Zhong, C.-W. Shu, and J. Qiu, “Runge–Kutta discontinuous Galerkin method using a new type of WENO limiters on unstructured meshes.,” Journal of Computational Physics, vol. 248, pp. 200–220, 2013.
  * [7] J. Zhu, X. Zhong, C.-W. Shu, and J. Qiu, “Runge-Kutta Discontinuous Galerkin Method with a Simple and Compact Hermite WENO Limiter on Unstructured Meshes.,” Commun. Comput. Phys., vol. 21, no. 3, pp. 623–649, 2017.
  * [8] H. Zhu, J. Qiu, and J. Zhu, “A simple, high-order and compact WENO limiter for RKDG method.,” Computers and Mathematics with Applications, vol. 79, pp. 317–336, 2020.
  * [9] M. Dumbser and R. Loubere, “A simple robust and accurate a posteriori sub-cell finite volume limiter for the discontinuous Galerkin method on unstructured meshes.,” Journal of Computational Physics, vol. 319, pp. 163–199, 2016.
  * [10] P. Giri and J. Qiu, “A high-order Runge-Kutta discontinuous Galerkin method with a subcell limiter on adaptive unstructured grids for two-dimensional compressible inviscid flows,” International Journal for Numerical Methods in Fluids, vol. 91, no. 8, pp. 367–394, 2019.
  * [11] M. Dumbser and M. Kaser, “Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems.,” Journal of Computational Physics, vol. 221, pp. 693–723, 2007.
  * [12] Y. Liu and Y.-T. Zhang, “A Robust Reconstruction for Unstructured WENO Schemes.,” J. Sci. Comput., vol. 54, pp. 603–621, 2013.
  * [13] J. Zhu and J. Qiu, “New finite volume weighted essentially nonoscillatory schemes on triangular meshes.,” SIAM J. Sci. Comput., vol. 40, no. 2, pp. A903–A928, 2018.
  * [14] C. Hu and C.-W. Shu, “Weighted essentially non-oscillatory schemes on triangular meshes.,” Journal of Computational Physics, vol. 150, pp. 97–127, 1999.
  * [15] J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications.  Springer New York, 2008.
  * [16] L. Krivodonova, J. Xin, J.-F. Remacle, N. Chevaugeon, and J. Flaherty, “Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws.,” Appl. Numer. Math., vol. 48, pp. 323–338, 2004.
  * [17] J. Qiu and C.-W. Shu, “A comparison of troubled-cell indicators for Runge-Kutta discontinuous Galerkin methods using weighted essentially nonoscillatory limiters.,” SIAM J. Sci. Comput., vol. 27, pp. 995–1013, 2005.
  * [18] D. M. Williams, L. Shunn, and A. Jameson, “Symmetric quadrature rules for simplexes based on sphere close packed lattice arrangements.,” Journal of Computational and Applied Mathematics, vol. 266, pp. 18–38, 2014.
  * [19] T. Sonar, “On the construction of essentially non-oscillatory finite volume approximations to hyperbolic conservation laws on general triangulations: polynomial recovery, accuracy and stencil selection.,” Comput. Methods Appl. Mech. Engrg., vol. 140, pp. 157–181, 1997.
  * [20] G. Jiang and C.-W. Shu, “Efficient implementation of weighted ENO schemes.,” Journal of Computational Physics, vol. 126, pp. 202–228, 1996\.
  * [21] M. Castro, B. Costa, and W. S. Don, “High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws.,” Journal of Computational Physics, vol. 230, pp. 1766–1792, 2011.
  * [22] C.-W. Shu, “TVD time discretizations.,” SIAM J. Sci. Stat. Comput., vol. 9, pp. 1073–1084, 1988.
  * [23] C. Geuzaine and J.-F. Remacle, “Gmsh: a three-dimensional finite element mesh generator withbuilt-in pre- and post-processing facilities.,” International Journal for Numerical Methods in Engineering, vol. 79, no. 11, pp. 1309–1331, 2009.
  * [24] C.-W. Shu, “Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws.,” Lecture Notes in Mathematics, Springer, vol. 1697, pp. 325–432, 1998.
  * [25] C.-W. Shu and S. Osher, “Effective implementation of essentially non-oscillatory shock-capturing schemes, II.,” Journal of Computational Physics, vol. 83, pp. 32–78, 1989.
  * [26] P. Woodward and P. Colella, “The numerical simulation of two-dimensional fluid flow with strong shocks.,” Journal of Computational Physics, vol. 54, pp. 115–173, 1984.
  * [27] P. D. Lax and X. D. Liu, “Solution of two-dimensional Riemann problem of gas dynamics by positive schemes.,” SIAM J. Sci. Comput., vol. 19, pp. 319–340, 1998.