A comparative Fourier analysis of discontinuous Galerkin schemes for advection–diffusion with respect to BR1, BR2, and local discontinuous Galerkin diffusion discretization

This work compares the wave propagation properties of discontinuous Galerkin (DG) schemes for advection–diffusion problems with respect to the behavior of classical discretizations of the diffusion terms, that is, two versions of the local discontinuous Galerkin (LDG) scheme as well as the BR1 and the BR2 scheme. The analysis highlights a significant difference between the two possible ways to choose the alternating LDG fluxes showing that the variant that is inconsistent with the upwind advective flux is more accurate in case of advection–diffusion discretizations. Furthermore, whereas for the BR1 scheme used within a third order DG scheme on Gauss‐Legendre nodes, a higher accuracy for well‐resolved problems has previously been observed in the literature, this work shows that higher accuracy of the BR1 discretization only holds for odd orders of the DG scheme. In addition, this higher accuracy is generally lost on Gauss–Legendre–Lobatto nodes.


INTRODUCTION
In this work, the wave propagation properties of discontinuous Galerkin (DG) schemes 1 for advection-diffusion problems are compared with respect to the behavior of several classical discretizations of the diffusion terms, that is, two versions of the local discontinuous Galerkin (LDG) scheme 2 as well as the first and second method of Bassi and Rebay, 3,4 termed BR1 and BR2, respectively. DG schemes are based on a variational formulation similar to classical finite element schemes but allow for discontinuous approximate solutions. Here, we will consider the simplified case that the approximate solution restricted to an arbitrary element is a polynomial of fixed degree N. Due to its flexibility and generality, the DG scheme is a popular numerical method in a variety of applications ranging from compressible fluid flow and aeroacoustics in [5][6][7] to electromagnetics in, 8,9 meteorology in, 10,11 and geophysics in. 12 The main advantages of the DG approach are its local conservation property, arbitrarily high order of accuracy, and superconvergence capabilities combined with compact stencils that facilitate hp-adaptivity as well as computation in parallel hardware environment. ORTLEB The investigation of wave propagation properties in terms of dispersion and diffusion errors depending on the wave number is of utmost importance for the analysis of accuracy and stability of any numerical scheme applied in the context of computational fluid dynamics. Especially in case of high order methods and for under-resolved turbulence simulation, a desired small numerical dissipation competes with robustness and thus has to be carefully analyzed. Therefore, dispersion and diffusion properties have been investigated for major classes of high order schemes such as the DG scheme in, 13 the spectral difference method in, 14 flux reconstruction schemes in, 15 and continuous Galerkin (CG) approximations in. 16 Dissipation and dispersion properties are usually inspected via Fourier analysis. For linear advection, Hu et al 13 have shown that the DG scheme admits one physical mode and N spurious modes that dissipate quickly for upwind fluxes but remain for central fluxes. Gassner and Kopriva 17 investigated the influence of Gauss-Legendre and Gauss-Legendre-Lobatto integration rules on the dispersion and dissipation characteristics of nodal DG schemes. Moura et al 18 observed that the spurious modes are in fact replications of the physical mode and contribute to the overall accuracy of the scheme. In particular, for higher wave numbers, the secondary eigenmodes may strongly influence the behavior of the scheme. Furthermore, based on the related flux reconstruction approach, Asthana and Jameson 19 derived a family of schemes with minimal dissipation and dispersion error for advection problems, and Vermeire and Vincent 20 investigated the behavior of fully discrete flux reconstruction schemes, which includes the influence of the chosen Runge-Kutta scheme for time integration.
In the above works, with the exception of Moura et al, 16 eigenanalysis of the numerical scheme is based on the linear advection equation. Furthermore, Manzanero et al 21 develop a generalized von Neumann technique to study the dispersion and dissipation properties of various DG schemes for nonconstant coefficient advection equations. However, parabolic equations have not been extensively studied in this context. In fact, the DG approach was originally introduced for the numerical approximation of first-order hyperbolic conservation laws as an alternative to widely used finite volume methods. Therefore, although quite natural for advection operators, the inherent discontinuity of the DG approximate solution at first does not offer any intrinsic way to discretize diffusion operators. When computing the diffusion flux at cell interfaces, neither a unique solution value nor unique derivatives are available. However, several techniques to compute diffusion fluxes within a DG framework have been successfully employed, either using specially designed penalty terms within a finite element approach as in 22,23 or by rewriting equations of advection-diffusion type into a system of first order equations using auxiliary variables for the solution derivatives as in. [2][3][4]24 Among these schemes, the first method by Bassi and Rebay 3 is the first extension of the DG scheme to the compressible Navier-Stokes equations and is usually termed the BR1 scheme. Hereby, after the reformulation of the viscous terms into a larger, extended first-order degenerate system of partial differential equations (PDEs) with the gradient as a new unknown, the standard DG approach is applied to the extended system, which necessitates to prescribe two types of numerical fluxes. The BR1 scheme uses the simplest approach, that is, arithmetic means for both types of fluxes. Motivated by the successful numerical results obtained with the BR1 scheme, Cockburn and Shu 2 analyzed various methods based on the reformulation into a first-order PDE and derived conditions on the numerical fluxes to guarantee stability, convergence, and a suboptimal error estimate of order N when using an approximation space of polynomial degree N. Their analysis shows suboptimal convergence of the BR1 scheme for odd N, whereas the choice of alternating numerical fluxes usually associated with the LDG scheme by Cockburn and Shu 2 lead to optimal convergence of order N + 1. Bassi and Rebay introduced a second approach in, 25 termed the BR2 scheme, which modifies the BR1 approach to yield a more compact stencil suitable for efficient implicit time integration. The BR2 diffusion discretization includes a penalty parameter that determines stability and accuracy of the scheme. Still, the BR1 scheme is considered an attractive approach because it is simple to code, parameter-free, and generic for nonlinear viscous fluxes and arbitrary grids. Therefore, it has recently been reinvestigated in Gassner et al, 26 where the neutral behavior of BR1 with respect to artificial dissipation over element interfaces and the resulting stability for the compressible Navier-Stokes equations has been proven. Considering the BR2 flux for the DG scheme on Gauss-Legendre-Lobatto nodes, the lifting operator may be calculated either by exact projection or using inexact numerical integration. Both versions are brought together in this work by extending recent results by Quaegebeur et al 27 on the equivalence of the BR2 scheme and the classical interior penalty formulation for linear diffusion in one space dimension.
Regarding wave propagation characteristics, the investigation of dissipation and dispersion properties of DG schemes applied to advection-diffusion problems is more recent than for pure advection. For pure diffusion problems, eigenanalysis for both well-resolved and under-resolved cases has been carried out, for example, by Huynh 28 for flux reconstruction schemes and by Alhawwary and Wang 29 for DG schemes. Furthermore, analytic Fourier eigenanalysis of DG diffusion discretizations for the well-resolved regime of wave numbers has been carried out for the DG scheme in Guo et al. 30 and Zhang and Shu. 31 For advection-diffusion problems, an eigenanalysis by Manzanero et al 32 for DG schemes considers the influence of a parameter-dependent Riemann solver for advective terms and the BR1 scheme for viscous terms. Their study investigates both the individual contribution of the dissipative mechanisms on the whole range of wave numbers and their combined effect. In addition, the authors correlate their findings to the results for 3D compressible Navier-Stokes flow. Furthermore, Watkins et al 33 carry out a von Neumann analysis of nodal DG schemes obtained via flux reconstruction, in order to investigate the stability, dissipation, and dispersion properties for advection-diffusion problems. In particular, their work analyzes the influence of different interface flux formulations. More precisely, for a DG scheme of polynomial degree N = 2 on Gauss-Legendre nodes, differences of its wave propagation properties are studied in case of either upwind or central flux for advection as well as either a particular alternate LDG scheme or the BR1 approach for diffusion. Their results show that the corresponding schemes with central flux discretization (central flux for advection and BR1 discretization for diffusion) produce smaller errors for well-resolved solutions whereas one-sided flux discretizations (upwind flux combined with LDG discretization) produce smaller errors in the under-resolved case. It is also worth mentioning that the analysis by Watkins et al 33 combines all eigenmodes into a wave number dependent error measure for the DG approximate solution instead of considering only the dispersion and dissipation properties of a single eigenmode which is regarded as the physical one.
Following the investigations by Zhang and Shu 31 and Guo et al, 30 a comparative Fourier analysis of the Gauss-Legendre and Gauss-Lobatto DG schemes has been carried out in Ortleb 34 for N = 1 using different viscous flux discretizations. This analysis confirmed the higher accuracy of Gauss-Legendre integration nodes in the well-resolved regime, that is, for small wave numbers or small cell sizes. In this work, the analysis carried out in Ortleb 34 is extended to the full range of wave numbers in order to gain additional insight into the numerical behavior for under-resolved waves. Here, we use the combined approach of Watkins et al to obtain a wave number dependent error bound. Extending the study by Watkins et al, 33 DG schemes on Gauss-Legendre-Lobatto nodes are considered as well and a larger range of interface fluxes for diffusion terms is investigated including the BR2 scheme. Whereas the focus of the analysis in Watkins et al 33 was put on the DG scheme for N = 2, in this work, the combined error measure is computed for a larger range of polynomial degrees to detect odd-even phenomena.
Regarding the wave number dependent error analysis, significant differences can be noted for the classical DG diffusion discretizations LDG and BR1/BR2. The original LDG scheme 2 represents a parameter-dependent family of diffusion discretizations. However, it commonly denotes the choice of alternating fluxes for the reformulated viscous term. For this choice, there are two alternatives, also considered with respect to their superconvergence properties by Cheng and Shu. 35 For advection-diffusion equations, this means that the diffusion fluxes are chosen either consistent with the convective flux, namely, the numerical flux for the unknown quantity is chosen equal to the upwind flux, whereas the numerical flux for the gradient is taken from the opposite direction, or the choice of diffusion fluxes is inconsistent with the convective flux when the directions are reversed. Whereas for pure diffusion problems, there is no preferred direction, the LDG scheme for advection-diffusion problems is often preferred in its consistent variant. Although Cheng and Shu 35 show superconvergence properties for both variants, they state a preference for the consistent variant in their conclusions. Furthermore, a recent investigation on IMEX time integration for LDG schemes by Wang et al 36 also uses the consistent variant. In this work, it is shown that the two choices of alternating fluxes lead to significant differences in the error versus wave number characteristics. More precisely, the variant that is inconsistent with the upwind advective flux is more accurate in case of advection-diffusion discretizations for a large range of wave numbers. Furthermore, the investigations in this work show that the higher accuracy of the BR1 flux for well-resolved solutions observed in Watkins et al 33 is restricted to the DG scheme on Gauss-Legendre nodes and to an even polynomial degree N. For the DG scheme on Gauss-Legendre nodes with odd N, both the BR2 flux for larger penalty parameters and the more accurate variant of LDG produce smaller errors for well-resolved problems, whereas for Gauss-Legendre-Lobatto nodes, the inconsistent LDG variant generally performs better as well.
This paper is organized as follows. Section 2 introduces the DG space discretization for linear advection-diffusion equations in one space dimension based on a reformulation of this equation as a system of first-order PDEs. For the BR2 diffusion flux introduced in Section 2, results on energy stability depending on the BR2 penalty parameter are reviewed and extended in Section 3. In Section 4, Fourier analysis is introduced to numerically compute the eigensolutions of the respective variants of the DG scheme for the advection-diffusion equation and derive the corresponding relative error versus nondimensional wave number. In Section 5, numerical experiments are carried out to verify the error analysis. Finally, a conclusion is given in Section 6.

THE DG SCHEME FOR ADVECTION-DIFFUSION EQUATIONS
In this work, we study the wave propagation properties of the DG-discretized linear advection-diffusion equation. Hence, we consider the problem with diffusion coefficient d > 0 and advective velocity a > 0, supplemented by the periodic initial condition U(x, 0) = U 0 (x) in L 2 ( ) and periodic boundary conditions. For space discretization, is partitioned into cells I = where  N (I ) denotes the space of polynomial functions on I j of degree at most N. As usual in DG schemes, the functions in V h may be discontinuous across element boundaries. At each element boundary, the left-hand side and right-hand side values of a piecewise continuous function v are denoted by v − and v + , respectively. The corresponding jump at element interfaces is denoted by [v] = v + − v − , and the arithmetic mean is given by {v} = 1 2 (v + + v − ). The BR1, BR2, and the LDG scheme are derived from the first-order reformulation of the advection-diffusion Equation (1) with an auxiliary variable Q. Denoting by (·, ·) j the usual inner product in L 2 (I j ), the element-wise DG space discretization to obtain the approximate solution u(t) ∈ V h as well as the auxiliary variable q(t) ∈ V h is derived from the variational formulation (q, r) = −(u, r x ) +û + 1 whereq andû represent suitable numerical fluxes determining the chosen DG diffusion scheme and upwind fluxes are applied to the advective term. Commonly, the semidiscretization (4) and (5) on a specific subinterval I j is transformed to the reference cell I = [−1, 1]. For this purpose, we use the transformation j defined by Λ ( ) = 1+ 2 x + 1 2 for ∈ I. Furthermore, in this work, numerical integration of the occurring integrals is carried out either exactly using Gauss-Legendre nodes or less accurately with Gauss-Legendre-Lobatto quadrature. Denoting by , = 1, … , N + 1, the set of quadrature nodes on I and by , = 1, … , N + 1, the corresponding quadrature weights, we hence replace the integrals occurring in Equations (4) and (5) by The formulation (4),(5) is thus replaced by < q, r> = − < u, r x > +û + 1 In the Equations 7 and (8), the simple choiceq BR1 represents the BR1 approach. Furthermore, the original LDG scheme 2 yields a parameter-dependent family of diffusion fluxesq LDG which, in one space dimension, contains the BR1 approach as a specific case with c 11 = c 12 = 0. Here, we use the common choice of alternating LDG fluxes with c 12 = ±1, c 11 = 0, which offers two different variants. One implementation is thus given byq LDG a which uses opposite wind direction compared with the upwind advective flux in Equation (7) and is therefore termed as inconsistent with the advective flux by Cheng and Shu. 35 The second variant is specified byq and termed consistent with the advective flux. The BR2 scheme is a modification of the BR1 approach, which was first suggested in Bassi et al 25 in order to obtain a more compact stencil suitable for efficient implicit time integration since the BR2 stencil contains only immediate neighbors in contrast to BR1. For the BR2 approach, the numerical flux in the auxiliary Equation (8) equals the choice for BR1, that is,û BR2 = {u}. Furthermore,q BR2 is determined by local lifting operators l + 1 2 which lift the jumps [u] + 1 2 into the DG approximation space. More precisely, we calculate l + 1 2 ([u]) ∈ V h based on the projection property where (·, ·) denotes the classical . The numerical fluxq in Equation (7) is then defined byq with a parameter e that was set to e = 1 in the original formulation by Bassi and Rebay 25 but is considered variable by Brezzi et al. 37 Replacing the exact integration in Equation (9) by the numerical quadrature rule used within the DG scheme does not change the lifting operator in case of Gauss-Legendre nodes. However, if the DG scheme is based on Gauss-Legendre-Lobatto integration, we obtain a second variant of the BR2 lifting operator due to inexact numerical integration. We denote by BR2 GLL the variant based on Gauss-Legendre-Lobatto nodes, whereas the lifting operator based on exact integration calculated, for example, by Gauss-Legendre quadrature is denoted by BR2 GL . For the BR2 GLL flux, the lifting operator is therefore computed as where < ·, · > GLL denotes the discrete inner product on V h computed via Gauss-Legendre-Lobatto quadrature, that is, Considering the penalty parameter e , it has been shown by Brezzi et al 37 using a coercivity condition, that the BR2 scheme is stable on triangular grids if e > 3. Because e is determined by the number of adjacent cells, this corresponds to e > 2 for the one-dimensional case. Quaegebeur et al 27 obtained a sharper bound on e via energy stability considerations, which yields e ≥ N N+1 for BR2 GL as provided in Theorem 2. However, the variant BR2 GLL is not considered in their analysis. Therefore, the results by Quaegebeur et al are extended in Section 3 to the calculation of the lifting operator via Gauss-Legendre-Lobatto quadrature, see Theorem 3. In summary, it is shown that the BR2 GLL scheme with parameter e is equivalent to the BR2 GL scheme with parameter̂e = N N+1 e . In addition, for the DG scheme on Gauss-Legendre-Lobatto nodes, the BR1 scheme is proven equivalent to BR2 GLL for e = 1.
Applying partial integration in space to the terms < u, v x > j , < q, v x > j in Equation (7) and < u, r x > j in Equation (8), which are exactly integrated also in case of Gauss-Legendre-Lobatto quadrature, we obtain the cell-wise strong form of the DG semidiscretization on the reference cell as follows, see, for example, Hesthaven et al. 38 We collect the nodal values of the approximate solution at the N + 1 quadrature points j ( ) within a DG cell (using either Gauss-Legendre or Gauss-Legendre-Lobatto nodes) into the solution vector u = ( . Furthermore, the Lagrange polynomials L k ( ) corresponding to the nodal set are used as the DG test and basis functions and collected into the vector valued function L given by L( ) = (L 1 ( ), … , L N+1 ( )) T . Defining the differentiation matrix D and discrete mass matrix M by their entries D k = L ′ k ( ) and M jk = jk j = M kj , the resulting DG formulation reads whereq andû again represent the numerical diffusion fluxes. We note that in case of Gauss-Legendre-Lobatto nodes, we have a lumped mass matrix M with M k ≈ ∫ 1 −1 L L k d , whereas for Gauss-Legendre nodes, integration is exact, that is,

ENERGY STABILITY OF THE BR2 SCHEME
Recently, Quaegebeur et al 27 investigated a class of energy stable flux reconstruction (ESFR) schemes for the one-dimensional linear diffusion equation using compact numerical fluxes. This class of ESFR schemes is based on the flux reconstruction framework developed by Huynh 28,39 and further advanced by Vincent et al, 40 which specifically contains the nodal DG schemes on Gauss-Legendre and on Gauss-Legendre-Lobatto nodes investigated in this work.
Given the first order equations and the approximation space V h defined in Equation (2) where we use the notation u , with the reference coordinate ∈ [−1, 1] and j denoting the transformation to a specific cell I j as defined in Section 2. Numerical diffusion fluxes are given byq| ±1 = q ± 1 2 ,û| ±1 =û ± 1 2 . The correction functions are specified by two parameters c and , such that where N denotes the Legendre polynomial of degree N and the parameters are given by N,c = . Considering the nodal DG schemes used in this work, the DG scheme on Gauss-Legendre nodes is given by chosing the left and right correction functions as the left and right Radau polynomials for both equations, that is, N, = N,c = 0, whereas the DG scheme on Gauss-Legendre-Lobatto nodes is obtained by chosing N, = N,c = N+1 N , see, for example, Huynh. 39 In particular, Quaegebeur et al 27 show that for ESFR schemes, the BR2 GL diffusion fluxes with exact projection property to calculate the lifting operators l + 1 2 are equivalent to interior penalty (IP) fluxesû IP if the penalty parameter is chosen appropriately. Although the analysis by Quaegebeur et al is given for general non-overlapping one-dimensional grids, for the eigensolution analysis in this work, the grid is assumed to be uniform with constant cell size Δx. From Quaegebeur et al, 27 ORTLEB Furthermore, Equations (3.6) and (3.7) in Quaegebeur et al 27 yield (21) Therefore, the minimum value of the above term is obviously attained for values of such that Inserting this into Equation (20) yields the assertion.
In addition, Quaegebeur et al provide More precisely, on uniform grids, the BR2 GL scheme is equivalent to an IP scheme with = e An extension of the analysis by Quaegebeur et al yields the following result for the BR2 GLL diffusion fluxes. Proof. In case of Gauss-Legendre-Lobatto quadrature, the cell boundaries are contained in the nodal set. Hence, we have where L 1 and L N+1 are the Lagrange polynomials corresponding to 1 = −1 and N+1 = 1 and 1 = N+1 = 2 N(N+1) are the corresponding quadrature weights, see, for example, Abramowitz and Stegun. 41 Inserting the weights and further simplification yields Hence, for ESFR schemes, the BR2 GLL diffusion discretization is equivalent to an IP scheme with = e N(N+1) 2Δx and therefore to a BR2 GL scheme witĥe = N N+1 e . If e ≥ 1, energy stability for BR2 GLL then follows from Theorem 1.
Furthermore, for DG schemes on Gauss-Legendre-Lobatto nodes, the BR1 scheme is again recovered by a particular choice of BR2, because we have Proof. Obviously, due to Theorem 3, proving the equivalence of the BR1 scheme to the BR2 GLL scheme with e = 1 is sufficient. For this purpose, partial integration is applied to the term < u, r x > j in Equation (8) and L 1 and L N+1 are inserted for the test function r.
for the BR1 scheme, this yields the following representation of the cell boundary values q + − 1 2 and q − + 1 2 on the cell I j , , the diffusion flux for the auxiliary variable q at a cell interface, for example, between cells I j and I j+1 , is then given bŷ q BR1

EIGENSOLUTION ANALYSIS
In the following, the accuracy of the DG schemes described in Section 2 is studied via Fourier analysis. The linear advection-diffusion Equation (1) (13), (14) yields with real (N + 1) × (N + 1) matrices A 0 , A −1 , B k , k = 0, ±1, ±2, depending only on the chosen nodal set and the diffusion fluxesû,q which characterize the respective DG scheme. Numerical solutions of the prescribed form can thus be found by solving the eigenvalue problem (27). We may furthermore reduce the set of parameters by defining the nondimensional wave number K = kΔx, the nondimensional numerical frequency Ω = Δxã , and the grid Peclet number Pe * = aΔx d to obtain the eigenvalue problem For a polynomial degree of N, for any K ∈ [0, (N + 1)], we obtain a set of N + 1 eigenvalues of Equation (28). In the literature, see, for example, Hu et al, 13 often only one of these is considered physical and the function assigning to each K its corresponding physical eigenvalue is termed the physical mode, whereas the remaining modes are named spurious. Figure 1A,B plots the real and imaginary part of the three modes for the DG (N = 2) scheme on Gauss-Legendre nodes with LDG diffusion flux for a grid Peclet number of Pe * = 20. For both variants of LDG, the eigenvalues coincide. Figure 1D,E shows the eigenvalue curves for the LDG diffusion flux on Gauss-Legendre-Lobatto nodes, whereas Figure 2 depicts the results for the BR1 and BR2 diffusion fluxes on Gauss-Legendre and Gauss-Legendre-Lobatto nodes for N = 2, respectively. For the BR2, scheme, the penalty parameter is hereby set to e = 3 and on Gauss-Legendre-Lobatto nodes the BR2 GLL variant is implemented. Comparing the choice of DG nodal set in Figures 1 and 2, we see that for the more accurate Gauss-Legendre quadrature rule, the physical mode stays close to the exact dispersion relation for a larger range of wave numbers in all cases. For the BR1 scheme and the BR2 scheme with e = 3, the numerical dispersion relations depicted in Figure 2 differ. In particular, the BR2 approach introduces more numerical dissipation to the physical mode for higher wave numbers as well as to the spurious modes due to the comparatively high value of e . Figure 3 shows the corresponding results when varying the penalty parameter e of the BR2 fluxes in case of the DG scheme on Gauss-Legendre nodes. Reducing e yields results closer to BR1. In fact, for e = N N+1 = 2 3 , that is, the smallest value of the penalty parameter yielding a provable energy stable scheme as described in Section 3, the eigenvalue curves are almost indistinguishable from those of the BR1 scheme. This behavior is not surprising, because the smallest choice of the BR2 penalty parameter in case of the DG scheme on Gauss-Legendre-Lobatto nodes is equivalent to BR1, see Section 3.

Energy distribution of the eigenmodes with respect to the exact solution
Solely regarding the physical mode does not yield a complete picture of the behavior of the numerical scheme. In fact, Moura et al 18 argue that the spurious modes contribute to the overall accuracy, in particular for higher wave numbers where the secondary eigenmodes may strongly influence the behavior of the scheme. Also deviating from the distinction of a physical mode, Watkins et al 33 combine all eigensolutions into an error bound depending on the wave number. We will follow their work, which also allows us to point out differences between the two variants of the LDG approach. For this purpose, we regard the complete set of N + 1 eigenvalues p , p = 1, … , N + 1, and corresponding normalized eigenvectors v p , p = 1, … , N + 1.
An initial wave with nondimensional wave number K on a cell I j , given by the initial nodal values , p = 1, … , N + 1. All eigenmodes are therefore present in the numerical solution of the DG scheme for the above initial condition. As also stated by Watkins et al, 33 the weights p determine how each mode contributes quantitatively to the numerical solution and the normalized weights Hereby, Figure 1C,F shows that although both variants of LDG have the same eigenmodes, the energy distribution among those modes is different for LDG a and LDG b . More precisely, for moderate wave numbers, the first eigenmode-which is often considered the physical one-has more influence on the numerical solution for LDG a than for LDG b . In addition, in this wave number regime, the spurious mode 3 has less influence for LDG a . Furthermore, Figure 2 depicts differences of both the eigenmodes and their energy distribution comparing the BR1 and BR2 diffusion discretization. Here, for low to moderate wave numbers, the BR2 scheme with e = 3 leads to a higher energy content of  the physical mode, whereas for high wave numbers, the energy content of the physical mode is higher for the BR1 scheme. In addition, the BR2 approach introduces a larger amount of numerical diffusion for the spurious mode 3 for this value of e as shown in Figure 2B,E. This mode still has a significant energy content as shown in Figure 2C,F and has therefore a considerable impact on the behavior of the numerical solution. The dependency of the numerical dissipation and the energy content of each mode on the BR2 penalty parameter e is indicated in Figure 3.
Increasing the polynomial degree, Figures 4 and 5 show the four modes for the DG(N = 3) scheme. Considering the two variants of alternating LDG fluxes, for moderate wave numbers the first eigenmode has again more influence on the numerical solution for LDG a than for LDG b , whereas the influence of the spurious mode 4 is significantly reduced for LDG a . Considering the comparison of the BR1 scheme with the BR2 scheme with parameter e = 3, differences regarding the DG nodal set can be perceived. On Gauss-Legendre nodes, the results are similar to the DG(N = 2) scheme relating among the eigenmodes for the discontinuous Galerkin (DG) scheme of order N = 3 using Gauss-Legendre (GL) nodes (first row) and Gauss-Legendre-Lobatto (GLL) nodes (second row). Diffusion term discretized by the two alternate variants of the local discontinuous Galerkin (LDG) diffusion, LDG a (solid lines), and LDG b (dashed lines) [Colour figure can be viewed at wileyonlinelibrary.com] to the higher energy content of the physical mode for the BR2 scheme in case of moderate wave numbers and of the BR1 scheme for high wave numbers. However, on Gauss-Legendre-Lobatto nodes, the BR2 scheme yields higher energy content of the physical mode for high wave numbers as well.

A wave number dependent error bound combining all eigenmodes
As in Watkins et al, 33 we now consider the error of the numerical solution depending on the nondimensional wave number. Hereby, we again consider the above initial condition (29). The numerical solution on cell I j corresponding to this initial condition is whereas the exact solution is given by For the absolute error of the numerical solution, we therefore have and using the norm of the relative error of the numerical solution is given by Whereas Watkins et al 33 use an upper bound of the above quantity (34), we will use the exact relative error with respect to the wave number in the following analysis. Particularly for well-resolved wave numbers, a combined error measure is advantageous because in the depiction of eigenmodes the variations with respect to different schemes are almost indistinguishable. In addition, the behavior of different schemes concerning each separate eigenmode is not as perceptible for higher polynomial degrees. Figures 6-8, and 9 show the relative error as defined in Equation (34) at a given time t versus the nondimensional wave number for the DG schemes up to fifth order using different diffusion fluxes. In particular, the behavior of the alternate LDG a and LDG b fluxes, the BR1 scheme and the BR2 scheme for e = 3 are compared. In addition, in Figure 10, a study with respect to the BR2 penalty parameter is carried out for the DG(N = 2) and DG(N = 3) schemes on Gauss-Legendre nodes.
Considering the two variants of LDG, since the energy distribution among modes differs, the two choices of alternating fluxes lead to significant differences with respect to the relative error. More precisely, the variant which is inconsistent with the upwind advective flux is more accurate in case of advection-diffusion discretizations for a large range of wave numbers.
Furthermore, although higher accuracy of the BR1 flux for well-resolved solutions has been observed in Watkins et al, 33 the plots in logarithmic scale at the right-hand side of each figure indicate that for very low wave numbers, this higher accuracy is restricted to the DG scheme on Gauss-Legendre nodes with even polynomial degree N. For the DG scheme on Gauss-Legendre-Lobatto nodes as well as for odd N, both the BR2 flux with penalty parameter e = 3 and the more accurate variant of LDG produce a smaller error for well-resolved solutions.
In  LDG a scheme is the most accurate approach, followed by the BR2 scheme with penalty paramter e = 3, the BR1 scheme and lastly the LDG b scheme. For N = 1, the error curves in Figure 6C are nearly indistinguishable but suggest a similar sequence for the schemes. Concerning the BR2 penalty parameter e , Figure 10 depicts the relative errors of the DG approximate solutions on Gauss-Legendre nodes for polynomial degrees of N = 2, 3 in case of the parameter choices e = N N+1 , 2, 3. For comparison, the results using the BR1 diffusion fluxes are included in the error plots. Obviously, the relative error of the BR2 scheme is closest to the relative error of the BR1 scheme for e = N N+1 . At this point, it should be remarked that even if the results for BR1 and BR2( e = N N+1 ) are similar, they do not coincide. In fact, on Gauss-Legendre nodes, the BR1 scheme is not equivalent to BR2 for any value of the penalty parameter e , since on this set of nodes, the BR1 diffusion discretization has a wider stencil than the BR2 and IP schemes. This does not contradict Theorem 4, because its assertions only apply to the DG scheme on Gauss-Legendre-Lobatto nodes.
Furthermore, with respect to the parameter study for e , it can be observed that BR2( e = N N+1 ) and BR2( e = 3) provide lower and upper bounds for the BR1 results as well as for BR2( e = 2). Along the whole range from low to high wave numbers, the BR2 schemes for e = N N+1 and e = 3 alternate to provide the most accurate solution among the schemes considered in Figure 10. Therefore, there is no optimal choice for e regarding accuracy for the whole range of low to high wave numbers. However, we observe that for very high wave numbers, the BR2( e = N N+1 ) and BR1 schemes yield the most accurate results for both values of N, whereas for well-resolved wave numbers, the most accurate scheme depends on the polynomial degree. For even polynomial degree N = 2, the parameter e = N N+1 yields the most accurate and e = 3 the less accurate result whereas for odd polynomial degree N = 3 the roles of the particular BR2 schemes are reversed and BR2( e = 3) is provides the best result. This odd-even phenomenon is investigated more thoroughly in Section 5.

An example regarding the influence of the grid Peclet number
In order to give an example of the influence of the grid Peclet number on the results of the eigensolution analysis, in particular for higher order DG schemes, Figure 11 shows the eigenmodes for the two alternate variants of the LDG diffusion fluxes within the DG(N = 5) scheme on Gauss-Legendre nodes for Pe * = 20 as well as Pe * = 100. In addition, Figure 12 gives a comparison of the numerical errors depending on the nondimensional wave number. Although the structure of the eigensolutions for the LDG fluxes differs for different Peclet numbers, with a high frequency mode for Pe * = 20 and replications of the physical mode for Pe * = 100, the basic observations regarding the differences with respect to diffusion discretizations still hold. Regarding the LDG variants, for moderate wave numbers, the first eigenmode has more influence on the numerical solution for LDG a than for LDG b . In addition, Figure 12 shows that LDG a is more accurate than LDG b . Regarding the results with respect to the BR1 scheme and the BR2 scheme for e = 3, which are also depicted in in Figure 12, the error plots in logarithmic scale indicate a smaller error for the BR2 scheme with penalty parameter e = 3 in comparison with the BR1 diffusion fluxes for both Peclet numbers in accordance to the observed reduction of accuracy for the BR1 scheme for odd polynomial degrees. A comparison of Figure 12A,C shows that for both Peclet numbers, the BR1 scheme is alternately the most or least accurate scheme depending on the wave number. A different among the eigenmodes for the discontinuous Galerkin (DG) scheme of order N = 5 using Gauss-Legendre (GL) nodes for Pe * = 20 (first row) and Pe * = 100 (second row). Diffusion term discretized by the two alternate variants LDG a (solid lines) and LDG b (dashed lines) [Colour figure can be viewed at wileyonlinelibrary.com] behavior depending on the grid Peclet number can be perceived for the LDG b scheme, which produces comparatively larger errors in the higher wave number range for Pe * = 100. Naturally, the differences between the diffusion schemes are generally smaller for higher Peclet numbers because then the influence of advection dominates regarding the underlying PDE. In that respect, the LDG a scheme and the BR2 scheme for e = 3 produce nearly identical results in Figure 12C,D for Pe * = 100.

NUMERICAL RESULTS
In this section, numerical experiments are carried out to solve the advection-diffusion problem (1) for a = 1 on the spatial domain = [−10, 10] discretized by E = 20 elements, that is, the grid size is Δx = 1. The diffusion coefficient d varies depending on the respective test case and can be obtained from the grid Peclet number, which in this case is We consider periodic initial conditions and supplement (1) with periodic boundary conditions. Now, the DG scheme (13), (14) is applied to this periodic advection-diffusion problem for polynomial degrees of N = 1, … , 5. For time integration, the classical explicit fourth-order Runge-Kutta scheme is used. Given the DG solution u(t) and the exact solution u ex (t), we measure the relative L 2 -error Hereby, the L 2 -norm of a nodal quantity v is approximated using the given quadrature rules as where , = 1, … , N + 1, denote either the Gauss-Legendre or the Gauss-Legendre-Lobatto weights corresponding to the DG nodal set employed.    1.17e+00 1.23e+00 1.22e+00 1.19e+00 1.18e+00 1.17e+00 9 ∕10 1.53e+00 1.54e+00 1.53e+00 1.52e+00 1.52e+00 1.52e+00 Abbreviation: LDG, local discontinuous Galerkin.

Evolution of a single wave
The purpose of this first test case is to verify the error analysis in Section 4. We consider initial solutions of the form for various sizes of the wave number K. Tables 1 and 2 list the L 2 -errors of the numerical solution versus the nondimensional wave number K N+1 for the DG schemes for N = 1, … , 5 on Gauss-Legendre and Gauss-Lobatto nodes, respectively. The simulations are run until final time t = 0.05, and the grid Peclet number is set to Pe * = 20, according to the eigenanalysis in Section 4. For this test case, the DG scheme is combined with both variants of the alternate LDG diffusion flux as well as the BR1 and BR2 diffusion discretization. The study in Section 4 regarding the BR2 penalty parameter is continued by chosing e = N N+1 , 2, 3 for the DG scheme on Gauss-Legendre nodes and e = 2, 3, 3(N+1) N for the DG scheme on Gauss-Legendre-Lobatto nodes using the BR2 GLL implementation described in Section 2. As shown in Section 3, the BR2 GLL ( e = 3 N+1 N ) scheme is equivalent to BR2 GL ( e = 3) and for the DG scheme on Gauss-Legendre-Lobatto nodes, BR2 GLL ( e = 1) is equivalent to the BR1 scheme. Regarding the numerical results listed in Tables 1 and 2, the following observations can be made: • In almost all set-ups, the LDG a variant performs better than LDG b . • In case of Gauss-Legendre nodes, for all N, there is some range of wave numbers where the BR1 scheme and the BR2 scheme for 3 = N N+1 beat the LDG a discretization, and the BR2

A well-resolved test case
As in Watkins et al, 33 we now consider the solution of the 1D advection-diffusion equation with a well-resolved approximate Gaussian as initial condition. Hereby, the initial condition u(0) as well as the exact solution u ex (t) can be computed from the analytical solution where is the -th spectral weight given below,k = 2 ∕L is the -th wave number associated with the domain length L = 20, and N̂k is the number of waves used, where N̂k is chosen the largest positive integer such thatk N̂k ≤ (N + 1) ∕Δx.
• for odd polynomial degrees N = 1, 3, 5, the LDG a diffusion discretization performs best for all investigated output times and both nodal sets except for the case N = 1, t = 1 on Gauss-Legendre-Lobatto nodes. Furthermore, the BR2 schemes with large penalty parameters e = 2, 3 and e = 3(N+1) N beat the BR1 scheme in all odd degree cases, with errors decreasing for increasing e , except for the DG(N = 1) scheme on Gauss-Legendre-Lobatto nodes and output time t = 10. In addition, on Gauss-Legendre nodes, the BR2 scheme with small penalty parameter e = N N+1 yields larger errors than BR1 for N = 3, 5.
• for even polynomial degrees N = 2, 4, the situation above is reversed, that is, the BR2 scheme for e = N N+1 performs best of all diffusion discretizations in case of Gauss-Legendre nodes followed by the BR1 scheme. Furthermore, regarding the BR2 penalty parameter, increasing e increases the error of the numerical solution. For the even degree case and Gauss-Legendre-Lobatto nodes, either the LDG a variant or the BR1 scheme yield the lowest numerical error, also depending on the output time t.

A low-resolution test case
For this test case, basically, the same set-up as in Section 5.2 is used. However, the standard deviation is reduced to = 1∕ √ 2 to produce a poorly resolved initial condition. The last five rows of Tables 3 and 4 contain the relative L 2 -errors at output times t = 0.1, 1, 10 for this low-resolution test case. The results show that • regarding the two alternative LDG variants, LDG a yields lower numerical errors in all set-ups except for the DG(N = 4) scheme on Gauss-Legendre nodes at t = 0.1 where both variants yield almost the same error. • On Gauss-Legendre-Lobatto nodes, the behavior of the schemes with respect to accuracy corresponds to the error plots in Section 4 for low output times t = 0.1 (for all polynomial degrees) and t = 1 (for N = 2, 4). More precisely, in these cases, the most accurate diffusion scheme is LDG a , followed by the BR2 schemes for e = 3(N+1) N , e = 3 and e = 2, in this order, whereas the second largest error is given by the BR1 scheme and the largest one by the LDG b approach.
• On Gauss-Legendre nodes, the most accurate approach is LDG a except for the three cases N = 1, t = 1; N = 3, t = 10 and N = 4, t = 0.1. Furthermore, considering the BR1 and BR2 approaches and even polynomial degrees N = 2, 4, the order of these schemes with respect to accuracy is BR2 for e = 3, 2, BR1, and lastly BR2 for e = N N+1 .

CONCLUSION
In this work, the dissipation and dispersion properties of DG schemes for the advection-diffusion equation in one space dimension have been compared with respect to different interface fluxes for the diffusion term. In particular, the consistent and the inconsistent variant of alternate LDG fluxes, the BR1 flux, and the BR2 scheme depending on a penalty parameter e have been employed. For the BR2 flux, recent results on energy stability by Quaegebeur et al 27 have been reviewed and extended, resulting in a simplified representation of the lower bound e = N N+1 in case of the calculation of the BR2 lifting operator by exact projection and the equivalency results in Theorems 3 and 4.
A comparison of the wave propagation properties of the DG scheme using BR2 diffusion fluxes with respect to e shows a similarity to the results of the BR1 flux for low values of e . In addition, there is no optimal choice of e which provides the most accurate result for all wave numbers and polynomial degrees. In fact, for N = 2, 3, an alternating behavior has been observed, where the results for e = N N+1 and e = 3 alternate to provide the smallest error. Considering the alternate LDG fluxes, the performance of LDG a is in general more favorable compared with LDG b , both for the well-resolved problem and the low-resolution test case, independent of the polynomial degree, and the nodal DG set.
For well-resolved wave numbers, the observation of higher accuracy of the BR1 scheme compared with LDG generally only holds in case of Gauss-Legendre integration nodes and even polynomial degree of the DG approximate solution. Furthermore, implemented in the DG scheme on Gauss-Legendre nodes, the BR2 schemes with large penalty parameters e = 2, 3 and e = 3(N+1) N beat the BR1 scheme in all odd degree cases for the well-resolved problem, with errors decreasing for increasing e . An odd-even phenomenon is observed in case of Gauss-Legendre nodes, that is, for even polynomial degree, the BR2 scheme for e = N N+1 performs best of all diffusion discretizations, followed by the BR1 scheme, and increasing the BR2 penalty parameter increases the error of the numerical solution.
For the low-resolution test case, the LDG a scheme yields the lowest error in most cases, which corresponds to the findings in the paper by Watkins et al. 33 Future work should also consider the Fourier analysis of the corresponding fully discrete schemes. Hereby, the dissipation and dispersion properties with respect to the diffusion discretization will surely depend on implicit or explicit time discretization of the diffusion terms.