Multi-element SIAC filter for shock capturing applied to high-order discontinuous Galerkin spectral element methods

TL;DR

This paper introduces a multi-element SIAC filter for high-order discontinuous Galerkin spectral element methods, enhancing shock capturing while maintaining accuracy in simulations of conservation laws like Euler and MHD equations.

Contribution

It develops a novel multi-element SIAC filtering technique that adaptively captures shocks in high-order DG methods without sacrificing accuracy.

Findings

Effective shock capturing in 2D Euler and MHD tests

Maintains high-order accuracy near discontinuities

Applicable to general systems of conservation laws

Abstract

We build a multi-element variant of the smoothness increasing accuracy conserving (SIAC) shock capturing technique proposed for single element spectral methods by Wissink et al. (B.W. Wissink, G.B. Jacobs, J.K. Ryan, W.S. Don, and E.T.A. van der Weide. Shock regularization with smoothness-increasing accuracy-conserving Dirac-delta polynomial kernels. Journal of Scientific Computing, 77:579--596, 2018). In particular, the baseline scheme of our method is the nodal discontinuous Galerkin spectral element method (DGSEM) for approximating the solution of systems of conservation laws. It is well known that high-order methods generate spurious oscillations near discontinuities which can develop in the solution for nonlinear problems, even when the initial data is smooth. We propose a novel multi-element SIAC filtering technique applied to the DGSEM as a shock capturing method. We design the SIAC filtering such that the numerical scheme remains high-order accurate and that the shock capturing is applied adaptively throughout the domain. The shock capturing method is derived for general systems of conservation laws. We apply the novel SIAC filter to the two-dimensional Euler and ideal magnetohydrodynamics (MHD) equations to several standard test problems with a variety of boundary conditions.

Figures

Tables

Table 1. Table 1. Euler convergence test for N = 7 , CFL = 0.1 formulae-sequence 𝑁 7 CFL 0.1 N=7,\texttt{CFL}=0.1 and T = 0.4 𝑇 0.4 T=0.4 without filtering.

$N_Q$	${ϵ}_{\infty }$	$EO{C}_{\infty }$	${ϵ}_{\mathrm{cons}}$
$1$	$5.72\cdot {10}^{-3}$	—	$4.4\cdot {10}^{-16}$
$2^2$	$4.56\cdot {10}^{-5}$	$6.97$	$1.3\cdot {10}^{-15}$
$4^2$	$1.74\cdot {10}^{-7}$	$8.03$	$3.1\cdot {10}^{-15}$
$8^2$	$4.82\cdot {10}^{-10}$	$8.50$	$8.9\cdot {10}^{-16}$
${16}^2$	$1.79\cdot {10}^{-12}$	$8.07$	$2.0\cdot {10}^{-14}$

Table 2. Table 2. Euler convergence test for N = 7 , CFL = 0.1 formulae-sequence 𝑁 7 CFL 0.1 N=7,\texttt{CFL}=0.1 and T = 0.4 𝑇 0.4 T=0.4 filtered with ( m , k ) = ( 1 , 6 ) , N d = 0.8 formulae-sequence 𝑚 𝑘 1 6 subscript 𝑁 𝑑 0.8 (m,k)=(1,6),N_{d}=0.8 .

$N_Q$	${ϵ}_{\infty }$	$EO{C}_{\infty }$	${ϵ}_{\mathrm{cons}}$
${10}^2$	$7.39\cdot {10}^{-2}$	—	$5.7\cdot {10}^{-6}$
${20}^2$	$3.97\cdot {10}^{-2}$	$0.90$	$4.6\cdot {10}^{-6}$
${40}^2$	$2.06\cdot {10}^{-2}$	$0.95$	$1.8\cdot {10}^{-6}$
${80}^2$	$1.05\cdot {10}^{-2}$	$0.97$	$5.4\cdot {10}^{-7}$

Table 3. Table 3. Euler convergence test for N = 7 , CFL = 0.1 formulae-sequence 𝑁 7 CFL 0.1 N=7,\texttt{CFL}=0.1 and T = 0.4 𝑇 0.4 T=0.4 filtered with ( m , k ) = ( 3 , 6 ) , N d = 2.5 formulae-sequence 𝑚 𝑘 3 6 subscript 𝑁 𝑑 2.5 (m,k)=(3,6),N_{d}=2.5 .

$N_Q$	${ϵ}_{\infty }$	$EO{C}_{\infty }$	${ϵ}_{\mathrm{cons}}$
$2^2$	$2.52\cdot {10}^{-2}$	—	$1.1\cdot {10}^{-3}$
$4^2$	$3.56\cdot {10}^{-3}$	$2.82$	$3.8\cdot {10}^{-5}$
$8^2$	$4.46\cdot {10}^{-4}$	$3.00$	$1.3\cdot {10}^{-6}$
${16}^2$	$5.56\cdot {10}^{-5}$	$3.00$	$3.6\cdot {10}^{-8}$

Table 4. Table 4. Euler convergence test for N = 7 , CFL = 0.1 formulae-sequence 𝑁 7 CFL 0.1 N=7,\texttt{CFL}=0.1 and T = 0.4 𝑇 0.4 T=0.4 filtered with ( m , k ) = ( 5 , 7 ) , N d = 4.5 formulae-sequence 𝑚 𝑘 5 7 subscript 𝑁 𝑑 4.5 (m,k)=(5,7),N_{d}=4.5 .

$N_Q$	${ϵ}_{\infty }$	$EO{C}_{\infty }$	${ϵ}_{\mathrm{cons}}$
$1$	$1.14\cdot {10}^{-3}$	—	$9.5\cdot {10}^{-5}$
$2^2$	$4.35\cdot {10}^{-5}$	$4.71$	$1.4\cdot {10}^{-6}$
$4^2$	$1.35\cdot {10}^{-6}$	$5.01$	$1.2\cdot {10}^{-8}$
$8^2$	$4.21\cdot {10}^{-8}$	$5.00$	$8.9\cdot {10}^{-11}$

Table 5. Table 5. Initial primitive states for the magnetic rotor test.

Location	$\varrho$	$v_1$	$v_2$	$p$	$B_1$
	$10$	$\frac{u_0}{r_0}\left(\frac{1}{2}-y\right)$	$\frac{u_0}{r_0}\left(x-\frac{1}{2}\right)$	$1$	$\frac{5}{\sqrt{4\pi }}$
$r_0\le r\le r_1$	$1+9s$	$\frac{s{u}_0}{r_0}\left(\frac{1}{2}-y\right)$	$\frac{s{u}_0}{r_0}\left(x-\frac{1}{2}\right)$	$1$	$\frac{5}{\sqrt{4\pi }}$
$r>r_1$	$1$	$0$	$0$	$1$	$\frac{5}{\sqrt{4\pi }}$

Equations

$$u_t+f(u{)}_x+g(u{)}_y=0,$$

$$\underset{\Omega }{\int }\left(u_t+f_x+g_y\right)\phi \mathrm{d}x\mathrm{d}y=0,$$

$$Q_n=\left[x_{n,1},x_{n,2}\right]\times [y_{n,1},y_{n,2}],n=1,\dots ,N_Q.$$

$$\underset{Q_n}{\int }\left(u_t+f_x+g_y\right)\phi \mathrm{d}x\mathrm{d}y=0,n=1,\dots N_Q,$$

$$X_n(\xi )=x_{n,1}+\frac{\xi +1}{2}\Delta x,Y_n(\eta )=y_{n,1}+\frac{\eta +1}{2}\Delta y,$$

$$J{u}_t+{\tilde{f}}_{\xi }+{\tilde{g}}_{\eta }=0,$$

$$J=\frac{\Delta x\Delta y}{4},\tilde{f}=\frac{\Delta y}{2}f,\tilde{g}=\frac{\Delta x}{2}g.$$

$${\phi }^{Q_n}=\sum _{i=0}^N\sum _{j=0}^N{\phi }_{ij}^{Q_n}{\psi }_i(\xi ){\psi }_j(\eta ),$$

$${\psi }_i(\xi )=\prod _{\stackrel{j=0}{j\ne i}}^N\frac{\xi -{\xi }_j}{{\xi }_i-{\xi }_j}\text{for}i=0,\dots ,N,$$

$$\underset{Q_0}{\int }\left(J{u}_t+{\tilde{f}}_{\xi }+{\tilde{g}}_{\eta }\right){\psi }_i(\xi ){\psi }_j(\eta )\mathrm{d}\xi \mathrm{d}\eta =0,$$

$$u(x,y,t){|}_{Q_n}$$

$$\tilde{f}(u(x,y,t)){|}_{Q_n}$$

$$u(x,y,t)\approx \bigcup _{n=1}^{N_Q}{u}^{Q_n}(\xi ,\eta ,t).$$

$$\begin{array}{rl}\underset{Q_0}{\int }J{u}_t^{Q_n}{\psi }_i(\xi ){\psi }_j(\eta )\mathrm{d}\xi \mathrm{d}\eta & \approx \sum _{n,m=0}^N\left(\sum _{p,q=0}^{N}J{\left(u_t^{Q_n}\right)}_{pq}{\psi }_p({\xi }_n){\psi }_q({\eta }_m)\right){\psi }_i({\xi }_n){\psi }_j({\eta }_m){\omega }_n{\omega }_m\\ & =J{\left(u_t^{Q_n}\right)}_{ij}{\omega }_i{\omega }_j,\end{array}$$

$$D_{ij}=\psi_{j}^{\prime}(\xi)\bigg{|}_{\xi=\xi_{i}},\quad i,j=0,\ldots,N.$$

$$J{\left(u_t^{Q_n}\right)}_{ij}=-$$

$$\frac{{\delta }_{Nj}}{{\omega }_N}\left[{\tilde{g}}^{\ast }({\xi }_i,1)-{\tilde{g}}_{ij}^{Q_n}\right]-\frac{{\delta }_{0j}}{{\omega }_0}\left[{\tilde{g}}^{\ast }({\xi }_i,-1)-{\tilde{g}}_{ij}^{Q_n}\right]+\sum _{m=0}^N{D}_{jm}{\tilde{g}}_{im}^{Q_n}),$$

$${\delta }_{\epsilon }^{m,k}(x)=\{\begin{array}{ll}\frac{1}{\epsilon }{P}^{m,k}\left(\frac{x}{\epsilon }\right)\text{leavevmode}\text{leavevmode}& \left|x\right|\le \epsilon \\ 0& \left|x\right|>\epsilon \end{array},$$

$$\underset{-1}{\overset{1}{\int }}{P}^{m,k}(\xi )\mathrm{d}\xi$$

$${\left(P^{m,k}\right)}^{(i)}(\pm 1)$$

$$\underset{-1}{\overset{1}{\int }}{\xi }^i{P}^{m,k}(\xi )\mathrm{d}\xi$$

$$\tilde{u}(x,t)=\underset{x-\epsilon }{\overset{x+\epsilon }{\int }}u(\tau ,t){\delta }_{\epsilon }^{m,k}(x-\tau )\mathrm{d}\tau \approx \underset{x-\epsilon }{\overset{x+\epsilon }{\int }}\left[\sum _{i=0}^N{u}_i(t){\psi }_i(\tau )\right]{\delta }_{\epsilon }^{m,k}(x-\tau )\mathrm{d}\tau =\sum _{i=0}^N{u}_i(t)\underset{x-\epsilon }{\overset{x+\epsilon }{\int }}{\psi }_i(\tau ){\delta }_{\epsilon }^{m,k}(x-\tau )\mathrm{d}\tau .$$

$${\Phi }_{ij}=\underset{x_i-\epsilon }{\overset{x_i+\epsilon }{\int }}{\psi }_j(\tau ){\delta }_{\epsilon }^{m,k}(x_i-\tau )\mathrm{d}\tau =\epsilon \underset{-1}{\overset{1}{\int }}{\psi }_j(\epsilon x+x_i){\delta }_{\epsilon }^{m,k}(\epsilon x)\mathrm{d}x\approx \epsilon \sum _{\nu =0}^{N^{\ast }}{\omega }_{\nu }{\psi }_j(\epsilon x_{\nu }+x_i){\delta }_{\epsilon }^{m,k}(\epsilon x_{\nu }),$$

$$\epsilon =\cos \left(\frac{\pi (\frac{N-N_d}{2})}{N}\right),$$

$$\tilde{\mathbf{u}}(t)=\mathbf{\Phi u}(t).$$

$$\tilde{u}(x,t)=\underset{x-\epsilon }{\overset{x+\epsilon }{\int }}u(\tau ,t){\delta }_{\epsilon }^{m,k}(x-\tau )\mathrm{d}\tau$$

$$=\sum _{\ell =1}^{N_Q}\underset{x-\epsilon }{\overset{x+\epsilon }{\int }}{u}^{Q_{\ell }}(\tau ,t){\delta }_{\epsilon }^{m,k}(x-\tau )\mathrm{d}\tau .$$

$$Q_{i,n}^{\epsilon }$$

$$Q_{i,n-1}^{\epsilon }$$

$$Q_{i,n+1}^{\epsilon }$$

$${\tilde{u}}_i^{Q_n}(t)$$

$$=\underset{Q_{i,n}^{\epsilon }}{\int }{u}^{Q_n}(\tau ,t){\delta }_{\epsilon }^{m,k}(x_i-\tau )\mathrm{d}\tau +\underset{Q_{i,n-1}^{\epsilon }}{\int }{u}^{Q_{n-1}}(\tau ,t){\delta }_{\epsilon }^{m,k}(x_i-\tau )\mathrm{d}\tau$$

$$+\underset{Q_{i,n+1}^{\epsilon }}{\int }{u}^{Q_{n+1}}(\tau ,t){\delta }_{\epsilon }^{m,k}(x_i-\tau )\mathrm{d}\tau .$$

$$E_i^{\epsilon }$$

$$E_{i,l}^{\epsilon }$$

Full text

Keywords: Discontinuous Galerkin, nonlinear hyperbolic conservation laws, SIAC filtering, Shock capturing

1. Introduction

Systems of nonlinear hyperbolic conservation laws cover a wide range of physical flow problems, e.g. modeled by the Euler or magneto-hydrodynamic (MHD) equations. Such flow phenomena are particularly interesting as they describe diverse physical processes like gas dynamics in chemical processes or plasma interactions in space, respectively [21, 23]. It is well-known that the solution of nonlinear hyperbolic systems can develop discontinuities, e.g. shocks, in finite time, regardless of the smoothness of the initial data [11]. Due to the complexity of nonlinear systems, we rely on numerical methods to approximate the solutions to such problems.

For low-order spatial approximations, like finite volume methods, discontinuous solutions are unproblematic because their inherent amount of numerical dissipation regularizes discontinuities naturally. However, for high-order numerical methods spurious oscillations near discontinuities, i.e. Gibbs phenomenon [13], arise. These unphysical overshoots might lead to unphysical solution states, e.g. negative density or pressure. Over the decades, many counter mechanisms have been developed to control overshoots and stabilize high-order approximations in shocked regions. Altogether, these methods can be subdivided into three main categories, which are slope limiters [3, 4, 12, 14, 15, 31], artificial dissipation techniques [7, 16, 17, 18, 28], and solution filters [2, 7, 27, 29].

In this work, we consider a shock capturing method that uses the global smoothness increasing accuracy conserving (SIAC) filter, a filter that has received much attention in the context of postprocessing data produced by discontinuous Galerkin approximations, see e.g. [27, 29, 35]. The SIAC filter increases smoothness by convolving the approximate solution with an appropriate smooth kernel function, e.g. B-Splines [19, 29] or Dirac-delta polynomial approximations, e.g. [34]. The accuracy conservation is more technical and is related to the fundamental building block of the discontinuous Galerkin solution space(s) and its solution ansatz [9]. The filter kernel is designed to reproduce certain polynomials orders by convolution, for example $m$. Consider a DG approximation that uses a polynomial solution space of $N$. If the filter is designed to recover a larger family of polynomial orders, i.e. $m>N$, then the filter conserves the accuracy of the approximation. If the filter recovers a smaller family of polynomial orders, $m<N$, then the accuracy of the overall approximation is bound by the filter accuracy. Typically, such SIAC filters were designed to obtain super-convergence in a post-processing step by a convolution of the numerical approximation against a specifically designed kernel function once at the final time [5, 9, 19, 26, 29, 33]. However, recent work has applied the SIAC filter as a shock capturing and/or regularization of general discontinuities strategy during the computation of the approximate solution for global spectral collocation methods [34, 37]. For such spectral methods, the SIAC filter is suitable to apply in shocked regions because said filter can recover full accuracy away from shocks, see e.g [19]. The filter for global spectral collocation methods is constructed with a Dirac-delta kernel sequence determined by two conditions that control the number of vanishing moments and the smoothness. With the discrete SIAC filter at hand, the shock regularization of a global approximation is then performed by a convolution of the solution with the high-order Dirac-delta kernels in every time step. In practice, the filtering procedure reduces to a simple matrix-vector multiplication and, thus, allows for an efficient and simple implementation.

The main contribution of this work is to extend the global filtering technique to element-wise approximations within a nodal discontinuous Galerkin (DG) method. We exploit the weak coupling of the discontinuous Galerkin spectral element method (DGSEM) at element interfaces to design a multi-element filter. Consequently, we convolve the polynomial approximations of one element and its nearest neighbor’s solutions with Dirac-delta kernels instead of the global representation of the solution. As we will point out in the derivation of this multi-element filter, it also recasts to locally applied matrix-vector formulations. We present the filtering matrices for one-dimensional and two-dimensional Cartesian DG discretizations. Moreover, we construct the multi-element SIAC filter such that the numerical scheme remains high-order accurate in smooth regions and that the shock capturing is applied adaptively throughout the domain. The latter we obtain by implementing a shock indicator, which is defined by the filtered solution itself. According to this indicator, we replace the DG approximation by the filtered one in oscillatory regions and, additionally, even introduce a smooth transition area, in which we use a convex combination of the filtered and unfiltered solutions.

The outline of this work is as follows: We begin in Sec. 2 with a construction of the DG method on two-dimensional Cartesian meshes. Next, in Sec. 3 we provide the SIAC filtering routines, where we first discuss the global filter, before we extend it to the multi-element variant as well as two spatial dimensions. Furthermore, we broach the issue of adaptive filtering and conservation properties in the same section. Finally, we provide several numerical benchmark tests in order to verify the applicability of the novel filter to shock problems for the two-dimensional Euler and ideal MHD equations in Sec. 4. Lastly, Sec. 5 gives concluding remarks and an outlook on possible further research projects.

2. Discontinuous Galerkin spectral element method

Throughout this work we consider the solution of hyperbolic conservation laws in two spatial dimensions which take the general form

[TABLE]

in a square domain $\Omega=\left[x_{L},x_{R}\right]\times\left[y_{B},y_{T}\right]\subset\mathbb{R}^{2}$ with appropriate boundary conditions and an initial solution $u(x,y,0)=u_{0}(x,y)$. Here $u$ is a conserved variable and $f,g$ are the nonlinear fluxes. We take (2.1) as the prototype equation for the conserved solution variables such as density or momentum. This simplifies the discussion and derivations for the SIAC filtering in Sec. 3. The discussion extends naturally to a system of hyperbolic conservations laws such as the Euler equations.

We first provide an overview for the derivation of the nodal DGSEM on Cartesian meshes. Complete details can be found in the book by Kopriva [20]. We derive the DGSEM from the weak form of the conservation law (2.1). As such, we multiply by an arbitrary discontinuous $L_{2}(\Omega)$ test function $\varphi$ and integrate over the domain

[TABLE]

where we suppress the ${u}$ dependence of the nonlinear fluxes. We subdivide the domain $\Omega$ into $N_{Q}$ non-overlapping quadratic elements

[TABLE]

For the present discussion we make the simplifying assumption that all elements have the same size, i.e. $\Delta x=x_{n,2}-x_{n,1}$and $\Delta y=y_{n,2}-y_{n,1}$ for all $n=1,\ldots,N_{Q}$. This divides the integral over the whole domain into the sum of the integrals over the elements. So, each element contributes

[TABLE]

to the total integral. Next, we create a transformation between the reference element $Q_{0}=[-1,1]^{2}$ and each element, $Q_{n}$. For rectangular meshes we create mappings $(X_{n},Y_{n}):Q_{0}\rightarrow Q_{n}$ such that $\left(X_{n}(\xi),Y_{n}(\eta)\right)=(x,y)$ are

[TABLE]

for $n=1,\ldots,N_{Q}$. Under the transformation (2.5) the conservation law in physical coordinates (2.1) becomes a conservation law in reference coordinates [20]

[TABLE]

where

[TABLE]

We select the test function $\varphi$ to be a piecewise polynomial of degree $N$ in each spatial direction

[TABLE]

on each spectral element $Q_{n}$, but do not enforce continuity at the element boundaries. The interpolating Lagrange basis functions are defined by

[TABLE]

with a similar definition in the $\eta$ direction. The values of $\varphi^{Q_{n}}_{ij}$ on each element $Q_{n}$ are arbitrary and linearly independent, therefore the formulation (2.4) is

[TABLE]

where $i,j=0,\ldots,N$.

For the DGSEM we approximate the conservative variable $u$ and the contravariant fluxes $\tilde{f},\tilde{g}$ with polynomial interpolants of degree $N$ in each spatial direction written in Lagrange form on each element $Q_{n}$, e.g.,

[TABLE]

This implies that the global representation of the solution $u$ is the union of these piecewise polynomials

[TABLE]

Next, we use integration-by-parts to move derivatives off the nonlinear fluxes and onto the test function, which generates surface and volume contributions. We resolve the discontinuities between elements at the surface by introducing Lax-Friedrichs numerical flux functions $f^{*}$ and $g^{*}$. We apply integration-by-parts a second time to move derivatives back onto the fluxes. For the nodal DGSEM any integrals present in the variational formulation are approximated with $N+1$ Legendre-Gauss-Lobatto (LGL) quadrature nodes and weights, e.g.,

[TABLE]

for each element $n=1,\ldots,N_{Q}$ and where $\left\{\xi_{i}\right\}_{i=0}^{N},\left\{\eta_{j}\right\}_{j=0}^{N}$ are the LGL quadrature nodes and $\left\{\omega_{i}\right\}_{i=0}^{N},\left\{\omega_{j}\right\}_{j=0}^{N}$ are the LGL quadrature weights. Further, we collocate the interpolation and quadrature nodes which enables us to exploit that the Lagrange basis functions (2.9) are discretely orthogonal and satisfy the Kronecker delta property, i.e., $\psi_{j}(\xi_{i})=\delta_{ij}$ with $\delta_{ij}=1$ for $i=j$ and $\delta_{ij}=0$ for $i\neq j$ to simplify (2.13). Due to the polynomial approximation (2.11) any derivatives fall on the Lagrange basis functions. These are approximated at high-order with the standard differentiation matrix [20]

[TABLE]

From these steps, we build the standard semi-discrete formulation of the strong-form DGSEM. We write the scheme in index notation as

[TABLE]

where $i,j=0,\ldots,N$ and $n=1,\ldots,N_{Q}$.

The semi-discrete formulation (2.15) on each element is integrated in time with an explicit five stage, fourth order Runge-Kutta method of Carpenter and Kennedy [6]. We select a stable explicit time step with an appropriate CFL condition which is equation and resolution dependent.

3. SIAC Filter

In this section we present the SIAC filter for a single domain and then discuss its extension and implementation into a multi-element DGSEM framework.

3.1. Single element filter

In [37] a SIAC filter was developed for a single domain spectral method. To define the global method we begin by introducing the delta sequence

[TABLE]

that is built from the polynomial $P^{m,k}(x)$, which is uniquely determined by the following conditions:

[TABLE]

In Fig. 1 we illustrate the polynomial approximation of the delta kernel (3.1) with $(m,k)=(3,8)$, $(m,k)=(5,8)$, and scaling parameter $\varepsilon=1$.

According to the SIAC filtering strategy [29, 34, 37] we regularize the solution produced by the numerical scheme with the manipulation

[TABLE]

For compact notation we introduce the filter matrix $\Phi$ and approximate its values with LGL quadrature by mapping the corresponding integration area $\left[x-\varepsilon,x+\varepsilon\right]$ into the reference element $E=\left[-1,1\right]$

[TABLE]

where $\left\{x_{\nu}\right\}_{\nu=0}^{N^{\ast}}$ and $\left\{\omega_{\nu}\right\}_{\nu=0}^{N^{\ast}}$ are the LGL quadrature points and weights for $N^{\ast}=2\left(\frac{m}{2}+k+1\right)$ to maintain the desired high-order accuracy of the approximation [34]. Moreover, we choose

[TABLE]

with $N_{d}$ determined empirically to ensure stable conserving results [34, 37]. The selection of the scaling parameter $\varepsilon$ is related to the quadrature accuracy of the integral (3.6). As was shown in [34], an arbitrary choice of $\varepsilon$ can lead to sub-optimal accuracy of the quadrature. However, depending on the number of LGL nodes used for the DG approximation the support of the kernel function (3.7) must be adjusted through the parameter $N_{d}$. In practice, for a fixed number of LGL nodes $N+1$, selecting different values of $\varepsilon$ changes the accuracy of the filter as well as the solution quality because of possible excessive smearing effects, see [34, 37] for details.

We can express the filtering process in terms of a matrix vector multiplication, i.e.

[TABLE]

For the global SIAC filtering technique we must address how the filter matrix is applied at the physical boundaries of the domain. However, at the physical boundaries no $\varepsilon$-stencils are defined. Thus, oscillations caused by shocks as well as by re-interpolation (Runge phenomena) cannot be smoothed in these areas. In the original approach for the global collocation the affected parts of the discretization are set to the analytical solution [37]. Using a local version of the filter we can avoid identifying interior points by an analytical reference solution, as discussed in the next section.

We also note that, by construction, the filter conserves mass solely for polynomial data of degree up to $m$, which especially in a global collocation method is difficult to realize. Thus, small conservation errors might be introduced by applying the filter matrix (3.6).

3.2. Multi-element filter

For the case with multiple Cartesian elements and DGSEM, we begin, again, with the one-dimensional case and then apply the tensor product decoupling of the DGSEM to move to higher spatial dimensions. In contrast to the previous section we now have DG solutions defined on multiple elements in a mesh ($N_{Q}>1$). However, we still want to apply the smoothing matrix $\Phi$ locally to the solution element-wise. It is, hence, necessary to determine how to couple the filtering across element interfaces to determine a multi-element SIAC technique.

To do so, we begin with (3.5), where, in one spatial dimension, we know that the solution $u$ is a union of piecewise polynomials over all elements $N_{Q}$

[TABLE]

Next, we focus on one physical node $x_{i}:=X_{n}(\xi_{i})$ within one element $Q_{n}$ and define the following sets

[TABLE]

Since $\varepsilon$ is sufficiently small, we assume that the $\varepsilon$-stencil is imbedded in these three sets and thus

[TABLE]

If we now define similar sets for the corresponding LGL node

[TABLE]

we can transform everything to reference space again and obtain

[TABLE]

Note, that we shift the arguments of the lagrange basis functions in the left and right elements by $\pm 2$ to guarantee the correct evaluation points. We can write (3.17) in compact notation by applying a modified $(N\!+\!1)\times 3(N\!+\!1)$ smoothing matrix to the solution, i.e.

[TABLE]

with the matrices

[TABLE]

Again, we evaluate these integrals by a Legendre-Gauss-Lobatto quadrature with $N^{\ast}=2\left(\frac{m}{2}+k+1\right)$ points as in (3.6). Because the two dimensional, multi-element filter is created through a tensor product of the one-dimensional filters, we select the same $\varepsilon$ (3.7) in each of the integrals (3.19)-(3.21).

Note, that since the neighboring elements only enter in the smoothing matrix at grid points near to the element boundaries, $\mathbf{\Phi}^{n-1}$ and $\mathbf{\Phi}^{n+1}$ are block matrices with mostly zero entries, especially when $N$ is large. In particular, only the first several rows of $\mathbf{\Phi}^{n-1}$ and the last several rows of $\mathbf{\Phi}^{n+1}$ are non-zero. We see that the multi-element SIAC filtering process is not entirely local to element $Q_{n}$; however, we only need solution information from its direct neighbors in the mesh.

An additional advantage in the design of this multi-element filtering technique is the treatment of element located at physical boundaries. We already noted that there is no $\varepsilon$-stencil defined in these boundary areas. Thus, we cannot apply the filter. However, from the multi-element technique we can introduce ghost elements, in which we can define a consistent solution depending on the physical boundary condition, e.g., to reflecting wall or Dirichlet values. This procedure removes Runge phenomena from the solution without the need to identify interior points by analytical values [37].

We can use the locally filtered solution as a shock detector to adaptively apply the multi-element filter only in elements where it is necessary. To do so, we define an indicator to measure the difference between the filtered and unfiltered solutions

[TABLE]

Next, we normalize this indicator with respect to the polynomial order and the number of elements and check in each element $n=1,\ldots,N_{Q}$, if

[TABLE]

for a given user defined tolerance TOL. If this condition is fulfilled, we replace the current element solution with the filtered solution. Otherwise the approximation is deemed to be sufficiently smooth and no filtering is applied. In order to extend this adaptation for systems of conservation laws we can use single variables to compute $\epsilon_{n}$, e.g., the density or pressure for the Euler equations.

Remark 1* (Adaptive filtering).*

The filtered solution acts as a self-contained shock detector because of the constraints used to construct the SIAC filter (3.2)-(3.4). In particular, the filter is designed to recover polynomial orders up to $m$. Therefore, in smooth regions of the flow the approximate solution and the filtered solution will be nearly the same, i.e., the indicator error (3.23) will be small. However, near discontinuities the approximate solution will contain large, spurious overshoots whereas the overshoots of the filtered solution will be considerably reduced, making (3.23) large.

Moreover, for convenience, we define

[TABLE]

and introduce a transition area between two tolerance levels, $\sigma_{\min}\leq\sigma_{n}\leq\sigma_{\max}$, to smoothly blend the filtered and unfiltered solutions. As such we define a parameter $0\leq\lambda\leq 1$ and then define the updated solution on a given element $Q_{n}$ to be a convex combination of the two solutions

[TABLE]

with

[TABLE]

A major concern for any shock capturing method is to maintain conservation, which ensures the correct shock speeds are maintained [25].

Remark 2* (Conservation).*

In its current incarnation the multi-element SIAC filter does not conserve the solution quantities, e.g., density, momentum and energy for the Euler equations. The unfiltered standard DGSEM conserves the solution variables up to machine precision, see e.g. [8]. However, the application of the local SIAC filter after each time step is no longer globally conservative because we re-distribute solution data, e.g. the mass, by the filtering process within each element. Whereas the conservation errors for the global filter are introduced by the necessary large interpolation order $N\gg m$, we can easily assure $N\leq m$ for the local element-wise filtering. However, we run into a different problem because our global approximation is no longer a polynomial, but solely built from piecewise polynomial data. Thus, again, we introduce conservation errors in our approximation, which are usually small. We examine the size of these conservation errors in Sec. 4.1.1 and show for the considered test cases that the conservation loss does not affect the solutions.

3.3. Two-dimensional filter

Next, we extend the one-dimensional local SIAC filter to higher spatial dimensions. For the filtering process we apply the same local smoothing matrix $\hat{\mathbf{\Phi}}$ as in the one-dimensional case in each spatial direction to the unfiltered element solutions. Conveniently, this is possible due to the tensor product ansatz of the DGSEM and the definition of the Dirac delta kernel (3.1). Just as in the previous section we first derive a global multi-dimensional filter and then modify it for the local multi-element case.

We begin again from the filtering assumption of (3.5), and find for the piecewise polynomial solution $u$ that

[TABLE]

where we define the multi-variable delta function to have the form

[TABLE]

We focus on one LGL node, transform into the reference space and use the tensor product property to split the integrand and obtain

[TABLE]

Here, the $\bar{\ell}$ points to the correct solution entry, which includes neighboring elements and is dependent on the storing data structure. The shifting of the evaluation points for the lagrange basis function is also hidden in the bar notation, i.e.

[TABLE]

From this definition of the filtering matrices it is possible to write the filtering matrix in a compact notation

[TABLE]

with

[TABLE]

provided the elements are labeled from bottom-left to top-right and $N_{Q,x}$ denotes the number of elements in the $x$-direction. In this case, we design the smoothing matrix $\mathbf{\Phi}=\begin{pmatrix}\mathbf{\Phi}^{n-1}&\mathbf{\Phi}^{n}&\mathbf{\Phi}^{n+1}\end{pmatrix}$ exactly as in one spatial dimension.

We want the resulting shock capturing DG scheme to be as local as possible and implement the 2D multi-element SIAC filter in a way to reflect this goal. First, we define

[TABLE]

which is nothing more than the solution vector of the three considered adjacent cells filtered in the $x$-direction. To filter in the $y$-direction, we just apply the smoothing matrix $\mathbf{\Phi}$ from the left hand side and obtain an overall filtered solution $\tilde{\mathbf{u}}^{n}$ from (3.31). The main advantage of this procedure is that the filtering procedure is done dimension by dimension. So for all elements $n=1,\ldots,N_{Q}$ we first filter in the $x$-direction to find $\hat{\mathbf{u}}^{n}$ with coupling only from the right and left neighbor cells. Next, we filter in the $y$-direction and compute the fully filtered solution $\tilde{\mathbf{u}}^{n}$ from the information stored in the intermediate array $\hat{\mathbf{u}}^{n}$ with coupling from the upper and lower neighbor elements. That is, we simply apply the one-dimensional filter twice for each grid point and, again, only need information from the direct neighbors. Additionally, we note, that this filtering procedure has no preferred direction such that the order of $x,y$ directions makes no difference.

4. Numerical tests

We verify the performance of the novel multi-element SIAC filter applied to a variety of two-dimensional shock tests for the Euler as well as ideal MHD equations. For all simulations, we consider two-dimensional Cartesian meshes discretized by uniform quadrilateral elements of equal size $\Delta x=\Delta y$. Further, we use the explicit five stage, fourth order Runge-Kutta method of Carpenter and Kennedy [6] with an stable explicit time step to advance the DG approximation in time. The (adaptive) filtering procedure (3.31) is performed for all element solutions after each time step. We begin with the numerical validation for the two-dimensional Euler equations in Sec. 4.1, where we also investigate the accuracy and conservation issues of the filter, before we apply it to more challenging shock tests for the ideal MHD equations in Sec. 4.2.

4.1. Euler tests

The two-dimensional Euler equations are described by a system of conservation laws, i.e.

[TABLE]

with

[TABLE]

Here, $\varrho$, $\vec{v}=(v_{1},v_{2})^{T}$, and $e$ denote the density, two-dimensional velocity field and inner specific energy, respectively. We close the system by the perfect gas equation, which relates the inner energy and pressure, i.e.

[TABLE]

where $\gamma$ denotes the adiabatic coefficient.

In the following, we apply the DGSEM with the multi-element SIAC filter derived herein to the two-dimensional Euler equations (4.1) to first show the high-order convergence and investigate the conservations properties of the scheme. Thereupon, we consider several benchmark problems for the two-dimensional Euler equations in order to verify the shock capturing capabilities of the novel filtering strategy.

4.1.1. Convergence and conservation studies

A substantial property of DG schemes is the high-order accuracy of the approximation for smooth solutions. Thus, we first consider an academic test case with a known analytical solution for the two-dimensional Euler equations, that allows us to compute numerical errors measured in a discrete $L_{\infty}$-norm. With the help of these error values for different discretization levels, we are able to compute the experimental order of convergence (EOC), which for the DGSEM is expected to agree with the theoretical order of $N\!+\!1$ as the mesh is refined.

The problem we use for convergence tests is defined in the domain $\Omega=[-1,1]^{2}$ with the initial conditions

[TABLE]

We use periodic boundary conditions and set $\gamma=5/3$. The analytical solution of the 2D Euler equations using these initial conditions is

[TABLE]

We run the simulation until the final time $T=0.4$ and intentionally choose a high polynomial degree of $N=7$. Consequently, we select a small explicit time step where $\texttt{CFL}=0.1$ to exclude errors from the time integration method. Further, we compute the errors of the approximation for different choices of $N_{Q}$ and calculate the EOC by the maximum error $\epsilon_{\infty}$ of the density. For the DGSEM approximation without filtering we obtain the convergence results illustrated in Table 1.

The results in Table 1 confirm the expected theoretical order of convergence $N+1$ for the DGSEM approximation without filtering. Moreover, we state the overall conservation error of the density computed by

[TABLE]

which regard to machine precision and, thus, agree with the desired conservative nature of the DGSEM [8].

As we are interested in the errors and convergence rates of the filtered solution as well, we turn off the adaptive filtering and investigate the convergence rates and conservation errors of the filtered solution as demonstrated in Table 2 for the SIAC filter using a Dirac-delta approximation with one vanishing moment.

Table 2 reveals that the order of convergence drops to one. Further, we lose conservation as pointed out in the previous section. We now repeat the convergence test with $m=3$ and $m=5$, see Tables 3 and 4.

As the convergence tests show, the experimental order of convergence depends on the number of vanishing moments in the Dirac-delta approximation, i.e. $EOC\approx{\min}\left(m,N+1\right)$. While one vanishing moment leads to a smearing effect and significantly drops the accuracy in smooth areas, increasing the number of vanishing moments improves the accuracy again. We observe that the conservation error decreases with an increasing number of vanishing moments as well as with finer grid resolutions.

In conclusion, we expect the Dirac-delta filter to effectively distinguish between shocks and smooth areas when using enough vanishing moments. Because of the large drop of accuracy for less vanishing moments, smooth areas will be mistaken for shocks more often which results in significantly less accuracy in those areas. We alleviate these issues with the adaptive filtering technique, so that small conservation errors are exclusively introduced in shocked regions.

4.1.2. Explosion problem

We now consider the first two-dimensional Euler test case that inherits shocks. The explosion problem is defined on the domain $\Omega=[-1,1]^{2}$, whereas its initial conditions consist of a region inside a circle with the radius $r=0.4$ centered at the origin and a region outside that circle, see e.g. [37]. The primitive variable initial conditions are then defined by two states, i.e.

[TABLE]

Here, the inner state applies for every $(x,y)\in\Omega$ such that $(x^{2}+y^{2})\leq r^{2}$ and the outside state applies otherwise. Further, the velocity vector is set to zero in the entire domain and we set $\gamma=5/3$.

The following plots in Figures 2 and 3 show the approximation of the density for this problem at the final time $T=0.25$ with $\texttt{CFL}=0.1$ and a polynomial degree of $N=7$ on $N_{Q}=80^{2}$ elements.

For the simulation result illustrated in Fig. 2, one vanishing moment ($m=1$) and six continuous derivatives at the endpoints ($k=6$) of the Dirac-delta approximation are chosen. The support width in (3.7) is calculated with $N_{d}=0.6$. Furthermore, the adaptive filtering with convex blending (3.25) is used with a minimum tolerance of $10^{-7}$ and a maximum tolerance of $10^{-3}$. Comparing the results to other configurations of the filter shows that the narrow domain width and the chosen tolerances effectively limit the smearing effect despite the small number of vanishing moments. The conservation error is $\epsilon_{\mathrm{cons}}=6.4\cdot 10^{-4}$.

In contrast, Fig. 3 uses three vanishing moments ($m=3$), which requires a larger support width defined by $N_{d}=2.5$ in (3.7). The corresponding conservation error is $\epsilon_{\mathrm{cons}}=2.0\cdot 10^{-5}$.

Both results reveal that the shocks are effectively regularized, even though small overshoots are observable. For this problem, we find that one vanishing moment is enough to obtain acceptable results if the tolerance of the adaptive filtering operation is adjusted. However, the three moment adaptive SIAC filter produces more accurate shock profiles, as can be seen in the slicing picture Fig. 4 that compares three filtered solutions against a reference solution [37]. The reference solution is computed by the publicly available high performance application code FLASH (http://flash.uchicago.edu/site/flashcode/) with a second order MUSCL-Hancock finite volume method (see e.g. [36]) on $2048\times 2048$ elements. We observe that the filter with one vanishing moment is quite dissipative, whereas the filtered approximations with higher moments produce small overshoots at the shock and rarefraction. Nonetheless, all three configurations are stable, nearly oscillation-free and close to the analytical profile.

4.1.3. Four state Riemann test

The next Euler shock test is described by the 2D-Riemann problem, see [22] and [24]. We consider the domain $\Omega=[0,1]^{2}$ with outflow boundary conditions. The initial conditions are defined by four different states in four quadrants:

**Top left:: **

$\Omega_{t\ell}=[0,0.5)\times(0.5,1]$.

**Bottom left:: **

$\Omega_{b\ell}=[0,0.5)\times(0,0.5]$.

**Top right:: **

$\Omega_{tr}=(0.5,1]\times(0.5,1]$.

**Bottom right:: **

$\Omega_{br}=(0.5,1]\times[0,0.5)$.

The four primitive variable states of the initial conditions are then assigned to these four quadrants:

[TABLE]

Depending on the particular choice of initial conditions, the simulation has to overcome at least one rarefaction wave, shock wave or contact wave. In [22], 19 particular configurations for initial values are given. Because configurations 17 and 19 deal with all of the mentioned waves, we apply the Dirac-delta filter to them. The results of the following configurations are obtained by using a polynomial degree of $N=7$ on a grid of $60^{2}$ elements with $\texttt{CFL}=0.1$.

Configuration 17:

The four initial conditions of configuration 17 are defined for the primitive variables as

[TABLE]

We run the approximation to $T=0.3$. The pseudo-color plot of the density in Fig. 5 is obtained by using the adaptive filter with a Dirac-delta approximation of $(m,k)=(5,7)$. The support width $\varepsilon$ is calculated by $N_{d}=4.5$ in (3.7), where the tolerances are set to $\sigma_{\text{min}}=-8$ and $\sigma_{\text{max}}=-3$. As illustrated in the plot, the filter ensures a successful regularization of shocks while keeping a high resolution of the vortex.

Configuration 19:

The initial conditions for the primitive variables of configuration 19 are defined as

[TABLE]

Again, we filter adaptively after every time step with a Dirac-delta approximation defined by $m=5,k=7$ and a support width calculated using $N_{d}=4.5$. Figure 6 demonstrates the approximation for the density and the pressure at the final time $T=0.3$, where the adaptive filtering again has minimum and maximum tolerances of $\sigma_{\text{min}}=-8$ and $\sigma_{\text{max}}=-3$. As before, the filter effectively regularizes shocked areas. Additionally, the adaptive filtering step ensures high resolution of the curly area in the center of the domain.

4.1.4. Double Mach reflection

As a final Euler test case, we consider the double Mach reflection problem, see e.g. [30]. This challenging problem is defined on the domain $\Omega=[0,3.25]\times[0,1]$ and involves both, strong shock interactions and different boundary conditions. An initial shock separates a left and right state defined by the primitive variables, i.e.

[TABLE]

In particular, the shock is initially located at $x=1/6,y=0$ along a linear slope, which includes an angle of $\alpha=\pi/3$ with the $x$-axis, defining the initial conditions

[TABLE]

Thus, the boundary conditions at the bottom for $x<1/6$ as well as at the left and upper domain boundaries are set to Dirichlet-values all through the simulation. Especially, the latter have to be implemented correctly by an analytical shock profile, since the shock moves forward and its position is determined by the following function

[TABLE]

Using this function, the values at the top boundary are set to

[TABLE]

Furthermore, the right boundary conditions are outflow and the bottom is considered a reflecting wall for $x>1/6$. Hence, in the filtering process the according ghost cell values are set constantly either to the analytical solution for Dirichlet or to the last interior value for outflow and reflection.

In our simulation of this problem, we use a high grid resolution of $325\times 100$ elements, a polynomial degree of $N=7$ and a CFL number of $\texttt{CFL}=0.1$. We provide a density plot at the final time $T=0.2$ in Figure 7, where we applied the adaptive multi-element SIAC filter with $m=3,k=6,N_{d}=2.5$ and the according tolerances $\sigma_{\text{min}}=-6,\sigma_{\text{max}}=-2$. Compared to the previous shock tests, the number of vanishing moments had to be reduced in order to get smoother results. With these particular configurations, the approximation shows the curved reflected shock as well as the formation of a turbulent vortex, which compare well to the simulation results in the literature, e.g. [32].

4.2. Ideal MHD tests

Next, we consider the ideal magneto-hydrodynamic (MHD) equations, which can be expressed in a compact form as a system of conservation laws

[TABLE]

with

[TABLE]

Here, $\vec{v}=(v_{1},v_{2},v_{3})^{T}$ denotes the velocity field and $\vec{B}=(B_{1},B_{2},B_{3})^{T}$ the magnetic field components. As for the Euler equations, the system is closed by the perfect gas equation, which for the ideal MHD equations reads

[TABLE]

where $\gamma$ again denotes the adiabatic coefficient.

For the ideal MHD equations it is well known, that magnetic monopoles are not observed in nature, but due to round-off errors, numerically $\nabla\cdot\vec{B}=0$ is not necessarily guaranteed. Hence, to maintain the divergence-free condition of the magnetic field $\nabla\cdot\vec{B}=0$ we implement a hyperbolic divergence cleaning method based on generalized Lagrangian multiplier [10].

In the following we demonstrate the performance of the local SIAC Dirac-delta filter for two common benchmark problems of the two-dimensional MHD simulations including (strong) shocks.

4.2.1. Orszag-Tang vortex

The first test case describes the evolution of a turbulent plasma cloud, the well-known Orszag-Tang vortex [1], which is initialized in the periodic domain $\Omega=\left[0,1\right]^{2}$ as

[TABLE]

The other variables are initially set to zero and $\gamma=5/3$. Furthermore, we use $\texttt{CFL}=0.5$, polynomials of degree $N\!=\!5$ and $40\times 40$ elements for the following simulations of this test case.

We show the evolution of the density in the plots below (Fig. 8) smoothed by the multi-element SIAC filter with $m=3,k=8$ and a fixed $\varepsilon=1.4$. As for the Euler tests, we apply the filtering adaptively as shown in (3.25) with the pressure as a shock indicator and $\sigma_{\min}=\sigma_{\max}=-8$. Further, we show the distribution of the cell-wise constant convex parameter $\lambda$ from (3.25) at the same stages in Fig. 9, which confirms the correct tracking of shocks as they evolve.

In order to assess the performance of the two-dimensional Dirac-delta filter, we compare the simulation results to a reference solution obtained by the publicly available high performance application code FLASH (http://flash.uchicago.edu/site/flashcode/). Particularly, this highly resolved reference solution is computed by a second order MUSCL-Hancock finite volume method (see e.g. [36]) on $1024\times 1024$ elements. Further, we apply the filter once more to the Orszag-Tang vortex with different parameters, i.e. we choose $m=1,k=6,\varepsilon=1.6,\sigma_{\min}=-9$ and $\sigma_{\max}=-6$. In Fig. 10 we provide two slices of the Orszag-Tang vortex at the final time $T=0.5$, in which we cut through the density distribution at $x=y$ (left) and $y=0.3$ (right) to compare the profiles obtained by the multi-element SIAC filter against a reference solution.

We see that the oscillations are smoothed out by the multi-element SIAC filter and the approximation matches the reference solution quite well, but the filtering technique still produces little overshoots at shocks. However, taking into account the coarse resolution of $40\times 40$ elements, the filter generates reasonable approximations for this shock test, as it stabilizes the approximation and regularizes it against spurious oscillations.

4.2.2. Magnetic rotor

The second test case describes a rotating dense circle in a static fluid, that generates strong circular shock waves [1]. In general this benchmark problem is defined in the same periodic domain $\Omega=\left[0,1\right]^{2}$ by the radius $r=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}$ and the slope $s=\frac{r_{1}-r}{r_{1}-r_{0}}$. The initial primitive variables for the magnetic rotor are stated in Table 5, where the unlisted quantities are initially zero in the entire domain and $\gamma=1.4$.

In our simulations we define $r_{0}=0.1,r_{1}=0.115$ and $u_{0}=2$. We use $\texttt{CFL}=0.5$, a polynomial degree of $N\!=\!4$ and $100\times 100$ elements. We show the density and pressure at $T=0.15$ in Figure 11. Due to the strong circular shocks combined with the high-order DG approximation this test case is extremely sensitive and unstable. Therefore, we apply the SIAC filtering matrix constructed by a Dirac-delta kernel with only one vanishing moment $m=1$ and $k=5$. Again, we smooth the approximation adaptively with the density as a shock indicator, $\sigma_{\min}=-9,\sigma_{\max}=-6$ and a fixed $\varepsilon=1.4$.

In Figure 11 we see, that the local SIAC filter performs well in terms of stabilizing the approximation and regularizing oscillatory regions, but is polluted by small mesh artifacts, which are particularly visible at the generated Alfvén waves in the pressure profile. Nonetheless, the novel filter produces reasonable approximations even for such a challenging test problem including strong circular shock waves.

5. Conclusion

We have presented a novel shock capturing technique for DG methods based on multi-element SIAC filtering. In particular, we have first introduced the DGSEM on two-dimensional Cartesian meshes, before we have derived the local filter matrix from the global SIAC filtering approach constructed by approximations of Dirac-delta kernels. Moreover, we have designed an adaptive filtering strategy and extended the overall method to higher spatial dimensions. Finally, we have verified the applicability of the filter to a variety of challenging shock problems for both, the two-dimensional Euler and the ideal MHD equations.

We have demonstrated that the multi-element SIAC filter indeed performs well in terms of regularizing oscillatory regions and stabilizing the approximation. However, the main issue of the constructed filter is the introduction of spurious, albeit small, conservation errors, which might be avoided by projections onto discontinuous basis functions across cell interfaces. This, together with the investigation of the entropic properties as well as the extension of the local SIAC filter to curvilinear elements are future research projects.

Acknowledgements

Gregor Gassner has been supported by the European Research Council (ERC) under the European Union’s Eights Framework Program Horizon 2020 with the research project Extreme, ERC grant agreement no. 714487. Professor Jacobs acknowledges funding from NSF-DMS 1115705.