1 Introduction
High frequency errors are always present in numerical simulations due to inaccuracies of the spatial derivatives at highwave numbers. These errors can be exacerbated due to aliasing errors or steep gradients in the solution which can exist (or arise) for solutions of partial differential equations (PDEs) with predominantly hyperbolic character. To combat such errors, one technique is to use a filter operator that removes high wave number oscillations, e.g.
chaudhur2017 ; Gassner:2013qf ; Hestahven:1008th ; meister2012 . The filtering procedure is separate from the scheme 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 timedependent PDEs. The family of discontinuous Galerkin (DG) spectral methods is wellsuited for hyperbolic problems due to its highorder nature and ability to propagate waves accurately over long times Airnsworth2004 . In particular, nodal collocation DG methods are attractive because of their computational efficiency Kopriva:2009nx . Additionally, if a nodal DG method is constructed on the LegendreGaussLobatto (LGL) nodes then the discrete differentiation and integration operators satisfy a summationbyparts (SBP) property kreiss1974finite ; kreiss1977 ; strand1994 ; svard2014 for any polynomial order gassner_skew_burgers
. This allows for semidiscrete stability estimates to be constructed for such highorder nodal DG methods, e.g
carpenter_esdg ; Gassner:2013ij .Recently, Lundquist and Nordström lundquist2020stable 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 reframing the filtering procedure as a transmission problem nordstrom2018well . Further, the work in lundquist2020stable 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 reinterpret the theoretical work from lundquist2020stable 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 don1994numerical ; hesthaven2008filtering ; vandeven1991 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 highaccuracy auxiliary filter matrix as well as a semidiscrete 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 lundquist2020stable 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.
2 Overview of nodal DG approximations
Discontinuous Galerkin methods are principally designed to approximate solutions of hyperbolic conservation laws Hestahven:1008th ; Kopriva:2009nx . Here, we consider such a conservation law in one spatial dimension
(2.1) 
where is the solution and is the flux function. The conservation law is then equipped with an initial condition and suitable boundary condition(s). Two prototypical examples of conservation laws are the linear advection equation and Burgers’ equation whose corresponding flux functions are
(2.2) 
Next, we provide the broad strokes to arrive at the nodal DG approximation of the conservation law (2.1). First, we transform the domain to the reference interval . To do so, we apply the mapping
(2.3) 
and rewrite the conservation law in the computational coordinate:
(2.4) 
2.1 The numerical scheme
The DG approximation is built from the weak form of the mapped equation (2.4). We list the nodal DG approximation steps below with full details given in Kopriva:2009nx :

Multiply by a test function and integrate over the reference domain.

Integratebyparts once and resolve discontinuities at the physical boundaries with a numerical flux function .

Integratebyparts again to obtain strong form DG.

Approximate the solution and flux with nodal polynomials of degree written in the Lagrange basis, e.g.
(2.5) where the interpolation nodes are taken to be the
LegendreGaussLobatto (LGL) nodes. 
Select the test function to be the Lagrange polynomials with .

Approximate integrals with LGL quadrature such that the DG scheme is collocated.

Arrive at the semidiscrete approximation that can be integrated in time.
The resulting strong form, nodal DG approximation is then
(2.6) 
where
(2.7) 
is the vector form for the degrees of freedom. The matrices in the semidiscrete form (
2.6) are the discrete derivative matrix Kopriva:2009nx(2.8) 
the matrix containing the LGL quadrature weights and the boundary matrix given by
(2.9) 
There exist several numerical flux functions depending on the particular mathematical flux function , see Toro toro2009 for details. A general (and simple) numerical flux which we use in the present work is the local LaxFriedrichs flux
(2.10) 
where 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 denotes a quadrature rule that is exact for polynomials up to degree . So, there is equality between the continuous and discrete integral
(2.11) 
provided the product of the functions and are polynomials of degree . From (2.11) the continuous and discrete norms are denoted
(2.12) 
The mass and derivative matrices of the LGL collocation DG scheme (2.6) form a summationbyparts (SBP) operator pair kreiss1974finite ; kreiss1977 ; strand1994 ; svard2014 for any nodal polynomial order gassner_skew_burgers , i.e.
(2.13) 
From the SBP property (2.13), stable versions of the nodal DG method can be constructed, e.g. carpenter_esdg ; chan2018 ; gassner_skew_burgers ; Gassner:2016ye .
2.2 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 is unique, and hence can be written in terms of other basis polynomial functions. Basically, the filtering procedure is:

Transform the coefficients of the nodal approximation into a modal set of basis functions, e.g., the orthogonal Legendre polynomials
(2.14) 
Because the modal basis is hierarchical, it is straightforward for one to perform a cutoff in modal space to filter higher order modes.

Transform the filtered modal solution coefficients back to the nodal Lagrange basis.
At present, we select the modal (normalized) Legendre basis polynomials to construct the filtering. It is straightforward to compute the Vandermonde matrix associated with the LGL nodal interpolation nodes
(2.15) 
which allows us to transform the nodal degrees of freedom, , to modal degrees of freedom, and vice versa
(2.16) 
The filtering matrix is then constructed as vandeven1991
(2.17) 
where is a diagonal modal cutoff matrix. The conditions the filter function , as defined by Vandeven vandeven1991 , are


must have continuous derivatives where
(2.18)
In the nodal DG community a typical choice is an exponential filter function, e.g. chaudhur2017 ; Gassner:2013qf ; Hestahven:1008th , to define the coefficients
(2.19) 
where , and are the filter parameters. The value indicates the number of the unaffected modes, is chosen such that is machine epsilon and 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 , and the filter strength is either “strong” with or “weak” with chaudhur2017 ; Gassner:2013qf ; Hestahven:1008th .
For any choice of the filter parameters the filter coefficients are constructed such that . 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 is below machine accuracy canuto2006 .
In summary, the filter matrix for the nodal DG approximation is given as
(2.20) 
The filter of the form (2.20) retains the highorder accuracy of the nodal DG approximation for smooth functions hesthaven2008filtering ; vandeven1991 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 kennedy1997 ; yee1999 as well as discontinuous Galerkin chaudhur2017 ; Gassner:2013qf ; Hestahven:1008th 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 cleanedup, possibly using filtering.
3 Stability
As previously mentioned, it is possible to develop semidiscrete 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 RungeKutta 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 reinterpret the work on provably stable filtering from lundquist2020stable into the nodal DG context. In a broad sense, with homogeneous boundary conditions, semidiscrete stability ensures that the discrete norm of the approximate solution is bounded by the discrete norm of the initial conditions, see Nordstrom:2016jk for complete details. For the nodal DG approximation such a stability statement takes the form
(3.1) 
where is the initial condition evaluated at the LGL nodes.
Pursuant to the work lundquist2020stable , we view the application of a filter matrix to a discrete solution at some intermediate time as a transmission problem nordstrom2018well :
(3.2)  
where the operator contains the derivative matrix 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.
(3.3) 
In turn, this contractive property guarantees that the filter procedure is stable because
(3.4) 
The contractivity property in the discrete norm (3.3) then implies that the following contractivity condition on the explicit filter matrix must hold
(3.5) 
This contractive matrix property was first identified in nordstrom2018well for stable transmission problems. The contractivity condition (3.5) expresses a precise interplay between the filter matrix and the mass matrix. As demonstrated in lundquist2020stable , 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 . The accuracy requirement on is necessary since otherwise (3.5) is provably indefinite lundquist2020stable .
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 indeed satisfies the contractivity condition (3.5). Both results require the following Lemma.
Lemma 1.
The matrix product is the LGL quadrature rule applied to the (normalized) Legendre polynomial functions and results in a diagonal matrix
(3.7) 
Proof.
The entries of this matrix product in terms of discrete inner products is
(3.8) 
From the accuracy of the LGL quadrature and the fact that are polynomials, we have equality between the discrete and continuous inner products
(3.9) 
provided . 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 with itself since . Thus,
(3.10) 
The discrete and continuous norms are equivalent canuto2006 . Provided that is a polynomial of degree , the discrete and continuous norms are related by for the LGL quadrature ISI:A1982NE30900005 . Using this fact, (3.10) becomes the diagonal matrix . ∎
We can now prove
Proposition 1.
The auxiliary filter is identical to the DG filter matrix, i.e. .
Proof.
We examine the difference between the two filter matrices
(3.11) 
Next, we factor out the matrix on the left and on the right to have
(3.12)  
Applying the result from Lemma 1 gives
(3.13) 
where we use that the matrices and are diagonal to obtain the desired result. ∎
Remark 2
The accuracy of the filter lies entirely in the filter function used to create the diagonal entries in the matrix , as shown by Vandeven vandeven1991 .
The following result is then selfevident.
Corollary 1.
The auxiliary filter exists and is as accurate as .
Further, the result from Lemma 1 allows us to prove
Proposition 2.
The nodal DG filter matrix 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, , grouped above is precisely that from Lemma 1, which gives
(3.15) 
Next, we recall that the modal cutoff matrix is diagonal with entries and, by construction, the term , see (2.18). Thus,
(3.16) 
We combine this fact with the diagonal quadrature matrix result (3.7) to simplify the middle term of (3.15) to be . The contractivity condition then becomes
(3.17) 
where we, again, apply the result from Lemma 1. From (3.7) and (3.16) we see that
(3.18) 
because for . Thus, (3.5) holds. ∎
Remark 3
The proof above holds provided the filter function is chosen such that and that it “clips” the highest mode, i.e. , then the nodal DG filter matrix satisfies the contractivity condition (3.5). Thus, other proposed filter functions like those found in Hussaini et al. hussaini1985spectral also produce a provably stable filter matrix.
4 Numerical results
Here we apply the nodal DG filter matrix in (2.20) to several test problems. For these tests we select the parameters for a “strong” version of the nodal DG filter from Section 2.2. We do so to demonstrate, in practice, the highorder accuracy of the filtered DG approximation for a smooth solution and how it performs for test cases that develops nonsmooth solutions over time. To integrate the semidiscrete DG approximation (2.6) in time, we use a thirdorder RungeKutta from Williamson williamson1980 .
4.1 Highorder convergence for linear advection
For the convergence test we consider the linear advection equation with flux function and take the wave speed to be the constant on the domain .
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 at values in the interval and integrate up as the final time. In Fig. 1 we present a semilog plot of the 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 thirdorder time integration technique.
This experimentally demonstrates that the nodal DG approximation filtered in every time step remains highorder accurate in space and time for a smooth solution.
4.2 Variable wave speed for linear advection
Next, we consider a more complicated test proposed in hesthaven2008filtering . For this case, the solution remains bounded, but develops steep gradients. Due to the highdegree polynomial approximation of the DG method, spurious oscillations can develop near these steep gradients and propagate throughout the domain, polluting the solution quality.
We consider the linear advection problem with a variable wave speed on the domain written in the form
(4.2) 
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 which gives the corresponding analytical solution gottlieb1981stability
(4.3) 
The solution in (4.3) develops a steep gradient around the point before it finally decays to a constant value .
We compare the analytical and approximate solutions at a final time on a single spectral element using polynomial order and as the explicit time step. The polynomial order and general setup are chosen such that a comparison can be made to the results in hesthaven2008filtering . 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 hesthaven2008filtering .
4.3 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
while the other skewsymmetric 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 is bounded for the nodal DG discretization constructed from the skewsymmetric formulation whereas no such bound exists for the discretization of the conservative form, see gassner_skew_burgers for complete details.
Therefore, the discretization of the skewsymmetric 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 with periodic boundary conditions and the initial condition
(4.5) 
We run the simulation with polynomial order up to 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:

[label=]

Conservative formulation; Unfiltered.

Conservative formulation; Filtered.

Skewsymmetric formulation; Unfiltered.

Skewsymmetric formulation; Filtered.
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 skewsymmetric 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 . 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 skewsymmetric 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 skewsymmetric scheme due to the formation of the shock jameson2008_energy . Further, the filtered skewsymmetric simulation is less energetic, as expected, because the act of filtering removes some solution energy.
5 Closing remarks
We proved that the commonly used nodal DG filter matrix satisfied a contractivity condition. Further, we proved that a highorder auxiliary filter matrix 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 LegendreGaussLobatto quadrature.
Numerical results were provided to demonstrate and verify that the filtering retained the highorder 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 semidiscrete 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
lundquist2020stable . 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.6 Declarations
6.1 Funding
This work was supported by Vetenskapsrådet, Sweden (award number 201805084 VR)
6.2 Conflicts of interest
The authors declare that they have no conflicts of interest in the present work.
6.3 Availability of data and material
Not applicable.
6.4 Code availability
The code used to generate the results in this work is available upon request with Andrew Winters (andrew.ross.winters@liu.se).
References
 (1) Ainsworth, M.: Dispersive and dissipative behaviour of high order discontinuous Galerkin finite element methods. Journal of Computational Physics 198, 106–130 (2004)
 (2) Canuto, C., Hussaini, M., Quarteroni, A., Zang, T.: Spectral Methods: Fundamentals in Single Domains. Springer, Berlin (2006)
 (3) Canuto, C., Quarteroni, A.: Approximation results for orthogonal polynomials in Sobolev spaces. Mathematics of Computation 38(157), 67–86 (1982)
 (4) Carpenter, M., Fisher, T., Nielsen, E., Frankel, S.: Entropy stable spectral collocation schemes for the NavierStokes equations: Discontinuous interfaces. SIAM Journal on Scientific Computing 36(5), B835–B867 (2014)
 (5) Chan, J.: On discretely entropy conservative and entropy stable discontinuous Galerkin methods. Journal of Computational Physics 362, 346–374 (2018)
 (6) Chaudhuri, A., Jacobs, G.B., Don, W.S., Abbassi, H., Mashayek, F.: Explicit discontinuous spectral element method with entropy generation based artificial viscosity for shocked viscous flows. Journal of Computational Physics 332, 99–117 (2017)
 (7) Don, W.S.: Numerical study of pseudospectral methods in shock wave applications. Journal of Computational Physics 110(1), 103–111 (1994)
 (8) Gassner, G.: A skewsymmetric discontinuous Galerkin spectral element discretization and its relation to SBPSAT finite difference methods. SIAM Journal on Scientific Computing 35(3), A1233–A1253 (2013)
 (9) Gassner, G.J., Beck, A.D.: On the accuracy of highorder discretizations for underresolved turbulence simulations. Theoretical and Computational Fluid Dynamics 27(3–4), 221–237 (2013)
 (10) Gassner, G.J., Winters, A.R., Kopriva, D.A.: Split form nodal discontinuous Galerkin schemes with summationbyparts property for the compressible Euler equations. Journal of Computational Physics 327, 39–66 (2016)
 (11) Gottlieb, D., Orszag, S.A., Turkel, E.: Stability of pseudospectral and finitedifference methods for variable coefficient problems. mathematics of computation 37(156), 293–305 (1981)
 (12) Hesthaven, J.S., Kirby, R.M.: Filtering in Legendre spectral methods. Mathematics of Computation 77(263), 1425–1452 (2008)
 (13) Hesthaven, J.S., Warburton, T.: Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. Springer (2008)
 (14) Hussaini, M., Kopriva, D., Salas, M., Zang, T.: Spectral methods for the Euler equations. IFourier methods and shock capturing. AIAA Journal 23(1), 64–70 (1985)
 (15) Jameson, A.: The construction od discretely conservative finite volume schemes that also globally conserve energy or entropy. Journal of Scientific Computing 34(2), 152–187 (2008)
 (16) Kennedy, C.A., Carpenter, M.H.: Comparison of several numerical methods for simulation of compressible shear layers, vol. 3484. NASA, Langley Research Center (1997)
 (17) Kopriva, D.A.: Implementing Spectral Methods for Partial Differential Equations. Scientific Computation. Springer (2009)
 (18) Kopriva, D.A., Gassner, G.: An energy stable discontinuous Galerkin spectral element discretization for variable coefficient advection problems. SIAM Journal on Scientific Computing 36(4), A2076–A2099 (2014)
 (19) Kreiss, H.O., Scherer, G.: Finite element and finite difference methods for hyperbolic partial differential equations. In: Mathematical aspects of finite elements in partial differential equations, pp. 195–212. Elsevier (1974)
 (20) Kreiss, H.O., Scherer, G.: On the existence of energy estimates for difference approximations for hyperbolic systems. Tech. rep., Deptpartment of Scientific Computing, Uppsala University (1977)
 (21) Lundquist, T., Nordström, J.: Stable and accurate filtering procedures. Journal of Scientific Computing 82(1), 1–21 (2020)
 (22) Meister, A., Ortleb, S., Sonar, T.: Application of spectral filtering to discontinuous Galerkin methods on triangulations. Numerical Methods for Partial Differential Equations 28(6), 1840–1868 (2012)
 (23) Nordström, J.: A roadmap to well posed and stable problems in computational physics. Journal of Scientific Computing 71(1), 365–385 (2016)
 (24) Nordström, J., Linders, V.: Wellposed and stable transmission problems. Journal of Computational Physics 364, 95–110 (2018)
 (25) Strand, B.: Summation by parts for finite difference approximations for . Journal of Computational Physics 110 (1994)
 (26) Svärd, M., Nordström, J.: Review of summationbyparts schemes for initialboundaryvalue problems. Journal of Computational Physics 268, 17–38 (2014)
 (27) Toro, E.F.: Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction. Springer (2009)
 (28) Vandeven, H.: Family of spectral filters for discontinuous problems. Journal of Scientific Computing 6(2), 159–192 (1991)
 (29) Williamson, J.H.: Lowstorage RungeKutta schemes. Journal of Computational Physics 35, 48–56 (1980)
 (30) Yee, H.C., Sandham, N.D., Djomehri, M.J.: Lowdissipative highorder shockcapturing methods using characteristicbased filters. Journal of Computational Physics 150(1), 199–238 (1999)
Comments
There are no comments yet.