High-resolution all-regime schemes Giacomo Dimarco * , Raphal Loubre - - PowerPoint PPT Presentation

High-resolution all-regime schemes

Giacomo Dimarco*, Raphaël Loubère†, Victor Michel-Dansac‡§, Andrea Thomann‖, Marie-Hélène Vignal¶ October 06 & 13, 2020 Séminaire MoCo

*Department of Mathematics and Computer Science, University of Ferrara, Italy †CNRS et Institut de Mathématiques de Bordeaux ‡Inria, centre Nancy - Grand Est, équipe-projet TONUS §Institut de Recherche Mathématique Avancée, Strasbourg ‖Dipartimento di Scienze e Alta Tecnologia, Università degli Studi dell’Insubria, Como, Italy ¶Institut de Mathématiques de Toulouse, Université Paul Sabatier Institut de Recherche Mathématique Avancée

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 “Second-order” L∞-stable schemes

[Dimarco, Loubère, M.-D., Vignal, ’18] [M.-D., Thomann, ’20]

“High-order” L∞-stable schemes

[M.-D., Thomann, in progress]

Conclusion and future work

Multi-scale models

Consider a multi-scale model Mε, depending on a parameter ε. In the space-time domain, ε can:

• be of the same order as the reference scale,
• be small compared to the reference scale,
• take intermediate values.

When ε → 0, we assume that the multi-scale model Mε goes towards a limit model M0: M0 = lim

ε→0 Mε.

Goal: Build a numerical scheme that “handles” (in a sense to be defjned) all scales of ε in Mε: that way, only one scheme is needed whatever the value of ε. Example of a multi-scale model: the isentropic Euler equations

1/57

The isentropic Euler equations: statement

Consider a compressible fmuid, whose fmow is described by the isentropic Euler equations for x ∈ Ω ⊂ Rd and t > 0: { ∂tρ + ∇ · (ρu) = 0, ∂t(ρu) + ∇ · (ρu ⊗ u) + ∇p(ρ) = 0, where ρ > 0 is the fmuid density, u ∈ Rd is its velocity, and p(ρ) = ργ is its pressure. Suitably rescaling the equations yields:

t

u 1

t

u u u 1 p 2 where M2 u p is the squared Mach number.

2/57

The isentropic Euler equations: statement

Consider a compressible fmuid, whose fmow is described by the isentropic Euler equations for x ∈ Ω ⊂ Rd and t > 0: { ∂tρ + ∇ · (ρu) = 0, ∂t(ρu) + ∇ · (ρu ⊗ u) + ∇p(ρ) = 0, where ρ > 0 is the fmuid density, u ∈ Rd is its velocity, and p(ρ) = ργ is its pressure. Suitably rescaling the equations yields: (Mε) { ∂tρ + ∇ · (ρu) = 0 (1ε) ∂t(ρu) + ∇ · (ρu ⊗ u) + 1 ε∇p(ρ) = 0 (2ε) where ε = M2 = ∥u∥/p′(ρ) is the squared Mach number.

2/57

The isentropic Euler equations: incompressible limit

(Mε) { ∂tρ + ∇ · (ρu) = 0 (1ε) ∂t(ρu) + ∇ · (ρu ⊗ u) + 1 ε∇p(ρ) = 0 (2ε) To write the formal limit when ε → 0, we take the following boundary conditions: u · n = 0 on ∂Ω. We now formally take ε → 0 in (Mε):                (2ε) = ⇒ ∇p(ρ) = 0 = ⇒ ∇ρ = 0 = ⇒ ρ(x, t) = ρ(t) ∫

Ω

(1ε) = ⇒ ∫

Ω

ρ′(t) + ρ(t) ∫

Ω

∇ · u = 0 = ⇒ |Ω|ρ′(t) + ρ(t) ∫

∂Ω

u · n = 0 = ⇒ ρ′(t) = 0 = ⇒ ρ(x, t) = cst = ⇒

(1ε) ∇ · u = 0 3/57

The isentropic Euler equations: incompressible limit

(Mε) { ∂tρ + ∇ · (ρu) = 0 (1ε) ∂t(ρu) + ∇ · (ρu ⊗ u) + 1 ε∇p(ρ) = 0 (2ε) The rigorous incompressible limit reads: (M0)        ρ = cst = ρ0, ∇ · u = 0, (10) ρ0 ∂tu + ρ0 ∇ · (u ⊗ u) + ∇π1 = 0, (20) where π1, the fjrst-order correction of the pressure, is given by: π1 = lim

ε→0

1 ε

• p(ρ) − p(ρ0)
• .

We get an explicit equation for π1: ∂t(10) − ∇ · (20) = ⇒ −∆π1 = ρ0 ∇2 :(u ⊗ u).

4/57

Asymptotic-preserving schemes

scheme for (Mε) ∆t, ∆x → 0 (Mε) model scheme for (M0) ∆t, ∆x → 0 (M0) model ε → 0

• unif. in ε

? ε → 0

• ∂tρ + ∇ · (ρu) = 0

∂t(ρu) + ∇ · (ρu ⊗ u) + 1 ε∇p(ρ) = 0 ε = M2; M = Mach number        ρ = ρ0 ∇ · u = 0 ρ0∂tu + ρ0∇ · (u ⊗ u) + ∇π1 = 0        ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0 (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0

5/57

Asymptotic-preserving schemes

scheme for (Mε) ∆t, ∆x → 0 (Mε) model scheme for (M0) ∆t, ∆x → 0 (M0) model ε → 0

• unif. in ε

? ε → 0

• ∂tρ + ∇ · (ρu) = 0

∂t(ρu) + ∇ · (ρu ⊗ u) + 1 ε∇p(ρ) = 0 ε = M2; M = Mach number        ρ = ρ0 ∇ · u = 0 ρ0∂tu + ρ0∇ · (u ⊗ u) + ∇π1 = 0        ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0 (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0

5/57

Asymptotic-preserving schemes

scheme for (Mε) ∆t, ∆x → 0 (Mε) model scheme for (M0) ∆t, ∆x → 0 (M0) model ε → 0

• unif. in ε

? ε → 0

• ∂tρ + ∇ · (ρu) = 0

∂t(ρu) + ∇ · (ρu ⊗ u) + 1 ε∇p(ρ) = 0 ε = M2; M = Mach number        ρ = ρ0 ∇ · u = 0 ρ0∂tu + ρ0∇ · (u ⊗ u) + ∇π1 = 0        ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0 (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0

5/57

Asymptotic-preserving schemes

scheme for (Mε) ∆t, ∆x → 0 (Mε) model scheme for (M0) ∆t, ∆x → 0 (M0) model ε → 0

• unif. in ε

? ε → 0

• ∂tρ + ∇ · (ρu) = 0

∂t(ρu) + ∇ · (ρu ⊗ u) + 1 ε∇p(ρ) = 0 ε = M2; M = Mach number        ρ = ρ0 ∇ · u = 0 ρ0∂tu + ρ0∇ · (u ⊗ u) + ∇π1 = 0        ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0 (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0 AP scheme        asymptotic consistency : the numerical solution of (Mε) goes towards a solution of (M0) as ε, ∆t and ∆x tend to 0 asymptotic stability : ∆t does not depend on ε

5/57

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 “Second-order” L∞-stable schemes

[Dimarco, Loubère, M.-D., Vignal, ’18] [M.-D., Thomann, ’20]

“High-order” L∞-stable schemes

[M.-D., Thomann, in progress]

Conclusion and future work

Barrier to the use of explicit schemes

(Mε) { ∂tρ + ∇ · (ρu) = 0 (1ε) ∂t(ρu) + ∇ · (ρu ⊗ u) + 1 ε∇p(ρ) = 0 (2ε) From (Mε), we write the pressure wave equation: ∂t(1ε) − ∇ · (2ε) = ⇒ ∂ttρ − 1 ε ∆p(ρ) = ∇2 : (ρ u ⊗ u). (3ε) From a numerical point of view, we remark:

• Explicit treatment of (3)ε =

⇒ conditional stability ∆t √ε ∆x

• Implicit treatment of (3)ε =

⇒ uniform stability with respect to ε To ensure the asymptotic stability, the discretization of (3)ε has to be implicit.

6/57

Barrier to the use of explicit schemes

Assume ρn and un are known at time tn: the classical explicit time semi-discretization reads        ρn+1 − ρn ∆t + ∇ · (ρu)n = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn ) = 0. (2) Set ρn = ρ0 and ∇ · un = 0. Does ρn+1 = ρ0 and ∇ · un+1 = 0 when ε → 0? 1

n 1 n

t u n t 0 un 2

0un 1 0un

t 0 u u n t1 p 2 un

1

un

2

u u n We do not recover the incompressible limit! Note that we have not actually used the limit 0.

7/57

Barrier to the use of explicit schemes

Assume ρn and un are known at time tn: the classical explicit time semi-discretization reads        ρn+1 − ρn ∆t + ∇ · (ρu)n = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn ) = 0. (2) Set ρn = ρ0 and ∇ · un = 0. Does ρn+1 = ρ0 and ∇ · un+1 = 0 when ε → 0? (1) = ⇒ ρn+1 = ρn − ∆t∇ · (ρu)n = ρ0 − ∆tρ0∇ · un = ρ0 2

0un 1 0un

t 0 u u n t1 p 2 un

1

un

2

u u n We do not recover the incompressible limit! Note that we have not actually used the limit 0.

7/57

Barrier to the use of explicit schemes

Assume ρn and un are known at time tn: the classical explicit time semi-discretization reads        ρn+1 − ρn ∆t + ∇ · (ρu)n = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn ) = 0. (2) Set ρn = ρ0 and ∇ · un = 0. Does ρn+1 = ρ0 and ∇ · un+1 = 0 when ε → 0? (1) = ⇒ ρn+1 = ρn − ∆t∇ · (ρu)n = ρ0 − ∆tρ0∇ · un = ρ0 (2) = ⇒ ρ0un+1 = ρ0un − ∆tρ0∇ · (u ⊗ u)n − ∆t1 ε∇p(ρ0) 2 un

1

un

2

u u n We do not recover the incompressible limit! Note that we have not actually used the limit 0.

7/57

Barrier to the use of explicit schemes

Assume ρn and un are known at time tn: the classical explicit time semi-discretization reads        ρn+1 − ρn ∆t + ∇ · (ρu)n = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn ) = 0. (2) Set ρn = ρ0 and ∇ · un = 0. Does ρn+1 = ρ0 and ∇ · un+1 = 0 when ε → 0? (1) = ⇒ ρn+1 = ρn − ∆t∇ · (ρu)n = ρ0 − ∆tρ0∇ · un = ρ0 (2) = ⇒ ρ0un+1 = ρ0un − ∆tρ0∇ · (u ⊗ u)n − ∆t1 ε∇p(ρ0) ∇ · (2) = ⇒ ∇ · un+1 = ∇ · un − ∇2 : (u ⊗ u)n ̸= 0 We do not recover the incompressible limit! Note that we have not actually used the limit ε → 0.

7/57

Asymptotic-preserving time semi-discretization

Assume ρn and un are known at time tn: an Implicit-Explicit (IMEX) time semi-discretization reads        ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0. (2) Set ρn = ρ0 and ∇ · un = 0. Does ρn+1 = ρ0 and ∇ · un+1 = 0 when ε → 0? 2 p

n 1 n 1 x n 1

1

n 1 n

t n

1

un

1 n n 1 n

1 t 0 un

1

un

1

We recover the incompressible limit! The time semi-discretization is asymptotically consistent.

8/57

Asymptotic-preserving time semi-discretization

Assume ρn and un are known at time tn: an Implicit-Explicit (IMEX) time semi-discretization reads        ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0. (2) Set ρn = ρ0 and ∇ · un = 0. Does ρn+1 = ρ0 and ∇ · un+1 = 0 when ε → 0? (2) = ⇒ ∇p(ρn+1) = 0 = ⇒ ρn+1(x) = ρn+1 1

n 1 n

t n

1

un

1 n n 1 n

1 t 0 un

1

un

1

We recover the incompressible limit! The time semi-discretization is asymptotically consistent.

8/57

Asymptotic-preserving time semi-discretization

Assume ρn and un are known at time tn: an Implicit-Explicit (IMEX) time semi-discretization reads        ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0. (2) Set ρn = ρ0 and ∇ · un = 0. Does ρn+1 = ρ0 and ∇ · un+1 = 0 when ε → 0? (2) = ⇒ ∇p(ρn+1) = 0 = ⇒ ρn+1(x) = ρn+1 ∫

Ω

(1ε) = ⇒ |Ω|ρn+1 = |Ω|ρn − ∆tρn+1 ∫

∂Ω

un+1 · n = ⇒ ρn+1 = ρn = ρ0 1 t 0 un

1

un

1

We recover the incompressible limit! The time semi-discretization is asymptotically consistent.

8/57

Asymptotic-preserving time semi-discretization

Assume ρn and un are known at time tn: an Implicit-Explicit (IMEX) time semi-discretization reads        ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0. (2) Set ρn = ρ0 and ∇ · un = 0. Does ρn+1 = ρ0 and ∇ · un+1 = 0 when ε → 0? (2) = ⇒ ∇p(ρn+1) = 0 = ⇒ ρn+1(x) = ρn+1 ∫

Ω

(1ε) = ⇒ |Ω|ρn+1 = |Ω|ρn − ∆tρn+1 ∫

∂Ω

un+1 · n = ⇒ ρn+1 = ρn = ρ0 (1) = ⇒ ρ0 = ρ0 + ∆tρ0∇ · un+1 = ⇒ ∇ · un+1 = 0 We recover the incompressible limit! The time semi-discretization is asymptotically consistent.

8/57

Asymptotic-preserving time semi-discretization

Assume ρn and un are known at time tn: an Implicit-Explicit (IMEX) time semi-discretization reads        ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0. (2) Is the pressure wave equation discretized implicitly? 1 n

1

1 n t 2

n 1

2 n

n 1

t2 1 p

n 1 2

u u n Yes: the time semi-discretization is asymptotically stable. We get an uncoupled formulation of this AP time semi-discretization: 2 into 1

n 1 n

t u n t p

n 1

t

2

u u n

9/57

Asymptotic-preserving time semi-discretization

Assume ρn and un are known at time tn: an Implicit-Explicit (IMEX) time semi-discretization reads        ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0. (2) Is the pressure wave equation discretized implicitly? (1)n+1 − (1)n ∆t − ∇ · (2) = ⇒ ρn+1 − 2ρn + ρn−1 ∆t2 − 1 ε∆p(ρn+1) = ∇2 : (ρu ⊗ u)n Yes: the time semi-discretization is asymptotically stable. We get an uncoupled formulation of this AP time semi-discretization: 2 into 1

n 1 n

t u n t p

n 1

t

2

u u n

9/57

Asymptotic-preserving time semi-discretization

Assume ρn and un are known at time tn: an Implicit-Explicit (IMEX) time semi-discretization reads        ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0. (2) Is the pressure wave equation discretized implicitly? (1)n+1 − (1)n ∆t − ∇ · (2) = ⇒ ρn+1 − 2ρn + ρn−1 ∆t2 − 1 ε∆p(ρn+1) = ∇2 : (ρu ⊗ u)n Yes: the time semi-discretization is asymptotically stable. We get an uncoupled formulation of this AP time semi-discretization: ∇ · (2) into (1) = ⇒ ρn+1 − ρn ∆t + ∇ · (ρu)n − ∆t ε ∆p(ρn+1) − ∆t∇2 : (ρu ⊗ u)n = 0.

9/57

IMEX formulation

       ρn+1 − ρn ∆t + ∇ · (ρu)n+1 = 0, (1) (ρu)n+1 − (ρu)n ∆t + ∇ · (ρu ⊗ u)n + 1 ε∇p(ρn+1) = 0. (2) Rewrite this time semi-discretization under the following IMEX (Implicit-Explicit) form: Wn+1 − Wn ∆t + ∇ · Fe(Wn) + ∇ · Fi(Wn+1) = 0, where we have set: W =

• ρ

ρu

• ,

Fe(W) =

• ρu ⊗ u
• ,

Fi(W) =

• ρu

1 εp(ρ)I

• .

10/57

CFL condition in 1D

Wn+1 − Wn ∆t + ∂xFe(Wn) + ∂xFi(Wn+1) = 0

• 1. The characteristic velocities coming from Fe and Fi are as follows:

     Fe : λ1

e = 0, λ2 e = 2u,

Fi : λ±

i = ±

• γργ−1

ε , and therefore both Fe and Fi are hyperbolic fmuxes.

• 2. The time step only depends on the explicit velocities:

∆tAP ∆x 2|u|.

• recall ∆tclass.

∆x√ε |u√ε ±

• γργ−1|

− − − →

ε→0 0

• Next step: space discretization of ∂xFe and ∂xFi.

11/57

Finite volume space discretization

We consider the fjnite volume framework:

• 1. the space domain is discretized with cells,
• 2. the cell values Wn

j are updated with a discretization of the fmuxes.

x × xj− 1

2

× xj+ 1

2

Wn

j−1

Wn

j

Wn

j+1

incoming fmux F(Wn

j−1, Wn j )

• utgoing fmux

F(Wn

j , Wn j+1)

We obtain the following fully discrete scheme:

Wn+1

j

= Wn

j + ∆t

∆x

• Fe(Wn

j−1, Wn j ) − Fe(Wn j , Wn j+1)

• + ∆t

∆x

• Fi(Wn+1

j−1 , Wn+1 j

) − Fi(Wn+1

j

, Wn+1

j+1 )

• ,

where the numerical fmuxes Fe and Fi have to be determined.

12/57

L∞-stable space discretization

Between cells WL and WR, we assume that the fmuxes are made of centred and viscosity parts:      Fe(WL, WR) = (Fe(WL) + Fe(WR)) 2 − De(WL, WR)(WR − WL), Fi(WL, WR) = (Fi(WL) + Fi(WR)) 2 − Di(WL, WR)(WR − WL).

• For De, any solver can be used (Rusanov, HLL, …).
• For Di, a linear stability analysis shows that:

◮ Di = 0 (i.e. a centred fmux) leads to a L2-stable scheme, ◮ L∞ stability is achieved by taking (for the Rusanov fmux):

Di(WL, WR) = max  

• γργ−1

L

ε ,

• γργ−1

R

ε  .

13/57

Importance of the upwind implicit viscosity

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).

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 0.14 0.16 0.18 0.2 0.22 0.24 0.26 Exact AP Di/=0 AP Di=0 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 0.14 0.16 0.18 0.2 0.22 0.24 0.26 Exact AP Di/=0 AP Di=0

×: centered implicit discretization = ⇒ L2 stability and less difgusive : upwind implicit discretization = ⇒ L∞ stability but more difgusive

14/57

Asymptotic preservation but difgusive results

ε = 0.99, 300 cells Class.: 273 loops CPU time 0.07 AP: 510 loops CPU time 1.46

0.01 0.02 0.03 0.04 0.05 0.06 1 2 3 4 5 6 x 10

−4

Time Time steps 1st-order AP

• Class. scheme

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8

Density 1st-order AP

• Class. scheme

x

ε = 10−4, 300 cells Class.: 4036 loops CPU time 0.82 AP: 57 loops CPU time 0.14

0.01 0.02 0.03 0.04 0.05 0.06 10

−5

10

−4

10

−3

Time Time steps 1st-order AP

• Class. scheme

0.2 0.4 0.6 0.8 1 0.9999 1 1

Density 1st-order AP

• Class. scheme

x

15/57

Asymptotic preservation but difgusive results

ε = 10−4 Underlying of the viscosity f

0.01 0.02 0.03 0.04 0.05 0.06 10

−6

10

−5

10

−4

Time Time steps AP scheme

• Class. scheme

0.2 0.4 0.6 0.8 1 0.9999 1 1

Density AP scheme

• Class. scheme

x

• It is necessary to use high-order schemes

for better results in all regimes of ε. Goal: build high-order schemes that

• remain in the IMEX framework, to conserve the AP properties;
• remain L∞-stable.

16/57

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 “Second-order” L∞-stable schemes

[Dimarco, Loubère, M.-D., Vignal, ’18] [M.-D., Thomann, ’20]

“High-order” L∞-stable schemes

[M.-D., Thomann, in progress]

Conclusion and future work

Derivation of high-order IMEX schemes

IMEX division: ∂tW + ∇ · Fe(W) + ∇ · Fi(W) = 0 General principle: Assume that Wn is known at time tn and compute the update Wn+1 by introducing s intermediate time steps. t

[

tn

]

tn+1

• t

[

tn

]

tn+1

| | | | |

. . . tn,k . . . s steps

• Quadratures between tn and tn+1, introducing intermediate values:

W(tn+1) = W(tn) − ∫ tn+1

tn

∇ · Fe(W(t)) dt − ∫ tn+1

tn

∇ · Fi(W(t)) dt Wn+1 = Wn − ∆t

s

∑

k=1

˜ bk ∇ · Fe(Wn,k) − ∆t

s

∑

k=1

bk ∇ · Fi(Wn,k)

17/57

Derivation of high-order IMEX schemes

• Expression of the intermediate values at times tn,k = tn + ck ∆t:

W(tn,k) = W(tn) + ∫ tn,k

tn

∂tW(t) dt

• Quadratures between tn and tn,k, ∀k ∈ 1, s:

W(tn,k) = W(tn) − ∫ tn,k

tn

∇ · Fe(W(t)) dt − ∫ tn,k

tn

∇ · Fi(W(t)) dt Wn,k = Wn − ∆t

k−1

∑

l=1

˜ ak,l ∇ · Fe(Wn,l) − ∆t

k

∑

l=1

ak,l ∇ · Fi(Wn,l) t

[

tn

]

tn,k

| | |

t (explicit part)

[

tn

]

tn,k-1 tn,k

| | |

t (implicit part)

[

tn

]

tn,k

| | |

18/57

Derivation of high-order IMEX schemes

The arbitrarily high-order IMEX time semi-discretization reads: ∀k ∈ 1, s, Wn,k = Wn − ∆t

k−1

∑

l=1

˜ ak,l ∇ · Fe(Wn,l) − ∆t

k

∑

l=1

ak,l ∇ · Fi(Wn,l), Wn+1 = Wn − ∆t

s

∑

k=1

˜ bk ∇ · Fe(Wn,k) − ∆t

s

∑

k=1

bk ∇ · Fi(Wn,k). Butcher tableaux (Runge-Kutta time discretizations): Explicit part Implicit part ˜ c1 · · · ˜ c2 ˜ a2,1 . . . . . . . . . ... ... . . . ˜ cs ˜ as,1 . . . ˜ as,s−1 ˜ b1 . . . . . . ˜ bs c1 a1,1 · · · c2 a2,1 a2,2 . . . . . . . . . ... ... . . . cs as,1 . . . as,s−1 as,s b1 . . . . . . bs

19/57

IMEX schemes: a fjrst-order example

We consider the following Butcher tableaux: ˜ b1 ˜ b2 explicit ˜ c1 ˜ c2 ˜ a2,1 b1 b2 implicit c1 c2 a1,1 a2,1 a2,2 To yield a fjrst-order accurate scheme, the s-stage Butcher tableaux have to satisfy the following consistency conditions: ∀l ∈ 1, s, ˜ ck =

s

∑

l=1

˜ ak,l, ck =

s

∑

l=1

ak,l,

s

∑

k=1

˜ bk = 1,

s

∑

k=1

bk = 1.

20/57

IMEX schemes: a fjrst-order example

With the compatibility conditions, the tableaux become: 1 − ˜ b2 ˜ b2 explicit ˜ c2 ˜ c2 1 − b2 b2 implicit c1 c2 c1 a2,2 − c2 a2,2 The simplest example (Forward-Backward Euler) is as follows: 1 explicit 1 1 1 implicit 1 1

21/57

IMEX schemes: a fjrst-order example

W(1) W(2) Wn+1 1 explicit 1 1 W(1) W(2) 1 implicit 1 1 W(1) W(2)

22/57

IMEX schemes: a fjrst-order example

W(1) W(2) Wn+1 1 explicit 1 1 W(1) W(2) 1 implicit 1 1 W(1) W(2)

W(1) = Wn − ∆t 0 ∇ · Fe(W(1)) − ∆t 0 ∇ · Fe(W(2)) − ∆t 0 ∇ · Fi(W(1)) − ∆t 0 ∇ · Fi(W(2))

22/57

IMEX schemes: a fjrst-order example

W(1) W(2) Wn+1 1 explicit 1 1 W(1) W(2) 1 implicit 1 1 W(1) W(2)

W(1) = Wn − ∆t 0 ∇ · Fe(W(1)) − ∆t 0 ∇ · Fe(W(2)) − ∆t 0 ∇ · Fi(W(1)) − ∆t 0 ∇ · Fi(W(2)) W(2) = Wn − ∆t 1 ∇ · Fe(W(1)) − ∆t 0 ∇ · Fe(W(2)) − ∆t 0 ∇ · Fi(W(1)) − ∆t 1 ∇ · Fi(W(2))

22/57

IMEX schemes: a fjrst-order example

W(1) W(2) Wn+1 1 explicit 1 1 W(1) W(2) 1 implicit 1 1 W(1) W(2)

W(1) = Wn − ∆t 0 ∇ · Fe(W(1)) − ∆t 0 ∇ · Fe(W(2)) − ∆t 0 ∇ · Fi(W(1)) − ∆t 0 ∇ · Fi(W(2)) W(2) = Wn − ∆t 1 ∇ · Fe(W(1)) − ∆t 0 ∇ · Fe(W(2)) − ∆t 0 ∇ · Fi(W(1)) − ∆t 1 ∇ · Fi(W(2)) Wn+1 = Wn − ∆t 1 ∇ · Fe(W(1)) − ∆t 0 ∇ · Fe(W(2)) − ∆t 0 ∇ · Fi(W(1)) − ∆t 1 ∇ · Fi(W(2))

22/57

IMEX schemes: a fjrst-order example

W(1) W(2) Wn+1 1 explicit 1 1 W(1) W(2) 1 implicit 1 1 W(1) W(2)

           W(1) = Wn W(2) = Wn − ∆t ∇ · Fe(W(1)) − ∆t ∇ · Fi(W(2)) Wn+1 = Wn − ∆t ∇ · Fe(W(1)) − ∆t ∇ · Fi(W(2))

Remark: Since the weights are equal to the last row, i.e. ˜ bl = ˜ a2,l and ˜ bl = ˜ a2,l for l ∈ {1, 2}, the fjnal update is the same as the last intermediate time step: Wn+1 = Wn − ∆t ∇ · Fe(Wn) − ∆t ∇ · Fi(Wn+1).

22/57

IMEX schemes: a second-order example

To get a second-order scheme, we require additional compatibility conditions on the Butcher tableaux coeffjcients:

s

∑

k=1

bk ck = 1 2,

s

∑

k=1

˜ bk ck = 1 2,

s

∑

k=1

bk ˜ ck = 1 2,

s

∑

k=1

˜ bk ˜ ck = 1 2. A simple example (Implicit-Explicit Midpoint) reads: 1 explicit

1

• 2

1

• 2 0

1 implicit

1

• 2

1

• 2

         W(1) = Wn W(2) = Wn − 1 2∆t ∇ · Fe(W(1)) − 1 2∆t ∇ · Fi(W(2)) Wn+1 = Wn − ∆t ∇ · Fe(W(2)) − ∆t ∇ · Fi(W(2))

23/57

IMEX schemes: a second-order example

To get a second-order scheme, we require additional compatibility conditions on the Butcher tableaux coeffjcients:

s

∑

k=1

bk ck = 1 2,

s

∑

k=1

˜ bk ck = 1 2,

s

∑

k=1

bk ˜ ck = 1 2,

s

∑

k=1

˜ bk ˜ ck = 1 2. A simple example (Implicit-Explicit Midpoint) reads: 1 explicit

1

• 2

1

• 2 0

1 implicit

1

• 2

1

• 2

     W∗ = Wn − 1 2∆t ∇ · Fe(Wn) − 1 2∆t ∇ · Fi(W∗) Wn+1 = Wn − ∆t ∇ · Fe(W∗) − ∆t ∇ · Fi(W∗)

This scheme is not L2-stable!

23/57

IMEX schemes: an L2-stable second-order example

To achieve L2 stability, we can use the ARS(2,2,2) scheme with β = 1 − √ 2 2 : 1 −

1 2β 1 2β

explicit β 1 β 1 −

1 2β 1 2β 1 2(1−β)

1 −

1 2(1−β)

implicit β 1 β

1 2(1−β)

1 −

1 2(1−β)

                 W∗ = Wn − β∆t ∇ · Fe(Wn) − β∆t ∇ · Fi(W∗) Wn+1 = Wn −

• 1 − 1

2β

• ∆t ∇ · Fe(Wn) − 1

2β∆t ∇ · Fe(W∗) − 1 2(1 − β)∆t ∇ · Fi(W∗) −

• 1 −

1 2(1 − β)

• ∆t ∇ · Fi(Wn+1)

24/57

Application to the isentropic Euler equations

We display the density ρ(x) for a rarefaction-shock Riemann problem with several values of ε.

0.2 0.4 0.6 0.8 1 1 1.05 1.1

ε = 10−1

0.2 0.4 0.6 0.8 1 1 1.005 1.01

ε = 10−2

0.2 0.4 0.6 0.8 1 1 1.0005 1.001

ε = 10−3

0.2 0.4 0.6 0.8 1 1 1.00005 1.0001

ε = 10−4

exact solution fjrst-order in time second-order in time

We observe large oscillations when ε goes to 0!

25/57

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 “Second-order” L∞-stable schemes

[Dimarco, Loubère, M.-D., Vignal, ’18] [M.-D., Thomann, ’20]

“High-order” L∞-stable schemes

[M.-D., Thomann, in progress]

Conclusion and future work

Better understand the oscillations: a model problem

To better understand the oscillations, we consider a model problem: ∂tw + ce∂xw + ci ε ∂xw = 0, ce > 0, ci > 0. It is a linear multi-scale advection equation, which we can write as: ∂tw + ∂xfe(w) + ∂xfi(w), fe(w) = cew, fi(w) = ci ε w. The 0 limit is obtained as follows:

xw tw

w cst and the fjrst-order AP time semi-discretization reads: wn

1

wn t ce

xwn

ci

xwn 1

The space discretization will be an upwind fmux.

26/57

Better understand the oscillations: a model problem

To better understand the oscillations, we consider a model problem: ∂tw + ce∂xw + ci ε ∂xw = 0, ce > 0, ci > 0. It is a linear multi-scale advection equation, which we can write as: ∂tw + ∂xfe(w) + ∂xfi(w), fe(w) = cew, fi(w) = ci ε w. The ε → 0 limit is obtained as follows: ε → 0 = ⇒ ∂xw = 0 = ⇒ ∂tw = 0 = ⇒ w = cst, and the fjrst-order AP time semi-discretization reads: wn+1 − wn ∆t + ce∂xwn + ci ε ∂xwn+1 = 0. The space discretization will be an upwind fmux.

26/57

Better understand the oscillations: L∞ stability failure

The oscillations are measured by the L∞ norm and the total variation: ∥wn∥∞ = max

j

|wn

j |

and TV(wn) = ∑

j

|wn

j+1 − wn j |.

TV =

At the continuous level, the L∞ norm and the total variation decrease over time; at the discrete level, we wish to get { L∞ stability: ∥wn+1∥∞ ∥wn∥∞, total variation diminishing: TV(wn+1) TV(wn). ⇐ ⇒ no oscillations First idea: Derive a L -stable second-order IMEX discretization. Obstruction: A theorem [Gottlieb, Shu, Tadmor, ’01] states that, 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.

27/57

Better understand the oscillations: L∞ stability failure

The oscillations are measured by the L∞ norm and the total variation: ∥wn∥∞ = max

j

|wn

j |

and TV(wn) = ∑

j

|wn

j+1 − wn j |.

TV =

At the continuous level, the L∞ norm and the total variation decrease over time; at the discrete level, we wish to get { L∞ stability: ∥wn+1∥∞ ∥wn∥∞, total variation diminishing: TV(wn+1) TV(wn). ⇐ ⇒ no oscillations First idea: Derive a L∞-stable second-order IMEX discretization. Obstruction: A theorem [Gottlieb, Shu, Tadmor, ’01] states that, 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.

27/57

Generic stability proof for an IMEX step

An idealized IMEX step for the model problem, with weights α and β, is: wn+1

j

− wn

j

∆t + βce wn

j − wn j−1

∆x + α ci ε wn+1

j

− wn+1

j−1

∆x = 0. Rearranging the terms, we get: wn+1

j

= wn

j − βce

∆t ∆x

• wn

j − wn j−1

• − α ci

ε ∆t ∆x

• wn+1

j

− wn+1

j−1

• .

Introducing µε = α ci

ε ∆t ∆x and λ = βce ∆t ∆x, the scheme reads:

wn+1

j

= (1 − λ)wn

j + λwn j−1 − µε(wn+1 j

− wn+1

j−1 ).

Next step: Derive conditions on µε and λ such that the IMEX step is L∞-stable, assuming periodic boundary conditions. Remark: λ 1 ⇐ ⇒ ∆t ∆x β 1 ce and µε 1 ⇐ ⇒ ∆t ∆x α ε ci .

28/57

Generic stability proof for an IMEX step

max

j

|wn

j | = (1 − λ) max j

|wn

j | + λ max j

|wn

j−1| j

1 wn

j j

wn

j 1 if 1

0 i.e. 1 and

j

1 wn

j

wn

j 1 j

1 wn

j

wn

j 1 since x

y x y

j

1 wn

1 j

wn

1 j 1 using the time update of the scheme j

1 wn

1 j

wn

1 j 1

since x y x y

j

1 wn

1 j

wn

1 j 1 j

1 wn

1 j j

wn

1 j 1

1

j

wn

1 j j

wn

1 j 1 if 1

0 and

j

wn

1 j 29/57

Generic stability proof for an IMEX step

max

j

|wn

j | = (1 − λ) max j

|wn

j | + λ max j

|wn

j−1|

= max

j

|(1 − λ)wn

j | + max j

|λwn

j−1| if 1 − λ 0 i.e. λ 1 and λ 0 j

1 wn

j

wn

j 1 j

1 wn

j

wn

j 1 since x

y x y

j

1 wn

1 j

wn

1 j 1 using the time update of the scheme j

1 wn

1 j

wn

1 j 1

since x y x y

j

1 wn

1 j

wn

1 j 1 j

1 wn

1 j j

wn

1 j 1

1

j

wn

1 j j

wn

1 j 1 if 1

0 and

j

wn

1 j 29/57

Generic stability proof for an IMEX step

max

j

|wn

j | = (1 − λ) max j

|wn

j | + λ max j

|wn

j−1|

= max

j

|(1 − λ)wn

j | + max j

|λwn

j−1| if 1 − λ 0 i.e. λ 1 and λ 0

= max

j (|(1 − λ)wn j | + |λwn j−1|) j

1 wn

j

wn

j 1 since x

y x y

j

1 wn

1 j

wn

1 j 1 using the time update of the scheme j

1 wn

1 j

wn

1 j 1

since x y x y

j

1 wn

1 j

wn

1 j 1 j

1 wn

1 j j

wn

1 j 1

1

j

wn

1 j j

wn

1 j 1 if 1

0 and

j

wn

1 j 29/57

Generic stability proof for an IMEX step

max

j

|wn

j | = (1 − λ) max j

|wn

j | + λ max j

|wn

j−1|

= max

j

|(1 − λ)wn

j | + max j

|λwn

j−1| if 1 − λ 0 i.e. λ 1 and λ 0

= max

j (|(1 − λ)wn j | + |λwn j−1|)

max

j

|(1 − λ)wn

j + λwn j−1| since |x| + |y| |x + y| j

1 wn

1 j

wn

1 j 1 using the time update of the scheme j

1 wn

1 j

wn

1 j 1

since x y x y

j

1 wn

1 j

wn

1 j 1 j

1 wn

1 j j

wn

1 j 1

1

j

wn

1 j j

wn

1 j 1 if 1

0 and

j

wn

1 j 29/57

Generic stability proof for an IMEX step

max

j

|wn

j | = (1 − λ) max j

|wn

j | + λ max j

|wn

j−1|

= max

j

|(1 − λ)wn

j | + max j

|λwn

j−1| if 1 − λ 0 i.e. λ 1 and λ 0

= max

j (|(1 − λ)wn j | + |λwn j−1|)

max

j

|(1 − λ)wn

j + λwn j−1| since |x| + |y| |x + y|

= max

j

|(1 + µε)wn+1

j

− µεwn+1

j−1 | using the time update of the scheme j

1 wn

1 j

wn

1 j 1

since x y x y

j

1 wn

1 j

wn

1 j 1 j

1 wn

1 j j

wn

1 j 1

1

j

wn

1 j j

wn

1 j 1 if 1

0 and

j

wn

1 j 29/57

Generic stability proof for an IMEX step

max

j

|wn

j | = (1 − λ) max j

|wn

j | + λ max j

|wn

j−1|

= max

j

|(1 − λ)wn

j | + max j

|λwn

j−1| if 1 − λ 0 i.e. λ 1 and λ 0

= max

j (|(1 − λ)wn j | + |λwn j−1|)

max

j

|(1 − λ)wn

j + λwn j−1| since |x| + |y| |x + y|

= max

j

|(1 + µε)wn+1

j

− µεwn+1

j−1 | using the time update of the scheme

max

j

||(1 + µε)wn+1

j

| − |µεwn+1

j−1 || since |x − y| ||x| − |y|| j

1 wn

1 j

wn

1 j 1 j

1 wn

1 j j

wn

1 j 1

1

j

wn

1 j j

wn

1 j 1 if 1

0 and

j

wn

1 j 29/57

Generic stability proof for an IMEX step

max

j

|wn

j | = (1 − λ) max j

|wn

j | + λ max j

|wn

j−1|

= max

j

|(1 − λ)wn

j | + max j

|λwn

j−1| if 1 − λ 0 i.e. λ 1 and λ 0

= max

j (|(1 − λ)wn j | + |λwn j−1|)

max

j

|(1 − λ)wn

j + λwn j−1| since |x| + |y| |x + y|

= max

j

|(1 + µε)wn+1

j

− µεwn+1

j−1 | using the time update of the scheme

max

j

||(1 + µε)wn+1

j

| − |µεwn+1

j−1 || since |x − y| ||x| − |y||

max

j (|(1 + µε)wn+1 j

| − |µεwn+1

j−1 |) j

1 wn

1 j j

wn

1 j 1

1

j

wn

1 j j

wn

1 j 1 if 1

0 and

j

wn

1 j 29/57

Generic stability proof for an IMEX step

max

j

|wn

j | = (1 − λ) max j

|wn

j | + λ max j

|wn

j−1|

= max

j

|(1 − λ)wn

j | + max j

|λwn

j−1| if 1 − λ 0 i.e. λ 1 and λ 0

= max

j (|(1 − λ)wn j | + |λwn j−1|)

max

j

|(1 − λ)wn

j + λwn j−1| since |x| + |y| |x + y|

= max

j

|(1 + µε)wn+1

j

− µεwn+1

j−1 | using the time update of the scheme

max

j

||(1 + µε)wn+1

j

| − |µεwn+1

j−1 || since |x − y| ||x| − |y||

max

j (|(1 + µε)wn+1 j

| − |µεwn+1

j−1 |)

max

j

|(1 + µε)wn+1

j

| − max

j

|µεwn+1

j−1 |

1

j

wn

1 j j

wn

1 j 1 if 1

0 and

j

wn

1 j 29/57

Generic stability proof for an IMEX step

max

j

|wn

j | = (1 − λ) max j

|wn

j | + λ max j

|wn

j−1|

= max

j

|(1 − λ)wn

j | + max j

|λwn

j−1| if 1 − λ 0 i.e. λ 1 and λ 0

= max

j (|(1 − λ)wn j | + |λwn j−1|)

max

j

|(1 − λ)wn

j + λwn j−1| since |x| + |y| |x + y|

= max

j

|(1 + µε)wn+1

j

− µεwn+1

j−1 | using the time update of the scheme

max

j

||(1 + µε)wn+1

j

| − |µεwn+1

j−1 || since |x − y| ||x| − |y||

max

j (|(1 + µε)wn+1 j

| − |µεwn+1

j−1 |)

max

j

|(1 + µε)wn+1

j

| − max

j

|µεwn+1

j−1 |

= (1 + µε) max

j

|wn+1

j

| − µε max

j

|wn+1

j−1 | if 1 + µε 0 and µε 0 j

wn

1 j 29/57

Generic stability proof for an IMEX step

max

j

|wn

j | = (1 − λ) max j

|wn

j | + λ max j

|wn

j−1|

= max

j

|(1 − λ)wn

j | + max j

|λwn

j−1| if 1 − λ 0 i.e. λ 1 and λ 0

= max

j (|(1 − λ)wn j | + |λwn j−1|)

max

j

|(1 − λ)wn

j + λwn j−1| since |x| + |y| |x + y|

= max

j

|(1 + µε)wn+1

j

− µεwn+1

j−1 | using the time update of the scheme

max

j

||(1 + µε)wn+1

j

| − |µεwn+1

j−1 || since |x − y| ||x| − |y||

max

j (|(1 + µε)wn+1 j

| − |µεwn+1

j−1 |)

max

j

|(1 + µε)wn+1

j

| − max

j

|µεwn+1

j−1 |

= (1 + µε) max

j

|wn+1

j

| − µε max

j

|wn+1

j−1 | if 1 + µε 0 and µε 0

= max

j

|wn+1

j

|

29/57

Instability of the ARS(2,2,2) scheme

Recall the ARS(2,2,2) discretization, with β = 1 −

√ 2 2 , λ = ∆t ∆xce and µε = ∆t ∆x ci ε :

                 w∗

j = wn j − βλ(wn j − wn j−1) − βµε(w∗ j − w∗ j−1),

wn+1

j

= wn

j −

• 1 − 1

2β

• λ(wn

j − wn j−1) − 1

2βλ(w∗

j − w∗ j−1)

− 1 2(1 − β)µε(w∗

j − w∗ j−1) −

• 1 −

1 2(1 − β)

• µε(wn+1

j

− wn+1

j−1 ).

• 1. The fjrst step is L -stable if

0 and 0 1

1 .

• 2. Inject

wj wj

1 from the fjrst step into the second one:

wn

1 j

1 1 2 1 1 wn

j

1 2 1 2 wj 2 wj

1

1 2 2 1 wn

1 j

wn

1 j 1

We have 1, so the coeffjcient of wn

j 1 is

loss of L stability!

30/57

Instability of the ARS(2,2,2) scheme

Recall the ARS(2,2,2) discretization, with β = 1 −

√ 2 2 , λ = ∆t ∆xce and µε = ∆t ∆x ci ε :

                 w∗

j = (1 − βλ)wn j + βλwn j−1 − βµε(w∗ j − w∗ j−1),

wn+1

j

= wn

j −

• 1 − 1

2β

• λ(wn

j − wn j−1) − 1

2βλ(w∗

j − w∗ j−1)

− 1 2(1 − β)µε(w∗

j − w∗ j−1) −

• 1 −

1 2(1 − β)

• µε(wn+1

j

− wn+1

j−1 ).

• 1. The fjrst step is L∞-stable if βµε > 0 and 0 βλ 1 =

⇒ λ 1

β.

• 2. Inject

wj wj

1 from the fjrst step into the second one:

wn

1 j

1 1 2 1 1 wn

j

1 2 1 2 wj 2 wj

1

1 2 2 1 wn

1 j

wn

1 j 1

We have 1, so the coeffjcient of wn

j 1 is

loss of L stability!

30/57

Instability of the ARS(2,2,2) scheme

Recall the ARS(2,2,2) discretization, with β = 1 −

√ 2 2 , λ = ∆t ∆xce and µε = ∆t ∆x ci ε :

                 w∗

j = (1 − βλ)wn j + βλwn j−1 − βµε(w∗ j − w∗ j−1),

wn+1

j

= wn

j −

• 1 − 1

2β

• λ(wn

j − wn j−1) − 1

2βλ(w∗

j − w∗ j−1)

− 1 2(1 − β)µε(w∗

j − w∗ j−1) −

• 1 −

1 2(1 − β)

• µε(wn+1

j

− wn+1

j−1 ).

• 1. The fjrst step is L∞-stable if βµε > 0 and 0 βλ 1 =

⇒ λ 1

β.

• 2. Inject µε(w∗

j − w∗ j−1) from the fjrst step into the second one:

wn+1

j

=

• 1 −

1 2β(1 − β) + λ1 − β β

• wn

j + β − 1

β λwn

j−1

+

• 1

2β(1 − β) − λ 2β

• w∗

j + λ

2βw∗

j−1 + 1 − 2β

2(1 − β)µε(wn+1

j

− wn+1

j−1 ).

We have 1, so the coeffjcient of wn

j 1 is

loss of L stability!

30/57

Instability of the ARS(2,2,2) scheme

Recall the ARS(2,2,2) discretization, with β = 1 −

√ 2 2 , λ = ∆t ∆xce and µε = ∆t ∆x ci ε :

                 w∗

j = (1 − βλ)wn j + βλwn j−1 − βµε(w∗ j − w∗ j−1),

wn+1

j

= wn

j −

• 1 − 1

2β

• λ(wn

j − wn j−1) − 1

2βλ(w∗

j − w∗ j−1)

− 1 2(1 − β)µε(w∗

j − w∗ j−1) −

• 1 −

1 2(1 − β)

• µε(wn+1

j

− wn+1

j−1 ).

• 1. The fjrst step is L∞-stable if βµε > 0 and 0 βλ 1 =

⇒ λ 1

β.

• 2. Inject µε(w∗

j − w∗ j−1) from the fjrst step into the second one:

wn+1

j

=

• 1 −

1 2β(1 − β) + λ1 − β β

• wn

j + β − 1

β λwn

j−1

+

• 1

2β(1 − β) − λ 2β

• w∗

j + λ

2βw∗

j−1 + 1 − 2β

2(1 − β)µε(wn+1

j

− wn+1

j−1 ).

We have β < 1, so the coeffjcient of wn

j−1 is 0 loss of L∞ stability! 30/57

A convex combination

The second step wn+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 wn+1,o1: wn+1 = θwn+1,o2 + (1 − θ)wn+1,o1 = ⇒ {

• rder 1, difgusive and stable, if θ = 0,
• rder 2, accurate and unstable, if θ = 1.

The full scheme reads as follows: wj wn

j

wn

j

wn

j 1

wj wj

1

wn

1 j

wn

j

1 1 2 wn

j

wn

j 1

1 2 wj wj

1

1 2 1 wj wj

1

1 1 2 1 wn

1 j

wn

1 j 1

1 wn

j

wn

j 1

1 wn

1 j

wn

1 j 1 31/57

A convex combination

The second step wn+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 wn+1,o1: wn+1 = θwn+1,o2 + (1 − θ)wn+1,o1 = ⇒ {

• rder 1, difgusive and stable, if θ = 0,
• rder 2, accurate and unstable, if θ = 1.

The full scheme reads as follows:                            w∗

j = wn j − βλ(wn j − wn j−1) − βµε(w∗ j − w∗ j−1),

wn+1

j

= wn

j − θ

• 1 − 1

2β

• λ(wn

j − wn j−1) − θ 1

2βλ(w∗

j − w∗ j−1)

− θ 1 2(1 − β)µε(w∗

j − w∗ j−1) − θ

• 1 −

1 2(1 − β)

• µε(wn+1

j

− wn+1

j−1 )

− (1 − θ)λ(wn

j − wn j−1) − (1 − θ)µε(wn+1 j

− wn+1

j−1 ). 31/57

Determination of θ

Next step: Determine the largest θ that ensures the L∞ stability. Applying the TVD proof to the convex combination scheme yields: Theorem: The convex combination scheme is L∞-stable if: 0 < β < 1, λ 1, θ 2β(1 − β). Remark: We are no longer constrained by the L2 condition β = 1 − √ 2 2 . Remark: λ 1 ⇐ ⇒ ∆t ∆xce 1 ⇐ ⇒ ∆t ∆x ce is a CFL condition Corollary: If

• pt

2 1 , the conditions are relaxed: 1 1 1 1

• pt

2 1

32/57

Determination of θ

Next step: Determine the largest θ that ensures the L∞ stability. Applying the TVD proof to the convex combination scheme yields: Theorem: The convex combination scheme is L∞-stable if: 0 < β < 1, λ 1, θ 2β(1 − β). Remark: We are no longer constrained by the L2 condition β = 1 − √ 2 2 . Remark: λ 1 ⇐ ⇒ ∆t ∆xce 1 ⇐ ⇒ ∆t ∆x ce is a CFL condition Corollary: If θ = θopt = 2β(1 − β), the conditions are relaxed: 0 < β < 1, λ min 1 β, 1 1 − β

• ,

θ = θopt = 2β(1 − β).

32/57

Numerical results with ARS(2,2,2)

We check the new scheme with the ARS(2,2,2) value of β: β = 1 −

√ 2 2 ,

λ =

1 1−β ≃ 1.414,

θ = 2β(1 − β) ≃ 0.414. For x ∈ (0, 1), we take a step function as initial condition: w0(x) = { 1 + ε if x ∈ (0.25, 0.75), 1

• therwise.

we set ce = ci = 1, tend = ε/(1 + ε), and we prescribe periodic BC.

0 25 0 5 0 75 1 1 1 05 1 1

x w

10

1

0 25 0 5 0 75 1 1 1 0005 1 001

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?

33/57

Numerical results with ARS(2,2,2)

We check the new scheme with the ARS(2,2,2) value of β: β = 1 −

√ 2 2 ,

λ =

1 1−β ≃ 1.414,

θ = 2β(1 − β) ≃ 0.414. For x ∈ (0, 1), we take a step function as initial condition: w0(x) = { 1 + ε if x ∈ (0.25, 0.75), 1

• therwise.

we set ce = ci = 1, tend = ε/(1 + ε), and we prescribe periodic BC.

0.25 0.5 0.75 1 1 1.05 1.1

x w

ε = 10−1

0.25 0.5 0.75 1 1 1.0005 1.001

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?

33/57

MOOD procedure

Remark : The convex combination removes the oscillations but the approximation remains difgusive. Increase the resolution using the MOOD framework (Multidimensional Optimal Order Detection).

wn compute a candidate approximation wc using the 2nd-order scheme detection of

• scillations

wn+1 = wc compute wn+1 with the convex combination wn+1 no yes

Figure above : MOOD method for one time step (compute wn+1 from wn)

34/57

Detection of oscillations for the model problem

How to detect oscillations?

• 1. First idea:

Compute wn+1 with the 2nd-order scheme and check whether ∥wn+1∥∞ ∥wn∥∞. Problem: If the wn is not smooth enough, the 2nd-order solution wn+1 will always violate the L∞ stability. Consequence: The 2nd-order scheme is never used, unless the solu- tion is smooth.

• 2. Second idea: Compute wn

1 with the 2nd-order scheme and check

whether wn

1

w0 . Why? We know that the numerical solution should not leave the bounds of the initial condition. Consequence: The 2nd-order solution is used as soon the initial bounds are not violated.

35/57

Detection of oscillations for the model problem

How to detect oscillations?

• 1. First idea:

Compute wn+1 with the 2nd-order scheme and check whether ∥wn+1∥∞ ∥wn∥∞. Problem: If the wn is not smooth enough, the 2nd-order solution wn+1 will always violate the L∞ stability. Consequence: The 2nd-order scheme is never used, unless the solu- tion is smooth.

• 2. Second idea: Compute wn+1 with the 2nd-order scheme and check

whether ∥wn+1∥∞ ∥w0∥∞. Why? We know that the numerical solution should not leave the bounds of the initial condition. Consequence: The 2nd-order solution is used as soon the initial bounds are not violated.

35/57

Numerical results with ARS(2,2,2) and the MOOD procedure

We check the new scheme with the ARS(2,2,2) value of β: β = 1 −

√ 2 2 ,

λ =

1 1−β ≃ 1.414,

θ = 2β(1 − β) ≃ 0.414. For x ∈ (0, 1), we take a step function as initial condition: w0(x) = { 1 if x ∈ (0.25, 0.75), 1 + ε

• therwise.

we set ce = ci = 1, tend = ε/(1 + ε), and we prescribe periodic BC.

0.25 0.5 0.75 1 1 1.05 1.1

x w

ε = 10−1

0.25 0.5 0.75 1 1 1.0005 1.001

x w

ε = 10−3

exact solution 1st-order 2nd-order convex combination MOOD

Question: Can we fjnd a “better” value of β ?

36/57

Optimizing the value of β

We know that the convex combination scheme is L∞-stable if 0 < β < 1, λ = λopt = min 1 β, 1 1 − β

• ,

θ = θopt = 2β(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 ↗

0.5 1 0.1 0.2 0.3 0.4 0.5

β θopt

0.5 1 1 1.5 2

β λopt

37/57

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 1 2 3 4

β CPU time (ms)

1st-order 2nd-order convex MOOD

• We note that, indeed, λopt ↗ =

⇒ CPU time ↘.

• Something happens around β = 1
• 2: let us look at the L∞ error.

38/57

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).

0.5 1 0.5 1 1.5 2

β L∞-error (%)

0.2 0.4 1.7 1.8 1.9 β L∞-error (%) 0.2 0.4 0.25 0.3 β L∞-error (%)

1st-order 2nd-order convex MOOD

• β > 1
• 2 =

⇒ the second-order scheme is not L2-stable anymore

• The error of the convex combination is actually minimal for β = 1 −

√ 2 2 !

ARS(2,2,2) is a good compromise between CPU time and error.

39/57

Second-order accuracy in space: smooth solution

The second-order accuracy in space is achieved in a classical way:

• for the second-order scheme:

◮ MUSCL reconstruction on the explicit fmux, ◮ BDF2 discretization on the implicit fmux;

• for the convex combination scheme:

◮ limited MUSCL reconstruction on the explicit fmux; ◮ fjrst-order on the implicit fmux (to conserve L∞ stability).

10 20 40 80 160 10−3 10−2 10−1 1 2 1

N

L∞-error, ε = 1

1000 2000 4000 8000 16000 10−5 10−4 1 2 1

N

L∞-error, ε = 10−3

1st-order 2nd-order convex MOOD

Figure above : error lines in L∞ norm for a smooth solution

40/57

Second-order accuracy in space: discontinuous solution

The error on a discontinuous solution is measured with: ∥wn∥L1

• = 1

∆x ∑

j

• |wn

j | + max mn

• max

j

wm

j − min j

wm

j

• −
• max

j

w0

j − min j

w0

j

• ,

constructed by adding overshoots and undershoots to the L1 norm.

10 20 40 80 160 0.1 0.2 0.4 0.8 0.5 1

N

L1-error, ε = 1

10 20 40 80 160 0.1 0.2 0.4 0.8 0.5 1

N

L1

• -error, ε = 1

1500 3000 6000 12000 24000 0.1 0.2 0.4 0.5 1

N

L1-error, ε = 10−3

1500 3000 6000 12000 24000 0.1 0.2 0.4 0.8 0.5 1

N

L1

• -error, ε = 10−3

1st-order 2nd-order convex MOOD

41/57

Application to the Euler equations

Recall the fjrst-order IMEX scheme for the Euler system:          ρn+1,o1 − ρn ∆t + ∇ · (ρ u)n+1,o1 = 0, (1) (ρ u)n+1,o1 − (ρ u)n ∆t + ∇ · (ρ u ⊗ u)n + 1 ε ∇p(ρn+1,o1) = 0. (2) We apply the same convex combination procedure: Wn+1,cc = θ Wn+1,o2 + (1 − θ) Wn+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?

42/57

Application to the Euler equations: detection of oscillations

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! We pick the Riemann invariants Φ± = u ∓ 2 γ − 1

• γργ−1

ε :

• ne of them satisfjes a maximum principle in a Riemann problem.

MOOD algorithm: to compute the time update of Wn,

• 1. compute the candidate order 2 approximation Wn+1,o2;
• 2. detect if both Riemann invariants break the maximum principle

at the same time;

• 3. if so, compute the convex combination Wn+1,cc.

43/57

Application to the Euler equations: numerical results in 1D

We consider a Riemann problem (left rarefaction wave and right shock).

0.2 0.4 0.6 0.8 1 1 1.5 2

x ρ

ε = 1

0.2 0.4 0.6 0.8 1 1 1.5

x ρu

ε = 1

0.2 0.4 0.6 0.8 1 1 1.00005 1.0001

x ρ

ε = 10−4

0.2 0.4 0.6 0.8 1 1 1.003 1.006 1.009

x ρu

ε = 10−4

exact

44/57

Application to the Euler equations: numerical results in 1D

We consider a Riemann problem (left rarefaction wave and right shock).

0.2 0.4 0.6 0.8 1 1 1.5 2

x ρ

ε = 1

0.2 0.4 0.6 0.8 1 1 1.5

x ρu

ε = 1

0.2 0.4 0.6 0.8 1 1 1.00005 1.0001

x ρ

ε = 10−4

0.2 0.4 0.6 0.8 1 1 1.003 1.006 1.009

x ρu

ε = 10−4

exact 1st order

44/57

Application to the Euler equations: numerical results in 1D

We consider a Riemann problem (left rarefaction wave and right shock).

0.2 0.4 0.6 0.8 1 1 1.5 2

x ρ

ε = 1

0.2 0.4 0.6 0.8 1 1 1.5

x ρu

ε = 1

0.2 0.4 0.6 0.8 1 1 1.00005 1.0001

x ρ

ε = 10−4

0.2 0.4 0.6 0.8 1 1 1.003 1.006 1.009

x ρu

ε = 10−4

exact 1st order 2nd order

44/57

Application to the Euler equations: numerical results in 1D

We consider a Riemann problem (left rarefaction wave and right shock).

0.2 0.4 0.6 0.8 1 1 1.5 2

x ρ

ε = 1

0.2 0.4 0.6 0.8 1 1 1.5

x ρu

ε = 1

0.2 0.4 0.6 0.8 1 1 1.00005 1.0001

x ρ

ε = 10−4

0.2 0.4 0.6 0.8 1 1 1.003 1.006 1.009

x ρu

ε = 10−4

exact 1st order 2nd order convex combination

44/57

Application to the Euler equations: numerical results in 1D

We consider a Riemann problem (left rarefaction wave and right shock).

0.2 0.4 0.6 0.8 1 1 1.5 2

x ρ

ε = 1

0.2 0.4 0.6 0.8 1 1 1.5

x ρu

ε = 1

0.2 0.4 0.6 0.8 1 1 1.00005 1.0001

x ρ

ε = 10−4

0.2 0.4 0.6 0.8 1 1 1.003 1.006 1.009

x ρu

ε = 10−4

exact 1st order 2nd order convex combination MOOD

44/57

Application to the Euler equations: numerical results in 2D

2 4 6 2 4 6

reference

−4 −2 4 2

Reference solution

• btained solving

the vorticity formulation ∂tω + u · ∇ω = 0, with ω = ∂xu1 − ∂yu2; we take ε = 10−5.

2 4 6 2 4 6

1st-order AP

2 4 6 2 4 6

2nd-order AP

2 4 6 2 4 6

TVD-AP

2 4 6 2 4 6

AP-MOOD

45/57

Application to the Euler equations: numerical results in 2D

2 4 6 2 4 6

reference

−4 −2 4 2

Reference solution

• btained solving

the vorticity formulation ∂tω + u · ∇ω = 0, with ω = ∂xu1 − ∂yu2; we take ε = 10−5.

2 4 6 2 4 6

1st-order AP

2 4 6 2 4 6

2nd-order AP

2 4 6 2 4 6

TVD-AP

2 4 6 2 4 6

AP-MOOD

46/57

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 “Second-order” L∞-stable schemes

[Dimarco, Loubère, M.-D., Vignal, ’18] [M.-D., Thomann, ’20]

“High-order” L∞-stable schemes

[M.-D., Thomann, in progress]

Conclusion and future work

Problem statement

Goal: Derive a convex combination based on higher order IMEX schemes, still considering the model problem with an upwind space discretization. The Butcher tableaux read: Explicit part Implicit part ˜ c1 · · · ˜ c2 ˜ a2,1 . . . . . . . . . ... ... . . . ˜ cs ˜ as,1 . . . ˜ as,s−1 ˜ b1 . . . . . . ˜ bs c1 a1,1 · · · c2 a2,1 a2,2 . . . . . . . . . ... ... . . . cs as,1 . . . as,s−1 as,s b1 . . . . . . bs Question: How to choose all the coeffjcients? We start with a third-order scheme.

47/57

A third-order example

Let us consider the following Butcher tableaux. b2 b3 explicit c2 c3 ˜ a2,1 ˜ a3,1 ˜ a3,2 b2 b3 implicit c2 c3 a2,2 a3,2 a3,3 To get relationships between the coeffjcients, we use:

• the fjrst-order consistency conditions,
• the second-order compatibility conditions,
• the third-order compatibility conditions (not detailed here).

48/57

A third-order example

We obtain the following Butcher tableaux, for γ / ∈ { 0, 1 3 } .

3γ2 3γ2 + 1 1 3γ2 + 1

explicit

3γ − 1 6γ γ + 1 2 3γ − 1 6γ −6γ3 − 3γ2 + 1 2(3γ − 1) γ(3γ2 + 1) 3γ − 1 3γ2 3γ2 + 1 1 3γ2 + 1

implicit

3γ − 1 6γ γ + 1 2 3γ − 1 6γ γ 1 − γ 2

Applying the L∞ stability proof, we show that the second step is not L∞-stable! Consequence: We need a convex combination at each sub-step.

49/57

The convex combination for an s-step scheme

The classical high-order IMEX scheme reads: ∀k ∈ 1, s, Wn,k = Wn − ∆t

k−1

∑

l=1

˜ ak,l ∇ · Fe(Wn,l) − ∆t

k

∑

l=1

ak,l ∇ · Fi(Wn,l), Wn+1 = Wn − ∆t

s

∑

k=1

˜ bk ∇ · Fe(Wn,k) − ∆t

s

∑

k=1

bk ∇ · Fi(Wn,k). With a convex combination at each sub-step, we get: k 1 s Wn k Wn

k t k 1 l 1

ak l Fe Wn l

k t k l 1

ak l Fi Wn l 1

k ck t

Fe Wn 1

k ck t

Fi Wn k Wn

1

Wn

s 1 t s k 1

bk Fe Wn k

s 1 t s k 1

bk Fi Wn k 1

s 1

t Fe Wn 1

s 1

t Fi Wn

1 50/57

The convex combination for an s-step scheme

The classical high-order IMEX scheme reads: ∀k ∈ 1, s, Wn,k = Wn − ∆t

k−1

∑

l=1

˜ ak,l ∇ · Fe(Wn,l) − ∆t

k

∑

l=1

ak,l ∇ · Fi(Wn,l), Wn+1 = Wn − ∆t

s

∑

k=1

˜ bk ∇ · Fe(Wn,k) − ∆t

s

∑

k=1

bk ∇ · Fi(Wn,k). With a convex combination at each sub-step, we get: ∀k ∈ 1, s, Wn,k = Wn − θk∆t

k−1

∑

l=1

˜ ak,l ∇ · Fe(Wn,l) − θk∆t

k

∑

l=1

ak,l ∇ · Fi(Wn,l) − (1 − θk)˜ ck∆t∇ · Fe(Wn) − (1 − θk)ck∆t∇ · Fi(Wn,k), Wn+1 = Wn − θs+1∆t

s

∑

k=1

˜ bk ∇ · Fe(Wn,k) − θs+1∆t

s

∑

k=1

bk ∇ · Fi(Wn,k) − (1 − θs+1)∆t∇ · Fe(Wn) − (1 − θs+1)∆t∇ · Fi(Wn+1).

50/57

L∞ stability of the third-order convex combination

Applying the proof to the third-order scheme with convex combination, we obtain the following result. Theorem: The third-order scheme with convex combination is L∞-stable provided that γ √ 3 3 , θ1 = 1, θ2 = 1, θ3 = 3γ − 1 6γ2 , θ4 < (3γ − 1)(3γ2 + 1) 18γ3 , and under the CFL condition λ 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 λ.

51/57

L∞ stability of the third-order convex combination

• 1. The value of θ3 is maximal for γ = 2

3, and we get θ3 = 3 8.

• 2. Then, the following bounds hold:

4

7 16 and 7 16 4 7 11 4

• 3. Introduce

0 1 and recast these inequalities as follows:

4

7 16 and 1 1

11 16

0 5 1 0 1 0 2 0 3 0 4

4

0 5 1 0 5 1

We get a tradeofg between CPU time ( ) and precision ( )!

52/57

L∞ stability of the third-order convex combination

• 1. The value of θ3 is maximal for γ = 2

3, and we get θ3 = 3 8.

• 2. Then, the following bounds hold:

0 < θ4 < 7 16 and 0 < λ < 7 − 16θ4 7 − 11θ4 .

• 3. Introduce

0 1 and recast these inequalities as follows:

4

7 16 and 1 1

11 16

0 5 1 0 1 0 2 0 3 0 4

4

0 5 1 0 5 1

We get a tradeofg between CPU time ( ) and precision ( )!

52/57

L∞ stability of the third-order convex combination

• 1. The value of θ3 is maximal for γ = 2

3, and we get θ3 = 3 8.

• 2. Then, the following bounds hold:

0 < θ4 < 7 16 and 0 < λ < 7 − 16θ4 7 − 11θ4 .

• 3. Introduce α ∈ (0, 1) and recast these inequalities as follows:

θ4 = 7 16α and λ = 1 − α 1 − 11

16α.

0 5 1 0 1 0 2 0 3 0 4

4

0 5 1 0 5 1

We get a tradeofg between CPU time ( ) and precision ( )!

52/57

L∞ stability of the third-order convex combination

• 1. The value of θ3 is maximal for γ = 2

3, and we get θ3 = 3 8.

• 2. Then, the following bounds hold:

0 < θ4 < 7 16 and 0 < λ < 7 − 16θ4 7 − 11θ4 .

• 3. Introduce α ∈ (0, 1) and recast these inequalities as follows:

θ4 = 7 16α and λ = 1 − α 1 − 11

16α.

0.5 1 0.1 0.2 0.3 0.4

α θ4

0.5 1 0.5 1

α λ = ⇒ We get a tradeofg between CPU time (λ) and precision (θ)!

52/57

L∞ stability of an s-step convex combination

Let us now consider a generic high-order IMEX discretization. After (very) lengthy computations, we obtain the following result: Theorem: For k ∈ 1, s and l ∈ 1, k − 1, we defjne Ak = θkakk + (1 − θk)ck, ˜ Ak = θkak1 + (1 − θk)˜ ck, Bkl = θkakl Al , ˜ Bkl = θk˜ akl. In addition, we recursively defjne the following expressions: Ck = ˜ Ak −

k−1

∑

l=2

BklCl, Ckl = ˜ Bkl −

k−1

∑

r=l+1

BkrCrl, Dk = 1 − λ ˜ Ak −

k−1

∑

l=2

BklDl, Dkl = Bkl − λ˜ Bkl −

k−1

∑

r=l+1

BkrDrl. Then, under the following restrictions for k ∈ 1, s and l ∈ 1, k − 1, Ak > 0, Ck 0, Dk 0, Ckl 0, Dkl 0, the scheme is L∞-stable.

53/57

Optimizing the Butcher tableaux coeffjcients

We wish to get:

• a large λ to minimize CPU time,
• large θk to maximize accuracy,

under the following constraints:

• inequality constraints on λ, θk and the Butcher coeffjcients for

the L∞ stability,

• equality constraints on the Butcher coeffjcients to ensure high-
• rder accuracy.

For an s-step scheme, we have:

• s

1 2 s 1 unknowns (the objective function depends on s 1 of them);

• 3s

s2 inequality constraints;

• 4 equality constraints for second-order accuracy, 18 for third-
• rder, 78 for fourth-order, …

54/57

Optimizing the Butcher tableaux coeffjcients

We wish to get:

• a large λ to minimize CPU time,
• large θk to maximize accuracy,

under the following constraints:

• inequality constraints on λ, θk and the Butcher coeffjcients for

the L∞ stability,

• equality constraints on the Butcher coeffjcients to ensure high-
• rder accuracy.

= ⇒ For an s-step scheme, we have:

• (s + 1)2 + s + 1 unknowns (the objective function depends on

s + 1 of them);

• 3s + s2 inequality constraints;
• 4 equality constraints for second-order accuracy, 18 for third-
• rder, 78 for fourth-order, …

54/57

Numerical results (pth-order, s-step schemes)

• “second-order” schemes (p = 2):

◮ two steps (s = 2): λ ≃ 1.41, 1

s

∑ θ ≃ 0.804

◮ three steps (s = 3): λ ≃ 3.00, 1

s

∑ θ ≃ 0.905

• “third-order” schemes (p = 3):

◮ three steps (s = 3): λ ≃ 0.87, 1

s

∑ θ ≃ 0.630

◮ four steps (s = 4): λ ≃ 0.55, 1

s

∑ θ ≃ 0.801

2.5 5 7.5 10 1 1.05 1.1

x w

ε = 10−1

250 500 750 1,000 1 1.0005 1.001

x w

ε = 10−3

s = 2 & p = 2 s = 3 & p = 2 exact solution s = 3 & p = 3 s = 4 & p = 3

55/57

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 “Second-order” L∞-stable schemes

[Dimarco, Loubère, M.-D., Vignal, ’18] [M.-D., Thomann, ’20]

“High-order” L∞-stable schemes

[M.-D., Thomann, in progress]

Conclusion and future work

Conclusion

• We have presented a technique to obtain high-resolution, L∞

∞ ∞-

stable IMEX Runge-Kutta discretizations:

◮ based on a convex combination between fjrst-order IMEX and

high-order IMEX,

◮ using the MOOD framework to enhance the resolution.

• We have applied this technique to the isentropic Euler equa-

tions, to get a high-resolution, L∞-stable and asymptotic- preserving scheme.

56/57

Future work

A few perspectives for this work include:

• the extension to the full Euler system (but obtaining an IMEX

form with hyperbolic fmuxes is challenging);

• the locality in space of θ (potentially thwarted by the non-

locality in space of the implicit discretization)

• enhancing the optimization of the Butcher coeffjcients, θ and λ;
• improving the oscillation detection in the Euler case to get

more robust MOOD approximations;

• in the long term, introducing a domain decomposition with re-

spect to ε, to use whichever scheme is more adapted at every scale.

57/57

Thank you for your attention !