u n Barrier to the use of explicit schemes 2 0. not actually used the limit We do not recover the incompressible limit! Note that we have 0 u 2 u n 1 u n 7/57 semi-discretization reads Assume ρ n and u n are known at time t n : the classical explicit time ρ n + 1 − ρ n + ∇ · ( ρ u ) n = 0 , ( 1 ) ∆ t ( ρ u ) n + 1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1 ε ∇ p ( ρ n ) = 0 . ( 2 ) ∆ t Set ρ n = ρ 0 and ∇ · u n = 0. Does ρ n + 1 = ρ 0 and ∇ · u n + 1 = 0 when ε → 0? ⇒ ρ n + 1 = ρ n − ∆ t ∇ · ( ρ u ) n = ρ 0 − ∆ t ρ 0 ∇ · u n = ρ 0 ( 1 ) = ⇒ ρ 0 u n + 1 = ρ 0 u n − ∆ t ρ 0 ∇ · ( u ⊗ u ) n − ∆ t 1 ( 2 ) = ε ∇ p ( ρ 0 )
Barrier to the use of explicit schemes 7/57 semi-discretization reads Assume ρ n and u n are known at time t n : the classical explicit time ρ n + 1 − ρ n + ∇ · ( ρ u ) n = 0 , ( 1 ) ∆ t ( ρ u ) n + 1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1 ε ∇ p ( ρ n ) = 0 . ( 2 ) ∆ t Set ρ n = ρ 0 and ∇ · u n = 0. Does ρ n + 1 = ρ 0 and ∇ · u n + 1 = 0 when ε → 0? ⇒ ρ n + 1 = ρ n − ∆ t ∇ · ( ρ u ) n = ρ 0 − ∆ t ρ 0 ∇ · u n = ρ 0 ( 1 ) = ⇒ ρ 0 u n + 1 = ρ 0 u n − ∆ t ρ 0 ∇ · ( u ⊗ u ) n − ∆ t 1 ( 2 ) = ε ∇ p ( ρ 0 ) ⇒ ∇ · u n + 1 = ∇ · u n − ∇ 2 : ( u ⊗ u ) n ̸ = 0 ∇ · ( 2 ) = � We do not recover the incompressible limit! Note that we have not actually used the limit ε → 0.
1 x t n 1 n t 0 Asymptotic-preserving time semi-discretization 1 n 1 0 n 1 0 u n n n 0 0 1 0 0 u n 1 u n 1 0 We recover the incompressible limit! The time semi-discretization is asymptotically consistent . 1 1 8/57 n n 0 1 n p 2 time semi-discretization reads Assume ρ n and u n are known at time t n : an Implicit-Explicit (IMEX) ρ n + 1 − ρ n + ∇ · ( ρ u ) n + 1 = 0 , ( 1 ) ∆ t ( ρ u ) n + 1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1 ε ∇ p ( ρ n + 1 ) = 0 . ( 2 ) ∆ t Set ρ n = ρ 0 and ∇ · u n = 0. Does ρ n + 1 = ρ 0 and ∇ · u n + 1 = 0 when ε → 0?
t n 1 n t 0 Asymptotic-preserving time semi-discretization 0 n 1 0 u n n 1 n 1 1 0 0 u n 1 u n 1 0 We recover the incompressible limit! The time semi-discretization is asymptotically consistent . 0 n 8/57 1 0 time semi-discretization reads Assume ρ n and u n are known at time t n : an Implicit-Explicit (IMEX) ρ n + 1 − ρ n + ∇ · ( ρ u ) n + 1 = 0 , ( 1 ) ∆ t ( ρ u ) n + 1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1 ε ∇ p ( ρ n + 1 ) = 0 . ( 2 ) ∆ t Set ρ n = ρ 0 and ∇ · u n = 0. Does ρ n + 1 = ρ 0 and ∇ · u n + 1 = 0 when ε → 0? ⇒ ∇ p ( ρ n + 1 ) = 0 = ⇒ ρ n + 1 ( x ) = ρ n + 1 ( 2 ) =
t 0 Asymptotic-preserving time semi-discretization 0 is asymptotically consistent . We recover the incompressible limit! The time semi-discretization 0 1 u n 1 u n 0 0 1 0 0 8/57 time semi-discretization reads Assume ρ n and u n are known at time t n : an Implicit-Explicit (IMEX) ρ n + 1 − ρ n + ∇ · ( ρ u ) n + 1 = 0 , ( 1 ) ∆ t ( ρ u ) n + 1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1 ε ∇ p ( ρ n + 1 ) = 0 . ( 2 ) ∆ t Set ρ n = ρ 0 and ∇ · u n = 0. Does ρ n + 1 = ρ 0 and ∇ · u n + 1 = 0 when ε → 0? ⇒ ∇ p ( ρ n + 1 ) = 0 = ⇒ ρ n + 1 ( x ) = ρ n + 1 ( 2 ) = ∫ ∫ = | Ω | ρ n − ∆ t ρ n + 1 u n + 1 · n = ⇒ ρ n + 1 = ρ n = ρ 0 ⇒ | Ω | ρ n + 1 ( 1 ε ) = Ω ∂Ω
Asymptotic-preserving time semi-discretization 0 is asymptotically consistent . 0 0 8/57 time semi-discretization reads Assume ρ n and u n are known at time t n : an Implicit-Explicit (IMEX) ρ n + 1 − ρ n + ∇ · ( ρ u ) n + 1 = 0 , ( 1 ) ∆ t ( ρ u ) n + 1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1 ε ∇ p ( ρ n + 1 ) = 0 . ( 2 ) ∆ t Set ρ n = ρ 0 and ∇ · u n = 0. Does ρ n + 1 = ρ 0 and ∇ · u n + 1 = 0 when ε → 0? ⇒ ∇ p ( ρ n + 1 ) = 0 = ⇒ ρ n + 1 ( x ) = ρ n + 1 ( 2 ) = ∫ ∫ = | Ω | ρ n − ∆ t ρ n + 1 u n + 1 · n = ⇒ ρ n + 1 = ρ n = ρ 0 ⇒ | Ω | ρ n + 1 ( 1 ε ) = Ω ∂Ω ⇒ ∇ · u n + 1 = 0 ⇒ ρ 0 = ρ 0 + ∆ t ρ 0 ∇ · u n + 1 = ( 1 ) = � We recover the incompressible limit! The time semi-discretization
1 n 1 n 2 n u n We get an uncoupled formulation of this AP time semi-discretization: u n u n Asymptotic-preserving time semi-discretization 2 into 1 p n 1 2 u Yes: the time semi-discretization is asymptotically stable . 1 n t 2 n t t p n 1 t 2 u 0 1 n 1 1 time semi-discretization reads Is the pressure wave equation discretized implicitly? 9/57 t n 1 2 Assume ρ n and u n are known at time t n : an Implicit-Explicit (IMEX) ρ n + 1 − ρ n + ∇ · ( ρ u ) n + 1 = 0 , ( 1 ) ∆ t ( ρ u ) n + 1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1 ε ∇ p ( ρ n + 1 ) = 0 . ( 2 ) ∆ t
We get an uncoupled formulation of this AP time semi-discretization: u n u n t 2 into 1 n 1 n t Asymptotic-preserving time semi-discretization Yes: the time semi-discretization is asymptotically stable . n 1 t 2 u 0 p 9/57 Is the pressure wave equation discretized implicitly? time semi-discretization reads Assume ρ n and u n are known at time t n : an Implicit-Explicit (IMEX) ρ n + 1 − ρ n + ∇ · ( ρ u ) n + 1 = 0 , ( 1 ) ∆ t ( ρ u ) n + 1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1 ε ∇ p ( ρ n + 1 ) = 0 . ( 2 ) ∆ t ( 1 ) n + 1 − ( 1 ) n ⇒ ρ n + 1 − 2 ρ n + ρ n − 1 ε∆ p ( ρ n + 1 ) = ∇ 2 : ( ρ u ⊗ u ) n − ∇ · ( 2 ) = − 1 ∆ t ∆ t 2
Asymptotic-preserving time semi-discretization Is the pressure wave equation discretized implicitly? Yes: the time semi-discretization is asymptotically stable . 9/57 time semi-discretization reads Assume ρ n and u n are known at time t n : an Implicit-Explicit (IMEX) ρ n + 1 − ρ n + ∇ · ( ρ u ) n + 1 = 0 , ( 1 ) ∆ t ( ρ u ) n + 1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1 ε ∇ p ( ρ n + 1 ) = 0 . ( 2 ) ∆ t ( 1 ) n + 1 − ( 1 ) n ⇒ ρ n + 1 − 2 ρ n + ρ n − 1 ε∆ p ( ρ n + 1 ) = ∇ 2 : ( ρ u ⊗ u ) n − ∇ · ( 2 ) = − 1 ∆ t ∆ t 2 We get an uncoupled formulation of this AP time semi-discretization: ⇒ ρ n + 1 − ρ n + ∇ · ( ρ u ) n − ∆ t ε ∆ p ( ρ n + 1 ) − ∆ t ∇ 2 : ( ρ u ⊗ u ) n = 0 . ∇ · ( 2 ) into ( 1 ) = ∆ t
IMEX formulation Rewrite this time semi-discretization under the following IMEX 1 0 where we have set: ( Implicit-Explicit ) form: 10/57 ρ n + 1 − ρ n + ∇ · ( ρ u ) n + 1 = 0 , ( 1 ) ∆ t ( ρ u ) n + 1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1 ε ∇ p ( ρ n + 1 ) = 0 . ( 2 ) ∆ t W n + 1 − W n + ∇ · F e ( W n ) + ∇ · F i ( W n + 1 ) = 0 , ∆ t � � � � � � ρ ρ u W = F e ( W ) = F i ( W ) = , , . ρ u ⊗ u ε p ( ρ ) I ρ u
CFL condition in 1D 2. The time step only depends on the explicit velocities: 11/57 W n + 1 − W n + ∂ x F e ( W n ) + ∂ x F i ( W n + 1 ) = 0 ∆ t 1. The characteristic velocities coming from F e and F i are as follows: F e : λ 1 e = 0 , λ 2 e = 2 u , � γρ γ − 1 λ ± F i : i = ± , ε and therefore both F e and F i are hyperbolic fmuxes. ∆ x √ ε � � ∆ t AP � ∆ x recall ∆ t class. � − − − → 2 | u | . | u √ ε ± � γρ γ − 1 | ε → 0 0 Next step : space discretization of ∂ x F e and ∂ x F i .
Finite volume space discretization W n j j j We consider the fjnite volume framework: outgoing fmux incoming fmux W n j W n We obtain the following fully discrete scheme: 2 x 1. the space domain is discretized with cells, 2. the cell values W n 12/57 2 j are updated with a discretization of the fmuxes. F ( W n j ) F ( W n j + 1 ) j − 1 , W n j , W n j − 1 j + 1 × × x j − 1 x j + 1 j + ∆ t + ∆ t � � W n + 1 F i ( W n + 1 j − 1 , W n + 1 ) − F i ( W n + 1 , W n + 1 � � = W n F e ( W n j − 1 , W n j ) − F e ( W n j , W n j + 1 ) j + 1 ) , ∆ x ∆ x where the numerical fmuxes F e and F i have to be determined.
13/57 2 centred and viscosity parts: R L 2 L ∞ -stable space discretization Between cells W L and W R , we assume that the fmuxes are made of F e ( W L , W R ) = ( F e ( W L ) + F e ( W R )) − D e ( W L , W R )( W R − W L ) , F i ( W L , W R ) = ( F i ( W L ) + F i ( W R )) − D i ( W L , W R )( W R − W L ) . • For D e , any solver can be used (Rusanov, HLL, …). • For D i , a linear stability analysis shows that: ◮ D i = 0 (i.e. a centred fmux) leads to a L 2 -stable scheme, ◮ L ∞ stability is achieved by taking (for the Rusanov fmux): � � γρ γ − 1 γρ γ − 1 D i ( W L , W R ) = max , . ε ε
Importance of the upwind implicit viscosity To highlight the relevance of upwinding the implicit viscosity, we 14/57 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 × : centered implicit discretization = ⇒ L 2 stability and less difgusive ⇒ L ∞ stability but more difgusive : upwind implicit discretization =
Asymptotic preservation but difgusive results Class.: CPU time 0.14 57 loops AP: CPU time 0.82 4036 loops 15/57 CPU time 1.46 Class.: 273 loops CPU time 0.07 AP: 510 loops − 4 6 x 10 ε = 0 . 99, 300 cells 1.8 Class. scheme Class. scheme 1.6 1st-order AP 5 1st-order AP 1.4 Density 4 Time steps 1.2 3 1 0.8 2 0.6 1 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 1 − 3 ε = 10 − 4 , 300 cells Class. scheme 10 1st-order AP Density Time steps Class. scheme 1 − 4 10 1st-order AP − 5 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
Asymptotic preservation but difgusive results Goal : build high-order schemes that • remain in the IMEX framework, to conserve the AP properties; 16/57 the viscosity Underlying of − 4 10 1 ε = 10 − 4 Class. scheme Class. scheme AP scheme AP scheme f Density Time steps − 5 10 1 − 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 for better results in all regimes of ε . • remain L ∞ -stable.
Multi-scale models and their numerical approximation A fjrst-order AP scheme for the isentropic Euler equations [Dimarco, Loubère, Vignal, ’17] High-order schemes in time and their failings [Dimarco, Loubère, M.-D., Vignal, ’18] [M.-D., Thomann, ’20] [M.-D., Thomann, in progress] Conclusion and future work “Second-order” L ∞ -stable schemes “High-order” L ∞ -stable schemes
Derivation of high-order IMEX schemes t s s W n t n t n s steps IMEX division : t n 17/57 t General principle : t n ∂ t W + ∇ · F e ( W ) + ∇ · F i ( W ) = 0 Assume that W n is known at time t n and compute the update W n + 1 by introducing s intermediate time steps. [ ] [ ] � | | | | | t n + 1 t n + 1 t n , k . . . . . . • Quadratures between t n and t n + 1 , introducing intermediate values: ∫ t n + 1 ∫ t n + 1 W ( t n + 1 ) = W ( t n ) − ∇ · F e ( W ( t )) dt − ∇ · F i ( W ( t )) dt ∑ ∑ ˜ W n + 1 b k ∇ · F e ( W n , k ) − ∆ t b k ∇ · F i ( W n , k ) = − ∆ t k = 1 k = 1
Derivation of high-order IMEX schemes t n t n (implicit part) t t n (explicit part) t t n t k W n 18/57 t n t n • Expression of the intermediate values at times t n , k = t n + c k ∆ t : ∫ t n , k W ( t n , k ) = W ( t n ) + ∂ t W ( t ) dt • Quadratures between t n and t n , k , ∀ k ∈ � 1 , s � : ∫ t n , k ∫ t n , k W ( t n , k ) = W ( t n ) − ∇ · F e ( W ( t )) dt − ∇ · F i ( W ( t )) dt k − 1 ∑ ∑ a k , l ∇ · F e ( W n , l ) − ∆ t a k , l ∇ · F i ( W n , l ) W n , k ˜ = − ∆ t l = 1 l = 1 ] [ | | | [ ] t n , k -1 t n , k | | | ] [ t n , k | | | t n , k
Derivation of high-order IMEX schemes . b s b 1 0 c s . The arbitrarily high-order IMEX time semi-discretization reads: ... 0 ... . . . . . c 1 0 0 ... b s b 1 c s . . . ... c 2 . . . . . . 0 . . 19/57 c 1 Butcher tableaux (Runge-Kutta time discretizations): Explicit part s k Implicit part 0 s 0 0 0 c 2 k − 1 ∑ ∑ ∀ k ∈ � 1 , s � , W n , k = W n − ∆ t a k , l ∇ · F e ( W n , l ) − ∆ t a k , l ∇ · F i ( W n , l ) , ˜ l = 1 l = 1 ∑ ∑ W n + 1 = W n − ∆ t ˜ b k ∇ · F e ( W n , k ) − ∆ t b k ∇ · F i ( W n , k ) . k = 1 k = 1 ˜ · · · · · · a 1 , 1 ˜ ˜ . . . . . . a 2 , 1 a 2 , 1 a 2 , 2 ˜ ˜ ˜ . . . . . . a s , 1 a s , s − 1 a s , 1 a s , s − 1 a s , s ˜ ˜ . . . . . . . . . . . .
IMEX schemes: a fjrst-order example 0 s s s s To yield a fjrst-order accurate scheme, the s -stage Butcher tableaux We consider the following Butcher tableaux: 0 c 2 c 1 implicit b 2 b 1 20/57 b 2 0 0 c 2 c 1 explicit b 1 ˜ a 1 , 1 ˜ ˜ a 2 , 1 a 2 , 2 a 2 , 1 ˜ ˜ have to satisfy the following consistency conditions: ∀ l ∈ � 1 , s � , ∑ ∑ ∑ ∑ ˜ ˜ c k = ˜ a k , l , c k = a k , l , b k = 1 , b k = 1 . l = 1 l = 1 k = 1 k = 1
IMEX schemes: a fjrst-order example 0 1 0 explicit 0 1 0 0 1 0 With the compatibility conditions, the tableaux become: 1 implicit 0 1 0 0 0 1 The simplest example ( Forward-Backward Euler ) is as follows: 21/57 0 c 2 b 2 b 2 explicit 0 c 2 0 0 0 b 2 implicit c 1 c 2 c 1 ˜ ˜ a 2 , 2 − c 2 a 2 , 2 1 − ˜ ˜ 1 − b 2
IMEX schemes: a fjrst-order example 0 1 0 0 0 1 0 implicit 1 0 1 0 0 1 0 explicit 0 1 22/57 W ( 1 ) W ( 2 ) W n + 1 W ( 1 ) W ( 2 ) W ( 1 ) W ( 2 )
IMEX schemes: a fjrst-order example 0 1 0 0 0 1 0 implicit 1 0 22/57 1 0 0 1 0 explicit 0 1 W ( 1 ) W ( 2 ) W n + 1 W ( 1 ) W ( 2 ) W ( 1 ) W ( 2 ) W ( 1 ) = W n − ∆ t 0 ∇ · F e ( W ( 1 ) ) − ∆ t 0 ∇ · F e ( W ( 2 ) ) − ∆ t 0 ∇ · F i ( W ( 1 ) ) − ∆ t 0 ∇ · F i ( W ( 2 ) )
IMEX schemes: a fjrst-order example 0 1 0 0 0 1 0 implicit 1 0 22/57 1 0 0 1 0 explicit 0 1 W ( 1 ) W ( 2 ) W n + 1 W ( 1 ) W ( 2 ) W ( 1 ) W ( 2 ) W ( 1 ) = W n − ∆ t 0 ∇ · F e ( W ( 1 ) ) − ∆ t 0 ∇ · F e ( W ( 2 ) ) − ∆ t 0 ∇ · F i ( W ( 1 ) ) − ∆ t 0 ∇ · F i ( W ( 2 ) ) W ( 2 ) = W n − ∆ t 1 ∇ · F e ( W ( 1 ) ) − ∆ t 0 ∇ · F e ( W ( 2 ) ) − ∆ t 0 ∇ · F i ( W ( 1 ) ) − ∆ t 1 ∇ · F i ( W ( 2 ) )
IMEX schemes: a fjrst-order example 0 1 0 0 0 1 0 implicit 1 0 1 0 1 0 explicit 0 1 0 22/57 W ( 1 ) W ( 2 ) W n + 1 W ( 1 ) W ( 2 ) W ( 1 ) W ( 2 ) W ( 1 ) = W n − ∆ t 0 ∇ · F e ( W ( 1 ) ) − ∆ t 0 ∇ · F e ( W ( 2 ) ) − ∆ t 0 ∇ · F i ( W ( 1 ) ) − ∆ t 0 ∇ · F i ( W ( 2 ) ) W ( 2 ) = W n − ∆ t 1 ∇ · F e ( W ( 1 ) ) − ∆ t 0 ∇ · F e ( W ( 2 ) ) − ∆ t 0 ∇ · F i ( W ( 1 ) ) − ∆ t 1 ∇ · F i ( W ( 2 ) ) W n + 1 = W n − ∆ t 1 ∇ · F e ( W ( 1 ) ) − ∆ t 0 ∇ · F e ( W ( 2 ) ) − ∆ t 0 ∇ · F i ( W ( 1 ) ) − ∆ t 1 ∇ · F i ( W ( 2 ) )
IMEX schemes: a fjrst-order example 0 the fjnal update is the same as the last intermediate time step: Remark : Since the weights are equal to the last row, i.e. 1 0 0 0 1 0 implicit 1 22/57 0 1 0 0 1 0 explicit 0 1 W ( 1 ) W ( 2 ) W n + 1 W ( 1 ) W ( 2 ) W ( 1 ) W ( 2 ) W ( 1 ) = W n W ( 2 ) = W n − ∆ t ∇ · F e ( W ( 1 ) ) − ∆ t ∇ · F i ( W ( 2 ) ) W n + 1 = W n − ∆ t ∇ · F e ( W ( 1 ) ) − ∆ t ∇ · F i ( W ( 2 ) ) ˜ a 2 , l and ˜ b l = ˜ b l = ˜ a 2 , l for l ∈ { 1 , 2 } , W n + 1 = W n − ∆ t ∇ · F e ( W n ) − ∆ t ∇ · F i ( W n + 1 ) .
IMEX schemes: a second-order example 0 1 explicit 0 1 To get a second-order scheme, we require additional compatibility 0 0 1 1 A simple example ( Implicit-Explicit Midpoint ) reads: implicit 0 1 2 0 0 0 1 2 0 2 23/57 s s s s conditions on the Butcher tableaux coeffjcients: ∑ ∑ ∑ ∑ ˜ ˜ b k ˜ b k ˜ b k c k = 1 2 , b k c k = 1 2 , c k = 1 2 , c k = 1 2 . k = 1 k = 1 k = 1 k = 1 W ( 1 ) = W n W ( 2 ) = W n − 1 2 ∆ t ∇ · F e ( W ( 1 ) ) − 1 2 ∆ t ∇ · F i ( W ( 2 ) ) � � � � 2 0 W n + 1 = W n − ∆ t ∇ · F e ( W ( 2 ) ) − ∆ t ∇ · F i ( W ( 2 ) )
IMEX schemes: a second-order example 0 1 explicit 0 To get a second-order scheme, we require additional compatibility 2 0 0 1 1 A simple example ( Implicit-Explicit Midpoint ) reads: implicit 0 1 2 0 0 0 1 2 0 1 23/57 s s s s conditions on the Butcher tableaux coeffjcients: ∑ ∑ ∑ ∑ ˜ ˜ b k ˜ b k ˜ b k c k = 1 2 , b k c k = 1 2 , c k = 1 2 , c k = 1 2 . k = 1 k = 1 k = 1 k = 1 W ∗ = W n − 1 2 ∆ t ∇ · F i ( W ∗ ) 2 ∆ t ∇ · F e ( W n ) − 1 � � � � W n + 1 = W n − ∆ t ∇ · F e ( W ∗ ) − ∆ t ∇ · F i ( W ∗ ) 2 0 This scheme is not L 2 -stable!
IMEX schemes: an L 2 -stable second-order example 0 0 0 1 1 implicit 0 1 0 1 0 0 0 1 1 1 1 1 0 0 1 2 1 1 0 explicit 0 24/57 0 0 0 0 √ To achieve L 2 stability, we can use the ARS(2,2,2) scheme with β = 1 − 2 : β β β β 1 − 1 − 2 ( 1 − β ) 2 ( 1 − β ) 2 β 2 β 1 − 1 − 2 β 2 β 2 ( 1 − β ) 2 ( 1 − β ) W ∗ = W n − β∆ t ∇ · F e ( W n ) − β∆ t ∇ · F i ( W ∗ ) � � W n + 1 = W n − 2 β∆ t ∇ · F e ( W ∗ ) 1 − 1 ∆ t ∇ · F e ( W n ) − 1 2 β � � 2 ( 1 − β ) ∆ t ∇ · F i ( W ∗ ) − ∆ t ∇ · F i ( W n + 1 ) − 1 − 2 ( 1 − β )
Application to the isentropic Euler equations 1 second-order in time fjrst-order in time exact solution 1 1 0 1 1 0 1 25/57 0 1 1 0 We display the density ρ ( x ) for a rarefaction-shock Riemann problem with several values of ε . ε = 10 − 1 ε = 10 − 2 1 . 1 1 . 01 1 . 05 1 . 005 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8 1 . 001 ε = 10 − 3 1 . 0001 ε = 10 − 4 1 . 0005 1 . 00005 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8 � We observe large oscillations when ε goes to 0!
Multi-scale models and their numerical approximation A fjrst-order AP scheme for the isentropic Euler equations [Dimarco, Loubère, Vignal, ’17] High-order schemes in time and their failings [Dimarco, Loubère, M.-D., Vignal, ’18] [M.-D., Thomann, ’20] [M.-D., Thomann, in progress] Conclusion and future work “Second-order” L ∞ -stable schemes “High-order” L ∞ -stable schemes
Better understand the oscillations: a model problem 0 The space discretization will be an upwind fmux. 0 1 x w n c i x w n c e t w n 1 w n and the fjrst-order AP time semi-discretization reads: cst w t w To better understand the oscillations, we consider a model problem: 0 x w 0 0 limit is obtained as follows: The 26/57 ∂ t w + c e ∂ x w + c i ε ∂ x w = 0 , c e > 0 , c i > 0 . It is a linear multi-scale advection equation, which we can write as: ∂ t w + ∂ x f e ( w ) + ∂ x f i ( w ) , f e ( w ) = c e w , f i ( w ) = c i ε w .
Better understand the oscillations: a model problem To better understand the oscillations, we consider a model problem: The space discretization will be an upwind fmux. and the fjrst-order AP time semi-discretization reads: 26/57 ∂ t w + c e ∂ x w + c i ε ∂ x w = 0 , c e > 0 , c i > 0 . It is a linear multi-scale advection equation, which we can write as: ∂ t w + ∂ x f e ( w ) + ∂ x f i ( w ) , f e ( w ) = c e w , f i ( w ) = c i ε w . The ε → 0 limit is obtained as follows: ε → 0 = ⇒ ∂ x w = 0 = ⇒ ∂ t w = 0 = ⇒ w = cst , w n + 1 − w n + c e ∂ x w n + c i ε ∂ x w n + 1 = 0 . ∆ t
t depends on , there are no implicit Runge-Kutta 27/57 no oscillations stability fails. Next step : Try to understand how the L 1 that are L -stable! discretizations of order unless Obstruction : A theorem [Gottlieb, Shu, Tadmor, ’01] states that, First idea : Derive a L -stable second-order IMEX discretization. decrease over time; at the discrete level, we wish to get j j and Better understand the oscillations: L ∞ stability failure The oscillations are measured by the L ∞ norm and the total variation: TV = � ∑ ∥ w n ∥ ∞ = max | w n j | TV ( w n ) = | w n j + 1 − w n j | . At the continuous level, the L ∞ norm and the total variation { L ∞ stability: ∥ w n + 1 ∥ ∞ � ∥ w n ∥ ∞ , ⇐ ⇒ total variation diminishing: TV ( w n + 1 ) � TV ( w n ) .
27/57 j Obstruction : A theorem [Gottlieb, Shu, Tadmor, ’01] states that, j no oscillations decrease over time; at the discrete level, we wish to get and Better understand the oscillations: L ∞ stability failure The oscillations are measured by the L ∞ norm and the total variation: TV = � ∑ ∥ w n ∥ ∞ = max | w n j | TV ( w n ) = | w n j + 1 − w n j | . At the continuous level, the L ∞ norm and the total variation { L ∞ stability: ∥ w n + 1 ∥ ∞ � ∥ w n ∥ ∞ , ⇐ ⇒ total variation diminishing: TV ( w n + 1 ) � TV ( w n ) . First idea : Derive a L ∞ -stable second-order IMEX discretization. unless ∆ t depends on ε , there are no implicit Runge-Kutta discretizations of order > 1 that are L ∞ -stable! Next step : Try to understand how the L ∞ stability fails.
Generic stability proof for an IMEX step j . c i c e 1 j j j w n j Rearranging the terms, we get: 28/57 w n j j An idealized IMEX step for the model problem, with weights α and β , is: w n + 1 w n + 1 − w n + 1 − w n j − w n j − 1 j − 1 + β c e + α c i = 0 . ∆ t ∆ x ε ∆ x ∆ t ∆ t � � � � w n + 1 w n + 1 − w n + 1 = w n j − β c e j − w n − α c i . j − 1 j − 1 ∆ x ε ∆ x ∆ t ∆ x and λ = β c e ∆ t Introducing µ ε = α c i ε ∆ x , the scheme reads: w n + 1 j − 1 − µ ε ( w n + 1 − w n + 1 = ( 1 − λ ) w n j + λ w n j − 1 ) . Next step : Derive conditions on µ ε and λ such that the IMEX step is L ∞ -stable, assuming periodic boundary conditions. ⇒ ∆ t � ∆ x ⇒ ∆ t � ∆ x ε Remark : λ � 1 ⇐ and µ ε � 1 ⇐ β α
1 if 1 1 since x 1 using the time update of the scheme 1 if 1 Generic stability proof for an IMEX step w n j 1 j 1 w n j 1 j 1 w n y x y since x 1 j 1 w n 1 j 1 j j 1 w n j 0 0 and j 1 w n j 1 1 w n j 1 1 j 1 w n j j j w n j 1 w n 1 j 0 1 and 0 i.e. w n j j j w n 1 j j j j w n j 1 j j 1 w n j j w n 1 1 y x y j w n j w n 1 29/57 max | w n j | = ( 1 − λ ) max | w n j | + λ max | w n j − 1 |
1 since x 1 using the time update of the scheme 1 if 1 Generic stability proof for an IMEX step w n j 1 j 1 w n j 1 j 1 w n y x y since x 1 j 1 1 j 1 j j 1 w n j 0 0 and j 1 w n j j 1 w n j 1 1 j 1 w n j w n w n 1 j j j j j j j 1 w n j w n 1 j 1 w n j w n j 1 j j 1 j w n 29/57 1 y 1 j x w n y max | w n j | = ( 1 − λ ) max | w n j | + λ max | w n j − 1 | = max | ( 1 − λ ) w n j | + max | λ w n j − 1 | if 1 − λ � 0 i.e. λ � 1 and λ � 0
1 since x 1 using the time update of the scheme 1 if 1 Generic stability proof for an IMEX step j 1 j 1 j 1 w n w n 1 1 1 j y x y since x 1 w n j j w n j 1 w n j 0 0 and j 1 j 1 j 1 w n j 1 1 j 1 w n j j w n x j j j j j 1 w n j w n j y j y j 1 w n 1 j j 1 w n 1 w n 1 j 29/57 max | w n j | = ( 1 − λ ) max | w n j | + λ max | w n j − 1 | = max | ( 1 − λ ) w n j | + max | λ w n j − 1 | if 1 − λ � 0 i.e. λ � 1 and λ � 0 = max j ( | ( 1 − λ ) w n j | + | λ w n j − 1 | )
1 using the time update of the scheme 1 if 1 Generic stability proof for an IMEX step w n 1 w n 1 j 1 j 1 1 j j w n 1 j y x y j w n 1 1 j 1 w n j 0 0 and j w n 1 j j 1 w n j 1 1 j since x 1 j 1 j j j j j j j 1 w n 29/57 j w n w n j 1 w n 1 j 1 j max | w n j | = ( 1 − λ ) max | w n j | + λ max | w n j − 1 | = max | ( 1 − λ ) w n j | + max | λ w n j − 1 | if 1 − λ � 0 i.e. λ � 1 and λ � 0 = max j ( | ( 1 − λ ) w n j | + | λ w n j − 1 | ) | ( 1 − λ ) w n j + λ w n j − 1 | since | x | + | y | � | x + y | � max
1 if 1 Generic stability proof for an IMEX step 1 j 1 w n 1 j 1 j w n w n j 1 w n 1 j y x j 1 since x 1 j 1 w n j 0 0 and j w n j j j 1 w n j 1 1 y 1 29/57 j j j j j j j j j 1 w n 1 j w n 1 j max | w n j | = ( 1 − λ ) max | w n j | + λ max | w n j − 1 | = max | ( 1 − λ ) w n j | + max | λ w n j − 1 | if 1 − λ � 0 i.e. λ � 1 and λ � 0 = max j ( | ( 1 − λ ) w n j | + | λ w n j − 1 | ) | ( 1 − λ ) w n j + λ w n j − 1 | since | x | + | y | � | x + y | � max | ( 1 + µ ε ) w n + 1 − µ ε w n + 1 = max j − 1 | using the time update of the scheme
1 if 1 Generic stability proof for an IMEX step j w n j j 1 w n 1 1 j j 1 w n j 1 w n 1 1 j j j 1 w n j 0 0 and 1 1 w n j j 1 w n j 1 29/57 j j j j j j j j j j max | w n j | = ( 1 − λ ) max | w n j | + λ max | w n j − 1 | = max | ( 1 − λ ) w n j | + max | λ w n j − 1 | if 1 − λ � 0 i.e. λ � 1 and λ � 0 = max j ( | ( 1 − λ ) w n j | + | λ w n j − 1 | ) | ( 1 − λ ) w n j + λ w n j − 1 | since | x | + | y | � | x + y | � max | ( 1 + µ ε ) w n + 1 − µ ε w n + 1 = max j − 1 | using the time update of the scheme || ( 1 + µ ε ) w n + 1 | − | µ ε w n + 1 � max j − 1 || since | x − y | � || x | − | y ||
1 if 1 Generic stability proof for an IMEX step j j j 1 w n 1 j j w n 1 j 1 1 w n j 1 j j w n 1 j 0 and 0 j w n 1 j j 29/57 j j j j j j j j max | w n j | = ( 1 − λ ) max | w n j | + λ max | w n j − 1 | = max | ( 1 − λ ) w n j | + max | λ w n j − 1 | if 1 − λ � 0 i.e. λ � 1 and λ � 0 = max j ( | ( 1 − λ ) w n j | + | λ w n j − 1 | ) | ( 1 − λ ) w n j + λ w n j − 1 | since | x | + | y | � | x + y | � max | ( 1 + µ ε ) w n + 1 − µ ε w n + 1 = max j − 1 | using the time update of the scheme || ( 1 + µ ε ) w n + 1 | − | µ ε w n + 1 � max j − 1 || since | x − y | � || x | − | y || j ( | ( 1 + µ ε ) w n + 1 | − | µ ε w n + 1 � max j − 1 | )
1 if 1 Generic stability proof for an IMEX step j j j j j j 1 j w n 1 j j w n 1 j 0 and 0 j w n 1 j j j 29/57 j j j j j j max | w n j | = ( 1 − λ ) max | w n j | + λ max | w n j − 1 | = max | ( 1 − λ ) w n j | + max | λ w n j − 1 | if 1 − λ � 0 i.e. λ � 1 and λ � 0 = max j ( | ( 1 − λ ) w n j | + | λ w n j − 1 | ) | ( 1 − λ ) w n j + λ w n j − 1 | since | x | + | y | � | x + y | � max | ( 1 + µ ε ) w n + 1 − µ ε w n + 1 = max j − 1 | using the time update of the scheme || ( 1 + µ ε ) w n + 1 | − | µ ε w n + 1 � max j − 1 || since | x − y | � || x | − | y || j ( | ( 1 + µ ε ) w n + 1 | − | µ ε w n + 1 � max j − 1 | ) | ( 1 + µ ε ) w n + 1 | µ ε w n + 1 � max | − max j − 1 |
Generic stability proof for an IMEX step j j j j j j j j j j j j w n 1 j j 29/57 j j j j j max | w n j | = ( 1 − λ ) max | w n j | + λ max | w n j − 1 | = max | ( 1 − λ ) w n j | + max | λ w n j − 1 | if 1 − λ � 0 i.e. λ � 1 and λ � 0 = max j ( | ( 1 − λ ) w n j | + | λ w n j − 1 | ) | ( 1 − λ ) w n j + λ w n j − 1 | since | x | + | y | � | x + y | � max | ( 1 + µ ε ) w n + 1 − µ ε w n + 1 = max j − 1 | using the time update of the scheme || ( 1 + µ ε ) w n + 1 | − | µ ε w n + 1 � max j − 1 || since | x − y | � || x | − | y || j ( | ( 1 + µ ε ) w n + 1 | − | µ ε w n + 1 � max j − 1 | ) | ( 1 + µ ε ) w n + 1 | µ ε w n + 1 � max | − max j − 1 | | w n + 1 | w n + 1 = ( 1 + µ ε ) max | − µ ε max j − 1 | if 1 + µ ε � 0 and µ ε � 0
Generic stability proof for an IMEX step j j j j j j j j j j j j j j 29/57 j j j j j max | w n j | = ( 1 − λ ) max | w n j | + λ max | w n j − 1 | = max | ( 1 − λ ) w n j | + max | λ w n j − 1 | if 1 − λ � 0 i.e. λ � 1 and λ � 0 = max j ( | ( 1 − λ ) w n j | + | λ w n j − 1 | ) | ( 1 − λ ) w n j + λ w n j − 1 | since | x | + | y | � | x + y | � max | ( 1 + µ ε ) w n + 1 − µ ε w n + 1 = max j − 1 | using the time update of the scheme || ( 1 + µ ε ) w n + 1 | − | µ ε w n + 1 � max j − 1 || since | x − y | � || x | − | y || j ( | ( 1 + µ ε ) w n + 1 | − | µ ε w n + 1 � max j − 1 | ) | ( 1 + µ ε ) w n + 1 | µ ε w n + 1 � max | − max j − 1 | | w n + 1 | w n + 1 = ( 1 + µ ε ) max | − µ ε max j − 1 | if 1 + µ ε � 0 and µ ε � 0 | w n + 1 = max |
1 is 1 . 1 from the fjrst step into the second one: 2 w j 1 w n 1 1 2 1 1 j Instability of the ARS(2,2,2) scheme w n 1 w j w j 2. Inject 1 0 and 0 1. The fjrst step is L -stable if j j 1 2 1 stability! loss of L 0 j 1, so the coeffjcient of w n We have 1 j w n 1 j 1 w n 2 1 2 1 1 w j 2 1 30/57 c i j 2 √ 2 , λ = ∆ t ∆ x c e and µ ε = ∆ t Recall the ARS(2,2,2) discretization, with β = 1 − ∆ x ε : w ∗ j − 1 ) − βµ ε ( w ∗ j − w ∗ j = w n j − βλ ( w n j − w n j − 1 ) , � � w n + 1 2 βλ ( w ∗ j − w ∗ = w n j − 1 − 1 λ ( w n j − w n j − 1 ) − 1 j − 1 ) 2 β � � 2 ( 1 − β ) µ ε ( w ∗ j − w ∗ µ ε ( w n + 1 − w n + 1 − j − 1 ) − 1 − j − 1 ) . 2 ( 1 − β )
1 is 1 from the fjrst step into the second one: 2 w j 1 2 1 j w n 1 1 2 Instability of the ARS(2,2,2) scheme 1 1 1 w n w j w j 2. Inject j j w j 2 j stability! loss of L 0 j 1, so the coeffjcient of w n We have 1 1 1 w n j 1 w n 2 1 2 1 1 1 30/57 j 2 c i √ 2 , λ = ∆ t ∆ x c e and µ ε = ∆ t Recall the ARS(2,2,2) discretization, with β = 1 − ∆ x ε : w ∗ j − 1 − βµ ε ( w ∗ j − w ∗ j = ( 1 − βλ ) w n j + βλ w n j − 1 ) , � � w n + 1 2 βλ ( w ∗ j − w ∗ = w n j − 1 − 1 λ ( w n j − w n j − 1 ) − 1 j − 1 ) 2 β � � 2 ( 1 − β ) µ ε ( w ∗ j − w ∗ µ ε ( w n + 1 − w n + 1 − j − 1 ) − 1 − j − 1 ) . 2 ( 1 − β ) ⇒ λ � 1 1. The fjrst step is L ∞ -stable if βµ ε > 0 and 0 � βλ � 1 = β .
1 is Instability of the ARS(2,2,2) scheme 1 1 j j 1 w n We have j j 1, so the coeffjcient of w n j 0 loss of L stability! 1 30/57 2 c i √ 2 , λ = ∆ t ∆ x c e and µ ε = ∆ t Recall the ARS(2,2,2) discretization, with β = 1 − ∆ x ε : w ∗ j − 1 − βµ ε ( w ∗ j − w ∗ j = ( 1 − βλ ) w n j + βλ w n j − 1 ) , � � w n + 1 2 βλ ( w ∗ j − w ∗ = w n j − 1 − 1 λ ( w n j − w n j − 1 ) − 1 j − 1 ) 2 β � � 2 ( 1 − β ) µ ε ( w ∗ j − w ∗ µ ε ( w n + 1 − w n + 1 − j − 1 ) − 1 − j − 1 ) . 2 ( 1 − β ) ⇒ λ � 1 1. The fjrst step is L ∞ -stable if βµ ε > 0 and 0 � βλ � 1 = β . 2. Inject µ ε ( w ∗ j − w ∗ j − 1 ) from the fjrst step into the second one: � � 2 β ( 1 − β ) + λ 1 − β j + β − 1 w n + 1 = 1 − λ w n j − 1 β β � � 2 β ( 1 − β ) − λ j + λ j − 1 + 1 − 2 β 2 ( 1 − β ) µ ε ( w n + 1 − w n + 1 w ∗ 2 β w ∗ + j − 1 ) . 2 β
Instability of the ARS(2,2,2) scheme j j 1 w n 1 j j 1 1 30/57 2 c i √ 2 , λ = ∆ t ∆ x c e and µ ε = ∆ t Recall the ARS(2,2,2) discretization, with β = 1 − ∆ x ε : w ∗ j − 1 − βµ ε ( w ∗ j − w ∗ j = ( 1 − βλ ) w n j + βλ w n j − 1 ) , � � w n + 1 2 βλ ( w ∗ j − w ∗ = w n j − 1 − 1 λ ( w n j − w n j − 1 ) − 1 j − 1 ) 2 β � � 2 ( 1 − β ) µ ε ( w ∗ j − w ∗ µ ε ( w n + 1 − w n + 1 − j − 1 ) − 1 − j − 1 ) . 2 ( 1 − β ) ⇒ λ � 1 1. The fjrst step is L ∞ -stable if βµ ε > 0 and 0 � βλ � 1 = β . 2. Inject µ ε ( w ∗ j − w ∗ j − 1 ) from the fjrst step into the second one: � � 2 β ( 1 − β ) + λ 1 − β j + β − 1 w n + 1 = 1 − λ w n j − 1 β β � � 2 β ( 1 − β ) − λ j + λ j − 1 + 1 − 2 β 2 ( 1 − β ) µ ε ( w n + 1 − w n + 1 w ∗ 2 β w ∗ + j − 1 ) . 2 β We have β < 1, so the coeffjcient of w n j − 1 is � 0 � loss of L ∞ stability!
A convex combination 1 w j 1 1 2 1 w j w j 1 1 1 2 1 w n 1 j w n j 2 w n 1 j 1 w n j 1 1 1 1 j w n j w n 1 w j 1 31/57 1 The full scheme reads as follows: w j w n j w n j w n j 1 w j w j 1 w n 2 j 1 j w n w n 1 1 j w n j The second step w n + 1 , o2 of the ARS(2,2,2) scheme is not L ∞ -stable. To address this issue, we introduce a convex combination between the second step and the fjrst-order scheme w n + 1 , o1 : { order 1, difgusive and stable, if θ = 0 , w n + 1 = θ w n + 1 , o2 + ( 1 − θ ) w n + 1 , o1 ⇒ = order 2, accurate and unstable, if θ = 1 .
A convex combination The full scheme reads as follows: j j 1 1 j 31/57 The second step w n + 1 , o2 of the ARS(2,2,2) scheme is not L ∞ -stable. To address this issue, we introduce a convex combination between the second step and the fjrst-order scheme w n + 1 , o1 : { order 1, difgusive and stable, if θ = 0 , w n + 1 = θ w n + 1 , o2 + ( 1 − θ ) w n + 1 , o1 ⇒ = order 2, accurate and unstable, if θ = 1 . w ∗ j − 1 ) − βµ ε ( w ∗ j − w ∗ j = w n j − βλ ( w n j − w n j − 1 ) , � � w n + 1 2 βλ ( w ∗ j − w ∗ = w n j − θ 1 − 1 λ ( w n j − w n j − 1 ) − θ 1 j − 1 ) 2 β � � 2 ( 1 − β ) µ ε ( w ∗ j − w ∗ µ ε ( w n + 1 − w n + 1 − θ j − 1 ) − θ 1 − j − 1 ) 2 ( 1 − β ) j − 1 ) − ( 1 − θ ) µ ε ( w n + 1 − w n + 1 − ( 1 − θ ) λ ( w n j − w n j − 1 ) .
32/57 Corollary : If 1 2 opt 1 1 1 1 0 , the conditions are relaxed: 1 2 opt c e Applying the TVD proof to the convex combination scheme yields: 2 Determination of θ Next step : Determine the largest θ that ensures the L ∞ stability. Theorem : The convex combination scheme is L ∞ -stable if: 0 < β < 1 , λ � 1 , θ � 2 β ( 1 − β ) . √ Remark : We are no longer constrained by the L 2 condition β = 1 − 2 . ⇒ ∆ t ⇒ ∆ t � ∆ x Remark : λ � 1 ⇐ ∆ xc e � 1 ⇐ is a CFL condition
32/57 2 Applying the TVD proof to the convex combination scheme yields: 1 c e Determination of θ Next step : Determine the largest θ that ensures the L ∞ stability. Theorem : The convex combination scheme is L ∞ -stable if: 0 < β < 1 , λ � 1 , θ � 2 β ( 1 − β ) . √ Remark : We are no longer constrained by the L 2 condition β = 1 − 2 . ⇒ ∆ t ⇒ ∆ t � ∆ x Remark : λ � 1 ⇐ ∆ xc e � 1 ⇐ is a CFL condition Corollary : If θ = θ opt = 2 β ( 1 − β ) , the conditions are relaxed: � 1 � θ = θ opt = 2 β ( 1 − β ) . 0 < β < 1 , λ � min β, , 1 − β
Numerical results with ARS(2,2,2) 1 0005 w 10 1 0 0 25 0 5 0 75 1 1 1 001 1 1 x w 10 3 exact solution fjrst-order second-order convex combination Question : Do we have to use the convex combination at every time step ? x 1 05 33/57 1 1 0 75 0 5 0 25 0 2 otherwise. 1 1 We check the new scheme with the ARS(2,2,2) value of β : √ β = 1 − λ = 1 − β ≃ 1 . 414, θ = 2 β ( 1 − β ) ≃ 0 . 414. 2 , For x ∈ ( 0 , 1 ) , we take a step function as initial condition: { 1 + ε if x ∈ ( 0 . 25 , 0 . 75 ) , w 0 ( x ) = we set c e = c i = 1, t end = ε/ ( 1 + ε ) , and we prescribe periodic BC.
Numerical results with ARS(2,2,2) 0 Question : Do we have to use the convex combination at every time step ? convex combination second-order fjrst-order exact solution w x 1 1 0 w x 1 1 otherwise. 1 2 33/57 1 We check the new scheme with the ARS(2,2,2) value of β : √ β = 1 − λ = 1 − β ≃ 1 . 414, θ = 2 β ( 1 − β ) ≃ 0 . 414. 2 , For x ∈ ( 0 , 1 ) , we take a step function as initial condition: { 1 + ε if x ∈ ( 0 . 25 , 0 . 75 ) , w 0 ( x ) = we set c e = c i = 1, t end = ε/ ( 1 + ε ) , and we prescribe periodic BC. 1 . 1 1 . 001 1 . 05 1 . 0005 ε = 10 − 1 ε = 10 − 3 0 . 25 0 . 5 0 . 75 0 . 25 0 . 5 0 . 75
MOOD procedure 2 nd -order scheme combination the convex compute Remark : The convex combination removes the oscillations but the detection of oscillations using the approximation w c compute a candidate w n (Multidimensional Optimal Order Detection). approximation remains difgusive. 34/57 � Increase the resolution using the MOOD framework w n + 1 = w c no w n + 1 w n + 1 with yes Figure above : MOOD method for one time step (compute w n + 1 from w n )
1 with the 2 nd -order scheme and check Detection of oscillations for the model problem 1 bounds are not violated. The 2 nd -order solution is used as soon the initial Consequence : bounds of the initial condition. We know that the numerical solution should not leave the Why? . w 0 whether w n How to detect oscillations? 2. Second idea : Compute w n tion is smooth. Consequence : The 2 nd -order scheme is never used, unless the solu- 1. First idea : 35/57 Compute w n + 1 with the 2 nd -order scheme and check whether ∥ w n + 1 ∥ ∞ � ∥ w n ∥ ∞ . Problem : If the w n is not smooth enough, the 2 nd -order solution w n + 1 will always violate the L ∞ stability.
Detection of oscillations for the model problem How to detect oscillations? bounds are not violated. The 2 nd -order solution is used as soon the initial Consequence : bounds of the initial condition. We know that the numerical solution should not leave the Why? 35/57 tion is smooth. Consequence : The 2 nd -order scheme is never used, unless the solu- 1. First idea : Compute w n + 1 with the 2 nd -order scheme and check whether ∥ w n + 1 ∥ ∞ � ∥ w n ∥ ∞ . Problem : If the w n is not smooth enough, the 2 nd -order solution w n + 1 will always violate the L ∞ stability. 2. Second idea : Compute w n + 1 with the 2 nd -order scheme and check whether ∥ w n + 1 ∥ ∞ � ∥ w 0 ∥ ∞ .
Numerical results with ARS(2,2,2) and the MOOD procedure 0 MOOD convex combination 2 nd -order 1 st -order exact solution w x 1 1 0 w x 1 1 otherwise. 1 2 1 36/57 We check the new scheme with the ARS(2,2,2) value of β : √ β = 1 − λ = 1 − β ≃ 1 . 414, θ = 2 β ( 1 − β ) ≃ 0 . 414. 2 , For x ∈ ( 0 , 1 ) , we take a step function as initial condition: { if x ∈ ( 0 . 25 , 0 . 75 ) , w 0 ( x ) = 1 + ε we set c e = c i = 1, t end = ε/ ( 1 + ε ) , and we prescribe periodic BC. 1 . 1 1 . 001 1 . 05 1 . 0005 ε = 10 − 1 ε = 10 − 3 0 . 25 0 . 5 0 . 75 0 . 25 0 . 5 0 . 75 Question : Can we fjnd a “better” value of β ?
37/57 , 2 1 1 0 0 1 0 1 Optimizing the value of β We know that the convex combination scheme is L ∞ -stable if � 1 � λ = λ opt = min θ = θ opt = 2 β ( 1 − β ) . 0 < β < 1, β, 1 − β � How to choose the “best” β ∈ ( 0 , 1 ) ? Remark on the relationship between λ opt , θ opt and β : • λ opt ↗ ⇐ ⇒ less restrictive CFL condition ⇐ ⇒ CPU time ↘ • θ opt ↗ ⇐ ⇒ scheme closer to second-order ⇐ ⇒ precision ↗ θ opt λ opt 0 . 5 0 . 4 0 . 3 1 . 5 0 . 2 0 . 1 β β 0 . 5 0 . 5
38/57 4 MOOD 0 convex 2 nd -order 1 st -order CPU time (ms) 1 2 3 Optimizing the value of β We plot the CPU time with respect to β ∈ ( 0 , 1 ) (when approximating the advection of a smooth initial condition with ε = 0 . 1). β 0 . 2 0 . 4 0 . 6 0 . 8 • We note that, indeed, λ opt ↗ = ⇒ CPU time ↘ . • Something happens around β = 1 2 : let us look at the L ∞ error. �
39/57 0 ARS(2,2,2) is a good compromise between CPU time and error. 2 MOOD convex 2 nd -order 1 st -order 0 2 1 0 1 Optimizing the value of β We plot the L ∞ error with respect to β ∈ ( 0 , 1 ) (when approximating the advection of a smooth initial condition with ε = 0 . 1). L ∞ -error (%) L ∞ -error (%) 1 . 9 1 . 8 1 . 7 β 1 . 5 0 . 2 0 . 4 L ∞ -error (%) 0 . 5 0 . 3 β 0 . 25 β 0 . 5 0 . 2 0 . 4 ⇒ the second-order scheme is not L 2 -stable anymore • β > 1 � 2 = √ • The error of the convex combination is actually minimal for β = 1 − 2 !
Second-order accuracy in space: smooth solution 1 MOOD convex 2 nd -order 1 st -order N 1 2 1 16000 8000 4000 2000 1000 The second-order accuracy in space is achieved in a classical way: N 2 40 • for the second-order scheme : • for the convex combination scheme : 10 1 20 80 160 40/57 ◮ MUSCL reconstruction on the explicit fmux, ◮ BDF2 discretization on the implicit fmux; ◮ limited MUSCL reconstruction on the explicit fmux; ◮ fjrst-order on the implicit fmux (to conserve L ∞ stability). L ∞ -error, ε = 1 L ∞ -error, ε = 10 − 3 10 − 1 10 − 4 10 − 2 10 − 5 10 − 3 Figure above : error lines in L ∞ norm for a smooth solution
Second-order accuracy in space: discontinuous solution 0.5 6000 3000 1500 L 1 N 1 160 24000 80 40 20 The error on a discontinuous solution is measured with: N 1 12000 0.5 160 1 MOOD convex 2 nd -order 1 st -order L 1 N 0.5 1 24000 12000 6000 3000 1500 N 0.5 10 80 j j j w m j w m j w 0 j w 0 j 41/57 10 40 20 ∑ � �� � � ��� ∥ w n ∥ L 1 o = 1 | w n j | + max max j − min − max j − min , ∆ x m � n constructed by adding overshoots and undershoots to the L 1 norm. L 1 -error, ε = 1 o -error, ε = 1 0 . 8 0 . 8 0 . 4 0 . 4 0 . 2 0 . 2 0 . 1 0 . 1 L 1 -error, ε = 10 − 3 o -error, ε = 10 − 3 0 . 8 0 . 4 0 . 4 0 . 2 0 . 2 0 . 1 0 . 1
Application to the Euler equations Recall the fjrst-order IMEX scheme for the Euler system: We apply the same convex combination procedure : 42/57 ρ n + 1 , o1 − ρ n + ∇ · ( ρ u ) n + 1 , o1 = 0 , ( 1 ) ∆ t ( ρ u ) n + 1 , o1 − ( ρ u ) n + ∇ · ( ρ u ⊗ u ) n + 1 ε ∇ p ( ρ n + 1 , o1 ) = 0 . ( 2 ) ∆ t W n + 1 , cc = θ W n + 1 , o2 + ( 1 − θ ) W n + 1 , o1 , with θ = 2 β ( 1 − β ) . � We use the value of θ given by the study of the model problem. � But how can we detect oscillations for the MOOD procedure?
Application to the Euler equations: detection of oscillations : at the same time; 2. detect if both Riemann invariants break the maximum principle MOOD algorithm: to compute the time update of W n , one of them satisfjes a maximum principle in a Riemann problem. 43/57 2 The previous detector ( L ∞ criterion on W ) is irrelevant for the Euler equations, since ρ and u do not satisfy a maximum principle. � We need another detection criterion! � γρ γ − 1 We pick the Riemann invariants Φ ± = u ∓ γ − 1 ε 1. compute the candidate order 2 approximation W n + 1 , o2 ; 3. if so, compute the convex combination W n + 1 , cc .
Application to the Euler equations: numerical results in 1D 1 exact x 1 1 0 x 1 1 We consider a Riemann problem (left rarefaction wave and right shock). 0 x 1 44/57 2 1 0 x 0 1 ρ ρ u 1 . 5 1 . 5 ε = 1 ε = 1 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8 ρ ρ u 1 . 009 1 . 0001 ε = 10 − 4 1 . 006 1 . 00005 ε = 10 − 4 1 . 003 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8
Application to the Euler equations: numerical results in 1D 1 1st order exact x 1 1 0 x 1 1 We consider a Riemann problem (left rarefaction wave and right shock). 0 x 1 44/57 1 2 0 0 1 x ρ ρ u 1 . 5 1 . 5 ε = 1 ε = 1 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8 ρ ρ u 1 . 009 1 . 0001 ε = 10 − 4 1 . 006 1 . 00005 ε = 10 − 4 1 . 003 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8
Application to the Euler equations: numerical results in 1D 1 2nd order 1st order exact x 1 1 0 x 1 1 We consider a Riemann problem (left rarefaction wave and right shock). 0 x 1 44/57 1 0 2 x 1 0 ρ ρ u 1 . 5 1 . 5 ε = 1 ε = 1 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8 ρ ρ u 1 . 009 1 . 0001 ε = 10 − 4 1 . 006 1 . 00005 ε = 10 − 4 1 . 003 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8
Application to the Euler equations: numerical results in 1D 0 x 0 We consider a Riemann problem (left rarefaction wave and right shock). 1 1 x 1 1 1 x exact 1st order 2nd order convex combination 1 44/57 1 0 2 x 1 0 ρ ρ u 1 . 5 1 . 5 ε = 1 ε = 1 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8 ρ ρ u 1 . 009 1 . 0001 ε = 10 − 4 1 . 006 1 . 00005 ε = 10 − 4 1 . 003 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8
Application to the Euler equations: numerical results in 1D 0 x 0 We consider a Riemann problem (left rarefaction wave and right shock). 1 1 x 1 1 1 x exact 1st order 2nd order convex combination MOOD 1 44/57 2 0 1 1 x 0 ρ ρ u 1 . 5 1 . 5 ε = 1 ε = 1 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8 ρ ρ u 1 . 009 1 . 0001 ε = 10 − 4 1 . 006 1 . 00005 ε = 10 − 4 1 . 003 0 . 2 0 . 4 0 . 6 0 . 8 0 . 2 0 . 4 0 . 6 0 . 8
Application to the Euler equations: numerical results in 2D 0 4 6 0 2 4 6 2nd-order AP 0 2 4 6 2 0 4 6 TVD-AP 0 2 4 6 0 2 4 6 AP-MOOD 2 1st-order AP 0 6 2 4 6 0 2 4 6 reference 0 4 2 Reference solution obtained solving the vorticity formulation 0 2 4 6 0 2 4 45/57 − 4 − 2 ∂ t ω + u · ∇ ω = 0, with ω = ∂ x u 1 − ∂ y u 2 ; we take ε = 10 − 5 .
Application to the Euler equations: numerical results in 2D 0 4 6 0 2 4 6 2nd-order AP 0 2 4 6 2 0 4 6 TVD-AP 0 2 4 6 0 2 4 6 AP-MOOD 2 1st-order AP 0 6 2 4 6 0 2 4 6 reference 0 4 2 Reference solution obtained solving the vorticity formulation 0 2 4 6 0 2 4 46/57 − 4 − 2 ∂ t ω + u · ∇ ω = 0, with ω = ∂ x u 1 − ∂ y u 2 ; we take ε = 10 − 5 .
Multi-scale models and their numerical approximation A fjrst-order AP scheme for the isentropic Euler equations [Dimarco, Loubère, Vignal, ’17] High-order schemes in time and their failings [Dimarco, Loubère, M.-D., Vignal, ’18] [M.-D., Thomann, ’20] [M.-D., Thomann, in progress] Conclusion and future work “Second-order” L ∞ -stable schemes “High-order” L ∞ -stable schemes
Problem statement . 0 Goal : Derive a convex combination based on higher order IMEX schemes, b 1 b s c 1 0 0 c 2 0 . . . . . . ... ... . . . c s b 1 b s Question : How to choose all the coeffjcients? We start with a third-order scheme. c s 47/57 . . still considering the model problem with an upwind space discretization. The Butcher tableaux read: Explicit part Implicit part c 1 0 0 0 c 2 0 . . 0 . . . ... ... . ˜ · · · · · · a 1 , 1 ˜ ˜ . . . . . . a 2 , 1 a 2 , 1 a 2 , 2 ˜ ˜ . . . ˜ . . . a s , 1 a s , s − 1 a s , 1 a s , s − 1 a s , s ˜ ˜ . . . . . . . . . . . .
A third-order example 0 b 2 b 3 implicit 0 c 2 c 3 0 0 0 0 0 0 To get relationships between the coeffjcients, we use: • the fjrst-order consistency conditions, Let us consider the following Butcher tableaux. 0 48/57 0 0 b 2 b 3 explicit 0 c 3 c 2 0 0 0 0 ˜ a 2 , 2 a 2 , 1 ˜ ˜ a 3 , 2 a 3 , 3 a 3 , 1 a 3 , 2 • the second-order compatibility conditions, • the third-order compatibility conditions (not detailed here).
A third-order example 2 0 0 0 0 1 implicit 0 0 0 0 0 0 0 0 2 Consequence : We need a convex combination at each sub-step. 0 49/57 0 1 3 . 0 0 2 explicit { } We obtain the following Butcher tableaux, for γ / ∈ 0 , 1 3 γ − 1 3 γ − 1 3 γ − 1 3 γ − 1 6 γ 6 γ 6 γ 6 γ − 6 γ 3 − 3 γ 2 + 1 γ ( 3 γ 2 + 1 ) γ + 1 γ + 1 1 − γ γ 2 ( 3 γ − 1 ) 3 γ − 1 3 γ 2 3 γ 2 3 γ 2 + 1 3 γ 2 + 1 3 γ 2 + 1 3 γ 2 + 1 Applying the L ∞ stability proof, we show that the second step is not L ∞ -stable!
With a convex combination at each sub-step, we get: k t F e W n l k t F i W n l k c k t F e W n k c k t F i W n k 1 t F e W n k 1 t F i W n k F e W n F i W n s 1 s W n W n 1 k The convex combination for an s -step scheme 1 b k s k 1 b k 1 s 1 t 1 s 1 t 1 s 1 1 W n k s a k l k 1 s W n k s k k The classical high-order IMEX scheme reads: 1 l 50/57 a k l l 1 k − 1 ∑ ∑ ∀ k ∈ � 1 , s � , W n , k = W n − ∆ t a k , l ∇ · F e ( W n , l ) − ∆ t a k , l ∇ · F i ( W n , l ) , ˜ l = 1 l = 1 ∑ ∑ W n + 1 = W n − ∆ t ˜ b k ∇ · F e ( W n , k ) − ∆ t b k ∇ · F i ( W n , k ) . k = 1 k = 1
The convex combination for an s -step scheme s s s k The classical high-order IMEX scheme reads: s 50/57 k k − 1 ∑ ∑ ∀ k ∈ � 1 , s � , W n , k = W n − ∆ t a k , l ∇ · F e ( W n , l ) − ∆ t a k , l ∇ · F i ( W n , l ) , ˜ l = 1 l = 1 ∑ ∑ W n + 1 = W n − ∆ t ˜ b k ∇ · F e ( W n , k ) − ∆ t b k ∇ · F i ( W n , k ) . k = 1 k = 1 With a convex combination at each sub-step, we get: k − 1 ∑ ∑ ∀ k ∈ � 1 , s � , W n , k = W n − θ k ∆ t a k , l ∇ · F e ( W n , l ) − θ k ∆ t a k , l ∇ · F i ( W n , l ) ˜ l = 1 l = 1 c k ∆ t ∇ · F e ( W n ) − ( 1 − θ k ) c k ∆ t ∇ · F i ( W n , k ) , − ( 1 − θ k ) ˜ ∑ ∑ W n + 1 = W n − θ s + 1 ∆ t ˜ b k ∇ · F e ( W n , k ) − θ s + 1 ∆ t b k ∇ · F i ( W n , k ) k = 1 k = 1 − ( 1 − θ s + 1 ) ∆ t ∇ · F e ( W n ) − ( 1 − θ s + 1 ) ∆ t ∇ · F i ( W n + 1 ) .
51/57 3 combination, we obtain the following result. Theorem : The third-order scheme with convex combination is and under the CFL condition Applying the proof to the third-order scheme with convex L ∞ stability of the third-order convex combination L ∞ -stable provided that √ θ 4 < ( 3 γ − 1 )( 3 γ 2 + 1 ) θ 3 = 3 γ − 1 γ � 3 , θ 1 = 1 , θ 2 = 1 , 6 γ 2 , , 18 γ 3 18 γ 3 θ 4 − ( 3 γ − 1 )( 3 γ 2 + 1 ) λ � ( 3 γ − 1 )(( 6 γ 2 + 1 ) θ 4 − ( 3 γ 2 + 1 )) . Next step : fjnd optimal values of γ , θ 4 and λ .
16 4 11 4 52/57 0 3 0 0 5 1 0 0 1 0 2 0 4 11 4 0 0 5 1 0 0 5 1 We get a tradeofg between CPU time ( ) and precision ( )! 16 1 1 and 8. 2. Then, the following bounds hold: 0 4 7 16 0 7 7 3. Introduce 0 1 and recast these inequalities as follows: 4 7 16 and L ∞ stability of the third-order convex combination 1. The value of θ 3 is maximal for γ = 2 3, and we get θ 3 = 3
52/57 0 4 0 5 1 0 0 1 0 2 0 3 4 16 0 0 5 1 0 0 5 1 We get a tradeofg between CPU time ( ) and precision ( )! 0 11 16 1 and 16 7 4 0 1 and recast these inequalities as follows: 3. Introduce 1 8. 2. Then, the following bounds hold: and L ∞ stability of the third-order convex combination 1. The value of θ 3 is maximal for γ = 2 3, and we get θ 3 = 3 0 < λ < 7 − 16 θ 4 0 < θ 4 < 7 . 7 − 11 θ 4
52/57 0 0 5 0 0 1 0 2 0 3 0 4 and 4 0 5 0 1 0 0 5 and 16 1 2. Then, the following bounds hold: 8. We get a tradeofg between CPU time ( ) and precision ( )! 1 L ∞ stability of the third-order convex combination 1. The value of θ 3 is maximal for γ = 2 3, and we get θ 3 = 3 0 < λ < 7 − 16 θ 4 0 < θ 4 < 7 . 7 − 11 θ 4 3. Introduce α ∈ ( 0 , 1 ) and recast these inequalities as follows: 1 − α θ 4 = 7 16 α λ = 16 α. 1 − 11
Recommend
More recommend
