I Lacroix-Violet

C Chainais-Hillairet

On the existence of solutions for a drift-diffusion system arising in corrosion modelling

In this paper, we consider a drift-diffusion system describing the corrosion of an iron based alloy in a nuclear waste repository. In comparison with the classical driftdiffusion system arising in the modeling of semiconductor devices, the originality of the corrosion model lies in the boundary conditions which are of Robin type and induce an additional coupling between the equations. We prove the existence of a weak solution by passing to the limit on a sequence of approximate solutions given by a semi-discretization in time.

Introduction

General framework of the study

In this paper, we consider a system of partial differential equations describing the corrosion of an iron based alloy in a nuclear waste repository. At the request of the French nuclear waste management agency ANDRA, investigations are conducted to evaluate the longterm safety assessment of the geological repository of high-level radioactive waste. The concept for the storage under study is the following: the waste shall be confined in a glass matrix, placed into cylindrical steel canisters and stored in a claystone layer at a depth of several hundred of meters. The long-term safety assessment of the geological repository has to take into account the degradation of the carbon steel used for the waste overpacks, which is mainly caused by generalized corrosion processes.

In this framework, the Diffusion Poisson Coupled Model (DPCM) has been developed by Bataillon et al [START_REF] Bataillon | Corrosion modelling of iron based alloy in nuclear waste repository[END_REF] in order to describe the corrosion processes at the surface at the steel canisters. It assumes that the metal is covered by a dense oxide layer which is in contact with the claystones. The model describes the development of the dense oxide layer. In most cases, the thickness of the oxide layer ranges from nanometers to micrometers and is much smaller than the size of the steel canister, which justifies a 1D-modelling. The model is made of electromigration-diffusion equations for the transport of charge carriers (three species are considered: electrons, Fe 3+ cations and oxygen vacancies), coupled with a Poisson equation for the electrostatic potential. The interaction between the oxide layer and the adjacent layers (metal and claystone) are described by Butler-Volmer laws for the electrochemical surface reactions and potential drops. The model includes moving boundary equations.

Numerical methods for the approximation of the DPCM model have been designed and studied by Bataillon et al in [START_REF] Bataillon | Numerical methods for the simulation of a corrosion model with moving oxide layer[END_REF]. Numerical experiments with real-life data shows the ability of model to reproduce expected physical behaviors. However, theoretical questions like existence of solutions and long-time behavior of the DPCM model have not been investigated yet. In this paper, we consider a simplified version of the DPCM model, already introduced in [START_REF] Bataillon | Numerical methods for the simulation of a corrosion model with moving oxide layer[END_REF]. In this case, the interfaces of the oxide layer are assumed to be fixed and only two charge carriers are taken into account: electrons and cations Fe 3+ . We will prove the existence of a global in time solution for this simplified model.

Presentation of the model

We recall here the simplified version of the DPCM model introduced in [START_REF] Bataillon | Numerical methods for the simulation of a corrosion model with moving oxide layer[END_REF]. It is a dimensionless model. The unknowns are the densities of electrons and cations, which are respectively denoted by N and P and the electrical potential denoted by Ψ. The current densities are denoted by J N and J P ; they contain both a drift part and a diffusion part. As we do not take into account the moving boundary equations in the simplified model, the domain describing the oxyde layer is the interval (0, 1). The equations describing the electrochemical behavior of the dense oxyde layer write, for t ≥ 0:

∂ t P + ∂ x J P = 0, J P = -∂ x P -3P ∂ x Ψ, in (0, 1), ( 1a
)
ε∂ t N + ∂ x J N = 0, J N = -∂ x N + N ∂ x Ψ, in (0, 1), ( 1b
) -λ 2 ∂ xx Ψ = 3P -N + ρ hl , in (0, 1), (1c)
where ρ hl is the net charge density of the ionic species in the host lattice (we will assume that it is a constant in all the sequel) and λ 2 and ε are positive dimensionless parameters arising in the scaling. Let us just mention that ε is the ratio of the mobility coefficients of cations and electrons, then ε ≪ 1.

As the equations for the carrier densities (1a) and (1b) have the same form, we will use the synthetical form :

ε u ∂ t u + ∂ x J u = 0, J u = -∂ x u -z u u∂ x Ψ, in (0, 1). ( 2
)
For u = P, N , the charge numbers of the carriers are respectively z u = 3, -1 and we respectively have ε u = 1, ε.

Let us now focus on the boundary conditions. Charge carriers are created and consumed at both interfaces x = 0 and x = 1. The kinetics of the electrochemical reactions at the interfaces are described by Butler-Volmer laws at interfaces. It leads to Robin boundary conditions on N and P . As in [START_REF] Bataillon | Numerical methods for the simulation of a corrosion model with moving oxide layer[END_REF], we assume that the boundary conditions for P and N have exactly the same form. Therefore, for u = P, N , they write:

-J u = r 0 u (u, Ψ) on x = 0, (3a)

J u = r 1 u (u, Ψ, V ) on x = 1, ( 3b
)
where V is a given applied potential (we just consider here the potentiostatic case) and r 0 u and r 1 u are linear and monotonically increasing functions with respect to their first argument. More precisely, due to the electrochemical reactions at the interfaces, r 0 u and r 1 u , for u = P, N , write:

r 0 u (s, x) = β 0 u (x)s -γ 0 u (x), (4a)
r 1 u (s, x, V ) = β 1 u (V -x)s -γ 1 u (V -x), ( 4b
)
where the functions (β i u ) i=0,1 , (γ i u ) i=0,1 are expressed with the interface kinetic coefficients (m i u , k i u ) i=0,1 , the positive transfer coefficients (a i u , b i u ) i=0,1 , the maximum occupancy for octahedral cations in the lattice P m and the electron density in the state of metal N m . They write, for u = P, N :

β i u (x) = m i u e -zub i u x + k i u e zua i u x , i = 0, 1, (5a)
γ 0 u (x) = m 0 u u m e -zub 0 u x , γ 1 u (x) = k 1 u u m e zua 1 u x . ( 5b
)
The boundary conditions for the Poisson equation take into account that the metal and the solution can be charged because they are respectively electronic and ionic conductors. Such an accumulation of charges induces a field given by the Gauss law. These accumulations of charges depend on the voltage drop at the interface by the usual Helmholtz law which links the charge to the voltage drop through a capacitance. The parameters ∆Ψ pzc 0 and ∆Ψ pzc 1 are the voltage drop corresponding to no accumulation of charges respectively in the metal and in the solution. Finally, the boundary conditions for the electric potential write:

Ψ -α 0 ∂ x Ψ = ∆Ψ pzc 0 , on x = 0, (6a)

Ψ + α 1 ∂ x Ψ = V -∆Ψ pzc 1 , on x = 1, (6b)
where α 0 and α 1 are positive dimensionless parameters arising in the scaling.

The system is supplemented with initial conditions, given in L ∞ (0, 1):

N (0, x) = N 0 (x); P (0, x) = P 0 (x), x ∈ (0, 1). (7)
In all the sequel, we will assume that the interface kinetic coefficients are given constants and satisfy:

m 0 u , k 0 u , m 1 u , k 1 u > 0, for u = P, N. ( 8
)
The positive transfer coefficients satisfy:

a 0 u , b 0 u , a 1 u , b 1 u ∈ [0, 1], for u = P, N, (9)
and we assume that

3P m -N m + ρ hl = 0. ( 10
)
Indeed, in the applications, we have ρ hl = -5, P m = 2 and N m = 1, so that the relation ( 10) is satisfied. Then, it is expected that the solution to the corrosion model verifies 0 ≤ P (t, x) ≤ P m and 0 ≤ N (t, x) ≤ N m almost everywhere. These estimates will be proved in Theorem 1.1; they are crucial for the proof of existence of solutions to the corrosion model.

Main result and outline of the paper

Let us first note that the system of equations (1a)-( 1c) is a simplified case (linear case) of the so-called drift-diffusion system. This model is currently used in the framework of semiconductors device modelling (see for instance [START_REF] Jüngel | Quasi-hydrodynamic semiconductor equations[END_REF][START_REF] Markowich | Semiconductor equations[END_REF][START_REF] Markowich | The stationary semiconductor device equations[END_REF]) or plasma physics (see [START_REF] Chen | Introduction to plasma physics and controlled fusion[END_REF]).

In this context the drift-diffusion model has been widely studied, from the analytical as from the numerical point of view. Let us refer to [START_REF] Da | On the semiconductor drift-diffusion equations[END_REF][START_REF] Fang | Global solutions of the time-dependent drift-diffusion semiconductor equations[END_REF][START_REF] Gajewski | On existence, uniqueness and asymptotic behaviour of solutions of the basic equations for carrier transport in semiconductors[END_REF][START_REF] Gajewski | On the uniqueness of solutions to the drift-diffusion model of semiconductors devices[END_REF][START_REF] Gasser | The initial time layer problem and the quasineutral limit in a nonlinear drift diffusion model for semiconductors[END_REF][START_REF] Jüngel | On the existence and uniqueness of transient solutions of a degenerate nonlinear drift-diffusion model for semiconductors[END_REF][START_REF] Jüngel | A nonlinear drift-diffusion system with electric convection arising in electrophoretic and semiconductor modeling[END_REF][START_REF] Jüngel | A hierarchy of hydrodynamic models for plasmas. Zeroelectron-mass limits in the drift-diffusion equations[END_REF][START_REF] Jüngel | Rigorous derivation of a hierarchy of macroscopic models for semiconductors and plasmas[END_REF][START_REF] Jüngel | The quasi-neutral limit in the quantum drift-diffusion equations[END_REF][START_REF] Markowich | Semiconductor equations[END_REF][START_REF] Schmeiser | A singular perturbation analysis of reverse biased pn-junctions[END_REF] for results on existence and uniqueness of solutions (using some fixed points theorem) and asymptotic limits of the model (as quasi-neutral limit, zero-relaxation time limit, ...).

From a numerical point of view, different methods have been proposed and studied (see for instance [START_REF] Brezzi | Numerical simulation of semiconductor devices[END_REF][START_REF] Chainais-Hillairet | Finite volume scheme for multidimensional drift-diffusion equations and convergence analysis[END_REF][START_REF] Chainais-Hillairet | Convergence of a finite-volume scheme for the drift-diffusion equations in 1D[END_REF][START_REF] Chen | Analysis of a finite element method for the drift-diffusion semiconductor device equations: the multidimensional case[END_REF][START_REF] Sacco | Mixed finite volume methods for semiconductor device simulation[END_REF][START_REF] Scharfetter | Large signal analysis of a silicon read diode oscillator[END_REF]).

In all these papers, the considered boundary conditions are Dirichlet, Neumann or mixed Dirichlet/Neumann boundary conditions (corresponding to the ohmic contacts and the insulated boundary segments of the semiconductor device). Then, the originality of the corrosion model described in this paper lies in the boundary conditions (3), ( 6) which are of Robin type, and induce an additional coupling between the equations. The authors have shown in [START_REF] Chainais-Hillairet | Existence of solutions for a steady state corrosion of steel model[END_REF] the existence of solution for such a model in the stationary case. Moreover in [START_REF] Bataillon | Numerical methods for the simulation of a corrosion model with moving oxide layer[END_REF], a numerical scheme for the model ( 1)-( 7) is given and analyzed.

Let us now state the main result of this paper: the existence and boundedness of a weak solution to system (1)-( 7).

Theorem 1.1. Assume (8), ( 9), [START_REF] Fang | Global solutions of the time-dependent drift-diffusion semiconductor equations[END_REF] and that N 0 , P 0 ∈ L ∞ (0, 1) satisfy

0 ≤ P 0 (x) ≤ P m , 0 ≤ N 0 (x) ≤ N m , for almost every x ∈ (0, 1). ( 11
)
Let us also assume that ∆Ψ pzc 0 and ∆Ψ pzc 1 verify

- 1 3a 0 P 1 + log a 0 P k 0 P α 0 ≤ ∆Ψ pzc 0 ≤ 1 a 0 N 1 + log a 0 N k 0 N α 0 , (12a) - 1 b 1 N 1 + log b 1 N m 1 N α 1 ≤ ∆Ψ pzc 1 ≤ 1 3b 1 P 1 + log b 1 P m 1 P α 1 . ( 12b
)
Then the problem (1)-( 7) admits a weak solution (P, N, Ψ) which is defined by

P, N ∈ L ∞ ([0, T ] × (0, 1)) ∩ L 2 (0, T ; H 1 (0, 1)), Ψ ∈ L 2 (0, T ; H 1 (0, 1)
) for all T > 0 and:

For u = N, P, ∀ϕ ∈ L 2 (0, T ; H 1 (0, 1)), -ε u T 0 1 0 u∂ t ϕdxdt + ε 1 0 u 0 (x)ϕ(0, x)dx - T 0 1 0 (-∂ x u -z u u∂ x Ψ) ∂ x ϕdxdt + T 0 β 1 u (V -Ψ(t, 1))u(t, 1) -γ 1 u (V -Ψ(t, 1)) ϕ(t, 1) + β 0 u (Ψ(t, 0))u(t, 0) -γ 0 u (Ψ(t, 0)) ϕ(t, 0) dt = 0, ( 13
)
and ∀ϕ ∈ L 2 (0, T ; H 1 (0, 1)),

λ 2 T 0 1 0 ∂ x Ψ∂ x ϕdxdt - T 0 λ 2 α 1 (V -Ψ(t, 1) -∆Ψ pzc 1 ) ϕ(t, 1)dt + T 0 λ 2 α 0 (Ψ(t, 0) -∆Ψ pzc 0 ) ϕ(t, 0)dt = T 0 1 0 (3P -N + ρ hl )ϕdxdt. ( 14
)
Moreover the electron and cation densities N and P satisfy for all t ∈ [0, T ] and x ∈ (0, 1)

0 ≤ P (t, x) ≤ P m , 0 ≤ N (t, x) ≤ N m . ( 15
)
In order to prove this Theorem, we consider a semi-discretization in time for (1)-( 7), as done by A. Jüngel and I. Violet in [START_REF] Jüngel | The quasi-neutral limit in the quantum drift-diffusion equations[END_REF] for the quantum drift-diffusion equations. We first show the existence of a solution to the semi-discrete problem. Then, we obtain some a priori estimates on the solution to the semi-discrete problem. These estimates will ensure the compactness of a sequence of approximate solutions. Finally by passing to the limit on the time step in the semi-discrete scheme, we obtain the existence of a weak solution for the problem ( 1)-( 7) satisfying [START_REF] Jüngel | A nonlinear drift-diffusion system with electric convection arising in electrophoretic and semiconductor modeling[END_REF]. The paper is organized as follows. In Section 2, we present the semi-discrete problem associated to (1)-( 7) and we prove the existence of a solution to this problem. In Section 3, we show some a priori estimates satisfied by the semi-discrete solution. Finally, in Section 4, we establish the convergence of a sequence of approximate solutions and passing to the limit in the semi-discrete scheme, we prove Theorem 1.1.

The semi discrete problem

In this Section, we propose a discretization in time of the system (1)- [START_REF] Chen | Introduction to plasma physics and controlled fusion[END_REF]. Therefore, we introduce a uniform subdivision (t k ) 0≤k≤K of the time interval [0, T ], for some T > 0. The time step is denoted by ∆t and we have t k = k∆t.

We set

P 0 = P 0 , N 0 = N 0 . ( 16
)
Starting from this initial condition, for all 0 ≤ k ≤ K, we consider the following system: [START_REF] Bataillon | Numerical methods for the simulation of a corrosion model with moving oxide layer[END_REF], for u = P, N, (18) supplemented with the following boundary conditions:

-λ 2 ∂ xx Ψ k = 3P k -N k + ρ hl , in (0,1), ( 17
)
ε u ∆t u k+1 -u k + ∂ x J k u = 0, J k u = -∂ x u k+1 -z u u k+1 ∂ x Ψ k , in (0,
Ψ k -α 0 ∂ x Ψ k = ∆Ψ pzc 0 , on x = 0, (19a)
Ψ k + α 1 ∂ x Ψ k = V -∆Ψ pzc 1 , on x = 1, (19b)
-J k u = β 0 u (Ψ k )u k+1 -γ 0 u (Ψ k ), on x = 0, for u = P, N, ( 20a
)
J k u = β 1 u (V -Ψ k )u k+1 -γ 1 u (V -Ψ k ), on x = 1, for u = P, N. ( 20b
)
Let us first note that, for each value of k, the system ( 17)-( 20) is decoupled. Indeed, knowing P k and N k , we define Ψ k as the solution to ( 17), [START_REF] Jüngel | The quasi-neutral limit in the quantum drift-diffusion equations[END_REF]. As ( 17) is a classical elliptic equation, if the right-hand side belong to L 2 (0, 1), Ψ k will be defined as the solution to a variational formulation (see below) of ( 17), [START_REF] Jüngel | The quasi-neutral limit in the quantum drift-diffusion equations[END_REF]. Therefore, Ψ k ∈ H 1 (0, 1) and as H 1 (0, 1) ⊂ C 0 ([0, 1]) the boundary conditions (20) have a sense. Then, knowing Ψ k , we define u k+1 = P k+1 or N k+1 as the solution to [START_REF] Jüngel | Rigorous derivation of a hierarchy of macroscopic models for semiconductors and plasmas[END_REF], [START_REF] Markowich | The stationary semiconductor device equations[END_REF]. The equation ( 18) on u k+1 is a linear convection-diffusion equation. Using the Slotboom change of variables [START_REF] Chainais-Hillairet | Existence of solutions for a steady state corrosion of steel model[END_REF][START_REF] Markowich | The stationary semiconductor device equations[END_REF] will lead to a classical elliptic equation for which we can also define a variational formulation and the associated weak solution.

Theorem 2.1. Under the assumptions of Theorem 1.1, the semi-discrete scheme (16)-(20) admits a unique weak solution

(P k , N k , Ψ k ) 0≤k≤K with Ψ k ∈ H 1 (0, 1) for all 0 ≤ k ≤ K and P k , N k ∈ H 1 (0, 1) for all 1 ≤ k ≤ K.
Proof. The proof will be done by induction. Let k ≥ 0, we assume that P k , N k ∈ L ∞ (0, 1) (a fortiori satisfied if P k , N k ∈ H 1 (0, 1) ⊂ C 0 ([0, 1])) are known. This property is satisfied for k = 0. Then, the variational formulation associated to the problem on Ψ k (17), [START_REF] Jüngel | The quasi-neutral limit in the quantum drift-diffusion equations[END_REF] writes:

λ 2 a Ψ (Ψ k , ϕ) = b k Ψ (ϕ), ∀ϕ ∈ H 1 (0, 1), (21)
with

a Ψ (Φ, ϕ) = 1 0 ∂ x Φ∂ x ϕdx + 1 α 1 Φ(1)ϕ(1) + 1 α 0 Φ(0)ϕ(0), b k Ψ (ϕ) = λ 2 α 1 (V -∆Ψ pzc 1 ) ϕ(1) + λ 2 α 0 ∆Ψ pzc 0 ϕ(0) + 1 0 (3P k -N k + ρ hl )ϕdx.
But, the functional space H 1 (0, 1) can be equipped either with its usual norm:

u H 1 = 1 0 (∂ x u) 2 dx + 1 0 u 2 dx,
or with the following equivalent norm (see [START_REF] Chainais-Hillairet | Existence of solutions for a steady state corrosion of steel model[END_REF]):

|||u||| H 1 = 1 0 (∂ x u) 2 dx + u(0) 2 + u(1) 2 .
Indeed, using the compact injection from H 1 (0, 1) to C([0, 1]) and an adaptation of Poincaré inequality, we obtain the existence of κ > 0 such that :

1 κ u H 1 ≤ |||u||| H 1 ≤ κ u H 1 ∀u ∈ H 1 (0, 1). ( 22
)
Using the norm ||| • ||| H 1 , it is immediately clear that a Ψ (•, •) is a bilinear, continuous and coercive form in H 1 (0, 1). Moreover, b k Ψ is a continuous linear form in H 1 (0, 1). Then, the Lax-Milgram Theorem ensures that the problem (21) admits a unique solution Ψ k ∈ H 1 (0, 1). Thanks to the embedding from H 1 (0, 1) to C 0 ([0, 1]), we also have Ψ k ∈ C 0 ([0, 1]).

Let us now focus on the problem on u k+1 = P k+1 or N k+1 (18), [START_REF] Markowich | The stationary semiconductor device equations[END_REF]. We first use the classical Slotboom change of variables introduced for instance in [START_REF] Markowich | The stationary semiconductor device equations[END_REF] :

u k+1 = e -zuΨ k ζ k+1 u . The current density J k u rewrites: J k u = -e -zuΨ k ∂ x ζ k+1 u
and we can consider the variational formulation:

a k u (ζ u , ϕ) = b k u (ϕ), ∀ϕ ∈ H 1 (0, 1), (23)
with

a k u (κ u , ϕ) = ε u ∆t 1 0 e -zuΨ k κ u ϕdx + 1 0 e -zuΨ k ∂ x κ u ∂ x ϕdx +β 1 u (V -Ψ k (1))e -zuΨ k (1) κ u (1)ϕ(1) + β 0 u (Ψ k (0))e -zuΨ k (0) κ u (0)ϕ(0), b u (ϕ) = ε u ∆t 1 0 u k ϕdx + γ 1 u (V -Ψ k (1))ϕ(1) + γ 0 u (Ψ k (0))ϕ(0).
As Ψ k is continuous in [0, 1], it is bounded and e -zuΨ k is also bounded and its lower bound is strictly positive. The hypothesis (8) also ensures that the functions β 0 u , β 1 u , γ 0 u , γ 1 u are strictly positive. Therefore a k u is a bilinear, continuous and coercive form in H 1 (0, 1) and b k u is a continuous linear form in H 1 (0, 1). Then, the Lax-Milgram Theorem answers that the problem ( 23) admits a unique solution denoted ζ k+1 u ∈ H 1 (0, 1) and u k+1 = e -zuΨ k ζ k+1 u ∈ H 1 (0, 1) is the unique weak solution to ( 18), (20a), (20b). This ends the proof of Theorem 2.1.

Thanks to the semi-discrete scheme ( 16)-( 20) and Theorem 2.1, we can define an approximate solution to the corrosion model ( 1)-( 7). This approximate solution is a piecewise constant function in time. It is denoted by (P ∆t , N ∆t , Ψ ∆t ) and defined by

P ∆t (t, x) = P k (x), N ∆t (t, x) = N k (x), Ψ ∆t (t, x) = Ψ k (x), ∀x ∈ (0, 1), ∀t ∈ (t k-1 , t k ], ∀1 ≤ k ≤ K. ( 24
)
Our objective is to prove the existence of a solution to the corrosion model by passing to the limit on a sequence of approximate solutions (P ∆t , N ∆t , Ψ ∆t ) ∆t when ∆t tends to 0. Therefore, we need to establish a priori estimates which will yield the compactness of the sequence of approximate solutions.

A priori estimates

In this Section we establish the a priori estimates satisfied by the approximate solution (P ∆t , N ∆t , Ψ ∆t ): L ∞ -estimates on N ∆t and P ∆t , an L 2 (0, T, H 1 (0, 1))-estimate on Ψ ∆t and L 2 (0, T, H 1 (0, 1))-estimates on N ∆t and P ∆t .

τ = λ 2 min 1 9P m , ε N m , ( 25
)
for all k ∈ {0, ..., K} we have

0 ≤ P k (x) ≤ P m , 0 ≤ N k (x) ≤ N m , ∀x ∈ (0, 1), ( 26
)
and therefore

0 ≤ P ∆t (t, x) ≤ P m , 0 ≤ N ∆t (t, x) ≤ N m , ∀x ∈ (0, 1), ∀t ∈ (0, T ]. ( 27
)
Proof. The proof will be done by induction. First, we note that (26) is satisfied for k = 0 thanks to ( 16) and [START_REF] Gajewski | On existence, uniqueness and asymptotic behaviour of solutions of the basic equations for carrier transport in semiconductors[END_REF]. Let us now assume that the property is satisfied for some k ∈ {0, ..., K -1} , we will establish that both inequalities in (26) remain true at step k + 1. Therefore, we adapt the technique used by A. Jüngel and Y. J. Peng in [START_REF] Jüngel | A hierarchy of hydrodynamic models for plasmas. Zeroelectron-mass limits in the drift-diffusion equations[END_REF].

Let u = P, N . We first prove the positivity of u k+1 . Multiplying ( 18) by (u k+1 ) -= min(u k+1 , 0) and integrating over (0,1), we obtain, after an integration by parts:

ε u ∆t 1 0 u k+1 -2 dx - ε u ∆t 1 0 u k u k+1 - dx + 1 0 ∂ x u k+1 -2 dx + z u 2λ 2 1 0 u k+1 -2 z P P k -P m + z N N k -N m dx = f 1 u + f 0 u , ( 28
)
3.2 L 2 (0, T, H 1 )-estimate on Ψ ∆t Proposition 3.2. Under the assumptions of Proposition 3.1, there exists a constant C depending only on the data of the problem (but not on ∆t) such that:

|||Ψ k ||| H 1 ((0,1)) ≤ C, ∀1 ≤ k ≤ K, ( 34
)
and

Ψ ∆t L 2 (0,T,H 1 (0,1)) ≤ C √ T . ( 35
) Proof. Applying (21) with ϕ = Ψ k , we get a Ψ (Ψ k , Ψ k ) = b k Ψ (Ψ k )/λ 2 . Setting η = min(1, 1/α 0 , 1/α 1 ), we have a Ψ (Ψ k , Ψ k ) ≥ η|||Ψ k ||| 2 . ( 36
)
Using Young and Cauchy-Schwarz inequalities, we get

1 λ 2 |b k Ψ (Ψ k )| ≤ η 2 |||Ψ k ||| 2 + 1 2ηλ 2 1 0 (3(P k -P m ) + (N k -N m )) 2 dx + 1 2η (V -∆Ψ pzc 1 ) 2 + 1 2η (∆Ψ pzc 0 ) 2 . ( 37
)
Then, using the L ∞ -estimate (26), we deduce (34) from ( 36) and (37). The obtention of (35) is straightforward.

Remark 3.1. Let us note that, due to the compact embedding from H 1 (0, 1) to C([0, 1]), the estimate (34) imply that Ψ k (for all k) is uniformly bounded on the interval [0, 1] by a constant depending only on the data of the problem.

3.3 L 2 (0, T, H 1 )-estimates on P ∆t , N ∆t Note that using (27), we easily obtain L 2 (0, T ; L 2 (0, 1)) estimates for P ∆t and N ∆t which imply the weak convergence of a sequence of approximate solutions in L 2 (0, T ; L 2 (0, 1)). But since we have to pass to the limit also in the space derivatives of P ∆t and N ∆t , it is not sufficient and we need estimates in L 2 (0, T ; H 1 (0, 1)). They are given in Proposition 3.3.

Proposition 3.3. Under the assumptions of Proposition 3.1, there exists a constant C, depending only on the data of the problem (but not on ∆t) such that

P ∆t L 2 (0,T ;H 1 (0,1)) ≤ C, N ∆t L 2 (0,T ;H 1 (0,1)) ≤ C. ( 38
)
Proof. Let us multiply (18) by u k+1 and integrate over (0, 1). Since

u k+1 u k+1 -u k ≥ 1 2 u k+1 2 - 1 2 u k 2 , we obtain ε u 2∆t 1 0 (u k+1 ) 2 -(u k ) 2 dx + 1 0 ∂ x J k u u k+1 dx ≤ 0.
not sufficient to pass to the limit in the convective terms T 0 1 0 N ∆t ∂ x Ψ ∆t ϕdx or in the boundary terms. We need some strong convergence. To prove this result, we will use a compactness argument due to Simon [START_REF] Simon | Compact sets in the space l p (0, t; b)[END_REF]. It is based on some estimates on the time translates of the approximate solutions. Let us introduce the shift operator. For u = P, N or Ψ, we set (σ ∆t u ∆t ) (x, t) = u ∆t (x, t + ∆t).

It means that (σ ∆t u ∆t ) (x, t) = u k+1 (x) for x ∈ (0, 1) and t ∈ (t k-1 , t k ]. Proposition 4.1. Under the assumptions of Proposition 3.1, there exists a constant C depending only on the data of the problem (and not on ∆t) such that:

σ ∆t P ∆t -P ∆t L 2 (0,T ;L 2 (0,1)) ≤ C∆t, σ ∆t N ∆t -N ∆t L 2 (0,T ;L 2 (0,1)) ≤ C∆t. ( 41
)
and

σ ∆t Ψ ∆t -Ψ ∆t L 2 (0,T ;L 2 (0,1)) ≤ C∆t. ( 42
)
Proof. Using equation ( 18), Minkowsky inequality and Hölder inequality, we obtain for u = N, P

σ ∆t u ∆t -u ∆t L 2 (0,T ;L 2 (0,1)) ≤ ∆t ε u ∂ x σ ∆t u ∆t L 2 (0,T ;L 2 (0,1)) +|z u | ∂ x Ψ ∆t L 2 (0,T ;L 2 (0,1)) σ ∆t u ∆t L ∞ (0,T ;L ∞ (0,1)) , ≤ ∆t ε u u ∆t L 2 (0,T ;H 1 (0,1)) +|z u | Ψ ∆t L 2 (0,T ;H 1 (0,1)) u ∆t L ∞ (0,T ;L ∞ (0,1)) . ( 43
)
Using estimates (27), ( 35) and ( 38) in (43), we obtain (41).

By definition (σ ∆t Ψ ∆t -Ψ ∆t ) (t, x) = Ψ k+1 (x) -Ψ k (x) for t ∈ (t k-1 , t k ]. Using ( 17), [START_REF] Jüngel | The quasi-neutral limit in the quantum drift-diffusion equations[END_REF], we have

-λ 2 ∂ xx Ψ k+1 -Ψ k = 3 P k+1 -P k -N k+1 -N k , x ∈ (0, 1), (44)
Ψ k+1 -Ψ k (0) -α 0 ∂ x Ψ k+1 -Ψ k (0) = 0, (45)
Ψ k+1 -Ψ k (1) + α 1 ∂ x Ψ k+1 -Ψ k (1) = 0. (46)
Multiplying (44) by Ψ k+1 -Ψ k and integrating over (0, 1), we obtain, after an integration by parts:

λ 2 1 0 ∂ x Ψ k+1 -Ψ k 2 dx + λ 2 α 1 Ψ k+1 -Ψ k (1) 2 + λ 2 α 0 Ψ k+1 -Ψ k (0) 2 = 3 1 0 P k+1 -P k Ψ k+1 -Ψ k dx + 1 0 N k+1 -N k Ψ k+1 -Ψ k dx. ( 47
)
First, we note that the left hand side in (47) is lower bounded by λ 2 η κ Ψ k+1 -Ψ k 2 H 1 with η = min(1, 1/α 0 , 1/α 1 ) and κ defined in [START_REF] Sacco | Mixed finite volume methods for semiconductor device simulation[END_REF]. Then, using Cauchy-Schwarz and Young inequalities, we obtain the following upper bound for the right hand side in (47):

λ 2 η 2κ Ψ k+1 -Ψ k 2 H 1 + 9κ 2λ 2 η 1 0 P k+1 -P k 2 dx + κ 2λ 2 η 1 0 N k+1 -N k 2 dx It yields λ 2 η 2κ Ψ k+1 -Ψ k 2 L 2 (0,1) ≤ 9κ 2λ 2 η 1 0 P k+1 -P k 2 dx + 1 0 N k+1 -N k 2 dx .
Summing over 1 ≤ k ≤ K and using the estimates (41), we deduce (42). It concludes the proof of Proposition 4.1.

Finally, thanks to the a priori estimates (Propositions 3.1, 3.2 and 3.3) and the estimates on the time translates (Propositions 4.1 ), we obtain the convergence of a subsequence of approximate solutions. Proposition 4.2. Under the assumptions of Theorem 1.1, we consider a sequence of approximate solutions (P ∆t , N ∆t , Ψ ∆t ) ∆t defined by the semi-discrete scheme ( 16)-( 20) and [START_REF] Schmeiser | A singular perturbation analysis of reverse biased pn-junctions[END_REF]. Then, there exists (P, N ) ∈ L ∞ ([0, T ] × (0, 1)) ∩ L 2 (0, T, H 1 (0, 1)) and Ψ ∈ L 2 (0, T, H 1 (0, 1)) such that, up to a subsequence, (u = P, N )

• Ψ ∆t → Ψ strongly in L 2 (0, T ; C([0, 1])), • u ∆t → u strongly in L 2 (0, T ; C([0, 1])), • ∂ x u ∆t ⇀ u weakly in L 2 (0, T ; L 2 (0, 1)), σ ∆t u ∆t ⇀ u weakly in L 2 (0, T ; L 2 (0, 1)), • ∂ x Ψ ∆t ⇀ ∂ x Ψ weakly in L 2 (0, T ; L 2 (0, 1)).
Proof. Firstly, due to Propositions 3.2 and 3.3, the sequences (Ψ ∆t ) ∆t and (u ∆t ) ∆t are bounded in L 2 (0, T ; H 1 (0, 1)). Then, there exist Ψ and u in L 2 (0, T ; H 1 (0, 1)) such that Ψ ∆t tends to Ψ and u ∆t tends to u in L 2 (0, T ; L 2 (0, 1)). Therefore, up to a subsequence, we obtain the weak convergence of (∂ x Ψ ∆t ) ∆t to ∂ x Ψ and of (∂ x u ∆t ) ∆t to ∂ x u in L 2 (0, T ; L 2 (0, 1)).

Secondly, to obtain the strong convergence in L 2 (0, T ; C([0, 1])) of sequences (Ψ ∆t ) ∆t and (u ∆t ) ∆t , we use a compactness argument. More precisely, we use Theorem 5 given by J. Simon in [START_REF] Simon | Compact sets in the space l p (0, t; b)[END_REF]. Indeed, we have:

H 1 (0, 1) ⊂ C([0, 1]) ⊂ L 2 (0, 1)
with a compact injection from H 1 (0, 1) to C([0, 1]) (due to the one dimension). Moreover due to Propositions 3.2, 3.3 and 4.1 the sequence (f ∆t ) ∆t (for f = Ψ, u) is bounded in L 2 (0, T ; H 1 (0, 1)) and σ ∆t f ∆t -f ∆t L 2 (0,T ;L 2 (0,1)) -→ 0 uniformly when ∆t → 0.

Then the sequence (f ∆t ) ∆t is relatively compact in L 2 (0, T ; C([0, 1])), which gives the strong convergence of (f ∆t ) ∆t in L 2 (0, T ; C([0, 1])).

Finally let us note that due to proposition 3.1, sequences (P ∆t ) ∆t and (N ∆t ) ∆t are bounded in L ∞ ([0, T ]×(0, 1)) and we obtain their ⋆-weak convergence in L ∞ ([0, T ]×(0, 1)). Then their limits belong to L ∞ ([0, T ] × (0, 1)). This completes the proof of proposition 4.2. for all ϕ ∈ L 2 (0, T ; H 1 (0, 1)). It is the expected result. Furthermore, passing to the limit in (26) gives for all t ∈ [0, T ] and all x ∈ (0, 1) 0 ≤ P (t, x) ≤ P m , 0 ≤ N (t, x) ≤ N m .

This ends the proof of Theorem 1.1.

Conclusion

In this paper, we prove the existence of a weak solution to a simplified corrosion model. The boundedness of the electronic and cationic densities is also established. Let us note that the dimensionless parameter ε in (1b) is in practice very small because it is the ratio of the mobility coefficients of cations and electrons. It could be set to 0 in the model. Then, the question of existence of a solution to the resulting model is still an open problem. Future work will focus on this question. The full understanding of the corrosion models reduced to a fixed domain will then permit to deal with the DPCM model on a moving domain (see [START_REF] Bataillon | Numerical methods for the simulation of a corrosion model with moving oxide layer[END_REF]).

3. 1 L 3 . 1 .

 131 ∞ -estimates on P ∆t , N ∆t Proposition Let us assume the hypothesis of Theorem 1.1 and that the time step ∆t verifies ∆t ≤ τ with

(

  3P -N + ρ hl )ϕdxdt,

But, using the induction hypothesis and the fact that z N z P ≤ 0, we have:

and (28) implies 1 ∆t

We remark that the condition on ∆t [START_REF] Simon | Compact sets in the space l p (0, t; b)[END_REF] and the induction hypothesis ensure that

Let us now study the sign of f 0 u and f 1 u . Therefore,we introduce, as in [START_REF] Bataillon | Numerical methods for the simulation of a corrosion model with moving oxide layer[END_REF], the following functions :

As shown in [START_REF] Bataillon | Numerical methods for the simulation of a corrosion model with moving oxide layer[END_REF], ξ 0 u and ξ 1 u are nonpositive functions under the hypothesis [START_REF] Gajewski | On the uniqueness of solutions to the drift-diffusion model of semiconductors devices[END_REF]. Using the boundary conditions [START_REF] Jüngel | The quasi-neutral limit in the quantum drift-diffusion equations[END_REF], [START_REF] Markowich | The stationary semiconductor device equations[END_REF], we have:

As the functions (β i u , γ i u ) i=0,1 are nonnegative on R and the functions (ξ i u ) i=0,1 are nonpositive on R for u = N, P , we deduce f 0 u ≤ 0 and f 1 u ≤ 0. Then, (29) and (30) imply that u k+1 - = 0 on (0, 1) and we obtain the positivity of u k+1 .

Let us now prove that u k+1 ≤ u m . Therefore, we multiply ( 18) by (u k+1 -u m ) + = max((u k+1 -u m ), 0) and we integrate over (0, 1). It yields:

, integrating by parts and using [START_REF] Jüngel | A hierarchy of hydrodynamic models for plasmas. Zeroelectron-mass limits in the drift-diffusion equations[END_REF], we get:

Again, using the induction hypothesis and the fact that z N z P < 0, we have

Combined with (31) and (32), it yields:

with

As we did for f 0 u and f 1 u , we use the boundary conditions [START_REF] Jüngel | The quasi-neutral limit in the quantum drift-diffusion equations[END_REF] and [START_REF] Markowich | The stationary semiconductor device equations[END_REF] and the definition of the functions (ξ i u ) i=0,1 in order to rewrite g 0 u and g 1 u :

It is now clear that g 0 u ≤ 0 and g 1 u ≤ 0 and therefore, thanks to the hypothesis on the time step [START_REF] Simon | Compact sets in the space l p (0, t; b)[END_REF], the inequality (33) implies (u k+1 -u m ) + = 0 and u k+1 ≤ u m for all x ∈ (0, 1). It ends the proof of Proposition 3.1.

But, integrating by parts the term 1 0 ∂ x J k u u k+1 dx and using [START_REF] Jüngel | A hierarchy of hydrodynamic models for plasmas. Zeroelectron-mass limits in the drift-diffusion equations[END_REF], we get

with

The sign properties of the functions ξ i u , β i u ,γ i u for i = 0, 1, ensure that

Then, we can multiply (39) by ∆t and sum it over k ∈ {0, . . . , K -1}. Using the L ∞estimate (26), we get

Thanks to [START_REF] Gajewski | On existence, uniqueness and asymptotic behaviour of solutions of the basic equations for carrier transport in semiconductors[END_REF], (26), the continuity of the functions (γ i u ) i=0,1 and Remark 3.1, we obtain

where C depends only on the data of the problem and not on ∆t. As the L 2 (0, T ; L 2 (0, 1))estimate of u ∆t is a straightforward consequence of the L ∞ -estimate (26), it concludes the proof of Proposition 3.3.

Proof of Theorem 1.1

In this Section, we first establish some compactness results and the convergence of a sequence of approximate solutions (N ∆t , P ∆t , Ψ ∆t ) ∆t when ∆t tends to 0. Then, by passing to the limit in the scheme ( 16)-( 20), we prove that the triple (N, P, Ψ) obtained as the limit of the sequence (N ∆t , P ∆t , Ψ ∆t ) ∆t is a weak solution to (1)-(7).

Compactness results

The estimates proven in Propositions 3.2 and 3.3 imply the weak convergence of a sequence of approximate solutions (N ∆t ) ∆t , (P ∆t ) ∆t and (Ψ ∆t ) ∆t in L 2 (0, T ; H 1 (0, 1)), when ∆t tends to zero. But, the weak convergence of (N ∆t ) ∆t and (Ψ ∆t ) ∆t in L 2 (0, T, H 1 (0, 1)) is

Passage to the limit

In order to conclude the proof of Theorem 1.1, it remains to prove that (P, N, Ψ) obtained as a limit of the semi-discrete scheme (see Proposition 4.2) satisfies ( 13) and ( 14).

Let us first recall the results of the previous sections. Theorem 2.1 gives the existence of a unique solution to the problem (u = N, P )

The variational formulation of ( 48)-( 53) writes : for all ϕ ∈ L 2 (0, T ; H 1 (0, 1)), + β 0 u (Ψ(0))u(0) -γ 0 u (Ψ(0)) ϕ(0) dt = 0, for u = N, P,


David Lopez

Bruno Lévy

Dynamic Mesh Optimization for Free Surfaces in Fluid Simulation

We propose a new free surface fluid front tracking algorithm. Based on Centroidal Voronoi Diagram optimizations, our method creates a serie of either isotropic or anisotropic meshes that confroms with an evolving surface.

Introduction

We consider free-surface fluid simulations that operate on a per time-step basis. Each step requires to compute the velocity field governing the fluid motion, to track its surface and to accurately represent it. We do not discuss how to simulate the physics and we assume that we have an everywhere and anytime known velocity field to move surface vertices according with an accurate advection method, eg Runge Kutta.

In this paper, we focus on accurately sampling and optimizing the advected surface regarding its geometry and combinatorics. Assuming curl-free velocity fields, like incompressible fluid flows, we will also minimize volume variations due to polygonal approximation. This, in a direct manner, without any other representation like tetrahedral mesh or sub-grid.

After a short review of previous work (section 2), we explain our remeshing operations (section 3) and discuss first results and perspectives (section 4).

previous work and foundations

Dynamic interfaces : Eulerian methods such as Level Set [START_REF] Enright | Using the particle level set method and a second order accurate pressure boundary condition for free surface flows[END_REF] or Volume of Fluid [START_REF] Hirt | Volume of fluid /vof/ method for the dynamics of free boundaries[END_REF] are populars because they implicitly handle topological changes like merging or breaking of droplets. On the other hand, Lagrangian methods like ours or "El Topo" [START_REF] Brochu | Robust topological operations for dynamic explicit surfaces[END_REF] provide accurate explicit surfaces that could be directly used (rendering...). We refer to the reviews in [START_REF] Brochu | Robust topological operations for dynamic explicit surfaces[END_REF] and [START_REF] Wojtan | Liquid simulation with mesh-based surface tracking[END_REF].

Evaluation : In order to evaluate the accuracy of our front tracker, we use velocity fields that are defined in closed form by an explicit formula. They give anywhere and anytime known velocity fields and do not suffer of physics simulation errors. Enright's [START_REF] Enright | Using the particle level set method and a second order accurate pressure boundary condition for free surface flows[END_REF] and Curl-noise [START_REF] Brochu | Robust topological operations for dynamic explicit surfaces[END_REF] tests are commonly used (see respectively figures 1 and 5).

Surface meshing : Our goal is to construct a mesh of good quality at each time step. We propose to optimize the placement of the mesh vertices.

Our optimization is based on the Restricted Voronoi Diagram, see figure 2 and [START_REF] Yan | Isotropic remeshing with fast and exact computation of restricted voronoi diagram[END_REF][START_REF] Lévy | Variational anisotropic surface meshing with voronoi parallel linear enumeration[END_REF] such distribution from a random sampling is to iteratively replace samples by their cell centroids as shown in figure 3 (Lloyd algorithm [START_REF] Stuart | Least squares quantization in PCM[END_REF]).

CV T -computing can also be viewed as an optimization problem [START_REF] Du | Centroidal voronoi tessellations: Applications and algorithms[END_REF] of a C 2 function [START_REF] Liu | On centroidal voronoi tessellation -energy smoothness and fast computation[END_REF] for which quasi-newton solvers could be used to converge faster.

Algorithm :

Each advected surface is remeshed with the following algorithm. Some details are given below. [START_REF] Edelsbrunner | Triangulating Topological Spaces[END_REF], some conditions are required to ensure that there exists an homeomorphism between S t+1 and the advected surface.

In addition, for free-surface fluid simulation, we need to detect merging and splitting events where the field is discontinuous and modify the surface topology accordingly. The management of the topology is realized on special Vor P |S cell configurations, with combinatorial corrections and/or local refinements (see [START_REF] Yan | Isotropic remeshing with fast and exact computation of restricted voronoi diagram[END_REF]). Surface volumetric optimization : For curl-free velocity fields, the fluid volume must remain constant throughout the simulation. We assume that S t+1 and S adv t are smooth and geometrically near to each other. As in 2D example given in fig. 3, for each voronoi cell we define a polygon (polyhedron in 3D) with Vor P |St+1 and Vor P |S adv t cell facets and Voronoi planes. Hence we can compute a local signed volume attached to each Voronoi seed. By minimizing the sum of squared local volumes (shaded triangles in figure 4), we improve the accuracy of the surface tracking. Anisotropic meshes are more suitable because they allow to represent very thin and sharp features with a small number of points. In addition, our local volumetric optimization further reduces physical errors involved by remeshing. Volume is well preserved throughout the simulation. Less than 0.2% of volume is lost during 3D Enright's test (same running conditions than in [START_REF] Brochu | Robust topological operations for dynamic explicit surfaces[END_REF]).

Merging and splitting should also be improved to complete intricate fluid simulations. We plan to couple our front tracker with a free surface fluid simulator [START_REF] Brochu | Matching fluid simulation elements to surface geometry and topology[END_REF].

Fig. 1 .

 1 Fig. 1. Enright's test at steps 0, 50, 75, 100, 150, 200, 225, 250 and 300 (top) and details of the isotropic and anisotropic 150th step mesh (resp. bottom right and left).

Fig. 2 .

 2 Fig. 2. from left to right : samples P on a surface S, Vor P |S , duals and Del P |S .

  for details. Given a 2-manifold surface S ∈ IR 3 and a sample set P ∈ S, the Voronoi Diagram of P restricted to S is defined as Vor P |S = {Ω i|S } with : Ω i|S = {p ∈ S | d(p, p i ) < d(p, p j ) ∀i, j ∈ P } where d(a, b) denotes the euclidean distance between a and b. A [restricted] Voronoi tessellation is said to be centroidal ([R]CV T for short) if each seed corresponds to it's cell centroid. A simple way to compute

Fig. 3 .

 3 Fig. 3. on the left, the Voronoi Tessellation of a point set (dots) restricted to a sphere, cells centroids are represented with squares. The three others images are respectively the result of 2, 5 and 20 Lloyd iterations.

Fig. 4 .

 4 Fig. 4. Local volume difference minimization : (a) a sampled surface (dotted line) and its remeshing (plain line), (b) surfaces with sample voronoi diagram superimposed, (c) highlighted volume loss and gains (resp. black and grey areas), (d) minimized volume differences, (e) original, remeshed and optimized surfaces (resp. dotted, grey and black plain lines)

4

  Results and future work Enright's and Curlnoise test screenshots are shown in figures 1 and 5 respectively. The mesh details in the first figure show isotropic and anisotropic meshes.

Fig. 5 .

 5 Fig. 5. Curl-noise test screenshots.

  Algorithm 1: simulation step Input: a 2-manifold surface St representing the fluid surface at time t Data: min lg, max lg : edge length bounds computed from S0, algo, nb iter, dim : RCVT parameters, nb volum optim iter Result: a 2-manifold surface St+1 representing the fluid surface at time t + 1 |S cell conformations : we need to build a valid surface (S t+1 ) from Vor P |S . As explained in

	// surface vertices advection
	S ← advect(St, velocity f ield)
	// surface adaptive sampling
	P ← V ertices(S)
	end
	end
	// sampling optimization
	P ← Compute RCV T (P, S, dim, algo, nb iter) // algo={Lloyd|q-newton}
	// surface building
	check and fix Vor P |S cell conformations	// detailed below
	St+1 ← duals(Vor P |S ) ;
	// surface optimization
	minimize local volume differences between S adv t	and St+1 // detailed below
	check and fix Vor P

for edge(pi, pj) ∈ edges(S) do if length(edge) < min length or length(edge) > max length then add (pi + pj) * 0.5 to P if length(edge) < min length then remove pi and pj from P

Acknowledgments

This work is partly supported by the European Research Council (project GOODSHAPE, ERC-StG-205693) and the ANR (project MORPHO).


Dang Van Nguyen

Denis Grebenkov

Cyril Poupon

Denis Le Bihan

Jing-Rebecca Li

Effective diffusion tensor computed by homogenization

Introduction

Diffusion MRI can give useful information on cellular structure and structural change (for a review see [1]). We show that the effective diffusion tensor obtained by mathematical homogenization theory (see e.g. [2,3]) is a good approximation to the long time apparent diffusion tensor under realistic DMR scanning conditions for both isotropic and anisotropic diffusion and general geometries. The homogenized diffusion tensor is obtained by solving three steady-state Laplace equations, which is a more computationally efficient approach than long time simulation in the time domain, either via Monte-Carlo simulation or numerical solution of the time-dependent Bloch-Torrey PDE.

Theory

In the two-compartment model, we consider the two compartments, i

Ω and e Ω , to be the ensemble of cells and the extra-cellular compartment, respectively. The two compartments have the same intrinsic diffusion coefficient D. The cell membrane is modeled by an infinitely thin permeable interface characterized by permeability κ. Given the diffusion gradient with profile f(t) and gradient strength γ

/ : q g r r =
, where γ is the gyro-magnetic ratio, the DMRI signal attenuation is ) , ( t q r Ψ , from which we define the apparent diffusion tensor D A from the Taylor expansion in q r :

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ - = Ψ ∫ ∫ 4 0 2 ) ( ) , ( log q O ds s f du q D q t q t u s A T r r r r
. The long time limit of D A can be well approximated by the effective diffusion tensor of homogenization,

Results and discussion

Simulated DMRI signals were obtained for the PGSE sequence with δ=10ms and Δ=10ms, 50ms and 190ms in two gradient directions

2 / q q r r =[1,0],
and [0,1] (2D case) and δ=10ms and Δ=10ms and 90ms in three gradient directions with v i =0,66, S/V=0,34μm -1 . The average radius of cells is 3.9μm for 2D and 8,8μm for 3D. We simulated the DMRI signal for D=2,8.10 -3 μm 2 /μs and κ=10 -5 μm/μs and κ=10 -4 μm/μs, giving computed ADCs of between 0,5.10 -3 and 2,2.10 -3 μm 2 /μs (3D case). In Table 1 and 2 we see the simulated apparent diffusion tensor D A approaches the steady-state tensor D eff computed by the homogenization method described above. The convergence of D A to D eff is faster at higher permeability. The calculation of D eff for 3D took between 10 minutes to two hours for spatial discretizations of between 80x80x80 and 160x160x160 on a Dell Latitiude E6410 laptop computer (Intel(R) Core(TM) i7 CPU M640 @ 2,8GHz).

q D q A T r r q D q eff T r r q D q A T r r q D q eff T
r r q r δ=10,Δ=10 δ=10,Δ=50 δ=10,Δ=190 δ=10,Δ=10 δ=10,Δ=50 δ=10,Δ=190

(1,0) 1,1.10 -3 8,3.10 -4 7,5.10 -4 7,1.10 -4 8,3.10 -4 4,5.10 -4 3,2. 10 -4 2,6.10 -4 (0,1) 1,1.10 -3 8,8.10 -4 7,9.10 -4 7,6.10 -4 8,9. 10 -4 4,8. 10 -4 3,3. 10 -4 2,7.10 -4 Table 1. The simulated apparent diffusion tensor D A approaches the steady-state tensor D eff computed by homogenization (2D case). κ=10 -4 μm/μs κ=10 -5 μm/μs

q D q A T r r q D q eff T r r q D q A T r r q D q eff T r r q r δ=10,Δ=10 δ=10,Δ=90 δ=10,Δ=10 δ=10, Δ=90
(1,0,0) 1,7.10 -3 1,6.10 -3 1,5.10 -3 1,4.10 -3 1,0.10 -3 0,8.10 -3 (0,1,0) 2,2.10 -3 2,1.10 -3 2,0.10 -3 2,0.10 -3 1,5.10 -3 1,2.10 -3 (0,0,1) 1,2.10 -3 1,0.10 -3 1,0.10 -3 0,8.10 -4 0,5.10 -3 0,4.10 -3 Table 2. The simulated apparent diffusion tensor D A approaches the steady-state tensor D eff computed by homogenization (3D case).

[ ] [ ] [ ] 3 2 1 , 0 , 0 , 0 l l l × × = Ω
which contains a representation sample of the cellular structure . In an nearly isotropic medium, D eff is nearly diagonal and all the diagonal entries are close to each other. But in the general anisotropic case, the eigenvalues of D eff are not equal to each other.

  0,0], [0,1,0] and [0,0,1] (3D case) by numerically solving the two-compartment Bloch-Torrey partial differential equation on a sample containing numerous convex cells in two dimensions with volume fraction v i =0,88, average surface to volume ratio S/V=0,51μm -1 (Fig 1) and convex-shaped cells in three dimensions (Fig 2)

Fig 1 .

 1 Fig 1. Convex cells with average radius 3.9μm, volume fraction v i =0,88

Fig 2 .

 2 Fig 2. Convex cells with average radius 8,8μm and volume fraction v i =0,66.

  [1] LeBihan (2007) Phys Med Bio 52. [2] Bensoussan et al. (1978) Asymptotic analysis for periodic structures, North-Holland, Amsterdam. [3] Cheng et al (1997) Proc Math Phy Engin Sci 453:145-161. [4] Novikov et al. (2011) Nature Physics 7, 508-514. unit vector in the kth coordinate direction and the unknown function v j can be found by solving the Laplace equation on the right over the box

Proc. Intl. Soc. Mag. Reson. Med. 20 (2012)


Jean-Antoine Désidéri
email: jean-antoine.desideri@inria.fr

MUTIPLE-GRADIENT DESCENT ALGORITHM FOR MULTIOBJECTIVE OPTIMIZATION

Keywords: multiobjective optimization, multidisciplinary, cooperative algorithms, Pareto optimality

The steepest-descent method is a well-known and effective single-objective descent algorithm when the gradient of the objective function is known. Here, we propose a particular generalization of this method to multi-objective optimization by considering the concurrent minimization of n smooth criteria {J i } (i = 1, . . . , n). The novel algorithm is based on the following observation: consider a finite set of vectors {u i } (u i ∈ R N , n ≤ N); in the convex hull of this family, there exists a unique element of minimal norm, say ω ∈ R N ; then, the scalar product of ω with any vector in the convex hull, and in particular, with any u i , is at least equal to ω 2 ≥ 0. Applying this to the objective-function gradients (u i = ∇J i ), we conclude that either ω = 0, and the current design point belongs to the Pareto set, or -ω is a descent direction common to all objective functions. We propose to construct a fixed-point iteration in which updates of the element ω are used as successive directions of search. This method converges to a point on the Pareto set. This result applies to both finite-dimensional and functional design spaces. Numerical illustrations have been provided in both cases using either analytical objective functions, or (discretized) functionals in [9] [5]. Here, following [6], a domain-decomposition framework is used to illustrate the necessity, in a (discretized) functional setting, to scale the gradients appropriately.

INTRODUCTION

Classically, in multi-objective optimization, several fundamental concepts are introduced: dominance in efficiency between design points, Pareto set, made of non-dominated solutions in design space, and Pareto front, its image in function space [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF]. The Pareto front provides the designer with the system maximum attainable performance. For complex systems, in particular those governed by partial-differential equations, a computational challenge is to devise algorithms permitting to identify numerically the Pareto set, or the most useful portions of it. In this respect, certain evolutionary strategies have been adapted to achieve this goal, and appear to provide the most robust algorithms. NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] is certainly one of the most widely-used methods for this purpose.

In the context of differentiable optimization, one would expect adequate strategies based on gradient evaluations to also be capable of capturing Pareto fronts, with less generality or robustness, but often far greater cost efficiency. However classical techniques, such as minimizing an agglomerated criterion, or one criterion at a time under the constraints of the others, are limited by hypotheses on the pattern of the Pareto front w.r.t. convexity and, or continuity. The Multiple-Gradient Descent Algorithm (MGDA ), originally introduced in [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF], and again formalized in [START_REF] Désidéri | Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF], is based on a very general principle permitting to define at each iteration, a descent direction common to all criteria. This direction is the support of the minimum-norm element in the convex hull of the local gradients. The efficacy of the algorithm to identify the Pareto front has been demonstrated in [START_REF] Zerbinati | Comparison between MGDA and PAES for Multi-Objective Optimization[END_REF] [START_REF] Désidéri | Multiple-gradient descent algorithm (MGDA) for multiobjective optimization[END_REF] in a test-case in which the Pareto front was non-convex. The method was compared in cost efficiency with an evolutionary strategy, and was found to offer very promising performance.

More recently, a variant, MGDA-II , has been proposed in which the descent direction is calculated by a direct procedure [START_REF] Désidéri | MGDA II: A direct method for calculating a descent direction common to several criteria[END_REF], which provides a valuable simplification of implementation.

Here, MGDA is tested in the fairly different context of a simulation by domain partitioning, as a technique to match the different interface components concurrently. For this, the very simple test-case of the finite-difference discretization of the Dirichlet problem over a square is considered. Full details have been provided in [START_REF] Désidéri | Application of MGDA to domain partitioning[END_REF]. The study aims at assessing the performance of MGDA in a discretized functional setting. One of the main teachings is the necessity, here found imperative, to scale the gradients appropriately.

DIRICHLET PROBLEM, DOMAIN PARTITIONING AND MATCHING DEFECTS

We consider the model problem consisting in solving Laplace's equation,

-∆u = f (Ω) (1)
over the square

Ω = [-1, 1] × [-1, 1] (2)
subject to homogeneous boundary conditions:

u = 0 (Γ = ∂Ω) (3)
For this, the domain Ω is partitioned in four disjoint sub-domains:

         Ω 1 = [0, 1] × [0, 1] Ω 2 = [-1, 0] × [0, 1] Ω 3 = [-1, 0] × [-1, 0] Ω 4 = [0, 1] × [-1, 0] (4)
with the following interfaces supported by the coordinate axes:

         γ 1 = { 0 ≤ x ≤ 1 ; y = 0 } γ 2 = { x = 0 ; 0 ≤ y ≤ 1} γ 3 = { -1 ≤ x ≤ 0 ; y = 0 } γ 4 = { x = 0 ; -1 ≤ y ≤ 0 } (5)
along which the following Dirichlet controls are applied:

         γ 1 : u = v 1 (x) γ 2 : u = v 2 (y) γ 3 : u = v 3 (x) γ 4 : u = v 4 (y) (6)
(see Fig. 1).

-

1 0 1 1 -1 Γ = ∂Ω : u = 0 γ 1 : u = v 1 (x) γ 3 : u = v 3 (x) γ 2 : u = v 2 (y) γ 4 : u = v 4 (y) Ω 1 Ω 2 Ω 4 Ω 3 x y
u = (u 1 , u 2 , u 3 , u 4 ) t (7)
(where the superscript t stands for transposition) is made based on a certain setting of the interface controls:

v = (v 1 , v 2 , v 3 , v 4 ) t (8)
by solving, possibly in parallel, the partial problems for i = 1, ..., 4:

         -∆u i = f u i = 0 u i = v i u i = v i+1 (Ω i ) (Γ ∩ ∂Ω i ) (γ i ) (γ i+1 ) (9)
In the above, and all throughout thereafter, by periodicity, the index i is understood modulo 4, i.e. u 5 = u 1 , γ 5 = γ 1 , etc. Since the interface controls are of Dirichlet type, the resulting compound solution u is continuous, and its derivative along each interface is also continuous. However, in general, unless the specified controls v i 's are equal to the restrictions of the global solution, the normal derivatives exhibit jump discontinuities, s i 's. Here, each interface is supported by a coordinate axis, and we adopt the following sign convention: on the interface γ i which is supported by the x (resp. y) axis for i = 1 and 3 (resp. 2 and 4), the jump, s i (x) (resp. s i (y)), is defined as the increment of the partial derivative ∂u/∂y (resp. ∂u/∂x) as y (resp. x) goes from 0 -to 0 + . Thus:

• over γ 1 (0 ≤ x ≤ 1 ; y = 0): s 1 (x) = ∂u ∂y (x, 0 + ) - ∂u ∂y (x, 0 -) = ∂u 1 ∂y - ∂u 4 ∂y (x, 0); • over γ 2 (x = 0 ; 0 ≤ y ≤ 1): s 2 (y) = ∂u ∂x (0 + , y) - ∂u ∂x (0 -, y) = ∂u 1 ∂x - ∂u 2 ∂x (0, y);
• over γ 3 (-1 ≤ x ≤ 0 ; y = 0): s 3 (x) = ∂u ∂y (x, 0 + ) -∂u ∂y (x, 0 -) = ∂u 2 ∂y -∂u 3 ∂y (x, 0);

• over γ 4 (x = 0 ; -1 ≤ y ≤ 0): s 4 (y) = ∂u ∂x (0 + , y) - ∂u ∂x (0 -, y) = ∂u 4 ∂x - ∂u 3 ∂x (0, y).
The above local measures of the defect in matching conditions can be associated with global measures defined as separate functionals defined over the interfaces:

J i = γ i 1 2 s 2 i w dγ i (10)
that is, explicitly:

J 1 = 1 0 1 2 s 1 (x) 2 w(x) dx J 2 = 1 0 1 2 s 2 (y) 2 w(y) dy (11)
J 3 = 0 -1 1 2 s 3 (x) 2 w(x) dx J 4 = 0 -1 1 2 s 4 (y) 2 w(y) dy (12)
Here, w(t) (t ∈ [0, 1]) is an optional weighting function, and w(-t) = w(t).

The jump s i depends on the partial solutions u i-1 and u i , which themselves, depend on (v i-1 , v i ) and (v i , v i+1 ) respectively. Hence, the integral J i depends on all four sub-controls except v i+2 . Nevertheless, these four integrals are thereafter considered as functionals of v.

The coordination problem is to devise a convergent iteration on the control v to satisfy in the limit the matching conditions

J 1 = J 2 = J 3 = J 4 = 0 (13)
To achieve this, the functional gradients are firstly established using the classical adjoint-equation approach, and several strategies are proposed and tested numerically.

ADJOINT PROBLEMS AND FUNCTIONAL GRADIENTS

A first calculation is made based on the four-component control v = (v 1 , v 2 , v 3 , v 4 ) t , resulting in the compound solution u = (u 1 , u 2 , u 3 , u 4 ) t , and the multi-component criterion J = (J 1 , J 2 , J 3 , J 4 ) t .

Then, one perturbs the control v of

v ′ = (v ′ 1 , v ′ 2 , v ′ 3 , v ′ 4 ) t = δv 1 (x), δv 2 (y), δv 3 (x), δv 4 (y) t (14)
Consequently, the compound solution u is perturbed of

u ′ = (u ′ 1 , u ′ 2 , u ′ 3 , u ′ 4 ) t = δu 1 (x, y), δu 2 (x, y), δu 3 (x, y), δu 4 (x, y) t ( 15
)
in which u ′ i is the solution of the following linearized system posed on sub-domain Ω i :

         ∆u ′ i = 0 u ′ i = 0 u ′ i = v ′ i u ′ i = v ′ i+1 (Ω i ) (Γ ∩ ∂Ω i ) (γ i ) (γ i+1 ) (16)
These state perturbations induce the following functional perturbations:

J ′ i = δJ i = γ i s i s ′ i w dγ i (17) in which s ′ i = δs i . But:                          s 1 s ′ 1 = ∂u 1 ∂y - ∂u 4 ∂y ∂u ′ 1 ∂y - ∂u ′ 4 ∂y (x, 0) s 2 s ′ 2 = ∂u 1 ∂x - ∂u 2 ∂x ∂u ′ 1 ∂x - ∂u ′ 2 ∂x (0, y) s 3 s ′ 3 = ∂u 2 ∂y - ∂u 3 ∂y ∂u ′ 2 ∂y - ∂u ′ 3 ∂y (x, 0) s 4 s ′ 4 = ∂u 4 ∂x - ∂u 3 ∂x ∂u ′ 4 ∂x - ∂u ′ 3 ∂x (0, y) (18)
We now recall Green's formula for two functions φ and ψ ∈ H 2 (̟), for a simply-connected planar domain ̟ with smooth enough boundary. Since

̟ φ ∆ψ = ̟ φ ∇.(∇ψ) = ∂̟ φ ψ n - ̟ ∇φ.∇ψ ( 19
)
where n is outward unit vector normal to the boundary ∂̟, and ψ n = ∂ψ/∂n the normal derivative, and symmetrically,

̟ ψ ∆φ = ̟ ψ ∇.(∇φ) = ∂̟ ψ φ n - ̟ ∇ψ.∇φ (20)
where φ n = ∂φ/∂n is the normal derivative of φ, one has:

̟ (φ ∆ψ -ψ ∆φ) = ∂̟ (φ ψ n -ψ φ n ) (21)
Consider the following eight adjoint systems (two per sub-domain):

     ∆p i = 0 p i = 0 p i = s i w (Ω i ) (∂Ω i \γ i ) (γ i )      ∆q i = 0 q i = 0 q i = s i+1 w (Ω i ) (∂Ω i \γ i+1 ) (γ i+1 ) (22)
Then apply Green's formula, (21), to the eight cases corresponding to

̟ = Ω i (i = 1, 2, 3, 4), φ = p i or q i , ψ = u ′ i ( 23
)
so that ∆φ = ∆ψ = 0 in ̟ giving:

∂Ω i φ ψ n = ∂Ω i ψ φ n (24)
On the boundary ∂̟ = ∂Ω i :

• φ = 0 except for φ = p i = s i w along γ i , and φ = q i = s i+1 w along γ i+1 ;

• ψ = u ′ i =      v ′ i along γ i v ′ i+1 along γ i+1 0 along Γ ∩ Ω i
Hence (24) reduces to:

γ i s i u ′ in w = γ i p in v ′ i + γ i+1 p in v ′ i+1 (25)
for φ = p i , and to:

γ i+1 s i+1 u ′ in w = γ i q i n v ′ i + γ i+1 q i n v ′ i+1 (26)
for φ = q i . These two equations are particularized to sub-domains Ω i as follows.

Sub-domain

Ω 1 : on γ 1 : n = n 14 = - ; on γ 2 , n = n 12 = -ı. Thus (25)-(26) write:          1 0 s 1 (x) - ∂u ′ 1 ∂y (x, 0) w(x) dx = 1 0 - ∂p 1 ∂y (x, 0) v ′ 1 (x) dx + 1 0 - ∂p 1 ∂x (0, y) v ′ 2 (y) dy 1 0 s 2 (y) - ∂u ′ 1 ∂x (0, y) w(y) dy = 1 0 - ∂q 1 ∂y (x, 0) v ′ 1 (x) dx + 1 0 - ∂q 1 ∂x (0, y) v ′ 2 (y) dy (27) Sub-domain Ω 2 : on γ 2 : n = n 21 = + ı; on γ 3 , n = n 23 = - . Thus (25)-(26) write:          1 0 s 2 (y) ∂u ′ 2 ∂x (0, y) w(y) dy = 1 0 ∂p 2 ∂x (0, y) v ′ 2 (y) dy + 0 -1 - ∂p 2 ∂y (x, 0) v ′ 3 (x) dx 0 -1 s 3 (x) - ∂u ′ 2 ∂y (x, 0) w(x) dx = 1 0 ∂q 2 ∂x (0, y) v ′ 2 (y) dy + 0 -1 - ∂q 2 ∂y (x, 0) v ′ 3 (x) dx (28) Sub-domain Ω 3 : on γ 3 : n = n 32 = +  ; on γ 4 , n = n 34 = + ı. Thus (25)-(26) write:          0 -1 s 3 (x) ∂u ′ 3 ∂y (x, 0) w(x) dx = 0 -1 ∂p 3 ∂y (x, 0) v ′ 3 (x) dx + 0 -1 ∂p 3 ∂x (0, y) v ′ 4 (y) dy 0 -1 s 4 (y) ∂u ′ 3 ∂x (0, y) w(y) dy = 0 -1 ∂q 3 ∂y (x, 0) v ′ 3 (x) dx + 0 -1 ∂q 3 ∂x (0, y) v ′ 4 (y) dy (
         0 -1 s 4 (y) - ∂u ′ 4 ∂x (0, y) w(y) dy = 0 -1 - ∂p 4 ∂x (0, y) v ′ 4 (y) dy + 1 0 ∂p 4 ∂y (x, 0) v ′ 1 (x) dx 1 0 s 1 (x) ∂u ′ 4 ∂y (x, 0) w(x) dx = 0 -1 - ∂q 4 ∂x (0, y) v ′ 4 (y) dy + 1 0 ∂q 4 ∂y (x, 0) v ′ 1 (x) dx (
                                                                 J ′ 1 = 1 0 s 1 (x) s ′ 1 (x) w(x) dx = 1 0 s 1 (x) ∂u ′ 1 ∂y - ∂u ′ 4 ∂y (x, 0) w(x) dx = 1 0 ∂(p 1 -q 4 ) ∂y (x, 0) v ′ 1 (x) dx + 1 0 ∂p 1 ∂x (0, y) v ′ 2 (y) dy + 0 -1 ∂q 4 ∂x (0, y) v ′ 4 (y) dy J ′ 2 = 1 0 s 2 (y) s ′ 2 (y) w(y) dy = 1 0 s 2 (y) ∂u ′ 1 ∂x - ∂u ′ 2 ∂x (0, y) w(y) dy = 1 0 ∂q 1 ∂y (x, 0) v ′ 1 (x) dx + 1 0 ∂(q 1 -p 2 ) ∂x (0, y) v ′ 2 (y) dy + 0 -1 ∂p 2 ∂y (x, 0) v ′ 3 (x) dx J ′ 3 = 0 -1 s 3 (x) s ′ 3 (x) w(x) dx = 0 -1 s 3 (x) ∂u ′ 2 ∂y - ∂u ′ 3 ∂y (x, 0) w(x) dx = - 1 0 ∂q 2 ∂x (0, y) v ′ 2 (y) dy + 0 -1 ∂(q 2 -p 3 ) ∂y (x, 0) v ′ 3 (x) dx - 0 -1 ∂p 3 ∂x (0, y) v ′ 4 (y) dy J ′ 4 = 0 -1 s 4 (y) s ′ 4 (y) w(y) dy = 0 -1 s 4 (y) ∂u ′ 4 ∂x - ∂u ′ 3 ∂x (0, y) w(y) dy = - 1 0 ∂p 4 ∂y (x, 0) v ′ 1 (x) dx - 0 -1 ∂q 3 ∂y (x, 0) v ′ 3 (x) dx + 0 -1 ∂(p 4 -q 3 ) ∂x (0, y) v ′ 4 (y) dy (31)
These formulas are of the form:

J ′ i = 4 j=1 γ j G i,j v ′ j dγ j (i = 1, ..., 4) (32)
in which the kernels, {G i,j }, are partial gradients given in terms of the partial derivatives of the eight adjoint states {p i , q i } (i=1,...,4) .

DISCRETIZATION

For purpose of numerical treatment, we assume that each sub-problem is discretized by standard centered finite-differences over a uniform (sub-)mesh of dimension N X × N Y rectangular cells. This permits a fast direct inversion by discrete separation of variables:

u h = (Ω X ⊗ Ω Y ) (Λ X ⊕ Λ Y ) -1 (Ω X ⊗ Ω Y ) f h (33)
Ω X and Ω Y are respectively the N X ×N X and N Y ×N Y orthogonal matrices associated with the discrete sine transform. For Dirichlet boundary conditions these matrices are also symmetric.

The matrices Λ X and Λ Y are diagonal matrices of the known eigenvalues of the second-order difference operators in x and y respectively (see [START_REF] Désidéri | Application of MGDA to domain partitioning[END_REF] for details).

For each sub-domain Ω i , derivatives normal to a given interface are calculated by one-sided second-order finite differences, and tangential derivatives, by central differencing.

The matching-defect integrals are then approximated by the trapezoidal rule:

J 1 . = h Z 2 N X -1 j=1 s 1 2 ,j w j , J 2 . = h Y 2 N Y -1 k=1 s 2 2 ,k w k , J 3 . = h X 2 N X -1 j=1 s 3 2 ,j w j , J 4 . = h Y 2 N Y -1 k=1 s 4 2 ,k w k .
(34) Two adjoint problems are solved on each sub-domain Ω i , again by direct inversions, to get the functions p i and q i whose derivatives are then approximated along interfaces and injected in (31). Then [START_REF] Désidéri | Application of MGDA to domain partitioning[END_REF]:

J ′ i . = 4 j=1 N Z k=1 G i,j .k v j ′ ,k h Z (i = 1, ..., 4) (35)
The discrete gradient of the criterion J i w.r.t. the nodal values of the control v j is given by:

∂J i ∂v j ,k . = G i,j .k h Z (36)
For each criterion J i , four such discrete gradients are calculated (one per control v j ), except that one of them is equal to 0. These four vectors are assembled in one, thereafter denoted ∇J i , of dimension 2(N X + N Y -2). Now, knowing (second-order approximations of) the criteria {J i } (i=1,...4) and their gradients {∇J i } (i=1,...,4) w.r.t. the 2(N X + N Y -2) nodal controls, we need to set up a strategy to iterate on these controls to satisfy the matching conditions at convergence.

GRADIENT-BASED COORDINATION ITERATIONS

Our main objective is to compare the standard steepest-descent method with the Multiple-Gradient Descent Algorithm (MGDA ) as potential iterative methods to satisfy the matching conditions by driving the defect functionals to 0. are known, one considers the convex hull U of these vectors, and identifies its minimum-norm element ω:

       ω = Argmin u∈U u 2 U = u ∈ R N / u = 4 i=1 α i u i ; α i ≥ 0 (∀i) ; 4 i=1 α i = 1 (44)
In our problem, the dimension N is the number of nodal controls:

N = 2(N X + N Y -2) (45)
A special parameterization of the convex hull is proposed to facilitate the determination of the element ω by a numerical optimization procedure (see Appendix A).

Then, once the element ω is determined, if ω = 01 , the current iterate is, or is treated as Pareto stationary. But here, the Pareto front is made of only one point corresponding to J i = 0 for all i. This situation corresponds to full convergence of the coordination algorithm. Otherwise (ω = 0), -ω is a descent direction for all criteria simultaneously. Thus, (39) is replaced by:

v (ℓ+1) = v (ℓ) -ρ ℓ ω (ℓ) (46)
Here again, we propose to adjust the step-size ρ ℓ according to (42). However here, it is not clear a priori that the proper scaling corresponds to ε ∼ 1.

NUMERICAL EXPERIMENTATION

Test-case

Let a and b be two adjustable constants, and:

                   ψ = (x + a) 2 + (y + b) 2 (-1 ≤ x ≤ 1) (-1 ≤ y ≤ 1) τ = 1 ln(a 2 + b 2 ) φ = τ ln ψ f (x) = 1 -x 2 f (y) = 1 -y 2 u e = f (x) f (y) φ (47)
As a result, φ is a harmonic function:

∆φ = 0 (48)
and this permits us to simplify somewhat the expression of the Laplacian of u e :

∆u e = ∆ f (x) f (y) φ + 2∇ f (x) f (y) .∇φ + f (x) f (y) ∆φ = -2 f (x) f (y) φ + 2∇ f (x) f (y) .∇φ (49) = -f e (50)
where

f e = 2 f (x) f (y) φ + 8τ ψ x(x + a)f (y) + y(y + b)f (x) (51)
Hence, for f = f e , the exact solution of the continuous problem is u = u e . The constants a and b have been introduced to destroy the symmetry in the solution. More specifically, the following settings were made: a = 5 4 and b = 3 4 . The corresponding problem has been discretized and solved using either one domain to establish a reference, or four to experiment multi-criterion optimization algorithms.

The single-domain discrete solution u h is depicted in Fig. 2 as a surface in 3D, and the corresponding contour map is given more precisely in Fig. 3.

Quasi-Newton steepest descent

In a first series of numerical experiments, the steepest-descent method was applied to drive to 0 the global criterion J. After a few trials, it appeared that best convergence was achieved by setting ε to 1 in (42), which corresponds to the quasi-Newton method.

Two experiments are reported presently. They differ in the setting of the initial interface conditions. In the first, the v i 's are initially assigned the restriction to the corresponding interface of the exact solution u e , which differs from the discrete solution by truncation errors. In this case, iterative errors are initially very small, which permits the asymptotic convergence to be assessed. In the second experiment, the controls are initially set to 0 in order to assess the global convergence.

Asymptotic convergence. The convergence history of the global criterion as well as its individual parts, {J i } (i=1,...,4) is represented in Fig. 4. The criterion J, in 20 iterations, goes from 2.9 × 10 -2 to a level below 10 -4 . Note that J 4 is somewhat smaller in magnitude and more subject to oscillations.

Global convergence.

The convergence history of the global criterion as well as its individual parts, {J i } (i=1,...,4) is represented in Fig. 5. The criterion J, in 200 iterations, is reduced by 8 orders of magnitude. The different criteria, apart from small oscillations, converge at essentially the same rate. In a linear convergence process, this rate is imposed by the most persistent mode, present in all criteria when the initial condition is arbitrary.

Convergence of the gradients. The evolution of the four gradients {∂J/∂v i } (i=1,...,4) over 200 iterations is given on Fig. 6-Fig. 9. They appear as high-frequency modes. Each one vanishes asymptotically over the corresponding interface. Discrete solution. The four-domain discrete solution is found perfectly smooth, in fact even smoother than the above single-domain discrete solution (Figs. 10 and11). This is due to a higher degree of iterative convergence.

Basic MGDA

Practical determination of the minimum-norm element ω. In the experiments of this section, at each iteration, the 3 parameters c 1 , c 2 and c 3 of (76) have been discretized uniformly by step of 0.01, and ω was set equal to the vector of minimum norm among the essentially 10 6 corresponding candidates. The resulting vector was therefore a very coarse approximation of the actual vector ω. Asymptotic convergence For this experiment, the discretized continuous solution is again imposed initially at the interfaces. The convergence history of the above basic algorithm is indicated in Fig. 12 for 20 iterations. After some initial adjustment, the trend is towards decaying, but at a deceiving slow rate. In an attempt to explain this poor convergence, the following observation was made: suppose the gradients of the individual criteria, {∂J i /∂v} (i=1,...,4) , are very different in magnitude. Remember that in a linear iterative process, unless initial conditions are very special, all quantities converge at the same rate, say Cρ (iter) , where ρ is the spectral radius, and C a constant which depends on the quantity considered. For example, in the previous experiment, J 4 itself was observed to be somewhat smaller than the other criteria, and so was its gradient. Then, the convex-hull minimum-norm element ω is paradoxically dominated by the gradient of smallest magnitude, since in the convex combination of the gradients, putting the largest weight on the smallest has the effect of reducing the norm of the combination. But this is not efficient, since this gradient corresponds to the already small criterion for which minimization is the least necessary. This observation has led us to calculate the direction ω as the minimum-norm element in the convex hull of normalized gradients. Here, the normalization was made by scaling each gradient to the corresponding value of the individual criterion. In other words, logarithmic gradients were considered. In fact, that is exactly what Newton's method does with the global criterion J.

The above experiment was then repeated, using logarithmic gradients to determine the vector ω. The corresponding convergence history is indicated in Fig. 13. This new result is now found very similar to the analogous result previously achieved by the quasi-Newton method (Fig. 4). The importance of scaling is therefore confirmed.

Global convergence. All interface controls v i 's are initially set to 0. The resulting convergence history over 200 iterations is indicated in Fig. 14 for the basic algorithm, and in Fig. 15 for the scaled algorithm based on logarithmic gradients. The first algorithm seems to convergence, but at a very slow rate. The second seems to subject to accidents and to experience a difficulty to enter the asymptotic convergence phase. The criteria stagnate. This may be caused by many factors on which current investigation is focused:

• the insufficiently accurate determination of ω;

• the non-optimality of the scaling of gradients;

• the non-optimality of the step-size, the parameter ε in (42) being maintained equal to 1 throughout;

• the large dimension of the design space, here 76 (4 interfaces associated with 19 d.o.f.'s).

MGDA-II

Recently, a variant, MGDA-II , has been proposed in which the descent direction is calculated by a direct procedure, which provides a valuable simplification of implementation, as well as more accurate definition of the direction of search [START_REF] Désidéri | MGDA II: A direct method for calculating a descent direction common to several criteria[END_REF]. This new algorithm is now presented again along with a new variant (MGDA-II b), and tested on the DDM problem.

Basic definition, scaling. In MGDA-II , the possibility to prescribe scales for the gradients, {S i } (i=1,...,n) (S i > 0 (∀i)) is offered. In the following experiments, at a given iteration, these scales are either set all equal to 1 ("no scaling prescribed"), or equal to the current values of the criteria ("prescribed scaling"):

S i = J i (i = 1, ..., 4) ; (52)
the latter implies that the descent direction is based on logarithmic gradients. These scales being supplied by the user, we define the following "scaled gradients":

J ′ i = ∇J i S i (53)
Automatic rescaling: MGDA-II b. The examination of the case n = 2 has led us to propose a slightly different handling of the scales. Here, one lets

A i = S i - i-1 k=1 c i,k (64)
only when this number is strictly-positive. Otherwise, the scale S i is redefined ("automatic rescaling") as follows:

S i = i-1 k=1 c i,k (65)
and one sets A i = ε i S i , for some small ε i . The result in (63) is still valid, and it now provides an information on gradients that have been weighted as prescribed whenever S i > i-1 k=1 c i,k , and otherwise by the procedure itself. This rescaling procedure is certainly perfectible.

Convergence experiments and discussion. MGDA-II has been tested on the partitioning problem in the four possible options corresponding to "no scaling prescribed" vs "prescribed scaling", and "automatic rescale off" vs "on". In MGDA-II b, when k<i c i,k was found greater or equal the prescribed S i (=1 or J i ), ε i was set to 0.01 (and maintained to 0 otherwise).

A first observation was made: the new procedure for determining ω is much faster, and MGDA-II seems to be less sensitive to round-off errors.

In Fig. 16 and Fig. 17, the automatic rescale is off, and the effect of scaling alone is evaluated. Over the first 200 iterations, the result is about the same. However the scaled version indicates a trend to convergence acceleration to be confirmed.

In Fig. 18 and Fig. 19, the automatic rescale is on, and the option of prescribed scaling is off/on. Again a better convergence is achieved when scales are prescribed.

In order to confirm these results, the best option "prescribed scales and automatic rescale" is compared with the basic method in Fig. 20 and Fig. 21 over 500 iterations. The trends indicate a linear convergence for the first method, and a seemingly-quadratic convergence for the second. Compared to the quasi-Newton method of Fig. 5, MGDA-II b is grossly-speaking twice slower, but it indicates a more definite trend to asymptotic convergence acceleration.

One last remark: in these experiments, we observe that scaling has the effect of making the convergence curves associated with the different criteria closer to one another. • the large dimension of the design space, here 76 (4 interfaces associated with 19 d.o.f.'s), was probably a handicap;

• a robust procedure to define the step-size should be devised; in our experiments, the parameter ε was not optimized but maintained equal to 1 throughout;

• the determination of ω in the basic method should be made more accurately by iterative refinement;

• the scaling of gradients was found important; alternatives to the logarithmic gradients should be analyzed and rationalized; more generally, preconditioning remains an open question;

• the MGDA-II variant was found faster and more robust;

• at present, our most sophisticated algorithm, MGDA-II b, also involves an automatic rescaling procedure; it indicates a definite trend to asymptotic convergence acceleration (quadratic convergence).

Figure 1 :

 1 Figure 1: Partition of a square in sub-domains {Ω i } (i=1,...,4) to solve the Dirichlet problem

  29) Sub-domain Ω 4 : on γ 4 : n = n 43 = -ı; on γ 1 , n = n 41 = +  . Thus (25)-(26) write:

  30) Then, (27)-(30) are injected in (18) and (17) to get:

- 1 - 1 -Figure 2 :Figure 3 :

 1123 Figure 2: Single-domain discrete solution u h

Figure 4 :Figure 5 :Figure 6 : 1 -Figure 7 : 2 -Figure 8 : 3 -Figure 9 :

 456172839 Figure 4: Quasi-Newton steepest descent -convergence history of criteria (discretized continuous solution imposed initially at interfaces)

Figure 11 :

 11 Figure 10: Four-domain discrete solution u h

Figure 13 :

 13 Figure 12: Basic MGDA -asymptotic convergence history of criteria

Figure 15 :

 15 Figure 14: Basic MGDA : global convergence

Figure 16 :Figure 17 :

 1617 Figure 16: MGDA-II , no scaling prescribed

In the numerical implementation, this condition is relaxed to be: ||ω|| < T OL, for a given tolerance T OL.

Conventional steepest-descent method

In the conventional approach, one considers a single matching defect measure, treating all interfaces as one:

The discrete gradient is then simply equal to the sum of the individual contributions of the interfaces:

The above global criterion can then be driven to 0 by the classical steepest-descent method [START_REF]WikipediA: The Free Encyclopedia[END_REF] [7]: at iteration ℓ, the control v is updated proportionally to (the opposite of) the discrete gradient:

for some appropriate positive step-size ρ ℓ (see below), and a new compound solution u (ℓ+1) is calculated, the defect-functional and its gradient reevaluated, and so on until a satisfactory convergence is achieved. Strictly speaking, in the standard steepest-descent method, once the direction of search is identified, by the calculation of the gradient ∇J (ℓ) , the step-size ρ ℓ is often defined via a onedimensional minimization:

This minimization is usually carried out by a numerical procedure. However here, we know of an additional information: the targeted value of J is known: J = 0. An estimation of the variation of J is given by the differential:

Hence, the step-size expected to diminish J (ℓ) of the amount δJ = -εJ (ℓ) is estimated to be:

In particular for ε = 1, we get the quasi-Newton method since the employed discrete gradient is only approximately equal to the gradient of the discrete J.

Multiple-gradient descent algorithm (MGDA )

In this subsection, we propose an alternative coordination algorithm in which the matching of the sub-solutions is treated as a multi-objective optimization problem, considering that all defect-functionals J i 's should be driven to 0 concurrently.

In the Multiple Gradient Descent Algorithm MGDA (see [START_REF] Désidéri | Multiple-Gradient Descent Algorithm (MGDA)[END_REF] for a detailed definition and convergence proof), once the individual discrete gradients,

Assuming that the gradients form a linearly-independent family, an assumption never contradicted in the numerical experiments, a family of orthogonal, but usually not orthonormal vectors are formed {u i } (i=1,...,n) according to the following:

and, for i = 2, 3, ..., n:

where:

and

for some arbitrary, but small ε i (0

The minimum-norm element in the convex of the family {u i } (i=1,...,n) is given by:

and one finds:

which confirms that ω does belong to the interior of the convex hull, so that:

and:

Convening that ε i = 0 in the regular case ( k<i c i,k = 1), and otherwise by modifying slightly the definition of the scaling factor to be

and redefining the "scaled gradient" accordingly (J ′ i = ∇J i /S ′ i ), the following holds:

that is, the same positive constant [START_REF] Désidéri | MGDA II: A direct method for calculating a descent direction common to several criteria[END_REF].

As a result of this direct, fast and accurate construction, the vector ω is usually different from the former definition, except in special cases, as for example, when n = 2 and the angle between the two gradients is obtuse. Nevertheless, the new ω also provides a descent direction common to all criteria, scaled essentially as initially prescribed.

CONCLUSION

In this study, various versions of the Multiple-Gradient Descent Algorithm (MGDA ) have been tested numerically over a domain-partitioning problem treated as a multi-objective problem in which matching defect integrals at the different interfaces are to be minimized concurrently, and in fact, all driven to 0.

The major objective of this experimentation was to assess the potential of MGDA to handle multi-objective problems in which the finite-dimensional setting was the result of discretization, thus approaching a more general functional setting. In this respect, the demonstration was made. Indeed convergence was achieved by MGDA . However the quasi-Newton method applied to the agglomerated criterion was globally found more efficient, but in the most sophisticated version (MGDA-II b), the algorithm seems to demonstrate a promising asymptotically-quadratic convergence.

Thus, if the convergence was not always found satisfactory, several observations should temper this somewhat deceiving conclusion, and many promising directions of improvement can be envisaged:

• in the problem under study, the Pareto set was reduced to the single point corresponding to all criteria equal to 0 associated with the unique solution of the discretized Poisson problem; this situation is really atypical of standard multi-objective problems; additionally, the criteria to be minimized were not really antagonistic, since they all converged at almost the same rate with the quasi-Newton method, leaving little possibility of improvement from the start; for these two reasons MGDA has been tested in a very straining situation for which it was not devised originally;

A GENERALITIES ABOUT CONVEX HULLS AND PARAMETERIZATION Notations

Parameterization -The convex hull may parameterized by identifying the set of allowable coefficients {α i } (i=1,...,n) . To satisfy the positivity condition automatically, one lets:

Then, the equality constraint

where S n is the unit sphere of R n , and precisely not R N . This sphere is easily parameterized using trigonometric functions of n -1 independent arcs φ 1 , φ 2 , ..., φ n-1 :

that is:

with φ 0 = π 2 . It is sufficient to consider the portion of the sphere corresponding to φ i ∈ [0, π 2 ] for all i ≥ 1 since the sign of the σ i 's makes no difference.

The usage of trigonometric functions is not really necessary, since one can let:

and get:

that is:

with c 0 = 0, and c i ∈ [0, 1] for all i ≥ 1.

In this way, the constraints on the coefficients {α i } have been replaced by the bounds 0 and 1 on the new parameters {c i }, independently of one another. However, the criterion to be minimized,

is now a polynomial of possibly large degree, namely 2(n -1), of the new parameters {c i }.

In the particular case of the coordination of 4 sub-domains by MGDA , and independently of the degree of refinement of the spatial discretization controlled by the integers N X and N Y , n = 4, and once the 10 scalar products (u i , u j ) (i, j = 1, ..., 4) calculated, the determination of the minimum-norm element ω is equivalent to minimizing a 6th-degree polynomial of (c 1 , c 2 , c 3 ) in the R 3 unit cube (limits included):


Dang Van Nguyen

Denis Grebenkov

Cyril Poupon

Denis Le Bihan

Jing-Rebecca Li

Dang Van Nguyen
email: vdnguyen@cmap.polytechnique.fr

Effective diffusion tensor computed by homogenization

Keywords: Effective diffusion tensor, apparent diffusion tensor, two-compartment model, Bloch-Torrey equation, Laplace equation

The convergence of the long-time apparent diffusion tensor under realistic DMR scanning conditions to the effective diffusion tensor obtained by homogenization theory was considered for a two-compartment model in both isotropic and anisotropic diffusion and general shapes. This gives a computationally efficient approach to get the homogenized diffusion tensor.

Introduction

Diffusion magnetic resonance imaging (DMRI) can give useful information on cellular structure and structural changes [1]. We show that the effective diffusion tensor obtained by mathematical homogenization theory [START_REF] Bensoussan | Asymptotic analysis for periodic structures[END_REF][START_REF] Cheng | [END_REF] is a good approximation to the long time apparent diffusion tensor under realistic DMR scanning conditions of a two-compartment model for both isotropic and anisotropic diffusion and general shapes. The homogenized diffusion tensor is obtained by solving d steady-state Laplace equations (where

3 , 2  d
is the problem dimension), which is a more computationally efficient approach than long time simulation in the time domain, either via Monte-Carlo simulation or numerical solution of the timedependent Bloch-Torrey PDE.

Theory

We consider two compartments, in  and ex  , to be an ensemble of cells and the extracellular space, respectively, with the same intrinsic diffusion coefficient D. The cell membrane is modeled by an infinitely thin permeable interface  characterized by permeability  . Fig. 1 illustrates a 2D sample in which the union of blue cells is considered as the interior compartment and the remaining part is the exterior one.

For a given diffusion gradient with the temporal profile

) (t f
and gradient strength

 / : q g   
(here  is the gyromagnetic ratio), the transverse magnetization is described by the two-compartment Bloch-Torrey partial differential equation (PDE) [4]  

) ( ) ( ) ( ) ( ) ( t , x M D t , x M x g t f i γ t t , x M         with the interface conditions at the inter-compartment interface  , Γ, x , t , x M t , x M κ x n t , x M D x n t , x M D in ex ex ex in in          )) ( ) ( ( ) ( ) ( ) ( ) (  
                                                  t k k b x k a x k t k k b x a x ds s f l e g ( i x t , x M x t , x M ds s f l e g ( i t , x M t , x M k k k k k k k k 0 0 ) ( ) exp ) ( ) ( ) ( ) exp ) ( ) (   where d k a b l k k k .. 1 ,    and ) ,.., ( 1 d x x x 
is the spatial position.

For the uniform initial condition, the DMRI signal is

x d t x M t q     ) , ( ) , ( 
. From the Taylor expansion of the DMRI signal in  , we define the apparent diffusion tensor D A which may depend on time t:

  4 0 2 0 ) ( ) , ( ln q O ds s f du q D q t q t u A T                  
The long time limit of D A can be approximated by the effective diffusion tensor obtained from homogenization theory [START_REF] Li | and analytical models of the long time apparent diffusion tensor[END_REF],

  d k d j eff jk eff D D .. 1 , .. 1    , where     Ω k j eff jk x d e v D D  
, k e  is the unit vector in the th k coordinate direction and the unknown functions v j can be found by solving the Laplace equation with the interface conditions at the inter-compartment interface  ,

  0 ) (     x v D j
Γ x , x v x v κ x n x v D x n x v D in j ex j ex ex j in in j          )) ( ) ( ( ) ( ) ( ) ( ) (  
, and the boundary conditions at the exterior boundaries of the computational box  are

    k k b x j a x j l e g x v x v k k k k ) ( ) ( ) (      , k k k k b x k j a x k j x x v x x v                        ) ( ) (

Results

We performed simulations on the computational domain   d m m, -  10 10 containing numerous convex cells . To construct these samples, we randomly generated a set of segments in 2D (Fig. 3) and a set of faces in 3D (Fig. 5). Since the boundaries should be periodic, we mirrored this set across all coordinate axes. Then, each segment (face) was thickened to create an extracellular compartment ex  . The interior compartment in  includes polygonal cells. The same intrinsic diffusion coefficient of water molecules

                    3 ) ( 2 2 2 0 2 0    g ds s f du q b t u T   .
Then, the entries of A D appear in front of the first-degree term of the polynomial fit of ) , ( ln t q    versus b-value. The Laplace equation was also solved over the same box to obtain eff D . All computations were carried out on a Dell Latitiude E6410 laptop computer.

2D results

The 2D sample was created by 20 thick segments that divide the computational domain into 116 convex cells (Fig. 3). The volume fraction i v is 0,77 , the average radius of cells is Results for two gradient directions with  =10ms, 50ms and 190ms are summarized in the following table:

10 -4 m/s 10 -5 m/s q D q A T   q D q eff T   q D q A T   q D q eff T   q q   / =10 =110 =190 =10 =110 =190 (

3D results

The 3D sample was created by 16 thick faces which cut the cube into 80 cells (Fig. 5). The average radius of cells is Results for three gradient directions with  =10ms, 50ms and 90ms are summarized in the following table: 10 -4 m/s 10 -5 m/s

q D q A T   q D q eff T   q D q A T   q D q eff T   q q   / =10 =50 =90 =10 =50 =90 (

Conclusions

The long-time apparent diffusion tensor D A approaches the steady-state tensor D eff computed by the homogenization in both cases of isotropic and anisotropic diffusion and the convergence is faster at higher permeability. This convergence seems to be linear in term of 1   that agrees with the result for 1D periodic structure in the long-time regime [START_REF] Dmitriy | [END_REF]. The computation of D eff is very much faster than that of D A . Since the samples were generated randomly, this conclusion is expected to be valid for general shapes.



  are normal vectors at  pointing outward ex  and in  , respectively), and the boundary conditions at the exterior boundaries of the computational box representative sample of the cellular structure:

Fig. 1 :

 1 Fig. 1: A two-compartment model with the extracellular

  both interior and exterior compartments. Two values of permeability, Torrey equation was solved for Stejskal -Tanner PGSE sequences (Fig.2) with two rectangular gradient pulses of duration

Fig. 2 :

 2 Fig.2: Stejskal-Tanner PGSE

  1,0) 5,33•10 -4 5,24•10 -4 5,23•10 -4 5,22•10 -4 3,84•10 -4 3,66•10 -4 3,65•10 -4 3,64•10 -4 (0,1) 5,69•10 -4 5,61•10 -4 5,60•10 -4 5,59•10 -4 4,06• 10 -4 3,95•10 -4 3,94•10 -4 3,93•10 -4

.

  performed for three gradient directions [1,0,0], [0,1,0] and [0,0,1]. Each A D was computed at 9 b-values changing from 0 The computation took between 30 minutes and 2 hours to compute each A D . On the contrary, it took less than one minute to compute eff D for the mesh size of 50476 vertices.

Fig. 5 FiveD A was showed in term of 1 

 51 Fig.5: 3D sample

  1,0,0) 1,24•10 -3 1,22•10 -3 1,22•10 -3 1,21•10 -3 9,65•10 -4 8.89•10 -4 8.78•10 -4 8,66•10 -4 (0,1,0) 1,73•10 -3 1,67•10 -3 1,67•10 -3 1,66•10 -3 1,26• 10 -3 1,11•10 -3 1,09•10 -3 1,07•10 -3 (0,0,1) 6,89•10 -4 6,87•10 -4 6,85•10 -4 6.85•10 -4 5,01•10 -4 4,93•10 -4 4,92•10 -4 4,91•10 -4

Fig. 6 :

 6 Fig. 6: The convergence of D A to D eff in 3D for the direction [1,0,0] and