Murman-Roe scheme The idea introduced by Roe for Euler equations was to solve the Riemann problem approximately. In case of scalar problem, the idea is to replace non-linear PDE u t + f ′ ( u ) u x = 0 with a linear PDE w t + a j + 1 2 w x = 0 2 ≈ f ′ ( u ( x j + 1 and solve this exactly. We must choose a j + 1 2 , t n )) � f j +1 − f j u j � = u j +1 u j +1 − u j a j + 1 2 = f ′ ( u j ) otherwise Solving the Riemann problem and evaluating the solution on ξ = 0 gives the Roe flux � f j a j + 1 2 > 0 2 = 1 2( f j + f j +1 ) − 1 g R 2 | a j + 1 2 | ( u j +1 − u j ) = j + 1 f j +1 otherwise 17 / 111
Murman-Roe scheme Being based on shocks only, the approximate solution is not good at modeling rarefactions. This can lead to entropy violating solutions. Consider the initial data � − 1 j ≤ − 1 u 0 j = +1 j ≥ 0 The Roe scheme gives the solution j = u 0 u n j which is a stationary shock and hence a weak solution. But the correct solution is a rarefaction. Note that the numerical viscosity vanishes at the initial discontinuity, which is the cause of the unphysical solution. 18 / 111
Monotone scheme The exact solutions have the property that if u ( x, 0) ≥ v ( x, 0) , then u ( x, t ) ≥ v ( x, t ) a.e in x and t . This motivates the notion of monotone schemes. Write a general finite volume scheme in the form u n +1 = H ( u n j − k , . . . , u n j + k ) j We say that the scheme is monotone if H is an increasing function of all its arguments . A 3-point scheme is of the form u n +1 = H ( u n j − 1 , u n j , u n j +1 ) = u n j − λ [ g ( u n j , u n j +1 ) − g ( u n j − 1 , u n j )] j This is monotone if g ( · , · ) is an increasing function of its first argument and a decreasing function of its second argument, and provided a CFL condition is satisfied. (Exercise: Check this for the Lax-Friedrich flux) A monotone scheme is consistent with any entropy condition . This would make monotone schemes to be the ideal choice but unfortunately, monotone schemes are atmost first order accurate . Hence the notion of monotone schemes is not very useful if we want to construct high order schemes. 19 / 111
TVD schemes The total variation of a grid function u h is defined as � TV ( u h ) = | u j +1 − u j | j The exact solutions of conservation laws have the property that their total variation does not increase with time. Hence we demand the same from the numerical solutions. A numerical scheme is said to be total variation diminishing if TV ( u n +1 ) ≤ TV ( u n h ) h Practically this helps to prevent the appearance of spurious oscillations in the case of discontinuous numerical solutions. To check if a scheme is TVD, we write it in incremental form u n +1 = u n j + C n 2 ∆ u n 2 − D n 2 ∆ u n 2 , ∆ u j + 1 2 = u j +1 − u j j + 1 j + 1 j − 1 j − 1 j 20 / 111
TVD schemes Theorem (Harten) If C j + 1 2 ≥ 0 , D j + 1 2 ≥ 0 , C j + 1 2 + D j + 1 2 ≤ 1 then the scheme is TVD. Unfortunately, TVD schemes need not be entropy consistent. So we have to additionally check this property. 21 / 111
Viscosity form We write the numerical flux in the viscosity form 2 = 1 2( f j + f j +1 ) − 1 g j + 1 2 λQ j + 1 2 ∆ u j + 1 2 where Q is called the viscosity coefficient 1 . If the finite volume scheme satisfies � � ∆ f j + 1 � � � � λ 2 � ≤ Q j + 1 2 ≤ 1 � � ∆ v j + 1 � 2 then it is TVD. If we further restrict the viscosity coefficient to � � ∆ f j + 1 � � 2 ≤ 1 � � λ 2 � ≤ Q j + 1 � � ∆ v j + 1 2 � 2 then the scheme is TVD and stable in maximum norm. 1 Such schemes are said to be essentially 3-point. 22 / 111
Example: Lax-Friedrichs scheme The numerical flux is given by g ( u, v ) = 1 2( f ( u ) + f ( v )) − 1 2 λ ( v − u ) for which the viscosity coefficient is Q LF 2 = 1 j + 1 This corresponds to the upper bound in the TVD condition. Hence the scheme is TVD provided the CFL condition � � ∆ f j + 1 � � � � λ max 2 � ≤ 1 � � ∆ v j + 1 j � 2 is satisfied. 23 / 111
Example: Murman-Roe scheme The numerical flux is given by g R ( u, v ) = 1 2( f ( u ) + f ( v )) − 1 2 | a ( u, v ) | ( v − u ) and its numerical viscosity coefficient is � � ∆ f j + 1 � � Q R � � 2 = λ | a ( v j , v j +1 ) | = λ 2 j + 1 � � ∆ v j + 1 � � 2 This corresponds to the lower bound in the TVD condition. The scheme is TVD provided the CFL condition � � ∆ f j + 1 � � � � λ max 2 � ≤ 1 � � ∆ v j + 1 j � 2 is satisfied. However, we know that this scheme admits entropy violating shocks. 24 / 111
Example: Lax-Wendroff scheme The numerical flux is given by g LW ( u, v ) = 1 2( f ( u ) + f ( v )) − 1 2 λa ( u, v )( f ( v ) − f ( u )) Hence its numerical viscosity coefficient is � � 2 ∆ f j + 1 2 ) 2 = λ 2 Q LW 2 = ( λa j + 1 2 j + 1 ∆ v j + 1 2 which does not satisfy the TVD condition . This also implies that it does not preserve monotonicity. 25 / 111
Entropy consistent scheme Consider a general finite volume scheme u n +1 = u n j − λ [ g n 2 − g n 2 ] , g j + 1 2 = g ( u j − k +1 , . . . , u j +1 ) j + 1 j − 1 j We say that the scheme is consistent with an entropy pair ( U, F ) if there exists a numerical entropy flux G j + 1 2 = G ( u j − k +1 , . . . , u j +1 ) which is consistent with the entropy flux F ( u ) in the sense G ( u, . . . , u ) = F ( u ) , ∀ u and the numerical solutions satisfy G n 2 − G n U ( u n +1 ) − U ( u n j ) j + 1 j − 1 j + 2 ≤ 0 ∆ t ∆ x This is a discrete approximation to the entropy inequality ∂U ∂t + ∂F ∂x ≤ 0 26 / 111
Definition (E-scheme) A consistent, conservative scheme is called an E-scheme if its numerical flux satisfies sign( v j +1 − v j )( g j + 1 2 − f ( u )) ≤ 0 for all u between v j and v j +1 . Remark Note that an E-scheme is essentially 3-point. Indeed letting v j +1 → v j with first v j +1 ≥ v j and then with v j +1 ≤ v j shows that g is essentially 3-point. Remark A 3-point monotone scheme is an E-scheme. Since g ( u, v ) is non-decreasing in u and non-increasing in v , we obtain g ( u, v ) ≤ g ( u, w ) ≤ g ( w, w ) = f ( w ) if u ≤ w ≤ v g ( u, v ) ≥ g ( w, v ) ≥ g ( w, w ) = f ( w ) if u ≥ w ≥ v and therefore sign( v − u )( g ( u, v ) − f ( w )) ≤ 0 , for all w between u and v In particular, the Godunov scheme is an E-scheme under CFL ≤ 1 . 27 / 111
Lemma Assume that CFL ≤ 1 . Then E-fluxes are characterized by � 2 ≤ g G g j + 1 if v j < v j +1 j + 1 2 2 ≥ g G g j + 1 if v j > v j +1 j + 1 2 where g G stands for Godunov numerical flux. Lemma Assume that CFL ≤ 1 . E-schemes are characterized by 0 ≤ Q G 2 ≤ Q j + 1 2 , ∀ j ∈ Z j + 1 28 / 111
Theorem (Viscous form and entropy condition) Assume that the CFL condition λ max | a ( u ) | ≤ 1 2 holds. An E-scheme whose coefficient of numerical viscosity satisfies 2 ≤ 1 Q G 2 ≤ Q j + 1 j + 1 2 is consistent with any entropy condition. The basic idea is to write any E-scheme as a convex combination of the Godunov scheme and a modified Lax-Friedrichs scheme, both of which satisfy entropy condition. 29 / 111
Godunov scheme The finite volume solution is made of piecewise constant states v ∆ ( x, t ) = v n j , x ∈ ( x j − 1 2 , x j + 1 2 ) , t ∈ [ t n , t n +1 ) which defines a Riemann problem at each cell face x = x j + 1 2 ∂w R + ∂ ∂xf ( w R ) = 0 , x ∈ ( x j , x j +1 ) , t ∈ [ t n , t n +1 ) ∂t � v n j , x < x j + 1 w R ( x, 0) = 2 v n j +1 , x > x j + 1 2 Under the CFL condition λ max | a ( u ) | ≤ 1 2 the solution at next time level is given by projecting the Riemann solution onto piecewise constant states � 0 � ∆ x 1 j ) d x + 1 2 v n +1 w R ( x/ ∆ t ; v n j − 1 , v n w R ( x/ ∆ t ; v n j , v n = j +1 ) d x j ∆ x ∆ x − ∆ x 0 2 30 / 111
Godunov scheme We can rewrite the above formula as a finite volume scheme with numerical flux g G 2 = f ( w R (0; v j , v j +1 )) j + 1 31 / 111
Approximate Riemann solver Let w ( x/t ; u l , u r ) be an approximation of the exact entropy solution w R ( x/t ; u l , u r ) of the Riemann problem ∂u ∂t + ∂ ∂xf ( u ) = 0 � u l , x < 0 u ( x, 0) = u r , x > 0 We will require that the approximate solution be consistent with the exact one in two respects Conservation : Integrate over rectangle ( − ∆ x 2 , + ∆ x 2 ) × (0 , ∆ t ) , and provided λ | a ( u ) | ≤ 1 2 , for all u between u l and u r we get � + ∆ x 1 w R ( x/ ∆ t ; u l , u r ) d x = 1 2 2( u l + u r ) + λ ( f ( u l ) − f ( u r )) ∆ x − ∆ x 2 32 / 111
Thus we require the approximate solution to satisfy � + ∆ x 1 w ( x/ ∆ t ; u l , u r ) d x = 1 2 2( u l + u r ) + λ ( f ( u l ) − f ( u r )) ∆ x − ∆ x 2 Entropy condition : Integrating the entropy inequality U t + F x ≤ 0 yields � + ∆ x 1 U ( w R ( x/ ∆ t ; u l , u r )) d x ≤ 1 2 2( U ( u l ) + U ( u r )) + λ ( F ( u l ) − F ( u r )) ∆ x − ∆ x 2 For consistency with the entropy condition, we require approximate solution to satisfy � + ∆ x 1 U ( w ( x/ ∆ t ; u l , u r )) d x ≤ 1 2 2( U ( u l ) + U ( u r )) + λ ( F ( u l ) − F ( u r )) ∆ x − ∆ x 2 Continuity : Finally, we require the solution to be continuous wrt the data w ( x/t ; u, u ) = u 33 / 111
Godunov-type scheme : With the help of such an approximate Riemann solution w , we define the Godunov-type scheme as follows � 0 ∆ x � 1 j ) d x + 1 2 v n +1 w ( x/ ∆ t ; v n j − 1 , v n w ( x/ ∆ t ; v n j , v n = j +1 ) d x j ∆ x ∆ x − ∆ x 0 2 Theorem Let w be the approximate Riemann solver which satisfies (1) conservation, (2) consistency with entropy condition for an entropy pair ( U, F ) and is (3) continuous. Then the Godunov-type scheme can be put in conservation form, is consistent with the conservation law and is consistent with the entropy condition associated with ( U, F ) under the CFL condition CFL ≤ 1 2 . 34 / 111
Roe scheme and its entropy modification The Roe scheme is an approximate Riemann solver with � u l , x/t < a ( u l , u r ) w ( x/t ; u l , u r ) = u r , x/t > a ( u l , u r ) with numerical flux g R ( u, v ) = 1 2( f ( u ) + f ( v )) − 1 2 | a ( u, v ) | ( v − u ) But we have seen that this admits entropy violating shocks. This solution has only shocks and hence there will be problem when the solution is a rarefaction. Harten and Hyman proposed the following approximate Riemann solver u l , x/t < a l u ∗ , w ( x/t ; u l , u r ) = a l < x/t < a r u r , x/t > a r Here the intermediate state u ∗ and a l , a r are yet to be specified (See GR1). 35 / 111
Entropy conservative schemes Let ( U, F ) be an entropy pair and define the entropy variable v = U ′ ( u ) and entropy potential ψ ( v ) = vf ( u ( v )) − F ( u ( v )) Consider the semi-discrete finite volume scheme g ∗ 2 − g ∗ d u j j + 1 j − 1 d t + 2 = 0 ∆ x Assume that the flux satisfies the condition (Tadmor) ( v j +1 − v j ) g ∗ 2 = ψ j +1 − ψ j j + 1 Now multiply semi-discrete scheme by v j = U ′ ( u j ) g ∗ 2 − g ∗ U ′ ( u j ) d u j j + 1 j − 1 d t + v j 2 = 0 ∆ x 36 / 111
Entropy conservative schemes We have the identity v j = 1 2( v j + v j +1 ) − 1 2 − 1 2( v j +1 − v j ) = { { v } } j + 1 2 � v � j + 1 2 so that 2 − 1 2 − 1 v j g ∗ 2 g ∗ 2 g ∗ 2 g ∗ 2 = { { v } } j + 1 2 � v � j + 1 2 = { { v } } j + 1 2 � ψ � j + 1 j + 1 j + 1 j + 1 j + 1 2 Similarly v j = 1 2( v j + v j − 1 ) − 1 2 + 1 2( v j − 1 − v j ) = { { v } } j − 1 2 � v � j − 1 2 so that 2 + 1 2 + 1 v j g ∗ 2 g ∗ 2 g ∗ 2 g ∗ 2 = { { v } } j − 1 2 � v � j − 1 2 = { { v } } j − 1 2 � ψ � j − 1 j − 1 j − 1 j − 1 j − 1 2 37 / 111
Entropy conservative schemes Hence 2 − 1 v j ( g ∗ 2 − g ∗ 2 g ∗ 2 ) = { { v } } j + 1 2( ψ j +1 − ψ j ) j + 1 j − 1 j + 1 2 − 1 2 g ∗ −{ { v } } j − 1 2( ψ j − ψ j − 1 ) j − 1 � � 2 − 1 2 g ∗ = { { v } } j + 1 2( ψ j +1 + ψ j ) j + 1 � � 2 − 1 2 g ∗ − { { v } } j − 1 2( ψ j + ψ j − 1 ) j − 1 � � 2 g ∗ = { { v } } j + 1 2 − { { ψ } } j + 1 j + 1 2 � � 2 g ∗ − { { v } } j − 1 2 − { { ψ } } j − 1 2 ) j − 1 Define the numerical entropy flux G ∗ 2 g ∗ 2 = { { v } } j + 1 2 − { { ψ } } j + 1 j + 1 j + 1 2 38 / 111
Entropy conservative schemes then we obtain the entropy equality G ∗ 2 − G ∗ d j + 1 j − 1 d tU ( u j ) + 2 = 0 ∆ x Now, let us check the consistency of the flux. The numerical flux is 2 = ψ j +1 − ψ j g ∗ j + 1 v j +1 − v j Firstly, this is a central flux since interchanging u j , u j +1 gives the same value. Secondly, if u j = u j +1 = u and correspondingly v j = v j +1 = v = U ′ ( u ) , the g ∗ ψ ′ ( v ) = f ( u ( v )) + vf ′ ( u ( v )) u ′ ( v ) − F ′ ( u ( v )) u ′ ( v ) = j + 1 2 f ( u ( v )) + [ U ′ ( u ) f ′ ( u ( v )) − F ′ ( u ( v ))] u ′ ( v ) = = f ( u ) The numerical entropy flux is also a central flux and if u j = u j +1 = u , we get G ∗ 2 = vf ( u ) − ψ ( v ) = F ( u ) j + 1 39 / 111
Entropy conservative schemes and hence is consistent with the entropy flux. Theorem If the numerical flux g ∗ 2 satisfies the condition j + 1 2 g ∗ � v � j + 1 2 = � ψ � j + 1 j + 1 2 then the semi-discrete finite volume scheme satisfies entropy conservation associated to the entropy U . 40 / 111
Example: Burger’s equation Take U ( u ) = 1 2 u 2 . Then � u � u ( s )( s ) d s = 1 U ′ ( s ) f ′ ( s ) d s = 3 s 3 F ( u ) = and ψ = ( u )(1 2 u 2 ) − 1 3 u 3 = 1 v = U ′ ( u ) = u, 6 u 3 The entropy conserving flux is 6 u 3 1 j +1 − 1 6 u 3 2 = ψ j +1 − ψ j = 1 j g ∗ 6( u 2 j + u j u j +1 + u 2 = j +1 ) j + 1 v j +1 − v j u j +1 − u j With this flux, the entropy equation is G ∗ 2 − G ∗ d t (1 d j + 1 j − 1 2 u 2 j ) + 2 = 0 ∆ x which implies that (assuming periodic bc) d 1 � 2 u 2 j = 0 d t j i.e., the energy is conserved. 41 / 111
Viscosity form: U ( u ) = 1 2 u 2 In this case v = U ′ ( u ) = u . Define the straight line path 2 ( ξ ) = 1 u j + 1 2( u j + u j +1 ) + ξ � u � j + 1 2 The entropy conservative flux is 1 1 � � 2 = ψ j +1 − ψ j 2 2 g ∗ ψ ′ ( u j + 1 = 2 ( ξ )) d ξ = f ( u j + 1 2 ( ξ )) d ξ j + 1 v j +1 − v j − 1 − 1 2 2 Equating this to the viscosity form � 1 2 ( ξ )) d ξ = 1 2( f j + f j +1 ) − 1 2 2 Q ∗ f ( u j + 1 2 � u � j + 1 j + 1 2 − 1 2 we can write the viscosity coefficient as (Exercise) 1 � 2 Q ∗ 2 ξf ′ ( u j + 1 2 = 2 ( ξ )) d ξ j + 1 − 1 2 42 / 111
Viscosity form: general case In this case v = U ′ ( u ) . Define the straight line path 2 ( ξ ) = 1 v j + 1 2( v j + v j +1 ) + ξ � v � j + 1 2 The entropy conservative flux is 1 1 � � 2 = ψ j +1 − ψ j 2 2 g ∗ ψ ′ ( v j + 1 = 2 ( ξ )) d ξ = f ( u ( v j + 1 2 ( ξ ))) d ξ j + 1 v j +1 − v j − 1 − 1 2 2 Equating this to the viscosity form � 1 2 ( ξ )) d ξ = 1 2( f j + f j +1 ) − 1 2 2 P ∗ f ( u j + 1 2 � v � j + 1 j + 1 2 − 1 2 we can write the viscosity coefficient as (Exercise) 1 � 2 P ∗ 2 ξf ′ ( u ( v j + 1 2 = 2 ( ξ ))) d ξ j + 1 − 1 2 43 / 111
Entropy consistent schemes We must produce entropy at shock waves, which means that we need an entropy inequality. Consider the semi-discrete finite volume scheme g j + 1 2 − g j − 1 d u j d t + 2 = 0 ∆ x Let us write the numerical flux in terms of viscosity form using jump in entropy variable 2 = 1 2( f j + f j +1 ) − 1 g j + 1 2 P j + 1 2 � v � j + 1 2 This can be re-written as 1 2( f j + f j +1 ) − 1 2 − 1 2 P ∗ 2 − P ∗ g j + 1 = 2 � v � j + 1 2( P j + 1 2 ) � v � j + 1 j + 1 j + 1 2 2 2 − 1 g ∗ 2 − P ∗ = 2( P j + 1 2 ) � v � j + 1 j + 1 j + 1 2 2 − 1 g ∗ 2 − P ∗ = 2 D j + 1 2 � v � j + 1 2 , D j + 1 2 = P j + 1 j + 1 j + 1 2 44 / 111
Entropy consistent schemes where g ∗ 2 is the entropy conservative flux. Let us derive the entropy equation. j + 1 2 ) − 1 2 + 1 2 ) = v j ( g ∗ 2 − g ∗ v j ( g j + 1 2 − g j − 1 2 v j D j + 1 2 � v � j + 1 2 v j D j − 1 2 � v � j − 1 j + 1 j − 1 2 Following steps as in the entropy conservative case 2 − 1 v j D j + 1 2 � v � j + 1 2 = { { v } } j + 1 2 D j + 1 2 � v � j + 1 2 � v � j + 1 2 D j + 1 2 � v � j + 1 2 2 + 1 v j D j − 1 2 � v � j − 1 2 = { { v } } j − 1 2 D j − 1 2 � v � j − 1 2 � v � j − 1 2 D j − 1 2 � v � j − 1 2 Hence � � 2 − 1 G ∗ v j ( g j + 1 2 − g j − 1 2 ) = 2 { { v } } j + 1 2 D j + 1 2 � v � j + 1 j + 1 2 � � 2 − 1 G ∗ − 2 { { v } } j − 1 2 D j − 1 2 � v � j − 1 j − 1 2 1 2 + 1 + 4 � v � j + 1 2 D j + 1 2 � v � j + 1 4 � v � j − 1 2 D j − 1 2 � v � j − 1 2 45 / 111
Entropy consistent schemes Define the numerical entropy flux 2 − 1 2 = G ∗ G j + 1 2 { { v } } j + 1 2 D j + 1 2 � v � j + 1 j + 1 2 then the entropy equation is G j + 1 2 − G j − 1 d = − 1 2 − 1 d tU ( u j ) + 2 4 � v � j + 1 2 D j + 1 2 � v � j + 1 4 � v � j − 1 2 D j − 1 2 � v � j − 1 2 ≤ 0 ∆ x where the inequality follows provided 2 ≥ P ∗ D j + 1 2 ≥ 0 ∀ j = ⇒ P j + 1 j + 1 2 i.e., the viscosity coefficient P must be larger than the viscosity coefficient P ∗ in the entropy conservative scheme . 46 / 111
Hyperbolic Systems Consider a system of hyperbolic conservation laws d ∂ u ∂ f j � ∂t + = 0 , f j = f j ( u ) ∂x j j =1 We say that ( U, F ) is an entropy pair for the above system if U ( u ) is strictly convex and F ′ j ( u ) = U ′ ( u ) f ′ j ( u ) , 1 ≤ j ≤ d Smooth solutions satisfy the additional equation d ∂U ∂F j � ∂t + = 0 ∂x j j =1 while in general we can only demand the inequality d ∂U ∂F j � ∂t + ≤ 0 ∂x j j =1 to hold in the sense of distributions which can be motivated from vanishing viscosity approach. The existence of such pairs is not guaranteed in general. 47 / 111
Strict convexity The function U ( u ) being strictly convex means that U ((1 − ξ ) u 1 + ξ u 2 ) < (1 − ξ ) U ( u 1 ) + ξU ( u 2 ) , 0 < ξ < 1 If U ( u ) is differentiable, this means the symmetric matrix U ′′ ( u ) ∂ 2 U [ U ′′ ( u )] ij = ∂u i ∂u j is positive definite, i.e., s ⊤ U ′′ ( u ) s > 0 , ∀ s � = 0 and hence all eigenvalues of U ′′ ( u ) are strictly positive. 48 / 111
Symmetric Hyperbolic Systems Suppose we do a change of variable from u to v d u ′ ( v ) ∂ v j ( u ( v )) u ′ ( v ) ∂ v � f ′ ∂t + = 0 ∂x j j =1 and we write this as d ∂ v ∂ v � A 0 ∂t + A j A 0 = 0 ∂x j j =1 where A 0 = u ′ ( v ) , A j = f ′ j ( u ) If A 0 is symmetric, positive definite and A j A 0 is symmetric, then we call this as a symmetric, hyperbolic form. The conservation law is said to be symmetrizable if such a change of variable exists. Theorem (Godunov, Mock), ([2], Thm 3.2, page 25) A necessary and sufficient condition for the conservation law to posses a strictly convex entropy U is that there exists a change of dependent variables u = u ( v ) that symmetrizes it. 49 / 111
Given a strictly convex function, the following theorem tell us when it will be an entropy function. Theorem ([2], Thm 3.1, page 24) Let U be a strictly convex function. A necessary and sufficient condition for U to be an entropy is that the matrices U ′′ ( u ) f ′ j ( u ) , 1 ≤ j ≤ d are symmetric. 50 / 111
Let U be a strictly convex entropy function. Define the entropy variables v ⊤ = U ′ ( u ) Since U ′′ ( u ) is positive definite, we can uniquely convert from u → v and v → u . If we transform the PDE from conserved variables to entropy variables, we obtain a symmetric, hyperbolic form. Given some conservation law, there is no general method to find an entropy function. Usually, for systems coming from Physics, we already know the existence of an entropy condition from the second law of Thermodynamics. For an entropy pair ( U, F ) , define the entropy potential ψ ( v ) = v · f ( u ( v )) − F ( u ( v )) 51 / 111
Entropy conservative scheme Consider the semi-discrete finite volume scheme f ∗ 2 − f ∗ ∂ u j j + 1 j − 1 ∂t + 2 = 0 ∆ x and assume that the numerical flux f ∗ 2 = f ∗ ( u j , u j +1 ) satisfies j + 1 2 · f ∗ � v � j + 1 2 = � ψ � j + 1 j + 1 2 Taking the dot product of the scheme with v j = U ( u j ) gives F ∗ 2 − F ∗ d j + 1 j − 1 d tU ( u j ) + 2 = 0 ∆ x where F ∗ 2 = F ∗ ( u j , u j +1 ) := { 2 · f ∗ { v } } j + 1 2 − { { ψ } } j + 1 j + 1 j + 1 2 is a consistent numerical entropy flux, i.e., F ∗ ( u , u ) = F ( u ) . 52 / 111
Higher order entropy conservative scheme The entropy conservative flux f ∗ ( u j , u j +1 ) leads to a second order scheme. This flux can be used as a building block to construct higher order schemes (LeFloch et al. [3]). Choose an integer p ≥ 1 and let α 1 , α 2 , . . . , α p be numbers which satisfy p p � � l 2 s − 1 α r = 0 , 2 rα r = 1 , s = 2 , 3 , . . . , p r =1 l =1 Define the numerical flux p r − 1 � � f ∗ , 2 p 2 = f ∗ , 2 p ( u j − p +1 , . . . , u j + p ) = f ∗ ( u j − s , u j − s + r ) α r j + 1 r =1 s =0 Then the semi-discrete FV scheme is 2 p ’th order accurate f ∗ , 2 p 2 − f ∗ , 2 p = ∂ f j + 1 j − 1 ∂x ( x j ) + O (∆ x ) 2 p 2 ∆ x 53 / 111
Higher order entropy conservative scheme and satisfies the entropy equation F ∗ , 2 p 2 − F ∗ , 2 p d U j j + 1 j − 1 d t + 2 = 0 ∆ x where p r − 1 � � F ∗ , 2 p F ∗ ( u j − s , u j − s + r ) 2 = α r j + 1 r =1 s =0 is a consistent entropy flux. Example : For p = 2 we get the fourth order accurate entropy conservative flux 2 = 4 3 f ∗ ( u j , u j +1 ) − 1 6 f ∗ ( u j − 1 , u j +1 ) − 1 f ∗ , 4 6 f ∗ ( u j , u j +2 ) j + 1 54 / 111
Entropy consistent scheme Consider the semi-discrete finite volume scheme 2 − f j − 1 f j + 1 ∂ u j ∂t + 2 = 0 ∆ x where the numerical flux is 2 − 1 2 = f ∗ 2 = D ⊤ f j + 1 2 D j + 1 2 � v � j + 1 2 , D j + 1 2 ≥ 0 j + 1 j + 1 Then we get the entropy inequality F j + 1 2 − F j − 1 d = − 1 2 − 1 4 � v � ⊤ 4 � v � ⊤ d tU ( u j ) + 2 2 D j − 1 2 � v � j − 1 2 D j + 1 2 � v � j + 1 2 ≤ 0 j − 1 j + 1 ∆ x where 2 − 1 2 = F ∗ } ⊤ F j + 1 2 { { v } 2 D j + 1 2 � v � j + 1 j + 1 j + 1 2 is a consistent numerical entropy flux. The flux f j + 1 2 is first order accurate since � v � j + 1 2 = O (∆ x ) . 55 / 111
Dissipation matrix Let us first write the entropy conservative flux in viscosity form 2 = 1 2( f j + f j +1 ) − 1 f ∗ 2 P ∗ 2 � v � j + 1 j + 1 j + 1 2 Define the symmetric matrix B ( v ) = ∂ ∂ v f ( u ( v )) then (Exercise) � 1 2 P ∗ 2 = 2 ξB ( v j + 1 2 ( ξ )) d ξ j + 1 − 1 2 where v j + 1 2 ( ξ ) is the linear path connecting v j , v j +1 2 + 1 v j + 1 2 ( ξ ) = { { v } } j + 1 2 ξ � v � j + 1 2 56 / 111
Dissipation matrix Let us write the entropy consistent flux also in viscosity form 2 = 1 2( f j + f j +1 ) − 1 f j + 1 2 P j + 1 2 � v � j + 1 2 This can be rewritten as 2 − 1 2 = f ∗ 2 − P ∗ f j + 1 2( P j + 1 ) � v � j + 1 j + 1 j + 1 2 2 � �� � D j + 1 2 Hence for entropy stability, we need to satisfy 2 ≥ P ∗ P j + 1 in the sense of SPD matrix ordering j + 1 2 Using a linearization � 1 2 2 = v ′ ( u ) j + 1 U ′′ ( u ( v j + 1 � v � j + 1 2 = H j + 1 2 � u � j + 1 2 , H j + 1 2 := 2 ( ξ ))) d ξ − 1 2 57 / 111
Dissipation matrix we can write the flux as 2 = 1 2( f j + f j +1 ) − 1 2 P j + 1 2 H j + 1 2 � u � j + 1 f j + 1 2 If we match the above flux to the Rusanov flux , 2 = λ m H − 1 P j + 1 2 H j + 1 2 = λ m I = ⇒ P j + 1 j + 1 2 where λ m is the maximum wavespeed at j + 1 2 . Tadmor shows that λ m H − 1 2 ≥ P ∗ (2) j + 1 j + 1 2 and hence the Rusanov flux is entropy consistent. We can write the Rusanov flux in terms of entropy variable jump as 2 = 1 2( f j + f j +1 ) − 1 f Rus 2 λ m H − 1 2 � v � j + 1 j + 1 j + 1 2 58 / 111
Dissipation matrix Remark : It is usually not possible to find explicit formula for H j + 1 2 . We can find an approximation to H j + 1 2 by performing a numerical integration, but we have to check that (2) is satisfied. Rusanov-type dissipation matrix : An easier scheme is to first approximate 2 ≈ u ′ ( v ) j + 1 u ′ ( v ) j + 1 2 = u ′ ( { � u � j + 1 2 � v � j + 1 2 , { v } } j + 1 2 ) and then add the Rusanov dissipation to the entropy conservative flux 2 − 1 2 = f ∗ 2 λ m u ′ ( v ) j + 1 2 � v � j + 1 f j + 1 j + 1 2 We discuss a similar approach later in the context of Euler equations, see also [4]. Roe-type dissipation matrix : Barth [5] shows that we can scale the eigenvectors in such a way that they satisfy the relation RR ⊤ = u ′ ( v ) 59 / 111
Dissipation matrix The Roe flux is given by 2 = 1 2( f j + f j +1 ) − 1 2 | R − 1 f Roe 2 R j + 1 2 | Λ j + 1 2 � u � j + 1 j + 1 j + 1 2 Converting jump in u to v 2 ≈ u ′ ( v ) j + 1 2 R ⊤ � u � j + 1 2 � v � j + 1 2 = R j + 1 2 � v � j + 1 j + 1 2 the dissipation in Roe flux can be written as 2 | R − 1 2 | R − 1 2 R ⊤ R j + 1 2 | Λ j + 1 2 � u � j + 1 = R j + 1 2 | Λ j + 1 2 R j + 1 2 � v � j + 1 j + 1 j + 1 j + 1 2 2 2 | R ⊤ = R j + 1 2 | Λ j + 1 2 � v � j + 1 j + 1 2 We have an SPD dissipation matrix 2 | R ⊤ D j + 1 2 = R j + 1 2 | Λ j + 1 j + 1 2 60 / 111
Dissipation matrix which can be used in combination with the entropy conservative flux 2 − 1 2 = f ∗ 2 | R ⊤ 2 R j + 1 2 | Λ j + 1 2 � v � j + 1 f j + 1 j + 1 j + 1 2 The above flux leads to an entropy consistent scheme (but the dissipation is probably not optimal). Remark : It is possible to carefully compute the dissipation matrix so that we get exact resolution of stationary contact waves, see [6]. 61 / 111
Higher order entropy consistent scheme To construct higher order scheme, we will follow the reconstruction approach. At each face j + 1 2 , we obtain a left and right value of entropy variables v L 2 , v R j + 1 j + 1 2 and define the numerical flux as 2 − 1 2 = f ∗ , 2 p 2 ( v R 2 − v L 2 D j + 1 2 ) f j + 1 j + 1 j + 1 j + 1 Note that we use the higher order entropy conservative flux, and use the reconstructed values to define the dissipative flux. Then we get the entropy equation F j + 1 2 − F j − 1 d − 1 4( v j − v j − 1 ) ⊤ D j − 1 2 ( v R 2 − v L d tU ( u j ) + 2 = 2 ) j − 1 j − 1 ∆ x − 1 4( v j +1 − v j ) ⊤ D j + 1 2 ( v R 2 − v L 2 ) j + 1 j + 1 where 2 − 1 2 = F ∗ , 2 p } ⊤ 2 ( v R 2 − v L F j + 1 2 { { v } 2 D j + 1 2 ) j + 1 j + 1 j + 1 j + 1 We do not know the sign of the right hand side. We have to design a reconstruction scheme that allows us to fix the sign of the terms on the right. 62 / 111
Sign preserving reconstruction Recall that the dissipation matrix can be written as 2 | R ⊤ D j + 1 2 = R j + 1 2 | Λ j + 1 j + 1 2 Define the new variables w k = R ⊤ 2 v k , k = . . . , j − 1 , j, j + 1 , . . . j + 1 We will use the ENO reconstruction scheme to obtain w L 2 , w R 2 . We can then j + 1 j + 1 compute 2 = R − 1 2 = R − 1 v L 2 w L v R 2 w R 2 , j + 1 j + 1 j + 1 j + 1 j + 1 j + 1 2 But it is not necessary to compute v L 2 , v R 2 since we can write the flux as j + 1 j + 1 2 − 1 2 = f ∗ , 2 p 2 | ( w R 2 − w L 2 R j + 1 2 | Λ j + 1 2 ) f j + 1 j + 1 j + 1 j + 1 63 / 111
Sign preserving reconstruction Fjordholm and Mishra [7] have shown that the ENO scheme preserves the sign of the jump, i.e., sign( w j +1 − w j ) = sign( w R 2 − w L 2 ) j + 1 j + 1 On the right of the entropy equation, we have terms like ( v j +1 − v j ) ⊤ D j + 1 2 ( v R 2 − v L 2 ) j + 1 j + 1 ( v j +1 − v j ) ⊤ R j + 1 2 | R ⊤ 2 ( v R 2 − v L = 2 | Λ j + 1 2 ) j + 1 j + 1 j + 1 ( R ⊤ 2 v j +1 − R ⊤ 2 v j ) ⊤ | Λ j + 1 2 | ( R ⊤ 2 v R 2 − R ⊤ 2 v L = 2 ) j + 1 j + 1 j + 1 j + 1 j + 1 j + 1 ( w j +1 − w j ) ⊤ | Λ j + 1 2 | ( w R 2 − w L = 2 ) j + 1 j + 1 ≥ 0 This shows that scheme satisfies entropy inequality. For second order scheme, the ENO scheme is same as reconstruction using minmod limiter. 64 / 111
Euler equations We will consider Euler equations in 1-D, which model friction-less fluids and can be written as ∂ u ∂t + ∂ f ∂x = 0 where ρ ρv , p + ρv 2 u = ρv f = ρe ( ρe + p ) v and ρ = mass density , ρv = momentum density , ρe = energy density and p is the pressure. The energy is made up of internal energy and kinetic energy ρe = ρε + 1 2 ρv 2 where ε = internal energy per unit mass. 65 / 111
Thermodynamic considerations We have four unknowns but only three equations. We need a way to relate the internal energy to the thermodynamics variable ρ, p, T where T is the absolute temperature. For a system in equilibrium, the thermodynamic variables satisfy an equation of state f ( ρ, p, T ) = 0 which may be solved for one of the variables, e.g., p = p ( ρ, T ) . The equation of state for an ideal gas can be written as p = ρRT where R is the gas constant . The specific heats at constant volume and constant pressure are defined as � ∂ε � ∂h � � τ = 1 C v = , C p = , where ρ, h = ε + pτ ∂T ∂T τ p 66 / 111
Thermodynamic considerations For an ideal gas ε = ε ( T ) , h = h ( T ) If we assume that C v , C p are independent of temperature (calorically and thermally perfect gas), which is a good assumption for air under normal conditions, we get ε = C v T, h = C p T, C p − C v = R Define the ratio of specific heats γ = C p R γR = ⇒ C v = γ − 1 , C p = C v γ − 1 Now we can close the model by using RT p ε = C v T = γ − 1 = ( γ − 1) ρ Some authors/books refer to this model as a polytropic ideal gas . 67 / 111
Thermodynamic considerations The second law of thermodynamics introduces a new state variable called entropy and denoted s . Under a quasi-static process T d s = d ε + p d τ For a polytropic ideal gas s = s 0 + C v ln( ε/ρ γ − 1 ) 68 / 111
Euler equation: Entropy function For our purpose, we can drop the constants and define the physical entropy function as s = ln( p/ρ γ ) From the Euler equations, we can derive an additional equation ∂t ( ρs ) + ∂ ∂ ∂x ( ρvs ) = 0 This motivates us to define the mathematical entropy function and entropy flux as U = − ρs F = − ρvs γ − 1 , (3) γ − 1 We can check that U ( u ) is a strictly convex function and F ′ ( u ) = U ′ ( u ) f ′ ( u ) is satisfied. In fact, we have many entropy functions of the form U = − ρη ( s ) F = − ρvη ( s ) γ − 1 , γ − 1 69 / 111
Euler equation: Entropy function where η ( s ) is any function that satisfies γη ′′ ( s ) < η ′ ( s ) The entropy condition for Navier-Stokes equations with Fourier law of heat condition requires that η ( s ) must be an affine function [8, 9]. Hence we will use the linear function η ( s ) = s from now onwards. For the entropy pair (3), the entropy variable is γ − s γ − 1 − βv 2 β = ρ 1 v = U ′ ( u ) = , 2 p = 2 βv 2 RT − 2 β and the entropy potential is ψ = v · f − F = ρv 70 / 111
Euler equation: Entropy conservative flux [6] An entropy conservative flux must satisfy the condition ( v j +1 − v j ) · f ∗ 2 = ψ j +1 − ψ j j + 1 This is one equation for the three components of the flux and is clearly under-determined. There are many possible solutions and we discuss one of them. Let f ∗ = [ f ∗ e ] ⊤ and let us choose ρ, v, β as independent variables. For any ρ , f ∗ m , f ∗ two grid functions a, b we have the identity 2 = a j +1 b j +1 − a j b j = a j + a j +1 ( b j +1 − b j ) + b j + b j +1 � ab � j + 1 ( a j +1 − a j ) 2 2 To keep the notation concise, we will drop the subscripts � ab � = { { a } } � b � + { { b } } � a � Using this relation, we can write � ψ � = � ρv � = { { ρ } } � v � + { { v } } � ρ � 71 / 111
Euler equation: Entropy conservative flux [6] The jump in entropy variables is − � s � { v 2 } } � v 2 � � v 1 � = γ − 1 − { } � β � − { { β } − � s � { v 2 } = γ − 1 − { } � β � − 2 { { β } }{ { v } } � v � � v 2 � = 2 { { β } } � v � + 2 { { v } } � β � � v 3 � = − 2 � β � We try to write � s � in terms of jumps in ρ, β ln( p/ρ γ ) = − ( γ − 1) ln ρ − ln β − ln 2 s = � s � = − ( γ − 1) � ln ρ � − � ln β � For a positive quantity, define the logarithmic average � a � � a � = � ln a � 72 / 111
Euler equation: Entropy conservative flux [6] Then � s � = − ( γ − 1) � ρ � � ρ � − � β � � β � and hence � � � v 1 � = � ρ � 1 { v 2 } � ρ � + ( γ − 1) � β � − { } � β � − 2 { { β } }{ { v } } � v � The condition for entropy conservative flux is � v 1 � f ∗ ρ + � v 2 � f ∗ m + � v 3 � f ∗ e = � ρv � Plugging in all the jump terms and collecting them together, we get f ∗ � � ρ } f ∗ } f ∗ � ρ � � ρ � + − 2 { { β } }{ { v } ρ + 2 { { β } � v � m �� � � 1 { v 2 } f ∗ } f ∗ m − 2 f ∗ + ( γ − 1) � β � − { } ρ + 2 { { v } � β � e = { { v } } � ρ � + { { ρ } } � v � 73 / 111
Euler equation: Entropy conservative flux [6] This equation must hold for all possible left and right states. Take the states such that � ρ � � = 0 , � v � = � β � = 0 which yields the mass flux f ∗ ρ = � ρ �{ { v } } Similarly, the other fluxes are obtained as m = { { ρ } } f ∗ } f ∗ } + { { v } ρ 2 { { β } and � � 2( γ − 1) � β � − 1 1 f ∗ { v 2 } f ∗ } f ∗ e = 2 { } ρ + { { v } m It is easy to check that these are consistent fluxes (Exercise). For some more numerical fluxes, see [10], [11]. 74 / 111
Kinetic energy consistency Kinetic energy per unit volume K = 1 2 ρv 2 satisfies the following equation � ∂v � 2 � � � � d p∂v ∂x d x − 4 p∂v K d x = µ d x ≤ ∂x d x (4) d t 3 ∂x Ω Ω Ω Ω Work done by pressure forces, absent in incompressible flows Irreversible destruction due to molecular diffusion Note: Convection contributes to only flux of KE across ∂ Ω . It does not change the total amount of KE inside the domain (except for boundary fluxes). The correct KE budget is important for simulation of turbulence since the KE cascade from large scales to small scales is a very important characteristic of turbulent flows. 75 / 111
KE preserving FVM ∂K − 1 2 v 2 ∂ρ ∂t + v ∂ ( ρv ) = ∂t ∂t � ∂v � 2 − ∂ ∂x ( p + ρv 2 / 2 − 4 3 µ∂v ∂x ) v + p∂v ∂x − 4 = 3 µ ∂x Centered numerical flux f ρ f ρ 0 f m } f ρ f j + 1 2 = = p + { ˜ { v } , g j + 1 2 = τ f e f e ˜ vτ − q j + 1 j + 1 j + 1 2 2 2 where 2 = 1 2 = 4 3 µv j +1 − v j 2 = − κT j +1 − T j { { v } } j + 1 2( v j + v j +1 ) , τ j + 1 , q j + 1 ∆ x ∆ x 76 / 111
KE preserving FVM Discrete KE equation � � 2 � � ∆ v j + 1 ∆ v j + 1 ∆ x d K j 2 − 4 � � = 2 p j + 1 ˜ 3 µ 2 ∆ x d t ∆ x ∆ x j j This is consistent with (4) and hence we refer to such a scheme as being KE consistent. Jameson’s KEP flux { { ρ } } { { u } } 2 = 1 } f ρ f j + 1 2 = { { p } } + { { v } , compare with f j + 1 2( f j + f j +1 ) } f ρ { { H } j + 1 2 We are free to choose ˜ p and the mass and energy fluxes in other ways. The entropy conservative flux also satisfies the KE consistency, since the momentum flux is of the correct form. 77 / 111
Euler equation: Rusanov-type dissipation The Rusanov flux is given by 2 = 1 2( f j + f j +1 ) − 1 f Rus 2 λ m � u � j + 1 j + 1 2 Derigs et al. [4] derive the exact relation ˜ � ρ � � ρ �{ { v } } E } 2 + ˜ ( ˜ � ρ �{ { v } } � ρ �{ { v } p E + ˜ p ) { { v } } � u � = H � v � , H = � E 2 � p 2 ˜ ( ˜ γ − 1 + ˜ 1 ˆ } 2 E E + ˜ p ) { { v } } + ˜ p { { v } � ρ � p = { { ρ } } p = � ρ � γ − 1 + 1 p ˆ ˜ { v 2 } ˜ } , ˆ 2 � β � , E = 2 � ρ �{ } 2 { { β } The matrix H is obviously symmetric and can also be shown to be positive definite. Then the entropy consistent flux can be taken as 2 − 1 2 = f ∗ f j + 1 2 λ m H j + 1 2 � v � j + 1 j + 1 2 78 / 111
Euler equation: Roe-type dissipation The Jacobian matrix is ρ ρv E H = h + 1 2 v 2 = ( ρe + p ) /ρ u ′ ( v ) = p + ρv 2 , ρv ρHv ( γ − 1) ρ + 1 γpE 4 ρv 4 E ρHv The eigenvectors are usually written as 1 1 1 a 2 = γp , R = v − a v v + a ρ 1 2 v 2 H − va H + va But this does not satisfy the condition RR ⊤ = u ′ ( v ) . Define the diagonal matrix � ρ � 2 γ , ( γ − 1) ρ , ρ S = diag γ 2 γ and the scaled eigenvector matrix ˜ 1 R = RS 2 79 / 111
Euler equation: Roe-type dissipation R ⊤ = u ′ ( v ) , so that the dissipation in Roe flux can be Then we can check that ˜ R ˜ written as − 1 R − 1 � u � ≈ − 1 R ⊤ � v � = − 1 R − 1 ˜ � u � ≈ ˜ R ˜ R ⊤ � v � , R | Λ | ˜ ˜ R | Λ | ˜ ˜ R ˜ R | Λ | ˜ ˜ R ⊤ � v � 2 2 2 and hence the entropy consistent flux can be taken as 2 − 1 2 = f ∗ ˜ 2 | ˜ R ⊤ f j + 1 R j + 1 2 | Λ j + 1 2 � v � j + 1 j + 1 j + 1 2 2 The dissipation matrix 2 = ˜ 2 | ˜ R ⊤ D j + 1 R j + 1 2 | Λ j + 1 j + 1 2 is symmetric and positive definite. 80 / 111
Entropy stable DG scheme Divide the domain into disjoint cells I j = ( x j − 1 2 , x j + 1 2 ) and inside each cell, we approximate the solution by a polynomial v h ( x ) = polynomial of degree k for x ∈ I j and v h ( x ) may be discontinuous at the cell faces x j + 1 2 . Note that we expand the entropy variable in terms of polynomials . The semi-discrete DG scheme is � � 2 · φ h ( x − 2 · φ h ( x + φ h · ∂ t u ( v h ) d x − f ( u ( v h )) · ∂ x φ h d x + f j + 1 2 ) − f j − 1 2 ) = 0 j + 1 j − 1 I j I j where the test functions φ h are also piecewise polynomials of degree k . The numerical flux is given by 2 − 1 2 = f ∗ 2 = v h ( x + 2 ) − v h ( x − 2 D j + 1 2 � v h � j + 1 2 , � v h � j + 1 2 ) f j + 1 j + 1 j + 1 j + 1 Let us take φ h = v h in the DG scheme. � � 2 · v h ( x − 2 · v h ( x + v h · ∂ t u ( v h ) d x − f ( u ( v h )) · ∂ x v h d x + f j + 1 2 ) − f j − 1 2 ) = 0 j + 1 j − 1 I j I j 81 / 111
Entropy stable DG scheme Since v h = U ′ ( u ( v h )) � � v h · ∂ t u ( v h ) d x = ∂ t U ( u ( v h )) d x I j I j The integrand in the second term can be written as f · ∂ x v = ∂ x ( f · v ) − v · ∂ x f ∂ x ( f · v ) − U ′ ( u ) f ′ ( u ) ∂ x u = ∂ x ( f · v ) − F ′ ( u ) ∂ x u = = ∂ x ( f · v ) − ∂ x F = ∂ x ( f · v − F ) = ∂ x ψ Hence the second term becomes � f ( u ( v h )) · ∂ x v h d x = ψ ( v h ( x + 2 )) − ψ ( v h ( x − − 2 )) j − 1 j + 1 I j 82 / 111
Entropy stable DG scheme The boundary flux terms can be written as 2 − 1 2 · v h ( x − f j + 1 2 ) = { { v h } } j + 1 2 · f j + 1 2 � v h � j + 1 2 · f j + 1 j + 1 2 where 2 = 1 2[ v h ( x − 2 ) + v h ( x + { { v h } } j + 1 2 )] j + 1 j + 1 Using the definition of the numerical flux 2 − 1 2 · v h ( x − 2 · f ∗ f j + 1 2 ) = { { v h } } j + 1 2 · f j + 1 2 � v h � j + 1 j + 1 j + 1 2 +1 4 � v h � ⊤ 2 D j + 1 2 � v h � j + 1 j + 1 2 2 − 1 2 )) + 1 2 ψ ( v h ( x + 2 ψ ( v h ( x − = { { v h } } j + 1 2 · f j + 1 2 )) j + 1 j + 1 +1 4 � v h � ⊤ 2 D j + 1 2 � v h � j + 1 j + 1 2 83 / 111
Entropy stable DG scheme Similarly 2 + 1 2 · v h ( x − 2 ) = { { v h } } j − 1 2 · f j − 1 2 � v h � j − 1 2 · f j − 1 f j − 1 j − 1 2 2 + 1 2 )) − 1 2 ψ ( v h ( x + 2 ψ ( v h ( x − = { { v h } } j − 1 2 · f j − 1 2 )) j − 1 j − 1 − 1 4 � v h � ⊤ 2 D j − 1 2 � v h � j − 1 j − 1 2 Adding all terms we get � � � 2 − 1 2 )) − 1 2 ψ ( v h ( x + 2 ψ ( v h ( x − ∂ t U ( v h ) d x + { { v h } } j + 1 2 · f j + 1 2 )) j + 1 j + 1 I j � � 2 − 1 2 )) − 1 2 ψ ( v h ( x + 2 ψ ( v h ( x − − { { v h } } j − 1 2 · f j − 1 2 )) j − 1 j − 1 = − 1 2 − 1 4 � v h � ⊤ 4 � v h � ⊤ 2 D j − 1 2 � v h � j − 1 2 D j + 1 2 � v h � j + 1 j − 1 j + 1 2 ≤ 0 84 / 111
Entropy stable DG scheme Defining the numerical entropy flux F j + 1 2 = { { v h } } j + 1 2 · f j + 1 2 − { { ψ ( v h ) } } j + 1 2 we get the cell entropy inequality � d U ( v h ) d x + F j + 1 2 − F j − 1 2 ≤ 0 d t I j Remark : We have shown entropy stability of DG scheme of arbitrary order of accuracy. Remark : We use integration by parts to prove the entropy stability. In practice we use quadrature to approximate the integrals but due to non-linearity, the quadratures are not exact. This inexact integration spoils the entropy stabililty property, which can lead to inaccurate solutions on coarse meshes. 85 / 111
DG scheme with SBP property The goal is to construct a DG scheme which satisfies the entropy condition even if the numerical integration is not exact. The two main ingredients are 1 summation by parts or SBP property 2 the use of entropy conservative/consistent fluxes For more details see [12], [13], [14], [15]. The DG scheme is based on nodal Lagrange basis functions. So we map each cell I i to the reference cell I = [ − 1 , +1] . We choose k + 1 Gauss-Lobatto-Legendre (GLL) quadrature nodes − 1 = ξ 0 < ξ 1 < . . . < ξ k = 1 Define the Lagrange basis polynomials each of degree k k ξ − ξ l � L j ( ξ ) = , L j ( ξ l ) = δ jl ξ j − ξ l l =0 ,l � = j 86 / 111
DG scheme with SBP property Define the continuous and discrete inner products � k � ( u, v ) := uv d ξ, ( u, v ) h := ω j u ( ξ j ) v ( ξ j ) I j =0 Recall that this GLL quadrature is exact for polynomials of degree upto 2 k − 1 . We want to approximate the first derivative of a grid function at the GLL nodes. The interpolation is given by k � u h ( ξ ) = L l ( ξ ) u l l =0 We differentiate the interpolant and evaluate at GLL nodes k k � � u ′ L ′ h ( ξ j ) = l ( ξ j ) u l = D jl u l , 0 ≤ j ≤ k l =0 l =0 87 / 111
DG scheme with SBP property where we defined the differentiation matrix D jl = L ′ l ( ξ j ) , 0 ≤ j, l ≤ k Define the mass and stiffness matrices S jl = ( L j , L ′ l ) h = ( L j , L ′ M jl = ( L j , L l ) h � = ( L j , L l ) , l ) The mass matrix is approximate and diagonal, k k � � M jl = ω r L j ( ξ r ) L l ( ξ r ) = ω r δ jr δ lr = ω j δ jl r =0 r =0 M = diag { ω 0 , ω 1 , . . . , ω k } also known as lumped mass matrix , while the stiffness matrix is exact. 88 / 111
Theorem (SBP property) Define B = diag {− 1 , 0 , . . . , +1 } Then MD + D ⊤ M = S + S ⊤ = B S = MD, which is a discrete analogue of integration by parts. Proof: By definition ( L j , L ′ S jl = l ) h k k � � ω r L j ( ξ r ) L ′ ω r δ jr L ′ = l ( ξ r ) = l ( ξ r ) r =0 r =0 ω j L ′ = l ( ξ j ) = M jj D jl k � = M jr D rl = ⇒ S = MD r =0 89 / 111
Next � � ( L j , L ′ L l , L ′ S jl + S lj = l ) + j � 1 ( L j L ′ l + L l L ′ = j ) d ξ − 1 � 1 d = d ξ ( L j L l ) d ξ − 1 = L j (1) L l (1) − L j ( − 1) L l ( − 1) = L j ( ξ k ) L l ( ξ k ) − L j ( ξ 0 ) L l ( ξ 0 ) = δ jk δ lk − δ j 0 δ l 0 = B jl and hence S + S ⊤ = B is proved. 90 / 111
Theorem For each j = 0 , 1 , . . . , k , we have − 1 j = 0 k k k � � � D jl = S jl = 0 , S lj = τ j = 0 1 ≤ j ≤ k − 1 l =0 l =0 l =0 +1 j = k Note : B = diag( τ 0 , τ 1 , . . . , τ k ) Proof: (1) The Lagrange polynomials form a partition of unity k � L l ( ξ ) = 1 ∀ ξ l =0 Hence k k � � L ′ D jl = l ( ξ j ) = 0 l =0 l =0 91 / 111
(2) Next, using previous Theorem and result from (1) k k � � S jl = ω j D jl = 0 l =0 l =0 (3) Finally, using S ⊤ = B − S k k k k � � � � S lj = B jl − S jl = B jl = B jj = τ j l =0 l =0 l =0 l =0 92 / 111
Semi-discrete DG scheme: Scalar case Find u h = � k l =0 u i l L l ( ξ ) such that for all test functions φ h of degree k � � ∂u h f ( u h ) ∂φ h 2 φ h ( x − 2 φ h ( x + ∂t φ h d x − ∂x d x + f i + 1 2 ) − f i − 1 2 ) = 0 i + 1 i − 1 I i I i The boundary flux is given by some numerical flux formula 2 = f ( u h ( x − 2 ) , u h ( x + k , u i +1 2 )) = f ( u i f i + 1 ) 0 i + 1 i + 1 and couples the solution in the two cells i and i + 1 . We map to reference cell I = [ − 1 , +1] x = x i ( ξ ) = 1 2 ) + 1 2( x i − 1 2 + x i + 1 2 ξ ∆ x i The DG scheme becomes � � ∆ x i ∂u h f ( u h ) ∂φ h ∂t φ h d ξ − ∂ξ d ξ + f i + 1 2 φ h (1) − f i − 1 2 φ h ( − 1) = 0 2 I I 93 / 111
Semi-discrete DG scheme: Scalar case We will approximate the integrals by GLL quadrature, and also approximate the flux k k � � f ( u i f i f ( u h ) ≈ f h ( ξ ) := j ) L j ( ξ ) = j L j ( ξ ) j =0 j =0 Taking the test functions φ h = L j ∆ x i d � � f h , L ′ d t ( u h , L j ) h − h + f i + 1 2 L j (1) − f i − 1 2 L j ( − 1) = 0 j 2 k k d u i ∆ x i � � � � L l , L ′ l f i d t ( L j , L l ) h − h + f i + 1 2 δ jk − f i − 1 2 δ j 0 = 0 l j 2 l =0 l =0 Define u i f i f i − 1 0 0 2 u i f i 0 1 1 . . . u i = f i = f i = ˆ . . . , , . . . u i f i 0 k − 1 k − 1 u i f i f i + 1 k k 2 94 / 111
Semi-discrete DG scheme: Scalar case Then we can write the DG scheme as 2 M d u i ∆ x i d t − S ⊤ f i = − B ˆ f i Using S ⊤ = B − S , we rewrite this as 2 M d u i ∆ x i d t + Sf i = B ( f i − ˆ f i ) or using S = MD d u i ∆ x i d t + Df i = M − 1 B ( f i − ˆ f i ) 2 Note that 2 ( Df i ) j = ∂f ∂x ( x j ) + O (∆ x ) k ∆ x i The last equation shows that we have a collocation method, i.e., we are collocating the PDE at the GLL nodes, together with a penalty term for the end nodes. 95 / 111
Semi-discrete DG scheme: System case We skip the detailed derivation in system case and directly give the equations. Consider a system of p equations and define the matrices and vectors M = M ⊗ I p , D = D ⊗ I p , S = S ⊗ I p , B = B ⊗ I p f i − 1 u i f i 2 0 0 0 u i f i u i = 1 f i = 1 f i = . ˆ , , . . . . . . . . 0 u i f i k k f i + 1 2 Then the semi-discrete scheme can be written as 2 M d u i ∆ x i d t − S ⊤ f i = − B ˆ f i d u i ∆ x i d t + Df i = M − 1 B ( f i − ˆ f i ) 2 96 / 111
Entropy conservative scheme We consider a single element, so we will omit the superscript i . The collocation scheme can be written as k ∆ x i d u j D jl f ( u l ) = τ j � ( f j − ˆ d t + f j ) 2 ω j l =0 We modify this scheme by introducing a symmetric flux f ∗ in the interior and boundary terms k ∆ x i d u j D jl f ∗ ( u j , u l ) = τ j � ( f j − ˆ f ∗ d t + 2 j ) (5) 2 ω j l =0 E.g., if f ∗ ( u j , u l ) = 1 2 ( f ( u j ) + f ( u l )) , we get k k k k � � � � D jl f ∗ ( u j , u l ) = 2 D jl f ( u j ) + D jl f ( u l ) = D jl f ( u l ) l =0 l =0 l =0 l =0 97 / 111
Theorem (Entropy conservative scheme) If f ∗ is consistent and symmetric, then (5) is conservative and high order accurate. If f ∗ is entropy conservative, then (5) is also entropy conservative. Proof: (1) Conservation property: The mean value changes as k k k k d ∆ x � � � � τ j ( f j − ˆ f ∗ S jl f ∗ ( u j , u l ) 2 ω j u j = j ) − 2 d t j =0 j =0 j =0 l =0 k k k � � � τ j ( f j − ˆ f ∗ ( S jl + S lj ) f ∗ ( u j , u l ) = j ) − j =0 j =0 l =0 k k k � � � τ j ( f j − ˆ f ∗ B jl f ∗ ( u j , u l ) = j ) − j =0 j =0 l =0 k k � � τ j ( f j − ˆ f ∗ = j ) − τ j f ( u j ) j =0 j =0 − ( f ∗ 2 − f ∗ = 2 ) i + 1 i − 1 The mean value in the cell changes only due to boundary fluxes. 98 / 111
(2) Accuracy property: Define f ∗ ( x, y ) = f ∗ ( u ( x ) , u ( y )) , f ( x ) = f ( u ( x )) Then f ∗ ( x, y ) is also symmetric and consistent f ∗ ( x, y ) = f ∗ ( y, x ) , f ∗ ( x, x ) = f ( x ) Hence ∂x ( x ) = ∂ f ∗ ∂x ( x, x ) + ∂ f ∗ ∂y ( x, x ) = 2 ∂ f ∗ ∂ f ∂y ( x, x ) The difference matrix D is exact for polynomials of degree upto k k 4 D jl f ∗ ( x ( ξ j ) , x ( ξ l )) = 2 ∂ f ∂y ( x ( ξ j ) , x ( ξ j ))+ O (∆ x k ) = ∂ f � ∂x ( x ( ξ j ))+ O (∆ x k ) ∆ x l =0 Hence the scheme has high order of accuracy. 99 / 111
(3) Entropy conservation: We compute the rate of change of total entropy in cell k k d ∆ x ∆ x d u j � � 2 ω j v ⊤ 2 ω j U j = j d t d t j =0 j =0 k k k � � � j ( f j − ˆ τ j v ⊤ f ∗ S jl v ⊤ j f ∗ ( u j , u l ) = j ) − 2 j =0 j =0 l =0 The second term can be written as k k � � ( B jl + S jl − S lj ) v ⊤ j f ∗ ( u j , u l ) j =0 l =0 k k k � � � τ j v ⊤ S jl ( v j − v l ) ⊤ f ∗ ( u j , u l ) = j f j + j =0 j =0 l =0 k k k � � � τ j v ⊤ = j f j + S jl ( ψ j − ψ l ) j =0 j =0 l =0 k � τ j ( v ⊤ = j f j − ψ j ) j =0 100 / 111
Recommend
More recommend
