Simon Klein

Mathieu Gibert
email: mathieu.gibert@grenoble.cnrs.fr

Antoine Bérut

Eberhard Bodenschatz

Simultaneous 3D measurement of the translation and rotation of finite size particles and the flow field in a fully developed turbulent water flow

Keywords: numbers: 47.27.Gs, 47.27.Jv, 47.80.Cb 3D Full Inertial particles, fully developed turbulence, Lagrangian Particle Tracking

We report a novel experimental technique that measures simultaneously in three dimensions the trajectories, the translation, and the rotation of finite size inertial particles together with the turbulent flow. The flow field is analyzed by tracking the temporal evolution of small fluorescent tracer particles. The inertial particles consist of a super-absorbent polymer that renders them index and density matched with water and thus invisible. The particles are marked by inserting at various locations tracer particles into the polymer. Translation and rotation, as well as the flow field around the particle are recovered dynamically from the analysis of the marker and tracer particle trajectories. We apply this technique to study the dynamics of inertial particles much larger in size (R p /η ≈ 100) than the Kolmogorov length scale η in a von Kármán swirling water flow (R λ ≈ 400). We show, using the mixed (particle/fluid) Eulerian second order velocity structure function, that the interaction zone between the particle and the flow develops in a spherical shell of width 2R p around the particle of radius R p . This we interpret as an indication of a wake induced by the particle. This measurement technique has many additional advantages that will make it useful to address other problems such as particle collisions, dynamics of non-spherical solid objects, or even of wet granular matter.

Introduction

Particle-laden flows are prevalent in natural and technological flows. For example they are relevant to warm cloud dynamics [START_REF] Falkovich | Acceleration of rain initiation by cloud turbulence[END_REF][START_REF] Shaw | Particle-turbulence interactions in atmospheric clouds[END_REF] with impact on the climate, or to dust and pollution transport in atmosphere and oceans, or can be found in most technological processes where matter is mixed in fluids. In the past few years, significant advances have been made in experimental approaches to particle dynamics in turbulent flows thanks to the rapid development of visualization based three dimensional measurement techniques such as Lagrangian Particle Tracking (LPT) [START_REF] Ouellette | A quantitative study of three-dimensional lagrangian particle tracking algorithms[END_REF][START_REF] Xu | Tracking lagrangian trajectories in position-velocity space[END_REF]. When studying the fluid's dynamics with LPT, one has to assure that the particles behave passively and do not perturb the flow. This goal can be considered reached when the particles are density matched to the fluid, which avoids buoyant and inertial forces, and when their sizes are much smaller than the smallest length scale of the velocity gradients in the flow. Such particles are then called tracer particles. These requirements are very stringent and can be hard to achieve in some flows, e.g. thermal turbulent convection, where fluid density changes strongly in the thermal boundary layer [START_REF] Sun | Experimental studies of the viscous boundary layer properties in turbulent Rayleigh-Bénard convection[END_REF][START_REF] Gibert | Heat convection in a vertical channel: Plumes versus turbulent diffusion[END_REF][START_REF] Gasteuil | Lagrangian temperature, velocity and local heat flux measurement in Rayleigh-Benard convection[END_REF]; highly turbulent flows, which exhibit very fast dynamics (fractions of milliseconds) on length scales that reach micrometers [START_REF] Voth | Measurement of particle accelerations in fully developed turbulence[END_REF] (indeed L/η ∝ Re 3/4 , where L is the large scale of the flow, η is the characteristic length-scale of the smallest eddy in the flow, also called the Kolmogorov length scale, and Re is the Reynolds number of the flow); or in the study of superfluid turbulence where the normal and superfluid coexist [START_REF] Bewley | Particles for tracing turbulent liquid helium[END_REF]. When the density of the particles is different from that of the fluid, or the particle size is large, the particles cannot be regarded as tracers of the fluid and are called inertial particles.

In this paper, we focus on the study of large particles in turbulent flows, i.e. particles much larger in size than the Kolmogorov length scale η of the flow. The particles that we consider here have a radius R p ≈ 100η ≈ 0.1L. Experiments on two-particle statistics of small and heavy particles can be found in [START_REF] Gibert | Inertial effects on two-particle separation in a turbulent flows[END_REF][START_REF] Gibert | Where do small, weakly inertial particles go in a turbulent flow?[END_REF]. Most theoretical and numerical simulations of the complex coupling between the motion of inertial particles and their carrier fluid rely on simplified equations such as the Maxey-Riley-Gatignol equation [START_REF] Maxey | Equation of motion for a small rigid sphere in a nonuniform flow[END_REF][START_REF] Gatignol | The faxen formulas for a rigid particle in an unsteady non-uniform stokes-flow[END_REF]. This equation is only valid for point-like particles and couples the motion of the particle to the fluid, without any back-reaction of the particle motion on the flow (one-way coupling). This approach has been proven to be sufficient to describe the rich dynamics of small and heavy particles [START_REF] Toschi | Lagrangian properties of particles in turbulence[END_REF], but has also found its limit [START_REF] Gibert | Inertial effects on two-particle separation in a turbulent flows[END_REF][START_REF] Qureshi | Turbulent transport of material particles: An experimental study of finite size effects[END_REF]16,[START_REF] Volk | Laser doppler measurement of inertial particle and bubble accelerations in turbulence[END_REF][START_REF] Calzavarini | Acceleration statistics of finite-sized particles in turbulent flow: the role of Faxen forces[END_REF][START_REF] Mercado | Lagrangian statistics of light particles in turbulence[END_REF] (even when corrected by the so called Faxén corrections) when considering particles larger than 4η. To go beyond this simplified model, Hohmann and Bec [START_REF] Homann | Finite-size effects in the dynamics of neutrally buoyant particles in turbulent flow[END_REF] developed a direct numerical simulation using a dynamical pseudopenalisation technique in order to satisfy the non-slip boundary conditions at the surface of a unique particle evolving in a turbulent flow. This strategy is very promising (as the one developed in [START_REF] Naso | The interaction between a solid particle and a turbulent flow[END_REF], in which the particle is static) but also comes with high computational cost that currently restricts it to very low Reynolds numbers (R max λ = 72 for [START_REF] Homann | Finite-size effects in the dynamics of neutrally buoyant particles in turbulent flow[END_REF] and R max λ = 20 for [START_REF] Naso | The interaction between a solid particle and a turbulent flow[END_REF], where R λ is the Reynolds number based on the Taylor micro-scale of the flow). Therefore, for the foreseeable future, the investigations on large particle dynamics in highly turbulent flows must rely heavily on experimental work. Few experiments such as [START_REF] N T Ouellette | Transport of Finite-Sized Particles in Chaotic Flow[END_REF] have been conducted in 2D flows, also 2D studies of 3D flows based on Particle Imaging Velocimetry can be found in [START_REF] Poelma | Particle-fluid interactions in grid-generated turbulence[END_REF], but a complete experimental 3D approach resolving the particles dynamics and the flow has to our knowledge not been reported. The recent work of Zimmermann et al. [START_REF] Zimmermann | Rotational Intermittency and Turbulence Induced Lift Experienced by Large Particles in a Turbulent Flow[END_REF][START_REF] Zimmermann | Tracking the dynamics of translation and absolute orientation of a sphere in a turbulent flow[END_REF] has been able to resolve the full motion (translation and rotation) of single finite-sized particles in a highly turbulent flow. Here, we present a technique that goes beyond and resolves simultaneously the full three dimensional motion of the particles together with the turbulent flow field around it.

Experimental setup

In order to generate a highly turbulent flow in a relatively small experiment, we used a von Kármán water flow, sometimes nick-named the "french washing machine" in the literature (see Fig. 1) [START_REF] Voth | Measurement of particle accelerations in fully developed turbulence[END_REF][START_REF] Marié | Experimental measurement of the scale-by-scale momentum transport budget in a turbulent shear flow[END_REF][START_REF] Cortet | Normalized kinetic energy as a hydrodynamical global quantity for inhomogeneous anisotropic turbulence[END_REF]. This flow is generated by two counter-rotating propellers of 28 cm in diameter facing each other. The propellers are driven at the same constant rotation rate. The turbulence chamber, shaped as a octagonal cylinder, measures 40 cm along the axis of the propellers and 38 cm in both height (vertical) and width (horizontal) in the cross-section. This apparatus is particularly well suited for LPT studies [START_REF] Voth | Measurement of particle accelerations in fully developed turbulence[END_REF] since at the geometric center the mean flow velocity is zero. Therefore, at the center the particle dynamics is mainly driven by the turbulent velocity fluctuations (the large

R λ =374 ± 8 u ′ ǫ L η τ η N f (m/s) (m 2 /s 3 ) (mm) (µm) (ms) (fr/τ η )
0.1 0.011 ± 4.10 -4 93 ± 4 98 ± 1 9.5 ± 0.2 27

Table 1. Parameters of the experiment. u ′ is the root-mean-square of the velocity. ǫ is the energy dissipation rate per unit mass. L ≡ u ′3 /ǫ is the integral length scale. η ≡ (ν 3 /ǫ) 1/4 and τ η ≡ (ν/ǫ) scale flow being deterministic ‡) [START_REF] Marié | Experimental measurement of the scale-by-scale momentum transport budget in a turbulent shear flow[END_REF][START_REF] Cortet | Normalized kinetic energy as a hydrodynamical global quantity for inhomogeneous anisotropic turbulence[END_REF].

To study the turbulent flow field, we seeded the water flow with red fluorescent polymer micro-spheres with a diameter of 107 µm and a density of ρ tr = 1.05 g/cm 3 §. These particles behave as tracer particles since in our flow the Kolmogorov length scale is η ≈ 100 µm (see Tab. 1) and the fluid density is ρ f = 1 g/cm 3 . We measured three-dimensional particle trajectories with high spatial and temporal resolutions using LPT [START_REF] Ouellette | A quantitative study of three-dimensional lagrangian particle tracking algorithms[END_REF][START_REF] Xu | Tracking lagrangian trajectories in position-velocity space[END_REF] with three high-speed CMOS cameras (Phantom V10, manufactured by Vision Research Inc., Wayne, USA). The particle velocities and accelerations were then obtained by smoothing and differentiating the trajectories [START_REF] Mordant | Experimental lagrangian acceleration probability density function measurement[END_REF]. The measurement volume of the LPT was determined as the largest sphere that fits inside the complex volume defined by the intersection of the fields of view from all three cameras and the expanded laser beam. In this experiment, the sphere was 7 cm in diameter slightly smaller than the integral length scale L = (9.3 ± 0.4)cm. The frame rate was set to 2.9 kHz with a resolution of 768 × 768 (each pixel corresponds roughly to 100µm in space). Given the massive amount of data required to conduct this experiment, we used the weighted averaging algorithm (simple and efficient) to locate the particles onto the 2D images. The number of particles tracked per frame was of the order of 200. In this configuration, as demonstrated in [START_REF] Ouellette | A quantitative study of three-dimensional lagrangian particle tracking algorithms[END_REF], the 2D particle finding algorithm and 3D stereoscopic reconstruction processes ensure that the particle positions are recovered with an accuracy of 0.1 pixels ≈ 10 µm. The relatively low seeding density further guaranteed that, using a "3 frames minimum acceleration" scheme [START_REF] Ouellette | A quantitative study of three-dimensional lagrangian particle tracking algorithms[END_REF], we obtained long and high quality particle tracks.

We measured u ′ the velocity fluctuation and determined the energy dissipation ‡ The ratio between the average velocity and the velocity fluctuations τ = u i / (u i -u i ) 2 does not exceed 25% over the entire measurement volume. § These fluorescent particles are hard-dyed (internally-dyed) polymer particles which utilise the Firefli process to incorporate the dye throughout the polymer matrix. This method produces bright fluorescent colors, minimises photo-bleaching, and prevents dye leaching into aqueous media. They are made of polystyrene and are sold by Thermo scientific. Their absorption (emission) wavelength peaks around 542 nm (612 nm). rate per unit mass ǫ from the inertial range scaling of the second-order longitudinal and transverse Eulerian velocity structure functions from tracking tracer particles. As explained in details in [START_REF] Gibert | Inertial effects on two-particle separation in a turbulent flows[END_REF], we also checked the consistency of the value using two exact inertial range relations: the Kolmogorov's "four-fifth law" and a theorem on the velocityacceleration mixed structure function [START_REF] Mann | Experimental study of relative, turbulent diffusion[END_REF][START_REF] Pumir | The lagrangian view of energy transfer in turbulent flow[END_REF][START_REF] Hill | Opportunities for use of exact statistical equations[END_REF]. The properties of the turbulent flow field are summarized in Table 1.

The aim of the investigation was to simultaneously measure the 3D trajectories of the tracer particles and the full dynamics and trajectories of finite size, inertial particles. The most straightforward way was for the particles to be invisible in the fluid and to mark them with tracer particles. This avoided optical distortions and shadowing. The spherical particles were made of a super-absorbent polymer (poly-acrylate) from Aqualinos, which in their dry form had a diameter between 1 and 2 mm; once immersed in water their average diameter grew to 1 cm ≈ L/10 ≈ 100η (see Fig. 2 (a)). Thus the gel particles were 99.9% water and had almost the same refractive index and density as water. As shown in Fig. 2 (b) they are indeed almost invisible in water and do not refract the light path. In order to be able to observe the big particle's motion with the cameras, we grafted their surface with the same fluorescent particles that we used to track the turbulent flow field. As we will see later, we needed at least 4 of these markers per inertial particle in order to recover the particle's center and radius. Therefore, we injected them with 6 to 8 markers, as shown by the red dots in Fig. 2 (a) . To achieve this we employed a lancet needle normally used by diabetics. We dipped it into the tracer particle powder to collect some of them on its tip and penetrated the surface of the big particle. After removing the lancet, tracer particles were left behind within the polymer gel. We quantitatively

The regular separation presented on Fig. 2 (a) is not necessary, therefore we do not enforce it.

checked that the markers were indeed grafted to the big particles by injecting some of them exactly as depicted in figure 2 (a) thanks to a specifically designed injection guide. No relative nor absolute displacement could be detected from this well defined injection pattern. We marked about 150 big particles, which we inserted into the turbulent flow. This particle number was necessary to achieve that statistically, at any given time, at least one (and no more than two) big particles were in the measurement volume.

The big solid gel-particles were elastically deformable, but the forces (shear and pressure) applied by the fluid onto them were not sufficient to deform them or trigger any internal flows (the internal effective viscosity of the gel being much greater than viscosity of the water). In fact, from measuring marker positions no deformation could be detected within the experimental uncertainties, while the particles were carried by the turbulent flow. Therefore, these gel-particle can, for the purpose of this study, be considered as solid. To prevent the particles from being damaged by the propellers, we placed grids (of 2 × 2 mm 2 mesh size and 75% aperture) in front of both of them (see Fig. 1). Measurements showed that the influence on the flow was negligible. We ensured that the measurement volume was sufficiently far from the grid so that the large particles' dynamics was not influenced by the meshes. This can be seen from the following analysis: the viscous relaxation time of the particle is,

τ ν ≡ 1 18 ρ p ρ f d 2 p ν (1)
where ρ p and ρ f are respectively the particle and fluid density, d p the particle diameter and ν the kinematic viscosity of the fluid. This applies only when the particle Reynolds number (Re p ≡ u ′ d p /ν) is very small. Following [START_REF] Clift | Bubbles, drops, and particles[END_REF][START_REF] Xu | Motion of inertial particles with size larger than kolmogorov scale in turbulent flows[END_REF], one can show that the relevant time scale that takes into account small but finite particle Reynolds number can be written as:

τ p = τ ν 1 + 0.132Re n p (2)
where the exponent n also depends on Re p : n = 0.82 -0.05 log 10 Re p . With ρ p /ρ f ≈ 1, d p = 0.01 m, ν = 10 -6 m/s 2 , and u ′ = 0.1 m/s this time scale is τ p = 0.38 s, which yields L p = u ′ τ p = 3.8 cm beyond which the big, inertial particle has reached the fluid velocity. The minimal distance between the grid and the measurement volume is ≈ 7 cm on both sides and thus we expect that the particles' dynamic is not impacted by the screens.

As pointed out before, we were seeding the turbulent flow with two types of particles: big particles marked by fluorescent tracers and the fluorescent tracers themselves. When tracking all fluorescent particle tracks where visible simultaneously; marker particles trapped in the big particles did not change their separation, while tracer particles would separate quickly.

Identifying the finite-size particles

One very significant advantage of the novel technique presented here is that it allows to follow two species of particles using only one LPT system. For comparison, Guala et al. [START_REF] Guala | Experimental study on clustering of large particles in homogeneous turbulent flow[END_REF] used two LPT systems, which makes such an experiment very demanding. To identify the two populations of particles we took advantage of the fact that the flow was highly turbulent. It is well known that tracer particles separate very quickly [START_REF] Bourgoin | The role of pair dispersion in turbulent flow[END_REF][START_REF] Gibert | Inertial effects on two-particle separation in a turbulent flows[END_REF]. In the contrary markers fixed to the surface of a big particle will stay at a constant distance. Specifically, two tracer particles initially separated by a distance R 0 ≡ |R(t = 0)| in the inertial range (i.e. η ≪ R 0 ≪ L) will separate on average following a ballistic regime predicted by Batchelor [START_REF] Batchelor | The application of the similarity theory of turbulence to atmospheric diffusion[END_REF] such as:

δR • δR = 11C 2 3 R 2 0 t t 0 2 (3)
for t ≪ t 0 = (R 2 0 /ǫ) 1/3 , where δR(t) ≡ R(t) -R(t = 0) is the vectorial separation increment and C 2 = 2.1 as suggested from a compilation of available data [START_REF] K R Sreenivasan | On the universality of the Kolmogorov constant[END_REF].

Finding pairs of tracks

The first step is to find pairs of fluorescent particles tracks whose distance to each other stays constant during the time of observation. As the fluorescent markers are embedded in the particles with radius R p , we scan only for tracks whose separations are 0.05 < R t /R p < 2(here R p = 8 mm). The lower bound comes from the fluorescent particles injection process that we used. Additionally, to limit possible false choices, we only conduct this analysis on pairs of tracks that coexist for longer than 2τ η . This is the characteristic time-scale after which the acceleration components of a fluid particle decorrelated entirely [START_REF] Voth | Measurement of particle accelerations in fully developed turbulence[END_REF]. To estimate the tolerance on the separation increments, we have used Equ. 3 with R 0 = 0.05R p , t = 2τ η and ǫ = 0.011 m 2 /s 3 which gave us δR 2 ≈ 0.1R p . In the analysis presented here we used a slightly tighter criterion:

δR 2 < 0.01R p , that corresponds approximately to our spatial resolution (100 µm/pixels).

Finding groups

To ease the explanations, two tracks that were paired by the algorithm described above (section 3.1) are now called "friends". We can now identify the fluorescent particle tracks that belong to a unique inertial particle. To do so we simply apply the adage "the friend of a friend is my friend". This allows to generate groups of tracks that trace the inertial particles' motion. By additionally imposing that at least four tracks coexist for more than 2τ η we exclude optical artifacts. Such a group is presented on Fig. 3 (b).

Track reconnection inside the groups

As one can see from Fig. 3 (b), the tracks that belong to a group occur in short segments. This can be attributed to limitations of the experiment, which include variations of the illumination intensity, particles blocking the line of sight of the cameras, particles blocking the illumination beam (shadowing), the presence of light-insensitive circuits on the CMOS sensors of the cameras, and the effect of thermal, electronic and environmental noise (see [START_REF] Xu | Tracking lagrangian trajectories in position-velocity space[END_REF]). To reconnect the tracks we use an algorithm similar to that developed by Haitao Xu at the Max Planck Institute for Dynamics and Self-Organization [START_REF] Xu | Tracking lagrangian trajectories in position-velocity space[END_REF]. Within a group, at the end of each track, we define a cylindrical search volume, parametrized by its length C L , its radius C R and its axis oriented along V 1 (the velocity at the end of the track-segment). In Fig. 3 (a) such a cylinder is shown exemplarily in green. Tracks of the group that have their starting point within this search volume are considered potential re-connection candidates (track 2 on Fig. 3 (a)). Then for each candidate, we compute the linear interpolation of the trajectory and compare the beginning and ending interpolated velocities (resp. V 3 and V 4 on Fig. 3 (a)) to the actual measured ones (resp. V 1 and V 2 on Fig. 3 (a)). If the beginning and ending relative velocity differences are smaller than a tolerance C T (e.g.

|V 3 -V 1 |/|V 1 | < C T )
we consider these as a possible reconnection. If there are several candidates, we reconnect the best match. We remind the reader that we are working within a group obtained by the previous step (section 3.2), that is to say a very small amount of tracks compared to the work in [START_REF] Xu | Tracking lagrangian trajectories in position-velocity space[END_REF]. In our experiments, we found that C R = 0.5 mm, C L = 30 mm and C T = 0.5 give satisfactory results. For example, we the result of this reconnection process on a particular group in Fig. 3 (c), and its statistical effect on several measurements in Fig. 3 (d) (please note that this more than doubles average track length inside the groups).

Cleaning the groups

The process described above is not error-proof. Indeed, while scanning examples of groups by eye we could identify wrong tracks. We therefore developed the "group cleaning" technique. In this step, we scan each group and we keep only tracks that have more than three direct "friends" as defined in section 3.1. We use the same parameters except that we relax the constraint on the maximal separation increment to δR 2 < 0.07R p . This we can do as the data was already connected by linear interpolation.

Concluding remarks on the particle tracking procedure

As shown by Fig. 4 we can now separate the two populations of fluorescent particles. In this process, as discussed above, we interpolate trajectories (section 3.3) solely to connect tracks. These interpolated data points need to be excluded when statistics of the fluid velocity, or vorticity along tracks is considered. We would also like to note that an inspection by eye of ≈ 5000 groups did show the procedure described above works only well, when only one inertial particle is in the measurement volume at any given time. We therefore discarded manually the rare cases where two or more particles (some colliding) were in view. This did not limit the statistics significantly. In the future it will be not hard to modify the algorithm to include this automatically.

Translation of the finite size particle

Here, we describe how we extract the center position x c and radius R p of the big inertial particles for each group identified by the procedure described in section 3. These groups have a minimum of four marker-trajectories. As described above, our injection technique ensured reproducible and well defined penetration depth into the inertial particle. The injection was such that the markers were inserted not more than 0.5 mm from the polymer surface. The center position x c and radius R p of the inertial particle is entirely defined by the sphere equation:

(x i -x c ) 2 + (y i -y c ) 2 + (z i -z c ) 2 = R 2 p ( 4
)
where (x i ; y i ; z i ) are the i-th marker's 3D-coordinates. For each permutation four of these markers, one gets the following linear system:

P • (A B C D) T = Tr(X • X T ) (5)
with

P ≡      -1 x 1 y 1 z 1 -1 x 2 y 2 z 2 -1 x 3 y 3 z 3 -1 x 4 y 4 z 4      X ≡      x 1 y 1 z 1 x 2 y 2 z 2 x 3 y 3 z 3 x 4 y 4 z 4      (6)
and

A ≡ x 2 c + y 2 c + z 2 c -R 2 , B ≡ -2x c , C ≡ -2y c and D ≡ -2z c .
Let us consider a group that contains n ≥ 4 markers trajectories at a given time t, the number of permutations of 4 trajectories out of n is M = n! 4!(n-4)! = n 4 . Therefore, at each time step, we solve M times the linear system given by Equ. 5. We define x c,m (t) the center coordinates obtained by the m-th permutation and d m (t) = 4 i=1 4 j=i+1 |x i,m (t)x j,m (t)| the sum of the distances between the markers. To define the actual center coordinates of the inertial particles, we perform an average of x c,m (t) weighted by d m (t):

x c (t) = M m=1 x c,m (t)d m (t)/ M m=1 d m (t). ( 7
)
We have made this choice as to correct for a non-homogeneous distribution of the markers on the inertial particle's surface. Indeed, four markers homogeneously distributed around the sphere define a center position and a radius more accurately than four agglomerated points. The result of the center finding process on a particular group is shown on Fig. 5 (a), together with the radius of the inertial particle R(t), and the distances from the center position to the markers of this group r i (t) = |x i (t)-x c (t)| in Fig. 5 (b). The latter shows that we are able to measure the particle radius with an uncertainty of ±200µm, which corresponds to the uncertainty of the fluorescent marker insertion. We also observed that R(t) exhibits unrealistic rapid variations when a marker track appear or disappear from the group (around frames 200, 350, 800 and 900). These kinks are non-physical and result from the uncertainty of the marker position inside the inertial particle. Just as for tracer particles [START_REF] Mordant | Experimental lagrangian acceleration probability density function measurement[END_REF], we measured the velocity of inertial particles from the center trajectory by convolving with a properly normalized, truncated, differentiating Gaussian smoothing kernel. The characteristic time scale of the filter was chosen as σ = 3.5 frames = 0.13τ η sufficient to minimize the effect of noise on the tracer particles' acceleration variance. At the kinks observed in the center position of the big particles x c (t) (Fig. 5 (b)) the values of the velocities were very large, discontinuous, and thus unphysical. We flagged and excluded those from the analysis together with 2σ = 7 frames before and after those discontinuities.

Rotation of the inertial particle

To measure the rotation of the inertial particles, we used the center position together with the trajectories of the markers around the center of the sphere. To do so, we took advantage of an algorithm first introduced by W. Kabsch [START_REF] Kabsch | A solution for the best rotation to relate two sets of vectors[END_REF][START_REF] Kabsch | A discussion of the solution for the best rotation to relate two sets of vectors[END_REF] that was developed to compare molecular conformations in chemistry. This algorithm seeks the optimal rotation matrix U between two sets of points by minimizing the root mean square of their separation. Let us consider two sets x and y of N paired points that have the same centroid at the origin. We are looking for the 3 × 3 unitary matrix U that aligns best x with y. This can be achieved by minimizing the root mean squared:

D = 1 N N i=1 (Ux i -y i ) 2 1/2 (8)
That can be rewritten as:

ND 2 = N i=1 [(Ux i ) 2 + y 2 i -2Ux i y i ] (9)
where one can see that minimizing the left hand side is equivalent to maximize the negative term on the right hand side. If we represent the points of set x (resp. y) by a N × 3 matrix X (resp. Y ), the quantity to maximize is:

N i=1 Ux i y i = Tr(Y T UX) = Tr((XY T )U) ( 10
)
where XY T is a square 3 × 3 matrix, that can be rewritten using its singular value decomposition (SVD) as XY T = V SW T . V and W T are orthogonal matrices of the left and right eigenvectors of XY T and S is a diagonal 3 × 3 matrix containing the eigenvalues s 1 ≥ s 2 ≥ s 3 . Using the commutation properties of the trace operator, Equ. 10 can be rewritten as:

Tr(V SW T U) = Tr(SW T UV ) = 3 i=1 s i T ii (11)
where T = W T UV is a product of orthogonal matrices and itself an orthogonal matrix with det(T ) = ±1. Therefore each elements of T are equal or smaller than 1. Thus to maximize Eq. 11, T has to be the identity matrix T = I.

The last constraint that we have to consider is that

U = W T V T (12)
has to be a proper rotation matrix with det(U) = 1 (rows/columns of U have to form a right handed system) [START_REF] Kavraki | Optimal alignment for lRMSD using rotation matrices[END_REF]. For det(XY T ) > 0, it follows that det(U) = 1, but when det(XY T ) < 0, det(U) = -1. In the later case, one has to settle for the second largest value of Eq. 11. This value is obtained when T 33 = -1 since s 1 ≥ s 2 ≥ s 3 . This finally allows us to write the optimal rotation matrix U as

U = W    1 0 0 0 1 0 0 0 d    V T , (13)
where d = sign(det(XY T )).

When applying this algorithm to our data, we constructed the N × 3 matrices X ni = x ni (t) -x c,i (t) and Y ni = x ni (t + ∆t) -x c,i (t + ∆t), where x ni are the position components of the n-th marker trajectory that existed at times t and t + ∆t (∆t is equal to the acquisition time). By construction for each group N ≥ 4 (see section 3) thus ensuring a proper rotation matrix. Then for each time steps, we computed the covariance matrix C = XY T , its SVD C = V SW T and the sign of its determinant d = sign(det(C)). Next, we applied Equ. 13 to obtain the optimal rotation matrix U. From the rotation matrix, one can extract the rotation vector or the Euler angles and compute the rotation rate of the inertial particle ω. One can estimate an uncertainty of about 5% on the components of ω from the analysis of the residuals of the least square optimisation described above. Figure 6 (a) shows the measured probability density function of the components for the rotation of the inertial particle. This distribution is non Gaussian and can be described by a stretched exponential function P (x) = exp(3s 2 /2) with our results.

As shown in Fig. 6 (b), we measured simultaneously the full motion of a finite size particle, translation and rotation, and the highly turbulent flow field from tracer particle tracking. This method thus allows us to study experimentally the interactions between the particle and the flow.

Particle flow correlation

This novel measurement technique is particularly well suited to study particle fluid correlation. The quantity that we want to study here is the second moment of the longitudinal velocity difference between the inertial particle and the flow, also called the mixed (particle/flow) Eulerian second order longitudinal velocity structure function:

D mix LL (r) = [(v p (0) -v f (r)) • r] 2 (14)
where v p and v f are the inertial particle and fluid velocity vectors, and r is the separation vector between the center of the inertial particle and the tracer particle considered (r = r/|r|). In general, this quantity depends on the position of the particle and on the separation vector r. Here however, we measured at the center of a von Kármán mixer, where the flow is fairly homogeneous and isotropic, therefore one can assume that D mix LL is only a function of r. Additionally, given the symmetries of the system, we averaged the data azimuthally around v p axis, where we expect rotational invariance.

In order to quantify the influence of the inertial particle on the turbulent flow field, we have normalized D mix LL by its expected value D f low LL in the case of a point-like particle (a tracer/fluid particle since ρ p /ρ f ≈ 1):

Dmix LL = D mix LL /D f low LL = D mix LL /C 2 (ǫr) 2/3 (15)
where D f low LL = C 2 (ǫr) 2/3 is the second order longitudinal velocity structure function of a homogeneous and isotropic turbulent flow (neglecting intermittency corrections) for separations r in the inertial range (η ≪ r ≪ L).

Dmix

LL is shown in Fig. 7. The statistical convergence of this quantity, obtained by binning the space around the big particle, was checked by splitting the dataset in two sub-samples, leading to the same conclusions detailed hereafter. First, one can see that both limits r = R p and r -→ ∞ are consistent with what we would expect. Indeed, Dmix LL (r = R p ) = 0 simply means that no fluid is entering or exiting the solid particle, and Dmix LL (r -→ ∞) = 1 shows that the particle/flow interaction occurs locally. We can identify a spherical shell of radius 1 < r/R p < 3, where the correlation between the inertial particle and the flow is greatly enhanced. The numerical results by Naso et al. [START_REF] Naso | The interaction between a solid particle and a turbulent flow[END_REF] are in good agreement with our experimental observations. Furthermore, beyond an isotropic spherical shell of distance R p from the surface of the inertial particle we observed a weak anisotropy with respect to the particle direction. This could indicate the presence of a wake induced by the particle. LL as a function of the relative separation distance r/R p (where R p is the particle's radius) and θ being the angle between the separation vector r and the velocity vector of the ball v p . The orientation of v p is shown by the grey arrow (with an arbitrary length). The grey circle in the middle of the map symbolizes the inertial particle. Assuming axis-symmetry around the velocity axis of the particle, the lower half is a reproduction of the upper half mirrored on the center line. The black dashed lines show the separations r = 2R p , 3R p , 4R p , 5R p .

Conclusion and perspectives

The main advantages of this new measurement technique are:

(i) it requires nothing more than the usual equipment used for "standard" LPT or Tomo-PIV (3 or 4 fast cameras and a light source);

(ii) it allows to follow in three dimensions the full dynamics (translation and rotation) of a solid body (here a spherical particle) together with the flow field carrying it;

(iii) the solid object, being made out of super-absorbent polymer, has the same index of refraction as the fluid which makes it invisible and therefore it does not block the fields of view of the cameras (no shadowing effects);

(iv) the geometry of the solid object is recovered during the acquisition (in this case the particle radius R p ) and could be studied dynamically;

(v) several solid objects can be followed simultaneously in the measurement volume.

This list demonstrates the capabilities of this novel measurement technique. Here we applied it to the dynamics of big inertial particles in turbulent flows (see section 6) and identified the zone of interaction between the particles and the flow. Other questions that can be addressed are, for example, the study of collisions of particles with equal and different particle sizes, and the dynamics of non-spherical objects in turbulent flows. At high particle number density this technique can also be used to investigate the dynamics of wet granular matter and of the flow between the grains.

Figure 1 .

 1 Figure 1. Picture of the Von Kármán water flow. The green sphere at the center of the apparatus represents the measurement volume.

Figure 2 .

 2 Figure 2. Super-absorbent polymer particles. (a) Once immersed in water and in their dry form. (b) Beaker full of these particles crossed by a laser beam that is not refracted. The red luminescent dot in the center of the picture is a marker particle grafted at the surface of the solid particle (symbolized by the red dots in (a)).

Figure 3 .

 3 Figure 3. reconnection inside the groups. (a) Definition of the search volume for possible track reconnection. (b) Group before reconnection. (c) Group after reconnection. (d) Probability density function of the track length within the groups, before and after reconnection.

Figure 4 .

 4 Figure 4. Result of the process described in section 3 on a 3 s movie. (a) All the trajectories of the fluorescent particles; (b) Same data set in light blue. The red trajectories correspond to the fluorescent particles bounded to the surface of the super-absorbent polymer spheres identified by the procedure described in section 3.

Figure 5 .

 5 Figure 5. Results of the center finding process for one example group. (a) Trajectories (thin lines) identified to belong to the particle and the center (thick green line). (b) Radius (bold black line) and the distances from the center to each trajectory (solid lines) for the tracking time of the particle. In this specific case, the particle was tracked for 1115 frames = 41.3 τ η ≈ 0.4s). The inset at the bottom shows the 3 coordinates of particle center x c subtracted by its filtered version (Gaussian kernel of σ = 3.5 frames) x c , in [mm].

Figure 6 .

 6 Figure 6. (a) Reduced probability density function of the rotation rate component of the inertial particle (the different symbols corresponds to the 3 components ω x , ω y , ω z ). The dashed line corresponds to a stretched exponential function. (b) 3D visualization of the inertial particle trajectory together with the tracer particles around it. The colors encode the velocities of the particles, warm colors indicate fast particles and cold ones low velocity particles.

Figure 7 .

 7 Figure 7. Normalized mixed longitudinal 2nd-order Eulerian velocity structure function DmixLL . Dmix LL as a function of the relative separation distance r/R p (where R p is the particle's radius) and θ being the angle between the separation vector r and the velocity vector of the ball v p . The orientation of v p is shown by the grey arrow (with an arbitrary length). The grey circle in the middle of the map symbolizes the inertial particle. Assuming axis-symmetry around the velocity axis of the particle, the lower half is a reproduction of the upper half mirrored on the center line. The black dashed lines show the separations r = 2R p , 3R p , 4R p , 5R p .

  1/2 are the Kolmogorov length and time scales, respectively, where ν is the kinematic viscosity of the fluid. N f is the frame rate of the camera, in frames per τ η and R λ = (15u ′ L/ν) 1/2 is the Taylor scale Reynolds number of the flow.

√ 3 1 -erf ln(|x/ √ 3|)+2s 2 √ 2swith a best fit parameter s = 0.55. Zimmermann et al.[START_REF] Zimmermann | Rotational Intermittency and Turbulence Induced Lift Experienced by Large Particles in a Turbulent Flow[END_REF] used the same functional form with s = 0.45 which is in good agreement

Acknowledgments

We are very grateful to Haitao Xu for his many important contributions to the techniques used in this investigation. We thank J. -F. Pinton, A. Pumir and R. Zimmermann for helpful discussions. Support from COST Action MP0806 is kindly acknowledged. This work was funded generously by the Max Planck Society and the Marie Curie Fellowship, Program PEOPLE -Call FP7-PEOPLE-IEF-2008 Proposal N o 237521.


M A T H I E U G I B E R T

¶, A N D E B E R H A R D B O D E N S C H A T Z

Where do small, weakly inertial particles go in a turbulent flow?

We report experimental results on the dynamics of heavy particles of the size of the Kolmogorov-scale in a fully developed turbulent flow. The mixed Eulerian structure function of two-particle velocity and acceleration difference vectors δ r v • δ r a p was observed to increase significantly with particle inertia for identical flow conditions. We show that this increase is related to a preferential alignment between these dynamical quantities. With increasing particle density the probability for those two vectors to be collinear was observed to grow. We show that these results are consistent with the preferential sampling of strain-dominated regions by inertial particles.

Introduction

Turbulence occurs whenever fluid viscous forces are small compared to the dominant driving forces of the flow. In practice this includes most macroscopic natural and technological flows. Very frequently, these turbulent flows are loaded with small, passive particles with non-negligible inertia, i.e. particles do not exactly follow fluid motion. Familiar examples of such flows include clouds in the atmosphere, sand storms in arid regions, flow driven transports of grains in pharmaceutical and food industries, fluidized beds in chemical engineering, and fuel sprays in combustion engines and turbines. The understanding of particle-turbulence interaction is therefore of central importance to our ability to make advances in key economical and societal issues like energy generation, climate change, and pollution control. Comparing to the fluid turbulence problem itself, which is well described by the non-linear Navier-Stokes equations, the particleturbulence interaction problem is on much weaker foundations. In principle, the particle dynamics can be obtained by solving particle equations coupled with the Navier-Stokes equations for the fluid with no-slip boundary conditions on particle surfaces. However, this approach is a formidable task for theoretical and numerical investigations [START_REF] Happel | Low Reynolds Number Hydrodynamics[END_REF][START_REF] Homann | Finite-size effects in the dynamics of neutrally buoyant particles in turbulent flow[END_REF]. Most theoretical analyses and numerical simulations therefore must rely on simplified equations, in which the particle-turbulence interactions are modeled. An interesting result obtained from these simplified equations is that in a turbulent flow particles with small to moderate inertia and sizes smaller than the Kolomogrov length scale tend to move to regions with high strain and low vorticity [START_REF] Maxey | Equation of motion for a small rigid sphere in a nonuniform flow[END_REF]. This so-called "preferential concentration" of inertial particles is under intensive investigation [START_REF] Falkovich | Intermittent distribution of heavy particles in a turbulent flow[END_REF][START_REF] Zaichik | Statistical models for predicting pair dispersion and particle clustering in isotropic turbulence and their applications[END_REF][START_REF] Chun | Clustering of aerosol particles in isotropic turbulence[END_REF][START_REF] Collins | Reynolds number scaling of particle clustering in turbulent aerosols[END_REF][START_REF] Sundaram | Collision statistics in an isotropic particle-laden turbulent suspension .1. direct numerical simulations[END_REF][START_REF] Ducasse | Intermittent particle distribution in synthetic free-surface turbulent flows[END_REF] and has been invoked when studying droplet collision rates in warm clouds [START_REF] Shaw | Particle-turbulence interactions in atmospheric clouds[END_REF][START_REF] Falkovich | Acceleration of rain initiation by cloud turbulence[END_REF]) and plankton encounter rates in the ocean [START_REF] Schmitt | Intermittent turbulence and copepod dynamics: Increase in encounter rates through preferential concentration[END_REF]. Inhomogeneous distributions of inertial particles in turbulent flows have also been observed in several experiments [START_REF] Aliseda | Effect of preferential concentration on the settling velocity of heavy particles in homogeneous isotropic turbulence[END_REF][START_REF] Wood | Preferential concentration of particles in homogeneous and isotropic turbulence[END_REF][START_REF] Salazar | Experimental and numerical investigation of inertial particle clustering in isotropic turbulence[END_REF][START_REF] Saw | Inertial clustering of particles in high-reynolds-number turbulence[END_REF][START_REF] Monchaux | Preferential concentration of heavy particles: A voronoï analysis[END_REF]). Here we show direct experimental evidence that regions of high particle concentrations are also regions of high strain.

Experimental method & observation

We conducted three-dimensional Lagrangian Particle Tracking measurements [START_REF] Ouellette | A quantitative study of three-dimensional lagrangian particle tracking algorithms[END_REF][START_REF] Xu | Tracking lagrangian trajectories in position-velocity space[END_REF]) and obtained the trajectories of solid particles with different densities in a turbulent water flow. Below we give a brief summary of the experiment (see also [START_REF] Gibert | Inertial effects on two-particle separation in a turbulent flows[END_REF]). The turbulent flow was generated by two counter-rotating baffled discs of 28 cm in diameter. The turbulence chamber, shaped as a hexagonal cylinder, measured 40 cm along the axis of the propellers and 38 cm in the cross-section, both in height (vertically) and width (horizontally). The axis of rotation of the discs was in the horizontal direction. The strong sweeping of the fluid near the bottom surface resuspended the heavy particles with the effect that the particle concentration remained sufficient for the measurement of multi-particle statistics. A spherically shaped measurement volume with a diameter of 22 mm was located at the center of the water chamber, where the mean flow was small compared to the velocity fluctuations. The spatial resolution of 45 µm per pixel and the temporal resolution of 31 frames per Kolmogorov time τ η guaranteed accurate measurements of particle velocity and acceleration [START_REF] Xu | Acceleration correlations and pressure structure functions in high-reynolds number turbulence[END_REF]. As shown in the Table 1 we used three different types of particles: polystyrene, glass, and stainless steel with densities between 1 ρ p /ρ f 8. For the Taylor microscale Reynolds number of the turbulent flow (R λ = 442), the particles had the size of the Kolmogorov scale (η = 74 µm). The resulting Stokes numbers, a relative measure of the particle inertia, were between 0.08 and 0.5 (see Table 1). Other relevant parameters of the turbulent flow are the fluctuating velocity u ′ = 0.15 m/s, the integral length-scale L = 87 mm and the Kolmogorov time-scale τ η = 5.2 ms. We measured 3 × 10 7 data points for the polystyrene particles, 8 × 10 7 for the glass particles and 3 × 10 6 for the steel particles.The fluctuations in Fig. 1 can be attributed to the statistical uncertainty. For these heavy particles, we observed that the magnitude of the mixed velocityacceleration Eulerian structure function δ r v•δ r a p increased with particle density, where δ r v and δ r a p are the relative velocity and the relative acceleration between two particles separated at distance r (e.g. δ r v = v(x + r, t)v(x, t)). For fluid tracers this mixed velocity-acceleration structure function in the inertial range equals a constant: δ r u • δ r a = -2ǫ, where u and a are fluid velocity and acceleration, respectively, and ǫ is energy dissipation rate per unit mass [START_REF] Mann | Experimental study of relative[END_REF][START_REF] Falkovich | Particles and fields in fluid turbulence[END_REF][START_REF] Pumir | The lagrangian view of energy transfer in turbulent flow[END_REF][START_REF] Hill | Opportunities for use of exact statistical equations[END_REF]. We also found this in our experiment, where the polystyrene particles behaved as tracers [START_REF] Gibert | Inertial effects on two-particle separation in a turbulent flows[END_REF]. As shown in Figure 1, for inertial particles, this mixed velocity-acceleration structure functions in the inertial range were still constant, but the values increased significantly (up to 30%) when the particle Stokes number increased from 0.08 to 0.44.

As reported earlier [START_REF] Gibert | Inertial effects on two-particle separation in a turbulent flows[END_REF] the relative velocities in this range changed by less than 10%. In addition, the RMS acceleration changed with particle density from 6.9, 6.6 and 8.5 m/s 2 for tracer, glass and steel particles, respectively. Therefore, the observed change in δ r v • δ r a p cannot be simply explained by changes in the magnitudes of δ r v and δ r a p , but must be associated with the alignment between the two vectors. Indeed, we observed a change of the cosines

cos θ = δ r v • δ r a p |δ r v||δ r a p | . (2.1)
of the angle between these two vectors with inertia of the particles. Considering the tracer particles (polystyrene), its average value cos θ monotonically increases from -0.052 to -0.01 when increasing the observation scale r (12 < r/η < 120). For glass and steel particles, the average cos θ increased in absolute value by 25% and 9%, respectively.

Alignment of relative acceleration and relative velocity

The effect of particle inertia on cos θ is best demonstrated by its probability density function (PDF). Figure 2a shows the PDFs of cos θ corresponding to δ r v and δ r a p at separation r = (28 ± 7)η for the three different types of particles (similar behaviour is observed for other separations in the inertial range). All PDFs are skewed to the side of cos θ < 0, indicating that the two vectors are preferentially oriented in opposite directions, which gives the negative values of δ r v • δ r a p , as shown in Figure 1.

A very interesting observation from Figure 2a is the increase of the probability of | cos θ| ≈ 1 with particle inertia. This can be seen even better in Figure 2b where the PDFs of the two types of inertial particles were normalized by the PDF of the fluid tracers. Clearly, the particle inertia enhanced the collinearity, either in the same direction or in the opposite direction, between the relative velocity and the relative acceleration vectors. We note that the PDFs are slightly skewed toward the anti-alignment side (cos θ < 0).

How can we understand this observation? For this range of Stokes numbers, as shown numerically and experimentally [START_REF] Bec | Acceleration statistics of heavy particles in turbulence[END_REF][START_REF] Salazar | Inertial particle relative velocity statistics in homogeneous isotropic turbulence[END_REF][START_REF] Gibert | Inertial effects on two-particle separation in a turbulent flows[END_REF]), the difference between particle and fluid velocity and acceleration at the particle position is very small. The same is true for the dynamical variables δ r a p and δ r v. Thus, if the inertial particles were homogeneously distributed in space, like fluid tracers, we should not observe any differences between the PDFs of cos θ for fluid tracers and inertial particles. This is not what we observed. Therefore our measurements strongly indicate that inertial particles preferentially explore certain regions of the flow, leading to the observed PDFs in Figure 2.

Implication to preferential sampling

To further elucidate which features of the turbulence are responsible for the observed behaviour, we investigate δ r u and δ r a of fluid tracers for simple flow fields. Let us consider a simple model, namely a steady, incompressible linear velocity field u i defined by :

u i = M ij x j (4.1)
where u i is the fluid velocity at position x i , M ij is the velocity gradient tensor. In our experiments, where the separation r are in the inertial range, we could treat M ij as an effective velocity gradient at scale r. It has been shown that the effective velocity gradient at scale r evolves with a time scale t 0 ≡ (r 2 /ǫ) 1/3 [START_REF] Xu | The pirouette effect in turbulent flows[END_REF][START_REF] Pumir | Tetrad deformation and alignment of vorticity and strain in a turbulent flow[END_REF], which is much larger than the particle relaxation time τ p for r in the inertial range, as it can be quantified by noticing that the ratio τ p /t 0 = St(η/r) 2/3 ≪ 1. Therefore, to study

-1 -0.5 0 0.5 1 0 0.5 1 1.5 Eq. (4.2) λ 2 / =0 λ 1 λ 2 / =0.1 λ 1 λ 2 / =0.2 λ 1 λ 2 / =0.3 λ 1
Figure 3. PDF of cos θ as defined in Eq. 2.1. The relative velocity and acceleration are obtained from two fluid tracers randomly distributed in a steady, purely straining linear velocity field.

The different symbols corresponds to different ratios of λ2/λ1. The dashed-line corresponds to Eq. 4.2.

the dynamics of particle relative velocity and relative acceleration, the flow field may be treated as frozen.

For a steady flow, the fluid acceleration field a i is determined by M ij : a i = u j ∂ui ∂xj = M ik M kj x j . The expressions for two-particle relative velocity and acceleration, together with the cosine of the angle between them, cos θ, can then be obtained in closed form. For a velocity field with pure rotation, M ij = -M ji , it can be easily shown that δ r a • δ r u = 0, i.e. the fluid acceleration is always normal to its velocity. The PDF of cos θ is a δ-function at cos θ = 0. This is not consistent with our observations. For a purely straining velocity field, the three eigenvalues of M ij are all real and can be arranged as λ 1 λ 2 λ 3 . Incompressibility requires that λ 1 + λ 2 + λ 3 = 0. Then it is straightforward to show that the alignment between relative velocity and relative acceleration vectors depends only on the orientation of the separation vector r between the two particles and the ratio of the eigenvalues, e.g., λ 2 /λ 1 . The PDF of cos θ at a given λ 2 /λ 1 can then be obtained by assuming that r distributes uniformly in all directions within this purely straining flow. When λ 2 = 0, i.e. the two-dimensional case, the PDF is:

PDF(cos θ) = 1 π(1 -cos 2 θ) 1/2 . (4.2)
For more general cases, λ 2 /λ 1 = 0, we obtained the PDFs numerically by choosing randomly 10 7 separation vectors r and computing cos θ (Eq. 2.1) in the linear velocity field u i (Eq. 4.1). Please note, numerical analysis of various test cases showed that the shape of the PDF was insensitive to specifics of the particle distribution. It is well known that in turbulent flows the average of the eigenvalue ratio of the velocity gradient is positive, with a value of λ 2 /λ 1 ≈ 0.2 [START_REF] Betchov | An inequality concerning the production of vorticity in isotropic turbulence[END_REF][START_REF] Ashurst | Alignment of vorticity and scalar gradient with strain rate in simulated navier-stokes turbulence[END_REF][START_REF] Kholmyansky | Velocity derivatives in the atmospheric surface layer at Re λ = 10 4[END_REF][START_REF] Lüthi | Lagrangian measurement of vorticity dynamics in turbulent flow[END_REF][START_REF] Berg | Backwards and forwards relative dispersion in turbulent flow: An experimental investigation[END_REF]. For the effective velocity gradient in the inertial range, this ratio decreases monotonically from this positive value to 0 when the scale r increases from η to the integral scale [START_REF] Pumir | Tetrad deformation and alignment of vorticity and strain in a turbulent flow[END_REF]. We therefore performed simulations with the range of λ 2 /λ 1 around these values. Figure 3 shows the PDFs of cos θ from the simulation results with 0 λ 2 /λ 1 0.3, together with the analytical solution Eq. 4.2. The PDFs for λ 2 /λ 1 > 0 in Fig. 3 are remarkably similar to the "relative PDFs" shown in Figure 2b, namely, the relative accelerations and the relative velocities prefer to be collinear, with slight skewness to the anti-alignment side. This observation, together with the fact that inertial particles have, to leading order, the velocities and accelerations of the fluid at the particle positions [START_REF] Bec | Acceleration statistics of heavy particles in turbulence[END_REF][START_REF] Salazar | Inertial particle relative velocity statistics in homogeneous isotropic turbulence[END_REF][START_REF] Gibert | Inertial effects on two-particle separation in a turbulent flows[END_REF], suggests that inertial particles sample preferentially strain-dominated regions in turbulence. Our approach should be compared with results from numerical simulations of simplified equations, for which the phenomenon of "preferential sampling of straindominated regions" has been observed [START_REF] Chun | Clustering of aerosol particles in isotropic turbulence[END_REF][START_REF] Collins | Reynolds number scaling of particle clustering in turbulent aerosols[END_REF].

There the velocity fields around the virtual inertial point particles was investigated. It will be interesting to see how well the PDFs of cos θ agree with the experimental results presented here. Moreover, studying the PDF of cos θ for fluid tracers might reveal new information about the structures of the turbulent flow itself, e.g. by conditioning of the PDF on local strain rates and vorticities.

Summary

In summary, we observed a significant increase of the mixed velocity-acceleration structure function δ r v • δ r a p when increasing the inertia of heavy particles of the size of the Kolmogorov-scale. We found that this increase could not be attributed to the change in the magnitudes of the velocity and/or acceleration increments (|δ r v| and/or |δ r a p |). A significant contribution to this effect is from the alignment or anti-alignment between the relative velocity and relative acceleration. Using a simple model, i.e. a linear velocity field, we show that in the inertial range heavy particles preferentially explore strain-dominated regions of a turbulent flow.

The tendency of heavy particles to preferentially explore strain-dominated regions may also explain our recent observations [START_REF] Xu | Motion of inertial particles with size larger than kolmogorov scale in turbulent flows[END_REF]) that the time-averaged radial distribution function of particle positions increases when the center of the von Kármán flow is approached. In that case, the long-time averaged mean flow imposes deterministic strain at the center of the apparatus. With our new findings reported here, one would expect that the long-time average would lead to a measurable preferential concentration towards the center of the apparatus.

Our analysis does not rely on the commonly used and hard to measure long time average of the relative particle positions. Instead it only considers dynamical properties and thus opens a novel way to study particle turbulence interaction. For higher Stokes number, the filtering induced by the particle response time [START_REF] Bec | Acceleration statistics of heavy particles in turbulence[END_REF][START_REF] Salazar | Inertial particle relative velocity statistics in homogeneous isotropic turbulence[END_REF] will have an additional effect on the analysis presented above.

We are grateful to R. J. Hill, A. Pumir and E. -W. Saw for many interesting discussions. Support from COST Action MP0806 is kindly acknowledged. This work was funded by the Max Planck Society, and the Marie Curie Fellowship, Programme PEOPLE -Call FP7-PEOPLE-IEF-2008 Proposal N o 237521.

Figure 1 .

 1 Figure 1. Mixed Eulerian velocity-acceleration structure function δrv • δra normalized by the expected values for the fluid particles -2ǫ. The values of the different plateaus represented by the dot-dashed lines are : blue 1, red 1.15, black 1.30.

Figure 2 .

 2 Figure 2. (a) Probability density function of the cosine of the angle θ between δrv and δrap for the three types of particles for a separation r/η = 28 ± 7. (b) Same PDFs as in (a) normalized by the PDF of the tracer particles. (b) shows the relative effect of the particle inertia on the alignment/anti-alignment of δrv with δrap.

Table 1 .

 1 Characteristics of the different particles used in the experiments. The Stokes number is defined as St ≡ d 2 p /12βη 2 , where dp is the particle diameter, β ≡ 3ρ f /(2ρp + ρ f ) is the modified density ratio, that takes into account the added mass effect, and η is the Kolmogorov length scale of the turbulence.

		ρp/ρ f dp (µm) dp/η St
	Polystyrene 1.06 74 ± 10 0.98 0.08
	Glass	4	75 ± 8 1.04 0.24
	Steel	7.8 75 ± 15 1.04 0.44

† Current address Institut N ÉEL CNRS/UJF (Grenoble France) ‡ Email address for correspondence: mathieu.gibert@grenoble.cnrs.fr ¶ Email address for correspondence: Haitao.Xu@ds.mpg.de Email address for correspondence: Eberhard.Bodenschatz@ds.mpg.de


Brice Goglin
email: brice.goglin@inria.fr

Managing the Topology of Heterogeneous Cluster Nodes with Hardware Locality (hwloc)

Keywords: topology, locality, affinities, I/O devices, clusters, hwloc

Modern computing platforms are increasingly complex, with multiple cores, shared caches, and NUMA architectures. Parallel applications developers have to take locality into account before they can expect good efficiency on these platforms. Thus there is a strong need for a portable tool gathering and exposing this information. The Hardware Locality project (hwloc) offers a tree representation of the hardware based on the inclusion and localities of the CPU and memory resources. It is already widely used for affinity-based task placement in high performance computing.

In this article we present how hwloc is extended to describe more than computing and memory resources. Indeed, I/O device locality is becoming another important aspect of locality since high performance GPUs, network or InfiniBand interfaces possess privileged access to some of the cores and memory banks. hwloc integrates this knowledge into its topology representation and offers an interoperability API to extend existing libraries such as CUDA with locality information. We also describe how hwloc now helps process managers and batch schedulers to deal with the topology of multiple cluster nodes, together with compression for better scalability up to thousands of nodes.

I. INTRODUCTION

High performance computing relies on powerful computing nodes made of tens of cores and accelerators such as GPUs or Xeon Phi. The architecture of these servers is increasingly complex because these resources are interconnected by multiple levels of hierarchical shared caches and a NUMA memory interconnect. Execution performance now significantly depends on locality, i.e. where a task runs with respect to its data allocation in memory, or with respect to the other tasks it communicates with.

Performance optimization of parallel applications require a thorough knowledge of the hardware, and many research projects aim to model the platform to tackle this challenge. Besides analytical performance models, one solution consists in static modeling of the hardware resource organization. Indeed, parallel developers need such information to properly use the platform. hwloc (Hardware Locality) is the de facto standard software for representing CPU and memory resources, and for binding software tasks in a portable and abstracted manner [START_REF] Broquedis | hwloc: a Generic Framework for Managing Hardware Affinities in HPC Applications[END_REF].

However, the locality importance has grown and it now applies to high-performance I/O devices such as accelerators or network interfaces. Moreover, several batch schedulers or process managers try to manage clusters of such heterogeneous nodes in a global manner, making locality an important aspect, outside of nodes as well.

We present how hwloc has evolved into a central place for gathering locality information about all hardware subsystems in HPC servers. It achieves this goal by combining topology information from many sources, including operating systems, domain-specific libraries and platform-specific instructions. It interoperates with these sources by extending their interfaces with locality information about the devices they manage. hwloc also offers ways to manipulate multiple nodes topologies with the ability to avoid duplication in case of nearlyidentical cluster nodes.

The remainder of this paper is organized as follows: Section II introduces the challenges and use cases for providing topology information. Section III then summarizes the hwloc model and describes how it manages all sources of information. I/O device locality within heterogeneous nodes is then presented in Section IV while the management of multiple nodes topologies is described in Section V.

II. CONTEXT AND STATE OF THE ART

A. Why Locality matters and Where

Locality has been cited as a critical aspect of performance of parallel applications for a long time, from distributed computing [START_REF] Szalay | The importance of data locality in distributed computing applications[END_REF] to single servers [START_REF] Steckermeier | Using locality information in userlevel scheduling[END_REF]. The complexity of modern computing platforms is increasing, even inside commodity nodes. Figure 1 shows the hierarchical organization of resources within a widespread type of servers where some physical devices have affinities for some cores and memory banks. Developers and users have to take the hardware topology into account when trying to optimize their codes.

We identify two major types of affinity. First, tasks have affinities for hardware resources they use. This includes memory banks, caches and TLBs that contain some of their data as well as I/O devices such as accelerators and network interfaces. Moving a task from one core to another (or worse, from one NUMA node to another) usually causes the execution to slow down because of cache affinities. Thus, it is wellknown that computing tasks should be bound to a single core to avoid such migration. Migration can also cause the performance to vary depending on the cores' locality with regard to the I/O devices used by the task [START_REF] Moreaud | Impact of NUMA Effects on High-Speed Networking with Multi-Opteron Machines[END_REF]. The second kind of affinity is between tasks. Indeed, parallel applications often involve communication, synchronization and/or sharing between some of the processes or threads. It usually means that related tasks should be placed on neighbor cores to optimize the communication/synchronization performance between them [START_REF] Song | Feedback-Directed Thread Scheduling with Memory Considerations[END_REF]. However, the affinity can also be reversed when single tasks have strong needs. For instance, memory-intensive applications may want to avoid sharing memory links or caches with others [START_REF] Kim | Fair Cache Sharing and Partitioning in a Chip Multiprocessor Architecture[END_REF].

Moreover, some energy-based affinities may also be involved. Technologies such as Intel TurboBoost can improve sequential performance on partially-idle multicore processors, while some processors can be shutdown completely when entirely idle.

Applications can have several of these types of affinities simultaneously, even with conflicting needs. We envision two ways to deal with these needs. First, tasks can be placed on the hardware resources according to their affinities. For instance, MPI process placement based on the communication scheme and on the platform topology is a very active area of research [START_REF] Hoefler | Generic Topology Mapping Strategies for Large-scale Parallel Architectures[END_REF], [START_REF] Subramoni | Design of a Scalable InfiniBand Topology Service to Enable Network-Topology-Aware Placement of Processes[END_REF], [START_REF] Jeannot | Process placement in multicore clusters: Algorithmic issues and practical techniques[END_REF]. Then, the actual communication between tasks can be adapted to the existing placement. For instance, the existence and the size of a shared cache between processes can be a reason to switch from one communication strategy to another [START_REF] Buntinas | Cache-Efficient, Intranode Large-Message MPI Communication with MPICH2-Nemesis[END_REF], [START_REF] Ma | Locality and Topology aware Intra-node Communication Among Multicore CPUs[END_REF]. The locality of I/O devices can also be used to better tune collective operations [START_REF] Moreaud | Adaptive MPI Multirail Tuning for Non-Uniform Input/Output Access[END_REF], [START_REF] Goglin | Dodging Non-Uniform I/O Access in Hierarchical Collective Operations for Multicore Clusters[END_REF].

B. Many Sources of Hardware Information

Tackling locality issues within parallel applications actually involves three steps: gathering the hardware topology, expressing the software affinities, and matching one with another. We focus on the former in this article: how to gather, abstract and expose useful hardware topology information? The importance of locality led many developers to retrieve topology information within their applications or libraries. Unfortunately, this work is difficult because of the amount and variety of the sources of locality information, ranging from operating system, to direct hardware query and highlevel tools.

Linux is widely used in HPC. Unfortunately, its ability to report topology information was designed over more than ten years and therefore suffers from a partial and non-uniform interface. Many hardware details are available from the sysfs virtual file system (/sys) but it misses processor details (only available in /proc/cpuinfo) and I/O information such as network connectivity. Moreover, some of these files are in humanreadable format, while some other pieces of information are split into many different machine-readable files. Extracting locality information from an application is therefore a lot of work. Some processors have dedicated instructions for retrieving topology information such as cpuid on x86. However, applications relying on this feature need to be updated for every new micro-architecture because special values with new meanings are often added and have to be supported. Tools such as the Intel compiler and LIKWID1 use this idea and end up failing to discover the right topology on some custom platforms. The operating system usually takes care of these cases, so these processor-specific instructions should not be needed in topology-aware applications, as long as the OS is recent enough.

Convenient topology discovery should be available in higher-level libraries that hide the difficulty of parsing lowlevel system files or architecture-specific registers. On Linux, numactl2 possesses knowledge of NUMA, CPU and I/O locality but lacks caches. Moreover, its programming interface is unstable, and it was designed for binding tasks only: it cannot be used for querying details about hardware characteristics.

As shown on Figure 2, many libraries exist for querying the topology of specific subsystems, for instance pciutils for PCI 3 , libibverbs for InfiniBand, CUDA for NVIDIA GPUs, etc. Unfortunately, there is almost no interoperability between these libraries and other topology-related tools. Therefore, it is for instance necessary to query sysfs, pciutils and CUDA when looking for the locality of a NVIDIA GPU with regard to host CPUs. Some higher level tools such as lscpu or lshw 4 merge the information from several sources but they lack a programming interface. In brief, all these sources of information still have to be used concurrently for a developer to gather locality information about all hardware resources. Some non-Linux operating systems may have better interfaces but they lack part of the information. For instance, Solaris does not report cache information. There is therefore a need for a portable, system-wide topology discovery tool.

C. Execution and Memory Binding

Besides consulting topology information, the other important technical requirement for tackling locality is binding. Applications need ways to specify that a task or memory buffer should be allocated to one (or some) hardware resources. Many command-line binding tools exist, including numactl, taskset and schedtool on Linux. But most of them may bind tasks only. Moreover, they only operate on sets of logical processors without any knowledge of processor sockets, caches, etc.

Manipulating sets of logical processors unfortunately raises the issue of resource numbering. The BIOS and operating system are indeed free to renumber hardware devices, especially processor cores, based on their expectations of what the best numbering is. Memory bandwidth needs usually lead to numbering by NUMA node first, while sequential performance would likely number by hyperthread first. It leads to cases where a standard dual-socket platform can have up to 8 different numbering schemes depending on the vendor or BIOS version (see examples on Figure 3). Applications cannot be portable anymore if they rely on physical resource numbers.

A higher-level approach, based on the hierarchical organization, is required to keep affinity information: two cores sharing a L2 cache are considered neighbors even if their OS indexes are 0 and 4 respectively. Discovery and binding interfaces must therefore be integrated so that the same objects are manipulated for querying information about the platform hierarchy and resource characteristics, and for binding on these resources, without exposing numbering issues. We now present Figure 3. Numbering of the processing units (PU) on dual-socket dualcore hyperthreaded platforms. Two inter-dependent tasks running on logical processors 0 and 1 are actually not close to each other on these platforms. The binding cannot be portable unless it is specified as positions within the hierarchy of resources instead of as PU numbers.

the Hardware Locality project that was notably designed to solve this particular problem.

III. HWLOC'S VIEW OF THE HARDWARE

The Hardware Locality project was announced in 2009 as the replacement and merger of former Open MPI PLPA 5 and Inria libtopology 6 projects. It quickly raised attention of HPC runtime developers as an easy way to discover the topology of servers and to bind tasks. hwloc is now used by most MPI implementations, many batch schedulers and parallel libraries 7 . We summarize in this section the early design choices that led to hwloc success before explaining its evolution into the central place for information and interoperability about the topology of multiple hardware subsystems.

A. Organizing the Information hwloc resource organization is based on the natural inclusive order of computing resources: every machine contains one or several sockets, that contain one or several cores. hwloc builds a Tree of Objects describing these computing resources organized just like they are physically packaged. hwloc cores can actually contain multiple Processing Unit objects (PU), defined as the smallest resource that can execute a thread or a process. PUs correspond to logical processors or hardware threads as found in technologies such as simultaneous multithreading or Intel hyperthreading.

Each hwloc object is characterized by a type, some hardware characteristics such as a socket number, and some optional parameters such as local cache or memory sizes. The inclusionordering is extended to memory objects by considering that cores sharing a cache or near a NUMA memory node are included in it. Thus, the tree is made of a mix of levels made of computing and memory resources, ordered by locality without depending on actual physical numbering. hwloc does not enforce the vertical ordering between these levels in the tree because some AMD platforms have two NUMA nodes per socket (see Figure 4) while some Itanium machines have multiple sockets per NUMA node. hwloc just moves larger objects above smaller ones depending on the architecture inclusion characteristics. Sections IV and V explain how hwloc was recently extended to more than computing and memory resources inside nodes. The hwloc programming interface allows walking the tree edges to find ancestors or children objects of a given type (e.g. when looking for the NUMA node close to a given core), iter-ate over objects of a same type (e.g. when binding processes on cores), etc. hwloc offers a convenient way to apply mapping or partitioning algorithms by matching applications affinity graphs onto the hwloc tree of hardware resources [START_REF] Jeannot | Process placement in multicore clusters: Algorithmic issues and practical techniques[END_REF]. More use cases and hwloc v1.0 early design details are presented in [START_REF] Broquedis | hwloc: a Generic Framework for Managing Hardware Affinities in HPC Applications[END_REF]. In the rest of the paper, we focus on major improvements in later releases.

One critique against the model is its lack of topology information within single levels of the tree. For instance, Xeon E5 and Xeon Phi processors assemble cores on a ring, and the NUMA memory interconnect is not always a complete graph. Both are ignored when objects are represented as an array of children. To workaround this constraint, hwloc now annotates the tree with distance matrices and creates additional hierarchical Groups of object close to each-other. Large SGI Altix UV platforms are therefore represented with multiple levels of Groups between the machine and NUMA nodes so that the physical organization as racks and blades is exposed.

B. Orchestrating multiple Sources of Information 1) Combining multiple Sources:

As explained in Section II-B, multiple sources have to be used to gather all topology information about the machine. On Linux, virtual files under both /sys and /proc have to be used. x86-specific instruction may also bring more precise information about the CPU type, especially for non-Linux operating systems. I/O information involves several specific libraries as well as other virtual files under /sys. The hwloc library is therefore organized as several backend components. Discovery usually uses several components to match the aforementioned combination of sources as show on Figure 5 8 . Information is first gathered from operating system components (most Unix systems, Windows and Mac OS) that use OS interfaces. It is then extended by platform-specific components (BlueGene/Q, x86, Xen). I/O discovery is finally performed using specific libraries such as pciutils and CUDA. Each component can specify conflicts with others, and priorities can be changed to avoid a component that would return wrong information on a given platform. Moreover, some intercomponent callbacks can be specified so that PCI discovery immediately checks whether a new PCI device corresponds to one of the CUDA devices.

2) Interoperating with external Libraries: hwloc components can be built either statically inside the main hwloc library or as separate plugins. This is necessary for a convenient distribution of binary packages to avoid strong dependencies on external libraries. Indeed, binary packages should support as many cases as possible, which means hwloc should be built using all aforementioned I/O libraries such as pciutils and CUDA, but such dependencies are not acceptable for administrators that do not have GPUs on their platform. Building as plugins is an easy way to make these dependencies optional: all plugins are installed by binary packages but hwloc only loads plugins whose dependency libraries are available on the system. Obviously, it is still possible to build a custom hwloc library from source and embed all components that are useful to a given platform. hwloc uses external libraries to gather topology information. However, it was not designed to replace these libraries that offer a lot more of features unrelated to topology. hwloc was rather designed as a central place of topology details that interoperates with existing libraries and extend them with locality information. It therefore offers several interoperability headers that let developers translate between external library objects and hwloc data structures. For instance, an application using CUDA or InfiniBand verbs can retrieve the locality of CUdevice or struct ibv_device. Applications can therefore keep using existing specific libraries for non-topology-related information and switch to hwloc for topology-related queries.

IV. MANAGING HETEROGENEOUS NODES

We explained in the previous section how the hwloc library combines multiple sources of topology information from all subsystems in the machine. We now detail how it actually manages heterogeneous servers combining CPUs, accelerators such as GPUs or Xeon Phi, and network or InfiniBand interfaces.

A. I/O Discovery

I/O controllers are often placed near one of the processor socket within servers. They are even integrated inside sockets on recent Intel processors. Hence, devices connected to these controllers have a privileged access to the local memory and cores. These I/O affinities actually matter to latency or throughput sensitive applications that use high performance GPUs or network interfaces. Thus, it is important to offer I/O locality information to applications so as to optimize their placement and use of I/O devices [START_REF] Moreaud | Adaptive MPI Multirail Tuning for Non-Uniform Input/Output Access[END_REF], [START_REF] Goglin | Dodging Non-Uniform I/O Access in Hierarchical Collective Operations for Multicore Clusters[END_REF].

We added the ability to expose I/O device affinity in hwloc. The inclusion-based tree has been extended to attach new I/O objects under hwloc computing and memory resources they are close to (usually a NUMA node). Since high performance I/O is only significant for PCI devices, PCI is the only I/O hierarchy that is currently discovered by hwloc, using either an external PCI library such as pciutils or Linux sysfs files. PCI bridges are also discovered (see Figure 6) in case applications need to know which devices share PCI links and the speed of these links. But the tree may also be simplified to only retain the actual locality of PCI devices. The main issue with I/O discovery is that applications do not manipulate PCI devices, they operate on software handles instead, such as network socket file descriptor, InfiniBand ibv_device, CUDA CUdevice, etc. When multiple similar devices are available in the system, finding which hardware device corresponds to the application software handle can be difficult. hwloc solves this issue by inserting OS device objects describing these software handles inside PCI device objects, and adding human readable types and names. For instance, a cuda0 device is inserted so that the locality of CUDA device #0 can be retrieved by walking up the tree across the PCI hierarchy up to NUMA node #0 on Figure 6.

B. Identifying Objects from Applications

One immediate advantage of this feature is for binding microbenchmarks. Instead of manually binding a network ping-pong to a core near the InfiniBand interface mlx4_0, binding can be performed automatically near a specific OS device with the hwloc-bind tool:

$ hwloc-bind os=mlx4_0 pingpong_benchmark $ hwloc-bind os=cuda1 cuda_benchmark $ hwloc-bind os=mic0 xeon_phi_benchmark

C. Identifying Objects from outside the Host

We described in the previous section how to identify I/O devices from host applications. We now look at identifying them from other points of view. The first use case is for matching Xeon Phi boards as viewed from the host and from inside the board. Indeed, one way to use Xeon Phi is to mix MPI ranks on the host CPUs and on the Phi. MPI communication have to be implemented depending on whether the Phi and the CPU are located inside the same server or not. One requirement from MPI implementers is therefore to identify which MPI ranks are inside the same server. We solved this issue by extracting the Phi serial number and making it available in the hwloc topologies of the host (in the Phi OS device) and of the Phi itself (in the root object). This is already in use in Open MPI 1.7.

Another case where devices have to be identified from outside the host is network and InfiniBand interfaces. MAC/IP addresses and InfiniBand GUID/LIDs respectively are the only way to refer to a specific remote host, especially when multiple interfaces or ports exist. We therefore added to hwloc I/O objects several additional attributes enabling such identification. This feature is already used by the netloc project as discussed in Section V-C.

Examples of InfiniBand and Xeon Phi I/O-specific attributes added to hwloc OS devices are presented on Figure 7.

Co-Processor L#5 (CoProcType=MIC MICFamily=x100

MICSerialNumber=ADKC32800176 MICActiveCores=61 MICMemorySize=16252928 ...) "mic0" OpenFabrics L#8 (NodeGUID=f452:1403:007a:7260 Port1GID0=fe80:0000:0000:0000:f452:1403:007a:7261 Port1State=4 Port1LID=0x1 Port1LMC=0) "mlx4_0"

D. Being Generic enough

hwloc represents CPU and memory objects using an exhaustive set of widespread resource types (PU, core, cache, socket, NUMA node, machine) as well as additional generic objects (such as Groups for describing intermediate affinity neighborhoods in the tree). I/O objects raised the need to provide even more convenient types. Unfortunately, there are many different types of I/O devices and hwloc cannot list all of them explicitly. Moreover, many of these objects have specific attributes such as the memory size, the cache type, the network address, etc. Therefore, there is a need for a generic way to annotate hwloc objects with custom attributes instead of adding many hardwired type-specific structures of attributes.

hwloc usually gathers about one hundred object attributes for an entire server. They are attached to the relevant objects within the tree or to the root object when the attribute applies to the entire topology. These generic attributes are stored as a pair of key and value strings, such as Ad-dress=00:11:22:33:44:55. This presents the drawback of requiring string manipulations for applications. As a consequence, widely-used attributes (such as cache sizes) are stored using explicit fields within the object structure instead of generic key/value string pairs. Each object needs a way to store these key/value pairs. The number may vary but it always remains small, usually lower than 10. So there is no need to optimize lookup with advanced data structures such as hash tables. hwloc just uses an array of key/value pairs in each object. This mechanism is already widely used in hwloc since many users requested new attributes to be added. For instance, Intel and Oraclespecific drivers consult hwloc CPU attributes to dynamically optimize the Open MPI implementations on the corresponding platforms.

V. MANAGING CLUSTERS OF NODES

We now look at managing with hwloc the topologies of multiple nodes, such as a cluster. This is used for batch schedulers such as Slurm or Torque and process launchers found in most MPI implementations [START_REF] Hursey | Advancing Application Process Affinity Experimentation: Open MPI's LAMA-Based Affinity Interface[END_REF]. They have to know how many cores each node features before deciding how many processes should run. Moreover knowing topology details let them allocate resources better. The topology of multiple nodes can then be combined to build a global cluster-wide topology, where placement algorithms can perform both inter-node and intra-node management simultaneously [START_REF] Subramoni | Design of a Scalable InfiniBand Topology Service to Enable Network-Topology-Aware Placement of Processes[END_REF], [START_REF] Jeannot | Process placement in multicore clusters: Algorithmic issues and practical techniques[END_REF].

A. Remote Node Topology

Managing the topology of multiple nodes requires a way to manipulate remote node topologies. Since most topology information pieces are read from virtual files on Linux, hwloc gained the ability to change the filesystem root path in order to use copies of the /sys and /proc files from another node. This feature is convenient for debugging the Linux discovery code without immediate access to a remote node but hwloc offers a more convenient solution: an API to import/export entire topologies to XML, either as a file, or as a memory buffer that can be transferred on the network. This is useful for developing topology-aware algorithms and testing on a variety of different platform topologies. But it is also already widely used by MPI process launchers: each compute node sends a XML copy of its local topology to the frontend node which implements the process placement algorithms cluster-wide, before actually starting processes on the compute nodes.

XML also has the advantage of being very easy to load, much easier than rereading topology information from the different sources as explained in Section II-B. Discovery the topology natively on Linux indeed reads information from several hundreds of files under /sys and /proc. A naive MPI implementation running one process per core would load the topology once per core, causing all these files to be accesses by all cores simultaneously. Table I shows that the native Linux discovery does not scale well with the number of cores working in parallel (contention in the Linux kernel filesystem locking code) while XML import scales well. It also shows that very large machines may benefit from always loading from XML (up to 70x faster) even when a single discovery is performed simultaneously.

B. Cluster Nodes are (almost) Identical

Once compute node topology has been retrieved on a master node, one may wonder if storing all of them locally scales to a high number of nodes. Moreover, cluster nodes are usually similar: clusters are made of a single (or few) types of nodes. Why storing the topologies of all nodes if most of them are identical? We identified three actual possible differences between cluster node topologies:

• different kinds of nodes (e.g. compute node vs fat node): topologies are very different;

• modified nodes (BIOS upgrade, software update, or hardware replacement): topologies can be different;

• similar nodes with different identification numbers such as MAC address, InfiniBand GUID, etc.

In the similar case, only for key/value pairs are modified. In other cases, the tree structure can be different. Therefore, we added to hwloc the ability to compute the difference between 2 similar nodes by recording which key/value pairs have been modified. This loss-less compression consist in identifying a few reference nodes whose topologies will be entirely stored (uncompressed). All other nodes are then compressed by only storing the difference between their topology and one of the references. This feature is already used in the netloc project.

Table II presents the memory occupancy improvement based on the compression of the topologies of two clusters of the University of Bordeaux9 . Each cluster is made of different kinds of nodes (6 for Plafrim and 5 for Avakas), but we observe more reference topologies (respectively 18 and 12) because of the modified case above. However, many topologies can indeed be reduced from several hundreds of kilobytes down to 1 or 2 kB in memory. Full topologies seem bigger in the compressed case because the share of fat nodes among reference topologies is higher.

Each difference is actually made of about 10 key/value pair changes. We could even improve compression further by ignoring keys that are not needed by the target application (for instance the platform serial number, or the MAC addresses if only InfiniBand is used).

C. Multiple Node Topology

Finally, we look at how to manage a full, cluster-wide topology. hwloc offers an API to assemble the topologies of multiple nodes into a global single one. However, the resulting topology must respect hwloc's tree model while networks interconnecting nodes can be random graphs. We explained in Section III-A that distance matrices can be used to annotate some levels of the tree but this idea is only satisfying for simple topologies such as NUMA interconnects or socket rings.

Moreover, cluster nodes are interconnected by NICs or InfiniBand HCAs at the bottom of the tree, not by the hwloc tree roots (the entire machine object). Therefore, assembling multiple nodes into a global hwloc topology does not seem convenient. That is why there is an ongoing work to develop a hwloc-companion called netloc 10 to combine hwloc node topologies with network graphs without enforcing a tree model [START_REF] Goglin | netloc: Towards a Comprehensive View of the HPC System Topology[END_REF].

VI. CONCLUSION AND FUTURE WORKS

The increasing complexity of computing platforms raises the need for developers to understand the hardware organization and adapt their application layout. As part of the overall optimization process, there is a strong need for a tool modeling the platform, and hwloc is the most popular software for exposing a static view of the topology of CPUs and memory.

We presented in this article how we have extended hwloc to more than these computing resources by also incorporating the topology of I/O devices and offering ways to manage multiple nodes. hwloc now integrates locality information from many sources and offers APIs to interoperate with these libraries and operating systems without replacing them. I/O locality information has been added to the hwloc tree representing the hardware as well as many attributes to help applications identify the resources they use, place tasks near them or adapt to their locality. An API to manipulate the topologies of remote hosts with compression for better scalability was also recently added.

All features listed in this paper are available in hwloc v1.9 (released in Spring 2014) 11 . On-going work is now focusing on improving topology detection on emerging ARM architectures for high-performance computing as well as automatic management of conflicts between redundant sources of information.

Figure 1 .

 1 Figure1. Topology of a dual-Xeon E5 host with GPUs (cuda0, cuda1), network (eth1), InfiniBand (mlx4 0) and disk (sda) connected to different sockets, simplified and reported by hwloc's lstopo tool.

Figure 2 .

 2 Figure 2. Overview of existing sources of locality information on Linux.

  4 http://ezix.org/project/wiki/HardwareLiSter PU numbering by NUMA node first, then by core, then by PU. PU numbering by core first, then by PU, then by NUMA node.

Figure 4 .

 4 Figure 4. AMD platform containing a single Opteron 6272 socket. hwloc's inclusion ordering is machine, socket, NUMA node, L3, L2, L1i, L1d, core, PU.

Figure 5 .

 5 Figure 5. hwloc's component-based organization.

Figure 6 .

 6 Figure 6. I/O device hierarchy connected to a NUMA node. Grey boxes inside dark green boxes on the right are OS device objects inside PCI devices. Small squares represent bridges, and decimal values are PCI link speeds in GB/s.

Figure 7 .

 7 Figure 7. Textual dump of some attributes gathered by hwloc for OS devices describing a Xeon Phi (mic0) and a InfiniBand HCA (mlx4 0).

Figure 8 .

 8 Figure 8. Object attributes include cache types (L1i, L3, etc), memory sizes, object numbers, PCI device and vendor numbers, PCI link speed, Xeon Phi memory and cores, CUDA memory and multiprocessors, as well as CPU vendor and model.

TABLE I HWLOC

 I TOPOLOGY DISCOVERY TIME DEPENDING ON THE SOURCE, EITHER NATIVE LINUX DISCOVERY, OR XML IMPORT. ON EACH HOST, WE MEASURE THE TIME FOR A SINGLE DISCOVERY AND FOR ALL CORES DISCOVERING SIMULTANEOUSLY.

	Host	16 cores	16 cores	160 cores
		without I/O	with 3 GPUs	SGI Altix UV
	# Processes	1	16	1	16	1	160
	Linux	26 ms	1 s	210 ms	6 s	390 ms 107 s
	XML	3 ms	7 ms	3 ms	7 ms	12 ms	22 ms

http://code.google.com/p/likwid/

http://oss.sgi.com/projects/libnuma/

http://mj.ucw.cz/sw/pciutils/

http://www.open-mpi.org/projects/plpa/

http://runtime.bordeaux.inria.fr/hwloc/

A non-exhaustive list of hwloc users is available on the project webpage http://www.open-mpi.org/projects/hwloc/.

The XML import backend cannot be combined, it is a global component that manages all objects, as explained in Section V-A).

The hwloc-compress-dir utility was used.

Available under the BSD license at http://www.open-mpi.org/projects/ netloc/.

ACKNOWLEDGMENTS

Some of cluster node topologies used in this study were provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour.

Some experiments presented in this paper were carried out using the PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from LABRI and IMB and other entities: Conseil Régional d'Aquitaine, FeDER, Université de Bordeaux and CNRS.