Carine Lucas
email: carine.lucas@univ-orleans.fr

Christine Kazantsev
email: christine.kazantsev@imag.fr

MULTISCALE APPROXIMATIONS FOR OCEAN EQUATIONS: THEORETICAL AND NUMERICAL RESULTS. *

Keywords: Multiscale, existence of solutions, numerical results, boundary layer. subject classifications. 34E13, 76D10, 35A05, 65N06

We study a stationary Quasi-Geostrophic type equation in one or two dimensional spaces, with a quickly varying topography. We consider an asymptotic expansion of this equation on several space and time scales. At each expansion's order, we split the approximated solution into an interior function, which represents the solution far from the western boundary, and a corrector function that takes into account the boundary layer. We derive the systems at each order for the two functions and prove mathematical properties on these systems. Then we present numerical tricks and results, with and without topography, in the one and two dimensional cases. The method is very efficient compared to classical ones (finite differences, finite elements) which are very expensive due to the quickly varying topography and thin boundary layer.

Introduction.

Introduced by Charney [START_REF] Charney | On the scale of atmospheric motions[END_REF] in 1948 for the atmospheric motion, used in 1950 by Charney et al. [START_REF] Charney | Numerical integration of the barotropic vorticity equation[END_REF] for the first computer weather forecasts on the ENIAC computer, the Quasi-Geostrophic potential vorticity equation (QG equation) is adapted to the motion of a wind-driven ocean by Bryan [START_REF] Bryan | A Numerical Investigation of a Nonlinear Model of a Wind-Driven Ocean[END_REF] in 1963. During the 70s, more realistic models are developed, like multilayer quasigeostrophic models (see [START_REF] Holland | The Role of Mesoscale Eddies in the General Circulation of the OceanNumerical Experiments Using a Wind-Driven Quasi-Geostrophic Model[END_REF]). Thanks to more powerful computers, primitive equations have replaced them for ocean forecasts. But QG models are still very useful simplified models, especially for preliminary studies. In this paper, we are interested in the stationary QG equation in a domain with an ε-periodic bathymetry. Numerically, two problems occur: the first one is to resolve correctly the boundary layer, that implies to work with a rather small space scale near the western boundary. The second difficulty is due to the ε-periodic bathymetry, which imposes a very small space scale everywhere. The cost of usual methods (finite difference, finite elements) explodes. In this work, we follow the multiscale approach of R. Klein, E. Mikusky, A. Owinoh [START_REF] Klein | Multiple Scales Asymptotics for Atmospheric Flows[END_REF], P. Ailliot, E. Frénot, V. Monbet [START_REF] Ailliot | Long term object drift forecast in the ocean with tide and wind[END_REF], D. Gérard-Varet [START_REF] Gérard-Varet | Highly rotating fluids in rough domains[END_REF] and D. Bresch, D. Gérard-Varet [START_REF] Bresch | Roughness-Induced Effects on the Quasi-Geostrophic Model[END_REF] and develop an asymptotic expansion of the QG equation. Writing the first order terms in ε, we get an approximate solution of this equation. But for each order, due to the presence of a boundary layer, we split the function in two parts: one deals with the interior of the domain, and the other one is a corrector which takes into account the western part of the domain. The systems obtained in this way are well-posed from a mathematical point of view, and easy to solve numerically. Notice that we only need to compute the first two orders of expansion to get a good approximation. The paper is organized as follows: in Section 2, we start with a derivation of the systems based on a multiple scale approximation of the QG equation. We also give mathematical results for both the one dimensional and two dimensional QG-type equations. Section 3 presents the numerical method and results for the stationary equation in dimension one or two, without and with the ε-periodic topography. In Section 4, we come back to the one dimensional equation with topography. We show that the effects of the topography can be computed separately as a corrector-like function we must add to the solution of the equation with flat bottom. Consequently, the cost of the ε-periodic topography is low.

Theoretical developments.

In this part, we derive the systems based on a multiscale approximation of a Quasi-Geostrophic type equation. First, we present the equation itself and its main characteristics, and then we introduce the multiscale asymptotic developments, starting with a simple example in one dimension. In this case, we also prove some properties of well-posedness before giving the equations satisfied at each order. At this point, we are able to come back to the complete two dimensional stationary Quasi-Geostrophic equation and express the systems which give the first orders of the approximation of the solution.

The Quasi-Geostrophic equation.

The main goal of this paper is the study of the following stationary Quasi-Geostrophic type equation:

-ε 2 ∆ 2 ψ + J ψ,ε 2 ∆ψ + b x ε ,y + ε -1 y = ε -1 f (x,y) in D, ψ = 0 on ∂D, ∆ψ = 0 on ∂D, (2.1)
where the unknown ψ is called the stream function, J is the jacobian defined by J(f,g) = ∂ x f∂ y g -∂ y f∂ x g, b is a periodic function, of null mean value (it represents the topography) and D = [0,1] × [0,1]. The parameter ε denotes a small number. We choose b such that, for x from 0 to 1, b(x/ε) has an entire number of periods. The function f will vanish on the boundaries corresponding to y = 0 and y = 1 and is at the main order (in ε). We recall that the stream function ψ is given by ψ = ∇ ⊥ u, where u is the velocity of the fluid.

Let us give a comment about this equation. When ε tends to zero, the type of the equation is modified: we do not need the boundary condition on the western part of the domain. So we have a boundary layer on the outflow boundary (for x = 0).

2.2.

The one dimensional case: mathematical results and multiscale system.

We start with the study of a restriction of Equation (2.1) to one dimension in space in order to introduce the method on a simple model (we add an evolution in time as it does not bring any real difficulty):

∂ t ψ(t,x) -∂ 2 x ψ(t,x) -ε -1 ∂ x ψ(t,x) = ε -1 f (t,x) in [0,T ] × D, ψ(t,0) = ψ(t,1) = 0 ∀t ∈ [0,T ], ψ(0,x) = 0 ∀x ∈ D, (2.2)
for

T ∈ R ⋆ + , D = ]0,1[ and f ∈ L 2 0,T ;H -1 (D) .
The parameter ε is supposed to be small.

Mathematical properties.

Let us give some mathematical properties of the solution of Equation (2.2). First, thanks to a classical analysis, if the source term f is in L 2 (0,T ;H -1 (D)), the solution ψ is in L ∞ 0,T ;L 2 (D) ∩ L 2 0,T ;H 1 0 (D) . We can also prove the following theorem using Lions theorem for parabolic equations (see [START_REF] Lucas | Effets de petites échelles, du tenseur des contraintes, des conditions au fond et à la surface sur les équations de Saint-Venant[END_REF]): Theorem 2.1. We denote by (.,.) the usual L 2 (D) scalar product and we define the bilinear form a ε by

a ε (u(t),v) = D ∂ x u(t) ∂ x v - 1 ε v dx.
Then, for all ε, there exists a unique solution

ψ ε ∈ C [0,T ];L 2 (D) ∩ L 2 0,T ;H 1 0 (D) such that, for all v ∈ H 1 0 (D) and all f ∈ L 2 0,T ;H -1 (D) , ∂ ∂t (ψ ε (t,x),v(x)) + a ε (ψ ε (t,x),v(x)) = 1 ε < f (t,x),v(x) > H -1 (D)×H 1 0 (D)
with ψ ε (0,x) = 0 for all x ∈ D.

Approximate system based on a multiscale analysis.

In this part, we propose a new approximation method to solve Equation (2.2). We introduce a new space scale and we perform an asymptotic development, for a small ε. This gives an approximate solution of our problem.

Construction of the approximate solution.

We seek an approximate solution ψ app of the problem (2.2) as the sum of an 'interior' term and a 'corrector' term that depends on the quick scale x/ε. We also decompose these two terms in powers of ε in order to write ψ app as:

ψ app (t,x) = ∞ i=0 ε i ψ int i (t,x) + ψ cor i t, x ε .
Let us suppose that the function f does not depend on the quick scale and write the equations satisfied by each ψ int i and ψ cor i . The interior function, defined by ψ interior app

(t,x) = ∞ i=0 ε i ψ int i (t,x)
, is solution of the equation:

∂ t ψ interior app (t,x) -∂ 2 x ψ interior app (t,x) - 1 ε ∂ x ψ interior app (t,x) = 1 ε f (t,x) in [0,T ] × D,
which reads:

∞ i=0 ε i ∂ t ψ int i (t,x) -∂ 2 x ψ int i (t,x) - 1 ε ∂ x ψ int i (t,x) = 1 ε f (t,x) in [0,T ] × D.
We assume that f does not contain any term of order ε j with j > 0 and we identify the powers of ε; we find the equations:

terms in 1 ε -∂ x ψ int 0 (t,x) = f (t,x) in [0,T ]× D, ψ int 0 (t,1) = 0 ∀t ∈ [0,T ],
terms of higher order (i ≥ 1)

∂ x ψ int i (t,x) = ∂ t ψ int i-1 (t,x) -∂ 2 x ψ int i-1 (t,x) in [0,T ]× D, ψ int i (t,1) = 0 ∀t ∈ [0,T ].
We perform the same reasoning for the corrector term defined by the relation ψ corrector app (t,y) = ∞ i=0 ε i ψ cor i (t,y) solution of the homogeneous equation. Here, the boundary conditions are such that the sum of the corrector and the interior terms satisfy Equation (2.2), that is:

   ψ corrector app (t,0) = -ψ interior app (t,0), ∀t ∈ [0,T ], ψ corrector app t,
x ε → 0 when x is far from the boundary layer (x ≫ ε).

Remark 2.1. These boundary conditions express the role of the corrector term. On the one hand, it must correct the interior solution such that their sum satisfies the conditions of the inital equation. On the other hand, far from the boundary layer, its influence should be very small. We define X = x/ε: for a fixed x, when ε → 0, X → +∞. Then we can see the corrector ψ corrector app as a function of X, which gives:

ε 2 ∂ t ψ corrector app (t,X) -∂ 2 X ψ corrector app (t,X) -∂ X ψ corrector app (t,X) = 0 in [0,T ]× R + , ψ corrector app (t,0) = -ψ interior app (t,0) ∀t ∈ [0,T ], lim X→+∞ ψ corrector app (t,X) = 0.
If we use the development of the corrector ψ corrector app , we get:

∞ i=0 ε i ε 2 ∂ t ψ cor i (t,X) -∂ 2 X ψ cor i (t,X) -∂ X ψ cor i (t,X) = 0 in [0,T ]× R + , ∞ i=0 ε i ψ cor i (t,0) = -∞ i=0 ε i ψ int i (t,0) ∀t ∈ [0,T ], lim X→+∞ ∞ i=0 ε i ψ cor i (t,X) = 0.
As before, we identify the powers of ε to write a system for each ψ cor i :

for i = 0 or 1        ∂ 2 X ψ cor i (t,X) + ∂ X ψ cor i (t,X) = 0 in [0,T ]× R + , ψ cor i (t,0) = -ψ int i (t,0) ∀t ∈ [0,T ], lim X→+∞ ψ cor i (t,X) = 0,
and for i ≥ 2

       ∂ 2 X ψ cor i (t,X) + ∂ X ψ cor i (t,X) = ∂ t ψ cor i-2 (t,X) in [0,T ]× R + , ψ cor i (t,0) = -ψ int i (t,0) ∀t ∈ [0,T ], lim X→+∞ ψ cor i (t,X) = 0.
The condition lim X→+∞ ψ cor i (t,X) = 0 will be replaced by ψ cor i (t,M ) = 0, for M large enough.

Existence of the ψ int i and ψ cor i . We start by studying the conditions for the existence of the interior functions. The function ψ int 0 is defined by ∂ x ψ int 0 (t,x) = -f (t,x) with ψ int 0 (t,1) = 0. If we want ψ int 0 in the space L ∞ 0,T ;L 2 (D) ∩ L 2 0,T ;H 1 (D) , we have to take the source term f in L ∞ 0,T ;L 1 (D) ∩ L 2 0,T ;L 2 (D) and then ψ int 0 is in L ∞ 0,T ;C( D) . If we also seek ψ int 0 as a continuous function in time, we must choose f in C [0,T ];L 1 (D) and consequently ψ int 0 is in C([0,T ]; D). For the following orders, we easily prove by recurrence (see [START_REF] Lucas | Effets de petites échelles, du tenseur des contraintes, des conditions au fond et à la surface sur les équations de Saint-Venant[END_REF]) that it is necessary that f belongs to:

W i,∞ 0,T ;L 1 (D) ∩ ••• ∩ W 1,∞ 0,T ;W i-1,1 (D) ∩ L ∞ 0,T ;W i,2 (D)
to define ψ int i and if we want the continuity in time, we impose f to be in:

C i [0,T ];L 1 (D) ∩ ••• ∩ C 1 [0,T ];W i-1,1 (D) ∩ C [0,T ];W i,2 (D) .
Looking carefully at the formulae for the corrector, we see that the existence of the ψ int i directly gives the existence of the ψ cor i . Remark 2.2. A priori we do not know that ψ app (0,x) = 0, for all x in D, that is:

∞ i=0 ε i ψ int i (0,x) + ψ cor i 0, x ε = 0, ∀x ∈ D.
So we impose f (0,x) = 0, for all x in D that ensures that the boundary conditions of (2.2) are satisfied when ε tends to zero.

Convergence when ε tends to zero.

We choose an integer N ≥ 0 and, for f regular enough, we define the approximation at the order N by:

ψN app (t,x) = N i=0 ε i ψ int i (t,x) + ψ cor i t, x ε .
With classical methods (see [START_REF] Lucas | Effets de petites échelles, du tenseur des contraintes, des conditions au fond et à la surface sur les équations de Saint-Venant[END_REF]), we prove that, when ε tends to zero, the approximate solution at the order N ψN app converges to the solution ψ of Equation (2.2).

Systems in two dimensions.

We now apply the previous method to the two dimensional case (2.1). The quick scale X is still defined by X = x/ε.

If we perform a classical asymptotic development in ψ, at the order 1/ε 2 , we get the following equation:

     -∂ 4 X ψ 0 (x,X,y) + ∂ X ψ 0 (x,X,y) = 0 in D, ∂ 2 X ψ 0 (x,X,y) = 0 for x = 0 and x = 1, ψ 0 (x,X,y) = 0 for x = 0 and x = 1.
Then we see that, as in one dimension, we can split ψ 0 into an 'interior' term ψ int 0 (x,X,y) and a 'corrector' term ψ cor 0 (X,y). The 'interior' term ψ int 0 (x,X,y), that depends on the X variable due to the X dependence of the topography function b, must vanish for x = 1 and is periodic in X. The 'corrector' term ψ cor 0 (X,y) has an influence only on the boundary layer (it vanishes far from the western coast) and must be equal to the opposite of ψ int 0 on the western boundary of the domain. Now we are able to write the two systems for the two terms at the first order:

   -∂ 4 X ψ int 0 (x,X,y) + ∂ X ψ int 0 (x,X,y) = 0 in D, ∂ 2 X ψ int 0 (0,0,y) = 0 and ∂ 2 X ψ int 0 (1, 1 ε ,y) = 0 for y ∈ [0,1], ψ int 0 (1,y) = 0 for y ∈ [0,1],
which means that ψ int 0 does not depend on X (due to the periodicity). Moreover, for the corrector, we have:

# " !      -∂ 4 X ψ cor 0 (X,y) + ∂ X ψ cor 0 (X,y) = 0 in D, ∂ 2 X ψ cor 0 (X,y) = 0 for X = 0 and X = 1 ε , ψ cor 0 (0,y) = -ψ int 0 (0,y), for y ∈ [0,1] (2.3)
and can be computed as soon as we know the interior term. Note that the dependence in the y variable is only related to ψ int 0 (0,y) and can be added afterwards. Up to now, we have written the first order term of Equation (2.1). We have expressed the corrector ψ cor 0 as a function of ψ int 0 and we know that ψ int 0 does not depend on X. But we do not have any equation for ψ int 0 . So we are led to study the next order. Equation (2.1) at the order 1/ε gives:

         -∂ 4 X ψ 1 -4∂ 3 X ∂ x ψ 0 + ∂ X ψ 0 ∂ y ∂ 2 X ψ 0 + ∂ x ψ 0 +∂ X ψ 1 -∂ y ψ 0 ∂ 3 X ψ 0 -∂ y ψ 0 ∂ X b = f in D, ∂ 2 X ψ 1 = 0 for x = 0 and x = 1, ψ 1 = 0 for x = 0 and x = 1.
Like before, we split the first order ψ 1 into ψ int 1 , such that ψ int 1 is periodic in X, and ψ cor 1 , equal to zero almost everywhere except near the boundary {0} × [0,1]. We have the equality:

-∂ 4 X ψ int 1 + ∂ x ψ int 0 + ∂ X ψ int 1 -∂ y ψ int 0 ∂ X b = f.
We compute the mean value in X, using the periodicity and the fact that ψ int 0 does not depend on X, and we get ∂ x ψ int 0 (x,y) = f (x,y), so we have an expression of ψ int 0 on D:

ψ int 0 (x,y) = x 1 f (s,y)ds.
Moreover, we have information on ψ int 1 : the X variation is given by: §

¦ ¤ ¥ -∂ 4 X ψ int 1 + ∂ X ψ int 1 = ∂ y ψ int 0 ∂ X b in D, with ∂ 2 X ψ int 1
equal to zero on the left and right boundaries of the domain D, and ψ int

1 vanishes for X = 1/ε.
Then, if we want to know ψ int 1 , we are led to identify the dependence of ψ int 1 in x, that is the function D(x,y) such that:

ψ int 1 (x,X,y) = X 0 ∂ e X ψ int 1 (x, X,y)d X + D(x,y),
and with D(1,y) given so that ψ int 1 vanishes on the right boundary of the domain. We will compute D thanks to Equation (2.1) at the order 1.

Before going to the next order, we can also express ψ cor 1 ; it must vanish far from the boundary layer and satisfy:

1 0 ( )          -∂ 4 X ψ cor 1 + ∂ X ψ cor 1 = -∂ X ψ cor 0 ∂ y ∂ 2 X ψ cor 0 +∂ y ψ cor 0 ∂ 3 X ψ cor 0 + ∂ y ψ int 0 ∂ 3 X ψ cor 0 + ∂ y ψ cor 0 ∂ X b in D, ∂ 2 X ψ cor 1 = 0 for X = 0 and X = 1 ε , ψ cor 1 (0,0,y) = -ψ int 1 (0,0,y) for y ∈ [0,1].
We just need the first order and ψ int 1 to get ψ cor 1 . So we have to express the function D to have the complete first order solution.

To do so, we write Equation (2.1) at the order 1:

               -∂ 4 X ψ 2 -4∂ 3 X ∂ x ψ 1 -6∂ 2 X ∂ 2 x ψ 0 -2∂ 2 X ∂ y ψ 0 + ∂ x ψ 0 ∂ y ∂ 2 X ψ 0 + ∂ X ψ 1 ∂ y ∂ 2 X ψ 0 +2∂ X ψ 0 ∂ y ∂ x ∂ X ψ 0 + ∂ X ψ 0 ∂ y ∂ 2 X ψ 1 + ∂ x ψ 0 ∂ y b + ∂ X ψ 1 ∂ y b + ∂ x ψ 1 +∂ X ψ 2 -3∂ y ψ 0 ∂ x ∂ 2 X ψ 0 -∂ y ψ 1 ∂ 3 X ψ 0 -∂ y ψ 0 ∂ 3 X ψ 1 -∂ y ψ 1 ∂ X b = 0 in D, ∂ 2 X ψ 2 = -∂ 2 x ψ 0 -∂ 2 y ψ 0 for x = 0 and x = 1, ψ 2 = 0 for x = 0 and x = 1.
Again, we split into interior and corrector functions and we study the former:

-∂ 4 X ψ int 2 -4∂ 3 X ∂ x ψ int 1 + ∂ x ψ int 0 ∂ y b + ∂ X ψ int 1 ∂ y b + ∂ x ψ int 1 + ∂ X ψ int 2 -∂ y ψ int 0 ∂ 3 X ψ int 1 -∂ y ψ int 1 ∂ X b = 0,
and we compute the mean value in X, with ψ int 2 periodic as for the previous orders. Then we have:

-4∂ 3 X ∂ x ψ int 1 X + ∂ x ψ int 0 ∂ y b X + ∂ X ψ int 1 ∂ y b X + ∂ x ψ int 1 X -∂ y ψ int 0 ∂ 3 X ψ int 1 X -∂ y ψ int 1 ∂ X b X = 0,
which gives ∂ x D(x,y). Since D(1,y) = 0 for all y in [0,1], a simple integration from x = 1 gives D (we give its explicit formulation in the following, in some particular cases). Remark 2.3. Note that it does not ensure that ψ vanishes on the boundaries [0,1] × {0} and [0,1] × {1}. In the following, we prove that, if we consider a topography that vanishes on these boundaries, this property is satisfied.

In the same way, Equation (2.1) at the order 1 gives the variation of ψ int 2 in X. The dependence in x will be expressed thanks to the average in X of Equation (2.1) at the order ε. For the corrector, we just need the previous terms.

Then, writing Equation (2.1) at each order, we can obtain each term of the asymptotic development.

2.3.1.

A particular case: b does not depend on y. We suppose that b does not depend on y and we study the first order interior equations. First, we define a function G(X) solution of:

     -G (4) (X) + G ′ (X) = b ′ (X) in [0,1/ε], G ′′ (X) = 0
for X = 0 and X = 1/ε, G(X) = 0 for X = 0 and X = 1/ε.

Then ψ int 1 reads ψ int 1 (x,X,y) = ∂ y ψ int 0 (x,y)G(X) + D(x,y) with the following equation for D:

∂ x D(x,y) + ∂ x ∂ y ψ int 0 (x,y)G(X) X -∂ 2 y ψ int 0 (x,y)G(X)∂ X b X = 0 in D, D(1,y) = 0 for y ∈ [0,1].
This system can be easily solved numerically.

2.3.2. Another particular case: b with separable variables.

We choose b as b(X,y) = b 1 (X)b 2 (y), with b 1 periodic and with a null mean value. As before, we define a function G solution of:

     -G (4) (X) + G ′ (X) = b ′ 1 (X) in [0,1/ε],
G ′′ (X) = 0 for X = 0 and X = 1/ε, G(X) = 0 for X = 0 and X = 1/ε, and

ψ int 1 is given by ψ int 1 (x,X,y) = ∂ y ψ int 0 (x,y)b 2 (y) G(X) + D(x,y).
The equation to solve to get D is much more intricate than in the previous case:

       ∂ y ψ int 0 (x,y)b 2 (y)b ′ 2 (y)b 1 (X) G ′ (X) X + ∂ x D(x,y) + ∂ x ∂ y ψ int 0 (x,y)b 2 (y) G(X) X -b 2 (y) ∂ 2 y ψ int 0 (x,y)b 2 (y) + ∂ y ψ int 0 (x,y)b ′ 2 (y) G(X)b ′ 1 (X) X = 0 in D, D(1,y) = 0 for y ∈ [0,1].
These relations can also be solved without any real difficulty.

Numerical results.

In this part, we give numerical results obtained thanks to the previous derivations. We study the first and the second orders of both the interior and corrector functions, in the one and two dimensional cases. We also explain how we changed the mathematical equalities into relations that can be implemented easily.

One dimensional case.

3.1.1. Implementation. The system have been implemented using Fortran. We define two different space steps (in the boundary layer and outside the boundary layer). We choose centered finite differences to compute the ψ cor i and the ψ int i , and for the former we use the LU decomposition with derivatives at the second order. We also approximate the solution of Equation (2.2) by centered finite differences and LU decomposition but without asymptotic development.

Stationary results.

We consider the case f (t,x) = x in order to be able to express the exact solution. One can also note that, in this case, the terms in the asymptotic development ψ i vanish from i = 2. We find that the approximated solution is very close to the exact solution and, contrarily to the "usual" solution, we do not have a strong stability condition and we are less restricted to choose the parameters.

Improvements for the evolution equation.

Then, we introduce the evolution in time (we just have to add some time derivatives). We need to reduce the computation time of the corrector: the condition lim X→+∞ ψ cor i (t,X) = 0 is replaced by ψ cor i (t,M ) = 0 for M large enough. The question is: how to choose M ? The first idea is to take M such that εM is outside of the domain. If we look at the shape of the corrector, we see that its contribution is significant only on a small part of the domain, and more precisely only on the boundary layer. Unfortunately, the computation time is not significantly reduced. Consequently, another idea is to compute the corrector only in the boundary layer. This method will be validated and used in what follows. In order to prove the accuracy of the reduction of the support length, we compare the exact corrector and the corrector that vanishes after a bound, chosen as one and a half times the boundary layer size. The results, presented in Figure 3.1, prove that, for both the stationary and the evolution problems, the truncated corrector is close to the exact one, even for long times, but regarding the computation time, our approximation is much better.

Then, we have an efficient method to compute the solution of Equation (2.2), based on asymptotic developments and a truncation of the boundary layer corrector term. In the following, we consider the two dimensional case and give some numerical results.

Numerical results

for the 2D case.

Programming issues.

The first difficulty due to this new equation is the complexity of the equation for the corrector, which is not of order two any more but of order 4. As noted above, we will add the y variation afterwards (we substitute ψ cor 0 (0,y) with 1). We use a centered scheme to express the fourth derivative. To get the second derivative on the boundaries, we introduce a fictive point on the exterior of the domain, which can be combined with the relation for the fourth derivative. The next problem occurs the computation of ψ cor 1 . We have an equation in the X variable whereas ψ int 0 depends on x and y and we recall that X = x/ε. We consider ψ int 1 as a piecewise constant function in X, i.e. for εX ∈ [x i ,x i+1 ], ψ int 1 (εX,y) = ψ int 1 (x i ,y). 3.2.2. Numerical results. We take f (x,y) = -sin(2πy) on a 500 × 500 grid. At the first order, the topography does not play any role, and, as for the one dimensional case, the corrector can be considered as equal to zero beyond a given bound (we chose here x = 0.2, see Figures 3.2 to 3.4).

In Figures 3.6 to 3.8, we chose the bottom profile function b(X,y) = 1/2cos(πX/10).

Here again, we can remark that the bottom has a non-negligible influence over the whole domain but also in the boundary layer. " . We have also tested the topography given by b(X,y) = 1/2cos(πX/10)sin(πy) and the results are plotted in Figures 3.10 to 3.12. As the topography vanishes on the boundaries y = 0 and y = 1, the boundary condition is satisfied by the sum ψ 0 + εψ 1 . Now we have a program that can compute quickly the solution of Equation (2.1), even for small ε (such as 0.001). However, we cannot assert that the boundary conditions are satisfied for y = 0 or 1. In the following section, we compare our results to the classical Quasi-Geostrophic model.

+ εψ 1 = ψ int 0 + ψ cor 0 + εψ int 1 + εψ cor 1 for ε = 0.01.

Comparison.

In order to compare our model to a classical Quasi-Geostrophic model, we consider the equation:

∂ t ∆ψ + σ∆ψ -∆ 2 ψ + J ψ,∆ψ + b x ε ,y + ε -3 y = f(x,y). (3.1)
Note that multiplying Equation (3.1) by ε 2 and assuming σ = 0, b x ε ,y = ε -2 b x ε ,y and f(x,y) = ε -3 f (x,y), the stationary regime satisfies Equation (2.1).

In Figure 3.13, we present the differences between the two programs, discretized with 200 points in both cases. We study separately the boundary layer. Outwards, the difference is very small, of order of 0.12%, see Figure 3.13, right. As ε is equal to 10 -2 , we have a good approximation. However, in the boundary layer, the approximation is not so accurate: let us study this phenomenon. The width of the boundary layer, following [START_REF] Kazantsev | Compact finite difference schemes for shallow-water ocean model[END_REF], is 2πε/ √ 3 ≈ 0.036 which means that the classical program have only 7 points in this region, which is not well solved. Our approximation gives a better result thanks to the new variable X: the computation uses one hundred times more points in the boundary layer.

We can also prove that our approximation is at the second order: for large enough ε (of order of 0.1 to have a good classical solution) we plot the infinite norm of the difference between the two results. In Figure 3.14, which represents the error as a function of ε in a log-log scale, we get a line with a slope equal to 2, so the error is in ε 2 . 4. Back to the one dimensional case: study of the topography. We can add the topography in Equation (2.2) and consider the following problem:

∂ t ψ(t,x) -∂ 2 x ψ(t,x) + b x ε ∂ x ψ(t,x) - 1 ε ∂ x ψ(t,x) = 1 ε f (t,x) in [0,T ] × D, ψ(t,0) = ψ(t,1) = 0 ∀t ∈ [0,T ], ψ(0,x) = 0 ∀x ∈ D, (4.1)
with

D = ]0,1[, T ∈ R ⋆ + , f ∈ L 2 0,T ;H -1 (D) .
The function b is periodic, its period is P , its mean value is zero, and, for x from 0 to 1, b(x/ε) has an entire number of periods. Note that we could assert the same a priori estimates and existence results as in Section 2.2. We are interested in the role of the topography in this equation: we write the asymptotic development as before and we isolate the bottom effect.

Asymptotic development.

As for a flat bottom, we split the unknown function into interior terms, that are periodic in x/ε, and corrector terms. The time-dependence does not really change the problem, so, for the sake of simplicity, we consider the stationary equation. Thanks to the periodicity, we get a system for the interior function. We will write it at each order in ε:

   -∂ 2 x ψ interior x, x ε + b x ε ∂ x ψ interior x, x ε - 1 ε ∂ x ψ interior x, x ε = 1 ε f (x), ψ interior (x = 1) = 0. (4.
2) The non-periodic part gives the system for the boundary layer corrector:

         -∂ 2 x ψ corrector x ε + b x ε ∂ x ψ corrector x ε - 1 ε ∂ x ψ corrector x ε = 0, ψ corrector (X = 0) = -ψ interior (x = 0), ψ corrector (X) -→ X→+∞ 0. (4.3)
Let us study separately these two systems at the first orders in ε.

The interior function.

The function ψ interior is periodic in x/ε and should satisfy Equation (4.2). We still denote by X the quick scale x/ε, and we perform an asymptotic development of our function, writing ψ interior as the sum of the ψ int i .

At the order 1/ε 2 , we have -∂ 2 X ψ int 0 (x,X) -∂ X ψ int 0 (x,X) = 0. The function ψ 0 is assumed to be periodic in the quick scale, so it does not depend on X and we have to study the next order to have more information.

At the order 1/ε, Equation (4.2) reads:

-∂ 2 X ψ int 1 (x,X) -∂ x ψ int 0 (x) -∂ X ψ int 1 (x,X) = f (x). (4.4)
We integrate in X and, with the periodicity, we find:

ψ int 0 (x) = 1 x f (s) ds.
Moreover, when we replace this value in Equation (4.4), we obtain (exactly as for ψ int 0 ) that ψ int 1 does not depend on X. Since the functions ψ int 0 and ψ int 1 do not depend on X, Equation (4.2) at the order 1 gives:

-∂ 2 x ψ int 0 (x) -∂ 2 X ψ int 2 (x,X) + b(X)∂ x ψ int 0 (x) -∂ x ψ int 1 (x) -∂ X ψ int 2 (x,X) = 0.
With the mean value in X, we find ψ int 1 (x) = f (x)f (1) and we have the X derivative of ψ int 2 : At the order ε, Equation (4.2) becomes:

∂ X ψ int 2 (x,X) = ∂ x ψ int 0 (x)
-∂ 2 x ψ int 1 (x) -2∂ x ∂ X ψ int 2 (x,X) -∂ 2 X ψ int 3 (x,X) + b(X)∂ x ψ int 1 (x) + b(X)∂ X ψ int 2 (x,X) -∂ x ψ int 2 (x,X) -∂ X ψ int 3 (x,X) = 0.
Its mean value can be written as

-∂ 2 x ψ int 1 (x) + b(X)∂ X ψ int 2 (x,X) X -∂ x ψ int 2 (x,X) X = 0,
where g X is the mean value of the function g in X.

Replacing the expression of ∂ X ψ int 2 in this equality, we obtain: As we have previously mentioned, at the order 1 the two terms differ in the boundary layer only. If we denote by φ 1 this difference, that is φ 1 (t,X) = ψ cor 1 (t,X)ψ cor 1 (t,X), it must satisfy the following equation:

-∂ x D D D(x) = N ∂ 2 x ψ int 0 (x) -M ∂ x ψ int 0 (x) + ∂ 2 x ψ int 1 (x),
∂ 2 X φ 1 (t,X) + ∂ X φ 1 (t,X) = b(X)∂ X ψ cor 0 (t,X) = b(X) 1 0 f (t,s) ds exp(-X), φ 1 (t,X = 0) = 0, φ 1 (t,X) -→ X→+∞ 0.
So if we want an approximation of the solution with topography ψ at the order 1, we just have to add the function φ 1 to the solution ψ for a flat bottom.

4.4.2. Topography effect at the order 2.

At the order 2, we have to consider the corrector terms but also the interior terms, as they both differ.

In the boundary layer, we define φ 2CL (t,X) = ψ cor 2 (t,X)ψ cor 2 (t,X) the solution of the system

∂ 2 X φ 2CL (t,X) + ∂ X φ 2CL (t,X) = b(X)∂ X ψ cor 1 (t,X), φ 2CL (t,X = 0) = ψ cor 2 (t,X = 0) -ψ cor 2 (t,X = 0), φ 2CL (t,X) -→ X→+∞ 0.
Since we are able to replace ∂ X ψ cor 1 (t,X) by its value, namely: we have an expression of the boundary layer term that controls the topography effect at the order 2.

In the interior of the domain, the difference between the two solutions, φ 2int , given by φ 2int (t,x,X) = ψ int 2 (t,x,X)ψ int 2 (t,x,X), is periodic in X. We define another function on [0,+∞[, denoted by G G G, periodic, with period P , solution of the following system:

∂ 2 X G G G(X) + ∂ X G G G(X) = b(X), G G G(0) = G G G(P ) = 0.
Then we can write φ 2int as:

φ 2int (t,x,X) = ∂ x ψ int 0 (t,x)G G G(X) + E(t,x) = -f (t,x)G G G(X) + E(t,x), with ∂ x E(t,x) = G G G(X) X ∂ x f (t,x) -b(X)∂ X G G G(X) X f (t,x), E(t,1) = f (t,1)G G G 1 ε .
Thus we can compute the term due to the topography that should be added to modify the solution in the interior of the domain at the order 2.

In the same way, we could express the topography terms for the upper orders but the results proposed here already give the first three orders of the approximation.

In order to prove the accuracy of this method, we performed the computations and we give the mains results in the following part.

Numerical experiments on the topography terms.

In this section, we compare the solution of the equation to the topography term to the solution of the equation for a flat bottom to which we add the topography terms we found in the previous section.

Computing details.

We consider the evolution equation but we choose f as a function of x only. We compute the solution with a flat bottom ψ and the solution with a varying bottom ψ using finite differences at the second order. As the stability condition (that reads dt ≤ 2dx 2 ε 2 /(4ε 2 + dx 2 )) is strong, the computations are time expensive. We stop when the solution is stationary. Then we compare the function ψ and the sum of ψ and the topography terms φ k (we only compute these terms for the final time). For the topography terms, we use a LU decomposition to get the function G G G, with derivatives at the second order. We take a constant space step as G G G is periodic (the values of G G G have been compared with its theoretical expression for various bottoms). We also need the expression of ψ cor 1 : here we use the handwritten expression but in the following we will get it thanks to ψ cor 1 and φ 1 . At last, we get the values of E with finite differences at the first order.

Numerical results.

In the following figures, we present the results for the bottom given by 1/2cos(πX/10) + 1/2cos(πX/5) (we studied two other choices and the conclusions are the same). We have plotted, for f (t,x) = x, the differences between:

• ψ and ψ,

• ψ + εφ 1 and ψ,

• and between ψ + εφ 1 + ε 2 φ 2CL + ε 2 φ 2int = ψ + εφ 1 + ε 2 φ 2 and ψ. We remark that the difference between the solution of the equation with the topography and the solution for a flat bottom to which we added the topography terms is of the predicted order in ε.

Contribution of the topography

Remark 4.3. The curve of the difference ψψ represents the effect of the topography. We have a logical oscillation similar to the bottom. However we can note a large influence in the boundary layer.

4.6. Complete method. The ultimate step is to replace the computation of the function ψ with finite differences by its approximation given thanks to the first orders of the asymptotic development. The results are plotted on Figure 4.2 (the bottom is still equal to 1/2cos(πX/10) + 1/2cos(πX/5)). We see that we are able to get an approximation of the oscillating solution at the order ε or ε 2 . Then we have an efficient method to compute the solution for a flat bottom. We can modify this function to get the solution of Equation (4.1) for any bottom. This program is rapid and can give an approximation of the solution at a determined order.

Approximations of the solution

Conclusion

In this paper, we studied the Quasi-Geostrophic equation from two different points of view: first, with an oscillating topography, we introduced a new variable and a corrector term to take into account the quick scale and the boundary layer respectively. The first orders of the asymptotic development give us systems that have been implemented to get approximations of the solution. This method is efficient as it gives results comparable to those obtained with usual technics but it is less expensive in time and computations. In the second part, we studied the role of the topography: we manage to find the terms that must be added to the solution with a flat bottom to take into account the topography. So we can solve the Quasi-Geostrophic equation without the topography term, and then add the contribution of the topography.

Fig. 3 . 1 .

 31 Fig. 3.1. Difference between the correctors for f (t,x) = x (on the left) and f (t,x) = x + sin(t) (on the right), for ε = 0.01 and a time step equal to 0.1. We define the boundary layer length as 0.1, and consider 31 points in the boundary layer, 51 points outside.

Fig. 3 . 2 .Fig. 3 . 3 .Fig. 3 . 4 .

 323334 Fig. 3.2. Plots of ψ int 0 and ψ cor 0 for ε = 0.01.

Fig. 3 . 5 .

 35 Fig. 3.5. Topography: b(X) = 1 2 cos " πX 10

Fig. 3 . 6 .of ψ int 1 and ψ cor 1 forFig. 3 . 7 .

 361137 Fig. 3.6. Plots of ψ int 1 and ψ cor 1 for ε = 0.01.

Fig. 3 . 8 .+ εψ cor 1 forFig. 3 . 9 .

 38139 Fig. 3.8. Level curves of ψ 0 + εψ 1 = ψ int 0 + ψ cor 0 + εψ int 1 + εψ cor 1 for ε = 0.01.

Fig. 3 . 10 . 1 and ψ cor 1 forFig. 3 .

 310113 Fig. 3.10. Plots of ψ int 1 and ψ cor 1 for ε = 0.01.

11Fig. 3 .

 3 Fig. 3.13. Difference between the two models: on the left, on the boundary layer and on the right, outside the boundary layer.

Fig. 3 .

 3 Fig. 3.14. Difference between the two programs in a log-log scale.

  )e s ds e -X .Then we know ψ int 2 up to a function D D D(x) that can be determined with the following order.

where M and N depend only on b.ψ int 2 is

 2 )e s ds e -u du X Then completely determined by the expression:ψ int 2 (x,X) = ∂ x ψ int 0 (x)X 0 u 0 b(s)e s ds + 1 e P -1 P 0 b(s)e s ds e -u du +D D D(x).

4. 4 . 1 .

 41 Topography effect on the boundary layer at the order 1.

  ds dy + e -X (f (t,1)f (t,0)),

Fig. 4 . 1 .

 41 Fig. 4.1. For ε = 0.01, dx = 10 -4 , dt = 10 -9 and b(X) = 1/2cos(πX/10) + 1/2cos(πX/5): difference between e ψ and the function ψ to which we added the topography terms.

Fig. 4 . 2 .

 42 Fig. 4.2. For ε = 0.01, dx = 10 -4 , dt = 10 -4 and b(X) = 1/2cos(πX/10) + 1/2cos(πX/5): difference between the function e ψ and the asymptotic approximation of ψ to which we added the topography terms.

Acknowledgement. The authors would like to thank D. Bresch for fruitful discussions and commentaries around this work.

We could carry on this process to get ψ int 3 and more generally each ψ int i but we do not give the details here. Remark 4.1. If we add the time dependence, ψ int 1 is solution of:

and the equation for D D D becomes

We have proved that we could give the expression of all the interior terms; let us consider now the boundary layer corrector terms.

Corrector terms for the boundary layer.

We develop the function ψ corrector as the sum of the ψ cor i and we isolate each order of Equation (4.3).

At the order 1, the equation and the boundary conditions give the expression of

Equation ( 4.3) at the order ε reads:

Consequently, ψ cor 1 is completely determined by the value of ψ int 1 for x = 0. The next orders lead us to the same conclusion: if we know all the (ψ cor i ) i≤j and (ψ int i ) i≤j+1 , then we can compute ψ cor j+1 . Remark 4.2. If we consider the evolution equation, from the third order, we must add the term ∂ t ψ cor i-2 in the right hand side of the equation for ψ cor i . We are now able to compute the first terms of the asymptotic development of the stationary solution of Equation (4.1) with the topography. In the following, we identify the topography effect through these equations and we show how to add this effect to the solution for a flat bottom.

Topography effect.

In this part, we denote by ψ the solution for a flat bottom, as in section 2.2.2, and by a tilde the solution of Equation (4.1) that, as in the previous sections, depends on an oscillating bottom.

Our goal is to express the effect of the topography, that is the term we could add to ψ to get ψ. To do so, we compare the first orders (we stop at the third one) of the solutions of the equations with and without topography, and we find the equation of a so-called 'topography effect'. We remark that the two developments are equal at the first order, and at the second order for the interior. So we focus on the differences between ψ int i and ψ int i for i ≥ 2 and between ψ cor i and ψ cor i for i ≥ 1.


Jacqueline Vaissière

Prediction of velum movement from phonological specifications

) that the strategies used by the two speakers are different with respect to control of movement velocity.

We propose a suprasegmental feature [+ strong] to account for observed movement patterns including timing relations between velum raising and lowering gestures and articulatory gestures of consonantal articulators such as the lips and the tongue tip.

Material and Method

.

----------------------

Insert Figure 2 here ---------------------------------------------Figure 2 represents the velum height time function in the same sentence and two other sentences uttered by the same speaker (i.e., speaker A). The speech signal envelope is shown below each velum time function. The speech signal could be played and listened to for any selected time intervals by an interactive analysis program, in order to verify the phonetic quality of the part or the material in question. Our goal is to determine the control variables that explain the behavior of the velum in such sentences. In particular, we attempt to predict velum height, velum movement velocity, and timing of velum movements relative to the speech waveform and to movements of other articulators.

In order to approach a solution to this complex problem, we followed a procedure of successive approximation, starting from careful but subjective interpretation of data, followed by hypothesizing a set of quantitative rules for predicting the time function of velum position for a given phonological representation of each utterance, then followed by revisions of the rule system and repeating the cycle. The rules were not implemented as computer programs, but were manually implemented to generate time functions (strictly according to the hypothesized rule system). Automatic interpretation of large varieties of utterance materials awaits future studies for testing the validity and completeness of our present tentative model.

Velum pellet height data were originally sampled at 8.5-ms intervals. They were subjected to minimal digital smoothing and plotted as time functions using the same position (ordinate) and time (abscissa) scales for all data studied, maintaining the same origin of the ordinate for the same speaker for comparison of different utterance samples. Time function curves from the same speaker (and the same session) were superimposed and manually copied on semi-transparent paper, with the time origin for each sample being adjusted appropriately for the purpose of comparison wherever applicable. Tracings of different data channels, such as velum height and acoustic signals, were obtained by maintaining the common time scale and origin. All figures illustrated in this paper were obtained in this manner, even though the magnification varies from figure to figure. Wherever appropriate, data from other subjects were consulted in interpreting results, but quantitative analyses pertain to the two speakers only.

Segmental and suprasegmental factors affecting velum height

Nasality Feature First, we should note that velum lowering as the phonetic manifestation of a nasal consonant does not necessarily result in a velum height valley in the time course of velum movement, as expected from a segment concatenation and smoothing point of view.

--------------------------------------------- Insert Figure 3 here ---------------------------------------------
Often, for example, the velum height minimum is found in the preceding acoustic vowel segment (i.e., before the articulatory closure for the nasal murmur) in a syllable containing a final nasal consonant. As expected, the velum was found in general to have a high position for consonants that are phonologically classified as oral ([-nasal]), and a low position for those consonants classified as nasal ([+ nasal]). Figure 3 illustrates the superimposed velum height time functions in utterances containing ([-nasal] consonants (fig. 3a), and in utterances containing [+ nasal] consonants (fig. 3b) for speaker B. As seen in figure 3, the stationary value of velum height for the [-nasal] consonants (fig. 3a) generally represents a peak of velum movement and is higher than that for the [+ nasal] consonants (fig. 3b). For a nasal consonant in the environment of [-nasal] segments (excluding vowels, see below), the velum time function has a local valley, which falls within the time interval of the acoustic nasal murmur. However, the [+ nasal] feature for a consonant is not sufficient to determine velum height; as can be seen from figure 3, both the maximum velum height achieved for [-nasal) consonants and the minimum velum height achieved for [+ nasal] consonants vary a great dear depending on the context. ---------------------------------------------Insert Figure 4 here

Nasality of Following Segmental Context

It has been discussed in the literature that the velum exhibits anticipatory lowering when the following context contains a nasal consonant [START_REF] Henke | Dynamic articulatory model of speech production using computer simulation[END_REF][START_REF] Moll | Investigation of the timing of velar movements during speech[END_REF][START_REF] Kent | Velar movement and timing: Evaluation of a model for binary control[END_REF][START_REF] Kent | Coarticulation in recent speech production models[END_REF]. Figure 4 illustrates with speaker B the influence of the following segmental context on a [-nasal] consonant, for consonant-vowel-consonant sequences CVC and CVN (in this paper we denote a nasal consonant by N and a non-nasal consonant by C). In figures 4a and b, the first [-nasal] consonant is followed by a vowel, and then by a [-nasal] consonant (a), or by a [+ nasal] consonant (b). In both Figures 4a andb, the first consonant C is not directly preceded by a nasal consonant.

There are three differences in velum movement pattern between CVC and CVN. First, in CVC sequences, the maximum velum height for the first C is usually higher than for the first C in the CVN sequences. Second, the maximum velum height for the first C in CVC sequences varies a great deal, in comparison with CVC (compare height values at vertical lines). Third, related to the first point, the velum exhibits anticipatory lowering for the [+ nasal] consonant in some of the CVN utterances, but not in others. The anticipatory lowering, when it occurs, seems to begin during or prior to the first consonant in the CVN sequence.

By careful examination of individual time functions of velum movement, we could not find any evidence for an inherent velum position for the vowels (V). In the non-nasal context (fig. 4a), the velum is high for the vowel ; in the nasal context (fig. 4b), the velum is lower for the vowel. But, as noted before, the actual valley often occurs within the acoustic vowel segment (according to the traditional definition) when it is followed by a tautosyllabic nasal consonant. In this study, we proceeded with the assumption that English vowels, unlike consonants, are unspecified for velum height (see below for further discussion).

Configurational Environment

The time functions for CVN sequences in figure 4b were further separated into two groups and are replotted in figures 4c and d. Figure 4c represents utterances satisfying one of three configurational conditions:

(1) C immediately precedes a stressed vowel.

(2) C is in word -initial position, or (3) C immediately precedes a major syntactic boundary. In figure 4d,C does not occur in any of the three configurational environments. To describe this grouping of data, we devised a binary feature [± strong], to be associated with each [-nasal] consonant. The C in any of the three configurational conditions (fig. 4c) was labeled [+ strong]; otherwise C was labeled [-strong]. An oral consonant with the feature [+ strong] is less likely to be influenced by anticipatory velum lowering, in comparison to an oral consonant with the feature [-strong].

In the set of utterance analyzed here, the consonant preceding a stressed vowel are syllable-initial consonants. Since (phonologically) word-initial consonants are also syllable-initial consonants, it is tempting to associate the behavior of [+ strong] with configurationally initial consonants, and [-strong] with non-initial consonant, perhaps referring to the syllable as the configurational domain [START_REF] Fujimura | Syllable-final nasal elements in English-an x-ray study of velum height[END_REF]. If so, we could dispense with the feature [+ strong]. However, consonants in condition 3 (before a major boundary) are clearly not syllable-initial. It may also be noted that the consonants in category (3) were typically syllable affixes [START_REF] Fujimura | An analysis of English syllables as cores and affixes[END_REF]Clements, in press). More data are needed to give ground to such an analysis based on syntagmatic structure.

Other Segmental Inherent Features

The feature [± strong] is not the only factor that determines anticipatory velum lowering, because the maximum velum height for the [-strong] oral consonants in figure 4d still exhibits a great deal of variation. The remaining differences in the amount of anticipatory velum lowering among the [-strong] samples and also among the [+ strong] samples where there is variation seem to be largely explainable by taking into account the inherent phonetic features of the consonant in question. First, velum height for the non-obstruent /r/ and /1/ or clusters containing /r/ and /1/ (dotted-line time functions in both figure 4c andd) is relatively more affected by the following N, for both [+ strong] and [strong] consonants than is the velum height for other consonants. Second, velum height for the velar consonants /k/ and /g/, as shown by the crossed-line time functions in figure 4d, is less affected by a following N, in comparison with other consonants. In the case of the velar stops, the high velum position is probably due to direct physical interference of the tongue surface pushing up the velum. This is also why the velar nasal consonant /ŋ/ has a higher velum position than other nasal consonants (fig. 15). ---------------------------------------------Insert Figure 5 here ---------------------------------------------Figure 5 illustrates velum movement in CVN utterances, contrasting stressed versus unstressed environments for both speakers. As can be seen in figure 5, for both speaker B and speaker A, the velum excursion from oral consonant C to nasal consonant N is generally more extensive in stressed syllables than in unstressed syllables. Both the peaks and valleys of the velum time functions are affected by stress; the velum height valley is not as deep for nasal consonants in unstressed syllables as it is for stressed syllable and the velum height peak is not as high for oral consonants in un tressed syllables as it is for stressed syllables.

Stress: Effects on Velum Position and Velum Velocity

However, the two speakers differ in terms of velum movement velocity in unstressed and stressed syllables. For speaker B the velocity of velum lowering seems to be different in stressed syllables and in un-stressed syllables; the maximum slope of the downward movement of the velum appears to be greater in the stressed syllable than in the unstressed syllables. or speaker A, in contrast, the velocity of velum lowering does not seem to depend on the stress status of the syllable ; the maximum slope of the downward movement of the velum is similar and consistent in stressed and un-stressed syllables.

Modeling Velum Height and Velum Velocity

For speaker A, the feature [± strong] was associated with a considerable difference in velum position, but did not have any direct influence on the velocity of velum lowering. The two facts may be considered consistent with [START_REF] Lindblom | Spectrographic study of vowel reduction[END_REF] well-known model of coarticulatory undershoot. In that model, articulatory movements between given target positions are determined by a filter, representing the characteristics of the articulator as a physical system. Undershoot of a phonologically specified target position occurs if the time interval specified for a transition from one target position to another is short relative to another is short relative to the time constant of the smoothing filter. Assuming that there is no direct control of velocity (i. e., assuming the input to the filter is always a step function) and that the target positions are the same (i. e., assuming the step size is the same), an articulator's velocity is determined by the filter. Then we should not expect to find any difference in movement velocity until the next movement is triggered. For speaker A, in most cases, as depicted in figure 5, we may assume that the starting position for each of the stressed and unstressed conditions is nearly the same, and the descent toward the next target position follows the same time course, displaying the characteristic response function of the physical system, which may be assumed to be the same for each stress condition. The extent of descent varies reflecting only the timing of the next input step function for an ascent. The unstressed syllables are expected to exhibit shorter excursions of descending movement because the syllable duration is shorter than for stressed syllables and because the input step may be smaller. If this interpretation is correct, the coarticulatory explanation assuming a binary specification of [±nasal] can be successful in explaining velum lowering patterns for speaker A.

For speaker B the velum movement pattern as a whole crucially depends on stress condition in way not predicted by a simple smoothing filter model. In other words, for speaker B, we cannot ascribe the difference in velum movement patterns to coarticulatory undershoot. Different velocities and multiple target positions for the velum for consonant are assumed below in modeling speaker B's movement patterns, depending on the suprasegmental context, particularly the feature [± strong].

Rules for modeling speaker A

For speaker A, the velum has a target position 'low' for nasal consonants 'high' for oral consonants. There is no target position for vowels.

Rules for Velum Lowering

The velocity of velum lowering is fixed for all stressed 'CVNs and for all unstressed CVNs (fig. 5). The timing of velum lowering, relative to crucial articulatory events such as stop closure is determined by the manner of articulation of the consonant preceding the nasal and by the suprasegmental configuration. In order to infer the rules of relative timing, we summarized the articulatory events, depicting the relevant movement patterns for different articulators, along with the acoustic signal envelope, as shown in figure 6. An excerpt of such depicted patterns is also shown in figure 10, illustrating the timing rule for velum raising.

Figure 7 shows a flowchart for selecting the applicable rule for the beginning time of velum lowering in different case . The rules are: L1) Begin to lower the velum at the release of the top closure (if is a stop in /CVN/). L2) Begin to lower the velum at the middle of frication (if C is a fricative in /CVN/). L3) Begin to lower the velum at the acoustic onset of the consonant (elsewhere in /CVN/s, and in /CN/s is [strong]. L4) Begin to lower the velum in the middle of the vowel (in /V=N/; i.e. if the nasal occurs at the beginning of a word and the preceding word end with a vowel). L5) Begin to lower the velum at the onset of the oral closure for the (in /C#N/, i.e., if a phrase-final C is strong and is followed by a phrase-initial N).

Figure 8 illustrates applications of these rules in three example sentences for this speaker. In this figure, the observed articulatory and acoustic events are provided as the input data, and the applications of rules (shown in circles) result in the triggering of each velum descent, as indicated by each straight downslope. The vertical dashed line relates this trigger time to the occurrence of the rule relevant articulatory/ acoustic event. The descent slope is assumed to be constant and was determined empirically to fit the observed data. ---------------------------------------------Insert Figure 7 here ------------------------------------------------------------------------------------------Insert Figure 8 here --------------------------------------------- ---------------------------------------------Insert Figure 9 here Figure 9 gives a flowchart for selecting the applicable rule in different cases. There are two velocities for the velum raising movement: high and low. The velum raising movement is timed according to the segmental composition and stress status of the tautosyllabic vowel. Applications of these rules are exemplified in figure 10. The rules are: R1) Begin to raise the velum at a constant lead time before the peak of the oral articulatory closure for the nasal consonant (i.e., the peak of the blade pellet time function for /n/, the lip pellet for the labial consonant /m/, as indicated by the down arrows in fig. 10) in /VNC/. The velum raising movement has high velocity.

--------------------------------------------- Insert Figure 6 here ---------------------------------------------

Rules for Velum Raising

R2) Begin to raise the velum at the closure peak for the nasal consonant in /NəC/. The velum raising movement has high velocity.

R3) Begin to raise the velum at the closure peak for the nasal consonant in /NVC/. The velum raising movement has low velocity.

Slow Component for Raising Movements (Carry-Over Effect)

We assumed two different constant slopes for velum raising, depending on context, in applying the rules stated above. However in both cases, the slope changes to a very low value (i.e., slow velocity movement). We shall call this phase of movement a slow component (or 'carry-over' effect). Apparently, the breakpoint for the onset or this very low velocity movement occurs at a fixed pellet height, h, as shown in figure 11. ---------------------------------------------Insert Figure 12 here ---------------------------------------------Figure 12 displays the time course of the velum predicted by these rules. The reconstructed time functions are superimposed on the actual velum time functions for the three sentences. As can be seen for this figure, the small set of rules can describe with fair accuracy the observed time functions of velum movement for speaker A. These rules are obviously incomplete to cover all possible phonological environments of a nasal consonant in English. Also, the rules may need revisions when they are tested against more extensive data.

Rules for Modeling Speaker B

For speaker B, only one segmental feature [± nasal] and two suprasegmental features [± strong] and [± stressed] are used in the rules for reconstructing velum height, velocity and timing. In contrast to speaker A, no mention of the manner features of the oral consonant is necessary.

- --------------------------------------------Insert Figure 13 here Figure 13 shows a flowchart for determining the velum target height and figure 14 shows a flowchart for determining velum velocity and timing. Figure 15 displays a reconstruction of the velum movement time course in sentences for speaker B.

Target Positions

For speaker B, several target positions for the velum were specified, as shown in figure 13: ---------------------------------------------Insert Figure 14 here ---------------------------------------------(1) 'Very high' velum position for an oral consonant if followed by an oral consonant (CVC sequence), or if it is [+ strong).

(2) 'High' velum position for an oral consonant that is [-strong].

(3) 'Low' velum position for a nasal con-sonant preceded by an oral consonant and an unstressed vowel, or (not indicated on the chart in fig. 13) after a pause. In our data, in these environments, the nasal consonant is in initial position. (4) 'Very low' velum position for a nasal consonant in two environments: following another nasal consonant (a case of assimilation), and following a stressed vowel (the /CVN/ or /CaN/ where C is [+ strong] nasal consonant is syllable-final).

(5) No velum target position is specified for the last consonant preceding the nasal.

Velum Lowering

The rules for timing and velocity, as illustrated in figure 14, are: LI) Begin to lower the velum at the articulatory release of the strong consonant C in /CVN/ or /C#N/ where C is [+ strong] i.e., where the [+ strong] consonant is the last consonant preceding the nasal. The velocity is high. That is, the velum starts to lower for the next nasal consonant immediately after the last [+ strong] oral consonant preceding the nasal, as soon as there are no stronger oral consonant ahead.

L2) Begin to lower the velum immediately after the release of the last [+ strong] oral consonant preceding the nasal when there are [-strong] consonants intervening between the [+ strong] consonant and the nasal ; the velocity is low. Velum height for the intermediate [-strong] consonants is affected (i. e., lowered) by such anticipatory lowering. There seems to be, however, a minimum velum height ('high') allowable for the [-strong] oral consonants. The velum does not lower beyond this height for [-strong] oral consonants.

The involvement of the tongue dorsum for oral or nasal palato-velar obstruents seems also to have a heightening influence on velum height (see above). For the sake of simplicity, the influences of the tongue position have not been taken into account in the reconstruction shown in figure 15.

Velum raising

R1) Begin to raise the velum in NC clusters following a stressed vowel (VNC sequences): the movement velocity is very high.

R2) Begin to raise the velum in other cases, t ol1 owing a nasal Consonant ; movement velocity is high.

There is no slow, carry-over component for velum raising for speaker B. The height of the velum in the oral consonant following a nasal consonant is not systematically lowered, as has been observed for speaker A.

Speech-Ready Gesture

No attempt was made to reconstruct velum behavior at the very beginning of sentence, because it is difficult to model by rules (fig. 12, 15). Tensing of all the muscles involved in speech production -those controlling the velum included -prior to speech is a common observation. The so-called 'speech-ready' gesture is highly variable in nature. The elevation of the velum at the beginning of the sentence does not seem to be fixed relative to either segmental events or supra segmental events. Instead, it varies apparently randomly or according to nonlinguistic factors, as confirmed by the study of velum behavior in different utterances of the same sentence by the same speaker. The striking variations in timing of the speech-ready gesture for the velum have been discussed in a previous paper [START_REF] Vaissière | Comment on Abbs' paper[END_REF].

Summary and Conclusion

For both speakers, the nasality target for a nasal consonant is an open velopharyngeal port and that for an oral consonant is a closed velopharyngeal port. However velum height as such is free to vary over a certain range, while the velopharyngeal port is kept closed or open correctly. (The X-ray data does not provide direct evidence that the nasal passage is actually open or closed.) But the variation follows a system of rules as discussed here; some of the characteristics may be linguistically meaningful. In this sense, the feature [±nasal] provide only a gross characterization in velum height control.

Other factors such as the suprasegmental feature [± strong], may influence the height of the velum, the velocity of velum movement, and the timing of velum gestures. This result illustrates that it is very important not to study the segmental and suprasegmental, or configurational and prosodic aspects of speech separately. Tongue height also plays a role in determining the amount of anticipatory velum lowering in oral consonants, probably as a consequence of direct physical interference between the tongue and the soft palate.

Contrary to our expectations [START_REF] Lubker | An electromygraphic-cineflourographic investigation of velar function during normal speech production[END_REF][START_REF] Fritzell | The velopharyngeal muscles in speech: an electromyographic and cineradiographic study[END_REF][START_REF] Bell-Berti | An electromyographic study of velopharyngeal function in speech[END_REF][START_REF] Bell-Berti | Coarticulatory effects of vowel quality on velar function[END_REF], our rules did not need to include specific targets for the velum for vowels. Perhaps the importance of other factors in our data obscured the often observed relationship between velum height and vowel height. The influence of vowel height on velum height was observed only in utterances of nonsense syllables in our study. In contrast with previous observations (for German stops, Kohler, 1980), there was no obvious difference in velum behavior between fortis and lenis consonants.

The relative importance of suprasegmenta1 features on velum movement patterns seems to be speakerdependent. As illustrated here, the behavior of the velum may differ, depending on whether the syllable is stressed or unstressed. It is possible that some of the differences between speakers A and B, such as the existence of the slow component, stem from different placements of the velum pellet. In any case, the extent of the influence of the suprasegmental factors relative to segmental factors seems to be of primary importance for understanding the differences between speakers. More data are clearly needed in order to optimize the feature set (involving [± stress] as well as [± strong], currently) and the rule systems or their common framework.

The approximation of movement patterns using a set of rules seems to be a valuable method for investigating complex behavior, such as the behavior of the velum in utterances of sentences. The report should be taken as an indication of possible rules that may be needed to account for some details of articulatory gestures.

Table 1. Sentences used for data analysis

Figure captions:

(1) Ned went through some or the assigned material an unnecessary number of times.

(2) Ben announced that an innocent-seeming infant had nimbly nabbed most of the bananas.

(3) The training of the remnants of the regiment neared completion. (4) A complement of militant nymphs promptly demanded amnesty.

(5) The ugly turtle collided with the collie. (6) The Boston steersman dismissed the atlas. 13) John made an attempt to remind Natasha of the fanbelt's condition. Fig. 1. Lateral trajectories of pellets placed on the velum, the rear, mid, and front parts of the tongue surface, the lower tooth, and the lower lip, as displayed by an interactive analysis program (Miller and Fujimura, 1952). Broken lines indicating the outlines of the speech organs were added to a computer display. The bottom part of the display shows the speech signal envelope. The dotted crosses represent the pellet positions for the point in time corresponding to vertical dotted line through the speech signal. The three numbers at the top identify, from left to right, the utterance record number, beginning and ending frame numbers of the current display. The frame number for the current cursor position is shown below the sentence transcription to the right. The time interval between the frames is about 8.5 ms. (1) A complement of militant nymphs promptly demanded amnesty (32359).

(2) Ben announced that an innocent-seeming infant had nimbly nabbed most of the bananas (32357).

(3) John made an attempt to remind Natasha of the fanbelt's condition (32356). Fig. 3. Velum height traces for utterances with utterances with non-nasal consonants (C) and utterances with nasal consonants (N) for speaker B. Time functions for different utterances are superimposed. The total duration of the utterances shown here is approximately 0.3 s. The vertical scale indicates velum height (relative to the maxilla) approximately in millimeter. In figures 3a and b, the time functions are lined up at the maximum and minimum velum height, respectively. Two cases of an oral consonant immediately surrounded by two nasal consonants (the final consonant of the first word in the sequences 'assigned material' and 'regiment neared') have not been plotted, since there is no peak in the velum tracings corresponding to the intermediate oral consonant. Fig. 6. Observed time functions (from top to bottom): speech signal lower lip height, tongue blade height, rear tongue height, and phonetic transcription for an utterance of the sentence 'A complement of militant nymphs promptly demanded amnesty'. Irrelevant portions of trace are suppressed to clarify the lip movement for the labial consonants (/p/ /m/, and /v/), the tongue tip movement for the consonant (/t/,/n/, /1/, and /r/), and the rear tongue movement for the consonant /k/. In the phonetic transcription, ≠ and # indicate word boundary and phrase boundary, respectively. The total duration from left to right is about 2.4 s. Fig. 7. A preliminary algorithm for predicting the starting of velum lowering for speaker A, based on the information depicted in figure 6. The event in the box at the bottom of the figure deter-mine the starting time of the velum movement (see text). For example, in /CVN/, the velum tarts to lower at the acoustic release of the stop consonant. Fig. 8. Velum movement predicted by the rules in figure 7. Solid lines enclose the speech waveform envelope for the three sentences shown in figure 2. Dotted lines, broken lines, and dot-dashed lines represent the blade, the lips, and the rear tongue pellets, respectively. The downward line above the speech envelope represents the lowering gesture whose onset has been determined using the rules described in figure 7. (1) Ned went through some of the assigned material an unnecessary number of times ( 12446).

(2) Ben announced that an innocent-seeming infant had nimbly nabbed most of the bananas (12447).

(3) The training of the remnants of the regiment neared completion (12448). (4) A complement of militant nymphs promptly demanded amnesty (12449).

(5) The ugly turtle collided with the collie (12451). (6) The Boston steersman dismissed the atlas (12452).

  Figure captions:(1) Ned went through some or the assigned material an unnecessary number of times.(2) Ben announced that an innocent-seeming infant had nimbly nabbed most of the bananas.(3) The training of the remnants of the regiment neared completion. (4) A complement of militant nymphs promptly demanded amnesty.(5) The ugly turtle collided with the collie. (6) The Boston steersman dismissed the atlas. (7) Don't disturb the anemic anemone after dinner. (6) We need a leader with courage in distant places. (9) Beb commands it: Beb commends it ; Beb commits it. (10) Beb compacts it : Beb compels it : Beb competes it. (11) If you displace the dummies, you'll have to disarrange the display. (12) If you display the puppies, they'll be hard to displace. (13) John made an attempt to remind Natasha of the fanbelt's condition. Fig.1. Lateral trajectories of pellets placed on the velum, the rear, mid, and front parts of the tongue surface, the lower tooth, and the lower lip, as displayed by an interactive analysis program(Miller and Fujimura, 1952). Broken lines indicating the outlines of the speech organs were added to a computer display. The bottom part of the display shows the speech signal envelope. The dotted crosses represent the pellet positions for the point in time corresponding to vertical dotted line through the speech signal. The three numbers at the top identify, from left to right, the utterance record number, beginning and ending frame numbers of the current display. The frame number for the current cursor position is shown below the sentence transcription to the right. The time interval between the frames is about 8.5 ms.

Fig. 2 .

 2 Fig. 2. Velum height (upper) and speech signal envelope (lower) time functions for the utterances given in phonetic transcription. The three pairs of traces represent three different utterances by the same speaker (speaker A):(1) A complement of militant nymphs promptly demanded amnesty (32359).(2) Ben announced that an innocent-seeming infant had nimbly nabbed most of the bananas (32357).(3) John made an attempt to remind Natasha of the fanbelt's condition (32356).

Fig. 4 .

 4 Fig. 4. Velum height time function for oral consonants (C) for speaker B. In all cases, the first oral consonant is not immediately preceded by a nasal consonant. Here, the time functions in figure 3b are separated into two groups : (a) C is followed by a vowel (V) and another consonant (C) (b) C is followed by a vowel (V) and a nasal consonant (N). The time functions in Figure 4h are further separated into two subgroups. (c. d) the described in the text 'Configurational Environment'). In figure 4c, C has the feature value ([+strong]: in figure 4d, C has the feature value [-strong].

Fig. 5 .

 5 Fig. 5. Effect of stress on the velum lowering velocity in stressed versus unstressed CV sequences. The figure represents the vertical position of the velum pellet (vertical scale) as a function of time. The same scale has been used for the two speakers.

Fig. 9 .

 9 Fig.9. An algorithm for predicting the timing and movement velocity of velum raising for speaker A.

Fig. 10 .

 10 Fig. 10. Upward velum movements predicted by the rules in figure 9.

Fig. 11 .

 11 Fig. 11. Predicted velum movement for the slow component (see text), for the sentences in figure 10. Compare with figure 8.

Fig. 12 .

 12 Fig. 12. Predicted velum movement (dotted lines) superimposed on the observed velum movement time functions (solid lines) for three sentences produced by speaker A. The partial absence of a solid line for the middle sentence is due to pellet mis-tracking.

Fig. 13 .

 13 Fig. 13. Prediction algorithm for velum target position for speaker B.

Fig. 14 .

 14 Fig. 14. Prediction algorithm for velum movement velocity and p timing for speaker B.

Fig. 15 .

 15 Fig.15. Predicted velum movement (dotted lines) and observed velum movement for six sentences for speaker B (table 1, sentences 1-6). The duration for each line is about 4.3 s.

  Figure 1

Figure

  Figure 13 :

Acknowledgements

The articulatory data used for this study were obtained in experiments conducted by Osamu Fujimura of AT&T Bell Laboratories with the cooperation of the University of Tokyo staff. We particularly thank Shigeru Kiritani, Research Institute of Logopedics and Pediatrics of the University of Tokyo, in this connection. Win Nelson and Joan Miller of AT&T Bell Laboratories developed the programs for analyzing the data, and I am very grateful to them for their help and kindness. And above all, I thank Osamu Fujimura for his constant help and attention during the course of this research. Marian Macchi of Bell Communications Research, in cooperation with him, assisted substantially in the preparation of this report.

Hirose at the University of Tokyo. The data acquisition was conducted by Drs. O. Fujimura, S. - -------------------------------------


TER Ilya Meignan--Masson

Toward the Mixing of Monitoring and Profiling Operations on HPC Platforms

System and job monitoring are two established way of measuring the utilization of HPC systems. Due to the scale and complexity of modern HPC systems, users also require to profile their application using the actual system. However, the usage of many tools usually requires more effort in support and can also lead to significant under-performance. We introduce Colmet, an existing system modified to combine monitoring and profiling operations. Users can dynamically reconfigure Colmet at runtime and set a different sampling period for each metric. We present the results of the performance analysis conducted.

Introduction

Modern HPC systems are composed of hundred of thousands if not millions of cores linked to nodes dedicated to storage with a very efficient network. They require large amounts of money and effort for their design, their building and for their operation. At this scale, their efficient usage is an essential goal.

The scale of those system makes tedious any manual monitoring and motivates dedicated monitoring systems. Some systems focus on presenting to system administrators an overview of the status of components in the system (Nagios [START_REF] Imamagic | Grid infrastructure monitoring system based on Nagios[END_REF]). In this paper, we define monitoring as the process of collecting, gathering and storing data about the execution of one or several applications. Some systems also include some way of displaying and presenting the collected data to the user. This data can be drawn from various sources and will be called metrics. Metric sources include hardware counters (e.g., using perf in Linux systems), operating system performance data (virtual memory or kernel related metrics) or network performance data. In the recent years, the issue of energy consumption has become a crucial matter as the biggest systems now require tremendous amounts of energy to operate, for their nodes and their cooling [1]. Metrics are collected periodically. We call sampling period the interval between two moments in time where the system collects data.

Designing monitoring systems for HPC platforms has the following challenges. First and foremost, the monitoring system will interfere with user applications running on the nodes. To avoid tampering as much as possible the data, the overhead of the system needs to be as small as possible. Another issue is the complexity of the hardware and software architecture in the system. Monitoring systems designers must choose between two directions. Either supporting and optimizing a specific architecture, which often results in a custom monitoring tool for each system. Or building a portable system which might not suit perfectly the monitored platform but in return will not require much effort to use in many places. Once the data is collected, another issue is to present it in a usable and meaningful way.

In a cluster, users are allocated only a subset of the nodes in what is called a job. In the development process, users of the HPC systems might find the need to perform a profiling of their job. This set of techniques, that aims at optimizing the performance of an application, is crucial in the context of HPC. Due to the complexity of the HPC system, executing the application on another system will only yield limited results.

To answer this specific need, distributed application profiling tools have been developed like the ones we present in Section 2. Designers of such tools face the same challenges that we described for monitoring systems. Moreover, a system where every user would be running its own profiling tool next to the global monitoring system would most likely suffer significant under-performance.

In this article, we propose to use the same tool for the system administrators and their monitoring needs as well as for the users and their profiling needs. We implemented this in Colmet [2], an existing monitoring system previously developed in the Datamove team. We modified the tool so that it could collect data at different sampling period for each metric. Our system limit the performance overhead both on the compute nodes and on the network with simple string substitution. It is also rather portable with ZeroMQ sockets for communication which is very common but also cgroups v1 which requires a specific set of Linux kernels to be running on the compute nodes. We believe our system will simplify the usage for both system administrators and users. Users will be provided a integrated and optimized way of profiling their application and system administrators will only have this unique tool to maintain and to fit to the system.

The paper is structured as follows. Section 2 gives an overview of multiple distributed monitoring and profiling systems. Section 3 details the design and the architecture of our system. The performance analysis that we conducted is presented in section 4 before the discussion about limitations and future work in section 5.

Related work

Monitoring tools are often customized to the system they are supporting for performance reasons if not built for it. However, modern systems have some usual properties that we describe in this section. We also present some profiling software designed for distributed systems.

To begin with, we note the existence of a recent review of different monitoring systems for HPC systems [START_REF] Stefanov | A Review of Supercomputer Performance Monitoring Systems[END_REF]. This review identifies one of the problems of monitoring large scale HPC system as requiring multiple tools with "custom scripts to get comprehensive monitoring of the HPC system" (section 3.1). It also mention that "there is no single tool to achieve all aspects of monitoring, analysis and alerting" (also section 3.1). On this part, the work presented here addresses part of those aspects while adding a profiling use-case.

The type of metrics that a system can collect is one the most important properties. Some systems will collect CPU performance counters values, virtual memory statistics and kernel related metrics (LIKWID [START_REF] Röhl | LIK-WID Monitoring Stack: A flexible framework enabling job specific performance monitoring for the masses[END_REF], DCDB [START_REF] Netti | From facility to application sensor data: modular, continuous and holistic monitoring with dcdb[END_REF], Colmet [3]), I/O or filesystem related metrics (Beacon [START_REF] Yang | End-to-end I/O monitoring on a leading supercomputer[END_REF], Colmet, DCDB), temperature (ExaMon [START_REF] Bartolini | Paving the way toward energy-aware and automated datacentre[END_REF], Colmet, DCDB) or Network performance (Colmet, DCDB). Some metrics are collected inband (i.e. on the compute nodes) and some must be accessed out-of-band (ie. on service nodes or on the routers).

A desirable property for a monitoring system is the ability to store data as a time series. As their name suggest, time series are collection of measures performed at successive moments in time. They are particularly important to give informations on the changes of a metric over time. Most modern systems (e.g., NWPerf [START_REF] Mooney | NWPerf: a system wide performance monitoring tool for large Linux clusters[END_REF], LIKWID, Colmet) are capable of storing metrics as time series often using dedicated databases.

These time series can contain the collected data at different scopes.

• The node-scope gives insights on the performance of the compute nodes. This level is particularly interesting for the system administrators to keep track of the performance of the nodes over a period of time or to correctly scale the system to the needs of the users. Systems like Ganglia [START_REF] Massie | The ganglia distributed monitoring system: design, implementation, and experience[END_REF] or DCDB are examples with this scope. • Consequently, job-level scope combines data collected on the compute nodes and job data given by the Resource and Job Management System (RJMS) to store the data corresponding to each job. This enables HPC system users to obtain insights on the performance of their application over time. NWPerf and Colmet are examples of system at this scope. Data gathered with this scope is sometimes used to perform analysis on the job performance and help the non-expert user to understand the data. LIKWID is an example of such tool providing simple analysis to the users. For this scope, a link between the data collected on a compute node and the job assignation data must be made. Collection agent on the compute nodes can be given this responsibility like in Colmet. Other systems like LIKWID rely on the routers to tag the information as it passes by them.

Even though such analysis can be helpful, HPC users are most of the time more knowledgeable about their application than the monitoring system designers. With an API to dynamically change the metrics collected and the sampling frequency at runtime, such users could tune the monitoring along their needs all along the execution of the application. To our knowledge, no system implements this feature.

Due to the volume of metrics, especially when the number of nodes in the cluster or the sampling frequency increases, some systems like Beacon or Colmet compress the data on the compute node before sending it on the network. This compression is a trade-off between overhead on the network and overhead on the compute node. Even if it does not compress the data, NWPerf designers included a network contention analysis in [START_REF] Mooney | NWPerf: a system wide performance monitoring tool for large Linux clusters[END_REF].

The criteria presented here are summarized in Table 1.

Profiling tools usually instrument the code or the executable of an application to trigger data collection at specific points of the execution or rather attach to the process and collect both the metric values as well as the instruction that the process is executing. Parallel profiling tools like TAU [START_REF] Shende | The Tau Parallel Performance System[END_REF] follow the former principle while other tools like STAT [START_REF] Arnold | Stack Trace Analysis for Large Scale Debugging[END_REF] use the latter. They share with monitoring systems some of their properties like the way metrics are collected (in-band or out-of-band), the time series storing, the per-job profiling.

Colmet : Design and Architecture

In this section, we present the design of Colmet with our modifications. Colmet was first designed after the publication of [START_REF] Emeras | Analysis of the jobs resource utilization on a production system[END_REF] in 2013. After a first version in Perl, a version has been written in Python [3] that is available in production on Grid5000 clusters. To improve the performance, a Rust re-writing had been started but was left unfinished. To mix operations of monitoring and of profiling, the tool needs to collect metrics at a small sampling period. It motivated the refactoring and improvement of the code to leverage the speed of Rust.

In terms of design, our version is close to the Python one on a lot of points. First and foremost, the tool uses some of the same backends (i.e., an abstraction of the protocol or algorithm needed to actually perform the data capture). Backends related to cgroups ( [START_REF] Flauzac | A review of native container security for running applications[END_REF]) metrics and to perf events are implemented but we intend to add the others like temperature, energy consumption (with the Running Average Power Limit -RAPL [START_REF] David | Rapl: Memory power estimation and capping[END_REF] metrics), or Infiniband related metrics. We note here that the Python version was using the taskstats interface, which is a kernel interface to access data about a process and we changed it to leverage cgroups created by the RJMS. It also reuses Elasticsearch as a time series database. The scope of the monitoring is at the job-level. It relies here on the RJMS to creates a cgroup for the job on the allocated compute nodes. Colmet is composed of two components : a node agent and a collector. Figure 1 presents the architecture that we describe in the following.The node agent runs on every node in the system. It stores the list of metrics to gather on this node as a tuple containing the metric, the job id and the sampling period. This relation allows to collect data at different sampling frequency for each metric and to link the data with the jobs. This means that some metrics can be set by system administrators on all the jobs for monitoring with a low sampling frequency. Typical values are at the order of the second. While other metrics can be set by the users for their job with possibly a higher sampling frequency, especially under the second. To this end, users can interact with a script to update the metrics to collect or the default sampling frequency on the compute nodes. Each time some metrics should be collected, the node agent queries the corresponding collection backends. The cgroup backend read the data respectively in the files cpu.stat and memory.stat created in the virtual directory of the cgroup. For perf event data, they are also read from a file created by a call to perf event open(). Data is then processed and sent to the collector. To avoid using too much the CPU, the processing simply substitutes an id in place of the metric name, sending shorter messages over the network. Unlike the node agent, the node collector is written in Python as we believe that the overhead on the collector node is not a significant performance issue. The collector listen on a port, receives the data and place it in the Elasticsearch database. This means that it requires a server to run. Communication between the node agents and the collector is done with ZeroMQ (ØMQ) sockets.

Ganglia

Performance analysis

An important part of this work was dedicated to the measurement of the overhead caused by the monitoring and profiling system. We focused on the execution time overhead. We studied the effects of the following parameters : the number of metrics, the sampling period and the version of Colmet (including a special value without Colmet launched) using 2 benchmarks from the NASA Parallel Benchmarks (NPB) suite as applications. In this section, we describe our experimental set-up, present the results we obtained and the analysis that guided our exploration.

Experiment setup

The first benchmark is the Embarrassingly Parallel (EP) benchmark. This benchmark accumulates statistics from pseudo-randomly generated numbers. It is compute-intensive with almost no communication and thus provides an estimate of the upper achievable limits for floating-point performance. The second one is the LU Simulated Computational Fluids Dynamic (CFD) application (LU) benchmark. It is intended to accurately represent the principal computational and data movement requirements of a real CFD application and thus is closer to a real application than the previous one (cf. the specification of the NPB [START_REF] Saini | NAS Parallel Benchmark Version 1[END_REF], in Sections 3.1.1 and 3.2.1). We used the MPI version of both benchmark. The size of the benchmark is quantified by the class. For each number of node, we chose a class that would be long enough to lower the influence of noise in the result.

All our experiments were carried out in Grid5000 [START_REF] Cappello | Grid'5000: a large scale and highly reconfigurable grid experimental testbed[END_REF], more precisely in the dahu cluster of the Grenoble site. Each node Figure 1: Architecture of Colmet in this cluster is provisioned with 2 Intel Xeon Gold 6130 CPUs with 16 cores each, 192GiB of memory and a 240GiB of SSD (Samsung MZ7KM240HMHQ0D3). All the nodes are also connected using a simple TCP network as well as a 100 Gbps Omni-Path network. Nodes are running Debian 11 with OpenMPI 4.1.0. Experiments were run using all the cores of the compute nodes but not all the hyperthreads because it would not improve the performance as HPC application computations mainly involve floating point operations and there is usually one Floating Point Unit (FPU) per core.

One technical point worth noting is the usage of the Nix package manager during the experiments to install the required software [START_REF] Dolstra | A safe and policy-free system for software deployment[END_REF]. This package manager provides strong guarantees on reproducibility. A package is specified by a definition containing the name and the version of all its dependencies. Nix then installs each package in its own environment ensuring that it can run properly as it was designed and tested. A system of linking helps installing each package only once to prevent from using more space than required on the disk.

Usual monitoring systems collect metrics with sampling periods in the order of second. For our use-case of profiling, the monitoring tool should be able to collect metris with a sampling period under the second while keeping the overhead as low as possible. We chose values for the sampling period between usual values like 5 seconds up to the millisecond because we believe that smaller values would not be of interest to the users.

We call configuration a set of values for the parameters. The experimental script [4] requests a job, install the required softwares and executes a set of configuration specified in a YAML file. Because a single experiment can require multiple version of Colmet, an experiment requires multiple run of the script, one for each version. To obtain results as accurate as possible, each configuration is replicated 10 times. We also randomize the order in which the configurations are executed to minimize the observational error. With the 10 values for each configuration, we compute the mean and the confidence intervals with a coefficient of 95%.

Experimental results and analysis

We present the result of 2 experiment that we performed to measure the impact of Colmet on the performance of the system.

The first one quantifies the impact of the number of metrics collected by the node agent. We ran the LU (class D) and the EP (class E) benchmark with a sampling period of 0.001, 0.01, 0.1, 1 and 5 seconds on 8 nodes (+1 for the collector). We then compared three version of Colmet : the Rust version with only 1 metric, the Rust version with 67 metrics (i.e., the greatest amount the system can gather) and the reference without Colmet. An ideal system would not be impacted by the number of metrics. Results are presented in Fig 2a The results for the first experiment clearly show that the number of metrics do not have a significant impact on the execution time of the application. The only configuration where the confidence interval are not intersecting is with the LU benchmark at 0.001 second in sampling period.

For the second experiment we do not have, at the moment, enough results to conclude on the hypothesis. Similarly to the first experiment, the execution time with the EP benchmark are spread on a large interval which leads to larger confidence intervals. An hypothesis for this behavior could be that the benchmark is compute-intensive and that Colmet then could be affected to a different core each time leading to nonuniform disruption of the execution which causes such difference between time of the same configuration. On the other hand, the results collected from the experiment with the LU benchmark correspond to the hypothesis. We also observe that, at 1 and 5 seconds, the Python version becomes significantly better than the Rust version. For the moment, we are not able to explain this part.

5 Future work

Limitations

Considering the problem, namely mixing monitoring and profiling in the same tool, the biggest limitation that this tool exhibits is the impossibility to instrument the code of the application. While this tool is perfectly capable of collecting metrics at a very high sampling frequency, correlating the gathered data with the instructions of the application will be hard if not impossible.

Even if we implemented the mechanism of dynamically changing the configuration of the node agent, it is not production ready especially in terms of permissions. The lack of proper permission system allows users to change the configuration of the node completely.

Another limitation on the design is the scalability. All our experiments were using a relatively small number of nodes with regard to the usual size of small to middle size clusters. At a larger scale, our choice of simple compression might not be enough to prevent an important overhead on the network because all the node agents are sending to the same collector.

Directions

First and foremost, future work should include a thorough and precise performance analysis especially on the improvement compared to the Python version. We believe that such analysis will require to carefully choose other application and to study their behavior. This analysis should also investigate the scalability consideration. To continue on the compression part, we believe that a design close to what has been done in Rezolus ([5], section Background) could improve the scalability of our design. The idea is to locally collect metrics at a high frequency and then accumulate the data over a moving window. The resulting data is then transmitted at a lower frequency to the collector. In our case, we could gather monitoring with this principle while keeping the current sending of all the data for the profiling.

In terms of permissions, we could have multiple group of metrics. Each group could be accessed by users with a certain permission. This could forbid users from modifying the metrics defined by the system administrators but also other users jobs. Building on the topic of security, we also note that it would be rather simple to add encryption of the data before sending on the network. But it would increase the overhead on the CPU of the compute nodes.

Another point that could be the object of future works is the integration with other monitoring or profiling tools. As detailed in Section 2, the subject of monitoring has been already extensively studied and many tools already implement many efficient way of collecting data. Future works could include finding a abstraction to include other tools as backends for Colmet. This would be a challenge in that other tools were not designed to support a different sampling period for each metric but it could allow Colmet users to benefit from works on other backends (not already implemented) or on data visualization and analysis (which could be useful for both monitoring and profiling).

Conclusion

In this paper, we presented the design and our implementation of Colmet, a tool capable of collecting metrics at different sampling period and to tag them with information on the job, allowing for job-scope monitoring and profiling with the same system. After giving design and implementation details, we presented the results of the performance analysis experiments performed. Experiments were carried out using multiple benchmarks from the NPB suite and studying the effect of the sampling period, of the number of metrics and of the version of Colmet on the execution time. The results show that the number of metrics does not cause a significant overhead on the execution of application. We also discussed some directions for future work.

  on the upper row. The second experiment compares the Rust version with the Python version. To do so, we used the same values as the first experiment and ran them with the Rust version with 67 metrics, the Python version and without Colmet as a reference. Here, our hypothesis is that the Rust version should be better than the Python version but worse than the reference. Results are presented in Fig 2a on the lower row.

Figure 2 :

 2 Figure 2: Results of the experiments The blacked dashed line is the reference without Colmet and the two red dashed line are the upper and lower value of the confidence interval.

Table 1 :

 1 Summary of the state of the artBlank means that there is no information on the related paper. Metrics : HC : Hardware Counters, OS : proc and sys virtual filesystem related, VMstats : virtual memory related, IO : filesystem related, network : high-performance network related

					NWPerf	ExaMon	Beacon	LIKWID	DCDB	Colmet	Colmet
										(Rust
										version)
	Metrics		VMstats,	HC/perf	HC/perf	IO	HC/perf,	VMstats,	HC/perf,	HC/perf
			OS	and	and VM-	and En-		VMstats,	OS,	IO,	VMstats,	and VM-
			custom		stats	ergy		network	and energy	network,	stats
								and IO		energy and
										IO
	Scope		node		job	node	node	job	node	job	job
	Dynamic								yes
	reconfigu-
	ration
	Features					Integrated	Integrated		Visualiza-	profiling
	for	the					analysis	online and		tion GUI	use-case
	users						and visu-	offline
							alization	analysis
	Performance	Hierarchical	performance		Compress-		performance	ID substi-	ID substi-
	considera-	collection	analysis		ion		analysis	tution	tution
	tions