Adam Larat

Notes on the Discontinuous Galerkin methods for the numerical simulation of hyperbolic equations

1 General Context

Bibliography

The roots of Discontinuous Galerkin (DG) methods is usually attributed to Reed and Hills in a paper published in 1973 on the numerical approximation of the neutron transport equation [START_REF] Reed | National topical meeting on mathematical models and computational techniques for analysis of nuclear systems[END_REF]. In fact, the adventure really started with a rather thoroughfull series of five papers by Cockburn and Shu in the late 80's [START_REF] Cockburn | The runge-kutta local projection p 1 -discontinuousgalerkin finite element method for scalar conservation laws[END_REF][START_REF] Cockburn | TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework[END_REF][START_REF] Cockburn | Tvb runge-kutta local projection discontinuous galerkin finite element method for conservation laws iii: Onedimensional systems[END_REF][START_REF] Cockburn | Tvb runge-kutta local projection discontinuous galerkin finite element method for conservation laws iv: the multidimensional case[END_REF][START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V -multidimensional systems[END_REF]. Then, the fame of the method, which could be seen as a compromise between Finite Elements (the center of the method being a weak formulation) and Finite Volumes (the basis functions are defined cell-wise, the cells being the elements of the primal mesh) increased and slowly investigated successfully all the domains of Partial Differential Equations numerical integration. In particular, one can cite the ground papers for the common treatment of convection-diffusion equations [START_REF] Bassi | High-order accurate discontinuous finite element solution of the 2d euler equations[END_REF][START_REF] Bassi | A high-order accurate discontinuous finite element method for the numerical solution of the compressible navier-stokes equations[END_REF] or the treatment of pure elliptic equations [START_REF] Arnold | Unified analysis of discontinuous galerkin methods for elliptic problems[END_REF][START_REF] Ortner | Discontinuous galerkin finite element approximation of nonlinear second-order elliptic and hyperbolic systems[END_REF]. For more information on the history of Discontinuous Galerkin method, please refer to section 1.1 of [START_REF] Hestaven | Nodal Discontinuous Galerkin Methods[END_REF].

Today, DG methods are widely used in all kind of manners and have applications in almost all fields of applied mathematics. (TODO: cite applications and structured/unstructured meshes, steady/unsteady, etc...). The methods is now mature enough to deserve entire text books, among which I cite a reference book on Nodal DG Methods by Henthaven and Warburton [START_REF] Hestaven | Nodal Discontinuous Galerkin Methods[END_REF] with the ground basis of DG integration, numerical analysis of its linear behavior and generalization to multiple dimensions.

Lately, since 2010, thanks to a ground work of Zhang and Shu [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF][START_REF] Zhang | On positivity-preserving high order discontinuous galerkin schemes for compressible euler equations on rectangular meshes[END_REF][START_REF] Zhang | Maximum-Principle-Satisfying and Positivity-Preserving High Order Schemes for Conservation Laws[END_REF][START_REF] Zhang | Positivity-preserving high order discontinuous galerkin schemes for compressible euler equations with source terms[END_REF][START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous galerkin schemes for conservation laws on triangular meshes[END_REF], Discontinuous Galerkin methods are eventually able to combine high order accuracy and certain preservation of convex constraints, such as the positivity of a given quantity, for example. These new steps forward are very promising since it brings us very close to the "Ultimate Conservative Scheme", [START_REF] Van Leer | Towards the ultimate conservative difference scheme V. A second order sequel to Godunov's method[END_REF][START_REF] Abgrall | Toward the ultimate conservative scheme: Following the quest[END_REF].

Hyperbolic conservation laws

This section is just a quick introduction to hyperbolic conservation laws in order to settle the notations.

Let d P N ˚be the number of spatial dimensions and Ω Ă R d be the domain of study. BΩ denotes the frontiers of the domain. To greatly simplify the following notes, we will ignore boundary conditions.

Let m P N ˚be the number of conserved variables and W the vector of these variables. Then, W is ruled by a Conservation Law if there exists a Flux Function " F :

S Ă R m ÝÑ pR m q d W Þ ÝÑ FpWq (1)
such that W verifies the PDE:

BW Bt `Ý Ñ ∇ ¨ F pWq " 0. (2)
Here, the set S denotes the possible physical constraints on the conserved variables (positivity of the density, realizability of the transported moments, etc. . . ). This conservation law is said to be hyperbolic, if the Jacobian

J pWq " B F BW (3)
is diagonalizable with real eigenvalues. In this case, we denote α the maximal absolute eigenvalue of J:

α pWq " max t|λ| , DV P R m such that J pWq .V " λVu . (4)
2 Runge-Kutta discontinuous Galerkin methods

Degrees of freedom

Let Ω Ă R d be the domain of study and

N ď k"1 T k " Ω such that @i, j P 1, N , dim pT i X T j q ă d, (5)
a tesselation of the domain.

In each element T k , we define a funcitonal basis `ψk l ˘l"1,...,P k and associate T k with the finite dimensional vector space:

V k h " # ϕ : Ω ÝÑ R, ϕ " ˜Pk ÿ l"1 a l ψ k l ¸χk + , (6)
where χ k denote the characteristic function of T k , taking value 1 in T k and 0 elsewhere.

Therefore, a numerical solution on M h " pT k q k"1,...,N has exactly P "

N ÿ k"1 P k Degrees of Freedom (DoF).

Spatial weak formulation

The continuous hyperbolic conservation law (2) is approximated on the high order functional vectorial space V h defined in [START_REF] Cockburn | Tvb runge-kutta local projection discontinuous galerkin finite element method for conservation laws iv: the multidimensional case[END_REF]. It reads:

Find W h pt, xq " N ÿ k"1 χ k pxq P k ÿ j"1
W j k ptqψ j k pxq a piecewise high order polynomial solution, such that:

@k 1 P 1, N , @i P 1, P k 1 , ż Ω ´Bt W h pt, xq `Ý Ñ ∇ ¨ F pW h pt, xqq ¯ψi k 1 pxqdx " 0. ( 7
)
Weak Formulation

Since the supports T k and T k 1 are disjointed, k 1 " k always and the problem comes to @k, i,

ż T k ÿ j d t W j k ptqψ j k pxqψ i k pxqdx `żT k Ý Ñ ∇ ¨ F pW h pt, xqq ψ i k pxqdx " 0.
By integration by parts, one easily gets the following three fold formulation:

ÿ j ˆżT k ψ j k ψ i k dx ˙dt W j k ptq loooooooooooooooomoooooooooooooooon 1 ´żT k F pW h pt, xqq ¨∇ψ i k pxqdx loooooooooooooooooomoooooooooooooooooon 2 `żBT k ψ i k psq F pW h pt, sqq ¨ n ds loooooooooooooooooomoooooooooooooooooon 3 " 0 (8)

Remark 1

One usually considers (but we don't have to)

F pW h pt, xqq " ÿ j F `Wj k ptq ˘ψj k pxq. ( 9
)
If not, 2 is integrated by means of an adequate volumic quadrature.

Remark 2

Along the boundaries BT k of T k , F pW h q does not have a mathematical sens other than the flux occuring locally between the two states on both sides. Then the flux appearing in 3 is replaced by the numerical flux:

F pW h pt, sqq ¨ n " F ˚pW ext pt, sq, W int pt, sq; nq . ( 10
)
In the case the flux function is also projected on the basis as in [START_REF] Cockburn | Tvb runge-kutta local projection discontinuous galerkin finite element method for conservation laws iii: Onedimensional systems[END_REF], the global formu-lation simply comes down to:

M ij d t W j ´ K t ij ¨ F pW j q `ÿ εPBT k ż ε ψ i k psqF ˚pW ext pt, sq, W int pt, sq; nq ds looooooooooooooooooooooooomooooooooooooooooooooooooon 4 " 0 (11)
Remark 3

• In general, last integral 4 is estimated by mean of an adequate quadrature of dimension d ´1. This implies a numerous amount of numerical flux computations.

• Luckily, since we transport and update all the derivatives of the solution within the mesh, we don't need to be high-order accurate on all these numerical fluxes.

A first order flux on all the degrees of freedom will provide a high order solution anyway. Therefore, the linear Rusanov flux (28) is very well suited for this purpose, since its evaluation is quite costless and it has good monotonicity properties (see section 2.5.1).

Basis functions

DG methods are split into two main families: modal and nodal DG. These prefixes refer to the choice on the basis functions ψ k l . In general, the modal basis function will be chosen so that the corresponding mass matrix

M k ij " ż T k ψ k i ψ k j ( 12
)
is everywhere diagonal, what greatly simplifies the updates. On the other hand, any access to the value of the solution at a given point (for example for the evaluation of a numerical flux somewhere) will require a combination of all the degrees of freedom in the element:

W h px k q " P k ÿ l"1 a l ψ k l px k q. ( 13
)
On the contrary, nodal basis function are attached to certain points of the considered element. In general, the points are chosen in a way it simplifies the volumic and edge terms evaluations. In particular, on the edges of the elements degrees of freedom should be attached to some known quadrature points. Then, the numerical flux integral evaluation has all the necessary data for the numerous numerical fluxes evaluation directly at hand.

For dimensions stricly higher than 1 and polynomial order at least quadratic, there is a conjecture that there is no orthogonal nodal basis1 . This might even be proven, alas I don't know yet about the demonstration.

Modal basis

• Legendre polynomials: this is the orthogonal interlocked polynomial basis for measure one. It is simply obtained by a Gram-Schmidt orthogonalization process starting from any polynomial basis of increasing order.

On r´1, 1s : 1, x, x 2 ´1{3, . . .

• Other orthogonal basis: Jacobi, Tchebychev, trigonometric, etc.

Nodal basis

Let us restrict to 1D for the moment. Givent a set px i q i"1,...,P k of points in T k , the associated Lagrangian basis is the only basis of polynomials of order P k ´1 such that L j px i q " δ ij . Then it is easy and classic to write

L i pxq " ź j‰i px ´xj q ź j‰i px i ´xj q . ( 15
)
In general, the set of nodes is:

• a regular distribution within the cell:

x j " x i´1{2 `j ∆x P k ´1, j " 0, . . . , P k ´1, (16)
• a set of smart quadrature points, like Gauss-Lobatto quadrature points since they contain the bounds of the interval.

In more dimensions, the generalization of the former discussion is more complex. Only the regular distribution can be algorithmically extended. In particular, the generalization of the Gauss-Lobatto Lagrange polynomials to 2 and 3 D is not obvious, see discution in [START_REF] Hestaven | Nodal Discontinuous Galerkin Methods[END_REF], chapter 6.

Time integration

Now, thanks to the spatial semi-integration, the numerical scheme [START_REF] Dumbser | A unified framework for the construction of one-step finite volume and discontinuous galerkin schemes on unstructured meshes[END_REF] comes down to an Ordinary Differential Equation. It can be integrated by any numerical procedure for ODEs. Here we focus on Runge-Kutta methods, since they offer a hierarchy of increasing order, which is suitable for the time integration of high order spatial weak formulations.

Runge-Kutta methods

Let

y 1 ptq " f pt, yptqq

be an arbitrary ordinary differential equation. An s-stages explicit Runge-Kutta method applied to this equation can be written in the form:

$ ' ' ' ' & ' ' ' ' % y n`1 " y n `∆t s ÿ i"1 b i k i , k i " f ˜tn `ci ∆t, y n `∆t ÿ jăi a ij k j ¸. ( 18
)
This writing can be summed up in the so-called Butcher Tableau:

0 0 c 2 a 21 . . . . . . . . . c s a s1 ¨¨¨a s,s´1 b 1 ¨¨¨¨¨¨b s (19)
These methods are consistent with equation ( 17) when

c i " ÿ j a ij . (20)
Then, under the additional constraint that ř b j " 1, the explicit numerical procedure is stable if the time step is smaller than a bound which depends on f Lipschitz constant.

In certain cases, especially when dealing with stiff multiscale problems, this stability constraint may be too harsh. One can overcome it by going implicit. The intermediate updates now become

k i " f ˜tn `ci ∆t, y n `∆t ÿ j a ij k j ¸, (21)
and the Butcher tableau now looks like

c 1 . . . A c s b 1 ¨¨¨b s ( 22
)
where A is a full matrix. The price to pay for the increased range of stability is a interdependance between all the intermediate steps, which is not the case in the explicit version, thanks to the triangular shape of A. The solution at time t `∆t is usually the solution of a big non-linear system. The accuracy study of Runge-Kutta methods is more complex. All I want to say here is that it is known that starting from s ě 5, there is no more RK method of order s.

Strong Stability Preserving integrators

In a fundamental paper published in 1988 [START_REF] Shu | Total-variation-diminishing time discretizations[END_REF], Chi-Wang Shu selects among all the RK methods a family of integrators he first calls "Total Variation Diminishing time discretizations", that will later be renamed as Strong Stability Preserving time discretizations in a review paper [START_REF] Gottlieb | Strong stability-preserving highorder time discretization methods[END_REF].

These methods can be seen as those which can be written as a convex combination of Explicit Euler (EE) time integrations:

y n`1 " s´1 ÿ i"0 a i y piq `∆tb i f `t `i∆t, y piq ˘. (23)
Since the stability constraint of the EE scheme is 1.0, these methods are stable under the fact that the time step ∆t is smaller than the Lipschitz constant of f time

CFL max " max i a i b i . (24)
Next, Shu and Osher [START_REF] Shu | Efficient implementation of essentially nonoscillatory shock-capturing schemes[END_REF][START_REF] Shu | Efficient implementation of essentially nonoscillatory shock-capturing schemes, {II}[END_REF] proved that up to fourth order, there exist an optimal SSP-RK integrator, meaning being both

• of the order of the number of stages,

• with the largest possible stability constraint: CFL max " 1.0.

Example 1 (Heun's Method )

This is the optimal second order SSP-RK integration :

# y n`1 2 " y n `∆tf pt, y n q , y n`1 " 1 2 " y n `´y n`1 2 `∆tf ´t `∆t, y n`1 2 ¯¯ı . (25)
2.5 Convex state preserving DG methods

Convex state preserving numerical flux

Let S Ă R m , convex, be the set of physical states.

Example 2 (Convex Constraints)

• Euler Equations: Density and pressure have to stay positive. If the Equation Of State (EOS) follows the Bethe-Weil conditions, these conditions imply the convexity of the physical states 2 .

• Moments of a Repartition Function: if the transported conserved quantities are the moments of a positive repartition function, then these moments need to stay moments of a positive distribution. In many cases, this implies convex constraints on the set of moments, [START_REF] Dette | The theory of canonical moments with applications in statistics, probability, and analysis[END_REF].

A numerical flux F ˚`W `, W ´; n ˘is said to be Convex State Preserving (CSP), when, for any states W n i´1 , W n i , W n i`1 P S, the Explicit Euler update belongs to S,

W n`1 i " W n i ´ν " F ˚`W n i`1 , W n i ; n ˘´F ˚`W n i , W n i´1 ; n ˘‰ P S, (26)
under a CFL constraint:

ν ď C. (27)

CSP Numerical Flux

Example 3

• Rusanov Flux:

F ˚`W `, W ´; n ˘" FpW `q ` FpW ´q 2 ¨ n ´α `W`´W´˘. (28)
α being defined by (4).

• Godunov Flux: since the Godunov flux is a convex combination of the physical states coming from the exact resolution of the Riemann problem at the interface, the update is physical and (26) comes true.

• HLLx Solvers: the same reasoning applies to HLL solvers, since the updated states are convex combinations of physical states, even though the resolution of the Riemann problem is approximated.

Corrolary 4 (SSP Methods)

This naturally extends to SSP integrators, since they are convex combination of Explicit Euler updates. Only the CFL condition must be multiplied by the additional constraint [START_REF] Zanotti | Space-time adaptive ader discontinuous galerkin finite element schemes with a posteriori sub-cell finite volume limiting[END_REF]. Hence the particular role of optimal SSP-RK integrators.

Application to RK-DG methods

Recall the general update of all the degrees of freedom of an element T k :

M ij d t W k j ´żT k F pWpt, xqq ¨∇ψ k i dx `żBT k ψ i k psqF ˚pW ext pt, sq, W int pt, sq; nq ds " 0. (11
) Once summed up on all the DoFs, we obtain the equation ruling the mean value W k :

|T k | dW k dt `żBT k F ˚¨ n ds " 0. (29)
If the ODE ( 11) is integrated by Explicit Euler, the update simply gives:

W k n`1 " W k n ´∆t |T k | ż BT k F ˚¨ n ds " 0. (30)
Note that by Corrolary 4, this discussion naturally extends to SSP integrators. Now, let's switch to 1D and assume that the polynomial order is P k ´1, so that there exists a number Q of Gauss-Lobatto quadrature points, p2Q ´3q ě pP k ´1q, such that the following equality is exact:

W k " Q ÿ q"1 ω q W k px q q. ( 31
)
Then, discrete update (30) can be rewritten into

W k n`1 " Q ÿ q"1 ω q W n k px q q ´∆t ∆x k " F ˚`W ḱ`1 , W k ˘´F ˚`W ḱ , W k´1 ˘*, (32)
" Q ÿ q"1 ω q " W n k px q q ´∆t ω q ∆x k " F ˚´W n k px q`1 q, W n k px q q ¯´F ˚´W n k px q q, W n k px q´1 q ¯* ,
with the convention that W n k px Q`1 q " W ḱ`1 and W n k px 0 q " W k´1 . In this form, we see that the DG update in the mean writes as a convex combination of abstract Euler Explicit updates at the Gauss-Lobatto quadrature points. So that

• if the numerical flux is CSP (26),
• if the solution at time t n belongs to S at all the Gauss-Lobatto quadrature points, @q P 0, Q `1 , W n k px q q P S,

the updated mean value W k n`1 will be in S, under a CFL constraint

ν ď ω 1 C, ( 34
)
since it is known that the smallest weights of the Gauss-Lobatto quadrature are always at the borders of the interval:

ω 1 " ω Q " min ω q .
Starting from that, Zhang and Shu [START_REF] Zhang | On maximum-principle-satisfying high order schemes for scalar conservation laws[END_REF] provided a limitation procedure which was proven to conserve accuracy. Given that the updated mean value W k n`1 belongs to S and that S is a convex set, we now look at the values of the update solution at the Gauss-Lobatto quadrature points. If any of these values W n`1 k px q q is outside S, there exists a unique value θ q Ps0, 1r, such that θ q W n`1 k px q q `p1 ´θq qW k n`1 is back on the frontier of S. By setting θ " min

q θ q , (35)
and

Ă W n`1 k " θW n`1 k `p1 ´θqW k n`1 , (36)
one obtains a new solution which is as accurate as W n`1 k but limited in a way that (32) will propagate the convex constraint preservation further.

ADER-DG integration

ADER-DG methods always start by a space-time integration of a weak formulation in space of equation [START_REF] Arnold | Unified analysis of discontinuous galerkin methods for elliptic problems[END_REF].

ż t n`1 t n "ż T k ˆBW Bt `Ý Ñ ∇ ¨ F pWq ˙¨ψ k i pXq dx  dt " 0, ( 37
) ô M ij ´Wk,n`1 j ´Wk,n j ¯´ż t n`1 t n ż T k F pWq ¨∇ψ k i dxdt looooooooooooooooomooooooooooooooooon 5 `ż t n`1 t n ż BT k ψ k i F ˚pt, nq ds dt looooooooooooooomooooooooooooooon 6 " 0.
In the two next sections, the goal is to set up a method to obtain the necessary data needed to integrate the two integrals 5 and 6 in last equation. In a first procedure, the time evolution is gotten from an arbitrary high order Taylor expansion in time of the conserved variables. Thanks to the Cauchy-Kowaleski procedure (also known as Lax-Wendroff procedure), the time derivatives of W are functions of the spatial derivatives of the fluxes, which are known at time t n . A Generalized Riemann Problem (GRP) is solved at each interface and the problem can be moved to next time step.

In a second procedure, a predictor step is ran in the form of a local space-time DG scheme. This allows to get a local information on the evolution of the data. The predicted solution is then used to compute the space-time flux and volume terms of equation (37).

Cauchy-Kowaleski ADER procedure

Generalized Riemann Problem

Most of this paragraph takes its roots in chapter 19 and 20 of Toro's book [START_REF] Toro | Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction[END_REF].

In a general context, the numerical flux needed to compute last integral 6 of equation (37) comes from the (possibly approximate) resolution of a Riemann problem:

$ ' ' & ' ' % B t W `B n F pWq " 0, t ą 0, x ¨ n P R, Wpt " 0, x ¨ n ă 0q " W ´, Wpt " 0, x ¨ n ą 0q " W `. (38)
By autosimilarity, this solution depends only on one variable ξ " x ¨ n{t and the numerical flux usually writes

F ˚pW `, W ´; nq " F pW ˚pξ " 0qq ¨ n, (39)
where W ˚pξ " 0q is the values taken by the (possibly approximate) solution of the Riemann Problem along the ordinate axis. At higher order in space, input values W `and W ´are fed with the limit values at the interface of the polynomial in each cell: W n k`1 px ì`1 2 q and W n k px í`1 2 q. Therefore a first order error in time comes from the fact that the Riemann solvers see a constant extrapolation of W `and W ´in the neighboring cells. This stimulated people to look at Generalized Riemann Problems (GRP) where the numerical flux would now be:

F ˚`t; W `px ą 0q, W ´px ă 0q; n ˘" F pW ˚pt, x " 0qq ¨ n, t Ps0, ∆tr. (40)
The accuracy in time would then be achieved by a Taylor expansion of W in time along the ordinate axis

W ˚pt, x " 0q " W ˚p0 `, x " 0q `N ÿ k"1 t k k! B k W Bt k p0 `, x " 0q `O `tN`1 ˘, (41)
and to look at an accurate enough way to solve all the unknowns (meaning the successive time derivatives at t " 0 `). It came out that the problem can be split into N sub-Riemann problems, one being the full space consistent Riemann problem (38), called the space high-order Riemann problem, the pN ´1q others being rather simple.

In fact, since W is the solution of ( 2), its k th -order time derivatives can be written as a function of its spatial derivatives up to order k: x Wpt " 0, x " 0 ´q if x ă 0, B pkq x Wpt " 0, x " 0 `q if x ą 0, (43) where W ˚pt " 0 `, x " 0q is the solution of the space high-order Riemann problem and J is the Jacobian of the flux at this state.

B pkq t Wpt, xq " G `Wpt,

Integration within DG formulation

Once this is done, the Taylor expansion (41) can be filled up to desired order N thanks to the Lax-Wendroff expansions (42), and the numerical flux is gained as the exact integration on t P r0, ∆ts of

F ´W˚p t, x " 0q ¯.
Next, the volume integral 5 needs to be evaluated. This is generally done by a numerical quadrature at a sufficient order of accuracy:

ż t n`1 t n ż T k F pWpt, xqq ¨∇ψ k i dxdt " Nt ÿ ts"1 Nx ÿ xq"1
ω s ω q ∆t∆x F pWpt s , x q qq ¨∇ψ k i pt s , x q q. (44)

But the solution is supposed to be regular within the cell (since it is locally approximated by a polynomial), so that the Taylor and Lax-Wendroff expansions (41) and ( 42) still apply at any spatial quadrature point x q . The same procedure can then be applied at every spatial quadrature point and the whole method is evaluated everywhere and each degree of freedom is updated thanks to (37). For a general algorithm to obtain the successive time derivative as in (42), in particular in the context of the Euler flux, see [START_REF] Dumbser | Building blocks for arbitrary high order discontinuous galerkin schemes[END_REF] and references therein.

Local space-time predictor ADER procedure

As we can see in the previous section, the Taylor expansion in time of the solution may lead to extremely complex moving from known spatial derivatives to time derivatives through Lax-Wendroff recursive expansion. Especially in the context of non-linear fluxes, or even worse, when the flux is an unknown black box. . . Anyway, another much more straightforward method is under development. As far as I know, this method really starts in [START_REF] Dumbser | Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws[END_REF] in the finite volume context and has recent applications in the DG context [START_REF] Zanotti | Space-time adaptive ader discontinuous galerkin finite element schemes with a posteriori sub-cell finite volume limiting[END_REF][START_REF] Dumbser | Building blocks for arbitrary high order discontinuous galerkin schemes[END_REF] (clearly not a throughful bibliography).

The idea here is that, in order to compute the time dependant numerical flux 6 in (37), only the outgoing information from each cell is needed. So, on a local space-time mesh, the solution can be evolved within each cell independantly. A predictor solution Q n k pt, xq is obtained and is next used to compute integrals 5 and 6 . Now, we restrict our study to the space-time slab T n k " T k ˆrt n , t n`1 s and suppose that two functional basis pψ t s q s"1,...,Nt and `ψx q ˘q"1,...,Nx are respectively defined on rt n , t n`1 s and T k , see sections 2.1 and 2.3, such that Q n k expands as

Q n k pt, xq " Nt ÿ s"1 Nx ÿ q"1 Q s,q k ψ t s ptqψ x q pxq. (45)
Then Q n k is supposed to locally verify the space-time problem

$ & % B t Q n k `Ý Ñ ∇ ¨ F pQ n k q " 0, in T k ˆst n , t n`1 s, Q n k pt " 0, xq " W n k pxq, in T k , Q n k pt, xq " Q n k pt, xq, on BT k . (46)
The last boundary condition is a little bit mysterious and corresponds to a certain personal point-of-view. . . Next, a local space-time DG formulation is led with an integration by part only in the time direction:

´żT n k Q n
k pt, xqB t ψ t s ptqψ x q pxq dxdt `żT k `Qn k pt n`1 , xqψ t s pt n`1 qψ x q pxq ´Qn k pt n , xqψ t s pt n qψ x q pxq ˘dx `żT n k Ý Ñ ∇ ¨ F pQ n k q ψ t s ptqψ x q pxq dxdt " 0, @s P 1, N t , q P 1, N x .

(47)

If Q n k and FpQ n k q are supposed to be spanned by the ψ t s ψ x q space-time basis, and the vector of local unknown Q s,q k is supposed to be ordered by time layers: Q " pQ s,q k q s"1,...,Nt, q"1,...,Nx " ´Q1,1 k , . . . , Q

where superscript p.q T stands for transposition, then equation (47) can be rewritten into

´

"`K t ˘T b M x ı Q `"1 t n`1 b M x ‰ Q ´"1 t n b M x ‰ Q `"M t b K x ı F pQq " 0, ô "`p K t q T ´1t n`1 `1t n ˘b M x ‰ Q " " M t b K x ı FpQq, (49)
with the following matrices:

M x i,j " ż T k ψ x i pxqψ x j pxq dx, (50)
We have got a local non-linear system which needs to be solved. Fortunately, according to [START_REF] Dumbser | A unified framework for the construction of one-step finite volume and discontinuous galerkin schemes on unstructured meshes[END_REF], the process would be always contractant and the solution can be reached within a few iterations per cells. Once this is done, the predictor Q n k is simply used in 5 and 6 to update equation (37): @T k P M h , @i P T k ,

  xq, B p1q x Wpt, xq, . . . , B pkqx Wpt, xq ˘. In fact, since only the solution at pt " 0 `, x " 0q is needed, the associated

			(42)
	Now, one can derivate problem (38) k times in space to get a Riemann problem on
	B
	Riemann problem can be simplified into
	$ ' ' ' ' & ' ' ' ' %	B t B pkq ´Bpkq x Wpt, xq ¯`J ´W˚p t " 0 `, x " 0q ¯Bx x Wpt " 0, xq " # B pkq	´Bpkq x Wpt, xq ¯" 0,

pkq

x Wpt, xq.

Personal conversation with V. Perrier

Remark 5

• In a paper published on Arxiv in 2012, [START_REF] Johnson | Outflow positivity limiting for hyperbolic conservation laws. part i: Framework and recipe[END_REF], Johnson and Rossmanith proved that such a limitation procedure on Gauss-Lobatto quadrature points is optimal in 1D.

• We have restrict the discussion to 1D because it is much simpler to explain in this context. However, the procedure generalizes to any dimension [START_REF] Zhang | On positivity-preserving high order discontinuous galerkin schemes for compressible euler equations on rectangular meshes[END_REF][START_REF] Zhang | Maximum-principle-satisfying and positivity-preserving high order discontinuous galerkin schemes for conservation laws on triangular meshes[END_REF][START_REF] Johnson | Outflow positivity limiting for hyperbolic conservation laws. part i: Framework and recipe[END_REF] but in a rather cumbersome manner. The interested reader is encouraged to read these articles and citation within.

Limits of such procedure

Looking for a Strong Stability Preserving strategy or not, generalization of Runge-Kutta integration to higher orders is rather complex, since an increasing number of intermediate stages is needed to obtain an additional order in time. Moreover, all these intermediate stages generally need to be stored for the final t n`1 update. This may be limiting in term of memory usage.

On the top of that, the global time step of the method is chosen a priori and kept during the whole multi-step Runge-Kutta process. Even if a security coefficient is applied globally on the time step, the local physics may evolve rapidly within the time step and the stability constraint on the time progress may be a posteriori violated. This risk increasing with the number of substeps.

These reasons explain why people start to turn to space-time formulation. A fully space-time formulation of DG methods is conceivable but finally looks like a huge implicit formulation: all the degrees of freedom of the mesh are coupled, especially with non-linear equations.

On the other hand, arbitrary high order one-step explicit space-time have appeared under the name of ADER (for Arbitrary high-order schemes using DERivatives). Even though these formulations are intrinsically space-time, the volume and flux terms 2 and 3 , in [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V -multidimensional systems[END_REF], are in fact somehow extrapolated from the available information at time t n and the method can be considered as one-step and explicit. This is the topic of the next section.