A square entropy stable flux limiter for P_NP_M schemes
aa r X i v : . [ m a t h . NA ] D ec A square entropy stable flux limiterfor P N P M schemes Claus R. Goetz and Michael Dumbser We study some theoretical aspects of P N P M schemes, which are a novelclass of high order accurate reconstruction based discontinuous Galerkin(DG) schemes for hyperbolic conservation laws. The P N P M schemes storeand evolve the discrete solution u h under the form of piecewise polynomialsof degree N , while piecewise polynomials w h of degree M ≥ N are used forthe computation of the volume and boundary fluxes. The piecewise poly-nomials w h are obtained from u h via a suitable reconstruction or recoveryoperator. The P N P M approach contains high order finite volume methods( N = 0 ) as well as classical DG schemes ( N = M ) as special cases of amore general framework. Furthermore, for N = M and N > , it leads to anew intermediate class of methods, which can be denoted either as Hermitefinite volume or as reconstructed DG methods. We show analytically why P N P M methods for N = M are, in general, not L -diminishing. To thisend, we extend the well-known cell entropy inequality and the following L stability result of Jiang and Shu for DG methods (i.e. for N = M ) to thegeneral P N P M case and identify which part in the reconstruction step maycause the instability. With this insight we design a flux limiter that en-forces a cell square entropy inequality and thus an L stability condition for P N P M schemes for scalar conservation laws in one space dimension. Fur-thermore, in this paper we prove existence and uniqueness of the solutionof the P N P M reconstruction operator.
1. Introduction
Two of the most popular families of high order schemes for hyperbolic conservation lawsare either reconstruction-based finite volume schemes, or discontinuous Galerkin (DG) methods that directly evolve the expansion coefficients of a piecewise polynomial datarepresentation, and thus in principle do not need any reconstruction stage. The P N P M philosophy introduced in [8, 7] provides a unified framework for the treatment of bothapproaches. At this point it is important to stress that reconstruction operators havealready been used in connection with the DG method before, but only in the contextof accuracy-enhancing postprocessors , which are applied to the discrete solution onlyonce at the end of the simulation, see [24, 5], or as nonlinear moment limiters , in Corresponding author,
University of Hamburg, Department of Mathematics, Bundesstr.55, D-20146 Hamburg, Germany, [email protected] University of Trento, Department of Civil, Environmental and Mechanical Engi-neering, Via Mesiano 77, I-38123 Trento, Italy, [email protected] , P N P M framework [8, 7]. Similar reconstruction-based DG schemes have subsequently beendeveloped and applied also for example in [19, 20, 29].In the context of finite volume schemes, high order of accuracy in space is achievedby reconstructing a high order piecewise polynomial representation of the solution(say, an ENO or WENO polynomial) from the known cell averages [16, 18, 3, 26].An example for a high order fully discrete one-step method following this approach isthe ADER scheme that is based on the approximate solution of generalized Riemannproblems arising from the piecewise polynomial reconstruction of the data, see e.g.[25, 27, 28, 1, 4, 14] for more details. Methods that only evolve cell averages in timehave the advantage that they can use large time-steps (CFL < N in space. At very high order, however, thisrequires very small time steps (CFL < / (2 N + 1) ). A natural approach for reachinghigher accuracy in the DG framework is then to also include a suitable reconstructionprocedure in the scheme. Since more degrees of freedom are available in each cell, wecan work with compact reconstruction stencils and still keep a larger CFL number, see[8]. This is the basic idea of the P N P M method: The solution is represented in a finiteelement space of piecewise polynomials of degree N ≥ (hence the P N in the nameof the method) and at each time step before the time evolution is carried out, a highorder reconstruction of piecewise polynomials of degree M ≥ N is computed (hencethe P M ). In this framework, the pure DG method can be viewed as a P N P N scheme,while the case P P M corresponds to the high order finite volume schemes. For N > and M > N a family of hybrid schemes emerges.The P N P M method has already been successfully applied to a variety of complex flowproblems, including viscous terms, non-conservative products and stiff source terms,see [8, 7, 9, 11], but our understanding of its analytical properties has not reached amature level yet. In particular, it was observed that even in the D case for scalar,linear problems, the method is in general not strictly L -diminishing [12, 6] and thefamous cell square entropy inequality of Jiang and Shu [17] for DG schemes is not validfor M > N .In this paper, we extend the technique of Jiang and Shu to the case
M > N andderive a modified sufficient criterion for entropy stability that reduces to the originalcriterion of Jiang and Shu for M = N . We demonstrate why this criterion is in generalnot satisfied for problems with large jumps in the solution. In the linear case the newcondition is a simple algebraic relation between the jump in the reconstruction andthe jump in the data. Those jumps can be expected to be of similar size for smoothsolutions, but can differ substantially in the vicinity of strong shocks. In order tostabilize the method, we propose a new flux limiter that strictly enforces the cellentropy inequality and thus L stability of general P N P M schemes.The rest of this paper is organized as follows: At first we describe the P N P M method2n Section 2. After that we discuss the reconstruction procedure in Section 3. Insection 4, we derive a condition on the reconstruction operator in order to satisfya cell entropy inequality. Our approach is analogous to the technique of Jiang andShu and we discuss which steps in their proof have to be modified when we include areconstruction operator in the scheme. We use the new insight where entropy stabilitymay be violated to develop a new flux limiter that ensures entropy stability in Section5. The new limiter is tested and validated in several numerical examples, which arepresented in Section 6. We summarize the work and draw conclusion in Section 7.Finally, in the appendix, we provide a proof that the P N P M reconstruction problemon the three cell central stencil has a unique solution for M = 3 N + 2 , which is themaximal possible M that can be reached in this case. It seems that such a result forarbitrary N has not been available in the literature so far.
2. The P N P M method In this section we describe the fundamental concepts of the P N P M method. Our focuslies on the analytical aspects of the scheme, so questions of practical implementationwill not be addressed in detail. For those we refer to [8]. We consider the followingframework:We solve the initial value problem for a scalar conservation law in one space dimen-sion, ∂∂t u ( x, t ) + ∂∂x f ( u ( x, t )) = 0 , ( x, t ) ∈ Ω × [0 , ∞ ) , (1) Ω ⊆ R , u ( x, t ) ∈ U ⊂ R , f : U → R , with the space of admissible states U , initial data u ( x,
0) = ¯ u ( x ) and, if necessary, suitable boundary conditions. Let T = { T ( i ) , i ∈ I } for some indexset I be a partition of Ω . In the following we assume for simplicity that we have auniform grid on the whole real line, Ω = R , T ( i ) = [ x i − , x i + ] , with x i + = (cid:18) i + 12 (cid:19) h, i ∈ Z , h > . Let u h ( x, t ) ∈ V h be a piecewise polynomial representation of the solution at time t ,such that in each cell T ( i ) = [ x i − , x i + ] the function u ( i ) = u h | T ( i ) is a polynomial ofdegree N . At time t = 0 , u h ( x, is obtained as the L projection of the initial data u ( x, onto the space V h of piecewise polynomials of degree N on T .We denote one-sided limits at the cell interfaces by u + i − ( t ) = lim x → x i − x>x i − u h ( x, t ) , u − i + ( t ) = lim x → x i + 12 x 3. Reconstruction Our goal in this section is to construct a function w h , such that for each cell T ( i ) the function w ( i ) = w h | T ( i ) is a polynomial of degree M > N . To this end, we firstchoose a reconstruction stencil S ( i ) ⊂ T for each T ( i ) and construct a polynomial w ( i ) s of degree M on the whole stencil S ( i ) . The function w ( i ) = w ( i ) s | T ( i ) is then givenas the restriction of that reconstruction polynomial w ( i ) s to the cell T ( i ) , while w h isgiven by the union of all the w ( i ) for all T ( i ) ∈ T . For simplicity, we only discuss thereconstruction procedure on the three cell central stencil S ( i ) = { T ( i − , T ( i ) , T ( i +1) } .Here, by a slight abuse of notation, we also write S ( i ) = T ( i − ∪ T ( i ) ∪ T ( i +1) to denotethe union of all cells in the stencil. Denote the space of piecewise polynomials of degree N on S ( i ) by V Nh (cid:16) S ( i ) (cid:17) = n p : S ( i ) → R : p | T ( i ) ∈ P N o , and the space of polynomials of degree M > N on S ( i ) by W Mh (cid:16) S ( i ) (cid:17) = n p : S ( i ) → R : p ∈ P M o . Note that while in the end we are looking for a reconstruction in the space of piecewise polynomials of degree M (i.e. a function in W h = V Mh ( T ) ), for the moment we considerthe space W M (cid:0) S ( i ) (cid:1) of polynomials that are continuously differentiable on the wholestencil S ( i ) .For each cell T ( j ) ∈ S ( i ) , let { Φ ( j ) ℓ , ℓ ∈ { , . . . , N }} be the shifted Legendre polyno-mials on T ( i ) , which form a basis of the polynomial space P N (cid:0) T ( j ) (cid:1) , that is orthogonalin the local L -sense: D Φ ( j ) m , Φ ( j ) n E T ( j ) = Z T ( j ) Φ ( j ) m ( x )Φ ( j ) n ( x ) dx = 0 , if n = m. For W Mh (cid:0) S ( i ) (cid:1) we choose a basis { Ψ k , k ∈ { , . . . , M }} such that h Ψ m , Ψ n i T ( i ) = 0 for n = m. Note that while all Ψ k are defined on the whole stencil S ( i ) , we require only theirrestrictions to the central cell T ( i ) to be orthogonal. To achieve this, we simply takethe Ψ k to be Legendre polynomials on the central cell T ( i ) and extend them to thewhole stencil. Note that this in particular means that Ψ k | T ( i ) = Φ ( i ) k , for k = 0 , . . . , N. 5n the following, we denote by w ( i ) s the continuous extension of the polynomial w ( i ) of degree M to the entire stencil S ( i ) , while w ( i ) s (cid:12)(cid:12)(cid:12) T ( i ) = w ( i ) . According to [8], werequire the reconstruction w ( i ) s ∈ W M (cid:0) S ( i ) (cid:1) from a given u h ∈ V Nh (cid:0) S ( i ) (cid:1) , to satisfy a generalized conservation principle or identity in the weak sense , i.e. we require (2) tohold for all cells in S ( i ) . That means D w ( i ) s , Φ ( j ) ℓ E T ( j ) = D u h , Φ ( j ) ℓ E T ( j ) for each ℓ ∈ { , . . . , N } , j ∈ { i − , i, i + 1 } , With w ( i ) s = M X k =0 ˆ w ( i ) k Ψ k , u ( j ) = N X k =0 ˆ u ( j ) k Φ ( j ) k , for each ℓ ∈ { , . . . , N } and each j ∈ { i − , i, i + 1 } this leads to the condition M X k =0 ˆ w ( i ) k h Ψ k , Φ ( j ) ℓ i T ( j ) = N X k =0 ˆ u ( j ) k h Φ ( j ) k , Φ ( j ) ℓ i T ( j ) . (5)Let’s define a basis of V Nh (cid:0) S ( i ) (cid:1) in the following way: At first, we extend each ofthe local basis functions Φ ( j ) ℓ to the whole stencil S ( i ) by setting ˜Φ ( j ) ℓ ( x ) = (cid:26) Φ ( j ) ℓ ( x ) , x ∈ T ( j ) , , else . Next, let k be an index k ≡ k ( ℓ, j ) , such that k (0 , i − 1) = 0 , . . . , k ( N, i − 1) = N,k (0 , i ) = N + 1 , . . . , k ( N, i ) = 2 N + 1 ,k (0 , i + 1) = 2 N + 2 , . . . , k ( N, i + 1) = 3 N + 2 . Finally, we set Θ k = ˜Φ ( j ) ℓ , ℓ ∈ { , . . . , N } , j ∈ { i − , i, i + 1 } , k = k ( ℓ, j ) . Then the Θ k form a basis of V Nh (cid:0) S ( i ) (cid:1) , that is orthogonal in the sense that h Θ m , Θ n i S ( i ) = Z S ( i ) Θ m ( x )Θ n ( x ) dx = 0 , if m = n. With this, we can rewrite (5) as M X ℓ =0 ˆ w ( i ) ℓ h Ψ ℓ , Θ k i S ( i ) = N +2 X m =0 ˆ u m h Θ m , Θ k i S ( i ) , k ∈ { , . . . , N + 2 } . (6)If M < N + 2 , this is an overdetermined system, which we solve at the aid of aconstrained least-squares method, see [10, 8]. In Appendix A we give an existence and6niqueness proof for M = 3 N + 2 , which, to our knowledge, has not been available inthe literature so far.In particular, by requiring conservation on the central cell of the stencil, we get ˆ w ( i ) ℓ = ˆ u ( i ) ℓ , for ℓ = 0 , . . . , N, which in turn implies that the residual r h = w h − u h , is orthogonal to V h and thus to u h on T ( i ) : Z T ( i ) r h ( x ) u h ( x ) dx = Z T ( i ) M X k = N +1 ˆ w ( i ) k Ψ k ( x ) ! N X k =0 ˆ u ( i ) k Φ ( i ) k ( x ) ! dx = 0 . 4. The cell entropy condition The focus of our analysis in this paper is the entropy stability (or lack thereof) of P N P M schemes. Recall that a convex function Q : U → R is called an entropy for theconservation law (1), if there exists a so-called entropy flux F with F ′ ( u ) = Q ′ ( u ) f ′ ( u ) .We say that a weak solution u of (1) satisfies the entropy inequality , if ∂∂t Q ( u ) + ∂∂x F ( u ) ≤ (7)in a distributional sense. For classical (i.e. C ) solutions, (7) is satisfied with equality.A discrete version of (7) is given for each cell T ( i ) ∈ T by Z T ( i ) ∂∂t Q ( u h ( x, t )) dx + F i + ( t ) − F i − ( t ) ≤ , (8)where F i + is a numerical entropy flux, consistent with F . For the special case of the square entropy Q ( u ) = u / , by summing over all T ( i ) ∈ T , condition (8) leads to the L -stability of the scheme: ddt k u h ( x, t ) k L (Ω) ≤ . Although we will mostly talk about entropy stability in this paper, it is in fact the L -stability of the scheme that we are most interested in.Jiang and Shu [17] have shown that the pure discontinuous Galerkin scheme, whichwe interpret as a P N P N scheme with M = N , i.e. with w h = u h in (4), satisfiesthe discrete cell entropy inequality (7) for the square entropy. This does, however,not hold in general for w h = u h . To our knowledge, this phenomenon has so far notyet been explained theoretically. In this section we follow the original constructionof Jiang and Shu and highlight which steps need to be modified to cover the general P N P M case. We show why the reconstruction procedure can lead to a violation of theentropy condition. 7e begin analogously to Jiang and Shu, i.e. in (4) take ϕ = u h ( x, t ) , so that onegets ddt Z T ( i ) ( u h ( x, t )) dx + f wi + ( t ) u − i + ( t ) − f wi − ( t ) u + i − ( t ) − Z T ( i ) f ( w h ( x, t )) ∂u h ∂x dx = 0 . (9)Let us denote E ( i ) ( t ) = Z T ( i ) ( u h ( x, t )) dx. Recall that in the P N P M method, contrary to the pure DG method, numerical fluxesand the physical flux in the volume integral are evaluated at the reconstructed poly-nomial w h , rather than at the DG polynomial u h . In general u h ∈ V h and w h ∈ W h are not from the same space. However, we formally interpret f wi + as a four-point flux f R ( u ) i + : f R ( u ) i + ( t ) = f R ( u ) i + (cid:16) u ( i − , . . . , u ( i +2) (cid:17) ( t )= ¯ f i + (cid:18) R (cid:16) u ( i − , u ( i ) , u ( i +1) (cid:17) − i + ( t ) , R (cid:16) u ( i ) , u ( i +1) , u ( i +2) (cid:17) + i + ( t ) (cid:19) . So we can formally do all analysis in the space of polynomials of degree N . Theconsistency condition f R ( u ) i + ( u, . . . , u ) = f ( u ) is satisfied if the reconstruction operator preserves constants, i.e. R ( u, . . . , u ) = u. The next step is to construct a numerical entropy flux. Recall that an entropy flux F for the square entropy is given by F ( u ) = f ( u ) u − g ( u ) where g is a primitive of f : g ( u ) = Z f ( u ) du. To simplify the notation, in the following let us drop the time-dependence in thenotation, but keep in mind that we are describing a semi-discrete method.We are looking for a numerical entropy flux F i + that is consistent with F . Take g as defined above and write f ( w ) = f ( u ) + ( f ( w ) − f ( u )) . Then we get from (9): ddt E ( i ) + f wi + u − i + − f wi − u + i − + g ( u + i − ) − g ( u − i + ) − Z T ( i ) ( f ( w h ) − f ( u h )) ∂u h ∂x dx = 0 . (10)8dding and subtracting f wi − u − i − − g ( u − i − ) from (10) and rearranging terms yields ddt E ( i ) + f wi + u − i + − g ( u − i + ) − f wi − u − i − + g ( u − i − ) − Z T ( i ) ( f ( w h ) − f ( u h )) ∂u h ∂x dx, − f wi − u + i − + f wi − u − i − − g ( u − i − ) + g ( u + i − ) = 0 . (11)Setting ˜ F i + = f wi + u − i + − g ( u − i + ) and V ( i ) = Z T ( i ) ( f ( w h ) − f ( u h )) ∂u h ∂x dxB ( i ) = − f wi − (cid:16) u + i − − u − i − (cid:17) + g ( u + i − ) − g ( u − i − ) − V ( i ) , eqn. (11) becomes ddt E ( i ) + ˜ F i + − ˜ F i − + B ( i ) = 0 . (12)The numerical entropy flux ˜ F i + is consistent with the entropy flux F . Therefore, if B ( i ) ≥ , we get the discrete entropy inequality ddt E ( i ) + ˜ F i + − ˜ F i − ≤ . (13)The rest of this section is devoted to analysing under which conditions we can guarantee B ( i ) ≥ and thereby the entropy inequality (7).By the mean-value theorem, there exists a ξ between u − i − and u + i − , such that g ( u + i − ) − g ( u − i − ) = Z u + i − u − i − f ( u ) du = (cid:16) u + i − − u − i − (cid:17) f ( ξ ) , and so a sufficient condition for B ( i ) ≥ is − (cid:16) u + i − − u − i − (cid:17) (cid:16) f wi − − f ( ξ ) (cid:17) − V ( i ) ≥ . (14)Note that up to here there our analysis is completely analogous to the Jiang andShu proof [17]. The only differences are that in our case f wi − is evaluated using thereconstructed solution w h and that we have to include the term V ( i ) .Assume that u − i − > u + i − , then condition (14) becomes ¯ f i − (cid:16) w − i − , w + i − (cid:17) ≥ u − i − − u + i − Z u − i − u + i − f ( u ) du − V ( i ) . (15)9bviously, this condition depends on the flux f , the numerical flux ¯ f i + and thereconstruction w h = R ( u h ) . In the following, we discuss the restrictions on w h imposedby (15) for given f and ¯ f i + .To better understand where the entropy condition may not be satisfied, we start withthe easiest case possible: We discuss the linear advection equation, f ( u ) = λu, λ > constant and take ¯ f i + to be the upwind flux. This allows us to find a simple algebraiccondition that relates the jump in the reconstruction to the jump in the data. It turnsout that already in this simple case, the reconstruction procedure can lead to a violationof the entropy condition.Assume that at the cell-interface i − we have a jump in the data with u − i − > u + i − .With the upwind flux, condition (15) reads λw − i − ≥ u − i − − u + i − Z u − i − u + i − f ( u ) du − V ( i ) . In the linear case, the flux-average f ( ξ ) in the jump is simply f ( ξ ) = 1 u + i − − u − i − Z u + i − u − i − λu du = 12 λ (cid:16) u − i − + u + i − (cid:17) , and the volume integral of the flux difference can be written in terms of the residual r h = w h − u h as V ( i ) = Z T ( i ) λ ( w h − u h ) ∂u h ∂x dx = Z T ( i ) λ r h ∂u h ∂x dx = 0 . (16)Since inside each cell T ( i ) the function ∂u h ∂x is a polynomial of degree k − and thereforecan be expressed as a linear combination of Φ ( i )0 , . . . , Φ ( i ) k − , i.e. we have ∂u h ∂x ∈ V h andthus the integral (16) vanishes according to the orthogonality of r h with respect to V h ,see also (3). In summary, for linear problems (14) reduces to λw − i − − f ( ξ ) ≥ . Thus, we get the pointwise condition w − i − ≥ (cid:16) u − i − + u + i − (cid:17) . (17)Analogously, if u − i − < u + i − , we get w − i − ≤ (cid:16) u − i − + u + i − (cid:17) . (18)It is easy to construct counterexamples, in which (17), (18) is not satisfied. However,numerical experiments suggest, that this occurs in practice only when there is a largejump in u h . 10 x u , w Figure 1: A P P reconstruction that does not satisfy condition (17), (18).Figure 1 shows an example where condition (17), (18) is not satisfied. We havepiecewise linear data (solid black lines) and compute two polynomials of degree five(red and blue). The leftmost there parts of the initial data are used to compute apolynomial w L (plotted in red), the parts to the right are used for the polynomial w R (plotted in blue). For the intercell flux between the cell ( − , and (1 , we considerthe interface i − / . Condition (18) states that in order to satisfy the entropyinequality we need w L (1) < ( u +1 + u − ) , which is marked as a magenta point in theplot. Clearly, this condition is not satisfied.In the next section, we describe a way to stabilize the method using a flux limiter approach. Our ansatz also covers the nonlinear case. 5. Flux limiting Now that we have seen where instabilities may occur, we develop a nonlinear fixfor that problem. In order to enforce a cell entropy condition, we employ a fluxlimiting approach. Denote by f ui + , f wi + numerical fluxes evaluated using u h and w h ,respectively. Our goal is to find a limiter θ i + ∈ [0 , such that the scheme with thelimited flux f i + = f ui + + θ i + f ri + , f ri + = f wi + − f ui + (19)11atisfies the entropy condition. We denote the terms that only include u h by F ui + = f ui + u − i + − g ( u − i + ) ,A ( i ) = − (cid:16) u + i − − u − i − (cid:17) f ui − + g ( u + i − ) − g ( u − i − ) . Note that these are exactly the terms that occur in the original work of Jiang and Shu[17] and we recall that the essential point of their proof is that A ( i ) ≥ since f ui − isa monotone E-flux. Moreover, let us use the standard notation for a jump at the cellinterface [[ u ]] i − = u + i − − u − i − . Inserting (19) in (12) and rearranging terms yields: dE ( i ) dt + (cid:16) F ui + + θ i + f ri + u − i + (cid:17) − (cid:16) F ui − + θ i − f ri − u − i − (cid:17) + A ( i ) − V ( i ) − θ i − [[ u ]] i − f ri − = 0 (20)Thus, defining a numerical entropy flux as ˆ F i + = F ui + + θ i + f ri + u − i + , the entropy condition dE ( i ) dt + ˆ F i + − ˆ F i − ≤ is satisfied if A ( i ) − V ( i ) − θ i − [[ u ]] i − f ri − ≥ . (21)In the linear case f ( u ) = λu we have V ( i ) = 0 due to the orthogonality of the residual r h w.r.t. V h , hence the condition (21) simplifies to . A ( i ) − θ i − [[ u ]] i − f ri − ≥ . (22)Note that A ( i ) ≥ and so in the linear case we can always find a θ i − ∈ [0 , thatsatisfies (22).In the nonlinear case, we have to account for the volume integral V ( i ) of the fluxdifference in (21). Since the limiter only acts on the cell-boundary , it is possible toinclude this new term, which only takes information from inside the cell . However, itmay happen that condition (21) cannot be satisfied with a θ i − ∈ [0 , . In this case,we found that good numerical results can be obtained by setting σ i − = A ( i ) − V ( i ) [[ u ]] i − f ri − (avoiding divisions by zero) and then compute the actual flux limiter from θ i − = max (cid:16) min( σ i − , , (cid:17) . 12n the update for the degrees of freedom inside each cell, we do not limit the poly-nomial w h inside the cell if θ i − ∈ (0 , . However, for θ i − = 0 , it can happenthat A ( i ) − V ( i ) < and in this case the flux limiter alone does not guarantee theentropy condition. If this occurs, additional limiting of the polynomial inside the cellis necessary. Numerical experiments suggest that in the case of entropy violation for θ i − = 0 , one can simply set w h | T ( i ) := u ( i ) + θ i ( w ( i ) − u ( i ) ) with θ i ∈ [0 , , so that A ( i ) + V ( i ) ≥ and thus entropy stability is guaranteed. Numerical evidence showsthat even setting θ i = 0 does not hurt the accuracy too much. Just as in the pure DGcase, it is important to stress that even with the limiter developed above, we can onlyenforce the entropy inequality and thereby the L -stability of the scheme. It does nottake into account requirements on the L ∞ norm. Moreover, the reconstruction methodwe use is linear. Due to Godunov’s theorem [13], the high order P N P M schemes withthis reconstruction cannot be monotone. Additional limiting or a different, nonlin-ear WENO/HWENO reconstruction procedure [22, 23, 2] are necessary to deal withspurious oscillations at shock waves. 6. Numerical results At first we test the effects of the flux limiter for a linear problem with smooth initialdata. In this case, even the unlimited P N P M scheme has little to no problems with L stability and so the main question will be whether the flux limiting affects the orderof accuracy. Consider the problem ∂u∂t + ∂u∂x = 0 , x ∈ ( − , with initial data u ( x, 0) = sin( πx ) and periodic boundary conditions. Errors are measured at t = 1 . Time integration isperformed with an ( M + 1) -stage linear Runge-Kutta method of order M + 1 [15] andthe Rusanov (local Lax-Friedrichs) numerical flux is used.We observe that on coarse grids the limited version of the P N P M scheme is only oforder N + 1 . On finer grids, however, the limiter is never active and the full order M +1 of the unlimited scheme is achieved. At which level of grid refinement the limiterbecomes inactive depends on N : For larger N , no limiting is needed on rather coarsegrids, e.g. for P P and P P the limiter is inactive on a grid with cells, whichcorresponds to a grid size of h = 1 / . Generally speaking, if the representation of thedata in the lower order polynomials is already sufficiently good, then the high orderrepresentation does not need to be limited.13 = 1 N = 2 Limiter off Limiter on Limiter off Limiter on M = 2 I E O E O E O E O 10 2 . E − 01 2 . E − . E − 02 2 . 67 6 . E − 02 1 . . E − 03 3 . 27 1 . E − 02 2 . . E − 04 2 . 91 5 . E − 03 1 . . E − 05 2 . 64 1 . E − 03 2 . M = 3 10 2 . E − 01 2 . E − 01 1 . E − 02 2 . E − . E − 02 3 . 19 5 . E − 02 2 . 14 6 . E − 04 4 . 79 4 . E − 03 2 . . E − 03 4 . 16 9 . E − 03 2 . 49 3 . E − 05 4 . 04 4 . E − 04 3 . . E − 05 4 . 05 2 . E − 03 1 . 65 2 . E − 06 3 . 67 3 . E − 05 3 . . E − 06 3 . 75 9 . E − 04 1 . 59 2 . E − 07 3 . 55 2 . E − 07 7 . M = 4 10 6 . E − 02 1 . E − 01 1 . E − 02 2 . E − . E − 03 5 . 15 5 . E − 02 1 . 84 5 . E − 04 4 . 98 3 . E − 03 2 . . E − 05 4 . 63 1 . E − 02 1 . 91 2 . E − 05 4 . 58 4 . E − 04 3 . . E − 06 3 . 78 3 . E − 03 2 . 06 9 . E − 07 4 . 49 3 . E − 05 3 . . E − 07 3 . 54 6 . E − 04 2 . 31 4 . E − 08 4 . 49 4 . E − 08 9 . M = 5 10 4 . E − 02 1 . E − 01 3 . E − 03 1 . E − . E − 04 6 . 55 5 . E − 02 1 . 76 3 . E − 05 6 . 49 4 . E − 03 1 . . E − 05 3 . 27 1 . E − 02 1 . 80 7 . E − 07 5 . 72 4 . E − 04 3 . . E − 06 3 . 15 4 . E − 03 1 . 87 1 . E − 08 5 . 54 3 . E − 05 3 . . E − 07 3 . 44 9 . E − 04 2 . 16 3 . E − 10 5 . 51 3 . E − 10 16 . M = 6 10 2 . E − 03 1 . E − . E − 05 7 . 07 4 . E − 03 2 . . E − 07 6 . 23 4 . E − 04 3 . . E − 09 5 . 89 3 . E − 05 3 . . E − 11 5 . 64 9 . E − 11 18 . N = 3 N = 4 Limiter off Limiter on Limiter off Limiter on M = 4 I E O E O E O E O 10 1 . E − 03 1 . E − . E − 05 5 . 13 1 . E − 04 2 . . E − 06 4 . 89 1 . E − 06 7 . . E − 08 4 . 67 4 . E − 08 4 . . E − 09 4 . 55 1 . E − 09 4 . M = 5 10 1 . E − 03 1 . E − 03 9 . E − 05 1 . E − . E − 05 5 . 38 1 . E − 04 2 . 95 1 . E − 06 6 . 05 1 . E − 06 6 . . E − 06 5 . 43 5 . E − 07 8 . 32 2 . E − 08 5 . 86 2 . E − 08 5 . . E − 08 5 . 49 1 . E − 08 5 . 49 5 . E − 10 5 . 63 5 . E − 10 5 . . E − 10 5 . 50 2 . E − 10 5 . 50 1 . E − 11 5 . 53 1 . E − 11 5 . M = 6 10 1 . E − 04 1 . E − 03 8 . E − 05 1 . E − . E − 06 6 . 09 1 . E − 04 3 . 13 1 . E − 06 6 . 33 1 . E − 06 7 . . E − 08 6 . 59 2 . E − 08 12 . 78 1 . E − 08 6 . 45 1 . E − 08 6 . . E − 10 7 . 46 1 . E − 10 7 . 46 1 . E − 10 6 . 48 1 . E − 10 6 . . E − 12 5 . 66 3 . E − 12 5 . 66 1 . E − 12 6 . 45 1 . E − 12 6 . Table 1: L errors for the linear advection equation with smooth initial data14 .2. Burgers equation Consider Burgers equation ∂u∂t + ∂∂x (cid:18) u (cid:19) = 0 , (23)with initial data u ( x, 0) = − − (cid:18) x − (cid:19) ! + 5 exp − (cid:18) x + 12 (cid:19) ! . (24)We solve (23), (24) on ( − , with transmissive boundary conditions. We use a fourthorder SSP Runge-Kutta method and denote by θ the average value of the limiter duringthe Runge-Kutta stages for one time-step. We present the results for P P on cellsand for P P at several time points during the simulation.Results are shown in Figure 2. We observe a similar behaviour as in the linearcase: In regions where the solution is smooth, the limiter is not active. In fact, mostlimiting is needed shortly before the shocks form. Again we can see that if the P N representation of the solution is sufficiently accurate, even on rather coarse grid the P M part requires no or only very little limiting. Consider the following Lighthill-Whitman type model for traffic flow: ∂ρ∂t + ∂∂x (cid:18) ρ exp (cid:18) − ρ (cid:19)(cid:19) = 0 , (25)where ρ ∈ (0 , describes the density of cars on a road. Note that in this case f is a strictly concave flux. We solve (25) on ( − , with periodic boundary conditions andinitial data given by ρ ( x, 0) = 12 + 14 sin( πx ) . The solution is given by a sinusodial wave travelling to the right that is deformed untila shock emerges. Note that in this model, lower values of ρ lead to a faster speed ofpropagation.Numerical results are shown in Figure 3. The left column depicts the numericalsolution on cells and the left column shows the results for P P on cells. Again,time integration is performed with a fourth order SSP Runge-Kutta method and θ denotes the average value of the limiter during the respective time-step. We observethat the limiter is only active on a very small number of cells, shortly before the shockis formed and close to the shock front after its formation, while the numerical solutionis not limited as long as the exact solution is smooth.15 x -505 u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ -1 -0.5 0 0.5 1 x -505 u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ -1 -0.5 0 0.5 1 x -50510 u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ -1 -0.5 0 0.5 1 x -505 u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ -1 -0.5 0 0.5 1 x -505 u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ -1 -0.5 0 0.5 1 x -505 u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ c) t = 0 . Figure 2: Burgers equation, left: P P solution on cells, right: P P solution on cells. Top row: t = 0 . , middle row: t = 0 . , bottom row: t = 0 . x u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ -1 -0.5 0 0.5 1 x u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ -1 -0.5 0 0.5 1 x u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ -1 -0.5 0 0.5 1 x u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ -1 -0.5 0 0.5 1 x u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ -1 -0.5 0 0.5 1 x u uLimiter onLimiter off -1 -0.5 0 0.5 1 x θ θ c) t = 0 . Figure 3: Traffic flow, left: P P solution on cells, right: P P solution on cells.Top row: t = 0 . , middle row: t = 0 . , bottom row: t = 0 . . Conclusion We have seen that the reconstruction step in P N P M schemes can cause unstable be-haviour with respect to the L norm of the numerical solution. We demonstratedanalytically how a numerical flux using reconstructed function values can lead to aviolation of a cell entropy condition. In particular for linear problems we were ableto derive a simple algebraic relation between the jump in the reconstruction and thejump in the data that determines whether the scheme is square entropy stable, or not.With this new condition, that also has a nonlinear analogon, we were able to con-struct a flux limiter that guarantees the entropy stability of the numerical solution.As expected, entropy stability only becomes an issue when the solution contains largejumps or very steep gradients, in which case the reconstruction can become oscillatory.Numerical experiments validate the limiter’s ability to stabilize the method while stillmaintaining high order accuracy in smooth regions for sufficiently fine meshes. A. Existence and uniqueness of the reconstruction inthe case M = 3 N + 2 . We show the existence and uniqueness of the solution to the reconstruction problem N +2 X m =0 ˆ w m h Ψ k , Φ ( j ) ℓ i T ( j ) = N X m =0 ˆ v ( j ) k h Φ ( j ) k , Φ ( j ) ℓ i T ( j ) , (26)for ℓ = 0 , . . . , N, j ∈ { i − , i, i +1 } . Recall that the functions { Φ ( j ) ℓ : ℓ = 0 , . . . , N } arethe shifted Legendre polynomials of degree ℓ on T ( j ) and the { Ψ m : m = 0 , . . . N +2 } are Legendre polynomials on the central cell T ( i ) , extended to the whole stencil.A basis { Θ k : k = 0 , . . . , N + 2 } of V Nh (cid:0) S ( i ) (cid:1) , that is orthogonal in the sense that h Θ m , Θ n i S ( i ) = Z S ( i ) Θ m ( x )Θ n ( x ) dx = 0 , if m = n. is given by ˜Φ ( j ) ℓ ( x ) = ( Φ ( j ) ℓ ( x ) , x ∈ T ( j ) , , else . and an an index k ≡ k ( ℓ, j ) , such that k (0 , i − 1) = 0 , . . . , k ( N, i − 1) = N,k (0 , i ) = N + 1 , . . . , k ( N, i ) = 2 N + 1 ,k (0 , i + 1) = 2 N + 2 , . . . , k ( N, i + 1) = 3 N + 2 . We let Θ k = ˜Φ ( j ) ℓ , ℓ ∈ { , . . . , N } , j ∈ { i − , i, i + 1 } , k = k ( ℓ, j ) . N +2 X k =0 ˆ w k h Ψ k , Θ m i S ( i ) = N +2 X k =0 ˆ v k h Θ k , Θ m i S ( i ) , m = 0 , . . . , N + 2 , (27)where we denote ˆ v k ≡ ˆ v k ( ℓ,j ) = ˆ v ( j ) ℓ . It is worth pointing out that condition (27)implies that v is the L -projection of w onto V Nh (cid:0) S ( i ) (cid:1) . Thus, our reconstructionproblem is an inverse projection problem, for which existence and uniqueness of thesolution are non-trivial.In order to guarantee existence and uniqueness for the solution of the reconstructionproblem, we have to check whether the matrix ˜ A = ( h Θ m , Ψ ℓ i S ( i ) ) m,ℓ ∈{ ,..., N +2 } is invertible.At first, note that the conditions for the reconstruction are invariant under linearcoordinate-transformations. So we can map T ( i − to the interval ( − , − , T ( i ) to ( − , , and T ( i +1) to (1 , . Denote the Legendre polynomial of degree k on ( − , by P k and the the shifted Legendre polynomials of degree k on ( − , and (1 , by P ( − k and P (+1) k , respectively.Then, after the coordinate change, the basis functions Φ ( i − k and Φ ( i +1) k become P ( − k and P (+1) k , respectively, while the Φ ( i ) k simply become P k . Moreover, the Ψ k also become P k . Thus, on ( − , − we get the matrix coefficient Z − − P k ( s ) P ( − ℓ ( s ) ds = Z − P k ( s + 2) P ( − ℓ ( s + 2) ds = Z − P k ( s + 2) P ℓ ( s ) ds, by shifting to ( − , and noting that by construction P ( − ℓ ( s + 2) = P ℓ ( s ) . Similarly,we have Z P ℓ ( s ) P (+1) k ( s ) ds = Z − P ℓ ( s − P k ( s ) ds. By a simple symmetry argument, it is then easy to see that Z − P ℓ ( s + 2) P k ( s ) ds = ( − k + ℓ Z − P ℓ ( s − P k ( s ) ds. Let a k,ℓ = Z − P k ( s ) P ℓ ( s + 2) ds, b k,ℓ = ( − k + ℓ a k,ℓ ˜ A is invertible, if and only if the matrix A = Z − P k ( s ) P ℓ ( s ) ds Z − P k ( s ) P ℓ ( s + 2) ds Z − P k ( s ) P ℓ ( s − ds k = 0 , . . . , Nℓ = 0 , . . . , N + 2 = k + 1 δ k,ℓ a k,ℓ b k,ℓ k = 0 , . . . , Nℓ = 0 , . . . , N + 2 is invertible. We have A = (cid:20) D A A (cid:21) with a diagonal part D ∈ R ( N +1) × ( N +1) and A = a k, . . . a k,N ( − k a k, . . . ( − k + N a k,N k =0 ,...,N ∈ R (2 N +2) × ( N +1) A = a k,N +1 . . . a k, N +2 ( − k + N +1 a k,N +1 . . . ( − k +3 N +2 a k, N +2 k =0 ,...,N ∈ R (2 N +2) × (2 N +2) Defining c k,ℓ = 12 ( a k,ℓ + b k,ℓ ) , it is then sufficient to show that B = a ,N +1 . . . a , N +3 ... ... a N,N +1 . . . a N, N +3 c ,N +1 . . . c , N +3 ... ... c N,N +1 . . . c N, N +3 ∈ R (2 N +2) × (2 N +2) is invertible. Before we show this, let us proof the following:20 emma. Let = Q be polynomial of degree N that does not vanish identically on ( − , , such that for some ℓ > N we have Z − P ℓ ( s + 2) Q ( s ) ds = 0 . Then Z − P ℓ +2 ( s + 2) Q ( s ) ds = 0 . Proof. By Rolle’s theorem and the well-known formula dds P ℓ +1 ( s ) = 2 P ℓ ( s ) k P ℓ k + 2 P ℓ − ( s ) k P ℓ − k + . . . , it is straightforward to check that for s ∈ ( − , and ℓ ≥ the following holds: • P ℓ ( s + 2) > , P ′ ℓ ( s + 2) > , P ′′ ℓ ( s + 2) > , • P ℓ +2 ( s + 2) > P ℓ ( s ) , P ′ ℓ +2 ( s + 2) > P ′ ℓ ( s ) , P ′′ ℓ +2 ( s + 2) > P ℓ ( s + 2) . Then, if for = Q and for some ℓ > N we have Z − P ℓ ( s + 2) Q ( s ) ds = 0 , it follows that Z − P ℓ +2 ( s + 2) Q ( s ) ds = 0 . (cid:3) Theorem. The matrix B is invertible. Therefore, the reconstruction problem has aunique solution. Proof. We proceed by showing that the rows of B are linearly independent. Let a k = [ a k,N +1 , . . . , a k, N +2 ] , c k = [ c k,N +1 , . . . , c k, N +2 ] and assume that λ k , µ k ∈ R , k = 0 , . . . , N , such that N X k =0 ( λ k a k + µ k c k ) = 0 . That means that for each ℓ = N + 1 , . . . , N + 3 we have N X k =0 ( λ k a k,ℓ + µ k c k,ℓ ) = N X k =0 (cid:18) λ k + 1 + ( − k + ℓ µ k (cid:19) a k,ℓ = 0 , Z − P ℓ ( s + 2) ( N X k =0 (cid:18) λ k + 1 + ( − k + ℓ µ k (cid:19) P k ( s ) ) ds = 0 . (28)Let Q ( s ) = N X k =0 λ k P k ( s ) + X k ≤ Nk even µ k P k ( s ) ,Q ( s ) = N X k =0 λ k P k ( s ) + X k ≤ Nk odd µ k P k ( s ) , then condition (28) reads for all even ℓ : Z − P ℓ ( s + 2) Q ( s ) ds = 0 , and for all odd ℓ : Z − P ℓ ( s + 2) Q ( s ) ds = 0 . By the above Lemma, that means that Q = Q = 0 . Therefore, λ k = µ k = 0 for all k and thus the rows of B are linearly independent. (cid:3) Acknowledgments The presented research has been financed by the European Research Council (ERC)under the European Union’s Seventh Framework Programme (FP7/2007-2013) withthe research project STiMulUs, ERC Grant agreement no. 278267. References [1] T. Aboiyar, E.H. Georgoulis, and A. Iske. Adaptive ADER Methods Using Kernel-Based Polyharmonic Spline WENO Reconstruction. SIAM Journal on ScientificComputing , 32:3251–3277, 2010.[2] D. Balsara, C. Altmann, C.D. Munz, and M. Dumbser. A sub-cell based indicatorfor troubled zones in RKDG schemes and a novel class of hybrid RKDG+HWENOschemes. Journal of Computational Physics , 226:586–620, 2007.[3] D.S. Balsera and C.W. Shu. Monotonicity perserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. J. Comput. Phys. ,160:405–452, 2000. 224] C. C. Castro and E. F. Toro. Solvers for the high-order riemann problem forhyperbolic balance laws. Journal of Computational Physics , 227:2481–2513, 2008.[5] B. Cockburn, M. Luskin, C. W. Shu, and E. Suli. Enhanced accuracy by post-processing for finite element methods for hyperbolic equations. Mathematics ofComputation , 72:577–606, 2003.[6] M. Dumbser. Arbitrary High Order Schemes for the Solution of Hyperbolic Con-servation Laws in Complex Domains . Shaker Verlag, Aachen, 2005.[7] M. Dumbser. Arbitrary high order PNPM schemes on unstructured meshes forthe compressible Navier–Stokes equations. Computers & Fluids , 39:60–76, 2010.[8] M. Dumbser, D. S. Balsara, E. F. Toro, and C.-D. Munz. A unified framework forthe construction of one-step finite volume and discontinuous Galerkin schemes onunstructured meshes. Journal of Computational Physics , 227:8209–8253, Septem-ber 2008.[9] M. Dumbser, A. Hidalgo, M. Castro, C. Parés, and E.F. Toro. FORCE schemes onunstructured meshes II: Non–conservative hyperbolic systems. Computer Methodsin Applied Mechanics and Engineering , 199:625–647, 2010.[10] M. Dumbser and M. Käser. Arbitrary high order non-oscillatory finite volumeschemes on unstructured meshes for linear hyperbolic systems. Journal of Com-putational Physics , 221:693–723, 2007.[11] M. Dumbser and O. Zanotti. Very high order PNPM schemes on unstructuredmeshes for the resistive relativistic MHD equations. Journal of ComputationalPhysics , 228:6991–7006, 2009.[12] G. Gassner. Neuartige Discontinuous Galerkin Verfahren für Advektions-Diffusionsgleichungen. Technical report, Institut für Aerodynamik und Gasdy-namik, Universität Stuttgar, 2004.[13] S.K. Godunov. Finite difference methods for the computation of discontinuoussolutions of the equations of fluid dynamics. Mathematics of the USSR - Sbornik ,47:271–306, 1959.[14] C. R. Goetz and A. Iske. Approximate solutions of generalized Riemann problemsfor nonlinear systems of hyperbolic conservation laws. Math. Comp. , 85:35–62,2016.[15] S. Gottlieb, C.-W. Shu, and E. Tadmor. Strong Stability-Preserving High-OrderTime Discretization Methods. SIAM Review , 43:89–112, January 2001.[16] A. Harten, B. Engquist, S. Osher, and S.R. Chakravarthy. Uniformly high or-der accurate essentially non-oscillatory schemes III. Journal of ComputationalPhysics , 71:231–303, 1987. 2317] G. Jiang and C.W. Shu. On a cell entropy inequality for discontinuous Galerkinmethods. Mathematics of Computation , 62:531–538, 1994.[18] G.S. Jiang and C.W. Shu. Efficient implementation of weighted ENO schemes. Journal of Computational Physics , 126:202–228, 1996.[19] H. Luo, L. Luo, R. Nourgaliev, V.A. Mousseau, and N. Dinh. A reconstructeddiscontinuous Galerkin method for the compressible Navier-Stokes equations onarbitrary grids. Journal of Computational Physics , 229:6961–6978, 2010.[20] H. Luo, Y. Xia, S. Spiegel, R. Nourgaliev, and Z. Jiang. A reconstructed discon-tinuous Galerkin method based on a Hierarchical WENO reconstruction for com-pressible flows on tetrahedral grids . Journal of Computational Physics , 236:477–492, 2013.[21] S. Osher. Riemann solvers, the entropy condition and difference approximations. SIAM Journal on Numerical Analysis , 21:217–235, 1984.[22] J. Qiu and C.W. Shu. Hermite WENO schemes and their application as limitersfor Runge-Kutta discontinuous Galerkin method: one-dimensional case. Journalof Computational Physics , 193:115–135, 2003.[23] J. Qiu and C.W. Shu. Hermite WENO schemes and their application as limiters forRunge-Kutta discontinuous Galerkin method II: two dimensional case. Computersand Fluids , 34:642–663, 2005.[24] J.K. Ryan, C.W. Shu, and H.L. Atkins. Extension of a post-processing techniquefor the discontinuous Galerkin method for hyperbolic equations with applicationsto an aeroacoustic problem. SIAM Journal on Scientific Computing , 26:821–843,2005.[25] V.A. Titarev and E.F. Toro. ADER: Arbitrary high order Godunov approach. Journal of Scientific Computing , 17(1-4):609–618, 2002.[26] V.A. Titarev and E.F. Toro. Finite–volume WENO schemes for three–dimensionalconservation laws. Journal of Computational Physics , 201:238–260, 2004.[27] V.A. Titarev and E.F. Toro. ADER schemes for three-dimensional nonlinearhyperbolic systems. Journal of Computational Physics , 204:715–736, 2005.[28] E. F. Toro and V. A. Titarev. Derivative Riemann solvers for systems of conserva-tion laws and ADER methods. Journal of Computational Physics , 212(1):150–165,2006.[29] L. Zhang, L. Wei, H. Lixin, D. Xiaogang, and Z. Hanxin. A class of hybrid DG/FVmethods for conservation laws I: Basic formulation and one-dimensional systems.