Houssaine El

Quenjel
email: quenjel@unice.fr

Analysis of accurate and stable nonlinear finite volume scheme for anisotropic diffusion equations with drift on simplicial meshes

This work addresses the development and analysis of a second-order accurate finite volume scheme for parabolic equations with anisotropy on general simplicial meshes. The discretization involves only vertex unknowns without processing additional ones. The scheme construction makes use of a nonlinear transformation of the linear elliptic term. Two propositions are mainly presented for the approximation of the mobility function at the interfaces. The existence of positive solutions for the discrete system is guaranteed thanks to the proved a priori estimates. The energy dissipation of the scheme is moreover ensured. The convergence of the approach is established. Numerical tests are given to show the efficiency, accuracy and robustness of the proposed approach, with respect to the anisotropy, while a particular emphasis is set on the effects of the approximate mobility. They also confirm the obtained theoretical results, especially the decay of the free energy when time grows.

Introduction

Anisotropic convection-diffusion equations are essential in modeling many real-world problems arising from complex reservoir models [START_REF] Chavent | Mathematical models and finite elements for reservoir simulation: single phase, multiphase and multicomponent flows through porous media[END_REF], chemotaxis processes [START_REF] Horstmann | From 1970 until present: the Keller-Segel model in chemotaxis and its consequences[END_REF] and systems of semiconductor devices [START_REF] Bank | Numerical methods for semiconductor device simulation[END_REF]. In such an application the numerical approximation of the model should ideally preserve the physical properties of the continuous problem. For instance, in the context of porous media flows the permeability can be strongly anisotropic and highly heterogeneous. Then, in these situations, the approximate solution accounting for the saturation or the concentration must be nonnegative without spurious oscillations in order to respect some standard laws of physics (e.g. Fick's law). This holds the name of a discrete maximum principle. Moreover, in the numerical modeling of semiconductors, the control of some dissipated entropies is very useful to understand the large time behavior of the physics being modeled. Hence, developing consistent, accurate and cheap numerical methods of type finite volumes on general meshes respecting the underlined properties is of great importance.

In view of its fundamental role, extensive research has been devoted to the discrete maximum principle for anisotropic linear diffusion problems. The pioneer two-point flux approximation (TPFA) scheme fulfills this property by construction [START_REF] Eymard | Finite volume methods[END_REF]. In the presence of convective effects, it is recommended to employ the Scharfetter-Gummel approach of the TPFA scheme [START_REF] Chainais-Hillairet | Finite-volume schemes for noncoercive elliptic problems with Neumann boundary conditions[END_REF][START_REF] Scharfetter | Large-signal analysis of a silicon read diode oscillator[END_REF] instead of the upwind version for better accuracy and stability. Unfortunately, these methods necessitates a stringent orthogonality condition on the mesh together with a scalar isotropic tensor. From a general viewpoint, it is well understood that no linear finite volume method for the linear diffusion preserves unconditionally the physical ranges of the approximate solution [START_REF] Droniou | Finite volume schemes for diffusion equations: introduction to and review of modern methods[END_REF][START_REF] Nordbotten | Monotonicity of control volume methods[END_REF] without additional assumptions. A tough attention was then given to nonlinear versions [START_REF] Camier | A monotone nonlinear finite volume method for approximating diffusion operators on general meshes[END_REF][START_REF] Gao | A second-order positivity-preserving finite volume scheme for diffusion equations on general meshes[END_REF][START_REF] Kapyrin | A family of monotone methods for the numerical solution of three-dimensional diffusion problems on unstructured tetrahedral meshes[END_REF][START_REF] Potier | Finite volume scheme satisfying maximum and minimum principles for anisotropic diffusion operators[END_REF][START_REF] Lipnikov | Interpolation-free monotone finite volume method for diffusion equations on polygonal meshes[END_REF][START_REF] Sheng | The finite volume scheme preserving extremum principle for diffusion equations on polygonal meshes[END_REF]. Most of these schemes make use of a convex splitting of the flux in the co-normal direction. Other unknowns, placed at the vertices or edges etc, are required to define accurately each component and they can be eliminated after that. Although they are accurate of second-order, these kinds of methods do not fulfill the coercivity of the scheme without further assumptions on the mesh and the anisotropic ratio. To ensure the latter, a general approach was addressed in [START_REF] Cancès | Monotone corrections for generic cell-centered finite volume approximations of anisotropic diffusion equations[END_REF]. It is applied to centered scheme or merely built on coercive centered finite volume methods. Yet, consistency or convergence is valid when a numerical assumption holds. Recently, a few convergent positive schemes were designed and analyzed for degenerate parabolic equations with anisotropy, see for instance [13,[START_REF] Ghilani | Positivity-preserving finite volume scheme for compressible two-phase flows in anisotropic porous media: The densities are depending on the physical pressures[END_REF][START_REF] Ghilani | Positive control volume finite element scheme for a degenerate compressible two-phase flow in anisotropic porous media[END_REF][START_REF] Quenjel | Enhanced positive vertex-centered finite volume scheme for anisotropic convection-diffusion equations[END_REF]. The methodology used in these references is accurate of first order for the diffusion since the idea is based on an upwinding strategy of Godunov type. In the work [14], the authors proposed and investigated an accurate vertex approximation gradient methodology. It uses vertex and cell unknowns where the latter are eliminated at the solver level. Also, the implementation of this methodology shows that the iterative solver is extremely sensible to highly anisotropic ratios of the permeability when the discrete mobility function is given by a centered arithmetic mean. Recall that the advantage of the methods [START_REF] Cancès | Numerical analysis of a nonlinear free-energy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations[END_REF]14] is that they can deal with more general grids and heterogeneities.

On the other hand, the considered model problem is dissipated and possesses an energy functional of Lyapunov type. The control of the latter allows to obtain relevant information on the asymptotic behavior of the solutions such as the exponential decay to the equilibrium [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 | Qualitative behavior of solutions of a degenerate nonlinear drift-diffusion model for semiconductors[END_REF]. In the discrete setting, the large time behavior has been addressed in [START_REF] Bessemoulin-Chatard | Exponential decay of a finite volume scheme to the thermal equilibrium for drift-diffusion systems[END_REF][START_REF] Bessemoulin-Chatard | A finite volume scheme for convection-diffusion equations with nonlinear diffusion derived from the Scharfetter-Gummel scheme[END_REF][START_REF] Chainais-Hillairet | Large-time behaviour of a family of finite volume schemes for boundary-driven convection-diffusion equations[END_REF][START_REF] Filbet | A finite volume scheme for boundary-driven convection-diffusion equations with relative entropy structure[END_REF] using a Two-Point finite volume discretizations together with the Scharfetter-Gummel strategy. The exponential decay to equilibrium states has been shown in [START_REF] Li | Large time behaviors of upwind schemes and B-schemes for Fokker-Planck equations on R by jump processes[END_REF] using jump processes. Establishing finite volume methods coupling between the energy diminishing, optimal accuracy, as well as positivity for such a model under general assumptions on the data and on the mesh is quite challenging. Only few papers have been devoted to such a topic in the context of multi-point framework [START_REF] Cancès | Numerical analysis of a nonlinear free-energy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations[END_REF]14].

In this paper, we propose a vertex-centered finite volume scheme on general simplicial meshes applied to a basic anisotropic linear diffusion equation with drift, due to its practical importance. The main objective of this contribution consists in extending the ideas of the scheme presented in [START_REF] Quenjel | Enhanced positive vertex-centered finite volume scheme for anisotropic convection-diffusion equations[END_REF]. Indeed, even if the discretization elaborated in [START_REF] Quenjel | Enhanced positive vertex-centered finite volume scheme for anisotropic convection-diffusion equations[END_REF] enjoys some advantages, it however suffers from several limitations. For instance, we reveal :

• the numerical convergence is only first order in space;

• it does not allow the control of the free energy;

• it is not suitable to preserve the long time behavior of the solutions since the methodology is based on the upwinding, see [START_REF] Bessemoulin-Chatard | Développement et analyse de schémas volumes finis motivés par la présentation de comportements asymptotiques[END_REF] for deep details.

In the present contribution, we aim to dispense with the previous items weakening the use of the scheme of [START_REF] Quenjel | Enhanced positive vertex-centered finite volume scheme for anisotropic convection-diffusion equations[END_REF] in studying the long time behavior. We below recall and state the assets of the new scheme.

• Unknowns are only placed at the vertices, which turns out the computational cost cheaper;

• accuracy of second order;

• the control of the free energy;

• preserving the asymptotic of solutions as fast as possible (numerical exponential decay);

• handling highly anisotropic and heterogeneous domains;

• honoring the physical bounds (positivity) on the approximate solution;

• theoretical convergence towards the continuous solution of the mathematical problem.

The task is first achieved by formulating the elliptic term into a very specific nonlinear manner allowing the dissipation of the entropy function. Next, we use the spatial discretization elaborated in [START_REF] Quenjel | Enhanced positive vertex-centered finite volume scheme for anisotropic convection-diffusion equations[END_REF]. Moreover, in [START_REF] Quenjel | Enhanced positive vertex-centered finite volume scheme for anisotropic convection-diffusion equations[END_REF], we employed an upwind approximation of the mobility-like function. In the current paper, we take advantage of a centered scheme to avoid the reduction in the numerical convergence rate. Two propositions on the centered approximation can be made. First, one has the possibility to take the arithmetic mean value as done for instance in [13]. Second, one considers a fractional logarithmic average. Both of them respect the above prescribed physical properties of the solutions and allow the convergence of the numerical scheme from a theoretical point of view.

In practice, the first choice does not support highly anisotropic tensors whereas the behavior of the logarithmic centered choice is independent of anisotropy. The content of this paper is outlined as follows. We next state the mathematical model we are interested in as well as the definition of weak solutions. Section 2 describes the finite volume discretization, namely the spatial and temporal meshes, the reconstruction operators and the discrete functional spaces. We next present the proposed numerical scheme under its pointwise and concise forms. In Section 3, we establish several properties of the scheme, especially the energy estimates, estimations on the free energy of the system and existence of discrete solutions. In Section 4, we study the convergence of these solutions towards a weak solution by means of a recent compactness argument. Finally, Section 5 shows numerical assessment of errors produced by our method and their associate convergence rates. The approach robustness with respect to the tensor anisotropy is also tested. It is illustrated that the proposed methodology is free diminishing in the sense that the energy dissipates in the course of time.

Let Ω be a bounded connected open polyhedral subset of R d (d ∈ N * ) with a Lipschitz boundary Γ. Let us set Q t f = Ω × (0, t f ) where t f > 0 is the final time. In this work, we are interested in the model example

s t -div D∇s + sD∇G(x) = 0 in Q t f , (1.1)
with the no-flux boundary condition

D∇s + sD∇G(x) • n = 0 on Γ × (0, t f ), (1.2)
and the initial state of the solution

s(0, •) = s 0 (x) in Ω. (1.
3)

The approach we adapt here is based on an equivalent formulation of the potential. Therefore, the equation (1.1) rewrites

s t -div sD∇(log(s) + G(x)) = 0 in Q t f . (1.4)
Now, let us state the main assumptions on the physical data of the problem. They are mandatory to give a proper sense to the sought solutions.

(H 1 ) The matrix D ∈ L ∞ (Ω) d×d is symmetric and uniformly elliptic in the sense that there exist α 0 and α 1 such that

α 0 |u| 2 ≤ D(x)u • u ≤ α 1 |u| 2 ∀ a.e. x ∈ Ω and ∀u ∈ R d . (H 2 )
The exterior potential G is a nonnegative function of C 1 (Ω, R).

(H 3 ) The initial datum s 0 is nonnegative and belongs to L 1 (Ω). We define the convex mapping

H : R + -→ R + by H(v) = v log(v) -v + 1. (1.5)
We require that this function possesses a finite entropy condition on s 0 i.e.

Ω

H(s 0 (x)) dx < +∞, and

Ω s 0 (x) dx > 0. (1.6)
The following statement surveys the definition of a weak solution.

(i) H(s) ∈ L ∞ (0, t f ; L 1 (Ω)) and ξ(s) := √ s ∈ L 2 (0, t f ; H 1 (Ω)).
(ii) The integral identity holds

- Qt f sψ t dx dt- Ω s 0 ψ(•, 0) dx+ Q tf D ∇s+s∇G(x) •∇ψ dx dt = 0, ∀ψ ∈ C ∞ c (Ω×[0, t f )).
(

The convergence (up to a subsequence) of the proposed numerical scheme (2.10)-(2.12) of the current work yields directly the existence of a weak solution to the model (1.1)-(1.3). More details are given in Propositions 4.1-4.2. A same strategy was employed in [START_REF] Cancès | Numerical analysis of a nonlinear free-energy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations[END_REF] in the context of discrete duality finite volume method. Even though no-flux Neumann boundary conditions are imposed in our case, one can mimic the compactness approach developed in [1.7. Existence Theorem.] [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF] to provide another alternative for the existence result. In this reference, the authors treated a general parabolic model with mixed Dirichlet-Neumann boundary conditions. Note that the current assumptions on the data (e.g. (H 3 )) and the definition of the weak solution (e.g. item (i) of Definition 1.1) are different from the ones prescribed in [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF]. As highlighted in [START_REF] Cancès | Numerical analysis of a nonlinear free-energy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations[END_REF], the uniqueness question of such a weak solution is still an open problem for s 0 ∈ L 2 (Ω). In the case where H(s 0 ) ∈ L 1 (Ω), the authors conjectured that the weak solutions in the sense of Definition 1.1 are renormalized ones, which yields the uniqueness. As a consequence, one could recover the convergence of the whole sequence of solutions.

We define the free energy E and the dissipation I for the evolution equation (1.1) by Thereby, the free energy is diminishing in time. This structural feature was intensively employed to study the large time behavior of solutions thanks to entropy's methods. For more details on this strategy, we refer to the article [START_REF] Carrillo | Entropy dissipation methods for degenerate parabolic problems and generalized Sobolev inequalities[END_REF]. In the recent years, a particular attention was paid to setting the entropy principles from a numerical point of view. For instance, using TPFA scheme, the numerical investigation of asymptotic of solutions to convection-diffusion equations based on the discrete version of the entropy methodology can be found in [START_REF] Bessemoulin-Chatard | A finite volume scheme for convection-diffusion equations with nonlinear diffusion derived from the Scharfetter-Gummel scheme[END_REF][START_REF] Chainais-Hillairet | Asymptotic behaviour of a finite-volume scheme for the transient drift-diffusion model[END_REF][START_REF] Glitzky | Exponential decay of the free energy for discretized electro-reaction-diffusion systems[END_REF][START_REF] Filbet | A finite volume scheme for boundary-driven convection-diffusion equations with relative entropy structure[END_REF][START_REF] Chainais-Hillairet | Large-time behaviour of a family of finite volume schemes for boundary-driven convection-diffusion equations[END_REF]. We also mention the work [START_REF] Cancès | Numerical analysis of a nonlinear free-energy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations[END_REF] where a discrete duality finite volume scheme enjoying (1.8) was studied.

E(t) = Ω H(s) + sG dx, I(t) = Ω sD∇(log(s) + G) • ∇(log(s) + G) dx.
We will see later on that our methodology respects the decay of the free energy even on anisotropic situations.

Finite volume discretization

In this section, we set the discrete tools that are necessary to define and study the proposed numerical scheme.

Meshes and reconstruction operators

Spatial discretization : the domain Ω is discretized with two different partitions. The first one consists in an admissible finite element triangulation referred to as the primal mesh. The second one or the dual mesh is built around the vertices of the elements. While the presentation of the terminologies below can be extended to the 3D setting, we rather expose the 2D version in order to keep the ease readability of the discrete framework.

Let T be a conforming triangulation of Ω composed of triangles satisfying T ∈T T = Ω. We define x T the mass center and h T the diameter of the simplex T ∈ T . Also, let ρ T be the largest diameter of the inscribed ball in T . The size and regularity of the primal mesh are respectively given by

h T = max T ∈T h T , and Θ T = max T ∈T h T ρ T .
When h T tends to zero, we assume that the refined meshes T h of the same triangulation have a uniformly bounded regularity Θ T h [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. We denote V the set of vertices of the mesh, and V T the nodes of the simplex T ∈ T . We can label each vertex of T thanks to the permutation τ acting on V T as depicted on the left side of Figure 1. We denote T i the set of elements having in common the vertex i ∈ V. We assign a unique dual cell or control volume K i centered at i for any given node i ∈ V. A simple way to construct K i is to connect the center of any triangle T ∈ T i with the midpoint of edges of T sharing the same vertex i. Given two finite volumes

K i and K τ (i) we set σ T iτ (i) = K i ∩ K τ (i) ∩ T . The notation σ T iτ (i) represents the (d -1)-dimensional measure of σ T iτ (i) . The unit normal on this interface oriented from K i to K τ (i) is denoted by n σ T iτ (i)
. The set E T K i refers to all the edges of K i contained in the triangle T ∈ T i . We designate by x i the center of the control volume K i , and by |K i | the d-dimensional measure of K i . The set of the control volumes K i will be referred to as M. Note that x i denotes also the coordinate of the node i.

Remark 2.1. In case of d ≥ 3 one always can label the faces sharing a vertex i ∈ V T for some element T thanks to a suitable permutation τ . A possible way to do so is given as follows. First, we fix the combination such that τ d+1 (i) = i. Then we apply τ d-times on iτ * (i) to obtain all faces sharing the node i.

τ * (i) = τ (i)τ • τ (i) • • • τ d (i), τ d = τ • • • • • τ d times ,
In this case, we would use the notation σ T iτ * (i) instead of σ T iτ (i) for the dual interfaces. Temporal discretization : the time interval (0, t f ) is discretized thanks to a uniform subdivision given by the increasing sequence (t k ) k=0,••• ,N such that t 0 = 0 and t N = t f . Then, the time step is fixed to t = t f /N . This uniform stepping is only taken to avoid heavy notations. The proofs carried out in this paper can be broadened to nonuniform time steps without any technical issues.

Reconstruction operators : first, let us define X T the standard P 1 finite elements space on T . We consider (ϕ i ) i∈V the classical basis of X T satisfying ϕ i (x j ) = 1 if j = i and ϕ i (x j ) = 0 else. We define the reconstruction operator π T : R #V -→ X T by setting

π T u = i∈V u i ϕ i , ∀u ∈ R #V .
Similarly, let X M be the space of piecewise constant functions on the dual mesh M. We define the piecewise constant operator π M : R #V -→ X M such that

π M u = i∈V u i 1 K i , ∀u ∈ R #V .
If g is a nonlinear function from R into R. We adopt the convention making g(u) = (g(u i )) i∈V for every u ∈ R #V . Then, one has

π T g(u) = i∈V g(u i )ϕ i , and π M g(u) = i∈V g(u i )1 K i , ∀u ∈ R #V . Let {π T u k } k=0,••• ,N be a sequence of functions of X T and let {π M u k } k=0,••• ,N
be a family of functions of the trial space X M . We define the time-dependent reconstruction functions π T , t , π M, t respectively by

π T , t u(•, t) = π T u k , and π M, t u(•, t) = π M u k , ∀t ∈ (t k-1 , t k ], ∀k = 1, • • • , N. (2.1)
This conduct us to define the time-dependent discrete gradient ∇π T , t u by

∇π T , t u(•, t) = ∇π T u k , ∀t ∈ (t k-1 , t k ], ∀k = 1, • • • , N.
The following result establishes the link between some reconstruction functions and the discrete gradient.

Lemma 2.1. (see [START_REF] Quenjel | Enhanced positive vertex-centered finite volume scheme for anisotropic convection-diffusion equations[END_REF]Lemma 2.4]) Let π T u be in X T . We consider the piecewise constant functions π T u and π T u defined respectively by

π T u(x) = max y∈T π T u(y) and π T u(x) = min y∈T π T u(y), ∀x ∈ T, ∀T ∈ T .
Then, one gets

π T u -π T u L 2 (Ω) ≤ #V T h T ∇π T u L 2 (Ω) d . (2.2) As a consequence π M u -π T u L 2 (Ω) ≤ #V T h T ∇π T u L 2 (Ω) d . (2.3)

Description of the numerical scheme

For the sake of simplicity we set Υ(s) = log(s) + G throughout the paper. Without loss of generality, we assume that the tensor D is constant per triangles. Before presenting the proposed numerical scheme, we would like to survey how we derive the main discrete terms. Let us select a control volume K i for some i ∈ V. We integrate on K i with ∂K i = T ∈T i σ T iτ (i) and make use of Green's formula to get

- ∂K i sD∇Υ(s) • n dσ = - T ∈T i σ T iτ (i) ∈E T K i σ T iτ (i) s∇Υ(s) • Dn σ T iτ (i) dσ. (2.4)
Being a finite volume method, our construction is based on the discretization of the local continuous fluxes. Following [START_REF] Quenjel | Enhanced positive vertex-centered finite volume scheme for anisotropic convection-diffusion equations[END_REF] we propose the following centered scheme to the integral flux across the interface σ T iτ (i)

- σ T iτ (i) s∇Υ(s) • Dn σ T iτ (i) dσ ≈ s iτ (i) G T iτ (i) (s T ), (2.5)
where s iτ (i) is given by combining the two different means

s c iτ (i) = s i + s τ (i) 2 , s iτ (i) =        2 √ s i - √ s τ (i) log(s i ) -log(s τ (i) ) 2 , if s i = s τ (i) , s i , otherwise , (2.6)
in the sense that

s iτ (i) = (1 -λ)s iτ (i) + λs c iτ (i) , (2.7)
where λ is a parameter belonging to (0, 1]. The fact that λ > 0 is required to reinforce the existence of λ > 0 such that

s iτ (i) ≥ λ max s i , s τ (i) . (2.8)
In our case λ = λ/2. This inequality is a key ingredient to prove the existence of positive solution to the proposed finite volume scheme. From a numerical viewpoint, we can consider a very small λ in order to show the efficiency of

s iτ (i) ≈ s iτ (i) compared to s c iτ (i) . The quantity G T iτ (i) (s T ) is a numerical flux-like function that approximates √ s∇Υ(s) along the co-normal vector Dn σ T iτ (i) G T iτ (i) (s T ) = j∈V T a T ij s jτ (j) Υ(s j ) -Υ(s τ (j) ) , (2.9)
where the weight a T ij is given by [40]

a T ij = 1 |T | σ T iτ (i) σ T jτ (j) Dn σ T iτ (i) • n σ T jτ (j)
.

By construction, notice that the numerical flux-like function is conservative. Indeed, the flux from τ (i) to i in the triangle T is given by G τ (i)i (s T ). Replacing i by τ (i) in (2.9), with τ

• τ (i) = i, implies G τ (i)i (s T ) = j∈V T a T τ (i)j s jτ (j) Υ(s j ) -Υ(s τ (j) ) . Now, observe that a T τ (i)j = 1 |T | σ T iτ (i) σ T jτ (j) Dn τ (i)i • n σ T jτ (j) = - 1 |T | σ T iτ (i) σ T jτ (j) Dn τ (i)i • n σ T jτ (j) = -a T ij .
This entails that G τ (i)i (s T ) = -G iτ (i) (s T ), which ensures the discrete local flux conservation.

To avoid a CFL constraint on the time steps in case of strong anisotropic ratios of the tensor D we consider a fully implicit (backward Euler) discretization in time. Then, the numerical scheme consists in finding a sequence (s k i ) i solving the following algebraic set of equations at each time

iteration k ∈ {1, • • • , N } , s 0 i = 1 |K i | K i s 0 (x) dx, ∀i ∈ V, (2.10)
|K i | s k i -s k-1 i t + T ∈T i σ T iτ (i) ∈E T K i s k iτ (i) G T iτ (i) (s k T ) = 0, ∀i ∈ V, (2.11)
s k i > 0, ∀i ∈ V. (2.12)
The last constraint gives sense to sought discrete solutions due to the use of the log function.

Remark 2.2. The nonlinear scheme takes into account a zero-flux boundary condition. The discretization can be applied to homogeneous Neumann and (positive) Dirichlet boundary conditions for instance. For this purpose, one first needs to specify the set of Dirichlet boundary vertices V D , see [START_REF] Quenjel | Enhanced positive vertex-centered finite volume scheme for anisotropic convection-diffusion equations[END_REF]. Then, the corresponding numerical scheme (2.10)-(2.12) remains the same for any degree of freedom i ∈ V \V D and we impose the considered (positive) Dirichlet function for i ∈ V D . The scope of this paper (positivity, a priori estimates, convergence) holds true for these kinds of boundary conditions.

The above system can be rewritten under a compact form. This formulation is of chief importance since it allows to straightforwardly conduct the analysis of the resulting finite volume scheme in a rigorous setting. To this purpose, we introduce the local diagonal matrix (S T ij ) 1≤i,j≤#V T whose entries are

S T ij = s k iτ (i) if j = τ (i) 0 if j = τ (i)
.

For every u ∈ R #V T , we denote δ T u the vector of R #V T such that the

i th component is (δ T u) i = u i -u τ (i) . We next define the discrete L 2 -norm on R #V by u T , v T T = #V i=1 |K i | u i v i .
Therefore, the equivalent form of the numerical scheme is obtained by performing a discrete integration by parts in space. It reads

s k T -s k-1 T t , v T T + T ∈T S T δ T Υ(s k T ) • A T S T δ T v T = 0, ∀v T ∈ X T , ∀k = 1, • • • , N. (2.

13)

The symmetric matrix A T is defined by A T = (a T ij ) 1≤i,j≤#V T . By taking v T = 1 T in (2.13) we recover the discrete conservation of mass claiming

s k T , 1 T T = Ω π M s k dx = Ω s 0 dx, ∀k = 1, • • • , N. (2.14)
The following lemma is one of the keystone results for the convergence analysis of our scheme. It shows in particular that the matrix A T is definite-positive.

Lemma 2.2. There exist two positive constant µ 1 and µ 2 that depend only on α 0 , α 1 and the mesh regularity so that for all T ∈ T , one has

µ 1 |T | h 2 T δ T u • δ T u ≤ δ T u • A T δ T u ≤ µ 2 |T | h 2 T δ T u • δ T u, ∀u ∈ R #V T . (2.15)
Proof. The proof is carried out in [START_REF] Quenjel | Enhanced positive vertex-centered finite volume scheme for anisotropic convection-diffusion equations[END_REF]Lemma 2.3]. It extends the ideas provided in [START_REF] Brenner | Convergence of a vertex centred discretization of two-phase Darcy flows on general meshes[END_REF]14].

Remark 2.3. From a theoretical viewpoint, s k iτ (i) can take other possible mean values. Indeed, the proposed scheme and the analysis are valid for any average s k iτ (i) satisfying

s ,k iτ (i) < s k iτ (i) ≤ s c,k iτ (i) , and s k iτ (i) ≥ λ max s k i , s k τ (i) ,
for some λ > 0. An example fulfilling these equalities is given by [37]

s k iτ (i) = (s k i ) r + (s k τ (i) ) r 2 1/r , 1 6 ≤ r ≤ 1.
In the sequel, unless specified, we denote by C generic positive constants that depend only on the physical data described in (H 1 )-(H 3 ) and possibly on the mesh regularity.

A priori estimates and existence of discrete solutions

A priori estimates serve for multiple purposes. They ensure the stability of the proposed approach, namely the coercivity of the scheme and the decay of the free energy. The coercivity means that we can get a control on the approximate gradient of ξ(s). This property allows in particular to prove the existence of the finite volume scheme. Therefore, one can derive a positive lower bound on any discrete solution. Contrary to [START_REF] Cancès | Numerical analysis of a nonlinear free-energy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations[END_REF] our energy estimates, namely the next result, are not based on the energy-dissipation estimations. Proposition 3.1. Let (s k i ) i∈V be a sequence such that the finite volume scheme (2.10)-(2.12) holds. Then, there exist positive constants C 1 , C 2 , and C 3 , depending only on the physical data and the regularity of the mesh such that

T ∈T |T | i∈V T s k iτ (i) ≤ C 1 , (3.1)
N k=1 t T ∈T |T | h 2 T i∈V T s k iτ (i) log(s k i ) -log(s k τ (i) ) 2 ≤ C 2 , (3.2)
N k=1 t ∇π T ξ(s k ) L 2 (Ω) d ≤ C 3 . (3.3)
Proof. First, we observe that

s k iτ (i) ≤ s k i + s k τ (i) and |T | ≤ C |T ∩ K i | for all i ∈ V T .
Summing by dual cells and using the mass conservation property (2.14) we directly check

T ∈T |T | i∈V T s k iτ (i) ≤ 2C i∈V |K i | s k i ≤ 2C s 0 L 1 (Ω) = C 1 .
Next, take v T = log(s k T ) in the formulation (2.13) and sum over k = 1, • • • , N . This gives

E 1 + E 2 + E 3 = 0, (3.4)
where each term writes

E 1 = N k=1 s k T -s k-1 T , log(s k T ) T , E 2 = N k=1 t T ∈T S T δ T log(s k T ) • A T S T δ T log(s k T ), E 3 = N k=1 t T ∈T S T δ T G T • A T S T δ T log(s k T ).
Thanks to the convexity of the function H defined in (1.5) we obtain

H(s N T ) -H(s 0 T ), 1 T T ≤ E 1 . (3.5)
According to the Cauchy-Schwarz inequality, the fact that

G(x i ) -G(x τ (i) ) ≤ ∇G ∞ h T for all i ∈ V T , Lemma 2.
2 and the estimation (3.1) we get

|E 3 | ≤ (E 2 ) 1/2 N k=1 t T ∈T S T δ T G T • A T S T δ T G T 1/2 ≤ C(E 2 ) 1/2   N k=1 t T ∈T |T | h 2 T i∈V T s k iτ (i) G(x i ) -G(x τ (i) ) 2   1/2 ≤ C ∇G ∞ (E 2 ) 1/2   N k=1 t T ∈T |T | i∈V T s k iτ (i)   1/2 ≤ C(E 2 ) 1/2 .
By virtue of Young's inequality we infer

|E 3 | ≤ C + 1 2 E 2 . (3.6)
In addition, we combine (3.4)-(3.6) and use the finite entropy condition (1.6) to ensure the uniform bound E 2 ≤ C. Owing to (2.15), it can be checked that

N k=1 t T ∈T |T | h 2 T i∈V T s k iτ (i) log(s k i ) -log(s k τ (i) ) 2 ≤ C 2 .
Applying once more Lemma 2.2 yields

∇π T ξ(s k ) 2 L 2 (Ω) d ≤ µ 1 T ∈T |T | h 2 T i∈V T ξ(s k i ) -ξ(s k τ (i) ) 2 .
By definition of s ,k iτ (i) (2.6) we automatically obtain

∇π T ξ(s k ) 2 L 2 (Ω) d ≤ µ 1 4 T ∈T |T | h 2 T i∈V T s ,k iτ (i) log(s k i ) -log(s k τ (i) ) 2 ≤ µ 1 4 T ∈T |T | h 2 T i∈V T s k iτ (i) log(s k i ) -log(s k τ (i) ) 2 .
We therefore derive the estimation (3.3).

Let us elaborate the discrete counterpart of the energy-dissipation relationship (1.8). To this end, we define the approximate free energy by

E k h = H(s k T ), 1 T T + G T , s k T T ≥ 0, ∀k = 0, • • • , N. (3.7)
The approximate dissipation is given by

J k h = T ∈T S T δ T Υ(s k T ) • A T S T δ T Υ(s k T ). (3.8)
Proposition 3.2. The finite volume scheme is free diminishing (2.10)-(2.12) in the sense that

E k h -E k-1 h t + J k h ≤ 0, ∀k = 1, • • • , N. (3.9)
Consequently, the discrete free energy is bounded and is decreasing in time

0 ≤ E k h ≤ E k-1 h ≤ H(s 0 ) L 1 (Ω) + G ∞ s 0 L 1 (Ω) , ∀k = 1, • • • , N. (3.10)
Furthermore, the time integral of the discrete dissipation is uniformly bounded

0 ≤ N k=1 tJ k h ≤ H(s 0 ) L 1 (Ω) + G ∞ s 0 L 1 (Ω) . (3.11)
Finally, there holds

π M, t H(s) L ∞ (0,t f ;L 1 (Ω)) ≤ H(s 0 ) L 1 (Ω) + G ∞ s 0 L 1 (Ω) . (3.12)
Proof. We proceed as in the proof of [START_REF] Cancès | Numerical analysis of a nonlinear free-energy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations[END_REF]Porposition 3.1]. By selecting v T = Υ(s k T ) in the formulation (2.13) we see

s k T -s k-1 T , Υ(s k T ) T + tJ k h = 0, ∀k = 1, • • • , N.
Thanks to the convexity of the function H(s) + sG and the nonnegativity of J k h we infer

E k h -E k-1 h ≤ E k h -E k-1 h + tJ k h ≤ s k T -s k-1 T , Υ(s k T ) T + tJ k h = 0.
This proves in particular (3.9). We now apply the Jensen inequality to obtain |K i | H(s 0 i ) ≤ K i H(s 0 ) dx for all i ∈ V. Whence, we show (3.10). The inequality (3.11) is a direct consequence of (3.9)-(3.10). The last inequality (3.12) results from (3.7), (3.10) and Assumption (H 2 ). This finishes up the proof. Lemma 3.1. Let (s k i ) i∈V be a family fulfilling the numerical scheme (2.10)-(2.12). Then, there exists a positive constant η h, t > 0 depending on the mesh parameters such that

s k i ≥ η h, t , ∀i ∈ V, ∀k = 1, • • • , N. (3.13)
Proof. By the mass conservation equality (2.14) and the assumption (1.6), one can find some i 0 ∈ V so that one has

0 < 1 |Ω| Ω s 0 (x) dx ≤ s k i 0 , and log(s k i 0 ) ≥ -C .
As a consequence of (2.7) -(2.8), we get

0 < s 0 i 0 := λ |Ω| Ω s 0 (x) dx ≤ λ s k i 0 ≤ s k i 0 j , ∀j ∈ V T \ {i 0 }.
Owing to (3.2) we claim for all T ∈ T i 0 , the set of elements sharing the vertex i 0 , that

s 0 i 0 log(s k i 0 ) -log(s k j ) 2 ≤ s k i 0 j log(s k i 0 ) -log(s k j ) 2 ≤ C h, t , ∀j ∈ V T \ {i 0 }. (3.14)
The role of λ is to obtain s 0 i 0 > 0 independently of s k j so that one can make use of [13, Lemma 3.10]. As a result, it is possible to estimate log(s k j ) in a straightforward way from (3.14). Note that the inequality (3.14) does not hold true if s k i 0 j = s ,k i 0 j since s ,k i 0 j is not far away from 0 with a positive constant independent of s k j . This is the reason why we perturb s k i 0 j to include a very small s c,k i 0 j satisfying the first inequality of (3.14). We deduce from (3.14) that, log(s k j ) ≥ -C h, t , ∀j ∈ V T \ {i 0 }, ∀T ∈ T i 0 . It is now sufficient to employ a similar procedure on the next vertex j and so on. This practical argument is therefore applied by induction to the set of mesh vertices, which is finite, see [13, Lemma 3.10]. Accordingly, one ends up with (3.13).

We now state without proof the existence of discrete positive solutions at each time level. The proof is based on the topological degree argument [START_REF] Deimling | Nonlinear Functional Analysis[END_REF]. The main idea consists in choosing a suitable homotopy function yielding a monotone scheme whose existence is shown by mimicking the guidelines of [START_REF] Eymard | Error estimates for the approximate solutions of a nonlinear hyperbolic equation given by finite volume schemes[END_REF]. Additionally, the energy estimate (3.3) and the positivity estimate (3.13) allow to prove the uniform continuity of the residual function, coming from the finite volume scheme (2.10)-(2.11), on a compact subset of (R + * ) #V . Then, no solution could exist on the boundary of this compact which directly implies the existence of at least one positive solution to the scheme in question. More details can be consulted in [14]. Proposition 3.3. For every k = 1, • • • , N , there exists at least a positive solution (s k i ) i∈V to the finite volume scheme (2.10)-(2.12).

Convergence of the finite volume scheme

In this section we are concerned with the convergence proof of the numerical scheme. To this end, we should establish some compactness arguments. A possible way to do that is to derive uniform estimations on the space and time translations [START_REF] Eymard | Finite volume methods[END_REF]Lemma 4.6]. Instead of using this traditional strategy, we here apply the practical time-compactness criterion recently elaborated in [START_REF] Andreianov | A nonlinear time compactness result and applications to discretization of degenerate parabolic-elliptic PDEs[END_REF]. This ensures the existence of a convergent subsequence of discrete solutions. Furthermore, the limit of the latter has to be identified to a weak solution. We first require the following results.

∈ R d , let us set Q t f ,a = Ω a × (0, t f ) with Ω a = {x ∈ Ω/[x, a] ⊂ Ω}.
Then, there exists a positive constant C independent of the discretization parameters such that the following integral satisfies

Qt f ,a |π M, t u(x + a, t) -π M, t u(x, t)| dx dt ≤ C |a| N k=1 t ∇π T u k 2 L 2 (Ω) d 1 2 , (4.1)
for all π M, t u defined by (2.1).

Let us now state and prove an estimate on the dual norm of the discrete time derivative. Let

ψ ∈ C ∞ c (Ω × [0, t f )).
In the sequel, we denote ψ k T the vector of

R #V such that (ψ k T ) i = ψ(x i , t k ) for all i = 1, • • • , #V. Lemma 4.2.
There exists a positive constant that depends only on the physical data and the regularity of the mesh such that

N k=1 s k T -s k-1 T , ψ k T T ≤ C ∇ψ ∞ , ∀ψ ∈ C ∞ c (Ω × [0, t f )). (4.2)
Proof. We consider v T = ψ k T in the concise form of the scheme (2.13). Then, we obtain the splitting X + Y = 0, where

X = N k=1 s k T -s k-1 T , ψ k T T , Y = N k=1 t T ∈T S T δ T Υ(s k T ) • A T S T δ T ψ k T .
Applying the Cauchy-Schwarz inequality on the last term, making use of 2.15, employing the smoothness of the test function ψ and (3.1) yield

|Y | ≤ N k=1 t T ∈T S T δ T Υ(s k T ) • A T S T δ T Υ(s k T ) 1/2 N k=1 t T ∈T S T δ T ψ k T • A T S T δ T ψ k T 1/2 ≤ C ∇ψ ∞   T ∈T |T | i∈V T s k iτ (i)   1/2 ≤ C ∇ψ ∞ .
This concludes the proof of the required inequality.

Q t f -→ (0, +∞), i.e. π Mm, t m s ----→ s a.e. in Q t f , as m -→ +∞. We recall that H(π Mm, t m s) L 1 (Qt f ) ≤ C.
Then, the sequence (π Mm, t m s) is equi-integrable thanks to De La Vallée Poussin criterion [START_REF] Bogachev | Measure theory[END_REF].

Consequently, Vitali's convergence theorem claims the strong convergence of the above sequence i.e.

π Mm, t m s ----→ s strongly in L 1 (Q t f ), as m -→ +∞.

In view of Lemma 2.1, we deduce π Tm, t m s = π Tm, t m ξ 2 (s) ----→ s strongly in L 1 (Q t f ) and a.e. in Q t f , as m -→ +∞.

To conclude the proof, it suffices to employ the standard identification of the limit process to check that ζ = ∇ξ(s) in the sense of distributions. Now, we are in a position to identify the limit function of the previous result as a weak solution. Proof. Following [START_REF] Cancès | Numerical analysis of a nonlinear free-energy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations[END_REF]Proposition 4.4], it can be checked that H(s) belongs to L ∞ (0, t f ; L 1 (Ω)). Let us now take v k Tm = ψ k Tm in the discrete weak formulation (2.13) and sum on k = 1, • • • , N . Then, the latter summation is decomposed as follows

X m + Y m + Z m = 0. (4.6) X m = N k=1 s k Tm -s k-1 Tm , ψ k Tm Tm , Y m = N k=1 t m T ∈Tm S T δ T log(s k Tm ) • A T S T δ T ψ k Tm , Z m = N k=1 t m T ∈Tm S T δ T G Tm • A T S T δ T ψ k Tm .
Notice that ψ N Tm = 0 R #Vm . Performing a discrete integration by parts in time implies

X m = -s 0 Tm , ψ 0 Tm Tm - N k=1 t m s k-1 Tm , ψ k Tm -ψ k-1 Tm t m Tm , = - Qt f π Mm, t m s(•, -t m ) Dπ Mm, t m ψ dx dt - Ω s 0 π Mm, t m ψ(x, 0) dx.
where Dπ Mm, t m ψ is the discrete time derivative given by

Dπ Mm, t m ψ(x, t) = ψ k i -ψ k-1 i t m , ∀(x, t) ∈ K i × (t k-1 , t k ], ∀i ∈ V, ∀k ≥ 1.
It follows from the smoothness of the test function ψ that (Dπ Mm, t m ψ) m (resp. (π Mm, t m ψ(•, 0)) m ) converges uniformly towards ψ t (resp. ψ(•, 0)). Thanks to (4.3) we get

X m = - Qt f s ψ dx dt - Ω s 0 ψ(x, 0) dx, as m -→ +∞.
Let us now study the convergence of the diffusive part Y m . To this purpose, we first consider

Y * m = 2 N k=1 t m T ∈Tm ξ(s k T,min )δ T ξ(s k Tm ) • A T δ T ψ k Tm , = 2 Qt f π Tm, t m ξ(s)∇π Tm, t m ξ(s) • ∇π Tm, t m ψ dx dt,
where s k T,min = min i∈V T {s k i } and we have defined

π Tm, t m ξ(s) = ξ(s k T,min ), ∀(x, t) ∈ T × (t k-1 , t k ], ∀T ∈ T m , ∀k ≥ 1.
Using Lemma 2.1 we estimate π Tm, t m ξ(s) -π Tm, t m ξ(s)

L 2 (Qt f ) ≤ π Tm, t m ξ(s) -π Tm, t m ξ(s) L 2 (Qt f ) ≤ h Tm C N k=1 t m ∇π Tm ξ(s k ) L 2 (Ω) d 1 2
≤ Ch Tm -→ 0, as m -→ +∞.

We deduce in particular that

π Tm, t m ξ(s) ----→ ξ(s) strongly in L 2 (Q t f ), as m -→ +∞. (4.7)
Thereby, we can pass to the limit in Y * m to obtain

Y * m -→ Y = 2 Qt f √ s∇ √ s • ∇ψ dx dt = Qt f ∇s • ∇ψ dx dt, as m -→ +∞.
We furthermore define the local diagonal matrix (S T ij ) 1≤i,j≤#V T such that

S T ij = s k T,min if j = τ (i) 0 if j = τ (i)
.

We denote S T,2 the diagonal matrix corresponding to the second choice of (2.6), namely the logarithmic mean. Note that Y * m can be rewritten under the form

Y * m = N k=1 t m T ∈Tm S T,2 δ T log(s k Tm ) • A T S T δ T ψ k Tm = N k=1 t m T ∈Tm δ T ξ(s k Tm ) • A T S T δ T ψ k Tm .
Let us moreover set

Y * * m = N k=1 t m T ∈Tm S T δ T log(s k Tm ) • A T S T δ T ψ k Tm .
Observe that s k T,min ≤ s k iτ (i) , for every i ∈ V T . Also, the regularity of ψ gives

|ψ k i -ψ k τ (i) | ≤ ∇ψ ∞ h T
for all i ∈ V T and k ≥ 0. We subtract the last two identities. We utilize once more the Cauchy-Schwarz inequality, the result (2.15), Proposition 3.1, and Lemma2.1 to find that

|Y * m -Y * * m | = N k=1 t m T ∈Tm S T δ T log(s k Tm ) • (S T,2 -S T )A T δ T ψ k Tm ≤ C ∇ψ ∞   N k=1 t m T ∈Tm |T | h 2 T i∈V T s k iτ (i) log(s k i ) -log(s k τ (i) ) 2   1 2 × N k=1 t m T ∈Tm |T | S T,2 -S T 2 ∞ 1 2 ≤ C ψ π Tm, t m ξ(s) -π Tm, t m ξ(s) L 2 (Qt f ) ≤ C ψ h Tm , as m -→ +∞.
Similarly, one shows that |Y m -Y * * m | -→ 0, as m -→ +∞. To conclude, it remains to treat the convective potential term Z m . Mimicking an analogous reasoning we consider

Z * m = N k=1 t m T ∈Tm S T δ T G Tm • A T S T δ T ψ k Tm = Qt f π Tm, t m ξ(s) 2 ∇π Tm G • ∇π Tm, t m ψ dx dt.
The regularity of the potential G entails the uniform convergence of the sequence (∇π Tm G) m towards ∇G in L 2 (Q t f ) d . Thanks to the strong convergence (4.7) we get

Z * m -→ Z = Qt f s∇G • ∇ψ dx dt, as m -→ +∞.
We continue in the same fashion by setting

Z * * m = N k=1 t m T ∈Tm S T δ T G Tm • A T S T δ T ψ k Tm .
Thanks to the smoothness of G we have G(x i ) -G(x τ (i) ) ≤ ∇G ∞ h T for all i ∈ V T and k ≥ 0. We follow the previous arguments together with (3.1) to discover

|Z * m -Z * * m | = N k=1 t m T ∈Tm (S T -S T )δ T G Tm • A T S T δ T ψ k Tm ≤ C ∇ψ ∞ N k=1 t m T ∈Tm |T | S T -S T 2 ∞ 1 2 ×   N k=1 t m T ∈Tm |T | h 2 T i∈V T s k T,min G(x i ) -G(x τ (i) ) 2   1 2 ≤ C ψ ∇G ∞ π Tm, t m ξ(s) -π Tm, t m ξ(s) L 2 (Qt f )   N k=1 t m T ∈Tm |T | i∈V T s k iτ (i)   1 2
≤ C ψ,G h Tm , as m -→ +∞.

We adapt the above procedure and the same steps to establish that |Z m -Z * * m | -→ 0, as m -→ +∞. Hence, the proof is complete.

Numerical tests

The goal of this section is to validate the efficiency and the robustness of our methodology through several examples including anisotropic and heterogeneous tensors. We solve the algebraic system (2.10)-(2.11) issuing from the numerical scheme by the Newton method that we implemented in Matlab. The Newton method computes a family of iterates(s k-1,η T ) η≥0 . To avoid the singularity of the logarithm, we reinforce the initial guess to be positive as follows s T = max(s T , 10 -10 ). In the case of its convergence, it tends to the sought solution s k T . Fixing the maximum number of iterations as 20, the stopping criterion is made on the difference of the successive iterates in the discrete 2 -norm i.e., the procedure stops when it yields

s k,η+1 T -s k,η T 2 ≤ ε,
where we set ε = 10 -10 .

Here we consider Ω = (0, 1) 2 the computational domain. It is covered by a family of refined triangular meshes coming from benchmarking problems [START_REF] Herbin | Benchmark on discretization schemes for anisotropic diffusion problems on general grids[END_REF]. We refer to Figure 2 for a depiction To evaluate the convergence speed of our approach, we calculate the numerical errors between the computed solution and the analytical one using the following norms

π M, t s -s e r = π M, t s -s e L r (Qt f ) , r = 2, ∞.
For clarity, in all the tests below, the logarithmic-mean (log-mean) corresponds to s iτ (i) defined in (2.7) when the parameter λ is assigned to the machine precision. That is, λ = 10 -16 so that s iτ (i) and s iτ (i) can be numerically identical. The arithmetic-mean is obtained for λ = 1.

Accuracy of the scheme

In this subsection we test the efficiency and the robustness of our scheme using a known exact solution to the continuous problem. We first consider a constant tensor D. Then, the analytical solution reads s e (x, y, t) = 1 + cos(πy)e -ct e x , where the final time is fixed to t f = 0.2 and c = π 2 D yy . The initial condition is then s 0 = s e (•, 0). This state degenerates on the segment [0, 1] × {1}. The exterior potential is set to G(x, y) = -x.

The second eigenvalue of D is taken as D yy (x, y) = 1. We also compare both formulas of the mobility function given in (2.6) in terms of the anisotropy ratio of the tensor D. The first expression i.e. s c iτ (i) will referred to as arithmetic-mean while the second one i.e. s iτ (i) will be called log-mean. In this example the time step is proportional to the size of the mesh i.e.

t/h 2 = 0.2. A summary of the time-space discretization is given in Table 1. We report in the following tables the obtained results. Each table displays, the errors as well as their convergence rates, the minimum (s min ) of the computed solution π M s k for k ≥ 1, and the number of iterations of the Newton solver, that we denote by Niter.

Mesh #V h T t
In Tables 234567, we observe that the second-order accuracy of the method is achieved for several anisotropic ratios less than 10 even on coarse meshes. This fact holds for both the log-mean and the arithmetic mean of the mobility-like function. We also record that the solution remains positive as established in Lemma 3.1. When the tensor becomes strongly anisotropic D xx ≥ 100, the underlined features are preserved for the log-mean case as provided in Tables 89. There is no huge difference in the number of Newton's iterations in the weak anisotropic case for the both averages. However, the arithmetic-mean choice turns out to break down the iterative solver. In order to converge, it necessitates a stringent constraint on the time step, that is t/h 2 should be less than 1/D xx , which allows extremely small time steps and makes the computational cost exorbitant.

The Newton solver corresponding to the arithmetic-mean choice fails to find the discrete solution when the anisotropic ratio is important. It produces spurious oscillations on the first iterates. This behavior is not recorded on the logarithmic-mean version of the solver. Reducing reasonably the time step or cutting the oscillations do not help to overpass the observed issue. In Figure 3, we plot the 2 -norm of the Newton residual function in terms of λ ∈ (0, 1) at the first time iteration (t 1 = t = 0.0031) on the second mesh and for D xx = 100, 1000. We observe that the norm augments as this parameter increases, meaning that we get close to the arithmetic-centered scheme. Through several experiments we also noticed that the norm of the residual function blows up when λ > 0.7.

Mesh

Large time behavior

This example focuses on the asymptotic behavior of our numerical scheme when time becomes large. It is inspired from [START_REF] Cancès | Numerical analysis of a nonlinear free-energy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations[END_REF]. We keep the same data of the first test except D and t f . According To illustrate the long time behavior of the proposed approach we look at the evolution of the discrete relative energy given as follows

E k h -E ∞ h = s k T log(s k T /s ∞ T ) -s k T + s ∞ T , 1 T T , ∀k = 0, • • • , N.
We represent on [14] C. Cancès and C. Guichard. Numerical analysis of a robust free energy diminishing finite volume scheme for parabolic equations with gradient structure. Foundations of Computational Mathematics, 17(6):1525-1584, 2017.

Formally, taking

  ψ = log(s) + G in the weak formulation (1.7) implies the energy-dissipation relationship dE dt + I = 0. Using the fact that I ≥ 0 we get dE dt = -I ≤ 0. (1.8)

Figure 1 :

 1 Figure 1: The spatial discretization of T ∈ T (left) and illustration of the dual mesh (right) given by dashed segments.

Lemma 4 . 1 .

 41 (see[START_REF] Eymard | Finite volume methods[END_REF] Lemma 4.2],[START_REF] Ghilani | Positive control volume finite element scheme for a degenerate compressible two-phase flow in anisotropic porous media[END_REF] Lemma 8.2]) For a given a

Proposition 4 . 2 .

 42 Keeping the statement of Proposition 4.1, the obtained function s is indeed a weak solution to the continuous model (1.1)-(1.3) in the sense of Definition 1.1.

Figure 2 :

 2 Figure 2: First two meshes of the considered triangulation.

Figure 3 :

 3 Figure 3: Impact of the parameter λ on the 2 -norm of the Newton residual function for D xx = 100 (left) and D xx = 1000 (right) using the second mesh at the first time iteration t 1 = t = 0.0031.

Figure 4 :

 4 Figure 4: Relative free energy with the log-mean choice of the mobility (left) and the arithmetic-mean (right) for D xx = D yy = 1, and t f = 4 on the first four meshes.

Figure 4 -

 4 5 this quantity as a function of time for the log-mean and arithmeticmean approximations of mobility in the isotropic and the anisotropic cases. The first plot addresses the isotropic case with D xx = D yy = 1 and t f = 4. It shows that the relative free energy is rapidly diminishing towards 0. The second plot concerns the anisotropic case D xx = 1 and D yy = 0.05 with t f = 55. We also observe the exponential decreasing of the free energy towards its stationary state. This requires more time, which is due to the impact of the tensor anisotropy.

Figure 5 :

 5 Figure 5: Relative free energy with the log-mean choice of the mobility (left) and the arithmetic-mean (right) for D xx = 1, D yy = 0.05, and t f = 55 on the first four meshes.

Figure 6 :

 6 Figure 6: From left-up to right-bottom, the log-mean approximate solution for different simulation times t ∈ {0.25, 1, 4, 10}.

Figure 7 :

 7 Figure 7: From left-up to right-bottom, the arithmetic-mean approximate solution for different simulation times t ∈ {0.25, 1, 4, 10}.

Figure 8 :

 8 Figure 8: Log-mean (left) and arithmetic-mean (right) cross section of the computed solution at the point (0, 0.4) for t ∈ {0.25, 1, 4, 10}.

[

  13] C. Cancès and C. Guichard. Convergence of a nonlinear entropy diminishing control volume finite element scheme for solving anisotropic degenerate parabolic equations. Mathematics of Computation, 85(298):549-580, 2016.

  Proposition 4.1. We assume that (H 1 )-(H 3 ) are fulfilled. Let (T m ) m be a sequence of refined meshes to T such that h Tm , t m go to 0 as m tends to infinity. Then, there exists a measurable function s : Q t f -→ (0, +∞) so that the following convergences hold up to a subsequenceπ Mm, t m s, π Tm, t m s ----→ s a.e. in Q t f , and strongly in L 1 (Q t f ),(4.3)Proof. According to the energy estimate (3.3) we have the uniform boundedness in L 2 (Q t f ) d of the sequence (∇π Tm, t m ξ(s)) m . Then, there exists a functionζ ∈ L 2 (Q t f ) d such that ∇π Tm, t m ξ(s) ----→ ζ, weakly in L 2 (Q t f ) d , as m -→ +∞. (4.5)Next, the mass conservation property gives a uniform estimation on the sequence (π Mm,t m ξ(s)) m in L 2 (Q t f )where we recall that ξ(s) = √ s. Thanks to Lemma 2.1, we infer that (π Tm, t m ξ(s)) m is uniformly bounded in the same space. Now, by virtue of Lemma 4.1 and Lemma 4.2 we apply the compactness argument [2, Theorem 3.9] to ensure the existence of a subsequence converging a.e. to a measurable function s :

	∇π Tm, t m ξ(s) ----→ ∇ξ(s)	weakly in L 2 (0, t f ; H 1 (Ω)).	(4.4)

Table 1 :

 1 Discretization data used to assess the accuracy of the scheme.

	Tri 1	37 0.250 0.125 E-01
	Tri 2	129 0.125 0.312 E-02
	Tri 3	481 0.063 0.793 E-03
	Tri 4 1857 0.031 0.192 E-03
	Tri 5 7297 0.016 0.512 E-04

Table 2 :

 2 π M, t s -s e 2 Rate π M, t s -s e ∞ Rates min Numerical errors with the log-mean choice of the mobility and D xx = 0.1.

	Niter

Mesh π M, t s -s e 2 Rate π M, t s -s e ∞ Rate

Table 3 :

 3 Numerical errors with the arithmetic-mean choice of the mobility and D xx = 0.1. Mesh π M, t s -s e 2 Rate π M, t s -s e ∞ Rate s min

	Niter

Table 4 :

 4 Numerical errors with the log-mean choice of the mobility and D xx = 1.

	Mesh	π M, t s -s e 2 Rate	π M, t s -s e ∞ Rate	s min	Niter
	Tri 1	0.148 E-01	-	0.143 E-00	-	0.112 E-00	45
	Tri 2	0.409 E-02	1.860	0.412 E-01	1.802 0.308 E-01 116
	Tri 3	0.107 E-02	1.956	0.108 E-01	1.950 0.798 E-02 395
	Tri 4	0.262 E-03	1.982	0.269 E-02	1.965 0.194 E-02 1309
	Tri 5	0.696 E-04	2.008	0.709 E-03	2.016 0.515 E-03 4037

Table 5 :

 5 Numerical errors with the arithmetic-mean choice of the mobility and D xx = 1.

	s min	Niter

Mesh

π M, t s -s e 2 Rate π M, t s -s e ∞ Rate

Table 6 :

 6 Numerical errors with the log-mean choice of the mobility and D xx = 10. Mesh π M, t s -s e 2 Rate π M, t s -s e ∞ Rate s min

	Niter

Table 7 :

 7 Numerical errors with the arithmetic-mean choice of the mobility and D xx = 10.

Table 8 :

 8 Mesh π M, t s -s e 2 Rate π M, t s -s e ∞ Rate s min Numerical errors with the log-mean choice of the mobility and D xx = 100. Mesh π M, t s -s e 2 Rate π M, t s -s e ∞ Rate s min

	Niter

Table 9 :

 9 Numerical errors with the log-mean choice of the mobility and D xx = 1000.

Heterogeneous five-spot problem

The objective of this last test is to exhibit the applicability of our scheme to a situation where the anisotropy is dependent on the position. This example is termed the five-spot problem in the context of porous media flows. Let us set Ω = Ω 1 ∪ Ω 2 such that Ω 1 = {(x, y) ∈ Ω/x < 0.5} and Ω 2 = Ω \ Ω 1 . Let us next consider the heterogeneous tensor

The domain is initially empty, that is s 0 = 0. The potential is chosen as G(x, y) = x. We impose a Dirichlet boundary condition s |∂Ω D = 1 on the left portion ∂Ω D = {x = 0} × {0 ≤ y ≤ 0.5} and an outflow condition on the right zone located at {x = 1} × {0.5 ≤ y ≤ 1}. The remainder border is assumed impervious. For the simulation, we use the third mesh of the family depicted in Figure 2. We consider the fixed time step t = 0.00089. The final simulation time is set to t f = 0.1. The obtained results are plotted on Figure 67. The first figure shows the solution computed accurately using the log-mean approximation of the mobility while the second one exhibits the solution of the scheme using the arithmetic mean. Both approaches produce the expected behavior. It is observed that the medium Ω 1 is more permeable than Ω 2 , which is due the the contrast of the intrinsic permeabilities. We moreover remark that our scheme generates no undershoots. These facts are finally confirmed by the cross section of the solution at the point (0, 0.4) illustrated in Figure 8. To sum up, both strategies behave similarly as the plotted solutions show, since the tensor is weakly homogeneous in each subdomain and the ratio is small.


C Nguyen

A Merle

Z O Guimaraes-Filho

X Garbet

R Sabot

J Decker

L.-G Eriksson

M Goniche

G Huysmans

H Lütjens

P Maget

Theoretical and Experimental Analysis of the Destabilization of Modes Driven by Fast Particles in Tore-Supra

A necessary step towards the control of burning plasmas is the understanding of instabilities driven by suprathermal particles, which are known to be dangerous for the plasma confinement. The mechanisms underlying their destabilization are studied for two modes identified in Tore-Supra. Linearly, experimental and analytical evidences are presented which show the strong sensitivity of the fast electron driven fishbone stability on the current profile, and the non-negligible role of diamagnetic effects and fast particle induced non-perturbative effects on non-chirping steady-state Beta Alfvén Eigenmodes. Nonlinearly, the main damping mechanism acting on Beta Alfvén Eigenmodes is shown to enable the existence of subcritical modes.

Introduction

On the way towards the understanding of burning plasmas, studies of heating induced suprathermal populations are of major importance. Such populations are known to excite Magneto-HydroDynamic (MHD) modes which can have deleterious effects on the plasma energy confinement. In Tore-Supra, Ion Cyclotron Resonance Heating (ICRH) induced suprathermal ions are found to drive acoustic modes called Beta Alfvén Eigenmodes (BAEs) [1,2]. Lower Hybrid (LH) heated electrons are also thought to excite meso-scale modes called fishbone modes albeit characterized by smaller radiogyradii, which is a good indicator that energy is the main parameter of the plasma-fast particle interaction, and that the dynamics studied in current devices is relevant to what should occur in larger reactor-size devices.

In the present work, we report on various experimental and theoretical studies carried out in Tore-Supra, to understand the mechanisms and conditions underlying the interaction of suprathermal particles with MHD modes. An upgrade of the reflectometry and Electron Cyclotron Emission diagnostics associated with correlation and fitting techniques makes possible to access the modes radial structure, which reveals noteworthy dynamics at the current diffusion time scale. Such observations motivate a detailed analysis of the two modes linear stability, and in particular of the role of diamagnetic effects, non-perturbative fast particle induced kinetic effects and current profile shaping. Such accurate analysis is necessary to enable a control of fast particle driven modes, as well as for MHD based plasma diagnostics. In section 2, Tore-Supra experimental observations and parameters are outlined. Next, the model and results of the linear analysis of the modes stability are given in section 3. Finally, on a more theoretical side, the role of nonlinear effects such as particle wave-trapping in the stability of the Beta-Alfvén Eigenmodes are investigated in section 4.

2 Experimental parameters and radial structure of fast particle driven modes Tore-Supra fast particle induced modes have been detected on the reflectometry and Electron Cyclotron Emission (ECE) diagnostics. Reflectometry is made of 2 channels, and ECE of 32 channels EXS/P7-02 which can be adjusted in one line of sight of the equatorial plane, hence allowing for a reconstruction of the modes radial structure and providing information on the mode poloidal mode number parity. For both BAEs and fishbones, a particular cyclical dynamics of the oscillation structures is found, which appears in close relation with the evolution of their frequency, as presented in Figs. 1,2.

Beta Alfvén Eigenmodes

BAEs have been found in sawtoothing high density discharges (∼ 4 × 10 1 9m -3 ) characterized by a core electron temperature of about 4-5 keV, and a plasma current of 0.9-0.95 MA, in the presence of ICRH of typical power 2MW to 4.6 MW [2]. For the shot #42039, shown in Fig. 1 and analyzed below, the simulation of the current profile by the CRONOS transport code [3], returns a monotonic q-profile with q=1 surface at a normalized radius of 0.32. At this point, the density is equal to n = 3.7 × 10 19 m -3 , the core ion temperature (obtained from the TPROF fit of neutron emission data [3]) to T i = 2.1 keV, the ratio of core electron to core ion temperatures to 1.3, the typical density and temperature scales are ∇n/n = -1.5m -1 and ∇T i /T i = -3.1m -1 .

As shown in Fig. 1b. and c., several modes are found to oscillate in the acoustic frequency range (around 50 kHz). Radially, the analysis of the three most intense modes shows that they are all localized inside the q = 1 surface, almost stable during the sawtooth period, yet with a slight radial slight [START_REF] Guimaraes-Filho | 37 th EPS Conference on Plasma Physics[END_REF]. Poloidally, phase-angle analysis of low-field and high-field side perturbations (Fig. 1c. ) suggests that the multiplicity of modes corresponds to modes with different poloidal mode numbers. Now looking at their frequency spectrum (Fig. 1b.), a particular dynamics comes forth. While quasi steady state at the millisecond time scale, the mode frequencies display a downwards-upwards cycle which matches the sawtooth period, and is more pronounced for the lower branches, leading to a splitting and merging of the modes. This dynamics is reminiscent of observations performed at Asdex-Upgrade [START_REF] Lauber | [END_REF], and suggests a role of temperature and mode number related physics, such as diamagnetic effects.

EXS/P7-02

Electron fishbone modes

Lower frequency modes with frequencies ranging from 8 to 14 kHz have been detected in discharges with typical density of 2-3e19m -3 , a central temperature of about 4 keV, a plasma current of 0.6MA and in the presence of 1 to 2MW of LH heating. They have been identified as fishbone modes, in their quasi steady-state (non-bursting) regime also identified in Ref. [6]. As was reported in Ref. [7], the frequency of these modes is characterized by jumps illustrated in Fig. 2b., associated to modulations of the central temperature, sometimes in the form of a neat oscillation (Fig. 2a.) with a time scale corresponding to a current diffusion time scale.

Using the matching techniques described in Ref. [START_REF] Guimaraes-Filho | 11th IAEA Conference on Energetic Particles[END_REF], the radial localization of the modes has been determined (Fig 2c .). It shows that the frequency jumps are associated to a continuous inward drift of a mode, and the apparition of a second mode at the first mode initial position, close to a q = 1 surface provided by the CRONOS transport code at the normalized radius 0.16. Remarkably, the discontinuous frequency jumps are associated to continuous evolutions of the mode position.

This evolution of the mode structure suggests a modication of equilibrium properties. More precisely, the correlation between the mode localization and a q=1 surface and the slow inward drift of one of the modes advocates for a cycle of the central q-profile, not accessible to a sufficient accuracy by the simulation. Moreover, the discontinous transition between modes indicates that this q-profile oscillation has a strong impact on stability. Such observations suggests the necessity of a detailed analysis of the role of the q-profile in the electron fishbone stability.

3 Linear analysis of the mode stability and oscillation frequency

Developed tools

As was shown in Ref. [START_REF] Zonca | 34th EPS Conference on Plasma Physics[END_REF], the linear dispersion relation of fast particle driven modes of the low shear Alfvén spectrum can be put in the same form

iΛ|s| = δW f + δW k (1)
where Λ represents the thermal plasma inertia including second order Finite Larmor Radius effects, δW f is the normalized plasma non-resonant fluid energy and δW k is the normalized kinetic resonant contribution of fast particles. Multiple derivations of Λ for BAEs [START_REF] Zonca | [END_REF]11,12,13,14] and for fishbone modes [15,16,17] exist, which can be consistently combined [START_REF] Zonca | [END_REF]. In this work, we considered Λ to be the square root of Λ 2 such that ∂Λ r /∂ω r > 0, where,

Λ 2 = 3 ω(ω -ω * i ) ω 2 A
for the lower frequency spectrum (Fishbone) (2)

Λ 2 = ω 2 ω 2 A 1 - ω m * pi ω + q 2 ωω ti ω A 1 - ω m * ni ω F ω ω ti - ω m * Ti ω G ω ω ti - N m (ω/ω ti ) 2 N m+1 (ω/ω ti ) D m+1 (ω/ω ti ) + N m-1 (ω/ω ti ) D m-1 (ω/ω ti ) for the higher spectrum (BAE) (3)
In these formulas, corresponding to the derivation of Ref. [START_REF] Zonca | [END_REF] with the poloidal non-degeneracy correction added in Ref. [12],

F(x) = x(x 2 + 3/2) + (x 4 + x 2 + 1/2)Z(x), G(x) = x(x 4 + x 2 + 2) + (x 6 + x 4 /2 + x 2 + 3/4)Z(x), N m (x) = (1 -ω m * ni /ω)[x + (1/2 + x 2 )Z(x)] -(ω m * Ti /ω)[x(1/2 + x 2 ) + (1/4 + x 4 )Z(x)], D m (x) = (1/x)(1 + 1/τ e ) + (1 -ω m * ni /ω)Z(x) -(ω m * Ti /ω)[x + (x 2 -1/2)Z(x)].
(m,n) represent the main poloidal and toroidal mode numbers, ω A = v A /qR with v A the Alfvén velocity, R the major radius and q the safety factor at the mode resonance surface r s , ω ti = v ti /qR with v ti = EXS/P7-02 2T i /m i the core ion thermal velocity, and τ e = T e /T i with T i and T e the equilibrium core electron and core ion temperatures. Diamagnetic effects are taken into account via the quantities

ω m * Ti = (m/r)(T i /e i B 0 )(∇ r T i /T i ), ω m * ni = (m/r)(T i /e i B 0 )(∇ r n i /n i ), ω m * pi = ω * Ti + ω * ni .
δ Ŵf depends on the mode global structure. For an |n| = |m| = 1 kink structure, expected for a traditional fishbone mode and possible for an |n| = |m| = 1 BAE [17,13], it is [18]:

δW f = 3π(1 -q min ) r 2 s R 2 13 144 -β 2 ps . (4)
where

β ps = (-R/r 2 s )
R r s 0 r 2 (dβ p /dr)dr and β p is the plasma poloidal normalized pressure. Finally, the hot particle contribution (referred to with the index h) included in δW k is derived from kinetic theory to include wave-particle resonances. For the modes considered in this paper and the energies expected for fast particles present in Tore-Supra experiments (about 300keV for fast ions and 50 to 80 keV for the suprathermal electrons), main wave-particle resonances are expected to take place with the particles drift precession frequency Ω d , such that,

δ Ŵk = 2 3/2 π 2 |m| qm 3/2 h µ 0 B 2 0 × Z a 0 1 r 2 s rdr| ξ| 2 Z 1+r/R 1-r/R dλ Ω2 d Ωb Z +∞ 0 E 5/2 dE QF h nΩ d -ω . ( 5
)
where integrations are taken over the radial and energy and anisotropy variables, (r, E, λ = µ/E) (µ is the magnetic invariant). Here, QF h is an operator acting on the hot particle distribution function

F h , such that QF h = ω∂ E F h + nω m * ph /T h , ξ = ξ/(k θ Ψ 0 /B 0 ωs)
is the mode MHD displacement, normalized to its asymptotic values (see [17,13]). Ωb and Ωd are the normalized particle bounce and drift frequencies (used with the convention of Ref. [2]).

A solver, MIKE (MHD Instability with Kinetic Electrons), has been developed to solve the dispersion relation (1) with the different types of inertia presented in Eqs.( 2) and (3). It is based on equilibrium data derived from the CRONOS transport code and has been coupled to the C3P0/LUKE kinetic code to precisely access the complete electron distribution function in the presence of LH heating [START_REF] Peysson | Joint Varenna-Lausanne International Workshop[END_REF]. Although initially developed for fast electrons driven modes, the solver has also been used for fast ions with a more simple ion distribution function, calculated with the PION Monte-Carlo code [START_REF] Eriksson | [END_REF] and fitted on an anisotropic Maxwellian of the form

F h = ñ(r)δ(λ -λ 0 )e -E/T h (r) with ñ(r) = [n h (r)/(2πm h T h ) 3/2 ][2 Ωb (λ 0 , r)]. Explicity, δW k assumed with such a distribution is δ Ŵk = πq(r s ) 2 |m| Z a 0 1 r 2 s rdr| ξ| 2 β h (r) Z 1+r/R 1-r/R dλ Ωb (λ 0 ) Ω2 d Ωb × - E res T h Z 5 (y 0 ) - E res T h Z 3 λδ (λ 0 )(6) - 1 Ωd R∂ r n h n h Z 5 (y 0 ) - R 0 ∂ r Ωb (r, λ 0 ) Ωb (r, λ 0 ) Z 5 (y 0 ) + R 0 ∂ r T h T h Z 7 (y 0 ) - 3 2 Z 5 (y 0 ) with β h = 2µ 0 n h T h /B 2 0 , E res = [r/nq(r)][e h B 0 Rω/ Ωd (r, λ 0 )], √ πZ p (y 0 ) = R +∞ -∞ dyy p+1 exp(-y 2 )/(y 2 - y 2 0 ), for y 2 0 = E res /T h .

Beta-Alfvén Eigenmodes

In this section, the roots of the dispersion relation are compared with Tore-Supra observations, in order to determine the relevant physical ingredients determining the BAE frequency and stability. In Ref. [START_REF] Lauber | [END_REF], diamagnetic effects were found to account for the BAE mode cycle during a sawtooth period at Asdex-Upgrade. Indeed, diamagnetic effects are directly related to the temperature gradient of the thermal ion population, which is known to be almost anihilated and flat in the plasma center after the sawtooth crash, and to get slowly restored during the sawtooth period. To mimic this dynamics, the dispersion relation with δW k = δW f = 0 is solved in Fig. 3a. for various values of the EXS/P7-02 rate are extremely sensitive to the shear. The growth rate in particular increases significantly close to its marginal point with a very small shear flattening : this is the indication that electrons characterized by a reversal of their precessional drift (whose number increases at lower shear), are fundamental for the mode excitation. This sensitivity suggests the possibility to reconstruct the q-profile from MHD observations, in particular in the presence of a cycle illustrated in Fig. 2.

Conditions for the existence of subcritical Beta Alfvén Eigenmodes

In this section, we study another feature of the fast particle driven stability, the existence of subcritical (linearly stable and nonlinearly unstable) modes. Finite frequency modes of the shear Alfvén spectrum which are preferentially driven by suprathermal particles, are mostly damped by resonant, Landau-like damping and continuum damping. In the nonlinear phase of a mode growth, both damping processes can be affected by the mode size or nonlinear frequency. Hence the possiblity to reach more unstable states nonlinearly.

Here, the focus is on modes which are simultaneously driven and damped resonantly. It is the case for BAEs which oscillate in a region where thermal ion Landau damping is strong, and can be resonantly driven by a population of fast ions. In order to prove the existence and triggering mechanisms of subcritical modes, a simple perturbative 1d model is used, which should be extended to non-perturbative systems as suggested in the previous section.

Model

We make use of a Poisson-Vlasov system with two species of distribution functions F 1 and F 2 , whose applicability to BAEs has been discussed in Ref. [START_REF] Nguyen | Submitted to Plasma Physics and Controlled Fusion[END_REF],

∂ t F i + v∂ x F i + e i E m i ∂ v F i = -ν i F i -F i eq (v) , i = 1, 2 (7)
-∂ x E = ∑ i=1,2 e i Z dv F i -F i eq (v) . (8)
excited by a initial perturbation of the form F 1 (x, v,t = 0) = F 1eq (1 + α cos(kx)). F 1 eq , F 2 eq are the species equilibrium distribution functions given by

F 1 eq = [n b / √ 2πv tb ] exp -1/2 [v/v tb ] 2 + [n e / √ 2πv te ] exp -1/2 [(v -v 0 )/v te ] 2 , F 2 eq = [n d / √ 2πv td ] exp -1/2 [v/v td ] 2 .
Hence, F 1 eq is divided into a thermal bulk (index b), and a bump on tail (index e) providing resonant excitation and is similar to the distribution used in traditional collisional bump-on-tail models [START_REF] Berk | [END_REF]24]. F 2 has a purely thermal equilibrium distribution, intended to provide self-consistent resonant damping (index d).

In the perturbative limit (n b n e , n d ) where the system is expected to oscillate approximately at the bulk plasma frequency, ω pb ≡ n b e 2 b /m b , and assuming a cold bulk ω pb kv b , a slighly dissipative plasma ω pb ν b , ν e , ν d and a monochromatic perturbation, the linear bump induced resonant excitation and damping are well assessed by the independent rates γ e and -γ d (γ d > 0) (explicitely

γ s=e,d = 0.5(ω pb /k 2 )(ω 2 ps /v 2 ts ) ImY (ω pb -δ se kv 0 )/ √ 2kv ts , with Y = (1 + Z)
and Z is the plasma dispersion function). When the mode amplitude is sufficiently large compared to restoring dissipative effects, ω Bs ν s (for s = e, d) and ω Bs ≡ |e s |Ek/m s , particles get trapped inside the mode structure, leading to a reduction of the resonant contribution. In a 1-d model with Krook collision operators, the nonlinear resonant contribution has been shown to be reduced by a factor cν * s = cν s /ω Bs 1 with α = 2.0 [START_REF] Berk | [END_REF]. Following this reduction, steady-state subcritical (or metastable) modes can take place due to a stronger nonlinear reduction of damping compared to drive. More precisely, using modified

Refinement of the low amplitude model

The nonlinear coefficients used in the results presented in the previous paragraph, c s = αν s /(αν s + ω Bs ) is a fit between the quasi-linear regime (ω Bs → 0) for which c s γ s → γ s , and the deeply trapped nonlinear regime (ω Bs ν s ), for which c s γ s → α(ν s /ω Bs )γ s . We now make use of a refined description of near-marginality (ω Bs → 0), derived from the model of Ref. [26] : c s γ s ∼ 1λ(ω Bs /ν s ) 4 γ s . Fitting this refined model with the deeply trapped regime, and assuming closeness to a saturation point (d/dt ω Bs , s = e, d), the evolution of the mode amplitude follows the equation [27], where α = 2.0 and λ = 1/8. Similarly as for the simple model [25] and for our set of parameters, the investigation of saturation points (d/dt = 0 which we solved numerically) returns one single positive stable solution in the linear unstable region, while giving zero or two solutions (a smaller unstable solution with d/dt < 0, and a larger stable solution with dt > 0) in the subcritical region. We compare the simulated amplitudes with the predicted one in A good agreement is found between the predicted and simulated mode amplitudes. The refined model is found to bring some improvement ( 17 % error, compared to 22 % for the simpler model). in particular for lower saturation amplitudes. The remaining discrepancy may be attributed to the breaking of the assumption of monochromaticity.

dω 2 Bb dt = ω 2 Bb 1 + αλ(ω Be /ν e ) 3 1 + αλ(ω Be /ν e ) 3 + λ(ω Be /ν e ) 4 γ e - 1 + αλ(ω Bd /ν d ) 3 1 + αλ(ω Bd /ν d ) 3 + λ(ω Bd /ν d ) 4 γ d -ν b . (9)

Triggering condition

The presented refinement is effective at low amplitudes, and in particular to assess the accuracy of the lower unstable root in the metastable region ω - Bb . From Eq. ( 9), a sufficient condition to obtain a metastable mode is the use of an initial perturbation ω Bb (0) ≥ ω - Bb . It is notable that with such a model, very small perturbation amplitudes (ω - Bb → 0 close to marginality) can lead to much higher stable saturation amplitudes (∝ γ NL close to marginality). In systems where metastability is expected, that is, when ν * e ν * d , it can be easily shown that the refined model has lower roots that the simpler model, which makes the triggering of a metastable mode even easier. However, the validity if Eq. ( 9) is reduced to systems with a reduced time evolution, which implies an additional condition on the inertial perturbation, of the form ω Bs > -γ LIN for all resonant species.

Conclusion

In this paper, we presented the analysis of two modes found in Tore-Supra. We put forward a particular dynamics of fast electron driven fishbones at the slow current diffusion time scale, characterized by a discontinuous frequency spectrum associated to a continuous structure evolution. Such an observation suggests an important effect of the current profile on the mode stability, which was confirmed by theory for a monotonous q-profile configuration. Next, the relevance of diamagnetic effects on Tore-Supra Beta Alfvén Eigenmodes was estimated and shown to be consistent with one part of the observed mode cycle during a sawtooth period. Indications were also given that fast particles should be considered non-perturbatively, even for non-chirping modes. Finally, we explained how subcritical Beta Alfvén Eigenmodes can take place, and presented necessary and sufficient conditions for the onset of such steady-state metastable modes. This work was supported by the European Commission.

Figure 1 :

 1 Figure 1: BAE cycle during a sawtooth period. a) is the core electron temperature, b) the reflectometry frequency spectrum and c) the phase angle between low and high field side perturbations for shot #42039. d) is the BAE radial structure obtained with ECE for an equivalent shot #44831.

Figure 2 :

 2 Figure 2: Fishbone cycle. a) is the core electron temperature b) the ECE frequency spectrum c) the mode radial structures obtained from a fit of the ECE signals.

Figure 5 :

 5 Figure 5: Simulated nonlinear states obtained, varying resonant damping and collisionality.reduction factors of the form αν s /(αν s +ω Bs ), we showed analytically and numerically in Ref.[25] that positive stationnary and stable solutions are possible : (i) in the complete linear region whereγ L ≡ γ eγ dν b > 0, (ii) in the subcritical region, if γ NL ≡ ν * e (γ eν b )ν * d (γ d + ν b ) > 0 and ∆ ≡ (γ NL ) 2 -4ν * e ν * d ν b |γ L | > 0.A numerical example of the regions determined by these conditions is shown in Fig.5along with the different nonlinear states obtained from our Poisson model for a particular set of parameters verifying ν * e ω Bb ∈ [0.001, 0.1] ν * d ω Bb = 0.003, and a large initial amplitude ω Bb (0) = 0.2 > ν e ∈ [0.001, 0.1], ν d = 0.005. A large number of stead-state subcritical states are found corresponding to the prediction.

Figure 6 :

 6 Figure 6: Comparizon between simulated and predicted saturation amplitudes.

Fig. 6

 6 for all simulated steady-state modes with the fixed parameters e s = 1.0 for s = b, e, d, m s = 1.0 s = b, e, n b = 1.0, v tb = 0.3, n e = 0.03, v te = 1.0 v 0 = 4.5, v td = 2.5 and k = 0.3, the variable parameters (m d /m l , ν d ) ∈ {(2.0, 0.04), (0.5, 0.005)} and full scans in n e ∈ [10 -3 , 1], ω 2 pd ∈ [0.0, 0.35].

normalized temperature gradient η = (∇T i /T i ) r s /(∇n i /n i ) r s ranging from 0 to 2.0 (that is, the average gradient estimated by the TPROF code), and for several mode numbers |n| = |m| ∈ [1,[START_REF] Lauber | [END_REF]. Although the assessed root significantly overestimates the mode experimental frequency (possibly due to a finite particle perpendicular energy [START_REF] Lauber | [END_REF] or to 4th order Finite Larmor Radius effects [21]), a downward shift of the order of 1 to 5 kHz is obtained in Fig. 3a. which is lower but of the same order as the experimental values in the first half of the sawtooth period (from 3 to 10 kHz). Similarly as for the experiment, a splitting occurs due the frequency |m|-dependence. Hence, diamagnetic effects appear as good candidates to explain the first half of the BAE cycle. Yet, from this test, they do not seem to provide a direct explanation for the second half of the sawtooth period, where the modes frequency increases significantly (up to a 30 kHz shift) while the core temperature gets stabilized.

Next, the non-perturbative role of the fast ion distribution is analyzed in Fig. 3c.d., using an artificial fraction ε, such that δW k = ε δW k experimental ( δW f is kept null) for the |n| = |m| = 1 mode (for which the use of a kink step-like MHD displacement for ξ in Eq. 5 makes some sense). As was already suggested by an earlier perturbative analysis of the BAE stability in Tore-Supra, nonperturbative effects such as a significant frequency shift (see the horizontal mismatch between theory and experiment in Fig. 12 of Ref. [2]) need to be assessed. In Fig. 3c., the frequency shift induced by the experimental distribution (ε = 1.0) is about 15 kHz, compared to the BAE accumulation point (given by ε = 0). This suggests that fast particles can have a role in the BAE frequency cycle, and that a non-perturbative treatment of stability is necessary (for such an ICRH power), in particular because such an upward shift can lead to stronger continuum damping. Nevertheless, a more refined model (for the radial structure, the role of |n|, details of the distribution function...) is still needed for this study (in particular, the mode excitation Fig. 3d. is too weak). 4), and a complete electron distribution simulated by LUKE for shot #40816. In the latter shot, an oscillation was detected at a normalized radius of 0.18, which matches a reasonable q=1 surface. Its frequency was found to be relatively stable, of about 8.5kHz and interrupted by very limited frequency jumps [7], which makes it appropriate for a stability analysis. It occured for an LH power of 1MW, at a density of 2.3e-19m -3 .

Electron fishbone modes

In Fig. 4, the role of the shape of the q-profile in the fishbone stability is analyzed. For this, the q-profile is enforced in δW k and chosen to be parabolic with a fixed q=1 surface at the normalized radius of 0.18, and a variable shear s = (r/q)dq/dr at this point. First, this analysis is in agreement with experimental observations. Indeed the q-profile corresponding to the frequency of 8.5kHz, the growth rate is found positive. Secondly, one can note that the mode frequency and growth


T Nicolas

R Sabot

X Garbet

H Lütjens

J.-F Luciani

R Guirlet

J Decker

A Sirinelli

Role of the sawtooth crash in the electron and impurity transport in the Tore Supra and JET tokamaks

Introduction

In the Tore Supra and JET tokamaks, the plasma density equatorial profile is measured with fast-sweeping reflectometry, in ohmic sawtoothing regime. The profiles are used to reconstruct the 2D density during an (m, n) = (1, 1) kink mode, following the tomography technique described in Ref. [1]. Crescent-shaped structures of amplitude 1% of the total density are observed inside the q = 1 surface immediately following the sawtooth crash. Despite the small amplitude of the structure, it constitutes a significant deviation with respect to the expected post-crash flat profiles. In Ref. [2], this behaviour was modelled using the 3D non-linear full MHD code XTOR-2F [3]. This code is capable of simulating several self-consistent sawtooth crashes [4].

Despite full reconnection occurring in the simulation, contrary to most experimental cases, the experimentally observed density structures were recovered and it was demonstrated that they arise because of the fast flows at the reconnection layer which first advect the density along the separatrix and then reinjects it inside the q = 1 surface. Thus the redistribution mechanism of the density is fundamentally different from that of the temperature, the evolution of which is governed by the fast diffusive transport along the reconnected field lines. The question of the particle transport induced by the sawtooth instability has been given more importance lately due to recent experiments on the JET tokamak with the ITER like wall [START_REF] Pütterich | In 24 th IAEA Fusion Energy Conf[END_REF], where Tungsten (W) is seen to accumulate in the core despite frequent sawtooth crashes.

The present contribution shows JET observations of the density in ohmic regime with the same diagnostic as in Tore Supra, confirming the results of this machine, and presents results obtained with a new module of the XTOR-2F code, allowing to simulate the density evolution of several impurities considered as passive scalars, evolving in the MHD flow of the main plasma.

The amount of particles reinjected inside the q = 1 surface, versus effectively expelled by the sawtooth crash, is determined, and conclusions are drawn.

Reflectometry results on the JET tokamak

The reflectometry technique and the subsequent tomographic inversion have recently been performed on the JET tokamak. As on Tore Supra, X-mode reflectometry is used, so that the density profile can be scanned up to the high field side, and profiles with dn e /dR > 0 can be observed if |dn e /dr| n e /a, where a is the minor radius. The results are presented in Fig. 1, which shows the density time trace and the corresponding tomographic 2D images of the density. The JET density profiles were reconstructed using the X-mode reflectometers [START_REF] Sirinelli | [END_REF]. JET reflectometers were swept simultaneously every 15 µs, this time covers the dwelt time between sweeps. 10000 profiles for a total duration of 150 ms were recorded in the ohmic phase before heating phase.

The impurity module in XTOR-2F

The XTOR-2F code solves the MHD equations with bifluid terms included in toroidal geometry. Bifluid diamagnetic terms are included in all the equations, in particular in Ohm law. The equilibrium is provided by the CHEASE code. [7]. The equations solved by the code for the following simulations are recalled in Ref. [2]. In order to simulate the density n Z of an impurity Z A X, a new equation (or as many as impurities) is added in the solver, with the following form: where V p = V p (r) e e e r is the pinch velocity, D Z (r) the diffusion coefficient, n Z0 (r) the initial density profile, and e e e r is the radial unit vector. Thus the only difference between the impurities lies in the choice of V p , D Z and n Z0 . The velocity v v v is the sum of the E × B drift and the parallel velocity v B/B of the plasma. Since ∇ • v v v = 0 the passive scalar assumption means there is no local conservation of particles. However (i) we checked that global particle conservation is not affected too much during our simulation (error is less than 1%) and (ii) a compressible version, including an additional equation for the parallel velocity of the impurity, is currently developed.

∂ n Z ∂t + v v v • ∇n Z + ∇ • (n Z -n Z0 ) V p = ∇ • D Z ∇ (n Z -n Z0 ) , k = 1, . . . , n imp , (1)
Thus current results are only preliminary.

The choice of the initial profile and of the transport coefficients is motivated by simplicity.

We aim at precisely assessing the proportion of impurities effectively ejected by the crash.

The easiest way to do this is by imposing a vanishing density outside the q = 1 surface, and monitoring the quantity

N Z (t) = V q=1 n Z (x x x,t)d 3 x x x V q=1 n Z (x x x,t 0 )d 3 x x x, (2)
where V q=1 is the volume contained inside the q = 1 surface, and t 0 is a time immediately before the crash. The results will be expressed in percentage for convenience. In the case where impurities are efficiently expelled from the q = 1 surface, N Z = 0 after the crash, while it has a finite value 0 < N Z < 1 in the opposite case. The impurity diffusion coefficient is ad hoc and chosen as to mimic the behaviour of critical gradient models, for example Ref. [8]. It increases by a factor of ten between the core and the edge, where it takes the value of the thermal diffusivity, D = χ ⊥ = 10 -6 in XTOR-2F units, normalized to a 2 /τ A , a the minor radius, τ A the Alfvén time.

40 th EPS Conference on Plasma Physics O2.108

Results

Fig. 2 shows the results obtained with this method. Panel (a) shows the initial and final fluxsurface averaged profiles, while panel (b) shows the evolution of N Z (t). The crash is seen to expel a significant part of the density, however the structures remaining inside the q = 1 surface immediately after the crash amount to approximately 15% of the total initial particle content.

The post-crash profile is not flat, with a quite large hump just outside the initial q = 1 surface, corresponding to the mixing radius of the sawtooth. At the end of reconnection (phase I), diffusion effects take over to increase N Z , and two slopes can be identified. In a first time (phase II) it increases quite rapidly and then more slowly (phase III). Indeed the hump diffuses both outward and inward, which contributes to increasing the particle content in the core and reducing the sawtooth efficiency. The smaller slope corresponds to the mere effect of the source.

The final value for N Z , once the redistribution due to the crash is over (end of phase II), is around 30%, meaning that 70% only of the particle content inside the q = 1 surface is expelled.

These results are in theory profile dependent, however similar behavior, both qualitatively and quantitatively, is recovered for all sorts of peaked profiles we have tested so far.

Conclusion

The confirmation of Tore Supra experimental results [2] on the JET tokamak and the physical analysis with XTOR-2F shows the relevance of the described structures in ohmic regime. The passive scalar study is very informative and shows that a significant part of the expelled particle content is reinjected just after the crash. This raises concern about the sawtooth flushing mechanism supposed to wash the core from impurities or He ash.

This work was supported by EURATOM and carried out within the framework of the European Fusion Development Agreement. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Figure 1 :

 1 Figure 1: Time trace of the density obtained in ohmic mode inside the q = 1 surface, JET pulse #83090 (a) and tomography images (b)

Figure 2 :

 2 Figure 2: Density before the crash and after the crash (a), and evolution of N Z , expressed as a percentage of the initial value (b).

th EPS Conference on Plasma Physics O2.108