E Prulière
email: pruliere@ujf-grenoble.fr

A Ammar

F Chinesta

Empirical Natural Closure Relation for Short Fiber Suspension Models

Keywords: Short fiber suspensions, Closure relations, Numerical modeling, Model reduction, Karhunen-Loève decomposition

This work focuses on the resolution of the Fokker-Planck equation that governs the evolution of the fibers orientation distribution. To reduce the computing time, that equation is solved along some flow trajectories in order to extract the significant information of the solution from the application of the Karhunen-Loève decomposition. Now, from this information one could solve the Fokker-Planck equation everywhere in the flow domain or simply adjust a closure relation that becomes optimal for such flow, solving the evolution of some orientation moments which require a less amount of computation. This paper focuses on this last strategy. For this purpose we start introducing the Karhunen-Loève decomposition that is applied later to automatically extract the main solution characteristics for adjusting empirically a natural closure relation.

Introduction

Short Fiber Suspensions flow models

Numerical modeling of non-Newtonian flows usually involves the coupling between equations of motion, which define an elliptic problem, and the fluid constitutive equation, which introduces an advection problem related to the fluid history. In short fiber suspensions (SFS) models, the extra-stress tensor depends on the fiber orientation whose evolution can be modeled from a transport problem. In all cases the flow kinematics and the fiber orientation are coupled: the kinematics of the flow governs the fiber orientation, and the presence and orientation of the fibers modify the flow kinematics. Thus, for example, in a contraction flow of a dilute suspension, large recirculating areas appear [START_REF] Lipscomb | The flow of fiber suspensions in complex geometries[END_REF]).

If one uses SFS flows in material forming processes, the final fiber orientation state depends on the process and exhibits flow-induced anisotropy. Thus, we need to compute the fiber orientation in order to predict the final mechanical properties of the composite parts, which depend strongly on the fiber orientation.

Moreover, the numerical simulation of such flows becomes interesting if one wants to identify their rheological parameters using some rheometric devices and an appropriate inverse technique.

The mechanical model governing the SFS flow is given by the following equations: [START_REF] Batchelor | Slender-body theory for particles of arbitrary cross-section in Stokes flow[END_REF], [START_REF] Hand | A theory of anisotropic fluids[END_REF], Hinch andLeal (1975, 1976))

• The momentum balance equation, when the inertia and mass terms are neglected, results

0 Divσ = (1)
where σ is the stress tensor.

• The mass balance equation for incompressible fluids 0 Divv =

(2) where v represents the velocity field.

• The constitutive equation for a dilute suspension of high aspect-ratio particles is given, with other simplifying assumptions [START_REF] Ch | Flow regimes for fiber suspensions in narrow gap[END_REF]), by where p denotes the pressure, I the unit tensor, η the viscosity which depends on the chosen model, D the str ai n ra te te ns or , p N a scalar parameter depending on both the fiber concentration and the fiber aspect ratio, " : " the tensorial product twice contracted (i.e.

kl ijkl ij D a D a =       :
) and a the fourth order orientation tensor defined by: ()

ad ρρρρ ψ ρρ =⊗ ⊗ ⊗ ∫ ( 4
)
where ρ is the unit vector aligned in the fiber axis direction, " ⊗ " denotes the tensorial product (i.e.

( ) We can also define the second order orientation tensor as:

j i ij ρ ρ = ρ ⊗ ρ ),
()

ad ρρ ψ ρρ =⊗ ∫ ( 6 )
It is easy to verify that if

) ( ) ( ρ - ρ δ = ρ ψ
, the fourth order orientation tensor can be written as aaa =⊗ ( 7 ) whose components are defined by kl ij ijkl a a a = .

For general expressions of ) (ρ ψ the previous relation is not exact, and equation ( 7) becomes a closure approximation known as the quadratic closure relation:

The 3D isotropic orientation state is defined by the uniform distribution function

() 1 4 ρ π Ψ= ( 1 2 )
and then, the second order orientation tensor related to that isotropic orientation state is 3

I a = ( 1 3 )
It is easy to verify that for isotropic orientation distributions (2D or 3D) the linear closure becomes exact.

• If we consider spheroidal fibers immersed in a dilute suspension, we can describe the orientation evolution by means of the Jeffery equation (Jeffery (1922))

( ) ( ) ( ) : d kD D dt ρ ρρρ ρ ρ =Ω + - ⊗ ( 14
)
where Ω is the vorticity tensor, and k is a constant that depends on the fiber aspect ratio r (fiber length to fiber diameter ratio):

( ) ( )

22 11 kr r =- + (15)
On the other hand the evolution of the fiber orientation distribution ψ is governed by the Fokker-Planck equation, (

() 0 dd dt dt ψρ ρ ψρ ρ  ∂ +=  ∂  (16) )
where the material derivative is given by:

d vG r a d dt t ψψ ψ ∂ =+ ∂ ( 1 7 )
Now, taking into account equations ( 6), ( 14) and ( 16), the equation that governs the evolution of the second order orientation tensor can be deduced

( ) ( ) 2 : da aa k D aa D aD dt =Ω -Ω+ + - (18)
A similar equation can be derived for the evolution of the fourth order orientation tensor, which in this case involves the sixth-order orientation tensor.

To take account of fiber interaction effects in semi-concentrated suspensions [START_REF] Folgar | Orientation behavior of fibers in concentrated suspensions[END_REF] pr opo se d the i ntro duc tio n of a di ffus io n ter m i n the F okker -P la nc k eq ua ti on, i .e .

() () ()

r dd D dt dt ψρ ρ ψρ ψρ ρ ρ ρ  ∂  ∂∂  +=    ∂∂ ∂     (19)
Fiber interaction being taken into account, the equation governing the evolution of a then yields:

( ) ( ) 2 : 4 r d da I aa k D aa D aD Da dt N  =Ω -Ω+ + - - -   (20)
The Fokker-Planck formalism circumvents the necessity of using closure relations, but it induces some difficulties related to its multidimensional character (the distribution function is defined in the physical and the configuration spaces) and moreover advection terms are defined in both spaces. By these reasons the number of works devoted to the treatment of the FP equation is relatively reduced (Lozinski and Chauvière (2003); Chauvière and [START_REF] Chauviere | Simulation of dilute polymer solutions using a Fokker-Planck equation[END_REF]). In these techniques, usually, to account for the multidimensional character of the FP equation, a time-splitting is often considered to decouple the advection problem in physical space and the advection-diffusion problem in the conformation space. The first problem can be solved by a numerical method appropriate for hyperbolic partial differential equations (discontinuous Galerkin, SUPG, …). Then, the advection-diffusion problem defined in the conformation space can be treated using different implicit techniques (SUPG, wavelets-Galerkin, spectral techniques, …) preserving stability, accounting for distribution relatively localized as well as periodic boundary conditions in the conformation space.

In general we can solve the FP equation from its associated Ito stochastic differential equation for a large set of realizations. The CONNFFESSIT method (Ottinger and Laso (1992)) was the first implementation of the stochastic approach. The Brownian Configuration Fields [START_REF] Hulsen | Simulation of viscoelastic flows using Brownian configuration fields[END_REF]) can be considered as an improvement of the CONNFFESSIT method. However, the control of the statistical noise is a major issue in stochastic micro-macro simulations, problem that does not arise in the deterministic Fokker-Planck approach.

This work focuses on the resolution of the Fokker-Planck equation that governs the evolution of the fibers orientation distribution. To reduce the computing time, that equation is solved along some flow trajectories in order to extract the significant information of the solution from the application of the Karhunen-Loève d ec omp osi tio n. N ow , fr om this info rma ti on o ne c ould s ol ve the F okker -P la nc k eq uati on e ver yw her e i n the flow domain or simply adjust a closure relation that becomes optimal for such flow, solving the evolution of some orientation moments which require a less amount of computation. Some antecedents exist concerning closure relation fitting; see for example [START_REF] Parsheh | Investigation of closure approximations for fiber orientation distribution in contracting turbulent flow[END_REF] or [START_REF] Gillet-Chaulet | Flow-induced anisotropy in polar ice and related ice-sheet flow modelling[END_REF].

The present paper focuses on the definition of an empirical natural closure approximation for each particular flow. For this purpose we start introducing the Karhunen-Loève decomposition [START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF]) that is applied later to automatically extract the main solution characteristics for adjusting empirically a natural closure relation.

The Karhunen-Loève decomposition

We assume that the evolution of a certain field ( ) t x u , is known. In practical applications, this field is expressed in a discrete form, that is, it is known at the nodes of a spatial mesh and for some times

( ) p i p i u t x u ≡ ,
. We can also write introducing a spatial interpolation ( ) ( )

, p ii ux u x tp t ≡= ∆ , [ ] 1, , pP ∀∈ A [ ] 1, , iN ∀∈ A
. The main idea of the Karhunen-Loève (KL) decomposition is how to obtain the most typical or characteristic structure () x φ among these

() p x u p ∀ ,
. This is equivalent to

obtaining a function () x φ that maximi ze s α . () () () ( ) 2 11 2 1 pP iN p ii pi iN i i xux x φ α φ = = == = =    = ∑∑ ∑ (21)
T he ma xi mi za ti on l ea ds to:

() () ( )( ) ()() 11 1 1 = ; pP jN iN pp ii j j pi j iN ii i xux x ux xx φ φ α φφ φ == = == = = =           =∀ ∑∑ ∑ ∑ # ## (22)
w hi ch c an be re w ri tte n in the fo rm

() ( ) ( ) () ()() 11 1 1 ; jN pP iN pp ij j i ijp iN ii i uxux x x xx φφ α φφ φ == = == = = =   =      =∀ ∑∑ ∑ ∑ # ## (23) Defining the vectors a such that its i-component is () i x a
, Eq. ( 23) takes the following matrix form

;

TT kk φφ α φφ φ φαφ =∀ ⇒ = ## # (24)
where the two points correlation matrix is given by

() ( ) ( ) 11 pP pP T pp pp ij ij pp ku x u xk u u == == =⇔ = ∑∑ (25)
which is symmetric and positive definite. If we define the matrix Q containing the discrete field history: 12 11 1 12 22 2 12

P P P NN N uu u uu u Q uu u    =    A A BB D B A (26)
is easy to verify that the matrix k in Eq. ( 25) results:

T kQ Q = (27)

"A posteriori" reduced modeling

If some direct simulations are carried out, we can determine ( )

[ ] [ ] P p N i u t x u p i p i , , 1 , , , 1 , , A A ∈ ∀ ∈ ∀ ≡ ,
and from these the n eigenvectors related to the n-highest eigenvalues

() [ ] [ ] n k N i x i k k , , 1 , , , 1 , A A ∈ ∀ ∈ ∀ φ = φ
that are expected to contain the most information about the problem solution. For this purpose we solve the eigenvalue problem defined by Eq. ( 24) retaining all the eigenavules belonging to the interval defined by the highest eigenvalue and that value divided by a large enough value ( 8 10 in our simulations). In practice n is much lower than N. Thus, we can try to use these n eigenfunctions for approximating the solution of a problem slightly different to the one that has served to define ( )

p i p i u t x u ≡ ,
. For this purpose we need to define the matrix B ( )

( ) ( ) () () () ( ) ( ) ( ) 11 1 12 22 2 12 12 n n NN N n xx x xx x B xx x φφ φ φφ φ φφ φ    =     A A BB D B A (28)
Now, if we consider the linear system of equations resulting from the discretization of a partial differential equation (PDE) in the form

() ( 1 ) mm KU F - = ( 2 9 )
where the superscript refers to the time step, then, assuming that the unknown vector contains the nodal degrees of freedom, it can be expressed as:

() () () 1 in mm m i i i UB ζφ ζ = = == ∑ (30)
Eq. ( 29) results

() ( 1 ) () ( 1 ) mm mm K U FK BF ζ -- = ⇒ = (31)
and by multiplying both terms by T B it results

() ( 1 ) Tm T m BKB BF ζ - = (32)
which proves that the final system of equations is of low order, i.e. the dimension of Remark 1. Equation ( 32) can be also derived introducing the approximation (30) into the PDE Galerkin form.

T B KB is n n × ,

. Simula ting c omple x flows

We consider the use of a decoupled strategy which solves the flow kinematics (from the orientation state known at the previous iteration or time step), and then the orientation state can be updated using the velocity field just computed. The flow kinematics can be solved accurately by using standard discretization techniques as for example the finite element method. This paper focuses on the resolution of the orientation problem which involves some specific difficulties.

Extracting the solution structure

One possibility for reducing the computing time associated with the resolution of the Fokker-Planck equation ( 18) consists of decoupling the resolution of the purely advection problem in the physical space and the advection-diffusion defined in the conformation space. Now, the first problem can be accurately integrated using the method of characteristics along some nodal trajectories, being the one defined in the conformation space discretized using some stabilized finite element strategy for accounting its advection character. The model reduction described later will operate only on the conformation space discretization.

We briefly summarize this strategy in the 2D case, the 3D case being a direct extension. In that case the Fokker-Planck equation can be written using polar coordinates as:

d dt DD k DD ϕ ϕ ϕϕ ϕ ϕ ϕϕ ϕ Ω   =- +   -Ω     -   -   (35)
where the components of both the vorticity and the strain rate tensor could depend on the physical coordinates x and time t.

Applying a stabilized finite element discretization of Eq. ( 33) in the conformation space ( ϕ ) it results (taking into account the solution periodicity in the angular coordinate):

22 2 * ** 00 0 2 * 0 r dd d dd d dt dt dt Dd π π π π ψϕ ψ ϕ ψ ϕ ψ ϕ ψψ ϕ ϕ ϕ ψψ ϕ ϕϕ ∂∂  ++ =  ∂∂  ∂∂ =- ∂∂ ∫∫ ∫ ∫ (36)
The weighting function * ψ is approximated as ψ whereas * ψ is approximated in a different way in order to account for the advection character of the differential operation. Different alternatives exist, out of the scope of the present work. In our simulations a SUPG (streamline upwind Petrov Galerkin) stabilization has been used. Finally, after introducing the different finite element interpolations and perform the numerical integration the following linear system is obtained:

11 12 1 21 KKF K KK F ψψ ψ ψ •• -- =+ ⇒ =+ (37)
In this equation vector F accounts for the normality condition of the orientation distribution that is e nfor ce d us ing a La gr ange multi pl ie r tec hni que.

Now, for integrating in the physical domain we use the method of characteristics by its simplicity and high accuracy. For this purpose, we assume that we are looking for the solution along the streamline related to point X where at the initial time t=0 the orientation distribution is known. Now, the orientation distribution updating is given by: Using this simple uncoupled strategy one can compute the evolution of the fiber orientation distribution along some nodal trajectories identified by the parameter l. Those distributions could be stored at some time steps (snapshots) leading to the following eigenvalue problems: Now, from the solution of Eq. ( 39) we can compute the most significant eignfunctions allowing to write:

ψ ψ ψ + ++ - == =+ = ∆   == = +    += + ∆   (38)
[ ] () ( l ) () ( l ) ; 1, , l l QQ l L φα φ =∈ A (39) with 12 () 1 () 1 () 1 12 () 2 () 2 () 2 () 12
() () () () () 1 () () () l in ll l l i li i tt B t ψζ φ ζ = = ≈= ∑ (41)
Vector ψ ha s c N components that are in fact the values of the distribution function at the nodes used to discretize the conformation domain (the unit circle or the unit surface in the 2D or 3D case respectively).

Now w e w ill define a global reduced approximation basis B from the different reduced approximation bases () l B . For this purpose, we proceed as follows:

• (1) BB = , 1 nn = and l=1 • While lL ≤ o 1 ll ←+ o 0 i = o While l in ≤ ̇ 1 ii ←+ ̇ If ( ) () () 1 n T ll ii j j j TOL φφ φ φ =  ->   ∑ Ü () l i BB φ  =  Ü 1 nn ←+
that is, w e update the reduced approximation basis introducing only vectors that are not too close to the e xi sti ng o ne s (for a gi ve n to le ra nc e TOL).

The final reduced approximation basis consists of the n vectors contained in matrix B , and it could be expected that that basis represents accurately all the orientation distributions in the w hole flow domain.

In order to accelerate the computation it could be interesting to solve the equation governing the evolution of the second order orientation equation instead of the resolution of the Fokker-Planck equation, along all the nodal trajectories, but as previously argued this approach involves the introduction of an appropriate closure relation allowing the expression of the fourth order orientation tensor from the expression of the second order one.

Remark 2. It is also possible to derive an evolution equation for the fourth order orientation tensor that involves the sixth order one, and also requires the definition of an appropriate closure relation. Despite that, this strategy could be more appealing because the constitutive equation directly uses the fourth order orientation tensor, and the procedure described in this work could be applied for deriving this closure relation, from now on we only consider the evolution problem related to the second order tensor and the associated closure problem. ((; , 0 ) ) () ()

in i i i xt x t t B t ψ ζ φ ζ = = =≈ = ∑ ( 42
)
where 0 (; , 0 ) xt x t= denotes the position at time t of a fluid particle located at 0 x at time 0 t = .

Taking into account the expressions of the second and fourth order orientation tensors, given by Eqs. ( 6) and ( 4) respectively, it results:

(2) 0 (2) (2) 1 ((; , 0 ) ) () () in i i i axtx t t a A t ζ ζ = = =≈ = ∑ (43)
and ( 4) 0 (4) (4) 1 ((; , 0 ) ) () ()

in i i i axtx t t a A t ζ ζ = = =≈ = ∑ (44)
where (2) a and (4) a denote the vector form of the second and fourth order orientation tensors respectively, being (2) i a and (4) i a the vector forms of the tensors defined by the integrals: We would mention that the value of () i φϕat point ϕ can be determined by standard interpolation from the nodal values of i φ contained in the vector i φ . Finally, matrix (2) A and (4) A are built from vectors

(2) i a and (4) i a w hi c h defi ne s thei r d iffe re nt co lumns .

Natural closure relation

Any closure approximation of the fourth order orientation tensor can be written using the Caley-Hamilton theorem as well as partial normalization and symmetry (see [START_REF] Dupret | Modelling the flow of fiber suspensions in narrows gaps[END_REF], in the general form:

6 (4) (4) i 1 3 (4) (4) i 1 in 3D in 2D i i i i i i aw aw β β = = = =  =     =   ∑ ∑ ( 4
=⊗   =⊗   =⊗   =⊗   =⊗   =⊗  ( 4 8 )
and where ( )

S
refers to the symmetric component.

All the coefficients i β are not independent due to the normalization conditions. Thus, we have in 3D From now on, and for the sake of generality, we consider the 3D case. Thus, Eq. ( 49) can be written in the matrix form:

DP P DP D D β β β β β β -               = - --+               -        (49) and in 2D 1 3 2 81 4 0 01 1 1 P β β β -     =- +        
ud fgh ββ =+ ( 5 1 )
Eq. ( 47) results:

1 1 6 (4) (4) i 1 0 =W i u i d i d d fg fh aw W W W W I W β ββ β β β - - = =     == = = + =         + ∑ &% (

Local natural closure relation

For computing the unknown beta coefficients we proceed by using a moving least square technique. For this purpose we define the approximation error

( ) 2 (2) (2) (2) (4) 1 (,) in ii i i i dd i Ja a a a W W ε βδ γ β = =  =- - -  ∑ & % (53)
where the weighting coefficient i γ scales with i λ r el ated to the ei ge nfunc tio n ()

( ) with

d d de d ε ε ε δ δ δ ε ε π  -   ==   (54)
Now, the optimality condition implies

( ) T (2) (2) (4) i 1 02 in ii i i i d i d J a aWa WW ε δ γ β β = = ∂  == - - -  ∂ ∑ & &% (55)
that is: = , it results a global least squares approximation that implies that closure coefficient does not depend anymore on (2) a .

( ) ( ) T (2) (2) i 1 T (2) (2) (4) i 1 in ii i d i in ii i i i aa W W aa W a W ε ε δ γ β δ γ = = = =  -=    =- -  ∑ ∑ & & & % (56) or ( ) ( ) 1 T (2) (2) i 1 T (2) (2) (4) i 1 in ii i d i in ii i i i aa W W aa W a W ε ε βδ γ δ γ - = = = =  =- ×    ×- -  ∑ ∑ && & % (

Numerical tests

Simple shear flow

We consider a 2D simple shear flow whose kinematics is defined by (, 0 )

T vy γ • =
being the shear rate 1 γ • = . The fiber suspension is characterized by fibers whose shape factor leads to 0.8 k = When the Fokker-Plank equation is integrated along a trajectory from the isotropic initial state until reaching the steady state characterized by a small enough variation between two consecutive time steps, that is:

6 11 0 () ( )1 0 ( ) ( ) nn tt t t ψψ ψ ψ - + -≤ -
and the Karhunen-Loève decomposition is applied to the snapshots selected during the orientation distribution evolution (around 20), the modes depicted in figure 1 are then obtained.

We can notice that the linear combination of only 5 orientation distributions (the ones depicted in figure 1) are enough to represent accurately the evolution of the orientation distribution along the considered streamline. It is also important to remark that these functions have not a full physical meaning, because, for example, they take negative values. However their linear combination can reproduce perfectly the evolution of the distribution function coming from the resolution of the Fokker-Planck equation. However, in the procedure described in the previous section the second moment of those distributions (second order orientation tensor) is calculated, and due to the characteristics of those distributions the derived tensors have not a unit trace. Despite this fact one could continue to apply the procedure proposed in the previous section ignoring those physical or conceptual difficulties.

Other possibility to avoid these modes without a full physical meaning, consists of assuming that the orientation distribution can be written everywhere on the considered streamline as a linear combination of some orientation distributions (snapshots), for example, those were used for applying the Karhunen-Loève decomposition. In order to avoid the consideration of two distributions too close, one could start with the initial distribution, adding a new distribution to the reduced approximation basis only when it is far enough from all the functions defining the reduced approximation basis. The advantage of this approach is that all the modes used for approximating the orientation distribution evolutions are real physical distributions and consequently the associated second order orientation tensors are perfectly defined from all points of view.

Figure 2 shows the different modes defined on the same streamline by using the just described strategy. We can notice the smoothness, periodicity, definite positive and normality of those distributions that define the reduced approximation basis which also contains 5 functions.

Finally, one could consider a last possibility related to the snapshot proper orthogonal decomposition. This technique is based on finding the significant modes from the application of the Karhunen-Loève decomposition, but assuming that those modes can be written as a linear combination of the snap N snapshot that were used to define the decomposition. The main advantage of this strategy derives from the fact that it is easy to prove that theses modes result from the eigenproblem defined by λ b ei ng the hi ghes t ei ge nv al ue . Fr om thos e e i genve ctor s the reduced approximation modes are computed using the fact that those modes are linear combination of the snapshots, i.e.:

[ ]

, i 1, , i i Qn φ =ϒ ∈ A ( 6 0 )
The most significant modes i φ related to the eigenvectors i

ϒ thro ugh i i Q φ =ϒ are depicted in figure 3.
In any case, these three numerical schemes offer three reduced approximation basis that can be used to define the natural closure according to the procedure defined in the previous section.

In order to compare these different strategies, we have computed the reference second order orientation tensor from the resolution of the Fokker-Planck equation, and then the different orientation tensors using the closures proposed in the previous section. The reference and approximated second order orientation tensors were compared at different times, for different values of the diffusion coefficients, using the following error measure:

ref aprox Ea a =- ( 6 1 )
Table 1 groups the different computed errors at the end of the evolution process. We can notice that when the diffusion coefficient is small enough the standard natural closure relation works perfectly because we are close to the hypothesis under which it was derived. However, when the diffusion coefficient is increased, the solution accuracy is degraded. On the other hand, the adaptive natural closure that w e have defined works perfectly for any diffusion coefficient, except for the too small values that induce numerical difficulties in the finite element discretization of the Fokker-Planck equation in the conformation space. We must mention that due to the fact that the Karhunen-Loève and the Snapshot Karhunen-Loève natural closures are built from modes (approximation functions) that have not a full physical meaning we have preferred using a global least square technique. Moreover, the global least squares fitting can be applied without a significant impact on the solution accuracy and with a considerable computing time reduction.

Extensional flow

We consider a 2D extensional flow w hose kinematics is defined by respectively, where only the snapshot natural closure is compared with the standard natural one. We can notice the superior accuracy of the snapshot natural closure for both diffusion coefficients.

Driven cavity flow model: an example of complex flow

In this section we consider the complex flow generated in a driven cavity . The velocity field is then solved by assuming a Newtonian behavior and by applying a standard mixed finite element formulation where the velocity and pressure approximations verify the LBB stability condition. The Fokker-Planck equation governing the evolution of the fiber orientation distribution function is then solved along some closed streamlines, where the periodicity condition of that distribution function was imposed as described in [START_REF] Ammar | A particle strategy for solving the Fokker-Planck equation governing the fiber orientation distribution in steady recirculating flows involving short fiber suspensions[END_REF]. From the computed orientation distribution function, the characteristic modes are extracted by using the technique previously described based on the application of the Karhunen-Loève decomposition, allowing to fit the empirical natural closure introduced in section 2.3 using the snapshot natural closure strategy previously introduced. Now, the evolution equation associated with the second order orientation tensor is solved by assuming different closure relations: linear, quadratic, hybrid, natural and the empirical natural one.

If the Fokker-Planck equation is solved in the whole domain considering a fine enough mesh, the fiber distribution function can be computed everywhere w ithout addressing any closure relation. Thus, the second order orientation tensor related to the just computed distribution function can be used as reference solution for comparison purposes.

The error between the computed second order orientation tensor and the reference solution coming from the solution of the Fokker-Planck equation is reported in Table 2.

From the previous results we can conclude that in the case of real complex flows it is difficult to conclude about the net superiority of a particular closure relation. In this kind of flow s there are regions in which advection terms dominate the diffusion ones (as for example in the neighborhood of the upper wall in the driven cavity flow) whereas in others regions diffusion terms dominate the advection ones. The contrast between these two competitive behaviors increases as the diffusion coefficient decreases because for high enough diffusion coefficients the spatial gradients of the orientation field are quite reduced.

Conclusions

This work focused on the resolution of the Fokker-Planck equation that governs the evolution of the fiber orientation distribution. To reduce the computing time, that equation is solved along some flow trajectories in order to extract the significant information of the solution from the application of the Karhunen-Loève decomposition. This extraction has been made using three different strategies: (i) the standard Karhunen-Loève decomposition; (ii) the simple use of some snapshots of the orientation distribution evolution (linear independence is checked during the basis construction) and (iii) applying the snapshot-Karhunen-Loève decomposition that makes optimal the snapshot basis. From this reduced approximation bases we adjust a closure relation that becomes optimal for such flow, that allows solving the evolution of some orientation moments which require a less amount of computation. The accuracy of the proposed approach has been illustrated in some numerical examples concerning some simple flows, and seems to be an appealing strategy when the salient diffusion effects have an unfavourable effect on the natural closure approximation.

Thus, a mixed strategy combining the natural closure in the regions where the orientation is dominated by advection terms and the adjustable one when the diffusion effects becomes important could be nice candidate to be considered when fast and accurate solutions are desired. [START_REF] Dupret | Modelling the flow of fiber suspensions in narrows gaps[END_REF] and the adaptive natural closure described in section 2 using the three different reduced approximation bases just defined: the first one using the modes coming from the direct application of the Karhunen-Loève decomposition, the second one using some appropriate snapshots and the last one based on the use of the snapshot Karhune-Loève decomposition. [START_REF] Dupret | Modelling the flow of fiber suspensions in narrows gaps[END_REF] and the snapshot natural closure described in section 2.

Diffusion

  's distribution, all the orientation probability is concentrated in the direction defined by ρ ˆ, and the corresponding orientation tensor results ρ

  a fluid volume located at point X at time t=0.

  denotes the nodal trajectory, l P the number of the selected distributions along each trajectory (snapshots), and L is the total number of nodal trajectories where the Fokker-Planck equation has been integrated, of course, L N << , with N the total number of nodes used to discretize the physical domain where the flow is defined and c N the number of nodes used in the conformation space discretization ( ϕ ).

Now

  the n approximation vectors contained in matrix B , that are expected to describe accurately the evolution of the orientation distribution solution along any flow streamline, i.e. 0 1

  () the tensor trace).

  most significant distribution functions (as explained in the previous section), one could be interested in determining the optimal value of d β , that is close to the procedure used by Dupret and Verleye to derive the natural closure relation, but they considered only the distribution related to an initial isotropic distribution of infinite slender fibers, without diffusion effects and in a simple shear flow.

  expression. In our simulations we have considered the following window function:

  57) which allows to compute at each (2) a the closure relation coefficients d β and then from Eq. (51)

.

  The orientation distribution (assumed 2D) is discretized using the SUPG finite element method on the mesh associated with the conformation space that consists of 200

  Figures 4 and 5 depict the time evolution of the error defined in Eq. (61) for

Figure 1 .

 1 Figure 1. Significant modes related to the evolution of the orientation distribution function along a streamline computed from the Karhunen-Loève decomposition.

Figure 2 .

 2 Figure 2. Reduced approximation modes related to the evolution of the orientation distribution function along a flow streamline.

Figure 3 .Figure 4 .Figure 5 .

 345 Figure 3. Significant modes related to the snapshot -Karhunen-Loève techniques.

  By multiplying the second equation by cosϕ , the first one by sinϕ and subtracting the first one from the second one, it results:

	(, ,) dx t dt ψϕ	,, xt dt ϕ d ϕ  (, ,) xt ψϕ ∂∂ ϕ ϕ  +=  ∂∂  	(, ,) x t ψϕ ϕ  r D ∂   ∂ 	(33)
	where	,, x t dt ϕ d ϕ	is computed from the Jeffery equation (13):
	-sin cos ϕ ϕ   	d dt ϕ	( 0 cos 12  12 0 sin k ϕϕ cos ) sin Ω  11 12 DD 12 11 11 12 DD 12 11 cos sin k DD DD ϕ ϕ ϕϕ cos sin cos sin ϕ ϕ ϕϕ   =+ -    -Ω -           -    -  	(34)
		(	sin			cos	)	0	12	12 0s i n cos
		+	(	sin	cos	)	11 12	12 11	cos sin

Table 1 .

 1 Error of the different closure approximations: standard natural closure

		Natural	KL -Natural	Snapshot	Snapshot-KL
		Closure	Closure	Natural Closure	Natural Closure
	0.1	0.0139	0.0095	0.0091	0.0355
	0.01	0.0046	0.0196	0.0113	0.0471

Table 2 .

 2 Error of the different closure approximations: linear, quadratic, hybrid, standard natural closure

	Diffusion	Linear	Quadratic	Hybrid	Natural	Empirical
						Natural
	1	0.0031	0.0704	0.0062	0.0019	0.0008
	0.1	0.1992	0.1516	0.0800	0.0353	0.0289
	0.01	1.1286	0.1452	0.1048	0.1216	0.1126


A Ammar

F Chinesta

Keywords: Model Reduction, Separated representation, Karhunen-Loève decomposition. Some Incipient Techniques For Improving Efficiency in Computational Mechanics

Some Incipient Techniques For Improving Efficiency in

Computational Mechanics

INTRODUCTION

Models encountered in computational physics and engineering, in particular computational mechanics, usually involve too many degrees of freedom, too many simulation time-steps, too many iterations (non-linear models, optimization or inverse identification), too many particles (in atomistic descriptions) or simply excessive simulation time (for example when simulation in real time is envisaged). Some examples of those scenarios are encountered in computational chemistry involving large time scale differences (reaction versus diffusion), direct numerical simulation of fluid flows (analysis of turbulence), kinetic theory descriptions of solid and fluids (including complex fluids) that usually involve multi-dimensional models with the associated curse of dimensionality problematic, virtual surgery, mechanical systems involving localized behavior, … In our former works (see [START_REF] Ammar | A New Family of Solvers for Some Classes of Multidimensional Partial Differential Equations Encountered in Kinetic Theory Modeling of Complex Fluids. Part II: Transient Simulation Using Space-Time Separated Representation[END_REF] and the references therein) some new strategies of model reduction and computational time saving have been proposed. In what follows the main ideas of those techniques are revisited.

GENERAL FORM OF THE PROBLEM

Let consider a general field Ψ involving two variables r and s (each one defined on a different space of dimension 1,2 or 3). Let consider two operators 1 and 2 respectively associated to the variables r and s.

(

) ( )

( , ) ( , ) 0 r s r s Ψ + Ψ = 1 2 (1)
For example, if we address the heat equation Ψ represents the temperature, r and s are respectively the time and the physical coordinates, 1 is the first order derivative with respect to time, and 2 is the Laplace operator. However, Eq. ( 1) can also represent more complex models. Thus, in steady state kinetic theory models Ψ represents the conformational distribution function, and r and s the physical and the conformational coordinates. In such models the first differential operator is defined in the physical space, accounting for the advective effects, and the second one is related to the induced microstructural evolution that sometimes is modeled from and advection-diffusion differential operator.

In other models such as the ones encountered in quantum mechanics, one find a Laplace operator defined in the conformation space that in this case is related to the physical space of each particle (electrons or nuclei). Thus, the resulting model results highly multidimensional, making difficult its numerical solution.

In the bead spring chain model (proposed to represent the mechanics of macromolecules) the variables involved are the connector vectors defining the contour length of the chain. The kinetic equation in this case involves, in addition to the general operators of equation ( 1), other "coupling" operators involving in a coupled way the derivatives with respect to different conformation coordinates.

For the sake of simplicity, from now on we are focusing on the solution of the generic equation ( 1).

MULTIDIMENSIONNAL APPROACH

The solution of equation ( 1) in the multidimensional case, in which a grid (or a mesh) of the domain Ω = Ω r UΩ s becomes too expensive or simply forbidden, needs for new advanced strategies. The existing strategies (sparse grids or stochastic techniques) allow addressing moderate multidimensional models. In the highly multidimensional case alternative strategies are needed. One possibility lies in defining separated representations and tensor product approximation basis:

( , ) ( ) ( )

i i i i r s F r G s α Ψ = 1 (2)
In order to build up this form, a series of projection and enrichment stages are considered.

The projection stage consists to find, for n given terms of the series the best set of the alphacoefficients associated to 1 1 , , , ,

n n F G F G 1
. The enrichment stage consists in finding the best functions ( ), ( ) R r S s with respect the problem variational formulation. The computer implementation of these steps needs for a discrete representation of all the functions: F i (r), G i (s), R(r) and S(s) using for example a 1D, 2D or 3D finite element interpolation in their associated domains Ω r and Ω s . We use for this purpose the vectors N and M containing the shape functions associated with each space. Finally F i , G i , R and S represent the nodal description of the associated functions. We define the following matrix related to the variational formulation equation ( 1

Ω Ω Ω Ω = Ω = Ω = Ω = Ω 2 2 2 2 MM NN M M N N 1 2 1 2 3 4 (3)
These integrals take into account the specific character of each operator. For example integration by parts is used in second order operators, upwinding for stabilizing advective terms, …

Projection stage

In this stage we consider the test functions given by * * ( , ) ( ) ( )

n i i i i r s F r G s α Ψ = 1 (4)
which writes in the discrete form:

* * ( , ) n T T i i i i r s α Ψ = 1 F MG N (5)
Using the discrete form of the Ψ function:

( , ) n T T i i i i r s α Ψ = 1 M F N G (6) the variational formulation of the equation (1) results * 1 1 0 n n i ij j i j H α α = = = 11 (7)
where

T T T T ij i j i j i j i j H = ⋅ + ⋅ F F G G F F G G 3 2 1 4
(8) This system must be solved taking into account some boundary or normality conditions. For example, we could consider a first function product verifying the boundary conditions, and then the other ones must vanish on the domain boundary (this condition is enforced in the enrichment stage). In this case we must enforce α 1 =1.

Enrichment stage

In this stage we consider the unknown field given by ( , ) ( ) ( ) ( ) ( )

i i i i r s F r G s R r S s α Ψ = + 1 (9)
being the test functions (12)

* * * ( , ) ( ) ( ) ( ) ( ) r s R r S s R r S s Ψ = + (
• For a given R (or for a test function in the form

of * * ( , ) ( ) ( ) r s R r S s Ψ = ) ( ) ( ) 1 T T n T T i j j j j i α = + = - + 1 R R R R S R F G R F G 3 2 1 4 3 2 1 4 (13)
It must be noticed that the number of the degrees of freedom involved in such non linear solution is the sum of the degrees of freedom involved in each discretization (instead the product resulting in mesh based discretization strategies). Finally, the new approximation functions F n+1 and G n+1 are obtained from R and S after normalization:

1 2 1 1 2 1 ( ) ( ) T n T n - + - + = = F R R R G S S S 1 2 (14)

ADAPTIVE MODEL REDUCTION

In this section the particular case of t

∂ = ∂ 1
, where t denotes the time, is considered. Let (t , x) be the couple of the two variables involved. In the finite element framework, Eq. ( 1) is usually solved within an incremental procedure. The unknown field at time t is noted by ( ) t x Ψ and is associated with the vector representation t 1 on Ω x . The weak form of the problem (1), using the notation introduced previously, writes

* * T T + = 1 1 1 1 2 2 2 4
(15) To obtain the evolution of the unknown field during the time an implicit or explicit time integration scheme could be used. But in general, the high number of degrees of freedom involved, makes difficult the problem solution at each time step when semi-implicit or fully-implicit schemes are considered. The main idea of the Karhunen-Loève decomposition is to represent the unknown field using a reduced number of functions expressed on the whole domain Ω x .

Karhunen-Loève décomposition

This decomposition lies in looking for the optimal functions φ(x) to represent the P discrete fields n 1 (n=1,.., P). These functions are obtained by maximizing 2 1

1 2

1 P N n i i n i N i i φ λ φ = = = 5 6 Ψ D E 9 A = 5 6 D E 9 A 1 1 1 ( 16
)
where N is the total number of nodes in the domain

Ω x .
The maximization of λ with respect to the functions φ is equivalent to enforce 0 δφ = , from which it results:

{ }

1 1 1 1
,

N N P N n n i j j i i i i j n i φ δφ λ φ δφ δφ = = = = F 5 6 Ψ Ψ = ∀ D E 9 A 1 1 1 1 whose matrix form writes , δφ φ λδφ φ δφ φ λφ = ∀ = T T k k (17) where 1 2 1 1 1 1 1 2 2 2 1 2
... ... , ...

3 4 Ψ Ψ Ψ 7 8 Ψ Ψ Ψ 7 8 = = 7 8 7 8 7 8 Ψ Ψ Ψ B C P P T P N N N k QQ Q 3 3 4 3 (18)

A posteriori model reduction

The time discretization of equation (15) writes:

1 n n + = 1 P 2 (19
) that in the case of using an explicit strategy ( )

n n = - P 1 2 4
(20) From the solutions n 1 one could extract the most significant eigenfunctions (the ones related to the highest eigenvalues): φ 1 , φ 2 ,... , φ r . These functions are stored into a matrix B which dimensions are (N,r) (in general r<<N). This matrix allow to describe the evolution of the unknown field using a reduced number of degrees of freedom (r instead N). Thus, writing

A priori model reduction

The main drawback of the a posteriori model reduction is the necessity to make at least one simulation using the fully finite element description. Besides, one could address the question related to the validity of the basis B. The purpose of the "a priori" model reduction [START_REF] Ammar | On the reduction of kinetic theory models related to finitely extensible dumbbells[END_REF] is the possibility of defining an adaptive strategy to build up the reduced basis during the model time evolution. Consequently, we must check the quality of the reduced solution at certain time steps. For this purpose we compute the residual, as proposed in [3]

1 1 n n n T n + + = -= - R 1 P Ba B P 2 2
(23) If the norm of this residual exceeds a given tolerance the basis must be enriched using some functions of the Krylov subspaces generated by the residual, and then the solution evolution is recomputed using the just updated approximation basis. This procedure continues until reaching convergence (that is, a small enough residual norm). It is obvious that successive enrichments could increase the number of the approximation functions. To avoid this unfavorable effect, a Karhunen-Loève decomposition is performed when the convergence of the enrichment procedure is reached.

TOWARDS A NEW HYBRID STRATEGY

This last section deals with models involving localized singularities or discontinuities, being the solution far from these singular regions smooth enough. We proved that evolving discontinuities do not accept efficient reduced descriptions. However, one expect that efficient simulations could be achieved by coupling reduced descriptions in the smooth regions and a fully finite element description in the zones involving singularities.

Let Ω 2 be the small region of interest where localized singularities exist and Ω 1 the remaining part of the domain, verifying

  weak form of the equation (1) becomes, in the framework of a fixed point strategy, in: • For a given S (or for a test function in the form of

A 1

 1 Each node of the discretization belongs to a single subdomain. The discrete finite element description + which are the respective degrees of freedom in the domains Ω 1 and Ω 2 . Thus the finite element problem could be the domain Ω 1 the "a priori" model reduction described in the previous section could be applied. A set of reduced basis function are incorporated in the matrix B. A model reduction over Ω 1 is made thanks to the


Cristian Ruiz

Olivier Richard

Oleg Iegorov

Brice Videau

Managing Large Scale Experiments in Distributed Testbeds

Keywords: Grid Computing, Distributed Application, Experiment conduction, Testbed, Experiment methodology Expérimentation, Expériences à grande échelle, Plates-formes d'expérimentation, Méthodologie d'expérimentation

Performing experiments that involve a large amount of resources or a complex conguration, proves to be a hard task. We present Expo, which is a tool for conducting experiments on distributed platforms. First, the tool is described along with the concepts of resource and task sets, which abstracts away some of the complexity in the manage of resources and software execution. Next, the tool is compared with other similar solutions based on some qualitative criteria, scalability and expressiveness tests as well as the feedback coming from using dedicated testbeds. The report nishes with the evaluation of Expo scalability and some uses cases in Grid'5000 and PlanetLab. Our experience showed that Expo is a promising tool to help the user with two primary concerns: how to perform a large scale experiment eciently and easily, together with its reproducilibity.

Introduction

As the software to perform simulations have improved in recent years, there is still the need to test and evaluate the software in real distributed infrastructures. Moreover, the option of experimental evaluation of an algorithm has been encouraged as an approach complementary to the theoretical evaluation [START_REF] Johnson | A theoretician's guide to the experimental analysis of algorithms[END_REF]. In order to addres limitations such as, software reconguration, lack of the control and access to monitoring systems testbeds were created [START_REF] Lanteri | Grid'5000: a large scale and highly recongurable grid experimental testbed[END_REF]. A testbed is a platform for experimentation of large distributed applications. It is sometimes shielded from the instabilities of production environments, and allows users to test particular modules of their applications in an isolated fashion. Some examples of testbeds are: PlanetLab [START_REF] Peterson | Experiences building planetlab[END_REF], Emulab [START_REF] White | An integrated experimental environment for distributed systems and networks[END_REF], GENI 1 , Grid5000 [START_REF] Bolze | Grid'5000: a large scale and highly recongurable experimental Grid testbed[END_REF] and ORBIT [START_REF] Ott | Orbit testbed software architecture: Supporting experiments as a service[END_REF]. While these platforms oer more stability and control over resources than grids, controlling, deploying and running applications on them is still a hard task. That is the reason why some tools have been developed to cope with the problems encountered when researches try to perform experiments involving a large amount of resources or a complex conguration. The main aspects those tools help the user with, are: (1) the description of the experiment, (2) the control and access to the resources, (3) task orchestration, (4) software deployment, [START_REF] Minyard | Gridshell: A login shell for orchestrating and coordinating applications in a grid enabled environment[END_REF] monitoring and the collection of results. In more detail a number of tasks must be completed before an experiment can be actually started. These tasks include resource discovery and acquisition as well as the deployment of the necessary software. Once the application is launched, its execution must be controlled, and once it is nished all the output must be collected. As a result, most of the applications on the testbeds are run in an ad-hoc, application-specic manner. This method may match the current requirements of experiments, but fails with the scale, heterogeneity, and dynamism of highly distributed systems. We can classify the concerns when performing experiments in distributed infrastructures in two categories: how to deploy software and manage the nodes and how to provide a mechanism to describe and log the experiment activity in order to render this process more reproducible, which is the cornerstone of the scientic method [START_REF] Dror | Experimental computer science: The need for a cultural change[END_REF]. Thereby the need of a generic tool to help users conduct their experiments in distributed environments is evident. In this paper we present Expo, already introduced in [START_REF] Videau | Toward an experiment engine for lightweight grids[END_REF][START_REF] Videau | Expo : un moteur de conduite d'experiences pour plates-forme dedies[END_REF], which is an experiment engine for distributed platforms. It helps users to conduct their experiments on any set of machines with the only requirement of remote shell access. We present its main advantages and compare it with others experiment control systems, as well as describing several uses cases and its scalability. The structure of this report is as follows: the next section shows the state of the art where related work is presented, regarding tools that an experimenter can use in order to conduct experiments in dierent testbeds. After in section 3 Expo is presented in depth with its principles and advantages, some cases of use are shown in Section 4 in two dierent testbeds. Results and comparison with other experiment tools are presented in Section 5 and nally section 6 presents the conclusions and future works.

State of the Art

There are many tools that researches can use in order to experiment with computational infrastructures. Here we recapitulate dierent complementary approaches passing from the testbeds managers, which are tools aimed to help the user with the experiment process, to other tools, which even though there are not directly related to the experimentation process, they share some similar concerns.

Testbed Managers

They are aimed at helping the user in the experimenting process, automating certain tasks. Sometimes they are tied to particular testbeds. Gush (GENI User Shell) -previously known as Plush [START_REF] Albrecht | Remote control: distributed application conguration, management, and visualization with plush[END_REF] (PlanetLab User Shell) -is a widely used tool for application management on PlanetLab and GENI testbeds. It helps users to nd desirable resources on the target platform, to prepare those resources with necessary software, to run the experiment, and to provide the maintenance while the experiment is running. With Gush users describe their experiments in a XML le, which consists of a set of building blocks. These blocks describe all the aspects of an experiment: software packages to install, resources to use and the main logic of an application. For partial tasks, such as resource discovery or remote software installation, Gush uses specic PlanetLab tools (SWORD, Stork, etc). As PlanetLab currently consists of 1133 nodes running at 544 sites all over the world, Gush's fault tolerance and scalability become extremely important features. Gush handles three types of failures: process failures, remote host failures, and controller failures (which is the responsible for managing the Gush client processes). OMF (cOntorl and Management Framework) [START_REF] Thierry Rakotoarivelo | Omf: a control and management framework for networking testbeds[END_REF] is a testbed Control, Measurement and Management Framework. This framework is used in dierent wireless testbeds around the world and also in Planetlab. Its architecture versatility aimed at federation of testbeds, conceived mainly for the testing of network protocols and algorithms in wireless infrastructures. The OMF architecture consist of 3 logical planes: Control, Measurement, and Management. Those planes provide the user with tools to develop, orchestrate, instrument and collect results as well as the tools to interact with the testbed services.

It uses a comprehensive domain specic language based on Ruby to provide experimentspecic commands and statements. This allows an experimenter to write an experiment description, which details the resource requirements, their conguration, the denition of the events to capture and the respective actions to trigger when performing the experiment. They propose an approach in which an experiment will be similar to a meta package. It will depend on many other packages, containing the experiment description, applications, topologies and measurements.

The aim of the system in [START_REF] Eide | An experimentation workbench for replayable networking research[END_REF] is to extend Emulab capacities to support replayable research. They propose a new model for experimenting with testbeds which overcomes several problems encountered when interacting with Emulab. This model is composed mainly of two parts: templates and records. A template describes all the testbed-based environment, comprises every conguration le, source code, input les and any data generated or needed by the application. These templates are versioned, therefore they can suer revision every time there is some change, which could be used for exploring dierent directions of research. A record is the persistent account of the activities and results that occurred within a run. This could be the data generated, all the raw intermediate data as well as the software, scripts, etc.

Additionally, there is a dynamic part that encompasses the activities of the experiments and how they are going to be scheduled when the experiment is executed as well as the actions triggered by the user or by special conditions. The user starts by describing a template, which denes the network topology and other experiment information. Then this template is commited to the system. The user next proceeds instanciating the template. When a template instance is created, testbed resources are allocated, congured, and booted. After the network and devices are up and running, the workbench automatically starts a run and starts any prescheduled activities withing that run.

2.2

Tools aimed to algorithm testing Splay [START_REF] Leonini | Splay: distributed systems evaluation made simple (or how to turn ideas into live systems in a breeze)[END_REF] is an integrated system, which covers the complete chain of distributed system design.

It is a unique tool allowing to develop, deploy and maintain an application in a distributed environment. Splay can be used on a classical testbed (PlanetLab, ModelNet, Grid5000), on a non-dedicated platform (such a network of idle workstations), in a normal cluster, and even on several testbeds simultaneously. Splay requires applications to be written in a Domain Specic Language based on Lua programming language. As Splay puts a particular emphasis on security, such applications are executed in sandboxed environments with minimal underlying operating system interaction. It uses a client/server architecture. The user can start his application on the server host or on some other machine, which can access the server. The server controls its execution and responds to failures. One of the most interesting features of Splay is its churn manager. It allows to reproduce a given live experiment with constant leaving and returning of participating nodes. The particular churn can be obtained from known analytical studies or publicly available churn scenarios.

Weevil [START_REF] Wang | Automating experimentation on distributed testbeds[END_REF] is a framework for automating experiments with distributed systems. It allows to generate application workloads based on a statistical model or on a specic system usage scenario programmed by the user. The latter feature makes Weevil a unique tool. After a specic workload is generated, the user can deploy his application on a predened set of remote hosts, run it and collect the results. Thus, the work with Weevil consists of two stages: Workload generation and application deployment.

In [START_REF] Broto | Autonomic management policy specication in tune[END_REF] a tool to manage software in distributed infrastructures is presented. It is based on the Concept of Autonomous computing in order to make the administration of an infrastructure as a component architecture. The main idea is then to automatically create a representation based on fractal components of the real system, with two main parts: Application components and platform components. All expressed with a subset of UML diagram. It has already been used in the installation of a cluster software and the deployment of a TLM electromagnetic simulation code in a grid infrastructure [START_REF] Fadi | Multi-scale modeling: from electromagnetism to grid[END_REF].

Continous integration

Encompases dierent principles aiming at the rapid test and deploy of changes in software. It could be cited tools such as capistrano2 , which automate the process of server administration, mainly the deployment of web application using the Ruby and rails framework. NMI [START_REF] Pavlo | The nmi build & test laboratory: continuous integration framework for distributed computing software[END_REF] is a framework to build and test software in a heterogeneous, multi-user, distributed computing environment. The principal aim is to oer to the user the continually testing of software changes. Because bugs in a software have to be xed early in the development process. The user describe the process of building and testing along with its external software dependencies by using a lightweight declarative syntax. It takes advantage of a batch system (Condor) in order to deploy the software to be tested into the distributed infrastructure. Working with condor gives some capabilities to the system such as: matchmaking, fault tolerance, grid resource access, resource control, authentication and le transfer. It works along with a versioning system, so as to log the results and changes as well as performing the tracking of all input so that to ensure the repeatability and reproducibility of tests.

Scientic Workows

These tools allow the experimenter to describe the complete ow of the experiment and the dependency between tasks and data generated. But the more important concept is the sharing of analysis and information through the composition and execution of workows. One example in this category is Taverna [START_REF] Oinn | Taverna: lessons in creating a workow environment for the life sciences[END_REF].

Parallel Executers

In this category we can cite projects like GridShell [START_REF] Minyard | Gridshell: A login shell for orchestrating and coordinating applications in a grid enabled environment[END_REF] and TakTuk [START_REF] Claudel | Taktuk, adaptive deployment of remote executions[END_REF], which help the user with the deployment of commands in large distributed infrastructures.

Expo

Expo [START_REF] Videau | Toward an experiment engine for lightweight grids[END_REF][START_REF] Videau | Expo : un moteur de conduite d'experiences pour plates-forme dedies[END_REF] 3 is a tool for conducting experiments on distributed platforms. It is written in Ruby and it has a modular structure, which allows to replace its components to add some new functionality. It is based on a client/server architecture. This architecture was designed having in mind that every time the lost of connection can occur, so that an asynchronous execution can be guaranteed. Communications between clients and servers are based on the SOAP protocol.

Currently, Expo has been tested in Grid5000 and PlanetLab testbeds and uses its tools to perform subtasks. The core functionality of Expo consists in taking commands from the user and execute them eciently on a large set of heterogeneous distributed resources. This eciency is achieved thanks to parallel execution using TakTuk [START_REF] Claudel | Taktuk, adaptive deployment of remote executions[END_REF], which spread itself and form a tree, improving the scalability and enabling the exploitation of the hierarchy structure of the platform. Expo users describe their experiments using Domain Specic Language (DSL) based on Ruby. As a result of the feedback obtained from Grid5000, we wanted to focus on the deployment and ecient execution of the experiment. That is why Expo oers the possibility of executing any software, as long as such software can be executed on the target machine. It interacts with the services of the Grid5000 testbed such as its API 4 in order to provide an easy interaction with Kadeploy [START_REF] Georgiou | A Tool for Environment Deployment in Clusters and light Grids[END_REF], which allows the user to deploy complete operating systems. It is not in any way tied to a specic platform, because of its modular architecture new modules can be added in order to interact with other testbeds or platforms. This was proved by developing the module for interacting with PlanetLab, which did not require any change in the subjacent components.

Expo components

Expo is the result of the interaction of dierent components that oer several functionalities. This interaction is shown in Figure 1. It was developed having in mind a modular architecture, each one of the component is described below:

Command Executer: It gives verbosity to the commands, getting all the command status and giving a low level control execution.

Experiment Controller: This part keeps track of the whole experiment run, controlling all the information related to an experiment. It logs all the commands executed, storing information such as: time of execution, output, errors etc. It enables to unify the tasks and resources into a logical unit: the experiment, making possible the addition of tasks or resources. The client part enables to save all the information into a structured data format for a postmortem study. Interface: This part could be seen as the set of interacting methods. An interactive mode driven by the following reasons: (1) an important part of the experiments are interactive, [START_REF] Bolze | Grid'5000: a large scale and highly recongurable experimental Grid testbed[END_REF] the writing of an experiment description le is a trial-and-error process which involve using dierent parameters, congurations and ows of control, (3) An interctive environment lets scientists look at data, test new ideas, combine algorithmic approaches, and evaluate their outcome directly [START_REF] Pérez | IPython: a System for Interactive Scientic Computing[END_REF]. This approach is already used by dierent scientic environments based on Python such as: IPython and Scipy [START_REF] Jones | SciPy: Open source scientic tools for Python[END_REF]. A Stand alone mode is also available, which execute the experiment description le without any user intervention. synchronous and asynchronous execution of commands. It allows the user to access the information about the nished commands or those that are running thanks to a unique identier. The information available are: Command status, start and end dates of the command, return status, stdout and stderr.

Later, it is shown the interface to all these components by presenting the language for describing the experiment. The features that control the ow of execution and the commands provided by Expo are described in [START_REF] Videau | Expo : un moteur de conduite d'experiences pour plates-forme dedies[END_REF].

Abstractions

In this subsection the abstractions are presented in depth, which is one of the more important features of Expo and some operations that have been added in order to ease and improve the exibility of the experiment description. These abstractions are implemented as objects in Ruby and have some operations, properties and iterators. They resemble other existing objects in Ruby.

Inria

Lyon

Resources

Expo introduces the notion of ResourceSet. These ResourceSet add resources into a logical unit and associate properties to them. For instance, we can gather together the nodes from a same cluster associating them the same frontend, as well as the same physical properties if the cluster is homogeneous. This information is actually used in order to perform the ecient deployment of commands. These sets of resources could suer diverse transformations using some operations and properties summarized in Table 1.

Figure 2 shows how the resourceSet object can be used to exploit the hierarchy of the infrastructure in Grid5000. We can divide the resources belonging to the same site as well as separate them per cluster. This can also be applied for the Planetlab testbed e.g. the resourceSet can have information about the location of the resources for the same country or site. In other cases it can be used to dene complex congurations as in the case we need to deploy an infrastructure where dierent nodes have dierent roles. and their interactions.

Tasks

A task associate a command with a set of resources. Therefore, dierent mappings between commands and resources can be expressed and easily managed. It can either encapsulates the parallel or the sequential execution of commands. Also it is possible to create a set of tasks, each one with dierent characteristics and resources. On those set of tasks, it is possible to control the way tasks are executed such as synchronously or asynchronously. This class encapsulates the results, allowing the user to the postmortem treatment of the execution of the experiment. Obviously the gathering of this information has a cost, but it does not inuence the self execution of the experiment. Among its features, it wraps up in ruby objects dierent functions normally used in an experiment, in order, for example, to measure time or getting some statistics. Making easier for the user the treatment of his/her experiment.

Some Uses cases

The aim of this section is to show how the experiments are described with Expo. The dierent functionalities that provide the user with the exibility for describing his/her experiment. Additionally, this section shows the interaction with two dierent testbeds, Grid5000 and PlanetLab.

In Figure 3 two examples of experiment description les are shown. One experiment description is an experiment performed in Grid5000, which consisted in the execution of an electromagnetic simulation over dierent sites using one node per site. We deployed an OS image with all the software already installed through the use of Grid5000 API that interacts with Kadeploy. The specication of the corresponding image to deploy is indicated as a parameter in the function that request the resources, which it is shown in the rst lines of the le. It is important to remark that the specication of the resources to use has a lot of parameters, because we needed this especial conguration but most of the parameters shown can be ignored and the default values will be used. Moreover, in the le we can see some Expo operators to ease the procedure of execution of a program using MPI such as generating a secure ssh communication and generating the correct hostle. Next all the inniban interfaces are deactivated, which is easily expressed with Expo. The application is executed and we got the results of the execution in two variables that we can treat with dierent operators and methods already dened in them. Finally we gather the les generated by the execution from the nodes with just one command. The other experiment performed over Planetlab consisted in executing the Linux command "hostname" in all the nodes of the slice, and to count how many of them reply. This information is put in a le that can be used to plot the availability of the nodes in the slice. It is a simple example that show how an experiment is carried out in Planetlab using Expo. More examples can be found on the Expo website 5 r e q u i r e ' g5k_api ' b a s e _ u r l = " h t t p : / / p u b l i c . g r e n o b l e . g r i d 5 0 0 0 . f r /~r u i z s a n / " environment = " t l m _ s i m u l a t i o n . env " g 5 k _ i n i t ( : s i t e => [ " nancy " , " r e n n e s " , " l i l l e " , " g r e n o b l e " , " s o p h i a " ] ,

: r e s o u r c e s => [ " nodes=1" ] , : environment => { b a s e _ u r l+environment => 5 } , : w a l l t i m e => 3 6 0 0 , : deployment_max_attempts => 3 , : name => "TLM_code" ) g5k_run $ a l l . gen_keys copy $ a l l . n o d e _ f i l e , $ a l l . f i r s t , : path=> " / r o o t / nodes . d e p l o y e d " ### d e a c t i v a t i n g t h e i b 0 i n t e r f a c e $ a l l . each { | node | t a s k node , " / s b i n / i f c o n f i g i b 0 down" } ### E x e c u t i o n o f t h e s i m m u l a t i o n with MPI id , r e s = t a s k $ a l l . f i r s t , " . / exec_tlm 1 369 192 510 250 1 s c " g e t _ r e s u l t s ( $ a l l , " /TLMME_multimode/ p r o f i l e . * " , "~/ p r o f i l e s " ) (a) TLM experiment description. r e q u i r e ' e x p o _ p l a n e t l a b ' ## g e t t i n g r e s o u r c e s form p l a n e t l a b S l i c e g e t _ r e s o u r c e s task_mon=Task : : new ( " hostname " , $ a l l , " M o n i t o r i n g " ) F i l e . open ( " P l a n e t l a b _ a v a i l .

Evaluation

As our platform of study is Grid5000, we took most of the feedback from there and the use cases. Some percentage of the experiments use jobs, which run interactive, some are batch and some need the deployment of an image as seen in Table 2. Some of the diculties that the user has to deal with are:

Managing heterogeneous clusters, varying the number of sites, clusters and computational nodes. r e q u i r e ' g5k_api ' r e q u i r e ' benchmark ' g 5 k _ i n i t ( : s i t e => [ " nancy " , " s o p h i a " ] , : r e s o u r c e s => [ " nodes =200" , " nodes =100" ] , : w a l l t i m e => 2 0 0 0 , : name => " Big_experiment " ) g5k_run # run t h e r e s e r v a t i o n s i z e s = [ 1 0 , 5 0 , 1 0 0 , 2 0 0 , 3 0 0 ] p u t s " nb_nodes time " ### i t e r a t e s o v e r t h e d i f f e r e n t s s e t s t o t r y s i z e s . each { | n | nodes=$ a l l . uniq [ 0 . . ( n -1)] task_mon=Task : : new ( " hostname " , nodes , " T e s t i n g " ) r e s =0 4: Comparison between Experiment description les: These les were used in the evaluation of the scalability of the two tools. It should be noticed here that the experiment description for Gush has to be changed every time we need to change the number of nodes to try with. Also that Gush needs a le for the resource description that is not shown. network experimentation. Unlike Expo which requires only the presence of Ruby and Perl and some modules in order to work. However, the language provided by OMF for describing the experiments oers a big exibility and uses the same approach as Expo. Splay was conceived to test algorithms and needs everything to be implemented in the Lua language, which makes it dicult to use it for testing and deploying existing software. Workbench is completely bounded to Emulab and does not support the manage of dierent testbeds. In contrast to almost all the tools, Gush shares many features with Expo as the ease of installation and the capacity of being adaptable to dierent testbeds and that is the reason why we evaluate both in more depth.

The evaluation consisted in the expressiveness of the language, the performance and scalability of the command execution. The comparison between both tools was done by carrying out an experiment, which involved a large amount of nodes. We dened an experiment that consisted in executing a command in a set of resources and measuring the time elapsed, while varying the number of nodes. Therefore, we compare the time to execute the commands and the exibility in the description of the experiment.

Figure 4 shows the descriptions of the experiment used for Gush and Expo. We can notice looking at the experiment description that for Gush either we have to change the le for each experiment so as to try dierent number of resources, or we can create a long description le with all the possibilities we want to try. This is not the case for Expo, which uses Ruby and provides a programmatically approach for describing the experiment allowing it to be exible enough to adapt to the normal activities or changes when we perform an experiment.

In Figure 5 the scalability of the mechanism for the execution of commands was evaluated. In this gure, we can see that Expo outperforms Gush due to the use of TakTuk parallel executer, also that Expo presents less variability in the time to execute the experiment, which is important to the reproducibility. It was noticed as well that when we tried to execute an experiment with more than 400 nodes, problems arise trying to perform it with Gush. [START_REF] Eide | An experimentation workbench for replayable networking research[END_REF] Conclusions and Future Works Experimentation in computer science and specially in distributed infrastructures has seen the emergence of dierent experiment control systems and from those we can draw as a conclusion that most of the tools distinguish almost the same phases in the experimenting process. The parts of the experiment process that a tool must control and help the user with are mainly three: (i) the control, (ii) the supervision and (iii) the management of the experiment. The rst part comprises the description of the experiment, the capture of data, the denition of the source of data and how to get it after the experiment has nished as well as the ow of control of the experiment. This is an important step for the reproducibility of the experiment. Second, the experiment supervision, which means the monitoring of the experiment. As a last phase the experiment management, i.e. the interaction with the platform, how to take advantage of the services provided by the infrastructure in order to carry out the experiment.

Expo oers a way to describe the experiment by using a programming language providing a lot of exibility and more importantly the abstractions that allows the user to express complex congurations. We had special attention in automating the typical tasks done when an experiment is performed. Because we think that automating the experimentation process is the way to go, being one of steps that will lead to the experiment reproducibility. Furthermore it is important to encourage the culture of experiment reproducibility, which is acknowledged to be a shortcoming in computer experimentation. The use of experiment tools will save user time, which can be spent in improving the software itself, it will save costs and allow others to reproduce the results more easily. It is important to integrate some features to Expo for the sake of reproducibility. We need to improve the part of the system that logs the experiment run in order to have detailed and easy to treat information that would enable a possible replay of the experiment. Incorporate some mechanisms to monitor and to generate a workload and more importantly to deal with fails.

Inria

Number of Nodes

Figure 1 :

 1 Figure 1: Expo Component Diagram

  t e => [ " l y o n " , " r e n n e s " ] , : r e s o u r c e s => [ " nodes=12" , " nodes=6" ] , : w a l l t i m e => 7200 ) g5k_run ## A l l The r e s o u r c e r e s e r v e d $ a l l $ a l l . gateway ## New r e s o u r c e s e t ##with o n l y r e s o u r c e s o f l y o n c l u s t e r lyon_nodes=$ a l l [ " l y o n " ] lyon_nodes . gateway (b) Expo Experiment description File.

Figure 2 :

 2 Figure 2: Example of exploiting the hierarchy in Grid5000 using Expo. Here is presented grapically the concept of resource Set and how it looks like in the code of the experiment description le.

  t x t " , 'w+ ' ) { | f | r e s=n i l f . p u t s " Date Time Num_Res" 2 4 0 . t i m e s { data_me=Time : : now . to_i id , r e s = task_mon . e x e c u t e time=r e s . t o t a l _ d u r a t i o n f . p u t s "#{data_me} #{time } #{ r e s . l e n g t h } " f . f l u s h s l e e p ( 6 0 ) } } (b) PlanetLab experiment.

Figure 3 :

 3 Figure 3: This show two examples of experiment description les used in dierent testbeds, Grid5000 and Planetlab. It is important to notice that we just need to load another module in order change the platform.

<

  ?xml v e r s i o n=" 1 . 0 " e n c o d i n g=" u t f -8"?> <gush> <p r o j e c t name=" T e s t i n g overhead "> <component name=" C l u s t e r 1 "> <r s p e c > <num_hosts>20</num_hosts> </r s p e c > <r e s o u r c e s > <r e s o u r c e type=" s s h " group=" l o c a l "/> </ r e s o u r c e s > </component> <e x p e r i m e n t name=" s i m p l e "> <e x e c u t i o n > <component_block name=" cb1 "> <component name=" C l u s t e r 1 "/> <p r o c e s s _ b l o c k name="p2"> <p r o c e s s name=" t e s t "> <path>hostname </path> <cmdline> <arg ></arg> </cmdline> </ p r o c e s s > </p r o c e s s _ b l o c k > </component_block> </e x e c u t i o n > </experiment > </ p r o j e c t > </gush> (a) Gush Experiment Description.

  2 0 . t i m e s { id , r e s = task_mon . e x e c u t e p u t s "#{ r e s . l e n g t h } #{ r e s . t o t a l _ d u r a t i o n } " } } (b) Expo Experiment Description.

Figure

  Figure 4: Comparison between Experiment description les: These les were used in the evaluation of the scalability of the two tools. It should be noticed here that the experiment description for Gush has to be changed every time we need to change the number of nodes to try with. Also that Gush needs a le for the resource description that is not shown.

TimeFigure 5 :

 5 Figure 5: Evaluation of the scalability of Gush and Expo when executing a command in a large set of resources. The red points represent outliers.

  Reservation System: The request of resources is perfom by this component, using the Command Executer and additionally, it associates all the reservation information to the experiment.

		Uses	Interface	Request Experiment Status	Experiment Controler
	Actor	Request	Hosts	cmd & Log	status
			API Testbed		Command Control
		Request	Hosts	cmd	status
	Resource Broker	Hosts Request	Reservation System		Command Executer

RR n°8106

http://www.geni.net Inria

https://github.com/capistrano/capistrano/wiki Inria

http://expo.gforge.inria.fr/

http://www.grid5000.fr/mediawiki/index.php/API RR n°8106

http://expo.gforge.inria.fr/ Inria

Cristian Ruiz, Olivier Richard, Oleg Iegorov, Brice Videau

Acknowledgments

Experiments presented in this paper were carried out using the Grid'5000 experimental testbed, being developed under the INRIA ALADDIN development action with support from CNRS, RE-NATER and several Universities as well as other funding bodies (see https://www.grid5000.fr).

Choosing and keeping track of the sites used.

Deploying software.

Easy control of nodes.

Running experiments at very large scale with more than 1000 nodes.

We used those diculties as criteria for comparing the dierent experiment control systems found and to justify our design choices in the development of Expo.

Evaluation of the experiment control systems

In section 2, we presented some of the tools that an experimenter can take advantage from. In this section, we focus on the tools that were conceived to help the user with the experiment cycle in distributed systems. The evaluation start by comparing those tools. In Table 3, we show the dierent criteria used to compare them. We can see that some of them favor the reproducibility of experiments like the capacity of workload generation, fault injection and the use of a versioning system. Others are platform dependent and we can see also some important characteristic such as the ease of use. The aim of this evaluation is to place our tool regarding all the existing experimental frameworks as well as looking for characteristic that are worth to be taken into account for Expo. Weevil regardless of its capacity of generating a workload and helping the user with the creation of scripts, it uses for describing the experiment a language used to create conguration les for complex software, which few users are aware of. OMF installation requires several nodes to perform a simple installation and it is more oriented to the Inria


Crustal structure of Guadeloupe Islands and the Lesser Antilles Arc from a new gravity and magnetic synthesis

Keywords: Lesser Antilles Arc, Guadeloupe Island, geophysical surveys, geophysical data synthesis, magnetic and gravity anomalies, modeling, seismic profiles, subsidence, deformation, crustal structures

Guadeloupe Island (West French Indies) is one of the twenty islands that compose the Lesser Antilles Arc, which results from the subduction of the Atlantic Ocean plate beneath the Caribbean one. The island lies in a complex volcano-tectonic system and the need to understand its geological context has led to numerous on-and offshore geophysical investigations. This work presents the compilation and processing of available, on-land, airborne and marine, gravity and magnetic data acquired during the last 40 years on Guadeloupe Islands and at the scale of the Lesser Antilles Arc. The overall dataset provides new Bouguer and reduced to the pole magnetic anomaly maps at the highest achievable resolution. Regionally, the main central negative gravity trend of the arc allows defining two subsident areas. The first one is parallel to the arc direction (~N160°E) to the north, whereas the second unexpected southern one is oriented parallel to oceanic ridges (N130°E). Along the Outer Arc, the long wavelength positive anomaly is interpreted, at least along the Karukera Spur, as an up-rise of the volcanic basement in agreement with the seismic studies. To the NE of Guadeloupe, the detailed analysis of the geophysical anomalies outlines a series of structural discontinuities consistent with the main bathymetric morphologies, and in continuity of the main fault systems already reported in this area. Based on geophysical evidences, this large scale deformation and faulting of the Outer Arc presumably primarily affects the Atlantic subducting plate and secondarily deforms the upper Caribbean plate and the accretion prism. At the scale of Guadeloupe Island, joined gravity and magnetic modeling has been initiated based on existing interpretation of old seismic refraction profiles, with a general structure in three main layers. According to our geophysical anomalies, additional local structures are also modeled in agreement with geological observations: i) the gravity and magnetic signals confirm an up-rise of the volcanic basement below the limestone platforms outcropping on Grande-Terre Island ; ii) the ancient volcanic complexes of Basse-Terre Island are modeled with high density and reverse magnetized formations; iii) the recent volcanic centre is associated with formations consistent with the low measured density and the underlying hydrothermal system. The EW models coherently image a NNW-SSE depression structure in half-graben beneath Basse-Terre Island, its western scarp following the arc direction in agreement with bathymetric and seismic studies to the north of the island. The so-defined depressed area, and particularly its opening in half-graben toward the SW, is interpreted as the present-day front of deformation of the upper plate associated with the recent volcanic activity on and around Guadeloupe. Based on this regional deformation model, perspectives are given for further integrated investigation of key targets to address the internal structure and evolution of the Lesser Antilles Arc and Guadeloupe volcanic system.