Superconvergence and the Numerical Flux: a Study Using the Upwind-Biased Flux in Discontinuous Galerkin Methods

One of the beneficial properties of the discontinuous Galerkin method is the accurate wave propagation properties. That is, the semi-discrete error has dissipation errors of order $$2k+1$$2k+1 ($$\le Ch^{2k+1}$$≤Ch2k+1) and order $$2k+2$$2k+2 for dispersion ($$\le Ch^{2k+2}$$≤Ch2k+2). Previous studies have concentrated on the order of accuracy, and neglected the important role that the error constant, C,  plays in these estimates. In this article, we show the important role of the error constant in the dispersion and dissipation error for discontinuous Galerkin approximation of polynomial degree k,  where $$k=0,1,2,3.$$k=0,1,2,3. This gives insight into why one may want a more centred flux for a piecewise constant or quadratic approximation than for a piecewise linear or cubic approximation. We provide an explicit formula for these error constants. This is illustrated through one particular flux, the upwind-biased flux introduced by Meng et al., as it is a convex combination of the upwind and downwind fluxes. The studies of wave propagation are typically done through a Fourier ansatz. This higher order Fourier information can be extracted using the smoothness-increasing accuracy-conserving (SIAC) filter. The SIAC filter ties the higher order Fourier information to the negative-order norm in physical space. We show that both the proofs of the ability of the SIAC filter to extract extra accuracy and numerical results are unaffected by the choice of flux.


Introduction
There have been many studies on the accurate wave propagation properties for discontinuous Galerkin methods [3,9,10,14,19,21,24,25]. These studies conclude that, in general, the semi-discrete error has dissipation errors of order 2k + 1 ( ≤ Ch 2k+1 ) and order 2k + 2 ( ≤ Ch 2k+2 ) for dispersion. These studies concentrate on the order of the error and ignore the role of the error constant. In this article, we give an explicit formula for how the error constant in the dispersion and dissipation error depends on the numerical flux as well as the polynomial degree. Specifically, when using a piecewise constant or quadratic approximation, the results suggest that one should use a more central flux to optimise the numerical wave propagation properties. For a piecewise linear or cubic approximation, the results confirm that an upwind flux is desirable. These results are illustrated through the use of the upwind-biased flux introduced by Meng et al. [15], which uses a convex combination of the upwind and downwind approximation values. In addition to this result, we introduce results for post-processing the numerical solution using the smoothness-increasing accuracy-conserving (SIAC) filter by Ryan and others [11-13, 18, 20]. This filter allows for extracting the higher order Fourier information that is used in the proof of the dissipation and dispersion errors and translates this information into physical space. To complete our discussion of superconvergence, we also review the pointwise superconvergence results for of Cao et al. [5] for the upwind-biased flux.
To summarize the contribution of this article, it is to (i) show the important role the error constant plays in the dissipation and dispersion error and its significance in tying the choice of flux to even versus odd polynomial degree approximations; (ii) show how this Fourier-type information can be extracted via SIAC filters; and (iii) tie the superconvergent information together through a review of pointwise superconvergent results.
The discontinuous Galerkin (DG) methods are one of the most extensively researched classes of numerical methods for solving partial differential equations that display convective or diffusive qualities and have been popularly adopted by the scientific and engineering communities as a method capable of achieving arbitrary orders of accuracy in space. The choice of the numerical flux function plays a pivotal role in the successful construction of DG methods and has an intrinsic effect on the superconvergence properties. The vast majority of theory for DG schemes for conservation laws has been developed with the (somewhat habituated) choice of a monotone numerical flux. Design and criteria for the selection of numerical flux functions is an area with a great deal of scope for further investigation. One approach is to not only explore the order of the dissipation and dispersion errors, but also the error constant to ensure that the errors from dispersion and dissipation are minimized.
For a Runge-Kutta (RK) DG solution to the linear advection equation, one can do better than the "natural" upwind flux. Recently, Meng et al. [15] introduced in the context of DG methods for linear hyperbolic equations a more general flux function: the upwind-biased flux. This function parametrises the ratio of information taken from the left compared to the right of cell interfaces. This choice avoids the requirement of exact knowledge of the eigenstructure of the Jacobian and may reduce numerical dissipation (yielding a better approximation in smooth regions) but it is made at the cost of the loss of monotonicity. These results provide the theoretical foundations for our investigations into the choice of flux function for RK-DG methods.
Recent interest [3,9,10,14,19,21,24,25] in analysis via a Fourier approach of DG solutions to the linear advection equation offers an alternative means by which to explore superconvergence. This analysis is limited to linear equations with periodic boundary conditions and a uniform mesh. However, as is justified by numerical examples [9], the results provide a guide for the behaviour of solutions in a more general 1 3 setting. Stability and (k + 1)-th order accuracy can be established via this approach. Furthermore, this analysis is useful for exploring the time propagation of errors. In [9], it was shown that a kth-order DG solution to the linear advection equation has one physical mode and k spurious ones which are damped exponentially fast over time, where C 1 , C 2 , C 3 , C ∈ ℝ + . The first term on the right-hand side of inequality (1), which dominates for T = O(1∕h −k−1 ), is attributed to the dispersion and dissipation errors of the physically relevant eigenvalues and grows linearly in time. The expected order of accuracy is = k + 1; however, a judicious choice of the numerical flux function in the semi-discrete DG scheme can yield a superconvergent order of accuracy as high as = 2k + 1. The third term, which decays exponentially fast over time with respect to h, accounts for dissipation of the spurious modes. The second term is due to projection of the initial condition and does not grow in time. Thus the error is on the order of 2k + 1 for long time integration but only k + 1 over short time. At certain special points (the superconvergent points which change with the choice of numerical flux) the accuracy of inequality (1) can be increased to O(h p+2 ) by carefully interpolating the initial projection.
The above superconvergence results can be extracted from the underlying approximation to achieve 2k + 1 in the L 2 -norm. This involves convolving the approximation against a specially designed kernel comprising a linear combination of B-splines [7,11,12,18,20], effecting increased smoothness by damping the non-physical eigenmodes of the DG operator and exploiting information concealed in the unwelcome fluctuations that characterise the numerical solution. We are able to extend the analysis of the SIAC filtered error, which is facilitated by a dual analysis in a similar fashion to [12] for the case of the upwind-biased flux and is limited to a contribution to the constant attached to the post-processed error term.
The outline of this paper is as follows: in Sect. 2, we discuss the preliminaries needed to perform both the pointwise and negative-order norm error analysis as well as reviewing the upwind-biased DG scheme. This includes a review of the pointwise superconvergent results for the upwind-biased DG scheme. In Sect. 3 we discuss the dispersion analysis related to this scheme and how the constants depend on the choice of the upwind biasing in the flux. In Sect. 4 we extend the negative-order norm error analysis to upwind-biased DG schemes and show that this simply affects the constant in the error. We end by supporting this analysis with numerical examples in Sect. 5 and conclusions in Sect. 6.

Preliminaries
We begin with some definitions used in the error estimates for discontinuous Galerkin solutions and SIAC filtering and complete the prefatory construction of the DG scheme.

Notation and Definitions
We first treat the discretisation of a multi-dimensional spatial domain and define the associated approximation spaces.
Consider the linear hyperbolic system for the conserved quantity u(x, t) , where x = (x 1 , ⋯ , x d ) ∈ Ω ⊂ ℝ d and where we assume that we have real constant coefficients a i ≥ 0 . We adopt throughout this paper the assumptions of a smooth initial condition and of periodicity in the boundary conditions: For much of the error analyses, it is required only that u 0 (x) ∈ H k+1 (Ω) but for the proof of Theorem 2.3, we require infinite differentiability to write the DG solution as a Maclaurin series.

Tessellation
Let Ω h be a tessellation of a d-dimensional bounded domain Ω into elements S of the regular quadrilateral-type shape. Denote by Ω h = ⋃ S∈Ω h S the union of boundary faces S of the elements S ∈ Ω h . A face e internal to the domain has associated with the "left" and "right" elements S L and S R and exterior-pointing normal vectors L = (n L 1 , ⋯ , n L d ) and R = (n R 1 , ⋯ , n R d ) , respectively, as described in [12]. Given a function v defined on neighbouring elements S L and S R which share a face e, we refer to its restriction in S L to the face e by writing v L ∶= (v| S L )| e and similarly for v R , the restriction of v to e in S R .
For clarity, we detail the discretisation of a two-dimensional domain Ω = [−1, 1] 2 into N x ⋅ N y rectangular cells. Elements take the form S = I i × J j where I i and J j are intervals given by for 1 ≤ i ≤ N x and 1 ≤ j ≤ N y respectively, where and Denoting by h , 1 ≤ j ≤ N y , we define and require regularity: h x,i , h y,j ≥ ch, c > 0. For simplicity of the filtered error estimates in Sect. 4, and by necessity in Sect. 3, we will adopt constant element lengths , y) , for example, is written as

Basis Polynomials
We use as basis functions the Legendre polynomials P n ( ) , which can be defined by the Rodrigues formula and which satisfy the orthogonality condition where nm is the Kronecker-delta function. We can then define the right and left Radau polynomials as respectively. It is known that the roots of R + k+1 ( ) and the roots of R − k+1 ( ) are real, distinct and lie in the interval [−1, 1].

Function Spaces
Due to the tensor product nature of the post-processing kernel, we will require the function space of tensor product polynomials Q k (S) of degree at most k in each variable. Thus, we base the proofs in this paper on the following finite element spaces: where L 2 (Ω) is the space of square-integrable functions on Ω . Nevertheless, we mention here that it has been observed ( [18]) that the filter also works for the standard polynomial space P k (S) . Note that for a one-dimensional domain Ω , these function spaces Q k (S) and

Operators on the Function Spaces
We list the following standard notations. The inner-product over Ω of two functions is defined as depending on whether the functions take scalar or vector values. We denote by ℙ h v and by Π h p the usual L 2 -projections of scalar and vector valued functions v and p , respectively. The L 2 -norm on the domain Ω and on the boundary Ω is defined as and the -norm and semi-norm of the Sobolev space H (Ω) are defined, respectively, as where is a d-dimensional multi-index of order | | and where D denotes multi-dimensional partial derivatives. The definitions for the above norms for vector-valued functions are analogous to the scalar case. In Sect. 4, we will utilise the negative-order Sobolev norm, defined as as a means of obtaining L 2 -error estimates for the filtered solution. Note that for all ≥ 1 , we have The negative order norm can be used to detect oscillations of a function [7] and is connected to the SIAC filter which smooths oscillations in the error. Finally, the difference quotients h,j v are given by where j is the jth component unit normal vector. For any multi-index = ( 1 , ⋯ , d ) , we define the th-order difference quotient by 1 3

Construction of the DG Scheme
Given a tessellation Ω h of the domain Ω , the method is facilitated by multiplying Eq. (2) by a test function v and integrating by parts over an arbitrary element S ∈ Ω h to obtain Next, we define the piecewise polynomial approximation space as Note that functions v ∈ V k h are allowed to be discontinuous across element boundaries. This is the distinguishing feature of DG schemes amongst finite element methods.
By replacing in Eq.
h , we obtain the DG method: find, for any v ∈ V k h and for all elements S, the unique where f is a single-valued numerical flux function used to enforce weak continuity at the cell interfaces. The initial condition u h (x, 0) ∈ V k h is usually taken to be the L 2 -projection ℙ h u 0 , although the analysis in Sect. 2.3 favours a function which interpolates u(x, 0) at the superconvergent points.
Summing Eq. (4) over the elements S, we get a compact expression for the global scheme: where we define for future use Remark 2. 1 We note that although this article considers the semi-discrete form of the DG method, the results of this article extend to fully discrete schemes such as the Lax-Wendroff DG method, provided the higher order derivatives of the fluxes are discretized using the LDG formulation as in.

Flux Function
To ensure stability of the scheme (4), it remains to define the numerical flux functions f featured in the cell boundary terms. In general, numerical solution from both sides of the cell interface. Traditionally [8], this function is chosen to be a so-called monotone flux, such as the Lax-Friedrichs flux, which satisfies the Lipschitz continuity, consistency ( f (u, u) = f (u) ) and monotonicity ( f (↑, ↓) ). For our test Eq. (2), where the linear flux f (u) = au determines a single wind direction, the usual choice in the literature is to satisfy the upwinding condition. In this paper, where f (u h ) = a� u h , we choose instead the upwind-biased flux defined here in the one-dimensional case with periodic boundary conditions, which was recently described in the context of DG methods by Meng et al. [15]. More information is taken from the left than from the right of cell boundaries and, when = 1 , the upwindbiased flux reduces to the purely upwind flux u − h . We do not allow = 1 2 , which gives a central flux, since then the scheme becomes unstable.
For clarity, we particularise the evaluation of this flux at cell boundary points. In two dimensions, the upwind-biased flux takes the form where 1 > 1 2 , and, similarly for 2 > 1 2 , Choosing, over the upwinding principle, the upwind biased flux offers several rewards [15]: a possibly reduced numerical viscosity and easier construction, for example. However, as a price paid for introducing the parameter , we sacrifice the established property of monotonicity of the flux function. In this paper, we consider in terms of superconvergence the severity of this loss of monotonicity.

Pointwise Superconvergence
In this subsection, we review current literature for the upwind-biased flux for DG methods [5]. In this case, the leading order term in the error is proportional to a sum, dependent upon , of left and right Radau polynomials. The main result, Theorem 2.3, is an extension of the observation, for example of Adjerid et al. [1,2], that the superconvergent points for the purely upwind DG scheme are generated by roots of right Radau polynomials. To be able to illustrate this idea, a "special" Radau polynomial can be defined as The main idea is that the roots of R ⋆ k+1 ( ) , which change with the value of , generating superconvergent points on the order of h k+2 in the element interior for the upwind-biased scheme.

Remark 2.3
Recall that when = 1 , one of the superconvergent points is the downwind end ⋆ k+1 = + k+1 = 1 . When k is even, the superconvergent points shift to the left with decreasing values of < 1 . On the other hand, when k is odd, the points shift to the right and ⋆ k+1 > 1 . For example, when k = 1 , the roots of R ⋆ k+1 ( ) are given by Shown in Table 1 are approximate values of the roots of R ⋆ k+1 ( ) for two values of . One can observe that for the lower value of , the roots shift to the left or right depending on k and that, for k = 1 and k = 3 , one of the roots is indeed greater than 1.
Following the lines of [1,5] interpolates the initial condition at roots of R ⋆ k+1 ( ) , where k is even. Lemma 2.1 dictates that it is not possible to obtain in the same way a kth degree polynomial interpolating u at all k + 1 roots of R ⋆ k+1 ( ) when k is odd; there are only k roots inside [−1, 1] . Instead, one can needs to define a global projection. Here, we provide a restatement of Lemma 3.5 in [5]. .

Then, the interpolation error satisfies
where Q ( ) is a polynomial of degree at most .
The interpolatory polynomials described in Lemma 2.2 are used as initial conditions in the proof of Theorem 2.3. Numerical results in Sect. 5 confirm that, in general, there are only k superconvergent points in each element when k is odd.
Here, we provide a restatement of Lemma 3.8 in [5] and Theorem 3.12. (2) with d = 1 obtained by a DG scheme (4) using kth order basis functions, a uniform mesh and the upwind-biased flux û h . Let the numerical initial condition be the interpolating polynomial ⋆ u(x, 0) described in Lemma 2.2.

Theorem 2.3 Consider the approximate solution u h to the one-dimensional linear hyperbolic conservation law
) be the scaling between the cell I j and the canonical element where c k+1 depends on u, h, k and t.

Remark 2.4
When u h (x, 0) = ⋆ u 0 (x) interpolates u 0 (x) at the roots of R ⋆ k+1 (x) , the coefficient of the term on the order of h k+1 in the series for the initial error satisfies In the proof of Theorem 2.3, this relation was extended to t > 0 . For odd k, In the next section, we will compliment the pointwise observations we have made with an analysis of the constants in the dispersion and dissipation error. As shown in the literature, these are of order 2k+ and 2k + 1 , respectively. We can then extract the full O(h 2k+1 ) superconvergence rate using the SIAC filtering. First, we review the construction of SIAC filtering kernels.

Smoothness-Increasing Accuracy-Conserving (SIAC) Post-processing
The hidden local accuracy of the DG solution, discussed in the previous section, may be extracted to a global measure by applying the SIAC filter introduced by [20]. In this section, we show that superconvergent accuracy on the order of h 2k+1 in the negative order norm, as is observed [11] for the upwind flux, still occurs when the upwind-biased DG method is used to solve linear hyperbolic conservation laws. To begin with, we observe that an error bound in the L 2 -norm follows from a negative order norm error estimate.
Here, we detail the component parts of the SIAC filter as defined in, for example, [11]. For a multi-index , we define and, given a point In this way, we construct a convolution kernel which comprises a linear combination of r + 1 B-splines ( ) ∈ C −2 of order such that has compact support and reproduces (by convolution) polynomials of degree strictly less than r. Typically, r = 2k and = k + 1 , where k is the degree of the polynomial basis. In Fourier space, the kernel can be written as as given in [13]. A further order of accuracy in Fourier space can be gained by requiring that the kernel is symmetric. This extra order of accuracy will only reveal itself in physical space if the negative-order norm estimates are of the same order. Notice, that the r + 1 order accuracy does not rely on the smoothness, , of the B-spline. The coefficients c are tensor products of the coefficients c found by requiring the reproduction of polynomials property K (r+1, ) h ⋆ x p = x p , p < r, in the one-dimensional case. It is important to note that derivatives of a convolution with this kernel maybe written in terms of difference quotients: Further properties of the kernel may be found in [12]. The SIAC filtered solution can be written as where the kernel has been applied by convolution at the final time. The filtered solution u ⋆ h (x, t) displays increased accuracy and oscillations in the error are reduced. The results in this paper treat only the symmetric kernel where the nodes are uniformly spaced. Extension to the one-sided filter given in [11] and [20] is a straight-forward task.

Dispersion Analysis
A further approach to the analysis of these methods which has proved fruitful [6,9,25] in recent years involves computing the eigenstructure of the amplification matrix. The choice of initial condition and set of basis functions can be crucial in obtaining optimal results. Recent work that demonstrates the importance of this choice includes [4,6,23]. In what follows, we analyse the eigenvalues, which are independent of the choice of basis.
Consider the local DG solution to Eq. (2) with d = 1 , periodic boundary conditions and a uniform mesh and let be the coefficient vector for basis polynomials j (x) on cell I j of degree at most k. Then, the DG scheme can be written as the following semi-discrete system of ODEs: where A = A 1 + A 2 , B and C are (k + 1) × (k + 1) matrices. Note that the term Cu j+1 from the right neighbour cell is a new contribution when < 1.
Continuing to follow the lines of [9,25], the coefficient vectors can be transformed to Fourier space via the assumption where is the wave number, i = ) is the element center, to obtain a global coefficient vector û . Substitution of the ansatz (13) into the scheme (12) gives a new ODE: where is called the amplification matrix. If G is diagonalisable, then it has a full set of eigenvalues 1 , ⋯ , k+1 and corresponding eigenvectors Λ 1 , ⋯ , Λ k+1 and the general solution of (14) takes the form where the constants C j , j = 1, ⋯ , k + 1 depend on the initial condition. Via this approach one can inspect the size of the error and its behaviour over time as in [9,25]. In this paper, we limit ourselves to an asymptotic expansion of the eigenvalues.
Using Mathematica to computationally perform an asymptotic analysis on = h = 0 , we can obtain the following sets of eigenvalues j of the amplification matrices G for the cases k = 0, 1, 2, 3. For For the case k = 3 , the algebra involved in computing the roots of R ⋆ k+1 and in computing the eigenvalues of the amplification matrix becomes prohibitively substantial and the need to evaluate components numerically makes it particularly difficult to obtain tidy closed expressions for the coefficients. For k = 3 , we have found that one eigenvalue satisfies while the leading order term of the other eigenvalues has a negative real part on the order of 1 h . For each value of k, the eigenvalue 1 has physical relevance, approximating −i with dispersion error on the order of h 2k+1 and dissipation error on the order of h 2k+2 . This is consistent with the previous findings of [3,9,10,19,21,25]. The coefficient of the leading order real term of the physically relevant eigenvalues 1 is negative. While for even k this coefficient vanishes in the limit → 1 2 , for odd k, due to the factor (2 − 1) −1 , it grows without bound with reducing values of . This (blow-up) behaviour is amplified in the coefficient of h 2k+2 . Note that such differences between odd and even k do not manifest when = 1 since then 2 − 1 = 1.
The remaining eigenvalues are non-physically relevant but have negative real part on the order of 1 h so that the corresponding eigenvectors in the solution (15) are sufficiently damped over time. This dampening is observed to be slower in the cases k = 1, 2 for lower values of [6]. The leading order terms vanish only when = 1 when the interpolation points coincide with the superconvergent points of the scheme. Numerical results suggest that for k = 2 , when we have k + 1 roots of R ⋆ k+1 ( ) , we are able to obtain the optimal O h 2k+1 accuracy using u h (x, 0) = ⋆ k+1 u 0 (x). , .

Extracting Superconvergence
To investigate the SIAC filtered error, denote by e h = u − u h the usual DG error and let be the DG solution to Eq. (2) post-processed by convolution kernel. We have the following estimate: Theorem 4.1 Let u h be the numerical solution to the linear hyperbolic conservation law (2) with the smooth initial condition obtained via a DG scheme (4) with upwind-biased flux. Then, for r = 2k and = k + 1.
The first term on the right-hand side of (16) is bounded by Ch r+1 from the integral form of Taylor's theorem and from the reproduction of polynomials property of the convolution (Lemma 5.1, [11]). Thus we need only to consider the second term for which by kernel properties of the th derivative D , the kernel's relation to the divided difference and by Young's inequality for convolutions. The tilde on K h in inequality (17) signals that the kernel uses ( − ) , which is a result of the property D ( ) = h ( − ) . Note that ��K h �� = ∑ r i=0 �c i � is just the sum of the kernel coefficients so we only need to show that || h e h || − ≤ Ch 2k+1 . Furthermore, the formulation of the DG scheme for the solution is similar to that for the divided differences and, as speculated in [7], This allows us to only have to consider the negative order norm of the solution itself; superconvergent accuracy in the negative order norm gives superconvergent accuracy in the L 2 -norm for the post-processed solution. The following result provides the required negative order norm error estimate.

Remark 4.1
Notice that the superconvergent points for the upwind-biased scheme, as described in the one-dimensional case in Lemma 2.2, change with the value of . However, the global superconvergence in the negative-order norm occurs regardless of the value of . (2) with the smooth initial condition obtained via a DG scheme (4) with the upwind-biased flux. Then,

Theorem 4.2 Let u h be the numerical solution to the linear hyperbolic conservation law
Proof For simplicity, we consider the case when = 0 . The case for > 0 is similar [22].
To extract information about the error at the final time, we work with the dual equation: find a continuous and analytic (x, t) such that If we test Eq. (20) against the solution u to the original system (2) and, similarly, multiply the original system by and integrate over the domain, integrating by parts and adding the two new equations yields A periodicity relation is obtained by integrating with respect to t and appealing to the Fundamental Theorem of Calculus: On the other hand, we can write instead Thus, we can estimate the term appearing in the definition of the negative order norm: Information about the model problem and the method can be incorporated into the second term on the right-hand side of Eq. (21) using the DG formulation and the dual equation: Hence we can rewrite Eq. (21) as where (20)

Remark 4.2
The effect of the introduction in the flux function of the parameter is limited to a contribution to the constant attached to the order term in the negative order norm error estimate. This is in contrast to the changing local behaviour seen in the pointwise analysis in Sect. 2.3. Despite the observations of Remark 2.4 regarding the local behaviour when k is odd, we can extract the same global order of accuracy, O(h 2k+1 ) , for any polynomial degree k.

Numerical Experiments
We present a numerical discussion for the test equation . The error curves cross the zero axis near these roots. In [2], Adjerid et al. commented that the intersection points align more closely as k increases and we observe that here as well.
For k = 2 we observe k + 1 superconvergent points while for k = 1 and k = 3 , in general, the error curves cross the zero axis only k times. Furthermore, as the value of reduces, we see an overall reduction in the magnitude of the errors for k = 2 . On the other hand, when k = 1 or k = 3 the magnitude of the errors in general increases for smaller values of . Inside certain anomalous elements, for example the fifth and tenth elements in Fig. 2, the curves miss the crosses or, as in the sixth element in Fig. 3, we observe an additional intersection and this may be due to the initial condition sin(x). Tables 2, 3 and 4 illustrate the order h k+1 accuracy of the DG solution in the L 2 -and L ∞ -norms. After post-processing by the SIAC filter, we observe the O h 2k+1 accuracy in the L 2 -norm described in Sect. 4 and we also see O h 2k+1 accuracy in the L ∞ -norm.
For odd k, convergence to the expected orders is slower for lower values of but is eventually achieved. Furthermore, if one compares the same degrees of mesh refinement for decreasing values of , one observes increasing errors for k = 1, 3 and reducing errors for k = 2 . For the post-processed solution, this is due in large part to the contribution of to the constant attached to the order term in the error estimate of Theorem 4.2.
The highly oscillatory nature of the DG solution, indicating the existence of the hidden superconvergent points, can be seen in Figs. 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 and 15 alongside the post-processed solutions which have increased smoothness and improved accuracy. The reduced numerical viscocity enforced by the upwind-biased flux is evident when comparing plots for = 1 and = 0.55.

Remark 5.1
For one-dimensional nonlinear hyperbolic scalar equations and systems, results using the Lax-Friedrichs flux can be found in [16,17]. The theoretical extension to multi-dimensional nonlinear hyperbolic equations is quite difficult to obtain as they rely on divided difference estimates.

Conclusions
In this article, we have discussed the effect of the choice of the flux on three types of superconvergence: pointwise, wave propagation properties, and the negative-order norm. We specifically provide an explict formula for the form of the error constant for the dispersion and dissipation errors. We illustrated this effect through the use of the  approximations. We also proved that the superconvergent extraction capabilities of the SIAC filter are uneffected. Numerical results were presented to confirm our findings.