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 . 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 , 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. . 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 . The filter kernel is designed to reproduce certain polynomials orders by convolution, for example . Consider a DG approximation that uses a polynomial solution space of . If the filter is designed to recover a larger family of polynomial orders, i.e. , then the filter conserves the accuracy of the approximation. If the filter recovers a smaller family of polynomial orders, , 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  . 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 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
in a square domain with appropriate boundary conditions and an initial solution . Here is a conserved variable and 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 . We derive the DGSEM from the weak form of the conservation law (2.1). As such, we multiply by an arbitrary discontinuous test function and integrate over the domain
where we suppress the dependence of the nonlinear fluxes. We subdivide the domain into non-overlapping quadratic elements
For the present discussion we make the simplifying assumption that all elements have the same size, i.e. and for all . This divides the integral over the whole domain into the sum of the integrals over the elements. So, each element contributes
to the total integral. Next, we create a transformation between the reference element and each element, . For rectangular meshes we create mappings such that are
We select the test function to be a piecewise polynomial of degree in each spatial direction
on each spectral element
, but do not enforce continuity at the element boundaries. The interpolating Lagrange basis functions are defined by
with a similar definition in the direction. The values of on each element are arbitrary and linearly independent, therefore the formulation (2.4) is
For the DGSEM we approximate the conservative variable and the contravariant fluxes with polynomial interpolants of degree in each spatial direction written in Lagrange form on each element , e.g.,
This implies that the global representation of the solution is the union of these piecewise polynomials
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 and . 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 Legendre-Gauss-Lobatto (LGL) quadrature nodes and weights, e.g.,
for each element and where are the LGL quadrature nodes and 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., with for and for 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 
From these steps, we build the standard semi-discrete formulation of the strong-form DGSEM. We write the scheme in index notation as
where and .
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  a SIAC filter was developed for a single domain spectral method. To define the global method we begin by introducing the delta sequence
that is built from the polynomial , which is uniquely determined by the following conditions:
For compact notation we introduce the filter matrix and approximate its values with LGL quadrature by mapping the corresponding integration area into the reference element
where and are the LGL quadrature points and weights for to maintain the desired high-order accuracy of the approximation . Moreover, we choose
with determined empirically to ensure stable conserving results [34, 37]. The selection of the scaling parameter is related to the quadrature accuracy of the integral (3.6). As was shown in , an arbitrary choice of 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 . In practice, for a fixed number of LGL nodes , selecting different values of 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.
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 -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 . 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 , 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 (). However, we still want to apply the smoothing matrix 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 is a union of piecewise polynomials over all elements
Next, we focus on one physical node within one element and define the following sets
Since is sufficiently small, we assume that the -stencil is imbedded in these three sets and thus
If we now define similar sets for the corresponding LGL node
we can transform everything to reference space again and obtain
Note, that we shift the arguments of the lagrange basis functions in the left and right elements by to guarantee the correct evaluation points. We can write (3.17) in compact notation by applying a modified smoothing matrix to the solution, i.e.
with the matrices
Again, we evaluate these integrals by a Legendre-Gauss-Lobatto quadrature with 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 (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, and are block matrices with mostly zero entries, especially when is large. In particular, only the first several rows of and the last several rows of are non-zero. We see that the multi-element SIAC filtering process is not entirely local to element ; 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 -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 .
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
Next, we normalize this indicator with respect to the polynomial order and the number of elements and check in each element , if
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 , 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 . 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
and introduce a transition area between two tolerance levels, , to smoothly blend the filtered and unfiltered solutions. As such we define a parameter and then define the updated solution on a given element to be a convex combination of the two solutions
A major concern for any shock capturing method is to maintain conservation, which ensures the correct shock speeds are maintained .
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. . 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 , we can easily assure 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 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 that
where we define the multi-variable delta function to have the form
We focus on one LGL node, transform into the reference space and use the tensor product property to split the integrand and obtain
Here, the 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.
From this definition of the filtering matrices it is possible to write the filtering matrix in a compact notation
provided the elements are labeled from bottom-left to top-right and denotes the number of elements in the -direction. In this case, we design the smoothing matrix 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
which is nothing more than the solution vector of the three considered adjacent cells filtered in the -direction. To filter in the -direction, we just apply the smoothing matrix from the left hand side and obtain an overall filtered solution from (3.31). The main advantage of this procedure is that the filtering procedure is done dimension by dimension. So for all elements we first filter in the -direction to find with coupling only from the right and left neighbor cells. Next, we filter in the -direction and compute the fully filtered solution from the information stored in the intermediate array 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 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 . Further, we use the explicit five stage, fourth order Runge-Kutta method of Carpenter and Kennedy  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.
Here, , , and 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.
where 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 -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 as the mesh is refined.
The problem we use for convergence tests is defined in the domain with the initial conditions
We use periodic boundary conditions and set . The analytical solution of the 2D Euler equations using these initial conditions is
We run the simulation until the final time and intentionally choose a high polynomial degree of . Consequently, we select a small explicit time step where to exclude errors from the time integration method. Further, we compute the errors of the approximation for different choices of and calculate the EOC by the maximum error 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 for the DGSEM approximation without filtering. Moreover, we state the overall conservation error of the density computed by
which regard to machine precision and, thus, agree with the desired conservative nature of the DGSEM .
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.
As the convergence tests show, the experimental order of convergence depends on the number of vanishing moments in the Dirac-delta approximation, i.e. . 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 , whereas its initial conditions consist of a region inside a circle with the radius centered at the origin and a region outside that circle, see e.g. . The primitive variable initial conditions are then defined by two states, i.e.
Here, the inner state applies for every such that and the outside state applies otherwise. Further, the velocity vector is set to zero in the entire domain and we set .
For the simulation result illustrated in Fig. 2, one vanishing moment () and six continuous derivatives at the endpoints () of the Dirac-delta approximation are chosen. The support width in (3.7) is calculated with . Furthermore, the adaptive filtering with convex blending (3.25) is used with a minimum tolerance of and a maximum tolerance of . 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 .