Stable filtering procedures for nodal discontinuous Galerkin methods

07/14/2020 ∙ by Jan Nordström, et al. ∙ Linköping University 0

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.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

High frequency errors are always present in numerical simulations due to inaccuracies of the spatial derivatives at high-wave 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 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 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 Legendre-Gauss-Lobatto (LGL) nodes then the discrete differentiation and integration operators satisfy a summation-by-parts (SBP) property kreiss1974finite ; kreiss1977 ; strand1994 ; svard2014 for any polynomial order gassner_skew_burgers

. This allows for semi-discrete stability estimates to be constructed for such high-order 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 re-framing 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 re-interpret 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 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 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


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


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


and rewrite the conservation law in the computational coordinate:


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 :

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

  2. Integrate-by-parts once and resolve discontinuities at the physical boundaries with a numerical flux function .

  3. Integrate-by-parts again to obtain strong form DG.

  4. Approximate the solution and flux with nodal polynomials of degree written in the Lagrange basis, e.g.


    where the interpolation nodes are taken to be the

    Legendre-Gauss-Lobatto (LGL) nodes.

  5. Select the test function to be the Lagrange polynomials with .

  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 Kopriva:2009nx


the matrix containing the LGL quadrature weights and the boundary matrix given by


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 Lax-Friedrichs flux


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


provided the product of the functions and are polynomials of degree . From (2.11) the continuous and discrete norms are denoted


The mass and derivative matrices of the LGL collocation DG scheme (2.6) form a summation-by-parts (SBP) operator pair kreiss1974finite ; kreiss1977 ; strand1994 ; svard2014 for any nodal polynomial order gassner_skew_burgers , i.e.


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:

  1. Transform the coefficients of the nodal approximation into a modal set of basis functions, e.g., the orthogonal Legendre polynomials

  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 to construct the filtering. It is straightforward to compute the Vandermonde matrix associated with the LGL nodal interpolation nodes


which allows us to transform the nodal degrees of freedom, , to modal degrees of freedom, and vice versa


The filtering matrix is then constructed as vandeven1991


where is a diagonal modal cutoff matrix. The conditions the filter function , as defined by Vandeven vandeven1991 , are

  • must have continuous derivatives where


In the nodal DG community a typical choice is an exponential filter function, e.g. chaudhur2017 ; Gassner:2013qf ; Hestahven:1008th , to define the coefficients


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


The filter of the form (2.20) retains the high-order 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 cleaned-up, possibly using filtering.

3 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 lundquist2020stable 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 Nordstrom:2016jk for complete details. For the nodal DG approximation such a stability statement takes the form


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 :


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.


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 must hold


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


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


The entries of this matrix product in terms of discrete inner products is


From the accuracy of the LGL quadrature and the fact that are polynomials, we have equality between the discrete and continuous inner products


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,


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. .


We examine the difference between the two filter matrices


Next, we factor out the matrix on the left and on the right to have


Applying the result from Lemma 1 gives


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 self-evident.

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).


We substitute the form of the filter matrix (2.20) into the contractivity condition (3.5) to obtain


The middle term, , grouped above is precisely that from Lemma 1, which gives


Next, we recall that the modal cutoff matrix is diagonal with entries and, by construction, the term , see (2.18). Thus,


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


where we, again, apply the result from Lemma 1. From (3.7) and (3.16) we see that


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 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 williamson1980 .

4.1 High-order 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


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 third-order time integration technique.

Figure 1: (left) Exact and approximate solution with for the linear advection with smooth Gauss pulse at . (right) Space and time convergence.

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.

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 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 consider the linear advection problem with a variable wave speed on the domain written in the form


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


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 .

Figure 2: Unfiltered (left) and filtered (right) nodal DG solution for the variable wave speed problem (4.2) at with polynomial order on a single spectral element.

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 skew-symmetric discretization writes the spatial derivative of the flux in a split formulation


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 skew-symmetric 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 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 with periodic boundary conditions and the initial condition


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.

  • Skew-symmetric formulation; Unfiltered.

  • Skew-symmetric 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 skew-symmetric DG formulation as well as a reference solution created with a standard finite volume scheme on 10000 grid cells.

Figure 3: (left) Solution energy evolution of four nodal DG variants for Burgers’ equation with the initial condition (4.5). (right) Plot of the filtered skew-symmetric solution () and a reference finite volume solution at .

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 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 jameson2008_energy . Further, the filtered skew-symmetric 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 high-order 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 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

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 2018-05084 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 (


  • (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 Navier-Stokes 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 skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods. SIAM Journal on Scientific Computing 35(3), A1233–A1253 (2013)
  • (9) Gassner, G.J., Beck, A.D.: On the accuracy of high-order 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 summation-by-parts 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 finite-difference 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. I-Fourier 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.: Well-posed 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 summation-by-parts schemes for initial-boundary-value 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.: Low-storage Runge-Kutta schemes. Journal of Computational Physics 35, 48–56 (1980)
  • (30) Yee, H.C., Sandham, N.D., Djomehri, M.J.: Low-dissipative high-order shock-capturing methods using characteristic-based filters. Journal of Computational Physics 150(1), 199–238 (1999)