Koen Hillewaert
email: koen.hillewaert@cenaero.be

Nicolas Chevaugeon

Philippe Geuzaine

Jean-Fran Cois Remacle

Cenaero

Hierarchic multigrid iteration strategy for the discontinuous Galerkin solution of the steady Euler equations

Keywords: discontinuous Galerkin, Euler equations, quadrature free, multigrid, Newton-Krylov method

We study the e cient use of the discontinuous Galerkin ÿnite element method for the computation of steady solutions of the Euler equations. In particular, we look into a few methods to enhance computational e ciency. In this context we discuss the applicability of two algorithmical simpliÿcations that decrease the computation time associated to quadrature. A simpliÿed version of the quadrature free implementation applicable to general equations of state, and a simpliÿed curved boundary treatment are investigated. We as well investigate two e cient iteration techniques, namely the classical Newton-Krylov method used in computational uid dynamics codes, and a variant of the multigrid method which uses interpolation orders rather than coarser tesselations to deÿne the auxiliary coarser levels.

INTRODUCTION

The discontinuous Galerkin ÿnite element method (DGFEM) has become a very popular method for the computation of unsteady ows, especially when accurate solutions are needed. This is due to the combination of an arbitrary order of accuracy with data and algorithmic locality. Since the interpolation functions are deÿned independently in each element the mass matrix has a block-diagonal structure and may be inverted directly. The method can be easily parallelized because DGFEM stencils do not grow in size with increasing order. The interpolation functions are not restrained by continuity requirements and can hence be chosen freely, allowing easy implementation of, e.g. spectral interpolation bases. Another consequence is the inherent capability of the method to handle h-and p-adaptation.

The method has not experienced a similar interest for the computation of steady ows. On one hand, the accuracy requirements are not as stringent for these computations. One does not care about temporal dispersion and dissipation, so one can use mesh reÿnement without too severe consequences on computational e ciency. On the other hand, DGFEM still has some evident drawbacks in terms of computational e ciency. First of all an e cient iteration method is still lacking. Furthermore, the method seems to be very expensive in terms of memory (more unknowns for the same formal order) and work load due to quadrature. Both are exacerbated by the apparent need for isoparametric elements near curved boundaries.

Fidkowski et al. have recently shown [START_REF] Fidkowski | Development of a higher-order solver for aerodynamic computations[END_REF] that in terms of CPU time for a given precision it is better to increase DGFEM solution order rather than mesh resolution, at least for smooth solutions. Hence it is likely that the high-order DGFEM method will grow to be more e cient than a state-of-the-art ÿnite volume solver for steady computations as well, especially if it is combined with h-and p-adaptation. Another advantage is that accurate results can be obtained with-by today's standards-very coarse meshes, thus relaxing the requirements for the mesh generators (which often are not parallelized and are hence limited in mesh size) and diminishing the overhead related to storage of mesh topology.

In this paper, we investigate methods for increasing computational e ciency in terms of computing time and memory requirements, both by algorithmical simpliÿcations and e cient iteration strategies.

The paper is organized as follows. First DGFEM is reviewed brie y. After that we discuss the application of two algorithmical simpliÿcations of DGFEM which decrease both computing time and memory requirements for the evaluation of the discretization residual. In the following section we discuss and compare several e cient iteration strategies for stationary DGFEM.

THE DISCONTINUOUS GALERKIN METHOD

We focus on the application of DGFEM to approximate the steady solutions u to the Euler equations of gas dynamics in a domain D. We note the Euler equations concisely as @u @t

+ ∇ • f(u) = 0 (1)
DGFEM approximates the solution u in the broken or discontinuous interpolation space U that is spanned by shape functions i . Each of these functions i is deÿned on one element T only of a tesselation T of domain D. The approximate solution ũ is as such discontinuous across element boundaries.

ũ = i u i i ; i ∈ U (2)
The expansion coe cients u i are found by requiring that the approximate solution ũ satisÿes the following weak formulation of (1):

@ ũ @t + ∇ • f( ũ); v = 0; ∀v ∈ U (3)
where the inner product in the broken space U has been deÿned as

(p; q) = D p q dV ∀p; q ∈ U (4)
Considering that any v ∈ U is a linear combination of shape functions i , and any shape function i is only supported on one element T we can rewrite (3) as

D j @ ũ @t dV + D j ∇ • f( ũ) dV = 0 ∀ j ∈ U = T j @ ũ @t dV + @T j f( ũ) •ñ dS - T ∇ j • f( ũ) dV ≈ T j @ ũ @t dV + @T j H( ũ+ ; ũ-;ñ) dS - T ∇ j • f( ũ) dV (5)
Since the solution is discontinuous across element boundaries, the boundary ux f( ũ) •ñ corresponds to the solution of a local Riemann problem. This ux is thus replaced by a numerical ux function H, depending on the solutions ũ+ and ũ-approaching the boundary from the inside and outside, respectively, and the boundary normal ñ. If this ux function is positive, it provides the correct amount of stabilization. The e ect of the stabilization is the damping of modes that are unresolved by the chosen discretization. For the Euler equations we typically use the ux di erence splitting scheme of Roe [START_REF] Roe | Approximate Riemann solvers, parameter vectors and di erence schemes[END_REF].

ALGORITHMIC IMPROVEMENTS

The DGFEM is a relatively expensive method since more unknowns are required for the same formal order than for instance high order continuous ÿnite elements. Moreover, exact quadrature involved in the correct integration of the weak form is even more penalizing. The number of term of each block of the system grows as p 2 in two dimensions (p 3 in 3). Each term is computed as the integral of a polynomial of order 2p +1 and the number of quadrature points for evaluating that term then grows like (2p + 1) 2 ((2p + 1) 3 for three dimensions). Furthermore, for each quadrature point on element-interfaces, an expensive ux function needs to be computed. The apparent necessity for (near)isoparametric elements adjacent to solid domain boundaries not only complicates the method but also increases the number of quadrature points. Hence, the simpliÿcation and acceleration of the quadrature method is an important step towards an e cient method.

Boundary representation

Bassi et al. [START_REF] Bassi | High-order accurate discontinuous Galerkin ÿnite element solution of the 2D Euler equations[END_REF] have shown that an accurate boundary representation is necessary to maintain the formal order of the scheme and to avoid spurious production of entropy at the intersection of boundary faces. They stated that an equal order interpolation of solution and coordinates for elements adjacent to the boundaries is unavoidable. This requires the use of a higher number of quadrature points to account for the variation of the metric. In a novel approach, suggested by Baumann et al. [START_REF] Baumann | An adaptive-order discontinuous Galerkin method for the solution of the Euler equations of gas dynamics[END_REF], the elements adjacent to the domain boundaries remain simplices but an accurate representation of the boundary normal is used in the quadrature points. Krivodonova et al. [START_REF] Krivodonova | High-order accurate implementation of solid wall boundary conditions in curved geometries[END_REF] subsequently implemented di erent variants of this approach and showed numerically that the formal order of the method is well conserved. Using only simplices instead of boundary ÿtted elements represents a huge simpliÿcation of the code, and tremendous savings in memory and computing time.

In our implementation the normal to be applied at the Gauss points is found by projecting the Gauss point on the true surface, which is described by an analytical expression. The projection uses an iterative procedure such that it's direction is aligned to the computed normal.

Quadrature free implementation

In Reference [START_REF] Lockard | E cient implementations of the quadrature-free discontinuous Galerkin method[END_REF] Lockard et al. further elaborate the quadrature free implementation of DGFEM as put forward by Atkins et al. [START_REF] Atkins | Quadrature-free implementation of discontinuous Galerkin method for hyperbolic equations[END_REF]. This method avoids explicit quadrature and results in a signiÿcant reduction of computational time. Since this method uses only expansion weights, and does not explicitly reconstruct values at the quadrature points, the number of ux evaluations is drastically reduced. Furthermore, all other operations may be cast in matrix-vector and matrix-matrix products for which highly optimized libraries exist.

The method starts from the expansion of the ux in terms of the base functions:

f( ũ) = i fi i (6)
If we assume the elements to be simplices, the metric of the element is constant. With these assumptions the volume term can then be rewritten as

T ∇ j f( ũ) dV ≈ fi T ∇ j i dV = fi • J -1 T • T ∇ j i d dÁ J T = fi • J -1 T • Cij J T (7)
where J T is the Jacobian of the mapping of the element T onto the reference element and J T is its determinant; ∇ deÿnes the gradient in parameter space:

J T = @(x; y) @( ; Á) (8) ∇a = @a @ k ẽ k (9)
The vector of matrices C is independent of the element, and may be precalculated for any order. For each element T , we need to precalculate the metric terms J -1 T and J T . Once the expansion coe cients fi have been determined, the evaluation of the volume term is reduced to matrix-matrix and matrix-vector products of relatively small dimensions. A similar development may be done for the surface term of the boundary ux. The weights fi are found trivially if the ux vector is a linear function of the approximate solution ũ. Atkins proposes an elegant method to project the uxes of the Euler equations for a perfect gas. Unfortunately, a similar method is di cult to ÿnd for generic ux functions, generic equations of state, or even for the expansion of the boundary ux H of the perfect gas Euler equations. Given the computational cost of the evaluation of the boundary function in case of the Euler equations, the quadrature free method loses much of its interest if it does not apply to the boundary terms as well.

If we use a nodal expansion basis we may however assume that a function which passes through the nodal values of the ux vector is a reasonable approximation of the ux expansion in the volume as well as on the boundary. The expansion is then trivial, since we compute the nodal weights fi pointwise from u i :

fi = f(u i ) (10)
Apart from the apparent deviation of the formal order of the scheme, a major problem of this approach is that the interpolation weight vectors u i are only weakly coupled numerically. This is because the projection of the ux [START_REF] Keyes | Newton-Krylov-Schwarz methods: interfacing sparse linear solvers with nonlinear applications[END_REF] uses ux values that are computed using only one expansion coe cient vector u i . When the ux becomes singular, which is the case at stagnation and sonic points, the weights u i may vary randomly without modifying the computed ux. Using normal quadrature eliminates this problem since any degree of freedom (DOF) is used for ux evaluations at di erent Gauss integration points. The instability is illustrated in Figure 1, which shows a zoom of the leading edge region of the NACA0012 airfoil. The free-stream conditions are M = 0:3 and angle of attack of 2 • and the solution has been approximated with second-order polynomials; a Newton-Krylov implicit iteration scheme has been used for the solution of the discretized equations. The quadrature free solution is obtained just before explosion of the computation.

A very simple ÿx consists in reverting to classical quadrature in elements where the equations are (nearly) singular, i.e. in which stagnation or sonic points occur. For simple subsonic ows, such as presented in this paper, stagnation points occur exclusively on solid boundaries. A ÿrst obvious step is hence to use full quadrature consistently on all solid boundaries. This ÿx does not penalize overall e ciency, since the number of a ected cells is small and we can easily a priori group all elements using quadrature free or full quadrature integration, respectively. Furthermore, boundary elements already need a special treatment to take boundary curvature into account into the integration. For more general ows, and especially for transonic and subsonic ows, an indicator based on the annihilation of one of the characteristic speeds should su ce. This strategy however introduces a switch-like behaviour, so one needs to look into a robust integration into the (implicit) iteration strategy. Furthermore, some e ciency loss will occur since the elements in which we need to apply full quadrature are not known a priori. These aspects have not yet been investigated since the goal here was to demonstrate the instability and identify its source.

ITERATION STRATEGIES

One of the main obstacles towards using DGFEM for stationary computations is the lack of an e cient iteration strategy. DGFEM has up to now been used primarily for the explicit and fully resolved computation of unsteady ows. To compute steady ows those explicit time-integrators are prohibitively slow, even when using locally deÿned optimal timesteps.

Newton-Krylov

The Newton-Krylov method is fairly often used in ÿnite volume methods [START_REF] Geuzaine | Newton-Krylov strategy for compressible turbulent ows on unstructured meshes[END_REF][START_REF] Brown | Hybrid Krylov methods for nonlinear systems of equations[END_REF][START_REF] Keyes | Newton-Krylov-Schwarz methods: interfacing sparse linear solvers with nonlinear applications[END_REF] and has only recently been applied to DGFEM [START_REF] Rasetarinera | An e cient implicit discontinuous spectral Galerkin method[END_REF][START_REF] Bassi | GMRES discontinuous Galerkin solution of the compressible Navier-Stokes equations[END_REF]. In this method the discretized non-linear set of equations is solved using Newton iterations. The linear system arising from Newton linearizations is iteratively solved with a Krylov subspace iteration method. In this case we use the classical combination of GMRES with ILU(0) as preconditioner. In this method we add a pseudo-time derivative to the original residual to improve the conditioning of the linear system:

ũn -ũn-1 n + ∇ • f( ũn ); v = 0 ( 11
)
The pseudo-timestep n is chosen locally to conform to a given CFL number, and thus provide the 'locally optimal' diagonal dominance. Since this amounts to solving a pseudounsteady problem the non-linear convergence is governed by the pseudo-timestep, and will be obviously faster as n becomes larger. Therefore, the CFL is updated following each Newton update, inversely proportional to a power of the current reduction of the L2-norm of the residual :

CFL n = max CFL 0 • 0 2 n 2 ; CFL max (12)
A major drawback of the Newton-Krylov-ILU method for high-order DGFEM is the large memory footprint. The elementary blocks composing the preconditioning matrix increase quadratically in size with the number of unknowns per element=face. If p is the polynomial order, then the number of unknowns per element grows like p 2 in 2 and p 3 in three dimensions. Then, the number of elements in each block of the ILU grows like p 4 resp. p 6 . Hence a matrix preconditioner becomes quite impractical even for relatively low orders.

A second problem concerns the application of the slope limiter. Classical techniques of ux limiting are not applicable for high order DGFEM because of the presence of volume terms in the formulation. Hence the slope limiter is not integrated in the computation of the residual, but e ectively acts as a post-processing ÿlter [START_REF] Cockburn | The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems[END_REF]. Such a ÿlter is easily integrated in an explicit, but not into an implicit iteration strategy, such as the Newton-Krylov method.

Finally, it is necessary to tune the evolution of the CFL number in order to maintain su cient diagonal dominance throughout the computation. This becomes particularly di cult for higher orders.

Multigrid

With most iterative methods more iterations are needed to converge as the number of degrees of freedom increases. This is due to either a poorer conditioning of the global system for implicit iterative solution methods, or smaller stable timesteps in the case of explicit time integration.

The multigrid method overcomes this problem by using a sequence of ever coarser representations of the discretized problem, in the hope of keeping the same convergence rate per iteration as for the coarsest representation. If the work needed to solve the problem is directly proportional to the number of degrees of freedom, one attains the textbook multigrid e ciency (TME). TME requires constant residual reduction per multigrid cycle, independent of the ÿnest representation of the solution. Multigrid methods succeed in doing this because

• High frequency errors tend-as a general trend-to be well noticed by a discretization method. Hence iterative solution methods tend to reduce (dissipate) those errors more e ciently. Multigrid uses successive coarser representations of the problem to convert low frequency errors to high frequency errors on those coarser levels, thus leading to a similarly fast elimination of the low frequency errors. • The relaxation parameters (e.g. a pseudo-timestep) have higher stable=convergent values with respect to the physical dimensions on the coarser levels, permitting faster evolution of the transient. • The overhead in terms of work and storage on the coarser levels is negligible in comparison to the ÿnest level.

These features allow multigrid to achieve competitive convergence rates with relatively simple iteration methods. Multigrid is as a consequence a very e cient method, with respect to both memory and computing time, especially when a large number of degrees of freedom is considered.

Classical multigrid methods use di erent, successively coarser tesselations of the domain. This tends to lead to complicated methods for the transfer of solutions between di erent representations and ad hoc deÿnitions of the defect correction equations. The complexity of the implementation is even more critical in parallel. Classical multigrid has recently been applied to DGFEM discretizations for the groundwater equations by Bastian et al. [START_REF] Bastian | Multigrid for higher order discontinuous Galerkin ÿnite elements applied to groundwater ow[END_REF]. In the context of hierarchic higher-order methods such as DGFEM, one may as well use lower order interpolation spaces based on the same tesselation to deÿne the coarser representations of the problem. The advantages are the ease of implementation and the possibility to deÿne a ÿrm theoretical basis for the deÿnition of transfer operators and the defect correction equations. A similar technology has been developed independently by Fidkowski et al. [START_REF] Fidkowski | Development of a higher-order solver for aerodynamic computations[END_REF]. The rationale for ÿnding the interlevel transfer operators presented in the current work is however more general and is applicable to both p-multigrid and classical multigrid for any weak formulation using discontinuous interpolation, even when interpolation spaces are not nested.

Full approximation storage.

In the basic form of multigrid, the two-level method, we transfer the current ÿne level iterant to a coarser level. On the coarser level a defect correction is computed, which is then transferred back to the ÿner representation. The full approximation storage (FAS) method is the variant of this scheme for non-linear equations. The description of the basic FAS cycle may be found in the classic textbook [START_REF] Wesseling | An Introduction to Multigrid Methods[END_REF].

We have redeÿned the FAS cycle for the case of hierarchic multiorder or p-multigrid iterations as follows. Suppose we want to ÿnd the approximated solution ũp to (1) in the space U p of order p, based on the tesselation T. First we deÿne the functional residual operator L() as

L( ũ) = ∇ • f(u) (13)
Consider then another interpolation space U q which is based on the same tesselation T, and where q¡p. The two-level FAS cycle then proceeds as follows:

1. Pre-smoothing: perform a number of iterations on level p, leading to solution ũp * . 2. Restriction: approximate ũp * by ũq * in U q . 3. Defect correction: solve exactly in U q :

(L( ũq ) -L( ũq * ) + L( ũp * ); v q ) = 0 ∀v q ∈ U q (14)

4. Prolongation: represent the correction ũq -ũq * in U p and add to ũp * : ũp = ũp * + P p q ( ũq -ũq * )

5. Post-smoothing: perform additional iterations on level p to smooth the corrected solution.

We can interpret the function of this implementation of the defect correction equation in the following manner: we eliminate the part of the residual which is representable in the lower order space. Suppose U q ⊂ U p , then we deÿne U p-q := U p \U q (16)

We may then decompose ũp ũp = ũq + ũp-q ; ũq ∈ U q ; ũp-q ∈ U p-q (17)

If the Galerkin weighted residuals are linear in the expansion coe cients, we may rewrite the defect correction weighted residual ( 14) as (L(u (p-q) * + ( ũq -ũq * )); v q ) = 0 ∀v q ∈ U q (18)

We see that the defect correction equation deÿnes the lower order correction to the solution that removes the corresponding order variations of the residual. Since both interpolation spaces use the same elements, the deÿnition of the defect residual and the prolongation=restriction are trivial. The prolongation and restriction of the solutions both use Galerkin or L 2 -projection. This proceeds as follows: let us consider a solution ũa ∈ U a and its Galerkin projection

I b a ũa = ũb ∈ U b ũa = i u a i a i ; a i ∈ U a I b a ũa = ũb = j u b j b j ; b j ∈ U a (19)
Orthogonalizing the di erence ũa -ũb to the space U b deÿnes the following set of equations for the expansion coe cients u b i :

i ( b k ; b j )u b j = j ( b k ; a i )u a i ∀ b k ∈ U b (20)
The solution transfer operator I b a has a discrete or matrix equivalent I b a deÿning the transfer between the expansion vectors u a = [u a 1 : : :

u a na ] T and u b = [u b 1 : : : u b n b ] T u b = I b a • u a I b a = (M b ) -1 • M ba (21)
where

M b ij = ( b i ; b j ) M ba ij = ( b i ; a j ) (22)
The 'restriction' of the residual vector follows directly from the weighted defect correction equation. Conventionally one goes the other way around: ÿrst a restriction operator is deÿned for the residual vector, and then-in the best of cases-the coarse grid operator is found by applying the discrete operators to the ÿne grid operator (Galerkin coarse grid approximation) or-more frequently-one uses the same discretization technique on the coarse representation (discrete coarse grid approximation). The deÿnition of the forcing term requires the computation of the 'restricted' residual (L( ũp ); v q ) (23)

To compute this term explicitly we would need to redeÿne routines for weighting all terms of the residual deÿned in U p with test functions in U q . To avoid this complication, we ÿrst expand L( ũp ) in U p :

L( ũp ) ≈ i r p i p i i r p i ( p i ; p j ) ≈ (L( ũp ); p j ) (24)
This projection requires the mass matrix and the Galerkin weighted residual deÿned on space U p , both of which are already available. The 'restricted residual' is then computed as

(L( ũp ); v q ) ≈ r p i p i ; v q (25)
We ÿnd the following matrix operator connecting the vector containing Galerkin weighted residuals for space U p to the restricted residuals in space U q :

Ĩ q p = M qp • (M p ) -1 (26)
The discrete form of the prolongation and restriction operators as given by Equations ( 21) and ( 26) are the same as proposed by Fidkowski et al. [START_REF] Fidkowski | Development of a higher-order solver for aerodynamic computations[END_REF] in the framework of p-multigrid DGFEM for nested interpolation spaces. The rationale of ÿnding the transfer operators indicated here is however more general. In the framework of discontinuous interpolation it is equally applicable to non-embedded spaces and even to classical multigrid methods with nested or non-nested grids. Hence this deÿnition may be exploited for hp-adaptive and spectral implementations. Theoretically one may apply it to continuous interpolation spaces as well. However, the transfer operators require the inversion of the mass matrix on both levels, which will probably be prohibitively costly.

Multigrid strategies.

The two-grid algorithm requires an accurate solve of the defect correction equation. This may still be untractable, and hence a number of strategies exploiting the two-grid cycle are deÿned (see Figure 2):

• Explicit v-and w-cycles: A ÿrst strategy consists in replacing the coarse exact solve by a recursive application of the two-grid cycle. Now we base the restricted residual on the defect correction equation. Depending on the number of two-cycle solves per defect correction we distinguish v-cycles and w-cycles. We deÿne our multigrid cycling strategy by the following parameters:

• number of coarser levels;

• number of pre-and post-smoothing steps (or accuracy) for each level;

• number of subcycles;

• smoother type for each level. Typically we use all available levels, and on each level-unless stated otherwise-we use 10 Runge-Kutta pseudo-time iterations for both the pre-and the post-smoothing step.

• Partially implicit v-and w-cycles: Multigrid convergence hinges upon the convergence of the lowest frequency (order) errors which are no longer converted to higher frequencies by coarsening. Returning to the Newton-Krylov implicit solver we observe that, for zeroth order interpolation, the storage implied in the ILU preconditioner is relatively small. The discretization is positive and hence no limiter is required to stabilize the scheme near solution discontinuities. So we have very convincing motivations to use the implicit solver on the coarsest levels.

Starting the Newton-Krylov iterations from the same CFL each time we get to the coarsest level is not a good idea. The non-linear convergence depends on it, and the required stabilization diminishes as the solution is further converged. Therefore, we update CFL 0 for the Newton-Krylov scheme following each multigrid cycle. We use the same rationale as in [START_REF] Bassi | GMRES discontinuous Galerkin solution of the compressible Navier-Stokes equations[END_REF], but now using the ratios of ÿne grid residual norms. Stable values for CFL 0 , CFL max and may be determined from a zeroth order run. • Full multigrid (FMG): Starting from the coarsest representation we successively reÿne the representation. Using the hitherto available coarser levels one performs a number of v-or w-cycles. Depending on the strategy a prespeciÿed number of cycles is done, or until a desired level of convergence has been reached. After that the solution is prolongated to the next ÿner level. This successive reÿning of the order is also called nested iteration. When ÿnally reaching the target order, one continues v-or w-cycling until full convergence. This combination is coined full multigrid.

APPLICATIONS

Grid convergence study: ow around a circular cylinder

To check the grid convergence for both the full quadrature and the quadrature free implementation we consider the computation of the inviscid ow around a circular cylinder at Mach 0.3. We use this testcase since it is easy to generate high quality nested meshes, and because the surface curvature is constant, so there is no discussion about grid clustering. The meshes are the same as those used by Krivodonova et al. [START_REF] Krivodonova | High-order accurate implementation of solid wall boundary conditions in curved geometries[END_REF]. The basic grid is fully regular and symmetric, the coarser grids are sequentially nested into the ÿnest mesh. Details of the meshes near the cylinder are displayed in Figure 3. The label of the meshes refers to the number of points on the cylinder times the number of parallel mesh layers between the cylinder and the free-stream boundary.

Orders ranging from 0 to 3 have been investigated. The numerical error has been measured by the L2-norm of the following quantity related to entropy (p ∞ and ∞ correspond to the free-stream values) [START_REF] Krivodonova | High-order accurate implementation of solid wall boundary conditions in curved geometries[END_REF]:

s = p p ∞ • ∞ -1 (27)
The evolution of the error norm is shown in Figure 4. The grid convergence rate r has been computed for each reÿnement separately and globally, and summarized together with the error norm = L 2 (s ) in Tables I andII. We see that the order of accuracy of the scheme in case of full quadrature corresponds to the results obtained by Krivodonova. The order of accuracy and solution accuracy are deteriorated by the quadrature free integration, but not in a very dramatic manner. Figures 5 and6 show comparisons of the computed Mach number isolines for both the full quadrature and the patched quadrature free implementation. The Mach number distribution is visualized using 21 isolines between 0 and 0.65.

Comparison of iterative strategies-subsonic ow around the NACA0012 airfoil

To compare the iteration strategies presented in this paper we compute the ow around the NACA0012 airfoil. The angle of attack is 2 • and the Mach number is 0.3. The geometry of the airfoil may be found in Reference [START_REF] Ladson | Computer program to obtain ordinates for NACA airfoils[END_REF]. The chord has been rescaled to avoid the truncation of the airfoil at the trailing edge. The free-stream boundary is located at 30 chords from the leading edge. The mesh size on the airfoil varies quadratically, from 0.01 at the leading edge, over 0.08 at midchord to 0.02 at the trailing edge. At the free-stream boundary the grid size is 20. The coarse mesh is obtained by doubling those sizes. The details of the meshes near the airfoil are shown in Figure 7. The coarse mesh contains in total 157 nodes, the ÿne mesh 622 nodes. Again the simpliÿed treatment for boundary curvature has been used.

Computed Mach number isolines for interpolation orders 1 and 4 on the coarse and on the ÿne mesh are compared in Figure 8. In Figure 9 we compare performance of three full multigrid strategies in terms of CPU time. All strategies use fully explicit w-cycles and are compared to the standard w-cycle. Residuals are always computed using the solution projected onto the highest order representation. The ÿrst strategy, labelled 'w-fmg1' only uses one cycle per level during the nested iteration. The second one ('w-fmg2') converges four orders of magnitude for each reÿnement. The last one ('w-fmg3') converges each level to machine accuracy. We see that w-fmg1 and w-fmg2 have more or less the same performance as the w-cycle. W-fmg3 performs much worse. This is due to the fact that from a given resolution onward the low order solution generates the same higher-order error. All the time spent in improving the low order solution with respect to this critical resolution is then lost. At least for this case, we have nothing to gain from FMG. The only advantage we may reasonably expect from FMG is an increased stability. Full resolution for each successive reÿnement is however a waste of CPU time.

In Figures 10 and11 we compare the convergence rates of the di erent iteration strategies in terms of number of multigrid cycles (if applicable) and CPU time for orders 2, 3 and 4 on the coarse and the ÿne mesh, respectively. The strategies included in this study are:

• fully explicit v-cycle (labelled 'v-cycle');

• fully explicit w-cycle (labelled 'w-cycle');

• v-cycle with implicit coarsest level (labelled 'v-cycle-i0');

• w-cycle with implicit coarsest level (labelled 'w-cycle-i0');

• Runge-Kutta pseudo-timestepping (labelled 'Runge-Kutta');

• Newton-Krylov-ILU (labelled 'Newton-Krylov'). All explicit levels use 10 Runge-Kutta pseudo-time iterations for both pre-and postsmoothing sweeps. On the implicit level we specify CFL 0 = 5 and = 0:5, whilst 30 Krylov vectors are used. The Runge-Kutta pseudo-time integration and the Newton-Krylov implicit scheme both use the same parameters as their pendants on the multigrid levels.

We see that the v-cycles have the same rate of convergence per cycle irrespective of the ÿnest-level interpolation order, i.e. we attain textbook multigrid e ciency. We see this behaviour only for v-cycles since the number of passes on each level of this strategy is independent of the ÿnest-level interpolation order. The number of cycles however is very dependent on the grid size. This is to be expected, since the maximum stable timestep is proportional to the grid size. If we had used a classical multigrid method, using coarser grids instead of decreasing order, we would have convergence rates per cycle independent of mesh size, provided the coarsest mesh remains the same. A combination of p-and classical multigrid could provide both order and grid size independent convergence rates. We note that in any conÿguration (fully or partly explicit) w-cycles have a higher asymptotic convergence rate per cycle, almost in a constant ratio to the number of coarser level sweeps (v-cycle orders 2 and 3: 1 pass, w-cycle order 2: 2 passes, w-cycle order 3: 4 passes). Since the extra work is conÿned to the coarser levels we can sometimes see a better performance in terms of CPU time.

The application of the Newton-Krylov strategy as the smoother on the coarsest level tends to speed up v-cycles considerably, both in number of cycles and CPU time. For w-cycles the e ect is most marked in the early stages, where the convergence rate seems to be dominated by low frequency errors. After an initial transient however we ÿnd asymptotic convergence rates that are similar to the ones for the fully explicit cycles. A possible explanation is that w-cycles, as they concentrate most of the workload on the coarser levels, smooth out the lowfrequency error more e ectively. Consequently this strategy, especially in the partly implicit case, would be dominated by the convergence of higher frequencies only in the last stage of convergence. v-cycles would tend to keep a balanced 'broad-band' error, and maintain e ciency on all levels up to the end. This e ect is probably only noticeable due to the high accuracy of the transfer operators.

Finally we also note that in all cases the implicit Newton-Krylov scheme converges ÿrst. However as order and number of DOFs increases, the multigrid scheme becomes a viable competitor, even when we use a poor smoother such as Runge-Kutta, and this for a fraction of the memory.

CONCLUSIONS

In order to enhance the computational e ciency of DGFEM two algorithmical simpliÿcations have been implemented. The ÿrst one is a simpliÿed representation of the boundary, thus avoiding the necessity of isoparametric elements. The second one is a simple implementation of the quadrature free integration method. Both simpliÿcations, although leading apparently to a decreased formal order of accuracy, retain su cient absolute precision. The quadrature free implementation has been shown to need stabilization, and a suitable strategy has been implemented for subsonic ows. A Newton-Krylov and di erent p-multigrid iteration strategies have been presented and compared. Even using very simple smoothers, the p-multigrid proves to be a performant option in terms of CPU time whilst drastically reducing memory footprint and code complexity.

Figure 1 .

 1 Figure 1. Singular point instability: (a) quadrature free; and (b) full quadrature.

Figure 2 .

 2 Figure 2. Multigrid strategies.

Figure 3 .

 3 Figure 3. Meshes for the cylinder testcase: (a) 16 × 4; (b) 32 × 8; (c) 64 × 16; and (d) 128 × 32.

Figure 4 .

 4 Figure 4. Grid convergence-cylinder testcase.

Figure 5 .

 5 Figure 5. Mach number distributions for the ow around the circular cylinder, 16 × 4 mesh: (a) order 1-full quadrature; (b) order 1-quadrature free; (c) order 3-full quadrature; and (d) order 3-quadrature free.

Figure 6 .

 6 Figure 6. Mach number distributions for the ow around the circular cylinder, 32 × 8 mesh: (a) order 1-full quadrature; (b) order 1-quadrature free; (c) order 3-full quadrature; and (d) order 3-quadrature free.

Figure 7 .

 7 Figure 7. Meshes for the ow around the NACA0012 airfoil: (a) coarse resolution; and (b) ÿne resolution.

Figure 8 .Figure 9 .

 89 Figure 8. Mach number distribution: (a) p = 1, coarse; (b) p = 4 coarse; (c) p = 1, ÿne; and (d) p = 4, ÿne.

Figure 10 .Figure 11 .

 1011 Figure 10. Comparison of iteration strategies, coarse mesh: (a) number of cycles, order = 2; (b) CPU time, order = 2; (c) number of cycles, order = 3; (d) CPU time, order = 3; (e) number of cycles, order = 4; and (f) CPU time, order = 4.

Table I .

 I Grid convergence-full quadrature-cylinder testcase.

		p = 1		p = 2		p = 3
			r		r		r
	16 × 4	6 :315e -02	-	8:228e -03	-	1:244e -03	-
	32 × 8	1 :169e -02	2.43	5:162e -04	3.99	6:425e -05	4.27
	64 × 16	1:812e -03	2.69	4:272e -05	3.59	4:552e -06	3.82
	128 × 32	2:676e -03	2.76	4:312e -06	3.31	5:876e -07	2.95
	Global	-	2.63	-	3.63	-	3.68

Table II .

 II Grid convergence-quadrature free-cylinder testcase.

		p = 1		p = 2		p = 3
			r		r		r
	16 × 4	6 :634e -02	-	9:216e -03	-	1:707e -03	-
	32 × 8	1 :251e -02	2.41	1:204e -03	2.94	1:695e -04	3.33
	64 × 16	2:226e -03	2.49	1:924e -04	2.65	1:343e -05	3.66
	128 × 32	5:710e -04	1.96	2:978e -05	2.69	1:132e -06	3.57
	Global	-	2.29	-	2.76	-	3.52

ACKNOWLEDGEMENTS

The visualizations and most of the meshes have been generated with the open source grid generator Gmsh (http://www.geuz.org/gmsh).

Contract=grant sponsor: CENAERO is funded by the Walloon Region and the structural funds ERDF and ESF; contract=grant number: EP1A122030000102


Slawinski Js

Billat Vl

Facsm

Slawinski Jean
email: jeanslawinski@yahoo.fr

Difference in mechanical and energy cost between highly, well and non trained runners

Keywords: oxygen uptake, biomechanics, training, running

Introduction:

Recently it has been shown that endurance training decreases the variability in stride rate. This decrease would lead to a reduction in the mechanical and the energy cost of running. Purpose: This study therefore aims to compare the mechanical and the energy cost of running according to the training status of the runner (highly, well and non trained endurance runners). Methods: The kinetic, potential and internal mechanical costs (C ke , C pe and C int ) were measured with a 3D motion analysis system (ANIMAN3D). The energy cost of running (C), was measured from pulmonary gas exchange using a breath by breath portable gas analyser (Cosmed K4b 2 , Roma, Italy). All the parameters were measured on track, for a speed of 4.84 ± 0.36 m.s -1 . Results: Highly trained runners did not exhibit significantly lower C compared to well or non trained runners (4.46 ± 0.38; 4.33 ± 0.32; 4.46 ± 0.46 J.kg - 1 .m -1 , respectively; p = 0.75). However, C pe was significantly lower in highly and well trained runners compared to non trained runners (0.43 ± 0.07; 0.45 ± 0.05; 0.54 ± 0.08 J.kg -1 .m -1 , respectively; p < 0.05). In contrast, C int was significantly higher in highly trained runners compared to well and non trained runners (respectively 0.80 ± 0.12; 0.60 ± 0.09; 0.59 ± 0.10 J.kg -1 .m -1 ; p < 0.05). Conclusion: Although there is a significant difference in C pe and in C int between runners of various training status, there is no difference in C. Differences in C pe and C int may be associated with the same self-optimizing mechanism which contributes to a reduction in the impact loads during the initial portion of the support phase of the stride.

Introduction

Paragraph Number 1

In middle and long distance events, maximal performance depends on the optimisation of aerobic and biomechanical factors. The energy cost of running, defined as the amount of energy spent per unit of distance [START_REF] Di Prampero | The energy cost of human locomotion on land and in water[END_REF], reflects the sum of both aerobic and biomechanical demands. In aerobic conditions, the metabolic power can be estimated by oxygen uptake (

2 O V &
). Among long distance runners with comparable maximal oxygen uptake (

max 2 O V &
), the energy cost of running (C) is a discriminating parameter in endurance performance [START_REF] Di Prampero | The energy cost of human locomotion on land and in water[END_REF]. While various factors have been found to affect C [START_REF] Morgan | Factors affecting running economy[END_REF], little is known about training effect on C. It seems that the use of interval training sessions over a long period (14 weeks to 5 years) brings about a decrease in C [START_REF] Billat | Interval training at 2max O V & effects on aerobic performance and overtraining markers[END_REF][START_REF] Petray | Running training instruction on running technique, and running economy in 10 year-old males[END_REF][START_REF] Smith | Effects of 4-week training using Vmax/Tmax on 2max O V & and performance in athletes[END_REF]. However, many studies did not find any link between these parameters, so the reasons for C improvement with training are not well known.

Paragraph Number 2

Biomechanical parameters have been identified as factors that could influence C.

Williams and Cavanagh [START_REF] Williams | Relationship between distance running mechanics, running economy, and performance[END_REF] demonstrated that 54% of the inter-individual variability of C could be explained by kinematic variables. These authors demonstrated that the more economical runners possess a characteristic running style. They present numerous differences in kinematic variables and some of them are correlated with C.

The mechanical cost of the movements of the centre of mass has also been identified as a potential biomechanical parameter that could explain the inter-individual variability in C [START_REF] Bourdin | Effect of vertical loading on energy cost and kinematics of running in trained male subjects[END_REF][START_REF] Candau | Energy cost and running mechanics during a treadmill run to voluntary exhaustion in humans[END_REF]. Indeed, C represents the energy used by the active muscle to produce mechanical work. Thus, the inter-individual differences in C have to be related to the variation in mechanical energies associated with the movement of the different segments of the body during running. Some studies have found a significant relationship between the inter-individual differences in C and the variation in mechanical energies of the centre of mass during running [START_REF] Bourdin | Effect of vertical loading on energy cost and kinematics of running in trained male subjects[END_REF][START_REF] Candau | Energy cost and running mechanics during a treadmill run to voluntary exhaustion in humans[END_REF]. Recently, Borrani et al. [START_REF] Borrani | Is the 2 O V & slow component dependent on progressive recruitment of fast-twitch fibers in trained runners?[END_REF] suggest that the modification of C could be related to a modification of an index of the internal cost. Therefore, in order to understand why there are significant differences in C among runners of different training status, there is a need to compare kinetic, potential, and internal energies indices in addition to the overall energy cost.

However, effects of endurance training on these biomechanical factors have been poorly investigated. Indeed, most of the current studies have focused on the effect of training on simple mechanical stride parameters such as stride rate or stride length [START_REF] Bailey | Variation in stride length and running economy in male novice runners subsequent to a seven-week training program[END_REF][START_REF] Williams | Relationship between distance running mechanics, running economy, and performance[END_REF]. More recently, a study performed in our laboratory [START_REF] Petray | Running training instruction on running technique, and running economy in 10 year-old males[END_REF], showed that endurance training leads to a decrease in the stride rate variability. It has been hypothesised that the decrease of the stride rate variability would induce a decrease in the mechanical cost of running by reducing the amplitude of movement of the centre of mass.

However, direct measurements of the effect of the training on kinetics, potential and internal mechanical cost (C ke , C pe and C int ) are not available.

Paragraph Number 3

Therefore, this study aims to examine the differences of mechanical and energy cost for groups of runners with different levels of training. The hypothesis being that the lower energy cost in highly trained runners would be associated with a lower mechanical cost. In order to test this hypothesis, C, C ke , C pe and C int were measured in three populations with different training status.

METHODS

Paragraph Number 4

Subjects

Three groups took part of the experimental protocol:

-Seven highly trained (3 men and 4 women) runners (170 ± 7 cm; 55.

O V & ).
All subjects were informed of risks and stress associated with the experimental protocol and gave a written voluntary informed consent in accordance with the guidelines of the hospital of Paris S t Louis.

Paragraph Number 5 Experimental design

The energy cost of running (C), the kinetic and potential mechanical cost associated with the movements of the centre of mass (C ke and C pe ), and the internal mechanical cost associated with the movements of the body's segments around the centre of mass (C int ) were measured during a constant load running exercise (V = 4.84 ± 0.36 m.s -1 ).

The energy cost of running (C), the kinetic and potential mechanical cost associated with the movements of the centre of mass (C ke and C pe ), and the internal mechanical cost associated with the movements of the body's segments around the centre of mass (C int ) were measured during a constant load running exercise (V = 4.84 ± 0.36 m.s -1 ).

Highly trained runners carried out a 30 minute run on a level road at a constant velocity (table 1). The velocity was set at the average velocity sustained during their best marathon race. Runners followed a pacing cyclist traveling at the required velocity. The cyclist was equipped with a speedometer with a speed precision of ± 0.

Paragraph Number 6

The well trained and the non trained subjects performed two exhaustive exercises. The first exercise was an incremental test allowing

Paragraph Number 7 Energy cost measurement (C).

Throughout the constant velocity exercises, the respiratory and pulmonary gas exchange variables were measured using a breath-by-breath portable gas analyzer (Cosmed K4b 2 , Roma, Italy ( 26)). Before each test, O 2 and CO 2 analyser were calibrated using ambient air and sample gas references. The flowmeter was calibrated with a 3-l syringe (Quinton instruments, Seattle, USA). Blood lactate concentration was measured for the three groups of runners, before and after the constant velocity exercises.

Paragraph Number 8

The breath-by-breath oxygen uptake data were reduced to 5-s stationary averages.

These data were then smoothed, using a 3-step average filter, to reduce the noise, so as to enhance the underlying characteristics. These data were finally fitted to mathematical exponential models (2) using an iterative nonlinear regression on Sigma Plot software (SPSS, Chicago, IL, USA). A single-exponential model (eq. 1) and a double-exponential models were used, comprising two exponential terms which start at two distinct time delays from the onset of exercise (eq. 2). The Fisher test, which was performed by the Sigma Plot software, was used to choose the model for which the fit was associated with the highest F value (2).

O V & (t) = A 0 + A1×(1-e -(t-TD1)/ τ 1 ) ×U 1 (linear and fast component) [1] + A2×(1-e -(t-TD2)/ τ 2 ) ×U 2 (slow component) [2]
With:

A 0 : the basal metabolic rate (ml.min -1 ).

A1 and A2 : the asymptotic amplitudes for the exponential terms (ml.min -1 ).

τ 1 and τ 2 : the time constants (s).

TD1 and TD2: the time delay from the onset of exercise (s).

C was expressed in J. kg -1 .m -1 and computed from the following equation:

C (J.kg -1 .m -1 ) = ) t ( O V 2 & × 2 O E × V -1 × M -1
To calculate the energy expenditure, an energy equivalent of oxygen (

F = 0.94 -10 -3 × T [3]
Where "T" is the best performance time (min) achieved in an endurance race. The intensity of the 30 min test was set at the average velocity sustained during the best marathon race of each runner. Therefore, F can be deduced from the best time achieved during the best marathon. From F and the oxygen uptake

O V & and v 2max O V & : 2max O V & (ml.kg -1 .min -1 ) = F ) t ( O V 2 & × M -1 [4] v 2max O V & (m.min -1 ) = 1 - 1 - 2 2max M V (t) O V O V × × & & [5]
Where "V" is the average velocity of the 30 min test, (expressed in m.min -1 ) and "M" is the subject's mass.

Paragraph Number 10

Mechanical cost measurement

During the constant load exercise, runners were filmed by a digital video camera (Sony, Japan, TRV 900) using a sampling frequency of 25 Hz. Two consecutive steps were analysed at the end of the exercise (figure 1). The video sequences were digitised on a computer using a video card (Perception Video Recorder, Digital System Inc.), and then transformed into a set of bitmap pictures. The parametric trajectories of the anatomical points of the runner were determined using a motion analysis system, ANIMAN3D [START_REF] Tavernier | A 3D motion analysis system using a numerical human model[END_REF]. This system uses a numerical manikin "MAN3D" (figure 2). The morphology of this numerical manikin integrates inertial properties of the different segments. The morphology of this numerical manikin could be adjusted to each runner from their height and weight. This model measures the 3D positions of the body, during running. In the present work, the lateral displacements of the different centres of masses have not been directly measured because only one camera was used. All the trajectories of the different centres of mass (segments and body centre of mass) were smoothed using a polynomial method in order to obtain by derivation the speed displacement of the different centres of mass. Variations of potential and kinetic energies (ΔE pe and ΔE ke respectively) were then calculated (in joules).

ΔE pe = ( ) min max H H g M - × × [6] ΔE ke = ( ) 2 min 2 max V V M 2 1 - × [7]
"M" is the body mass (kg), "g" is the gravitational acceleration (9.81 m.s -2 ), "H max and H min " are the maximal and minimal heights of the body centre of mass (CM) during one step (m). "V max and V min " are the maximal and minimal horizontal velocities of the CM during one step (m.s -1 ).

Paragraph Number 11

The variation of the internal energy (ΔE int ) represents the variations in mechanical energies associated with the movement of the different segments of the body around the CM. The present model used in this work (figure 2) is composed of 17 different segments. The variation of internal energy (Δ int E ) is calculated,

( ) ∑ = + = 17 1 i 2 s s 2 s s int ω I V m 2 1 Δ ΔE [8]
Where "s" represents a segment, "m s " is the mass of the segment (kg) considered and "I s " is the moment of inertia of the segment (kg.m 2 ). " s V " is the velocity of the centre of mass of the segment (m.s -1 ) with respect to the referential linked to the CM. " s ω "

is the angular velocities of the centre of mass of the segment (rad.s -1 ) with respect to the referential linked to the CM of the segment. Kinetic, potential and internal mechanical cost expressed in J.kg -1 .m -1 (C ke , C pe , C int ) are equal to the positive variation of ΔE ke , ΔE pe and ΔE int divided by the mass of the subject and his step amplitude. These parameters were measured during the steady state of

RESULTS

Paragraph Number 10

Characteristics of the highly, well and non trained runners

2max O V & and v 2max O V &
were significantly greater in trained subjects. Indeed, both parameters were greater for highly trained runners compared to well trained runners or compared to non-trained runners (table 1).

The speed of the constant load exercise was not significantly different between the groups (p=0.33). The different measures obtained from the different populations, during the constant load exercise (C pe , C ke , C int and C), were performed for similar running speeds (table 1).

Paragraph Number 13

Effect of the endurance training status on C, C pe , C ke , C int and the step rate.

The metabolic cost of running reported in table 2 is the "gross" value, including the O 2 consumption at rest. The less trained runners do not consume more energy than the highly trained runners. Indeed, in order to move forward at the same running speed (table 2) the highly trained runners present a same C as well-trained, or non-trained runners. Moreover, when the energetic equivalent of lactate (EEL) is added to C there is also no difference between the groups.

Paragraph Number 14

C pe was significantly lower in highly and well trained runners compared to non trained runners (20%). Inversely, the step rate (SR) was significantly higher in highly trained runners compared to well and non trained runners. Moreover, there was an inverse correlation between SR and C pe or the vertical displacement of the centre of mass (CM)(respectively r = -0.63; p = 0.001 and r = -0.81; p ≤ 0.0001). Figure 3 shows that runners who present the greatest vertical displacement of the CM also possess the weakest limb displacement speed around the CM. Indeed, the vertical displacement of the CM for the different group is 6.7 ± 0.9 cm for the highly trained runners compared to 7.7 ± 0.7 cm for the well trained runners, and 8.3 ± 1.3 cm for the non trained runners (p = 0.02).

Paragraph Number 15

C int was also significantly higher in highly trained runners: 30% greater in comparison with well trained runners and 37% greater in comparison with non trained runners. C int was correlated with the vertical displacement of the CM (r = -0.46; p = 0.03), but not with SR (r = 0.41; p = 0.07) or with C pe (r = -0.32; p = 0.16).

Paragraph Number 16

The running speed was similar in all the runners, therefore C ke was not significantly different in the three populations of different training status (p = 0.89).

Paragraph Number 17

Relationship between C, C ke , C pe , C int , SR and

DISCUSSION

Paragraph Number 18

The purpose of this work was to test the hypothesis that the lower energy cost in highly trained runners is associated with a lower mechanical cost. To test this hypothesis, C, C ke , C pe and C int were measured in runners who differ in training status. The results showed that contrary to our hypothesis, the energy cost of running was not different for runners with different training status. However, the mechanical cost is different for runners with different training status. Indeed, the mechanical cost associated to the vertical movements of the centre of mass (C pe ) was smaller in highly trained runners. Inversely, the mechanical cost associated with the movement of the segments around the centre of mass (C int ) and the step rate (SR) were greater.

Paragraph Number 19

There was no difference in the energy cost of running between the three groups, whereas such a difference would be expected and has been demonstrated in previous studies [START_REF] Billat | Interval training at 2max O V & effects on aerobic performance and overtraining markers[END_REF][START_REF] Petray | Running training instruction on running technique, and running economy in 10 year-old males[END_REF][START_REF] Smith | Effects of 4-week training using Vmax/Tmax on 2max O V & and performance in athletes[END_REF]. Numerous mechanisms such as a decrease in the ventilation, an optimisation of the neuromuscular recruitment, an improvement in the transport of oxygen by the blood, the use of this oxygen by the active muscles or the increase of the muscular elasticity [START_REF] Daniels | Physiological characteristics of champion male athletes[END_REF][START_REF] Paavolainen | Explosive strength training improves 5-km running time by improving running economy and muscle power[END_REF]) could be at the origin of the improvement in the energy cost with training. However, mechanisms associated with this decrease are not well known and numerous studies do not report any change after a period of training [START_REF] Daniels | Changes in VO2 max and running performance with training[END_REF][START_REF] Lake | Six weeks of training does not change running mechanics or improve running economy[END_REF][START_REF] Petray | Running training instruction on running technique, and running economy in 10 year-old males[END_REF][START_REF] Slawinski | Effect of supra-lactate threshold training on the relationship between mechanical stride descriptors and aerobic energy cost in trained runners[END_REF][START_REF] Wilcox | Change in running economy relative to 2max O V & during a cross-country season[END_REF]. Daniels and Daniels [START_REF] Daniels | Physiological characteristics of champion male athletes[END_REF] have shown that for intensities typically used in races up through the marathon, no differences exist in C within the range of most intensities (80-100% of max 2

O V &

). In the present study, the range of intensities was comprised within 80% and 97% of

v max 2 O V &
. Therefore, the lack of difference in C among runners of different training levels may be associated with the fact that C is measured for speeds upper to speeds typically used during a marathon race.

Paragraph Number 20

Following this hypothesis, the present results show that the mechanical cost of the vertical movements of the CM (C pe ) is smaller in highly trained runners. This result is in accordance with the generally acknowledged principle which states that smaller verticals oscillations of the body centre of mass are associated with a high level of training [START_REF] Morgan | Factors affecting running economy[END_REF][START_REF] Williams | Relationship between distance running mechanics, running economy, and performance[END_REF]. However, there was no relationship observed between C pe and C (figure 4). According to the training status, the variation in C pe is not associated with the variation in C. Similarly, the kinetic mechanical cost (C ke ) or the mechanical cost of the movements of the segments around the centre of mass (C int ) were not correlated with C (figure 4). These results agree with numerous studies which have demonstrated that the relationships between C and the mechanical cost are both weak and inconsistent [START_REF] Dalleau | The spring-mass model and the energy cost of treadmill running[END_REF][START_REF] Morgan | Factors affecting running economy[END_REF][START_REF] Williams | Relationship between distance running mechanics, running economy, and performance[END_REF]. The mechanical cost reflects a global and indirect expression of the muscular effort which explains why a low correlation is observed between C and the mechanical cost. Other mechanical parameters should be correlated with C.

Indeed, it has been shown that less economical runners possess a more compliant running style during ground contact [START_REF] Dalleau | The spring-mass model and the energy cost of treadmill running[END_REF][START_REF] Heise | Leg spring" characteristics and the aerobic demand of running[END_REF].

Paragraph Number 21

Collectively, the results reported here, do not confirm our primary hypothesis and suggest that the mechanical cost (C pe and C int ) is not correlate with C but with 2max O V & . The mechanism that explains the modifications of the mechanical parameters (SR, C pe and C int ) is in relation with the level of training of the runners rather than their energy cost. However, this mechanism has to be elucidated.

Paragraph Number 22

If humans are self-optimizing machines, the minimal cost, being an optimally criterion, may be identified that governs the kinematic and kinetic detail of the performance. For example, it has been shown that adult and children tend to walk or run at frequencies that are determined by the oxygen cost [START_REF] Cavanagh | The effect of stride length variation on oxygen uptake during distance running[END_REF]. However, it is unlikely that the metabolic cost is the only optimality criteria adopted for human activity. The potential for injury may also result in the development of optimal criteria. Farley and Taylor [START_REF] Farley | A mechanical trigger for the trot-gallop transition in horse[END_REF] reported that horses naturally switch from a trot to a gallop actually increasing their metabolic cost but reducing the peak of forces on the muscles, tendons and bones. They suggest that this mechanism reduces the chance of injury for the horse. During human running, the body is subjected to high impact loads during the initial portion of the support phase of the stride. These impacts have magnitudes up to 2.3 times body weight (BW) with an impact load rate of 113 BW.s -1 [START_REF] Munro | Ground reaction forces in running: a reexamination[END_REF]. Increases in impact shock can result from an increase in running speed [START_REF] Hamill | Variation in ground reaction force parameters at different running speeds[END_REF], from running downhill [START_REF] Hamill | Effect of grade running on kinematics and impact force[END_REF], from an increase in stride length [START_REF] Derrick | Energy absorption of impacts during running at various stride lengths[END_REF][START_REF] Hamill | Shock attenuation and stride frequency during running[END_REF] or fatigue [START_REF] Verbitsky | Shock transmission and fatigue in human running[END_REF]. It would appear therefore, that impact shock attenuation may be an important factor on which individuals optimise. In a same way, Zamparo et al. [START_REF] Zamparo | The self selected speed of running in recreational long distance runners[END_REF] have recently demonstrated that the self selected speed of running depends not only on cardiovascular factors, but also on biomechanical factors.

Paragraph Number 23

The running style of the highly trained runners may be associated with the same selfoptimizing mechanism which contributes to reducing the impact loads during the initial portion of the support phase of the stride. Indeed, this style is characterised by a high SR and a low C pe and a low vertical displacement of the CM. Moreover, the figure 3 shows that the vertical displacement of the CM is inversely correlated to SR.

This inverse relationship shows that a high stride rate is linked to a smaller vertical displacement of the CM. During the flight phase, the CM of highly trained runners may reach a lower height (table 2), which decreases the impact shock when the foot hits the ground. This result is in accordance with the results of Farley and Gonzales [START_REF] Farley | Leg stiffness and stride frequency in human running[END_REF] which show that at higher step rate, the initial vertical force peak is absent and that the vertical force peak is decreased. This style of running may be associated with 6)) compared to the other populations studied. For example, running 32 km.week -1 will produce over 1.3 million impacts to the body over a period of one year [START_REF] Derrick | Impacts and kinematic adjustments during an exhaustive run[END_REF]. More frequent impacts will place greater stress on the muscles, tendons and bones resulting in an increased risk of injury and degenerative disease [START_REF] Derrick | Energy absorption of impacts during running at various stride lengths[END_REF]. Therefore, to attenuate this risk, mechanical adjustments such as a decrease in the vertical oscillation of the CM may take place with training. Following this hypothesis, the modification of the mechanics of running with the training status does not have any relationship with the energy cost of running (C).

Paragraph Number 24

To conclude, highly trained runners did not display a lower C. However, C pe was 20% lower and C int was 30% greater in highly trained runners. The endurance training leads to an increase of C int and a decrease in C pe without any modification of C. These mechanical adjustments may be associated with the same self-optimizing mechanisms which contribute to reduce the impact loads during the initial portion of the support phase of the stride.

The authors gratefully acknowledge the Fédération Française de Ski and Michel

Tavernier for their assistance during the experiment.

  and the velocity at the lactate threshold (vLT) to be determined. The running speed of this test was progressively increased from 12 km.h -1 to exhaustion. The velocity increments between the stages (3 min duration) were set at 1 km.h -1 . All stages were followed by a 30-s period of rest. During this period, a fingertip capillary blood sample was collected in order to calculate the energetic equivalent of lactate (EEL). The second exercise was a test to exhaustion at the velocity corresponding to 50% of the work rate difference between vLT and v max 2 O V & (vΔ50) for the well trained runners. For the non trained runner, the velocity of the exercise was set at 95% of v 2max O V & . C, C ke , C pe and C int were determined during the last minute of the exercise. All test sessions were completed on a 400-m covered synthetic track. Throughout the tests, the subjects adopted the required velocity using an audio rhythm which gave the time to cover 25m. Visual marks were set at 25-m intervals along the track (inside the first lane).

2

 2

  during the 30 min test, it was easy to calculate 2max

Fig. 1 :

 1 Fig. 1: Representation of space calibration and of the field of the camera

Fig. 3 :

 3 Fig. 3 : Correlation between the vertical displacement of the centre of mass (VD of the CM) and the step rate (SR). The verticals bars represent the standard error of measure.

  and C int were significantly different according to the training status of the runners, there was no correlation between C and C pe or C int (figure4). The correlation coefficients were near r = 0.05 and p = 0.83 for C ke , r = -0.13 and p = 0.58 for C pe and r = 0.41, and p = 0.06 for C int . Furthermore, there was no correlation between the SR and C (r = 0.15; p = 0.53). These results showed that the interindividual variability of C can not be explained by biomechanical parameters such as SR, C pe , C ke or C int . However, SR, C pe and C int were significantly correlated with 2max O V & (respectively r = 0.58; p = 0.007; r = -0.73; p = 0.0001; r = 0.60; p = 0.004).

Fig. 4 :

 4 Fig. 4: Relationship between the energy cost of running (C), the kinetic mechanical cost (C ke ), the potential mechanical cost (C pe ) and the internal mechanical cost (C int ).

  the training mode of the runners. Indeed, the nature of the adaptive responses to training is specific to the training stimulus. Highly trained runners cover very long distances during training session, and perform a very high number of strides (more than 20 km.day -1 ; (

Table 2 :

 2 Effects of the endurance training status on the gross energetic cost (C), the gross energetic cost plus the energetic equivalent of lactate (C + EEL), the step rate (SR), the step length (SL), the kinetic, potential and internal mechanical costs (C ke , C pe , C int , respectively).

		Highly trained (±SD)	Well trained (±SD)	Non trained (±SD)	ANOVA
	C ke (J.kg -1 .m -1 )	1.54±0.70	1.41±0.40	1.43 ±0.33	0.89
	C pe (J.kg -1 .m -1 )	0.43±0.07	0.45±0.05	0.54 ±0.08 † §	0.02
	C int (J.kg -1 .m -1 ) 0.80±0.12	0.60±0.09 †	0.59 ±0.10 †	≤0.01
	C (J.kg -1 .m -1 )	4.46±0.38	4.33±0.32	4.46 ±0.46	0.75
	C + EEL (J.kg -1 .m -1 )	4.56±0.39	4.75±0.32	4.89 ±0.44	0.30
	R	0.94±0.08	1.00±0.03	1.06 ±0.05 †	<0.01
	SR (Hz)	3.21±0.22	3.01±0.11 †	2.96 ±0.21 †	0.05
	SL (m)	1.52±0.14	1.66±0.08 †	1.58 ±0.12	0.09
	† Significantly different from highly trained runner (p≤0.05).

§

Significantly different from well trained runner (p≤0.05).


Eric Lao
email: eric.lao@lip6.fr

Marie-Minerve Louërat
email: marie-minerve.louerat@lip6.fr

Jean-Paul Chaput

Highly configurable place and route for analog and mixed-signal circuits

Analog design remains a manual task because of the complexity of the interations among devices. Automation tools dedicated to analog circuits are not as mature as the digital automation tool but have been improved a lot, at a point that they can help at individual steps in the analog design flow. This Ph. D. thesis is oriented toward creating a semi-automated mixed design flow controlled by the designer with a particular focus on an analog and mixed signal placer and router. The results show the ability of our tool at generating multiple placed and routed layouts respecting designer's constraints.

I. CONTEXT AND MOTIVATIONS

An automatic place and route tool should produce analog device-level layouts similar in density and performance to the high-quality manual layouts. To achieve this task, the capability to deal with layout constraints, in order to decrease unwanted parasitics due to the process variations, is mandatory [START_REF] Helmut | Analog Layout Synthesis -A Survey of Topological Approaches[END_REF].

Recent research focused on using simulated annealing algorithm [START_REF] Kirkpatrick | Optimization by Simulated Annealing[END_REF] in combination with topological representations to respect these analog constraints. Topological representations encode the positioning relations between devices and the optimizer alters their relative positions. The most popular representations are Sequence-Pair [START_REF] Ma | Simultaneous Handling of Symmetry, Common Centroid, and General Placement Constraints[END_REF], B*-Tree [START_REF] Lourenc ¸o | Floorplanaware analog IC sizing and optimization based on topological constraints[END_REF], Transitive Closure Graph [START_REF] He | Symmetry-aware TCG-based placement design under complex multi-group constraints for analog circuit layouts[END_REF], Ordered-Tree [START_REF] Xiao | Analog placement with common centroid and 1-D symmetry constraints[END_REF] and slicing floorplans [START_REF] Unutulmaz | On the convex formulation of area for slicing floorplans[END_REF] and they are coupled with some constraints to respect at the same time.

Although most of the recent works focus on highly automated approaches based on optimization processes [START_REF] -H. Lin | Recent Research and Development and New and Challenges and in Analog and Layout Synthesis[END_REF], we believe that giving more control to designers and using their interventions to set some constraints is mandatory to get resulting layouts accepted by designers. Our semi-automatic approach also helps designers to debug more efficiently and makes adjustments easier while some tiresome and error-prone tasks are automated.

II. ANALOG PLACEMENT

Digital and analog circuits have a dedicated area on a system-on-chip circuit so they can be independently designed within a specific space (fig. 1). It is also common to design analog circuit in rows of devices where the height of each row of devices should be adjustable so it can match its dedicated area. In order to organize devices in row, we use the slicing tree structure which will be specified by designers. Having designers defining the slicing tree, it means they will have control over the global topological To set rows of devices, the height of each device needs to be similar and to do so, designers have to consider possible aspect ratios for each device by varying their number of fingers. In a bottom-up propagation manner, designers will have choices of multiple placements' results ordered according to their total aspect ratio. Designers also have control over margin tolerances which can eliminate some placements if devices' height in a row are too different. Once designers define their slicing tree and the previously mentioned constraints, our placer will process the accepted placements in a few seconds and present them on an interactive graph (fig. 2). Each point of the graph represents a valid placement and its coordinates are its width and height. A point can be selected and our placer will generate, in less than a second, the according placement. Designers can pick the most appropriate placement based on his experiences and the aspect ratio of the floorplan.

III. ANALOG AND MIXED SIGNAL ROUTING

Once the placement phase is performed with the slicing tree, it is then used for the routing phase. The purpose of this work is to place and route mixed signal circuits, this is the reason why our router has been built through a common design flow for routing digital and analog parts. It is divided into two steps: a global routing step and a detailed routing step.

1) Global Routing:

To iniate the global routing step, the slicing tree used for the placement is going to be converted into a graph of relation between rectangular areas (fig. 3&4). Each area represents either a device or a channel where wire will be routed. The goal is to establish by which areas each wire is going to pass by. The search of the shortest paths is performed using the Dijkstra algorithm because it is easy to configure constraints. Designers can set restrictions concerning devices which must avoid routing over themselves. Designers can also set routing constraints such as symmetry, wirelength and topology constraints for some specific signal paths that require special attention. Based on the routing estimation, routing channel are then enlarged in order to leave enough space for wires to be placed by the detailed routing phase.

2) Detailed Routing: The goal of detailed routing is to assign route segments of signal nets to specific routing tracks, vias, and metal layers with given global routes of those nets. The detailed router works as a channel router manner. Each cut of the slicing tree is considered as a routing channel that can be expanded based on a routing estimation leaving enough space for wires. Constraints specified by designers need to be respected in this step. Routing channel is performed by respecting symmetric nets, topology-matching constraints and wirelength constraints.

IV. RESULTS

Our tool was implemented in C++ programming language on a Intel(R) Core(TM) i5-4590S CPU @ 3.00 GHz workstation with 6 GB RAM. To illustrate the capability of our tool, we experiment it on a fully differential transconductor [START_REF] Iskander | Hierarchical sizing and biasing of analog firm intellectual properties[END_REF], designed in 130nm CMOS technology process (fig. 5).

The fully differential transconductor is composed of a total of 32 devices and we consider 2 possible variations for each device. The slicing tree takes into account 11 symmetries Fig. 5. Placed and routed fully differential transconductor [START_REF] Iskander | Hierarchical sizing and biasing of analog firm intellectual properties[END_REF] in 130nm CMOS technology process for this circuit and tolerance margins are set in a way to have reasonnable amount of accepted possibilities. Designers can choose their final placement based on their experiences and preferences. After choosing a placement, we perform the routing step that is performed in less than 5 seconds to obtain the results on fig. 5. In case the final result does not satisfy designers, they can try another placement to be routed easily and in matters of seconds.

V. CONCLUSION This Ph. D. thesis is about creating a highly configurable analog and mixed signal placer and router. Placement is realized using the slicing tree as floorplan representation to organize circuits in rows of devices for regularity. Designers choose their desire placement among multiple placement results that can be selected and generated in a few seconds. Their choice is based on their experience and the aspect ratio of the placement. The mixed signal router is able to handle digital and analog router through a common software structure. Analog circuits are routed with particular constraints such as symmetries, topology and wirelength constraints.

Fig. 1 .

 1 Fig. 1. Mixed signal circuit floorplanning and analog placement in 3 rows

Fig. 2 .

 2 Fig. 2. Circuits possible dimensions graph

Fig. 3 .Fig. 4 .

 34 Fig. 3. Global routing structure used for the Dijkstra algorithm