L Isoardi

G Chiavassa

G Ciraolo

P Haldenwang

E Serre

Ph Ghendrih

Y Sarazin

F Schwander

P Tamain

Penalization modeling of a limiter in the tokamak edge plasma

Keywords: penalization, plasma, limiter, tokamak, Scrape Off Layer

An original penalization method is applied to model the interaction of magnetically confined plasma with limiter in the frame of a minimal transport model for ionic density and parallel momentum. The limiter is considered as a pure particle sink for the plasma and consequently the density and the momentum are enforced to be zero inside. Comparisons of the numerical results with one dimensional analytical solutions show a very good agreement. In particular, the penalization scheme followed in this paper tends to ensure an almost sonic plasma condition at the plasma-obstacle interface, Bohm-like criterion, with relatively weak dependence on the target Mach number profile within the obstacle. The new system being solved in a periodic obstacle free domain, an efficient pseudo-spectral algorithm based on a Fast Fourier transform is also proposed, and associated with an exponential filtering of the unphysical oscillations due to Gibbs phenomenon. Finally, the efficiency of the method is illustrated by investigating the flow spreading from the plasma core to the Scrape Off Layer at the wall in a two-dimensional system with one, then two neighboring limiters.

Introduction

Nuclear fusion occurs in magnetic confinement devices when hydrogen isotopes, namely deuterium and tritium, combine to form helium in a very hot (10 8 K) and low density (10 20 m -3 ) ionized gas or plasma. Progress toward the development of fusion as a secure and reliable energy source requires the accelerated development of computational tools and techniques that aid the acquisition of the scientific understanding needed to develop predictive and interpretative models [START_REF] Rosenbluth | Physics fundamentals for ITER[END_REF]. These are mandatory to address the physics of the plasma in particular near the boundary layers, the edge plasma, where the empirical scaling laws appear to be rather uncertain [START_REF]ITER Physics Expert Group on Divertor Modelling and Database[END_REF]. This might be related to the large number of dimensionless parameters that could govern the physics of this specific region, such as those related to atomic processes, or to the geometrical features of the plasma facing components.

The issue of plasma-wall interaction is a matter of growing concern since it associates a broad range of issues: safety, engineering of the wall components and reactor performance [START_REF] Loarte | Chapter 4: Power and particle control[END_REF]. Most important among the latter is the ability of future devices like the tokamak ITER to reach operation regimes with burning plasmas. Regarding the edge plasma, the issue is threefold; first, one must ensure the build-up of the H-mode transport barrier located in this edge region, second, one must control the energy out-flux from the plasma to the wall components so that it remains below the prescribed technological limit, and, third, one must comply with safety regulations regarding in particular the tritium wall inventory, the latter stemming from the plasmawall interaction. All these issues are closely related to the wall boundary conditions.

Conversely to neutral fluid flows, solid surfaces are sinks for charged particles so that the plasma boundary condition on the surface is a flow directed into the solid surface [START_REF] Tonks | A General Theory of the Plasma of an Arc[END_REF]. Indeed, the very large density within the solid compared to that of the plasma ensures the slowing down of the impinging ions within 10 -8 m as well as their prompt recombination into neutral atoms. This important difference between neutral fluids and plasmas has several key consequences. Regarding the plasma facing components, these are submitted to an intense particle and energy flux that must be controlled to ensure a long lifetime of the components [START_REF] Jacquinot | Steady-state operation of tokamaks: key experiments, integrated modelling and technology developments on[END_REF]. Regarding the plasma, a specific plasma boundary layer, the sheath, builds at the solid interface [START_REF] Riemann | Kinetic theory of the plasma sheath transition in weakly ionized plasma[END_REF]. The sheath extends over several Debye lengths, typically 10 -5 m, a scale much smaller than most other scales of the edge plasma. The sheath physics is governed by the departure from quasi-neutrality leading to a supersonic constraint and a positively charged plasma in the sheath. Asymptotic matching allows one to define boundary conditions at the sheath entrance, usually referred to as the Bohm boundary conditions . These constrain the particle flux to a sonic flow and determine the plasma energy flux lost by the plasma. The losses must be balanced by a particle and energy source in order to reach a dynamical steady state. Part of the particle refueling process is governed by neutral particle desorption from the plasma facing components providing therefore a volumetric plasma particle source via ionization. This source or recycling process further enforces the impact of the wall components geometry on the overall behavior of the edge plasma [START_REF] Tsitrone | Role of wall implantation of charge exchange neutrals in the deuterium retention for Tore Supra long discharges[END_REF].

The tokamak geometry, as exemplified by ITER, exhibits a toroidal geometry where the magnetic surfaces (surface tangent to the magnetic field) are organized as nested tori. The free and hence fast parallel motion along the magnetic field tends to ensure that most relevant physical quantities are nearly constant on the magnetic surfaces. This effect is less pronounced toward the edge due in particular to the particle sources and large scale plasma flows. Further in the edge, in the so-called Scrape-Off Layer (SOL), the wall components intercept the magnetic field lines. The geometry of the wall components and that of the magnetic field lines are then organized to control as best possible the plasma-wall interaction and in particular the plasma out-flux along the magnetic field lines into the wall components, as well as the neutral particle influx from these wall components. Two main configurations are implemented in present devices, the divertor geometry and the limiter geometry [START_REF] Stangeby | The plasma boundary of magnetic fusion devices[END_REF]. In both cases, the double periodicity that characterize the toroidal geometry is lost. Further complexity is introduced in these two configurations by secondary limiters. Limiters are protruding structures attached to the main vessel with a relatively small extent along the field lines, from 0.1 m to 1 m, compared to a standard field line length, from 10 m to 50 m. Properly taking into account such components is a challenge in the modeling effort. To address this issue we shall consider the case of a toroidally symmetric limiter, see a schematic view on Fig. 1a. This simplified geometry allows one to test and validate means to implement the appropriate Bohm boundary conditions without requiring surface boundary conditions.

Given the complexity of the physics of plasma-wall interaction, the macroscopic description of the plasma based on the three-dimensional fluid-like equations appears to be much more tractable numerically than the kinetic one. Simulations of large-scale motions and instabilities become feasible and one should be able to perform routine simulations prior to ITER start-up. Fluid-like equations are obtained directly by advancing velocity moments of the kinetic equation in time but require closures assuming high collisionality, strong fields and long times [START_REF] Braginskii | Transport processes in a plasma, Reviews in Plasma Physics[END_REF]. Usually, limiter or divertor geometries are taken into account in numerical codes by implementing the Bohm boundary conditions on the obstacle surface. Due to the very fast dynamics of the plasma at the limiter interface, such boundary conditions on the surface of the obstacle can require body-fitted unstructured meshes or remeshing strategies which can be time consuming.

In fluid mechanics, the volume penalization is an elegant and efficient way of introducing obstacles keeping cartesian grids. In the volume-penalization approach introduced by Arquis and Caltagirone [START_REF] Arquis | Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide-milieu poreux: application à la convection naturelle[END_REF] for incompressible flows, a Darcy drag term is added to the Navier-Stokes equations such that the velocity is penalized toward zero inside the obstacle in order to mimic the noslip boundary condition. The new system is handled now in an obstacle free domain for which fast solvers are available. Moreover, this technique allows an easy modification of the obstacle geometry and the associated boundary conditions.

Many different numerical simulations using penalization technique have put in evidence its efficiency for incompressible flows governed by Navier-Stokes equations [START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF][START_REF] Kevlahan | Computation of turbulent flow past an array of cylinders using a spectral method with Brinkmann penalization[END_REF][START_REF] Kevlahan | An Adaptive Wavelet collocation Method for Fluid-Structure Interaction at High Reynolds Number[END_REF][START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF][START_REF] Schneider | Adaptive wavelet simulation of a flow around an impulsively started cylinder using penalization[END_REF][START_REF] Keetels | Fourier spectral and wavelet solvers for the incompressible Navier-Stokes equations with volume penalization: Convergence of a dipole-wall collision[END_REF][START_REF] Vasilyev | Hybrid wavelet collocation-Brinkman penalization method for complex geometry flows[END_REF]. In this context, different discretization methods, such as wavelets [START_REF] Schneider | Adaptive wavelet simulation of a flow around an impulsively started cylinder using penalization[END_REF][START_REF] Kevlahan | An Adaptive Wavelet collocation Method for Fluid-Structure Interaction at High Reynolds Number[END_REF][START_REF] Vasilyev | Hybrid wavelet collocation-Brinkman penalization method for complex geometry flows[END_REF], pseudospectral [START_REF] Kevlahan | Computation of turbulent flow past an array of cylinders using a spectral method with Brinkmann penalization[END_REF][START_REF] Keetels | Fourier spectral and wavelet solvers for the incompressible Navier-Stokes equations with volume penalization: Convergence of a dipole-wall collision[END_REF], or finite difference/volume [START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF][START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF], have been successfully tested and the penalized system has been shown to converge toward the Navier-Stokes equations with no-slip boundary conditions when increasing the penalization term [START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF]. In [START_REF] Paccou | A penalization method applied to the wave equation[END_REF], this technique has been applied and analyzed for the one dimensional linear wave equation and recently a different approach has been proposed for linear and semi-linear symmetric hyperbolic problems [START_REF] Fornet | A Kreiss' symmetrizer domain penalization approach[END_REF][START_REF] Fornet | Penalization approach to semi-linear symmetric hyperbolic problems with dissipative boundary conditions[END_REF]. Penalization methods have also been developed for the compressible Navier-Stokes system, for low Mach numbers flows [START_REF] Cho | Immersed Boundary Method for Viscous Flow Around Moving Bodies[END_REF] acoustics [START_REF] Liu | A Brinkman penalization method for compressible flow in complex geometries[END_REF] and for supersonic shocked flows where convective terms are dominant in [START_REF] Boiron | A High-Resolution Penalization Method for large Mach number Flows in the presence of obstacles[END_REF].

In the present work, the limiter within a tokamak with circular crosssection, Fig. 1a is modeled using an original penalization technique. The method is implemented in a reduced model of transport equations for the ion density n i and the particle flux Γ // = n i v // in the direction parallel to the magnetic field. The plasma region includes both the edge region with closed magnetic field lines and the SOL region, where field lines are intercepted by wall components as shown on Fig. 1a. In this model, the convective parallel transport, governed by nonlinear hyperbolic equations is dominant while the perpendicular transport is only diffusive.

The paper is organized as follows. The physical setting including description of the geometry of the problem and governing equations are presented in Section 2. In Section 3, the penalization method is described and analyzed (part 3.1). The penalized one dimensional nonlinear hyperbolic system is discretized and solved in Section 3.2 with a high order shock capturing finite difference scheme. Numerical results are compared to a reference exact solution in order to validate the proposed penalization method. The following Section, Section 4, addresses the suitability of a Fourier pseudo spectral method for the one-dimensional penalized system. In Section 5, two dimensional simulations are presented and compared to reference analytical and numerical solutions of the tokamak edge plasma. Conclusions and perspectives are drawn in Section 6.

Physical setting: Geometry and governing equations

The two-dimensional geometric model Fig. 1c represents the transition region from the periodic core to the non periodic Scrape Off Layer (SOL) of a tokamak with toroidal limiter. While the limiter geometry used here is a sketch of actual limiter geometries [START_REF] Jacquinot | Steady-state operation of tokamaks: key experiments, integrated modelling and technology developments on[END_REF][START_REF] Stangeby | The plasma boundary of magnetic fusion devices[END_REF], this geometry is sufficient to address the main issues related to the insertion of a limiter in the plasma, with the exception of the specific issues raised by shallow incidence angles. The dimensions of the domain are ∆ r and L , respectively in the radial and parallel directions, that are used to normalize the two directions of the model. r stands for the radial direction normalized by ∆ r in a case where the poloidal section of the magnetic surfaces are circular (Fig. 1a). The parallel direction s is aligned on the magnetic field line and normalized by L , namely the connection length of the magnetic field line between two intersections with the limiter. The aspect ratio of the simulation domain is defined by A = L /∆ r . A limiter of parallel extension ℓ yields two further dimensionless parameters, the normalized parallel extension ∆ s = ℓ /L and the limiter height ∆ r /2. The latter defines the relative weight of the SOL and edge plasmas in the simulation domain. For symmetry reasons, this limiter is located at the middle of domain in the parallel direction at s = 0.5. Due to the very strong recombination rate of the plasma within this solid state component, this limiter causes loss of periodicity that characterizes the SOL plasmas. A minimum system, involving the ion density n i and the parallel particle flux (or a parallel momentum) Γ = n i u has been considered to investigate the interaction between the SOL and the edge. The equations deduced from the balance equations for compressible flows will represent the plasma equations when adding electrons and ions in the quasineutral limit and neglecting the electron inertia terms :

   ∂ τ i ∂t n i + ∂ L ∂s Γ -∂ ∆r∂r D ∂ ∆r∂r n i = 0 ∂ τ i ∂t Γ + ∂ L ∂s Π -∂ ∆r∂r ν ∂ ∆r∂r Γ = 0. (1)
In this set of equation the time is normalized to an ion characteristic time τ i since the electron inertia terms are not taken into account. The flow is assumed isothermal with constant ion and electron temperatures T i and T e involving a constant sound speed c s , such that m i c 2 s = (T e + T i ) where m i is the ion mass. The quantity Π = m i (Γ 2 /n i + c 2 s n i ) is the total plasma pressure (kinetic plus thermodynamic and electrons plus ions). A cross-field particle diffusion term D has been added in the density equation. Similarly, the momentum diffusion term, with coefficient ν, is also anisotropic and only taken into account in the radial direction. These diffusive transport mechanisms are introduced to take into account the effect of turbulent fluctuations in the cross-field direction, this turbulent transport being significantly larger than the collisional one. In the present work, D and ν will be assumed time-independent and spatially uniform.

Boundary conditions are periodic in the parallel direction s. In the radial direction r, the particle flux coming from the core, at r = 0, has been imposed through a constant density gradient while the parallel momentum is assumed equal to zero at the core interface. On the tokamak wall, at r = 1, radial absorption conditions (zero fluxes) for both density and parallel momentum have been imposed. On the limiter, the physics of the sheath is taken into account in the parallel direction via the Bohm boundary conditions [START_REF] Stangeby | The plasma boundary of magnetic fusion devices[END_REF]. In the simplest and usual formulation, these conditions are met by imposing a constant velocity equal to ±c s on the two opposite sides of the limiter. The acoustic velocity c s is the reference velocity so that the normalized velocity is given by the Mach number M. The density n i is normalized by a constant and a-priori arbitrary density n 0 , hence defining the normalized density N = n i /n 0 . Finally, the parallel momentum is normalized by the product n 0 c s , leading one to define Γ = Γ /(n 0 c s ). The parallel transport being convective tends to generate small time scales. It is therefore convenient to use τ = L /c s to normalize time, hence setting τ i = τ . The density and momentum evolution equations, Eq.( 1), then read for N and Γ :

∂ t N + ∂ s Γ = A P e ∂ 2 r N ∂ t Γ + ∂ s Γ 2 N + N = A Sc P e ∂ 2 r Γ, (2)
where (r, s) ∈ [0, 1] × [0, 1]. The flow is controlled by the aspect ratio A = L /∆ r , the ratio between the momentum and particle diffusivities, ν/D, and the ratio between the convective and diffusive transports, ∆ r c s /D. By analogy with fluid mechanics the two latter control parameters are named the Schmidt number Sc and the Peclet number P e respectively. Eq.( 2) is characterized by a strong anisotropy in the flow dynamics with a rapid transport by convection in the parallel direction and a slow transport by diffusion in the radial direction. Key features of the rapid transport in the parallel direction, and in particular the penalization issues addressed in this paper, will be identified by analyzing the following one-dimensional problem:

∂ t N + ∂ s Γ = S ∂ t Γ + ∂ s Γ 2 N + N = 0. ( 3
)
where S is a time-independent source term added to ensure the existence of non trivial steady state solutions, and the momentum Γ = NM. This system can be deduced from Eq.( 2) by setting to zero the diffusion coefficients. It defines a nonlinear, 1 dimension, hyperbolic problem.

Steady state solutions must verify the following equations for N and M.

∂ s N = - 2M 1 -M 2 S ; ∂ s M = (1 + M 2 ) (1 -M 2 ) S N . (4)
As for any fluid, this equation states that in a generic fashion an infinite gradients build-up at M 2 = 1, and prevent the transition from subsonic to supersonic flows. When the gradient lengths reach values comparable to λ D /L ≈ 10 -6 , the ratio of the Debye length by the connection length , the scale separation allowing one to consider a quasi-neutral plasma is not met anymore and the full Poisson equation should be considered. However, this only occurs when M 2 ≈ 1 -λ D /L . Prior to this sheath condition, hence for quasi-neutral plasmas, strong parallel gradients will appear at a distance from the wall that can be significantly larger than the Debye length. When addressing simulations in the vicinity of the boundary, the very sharp density gradients will thus require special numerical treatment to avoid spurious oscillations. It is important to stress that while a precise computation of the density N and Mach number M near the boundary will be demanding, the calculation of either the plasma particle outflux Γ or the total plasma pressure Π do not exhibit these strong gradients and should be more accurately computed.

When the source term S is constant along the field line, Eq.( 3) with Bohm conditions at the limiter (M 2 = 1) admits a steady state solution for the density and parallel momentum. This solution is given in appendix A and will be used in the numerical test of the penalization method.

The penalization method

The basic idea of all penalization techniques for incompressible flow was to enforce the no-slip boundary conditions V = 0, where V is the velocity field, by adding a penalized velocity term in the momentum equation. In [START_REF] Arquis | Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide-milieu poreux: application à la convection naturelle[END_REF][START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF], the penalization term 1 η χ * (V ) T , where χ is the characteristic (or mask) function of the solid domain Ω and η the penalization parameter, is added to the momentum equation of the incompressible Navier-Stokes system. This technique corresponds to a Brinkman-type porous media model with variable permeability, where the fluid domain has a very large permeability in front of the obstacle. The analysis carried out in [START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF] demonstrates that the solution of the penalized system converges to the solution of the incompressible Navier-Stokes system in the fluid domain, and that the velocity converges to zero in the solid domain at a theoretical rate of O(η 3/4 ). Extensions to compressible flows have been proposed in [START_REF] Cho | Immersed Boundary Method for Viscous Flow Around Moving Bodies[END_REF][START_REF] Liu | A Brinkman penalization method for compressible flow in complex geometries[END_REF][START_REF] Boiron | A High-Resolution Penalization Method for large Mach number Flows in the presence of obstacles[END_REF] using the same underlying idea.

In the case of purely hyperbolic problems, penalization techniques have only been developed for linear one dimensional equation [START_REF] Paccou | A penalization method applied to the wave equation[END_REF] or linear and semi-linear systems [START_REF] Fornet | A Kreiss' symmetrizer domain penalization approach[END_REF][START_REF] Fornet | Penalization approach to semi-linear symmetric hyperbolic problems with dissipative boundary conditions[END_REF]. To our knowledge, no such technique has yet been defined for nonlinear equations, and we propose in section 3.1 an original volume penalization for the system (1).

Description and analysis

The limiter is modeled as a particles sink using a penalized system for Eqs. (2). In such system, the plasma-limiter interaction is no longer described by demanding to density and velocity to satisfy Bohm conditions. The limiter (or obstacle) Ω is represented by a characteristic function χ(s, r) satisfying χ(s, r) = 1 if (s, r) ∈ Ω and χ(s, r) = 0 if (s, r) ∈ Ω. Density and momentum depend now on the penalization parameter η which must be small (η << 1) in the following system :

∂ t N η + ∂ s Γ η + χ η N η = A P e ∂ 2 r N η ∂ t Γ η + (1 -χ) ∂ s Γ 2 η Nη + N η + χ η (Γ η -Γ Ω ) = ASc P e ∂ 2 r Γ η . (5)
Γ Ω is the 'momentum in the obstacle'. Its value is given in term of a fixed Mach number

M Ω by Γ Ω = N η M Ω . (6)
M Ω (s, r) is imposed inside the limiter (non identically zero) and is zero outside. More details will be given on practical choices of M Ω in section 3.2. The pressure force is canceled within the limiter in order to avoid a divergence of the Mach number induced by the vanishing density.

It can be shown, as is done in [START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF], that the penalized terms added in Eq.( 5) enforce the initial boundary conditions on density and momentum by writing formally that

N η = N + η Ñ, M η = M + η M . (7)
It then follows that Γ η = Γ + η Γ to first order in η. By introducing these relations in Eq.( 5)and by identifying the terms of order η -1 one finds:

χN = 0, χ(Γ -NM Ω ) = 0, (8)
while for the η 0 terms one obtains:

∂ t N + ∂ s Γ + χ Ñ = A P e ∂ 2 r N ∂ t Γ + (1 -χ) ∂ s Γ 2 N + N + χ( Γ -Ñ M Ω ) = ASc P e ∂ 2 r Γ. (9)
In the plasma domain, the mask function χ is zero and Eq.( 9) reduces to the initial system of Eq.( 2). Inside the limiter, χ = 1, and hence Eq.( 8) becomes:

N = 0, Γ = NM Ω , (s, r) ∈ Ω ( 10
)
which yields immediately N = 0 and Γ = 0 in the limiter, i.e, the physical plasma conditions in the limiter. Thus, the penalized system Eq.( 5) enforces, at least formally, the physical boundary conditions on density and momentum in the limiter, while solving the initial system of equations in the plasma region. These conditions could also be obtained in a simpler way using directly Γ Ω = 0 in Eq.( 5). But, as proven in Section 3.2, it is essential for computational reasons to control the ratio Γ/N inside the limiter, since this quantity will be directly involved in the time step stability criterion. This is the role played by the fixed function M Ω since inside the limiter N and Γ will tend to zero but with a bounded ratio:

Γ N = M Ω (s, r) ∈ Ω. ( 11
)
Convergence of the penalized solution N η , Γ η toward the solution of Eq.( 2) when η tends to zero is examined in numerical simulations in the following section.

Numerical validation

The method has been validated in the one-dimensional case by comparing numerical results to the exact analytical solution developped in Appendix (A). The penalized system is deduced from Eq.( 3) and reads:

     ∂ t N + ∂ s Γ + χ η N = (1 -χ)S ∂ t Γ + (1 -χ) ∂ s Γ 2 N + N + χ η (Γ -Γ Ω ) = 0 s ∈ [0, 1], t ∈ ]0, +∞[. (12)
The subscript η has been removed for simplicity. The mask function χ is defined such as χ = 1 for s ∈]0.5(1 -∆ s ), 0.5(1 + ∆ s )[. Thus, the values of the unknowns at the boundaries s = 0.5(1 ± ∆ s ) are computed with the equations and not imposed directly by the penalization method. Periodic boundary conditions are considered at s = 0 and s = 1. Since the steady state solution of Eq.( 1) exhibits strong gradients in the vicinity of the obstacle, the same behavior must be recovered with Eq.( 12) so that we are led to use a robust shock-capturing finite difference method in order to avoid any spurious numerical oscillations. Considering that the equations are now written in an obstacle free domain [0, 1], a uniform discretization {s i = i∆s, i = 0, ..., n} of this interval is used. In what follows, N i and Γ i stand for the grid point computed values N(s i ) and Γ(s i ), respectively. Eq.( 12) has been discretized by a classical method of lines, which leads one to the following system of ODEs at each point of the computational grid:

dN i dt + f i+1/2 -f i-1/2 ∆s + χ(s i ) η N i = (1 -χ(s i ))S i dΓ i dt + (1 -χ(s i )) g i+1/2 -g i-1/2 ∆s + χ(s i ) η (Γ i -Γ Ω,i ) = 0, (13)
f j , g j are the numerical flux functions, consistent with the convective physical flux functions of Eq.( 12). Evaluation of these numerical fluxes is obtained using a shock capturing scheme, based on the Roe approximate Riemann solver [START_REF] Leveque | Finite Volume Methods for Hyperbolic Problems[END_REF], associated with a robust third order ENO spatial reconstruction [START_REF] Shu | Efficient Implementation of Essentially Non-Oscillatory Shock Capturing Schemes II[END_REF].

The time integration of Eq.( 13) is carried out by using a second order TVD Runge-Kutta scheme described in [START_REF] Shu | Efficient Implementation of Essentially Non-Oscillatory Shock Capturing Schemes II[END_REF]. The convective terms are treated in an explicit fashion, while the penalization terms are treated implicitly. In that way the time step restriction is independent of the penalization parameter η and depends only on the eigenvalues of the Jacobian matrix of the convective flux F = (Γ, Γ 2 N + N) T . Since eigenvalues are equals to M ± 1, the stability time step condition for the full discretized scheme reads:

∆t = C ∆s max s∈[0,1] (|M(s)| + 1) . ( 14
)
The CFL constant is set in the program to C = 0.8 and the time step is adapted automatically at each time step using relation [START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF]. In this relation, we remark that the time step depends on the Mach number M in the obstacle. Controlling the latter then provides a means to optimize the numerical scheme. This justifies the introduction of the function M Ω in [START_REF] Riemann | Kinetic theory of the plasma sheath transition in weakly ionized plasma[END_REF]. It is also important to notice that the implicit treatment of the penalization terms will not introduce particular numerical issues, since their specific forms will make the linear system diagonal. In all numerical results, the computations has been started from the initial conditions N(s, t = 0) = 1 and Γ(s, t = 0) = 0, and the steady state has been reached in order to compare the exact solution given in Eqs. [START_REF] Shu | Efficient Implementation of Essentially Non-Oscillatory Shock Capturing Schemes II[END_REF][START_REF] Serre | 3D pseudo-spectral method for convection in rotating cylinder[END_REF] to the numerical solution. The limiter size in the parallel direction is set at ∆ s = 0.1. The particle source term S is kept constant and equal to 2. The number of grid points n is fixed, n = 100.

Numerical results for density and momentum are presented on Fig. 2. Excellent agreement with the analytical solution can be observed. Inside the limiter, density and momentum are driven to zero, thus following the prescribed behavior (see Fig. 3a). The ratio M = Γ/N, displayed on Fig. 3b perfectly matches M Ω , as predicted by the analysis in Section 3.1, Eq.( 11). The discrete L 1 errors are defined as:

E 1 = N -N exact 1 = 1 n n-1 i=0 |N(s i ) -N exact (s i )|, (15)
E 1,Ω = N -N exact 1,Ω = 1 n Ω s i ∈Ω |N(s i ) -N exact (s i )|. ( 16
)
n Ω is the number of discretization points within the limiter and N exact the exact density computed by Eq.( 26).

Due to the presence of strong gradients, the measured order of convergence for the global error (Eq.( 15)) is 0.94. Nevertheless, our numerical experiments show that the present high order scheme is much more accurate in smooth zones of the solution and less dissipative near discontinuities than a first order scheme.

The influence of the penalization parameter η on the accuracy of the numerical scheme, Eq.( 13), is illustrated on Fig. 4. Inside the limiter, the local error defined in Eq.( 16) decays linearly with η as reported in [START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF][START_REF] Boiron | A High-Resolution Penalization Method for large Mach number Flows in the presence of obstacles[END_REF][START_REF] Kevlahan | Computation of turbulent flow past an array of cylinders using a spectral method with Brinkmann penalization[END_REF][START_REF] Kevlahan | An Adaptive Wavelet collocation Method for Fluid-Structure Interaction at High Reynolds Number[END_REF][START_REF] Paccou | A penalization method applied to the wave equation[END_REF][START_REF] Schneider | Adaptive wavelet simulation of a flow around an impulsively started cylinder using penalization[END_REF]. Such behavior gives a real control on the way the boundary conditions are imposed numerically. The global error (Eq.( 15)) also decays and levels off as it reaches the accuracy of the numerical scheme that means that the global error becomes independent of η as this parameter decreases. This interesting behavior is typical of penalization methods and is due to the presence of a thin boundary layer whose size is proportional to √ η (see references [START_REF] Kevlahan | Computation of turbulent flow past an array of cylinders using a spectral method with Brinkmann penalization[END_REF][START_REF] Paccou | A penalization method applied to the wave equation[END_REF][START_REF] Fornet | Penalization approach to semi-linear symmetric hyperbolic problems with dissipative boundary conditions[END_REF][START_REF] Pasquetti | A pseudo-penalization method for high Reynolds number unsteady flows[END_REF]). Once this value is small in comparison to the grid spacing ∆s, the boundary layer is not resolved and the expected order of the numerical scheme is recovered. For the present non linear system [START_REF] Kevlahan | Computation of turbulent flow past an array of cylinders using a spectral method with Brinkmann penalization[END_REF], no theoretical result exists but our numerical experiments confirm the same tendency.

These results indicate that the choice of the function M Ω does not appear to influence significantly the behavior of the solution in the plasma. In particular, the Mach number at the obstacle boundaries still tends to ±1 although one has enforced unrealistic values M Ω (s ∓ ∆ s /2) = ±2 or M Ω (s ∓ ∆ s /2) = ±0.5 within the limiter volume, as readily observed on Fig. 5. The boundary conditions imposed by the penalization method (N = 0 and Γ = 0) thus allow one to recover Bohm-like conditions with a relatively reduced sensitivity to the chosen function M Ω . This feature is presently being analyzed and appears to be related to a minimization of the plasma pressure governed by the particle sink at the obstacle. The issue at hand is not directly related to the penalization technique presented here and will be fully analyzed in another paper. Furthermore, the stability time step (Eq. ( 14)) depends on the maximum Mach number over the whole domain, expressions of M Ω with maximum values close to one are preferable.

One dimensional pseudo-spectral method

Periodicity in the parallel direction is recovered in the SOL with our penalization method so that one is naturally led to question the suitability of the Fourier approximation in this renewed framework. In many situations (see for example [START_REF] Serre | 3D pseudo-spectral method for convection in rotating cylinder[END_REF][START_REF] Severac | A Spectral Viscosity LES for the simulation of turbulent flows within rotating cavities[END_REF]), Fourier series are well known for their efficiency both in terms of accuracy and in terms of performances. They are based on efficient FFT (Fast Fourier Transform) algorithms which can be easily parallelized and efficiently run on supercomputers. Spectral methods nondispersiveness and non-dissipativity should make them well fitted to computations of plasma turbulence requiring fine scales resolution and long time integration. However, their implementation to solve Eq.( 12) is not straightforward. The uniform convergence is no longer exponential and localized discontinuities make the Fourier projection suffer from oscillatory behavior known as the Gibbs phenomenon. Discontinuities are introduced by the mask function χ and steep gradients that occur in the vicinity of the limiter where the parallel density gradient goes to infinity (Eq.( 4)) but also inside the limiter where density and momentum are driven to zero on a very short distance by the penalization technique. Moreover, the hyperbolic character of the equations in the parallel direction can propagate these oscillations. Thus, from a practical point of view, we have to address here the problem of the computation stability, i.e. if these Gibbs oscillations remain bounded during the calculations, and the question of the global accuracy, i.e. knowing that these Gibbs oscillations are first order only, is there a means to recover high-order information. For incompressible flows, Kevlahan and Ghidaglia [START_REF] Kevlahan | Computation of turbulent flow past an array of cylinders using a spectral method with Brinkmann penalization[END_REF] showed that Gibbs oscillations remained bounded and that the computation remains stable during the flow evolution. Moreover, recent advances in the theory (Gelb and Gottlieb [START_REF] Gelb | The Resolution of the Gibbs Phenomenon for Fourier Spectral Methods[END_REF]) allow one to perform stable computations of discontinuous solutions (shocks waves) of hyperbolic partial differential equation and indicate that high-order information can be recovered. The 1D-system Eq.( 12) where the periodicity has been recovered by using the penalization technique is considered. The time integration is computed using a classical second order Runge-Kutta scheme where the convective terms are treated explicitly and the penalization terms treated implicitly as previously. The time step restriction is independent of the penalization parameter η and is given by Eq.( 14). The CFL constant is also set equal to C = 0.8, the source term S, the velocity M Ω inside the limiter and the limiter extent ∆ s are left unchanged.

The system is discretized in physical space and the nonlinear convective term is evaluated using the classical pseudo-spectral technique which consists of performing differentiations in spectral space and the product in the physical space [START_REF] Peyret | Spectral Methods for Incompressible Viscous Flow[END_REF]. The link between the two spaces is made by the FFT algorithm from the JMFFT library developed by the CNRS computing center IDRIS.

As expected by theory, unbounded Gibbs oscillations quickly appear inside and around the limiter and propagate in the whole computational domain leading to unstable computations regardless of the grid resolution up to n = 4096 Fourier modes. Following the work for hyperbolic problems of Gottlieb and Gottlieb in [START_REF] Gottlieb | Spectral methods for compressible reactive flows[END_REF], an exponential filtering has been implemented. If {a k , b k , -n/2 ≤ k ≤ n/2} are respectively the Fourier coefficients of U = (N, Γ) T and F = (Γ, Γ 2 /N + N) T , we write at each time step :

a tn+1 k = exp -α( k N ) 2p (a tn k -∆t b tn k ). ( 17
)
The filter is characterized by its order, here p = 2, and by a coefficient α which is tuned in order to provide the best compromise between stability and accuracy of the computations.

Numerical results show a global agreement with the analytical solution over the whole domain as presented in Fig. 6, density and momentum being driven to zero within the limiter. As previously, with a finite difference approximation, Bohm conditions M 2 = 1 are obtained at the limiter.

Nevertheless, weak Gibbs oscillations in the Γ profile remain, localized around the limiter despite the filtering (Fig. 6b). These oscillations remain bounded during the computation but the solution remains too far from the analytical solution. The use of a stronger filtering reduces the amplitude of the oscillations but a significant fraction of the wavenumber spectrum of the steady state solution is also damped, modifying in this way the general behavior of the solution. In particular, too large a filtering modifies the numerical representation of the limiter geometry. Even with a procedure for reconstruction in postprocessing to recover the pointwise convergence for the Fourier projections of piecewise continuous data, as proposed by [START_REF] Keetels | Fourier spectral and wavelet solvers for the incompressible Navier-Stokes equations with volume penalization: Convergence of a dipole-wall collision[END_REF], one is not able to obtain a better accuracy than the finite difference scheme (Section 3) in this situation. Moreover, the sensitivity of the solution to the tuning parameter α is very high and the choice of an optimal value may be tricky, also depending on the grid resolution. In conclusion, spectral method cannot be retained using the procedure outlined above.

The choice was therefore made to consider a smoothed mask function χ for the limiter, Fig. 7. This procedure is classical when combining spectral and penalization techniques [START_REF] Minguez | High-order LES of flow over the Ahmed reference body[END_REF][START_REF] Keetels | Fourier spectral and wavelet solvers for the incompressible Navier-Stokes equations with volume penalization: Convergence of a dipole-wall collision[END_REF][START_REF] Pasquetti | A pseudo-penalization method for high Reynolds number unsteady flows[END_REF]. Using the filtering procedure of Eq.( 17), the results exhibit similar features as that of the previous solutions obtained with the discontinuous mask function, but gradients are now smaller inside the limiter and the solutions are free of Gibbs oscillations, Fig. 8. The solution still satisfies the Bohm condition despite the use of this smooth shape limiter. Note that the location of the near sonic flow is known within the uncertainty of the smooth limiter location, of the order of several grid points, to be compared to the expected uncertainty stemming from theoretical considerations and estimated to be of the order of several Debye lengths. In practice such an error will govern an error on the particle outflux that will scale like 0.1ρ * where ρ * is the ratio between the characteristic scale of turbulence and the size of the device. For simulations of interest this error is then both small, of the order of 10 -3 or smaller, and small regarding the expected errors stemming from the magnetic equilibrium.

Using Eqs. [START_REF] Schneider | Adaptive wavelet simulation of a flow around an impulsively started cylinder using penalization[END_REF][START_REF] Keetels | Fourier spectral and wavelet solvers for the incompressible Navier-Stokes equations with volume penalization: Convergence of a dipole-wall collision[END_REF] expected by theory, and previously observed using a finite difference approximation, Fig. 4, the decrease of the magnitude of the error is linear with η inside the limiter. Using such a smoothed characteristic function for the limiter, one cannot use the analytical solution calculated in Appendix(A) for validation and consequently it will be impossible to determine the exact order of the scheme. The quality of the Fourier computations have been then investigated by calculating the truncation error with respect to the highest available resolution n max = 4096 considered as a benchmark computation. For all the computations, the penalization parameter and the filtering coefficient are fixed at η = 10 -7 and α = 10, respectively. The absolute value of the difference between the benchmark computation and computed solutions at lower resolution has been plotted on Fig. 10. As expected, the results are characterized by a decrease of the error when increasing the resolution and when considering the solution further away from the limiter even though the error behaves slightly differently for density and parallel momentum. For both fields, the maximum error is reached close to the limiter where the gradients are the steepest. The following global relative error for density and parallel momentum have been calculated over the whole domain Ω and plotted on Fig. 9b In all computations, the penalization parameter η = 10 -7 and the filtering coefficient α = 10 are used.

δ N (Ω) = N n -N nmax L 2 (Ω) N nmax L ∞ (Ω) (18)
δ Γ (Ω) = Γ n -Γ nmax L 2 (Ω) Γ nmax L∞(Ω) . (19)
As already noticed in section 2, the momentum is more accurately computed than the density due to the very sharp density gradients in the vicinity of the boundaries. The error decays faster than n -1 and n -3 for the density and the parallel momentum, respectively.

Two-dimensional results

In the previous section we have shown the suitability of the Fourier approximation in the case of the one-dimensional system [START_REF] Kevlahan | Computation of turbulent flow past an array of cylinders using a spectral method with Brinkmann penalization[END_REF]. Here, the approximation is extended to the two-dimensional system (5) in a computational domain of aspect ratio A = 50. The parallel s and radial r directions are discretized with n = 128 Fourier modes and with a second-order finitedifference scheme with n r = 50 grid points, respectively. A second-order Runge-Kutta scheme is used for the time integration and the diffusive radial terms are discretized implicitly for numerical stability reasons. For all computations, the penalization terms are still treated implicitly and the time step is defined as in the one-dimensional study. Boundary conditions at the tokamak wall and at the core have been defined in section 2. Numerical results are first compared with our former results [START_REF] Isoardi | Modelling SOL flow pattern spreading in the edge plasma[END_REF] obtained in the same geometry, Fig. 1b, with a limiter of zero thickness modeled by the Bohm conditions. Then, results are presented in a geometry with two neighboring limiters illustrating the potential of the penalization technique to easily cope with changes in the geometry.

Comparison with results obtained with a limiter modeled by surface boundary conditions

Former results [START_REF] Isoardi | Modelling SOL flow pattern spreading in the edge plasma[END_REF] were obtained by discretizing Eq.( 2) with a TVD second-order finite difference scheme. The same physics is captured here, Fig. 11. Two transport regimes are found in steady state depending on the relative magnitude of the diffusive transport and parallel momentum. These indicate that the SOL and the edge region may be more or less connected. At large diffusivity (Figs. (11a,b)), i.e. small Peclet number, the edge region and the SOL are strongly connected. Density from the edge can reach out to the tokamak wall. The parallel momentum increases toward the limiter up to a value corresponding to M = ±1 at the limiter as expected from the Bohm conditions. At low diffusivity (Figs. (11c,d)), i.e. large Peclet number, the solution is characterized by a strong parallel dynamics so that the transition from the close field lines to the SOL is rather abrupt. As in reference [START_REF] Isoardi | Modelling SOL flow pattern spreading in the edge plasma[END_REF], the parallel momentum can nevertheless spread from the SOL to the edge plasma leading to imprints of non-zero velocities in the closed field line region with periodic conditions. In the regime at large Peclet number, which corresponds to most situations met in tokamak plasmas, systematic quantitative comparisons with reference [START_REF] Isoardi | Modelling SOL flow pattern spreading in the edge plasma[END_REF] have been made. An important point for comparison is the steady state profile of the density averaged in the parallel direction < N > . In the edge plasma, the equation is a diffusive equation without loss of thus yielding a linear profile, the slope being imposed by the core boundary condition. Assuming that the density at the sheath is proportional to the averaged density, the expected profile in the SOL is an exponential with e-folding length λ n . In the limit Sc → 0, and provided the mean density < N > is equal to the density of the stagnation point (M = 0, s = 0 in the present simulation), one readily recovers the standard SOL width for λ n , namely λ n = (2A SOL /P e) 1/2 ∆ r , where A SOL = A(1 -∆ s ). In this expression one thus takes into account the limiter parallel extension. Given the expression of P e one finds the usual expression λ n = (D L SOL /c s ) 1/2 , with L SOL = L -ℓ . The comparison is illustrated on Fig. 12a. The result from [START_REF] Isoardi | Modelling SOL flow pattern spreading in the edge plasma[END_REF] with surface boundary condition (plain line), λ n /∆ r = 0.141, exhibits the appropriate slope and the result using the penalization technique (dashed line) allows one to recover this basic behavior, λ n /∆ r = 0.134. The small difference in λ n is governed by the limiter size. A slight deviation can however be observed between the two numerical solutions. It appears to be mainly governed by specific physics induced by the parallel extent of the limiter, vanishing limit ∆ s → 0 in [START_REF] Isoardi | Modelling SOL flow pattern spreading in the edge plasma[END_REF], ∆ s = 0.1 for the solution with penalization. This effect is more pronounced when considering the parallel variation of the total pressure, i.e. ∆Π/Π = Π lim. /Π stag. -1, where Π lim. is the total pressure at the limiter and Π stag. is the total pressure at the stagnation point (s = 0, M = 0), Fig. 12b. This quantity illustrates the pressure variation along a field line and thus accounts for the effect of the diffusive term in the momentum balance equation. The overall radial profile of this pressure variation is the same for both simulations. The momentum source, region with negative pressure variation is localized close to the limiter head, both in the SOL and edge plasma. In the edge plasma, toward the core plasma, the pressure variation is positive indicating a gain in momentum that results from a transfer from the SOL flow build-up. This momentum sink region rapidly decays away, essentially because this region has a higher pressure (higher density) and small Mach number. Conversely, toward the outer SOL the pressure variation is positive and increases because this region is a low pressure region. The diffusive momentum transfer from the higher density SOL region thus has more impact. It is important to underline here that the overall pressure variation is small, typically 5%, indicating that the pressure redistribution, proportional to the Schmidt number, is a weak effect in this diffusive model. The model with boundary conditions (plain line) and penalization (dashed line) shown on Fig. 12b also exhibit slight differences. First, the region of momentum source is more clearly located in the SOL for the surface boundary condition. This effect stems from the parallel extent of the limiter. Indeed, for a limiter with vanishing extent, the edge plasma close to the limiter head must exhibit a cross-over from a near M = 1 branch to a near M = -1 solution due to coupling to the neighboring SOL conditions at the limiter head. This hinders the development of the region with significant flow in the plasma edge. Also, when the limiter has finite parallel extent, there is a density loss term governed by the transverse diffusion to the limiter head that must be accounted for. This localized particle loss then governs a parallel flow to compensate it. The overall effect is to peak the region of momentum loss at the separatrix between the SOL and edge plasma. For the same reasons, the region with momentum sink in the edge plasma extends slightly more toward the core. The outer SOL region is weakly affected. On the overall, although some differences can be observed, they remain small and essentially governed by the difference in the models, thus validating the penalization method. [START_REF] Isoardi | Modelling SOL flow pattern spreading in the edge plasma[END_REF] with a limiter with zero thickness imposed by surface Bohm boundary conditions (solid line), the present penalization (dashed line) and the analytical slope [START_REF] Isoardi | Modelling SOL flow pattern spreading in the edge plasma[END_REF] (dotted line). The interaction between the edge and the SOL is for r = 0.5. Fourier modes n = 128 and n r = 50, penalization parameter η = 10 -7 and filtering coefficient α = 10 are used.

Flow features with two neighboring limiters

One of the advantages of penalization methods is their ability to deal with geometries of arbitrary complexity. This is achieved by changing the characteristic χ function Eq.( 5). In order to illustrate this point, a geometry with two limiters in the SOL has been considered. A secondary limiter with a slightly smaller parallel extent ∆ ′ s = 0.06 has been added at s = 0.75 in the parallel direction (s = 0.5 for the main limiter) and r = 0.75 in the radial direction (r = 0.5 for the main limiter). Otherwise its shape is similar to that of the main limiter Fig. 1b. The main and secondary limiters generate the density and parallel flow patterns shown on Fig. 13. As expected, the symmetry with respect to the main limiter is now broken. The secondary SOL, with short connection length between the limiters exhibits, a lower density since the parallel life time in this region is typically a factor 4 shorter than for the single limiter SOL addressed in the previous Section. This directly translates into a factor 2 decrease of the e-folding length with respect to the single limiter case, λ n /∆ r = 0.058, hence governing a low secondary SOL density. The e-folding length of the primary SOL is also slightly modified, λ n /∆ r = 0.116. The profile of the parallel averaged density exhibits these two radial regions with the two different e-folding rates one from r ≈ 0.5 to r ≈ 0.75, the other from r ≈ 0.75 to r ≈ 1.0. Interestingly enough the average density still appears to follow the behavior that one expects from the simple theoretical analysis, Fig. 14. Distortions with respect to such results are likely to mainly occur in the vicinity of the limiters as reported in [START_REF] Isoardi | Modelling SOL flow pattern spreading in the edge plasma[END_REF], this region being of immediate interest for the particle and energy handling issues. The flux pattern is more strongly modified with a straightforward modification of the SOL flow pattern. This effect also has an imprint into the edge region although the secondary limiter is localized in the far SOL at typically 2 e-folding lengths from the separatrix determined by the main limiter head. The secondary limiter, thus introduces further dimensionless parameters (3 with this simplified geometry) that modify the overall solution, including the edge plasma, and could therefore play a role in the scaling laws. This would enforce the effect of numerous parameters related to the detailed geometry of the far SOL limiters and offer a possible explanation for the observed dispersion in the SOL data base.

Discussion and concluding remarks

The present work shows that a volume penalization technique can be used efficiently to model the plasma/limiter interaction in the frame of a two-dimensional minimal transport model based on plasma density and parallel momentum. The penalization issues have been addressed by analyzing a one-dimensional nonlinear hyperbolic problem corresponding to the rapid transport in the parallel direction. In that way, this numerical work extends former mathematical studies restricted to linear equations. The penalization error corresponding to the difference between the solution of the penalized equations and the exact solution of the equations with the surface Bohm boundary condition can be controlled by choosing appropriate values of the penalization parameter η. As reported in the literature, the local error defined inside the limiter is shown to decay linearly with η, that gives a real control on the way the boundary conditions are imposed numerically. Since the periodicity of the domain is recovered, present results have also shown the suitability of Fourier spectral methods when the limiter is modeled by a smoothed mask function. As expected the convergence of the spectral method is severely reduced by the strong gradients but the solution can be freed of Gibbs oscillations. The major interest of spectral method will appear further for three-dimensional computations. Indeed, elliptic problems spectrally discretized in toroidal (naturally periodic) and poloidal directions will be reduce to one dimensional systems more quickly resolved than the 2D problems provided by the finite difference approximation. The present method appears to provide a plasma velocity which is almost sonic at the interface between the plasma and the obstacles as expected from the sheath conditions through the Bohm criterion. This indicates that for reasonable choices of the constrained Mach number profile within the obstacles, the numerical solutions will exhibit close to sonic boundary conditions, and, most important, the proper value of the particle outflux from the plasma. Moreover, two-dimensional simulations have shown that the physics of the SOL is recovered, in particular regarding the radial profile of the density averaged in the parallel direction. This approach also allows one to readily take into account the radial diffusive transport of particles to the limiter head. This is particularly important when the limiter has a finite parallel extension. A consequence of this new physics is to slightly modify the flow pattern, in particular by enforcing a flow in the edge plasma to compensate for this specific sink. As a consequence, the total plasma pressure variation, which characterizes the diffusive momentum transfer, is slightly modified although the overall shape is unaffected. It exhibits a momentum source region localised in the edge and SOL plasma near the separatrix coupled to momentum sinks both in the plasma core and in the outer SOL. Since the total pressure is larger in the edge plasma than in the outer SOL, the momentum spreading is more effective toward the outer SOL.

Having validated the penalization method for the plasma/limiter interaction, it is straightforward to address the multi-limiter case. We have considered a geometry with a primary limiter together with a secondary limiter typically at 2 e-folding lengths in the far SOL. This simulation has allowed us to analyse the properties of the secondary SOL and in particular its small e-folding length, λ n . One has also found that the plasma flow pattern was modified, including the flow pattern in the edge plasma. One therefore finds that the transport properties are changed by the secondary distant limiter. Consequently, one must introduce 3 further dimensionless parameters in order to account for the simplified geometry of the secondary limiter. Potentially, this can give a first principle understanding of the scatter reported in the empirical data. This would also hinder the capability of extracting edge plasma scaling laws without taking into account the features of the components that modify the far SOL transport.

The complexity introduced by the Bohm boundary condition for the plasma-wall interaction has strongly restrained the simulations of the edge and SOL plasmas. The penalization technique presented in this paper alleviates these difficulties and allows one to implement in a rather straightforward fashion realistic plasma facing components. Of course, introducing small objects will have implication in the required code resolution. However, for a code based on a proper description of the turbulent transport, hence with a typical mesh size given by the ion Larmor radius, this issue will be addressed anyhow. As shown by the few results proposed in this paper, the capability to properly describe these components can affect the overall transport pattern up to the edge plasma. The penalization technique thus appears as a powerful tool to investigate the SOL and edge plasma transport properties, and in particular assess the geometrical effects first addressed in the examples chosen to illustrate the present work.

Acknowledgement

This work was carried out within the framework of the 'Fédération de Recherche Fusion par Confinement Magnétique -ITER'. This work has been supported by the ANR ESPOIR (Edge Simulation of the Physics Of Iter Relevant turbulent transport). This work supported by the European Communities under the contract of Association between EURATOM and CEA was 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.

A. Exact solution for the 1D case

We consider the following system of stationary equations

By symmetry, one can assume N(1) = N(0) = S 2 , which yields: Γ(s) = S(s -0.5).

(

) 22
The second equation leads one to:

Γ 2 + N 2 = SN. ( 23
)
Using the previous expression of Γ one then obtains the density on for the subsonic branch (the other solution corresponds to the supersonic branch):

N(s) = S 0.5 + s(1 -s)

0.5 . (24)
By a translation of s, one can derive the solution for the 1D configuration described in section 3.2, where the limiter is located between s L = 0.5-d and s R = 0.5 + d. Introducing either the variable z L = s/s L or z R = -(1 -s)/s L , the solution of interest then reads:

• if s ∈ [0, s L ], z L ∈ [0, 1]:    N exact (s) = S s L 1 + 1 -z 2 L 0.5 Γ exact (s) = S s L z L . (25)
• if s ∈ [s R , 1], z R ∈ [-1, 0]    N exact (s) = S s L 1 + 1 -z 2 R ) 0.5 Γ exact (s) = S s L z R . (26)

Figure 1 :

 1 Figure 1: Geometry of the model. (a) Schematic representation of a toroidal tokamak with a circular poloidal section with a toroidal limiter. The white line shows an open magnetic field line of length L // . (b) Poloidal cut of the vacuum chamber. The inner white annulus shows the edge region not taken into account. (c) Two-dimensional computational domain including the transition between the edge and the SOL region up to the wall.

Figure 2 :Figure 3 :

 23 Figure 2: One-dimensional profiles of the numerical (circles) and analytical (solid line) solutions for density N (a) and parallel momentum Γ (b). Vertical dotted lines indicate the limiter location. Resolution n = 100, penalization parameter η = 10 -10 and function M Ω (s) = 2χ(s)(0.5 -s)/∆ s are used.

Figure 4 :

 4 Figure 4: Local ( * ) and global (o) L 1 errors computed from Eqs. (15-16) versus the penalization parameter η. Dashed line represents the function η.

Figure 5 :

 5 Figure 5: Mach number profiles around the limiter (o) for two functions M Ω . (a) M Ω = 4χ(s)(0.5 -s)/∆ s (line). (b) M Ω = χ(s)(0.5 -s)/∆ s . Vertical dotted indicate the limiter location. Resolution n = 100 and penalization parameter η = 10 -10 are used.

Figure 6 :

 6 Figure 6: One dimensional spectral solution profiles (circles) and analytical solution (solid line) for density N (a) and parallel momentum Γ (b), with the magnification of the momentum profile closed to the limiter. Vertical dotted lines indicate the limiter location. Fourier mode n = 128, penalization parameter η = 10 -7 and filtering coefficient α = 10 are used.

Figure 7 :

 7 Figure 7: Mask functions defining the limiter. Smoothed (o) and step (+) χ-functions. Fourier mode n = 256 is used.

Figure 8 :

 8 Figure 8: One dimensional spectral solution profiles using the smoothed characteristic χ-function . Density N (a) and parallel momentum Γ (b). Fourier mode n = 128, penalization parameter η = 10 -7 and filtering coefficient α = 10 are used.

Figure 9 :

 9 Figure 9: (a) Local L 1 error ( * ) from Eqs.(16) versus the penalization parameter η. The dashed line represents the function η. (b) Global L 2 error from Eq.(19) for the density (o) and the momentum (*). Penalization parameter η = 10 -7 and filtering coefficient α = 10 are used.

Figure 10 :

 10 Figure 10: Absolute value of the differences between benchmark solution computed with n = 4096 Fourier modes and solutions at lower resolutions. (a) for the density and (b) for the particle flux. n = 128 (dashed line), n = 1024 (dotted line) and n = 2048 (solid line).In all computations, the penalization parameter η = 10 -7 and the filtering coefficient α = 10 are used.

Figure 11 :

 11 Figure 11: Contours plot of density (a,c) and particle flux (b,d) in steady state conditions given by Sc = 1 and A = 50. The radial extension of the edge region is between r = 0 and r = 0.5. The SOL is defined between r = 0.5 and r = 1. Strong diffusive transport at P e = 100 (a,b) and weak diffusive transport at P e = 5000 (c,d). Fourier modes n = 128 and n r = 50, penalization parameter η = 10 -7 and filtering coefficient α = 10 are used.

Figure 12 :

 12 Figure12: Mean density (a) and total pressure variation from the stagnation point close to the limiter (b) profiles in the radial direction for Sc = 1, P e = 5000, A = 50. Solutions from reference[START_REF] Isoardi | Modelling SOL flow pattern spreading in the edge plasma[END_REF] with a limiter with zero thickness imposed by surface Bohm boundary conditions (solid line), the present penalization (dashed line) and the analytical slope[START_REF] Isoardi | Modelling SOL flow pattern spreading in the edge plasma[END_REF] (dotted line). The interaction between the edge and the SOL is for r = 0.5. Fourier modes n = 128 and n r = 50, penalization parameter η = 10 -7 and filtering coefficient α = 10 are used.

Figure 13 :Figure 14 :

 1314 Figure 13: Contours plot of density (a) and particle flux (b) for two limiters in steady state conditions given by Sc = 1, P e = 1000 and A = 50. The radial extension of the edge region is between r = 0 and r = 0.5. The SOL is defined between r = 0.5 and r = 1. Fourier modes n = 128 and n r = 50, penalization parameter η = 10 -7 and filtering coefficient α = 10 are used.

2 NS = 1 0∂

 21 + N = 0 s ∈ ]0, 1[, M(0) = Γ(0) N (0) = -1, M(1) = Γ(1) N (1) = 1. (20)Assuming S is constant, by integration of the first equation, one readily finds that: s Γ ds = Γ(1) -Γ(0) = N(1) + N(0).


Guo Xian
email: tanguoxian@pmail.ntu.edu.sg

Tan Christian

Viard-Gaudin Alex
email: christian.viard-gaudin@univ-nantes.fr

C Kot

Online Writer Identification Using Fuzzy C-means Clustering of Character Prototypes

Keywords: Writer identification, information retrieval, online handwriting, fuzzy c-means

New kinds of documents such as handwritten online documents are emerging, which are produced by digital devices such as Tablet PC, personal handheld devices or digital paper coupled with digital pens. The rapid increase in the number of such handwritten online documents leads to mounting pressure on finding innovative solutions towards faster processing, indexing and retrieval of the documents from databases. One such method is to extract writer information derived from the raw ink signal for indexing and retrieval of the documents. This paper proposes a text independent method that does not place any constraints on the content being written or writing styles of the writers. We subsequently extract writer information at the character level from online handwritten documents and present a fuzzy c-means approach to cluster and classify the character prototypes for writer identification. The proposed system attained an accuracy of 97.6% on 82 writers and an accuracy of 98.3% when retrieved from a scaled up larger database of 120 writers.

Introduction

Technology has become an integral part of modern lifestyles that our lives are becoming intertwined with technology itself. Numerous initiatives have been funded to research and develop more efficient algorithms, software and computing platforms to handle the surge in demand for seamless interaction and to deliver interactive environments with a new level of intelligence [START_REF]Moblile PC User Experience Guidelines for Developers[END_REF][START_REF] Oviatt | Designing the User Interface for Multimodal Speech and Pen-Based Gesture Applications: State-of-the-Art Systems and Future Research Directions[END_REF]. All these led to the emergence and proliferation of a kind of document: handwritten on-line documents. They are produced by state-of-the-art devices such as Tablet PC, personal handheld devices or digital paper coupled with digital stylus pens [START_REF]Mobility and Tablet PC Road Map[END_REF]. The rapid increase in the number of such documents requires efficient management tools for proper indexing and retrieving from databases.

Online handwritten digital documents are defined as those digital documents that not only provide information obtainable from offline digital documents, but also contain temporal information of the handwriting process [START_REF] Jain | Indexing and Retrieval of On-line Handwritten Documents[END_REF]. Such additional information provides vital clues as to the identities of the writer. Writer identification systems must be clearly distinguished from writer verification systems.

Writer verification performs a one-to-one matching between a test writer and a database of writers and attempts to ascertain the authenticity of the test writer. On the other hand, writer identification involves executing a one-to-many match and returns a ranked list of results for the search. The difference, though subtle, lies in the applications in which they can be utilized in.

Online document indexing using writer information provides two-fold distinct advantages. Firstly, from information security's point of view, writer identification has ubiquitous applications in digital rights management and forensic analysis in the prevention of fraud and identity theft cases. Secondly, in environments where large amounts of documents, forms, notes and meeting minutes are constantly being processed and managed, knowing the identity of the writer would provide an additional value. One such application is to process and retrieve the identities of students for subsequent verification purposes.

Online writer recognition systems can make use of global features such as texture, curvature and slant features [START_REF] Busch | Texture for Script Identification[END_REF][START_REF] Hochberg | Automatic script identification from document images using clusterbased templates[END_REF] as well as a combination of local features such as graphemes, allographs 1 and connected components [START_REF] Schomaker | Automatic Writer Identification Using Connected-Component Contours and Edge-Based Features of Uppercase Western Script[END_REF][START_REF] Srinivasan | On Computing Strength of Evidence for Writer Verification[END_REF]. They can be generally classified into text-dependent or text-independent techniques. Previous works for online writer recognition such as the method proposed by Pitak et al. [START_REF] Pitak | On-line Writer Recognition for Thai Based on Velocity of Bary center of Pen-point Movement[END_REF] adopted a Fourier transformation approach. The extracted features are the velocities of the barycenter of the pen movements and they are transformed into the frequency domain using Fourier transform. The advantage in adopting such a model is that it is text-independent, but at the expense of a lower noise tolerance. The noise must be filtered out as much as possible in the pre-processing stage, otherwise the noise might be mistaken for high velocity components once the features are transformed into the frequency domain. Bensefia et al. [START_REF] Bensefia | A Writer Identification and Verification System[END_REF][START_REF] Bensefia | Information Retrieval Based Writer Identification[END_REF] proposed using a sequential clustering approach at the grapheme level to categorize different writers for their writer identification system. This approach attained an identification rate of 86% on a test set of 150 writers. The advantage of this method is that it does not depend on any lexicon and is therefore language independent. Works using allograph based methods are also popular in writer identification. Niels et al. used a dynamic time warping approach [START_REF] Jain | Indexing and Retrieval of On-line Handwritten Documents[END_REF][START_REF] Niels | Generating Copybooks From Consistent Handwriting Styles[END_REF][START_REF] Niels | Automatic Allograph Matching In Forensic Writer Identification[END_REF] to hierarchically cluster allographs and build a set of membership vectors, which contains the frequency of occurrence of each allograph for each character. This prototypic template of membership vectors then represents the handwriting styles of the different writers. Niels et al. reported a top1 accuracy of 89% based on this method [START_REF] Niels | Automatic Allograph Matching In Forensic Writer Identification[END_REF]. Chan et al. [START_REF] Chan | Online Text Independent Writer Identification Using Character Prototypes Distribution[END_REF] also made use of a character prototype distribution to model the specific allographs used by a given writer. They managed to achieve a top1 accuracy of 95% based on this textindependent approach. Even though working at the character level as opposed to using the grapheme or word level appears to be quite challenging, character prototype approaches are able to produce a more consistent set of templates for writer identification. Our proposed work further improves upon the work done by Chan et al. by adopting a fuzzy c-means algorithm, resulting in significant improvements in the performance.

The remainder of this paper is organized as follows: Section 2 describes the proposed methodology and experimental setup. Preprocessing such as normalization and resampling and our proposed fuzzy c-means algorithm is discussed here. Section 3 then presents the experimental results. Finally, discussions and future areas to explore are given in section 4.

Proposed Writer Identification Method

The writer identification can be divided into three stages, namely the prototype training stage, the reference and test document labeling stage and finally, the classification stage. During the prototype training stage, prototypes built at the character level are trained using the IRONOFF database [START_REF] Viard-Gaudin | The IRESTE On/Off (IRONOFF) Dual Handwriting Database[END_REF] of 16585 isolated French words written by 373 subjects. The purpose of this stage is to build a set of character prototypes using the 16585 isolated words to model the different allographs of the 26 Latin alphabets ('a' to 'z'). Following this, in the document labeling stage, automatically segmented characters from a set of 82 French reference and 82 test documents (extended later to 120 reference and 120 test documents) are then mapped to the set of prototypes built previously and transformed into a distribution of frequency vectors. These reference and test documents that we collected belong to a separate dataset from the IRONOFF dataset. A separate dataset needed to be collected because the IRONOFF dataset contains only isolated words and hence are not representative of actual online documents. Finally, the frequency vectors are used for classification in order to identify the writer corresponding to the test document. A detailed account of each of the three stages is given in [START_REF] Chan | Online Text Independent Writer Identification Using Character Prototypes Distribution[END_REF]. Figure 1 shows the block diagram of the entire writer identification process.

FIG. 1 -Block diagram for proposed methodology

One key point of the method relies on the automatic segmentation and labeling of the text, which is performed at the character level. This has to be done to extract and define the allographic prototypes from an independent word database (IRONOFF) and subsequently on all the reference documents and test documents to assess the writer retrieval performance of the proposed method. An industrial character segmentation and recognition engine, "MyScript SDK" [START_REF]Vision Objects Industrial Text Recogniser SDK[END_REF], with the French linguistic resource attached for increased accuracy, has been used for this purpose. After the characters are segmented, the segmented characters then underwent further preprocessing where the size of each segmented character is normalized and resampled to 30 points. A process of feature extraction on each of the resampled points is then carried out. The features being used are the x and y co- Segmentation ordinates, the directions of x and y co-ordinates, the curvatures of x and y co-ordinates and the Pen-up or Pendown information [START_REF] Chan | Online Text Independent Writer Identification Using Character Prototypes Distribution[END_REF]. Accordingly, the matching process between a character extracted from a document and a prototype is carried in a feature space of dimension: 30 points × 7 features = 210.

In this paper, we propose a fuzzy c-means algorithm in the document labeling stage. The prototype training stage serves to identify common individual handwriting styles into individual prototypes at the character level. Following this, the document labeling stage then utilizes these prototypes to create individual distributions of handwriting styles for each of the test and reference documents in the database. Based on the results of the distributions, they provide statistical information about the handwriting styles of each writer. Our proposed method adopts a fuzzy cmeans algorithm [START_REF] Han | Data Mining: Concepts and Techniques[END_REF] which uses an exponential kernel function as described in Eq. 1 to create these individual distributions of handwriting styles.

( ) (

)

1 1 exp ( , ) exp ( , ) p k k i p k M N p k dist x x C dist x x α α α α α β β = = -× = -× ∑ ∑ (1)
where

k i C α is the total amount of characters from alphabet α, α { } ' ', ' ',..., ' ' a b z ∈
, that is assigned to prototype k from writer i. In Eq. 1, p represents a given segmented character that has been recognized as of the alphabet α,

p { } M ,..., 3 , 2 , 1 ∈
, with M being the number of characters corresponding to alphabet α.

( , )

p k dist x x α α
is the Euclidean distance between the feature vector x for point p of alphabet α and the feature vector x for prototype k of alphabet α. In Eq. 1, β is a tuning parameter which is set to be 0.01 in our experiments.

Characters from the reference and test documents are then assigned a partial membership to the prototypes based on their distance metric to the prototypes. Therefore, characters which lie further away from certain prototypes are assigned a lesser degree to that particular prototype. k i C α is then used to calculate the distribution of frequency vectors [START_REF] Chan | Online Text Independent Writer Identification Using Character Prototypes Distribution[END_REF] to be used during classification.

It can be seen that a character does not need to be labeled to one particular prototype. There should be a fuzzy logic behind labeling the character to all the prototypes that are present and then building a distribution based on this fuzzy logic. The experimental results obtained serves as strong evidence to attest that our proposed fuzzy c-means approach for clustering into prototypes during the document labeling stage yielded a higher level of accuracy.

Experimental Results

The first set of experiments was initially conducted on a smaller set of 82 test and 82 reference documents from 82 different writers. As presented in table 1, the proposed methodology using fuzzy c-means to label the test and reference documents to the prototypes resulted in a high accuracy of 97.6% for writers that are ranked correctly in the top 1 position. This translates into a misclassification error of 2 misclassified writers out of 82, with both of them being misclassified in the top 2 position. This indicates that the writer identification system has confused the misclassified documents with only 1 other document from a different writer. Comparisons with previous results obtained by Siew et al. [START_REF] Chan | Online Text Independent Writer Identification Using Character Prototypes Distribution[END_REF], who also performed writer identification based on character prototyping, show a significant improvement over their proposed methodology. An accuracy of 95.1% (4 misclassified writers out of 82) was reported in their method where the four writers were wrongly identified ranked at top 2, 4, 9, and 12 positions. This means that their writer identification system has confused the misclassified documents with up to 11 other documents. Therefore, our proposed methodology is able to perform with significantly higher accuracies. This improvement over Siew et al.'s results can be explained as follows. Their methodology hinges on the concept that each character can only be assigned to one particular prototype for which a distribution of handwriting styles is built. This is flawed in reality because there often exist overlapping handwriting styles for different writers. Our observations reveal that there are numerous instances when the characters are close to more than one prototype in the vector space. This can be explained by the fact that a writer can have strong, dominant handwriting style and weak handwriting styles. Weak handwriting styles change according to various circumstantial and temporal states [START_REF] Huber | Handwriting Identification: Facts and Fundamentals[END_REF], which can affect the strong dominant handwriting style and lead to reminiscence of multiple overlapping handwriting styles. Therefore, in our proposed methodology, each character is not just assigned to one particular prototype, but rather, each character is assigned a certain degree of all the prototypes depending on how close they are to that prototype. The more similar the character is to a certain prototype, the greater the degree that prototype has on the character. For example, one character might be attributed among three different prototypes if the character shares common traits of all three handwriting styles.

Experiments were also performed using the Kullback-Leibler (KL) divergence [START_REF] Cover | Elements of Information Theory[END_REF] as a different metric for the fuzzy c-means algorithm to determine the best performing metric for our writer identification system. We can observe from Table 1 that using KL divergence resulted in a low 87.8% identification rate. This can be attributed to the asymmetric nature of the KL divergence. Therefore, in our system, the better performing metric to use for our fuzzy cmeans algorithm is the Euclidean distance metric.

The database was subsequently enlarged to 120 test and 120 reference documents from 120 different writers. The purpose of this experiment was to investigate the scalability of the system as to whether it can handle larger loads on an enlarged database. Table 1 show that a top 1 accuracy of 98.3% (2 misclassified writers out of 120) was obtained with both the remaining two wrongly classified writers at the top 2 positions. This indicates that the proposed methodology is indeed scalable and can handle applications where scalability is critical in designing the system.

Fuzzy C-means Kernel Design

Experiments were also conducted using different kernel functions [START_REF] Hoppner | Fuzzy Cluster Analysis: Methods for Classification, Data Analysis and Image Recognition[END_REF] for the fuzzy c-means algorithm to determine the kernel function that can perform best in our writer identification system.

Gaussian kernel function:

The distribution of the feature vectors was assumed to be Gaussian with zero mean and unit variance.

2. Inverse kernel function: Eq.2 describes the formulation for the inverse kernel function where the notations have the same meaning as that shown in Eq. 1.

1 1 1 ( , ) 1 ( , ) p k k i p k M N p k dist x x C dist x x α α α α α = = = ∑ ∑ (2)
Table 2 shows the comparison in performance of the writer identification among using the three different kernel functions. It can be seen that the inverse kernel function performs poorly, which can be explained by the poor behavior of such kernel functions as it approaches the centroids of the prototypes. The results obtained using the Gaussian kernel function is similar to that for the exponential kernel function. We believe that adjustments in the mean and variance parameters will allow the Gaussian kernel function to achieve the same results as the exponential kernel function since Gaussian functions are essentially exponential in nature. Identification rate using exponential kernel function (Eq. [START_REF]Moblile PC User Experience Guidelines for Developers[END_REF] 98.3%

Identification rate using Gaussian kernel function 97.5%

Identification rate using inverse kernel function (Eq. 2) 96.7%

Effect of number of Character Prototypes on Accuracy

It is hypothesized that different alphabets require different number of prototypes to effectively model all the possible writing styles of that character. For example, there are more ways and styles to write the alphabet 'f' than to write the alphabet 'c'. A preliminary level of analysis has been performed to find a global optimal value for the number of prototypes needed. Figure 2 shows the achievable by varying the number of global prototypes used for every letter of the alphabet. As seen from figure 2, the identification rate is highest when the number of prototypes varies from 10 to 30.

Additional number of prototypes beyond the 30 prototypes will result in a drop in the performance of the identification system. This can be explained by the principle of Occam's razor. A large number of prototypes create sparse dimensionality which deteriorates the performance of the classification. Likewise, insufficient number of prototypes will be unable to effectively separate between intra-class variations. Based on the above analysis, the optimum number of global prototypes is taken to be 10 in the experiment.

Discussions

From the experimental results, the proposed methodology is able to generate high accuracies of 97.6% (two misclassified writers out of 82) for the identification of 82 different writers, with both the misclassified writers being identified correctly in the rank 2 nd position. This is a remarkable improvement over Siew et al.'s [START_REF] Chan | Online Text Independent Writer Identification Using Character Prototypes Distribution[END_REF] methodology which only attained a top 1 accuracy of 95.12% (four misclassified writers out of 82), with all the misclassified writers only being correctly identified in the rank 11 th position. This can be explained by the fact that our proposed methodology is based on fuzzy logic that mimics the fact that a writer's handwriting style has reminiscence of multiple overlapping handwriting styles. Furthermore, when we further increase the size of the database to 120 different writers, a similar result of having only two misclassified writers out of 120 was attained, leading to a top 1 accuracy of 98.3%. This concludes that the system can provide high accuracies in identifying writers and is highly scalable as well. Our results also indicated that KL divergence performed poorly in our writer identification system. This might be due to the fact that KL divergence itself is asymmetric. Furthermore, we have also showed that the optimum number of prototypes to use for our writer identification system is 10.

One main area of improvement to be explored is to investigate how different alphabets affect the character prototypes. Certain alphabets have more writing styles than others and are more discriminatory in writer identification. For instance, character shapes like 'f', 'h' or 't' have more variability than others like 'e' or 'c'. In this paper, we have presented the global optimum number of prototypes to use. However, we can go one step further in future to investigate a variable number of prototypes required to effectively model each alphabet, which may yield even better performance. Furthermore, more work needs to be done to investigate the extent of impact that different alphabets have on the accuracy of writer identification. In this way, a dynamic algorithm can then be proposed to handle and create an alphabet matrix for each writer. The alphabet matrix will be constructed based on the discriminative power of different alphabets for each writer and will then determine if certain alphabets can be ignored for writer identification. This will result in a smaller feature vector set being used, which will thus reduce the computational complexity of the system.

Figure 2 .

 2 Figure 2. Graph of Identification Rate against Number of Clusters

Table 1 .

 1 Performance of writer identification using different distance metrics

	Size of		Fuzzy C-means
	reference document	1NN 1	Euclidean as	KL divergence 2 as
	database		distance metric	distance metric
	82	95.1%	97.6%	87.8%
	120	96.7%	98.3%	91.7%

Table 2 .

 2 Performance of the Fuzzy c-means algorithm using different kernel functions.

Allographs are different shapes and forms of the same alphabet.

1-Nearest Neighbor algorithm adopted by Siew et al.[START_REF] Chan | Online Text Independent Writer Identification Using Character Prototypes Distribution[END_REF]

Kullback-Leibler divergence

Acknowledgements

This research is jointly supported by Nanyang Technological University of Singapore, the French Merlion Scholarship and the ANR grant CIEL 06-TLOG-009.