The jump filter in the discontinuous Galerkin method for hyperbolic conservation laws

Consider the scalar conservation laws in one-dimensional form

$\displaystyle\begin{cases}u_{t}+f(u)_{x}=0,\;x\in\Omega=[a,b],\\ u(x,0)=u_{0}(x),\end{cases}$

(2.1)

with periodic or compactly supported boundary conditions. Let $\mathcal{T}_{h}$ denote the partition of $\Omega$ with the mesh denoted by $K_{j}=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$ for $j=1,\cdots,N$. The center of the cell is $x_{j}=\frac{1}{2}(x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}})$ and the mesh size is denoted by $h_{j}=x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}$ with $h=\max\limits_{1\leq j\leq N}h_{j}$ being the maximum cell size. The mesh is assumed to be regular, which means the ratio between the maximum and minimum mesh sizes keeps bounded in the mesh refinements.

The DG approximation space $V_{h}^{k}$ is the space of all piecewise polynomial functions defined on $\Omega$:

$\displaystyle V_{h}^{k}=\{v\in L^{2}([a,b]):\;v|_{K_{j}}\in P^{k}(K_{j}),\;j=1,\cdots,N\}.$

(2.2)

We also denote the inner product over the interval $K_{j}$ and the associated norm by

$\displaystyle(w,r)_{K_{j}}=\int_{K_{j}}\,w\,r\,\mathrm{d}x,\quad\left\|w\right\|_{K_{j}}=\sqrt{(w,w)_{K_{j}}}.$

(2.3)

It results in the semi-discrete DG scheme as follows: Find $u_{h}\in V_{h}^{k}$, such that

$\displaystyle(\partial_{t}u_{h},v)_{K_{j}}+<\widehat{f(u_{h})},v>_{\partial K_{j}}-(f(u_{h}),v_{x})_{K_{j}}=0,\;\forall v\in V_{h}^{k}\mbox{ and }j=1,\cdots,N,$

(2.4)

where $<\widehat{f(u_{h})},v>_{\partial K_{j}}=\hat{f}_{j+\frac{1}{2}}v(x_{j+\frac{1}{2}}^{-})-\hat{f}_{j-\frac{1}{2}}v(x_{j-\frac{1}{2}}^{+})$. Here, $\hat{f}_{j+\frac{1}{2}}$ is the numerical flux defined on the cell interface $\partial K_{j}$ and in general depends on the values of $u_{h}$ from both sides of the interface $\hat{f}_{j+\frac{1}{2}}=\hat{f}(u_{h}(x_{j+\frac{1}{2}}^{-}),u_{h}(x_{j+\frac{1}{2}}^{+}))$. The Godunov flux or the Lax-Friedrichs flux is often employed when solving scalar equations.

In the linear case of Eqn. (2.1), the smoothness of the solution depends only on the initial conditions. However, for general nonlinear cases, the discontinuous solutions, known as shocks, may develop even if the initial data are smooth. The high order DG scheme would generate spurious oscillations when the solution contains shocks, so-called the Gibbs phenomenon. These spurious oscillations may not only generate some artificial structures, but also cause many overshoots and undershoots that make the numerical scheme less robust. To overcome this issue, the limiter approach was proposed, such as the TVB limiters [12, 11, 13], the WENO limiters [65, 85], and the vertex-based limiters [44, 45, 46], which is processed after each time stage. Another strategy, known as the entropy viscosity approach, is a predominant approach in the continuous Galerkin (CG) community for solving hyperbolic conservation laws [34, 23, 60, 22], which is based on the viscosity solution [47]. The viscosity solution of the Eqn. (2.1) is defined as the solution to the equation

$\displaystyle u^{\varepsilon}_{t}+f(u^{\varepsilon})_{x}=\epsilon u^{\varepsilon}_{xx},$

(2.5)

with $u^{\varepsilon}(x,0)=u(x,0)$. The physical solution $u(x,t)$ is obtained as the limit viscosity solution $u(x,t)=\lim\limits_{\varepsilon\rightarrow 0^{+}}u^{\varepsilon}(x,t)$.

We propose the following DG scheme with local viscosity: Find $u_{h}\in V_{h}^{k}$, such that

$\displaystyle(\partial_{t}u_{h},v)_{K_{j}}+<\widehat{f(u_{h})},v>_{\partial K_{j}}-(f(u_{h}),v_{x})_{K_{j}}+(\nu_{j}(u_{h})_{x},v_{x})_{K_{j}}=0,\;\forall v\in V_{h}^{k},$

(2.6)

where $\nu_{j}$ is the local viscosity coefficient and will be determined later. In the DG scheme (2.6), the introduction of this local viscosity term aims to mitigate spurious oscillations in shock solutions without compromising the high-order accuracy of smooth solutions. To achieve this, we define the local viscosity function as follows:

$\displaystyle\nu_{j}=\sigma_{j}(1-\xi^{2}),$

(2.7)

where $\xi=\frac{2(x-x_{j})}{h_{j}}$ and the jump viscosity coefficient $\sigma_{j}$ is given by

$\displaystyle\sigma_{j}=c_{f}(h_{j}\left\|[\mspace{-2.5mu}[u_{h}]\mspace{-2.5mu}]\right\|_{\partial K_{j}}+\sum\limits_{l=1}^{k}l(l+1)h_{j}^{l+1}\lVert[\mspace{-2.5mu}[\partial_{x}^{l}u_{h}]\mspace{-2.5mu}]\rVert_{\partial K_{j}}),$

(2.8)

with $c_{f}\geq 0$ as a carefully designed free parameter. Here, $\left\|[\mspace{-2.5mu}[w]\mspace{-2.5mu}]\right\|_{\partial K_{j}}=|w(x_{j-\frac{1}{2}}^{+})-w(x_{j-\frac{1}{2}}^{-})|+|w(x_{j+\frac{1}{2}}^{+})-w(x_{j+\frac{1}{2}}^{-})|$ denotes the jump of $w$ at cell interface $\partial K_{j}$, and $\partial^{l}_{x}u_{h}$ represents the $l$-th derivative of $u_{h}$ with respect to $x$.

Assuming $\sigma_{j}$ is constant in $K_{j}$, which can be estimated by the numerical solution at the current time stage and then through integration by parts, the local viscosity term is equivalent to:

	$\displaystyle(\nu_{j}(u_{h})_{x},v_{x})_{K_{j}}$	$\displaystyle=-\sigma_{j}(\partial_{x}((1-\xi^{2})\partial_{x}u_{h}),v)_{K_{j}},$		(2.9)
		$\displaystyle=-\sigma_{j}\frac{2}{h_{j}}(\partial_{\xi}((1-\xi^{2})\partial_{\xi}\tilde{u}_{h}),\tilde{v})_{\widehat{K}},$		(2.10)

where the reference cell $\widehat{K}=[-1,1]$. The numerical solution on $K_{j}$ can be expressed by the Gauss-Legendre basis

$\displaystyle u_{h}|_{K_{j}}=\sum\limits_{l=0}^{k}u_{j}^{l}(t)P_{l}(x)=\sum\limits_{l=0}^{k}u_{j}^{l}(t)\tilde{P}_{l}(\xi).$

(2.11)

Since $\tilde{P}_{l}(\xi)$ satisfies the Sturm-Liouville equation

$\displaystyle\frac{d}{d\xi}((1-\xi^{2})\frac{d}{d\xi}\tilde{P}_{l}(\xi))+l(l+1)\tilde{P}_{l}(\xi)=0,$

(2.12)

we obtain

$\displaystyle(\nu_{j}(u_{h})_{x},v_{x})_{K_{j}}$

$\displaystyle=\sigma_{j}\frac{2}{h_{j}}\sum\limits_{l=0}^{k}u_{j}^{l}(l(l+1)\tilde{P}_{l}(\xi),\tilde{v})_{\widehat{K}}.$

(2.13)

If we were to directly implement this term in (2.6), it would affect the maximum allowable time step when using an explicit time discretization. Alternatively, this term can be included in a time-splitting fashion [20, 31]. In this approach, we first advance the original DG scheme (2.4) by one time stage, followed by incorporating the diffusion term

$\displaystyle(\partial_{t}u_{h},v)_{K_{j}}+(\nu_{j}(u_{h})_{x},v_{x})_{K_{j}}=0,$

(2.14)

which is equivalent to

$\displaystyle(\partial_{t}\tilde{u}_{h},\tilde{v})_{\widehat{K}}+\sigma_{j}(\frac{2}{h_{j}})^{2}\sum\limits_{l=0}^{k}u_{j}^{l}(l(l+1)\tilde{P}_{l}(\xi),\tilde{v})_{\widehat{K}}=0.$

(2.15)

By taking $\tilde{v}=\tilde{P}_{l}(\xi)$, it results that

$\displaystyle\frac{d}{dt}u_{j}^{l}+\sigma_{j}^{l}u_{j}^{l}=0,\;l=0,\cdots,k,$

(2.16)

where

$\displaystyle\sigma_{j}^{l}=\sigma_{j}(\frac{2}{h_{j}})^{2}l(l+1).$

(2.17)

For one time step with size $\Delta t$, this ODE system can be solved as follows:

$\displaystyle u_{j}^{l}(t+\Delta t)=\exp(-\sigma_{j}^{l}\Delta t)u_{j}^{l}(t),\;l=0,\cdots,k.$

(2.18)

Notice that for $l=0$, $u_{j}^{0}$ is the cell average and $\sigma_{j}^{0}=0$, thus the local cell conservation of the DG scheme is preserved.

In (2.17), to further minimize numerical dissipation, we opt for a smaller mode damping parameter $\sigma_{j}^{l}$ defined as follows:

	$\displaystyle\sigma_{j}^{l}=\sigma_{j,l}(\frac{2}{h_{j}})^{2}l(l+1),$		(2.19)
	$\displaystyle\sigma_{j,l}=c_{f}^{l}(h_{j}\left\|[\mspace{-2.5mu}[u_{h}]\mspace{-2.5mu}]\right\|_{\partial K_{j}}+\sum\limits_{m=1}^{l}m(m+1)h_{j}^{m+1}\lVert[\mspace{-2.5mu}[\partial_{x}^{m}u_{h}]\mspace{-2.5mu}]\rVert_{\partial K_{j}}).$		(2.20)

Here, $\sigma_{j,l}$ replaces $\sigma_{j}$ in equation (2.17), focusing solely on jumps within the mode number $l$. Compared to equation (2.17), this damping technique eliminates the need to consider jump information from higher-order derivatives when filtering low-level moments. Consequently, it enhances the preservation of original polynomial data, thereby maintaining the accuracy of the standard DG scheme more effectively.

Similar ideas are explored in [52, 49], wherein an additional damping source term is introduced into the DG scheme:

$\displaystyle(\partial_{t}u_{h},v)_{K_{j}}+<\widehat{f(u_{h})},v>_{\partial K_{j}}-(f(u_{h}),v_{x})_{K_{j}}-\sum\limits_{l=0}^{k}\frac{\tilde{\sigma}_{j}^{l}}{h_{j}}(u_{h}-P_{h}^{l-1}u_{h},v)_{K_{j}}=0,$

(2.21)

where

$\displaystyle\tilde{\sigma}_{j}^{\ell}=\frac{2(2l+1)}{2k-1}\frac{(h_{j})^{l}}{l!}\left([\mspace{-2.5mu}[\partial_{x}^{\ell}u_{h}]\mspace{-2.5mu}]_{j+\frac{1}{2}}^{2}+[\mspace{-2.5mu}[\partial_{x}^{\ell}u_{h}]\mspace{-2.5mu}]_{j-\frac{1}{2}}^{2}\right)^{\frac{1}{2}},$

(2.22)

and $P_{h}^{l}$ denotes the standard $L^{2}$ projection into $P^{l}$ polynomial space and $P_{h}^{-1}=P_{h}^{0}$. Here, $[\mspace{-2.5mu}[w]\mspace{-2.5mu}]_{j+\frac{1}{2}}=w(x_{j+\frac{1}{2}}^{+})-w(x_{j+\frac{1}{2}}^{-})$.

The conservation property of the scheme can be readily obtained by setting $v=1$ in (2.6), leading to

$\displaystyle\dfrac{\mathrm{d}}{\mathrm{d}t}\int_{K_{j}}u_{h}\mathrm{d}x=-\hat{f}_{j+\frac{1}{2}}+\hat{f}_{j-\frac{1}{2}}.$

(2.23)

Summing over $j$, with either periodic or compactly supported boundary conditions, we achieve global conservation. The $L^{2}$ stability of the semi-discrete scheme is achieved by taking $v=u_{h}$ in (2.6). Summing it over $j$, we obtain

$\displaystyle\dfrac{1}{2}\dfrac{\mathrm{d}}{\mathrm{d}t}\left\|u_{h}\right\|^{2}=-\sum_{j}\Theta_{j+\frac{1}{2}}-\sum_{j}\left\|\sqrt{\nu_{j}}u_{h}\right\|^{2}_{K_{j}},$

(2.24)

where $\Theta_{j+\frac{1}{2}}=\int_{u_{h}(x_{j+\frac{1}{2}}^{-})}^{u_{h}(x_{j+\frac{1}{2}}^{+})}\left(f(y)-\hat{f}(u_{h}(x_{j+\frac{1}{2}}^{-}),u_{h}(x_{j+\frac{1}{2}}^{+}))\right)\,\mathrm{d}y\geq 0$ as for the cell entropy inequality for (2.4) in [38], and $\left\|u_{h}(\cdot,t)\right\|$ denotes the global $L^{2}$ norm. Consequently, we have the $L^{2}$ stability $\left\|u_{h}(\cdot,t)\right\|\leq\left\|u_{h}(\cdot,0)\right\|$.

In the two-dimensional case, the situation is similar. Consider the two-dimensional scalar conservation laws

$\displaystyle\begin{cases}u_{t}+\nabla\cdot\boldsymbol{F}(u)=0,\;\\ u(\boldsymbol{x},0)=u_{0}(\boldsymbol{x}),\end{cases}$

(2.29)

with suitable boundary conditions on the rectangular domain $\Omega$. After we have the regular Cartesian grid $\mathcal{T}_{h}$ of $\Omega$, the cell $K_{i,j}\in\mathcal{T}_{h}$ is denoted by $K_{i,j}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times[y_{j-\frac{1}{2}},y_{j+\frac{1}{2}}]$ for $i=1,\cdots,N_{x}$; $j=1,\cdots,N_{y}$. The center of the cell $K_{i,j}$ is $(x_{i},y_{j})=(\frac{1}{2}(x_{i-\frac{1}{2}}+x_{i+\frac{1}{2}}),\frac{1}{2}(y_{j-\frac{1}{2}}+y_{j+\frac{1}{2}}))$ and the mesh size is denoted by $h_{x_{i}}=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$, $h_{y_{j}}=y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}}$.

The DG scheme with local viscosity on the rectangular mesh is: Find $u_{h}\in V_{h}^{k}=\{v\in L^{2}(\Omega):\;v|_{K_{i,j}}\in P^{k}(K_{i,j}),\;i=1,\cdots,N_{x};\;j=1,\cdots,N_{y}\}$, such that for any $v\in V_{h}^{k}$, we have

$$(\partial_{t}u_{h},v)_{{K_{i,j}}}-(\boldsymbol{F}(u_{h}),\nabla v)_{{K_{i,j}}}+<\hat{\boldsymbol{F}}_{h}^{\boldsymbol{n}},v>_{\partial K_{i,j}}+(\boldsymbol{\nu}_{i,j}\odot\nabla u_{h},\nabla v)_{K_{i,j}}=0,\;$$

(2.30)

where $<\hat{\boldsymbol{F}}_{h}^{\boldsymbol{n}},v>_{\partial K_{i,j}}=\int_{\partial K_{i,j}}\hat{\boldsymbol{F}}_{h}\cdot\boldsymbol{n}v\mathrm{d}s$ and $\odot$ denotes Hadamard product. Here, $\boldsymbol{n}=(n_{x},n_{y})^{T}$ is the outward unit normal of the cell boundary $\partial_{K_{i,j}}$ and $\hat{\boldsymbol{F}}_{h}$ is a two point numerical flux. Following the one-dimensional case, we can choose the local viscosity function as

$\displaystyle\boldsymbol{\nu}_{i,j}=(\boldsymbol{\nu}_{i,j}^{(1)},\boldsymbol{\nu}_{i,j}^{(2)})=\sigma_{i,j}((1-\xi^{2}),(1-\eta^{2})),$

(2.31)

where $\xi=\frac{2(x-x_{i})}{h_{x_{i}}}$, $\eta=\frac{2(y-y_{j})}{h_{y_{j}}}$.

By integration by parts, the local viscosity term is equivalent to

$$(\boldsymbol{\nu}_{i,j}\odot\nabla u_{h},\nabla v)_{K_{i,j}}=-\sigma_{i,j}(\partial_{\xi}((1-\xi^{2})\partial_{\xi}\tilde{u}_{h})+\partial_{\eta}((1-\eta^{2})\partial_{\eta}\tilde{u}_{h}),\tilde{v})_{\widehat{K}},$$

(2.32)

where the reference cell $\widehat{K}=[-1,1]^{2}$.

The numerical solution on $K_{i,j}$ can be expressed by the Gauss-Legendre basis

$\displaystyle u_{h}|_{K_{i,j}}=\sum\limits_{l=0}^{k}\sum\limits_{p+q=l}u_{i,j}^{p,q}(t)P_{p,q}(x,y)=\sum\limits_{l=0}^{k}\sum\limits_{p+q=l}u_{i,j}^{p,q}(t)\tilde{P}_{p,q}(\xi,\eta).$

(2.33)

Since $\tilde{P}_{p,q}(\xi,\eta)$ satisfies the Sturm-Liouville equation

$\displaystyle\frac{\partial}{\partial\xi}((1-\xi^{2})\frac{\partial}{\partial\xi}\tilde{P}_{p,q}(\xi,\eta))+\frac{\partial}{\partial\eta}((1-\eta^{2})\frac{\partial}{\partial\eta}\tilde{P}_{p,q}(\xi,\eta))+\lambda_{p,q}\tilde{P}_{p,q}(\xi,\eta)=0,$

(2.34)

where

$\displaystyle\lambda_{p,q}=p(p+1)+q(q+1).$

(2.35)

Then we have

$\displaystyle(\boldsymbol{\nu}_{i,j}\odot\nabla u_{h},\nabla v)_{K_{i,j}}=\sigma_{i,j}\sum\limits_{l=0}^{k}\sum\limits_{p+q=l}u_{i,j}^{p,q}(t)(\lambda_{p,q}\tilde{P}_{p,q},\tilde{v})_{\widehat{K}}.$

(2.36)

We also use the time splitting fashion to deal with this term as in the one-dimensional case. By taking $\tilde{v}=\tilde{P}_{p,q}(\xi,\eta)$, we obtain the following ODE system

$\displaystyle\frac{d}{dt}u_{i,j}^{p,q}+\sigma_{i,j}^{l}u_{i,j}^{p,q}=0,\;p+q=l=0,\cdots,k,$

(2.37)

where

$\displaystyle\sigma_{i,j}^{l}=\sigma_{i,j}\frac{2}{h_{x_{i}}}\frac{2}{h_{y_{j}}}\lambda_{p,q}.$

(2.38)

For one time step $\Delta t$, it can be solved as follows:

$\displaystyle u_{i,j}^{p,q}(t+\Delta t)=\exp(-\sigma_{i,j}^{l}\Delta t)u_{i,j}^{p,q}(t),\;p+q=l=0,\cdots,k.$

(2.39)

The jump viscosity coefficient $\sigma_{i,j}$ is given by

	$\displaystyle\sigma_{i,j}=c_{f,x}^{p,q}(h_{y_{j}}\left\|[\mspace{-2.5mu}[u_{h}]\mspace{-2.5mu}]\right\|_{\partial K_{i,j}^{\text{vertical}}}+\sum\limits_{l=1}^{k}l(l+1)h_{y_{j}}^{l+1}\sum_{|\boldsymbol{\alpha}|=l}\left\|[\mspace{-2.5mu}[\partial^{\boldsymbol{\alpha}}u_{h}]\mspace{-2.5mu}]\right\|_{\partial K_{i,j}^{\text{vertical}}})$
	$\displaystyle+c_{f,y}^{p,q}(h_{x_{i}}\left\|[\mspace{-2.5mu}[u_{h}]\mspace{-2.5mu}]\right\|_{\partial K_{i,j}^{\text{horizontal}}}+\sum\limits_{l=1}^{k}l(l+1)h_{x_{i}}^{l+1}\sum_{|\boldsymbol{\alpha}|=l}\left\|[\mspace{-2.5mu}[\partial^{\boldsymbol{\alpha}}u_{h}]\mspace{-2.5mu}]\right\|_{\partial K_{i,j}^{\text{horizontal}}}),$		(2.40)

where $c_{f,x}^{p,q}$ and $c_{f,y}^{p,q}$ are a free parameter. Here, the vector $\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2})$ is the multi-index with $|\boldsymbol{\alpha}|=\alpha_{1}+\alpha_{2}$, and $\partial^{\boldsymbol{\alpha}}w$ is defined as $\partial^{\boldsymbol{\alpha}}{w}=\partial^{\alpha_{1}}_{x}\partial^{\alpha_{2}}_{y}{w}$. Besides, $\left\|[\mspace{-2.5mu}[w]\mspace{-2.5mu}]\right\|_{\partial K_{i,j}^{\text{vertical}}}=\frac{1}{h_{y_{j}}}\int_{y_{j-\frac{1}{2}}}^{y_{j+\frac{1}{2}}}\lVert[\mspace{-2.5mu}[w]\mspace{-2.5mu}]_{i+\frac{1}{2},j}\rVert+\lVert[\mspace{-2.5mu}[w]\mspace{-2.5mu}]_{i-\frac{1}{2},j}\rVert\mathrm{d}s$, where $\lVert[\mspace{-2.5mu}[w]\mspace{-2.5mu}]_{i+\frac{1}{2},j}\rVert=|w(x_{i+\frac{1}{2},j}^{-},y)-w(x_{i+\frac{1}{2},j}^{+},y)|$ represents the absolute value of the jump of $w$ across the right interface of an element $K_{i,j}$. A similar formula can be used to express $\left\|[\mspace{-2.5mu}[\partial^{\boldsymbol{\alpha}}u_{h}]\mspace{-2.5mu}]\right\|_{\partial K_{i,j}^{\text{horizontal}}}$. The trapezoidal rule is employed to approximate this damping term along each edge of $K_{i,j}$ in all our two-dimensional tests. Just as in the one-dimensional case, to further minimize numerical dissipation, we can also adjust the mode damping parameter $\sigma_{i,j}^{l}$ as follows:

	$\displaystyle\sigma_{i,j}^{l}=\sigma_{i,j,l}(\frac{2}{h_{x_{i}}})(\frac{2}{h_{y_{j}}})\lambda_{p,q},$		(2.41)
	$\displaystyle\sigma_{i,j,l}=c_{f,x}^{p,q}(h_{y_{j}}\left\|[\mspace{-2.5mu}[u_{h}]\mspace{-2.5mu}]\right\|_{\partial K_{i,j}^{\text{vertical}}}+\sum\limits_{m=1}^{l}m(m+1)h_{y_{j}}^{m+1}\sum_{|\boldsymbol{\alpha}|=m}\left\|[\mspace{-2.5mu}[\partial^{\boldsymbol{\alpha}}u_{h}]\mspace{-2.5mu}]\right\|_{\partial K_{i,j}^{\text{vertical}}})$
	$\displaystyle+c_{f,y}^{p,q}(h_{x_{i}}\left\|[\mspace{-2.5mu}[u_{h}]\mspace{-2.5mu}]\right\|_{\partial K_{i,j}^{\text{horizontal}}}+\sum\limits_{m=1}^{l}m(m+1)h_{x_{i}}^{m+1}\sum_{|\boldsymbol{\alpha}|=m}\left\|[\mspace{-2.5mu}[\partial^{\boldsymbol{\alpha}}u_{h}]\mspace{-2.5mu}]\right\|_{\partial K_{i,j}^{\text{horizontal}}}).$		(2.42)

We replace $\sigma_{i,j}$ by $\sigma_{i,j,l}$ in (2.38), which considers jumps only within the mode number $l$. Additionally, this damping term (2.41) serves to determine the magnitude of the viscosity term throughout this paper, allowing for adaptive damping of various moments based on the jump size. Again, it leads to reduced numerical dissipation compared with (2.40).

Finally, it’s worth noting that there are no fundamental obstacles to extending the jump filter to Cartesian meshes in higher dimensions.

Since the system of conservation laws can be regarded as multiple scalar equations coupled together, extending the jump filter to systems does not pose any fundamental difficulty. Considering one-dimensional system $\boldsymbol{u}_{t}+\boldsymbol{f}(\boldsymbol{u})_{x}=\boldsymbol{0}$ with suitable boundary conditions, following the analogous approach as in the one-dimensional scalar case, we introduce appropriate viscous terms on each component. By employing the time-splitting method, we derive the following system of ODEs:

$\displaystyle\frac{d}{dt}\boldsymbol{u}_{j}^{l}+\sigma_{j}^{l}\boldsymbol{u}_{j}^{l}=\boldsymbol{0},\;l=0,\cdots,k,$

(2.43)

where

$\displaystyle\sigma_{j}^{l}=\max\limits_{1\leq s\leq\mathcal{N}}\sigma_{j}^{s}(\frac{2}{h_{j}})^{2}l(l+1),$

(2.44)

where $\mathcal{N}$ represents the number of equations in the system, and $\sigma_{j}^{s}$ denotes the jump viscosity coefficient based on the $s$-th variable. For the one time step size $\Delta t$, this ODE system can be solved as follows:

$\displaystyle\boldsymbol{u}_{j}^{l}(t+\Delta t)=\exp(-\sigma_{j}^{l}\Delta t)\boldsymbol{u}_{j}^{l}(t),\;l=0,\cdots,k.$

(2.45)

When dealing with the system of conservation laws, it suffices to compute the jump viscosity coefficient for each component of the system and select the maximum one which is aimed to ensure stability of the scheme when solving problems involving strong shock waves.

We are aware that classical limiters like the TVB limiter or the WENO-type limiter often yield satisfactory outcomes when applied in the characteristic space. In our case, we compute the jump and then apply damping to various moments based on the jump magnitude. However, it’s worth noting that the eigenvectors of a matrix are not unique, and selecting different eigenvectors can result in different damping values if we calculate jump in characteristic space, potentially resulting in varying computed results. To address this challenge, we directly compute the jumps in the conservative space and select an appropriate $c_{f}^{l}$ in (2.20), or $c_{f,x}^{p,q}$ and $c_{f,y}^{p,q}$ in (2.42). Through numerical experiments, we have demonstrated that this approach effectively mitigates spurious oscillations even in the presence of strong shocks.

We perform the limiting procedure with (2.19) with the parameter $c_{f}^{l}=\frac{\beta_{j}}{4l(l+1)}$ in (2.20) where $\beta_{j}$ is the spectral radius of the Jacobian matrix $\frac{\partial\boldsymbol{f}}{\partial\boldsymbol{u}}(\bar{\boldsymbol{u}}_{j})$ which contains the information of wave speed and $\bar{\boldsymbol{u}}_{j}$ denotes the averages of the conservative variables of cell $K_{j}$.

Next, we delve into the Euler equations:

$$\boldsymbol{U}_{t}+\boldsymbol{F}(\boldsymbol{U})_{x}=0,$$

(3.47)

with $\boldsymbol{U}=(\rho,m,E)^{T}$, $\boldsymbol{F}(\boldsymbol{U})=(m,\frac{m^{2}}{\rho}+p,\frac{m}{\rho}(E+p))^{T}$, where $m=\rho u$, $E=\frac{1}{2}\rho u^{2}+\rho e$, and $p=(\gamma-1)\rho e$. Here, $\rho$ denotes density, and $u$ and $m$ represent velocity and momentum, respectively. $E$ stands for the total energy, $p$ denotes pressure, and $e$ signifies internal energy. $\gamma$ represents the ratio of specific heats, and $H=(E+p)/\rho$ is the enthalpy. For the Euler equations in 1D, we set $c_{f}^{l}=\frac{\beta_{j}}{4l(l+1)}\frac{1}{H}$ in (2.20). Notably, in Euler equations, we observed that the factor $\frac{1}{H}$ enhances the scheme’s capability to capture strong shock waves. In this section, we present plots of the full polynomials of DG solutions.