Invariant domains and first-order continuous finite element approximation for hyperbolic systems
IINVARIANT DOMAINS AND FIRST-ORDERCONTINUOUS FINITE ELEMENT APPROXIMATIONFOR HYPERBOLIC SYSTEMS ∗ JEAN-LUC GUERMOND † AND
BOJAN POPOV † Abstract.
We propose a numerical method to solve general hyperbolic systems in any spacedimension using forward Euler time stepping and continuous finite elements on non-uniform grids.The properties of the method are based on the introduction of an artificial dissipation that is definedso that any convex invariant sets containing the initial data is an invariant domain for the method.The invariant domain property is proved for any hyperbolic system provided a CFL condition holds.The solution is also shown to satisfy a discrete entropy inequality for every admissible entropy of thesystem. The method is formally first-order accurate in space and can be made high-order in time byusing Strong Stability Preserving algorithms. This technique extends to continuous finite elementsthe work of [Hoff(1979), Hoff(1985)], and [Frid(2001)].
Key words.
Conservation equations, hyperbolic systems, parabolic regularization, invariantdomain, first-order method, finite element method.
AMS subject classifications.
1. Introduction.
The objective of this paper is to investigate a first-order ap-proximation technique for nonlinear hyperbolic systems using continuous finite el-ements and explicit time stepping on non-uniform meshes. Consider the followinghyperbolic system in conservation form(1.1) (cid:40) ∂ t u + ∇· f ( u ) = 0 , for ( x , t ) ∈ R d × R + . u ( x ,
0) = u ( x ) , for x ∈ R d . where the dependent variable u takes values in R m and the flux f takes values in( R m ) d . In this paper u is considered as a column vector u = ( u , . . . , u m ) T . The fluxis a matrix with entries f ij ( u ), 1 ≤ i ≤ m , 1 ≤ j ≤ d and ∇· f is a column vector withentries ( ∇· f ) i = (cid:80) ≤ j ≤ d ∂ x j f ij . For any n = ( n . . . , n d ) T ∈ R d , we denote f ( u ) · n the column vector with entries (cid:80) ≤ l ≤ d n l f il ( u ), where i ∈ { m } . The unit spherein R d centered at 0 is denoted by S d − ( , D the spatial domainwhere the approximation is constructed. The domain D is the d -torus in the case ofperiodic boundary conditions. In the case of the Cauchy problem, D is a compact,polygonal portion of R d large enough so that the domain of influence of u is alwaysincluded in D over the entire duration of the simulation.The method that we propose is explicit in time and uses continuous finite elementson non-uniform grids in any space dimension. The algorithm is described in § ∗ This material is based upon work supported in part by the National Science Foundation grantsDMS-1217262, by the Air Force Office of Scientific Research, USAF, under grant/contract numberFA99550-12-0358, and by the Army Research Office under grant/contract number W911NF-15-1-0517. Draft version, October 15, 2018 † Department of Mathematics, Texas A&M University 3368 TAMU, College Station, TX 77843,USA. 1 a r X i v : . [ m a t h . NA ] S e p J.L. GUERMOND, B. POPOV of the paper are Theorem 4.1 and Theorem 4.2. It is shown in Theorem 4.1 that theproposed scheme preserves all the convex invariant sets as defined in Definition 2.3and it is shown in Theorem 4.2 that the approximate solution satisfies a discreteentropy inequality for every entropy pair of the hyperbolic system. Similar resultshave been established for various finite volumes schemes by [Hoff(1979), Hoff(1985)],[Perthame and Shu(1996)], [Frid(2001)] for the compressible Euler equations and thep-system. Our scheme has no restriction on the nature of the hyperbolic system,besides the speed of propagation being finite. To the best of our knowledge, we arenot aware of any similar scheme in the continuous finite element literature.The paper is organized as follows. The notions of invariant sets and invariantdomains with various examples and other preliminaries are introduced in Section 2.The method is introduces in Section 3. Stability properties of the algorithm areanalyzed in Section 4. Numerical illustrations and comparisons with existing first-order methods are presented in Section 5.
2. Preliminaries.
The objective of this section is to introduce notation andpreliminary results that will be useful in the rest of the paper. We mostly usethe notation and the terminology of [Chueh et al.(1977)Chueh, Conley, and Smoller,Hoff(1979), Hoff(1985), Frid(2001)]. The reader who is familiar with the notions ofinvariant domains and Riemann problems may skip this section and go directly to § We assume that (1.1) is such that there is a clearnotion for the solution of the Riemann problem. That is to say there exists an(nonempty) admissible set
A ⊂ R m such that for any pair of states ( u L , u R ) ∈ A×A and any unit vector n ∈ S d − ( , ∂ t u + ∂ x ( f ( u ) · n ) = 0 , ( x, t ) ∈ R × R + , u ( x,
0) = (cid:40) u L , if x < u R , if x > , has a unique (physical) solution, which we henceforth denote u ( n , u L , u R ).The theory of the Riemann problem for general nonlinear hyperbolic systems withdata far apart is an open problem. Moreover, it is unrealistic to expect a general the-ory for any system with arbitrary initial data. However, when the system is strictlyhyperbolic with smooth flux and all the characteristic fields are either genuinely non-linear or linearly degenerate, it is possible to show that there exists δ > (cid:107) u L − u R (cid:107) (cid:96) ≤ δ , see [Lax(1957)] and [Bressan(2000), Thm 5.3].In particular there are 2 m numbers(2.2) λ − ≤ λ +1 ≤ λ − ≤ λ +2 ≤ . . . ≤ λ − m ≤ λ + m defining up to 2 m + 1 sectors (some could be empty) in the ( x, t ) plane:(2.3) xt ∈ ( −∞ , λ − ) , xt ∈ ( λ − , λ +1 ) , . . . , xt ∈ ( λ − m , λ + m ) , xt ∈ ( λ + m , ∞ ) . The Riemann solution is u L in the sector xt ∈ ( −∞ , λ − ) and u R in the last sector xt ∈ ( λ + m , ∞ ). The solution in the other sectors is either a constant state or anexpansion, see [Bressan(2000), Chap. 5]. The sector λ − t < x < λ + m t , 0 < t , is nvariant domains and C finite element approximation of hyperbolic systems λ max ( n , u L , u R ) := max( | λ − | , | λ + m | )such that for t ≥ u ( x, t ) = (cid:40) u L , if x ≤ − tλ max ( n , u L , u R ) u R , if x ≥ tλ max ( n , u L , u R ) . Actually, even if the above structure of the Riemann solution is not available orvalid, we henceforth make the following assumption:(2.5) The unique solution of (2.1) has a finite speed of propagation for any n ,i.e., there is λ max ( n , u L , u R ) such that (2.4) holds.For instance, this is the case for strictly hyperbolic systems that may have char-acteristic families that are either not genuinely nonlinear or not linearly degener-ate, see e.g., [Liu(1975), Thm.1.2] and [Dafermos(2000), Thm. 9.5.1]. We refer to[Osher(1983), Thm. 1] for the theory of the Riemann problem for scalar conservationequations with nonconvex fluxes. In the case of general hyperbolic systems, we re-fer to [Bianchini and Bressan(2005), Section 14] for characterizations of the Riemannsolution using viscosity regularization. We also refer to [Young(2002), Thm. 2] forthe theory of the Riemann problem for the p -system with arbitrary data (i.e., withpossible formation of vacuum).The following elementary result is an important, well-known, consequence of (2.4),i.e., the Riemann solution is equal to u L for x ∈ ( −∞ , λ − t ) and equal u R for x ∈ ( λ + m t, ∞ ): Lemma 2.1.
Let u L , u R ∈ A , let u ( n , u L , u R ) be the Riemann solution to (2.1) ,let u ( t, n , u L , u R ) := (cid:82) − u ( n , u L , u R )( x, t ) d x and assume that t λ max ( n , u L , u R ) ≤ , then (2.6) u ( t, n , u L , u R ) = 12 ( u L + u R ) − t (cid:0) f ( u R ) · n − f ( u L ) · n (cid:1) . If the system (1.1) has an entropy pair ( η, q ), and if the Riemann solution isdefined to be entropy satisfying, i.e., if the following holds(2.7) ∂ t η ( u ( n , u L , u R )) + ∂ x (cid:0) q ( u ( n , u L , u R )) · n (cid:1) ≤ , in some appropriate sense (distribution sense, measure sense, etc.), then we have thefollowing additional result. Lemma 2.2.
Let ( η, q ) be an entropy pair for (1.1) and assume that (2.7) holds.Let u L , u R ∈ A and let u ( n , u L , u R ) be the Riemann solution to (2.1) . Assume that t λ max ( n , u L , u R ) ≤ , Then (2.8) η ( u ( t, n , u L , u R )) ≤ ( η ( u L ) + η ( u R )) − t ( q ( u R ) · n − q ( u L ) · n ) . Proof . Under the CFL assumption t λ max ( n , u L , u R ) ≤ , the inequality (2.7)implies that(2.9) (cid:90) − η ( u ( n , u L , u R ))( x, t ) d x ≤ ( η ( u L ) + η ( u R )) − t ( q ( u R ) · n − q ( u L ) · n ) . Jensen’s inequality η ( u ( t, n , u L , u R )) ≤ (cid:82) − η ( u ( n , u L , u R )( x, t )) d x then implies thedesired result. J.L. GUERMOND, B. POPOV
We introduce in this section the notionsof invariant sets and invariant domains. Our definitions are slightly different fromthose in [Chueh et al.(1977)Chueh, Conley, and Smoller, Hoff(1985), Smoller(1983),Frid(2001)]. We will associate invariant sets only with solutions of Riemann problemsand define invariant domains only for an approximation process.
Definition 2.3 (
Invariant set ). We say that a set A ⊂ A ⊂ R m is invariant for (1.1) if for any pair ( u L , u R ) ∈ A × A , any unit vector n ∈ S d − ( , , and any t > ,the average of the entropy solution of the Riemann problem (2.1) over the Riemannfan, say, t ( λ + m − λ − ) (cid:82) λ + m tλ − t u ( n , u L , u R )( x, t ) d x , remains in A . Note that, the above definition implies that given t > I suchthat ( λ − t, λ + m t ) ⊂ I , we have that I (cid:82) I u ( n , u L , u R )( x, t ) d x ∈ A . Note also that mostof the time expansion waves and shocks are not invariant sets.We now introduce the notion of invariant domain for an approximation process.Let X h ⊂ L ( R d ; R m ) be a finite-dimensional approximation space and let S h : X h (cid:51) u h (cid:55)−→ S h ( u h ) ∈ X h be a discrete process over X h . Henceforth we abuse the languageby saying that a member of X h , say u h , is in the set A ⊂ R m when actually we meanthat { u h ( x ) | x ∈ R } ⊂ A . Definition 2.4 (
Invariant domain ). A convex invariant set A ⊂ A ⊂ R m is saidto be an invariant domain for the process S h if and only if for any state u h in A , thestate S h ( u h ) is also in A . For scalar conservation equations the notions of invariant sets and invariantdomains are closely related to the maximum principle, see Example 2.3. In thecase of nonlinear systems, the notion of maximum principle does not apply andmust be replaced by the notion of invariant domain. To the best of our knowl-edge, the definition of invariant sets for the Riemann problem was introduced in[Nishida(1968)], and the general theory of positively invariant regions was developedin [Chueh et al.(1977)Chueh, Conley, and Smoller]. Applications and extensions tonumerical methods were developed in [Hoff(1979), Hoff(1985)] and [Frid(2001)].The invariant domain theory when m = 2 and d = 1 relies on the existenceof global Riemann invariants; the best known examples are the hyperbolic systemsof isentropic gas dynamics in Eulerian and Lagrangian form, see Example 2.4 and[Lions et al.(1996)Lions, Perthame, and Souganidis]. For results on general hyper-bolic systems, we refer to [Frid(2001)], where a characterization of invariant domainsfor the Lax-Friedrichs scheme and some flux splitting schemes is given. In particularthe existence of invariant domains is established for the above mentioned schemes forthe compressible Euler equations in the general case m = d +2 (positive density, inter-nal energy, and minimum principle on the specific entropy), see [Frid(2001), Thm. 7and Thm. 8]. Similar results have been established for various finite volume schemesin two-space dimension for the Euler equations in [Perthame and Shu(1996), Thm. 3].The objective of this paper is to propose an explicit numerical method based oncontinuous finite elements to approximate (1.1) such that any convex invariant set of(1.1) is an invariant domain for the process generated by the said numerical method.To facilitate the reading of the paper we now illustrate the abstract notions ofinvariant sets and invariant domains with some examples. Assume that m = 1 and d is arbitrary,i.e., (1.1) is a scalar conservation equation. Provided f ∈ Lip( R ; R d ), any boundedinterval is an admissible set for (1.1). For any Riemann data u L , u R , the maxi-mum speed of propagation in (2.4) is bounded by λ max ( u L , u R ) := (cid:107) f · n (cid:107) Lip( u min ,u max ) nvariant domains and C finite element approximation of hyperbolic systems u min = min( u L , u R ), u max = max( u L , u R ). If f is convex and is of class C , we have λ max ( u L , u R ) = max( | n · f (cid:48) ( u L ) | , | n · f (cid:48) ( u R ) | ) if n · f (cid:48) ( u L ) ≤ n · f (cid:48) ( u R ) and λ max ( u L , u R ) = n · ( f ( u L ) − f ( u R )) / ( u L − u R ) otherwise. Any interval [ a, b ] ⊂ R is ad-missible and is an invariant set for (1.1), i.e., if u R , u L ∈ [ a, b ], then a ≤ u ( n , u L , u R ) ≤ b for all times; this is the maximum principle. For any a ≤ b ∈ R , the interval [ a, b ]is an invariant domain for any maximum principle satisfying numerical scheme. Notethat the maximum principle can be established for a large number of numerical meth-ods (whether monotone or not), see for example [Crandall and Majda(1980)]. The one-dimensional motion of an isentropic gasis modeled by the so-called p -system, and in Lagrangian coordinates the system iswritten as follows:(2.10) (cid:40) ∂ t v + ∂ x u = 0 ,∂ t u + ∂ x p ( v ) = 0 , for ( x, t ) ∈ R × R + . Here d = 1 and m = 2. The dependent variables are the velocity u and the specificvolume v , i.e., the reciprocal of density. The mapping v (cid:55)→ p ( v ) is the pressure and isassumed to be of class C ( R + ; R ) and to satisfy(2.11) p (cid:48) < , < p (cid:48)(cid:48) . A typical example is the so-called gamma-law, p ( v ) = rv − γ , where r > γ ≥ u = ( v, u ) T , any set A in (0 , ∞ ) × R is admissible.Using the notation d µ := (cid:112) − p (cid:48) ( s ) d s , and assuming (cid:82) ∞ d µ < ∞ , the systemhas two families of global Riemann invariants:(2.12) w ( u ) = u + (cid:90) ∞ v d µ, and w ( u ) = u − (cid:90) ∞ v d µ. Note that (cid:82) ∞ d µ < ∞ if γ >
1. If γ = 1 we can use w ( u ) = u − √ r log v and w ( u ) = u + √ r log v . Let a, b ∈ R , then it can be shown that any set A ab ∈ R + × R of the form(2.13) A ab := { u ∈ R + × R | a ≤ w ( u ) , w ( u ) ≤ b } is an invariant set for the system (2.10) for γ ≥
1, see [Hoff(1985), Exp. 3.5, p. 597] fora proof in the context of parabolic regularization, or use the results from [Young(2002)]for a direct proof. Moreover, A ab is an invariant domain for the Lax-Friedrichs scheme,see [Hoff(1979), Thm. 2.1] and [Hoff(1985), Thm. 4.1].Since in the rest of the paper the maximum wave speed is the only informationwe are going to need from the Riemann solution, we give the following result. Lemma 2.5.
Let ( v L , u L ) , ( v R , u R ) ∈ R + × R with v R , v L < ∞ . Then λ max ( u L , u R ) = (cid:40)(cid:112) − p (cid:48) (min( v L , v R )) , if u L − u R > (cid:112) ( v L − v R )( p ( v R ) − p ( v L )) , (cid:112) − p (cid:48) ( v ∗ ) , otherwise,where v ∗ is the unique solution of φ ( v ) := f L ( v ) + f R ( v ) + u L − u R = 0 and f Z ( v ) := − (cid:112) ( p ( v ) − p ( v Z )( v Z − v ) , if v ≤ v Z (cid:90) vv Z d µ, if v > v Z . J.L. GUERMOND, B. POPOV
Upon setting w max1 := max( w ( u L ) , w ( u R )) and w min2 := min( w ( u L ) , w ( u R )) wehave also have v ≤ min( v L , v R , v ∗ ) , i.e., λ max ( u L , u R ) ≤ (cid:112) − p (cid:48) ( v ) , where v := ( γr ) γ − (cid:18) γ − w max1 − w min2 ) (cid:19) γ − . Proof . It is well know that the solution of the Riemann problem consists of threeconstant states u L , u ∗ , and u R connected by two waves: a 1-wave connects u L and u ∗ , and a 2-wave connects u ∗ and u R . Moreover, a vacuum forms if and only iflim v → + ∞ φ ( v ) ≥
0, see [Young(2002)] for details. In the presence of vacuum theequation φ ( v ) = 0 has no solutions and in this case we conventionally set v ∗ := + ∞ and (cid:112) − p (cid:48) ( v ∗ ) := 0. Note that since φ is an increasing and concave up function withlim v → φ ( v ) = −∞ , the solution v ∗ is unique. We also have that the maximumspeed of the exact solution is λ max ( u L , u R ) = max( (cid:112) − p (cid:48) ( v L ) , (cid:112) − p (cid:48) ( v ∗ ) , (cid:112) − p (cid:48) ( v R )).The only possibility for λ max ( u L , u R ) = (cid:112) − p (cid:48) ( v ∗ ) is if v ∗ ≤ min( v L , v R ), i.e., thesolution contains two shock waves which is equivalent to φ (min( v L , v R )) ≥
0. Us-ing the definition of φ we derive that λ max ( u L , u R ) = (cid:112) − p (cid:48) ( v ∗ ) if and only if φ (min( v L , v R )) = u L − u R − (cid:112) ( v L − v R )( p ( v R ) − p ( v L )) ≥
0. This finishes the proofof the first part of the lemma.The exact value of v ∗ can be found using Newton’s method starting with a guess v ≤ v ∗ . This guarantees that at each step of Newton’s method the estimated max-imum speed is an upper bound for the exact maximum speed. One can obtain sucha guess v by using the invariant domain property (2.13), i.e., we define the state u := ( v , u ) by w max1 = w ( u ) and w min2 = w ( u ) thereby giving v = ( γr ) γ − (cid:18) γ − w max1 − w min2 ) (cid:19) γ − . The invariant domain property guarantees that v ≤ v ∗ . Hence, the result is estab-lished. Remark 2.1.
Note that the estimate on λ max ( u L , u R ) given in Lemma 2.5 is validwhether vacuum is created or not in the Riemann solution. Remark 2.2.
We only consider the case where both u L and u R , are not vacuumstates in Lemma 2.5, since the algorithm that we propose in this paper never producesvacuum states if vacuum is not present in the initial data. Consider the compressible Euler equations(2.14) ∂ t c + ∇· ( f ( c )) = 0 , c = ρ m E , f ( c ) = mm ⊗ m ρ + p I m ρ ( E + p ) , where the independent variables are the density ρ , the momentum vector field m andthe total energy E . The velocity vector field u is defined by u := m /ρ and the internalenergy density e by e := E − | u | . The quantity p is the pressure. The symbol I denotes the identity matrix in R d . Let s be the specific entropy of the system, andassume that − s ( e, ρ − ) is strictly convex. It is known that(2.15) A r := { ( ρ, m , ρE ) | ρ ≥ , e ≥ , s ≥ r } nvariant domains and C finite element approximation of hyperbolic systems r ∈ R . It is shown in [Frid(2001), Thm.7 and 8] that the set A r is convex and is an invariant domain for the Lax-Friedrichsscheme.Let n ∈ S d − ( ,
1) and let us formulate the Riemann problem (2.1) for theEuler equations. This problem was first described in the context of dimension split-ting schemes with d = 2 in [Chorin(1976), p. 526]. The general case is treated in[Colella(1990), p. 188], see also [Toro(2009), Chapter 4.8]. We make a change of ba-sis and introduce t , . . . , t d − so that { n , t , . . . , t d − } forms an orthonormal basisof R d . With this new basis we have m = ( m, m ⊥ ) T , where m := ρu , u := u · n , m ⊥ := ρ ( u · t , . . . , u · t d − ) := ρ u ⊥ . The projected equations are(2.16) ∂ t c + ∂ x ( n · f ( c )) = , c = ρm m ⊥ E , n · f ( c ) = m ρ m + pu m ⊥ u ( E + p ) . Using the density ρ and the specific entropy s as dependent variables for the pressure, p ( ρ, s ), the linearized Jacobian is u ρ T ρ − ∂ ρ p u T ρ − ∂ s p u I T u . The eigenvalues are u , with multiplicity d , u + (cid:112) ∂ ρ p ( ρ, s ), with multiplicity 1, and u − (cid:112) ∂ ρ p ( ρ, s ), with multiplicity 1. One key observation is that the Jacobian doesnot depend on m ⊥ , see [Toro(2009), p. 150]. As a consequence the solution of theRiemann problem with data ( c L , c R ), is such that ( ρ, u, p ) is obtained as the solutionto the one-dimensional Riemann problem(2.17) ∂ t ρm E + ∂ x m ρ m + pu ( E + p ) = 0 , with e = E − m ρ with data c n L := ( ρ L , m L · n , E L ), c n R := ( ρ R , m R · n , E R ), where E Z = E Z − (cid:107) m ⊥ Z (cid:107) (cid:96) ρ Z , Z ∈ { L, R } . Moreover, for an ideal gas obeying the caloric equation of state p =( γ − ρe , it can be shown (see [Toro(2009), p. 150]) that m ⊥ is the solution of thetransport problem ∂ t m ⊥ + ∂ x ( u m ) = 0. The bottom line of this argumentation isthat the maximum wave speed in (2.16) is λ max ( c L , c R ) = max( | λ − ( c n L , c n R ) | , | λ +3 ( c n L , c n R ) | ) . where λ − ( c n L , c n R ) and λ +3 ( c n L , c n R ) are the two extreme wave speeds in the Riemannproblem (2.17) with data ( c n L , c n R ).We now determine the values of λ − ( c n L , c n R ) and λ +3 ( c n L , c n R ). We only considerthe case where both states, c L and c R , are not vacuum states, since the algorithmthat we are proposing in this paper never produces vacuum states if vacuum is notpresent in the initial data. That is, we assume ρ L , ρ R > p L , p R ≥
0. Then thelocal sound speed is given by a Z = (cid:113) γp Z ρ Z where Z is either L or R . We introduce the J.L. GUERMOND, B. POPOV following notations A Z := γ +1) ρ Z , B Z := γ − γ +1 p Z and the functions φ ( p ) := f ( p, L ) + f ( p, R ) + u R − u L (2.18) f ( p, Z ) := ( p − p Z ) (cid:16) A Z p + B Z (cid:17) if p ≥ p Z , a Z γ − (cid:18)(cid:16) pp Z (cid:17) γ − γ − (cid:19) if p < p Z , (2.19)where again Z is either L or R . It is shown in [Toro(2009), Chapter 4.3.1] that thefunction φ ( p ) ∈ C ( R + ; R ) is monotone increasing and concave down. Observe that φ (0) = u R − u L − a L γ − − a R γ − . Therefore, φ has a unique positive root if and only ifthe non-vacuum condition(2.20) u R − u L < a L γ − a R γ − p ∗ , i.e., φ ( p ∗ ) = 0 and p ∗ can be found via Newton’s method. If (2.20) does not hold we set p ∗ = 0. Then itcan be shown that, whether there is formation of vacuum or not, we have λ − ( c n L , c n R ) = u L − a L (cid:32) γ + 12 γ (cid:18) p ∗ − p L p L (cid:19) + (cid:33) , (2.21) λ +3 ( c n L , c n R ) = u R + a R (cid:32) γ + 12 γ (cid:18) p ∗ − p R p R (cid:19) + (cid:33) , (2.22)where z + := max(0 , z ). Remark 2.3. (Fast algorithm)
Note that if both φ ( p L ) > φ ( p R ) >
0, thereis no need to compute p ∗ , since in this case λ − ( u L , u R ) = u L − a L and λ +3 ( u L , u R ) = u R + a R , i.e., two rarefaction waves are present in the solution with a possible formationof vacuum. This observation is important since traditional techniques to compute p ∗ may require a large number of iterations in this situation, see [Toro(2009), p. 128].Note finally that there is no need to compute p ∗ exactly since one needs only an upperbound on λ max . A very fast algorithm, with guaranteed upper bound on λ max up toany prescribed accuracy (cid:15) of the type λ max ≤ ˜ λ max ≤ (1 + (cid:15) ) λ max , is described in[Guermond and Popov(2015b)].
3. First order method.
We describe in this section an explicit first-order finiteelement technique that, up to a CFL restriction, preserves all convex invariant sets of(1.1) that contain reasonable approximations of u . Although most of the argumentsinvoked in this section are quite standard and mimic Lax’s one-dimensional finitevolume scheme, we are not aware of the existence of such a finite-element-based schemein the literature. We want to approximate the solution of (1.1)with continuous finite elements. Let ( T h ) h> be a shape-regular sequence of affinematching meshes. The elements in the mesh sequence are assumed to be generatedfrom a finite number of reference elements denoted (cid:98) K , . . . , (cid:98) K (cid:36) . For example, themesh T h could be composed of a combination of triangles and parallelograms in twospace dimensions ( (cid:36) = 2 in this case); it could also be composed of a combination oftetrahedra, parallelepipeds, and triangular prisms in three space dimensions ( (cid:36) = 3 nvariant domains and C finite element approximation of hyperbolic systems (cid:98) K r to an arbitrary element K ∈ T h is denoted T K : (cid:98) K r −→ K and its Jacobian matrix is denoted J K , 1 ≤ r ≤ (cid:36) . Wenow introduce a set of reference Lagrange finite elements { ( (cid:98) K r , (cid:98) P r , (cid:98) Σ r ) } ≤ r ≤ (cid:36) (theindex r ∈ { (cid:36) } will be omitted in the rest of the paper to alleviate the notation).Then we define the scalar-valued and vector-valued Lagrange finite element spaces P ( T h ) = { v ∈ C ( D ; R ) | v | K ◦ T K ∈ (cid:98) P , ∀ K ∈ T h } , P ( T h ) = [ P ( T h )] m . (3.1)where (cid:98) P is the reference polynomial space defined on (cid:98) K (note that the index r hasbeen omitted). Denoting n sh := dim (cid:98) P and denoting by { (cid:98) a i } i ∈{ n sh } the Lagrangenodes of (cid:98) K , we assume that the space (cid:98) P is such that(3.2) min ≤ (cid:96) ≤ n sh (cid:98) v ( (cid:98) a (cid:96) ) ≤ (cid:98) v ( (cid:98) x ) ≤ max ≤ (cid:96) ≤ n sh (cid:98) v ( (cid:98) a (cid:96) ) , ∀ (cid:98) v ∈ (cid:98) P , ∀ (cid:98) x ∈ (cid:98) K. Denoting by P and Q the set of multivariate polynomials of total and partial degreeat most 1, respectively; the above assumption holds for (cid:98) P = P when K is a simplexand (cid:98) P = Q when K is a parallelogram or a cuboid. This assumption holds also forfirst-order prismatic elements in three space dimensions.Let { a i } i ∈{ I } be the collection of all the Lagrange nodes in the mesh T h , andlet { ϕ i } i ∈{ I } be the corresponding global shape functions. Recall that { ϕ i } i ∈{ I } forms a basis of P ( T h ) and ϕ i ( a j ) = δ ij . The Lagrange interpolation operator in P ( T h ) is denoted Π h : C ( D ) −→ P ( T h ). Recall that Π h ( v ) = (cid:80) ≤ i ≤ I v ( a i ) ϕ i . Wedenote by S i the support of ϕ i and by | S i | the measure of S i , i ∈ { I } . We alsodefine S ij := S i ∩ S j the intersection of the two supports S i and S j . Let E be a unionof cells in T h ; we define I ( E ) := { j ∈ { I } | | S j ∩ E | (cid:54) = 0 } the set that containsthe indices of all the shape functions whose support on E is of nonzero measure. Weare going to regularly invoke I ( K ) and I ( S i ) and the partition of unity property: (cid:80) i ∈I ( K ) ϕ i ( x ) = 1 for all x ∈ K .We define the operator C : P ( T h ) −→ R I so that C ( v h ) is the coordinate vectorof v h in the basis { ϕ i } i ∈{ I } , i.e., v h = (cid:80) Ii =1 C ( v h ) i ϕ i . Note that C ( v h ) i = v h ( a i ).We are also going to use capital letters for the coordinate vectors to alleviate thenotation; for instance we shall write V = C ( v h ) when the context is unambiguous.Note finally that the above assumptions on the mesh and the reference elementsimply the following property: for all x ∈ K and all K ∈ T h ,(3.3) min (cid:96) ∈I ( K ) C ( v h ) (cid:96) ≤ v h ( x ) ≤ max (cid:96) ∈I ( K ) C ( v h ) (cid:96) , ∀ v h ∈ P ( T h ) . We define similarly the C : P ( T h ) −→ ( R m ) I , i.e., v h = (cid:80) Ii =1 C ( v h ) i ϕ i , or equiva-lently C ( v h ) i = v h ( a i ).Let M ∈ R I × I be the consistent mass matrix with entries (cid:82) S ij ϕ i ( x ) ϕ j ( x ) d x ,and let M L be the diagonal lumped mass matrix with entries(3.4) m i := (cid:90) S i ϕ i ( x ) d x. The partition of unity property implies that m i = (cid:80) j ∈I ( S i ) (cid:82) ϕ j ( x ) ϕ i ( x ) d x , i.e., theentries of M L are obtained by summing the rows of M .0 J.L. GUERMOND, B. POPOV
Let u h ∈ P ( T h ) be a reasonable approximation of u (weshall be more precise in the following sections). Let n ∈ N , τ be the time step, t n be the current time, and let us set t n +1 = t n + τ . Let u nh ∈ P ( T h ) be the spaceapproximation of u at time t n and set U n = C ( u nh ). We propose to compute u n +1 h by(3.5) m i U n +1 i − U ni τ + (cid:90) D ∇· (Π h f ( u nh )) ϕ i d x − (cid:88) j ∈I ( S i ) d ij U nj = 0 , where U n +1 = C ( u n +1 h ) and the lumped mass matrix is used for the approximationof the time derivative. The coefficient d ij is an artificial viscosity for the pair ( i, j )that has yet to be clearly identified. For the time being we assume that(3.6) d ij ≥ , if i (cid:54) = j, d ij = d ji , and d ii := (cid:88) i (cid:54) = j ∈I ( S i ) − d ji . Using that ∇· (Π h f ( u nh )) = (cid:80) j f ( U nj ) ·∇ ϕ j , the above equation simplifies into(3.7) m i U n +1 i − U ni τ + (cid:88) j ∈I ( S i ) f ( U nj ) · c ij − U nj d ij = 0 , where the coefficients c ij ∈ R d are defined by(3.8) c ij = (cid:90) D ϕ i ∇ ϕ j d x, Remark 3.1. (Conservation)
The definition d ii := (cid:80) i (cid:54) = j ∈I ( S i ) − d ji implies that (cid:80) j ∈I ( S i ) d ji = 0, which in turn implies conservation, i.e., (cid:82) D u n +1 h d x = (cid:82) D u nh d x + (cid:82) D ∇· (Π h f ( u nh )) d x . Note also that the symmetry assumption in (3.6) implies d ii := (cid:80) i (cid:54) = j ∈I ( S i ) − d ij , which is often easier to compute. We motivate the choice of theartificial viscosity coefficients d ij in this section. Observing that the partition ofunity property (cid:80) j ∈I ( S i ) ϕ j = 1 and (3.6) imply conservation, i.e.,(3.9) (cid:88) j ∈I ( S i ) c ij = 0 , (cid:88) j ∈I ( S i ) d ij = 0 . we re-write (3.7) as follows:(3.10) m i U n +1 i − U ni τ = − (cid:88) j ∈I ( S i ) ( f ( U nj ) − f ( U ni )) · c ij + d ij ( U nj + U ni ) . Using again conservation, i.e., d ii = − (cid:80) i (cid:54) = j ∈I ( S i ) d ij , we finally arrive at(3.11) U n +1 i = U ni (cid:16) − (cid:88) i (cid:54) = j ∈I ( S i ) τ d ij m i (cid:17) + (cid:88) i (cid:54) = j ∈I ( S i ) τ d ij m i U n +1 ij . where we have introduced the auxiliary quantities(3.12) U n +1 ij := 12 ( U nj + U ni ) − ( f ( U nj ) − f ( U ni )) · c ij d ij . nvariant domains and C finite element approximation of hyperbolic systems
11A first key observation is that (3.11) is a convex combination provided τ is smallenough. A second key observation at this point is that upon setting n ij := c ij / (cid:107) c ij (cid:107) (cid:96) , U n +1 ij is exactly of the form u ( t, n ij , U i , U j ) as defined in (2.6) with a fake time t = (cid:107) c ij (cid:107) (cid:96) / d ij . The CFL condition tλ max ( n ij , u L , u R ) ≤ in Lemma 2.1 motivatesthe following definition for the viscosity coefficients d ij (3.13) d ij := max( λ max ( n ij , U ni , U nj ) (cid:107) c ij (cid:107) (cid:96) , λ max ( n ji , U nj , U ni ) (cid:107) c ji (cid:107) (cid:96) ) , where recall that λ max ( n ij , U i , U j ) is defined in the assumption (2.5). Remark 3.2. (Symmetry)
If either a i or a j is an interior node in the mesh, oneintegration by parts implies that c ij = − c ji , which in turn implies λ max ( n ij , U i , U j ) = λ max ( n ji , U j , U i ). In conclusion λ max ( n ij , U i , U j ) (cid:107) c ij (cid:107) (cid:96) = λ max ( n ji , U j , U i ) (cid:107) c ji (cid:107) (cid:96) if either a i or a j is an interior node. Remark 3.3. (Upwinding)
Note that in the scalar one-dimensional case when theflux f is linear, (3.5) gives the usual upwinding first-order method.
4. Stability analysis.
We analyze the stability properties of the scheme (3.5)with the viscosity defined in (3.13).
Upon defining h K := diam( K ), the globalmaximum mesh size is denoted h = max K ∈T h h K . The local minimum mesh size, h K ,for any K ∈ T h is defined as follows:(4.1) h K := 1max i (cid:54) = j ∈I ( K ) (cid:107)∇ ϕ i (cid:107) L ∞ ( S ij ) , and the global minimum mesh size is h := min K ∈T h h K . Due to the shape regularityassumption, the quantities h K and h K are uniformly equivalent, but it will turnout that using h K instead of h K gives a sharper estimate of the CFL number. Let n sh := card( I ( K )) and let us define ϑ K := n sh − . Note that(4.2) 0 < ϑ min := min ( T h ) h> min K ∈T h ϑ K < + ∞ , since there are at most (cid:36) reference elements defining the mesh sequence. We alsointroduce the mesh-dependent quantities(4.3) µ min := min K ∈T h min i ∈I ( K ) | K | (cid:90) K ϕ i ( x ) d x, µ max := max K ∈T h max i ∈I ( K ) | K | (cid:90) K ϕ i ( x ) d x. Note that µ min = µ max = n sh = d +1 for meshes uniquely composed of simplices and µ min = µ max = 2 − d for meshes uniquely composed of parallelograms and cuboids. Wenow prove the main result of the paper. Theorem 4.1.
Let A ⊂ A be an invariant set for (1.1) in the sense of Defini-tion 2.3. Assume that A is convex and (4.4) λ max ( A ) := max n ∈ S d − ( , max u L , u R ∈ A λ max ( n , u L , u R ) < ∞ , where S d − ( , is the unit sphere in R d . Assume that u h ∈ A and τ is such that (4.5) 2 τ λ max ( A ) h µ max µ min ϑ min ≤ . Then J.L. GUERMOND, B. POPOV (i) A is an invariant domain for the solution process u nh (cid:55)−→ u n +1 h for all n ≥ .(ii) Given n ≥ and i ∈ { I } , let B ⊂ A be a convex invariant set such that U nl ∈ B for all l ∈ I ( S i ) , then U n +1 i ∈ B .Proof . We prove the statement (i) by induction. Assume that u nh ∈ A for some n ≥
0; we are going to prove that u n +1 h ∈ A . Note that u h ∈ A by assumption. Let i ∈ { I } and consider the update (3.5) rewritten in the form (3.11). Observe thatupon defining n ij := c ij / (cid:107) c ij (cid:107) (cid:96) , the quantity U n +1 ij defined in (3.12) is exactly ofthe form u ( t, n ij , U ni , U nj ) as defined in (2.6) with the flux f · n ij and the fake time t = (cid:107) c ij (cid:107) (cid:96) / d ij . The definition d ij ≥ λ max ( n ij , U ni , U nj ) (cid:107) c ij (cid:107) (cid:96) , (4.6)is the CFL condition for the conclusions of Lemma 2.1 to hold with fake time t = (cid:107) c ij (cid:107) (cid:96) / d ij . Since A is a convex invariant set we have U n +1 ij := u ( t, n ij , U ni , U nj ) ∈ A for all j ∈ I ( S i ). Let us now prove that (3.11) is indeed a convex combination byproving that 1 − (cid:80) i (cid:54) = j ∈I ( S i ) 2 τd ij m i = 1 + τd ii m i ≥
0. Note first that (cid:107) c ij (cid:107) (cid:96) ≤ (cid:90) S ij (cid:107)∇ ϕ j (cid:107) (cid:96) ϕ i d x ≤ h − (cid:90) S ij ϕ i d x ≤ h − µ max | S ij | . The definition of d ii implies that − d ii ≤ λ max ( A ) h µ max (cid:88) i (cid:54) = j ∈I ( S i ) | S ij | ≤ λ max ( A ) h µ max ϑ min | S i | . The using that µ min | S i | ≤ m i , we infer that − τ d ii m i ≤ τ λ max ( A ) h µ max µ min ϑ min ≤ , which proves the result owing to the CFL assumption (4.5). Hence (3.11) defines U n +1 i as a convex combination between U ni and the collection of states { U n +1 ij } j ∈I ( S i ) . Theconvexity of A implies that U n +1 i ∈ A , since U ni ∈ A by assumption and we haveestablished above that U n +1 ij ∈ A for all j ∈ I ( S i ). The space approximation beingpiecewise linear, a convexity argument implies again that u n +1 h ∈ A , which proves theinduction assumption.Note in passing that we have also proved the following local invariance property:given any convex invariant set B ⊂ A that contains { U nl } l ∈I ( S i ) , U n +1 i is also in B ,i.e., the local statement (ii) holds. This completes the proof. Remark 4.1.
The arguments invoking the convex combination (3.11) and the one-dimensional Riemann averages (3.12) are similar in spirit to those used in the proofof Theorem 3 in [Perthame and Shu(1996)].
We now derive a local entropy inequality.
Theorem 4.2.
Let A ⊂ A be a convex invariant set for (1.1) . Let ( η, q ) be anentropy pair for (1.1) . Assume that (2.7) holds for any Riemann data ( u L , u R ) in A , and any n ∈ S d − ( , . Assume also that (4.4) and (4.5) hold, then we have thefollowing for any n ≥ and any i ∈ { I } : m i τ ( η ( U n +1 i ) − η ( U ni )) + (cid:90) D ∇· (Π h q ( u nh )) ϕ i d x + (cid:88) i (cid:54) = j ∈I ( S i ) d ij η ( U nj ) ≤ . (4.7) nvariant domains and C finite element approximation of hyperbolic systems Proof . Let ( η, q ) be an entropy pair for the system (1.1). Let i ∈ { I } , thenrecalling (3.11), the CFL condition and the convexity of η imply that η ( U n +1 i ) ≤ (cid:16) − (cid:88) i (cid:54) = j ∈I ( S i ) τ d ij m i (cid:17) η ( U ni ) + (cid:88) i (cid:54) = j ∈I ( S i ) τ d ij m i η ( U n +1 ij ) . Owing to Lemma 2.2 we have η ( U n +1 ij ) ≤ ( η ( U ni ) + η ( U nj )) − t ( q ( U nj ) · n ij − q ( U ni ) · n ij ) . with t = (cid:107) c ij (cid:107) (cid:96) / d ij ; hence, m i τ ( η ( U n +1 i ) − η ( U ni )) ≤ (cid:88) i (cid:54) = j ∈I ( S i ) d ij ( η ( U n +1 ij ) − η ( U ni )) ≤ (cid:88) i (cid:54) = j ∈I ( S i ) d ij ( η ( U nj ) − η ( U ni )) − (cid:107) c ij (cid:107) (cid:96) ( q ( U nj ) · n ij − q ( U ni ) · n ij ) . The conclusion follows from the definitions of n ij , c ij and d ij . Remark 4.2.
One recovers the equation (3.5) from (4.7) with η ( v ) = v . Note alsothat (4.7) gives the global entropy inequality (cid:80) ≤ i ≤ I m i η ( U n +1 i ) ≤ (cid:80) ≤ i ≤ I m i η ( U ni ). Remark 4.3.
The meaning of the entropy inequality (2.7) might be somewhatambiguous in some cases, especially when u is a measure. Since it is only the in-equality (2.8) that is really needed in the proof of Theorem 4.2, we could replacethe assumption (2.7) by (2.8). This would avoid having to invoke measure solutionssince u ( t, n , u L , u R ) should always be finite for the Riemann problem (2.1) to have areasonable (physical) meaning. In the formulation (3.5) theterm (cid:80) j ∈I ( S i ) d ij U j models some edge-based dissipation, i.e., d ij is a dissipationcoefficient associated with the pair of degrees of freedom of indices ( i, j ). This for-mulation is related in spirit to that of local extremum diminishing (LED) schemesdeveloped for scalar conservation equations in [Kuzmin and Turek(2002), Eq. (32)-(33)], see also [Jameson(1995), § −∇· ( ν ∇ ψ ). For instance, assuming that theviscosity field ν is piecewise constant over each mesh cell K ∈ T h , we write:(4.8) (cid:90) D −∇· ( ν ∇ ψ ) ϕ i d x = (cid:88) K ⊂ S i ν K (cid:90) K ∇ ψ ·∇ ϕ i d x. Unfortunately, it has been shown in [Guermond and Nazarov(2013)] that the bilinearform ( ψ, ϕ ) −→ (cid:82) K ∇ ψ ·∇ ϕ d x is not robust with respect to the shape of the cells.More specifically, the convex combination argument, which is essential to prove themaximum principle for scalar conservation equations in arbitrary space dimensionwith continuous finite elements, can be made to work only if (cid:82) S ij ∇ ϕ i ·∇ ϕ j d x < ϕ i , ϕ j , with common support of nonzero measure. Thisis the well-known acute angle condition assumption, which a priori excludes a lot of4 J.L. GUERMOND, B. POPOV meshes in particular in three space dimensions. To avoid this difficulty, it is proposedin [Guermond and Nazarov(2013)] to replace (4.8) by (cid:80) K ⊂ S i ν K b K ( ψ, ϕ i ), where(4.9) b K ( ϕ j , ϕ i ) = − ϑ K | K | if i (cid:54) = j, i, j ∈ I ( K ) , | K | if i = j, i, j ∈ I ( K ) , i (cid:54)∈ I ( K ) or j (cid:54)∈ I ( K ) . The essential properties of b K can be summarized as follows: Lemma 4.3.
There is c > depending only on the collection { ( (cid:98) K r , (cid:98) P r , (cid:98) Σ r ) } ≤ r ≤ (cid:36) and the shape-regularity, such that the following identities hold for all K ∈ T h and all u h , v h ∈ P ( T h ) : b K ( ϕ i , ϕ j ) = b K ( ϕ j , ϕ i ) , b K ( ϕ i , (cid:88) j ∈I ( K ) ϕ j ) = 0 . (4.10) b K ( u h , v h ) = ϑ K | K | (cid:88) i ∈I ( K ) (cid:88) I ( K ) (cid:51) j
Let { ν nK } K ∈T h be defined by (4.14) ν K := max i (cid:54) = j ∈I ( K ) λ max ( n ij , U i , U j ) (cid:107) c ij (cid:107) (cid:96) (cid:80) T ⊂ S ij − b T ( ϕ j , ϕ i ) . Then the conclusions of Theorem 4.1 and Theorem 4.2 hold under the assumptions (4.4) and (4.5) and with the solution process u nh −→ u n +1 h , n ≥ , defined by (4.13) .Proof . Let us denote ˜ d ij := − (cid:80) K ∈ S ij ν nK b K ( ϕ j , ϕ i ), then (4.13) can be recast asfollows: m i U n +1 i − U ni τ + (cid:88) j ∈I ( S i ) f ( U nj ) · c ij − U nj ˜ d ij = 0 , which in turn implies that U n +1 i = U ni (cid:16) − (cid:88) i (cid:54) = j ∈I ( S i ) τ ˜ d ij m i (cid:17) + (cid:88) i (cid:54) = j ∈I ( S i ) τ ˜ d ij m i U n +1 ij . nvariant domains and C finite element approximation of hyperbolic systems U n +1 ij := 12 ( U nj + U ni ) − ( n ij · f ( U nj ) − n ij · f ( U ni )) (cid:107) c ij (cid:107) (cid:96) d ij . Here again U n +1 ij is of the form u ( t, n ij , U i , U j ) as defined in (2.6) with the fake time t = (cid:107) c ij (cid:107) (cid:96) / d ij , hence we need to make sure that λ max ( n ij , u L , u R ) (cid:107) c ij (cid:107) (cid:96) / d ij ≤ to preserve the invariant domain property. Recalling that d ij has been defined by d ij := λ max ( n ij , u L , u R ) (cid:107) c ij (cid:107) (cid:96) (see (3.6)), the above condition reduces to showingthat d ij ≤ ˜ d ij . The definitions of ν K and ˜ d ij implies that˜ d ij = − (cid:88) K ∈ S ij ν nK b K ( ϕ j , ϕ i ) ≥ − (cid:88) K ∈ S ij λ max ( n ij , U i , U j ) (cid:107) c ij (cid:107) (cid:96) (cid:80) T ⊂ S ij − b T ( ϕ j , ϕ i ) b K ( ϕ j , ϕ i ) ≥ − (cid:88) K ∈ S ij d ij (cid:80) T ⊂ S ij − b T ( ϕ j , ϕ i ) b K ( ϕ j , ϕ i ) = d ij , whence the desired result. We now prove that 1 − (cid:80) i (cid:54) = j ∈I ( S i ) 2 τ ˜ d ij m i ≥ d ij ≤ λ max ( A ) h − µ max | S ij | ,hence ν K ≤ λ max ( A ) µ max h max k (cid:54) = l ∈ I ( K ) | S kl | (cid:80) T ⊂ S kl − b T ( ϕ k , ϕ l ) , which in turn implies that˜ d ij ≤ λ max ( A ) µ max h (cid:88) K ⊂ S ij − b K ( ϕ i , ϕ j ) max k (cid:54) = l ∈ I ( K ) | S kl | (cid:80) T ⊂ S kl − b T ( ϕ k , ϕ l )Recalling the definition of b T (( ϕ k , ϕ l ) we have (cid:80) T ⊂ S kl − b T ( ϕ k , ϕ l ) ≥ ϑ min | T | = ϑ min | S kl | ; hence˜ d ij ≤ λ max ( A ) µ max ϑ min h (cid:88) K ⊂ S ij − b K ( ϕ i , ϕ j ) = λ max ( A ) µ max ϑ min h (cid:88) K ⊂ S ij ϑ K | K | . Finally we have − ˜ d ii := (cid:88) i (cid:54) = j ∈I ( S i ) ˜ d ij ≤ λ max ( A ) µ max ϑ min h (cid:88) i (cid:54) = j ∈I ( S i ) (cid:88) K ⊂ S ij ϑ K | K | = λ max ( A ) µ max ϑ min h | S i | . This means that the bound on − ˜ d ii is the same as that on − d ii in the proof ofTheorem 4.1. This concludes the proof.
5. Numerical illustrations.
We illustrate in this section the method describedin the paper, i.e., (3.5)-(3.13), and discuss possible variants.
We give in thissection a counter-example showing that a method that is formally first-order consistentand satisfies the invariant domain property may not necessarily be convergent.To illustrate or point, let us focus our attention on scalar conservation equationsand let us consider an algebraic approach that is sometimes used in the literature, see6
J.L. GUERMOND, B. POPOV e.g., [Kuzmin et al.(2005)Kuzmin, L¨ohner, and Turek, p. 163], [Kuzmin and Turek(2002),Eq. (32)-(33)]. Instead of constructing a convex combination involving (entropy sat-isfying) intermediate states like in (3.11), we re-write (3.10) as follows:(5.1) m i U n +1 i − U ni τ = − (cid:88) i (cid:54) = j ∈I ( S i ) ( f ( U nj ) − f ( U ni )) · c ij + (cid:88) j ∈I ( S i ) d ij U nj . Or, equivalently(5.2) m i U n +1 i − U ni τ = − (cid:88) i (cid:54) = j ∈I ( S i ) f ( U nj ) − f ( U ni ) U nj − U ni · c ij ( U nj − U ni ) + (cid:88) j ∈I ( S i ) d ij U nj . Let us set k ij := f ( U nj ) − f ( U ni ) U nj − U ni · c ij , (with k ij := 0 if U nj = U ni ), then(5.3) U n +1 i = U ni (cid:16) − τm i (cid:88) i (cid:54) = j ∈I ( S i ) ( − k ij + d ij ) (cid:17) + (cid:88) i (cid:54) = j ∈I ( S i ) τm i ( − k ij + d ij ) U nj . Let us finally set(5.4) d ij := max(0 , k ij , k ji ) , i (cid:54) = j, and d ii := − (cid:88) i (cid:54) = j ∈I ( S i ) d ij . This choice implies that − k ij + d ij ≥ i ∈ { N } , j ∈ I ( S i ). As a result, U n +1 i ∈ conv { U nj , j ∈ I ( S i ) } under the appropriate CFL condition; hence, the solu-tion process u nh (cid:55)−→ u n +1 h described above in (5.3)-(5.4) satisfies the maximum princi-ple. Although, this technique looks reasonable a priori, it turns out that it is not dif-fusive enough to handle general fluxes as discussed in [Guermond and Popov(2015a), § n ij · f ( U nj ) − f ( U ni ) U nj − U ni , which is invoked in the above definition. This definition of the wavespeed is correct in shocks, i.e., if the Riemann problem with data ( U i , U j ) is a sim-ple shock; but it may not be sufficient if the Riemann solution is an expansion or acomposite wave, which is likely to be the case if f is not convex.We now illustrate numerically the observation made above. We consider the so-called KPP problem proposed in [Kurganov et al.(2007)Kurganov, Petrova, and Popov].It is a two-dimensional scalar conservation equation with a non-convex flux:(5.5) ∂ t u + ∇ · f ( u ) = 0 , u ( x ,
0) = u ( x ) = (cid:40) π , if (cid:112) x + y ≤ , π , otherwise . , where f ( u ) = (sin u, cos u ). This is a challenging test case for many high-order numeri-cal schemes because the solution has a two-dimensional composite wave structure. Forexample, it has been shown in [Kurganov et al.(2007)Kurganov, Petrova, and Popov]that some central-upwind schemes based on WENO5, Minmod 2 and SuperBee re-constructions converge to non-entropic solutions.The computational domain [ − , × [ − . , .
5] is triangulated using non-uniformmeshes and the solution is approximated up to t = 1 using continuous P finite ele-ments (29871 nodes, 59100 triangles). The time stepping is done with SSP RK3. Thesolution shown in the left panel of Figure 5.1 is obtained using (5.4) for the definition nvariant domains and C finite element approximation of hyperbolic systems d ij . The numerical solution produces very sharp, non-oscillating, entropy violatingshocks, the reason being that the artificial viscosity is not large enough. Note that thesolution is maximum principle satisfying (the local maximum principle is satisfied atevery grid point and every time step) and no spurious oscillations are visible. The nu-merical process converges to a nice-looking (wrong) piecewise smooth weak solution.The numerical solution shown in the right panel of Figure 5.1 is obtained by using ourdefinition of d ij , (3.13) (note in passing that the results obtained with (4.13)-(4.14)together with (3.13) are indistinguishable from this solution). The expected helicoidalcomposite wave is clearly visible; this is the unique entropy satisfying solution. Fig. 5.1 . KPP solution with continuous P elements (29871 nodes, 59100 triangles). Left:entropy violating solution using (3.5) - (5.4) ; Right: entropy satisfying solution using (3.5) - (3.13) . In conclusion, the above counter-example shows that satisfying the invariant do-main property/maximum principle does not imply convergence, even for a first-ordermethod. It is also essential that the method satisfies local entropy inequalities to beconvergent; this is the case of our method (3.5)-(3.13) (see Theorem 4.2), but it is notthe case of the algebraic method (5.3)-(5.4).
Remark 5.1.
The reader should be aware that we are citing [Kuzmin et al.(2005)Kuzmin, L¨ohner, and Turek,p. 163], [Kuzmin and Turek(2002), Eq. (32)-(33)] a little bit out of context. Thescheme as originally presented in the above references was only meant to solve thelinear transport equation, and as such it is a perfectly good method. Problems arisewith (5.4) only when one extends the methodology to nonlinear nonconvex fluxes, aswe did in (5.2).
The construction of the intermediate states in (3.10) isnot unique. For instance we can extend a construction used by [Hoff(1979), Cor.1] in one space dimension for the p -system. Let us assume that i ∈ { , . . . , N } issuch that every j ∈ I ( S i ) \ { i } , there is a unique σ i ( j ) ∈ I ( S i ) \ { i, j } such that c ij := (cid:82) S i φ i ∇ φ j d x = − (cid:82) S i φ i ∇ φ σ i ( j ) d x =: − c iσ i ( j ) . This property holds in onespace dimension for any mesh if a i is an interior node. It holds in higher spacedimension provided the mesh has symmetry properties and a i is an interior node; forinstance it holds if the mesh is centrosymmetric, i.e., the support of φ i is symmetricwith respect to the node a i for any i ∈ { , . . . , N } . Then we can re-write (3.7) asfollows:(5.6) m i U n +1 i − U ni τ = d ii U ni − (cid:88) j ∈J ( S i ) ( f ( U nj ) − f ( U nσ i ( j ) )) · c ij + d ij U nj + d iσ i ( j ) U nσ i ( j ) . J.L. GUERMOND, B. POPOV where the set J ( S i ) ⊂ I ( S i ) is such that σ i : J ( S i ) −→ σ i ( J ( S i )) is bijective and J ( S i ) ∪ σ i ( J ( S i )) = I ( S i ) \ { i } . Then upon recalling that d ii := − (cid:80) j ∈J ( S i ) ( d ij + d iσ i ( j ) ), we have(5.7) U n +1 i = U ni (cid:18) − (cid:88) j ∈J ( S i ) τm i ( d ij + d iσ i ( j ) ) (cid:19) + (cid:88) j ∈J ( S i ) τ ( d ij + d iσ i ( j ) ) m i U n +1 ij , where we have defined the intermediate state U n +1 ij by(5.8) U n +1 ij = d iσ i ( j ) d ij + d iσ i ( j ) U nσ i ( j ) + d ij d ij + d iσ i ( j ) U nj − ( f ( U nj ) − f ( U nσ i ( j ) )) · c ij d ij + d iσ i ( j ) . The state U n +1 ij is of the form u ( t, n ij , U nσ i ( j ) , U nj ) := (cid:82) α R α L u ( n ij , U nσ i ( j ) , U nj )( x, t ) d x ,where α L = − d iσi ( j ) d ij + d iσi ( j ) , α R = d ij d ij + d iσi ( j ) and t := (cid:107) c ij (cid:107) (cid:96) d ij + d iσi ( j ) , provided d iσ i ( j ) ≥ ( λ − ) − ( n ij , U nσ i ( j ) , U nj ) (cid:107) c ij (cid:107) (cid:96) , ∀ j ∈ J ( S i ) , (5.9) d ij ≥ ( λ + m ) + ( n ij , U nσ i ( j ) , U nj ) (cid:107) c ij (cid:107) (cid:96) , ∀ j ∈ J ( S i ) , (5.10)where we defined x + = max( x,
0) and x − = − min( x, J i ( S i )is(5.11) min( d ij , d iσ i ( j ) ) ≥ λ max ( n ij , U nσ i ( j ) , U nj ) (cid:107) c ij (cid:107) (cid:96) , j ∈ J ( S i ) . Note that the above argument holds only if a i is an interior node satisfying thesymmetry property c ij = − c iσ i ( j ) . If this is not the case, then we can always use thelower bound (4.6), i.e., d ij ≥ λ max ( n ij , U ni , U nj ) (cid:107) c ij (cid:107) (cid:96) .In conclusion the diffusion matrix ( d ij ) ≤ i,j ≤ N can be constructed as follows:(1) For every node i satisfying the symmetry property c ij = − c iσ i ( j ) for every j ∈ J ( S i ), we define (cid:101) d ij = (cid:101) d iσ i ( j ) = λ max ( n ij , U nσ i ( j ) , U nj ) (cid:107) c ij (cid:107) (cid:96) ; (2) For everyother index i not satisfying the symmetry property mentioned above, we define (cid:101) d ij = λ max ( n ij , U ni , U nj ) (cid:107) c ij (cid:107) (cid:96) ; (3) We construct the diffusion matrix by setting d ij :=max( (cid:101) d ij , (cid:101) d ji ) for j (cid:54) = i and d ii := − (cid:80) i (cid:54) = j ∈I ( S i ) d ij . This construction guarantees con-servation, i.e., (cid:80) i ∈I ( S j ) d ij = 0 and first-order consistency, i.e., (cid:80) j ∈I ( S i ) d ij = 0. Remark 5.2.
Quite surprisingly, in the case of scalar linear transport the aboveconstruction and the construction done in § We show in this sec-tion that the invariance property and what is usually understood in the literatureas monotonicity are two different concepts and just looking at monotonicity may bemisleading.
We consider the p-system and solve the Riemann problem cor-responding to the initial data ( v L , u L ) = (1 , v R , u R ) = (2 γ − , γ − ). The compu-tational domain is the segment [0 ,
1] and the separation between the left and rightstates is set at x = 0 .
75. The solution is a single rarefaction wave from the firstfamily (i.e., w ( v L , u L ) = w ( v R , u R )):(5.12) v ( x, t ) = x − x t ≤ − x − xt ) − γ +1 if − ≤ x − x t ≤ − − γ +1 γ − γ − otherwise nvariant domains and C finite element approximation of hyperbolic systems u ( x, t ) = x − x t ≤ − γ − (cid:16) − ( x − xt ) γ − γ +1 (cid:17) if − ≤ x − x t ≤ − − γ +1 γ − γ − otherwiseThis case is such that ( v ∗ , u ∗ ) = ( v R , u R ), hence the second wave corresponding tothe eigenvalues λ ± is not present. We use continuous piecewise linear finite elementswith the algorithm (3.5)-(3.13). The time stepping is done with the SSP RK3 tech-nique. We show the profile of v at t = 0 .
75 in Figure 5.2 for meshes composed of10 , × , × , , × , × , 10 , 2 × cells. We observe that the profile Fig. 5.2 . Left: v -profile for the p-system at t = 0 . , grid points. Right: close up view ofthe v -profile for various grid sizes: , × , × , , × , × , grid points. is not monotone. There is an overshoot at the right of the foot of the (left-going)wave. Actually this overshoot does not violate the invariant domain property; wehave verified numerically that, at every time step and for every grid point in eachmesh, the numerical solution is in the smallest invariant domain of type (2.13) thatcontains the piecewise linear approximation of the initial data. This result seems abit surprising, but it is perfectly compatible with Theorem 4.1. Since the numericalsolution cannot stay on the exact rarefaction wave (green line connecting U L and U L in Figure 5.3), the second wave reappears in the form of an overshoot at the end ofthe rarefaction wave (see right panel of the Figure 5.2). Fig. 5.3 . The overshooting mechanism for a single rarefaction wave in the phase space for thep-system. Initial data in black; additional points after one time step in red; after two time steps inblue. Observe the position of U . Let ( U L . . . , U L , U R . . . , U R ) be the initial sequence of degrees of freedom. Af-ter one time step two additional points appear in the phase space, denoted on Fig-0 J.L. GUERMOND, B. POPOV ure 5.3 by U and U . Because of the invariant domain property, these pointsare under the rarefaction wave. Then the sequence of degrees of freedom at time t = τ is ( U L . . . , U L , U , U , U R . . . , U R ). Six additional points U , . . . , U appearafter two time steps and the sequence of degrees of freedom at time t = 2 τ is( U L . . . , U L , U , . . . , U , U R . . . , U R ). The point U is the one whose v -componentmay overshoot because the exact solution of the Riemann problem with the left state U and the right state U R is composed of two rarefaction waves and the maximumvalue of v on these rarefactions is necessarily larger than v R (see red line in Fig-ure 5.3). Note that this is not a Gibbs phenomenon at all; in particular the amplitudeof the overshoot decreases as the mesh is refined as shown in the close up view inthe right panel of the Figure 5.2. This phenomenon is actually very common innumerical simulations of hyperbolic systems but is rarely discussed; it is sometimescalled ”start up error” in the literature, see for example the comments on page 592in [Kurganov and Tadmor(2002)] and the comments at the bottom of page 1005 in[Liska and Wendroff(2003)]. The (relative) L -norm of the error on both v and u at t = 0 .
75 is shown in Table 5.1. The method converges with an order close to 0 . /h v rate u rate10 × × × × × × Table 5.1
Convergence rates for the p-system
We consider now the compressibleEuler equations. We solve the Riemann problem also known in the literature as theLeblanc Shocktube. The data are as follows: γ = and ρ L = 1 . , u L = 0 . , p L = 0 . ρ R = 0 . , u R = 0 . , p R = 10 − . The structure of the solution is standard; it consists of a rarefaction wave moving tothe left, a contact discontinuity in the middle and a shock moving to the right. Thedensity profile is monotone. We solve this problem with the algorithm (3.5)-(3.13)using piecewise linear finite elements. The density profile computed with 50 , , , ,
000 and 800 ,
000 grid points is shown in the left panel of Fig-ure 5.4. The right panel in the figure shows a close up view of the region at the footof the expansion wave. Of course the scheme does not have any problem with thepositivity of the density and the internal energy, but we observe that the numericalprofile is not monotone; there is a small dip at the foot of the expansion. There isnothing wrong here, since, for each mesh, the numerical solution is guaranteed byTheorem 4.1 to be in the smallest convex invariant set that contains the Riemanndata. This phenomenon is similar to what has been observed for the p-system in theprevious section. This example shows again that the invariant domain property is a nvariant domains and C finite element approximation of hyperbolic systems Fig. 5.4 . Left: Density profile for the Leblanc Shocktube at t = 0 . . Right: close up view of thedensity profile at the foot of the rarefaction wave. different concept than monotonicity, and just looking at monotonicity is not enoughto understand hyperbolic systems.
6. Concluding remarks.
We have proposed a numerical method to solve hy-perbolic systems using continuous finite elements and forward Euler time stepping.The properties of the method are based on the introduction of an artificial dissipa-tion that is defined so that any convex invariant sets is an invariant domain for themethod. The main result of the paper are Theorem 4.1 and Theorem 4.2. The methodis formally first-order accurate with respect to space and can be made higher-orderwith respect to the time step by using any explicit Strong Stability Preserving timestepping technique. Although, the argumentation of the proof of Theorem 4.1 re-lies on the notion of Riemann problems, the algorithm does not require to solve anyRiemann problem. The only information needed is an upper bound on the local max-imum speed. Our next objective is to work on a generalization of the FCT technique(see [Kuzmin et al.(2005)Kuzmin, L¨ohner, and Turek]) to make the method at leastformally second-order accurate in space and still be domain invariant.
REFERENCES[Bianchini and Bressan(2005)] S. Bianchini and A. Bressan. Vanishing viscosity solutions of nonlin-ear hyperbolic systems.
Ann. of Math. (2) , 161(1):223–342, 2005.[Bressan(2000)] A. Bressan.
Hyperbolic systems of conservation laws , volume 20 of
Oxford LectureSeries in Mathematics and its Applications . Oxford University Press, Oxford, 2000. Theone-dimensional Cauchy problem.[Chorin(1976)] A. J. Chorin. Random choice solution of hyperbolic systems.
J. ComputationalPhys. , 22(4):517–533, 1976.[Chueh et al.(1977)Chueh, Conley, and Smoller] K. N. Chueh, C. C. Conley, and J. A. Smoller. Pos-itively invariant regions for systems of nonlinear diffusion equations.
Indiana Univ. Math.J. , 26(2):373–392, 1977.[Colella(1990)] P. Colella. Multidimensional upwind methods for hyperbolic conservation laws.
J.Comput. Phys. , 87(1):171–200, 1990.[Crandall and Majda(1980)] M. G. Crandall and A. Majda. Monotone difference approximations forscalar conservation laws.
Math. Comp. , 34(149):1–21, 1980.[Dafermos(2000)] C. M. Dafermos.
Hyperbolic conservation laws in continuum physics , volume 325of
Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathe-matical Sciences] . Springer-Verlag, Berlin, 2000.[Frid(2001)] H. Frid. Maps of convex sets and invariant regions for finite-difference systems ofconservation laws.
Arch. Ration. Mech. Anal. , 160(3):245–269, 2001.[Guermond and Nazarov(2013)] J.-L. Guermond and M. Nazarov. A maximum-principle preserving J.L. GUERMOND, B. POPOV C finite element method for scalar conservation equations. Comput. Methods Appl. Mech.Engrg. , 272:198–213, 2013.[Guermond and Popov(2015a)] J.-L. Guermond and B. Popov. Error estimates of a first-order la-grange finite element technique for nonlinear scalar conservation equations.
SIAM J. Nu-mer. Anal. , 2015a. in review.[Guermond and Popov(2015b)] J.-L. Guermond and B. Popov. Fast estimation of the maximumspeed in a riemann problem. 2015b. Submitted.[Hoff(1979)] D. Hoff. A finite difference scheme for a system of two conservation laws with artificialviscosity.
Math. Comp. , 33(148):1171–1193, 1979.[Hoff(1985)] D. Hoff. Invariant regions for systems of conservation laws.
Trans. Amer. Math. Soc. ,289(2):591–610, 1985.[Jameson(1995)] A. Jameson. Positive schemes and shock modelling for compressible flows.
Internat.J. Numer. Methods Fluids , 20(8-9):743–776, 1995. Finite elements in fluids—new trendsand applications (Barcelona, 1993).[Kurganov and Tadmor(2002)] A. Kurganov and E. Tadmor. Solution of two-dimensional Riemannproblems for gas dynamics without Riemann problem solvers.
Numer. Methods PartialDifferential Equations , 18(5):584–608, 2002.[Kurganov et al.(2007)Kurganov, Petrova, and Popov] A. Kurganov, G. Petrova, and B. Popov.Adaptive semidiscrete central-upwind schemes for nonconvex hyperbolic conservation laws.
SIAM Journal on Scientific Computing , 29(6):2381–2401, 2007.[Kuzmin and Turek(2002)] D. Kuzmin and S. Turek. Flux correction tools for finite elements.
Jour-nal of Computational Physics , 175(2):525–558, 2002.[Kuzmin et al.(2005)Kuzmin, L¨ohner, and Turek] D. Kuzmin, R. L¨ohner, and S. Turek.
Flux–Corrected Transport . Scientific Computation. Springer, 2005. 3-540-23730-5.[Lax(1957)] P. D. Lax. Hyperbolic systems of conservation laws. II.
Comm. Pure Appl. Math. , 10:537–566, 1957.[Lions et al.(1996)Lions, Perthame, and Souganidis] P.-L. Lions, B. Perthame, and P. E. Souganidis.Existence and stability of entropy solutions for the hyperbolic systems of isentropic gasdynamics in Eulerian and Lagrangian coordinates.
Comm. Pure Appl. Math. , 49(6):599–638, 1996.[Liska and Wendroff(2003)] R. Liska and B. Wendroff. Comparison of several difference schemes on1D and 2D test problems for the Euler equations.
SIAM J. Sci. Comput. , 25(3):995–1017(electronic), 2003.[Liu(1975)] T. P. Liu. The Riemann problem for general systems of conservation laws.
J. DifferentialEquations , 18:218–234, 1975.[Nishida(1968)] T. Nishida. Global solution for an initial boundary value problem of a quasilinearhyperbolic system.
Proc. Japan Acad. , 44:642–646, 1968.[Osher(1983)] S. Osher. The Riemann problem for nonconvex scalar conservation laws and Hamilton-Jacobi equations.
Proc. Amer. Math. Soc. , 89(4):641–646, 1983.[Perthame and Shu(1996)] B. Perthame and C.-W. Shu. On positivity preserving finite volumeschemes for Euler equations.
Numer. Math. , 73(1):119–130, 1996.[Smoller(1983)] J. Smoller.
Shock waves and reaction-diffusion equations , volume 258 of
Grundlehrender Mathematischen Wissenschaften [Fundamental Principles of Mathematical Science] .Springer-Verlag, New York-Berlin, 1983.[Toro(2009)] E. F. Toro.
Riemann solvers and numerical methods for fluid dynamics . Springer-Verlag, Berlin, third edition, 2009. A practical introduction.[Young(2002)] R. Young. The p -system. I. The Riemann problem. In The legacy of the inversescattering transform in applied mathematics (South Hadley, MA, 2001) , volume 301 of