El-Haj Laamri
email: elhaj.laamri@iecn.u-nancy.fr

Global Existence of classical solutions for a class of reaction-diffusion systems

Global existence of classical solutions for a class of reaction-diffusion systems

El Haj Laamri

Introduction

This paper is motivated by the general question of global existence in time of solutions to the following reaction-diffusion system (S)

                       u t -d 1 ∆u = w γ -u α v β (0, +∞) × Ω, (E 1 ) v t -d 2 ∆v = w γ -u α v β (0, +∞) × Ω, (E 2 ) w t -d 3 ∆w = -w γ + u α v β (0, +∞) × Ω, (E 3 ) ∂u ∂n (t, x) = ∂v ∂n (t, x) = ∂w ∂n (t, x) = 0 (0, +∞) × ∂Ω, u(0, x) = u 0 (x) ≥ 0 x ∈ Ω, v(0, x) = v 0 (x) ≥ 0 x ∈ Ω, w(0, x) = w 0 (x) ≥ 0 x ∈ Ω,
where Ω is a bounded regular open subset of R N , (d 1 , d 2 , d 3 , α, β, γ) ∈ (0, +∞) 3 × [1, +∞) 3 . Note that the system (S) satisfies two main properties, namely : (P ) the nonnegativity of solutions of (S) is preserved for all time ; (M ) the total mass of the components u, v, w is a priori bounded on all finite intervals (0, t).

If α, β and γ are positive integers, system (S) is intended to describe for example the evolution of a reversible chemical reaction of type αU + βV ⇋ γW where u, v, w stand for the density of U , V and W respectively. This chemical reaction is typical of general reversible reactions and contains the major difficulties encountered in a large class of similar problems as regards global existence of solutions.

Let us make precise what we mean by solution.

By classical solution to (S) on Q T = (0, T ) × Ω, we mean that, at least (i) (u, v, w) ∈ C([0, T ); L 1 (Ω) 3 ) ∩ L ∞ ([0, τ ] × Ω) 3 , ∀τ ∈ (0, T ) ; (ii) ∀k, ℓ = 1 . . . N , ∀p ∈ (1, +∞)

∂ t u, ∂ t v, ∂ t w, ∂ x k u, ∂ x k v, ∂ x k w, ∂ x k x ℓ u, ∂ x k x ℓ v, ∂ x k x ℓ w, u, v, w ∈ L p ((0, T ) × Ω) ;
(iii) equations in (S) are satisfied a.e (almost everywhere).

By weak solution to (S) on Q T = (0, T ) × Ω, we essentially mean solution in the sense of distributions or, equivalently here, solution in the sens of the variation of constants formula with the corresponding semigroups. More precisely

u(t) = S d 1 (t)u 0 + 0 S d 1 (t -s)(w γ (s) -u α (s)v β (s)) ds v(t) = S d 2 (t)v 0 + 0 S d 2 (t -s)(w γ (s) -u α (s)v β (s)) ds w(t) = S d 3 (t)u 0 + 0 S d 3 (t -s)(-w γ (s) + u α (s)v β (s)) ds where S d i (.) is the semigroup generated in L 1 (Ω) by -d i ∆ with homogeneous Neumann boun- dary condition, 1 ≤ i ≤ 3.
By just integrating the sum (E 1 ) + (E 2 ) + 2(E 3 ) in space and time, and taking into account the boundary conditions

Ω ∆(d 1 u + d 2 v + d 3 w) = 0 , we obtain Ω u(t) + v(t) + 2w(t) = Ω u 0 + v 0 + 2w 0 t ≥ 0. ( 1
)
Together with the nonnegativity of u, v and w, estimate (1) implies that

∀t ≥ 0 , u(t) L 1 (Ω) , v(t) L 1 (Ω) , w(t) L 1 (Ω) ≤ u 0 + v 0 + 2w 0 L 1 (Ω) . (2)
In other words, the total mass of three components does not blow up ; u(t), v(t) and w(t) rest bounded in L 1 (Ω) uniformly in time.

Although one has uniform L 1 -bound in time, classical solutions may not globally exist for diffusion coefficients d 1 , d 2 , d 3 which are not equal (global existence obviously holds if d 1 = d 2 = d 3 ). As surprisingly proved in [START_REF] Pierre | Blow-up in reaction-diffusion systems with dissipation of mass[END_REF] and [START_REF] Schmitt | Existence globale ou explosion pour des systèmes de réaction-diffusion avec contrôle de masse[END_REF], it may indeed happen that, under assumptions (P ) and (M ), solutions blow up in finite time in L ∞ ! In particular, classical bounded solutions do not exist globally in time.

If u 0 , v 0 , w 0 ∈ L ∞ (Ω), local existence and uniqueness of nonnegative and uniformly bounded solution to (S) are known (see e.g. [START_REF] Rothe | Global solutions of reaction-diffusion systems[END_REF]). More precisely, there exists T > 0 and a unique classical solution (u, v, w) of (S) on [0, T ). If T max denotes the greatest of these T's, then

T max < +∞ =⇒ lim tրTmax u(t) L ∞ (Ω) + v(t) L ∞ (Ω) + w(t) L ∞ (Ω) = +∞. ( 3
)
To prove global existence (i.e. T max = +∞), it is sufficient to obtain an a priori estimate of the form ∀t ∈ [0, T max ),

u(t) L ∞ (Ω) + v(t) L ∞ (Ω) + w(t) L ∞ (Ω) ≤ H(t), (4)
where H : [0, +∞) → [0, +∞) is a nondecreasing and continuous function.

This type of estimates is far of being obvious for our system except the case where diffusion coefficients

d 1 , d 2 , d 3 are equal i.e d 1 = d 2 = d 3 = d. Indeed, Z = u + v + 2w satisfies (E)      Z t -d∆Z = 0 (0, +∞) × Ω, ∂Z ∂n = 0 (0, +∞) × ∂Ω, Z(0, x) = Z 0 (x) x ∈ Ω,
where Z 0 (x) = u 0 (x) + v 0 (x) + 2w 0 (x).

In particular, we deduce by maximum principle that

u(t) + v(t) + 2w(t) L ∞ (Ω) ≤ u 0 + v 0 + 2w 0 L ∞ (Ω) , t ≥ 0.
Together with nonnegativity, this implies

u(t) L ∞ (Ω) + v(t) L ∞ (Ω) + w(t) L ∞ (Ω) ≤ u 0 + v 0 + 2w 0 L ∞ (Ω) , t ≥ 0.
In other words, u(t), v(t) and w(t) stay uniformly bounded in L ∞ (Ω) and therefore T max = +∞.

In the case where the diffusion coefficients are different from each other, global existence is considerably more complicated. It has been studied by several authors in the following cases.

First case α = β = γ = 1.
In this case, global existence of classical solutions has been obtained by Rothe [START_REF] Rothe | Global solutions of reaction-diffusion systems[END_REF] for dimension N ≤ 5. Later, it has first been proved by Pierre [10] for all dimensions N and then by Morgan [START_REF] Morgan | Global existence for semilinear parabolic systems[END_REF].

The exponentiel decay towards equilibrium has been studied by Desvillettes-Fellner [START_REF] Desvillettes | Exponential decay toward equilibrium via entropy methods for reaction-diffusion equations[END_REF] in the case of one space dimension.

The global existence of weak solutions has been proved by Laamri [START_REF] Laamri | Existence globale pour des systèmes de réaction-diffusion dans L 1[END_REF] for initial data u 0 , v 0 and w 0 only in L 1 (Ω). Second case γ = 1 regardless of α and β.

In this case, global existence of classical solutions has been obtained by Feng [START_REF] Feng | Coupled system of reaction-diffusion equations and Applications in carrier facilitated diffusion[END_REF] in all dimensions N and more general boundary conditions.

Third case α + β ≤ 2 or γ ≤ 2.
In this case, Pierre [START_REF] Pierre | Global Existence in Reaction-Diffusion Systems with Dissipation of Mass : a Survey[END_REF] has proved global existence of weak solutions for initial data u 0 , v 0 and w 0 only in L 2 (Ω).

Our paper mainly completes the investigations of [ [START_REF] Feng | Coupled system of reaction-diffusion equations and Applications in carrier facilitated diffusion[END_REF], [START_REF] Morgan | Global existence for semilinear parabolic systems[END_REF], [10], [START_REF] Rothe | Global solutions of reaction-diffusion systems[END_REF]] and [ [START_REF] Laamri | Existence globale pour des systèmes de réaction-diffusion dans L 1[END_REF], [START_REF] Pierre | Global Existence in Reaction-Diffusion Systems with Dissipation of Mass : a Survey[END_REF]]. As far as we know, our results are new either when α + β < γ, or when 1 < γ < N + 6 N + 2 regardless of α and For the sake of completeness and for the reader's convenience, we shall also give a direct proof different from that of Feng [START_REF] Feng | Coupled system of reaction-diffusion equations and Applications in carrier facilitated diffusion[END_REF] in the special case γ = 1.

Notation : Throughout this study, we denote by C i 's various positive numbers depending only on the data and for p ∈ [1, +∞[

u(t) p = Ω |u(t, x)| p dx 1/p , u L p (Q T ) = T 0 Ω |u(t, x)| p dtdx 1/p , u(t) ∞ = esse sup x∈Ω |u(t, x)|, u L ∞ (Q T ) = esse sup (t,x)∈Q T |u(t, x)|.

The main results

One of the main ingredients for the proof of our results is the following lemma which is based on the regularizing effects of the heat equation. This lemma has been introduced by Hollis-Martin-Pierre in [START_REF] Hollis | Global existence and boundedness in reaction-diffusion systems[END_REF].

Lemma 1 Let T > 0 and (φ, ψ) the classical solution of                      φ t -d 1 ∆φ = f (φ, ψ) (t, x) ∈ (0, T ) × Ω ψ t -d 2 ∆ψ = g(φ, ψ) (t, x) ∈ (0, T ) × Ω ∂φ ∂n (t, x) = 0 (t, x) ∈ (0, T ) × ∂Ω ∂ψ ∂n (t, x) = 0 (t, x) ∈ (0, T ) × ∂Ω φ(0, x) = φ 0 (x) x ∈ Ω ψ(0, x) = ψ 0 (x) x ∈ Ω.
Assume that f + g = 0, then for each p ∈ (1, +∞), there exists C such that for all t ∈ (0, T )

ψ L p (Qt) ≤ C φ L p (Qt) + 1 . (5)
A more general version of this lemma can be founded in [11, lemma 3.4].

The case α + β < γ

Theorem 1 Assume that 0 ≤ u 0 , v 0 , w 0 ≤ M where M is a positive real. If α + β < γ, then the system (S) admits a global classical solution.

Proof :

• Let T ∈ (0, T max ) and let t ∈ (0, T ]. Thanks to the nonnegativity of u, v and w, we deduce from the equation (E 1 ) that u is bounded from above by the solution U of (P 1 )

     U t -d 1 ∆U = w γ (t, x) ∈ (0, T ) × Ω ∂U ∂n (t, x) = 0 (t, x) ∈ (0, T ) × ∂Ω U (0, x) = u 0 (x) x ∈ Ω,
and we deduce from the equation (E 2 ) that v is bounded from above by the solution V of

(P 2 )      V t -d 2 ∆V = w γ (t, x) ∈ (0, T ) × Ω ∂V ∂n (t, x) = 0 (t, x) ∈ (0, T ) × ∂Ω V (0, x) = v 0 (x) x ∈ Ω.
Therefore it is sufficient to show that w ∈ L p (Q T ) for p large enough.

• Let q > 1. Multiplying the equation (E 3 ) by w q and integrating over Q T , we get

1 q + 1 Ω w q+1 (T ) + qd 3 Q T |∇w| 2 w q-1 + Q T w q+γ = Q T u α v β w q + K 0 (6)
where

K 0 = 1 q + 1 Ω w q+1 0 .
Thanks to Hölder's inequality, we have

Q T u α v β w q ≤ u α L αr (Q T ) v β L βs (Q T ) w q L γ+q (Q T ) (7)
where

1 r + 1 s + q q + γ = 1.
Since α + β < γ, we can choose r such that rα ≤ q + γ and s such that sβ ≤ q + γ. To convince oneself, it is enough to draw the straight line with cartesian equation x + y = γ q + γ and to identify the points with coordinates ( α q + γ , 0) and (0,

β q + γ ).
Then

L q+γ (Q T ) ⊂ L αr (Q T ) and L q+γ (Q T ) ⊂ L βs (Q T ). Consequently, there exists C 1 such that Q T u α v β w q ≤ C 1 u α L γ+q (Q T ) v β L γ+q (Q T ) w q L γ+q (Q T ) . (8)
By virtue of lemma 1, there exists C 2 such that

u L γ+q (Q T ) ≤ C 2 (1 + w L γ+q (Q T ) ) (9)
and there exists C 3 such that

v L γ+q (Q T ) ≤ C 3 (1 + w L γ+q (Q T ) ). ( 10
)
Thanks to ( 9) and (10), estimate (8) can be written

Q T u α v β w q ≤ C 4 1 + w L γ+q (Q T ) α 1 + w L γ+q (Q T ) β 1 + w L γ+q (Q T ) q . ( 11
)
If w L γ+q (Q T ) ≤ 1 then the proof ends up. Otherwise, there exists C 5 such that

Q T u α v β w q ≤ C 5 w q+α+β L γ+q (Q T ) . (12)
So we deduce from (6)

Q T w q+γ ≤ C 5 w q+α+β L γ+q (Q T ) + K 0 . (13)
With the notation R := Q T w q+γ , estimate (13) can be written

R ≤ C 5 R q+α+β q+γ + K 0 . (14)
Since q + α + β < q + γ, by applying Young's inequality to [START_REF] Quittner | Souplet : Superlinear Parabolic Problems : Blow-up, Global Existence and Steady States, Advanced Texts[END_REF], we obtain

(1 -ε)R ≤ K 0 + C 6 . (15)
Then, for ε ∈ (0, 1), we have the desired estimate

w L q+γ (Q T ) ≤ C 7 . (16)
Going back to (P 1 ) and (P 2 ), we have, by choosing q such that q + γ γ > N + 2 2 and thanks to the L p -regularity theory for the heat operator (see [START_REF] Ladyzenskaya | Linear and quasilinear equations of parabolic type[END_REF]),

u L ∞ (Q T ) ≤ C 8 (17) v L ∞ (Q T ) ≤ C 9 . (18)
Now going back to (E 3 ), we deduce from ( 17) and ( 18) that there exists C 10 such that

w L ∞ (Q T ) ≤ C 10 . (19)
This implies that T max = +∞.

Remark This method seems to be specific to the case α + β < γ. It fails when α + β ≥ γ since some restrictions on the parameters α, β, γ and on the diffusion coefficients will appear.

Case where

d 1 = d 3 or d 2 = d 3 or d 1 = d 2 .
Theorem 2 Assume that 0 ≤ u 0 , v 0 , w 0 ≤ M . We deduce by maximum principle

u(t) + w(t) ∞ ≤ u 0 + w 0 ∞ . (20)
Together with the nonnegativity of u et w, this implies that u(t) and w(t) are uniformly bounded in L ∞ (Ω).

By going back to (E 2 ) and thanks to the L p -regularity theory for the heat operator (see [START_REF] Ladyzenskaya | Linear and quasilinear equations of parabolic type[END_REF]), we conclude that v(t) ∞ is uniformly bounded in L ∞ (Ω) on all interval [0, T ] so that T max = +∞.

(ii) Assume that d 1 = d 2 = d. The case α + β < γ was already handled in the theorem 1, so it remains only to tackle the case γ < α + β. Moreover, one can assume that u 0 = v 0 since if u 0 = v 0 the result is obvious.

Since

d 1 = d 2 = d, we have (u -v) t -d∆(u -v) = 0 ; ∂(u -v) ∂n = 0 ; (u -v)(0, x) = u 0 (x) -v 0 (x).
The maximum principle then implies u(t) -v(t) ∞ ≤ u 0 -v 0 ∞ = C. Hence we have

u α+β = u α v β + u α (u β -v β ) = u α v β + u α β(θu + (1 -θ)v) β-1 (u -v) where θ ∈]0, 1[ ≤ u α v β + u α β2 β-1 C(u β-1 + v β-1 ).
Thanks to Young's inequality, there exists C 11 > 0 and C 12 > 0 such that

C 11 u α+β ≤ u α v β + C 12 . ( 21
)
By virtue of (21), equation (E 1 ) implies that

u t -d 1 ∆u + C 11 u α+β ≤ w γ + C 12 . ( 22
)
Let q > 1. Multiplying ( 22) by u q and integrating over Q T , we obtain

1 q + 1 Ω u q+1 (T )+qd 2 Q T |∇u| 2 u q-1 +C 11 Q T u q+α+β ≤ Q T w γ u q +C 12 Q T u q +K 1 (23) where K 1 = 1 q + 1 Ω u q+1 0 .
Thanks to Hölder's inequality, we have

Q T w γ u q ≤ Q T w γr 1/r Q T u qs 1/s (24)
where r = α + β + q γ and s = α + β + q q + α + β -γ .

Lemma 1 implies that there exists C 13 such that

Q T w γr 1/r = w γ L q+α+β (Q T ) ≤ C γ 13 1 + u L q+α+β (Q T ) γ .
If u L q+α+β (Q T ) ≤ 1 then the proof ends up. Otherwise, there exists C 14 such that

Q T w γr 1/r ≤ C 14 u γ L q+α+β (Q T ) . (25)
Since qs < q + α + β, we have L q+α+β (Q T ) ⊂ L qs (Q T ), then there exists C 15 such that

Q T u qs 1/s ≤ C 15 u q L q+α+β (Q T ) . (26)
Denote S := Q T u q+α+β . Estimates (25) and (26) imply that

Q T w γ u q ≤ C 16 S q+γ q+α+β . (27)
Moreover, since L q+α+β (Q T ) ⊂ L q (Q T ), there exists C 17 such that

C 12 Q T u q ≤ C 17 S q q+α+β . (28)
Since γ < α + β, by applying Young's inequality, there exists C 18 such that

C 16 S q+γ q+α+β ≤ ε 2 S + C 18 . (29)
Applying again Young's inequality, there exists C 19 such that

C 17 S q q+α+β ≤ ε 2 S + C 19 . (30)
Consequently, estimate (23) implies

(C 11 -ε)S ≤ C 18 + C 19 + K 1 . (31)
By choosing ε < C 11 in (31), there exists C 20 such that

u L q+α+β (Q T ) ≤ C 20 . (32)
Thanks to lemma 1 and estimate (32) there exists C 21 such that

w L q+α+β (Q T ) ≤ C 21 . (33)
By going back to (P 1 ) and (P 2 ), we have by choosing q such that q + α + β γ > N + 2 2 and thanks to the L p -regularity theory for the heat operator (see [START_REF] Ladyzenskaya | Linear and quasilinear equations of parabolic type[END_REF])

u L ∞ (Q T ) ≤ C 22 (34) v L ∞ (Q T ) ≤ C 23 . (35)
Now let's go back to (E 3 ), we deduce from (34) and (35) that there exists C 24 such that

w L ∞ (Q T ) ≤ C 24 . (36)
This implies that T max = +∞. Remark Even in the last case i.e d 1 = d 2 , global existence or blow-up in the limit case α+β = γ remain an open problem.

Case

1 ≤ γ < N + 6 N + 2
regardless of α and β. Proof : Let T ∈ (0, T max ) and let t ∈ (0, T ]. Thanks to the nonnegativity of u, v and w, we deduce from the equation (E 1 ) that u is bounded from above by the solution U of (P 1 )

     U t -d 1 ∆U = w γ (t, x) ∈ (0, T ) × Ω ∂U ∂n (t, x) = 0 (t, x) ∈ (0, T ) × ∂Ω U (0, x) = u 0 (x) x ∈ Ω.
Therefore it is sufficient to show that w ∈ L p (Q T ) for p large enough. For this we have to distinguish the case γ = 1 and the case γ > 1.

• Case γ = 1 and α, β ≥ 1.

Let us recall that global existence of classical solutions for (S) when α = β = γ = 1 has been studied by several authors. It has been obtained by Rothe [START_REF] Rothe | Global solutions of reaction-diffusion systems[END_REF] for dimension N ≤ 5. Later, it has first been proved by Pierre [10] for all dimensions N and then by Morgan [START_REF] Morgan | Global existence for semilinear parabolic systems[END_REF]. Independantly, Feng [START_REF] Feng | Coupled system of reaction-diffusion equations and Applications in carrier facilitated diffusion[END_REF] has proved global existence in the case γ = 1 regardless of α and β and more general boundary conditions. For the sake of completeness and for the reader's convenience, we give here a simple and direct proof in the last case (γ = 1 regardless of α and β). In our proof, we use an idea introduced by Pierre in [10] and applied in [START_REF] Martin | Nonlinear reaction-diffusion systems[END_REF].

For any p ≥ 1, we deduce from (P 1 ) and the semigroup property

u(t) p ≤ u 0 p + t 0 w(s) p ds. (37)
By applying Hölder's inequality for p > 1 and thanks to (5), we obtain

Taking the p th power of (39) we obtain

h(t) ≤ 2 p-1 C p 26 + 2 p-1 C p 27 t 0 h(s) ds. (40)
But, inequality (40) is a linear Gronwall's inequality, then

u L p (Q T ) ≤ C 28 . (41)
Repeating the method above with v instead of u, we obtain

v L p (Q T ) ≤ C 29 . (42)
Estimates ( 41) and (42) imply that for some q > N + 2 2

u α v β L q (Q T ) ≤ C 30 . (43)
Going back to equation (E 3 ) we have, thanks to the L q -regularity theory for the heat operator,

w L ∞ (Q T ) ≤ C 31 . (44)
This concludes the proof for the case γ = 1 regardless of α and β.

• Case 1 < γ < N + 6 N + 2 .
The proof in this case is based on lemma 1 and these two following lemmas.

Lemma 2 (Michel Pierre) Let T > 0 and let Z the solution of

     Z t -∆(A(t, x)Z) ≤ 0 (t, x) ∈ (0, T ) × Ω, ∂Z ∂n (t, x) = 0 (t, x) ∈ (0, T ) × ∂Ω, Z(0, x) = Z 0 (x) x ∈ Ω. Assume that 0 < d < A(t, x) < D where (d, D) ∈ (0, +∞) 2 . Then, there exists C = C(T, d, D, Ω) such that Z L 2 (Q T ) ≤ C Z 0 L 2 (Ω) .
For a general version of this lemma, see [11, proposition 6.1] or [3, theorem 3.1].

Lemma 3 Let (p, q) such that 1 ≤ p ≤ q ≤ +∞, d > 0 and S d (t) the semigroup generated in L p (Ω) by -d∆ with homogeneous Neumann boundary condition. Then

S d (t)Y q ≤ (C(Ω)m(t)) -N 2 ( 1 p -1 q ) Y p , for all Y ∈ L p (Ω), t > 0 ( 45
)
where m(t) = min(1, t).

For a proof of this lemma see for instance [START_REF] Rothe | Global solutions of reaction-diffusion systems[END_REF]Lemma 3,p. 25] or [1, Theorem 3.2.9, p. 90].

We now go back to the proof of theorem 3.

By applying lemma 2 to the system (S) where Z = u + v + 2w and A = d 1 u + d 2 v + 2d 3 w u + v + 2w , we have u, v, w ∈ L 2 (Q T ). More precisely, there exists C 32 such that • In the case where the diffusion coefficients are not equal (i.e. d i = d j for all 1 ≤ i = j ≤ 3), global existence of classical solutions for (S) or blow-up is still an open question when

u L 2 (Q T ) , v L 2 (Q T ) , w L 2 (Q T ) ≤ C 32 . (46)
N + 6 N + 2 ≤ γ ≤ α + β.
Our guess is that system (S) admits a classical global solution for all N + 6 N + 2 ≤ γ < α + β and that there is a finite time blow-up when γ = α + β and the dimension N is large.

β.

  For the sake of clarity, we decided to focus in this work on the question of global existence in time of solutions in the case of homogeneous Neumann boundary conditions. So, we shall prove global existence of classical solutions to system (S) in the following cases : * α + β < γ ; * (d 1 = d 3 or d 2 = d 3 ) and for any (α, β, γ) ; * d 1 = d 2 and for any (α, β, γ) such that α + β = γ ; * 1 < γ < N + 6 N + 2 and for any (α, β).

  (i) If d 1 = d 3 or d 2 = d 3 , then system (S) admits a global classical solution for any (α, β, γ). (ii) If d 1 = d 2 , then the system (S) admits a global classical solution for any (α, β, γ) such that α + β = γ. Proof : (i) Assume that d 1 = d 3 = d, we have (u + w) t -d∆(u + w) = 0 ; ∂(u + w) ∂n = 0 ; (u + w)(0, x) = u 0 (x) + w 0 (x).

Theorem 3

 3 Assume that 0 ≤ u 0 , v 0 , w 0 ≤ M where M > 0. If 1 ≤ γ < N + 6 N + 2, then the system (S) admits a global classical solution for any (α, β) ∈ [1, +∞) 2 .

  Now, we have thanks to the estimate (45) with p > 1 and q = +∞ u(t) ∞ ≤ u 0 ∞ + C 33 All our results are still true if we replace homogeneous Neumann boundary conditions by homogeneous Dirichlet boundary conditions, it suffices to replace lemma 3 by the following one.Lemma 4 Let (p, q) such that 1 ≤ p ≤ q ≤ +∞, d > 0 and S d (t) the semigroup generated in L p (Ω) by -d∆ with homogeneous Dirichlet boundary. ThenS d (t)Y q ≤ (4πt) ) Y p , for all Y ∈ L p (Ω), t > 0. (55)For a proof of this lemma, see for instance[START_REF] Quittner | Souplet : Superlinear Parabolic Problems : Blow-up, Global Existence and Steady States, Advanced Texts[END_REF] Proposition 48.4, p. 441].

	3 Conclusion
	• -N 2 ( 1 p -1
							t	(t -s)	-N 2p	w γ (s) p ds.			(47)
							0
	By applying Hölder's inequality, we obtain
	0	t	(t -s)	-N 2p	w(s) γ	p ds ≤	-Np ′ 2p	ds	1/p ′	0	t	w γ (s) p p ds	1/p	.	(48)

t 0 (t -s) q

Acknowledgments : I am indebted to Michel Pierre for numerous and helpful discussions concerning this work. I would like to thank warmly Philippe Souplet for his advice which improved the second part of the proof of theorem 3. I am sincerely grateful to Laurent Desvillettes for his careful reading of a initial version of this paper and for his useful comments. I thank Didier Schmitt for many conversations. Finally, I also thank the referee for valuable comments and remarks which allow me to improve the writing of this paper.

We first remark that the integral t 0 (t -s) -N 2(p-1) ds converges when p > N + 2 2 and we have 1) dy.

On the other hand, lemma 1 implies that t 0 w γ (s) p p ds

If u L pγ (Qt) ≤ 1 then the proof ends up. Otherwise there exists C 35 such that

Since

, it follows that (47) can be written

If p(γ -1) < 2, by choosing ε ∈ (0, min(p, 2 -p(γ -1)), we deduce from ( 46) and (51) that there exists

Note that the above condition p(γ -1)

We establish in the same way that there exists C 38 such that

Finally, for (E 3 ), we deduce from ( 52) and ( 53) that there exists C 39 such that

This concludes the proof in the case 1 < γ < N + 6 N + 2 .

Remark : Our conjecture is that γ * = N + 6 N + 2 is not optimal. In fact, when N = 1 one can prove that the result of theorem 3 still holds for γ * = 7/2.


Christian Obrecht

Frédéric Kuznik

Bernard Tourancheau

Jean-Jacques Roux

A new approach to the lattice Boltzmann method for graphics processing units

Keywords: GPU programming, CUDA, Lattice Boltzmann method, Parallel computing

Emerging many-core processors, like CUDA capable nVidia GPUs, are promising platforms for regular parallel algorithms such as the Lattice Boltzmann Method (LBM). Since global memory on graphic devices shows high latency and LBM is data intensive, memory access pattern is an important issue to achieve good performances. Whenever possible, global memory loads and stores should be coalescent and aligned, but the propagation phase in LBM can lead to frequent misaligned memory accesses. Most previous CUDA implementations of 3D LBM faced this problem using low latency on chip shared memory. Instead of this, our CUDA implementation of LBM follows carefully chosen data transfer schemes in global memory. On the 3D lid-driven cavity test case, we obtained 2× to 3× speed-up over previously published performances, achieving up to 86% of the global memory maximal throughput. As a consequence we show that highly efficient implementations of LBM on GPU are possible, even for complex models.

Introduction

During the last decade, the computational power of commodity graphics hardware has dramatically increased, as shown in figure 1, nearing 1 GFlop/s with nVidia's latest GT200. Yet, one should be aware that this performance is attainable only for single precision computations, which are not fully IEEE-754 compliant. Nonetheless, due to their low cost, GPUs become more and more popular for scientific computations (see [START_REF] Dongarra | Exploring new architectures in accelerating CFD for Air Force applications[END_REF][START_REF] Tomov | Towards dense linear algebra for hybrid GPU accelerated manycore systems[END_REF] ).

Lattice Boltzmann method, which originates from the lattice gas automata methods, is an efficient alternative to the numerical solving of Navier-Stokes equations for simulations of complex fluid systems. Besides its numerical stability and accuracy, one of the major advantage of LBM is its data parallel nature. Nevertheless, using LBM for practical purposes requires large computational power. Thus, several attempts to implement LBM on GPUs were made recently.

In this paper, we intend to present some optimisation principles for CUDA programming. These * Corresponding author: christian.obrecht@insa-lyon.fr principles led us to a GPU implementation of 3D LBM which appears to be more efficient than the previously published ones.

CUDA

The Compute Unified Device Architecture (CUDA), released by nVidia in early 2007, is up to now the leading technology for general purpose GPU programming (see [START_REF] Halfhill | Parallel processing with CUDA[END_REF] ). It consists of hardware specifications, a specific programming model, and a programming environment (API and SDK).

Architecture

General purpose GPU programming usually requires to take some architectural aspects into consideration. CUDA hardware specifications make the optimisation process easier by providing a general model for the nVidia GPUs architecture from the G80 generation on. Figure 2 shows the main aspects of the CUDA hardware specifications. A GPU consists in several Streaming Multiprocessors (SMs). Each SM contains Scalar Processors (SPs), an instruction unit, and a shared memory, concurrently accessible by the SPs through 16 memory banks. Two cached, read-only memories for constants and textures are also available. The device memory, usually named global memory is accessible by both the GPU and the CPU. Table 1 specifies some of the features of the GT200 processor on which our implementations were tested.

SPs are only able to perform single precision computations. From compute capability 1.3 on, CUDA supports double precision. On this kind of hardware, each SM is linked to a double precision computation unit. Both single and double precision calculations are mostly IEEE-754 compliant. Divergences from the standard are mainly: • Partial support of rounding modes.

• No floating point exception mechanism.

• Multiply-add operations with truncated intermediate results.

• Non compliant implementations of some operations like division or square root.

Programming

CUDA programming model (see [START_REF] Nvidia | Compute Unified Device Architecture Programming Guide version 2[END_REF] ) relies on the concept of kernel. A kernel is a function that is executed in concurrent threads on the GPU. Threads are grouped into blocks which in turn form the execution grid (see figure 3). The CUDA technology makes use of a slightly modified version of the C (or C++) language as a 2 programming language. The code of a CUDA application consists in functions which can be classified in four categories:

1. Sequential functions run by the CPU. 2. Launching functions allowing the CPU to start a kernel. 3. Kernels run by the GPU. [START_REF] Nvidia | Compute Unified Device Architecture Programming Guide version 2[END_REF]. Auxiliary functions which are inlined into the kernels at compile time.

The execution grid's layout is specified at run time. A grid may have one or two dimensions. The blocks of threads within a grid must be identical and may have up to three dimensions. A thread is identified in respect of the grid using the two structures threadIdx and blockIdx, containing the three fields x, y, and z.

A block may only be executed on a single SM, which yields to an upper bound of the number of threads within a block 1 . Scheduling is carried out at hardware level and may not be adjusted. It is yet possible to place synchronisation barriers, but their scope is limited to blocks. The only way to ensure global synchronisation is to use a kernel for each step.

The local variables of a kernel are stored in the registers of the SMs. Their scope is limited to threads and they cannot be used for communication purposes. Data exchanges between threads require the use of shared memory. The management of these exchanges is left to the programmer. It is worth noting that no protection mechanism is available, hence concurrent writes at the same memory location yield unpredictable results. The shared memory's scope is limited to blocks. Communication between threads belonging to different blocks requires the use of global memory.

Optimisation principles

Computational aspect

Generally speaking, the occupancy rate of the SPs, i.e. the ratio between the number of threads run and the maximum number of executable threads, is an important aspect to take into consideration for the optimisation of a CUDA kernel. Even though a block may only be run on a single SM, it is possible however to execute several 1 As of compute capability 1.3, this maximum is 1,024.

blocks concurrently on the same SM. Hence tuning the execution grid's layout allows to increase the occupancy rate. Nevertheless, reaching the maximal occupancy is usually not possible: the threads executed in parallel on one SM have to share the available registers. On compute capability 1.3 architectures, for instance, maximal occupancy is achieved only for kernels using at most 16 registers, that is to say only the simplest ones. It should be noted that shared memory, which is rather scarce too, may also be a limiting factor for the occupancy.

The rather elementary optimisation technique consisting in common sub-expression elimination should be used with care. As a matter of fact, this method implies to store the values of these subexpressions in temporary variables, thence increasing the use of registers, which in turn may lead to lower occupancy. In some cases, this common sense technique has negative effects, and it may be better to recompute some values than to store them. Anyway, general principles regarding this topic are not relevant. Since the compiler performs aggressive optimisation, the number of register needed for a given kernel is scarcely predictable.

The hardware scheduler groups threads in warps of 32 threads. Though not mandatory, the number of threads in a block should be a multiple of the warp size. Whenever a warp is running, all the corresponding threads are executed concurrently by the SPs, except when conditional branching occurs. Divergent branches are executed sequentially by the SM. Even though serialisation only happens at warp level, conditional structures should be avoided as much as possible, being likely to have a major impact on actual parallelism.

Regarding optimisation, the cost of arithmetic operations (in clock cycles) must also be taken in consideration. Table 2 It should be noted that, since addition and mul-tiplication are merged in one single axpy operation whenever possible, evaluating the actual algorithmic complexity of a computation is not straightforward. It's also worth noting that division is rather expensive and should be used parsimoniously.

Data transfer aspect

For many applications, memory transactions optimisation appears to be even more important than computations optimisation. Registers do not arise any specific problem apart from their limited amount. Shared memory is in terms of speed similar to register but is accessed by the SPs through 16 memory banks. For efficient accesses, each thread in a half-warp must use a different bank. When this condition is not met, the transaction is repeated as much as necessary.

Global memory, being the only one accessible by both the CPU and the GPU, is a critical path for CUDA applications. Unlike registers and shared memory, global memory suffers high latency ranging from 400 to 600 clock cycles. Nonetheless, this latency can be mostly hidden by the scheduler which stalls inactive warps until data is available. Furthermore, global memory throughput is significantly less than register throughput. For data intensive applications like LBM, this aspect is generally the limiting factor.

Global memory accesses are performed by halfwarp on 32, 64, or 128 bytes segments whose start addresses are multiple of the segment's size. To optimise global memory transactions, memory accesses should be coalesced and aligned whenever possible. To achieve coalescence, threads within a half-warp must access contiguous memory locations.

Data transfer modelling

In CUDA applications, the execution of a kernel can generally be split into three steps:

1. Reading data from global memory. 2. Processing data using registers (and possibly shared memory). 3. Writing processed data to global memory. Code 1 follows this scheme in the case where the amount of data read and written are equal. Function launch_kernel calls function kernel with an execution grid containing L 3 threads. One may notice some syntactic specificities of the CUDA programming language: the use of the tripled angle brackets for kernel invocation, the __global__ keyword for kernel definition, the __device__ keyword for auxiliary functions. The kernel performs the reading and writing of N 32-bit words. The second step is simplistic, though not suppressed in order to ensure actual data transfer to the GPU. Global memory accesses are optimal, provided L is set to an appropriate value. In the present study, L = 128 was chosen. Measuring the execution time of the launch_kernel function enables to estimate the average time T for data transfer between GPU and global memory relatively to the amount N of data exchanged. Figure 4 shows the results for one warp with T in nanoseconds and N in 32-bit words, obtained using a GeForce GTX 295 graphics board.

The quasi-linear aspect of these measurements reveals that, in ideal cases, the hardware scheduler is able to hide the latency of global memory. Numerically, we obtain:

T ≈ 2.78 × N + 0.99 (1)
The average throughput is almost constant relatively to N and is about 90.7 GB/s for the GeForce GTX 295. Reckoning the characteristics of the benchmark program, we consider the obtained value as the effective maximal throughput for data transfer between GPU and global memory. This upper bound is useful in evaluating the performances of a CUDA application leaving aside the hardware in use. The obtained value is about 81% of the theoretical maximal throughput, which is comparable to the result found in [START_REF] Tölke | Implementation of a Lattice Boltzmann kernel using the Compute Unified Device Architecture developed by nVIDIA[END_REF].

On the same hardware, the bandwidthTest program from the CUDA SDK gives 91.3 GB/s, which is rather close to the value we obtained. Yet, this program uses only memory copy functions instead of a kernel, hence yielding less relevant results from a practical standpoint.

Lattice Boltzmann Method

The lattice Boltzmann method is based on a threefold discretisation of the Boltzmann equation: time, space and velocity (see [START_REF] Mcnamara | Use of the Boltzmann Equation to Simulate Lattice-Gas Automata[END_REF] ). Velocity space reduces to a finite set of well chosen velocities {e i | i = 0, . . . N} where e 0 = 0. Figure 5 illustrates the D3Q19 stencil we used.

Instead of reviewing the well-known Lattice Bhatnagar-Gross-Krook (LBGK) model (see [START_REF] Qian | Lattice BGK models for Navier-Stokes equation[END_REF] ), we will outline the Multiple Relaxation Time model presented in [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]. The analogous of the one-particule

f (x, t) = (f 0 (x, t), . . . f N (x, t))
T for given lattice node x and time t, T being the transpose operator. The mass fractions can be mapped to a set of moments {m i | i = 0, . . . N} by an invertible matrix M such as:

f (x, t) = M -1 m(x, t) (2)
where m(x, t) is the moment vector. With the D3Q19 stencil, the density is ρ = m 0 , the momentum is j = (m 3 , m 5 , m 7 ). Higher order moments as well as matrix M are given in detail in [8, app. A]. Using these notations, the lattice Boltzmann equation can be written as:

f (x + δte i , t + δt) -f (x, t) = M -1 S m(x, t) -m (eq) (x, t) (3)
where m (eq) (x, t) is the equilibrium-moment vector and:

S = diag(s 0 , . . . s N )
is the relaxation rates matrix. The LBGK model is a special case of MRT where all relaxation rates s i = 1/τ . In a numerical point of view, MRT should be preferred to LBGK, being more stable and accurate.

Previous Work

Data organisation schemes for LBM are mainly of two kinds. First, the Array of Structures (AoS) type, which for D3Q19 is equivalent to a L 3 × 19 array. Second, the Structure of Arrays (SoA) type, which for D3Q19 is equivalent to a 19 × L 3 array. For CPU implementations of the LBM, the AoS is relevant insofar as it improves the locality of mass fractions associated to a same node (see [START_REF] Pohl | Optimization and Profiling of the Cache Performance of Parallel Lattice Boltzmann Codes[END_REF] ). Up to now, for all GPU implementations of LBM, a thread is allocated to each lattice node, which is probably the simplest way to take advantage of the massively parallel structure of the GPU. With this approach, ensuring coalescence of global memory accesses requires to use a SoA kind of organisation.

With values of L divisible by 16, every mass fractions associated to a half-wrap lay in a same segment of global memory. Yet, this is not sufficient to ensure optimal memory transaction. Indeed, for the minor spatial dimension, propagation corresponds to one unit shifts of memory addresses. In other words, for most mass fractions, propagation phase leads to misalignments. Getting round this problem was up to now the main issue regarding GPU implementations of the LBM.

The first attempt of implementing a D3Q19 model using CUDA is due to Ryoo et al. (see [START_REF] Ryoo | Optimization principles and application performance evaluation of a multithreaded GPU using CUDA[END_REF] ). It consists mainly in a port of the 470.lbm code from the SPEC CPU2006 benchmark (see [START_REF] Henning | SPEC CPU2006 Benchmark Descriptions[END_REF] ). In terms of optimisation, switching from AoS to SoA is the only important modification undertaken. To the best of our knowledge, misalignment problems caused by propagation are not taken into consideration. The announced speed-up factor of 12.3 is rather low compared to subsequent results.

The two-dimensional D2Q9 implantation submitted by Tölke in [START_REF] Tölke | Implementation of a Lattice Boltzmann kernel using the Compute Unified Device Architecture developed by nVIDIA[END_REF] solves the misalignment problems using one-dimensional blocks and shared memory. More precisely, propagation within one block is split in two steps: a longitudinal shift in shared memory followed by a lateral shift in global memory. This approach is outlined in figure 6.

Since blocks follow the minor dimension, no more misalignment arises. Nevertheless, because of the limited scope of shared memory, mass fractions leaving or entering a block require specific handling. The retained solution is to store out-coming mass fractions in places temporarily left vacant by in coming mass fractions (see figure 7).

A drawback of the shared memory approach is consequently the need for a second kernel exchang- ing data in order to place properly mass fractions located at the blocks' boundaries. Obviously, this further processing has a non negligible cost.

Following the same method than Tölke, Habich in [START_REF] Habich | Performance Evaluation of Numeric Compute Kernels on nVIDIA GPUs[END_REF] describes an implementation of a D3Q19 model. The transition from D2Q9 to D3Q19 leads to lower performances, achieving only 51% of the effective maximal throughput. Habich assumes this decrease is due to the low occupancy rate. As a matter of fact, given the limited amount of registers, the number of threads run in parallel on a SP cannot exceed one or two warps.

This point of view is probably erroneous as we shall see subsequently. The lower performances seem more likely due to the increase of the execution time of the second kernel, since there is dramatically more data to exchange than for twodimensional LBM. As an example, for D2Q9 on a 2,048 2 grid with 128 threads per block, there are 2,048×16×6 = 196,608 mass fractions to move. For D3Q19, on a 160 3 grid with 32 threads per block, there are 160 2 × 5 × 14 = 1,792,000 mass fractions to move. Relatively to the number of nodes, the ratio is about 9.3.

A way to obtain better performances for threedimensional LBM consists in using stencils contain-6 ing less mass fractions, like D3Q13. This approach was studied by Tölke and Krafczyk in [START_REF] Tölke | TeraFLOP computing on a desktop PC with GPUs for 3D CFD[END_REF], obtaining 61% of the effective maximal throughput. The D3Q13 stencil, which corresponds to the points of contact in a close-packing of spheres, is the simplest three-dimensional structure sufficiently isotropic for LBM. Yet, due to the lesser amount of information processed, D3Q13 is less accurate than D3Q19. Furthermore, node addressing becomes quite complex.

Bailey et al. in [START_REF] Bailey | Accelerating Lattice Boltzmann Fluid Flow Simulations Using Graphics Processors[END_REF] announce a 20% improvement of maximal performances for their implementation of D3Q19 compared to those published in [START_REF] Habich | Performance Evaluation of Numeric Compute Kernels on nVIDIA GPUs[END_REF]. The description of the tested optimisations is not very explicit, but it seems that the main intention was to increase occupancy. One of the proposed technique consists in imposing at compile time an upper bound to the number of registers used by the computation kernel. Of course this directive causes the compiler to fall back on register spilling. Taking the cost of global memory accesses into account, we consider this approach as not relevant.

Proposed Implementations

All but one CUDA implementations of LBM mentioned in the former section use shared memory for propagation. As formerly outlined, this approach imposes the use of a second kernel taking care of the mass fractions crossing the blocks' boundaries. Though rather basic, the CUDA profiler allows to gather some informations during kernel run time (see [START_REF] Nvidia | [END_REF] ). Concerning LBM, this tool led us to make two assumptions:

1. The additional cost caused by misalignment has the same order of magnitude than the one caused by the exchange kernel. 2. The cost of a misaligned read is less than the cost of a misaligned write.

Hence we adopted the following approach for our implementations of D3Q19:

• SoA type of data organisation.

• One-dimensional blocks following the minor dimension.

• Propagation performed by global memory transactions.

• Deferment of misalignment on reading.

We experimented two propagation schemes: a split scheme and a reversed scheme. The split scheme was tested with a LBGK model and ongrid boundary conditions. The reversed scheme was tested with a MRT model and mid-grid boundary conditions. To ease cross platform development, we employed the CMake build system (see [START_REF] Martin | Mastering CMake, A Cross-Platform Build System[END_REF] ). Moreover, we used the XML based VTK format for output (see [START_REF] Schroeder | The VTK User's Guide[END_REF] ).

Split scheme

With the split scheme, propagation is parted in two components: shifts that induce misalignment are performed at reading, the others are performed at writing, as outlined in figure 8. For the sake of simplicity, the diagram shows the two-dimensional case. Boundary conditions are implemented using ongrid bounce back: nodes at the cavity's borders, except the lid, are considered as solid and simply return the in-coming mass fractions in the opposite direction.

To summarise, the corresponding kernel breaks up into:

1. Reading along with propagation in minor dimension. 2. On-grid bounce back boundary conditions. 3. Computations using LBGK model. 4. Writing along with propagation in major dimensions.

Reversed scheme

With reversed scheme, propagation is entirely performed at reading, as outlined in figure 9. Again, the diagram shows the two-dimensional case only. Boundary conditions are implemented using midgrid bounce back: nodes at the cavity's borders, except the lid, are considered as fluid with null velocity. Unknown mass fractions are determined using:

f i -f eq i = f j -f eq j ( 4
)
whith i and j such as e i = -e j (see [START_REF] Zou | On pressure and velocity flow boundary conditions and bounceback for the lattice Boltzmann BGK model[END_REF] ). For null velocity, f eq i = f eq j . Hence the former equation yields f i = f j .

To summarise, the corresponding kernel breaks up into:

1. Reading carrying out propagation. 2. Mid-grid bounce back boundary conditions. 3. Computations using MRT model. 4. Writing without propagation.

Results

Validation

Numerical validation is an important issue in GPU computing, since most calculations are performed using single precision. This topic being thoroughly studied for GPU implementations of LBM in [START_REF] Kuznik | LBM Based Flow Simulation Using GPU Computing Processor[END_REF], we will rather focus on physical validation. We used the well-known lid-driven cavity test case, comparing velocity coordinates with results published by Albensoeder and Kuhlmann in [START_REF] Albensoeder | Accurate threedimensional lid-driven cavity flow[END_REF]. More precisely, for velocity u(u, v, w), we gathered u on line x = L/2 and y = L/2, as well as w on line z = L/2 and y = L/2. For Reynolds number Re = 1,000, our LBGK code outcomes are in quite good accordance with the reference values. Not surprisingly, the MRT implementation achieves almost perfect correspondence as shown in figure 10.

Performances

Binaries for nVidia GPUs are generated through a two stages process (see [START_REF] Murphy | NVIDIA's Experience with Open64[END_REF] ). First, the nvopencc program compiles CUDA code into Parallel Thread eXecution (PTX) pseudo assembly language. Second, the ocg assembler translates the PTX code into actual binary. Analysing PTX outputs allows to enumerate the floating point operations in a kernel and thence to evaluate the actual algorithmic complexity of the computations. Table 3 Though being more complex than LBGK, MRT has almost the same computational cost. It is worth noting that this cost is of the same order of magnitude than one single global memory transaction, that is to say 400 to 600 cycles. Thus, taking the hardware scheduler into account, the impact of computations on global processing time is negligible. Most of the execution time of our kernels is consumed by data transfer, the remaining being probably induced by scheduling. In terms of optimisation, increasing the occupancy rate of the SMs is not especially crucial.

The former opinion is supported by the analysis of the performances of our implementations. Million Lattice node Updates Per Second (MLUPS) is the usual unit for performance measurement in LBM. For both implementations, memory addressing is analogous to the one used in code 1. Therefore, the size of the blocks corresponds to the size of the cavity. Tables 4 and5 show the obtained performances on a Debian GNU/Linux 5.0 workstation fitted with a GeForce GTX 295 graphics card. One may notice that the data transfer rate is rather close to maximum. Global memory throughput is presently the limiting factor for LBM computations on GPU. Moreover, it is worth mentioning that these satisfactory performances are achieved with quite low SM occupancy.

Confronting the obtained performances to published results corroborates our approach. Depending on the size of the cavity, we observe 2× to 3× speed-up compared to the performances mentioned in [START_REF] Habich | Performance Evaluation of Numeric Compute Kernels on nVIDIA GPUs[END_REF][START_REF] Bailey | Accelerating Lattice Boltzmann Fluid Flow Simulations Using Graphics Processors[END_REF]. Yet these studies were led on GeForce 8800 GTX graphics cards, which belong to the previous generation though being comparable to the hardware we used in terms of memory throughput. Therefore, these comparisons should be considered with care, and we additionally compared a D2Q9 version of our code to the one published in [START_REF] Tölke | Implementation of a Lattice Boltzmann kernel using the Compute Unified Device Architecture developed by nVIDIA[END_REF] on the GTX 295 obtaining a 15% betterment of the performances.

Summary

The present study proposes a model for data transfer on the latest generation of nVidia GPUs. Optimisation principles, leading to efficient implementations of 3D LBM on GPUs, are drawn as well. We state the impact of global memory transfer as the main limiting factor for now. Our implementations achieved up to 86% of the effective maximal throughput of global memory. On the 3D lid-driven cavity test case, we obtained 2× to 3× speed-up over previously published implementations. Moreover, we show that, compared to LBGK, the more stable and accurate MRT, despite its higher computational cost, yields equivalent performances on GPUs. Our approach, being simpler than the previous ones, exerts less pressure on hardware. Hence, our method will allow to implement more complex models in the near future.

Figure 1 :

 1 Figure 1: Peak performances GPU vs CPU (source nVidia)

Figure 2 :

 2 Figure 2: CUDA hardware (source nVidia)

Figure 3 :

 3 Figure 3: CUDA programming model (source nVidia)

#

  define id(j, k) k + SIZE*(j) __device__ int index(void) { int x = threadIdx.x; int y = blockIdx.x; int z = blockIdx.y; return x + y*L + z*L2; } __global__ void kernel(int N, float* t) { int k = index(); for (int j = 0; j <= N; j++) { t[id(j+1, k)] = t[id(j, k)]*0.5; } } extern "C" void launch_kernel(int N, float* t) { dim3 grid(L, L); dim3 block(L); kernel<<<grid, block>>>(N, t); cudaThreadSynchronize(); } Code 1: Data transfer benchmark kernel

Figure 4 :

 4 Figure 4: Average transfer time

Figure 5 :

 5 Figure 5: The D3Q19 stencil

Figure 6 :Figure 7 :

 67 Figure 6: Propagation using shared memory

Figure 8 :

 8 Figure 8: Split propagation scheme

Figure 9 :

 9 Figure 9: Reversed propagation scheme

Figure 10 :

 10 Figure 10: Validation for Re = 1,000

  64

Table 1 :

 1 Features of the GT200

	Number of SMs	30
	Number of SPs per SM	8
	Registers per SM	16,384
	Shared Memory	16 KB
	Constant Cache	8 KB
	Texture Cache	8 KB
	Global Memory	896 MB or 1 GB

• No denormalized numbers. Numbers with null exponent are considered as zero.

Table 2 :

 2 displays the time needed for a warp to perform the most common single precision floating point operations : Cost of floating point operations

	Operation	Cycles
	Add, multiply, multiply-add	4
	Reciprocal, logarithm	16
	Sine, cosine, exponential	32
	Divide	36

Table 3 :

 3 assembles the obtained results for both the LBGK and the MRT kernels (rcp stands for reciprocal, mad for multiply-add). Algorithmic complexity of LBGK and MRT kernels

		add	sub	mul	div	rcp	mad	cycles
	LBGK	63	30	48	0	1	34	716
	MRT	76	80	51	0	0	18	900

Table 4 :

 4 Performances for LBGK

		3	96 3	128 3	160 3
	Performance (MLUPS)	471	512	481	482
	Ratio to max. throughput	79%	86%	81%	81%
	Occupancy rate	31%	19%	25%	16%
		64 3	96 3	128 3	160 3
	Performance (MLUPS)	484	513	516	503
	Ratio to max. throughput	81%	86%	86%	84%
	Occupancy rate	25%	19%	25%	16%

Table 5 :

 5 Performances for MRT


Nadia Maïzi

Mathilde Drouineau
email: m.drouineau@mines-paristech.fr

Edi Assoumou

Vincent Mazauric

Flexibility and reliability in long-term planning exercises dedicated to the electricity sector

Keywords: Long-term energy planning, Flexibility, Reliability of supply -

Long-term planning models are useful to build plausible options for future energy systems and must consequently address the technological feasibility and associated cost of these options. This paper focuses on the electricity sector and on problems of flexibility and reliability in power systems in order to improve results provided by long-term planning exercises: flexibility needs are integrated as an additional criterion for new investment decisions and, reliability requirements are assessed through the level of electrical losses they induced and a related cost. These approaches are implemented in a long-term planning model and demonstrated through a study of the Reunion Island.

Introduction

Electricity is a very convenient way to deliver huge amounts of power in areas where demand is concentrated. Due to the predicted population densification, electricity consumption is set to significantly increase over the next decades. Besides, the International Energy Agency estimates that during the next thirty years US$10 trillion will be spent on generation, transmission and distribution of electricity [START_REF]Organisation for Economic COopération and Developement[END_REF], in order to replace existing capacities in developed countries, to accompany the development of energy markets in developing countries or to substitute energy vectors that are less clean or growing scarce.

In this context, two critical features challenge forthcoming changes in power systems:

• the emergence of different paradigms for serving electricity than those for which the system was designed [START_REF] Ilic | From hierarchical to open access electric power systems[END_REF]. For instance, high shares of renewable energy sources may become a critical aspect of future energy systems, both for centralized scheme, with the Desertec concept [START_REF] Desertec | [END_REF], and for distributed architecture, with the smartgrid concepts [START_REF] Siddiqui | The green grid: Energy savings and carbon emissions reductions enabled by a smart grid[END_REF].

• the will to improve the energy efficiency of electrical power transmission, given that the current system clearly lacks efficiency. For instance, since 1950 in the US, electrical losses have been twice as high as electricity consumption as pointed out in [START_REF]Annual Energy Review 2007[END_REF] and illustrated in figure 1 for the industrial sector (same trends for the commercial and residential sectors). Moreover, the world net generation of electricity increasingly relies on fossil fuels as shown in figure 2, which implies a higher impact on the environment due to the level of greenhouse gas emitted by the electrical sector. [START_REF]Key world energy statistics from the IEA[END_REF]. Other includes geothermal, solar, wind, combustible renewables and waste, and heat.

These issues stress the need for assessing future electrical power systems and the MARKAL (MARKet ALlocation) type of technology-rich models provide a partial solution to this problem [START_REF] Fishbone | MARKAL, a linear-programming model for energy systems analysis: technical description of the BNL version[END_REF][START_REF] Etsap | Energy Technology Systems Analysis Programme[END_REF][START_REF] Assoumou | La famille MARKAL de modèles de planification énergétique: un complement aux exercices de modélisation dans le contexte français[END_REF]. Actually, models for energy planning -like the MARKAL family of models -have proven useful to determine plausible evolutions of the energy sector in the mid-to long-term when facing strong environmental pressures, such as carbon mitigation or fossil energy depletion. Indeed, this kind of models offers substitution possibilities throughout the whole energy chain and subsequently enables to address different issues such as pointing out the main drivers of the energy system at a given regional scale, anticipating changes in and impacts of energy prices, or estimating pollutant emissions.

However, the feasibility of the results provided by energy planning models -namely their sustainability and robustness -depends on the description of spatiality and time constraints of power systems:

• Spatiality constraints are related to the geographical distribution of power plants, the structure and availability of the transmission network, and the location of the consumers.

• Time constraints are related to the electric current and the management of the network in order to minimize the variations in tension and frequency at any time, while keeping the system within safety limits. After recovering a stationary behavior, the total electric power supplied by the stations must be equal to the network demand. In order to ensure this power balance in real time, electricity production modes will be chosen for their dynamic nature and for their location in the network.

These spatiality and time constraints come from the needs for flexibility and reliability when operating power systems. In this paper, we choose to improve the description of flexibility and reliability requirements in long-term energy planning models in order to provide more plausible results for the electricity sector. We quantitatively address these issues in section 1.2 and 2.2, respectively for the flexibility and reliability needs.

We start with a preliminary presentation of the current MARKAL modeling approach for the electricity sector. The paper is then divided into three contributions:

1. We present an augmented MARKAL model that introduces flexibility as an additional criterion for electricity generation investments: the choice of power stations relying on a cost minimization is then tempered by specific technical constraints (section 1).

2. We describe reliability requirements on power systems and we present a methodology that provides the cost of reliability of supply in future power systems. This analysis is based on a thermodynamic approach, which enables a synthetic description of power transmission (section 2).

3. Finally, we study the case of the Reunion Island for which we conduct a longterm analysis of the electricity sector, and for which we aim at implementing the approaches presented in sections 1 and 2. The Reunion Island is interesting because it relies on a small and isolated power system, which exarcerbates the needs for flexibility and reliability; and because they aimed to have in 2030 an energy consumption based to 100% on renewable energy sources [START_REF] Arer | PETREL -Ile de la Réunion Plan Economique de Transition et de Relance via des Energies 100% Locales à l'île de la Réunion[END_REF] (section 3).

The last section gives some concluding remarks.

Preliminary presentation of the MARKAL models

To start with, we present the main features of the MARKAL long-term planning models and underline some strengths and weaknesses concerning the representation of the electricity sector.

Long-term global prospective models permit the assessment of multi-sectoral energy policies and are therefore persuasive. Among these models, the MARKAL family of models are technological models developed since the mid-eighties [START_REF] Fishbone | MARKAL, a linear-programming model for energy systems analysis: technical description of the BNL version[END_REF] under the auspices of the International Energy Agency [START_REF] Etsap | Energy Technology Systems Analysis Programme[END_REF]. MARKAL, in its basic version, is a technically optimum model. It relies on an explicit formulation of the input/output relationships for each technology and minimizes -over the chosen time horizon and for a given final outcome -the actualized global cost. The optimization is subject to constraints such as energy management features, caps for CO 2eq emissions, limitations on fuel shares in electricity generation, etc. The decision variables depend on the choice of the activity level of technologies, and of capacity investments. The equilibria of energy flows are generally expressed over the year and evaluated on total energy rather than on hourly power demand. A synthetic description of the input/output relationships is given in the scheme covering the whole energy chain depicted in the figure 3, usually called the Reference Energy System. MARKAL offers a more detailed description for the electricity sector. Specific technical constraints [START_REF]System for the analysis of global energy market, model documentation[END_REF] are represented in the model:

• Flow equilibrium constraints:

Electricity and heat are represented in more detail in the model. The time divisions applied to these two energy vectors are shorter and each period is broken down into six sub-periods showing the combinations between, on the one hand, three seasons (summer, winter, intermediate), and on the other hand, day and night. The flux equilibrium equations are then published separately for each of these sub-periods.

• Peak reserve capacity constraints: A peak reserve constraint guarantees the setting-up of a supplementary capacity reserve to cope with high demand periods. A peak equation stipulates that the total production capacity must be oversized by a certain percentage1 to satisfy the demand (for exports, processes and demand technologies) and to insure against several contingencies. It forces to increase the production capacity by the chosen level of reserve. The user specifies two parameters: a global electricity or heat reserve factor, and the contribution2 of each electricity or heat production technology to the reserve factor.

1 Flexibility desciption in MARKAL

Description of the need for flexibility

The need for flexibility is related to the structure of the load curve: installations that run for a relatively long time throughout the year must be installed and actually used. Therefore, to satisfy the electricity demand, all functioning electricity systems require, a priori, the installation and effective use of power stations devoted to the baseload supply, as well as others devoted to peakload supply (few hours per year). There is a strong link between installation decisions motivated by power dimensioning, and the effective use of these installations for energy supply. Thus, the need for flexibility in time management of electricity systems corresponds to a challenging technical constraint, making it necessary to use more expensive technologies that can rapidly adapt their production level in order to meet the consumption level. These flexibility constraints are on the whole largely ignored in MARKAL models, whereas in reality requirements for flexibility call for long-term investments. Indeed, least expensive ways of production are favored, as the MARKAL approach relies on the minimization of economic criteria [START_REF] Berger | Canadian markal: An advanced linear programming system for energy and environmental modelling[END_REF] -which aims at choosing the electricity system that is cost-efficient. Ignoring flexibility presents several drawbacks for the MARKAL results in the electricity sector. They are listed below:

1. It leads to an over-estimation of the importance of power stations that have low production costs but may also have low reactivity on the load curve: this impacts the capacity needs, the type of technologies chosen, and their level of use.

2. The system's running costs can be miscalculated since, even with identical technologies, the constraints of operating at partial load increase the costs of meeting with demand (slower amortization, the need to use a greater number of power stations).

3. The consequences of the increased share of intermittent renewable sources in the power generation mix are not fully handled. In France, the theoretical potential is estimated at 66 TWh for 30 GW of onshore wind power and 97 TWh for 30 GW of offshore wind power [START_REF]Programmation pluriannuelle des investissements de production électrique[END_REF]. As development of these power plants increases, there will be a greater need for more flexibility. The feasibility of scenarios figuring a low share of thermal production and a large share of wind or solar power is questionable. On top of that, intermittent sources call for thermal production capacities. As these plants are the only direct sources of greenhouse gas emissions in the electricity sector, they counterbalance the expected potential for reducing emissions attributed to renewable energy sources. This should be revised in order to better assess greenhouse gas emissions.

For these reasons, flexibility has to be integrated as a constraint within the model's framework.

An augmented flexible MARKAL model

This approach quantifies the average production needs for several specific operating modes, such as semi-base load or peak, and then introduces flexibility as one of the model's selection criteria. It was applied to the French electricity sector in [START_REF] Maïzi | Key features of the electricity production sector through long-term planning : the french case[END_REF], which is dominated by nuclear power. Technological choices are differentiated for each mode, and the flexibility criterion, defined globally, does not depend on short time divisions for the model's equations. Its quantitative assessment, relies on a time analysis that is both shorter and easier to put into place outside the model. The pre-existent mechanisms specific to the electricity sector are conserved.

Results for the standard and for the augmented MARKAL model -taking flexibility into account -are only given for fossil plants regarding that these power plants are sensitive to dynamical features (figures 4 and 5). The figures compare fossil capacities effectively used for electricity generation.

With the standard MARKAL model, gas, fuel or oil power plants are installed to satisfy the peak equation, but the figure 4 shows that they are not used for electricity production. Actually, power stations are used over long periods of time leading preferentially to baseload profiles of power plants running with cheaper fuels.

With the augmented approach, the figure 5 shows on the contrary that oil power plants appear in the electricity generation. Indeed, even if they are more expensive, they are essential plants to meet the peak demand. It confirms that the optimization process with the augmented MARKAL model provides results more plausible concerning the electricity sector. With this model, we also observe that a greater amount of energy is produced from fossil fuels (200 PJ in 2030 compared to 150 PJ), which corresponds to a best representation of the need for flexible technologies.

When flexibility is not properly handled, the results issued from MARKAL reflect a minimal (under-constrained) condition to satisfy global electricity demand. Thus it fails to address properly the questions related to environmental issues, and namely the level of carbon emissions -except for regions where fossil fuels dominate the baseload mix. In France, the level of greenhouse gas emissions depends on the precision of the thermal production assessment [START_REF] Maïzi | Key features of the electricity production sector through long-term planning : the french case[END_REF].

Regarding these results, flexibility appears to be a key parameter for reaching a better representation of the electricity generation system and consequently achieving sustainable power systems. It has also appeared that providing power systems with enough installed capacities -both satisfying the electric demand and following the load curve -is difficult when ignoring the need for flexibility. Flexibility is a step towards planning long-term generation shares where installed capacities are well-anticipated.

The need for reliability of supply

At the side of the need for flexibility also stands the question of the reliability of electricity supply. Reliability of supply relies on the ability of the electric system to withstand sudden disturbances [START_REF]European Network of Transmission Systems Operators for Electricity[END_REF] and is ensured by technical properties of the production means.

Currently, power systems clearly lack of efficiency because this sector is severely disadvantaged by the efficiency of the Carnot cycle of the thermal units, which are about one order of magnitude greater than electrical losses at the transmission, distribution and consumption levels (figure 6). In this context, renewable and distributed energy sources seem to be attractive alternatives for cleaner and unlimited power generation. Furthermore, a wide integration of renewable energy sources is expected to improve the overall efficiency of power systems. Renewable energy sources decrease electrical losses at the production level, whereas distributed energy sources decrease electrical losses at the transmission and distribution levels. However, such a gain in efficiency may be counterbalanced by a decrease in the reliability of electricity supply. The reliability of electricity supply is defined as the capacity of a power system to handle safely sudden disturbances, such as production or load fluctuations. To maintain an appropriate reliability of supply, power systems must provide amounts of energy dedicated to the reliability of supply, which induces additional investments and extra-losses. Consequently, the need for reliability of supply may decrease the expected gain in efficiency with renewable and distributed energy sources. In this section, we introduce the reliability of electricity supply within an optimal discussion and aim at assessing the cost of reliability of supply.

Description of reliability requirements

Firstly, in order to provide relevant comparisons between future electrical systems, we must define accurately the three main kind of losses on power systems:

Losses induced by the Carnot cycles The electricity industry is severely disadvantaged by the efficiency of the Carnot cycle of the thermal units, which are about one order of magnitude greater than electrical losses at the transmission, distribution and consumption levels (figure 6).

Conveyance losses occur during power transmission through the network. They mainly depend on whether or not the transmission grid is congested, on the voltage level, or on the network architecture. They can be assessed from the duration of peak, semi-base or base loads. When production capacities are centralized, transmission takes place through longer distances, and conveyance losses may increase, despite high voltage lines. In fact, for a given geographical distribution of loads and generators, the more the meshing of the grid increases, the more the Joule losses decrease, the voltage profile improves and the system is more stable.

Reliability-induced losses are linked to the needed level of reliability. The level of losses corresponds to the stocks of electromagnetic coupling energy and of kinetic and spinning energy required to handle load fluctuations. A part of the electricity production must be dedicated to maintain the appropriate stocks of electromagnetic, kinetic and spinning energy. Consequently, the electricity industry must be (over-)sized to produce enough electricity for the consumers and also for providing the stocks. The part of the electricity dedicated to the stocks induces investments in dedicated capacities and additional electrical losses.

With a wide integration of renewable energy sources on power systems, the losses induced by the Carnot cycles are expected to decrease, as energy production will rely on fewer thermal units. With a wide development of distributed energy sources, the conveyance losses are expected to have similar benefits, as the average distance between production and consumption decrease. Conversely, reliability-induced losses are expected to increase, as renewable and distributed energy sources do not maintain sufficient stocks of electromagnetic coupling energy and of kinetic and spinning reserves.

Secondly, we describe more precisely reliability of electricity supply. According to the UCTE handbook [START_REF]European Network of Transmission Systems Operators for Electricity[END_REF], the reliability of an electric system is addressed by considering both the adequacy of the system, i.e. the ability to supply the aggregate electrical demand and energy requirements of the customers at all times, and the security of the system, i.e. the ability to withstand sudden disturbances such as unanticipated loss of system elements.

To ensure these two elements, power systems rely on the dynamic management of frequency and voltage [START_REF] Bergen | Power System Analysis[END_REF]. Frequency and voltage are crucial quantities in power systems, whose great deviations can lead to brownouts or power outages when the system recovers from production or load fluctuations or when it experiences transient states (e.g. lightnings). Basically, maintaining appropriate variations of frequency depends on the levels of kinetic and spinning reserves on the system, and maintaining appropriate variations of voltage depends on the electromagnetic energy on the system. Levels of kinetic and spinning reserves are ensured by the mechanical inertia of the rotating generators on the system and can be adjusted by investing in weighing generation machines or flywheel. The electromagnetic coupling energy is the energy involved in the electromechanical and electromagnetic power conversions. It is absolutely necessary to enable power transmission and it is stocked in electromagnetic materials, namely rotating power machines and transformers. The level of electromagnetic coupling energy can be adjusted by investing preferentially in rotating power machines and transformers.

Furthermore, with high shares of renewable energy sources on power system, levels of electromagnetic energy and of kinetic and spinning reserves decrease, whereas intermittent energy sources induce abundant production fluctuations, with higher magnitude. High shares of renewable energy sources in electricity mix may decrease the reliability of electricity supply.

Assessing the cost of reliability

In the following, we intend to compare the level of network losses for power systems with different production shares of renewable energy sources -in particular intermittent ones. To do so, we apply variational principles deduced from thermodynamics to achieve a global description of power systems and of their subsequent level of losses [START_REF] Mazauric | From thermostatistics to Maxwell's equations: a variational approach of electromagnetism[END_REF][START_REF] Mazauric | A global approach of electromagnetism dedicated to further long-term planning[END_REF].

The system on which we apply the variational principles is the electromagnetic field, conveying electricity from the generation units to the consumers:

1. In the thermodynamic framework, the deviation between the mechanical power flowing through the network and the variation with time of the Helmoltz freeenergy F gives an evaluation of the Joule losses P Joule , at the transmission, distribution and consumption levels. P Joule gives an evaluation of the irreversibility experienced by the system. According to the thermodynamics, Joule losses are always positive and the lower they are, the more reversible is the evolution of the system. The latter statement can be expressed with an optimal condition, which means that the system always tends to minimize the Joule losses during its evolution. Ideally, a centralized power system can be described from the Helmoltz free-energy, and this reversibility condition provides the very minimum amount of Joule losses during the power transaction and matches the conveyance losses.

Added to the losses induced by the Carnot cycles, it gives the total amount of primary energy lowered in a centralized system.

2. Conversely, in order to explicitly take reliability into account, the description should also include the inertial behavior of the electromagnetic coupling. Hence, the actual level of Joule losses should be assessed from another reversibility condition obtained with the Gibbs free-energy G = F -ΦI exc , where ΦI exc is the electromagnetic coupling and I exc the excitation provided to the generator. As a physical result, the latter description provides the level of both conveyance and reliability losses. Compared to the case (1), the argument of the optimum is drastically changed. This leads to new stability requirements for dynamic management.

Added to the losses induced by the residual Carnot cycles, it gives the total amount of primary energy lowered in a distributed system.

This analysis provides the level of Joule losses P Joule for operating electric power systems in a reliable way. It will enable to compare the overall amounts of losses for a centralized power system and for a distributed one, when the electrical demand and the level of reliability of supply are the same. The cost of electricity production in centralized systems is severely disadvantaged by the efficiency of the Carnot cycle of the thermal units, whereas in distributed system it is more penalized by the cost of reliability of electricity supply. We now focus on a TIMES model dedicated to the supply and power sector of the Reunion Island (figure 7), where TIMES is the latest evolution of the MARKAL family of models. The Reunion Island is blessed with high potentials of renewable energy sources, which may decrease the reliability of electricity supply, its power system is small, weakly-meshed and isolated, and its regional government has set the binding target of an energy consumption based to 100% on renewable energy sources in 2030 [START_REF] Arer | PETREL -Ile de la Réunion Plan Economique de Transition et de Relance via des Energies 100% Locales à l'île de la Réunion[END_REF]. The current use of renewable energy sources in the electricity sector is 36% [START_REF] Edf | Bilan Prévisionnel Pluriannuel : investissements en production (La Réunion)[END_REF], so the energy system has to change substantially for reaching the target. In this section we present the model development and the results obtained with the TIMES-Reunion model. The case study of the Reunion Island finally poses the question of the technological feasibility of the results proposed by the model. The technological feasibility may be studied in the light of the flexibility and reliability requirements presented previously.

The case of the Reunion Island

Model development

Resources supplies

The resources in the Reunion Island are domestic sugarcane bagasse and importations of coal, and fuel oils (heavy and distillate). The fossil energy import prices are based on the projections of the World Energy Outlook [21]. We assume that fuel oil prices follow the projections of crude oil prices. Fossil-fuel price assumptions are listed in table 1. The Reunion Island also produces around 10% of its annual electricity consumption with the combustion of the sugarcane bagasse. The cost of the bagasse sugarcane is set to zero as the bagasse is a co-product of the sugar factories and that these factories are on the same production areas than the thermal power plants using sugarcane bagasse. Electricity production from bagasse takes place in the power plants of Le Gol (111.5 MW) and Bois-Rouge (100 MW), which also work with coal, thus producing electricity apart from the season of sugar production.

Electricity demand

In 2008, the general features of the electricity sector were the followings [START_REF] Edf | Bilan Prévisionnel Pluriannuel : investissements en production (La Réunion)[END_REF][START_REF]Bilan énergétique 2009 de l'île de la Réunion[END_REF]: Electricity consumption rose up to 2546 GWh, divided into 50% coal, 14% other fossil fuels, 25% hydroelectricity, 10% sugarcane bagasse and 1% others; the electricity peak demand was 408 MW; and the total installed capacities were slightly less than 650 MW.

Since 1995, the growth rate of electricity demand has decreased from 6.7% to 2.8%, and it is expected to continue decreasing and reach a value between 1 and 2% in 2025 [START_REF] Edf | Bilan Prévisionnel Pluriannuel : investissements en production (La Réunion)[END_REF]. A projection for electricity consumption growth until 2025 was provided by Électricité de France (EDF) [START_REF] Edf | Bilan Prévisionnel Pluriannuel : investissements en production (La Réunion)[END_REF] (see table 2).

Existing and coming power plants

Data for existing capacity, capacity factors, and efficiency are derived from reports on existing power plants by EDF and the Regional Agency for Energy in the Reunion Island (ARER) [START_REF] Edf | Bilan Prévisionnel Pluriannuel : investissements en production (La Réunion)[END_REF][START_REF]Bilan énergétique 2009 de l'île de la Réunion[END_REF]. Following discussion with experts, some of the technico-economical data have been revised to correspond more accurately to the electricity mix, in particular with the current spread of renewable energy sources. Costs and performance characteristics of new power plants are derived from the database of the European RES20203 project.

Renewable potentials

According to experts and to the literature, the available renewable sources for the Reunion Island and their technical and economic resource potential are those presented in the table 3.

Scenarios specification

Various scenarios are investigated wherein the electricity system is required to simultaneously meet a given level of electricity demand (table 2) and reach the 2030 target with an electricity consumption based to 100% on renewable energy sources. Height scenarios are built around three main assuptions concerning levels of fossil fuels imports, electricity demand, and sugarcane bagasse potential. They are summarized in the table 4:

• An upper limit for the fossil fuels imports is set in 2008 and linearly decreased to 0 in 2030. The objective is to study the evolution of the electricity production without fossil fuels. We propose an alternative scenario with a softened constraint, where the limit is only set on coal imports. These scenarios are compared with a business-as-usual with no limits on importations.

• We consider scenarios with lower electricity consumption to study how lowering the demand helps reaching the 100% target.

• For all these scenarios, we set the potential for the renewable energy sources at the maximum rates described previously, except for sugarcane bagasse. We finally consider the option with a higher available potential, where the sugarcane industry is only dedicated to energy production. The scenario MedDEM corresponds to the business-as-usual scenario.

Results

Electricity shares of 2008

In table 5, we compare the 2008 electricity production given by the model with actual data [START_REF] Edf | Bilan Prévisionnel Pluriannuel : investissements en production (La Réunion)[END_REF] and the similarity of the results is encouraging. In particular, the proportion between production based on fossil fuels and production based on renewable sources are the same. The main difference concerns the production based on fossil fuels. With the model, coal participates to a higher share of electricity production. This is explained by the fact that in this TIMES-Reunion model, load curves and peak demands are not fully described. Thus, the modeled electricity system requires in a lesser extent power plants dedicated to peak loads such as fuel oil turbines.

Towards an electricity consumption based to 100% on renewable energy sources

We then present the evolution of the electricity mix when the limitation of coal or all fossil fuel imports is set on the standard level of electricity demand. In the business as usual scenario (figure 8(a)), the actual shares of production roughly follow the increase in the demand. However, there is no development of renewable energy sources despite their potentials. This is explained by the cost of power plants relying on these sources, compared to those relying on coal. The model may be improved with a precise description of the incentives' system favoring renewable energy sources or with another constraint on the greenhouse gas emissions in the Reunion Island. Furthermore, in this scenario, power plants based on fuel oils are still hardly used, which corroborates the need to describe the time characteristics of the electricity demand.

The two other scenarios (figures 8(b) and 8(c)) show that there is a large room for the development of renewable energy sources when a constraint is set on coal and all fossil fuels import. In both scenarios, new renewable energy sources appear such as geothermy, ocean thermal energy, or wave energy. Hydroelectricity increases sharply. Nevertheless, with the constraint on both coal and petroleum products, the figure 8(c) shows that around 500 GWh of electricity demand can not be met in 2030.

Lower growth of electricity demand

The decrease in the electricity demand may be due to two main features: a slower growth rate of the demand, and a deep improvements of energy efficiency of end-use devices. If electricity consumption goes in this direction, it will consequently ease the Reunion Island to reach its ambitious target of 100% of renewable energy sources by 2030. The new electricity consumption growth is taken from a lower scenario also provided by EDF [START_REF] Edf | Bilan Prévisionnel Pluriannuel : investissements en production (La Réunion)[END_REF]. The results with a lower electricity demand are depicted in the figure 9, which presents the changes induced on the electricity mix.

Higher potential for sugarcane bagasse production

Finally, we consider a radical change of the sugarcane industry, where sugarcane production is only dedicated to electricity production, as suggested in [START_REF] Arer | PETREL -Ile de la Réunion Plan Economique de Transition et de Relance via des Energies 100% Locales à l'île de la Réunion[END_REF]. In this case, the available potential of sugarcane bagasse is increased as the whole sugarcane is used for energy purposes. Interestingly, it seems that the whole available potentials for renewable energy sources can meet the lower electricity demand (figure 10).

Concluding remarks

In this paper, we have studied spatiality and time constraints for power generation system with the flexibility and reliability requirements. This work allows the following understandings:

• On the one hand, flexible features of electricity are handled using an augmented version of the MARKAL model. Regarding its results, flexibility appears to be a key parameter for reaching a better representation of the electricity generation system and consequently achieving sustainable systems. Taking flexibility into account has made it possible to provide electric systems with enough installed capacities, both satisfying the electric demand and following the load curve. It is a step towards long-term planning models where electrical installed capacities are well-anticipated.

• On the other hand, for reliability requirements, we have proposed a method that provides the cost of reliability of supply. This method is based on a thermodynamic approach applied to power systems. It emphasizes the need to take reliability requirements into account for the design of future power systems. As a result, the total cost of future energy systems may increase.

• Finally, we propose to implement the approaches dealing with the flexibility and reliability needs of power systems in a TIMES model dedicated to the Reunion Island. This case study is a good illustration of the previous results: its small island grid without interconnection and the proposed shares for renewable energy sources raise the two challenging issues -flexibility and reliability requirements -for future power systems.

Figure 1 :

 1 Figure 1: Industrial Total Energy Consumption, Major Sources [5]: electricity use (and related losses) expanded dramatically.

Figure 2 :

 2 Figure 2: Electricity generation in TWh by fuel since 1971 [6]. Other includes geothermal, solar, wind, combustible renewables and waste, and heat.

Figure 3 :

 3 Figure 3: Synthetic view of the Reference Energy System, issued from [11].

Figure 4 :

 4 Figure 4: Fossil plant production 2000-2030: a standard MARKAL approach.

Figure 5 :

 5 Figure 5: Fossil plant production 2000-2030: an augmented MARKAL approach.

Figure 6 :

 6 Figure 6: Energy supply-chain in Mtoe in 2007 (e!Sankey diagram). The electricity sector is severely disadvantaged by the efficiency of the Carnot cycle of the thermal units, and by the electrical losses at the transmission, distribution and consumption levels [6].

Figure 7 :

 7 Figure 7: Electricity generation on the Reunion Island [19].

  (a) Business as Usual. (b) Without importation of coal in 2030. (c) Without importation of fossil fuels in 2030.

Figure 8 :

 8 Figure 8: Shares of the electricity production for the electricity demand proposed in the table 2. The importations of coal and petroleum products are constrained over the time horizon: from their levels in 2008 to 0 in 2030.

Figure 9 :

 9 Figure 9: Shares of the electricity production with a lower electricity demand.

Figure 10 :

 10 Figure 10: Shares of the electricity production with a higher sugarcane bagasse potenatial and a lower electricity demand.

Table 1 :

 1 Fossil-fuel price assumptions [21].

	Real terms (2008 prices)	Unit	2000	2008	2015	2020	2025	2030
	OECD steam coal imports	$2008/ton	41.22 120.59 91.05 104.16 107.12	109.4
	IEA crude oil imports	$2008/barrel	34.3	97.19	86.67 100.00 107.50 115.00
	Heavy fuel oil	e2008/ton	-	196	174	201	216	231
	Distillate fuel oil	e2008/hl	-	47	42	48	51	55

Table 2 :

 2 Electricity consumption growth in the Reunion Island from the medium scenario of EDF[START_REF] Edf | Bilan Prévisionnel Pluriannuel : investissements en production (La Réunion)[END_REF]. This scenario is extended from 2025 to 2030 with a growth rate of 1.5%.

		Unit	2008	2010	2015	2020	2025	2030
	Electricity consumption GWh 2 546 2 710 3 110 3 500 3 805 4100
	Growth rate	%	3.4	3.2	2.6	2.4	1.5	1.5
	Power	MW	408	445	520	595	670	720
	Energy sources		Current levels		Potentials
	Biomass		260 GWh		400 GWh
	Hydropower	121 MW (553 GWh) 177 MW until 2012,
					then 268 MW
	Wind		16,8 MW		50 MW
	Solar PV		10 MW			160 MW
	Ocean Thermal		-		10 MW in 2020,
	Energy Conversion				100 MW in 2030
	Wave Energy		-		30 MW (by 2014)
	Geothermy		-			30 MW
	Storage Capacities		-		1 MW in 2009,
						then 10 MW

Table 3 :

 3 Renewable energy potentials in the Reunion Island.

Table 4 :

 4 Scenarios specification. The scenarios are built around three main assuptions concerning electricity demand, fossil fuel imports and sugarcane bagasse potential.

Table 5 :

 5 The 2008 electricity shares given by the TIMES-Reunion model compared to the actual values of EDF[START_REF] Edf | Bilan Prévisionnel Pluriannuel : investissements en production (La Réunion)[END_REF].

	Energy sources	Model (%) EDF (%)
	Coal	56.90	50.55
	Fuel Oils (Distillate and Heavy)	9.06	13.30
	Sugarcane bagasse	10.21	10.31
	Hydroelectricity	21.71	24.86
	Wind energy	1.19	0.53
	Solar energy	0.41	0.42
	Municipal waste	0.52	0.03
	Production	2 547 GWh 2 546 GWh

For instance, in Metropolitan France, the value of the peaking factor can be set to 1.6.

The coefficient (from 0 to 1) specifies the fraction of the technology's capacity that is allowed to contribute to the peak load and makes it possible to differentiate between the contributions of different power plants. Nuclear power plant has a peak coefficient of 1, whereas wind farm has a peak coefficient of 0.2 or 0.3.

RES2020 is a European project which aims at monitoring and assessing the implementation of the directives on Renewable Energy Sources and the policy recommendations for 2020 in the EU-27. With this project, a number of future options for policies and measures are defined and studied with the use of TIMES.

Acknowledgment

This research was supported by the Chair Modeling for sustainable development, driven by MINES ParisTech, Ecole des Ponts ParisTech, AgroParisTech, and ParisTech; and supported by ADEME, EDF, RENAULT, SCHNEIDER ELECTRIC and TOTAL.