Jean-François Chaix
email: jean-francois.chaix@univmed.fr

Mathieu Rossat

Vincent Garnier

An experimental evaluation of two effective medium theories for ultrasonic wave propagation in concrete

Keywords: numbers: 43.35.Cg, 43.20.Gp, 43.20.Jr

This study compares ultrasonic wave propagation modeling and experimental data in concrete. As a consequence of its composition and manufacturing process, this material has a high elastic scattering (sand and aggregates), and air (microcracks and porosities) content.

The behavior of the "Waterman-Truell" and "Generalized Self Consistent Method" dynamic homogenization models are analyzed in the context of an application for strongly heterogeneous solid materials, in which the scatterers are of various concentrations and types. The experimental validations of results predicted by the models are carried out by making use of the phase velocity and the attenuation of longitudinal waves, as measured by an immersed transmission setup. The test specimen material has a cement-like matrix containing spherical inclusions of air or glass, with radius close to the ultrasonic wavelength. The models are adapted to the case of materials presenting several types of scattering particle, and allow the propagation of longitudinal waves to be described at the scale of materials such as concrete. The validity limits for frequency and for particle volume ratio, can be approached through the comparison with experimental data. The potential of these homogenization models, for the prediction of phase velocity and attenuation in strongly heterogeneous solids, is demonstrated.

I. INTRODUCTION

The general context of this study is the evaluation of damage in concrete structures using ultrasonic methods. Damage often leads to the appearance and then to the development of microcracks and/or microporosities in the medium [START_REF] Grattan-Bellew | Micro-structural investigation of deteriorated Portland cement concretes[END_REF][START_REF] Handoo | Physicochemical, mineralogical, and morphological characteristics of concrete exposed to elevated temperatures[END_REF] . This type of damaged medium can be described as a heterogeneous medium composed of an elastic matrix (the cement) containing elastic inclusions (the aggregates) and air inclusions (the microcracks).

Interaction between the wave and the obstacles leads to scattering of the wave. In the case of concrete, in which there are numerous obstacles, it is necessary to take multiple scattering effects into account [START_REF] Chaix | Ultrasonic wave propagation in heterogeneous solid media : Theoretical analysis and experimental validation[END_REF][START_REF] Philippidis | Experimental study of waves dispersion and attenuation in concrete[END_REF] .

The observable parameters, which are the velocity and attenuation of ultrasonic waves, are difficult to relate to the size of the aggregates, the volume fraction of porosities, the elastic constants, and damage in the medium. The understanding of the direct problem in concrete requires the propagation phenomena of ultrasonic waves, in heterogeneous media containing high densities of different types and shapes of scattering particle, to be modeled. The experimental validation phase of the models not only demonstrates the potential and limitations of this type of model, but also highlights the work which still remains to be accomplished [START_REF] Chaix | Ultrasonic wave propagation in heterogeneous solid media : Theoretical analysis and experimental validation[END_REF][START_REF] Molero | Characterization of mortar samples using ultrasonic scattering attenuation[END_REF][START_REF] Kim | Characterization of ultrasonic properties of concrete[END_REF] .

The prediction of propagation in a multi-scattering particle medium can be made using the exact calculation of the field scattered by all of the obstacles. However, the large number of scattering particles and the need to know the position, type, size and shape of each of them rapidly lead this technique to become costly in terms of computation time, and very difficult to implement whenever some of the parameters are random, such as the exact positions and shapes of the scatterers. Medium homogenization methods, which integrate the multiple scattering equations, can thus be a good modeling solution for a wave propagating through an acoustically equivalent homogeneous medium. There are numerous medium homogenization models. The simplest ones are independent of frequency, take only simple scattering into account, and are of the static or quasi static type: Kuster and Toksöz 7 , Berryman [START_REF] Berryman | Long-wavelength propagation in composite elastic media: Spherical inclusions[END_REF] for example. These models are limited to use in the Rayleigh domain, in which the wavelength is much greater than the size of the obstacle, and for scattering particle concentrations which are often low.

If the size of the inclusions is of the order of one wavelength, the wave propagation characteristics such as phase velocity and attenuation become frequency dependent [START_REF] Chaix | Ultrasonic wave propagation in heterogeneous solid media : Theoretical analysis and experimental validation[END_REF][START_REF] Philippidis | Experimental study of waves dispersion and attenuation in concrete[END_REF][START_REF] Molero | Characterization of mortar samples using ultrasonic scattering attenuation[END_REF][START_REF] Kim | Characterization of ultrasonic properties of concrete[END_REF] . In the stochastic domain, many dynamic models provide a good description of propagation.

The integration of multiple scattering effects combined with the application of a statistical mean to the positions of isotropic scattering particles was initiated by Foldy [START_REF] Foldy | The multiple scattering of waves[END_REF] . The extension of this model to non-isotropic scattering particles was proposed by Waterman and Truell [START_REF] Waterman | Multiple scattering of waves[END_REF] and Twersky [START_REF] Twersky | On propagation in random media of discrete scatterers[END_REF] , and modified versions of these methods were proposed more recently by Linton and Martin [START_REF] Linton | Multiple scattering by multiple spheres: A new proof of the Lloyd-Berry formula for the effective wavenumber[END_REF] . In parallel with these developments, hypotheses concerning the relative positions of the scattering particles were formulated by Lax [START_REF] Lax | Multiple scattering of waves[END_REF] , who proposed the quasi-crystalline approximation (QCA). This approximation was coupled to different correlation functions, which led to various models such as those of Tsang 14 , Varadan [START_REF] Varadan | A multiple scattering theory for elastic wave propagation in discrete random media[END_REF] , and Sabina and Willis [START_REF] Sabina | A simple self-consistent analysis of wave propagation in particulate composites[END_REF] . For these models, the choice of inter-correlation functions is important, and solutions are difficult to obtain in the short wavelength domain. The numerical computations thus become complex and their convergence towards a solution is no longer certain.

Finally, other methods rely on one or another of the two preceding families, but use an iterative process which makes it possible to converge towards the wavenumber of the equivalent medium: Kim [START_REF] Kim | Dispersion of elastic waves in random particulate composites[END_REF] proposed an initial auto-coherent scenario, which he called the Coherent Potential Approximation (CPA). Another approach based on the same principle was proposed by Yang [START_REF] Yang | A Dynamic Generalized Self-Consistent Model for Wave Propagation in Particulate Composites[END_REF] : the Generalized Self Consistent Model (GSCM). Its results are close to those of Waterman and Truell's model [START_REF] Waterman | Multiple scattering of waves[END_REF] and to the experimental comparisons with lead/epoxy composites carried out by Kinra and Rousseau [START_REF] Kinra | Acoustical and Optical Branches of Wave Propagation[END_REF] . These auto-coherent methods are described as having potentially good performance in strongly heterogeneous media.

In the context of an application in solid media with high scattering particle concentrations, this study is based on the Waterman and Truell model (WT), the formalism of which is by far the simplest, and which has already demonstrated promising results [START_REF] Chaix | Ultrasonic wave propagation in heterogeneous solid media : Theoretical analysis and experimental validation[END_REF] , and the Generalized Self Consistent Model (GSCM) proposed by Yang which has a strong potential for the study of highly charged media.

More generally, the experimental validations related to these models are numerous in the case of mediums with a fluid matrix, but more rare in the case of solid matrices. In this case, the models' limitations, in terms of the frequency range or the volume fraction of the obstacles in solids, need to be defined.

The application of these models to the case of concrete is envisaged. As a consequence of the shape and the nature of the obstacles, the problem is analyzed for three dimensions objects over a wide range of frequencies, in order to cover the potential domain of ultrasonic testing applications. The cases of elastic inclusions, as well as that of air inclusions and the combination of both types of inclusions, are studied. This paper proposes to compare the modeling results of the two effective theories with experimental data obtained in concrete media. The long-term objective is to be able to predict the acoustic behavior of waves in damaged concrete.

In Section II, we describe the models we used, their extension, allowing several types of inclusions to be accounted for, and their behavior with respect to the media analyzed in the present study. The experimental study, described in Section III, presents the test specimens, the acquisition protocol and the measurement principle. In Section IV, we compare the results provided by the models and the experiments. Finally, our conclusions are given in Section V.

II. DYNAMIC HOMOGENIZATION MODEL

A. Homogenization principle

In order to compute the field scattered by the full set of obstacles, the theory of equivalent media, or homogenization, allows a heterogeneous medium to be replaced by an acoustically equivalent homogeneous medium (Fig. 1).

These propagation models are based on the statistical average of multiple scattering equations, relying on the acoustic scattering behavior of a single object. Spherical geometries are considered in this study, because of their similarity to the shape of real objects present in concrete, such as aggregates or porosities. Although the microcracks are less similar to spheres, nevertheless we assume the hypothesis that the scattering effects arising from a set of randomly oriented microcracks can have similar behaviours than a set of randomly distributed spheres. So we describe the scatterers by a distribution of spheres of rocks or air.

The propagation characteristics of the longitudinal wave, i.e. the phase velocity * c and the attenuation * α are respectively the real and imaginary parts of the complex wavenumber * k (the asterisk indicates the equivalent medium):

* * * . i c k α + ω = (1) where f . . 2 π = ω
is an angular velocity and f is the frequency. This wavenumber for the equivalent homogeneous medium, is this one defined using the dynamic homogenization models, and depends on the wavenumber in the matrix, the ones in the inclusions, and the sizes and densities of the inclusions.

B. The Waterman-Truell model

In the initial approach (Fig. 2) the scattering particles are all identical, with a radius a, a volume fraction inside the medium which is assumed to be constant and equal to τ, and positions which are statistically independent. In the case of spherical scattering particles, the results for a longitudinal wave propagating in the equivalent medium lead to the following relationship for the complex wavenumber 10 :

( ) ( ) 2 2 1 0 2 2 1 0 2 1 * k f . n . . 2 k 0 f . n . . 2 1 k k π π - π + = (2) where 3 0 a . . 4 . 3 n π τ =
is the number of scattering particles per unit volume, K 1 is the wavenumber in the matrix and f(0) and f(π)are the forward and backward form function for a single scattering particle. The scattering properties of a single object are characterized by the form function. This gives the angular distribution of the amplitude of the scattered wave surrounding the obstacle, in particular in the forward and backward directions. In the case of longitudinal waves the form function, for an observation angle θ is given by the following expression 20 :

( ) ( )( ) ( ) ( ) ∞ = θ + = θ 0 n 0 n 0 n 0 n 11 1 cos P . T . 1 n . 2 . k . i 1 f (3)
This far-field expression is based on the decomposition of potentials using the Legendre polynomials ( 0 n P ), and on the use of the T-matrix ( T ) [START_REF] Waterman | Matrix theory of elastic wave scattering[END_REF] in the scattering calculations. ( )

0 n 0 n 11 T
is the coefficient, of n order, derived from this matrix, corresponding to the longitudinal wave scattering of a longitudinal wave. The expression for this coefficient in the case of a spherical, fluid or elastic object in an elastic medium is computed from the ultrasound parameters of the matrix and the scattering particles using a resonant scattering theory [START_REF] Brill | Resonance theory of elastic waves ultrasonically scattered from an elastic sphere[END_REF] . This technique is widely used in the scattering calculations and had been experimentally validated for the case of elastic sphere in elastic medium [START_REF] Sessarego | Scattering by an elastic sphere embedded in an elastic isotropic medium[END_REF] with very good agreement.

C. The "Generalized Self Consistent Method" model

Yang [START_REF] Yang | A Dynamic Generalized Self-Consistent Model for Wave Propagation in Particulate Composites[END_REF] considered a scattering medium composed of an equivalent matrix in which scattering particles are included (Fig. 3). As opposed to the previous model, the equivalent matrix corresponds to a homogenized medium instead of the real matrix. This is the hypothesis which makes the model self-coherent. It is an iterative approach which converges, tending towards the characteristics of the equivalent medium.

The scattering particles are all identical, with a radius a, and have a volume fraction in the medium assumed to be constant and equal to τ, and statistically independent positions. These scattering articles are thus considered to be inside a shell of characteristic size b having the same acoustic properties as the real matrix.

In the case of spherical diffusers, the radius b of the shell is related to the radius a of the obstacle by the obstacle volume fraction in the medium, written as:

3 3 b a = τ (4)
The result of the "Generalized Self Consistent Method" is obtained from equation ( 2), with the wavenumber of the matrix k 1 being replaced by that of the homogenized medium k*. For the longitudinal wave propagating in the equivalent medium, we thus obtain the following expression:

( ) ( ) 2 2 * gscm 0 2 2 * gscm 0 k f . n . . 2 k 0 f . n . . 2 1 1 π π - π + = (5) where ( ) ( ) ∞ = - = 0 n n 0 n gscm A . i 0 f and ( ) ( ) ∞ = = π 0 n n 0 n gscm A . i f (6)
These are the forward and backward scattering functions of a longitudinal wave impinging on a spherical obstacle of radius a, included inside a shell of radius b, which itself is embedded in an equivalent matrix. These functions are calculated starting from the scattering coefficient A 0n defined by Yang [START_REF] Yang | A Dynamic Generalized Self-Consistent Model for Wave Propagation in Particulate Composites[END_REF] . Taking into account inclusions inside a shell of radius b (eq. 4) depending on volume fraction of inclusion induces modifications of scattering functions. This way to model inclusions in 2D effective medium model coupled with the self-consistent scheme provided good agreements with experiment in the SiC/Ti fiber/matrix composite [START_REF] Huang | Frequency Dependencies on Ultrasonic Wave Velocity and Attenuation in Fiber Composites: Theory and Experiments[END_REF] and in particular at high concentration.

Although equation ( 5) cannot be solved directly, an iterative pattern initialized by the result from equation ( 3) can be used to find the solution. The iterative pattern is then stopped when the results converge.

D. Taking different obstacles into account

In the case where there are several sizes and types of scattering particle, it is possible, for both models, to introduce (in equations 3 and 6) the influence of these scattering particles, based on the average of the scattering functions, such that:

( ) ( ) ( ) θ = θ i 0 i , f . n i n f (7)
where i is the index of the different scattering particles, n(i) is the number of particles (i) per unit volume, and ( ) i , f θ is the scattering function of the obstacle i, at an angle θ .

E. Modeled outcome for wave propagation in concretes

In order to observe the evolution dynamics of the two propagation models, the responses of the models are studied, in terms of phase velocity and attenuation, as a function of the longitudinal wave's frequency in media similar to concrete. The scattering particles related to the composition and to the damage are introduced separately, then simultaneously. The maximum order used in scattering calculations (eq. 3 and 6) for n is 20, at this value we didn't observed significant variations in function forms for the studied frequency domain. Variations in the volume fraction of the scattering particle are proposed. A second order polynomial approximation describes the attenuation as a function of frequency (Fig. 4). The measurement protocol used to obtain these values is that shown in Section III.B.

The variation of the attenuation is approximated by a second order polynomial function, which corresponds to the best approximation found for all orders between 1 and 4.

Although this approximation is valid for the domains corresponding to the passband of the sensors (from 200 kHz to 1.2 MHz), the classical range of test frequencies rarely exceeds the value of approximately 500 kHz, which is included by the sensor passband.

Table 1 presents the materials used for the scattering particles into the mortar matrix. Glass density was provided by glass beads manufacturer and glass ultrasonic velocities were chosen near to literature values [START_REF] Gür | Determination of the influence of TiO2 on the elastic properties of a mica based glass ceramic by ultrasonic velocity measurements[END_REF] . The value for the longitudinal velocity has been verified, on several beads of 10 mm diameters, by contact transmission measurement. Expanded polystyrene parameters included in specimens were considered equal to the ones of air according to the manufacturer which ensures that the material contains 98 percent air.

In order to tend towards the real composition of concretes, glass beads were introduced into the mortar matrix. Glass was chosen for its physical, mechanical and acoustic properties which are close to those of the rocks normally found in concrete. The simulation of damage in concrete is implemented by introducing spherical cavities. This damage often occurs through the appearance or evolution of porosities or microcracks in the mortar matrix.

Modeling concrete composition

The numerical application relates to the two studied models, which we note as WT for the Waterman-Truell's model [START_REF] Waterman | Multiple scattering of waves[END_REF] , and GSCM for that of the "Generalized Self Consistent Method" [START_REF] Yang | A Dynamic Generalized Self-Consistent Model for Wave Propagation in Particulate Composites[END_REF] . The introduced media are identical for both models, and correspond to 10%, then 30% of 3 mm diameter spherical glass diffusers in the fine mortar matrix. The maximum chosen values of volume fraction correspond to values similar to those encountered in concrete, for a given obstacle dimension. The dimensions of aggregates in media such as concrete can vary between 0 and 20 mm; the simulation relates to beads with a diameter of 3 mm. The results obtained in terms of phase velocity and attenuation as a function of frequency are shown in Fig. 5.

The introduction of solid inclusions whose acoustic impedance is greater than that of the matrix leads, whatever the model, to an increase in phase velocity and attenuation relative to those obtained for the matrix alone (Fig. 4). The increase in velocity results from the higher velocity in the inclusions, and the increase in attenuation results from the scattering produced by the introduced objects. The greater the volume fraction of scattering objects, the greater the velocity and attenuation values over the studied frequency domain.

We observe different behaviors for the two models, in terms of their dependence on frequency. When the scattering particle volume fraction increases, the differences between the two models become stronger. In terms of velocity analysis, we observe, for both models, a velocity minimum at a low frequency (around 450 to 500 kHz, which corresponds to values of k.a close to 1) for 10% of scattering particles, and a slightly higher value (nearer to 500 kHz, with k.a close to 1) for 30% of scattering particles.

The evolution of the attenuation also reveals different behaviors for the two models, between which the greatest differences are observed at high frequencies. Over the first part of the frequency domain (< 500 kHz, k.a < 1), the models are relatively close to one another. Strong differences are observed between the two models at higher frequencies, with diverging values of attenuation. More the volume fraction of glass inclusions are high, more the observed differences are important. These behaviors confirm those observed by Yang [START_REF] Yang | A Dynamic Generalized Self-Consistent Model for Wave Propagation in Particulate Composites[END_REF] in glass/epoxy composites.

Modeling damage in concretes

In order to model damage in a medium, the 3 mm diameter spherical cavities are introduced into the same fine mortar matrix, with volume densities of 10%, and then 30%.

The value of 30% of air inclusions in a granular medium is high in comparison with the evolution in volume expected when such media are damaged. The results obtained in terms of phase velocity and attenuation, as a function of frequency, are shown in Fig. 6.

Whatever model is used, the introduction of air inclusions with an acoustic impedance lower than that of the matrix leads to a fall in velocity and an increase in attenuation, in comparison with the values obtained for the matrix alone (Fig. 4). The fall in velocity is caused by the slower velocity in the inclusions, and the increase in attenuation is caused by the scattering produced by the introduced obstacles.

The greater the volume fraction of the scattering objects, the greater the changes in velocity and attenuation over the full domain of studied frequencies. Similar curves are obtained for both models. However, the differences between the models increase with increasing scattering particle volume fraction.

We systematically observe a fall in velocity at low frequencies, with a local minimum in the vicinity of 450 to 500 kHz, for 10% of scattering particles (k.a close to 1) and slightly greater (closer to 500 kHz, k.a close to 1), for 30% of scattering particles. The strongest differences between the WT and GSCM models are observed at low frequencies.

For the attenuation, we also note similar behaviors, with transitions (sudden changes in slope) occurring at different frequencies. Like for the case of composition modeling, one can observe increasing differences between the two models with increasing scatterers volume fraction.

Coupling between composition and damage

The coupling of two types of scattering particle is produced by simultaneously introducing 10% of spherical glass scattering particles with a diameter of 3 mm, and 10% of spherical cavities with a diameter of 3 mm, into a fine mortar matrix. This operation allows the coherency of the results given by this simulation, in comparison with the results of the preceding simulations, and the models' respective abilities to simultaneously integrate several types of diffuser, to be evaluated. The results obtained in terms of phase velocity and attenuation as a function of frequency are shown in Fig. 7.

The results obtained clearly show the combined effects of the different separately observed behaviors. The differences between the models are less significant than in the two preceding cases, which can be explained by the behavior inversions depending on the model, with respect to the different types of scattering particle.

As the glass and cavities in this simulation identical in size, the strong evolutions obtained at 500 kHz (k.a close to 1 in both cases) are quite noticeable for both models. The velocities are distributed around that found for the matrix alone. The greatest divergences between the two models are observed in the high frequency domain.

The variations in simulated attenuation are very similar, which can be explained by the inversion of the attenuation behaviors between the models, depending on the types of scattering particle. The attenuation values are of course greater than when each case is treated separately.

III. EXPERIMENTAL VALIDATION STUDY

The experiments we carried out were designed to verify the validity of the wave propagation models, and to evaluate the limits of each model with respect to the types, sizes and densities of scattering particles, as well as in terms of the frequency of the acoustic waves used for the tests.

A. Definition of the test specimens

All of the test specimens we used were made from the same, very fine mortar matrix obtained with a mixture of cement, water and fine sand (grain size < 0.5 mm). This matrix is defined in Table 2, and is referred to as FM (Fine Mortar). Glass beads of different sizes (diameter ∅) were then added to the FM in order to simulate the composition of concretes (FM-G3 and FM-G10), and expanded polystyrene balls (98% air) allowed damage to be simulated by increasing the air scattering volume fraction inside the medium (FM-P3). The 3 mm beads allowed sand grains to be accounted for, and the 10 mm beads were used to represent aggregates. The air beads can simulate damage by increasing the porosity of concretes and/or damage caused by the appearance and evolution of randomly oriented microcracks. The coupling between composition and damage is achieved by using a test specimen containing both types (glass and air) of diffuser (FM-G3-P3).

The fine mortar test specimen (FM) allows the characteristics of the matrix, into which the scattering particles are later introduced, to be determined. The matrix starting values for the models are those measured and presented in Fig. 4. The data concerning the glass beads and expanded polystyrene balls included into the test specimens are shown in Table 1. The expanded polystyrene balls are composed of 98% air, and are considered as equivalent to air inclusions. The test specimens have the shape of a disk, with a diameter of 250 mm and thickness of 50 mm.

B. Experimental measurement protocol

The ultrasound measurement protocol we developed was adapted to the characterization of heterogeneous media. The measurements were carried out using the immersion transmission of a longitudinal wave, which allowed reproducible coupling levels to be obtained between the specimens and the transducers. Three pairs of transducers were used and the total passband found at -12 dB ranged from 200 kHz to 1.2 MHz. The measurement was based on a comparative method between the direct signal transmitted through water and that corresponding to a trajectory through water and the tested specimen (Fig. 8).

The distance between the sensors was the same in both cases, with the specimen being placed in the far field of the emitting transducer. The quantities averaged over the scatterers positions are obtained from averages, for each specimen, on uncorrelated transducers positions. This operation allows the amplitude of the incoherent waves to be reduced, thus allowing the coherent transmitted wave to be separated and analyzed. The broad-band temporal signals obtained are reduced by Fourier analysis, in order to extract their phase and amplitude.

The phase velocity and attenuation of the ultrasonic sensor are then determined, from the thickness of the specimen (e) and the following expressions:

( ) ( ) ( ) ω φ ∆ + ω + ω φ - ω φ ω =

IV. RESULTS AND DISCUSSION

In order to validate the behavior of the studied models during the previous numerical simulations, the present section proposes a comparison between theory and experiment, based on the specimen compositions described in Table 2. Firstly, the results for test specimens with compositions simulating similar concrete formulations (FM-G3 and FM-G10) were studied. Secondly, the results for test specimens simulating damaged concrete (FM-P3), and thirdly results from coupling between two types of scattering particle, based on the (FM-G3 P3) test specimen, are presented. Fig. 9 shows the velocity and attenuation results found for the FM-G3 test specimen made of fine mortar including 34.1% of 3 mm diameter glass beads, with the corresponding models.

The phase velocity results are found to be in nearly perfect agreement between the experiment and the GSCM model, in the observed experimental passband. The maximum differences between the two curves remain within the measurement uncertainties (dotted curves). Dispersion on measured phase velocity at low frequency can be attributed to electrical noise that can occurred for small signal to noise ratio [START_REF] Kalashnikov | Errors and Uncertainties in the Measurement of Ultrasonic Wave Attenuation and Phase Velocity[END_REF] especially at the limit of the passband.

As far as the attenuation is concerned, the general shape of the curves is close to that of the WT model, although the experimental values obtained are quite different to those given by this model. They are close to the GSCM model values in the low frequency part of the passband, and are removed from them at frequencies greater than 700 or 800 kHz (k.a is approximately equal to 1.5). The overestimating behavior of the WT model is classically observed in literature for acrylic spheres in ethylene glycol [START_REF] Vander Meulen | Theoretical and experimental study of the influence of the particle size distribution on acoustic wave properties of strongly inhomogeneous media[END_REF] and for glass spheres in epoxy matrix [START_REF] Layman | The interaction of ultrasound with particulate composites[END_REF] . The good agreement of the GSCM model for high frequencies is also obtained the solid heterogeneous medium [START_REF] Layman | The interaction of ultrasound with particulate composites[END_REF] . Fig. 10 shows the velocity and attenuation results relative to the FM-G10 test specimen composed of fine mortar including 32.1% of 10 mm diameter glass beads, with the corresponding modeled curves. This allows differing behaviors to be observed, as a function of changes in size, and possible limits in useful range of the models, in terms of frequency or scattering domain, to be observed.

The results for phase velocity show that there is a good agreement between the experimental and GSCM model values, between 200 kHz and 700 kHz (k.a is approximately 5), the latter frequency being that beyond which the measured velocity begins to fall. This decrease in velocity then brings the experimental behavior closer to that predicted by the WT model at high frequencies. Theses behaviors can be related to the one observed by Layman and al [START_REF] Layman | The interaction of ultrasound with particulate composites[END_REF] in case of glass spheres in epoxy matrix in which no conclusion can be done between WT and GSCM models for the phase velocity dispersion.

The obtained results for the velocity in the present study show better agreements for GSCM model than WT one when elastic spheres are introduced in elastic matrix.

Concerning the attenuation, the general shape of the curves is close to that of the GSCM model and the values obtained are generally very close to those given by the model.

As previously shown (for 3 mm beads), the experimental values depart from the model at frequencies higher than 400 kHz (k.a values approximately equal to 2.5). The greatest differences observed are 10 Np/m at the highest frequencies (1.2 MHz).

The GSCM model thus appears to more accurately describe the behavior of ultrasonic waves, when the scatterers introduced into the elastic matrices are elastic spherical particles.

In order to evaluate the influence of introducing cavities into the medium, Fig. 11 shows the velocity and attenuation results for the FM-P3 test specimen. It is made of fine mortar, with 30% of 2.84 mm diameter expanded polystyrene balls, together with the results of the corresponding model.

The phase velocity results show good agreement between the experiment and the WT model, over the full frequency domain. The differences between experiment and WT model remain very near to the measurement uncertainties limits (dotted curves). The two local minima predicted by the GSCM model do not appear in the experiment, for which just one minimum can be observed.

The attenuation curves show that the measured values are very close to the WT model, up until 700 kHz (k.a = 1.6), and then become greater than those of the model beyond this frequency. The general shape of the experimental and WT modeled curves are very similar. The attenuation foreseen by the GSCM model is quite different to the measured values. This behavior should be complemented by Fig. 6, which shows that beyond the observed frequencies the curve should continue to increase, up until a turning point at a much higher value. Measured attenuation at high frequency reaches important values and has a dispersive behavior that can be related to the one obtained at low frequency on velocity on Fig. 9.

The WT model appears to have the best performance in describing the behavior of ultrasonic waves when spherical cavities are introduced into elastic matrices.

Finally, Fig. 12 presents the phase velocity and attenuation results corresponding to the FM-G3-P3 test specimen combining two types of diffuser. This specimen comprises fine mortar including 34.1% of 3 mm diameter glass beads and 17.1% of 2.84 mm diameter expanded polystyrene balls.

The results for phase velocity show that there is a good agreement between the shape of the experimental and studied model curves, over the full frequency range.

However, there are significant differences between the theory and the models, in particular for the WT model. The GSCM model is in agreement with the experimental values beyond 700 to 800 kHz (k.a is of the order of 1.5). These observations should be compared with the velocity curves shown in Figs. 9 and11, where it can be seen that the WT model does not fully represent the solid scattering particles, and the GSCM has evolution dynamics which are too strong (in the vicinity of k.a values close to 1), in the case of the introduction of the cavities.

The attenuation computed by the GSCM model is close to the experimental values.

The WT model predicts an attenuation curve similar in shape to that found experimentally, but with higher values than in the experiment. Neither of these two models allows a quantitatively correct description to be found, in agreement with the experimental measurements. Their good qualitative agreement nevertheless confirms the potential of these models in this type of medium.

The different results show the good agreements of GSCM model for the case of elastic inclusions in elastic matrix and good agreements of WT model for the case of cavities in elastic matrix up to about k.a=1.5 and when volume fraction of particles reaches 34%. The frequency limit is sufficient to cover the typical frequency range of ultrasonic testing methods for concrete and the volume fractions of obstacle are representatives of the ones obtained in concrete material.

Combination of the two types of scatterers provides relative agreements but not quantitatively perfect in particular for the low frequencies range. This case highlights the potential of effective medium modeling to reach high volume fraction but also improvements we need to do. The observed differences in results between model and experiment can be due to different sources.

First, small errors in measurement at the limits of the total passband 25 (see Fig. 9 for velocities at low frequencies or Fig. 11 for attenuation at high frequencies) can appear on measured velocities and attenuations. Another parameter is the averaging process on scatterers positions done by uncorrelated spatial averaging which can be quite different or particular because of the nature of statistic observations. Second, all parameters used for the input data (Tables 1 and2) are known with unknown uncertainties and may be inaccurate due to difficulties to obtain representative values in literature (properties of glass and expanded polystyrene), to control manufacture process (real values of volume fractions of inclusions, similar cement matrix for all specimen, …) or to implement good evaluation of this parameters by measurement systems. Third, approximations done in the construction of the different models can be definitely the largest part of the observed disagreements. To comment this point we can focuses on each case of inclusion type (see Fig 9 and Fig. 11).

In the two models, differences can appear due to mismatch of approximation in calculation for the scattering by one obstacle or approximation in the construction scheme of effective medium model. If the way to model scattering by sphere used in WT model is classically employed and is validated with good accuracy, the one used in GSCM model, with transition medium around the sphere, is quite different and cannot be validated separately. This method is coupled with the self-consistent scheme and shows, in past study [START_REF] Huang | Frequency Dependencies on Ultrasonic Wave Velocity and Attenuation in Fiber Composites: Theory and Experiments[END_REF] , interesting improvement for the elastic two dimensional case (scattering by elastic cylinders). For the case of elastic three dimensional spheres, we also observed improvement of results provided by GSCM model in comparison to WT model. This is not true for the case of cavities in elastic matrix. The way used in GSCM model seems to be adapted for elastic inclusions and to be failed for cavities. Maybe the strong acoustic contrast between inclusions and matrix induces too high variations in scattering calculations done in GSCM model.

The two employed models are originally constructed on the basis of fluid media including scatterers with existence of an only one longitudinal wave in the effective medium and the conversion modes are neglected in this multiple scattering theory [START_REF] Conoir | Effective wavenumbers and reflection coefficients for an elastic medium containing random configurations of cylindrical scatterers[END_REF] . So, the energy involved in mode conversions is considered as lost that can induces overestimating of attenuation and also dispersion in phase velocity. If the self-consistent scheme provides new solutions in case of high concentrated elastic scatterers in elastics matrix, taking into account mode conversions using Quasi Crystalline Approximation [START_REF] Lax | Multiple scattering of waves[END_REF][START_REF] Linton | Multiple scattering by multiple spheres: A new proof of the Lloyd-Berry formula for the effective wavenumber[END_REF][START_REF] Conoir | Effective wavenumbers and reflection coefficients for an elastic medium containing random configurations of cylindrical scatterers[END_REF] will surely improve modeling results, like it is already done for two dimensional random composites 30 .

V. CONCLUSION

We have analyzed the longitudinal wave propagation in a heterogeneous medium like concrete including elastic or fluid scattering particles. The evolution of the phase velocity and the attenuation have been computed with two homogenization models (Waterman and Truell -WT, and "Generalized Self Consistent Method" -GSCM) over a wide range of frequencies, as a function of the volume ratio and size of spherical scattering particles. An experimental validation study has allowed each of these models to be qualified with respect to changes in concrete. Both models have shown different potential, with the GSCM appearing to be the most suitable for describing elastic solids containing elastic scattering particles, whereas the WT model allows cavities to be taken into account.

The observed frequency limits are found to be in good agreement, up until k.a values of the minimum order of 1.5 for the velocity and for the attenuation. Both models have shown good potential in representing the two types of scatterers. The key to modeling improvements may lie in the coupling of these two models, by way of an incremental approach to their representation for different types of scattering particles. A first step can be the modeling of concrete composition by elastic inclusions in a cement matrix with GSCM model to obtain equivalent homogeneous medium. Then this equivalent medium can be used like matrix for the second step including spherical cavities with WT model to take into account damage in concrete. Thus, with these two steps, one can reach high volume fraction of obstacles in elastic matrix. This outcome allows the use of these models to be validated in media with very high scattering particles volume fraction, for the prediction of phase velocity and attenuation but neglecting multiple scattering effects between the two types of scatterers. Another way to improve ultrasonic non destructive characterization of heterogeneous media with high particles concentration can be conducted using model which integrates modes conversion. TABLE 1. Properties of the media under consideration

Tables

  The model is based on the knowledge of a mortar matrix ( ) scattering particles have been introduced. The longitudinal ( )propagation velocity characteristics of the matrix are constant.

  the phases and amplitudes of the signals obtained through water, and those obtained through water and the sample. The coefficients coefficients at the water-specimen and specimen-water interfaces respectively. The terms coefficients for phase and amplitude, respectively, of the beam divergence in the far field during transmission of the acoustic waves. These terms are computed from the Thompson model 25 . The measurement uncertainties were evaluated for the full bandwidth, using repeatability tests. This uncertainty varies only slightly as a function of frequency and the type of test specimen. The values retained were the maxima, i.e. ±45 m/s for the phase velocities and ±3 Np/m (values corresponding to a 95.4 percent confidence interval). These values are given on the following figures by dotted curves.

FIG. 1 .FIG. 2 .FIG. 3 .

 123 Figure captionsFIG.1. Real medium, geometric model and equivalent homogeneous medium.

FIG. 5 .FIG. 8 .

 58 FIG. 5. (Color online) Phase velocity and attenuation modeled versus frequency in fine

J.Y. Kim, "Models for wave propagation in two-dimensional random composites: A comparative study", J. Acoust. Soc. Am. 127, 2201-2209, (2010).


Stéphane Gauthier
email: stephane.gauthier@univ-paris1.fr

Optimal Tax Base with Administrative Fixed Costs

Keywords: JEL classification numbers H21 indirect taxation, VAT, tax base, administrative costs

This note characterizes the optimal base for commodity taxation in the presence of administrative fixed costs varying across goods. For low tax rates, the optimal base only comprises commodities whose discouragement index is greater than the ratio of their administrative costs to the tax they yield. An illustration with UK data shows that a category of goods should be taxed only if the revenue generated on this category is at least ten times greater than its administrative fixed cost. The cost imputable to the category of goods taxed at the standard rate would be at most 6 percent of total VAT revenue. The administration cost associated with categories of goods currently tax free could justify exemption.

Introduction

The theory of optimal taxation focuses on the excess burden of taxation as the main source of the social loss caused by taxation. The empirical evidence however suggests that the administrative costs required to collect taxes may be substantial [START_REF] Slemrod | Optimal taxation and optimal tax systems[END_REF]). Most measures of such costs are based on staff salary or equipment costs, but little is known about their true underlying determinants and precise shape [START_REF] Shaw | Administration and compliance[END_REF]). [START_REF] Polinsky | Pigouvian taxation with administrative costs[END_REF], [START_REF] Kaplow | Optimal taxation with costly enforcement and evasion[END_REF] and [START_REF] Mayshar | Taxation with costly administration[END_REF] have derived optimal tax rules when the administrative cost function displays usual continuity and convexity properties with respect to the level of tax rates. Still, as argued by [START_REF] Slemrod | The costs of taxation and the marginal efficiency cost of funds[END_REF] and [START_REF] Alm | What is an optimal tax system?[END_REF], it is likely that this function exhibits significant discontinuities and/or nonconvexities. [START_REF] Yitzhaki | A note on optimal taxation and administrative costs[END_REF] and [START_REF] Wilson | On the optimal tax base for commodity taxation[END_REF] have given early insights into the optimal indirect taxes when these discontinuities are due to idiosyncratic fixed costs varying across commodities. A minimum number of employees may be required for performing the administration of taxed goods. Otherwise, no tax would be recovered because of, e.g., tax evasion and black market operations. The tax authority then faces two alternatives: a category of goods can be either taxed, in which case the tax authority must bear some specific administration fixed cost, or this category can remain tax free. The characterization of the tax base therefore involves a discrete choice between taxation and exemption. In a simple Cobb-Douglas partial equilibrium economy with one representative consumer and uniform commodity taxation, [START_REF] Yitzhaki | A note on optimal taxation and administrative costs[END_REF] has shown that this choice relies on price sensitivity of demand, as is usual in a Ramsey setup, and also on the level of the demand. The tax base actually comprises all the goods with a low enough ratio of administrative cost to the amount of tax they yield. Hereafter this ratio will be referred to as the 'Yitzhaki ratio.' This note generalizes the papers of [START_REF] Yitzhaki | A note on optimal taxation and administrative costs[END_REF] and [START_REF] Wilson | On the optimal tax base for commodity taxation[END_REF] in three directions: (i) optimal tax rates possibly differ across goods; (ii) households have heterogeneous preferences with respect to consumption and labor; (iii) the social planner has a redistributive objective. This generalization is made possible by considering a new problem which does not directly refer to the discrete choice between taxation and exemption, but instead the continuous choice of the probability that a category of similar goods will be subject to taxation. The optimal choice of this probability can be treated by appealing to the usual Lagrangian method, which enables us to compare solutions where the category is either taxed or exempted. This construction yields a general rule for including a good in the tax base. This rule takes a very simple form in the particular case of low tax rates: a good must be included in the tax base only if its associated [START_REF] Mirrlees | Optimal tax theory : A synthesis[END_REF] discouragement index is higher than its Yitzhaki ratio. Thus the influence of equity on the decision to tax only transits through the optimal levels of tax rates: a good whose demand should not be strongly discouraged, because it is consumed by agents whose social value is high, is less likely to be taxed.

An empirical illustration on UK data suggests that a category of goods should be taxed only if the revenue generated on this category is at least ten times greater than its administrative fixed cost. Currently untaxed categories of goods may fail to satisfy this requirement, which could justify exemption. The illustration also gives upper bounds for the ratio of administrative fixed costs to total VAT revenue. The highest bounds apply to the categories 'Household Goods and Services', 'Petrol and Diesel', 'Tobacco' and 'Leisure Services'. The administrative cost imputable to the whole category of goods taxed at the standard rate would be at most 6 per cent of total VAT revenue.

The setup

The economy consists of a tax authority, a continuum of heterogeneous households h ∈ H, and a continuum of categories of consumption goods i ∈ C produced from labor provided by households. The tax authority raises a given amount of revenue R by using linear taxes on consumption goods. In the Ramsey approach the tax authority first sets the tax rates t i applying to goods in category i, and then households choose how much to consume and how much to work, considering as given the various tax rates.

We proceed backward and begin with the study of households' behavior. This behavior will act as a constraint in the problem of the tax authority. The preferences of household h are represented by a utility function which is separable across consumption goods and labor,

C u h i (x h i )dµ i -h , (1)
where µ is a non-atomic measure which captures the importance of the different categories of goods.

The problem of household h consists in choosing a bundle (x h i ) and a labor supply h which maximize her utility (1) subject to the budget constraint

C (1 + t i )x h i dµ i ≤ h . ( 2
)
Consumption goods are produced from labor according to a linear technology normalized so that producer prices equal 1, and thus 1 + t i is the consumer price of goods in category i. The demand function solution to the problem of maximizing (1) subject to (2) is

x h i = ξ h i (t i ) for every i ∈ C. Indirect utility is C u h i (ξ h i (t i )) -(1 + t i )ξ h i (t i ) dµ i ≡ C v h i (t i )dµ i .
The contribution of category i goods to the welfare of household h can therefore be measured by v h i (t i ).

Some of the arguments used in section 3 to characterize the composition of the tax base rely on the assumption that utility is derived from a continuum of goods (see Remark 1). The modeling assumption that goods differ continuously with some characteristics may not fully accord with the formulation (1), since a tax rate placed on good i is likely to cause substitution into almost identical neighboring goods. Section 4 considers a general utility function, with a continuum of consumption goods, allowing for cross-price and income effects. It shows that the simple formulation (1) actually provides the main insights.

The decision to tax

Following [START_REF] Yitzhaki | A note on optimal taxation and administrative costs[END_REF], the tax authority must pay a fixed cost c i when it decides to tax (or subsidize) commodity i. Otherwise, this commodity remains tax free. The tax authority thus chooses both the tax rates applying to taxed commodities and the composition of the tax base.

Whether a commodity should be taxed or exempted is a discrete decision to which the standard Lagrangian method does not directly apply. It can be treated as a continuous decision by proceeding as if it were possible to tax commodities randomly. Suppose that commodity i is taxed at rate t i with probability π i (0 ≤ π i ≤ 1) and exempted otherwise. When this commodity is taxed, household h gets indirect utility v h i (t i ) and pays the tax t i ξ h i (t i ) while the authority bears the administration cost c i . When it is tax free, household h indirect utility is v h i (0). The problem of the tax authority is to select t i and π i for every i ∈ C. Let the measure ν stand for the distribution of households, and let γ h be the social valuation of the welfare of household h. At the optimum, the profile (t i , π i ) maximizes

H γ h   C π i v h i (t i )dµ i + C (1 -π i )v h i (0)dµ i   dν h (3) subject to the budget constraint C π i   t i H ξ h i (t i ) dν h -c i   dµ i ≥ R (λ) (4)
and for every i ∈ C,

π i ≥ 0, (ρ i ) π i ≤ 1. (σ i )
The variables in brackets are the associated Lagrange multipliers.

Using the Lagrangian approach, a first-order (necessary) condition for π i to be a maximum is

H γ h v h i (t i ) -v h i (0) + λ t i ξ h i (t i ) -c i dν h + ρ i -σ i = 0. ( 5
)
In addition, the Kuhn and Tucker exclusion conditions must be satisfied,

ρ i ≥ 0, ρ i π i = 0, (6)
σ i ≥ 0, σ i (1 -π i ) = 0. (7) Let L i (t i , λ) ≡ H γ h v h i (t i )dν h + λt i H ξ h i (t i )dν h
stand for the contribution of commodity i to social welfare (net of its associated administrative costs). The first-order conditions ( 5), ( 6) and ( 7) directly yield the following result:

Proposition 1 Assume that commodity i is taxed (or subsidized) at some arbitrary rate t i when it belongs to the tax base. A necessary condition for commodity i to belong to the tax base is

L i (t i , λ) -L i (0, λ) ≥ λc i , or equiva- lently H β h v h i (t i ) -v h i (0) dν h + t i H ξ h i (t i )dν h ≥ c i , (8)
where β h ≡ γ h /λ is the marginal social valuation of the income of household h.

Similarly, a necessary condition for commodity i to be tax free is

L i (t i , λ)- L i (0, λ) ≤ λc i , or equivalently H β h v h i (t i ) -v h i (0) dν h + t i H ξ h i (t i )dν h ≤ c i . (9)
In the second stage described in Section 2 households will be faced with the situation where commodity i is either exempted (π i = 0) or taxed (π i = 1). The relevant solutions must therefore be such that π i ∈ {0, 1}. There is in fact no loss in focusing attention on this type of solutions in Proposition 1. This is clear when inequality (8) or ( 9) is strict. This is also clear when there is only one non-atomic commodity for which (8) holds at equality, since setting the probability on that commodity to either 0 or 1 can then be done without any further implications.

In the remaining case where there is an atomic group of commodities for which (8) holds at equality, and π i / ∈ {0, 1} for some of these commodities, it is always possible to split this group into two (atomic) subgroups, one consisting of taxed commodities and the other consisting of exempted commodities, such that both the budget constraint and the social objective are unaffected. To see this, suppose that (8) holds at equality for every commodity in an atomic group I ⊆ C, i.e.,

H γ h v h i (t i ) -v h i (0) dν h = -λ   t i H ξ h i (t i )dν h -c i   (10)
for all i ∈ I, µ(I) > 0. Consider any interval I = i inf , i sup ⊆ I, and suppose that π i = {0, 1} for some i ∈ I. Let R I be the amount collected from I. There is then a new profile (π * i ), π * i ∈ {0, 1} for every i ∈ I, such that the profiles (π i ) and (π * i ) yield the same amount R I of collected tax,

I π i   t i H ξ h i (t i ) dν h -c i   dµ i = I π * i   t i H ξ h i (t i ) dν h -c i   dµ i . (11)
Indeed, the function

G(i) = i sup i   t i H ξ h i (t i ) dν h -c i   dµ i is continuous with i, i ∈ i inf , i sup ,
when the measure µ is non-atomic. Since G(i inf ) ≥ R I and G(i sup ) = 0, it follows from the intermediate value theorem that there is at least one i ∈ I such that G(i) = R I . By (10), G(i) is decreasing, so that such a i is in fact unique. Proceeding similarly for every interval in I yields a profile (π * i ), π * i ∈ {0, 1} for every i ∈ I, allowing to collect the same amount of tax from I as the initial profile (π i ). Finally, using ( 10), ( 11), and again (10) shows that the value of the social objective is also the same with both profiles,

I   H γ h v h i (0)dν h + π * i H γ h v h i (t i ) -v h i (0) dν h   dµ i = I   H γ h v h i (0)dν h + π i H γ h v h i (t i ) -v h i (0) dν h   dµ i .
Remark 1 The existence of a solution where π i ∈ {0, 1} for every good relies on the measure µ. When this measure is non-atomic, the previous argument shows that there may exist several solutions where π i ∈ [0, 1] for every good, but there is always one solution where every good is either taxed or exempted.

When the measure µ has atoms, e.g., in the discrete good version of this model, the function G used above is no longer continuous. Then it may not be possible to find a solution π i ∈ {0, 1} for every good.

Proposition 1 characterizes the optimal tax base associated with arbitrary tax rates, not necessarily those chosen by the government to maximize social welfare. An optimal tax rate t * i on commodity i is an extremum of L i (t i , λ). It satisfies the first-order condition

-d * i ≡ H t * i ξ i (t * i ) ∂ξ h i ∂t i (t * i )dν h = -((1 -β) -βφ * i ) , ( 12
)
Fig. 1 The tax base in the Cobb-Douglas case

with β ≡ H β h dν h , ξ i (t * i ) ≡ H ξ h i (t * i )dν h , φ * i ≡ cov β h β , ξ h i (t * i ) ξ i (t * i )
. This is the familiar many-person Ramsey formula. At the optimum, the Mirrlees' (1976) discouragement index d * i equals a constant plus an equity correction varying with the commodity.

In an optimal indirect tax structure, the tax base comprises goods satisfying ( 8), the tax rates (t * i ) satisfy ( 12) and the marginal social cost of public funds λ is determined by the budget constraint (4) of the tax authority.

In order to describe how the composition of the tax base described in Proposition 1 and the Ramsey rule (12) interact, consider for instance the particular case where preferences over consumption goods are represented by a Cobb-Douglas utility function, u h i (x h i ) = α h i log x h i for every h, i (whereas the contribution of leisure to utility remains a linear function of its consumption). Appealing to the Ramsey rule,

H α h i dν h = (1 + t * i ) H β h α h i dν h , inequality (8) becomes t * i -log (1 + t * i ) 1 + t * i ≥ c i H α h i dν h . ( 13
)
In the representative agent case considered by [START_REF] Yitzhaki | A note on optimal taxation and administrative costs[END_REF], α h i = α i whatever h and i are, so that the Ramsey tax rate is uniform and the tax base consists of all the goods whose c i /α i ratio is low enough.

Figure 1 illustrates how the optimal tax base is determined in the presence of heterogeneity across households. The LHS of ( 13 Inequality (13) shows that, in the Cobb-Douglas case, heterogeneity across households and equity considerations influence the composition of the tax base through the Ramsey tax rates only. This property extends to more general preferences satisfying (1) in the empirically plausible configuration of low rates of tax, i.e., when t * i is close enough to 0 for every i ∈ C.

Proposition 2 Assume that commodity i is taxed at a low Ramsey tax rate t * i when it belongs to the tax base. It should be included in the tax base only if its associated discouragement index is greater than its Yitzhaki ratio, i.e.,

|d * i | ≥ c i |t * i | ξ i (0)
.

Proof Appealing to Roy's identity, the first-order Taylor expansion of

L i (t * i , λ) at t * i = 0 yields L i (t * i , λ) L i (0, λ) + λt * i ξ i (0)   1 - H β h ξ h i (0) ξ i (0) dν h   . (14)
The tax rate t * i is given by ( 12), which can be rewritten as

d * i = 1 - H β h ξ h i (t * i ) ξ i (t * i ) dν h . ( 15
)
For low tax rates, t * i 0 for every i ∈ C,

H β h ξ h i (t * i ) ξ i (t * i ) dν h H β h ξ h i (0) ξ i (0) dν h + t * i H ∂ ∂t * i β h ξ h i (t * i ) ξ i (t * i ) t * i =0 dν h , so that (15) yields 1 - H β h ξ h i (0) ξ i (0) dν h d * i + t * i H ∂ ∂t * i β h ξ h i (t * i ) ξ i (t * i ) t * i =0 dν h .
Reintroducing this expression into (14), and neglecting the higher order term in t * 2 i , one gets

L i (t * i , λ) L i (0, λ) + λt * i ξ i (0)d * i . ( 16
)
The result then follows from Proposition 1.

Proposition 2 gives a simple picture of the optimal tax base. Assume for instance that the Ramsey tax rate t * i is positive. Then, commodity i is more likely to be exempted whenever (a) it is costly to administrate (c i is high), (b) it yields a low amount of taxes (t * i ξ i (0) is low), and (c) its demand should not be strongly discouraged (d * i > 0 is low), i.e., the efficiency cost induced by taxation of commodity i is high (β is high) and this good is consumed by households with high social value (φ * i is positive).

An empirical illustration

The method used in Sections 2 and 3 can be adapted to characterize the tax base in the presence of cross-price and income effects. Let the preferences of household h be represented by a general utility function u h (x, ) with usual monotony and convexity properties. When this household is faced with tax rate t i = {0, t * i } on commodity i, her budget constraint is

C (1 + t i )x h i dµ i ≤ h .
The demand function of household h for good i is ξ h i (t) and her labor supply is h (t), where t = (t i ) denotes the whole vector of tax rates. Her indirect utility is v h (t).

As in Section 3 we assume that the tax rate t i is set according to a Bernoulli distribution with parameter π i (0 ≤ π i ≤ 1). Social welfare then can be expressed as

H γ h E t v h (t) dν h (17)
and the social budget constraint is

H E t   C t i ξ h i (t)dµ i   dν h - C π i c i dµ i ≥ R. ( 18
)
The optimal indirect tax structure is a profile (t * i , π * i ) which maximizes (17) subject to (18) and the constraints 0 ≤ π i ≤ 1 for every i ∈ C. When the optimal tax rates are low enough, it is shown in appendix that commodity i should be included in the tax base, i.e., π *

i = 1, only if   1 - H β h ξ h i (0) ξ i (0) dν h   t * i d * i t * i ≥ c i ξ i (0) , ( 19
)
where β h = α h γ h /λ, and α h is the marginal utility of income of household h. This generalizes the inequality given in Proposition 2 for a utility function for which the derived demands are income sensitive and for which there are cross-price effects.

Inequality ( 19) can be used to get bounds for actual administrative costs. Appealing to UK data from the Institute of Fiscal Studies, goods are grouped in twenty homogeneous discrete categories of goods reported in Table 1. Approximating integration to summation, the necessary condition (19) will be applied to these categories of goods. All the goods in a category are taxed at the same rate (except for goods subject to excises). [START_REF] Belan | Optimal grouping of commodities for indirect taxation[END_REF] recovered from these data a marginal social cost of public funds λ equal to 1.11, and a poorly redistributive VAT with α 1 γ 1 = 0.03, α 4 γ 4 = 0.54 and α 5 γ 5 = 0.43, while α h γ h = 0 for all the other deciles of consumption expenditures. The budget shares by decile for these categories then enable us to compute the actual discouragement indices (applying to uncompensated demand) which appear in the LHS of inequality ( 19). They are reported in the second column of Table 1, with the shares of VAT from each category in the total revenue (in the third column).

The quantities in ( 19) are evaluated at t = 0 since the rule ( 19) has been obtained under the assumption that the tax rates are close to 0. On the other hand, the computation of the discouragement indices is made using the actual (observed) VAT structure, with some goods being taxed. This computation is therefore valid provided that the quantities in ( 19) are not too sensitive to changes in the tax rates. With this caveat in mind, each index is found positive and close to 0.1. Hence, by ( 19), a category will be taxed only if the tax it yields is at least ten times greater than its fixed administrative cost.

Existing empirical evidence suggests that total VAT administrative costs, i.e., including variable costs but excluding possible compliance costs, would be around 1 or 2 percent of total VAT revenue in the UK [START_REF] Sandford | Administrative and Compliance Costs of Taxation[END_REF]; [START_REF] Bickley | Value Added Tax: Concepts, Policy Issues, and OECD Experiences[END_REF]; [START_REF] Evans | Studying the studies: An overview of recent research into taxation operating costs[END_REF]). In order to check whether these estimates are consistent with those obtained from (19), inequality ( 19) is expressed in share of total VAT revenue,

t * k ξ k (0) C t * i ξ i (0)dµ i d * k ≥ c k C t * i ξ i (0)dµ i (20)
for every category k. Assuming that the actual (observed) demand is given by ξ i (t * i ) and can be evaluated as if the various tax rates were zero, the observed share of VAT from every taxed category in total VAT revenue appears in the LHS of inequality (20), while the RHS gives the actual ratio of administrative fixed costs to total VAT revenue. The product of the share of VAT from every taxed category in total revenue and the discouragement index of the category thus provides an upper bound for the ratio of administrative fixed costs to total VAT revenue. These upper bounds for the ratio of administrative fixed costs to total VAT revenue are reported in the fourth column of Table 1. Their highest values stand between 1% and 2% of total VAT revenue. They apply to four categories ('Tobacco', 'Petrol and Diesel', 'Household Goods and Services', and 'Leisure Services'). Otherwise the administrative fixed costs seem low, e.g., those relative to 'Public Transport' do not exceed 0.2% of total VAT revenue,

that is 0.002 ≥ c k C t * i ξ i (0)dµ i
for k = 'Public Transport.' An upper bound for the ratio of total administrative fixed cost to total VAT revenue

T c k C t * i ξ i (0)dµ i
can be found by summation of the bounds for the individual taxed goods. The fixed cost imputable to the whole category of goods taxed at the standard rate would be at most 5% and 6% of total revenue. These upper bounds on administrative fixed costs are above the existing estimates of the total administrative costs, and they have approximately the same order of magnitude as these estimates. One may therefore conclude that the results derived from ( 19) are roughly consistent with the existing UK estimates.

Finally one can also get an idea about the costs attached to untaxed categories by considering what would happen if the optimal rate on these categories were the reduced rate of VAT of 5%. The fifth column of Table 1 reports the corresponding bounds for administrative fixed costs, postulating no reaction of demand. These bounds can be viewed as lower bounds for the ratio of administrative costs of untaxed categories to total VAT revenue. These bounds seem low, about 0.1 or 0.2% of total revenue. In comparison with the actual administrative costs found in the literature, it seems therefore possible that these categories should remain tax free, partly because of relatively high fixed administrative costs.

t k d k t k ξ k C t i ξ i dµ i t k ξ k C t i ξ i dµ i d k tk ξ k C ti ξ i dµ i d k Untaxed Meat & Fish 0 0.

Appendix

This appendix provides an explicit derivation of inequality ( 19). Let T be the set of all taxed commodities, i.e., t i = 0 for every i ∈ T , and t i = 0 otherwise (i ∈ C\T ). The optimal tax rate t * k on commodity k is an extremum of the Lagrangian function

H γ h v h dν h + λ H T t i ξ h i -c i dµ i dν h -λR,
where v h and ξ h i are evaluated at (t). The first-order condition in t k writes

d * k ≡ - H T t i ξ k ∂ξ h i ∂t k dµ i dν h = H 1- γ h α h λ ξ h k ξ k dν h
where ξ k is the aggregate consumption of good k, α h is the marginal utility of income for household h, and d * k is the discouragement index for category k. All these quantities are evaluated at (t * ). Note that, unlike the usual formulation, the discouragement indices refer to uncompensated demand (and not compensated demand) in order to be used in the empirical illustration. The optimal tax rate on good k satisfies

d * k = 1 - H β h ξ h k (t) ξ k (t) dν h 1- H β h ξ h k (0) ξ k (0) dν h + T t * i ∂d * k ∂t i t * i =0 dµ i , (21)
where the approximation is obtained by assuming low tax rates. Let t i = {0, t * i } and assume that t i ∼ t * i × B(π i ), where B(π i ) stands for the Bernoulli distribution with parameter π i . In order to derive the first-order condition associated with the inclusion of commodity k in the tax base, we let t = (t k , t -k ) where t -k is the set of all the tax rates different from t k . The social welfare becomes

π k E t -k   H γ h v h (t * k , t -k )dν h   + (1 -π k )E t -k   H γ h v h (0, t -k )dν h   ,
while the total receipt from commodity taxation is

π k E t -k   t * k H ξ h k (t * k , t -k )dµ k dν h + H C\{k} t i ξ h i (t * k , t -k )dµ i dν h    + (1 -π k )E t -k    H C\{k} t i ξ h i (0, t -k )dµ i dν h    .
Let also

L k (t * k , λ) = H γ h E t -k v h (t * k , t -k ) dν h + λE t -k    H t * k ξ h k (t * k , t -k )dµ k dν h + H C\{k} t i ξ h i (t * k , t -k )dµ i dν h    .
By the same argument as the one used to get Proposition 1, commodity k should be included in the tax base only if

L k (t * k , λ) -L k (0, λ) -λc k ≥ 0. ( 22
)
For low tax rates, i.e., t * i close to 0, Roy's identity yields

v h (t * k , t -k ) v h - T t * i α h ξ h i dµ i
where quantities are evaluated at (t * k , t -k ) = (0, 0). Therefore, at this point,

L k (t * k , λ) H γ h   v h - T t * i α h ξ h i dµ i   dν h + λ H T t * i ξ h i dµ i dν h .
Similarly,

L k (0, λ) H γ h   v h - T \{k} t * i α h ξ h i dµ i    dν h + λ H T \{k} t * i ξ h i dµ i dν h .
The last two expressions only differ with respect to the inclusion of commodity k in the tax base. This commodity belongs to the tax base T in L k (t * k , λ), but it is not in T (in fact it belongs to T \{k}) in L k (0, λ). It follows that

L k (t * k , λ) -L k (0, λ) H λ -γ h α h t * k ξ h k dν h = λt * k ξ k   1 - H β h ξ h k ξ k dν h   .
Appealing to (21), this approximation rewrites

L k (t * k , λ) -L k (0, λ) λt * k ξ k   d * k - T t * i ∂d * k ∂t i dµ i   λt * k ξ k d * k (23)
where the last approximation is obtained by neglecting the higher order terms in t * k t * i . Inequality (19) then comes from ( 22) and ( 23).

  ) is the U-shaped function of the tax rate, with a global minimum at t * i = 0. When satisfied at equality, (13) typically yields two thresholds t inf

	i i ≤ t sup and t sup i ≤ t * i , and it should enter the tax base (t inf i < 0 < t sup i ). Commodity i ≤ t inf i should be tax free only if t inf i only if either t * i ≥ t sup i or t * i .

Table 1

 1 UK VAT base with administrative fixed costs

			0899	0	-	0.0023
	Bread & Cereals	0	0.0901	0	-	0.0007
	Dairy	0	0.0900	0	-	0.0015
	Tea & Coffee	0	0.0901	0	-	0.0004
	Fruits and Vegetables	0	0.0902	0	-	0.0014
	Other Untax. Foods	0	0.0903	0	-	0.001
	Books & Newspapers	0	0.0905	0	-	0.0006
	Children's Clothing	0	0.0907	0	-	0.0003
	Reduced rate
	Domestic Fuels	0.05	0.0900	0.0356	0.0032	-
	Excises					-
	Tobacco	0.791 0.0890	0.1379	0.012	-
	Beer	0.28	0.0904	0.0634	0.0057	-
	Petrol and Diesel	0.732 0.0908	0.1465	0.0133	-
	Wine & Spirits	0.558 0.0913	0.0538	0.0049	-
	Standard rate					-
	Public Transport	0.175 0.0903	0.0211	0.0019	-
	Standard VAT Food	0.175 0.0904	0.0381	0.0034	-
	Food out	0.175 0.0919	0.0590	0.0054	-
	Household G & S	0.175 0.0912	0.2005	0.0182	-
	Adult Clothing	0.175 0.0914	0.0806	0.0073	-
	Leisure Goods	0.175 0.0926	0.0508	0.0047	-
	Leisure Services	0.175 0.0931	0.1101	0.010	-

Acknowledgements I thank Gaël Giraud, Guy Laroque, Bertrand Wigniolle, and Shlomo Yitzhaki for thoughtful comments on the paper. I also thank Eckhard Janeba, editor of this Journal, and two anonymous referees. I am particularly grateful to one reviewer for carefully reading my manuscript, providing guidance and constructive comments which helped me to improve the paper. The usual disclaimers apply.


Christian Obrecht
email: christian.obrecht@insa-lyon.fr

Frédéric Kuznik

Bernard Tourancheau

Jean-Jacques Roux

The TheLMA project: Multi-GPU Implementation of the Lattice Boltzmann Method

In this paper, we describe the implementation of a multi-GPU fluid flow solver based on the lattice Boltzmann method (LBM). The LBM is a novel approach in computational fluid dynamics, with numerous interesting features from a computational, numerical, and physical standpoint. Our program is based on CUDA and uses POSIX threads to manage multiple computation devices. Using recently released hardware, our solver may therefore run eight GPUs in parallel, which allows to perform simulations at a rather large scale. Performance and scalability are excellent, the speedup over sequential implementations being at least of two orders of magnitude. In addition, we discuss tiling and communication issues for present and forthcoming implementations.

Introduction

First introduced by McNamara and Zanetti in 1988 [START_REF] Mcnamara | Use of the Boltzmann Equation to Simulate Lattice-Gas Automata[END_REF], the lattice Boltzmann method (LBM) is an increasingly popular alternative to classic CFD methods. From a numerical and physical standpoint, the LBM possesses numerous convenient features, such as explicitness or ability to deal with complex geometries. Moreover, this novel approach is especially well-fitted for high performance CFD applications because of its inherent parallelism.

Since the advent of the CUDA technology, the use of GPUs in scientific computing becomes more and more widespread [START_REF] Agullo | Numerical linear algebra on emerging architectures: The PLASMA and MAGMA projects[END_REF][START_REF] Dongarra | Exploring new architectures in accelerating CFD for Air Force applications[END_REF]. Several successful single-GPU implementations of the LBM using CUDA were reported during the past years [START_REF] Kuznik | LBM Based Flow Simulation Using GPU Computing Processor[END_REF][START_REF] Tölke | TeraFLOP computing on a desktop PC with GPUs for 3D CFD[END_REF]. Nevertheless, these works are of moderate practical impact since memory in existing computing devices is too scarce for large scale simulations. Multi-GPU implementations of the LBM are still at an early stage of development and reported performance is below what is expected from such hardware [START_REF] Riegel | Implementation of a Lattice-Boltzmann method for numerical fluid mechanics using the nVIDIA CUDA technology[END_REF].

Recently released motherboards are able to handle up to eight computing devices. We therefore chose to develop a POSIX thread based multi-GPU LBM solver, as a first step towards a distributed version. The structure of the paper is as follows. We first briefly present the LBM from an algorithmic perspective, and summarize the principles of LBM implementation using CUDA. Then, we describe the implementation of our solver and discuss performance, scalability, and tiling issues.

To conclude, we study inter-GPU communication and propose a performance model.

LBM Flow Solvers

Although first introduced as an extension to the lattice gas automata method [START_REF] Frisch | Lattice-Gas Automata for the Navier-Stokes Equation[END_REF], the LBM is nowadays more commonly interpreted as a discrete numerical procedure to solve the Boltzmann transport equation:

∂ t f + ξ • ∇ x f + F • ∇ ξ f = Ω(f ) (1)
where f (x, ξ, t) describes the evolution in time t of the distribution of one particle in phase space (i.e. position x and particle velocity ξ), F is the external force field, and Ω is the collision operator.

The distribution function f is linked to the density ρ and the velocity u of the fluid by:

ρ = f dξ (2) ρu = f ξdξ (3)
For LBM, Eq. 1 is usually discretized on a uniform orthogonal grid of mesh size δx with constant time steps δt. The particle velocity space is likewise replaced by a set of N + 1 velocities ξ i with ξ 0 = 0, chosen such as vectors δt ξ i link any given node to some of its neighbors. 1 Figure 1 shows the D3Q19 stencil we chose for our implementation. This stencil is a good trade-off between geometric isotropy and complexity of the numerical scheme for the three-dimensional case. The discrete counterpart of the distribution function f is a set of N + 1 particle population functions f i corresponding to the finite set of particle velocities. In a null external force field, Eq. 1 writes:

f i (x + δt ξ i , t + δt) -f i (x, t) = Ω f i (x, t) (4)
where a i denotes the transpose of (a 0 , . . . , a N ). Density and velocity of the fluid are given by:

ρ = N i=0 f i (5) ρu = N i=0 f i ξ i (6)
The description of possible collision operators is beyond the scope of this contribution. Nonetheless, it should be noted from a mathematical perspective that the collision operators are usually linear in f i and quadratic in ρ and u, as for the multiple-relaxation-time model we chose to implement [START_REF] Humières | Generalized lattice-Boltzmann equations[END_REF][START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]. Equation 4 naturally breaks in two elementary steps:

fi (x, t) = f i (x, t) + Ω f i (x, t) (7)
f i (x + δt ξ i , t + δt) = fi (x, t) (8)
Equation 7 describes the collision step in which an updated particle distribution is computed, and Eq. 8 describes the propagation step in which the updated particle populations are transferred to the neighboring nodes. From an algorithmic standpoint, a LBM flow solvers is therefore outlined by the following pseudo-code:

end for 10. end for

To conclude, it should be noted that parallelization of the LBM is rather straightforward, ensuring global synchronization at each time step being the only constraint.

CUDA Implementation Guidelines

When using CUDA computing devices for LBM simulations, the particle distribution has to be stored in device memory, for obvious performance reason. The simplest strategy to ensure global synchronization is to launch one kernel at each time step. Moreover, due to the lack of control over memory transaction scheduling, two consecutive particle distributions have to be retained at once. On a Tesla C1060, with a D3Q19 stencil in single precision, the maximum number of handleable nodes is therefore about 2.83 × 10 7 .

The LBM being a data-parallel procedure, CUDA implementations usually assign one thread per node. This option allows to take advantage of the massive parallelism of the targeted architecture. The execution grid has to obey the following constraints and limitations:

• The grid may only have one or two dimensions.

• Blocks may have up to three dimensions but the size of the blocks is limited by the amount of available registers per streaming multiprocessor (SM).

• The size of the blocks should be a multiple of the warp size, i.e. 32 on existing hardware, the warp being the minimum execution unit.

A simple and efficient execution grid layout consists in a two-dimensional grid of one-dimensional blocks. For cubic cavities, taking the previously mentioned device memory limitation into account, this solution leads to blocks of size up to 256. Thus, there are at least 64 available registers per node, which proves to be sufficient for most three-dimensional collision models.

As for CPU implementation of the LBM, the particle distribution is stored in a multi-dimensional array. Global memory transactions are issued by half-warps [START_REF]Compute Unified Device Architecture Programming Guide version 3[END_REF]. Hence, in order to enable coalesced data transfer when using the former execution grid layout, the minor dimension of the array should be the spatial dimension associated to the blocks. The disposition of the remaining dimensions, including the velocity index may be tuned so as to minimize TLB misses [START_REF] Papadopoulou | Micro-benchmarking the GT200 GPU[END_REF]. However, this data layout does not ensure optimal global memory transactions. Propagation of particle populations along the blocks' direction yields one-unit shifts in addresses and therefore causes misaligned memory accesses as illustrated in Fig. 2 (in the two-dimensional case for the sake of simplicity).

Figure 2: Misalignment issue

For CUDA devices of compute capability 1.0 or 1.1, misaligned memory transactions may have dramatic impact on performance, since sixteen accesses are issued in lieu of a single one. Performing propagation within the blocks using the shared memory as described in [START_REF] Tölke | Implementation of a Lattice Boltzmann kernel using the Compute Unified Device Architecture developed by nVIDIA[END_REF] is an efficient way to avoid this issue. For compute capability 1.2 or 1.3, however, a misaligned memory access is carried out in as few 128 B, 64 B, or 32 B transactions as possible. An alternative approach for this kind of hardware consists in using in-place propagation instead of out-of-place propagation. Figure 3 illustrates the two propagation scheme. With the GT200, misaligned reads are far less expensive than misaligned writes [START_REF] Obrecht | Global Memory Access Modelling for Efficient Implementation of the LBM on GPUs[END_REF]. Thus, the in-place propagation approach in only slightly less efficient than the shared memory method for simulations in three dimensions [START_REF] Obrecht | A new approach to the lattice Boltzmann method for graphics processing units[END_REF], while being simpler to implement and exerting less pressure on hardware.

Multi-GPU Implementation of the LBM

As usual in parallel computing, optimal performance for multi-GPU applications requires efficient overlapping of computation and communication. In our implementation of the LBM, we used the zero-copy feature provided by the CUDA technology which allows GPUs to directly access locked CPU memory pages. As outlined by Fig. 4, each interface between sub-domains is associated to four buffers: two for in-coming populations and two for out-going populations. Pointers are switched at each time step.

Figure 4: Inter-GPU communication scheme Using CUDA streams is another possible way to perform inter-GPU communication. Yet, the zero-copy feature proves to be the most convenient solution since communication is taken care of in GPU code instead of CPU code. As for global memory accesses, zero-copy transactions should be coalescent in order to preserve performance. To enable coalescence, the blocks have to be parallel to the interfaces of the sub-domains. For the sake of simplicity, we decided to use one-dimensional blocks for our multi-GPU implementation as in our single-GPU implementation. Because of the targeted hardware, our solver handles at most eight sub-domains. We therefore chose to split the cavity in rectangular cuboids along one direction orthogonal to the blocks. This direction is selected such as to minimize the amount of data to exchange. Nevertheless, for large cavities, the dimension in the direction of the blocks may exceed the maximum block size. Thus, the simple grid layout described in section 3 may not be used in general. We propose instead, to use blocks of size 2 n with a two-dimensional grid of size (ℓ

x × ℓ y × ℓ z × 2 -m ) × (2 m-n ),
where ℓ x , ℓ y , and ℓ z are the dimensions of the cavity, n and m are adjustable parameters. The retrieval of the coordinates is done using a code equivalent to the following:

w = blockIdx.x<<m | blockIdx.y<<n | threadIdx.x; x = w % lX; y = (w/lX) % lY; z = w/(lX*lY);
The proposed grid layout leads to shuffle the blocks' schedule which tends to reduce partition camping [START_REF] Ruetsch | Optimizing matrix transpose in CUDA[END_REF] and allows efficient communication-computation overlapping as shown in section 5. The optimal values for n and m, which were determined empirically, are n = 7 and m = 15.

From a software engineering perspective, the CUDA technology is a great improvement over the early days of GPGPU [START_REF] Fan | GPU cluster for high performance computing[END_REF]. Yet, there are still some important limitations compared to usual software development. The inability of the CUDA compilation tool-chain to link several GPU binaries for instance, makes difficult to follow an incremental, library oriented approach. To improve code reusability, we therefore developed the TheLMA framework [START_REF]Thermal LBM on Many-core Architectures[END_REF]. TheLMA stands for Thermal LBM on Many-core Architectures, thermal flow simulations being our main topic of interest. Figure 5 outlines the structure of the framework. The TheLMA framework mainly consists in two parts. The first part is a set of C source files which provides various utility functions such as command line parsing, POSIX threads setup, simulation statistics, and data output in various formats. The second part is a set of CUDA source files which are included in the thelma.cu file at compile time. The later is basically a container providing some general macro definitions.

The initialisation module mainly creates CUDA contexts in order to assign a GPU to each POSIX thread. The geometry module is responsible for setting up the boundary conditions and any obstacle that may lay inside the cavity. Permitted velocity directions for each node, as well as specific boundary conditions, are encoded using bit fields. The computation module contains the core kernel of the simulation which performs both collision and propagation. Last, the result module retrieves the macroscopic variables of the fluid.

The source organisation we describe is meant to ease subsequent code changing. For instance, the implementation of a different collision model or propagation scheme would mainly require the modification of the computation module. Likewise, the setting up of a given simulation configuration principally involves changes in the geometry module.

Performance and Scalability

To validate our code, we implemented the lid-driven cubic cavity test case in which a constant velocity is imposed at the top lid, whereas the other walls have null velocity boundary condition. We chose to assign the x direction to the blocks and to split the cavity along the z direction. Performance for LBM solvers is usually given in MLUPS (Million Lattice node Updates Per Second). We measured performance in single precision for a 192 3 and a 256 3 lattice on a Tyan B7015 server with 8 Tesla C1060 computing devices. The former lattice size, being a multiple of 24, allows to run the solver using 2, 3, 4, 6, or 8 GPUs with balanced sub-domains. Figure 6 shows our measurement for both lattices.

The maximum achieved performance is about 2,150 MLUPS obtained on the 192 3 lattice using 6 GPUs. To set a comparison, this result is comparable to the one obtained with the acknowledged Palabos code in double precision on a Blue Gene/P computer using 4,096 cores [START_REF] Latt | Palabos Benchmarks (3D Lid-driven Cavity on Blue Gene/P)[END_REF]. Scalability is excellent in all cases but one, with no less than 91% parallelization efficiency. Tables 1 and2 give the inter-GPU data exchange throughput required to achieve optimal scalability.

With the Tyan S7015 motherboard of our server, the bandwidthTest program that comes with the CUDA development kit gives 5.72 GB/s host to device and 3.44 GB/s device to host maximum cumulative throughput. These values are obtained using the copy functions provided by the CUDA API with pinned memory and not zero-copy transactions as in our program. The communications being symmetric in our case, we may use the mean of these two values, i.e. 4.58 GB/s, as a lower bound for sustained data exchange between GPUs and main memory. Yet, we see that, except in the case of 8 GPUs on a 192 3 lattice, the required throughput is at most comparable to this bound and is not likely to overflow the capacity of the PCI-E links. We can furthermore conclude that our data access pattern generally enables excellent communication-computation overlapping.

Considering the unfavorable case, we deduce that the number of interfaces with respect to the size of the cavity is too large to take advantage of surface to volume effects. This naturally leads to the question of a less simplistic tiling of the cavity than the one we adopted. Yet, dividing a cubic cavity in eight identical cubic tiles would lead to only three interfaces, but would yield non-coalesced, i.e. serialized, zero-copy transactions for one of the interfaces. It is possible with simple code modifications, like splitting the cavity along the x direction instead of the z direction, to evaluate the impact of serialized inter-GPU communication. Figure 7 displays the performance of our code on a 192 3 lattice after such a transformation. We can see the overwhelming impact of serialized communication on performance, which is obviously communication bound. We may therefore conclude that gaining flexibility in tiling would most likely require to adopt a more elaborate execution grid layout.

Inter-GPU Communication

In order to get better insight into inter-GPU communication, we implemented a benchmark program based on a stripped-off version of our LBM solver. The purpose of this software is to evaluate the achievable sustained throughput when using zero-copy transactions over several GPUs. Data to exchange is formed of k two-dimensional L × L arrays of 32-bit words. The communication graph may be either circular or linear. Whereas linear graphs correspond to our solver, circular graph are useful to simulate balanced communication which is likely to occur with multi-node implementations.

The communication graph is specified by an ordered list of involved devices. In our tests, the S7015 being a two-socket motherboard, we used this parameter to obtain optimal balancing between the two northbridges. In addition, the data arrays may optionally be transposed at each access which causes serialization of the zero-copy transactions. Figure 8 shows the obtained throughput averaged over 50,000 iterations for circular graphs, linear graphs, and linear graphs with transposition. In order for our benchmark program to exchange an amount of data comparable to our solver, the values of k and L were set to k = 5 and L = 192. Yet, extensive tests have shown that these parameters are of negligible impact over the measurements.

The variations of throughput for circular and linear graphs are similar, which is natural insofar the configurations with respect of the number of devices are identical. The throughput in the

P = ℓ 3 × T i 2(i -1)D × ℓ 2 = T i × ℓ 2(i -1)D ( 9
)
where i is the number of devices, T i is the corresponding throughput (in MB/s), D is the amount of out-going and in-coming data by node (in bytes), and ℓ is the dimension of the cavity (in nodes).

In the case of 8 GPUs on a 192 3 lattice, the estimated performance is 1,977 MLUPS instead of 1,870 MLUPS, i.e. a 6% relative error. For the serialized communication version, Fig. 9 displays the comparison between the corresponding estimation and the actual performance. The estimation is in good accordance with the measured values, the relative error being at most 15%, which corroborates our assumption.

Conclusion

In this contribution, we describe a multi-GPU implementation of the LBM capable of handling up to eight CUDA enabled computing devices. Performance is comparable with efficient parallel implementations on up-to-date homogeneous supercomputers and is at least two orders of magnitude higher than with optimized sequential code. We propose an execution grid layout providing excellent computation-communication overlapping and an efficient inter-GPU communication scheme. Though simple, this scheme yields excellent scalability with nearly optimal parallelization efficiency in most cases.

Our solver is implemented over the TheLMA framework, which aims at improving code reusability. A thermal LBM solver also based on TheLMA is currently under development. In the present version, management of multiple CUDA contexts is based on POSIX threads. Yet, we plan to extend our framework to more generic parallelization interfaces (e.g. MPI) in order to make possible distributed implementations. This extension requires in conjunction more elaborate data exchange procedures so as to gain flexibility in domain decomposition.

In addition, we study inter-GPU communication using a specific benchmark program. The obtained results allowed us to express a rather accurate performance model for the cases where our application is communication bound. This tool reveals moreover that data throughput depends to a certain extent on hardware locality. An extended version, using for instance the Portable Hardware Locality (hwloc) API [START_REF] Broquedis | hwloc: a generic framework for managing hardware affinities in HPC applications[END_REF], would help performing further investigation and could be a first step for adding dynamic auto-tuning in our framework.

Figure 1 :

 1 Figure 1: D3Q19 stencil

  (a) Out-of-place propagation (b) In-place propagation

Figure 3 :

 3 Figure 3: Propagation schemes

Figure 5 :

 5 Figure 5: The TheLMA framework

Figure 6 :

 6 Figure 6: Performance and scalability on a cubic lattice

Figure 7 :

 7 Figure 7: Performance for coalesced and serialized communication

Figure 8 :

 8 Figure 8: Inter-GPU communication throughput

Figure 9 :

 9 Figure 9: Performance model for serialized communication

Table 1 :

 1 Required throughput at 100% efficiency on a 192 3 lattice

	Number of GPUs	1	2	3	4	6	8
	Kernel duration (ms) 18.29 9.14 6.10 4.57 3.05	2.29
	Exchanged data (MB)	0	2.95 5.90 8.85 14.75 20.64
	Required throughput (GB/s)	0	0.32 0.97 1.93 4.84	9.03
	Number of GPUs	1	2	4	8
	Kernel duration (ms) 84.31 42.15 21.08 10.54
	Exchanged data (MB)	0	5.24 15.73 36.70
	Required throughput (GB/s)	0	0.12	0.75	3.48

Table 2 :

 2 Required throughput at 100% efficiency on a 256 3 lattice

Commonly, only nearest neighbors are linked.