Second order Implicit-Explicit Total Variation Diminishing schemes - PowerPoint PPT Presentation

Second order Implicit-Explicit Total Variation Diminishing schemes for the Euler system in the low Mach regime Victor Michel-Dansac Séminaire ANEDP , Lille, 22/03/2018 Giacomo Dimarco, Univ. of Ferrara, Italy Raphaël Loubère, Univ. of Bordeaux, CNRS, France Marie-Hélène Vignal, Univ. of Toulouse, France Funding: ANR MOONRISE, SHOM

Outline General context: multi-scale models and principle of AP schemes 1 An order 1 AP scheme for the Euler system in the low Mach limit 2 Second-order schemes in time 3 Second-order schemes in space and application to Euler 4 Work in progress and perspectives 5

General context 1/26 electron density (log scale) : t=0.1 Multiscale model M ε , depending on a parameter ε 0 2.4 − 1 2.3 − 2 2.2 − 3 In the (space-time) domain, ε can 2.1 − 4 2 − 5 1.9 − 6 1.8 be of same order as the reference scale; 1.7 − 7 1.6 − 8 1.5 − 9 be small compared to the reference scale; 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 take intermediate values. epsilon 0 10 −5 10 When ε is small: M 0 = lim ε → 0 M ε asympt. model −10 10 0 0.2 0.4 0.6 0.8 1 x Difficulties: Classical explicit schemes for M ε : they are stable and consistent if the mesh resolves all the scales of ε . = ⇒ very costly when ε → 0 Schemes for M 0 = ⇒ the mesh is independent of ε ➠ M 0 is not valid everywhere, it needs ε ≪ 1 But: ➠ the interface may be moving: how to locate it?

Principle of AP schemes 2/26 A possible solution: Asymptotic Preserving (AP) schemes Use the multi-scale model M ε even for small ε . Discretize M ε with a scheme preserving the limit ε → 0. ➠ The mesh is independent of ε : Asymptotic stability . ➠ Recovery of an approximate solution of M 0 when ε → 0: Asymptotic consistency . ➠ Asymptotically stable and consistent scheme = ⇒ Asymptotic preserving scheme (AP) . ([Jin, ’99] kinetic → hydro) The AP scheme may be used only to reconnect M ε and M 0 . ε 0 M 0 class. scheme M ε intermediate class. scheme ε = O ( 1 ) ε � 1 zone M ε AP scheme

Outline General context: multi-scale models and principle of AP schemes 1 An order 1 AP scheme for the Euler system in the low Mach limit 2 Second-order schemes in time 3 Second-order schemes in space and application to Euler 4 Work in progress and perspectives 5

The multi-scale model and its asymptotic limit 3/26 ➠ Isentropic Euler system in scaled variables: x ∈ Ω ⊂ IR d , t ≥ 0 � ∂ t ρ + ∇ · ( ρ u ) = 0 ( 1 ) ε (with p ( ρ ) = ρ γ ) ( M ε ) ∂ t ( ρ u )+ ∇ · ( ρ u ⊗ u )+ 1 ε ∇ p ( ρ ) = 0 ( 2 ) ε Parameter: ε = M 2 = | u | 2 / ( γ p ( ρ ) / ρ ) , M = Mach number Boundary and initial conditions: � ρ ( x , 0 ) = ρ 0 + ε ˜ ρ 0 ( x ) u · n = 0 on ∂ Ω and u ( x , 0 ) = u 0 ( x )+ ε ˜ u 0 ( x ) , with ∇ · u 0 = 0 The formal low Mach number limit ε → 0 : ( 2 ) ε ⇒ ∇ p ( ρ ) = 0 ⇒ ρ ( x , t ) = ρ ( t ) � ( 1 ) ε ⇒ | Ω | ρ ′ ( t )+ ρ ( t ) u · n = 0 ⇒ ρ ( t ) = ρ ( 0 ) = ρ 0 ⇒ ∇ · u = 0 ∂ Ω

The multi-scale model and its asymptotic limit 4/26 The asymptotic model: Rigorous limit [Klainerman & Majda, ’81]:  ρ = cst = ρ 0 ,   ( M 0 ) ρ 0 ∇ · u = 0 , ( 1 ) 0   ρ 0 ∂ t u + ρ 0 ∇ · ( u ⊗ u )+ ∇π 1 = 0 , ( 2 ) 0 � � where 1 π 1 = lim p ( ρ ) − p ( ρ 0 ) . ε ε → 0 − ∆ π 1 = ρ 0 ∇ 2 :( u ⊗ u ) Explicit eq. for π 1 : ∂ t ( 1 ) 0 − ∇ · ( 2 ) 0 = ⇒ The pressure wave equation from M ε : ∂ tt ρ − 1 ε ∆ p ( ρ ) = ∇ 2 : ( ρ u ⊗ u ) ∂ t ( 1 ) ε − ∇ · ( 2 ) ε = ⇒ ( 3 ) ε From a numerical point of view ⇒ conditional stability ∆ t ≤ √ ε ∆ x Explicit treatment of ( 3 ) ε = Implicit treatment of ( 3 ) ε = ⇒ uniform stability with respect to ε

An order 1 AP scheme in the low Mach limit 5/26 Time semi-discretization: [Degond, Deluzet, Sangam & Vignal, ’09], [Degond & Tang, ’11], [Chalons, Girardin & Kokh, ’15] If ρ n and u n are known at time t n :  ρ n + 1 − ρ n + ∇ · ( ρ u ) n + 1 = 0 ,  ( 1 )  (AS)  ∆ t ( ρ u ) n + 1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1  ε ∇ p ( ρ n + 1 ) = 0 .  ( 2 )  (AC) ∆ t ∇ p ( ρ n + 1 ) = 0 ε → 0 = ⇒ gives consistency at the limit implicit treatment of the pressure wave eq. = ⇒ uniform stability in ε ρ n + 1 − 2 ρ n + ρ n − 1 − 1 ε ∆ p ( ρ n + 1 ) = ∇ 2 : ( ρ U ⊗ U ) n ∆ t 2 ∇ · ( 2 ) inserted into ( 1 ) : gives an uncoupled formulation ρ n + 1 − ρ n + ∇ · ( ρ u ) n − ∆ t ε ∆ p ( ρ n + 1 ) − ∆ t ∇ 2 : ( ρ u ⊗ u ) n = 0 ∆ t

An order 1 AP scheme in the low Mach limit 6/26 The scheme proposed in [Dimarco, Loubère & Vignal, ’17]: ➠ Framework of IMEX (IMplicit-EXplicit) schemes: � ρ � ρ u � � � � 0 ∂ t + ∇ · + ∇ · = 0 . p ( ρ ) ρ u ρ u ⊗ u ε Id � �� � � �� � � �� � W F e ( W ) F i ( W ) ➠ The CFL condition comes from the explicit flux F e ( W ) : in 1D, we have ∆ x √ ε � � ∆ t AP ≤ ∆ x = ∆ x recall ∆ t class. ≤ j | , � λ n 2 | u n | u n γρ γ − 1 | i ± j where λ n j are the eigenvalues of the explicit Jacobian matrix DF e ( W n j ) . ➠ A linear stability analysis yields: if the implicit part is ⇒ L 2 stability; • centered = ⇒ TVD and L ∞ stability. • upwind = SSP Strong Stability Preserving, [Gottlieb, Shu & Tadmor, ’01]

Importance of the upwind implicit viscosity 7/26 To highlight the relevance of upwinding the implicit viscosity, we display the density ρ in the vicinity of a shock wave and a rarefaction wave ( ε = 0 . 99, 45 cells in the left panel, 150 cells in the right panel). 2 2 Exact Exact AP Di/=0 AP Di/=0 AP Di=0 AP Di=0 1.9 1.9 1.8 1.8 1.7 1.7 1.6 1.6 1.5 1.5 1.4 1.4 1.3 1.3 1.2 1.2 1.1 1.1 1 1 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.14 0.16 0.18 0.2 0.22 0.24 0.26 ⇒ L 2 stability and less diffusive × : centered implicit discretization = ⇒ L ∞ stability but more diffusive : upwind implicit discretization =

AP but diffusive results, 1D test case 8/26 ε = 0 . 99, 300 cells − 4 6 x 10 1.8 Class. scheme Class. scheme 1.6 1st-order AP 5 1st-order AP Class: 273 loops 1.4 Density 4 Time steps 1.2 7 CPU time 0.07 3 1 0.8 AP: 510 loops 2 0.6 1 CPU time 1.46 0.4 0 0.2 0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.2 0.4 0.6 0.8 1 x Time ε = 10 − 4 , 300 cells 1 − 3 Class. scheme 10 1st-order AP Density Class: 4036 loops Time steps Class. scheme 1 − 4 10 CPU time 0.82 1st-order AP AP: 57 loops − 5 CPU time 0.14 0.9999 10 0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.2 0.4 0.6 0.8 1 x Time

AP but diffusive results, 1D test case 9/26 − 4 10 1 Class. scheme Class. scheme AP scheme AP scheme ε = 10 − 4 Density of Time steps − 5 10 1 Underlying of the viscosity − 6 10 0.9999 0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.2 0.4 0.6 0.8 1 x Time ⇓ It is necessary to use high order schemes But they must respect the AP properties we also wish to retain the L ∞ stability

Outline General context: multi-scale models and principle of AP schemes 1 An order 1 AP scheme for the Euler system in the low Mach limit 2 Second-order schemes in time 3 Second-order schemes in space and application to Euler 4 Work in progress and perspectives 5

Principle of IMEX schemes 10/26 Bibliography for stiff source terms or ODE problems: Ascher, Boscarino, Cafflish, Dimarco, Filbet, Gottlieb, Happenhofer, Higueras, Jin, Koch, Kupka, LeFloch, Pareschi, Russo, Ruuth, Shu, Spiteri, Tadmor... ∂ t W + ∇ · F e ( W )+ ∇ · F i ( W ) = 0 . IMEX division: Step n: W n is known General principle: Quadrature formula with intermediate values: � t n + 1 � t n + 1 W ( t n + 1 ) = W ( t n ) − ∆ t ∇ · F e ( W ( t )) dt − ∆ t ∇ · F i ( W ( t )) dt t n t n � �� � � �� � s s W n + 1 = W n ˜ b j ∇ · F e ( W n , j ) − ∆ t b j ∇ · F i ( W n , j ) ∑ ∑ − ∆ t j = 1 j = 1 Quadratures exact on the constants: ∑ s j = 1 ˜ b j = ∑ s j = 1 b j = 1 Intermediate values at times t n , j = t n + c j ∆ t : � t n , j � c j W n , j ≈ W ( t n , j ) = W ( t n )+ ∂ t W ( t ) dt = W n +∆ t 0 ∂ t W ( t n + s ∆ t ) ds t n

Principle of IMEX schemes 11/26 Quadrature formula for intermediate values: i = 1 , ··· , s W n , j = W n − ∆ t ∑ a j , k ∇ · F e ( W n , k ) − ∆ t ∑ a j , k ∇ · F i ( W n , k ) , ˜ k < j k ≤ j s s ∑ ∑ Quadratures exact on the constants: ˜ a j , k = ˜ c j , a j , k = c j k = 1 k = 1 s s W n + 1 = W n − ∆ t b j ∇ · F e ( W n , j ) − ∆ t ˜ b j ∇ · F i ( W n , j ) ∑ ∑ j = 1 j = 1 Butcher tableaux: Explicit part Implicit part ··· ··· 0 0 0 0 c 1 a 1 , 1 0 0 ˜ ... ... c 2 a 2 , 1 0 0 c 2 a 2 , 1 a 2 , 2 0 . . . . . . ... ... ... ... . . . . . . . . . . . . ˜ ˜ ... ... c s a s , 1 a s , s − 1 0 c s a s , 1 a s , s − 1 a s , s ˜ ˜ ... ... ... ... b 1 b s b 1 b s Conditions for 2nd order: ∑ b j c j = ∑ b j ˜ c j = ∑ ˜ b j c j = ∑ ˜ b j ˜ c j = 1 / 2