Stable filtering procedures for nodal discontinuous Galerkin methods

We prove that the most common filtering procedure for nodal discontinuous Galerkin (DG) methods is stable. The proof exploits that the DG approximation is constructed from polynomial basis functions and that integrals are approximated with high-order accurate Legendre-Gauss-Lobatto quadrature. The theoretical discussion serves to re-contextualize stable filtering results for finite difference methods into the DG setting. It is shown that the stability of the filtering is equivalent to a particular contractivity condition borrowed from the analysis of so-called transmission problems. As such, the temporal stability proof relies on the fact that the underlying spatial discretization of the problem possesses a semi-discrete bound on the solution. Numerical tests are provided to verify and validate the underlying theoretical results.

itself, and is often done in an ad hoc fashion (filtering is applied "as little as possible, but as much as needed").
There exist a multitude of numerical methods to approximate the solution of hyperbolic time-dependent PDEs. The family of discontinuous Galerkin (DG) spectral methods is well-suited for hyperbolic problems due to its high-order nature and ability to propagate waves accurately over long times [1]. In particular, nodal collocation DG methods are attractive because of their computational efficiency [17]. Additionally, if a nodal DG method is constructed on the Legendre-Gauss-Lobatto (LGL) nodes then the discrete differentiation and integration operators satisfy a summation-by-parts (SBP) property [19,20,25,26] for any polynomial order [8]. This allows for semi-discrete stability estimates to be constructed for such high-order nodal DG methods, e.g [4,18].
Recently, Lundquist and Nordström [21] removed the ad hoc nature of filtering in the context of finite difference methods. Therein, they discuss a contractivity condition on the explicit filter matrix by re-framing the filtering procedure as a transmission problem [24]. Further, the work in [21] develops a necessary condition on the existence of an auxiliary filter matrix and how it must be related to the particular discrete integration (quadrature). Also, an implicit implementation of filtering was proved to be stable, i.e. sufficiency was obtained. The explicit implementation was more easily implemented, slightly more accurate and numerically shown to be stable but a proof was not obtained.
The goal of the present work is to re-interpret the theoretical work from [21] into the nodal DG context. In doing so, we will remove the ad hoc nature of the nodal DG filtering procedure and prove that the commonly used explicit filter technique from the spectral community [7,12,28] is stable in time. Just as in the case of finite difference methods, the temporal stability of the filtering for nodal DG relies upon the existence of a high-accuracy auxiliary filter matrix as well as a semi-discrete bound on the solution.
The remainder of the paper is organized as follows: Section 2 provides an overview of the nodal DG method and the commonly used filtering procedure. Then, Section 3 generalizes the theoretical time stability results from [21] into the DG context. Numerical results that support and verify the theoretical findings are given in Section 4. Our concluding remarks and outlook are given in the final section.

Overview of nodal DG approximations
Discontinuous Galerkin methods are principally designed to approximate solutions of hyperbolic conservation laws [13,17]. Here, we consider such a conservation law in one spatial dimension ∂u ∂t is the solution and f ≡ f (u) is the flux function. The conservation law is then equipped with an initial condition u(x, 0) ≡ u ini (x) and suitable boundary condition(s where the interpolation nodes are taken to be the N + 1 Legendre-Gauss-Lobatto (LGL) nodes.
6. Approximate integrals with LGL quadrature such that the DG scheme is collocated. 7. Arrive at the semi-discrete approximation that can be integrated in time.
The resulting strong form, nodal DG approximation is then is the vector form for the degrees of freedom. The matrices in the semi-discrete form (2.6) are the discrete derivative matrix [17] D ij = j (ξ i ), i, j = 0, . . . , N, where λmax is the maximum wave speed of the flux Jacobian. Two features to note for this flavour of nodal DG approximation are: The diagonal mass matrix M denotes a quadrature rule that is exact for polynomials up to degree 2N − 1. So, there is equality between the continuous and discrete integral provided the product of the functions v and w are polynomials of degree ≤ 2N − 1. The mass and derivative matrices of the LGL collocation DG scheme (2.6) form a summation-by-parts (SBP) operator pair [19,20,25,26] for any nodal polynomial order (2.13) From the SBP property (2.13), stable versions of the nodal DG method can be constructed, e.g. [4,5,8,10].

Construction of filtering for nodal DG
In the context of DG methods, the general idea of filtering exploits that the polynomial representation of the function U is unique, and hence can be written in terms of other basis polynomial functions. Basically, the filtering procedure is: 1. Transform the coefficients of the nodal approximation into a modal set of basis functions, e.g., the orthogonal Legendre polynomials (2.14) 2. Because the modal basis is hierarchical, it is straightforward for one to perform a cutoff in modal space to filter higher order modes. 3. Transform the filtered modal solution coefficients back to the nodal Lagrange basis.
At present, we select the modal (normalized) Legendre basis polynomials {L j (ξ)} N j=0 to construct the filtering. It is straightforward to compute the Vandermonde matrix V associated with the LGL nodal interpolation nodes which allows us to transform the nodal degrees of freedom, The filtering matrix is then constructed as [28] where C is a diagonal modal cutoff matrix. The conditions the filter function σ(η), as defined by Vandeven [28], are In the nodal DG community a typical choice is an exponential filter function, e.g. [6,9,13], to define the coefficients where α, s and Nc are the filter parameters. The value Nc indicates the number of the unaffected modes, α is chosen such that exp(−α) is machine epsilon and s is an even number determining the order (sometimes referred to as the strength) of the filter. Two common choices for the filtering parameters are to take α = 36, Nc = 4 and the filter strength is either "strong" with s = 16 or "weak" with s = 32 [6,9,13].
For any choice of the filter parameters the filter coefficients are constructed such that 0 ≤ σ i ≤ 1. Note, this exponential filter does not strictly adhere to Vandeven's definition of the filter function (2.19), but it does so in practice by choosing α such that σ(1) is below machine accuracy [2]. In summary, the filter matrix for the nodal DG approximation is given as (2.20) The filter of the form (2.20) retains the high-order accuracy of the nodal DG approximation for smooth functions [12,28] as shown numerically in Section 4.1.

Remark 1
Filtering have often been used as a stabilization technique for numerical methods such as in finite difference [16,30] as well as discontinuous Galerkin [6,9,13] methods. We strongly advise against such use. Instead, one should first construct construct a (provably) stable numerical scheme. After this, the solution quality can be addressed and cleaned-up, possibly using filtering.

Stability
As previously mentioned, it is possible to develop semi-discrete stability estimates for the nodal DG approximation via the SBP property (2.13). The filtering is a separate procedure which changes the approximate solution during the time integration procedure, for example after each explicit time step or even after each stage of a Runge-Kutta method. Here, we explore what influence this filtering step has on the stability estimate for the nodal DG approximation.
To discuss the filtering procedure and its affect on stability we re-interpret the work on provably stable filtering from [21] into the nodal DG context. In a broad sense, with homogeneous boundary conditions, semi-discrete stability ensures that the discrete norm of the approximate solution is bounded by the discrete norm of the initial conditions, see [23] for complete details. For the nodal DG approximation such a stability statement takes the form where U ini is the initial condition evaluated at the LGL nodes.
Pursuant to the work [21], we view the application of a filter matrix to a discrete solution at some intermediate time t 1 as a transmission problem [24]: where the operator D contains the derivative matrix D as well as the boundary conditions. For the present discussion, the filtering stated in the final line of (3.2) is performed in an explicit fashion. For stability it must hold that the filter is contractive, i.e.
In turn, this contractive property guarantees that the filter procedure is stable because The contractivity property in the discrete norm (3.3) then implies that the following contractivity condition on the explicit filter matrix F must hold  This contractive matrix property was first identified in [24] for stable transmission problems. The contractivity condition (3.5) expresses a precise interplay between the filter matrix and the mass matrix. As demonstrated in [21], a necessary condition for the explicit filter matrix to satisfy (3.5) is that an auxiliary filter matrix (3.6) exists and possesses the same accuracy as F. The accuracy requirement on Faux is necessary since otherwise (3.5) is provably indefinite [21].
We will show that the auxiliary filter matrix (3.6) is identical to the original filter matrix (2.20) for the LGL nodal DG approximation. Furthermore, the filter matrix F indeed satisfies the contractivity condition (3.5). Both results require the following Lemma.  Proof. The entries of this matrix product in terms of discrete inner products is From the accuracy of the LGL quadrature and the fact that {L j (ξ)} N j=0 are polynomials, we have equality between the discrete and continuous inner products provided j + k ≤ 2N − 1. The result above utilizes that the Legendre basis is orthonormal. The quadratures in (3.8) are therefore exact for all inner products except the one related to L N with itself since 2N > 2N − 1. Thus, (3.10) The discrete and continuous norms are equivalent [2]. Provided that φ is a polynomial of degree N , the discrete and continuous L 2 norms are related by ||φ|| N = 2 + 1/N ||φ|| for the LGL quadrature [3] . Using this fact, (3.10) becomes the diagonal matrix K.
We can now prove

Proposition 1
The auxiliary filter is identical to the DG filter matrix, i.e. Faux = F.
Proof. We examine the difference between the two filter matrices (3.11) Next, we factor out the matrix V on the left and V −1 on the right to have (3.12) Applying the result from Lemma 1 gives where we use that the matrices C and K are diagonal to obtain the desired result.

Remark 2
The accuracy of the filter F lies entirely in the filter function σ used to create the diagonal entries in the matrix C, as shown by Vandeven [28].
The following result is then self-evident.

Corollary 1 The auxiliary filter Faux exists and is as accurate as F.
Further, the result from Lemma 1 allows us to prove Proposition 2 The nodal DG filter matrix F is contractive in the sense of (3.5).
Proof. We substitute the form of the filter matrix (2.20) into the contractivity condition (3.5) to obtain (3.14) The middle term, V T M V, grouped above is precisely that from Lemma 1, which gives Next, we recall that the modal cutoff matrix is diagonal with entries {σ i } N i=0 and, by construction, the term σ N = 0, see (2.18). Thus, 0). (3.16) We combine this fact with the diagonal quadrature matrix result (3.7) to simplify the middle term of (3.15) to be C K C = C 2 . The contractivity condition then becomes where we, again, apply the result from Lemma 1. From (3.7) and (3.16) we see that

Remark 3 The proof above holds provided the filter function is chosen such that σ(η) ∈
[0, 1] and that it "clips" the highest mode, i.e. σ(N ) = 0, then the nodal DG filter matrix F satisfies the contractivity condition (3.5). Thus, other proposed filter functions like those found in Hussaini et al. [14] also produce a provably stable filter matrix.

Numerical results
Here we apply the nodal DG filter matrix in F (2.20) to several test problems. For these tests we select the parameters for a "strong" version of the nodal DG filter F from Section 2.2. We do so to demonstrate, in practice, the high-order accuracy of the filtered DG approximation for a smooth solution and how it performs for test cases that develops non-smooth solutions over time. To integrate the semi-discrete DG approximation (2.6) in time, we use a third-order Runge-Kutta from Williamson [29].

High-order convergence for linear advection
For the convergence test we consider the linear advection equation with flux function f (u) = au and take the wave speed to be the constant a = 1 on the domain Ω = [0, 1].
We use a smooth Gaussian pulse to set the initial and boundary conditions as well as compute errors (4.1) We vary the polynomial degree N at values in the interval [7,64] and integrate up T = 0.5 as the final time. In Fig. 1 we present a semilog plot of the L∞ error versus the polynomial order. We see that the error decays exponentially until the errors in the approximation are dominated by the time integration. Thereafter, we halve the time step size and see that the stagnation point of the error drops by a factor of eight, as expected for a third-order time integration technique. This experimentally demonstrates that the nodal DG approximation filtered in every time step remains high-order accurate in space and time for a smooth solution.

Variable wave speed for linear advection
Next, we consider a more complicated test proposed in [12]. For this case, the solution remains bounded, but develops steep gradients. Due to the high-degree polynomial approximation of the DG method, spurious oscillations can develop near these steep gradients and propagate throughout the domain, polluting the solution quality.
We  This wave speed remains positive at the boundaries of the domain, but it can change sign within the domain. We take the initial condition to be u ini (x) = sin(πx) which gives the corresponding analytical solution [11] u(x, t) = sin 2 tan −1 e −t tan πx − 1 2 + 1 . The solution in (4.3) develops a steep gradient around the point x = (1−π)/π ≈ −0.68 before it finally decays to a constant value lim t→∞ u(x, t) = sin(1). We compare the analytical and approximate solutions at a final time T = 4 on a single spectral element using polynomial order N = 256 and ∆t = 1/2000 as the explicit time step. The polynomial order and general setup are chosen such that a comparison can be made to the results in [12]. In Figure 2 we present the unfiltered approximation on the left and the filtered approximation on the right. Clearly, the unfiltered DG scheme contains spurious oscillations whereas the filtered solution suppresses such behaviour. Also, qualitatively, the results of the new DG filter scheme are very similar to those from [12].

Stability demonstration for Burgers'
This example is designed to illustrate the importance of the stability of the underlying spatial discretization and how the filtering can influence the behavior of the solution in time. For this we consider Burgers' equation and two forms of the nodal DG spatial discretization. One discretizes the conservative form of the PDE where the nonlinear Burgers' flux is simply f (u) = u 2 /2 while the other skew-symmetric discretization writes the spatial derivative of the flux in a split formulation

(4.4)
On the continuous level these two forms of Burgers' equation are equivalent; however, on the discrete level they exhibit different behavior. Most notably, the solution energy u 2 /2 is bounded for the nodal DG discretization constructed from the skew-symmetric formulation whereas no such bound exists for the discretization of the conservative form, see [8] for complete details.
Therefore, the discretization of the skew-symmetric form is provably stable and the conservative form is not. As discussed in the previous Sections, the DG filtering is a procedure divorced from the spatial discretization. If the underlying numerical scheme is provably stable, and the filtering is contractive, it will remove energy from the solution in a stable way. However, if the underlying scheme is unstable the filtering still removes energy and the approximation might be stable, but no further conclusions regarding the solution energy can be drawn.
To illustrate this we consider the domain Ω = [0, 2] with periodic boundary conditions and the initial condition (4.5) We run the simulation with polynomial order N = 128 up to T = 2.25 as a final time.
Further, the solution is filtered at 16 equally spaced times during the simulation. Due to the nonlinear nature of Burgers' equation the initial conditions will steepen and eventually a shock will form. We run four variants of the nodal DG scheme relevant to the present discussion: • Conservative formulation; Unfiltered.
On the left in Figure 3 we present the evolution of the solution energy, normalized with its initial value, over time. On the right in Figure 3 we give the approximate solution at the final time produced by the filtered skew-symmetric DG formulation as well as a reference solution created with a standard finite volume scheme on 10000 grid cells.  Due to the high polynomial order, the simulation is well resolved and the four variants are nearly indistinguishable for most of the simulation time. However, as the gradients steepen we see that the conservative formulation, for which no energy stability statement exists, behave erratically. The unfiltered conservative form simulation crashes at T ≈ 1.8. The filtered conservative form simulation successfully runs as the filtering keeps the solution energy "under control." This is illustrated in Figure 3 where we observe growth in the solution energy between the filter applications because the underlying spatial discretization is unstable. The solution energy of unfiltered and filtered skew-symmetric simulations both remain bounded because the underlying spatial discretization possesses an energy bound. Note, there is a small amount of dissipation in the solution energy for the unfiltered skew-symmetric scheme due to the formation of the shock [15]. Further, the filtered skew-symmetric simulation is less energetic, as expected, because the act of filtering removes some solution energy.

Closing remarks
We proved that the commonly used nodal DG filter matrix F satisfied a contractivity condition. Further, we proved that a high-order auxiliary filter matrix Faux exists for the nodal DG approximation which is necessary for the contractivity condition to be satisfied. Together, these results implied that the explicit filtering procedure in the context of nodal DG methods "removes" information in a stable way when measured in the norm induced by the Legendre-Gauss-Lobatto quadrature.
Numerical results were provided to demonstrate and verify that the filtering retained the high-order accuracy of the nodal DG approximation, that it suppressed spurious oscillations near steep gradients and was stable provided the underlying spatial discretization of the method had a semi-discrete bound.
The generalization of the results described in this work to multiple spatial dimensions is straightforward due to the tensor product nature of the nodal DG method. The same is true for problems involving curved physical boundaries, similar to those discussed in [21]. An interesting open question for future work is the extension of provably stable filtering to problems involving multiple DG elements where the interface coupling will play a key role.

Funding
This work was supported by Vetenskapsrådet, Sweden (award number 2018-05084 VR)

Conflicts of interest
The authors declare that they have no conflicts of interest in the present work.

Availability of data and material
Not applicable.