P Xavier

A Ghani

D Mejia

M Miguel-Brebion

M Bauerheim

L Selle

T Poinsot

Experimental and numerical investigation of flames stabilised behind rotating cylinders: interaction of flames with a moving wall

Keywords: combustion, flames, laminar reacting flows Mirror Cylinder Glass balls Seeding

Steady methane/air laminar premixed flames stabilised on a cylindrical bluff body subjected to a continuous rotation are analysed using joint direct numerical simulations (DNS) and experiments. DNS are carried out using a 19 species scheme for methane/air combustion and a lumped model to predict the cylinder temperature. Rotation of the cylinder induces a symmetry breaking of the flow, and leads to two distinct flame branches in the wake of the cylinder. DNS are validated against experiments in terms of flame topologies and velocity fields. DNS are then used to analyse flame structures and thermal effects. The location and structure of the two flames are differently modified by rotation and heat transfer: a superadiabatic flame branch stabilises close to the hot cylinder and burns preheated fresh gases while a subadiabatic branch is quenched over a large zone and anchors far downstream of the cylinder. Local flame structures are shown to be controlled to first order by the local enthalpy defect or excess due to heat transfer between the cylinder and the flow. An analysis of the local wall heat flux around the cylinder shows that, for low rotation speeds, the superadiabatic flame branch contributes to wall heat fluxes that considerably exceed typical values found for classical flame/wall interactions. However, for high rotation speeds, fluxes decrease because the cylinder is surrounded by a layer of burned gases that dilute incoming reactants and shield it from the flame.

Introduction

The stabilisation of premixed flames in the wake of obstacles is a key phenomenon in the field of combustion. Bluff bodies are commonly used as 'flame holders' to create a low-pressure recirculation zone that allows the flame to anchor in a low-speed region [START_REF] Williams | Flame stabilization and propagation in high velocity gas streams[END_REF][START_REF] Longwell | Flame stabilization by bluff bodies and turbulent flames in ducts[END_REF][START_REF] Williams | Some properties of rod stabilized flames of homogeneous gas mixtures[END_REF][START_REF] Penner | Recent studies on flame stabilization of premixed turbulent gases[END_REF] where incoming fresh gases mix with hot recirculated gases and ignite [START_REF] Zukoski | Experiments concerning the mechanism of flame blowoff from bluff bodies[END_REF][START_REF] Chen | A comparison of bluff-body and swirl-stabilized flames[END_REF][START_REF] Correa | Measurements and modeling of a bluff body stabilized flame[END_REF][START_REF] Masri | The structure of the recirculation zone of a bluff-body combustor[END_REF][START_REF] Smith | LES blowout analysis of premixed flow past V-gutter flameholder[END_REF][START_REF] Nair | Near-blowoff dynamics of a bluff-body stabilized flames[END_REF].

Various bluff-body geometries can be used: cylinder, triangle, square, closed cone, closed and open V-gutter, plate, etc. [START_REF] Zukoski | The role of wake transition in the process of flame stabilization on bluff bodies[END_REF][START_REF] Ballal | Weak extinction limits of turbulent flowing mixtures[END_REF][START_REF] Sanquer | Some specific characteristics of turbulence in the reactive wakes of bluff bodies[END_REF]. The fundamental mechanisms involved in flame stabilisation are numerous. Heat transfer between flame holders and the flow, low-temperature chemistry and recirculating gases influence the flame root location (where chemical reactions begin) and determine the ability of a flame to be stabilised [START_REF] Plee | Characteristic time correlation for lean blow off of bluffbody stabilized flames[END_REF]Glassman 1996;[START_REF] Shanbhogue | Lean blowoff of bluff body stabilized flames: scaling and dynamics[END_REF][START_REF] Chen | Flame lift-off and stabilization mechanisms of nonpremixed jet flames on a bluff-body burner[END_REF][START_REF] Kedia | Impact of flame-wall interaction on premixed flame dynamics and transfer function characteristics[END_REF][START_REF] Mejia | Wall-temperature effects on flame response to acoustic oscillations[END_REF].

Since flame holders are usually obstacles placed in the flow, they also trigger hydrodynamic instabilities that can lead to situations where the flame is 'stabilised' (i.e. anchored on the flame holder) but 'unstable' (i.e. showing strong fluctuations of heat-release rate coupled with the unsteady flow created by the flame holder). One specific example of such instabilities is thermoacoustic modes where fluctuations in heat-release rate lead to self-excited instabilities due to a constructive coupling between acoustics and combustion [START_REF] Rayleigh | The explanation of certain acoustic phenomena[END_REF][START_REF] Culick | Combustion instabilities in liquid-fueled propulsion systems[END_REF][START_REF] Candel | Combustion dynamics and control: progress and challenges[END_REF][START_REF] Lieuwen | Combustion instabilities in gas turbine engines. Operational experience, fundamental mechanisms and modeling[END_REF]. These instabilities can restrict the range of operating conditions and yield unacceptable pollutant emissions [START_REF] Rhee | Dynamic behavior of premixed turbulent v-flame[END_REF][START_REF] Huang | Dynamics and stability of lean-premixed swirl-stabilized combustion[END_REF][START_REF] Stohr | Experimental study of vortex-flame interaction in a gas turbine model combustor[END_REF][START_REF] Kwong | Coupled dynamics of lift-off and precessing vortex core formation in swirl flames[END_REF]. Even if hydrodynamic modes are not the sole driver of thermoacoustic instabilities, their role in these instability loops has been largely documented [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF][START_REF] Lieuwen | Unsteady Combustor Physics[END_REF][START_REF] Terhaar | Key parameters governing the precessing vortex core in reacting flows: an experimental and analytical study[END_REF]. For example, vortex generation behind obstacles is one typical mechanism leading to thermoacoustic activity [START_REF] Poinsot | Vortexdriven acoustically coupled combustion instabilities[END_REF][START_REF] Ducruix | Combustion dynamics and instabilities: elementary coupling and driving mechanisms[END_REF][START_REF] Ghani | LES of longitudinal and transverse self-excited combustion instabilities in a bluff-body stabilized turbulent premixed flame[END_REF] which has been identified for a very long time [START_REF] Reynst | Pulsating Combustion[END_REF]. Therefore, our capacity to mitigate these hydrodynamic instabilities is of great importance as it could improve flame stabilisation.

One method to suppress or control these hydrodynamic instabilities is to implement an active or passive control to modify the flow field: micro-jets, vibrating/moving parts, plasma or morphing actuators have been reported [START_REF] Viets | Boundary layer control by unsteady vortex generation[END_REF][START_REF] Gelzer | synthetic jets[END_REF][START_REF] Moreau | Airflow control by non-thermal plasma actuators[END_REF][START_REF] Dong | Elimination of vortex streets in bluff-body flows[END_REF][START_REF] Godoy-Diana | A model for the symmetry breaking of the reverse Benard von Karman vortex street produced by a flapping foil[END_REF][START_REF] Cattafesta | Actuators for active flow control[END_REF][START_REF] Chinaud | Trailing-edge dynamics and morphing of a deformable flat plate at high Reynolds number by time-resolved PIV[END_REF]. Controlling flows around objects has been extensively used in the aerodynamic community to reduce drag or to modify the global flow topology [START_REF] Roshko | Perspectives on bluff body aerodynamics[END_REF][START_REF] Leweke | Stability of vortex shedding modes in the wake of a ring at low Reynolds numbers[END_REF][START_REF] Mittal | Generation of streamwise vortical structures in bluff body wakes[END_REF][START_REF] Cimbala | Large structure in the far wakes of twodimensional bluff bodies[END_REF][START_REF] Prasad | The instability of the shear layer separating from a bluff body[END_REF].

A canonical example to illustrate these hydrodynamic flow instabilities is a laminar flow around a cylinder (Re 2000, based on the cylinder diameter). Kelvin-Helmholtz (KH) instabilities develop in shear layers where velocity gradients induce coherent structures whereas Bénard-von Kármán (BVK) instabilities lead to an asymmetric vortex shedding in the wake of the obstacle [START_REF] Cantwell | An experimental study of entrainment and transport in the turbulent near wake of a circular cylinder[END_REF][START_REF] Monkewitz | A note on vortex shedding from axisymmetric bluff bodies[END_REF][START_REF] Plaschko | Periodic flow in the near wake of straight circular cylinders[END_REF][START_REF] Kelso | An experimental study of round jets in cross-flow[END_REF]. The transition to observe such instabilities does not depend only on the Reynolds number Re [START_REF] Rao | A universal three-dimensional instability of the wakes of two-dimensional bluff bodies[END_REF], but also on the blockage ratio determined by the ratio of the cylinder diameter to the width of the flow passage [START_REF] Sahin | A numerical investigation of wall effects up to high blockage ratios on two-dimensional flow past a confined circular cylinder[END_REF], or the position of the cylinder with respect to walls [START_REF] Rao | Three-dimensionality in the wake of a rotating cylinder in a uniform flow[END_REF]. One method to cancel the KH or BVK instabilities is to rotate the cylinder along its axis. Rotation can be used for symmetry breaking and wake control. Depending on the rotation rate (defined as the ratio of the cylinder boundary velocity to the incoming velocity), the usual von Kármán vortex street observed for a fixed cylinder can be suppressed [START_REF] Badr | Unsteady flow past a rotating circular cylinder at Reynolds numbers 10 3 and 10 4[END_REF][START_REF] Schumm | Self-excited oscillations in the wake of two-dimensional bluff bodies and their control[END_REF][START_REF] Modi | Moving surface boundary-layer control: a review[END_REF]; Coutanceau & Menard 1998; [START_REF] Rao | Three-dimensionality in the wake of a rotating cylinder in a uniform flow[END_REF][START_REF] Bourguet | Flow-induced vibrations of a rotating cylinder[END_REF]. An example of the vorticity field obtained by direct numerical simulation (DNS) for a non-rotating and a rotating cylinder (figure 1) shows that rotation can be used to control vortex shedding by creating an asymmetric flow field and when the rotation rate is high enough (figure 1b), the vortex street observed for a fixed cylinder totally disappears. Flow control around an obstacle in reactive conditions has received less attention [START_REF] Mejia | Threedimensionality in the wake of a rotating cylinder in a uniform flow[END_REF]) and many questions arise regarding flame stabilisation in these conditions.

Therefore, the general aim of the present work is to better understand the stabilisation mechanisms of laminar methane/air premixed flames in the wake of a rotating cylinder. Based on existing non-reacting flow studies, it is expected that rotation will change the flow field, the flame roots location and thus the mean flame topology. Note that this study only deals with steady flames even though the cylinder is rotating. Another goal of this study is that it is a novel technique to create original steady laminar flame structures that differ from the usual laminar premixed flame configurations: planar, cylindrical [START_REF] Kitano | A study of cylindrical premixed flames with heat loss[END_REF][START_REF] Groot | A computational study on propagating spherical and cylindrical premixed flames[END_REF] or spherical [START_REF] Dowdy | The use of expanding spherical flames to determine burning velocities and stretch effects in hydrogen/air mixtures[END_REF][START_REF] Chen | On the extraction of laminar flame speed and Markstein length from outwardly propagating spherical flames[END_REF][START_REF] Bonhomme | Curvature and confinement effects for flame speed measurements in laminar spherical and cylindrical flames[END_REF].

The specific objectives of the paper are to: (i) analyse the stabilisation of a laminar premixed flame behind a rotating flame holder and compare experimental to numerical results; (ii) analyse the different flame structures in order to understand the key mechanisms leading to the observed flames; and (iii) determine the effects of heat transfer between the flow and the cylinder. This investigation is conducted with both experimental and numerical approaches. The paper first describes the experimental and numerical methods used to study flame stabilised on rotating cylinders ( § 2). Next, flame topologies are presented as a function of the rate of rotation of the cylinder and experimental results are compared to DNS data ( § 3). Third, the structure of each flame branch is analysed using DNS data, taking into account non-adiabadicity in the vicinity of the cylinder: effects of stretch, dilution and cooling will be studied ( § 4). Finally, the distribution of heat transfer on the cylinder wall is studied to examine the interaction between the flame fronts and the rotating cylinder ( § 5).

Experimental and numerical methods

Experimental set-up

The experimental set-up consists of a rectangular cross-section combustion chamber (h = 34 mm by l = 94 mm), operated at atmospheric pressure, in which a stainless steel cylinder with a diameter d = 8 mm, serves as the flame holder (figure 2a). The cylinder is connected to a brushless direct current electric motor with rotational speed ranging from 600 to 20 000 revolutions per minute (RPM). During experiments, the motor experienced low frequency rotation fluctuations so that the uncertainty on the rotation speed was ±100 RPM. The flame shape is however weakly sensitive to these fluctuations so that it remains steady. The combustion chamber is equipped with three planar quartz windows allowing optical access. The lean air-methane mixture is injected in a plenum through six injectors. The flow is then laminarised by a bed of 1 mm glass balls and two honeycombs. The upper part of the plenum is water cooled to ensure a constant temperature for the fresh mixture (T in = 288.15 K). The velocity field is measured with particle image velocimetry (PIV). A double cavity Nd:YAG laser (Quantel Big Sky), operating at 532 nm, fires two laser beams, with a delay of 200 μs. The laser beam is expanded through a set of fused silica lenses (spherical and diverging). Since the flow is not symmetric when the cylinder is rotating, a tilted mirror is inserted to lighten the cylinder shadow region with a reflection of the incoming laser sheet (figure 2a). The laser sheet thickness is measured to 500 μm. Olive oil particles of 1 μm are seeded through two injection systems located just before the glass balls bed (venturi seeder). Mie scattering is collected on a Imager Intense camera (Lavision), operating at a repetition rate of 1 Hz and a resolution of 1376 × 1040 pixels. A f/16 182 mm telecentric lens (TC4M64, Opto-engineering) is used to reduce parallax displacements occurring with classical lenses and obtain more accurate velocity vector fields in the vicinity of the cylinder. A narrow band-pass optical filter, centred at λ c = 532 ± 10 nm is inserted in front of the camera sensor to record only Mie scattering of olive oil particles. PIV images are processed with a cross-correlation multi-pass algorithm (Davis 8.2.3), resulting in a final window of 16 × 16 pix 2 and a 50 % overlap. Ninety images are collected for each operating condition, with a spatial vector spacing of 0.2 mm.

The flame is also visualised thanks to a second Imager Intense camera (Lavision), equipped with a similar telecentric lens as the PIV system. The exposure time is set to 500 ms and a band-pass optical filter, centred at λ c = 430 ± 10 nm, is placed in front of the camera to record spontaneous CH * emission.

A 50/50 beam splitter is inserted between the two cameras to simultaneously record PIV and CH * signals, with only one optical access for visualisation. The experimental field of view (EFOV) shown in figure 2(b) captures a part of the combustion chamber section (36 mm × 26 mm) and the origin is taken at the centre of the cylinder.

The inlet temperature of fresh gases T in and the ambient pressure P 0 are systematically checked prior to measurements, and images are recorded when thermal equilibrium of the combustion chamber is reached (measured with a K-thermocouple fixed on the burner structure). During experiments, the flame started to slightly flicker for certain operating conditions (pressure fluctuations less than 2 Pa). These effects contribute to the slight flame thickening observed in the experiments.

Direct numerical simulations

The reactive multi-species Navier-Stokes equations are solved using a fully compressible unstructured solver called AVBP [START_REF] Schonfeld | Steady and unsteady flows simulations using the hybrid flow solver avbp[END_REF][START_REF] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF]. A major advantage of the configuration of figure 2 is that it is almost perfectly two-dimensional so that a two-dimensional DNS with detailed chemistry can be used. Chemical kinetics are described by the 'Lu 19' analytical mechanism of [START_REF] Lu | A criterion based on computational singular perturbation for the identification of quasi steady state species: a reduced mechanism for methane oxidation with no chemistry[END_REF] for methane/air combustion. This mechanism involves 19 transported species and 11 quasi-steady state (QSS) species. The viscosity is modelled with a power law, each species has its own Lewis number and a constant Prandtl number is assumed. It was validated for a large range of equivalence ratio (φ = 0.5-1.5), pressure (P 0 = 1-30 atm) and correctly predicts auto-ignition times as well as CO levels. It was implemented in the DNS code and validated by comparison with GRI-MECH3.0 computations using the Cantera open source solver [START_REF] Smith | [END_REF][START_REF] Goodwin | Cantera C++ users guide[END_REF]. Laminar flame speeds and adiabatic temperatures are correctly reproduced (figure 3b). In addition, properties of flame/wall interactions were validated by running one-dimensional flames. For the canonical head-on quenching scenario, the maximum normalised wall heat flux is 0.35 and the corresponding Peclet number is 3, being in agreement with the literature [START_REF] Huang | Heat transfer during laminar flame quenching: effect of fuels[END_REF][START_REF] Jarosinski | A survey of recent studies on flame extinction[END_REF][START_REF] Wichman | Head-on quenching of a premixed flame by a cold wall[END_REF]. The laminar flame thickness for the φ = 0.7 equivalence ratio flame studied here is 0.680 mm while the mesh resolution in the flame zone is 0.06 mm so that at least 12 points are used to resolve the flame front. Mesh sensitivity was verified by testing a finer mesh, with a cell resolution of 25 μm (instead of 60 μm in this study), leading to negligible modifications in flame position and velocity field.

The two-step Taylor-Galerkin convection scheme [START_REF] Colin | Development of high-order Taylor-Galerkin schemes for unsteady calculations[END_REF] used for DNS provides third-order accuracy in time and space. This scheme is characterised by a high spectral resolution, excellent dispersion and dissipation properties. Navier-Stokes characteristic boundary conditions (NSCBC) [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF][START_REF] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF]) are applied at inlet and outlet. Lateral walls were modelled as no-slip and isothermal (T w = 288 K). In order to account for the rotation of the cylinder, the velocity of its wall is set at u θ = ω/2 in the frame of the burner, where ω is the angular rotation speed. The temperature T c is imposed at the cylinder external face (see § 2.2.1). The computational domain extends from z = -35 to z = 45 mm, which corresponds to the shaded zone in figure 2(b). The velocity profile at the inlet (figure 3a) is imposed by fitting experimental velocity data.

Cylinder temperature

Many recent studies emphasise the importance of the flame holder temperatures on the flame structure and dynamics [START_REF] Duchaine | Sensitivity analysis of transfer functions of laminar flames[END_REF][START_REF] Kedia | An analytical model for the prediction of the dynamic response of premixed flames stabilized on a heat-conducting perforated plate[END_REF], 2015;[START_REF] Mejia | Joint experimental and numerical study of the influence of flame holder temperature on the stabilization of a laminar methane flame on a cylinder[END_REF]. The temperature of the flame holder modifies the quenching distance which is the minimum distance between flame front and wall. This distance then controls the flame dynamics so that an accurate DNS of this flame requires a precise determination of the cylinder temperature. This does not necessarily require a full coupled simulation of heat transfer between the cylinder and the flow. An estimation of the Biot number Bi is 0.01, so that very small temperature gradients are expected in the cylinder. Moreover, a one-dimensional model of the wall temperature, considering a periodic and unsteady wall heat flux (related to the rotation speed) gives penetration depths of around 1 mm with negligible temperature fluctuations in the bluff body itself. This proves that the cylinder responds to an average heat flux and thus its temperature is constant when steady state is reached. The temperature within the cylinder can be assumed to be homogeneous and constant so that a simplified lumped model is used to evaluate its temperature, starting from the unsteady heat equation, integrated over the cylinder:

ρC p V dT c dt = S q • n dS, (2.1)
where T c is the cylinder temperature, S is its surface, V its volume and q the energy flux at the fluid-solid interface. Thus the cylinder temperature T c is:

dT c dt = S ρC p V (Φ cond -Φ rad ), (2.2)
where radiative Φ rad and convective Φ cond fluxes are considered at the cylinder boundary. The radiative flux Φ rad = σ (T 4 c -T 4 ext ) is assumed to be absorbed by the combustion chamber walls at temperature T ext = T w . Additionally, gases are considered as transparent and do not radiate to the cylinder. In the present configuration,

α u b (m s -1 ) φ T in (K) P 0 (bar) T c (K) (ωr/u b )

2.3.

Operating conditions A large range of operating conditions have been tested in order to determine the stability map of the burner. Based on these observations, the bulk velocity taken in the plenum was set to u b = 1.07 m s -1 and the equivalence ratio of the air-methane mixture to φ = 0.70. The Reynolds number Re of the incoming flow is 580 (based on the cylinder diameter d). To characterise the rate of rotation of the cylinder, the normalised parameter α is introduced as:

α = ωd 2u b . (2.3)
This parameter compares the radial velocity at the cylinder to the bulk flow velocity u b , and is the control parameter for non-reactive flows around rotating cylinders. The fresh gas temperature is T in = 288.15 K and the outlet pressure P 0 = 0.99 bar. In this study, all parameters are kept constant except for α whose variations are reported in table 1 together with the corresponding cylinder temperature.

Experimental and numerical results: flame classification

Figure 4 shows the normalised time averaged CH* signal (grey scale) and an iso-contour of the heat-release rate extracted from the DNS (20 % of the maximum, dashed lines) when the flames have reached steady state. Even though CH* traces do not match exactly reaction rate zones, they are a good indicator of the flame position. The flame roots location obtained in the experiments is well reproduced with the DNS. However, some differences on the flame angle are observed for z/d > 2. This latter effect is caused by the chemical scheme that overestimates the laminar burning velocity s 0 l (figure 3b), resulting in a slightly wider flame angle. A V-shape symmetric flame is observed when there is no rotation, i.e. α = 0.00 (figure 4a).

The flame intensity decreases when approaching the cylinder due to the flame/wall interaction. During the transient period, the cylinder temperature increases due to the presence of burned gases in the wake of the cylinder and flame roots move closer to the cylinder. When thermal equilibrium is reached, the flame roots stabilise at z/d = 0.66. Profiles of axial velocity u (normalised by the bulk velocity upstream of the cylinder, u b = 1.07 m s -1 ) without rotation (α = 0, figure 5, first row) show that the flow speed upstream of the cylinder (z/d = -0.9) is affected by the presence of the cylinder. At z/d = 0.0, the axial velocity exceeds the bulk velocity u b by a factor of ≈1.7 because of the area restriction (figure 5b, top). Negative u values at z/d = 0.9 (figure 5c, top) reveal a weak symmetric recirculation zone in the wake of the cylinder.

-2 When the cylinder has a moderate rotation rate, the flame becomes asymmetric (figure 4b: α = 1.16 and c: α = 2.30). Indeed, the lower branch stabilises close to the cylinder and the flame root is more curved when the rotation rate is increased, thus leading to a L-shape flame (figure 4c). This branch exhibits a maximum intensity near the cylinder, and decreases along the flame. In contrast, the upper branch is not anchored to the cylinder any more and it is located far downstream of the cylinder. Its intensity is very low and this flame branch is virtually quenched over a large region downstream of the flame holder. Velocity profiles (figure 5, middle row) show an asymmetric flow at the cylinder and downstream. The cylinder rotation locally increases the flow velocity by a factor of 2, and the recirculation region is affected, being shifted in the direction of the cylinder rotation.

-1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -4 0 2 4 1 -1 0 2 1 -1 0 2 (a) ( b) (c) Cylinder Cylinder Cylinder
When α is increased to higher values (figure 4d: α = 3.07 and 4e: α = 4.10), the lower branch is distorted in the direction of rotation by the high velocities induced by the cylinder. A clear L-shape is observed and the lower flame root location is still close to the cylinder while its intensity remains high. The upper branch is now located closer to the cylinder with higher reactivity, almost leading to an asymmetric U-shape flame.

The velocity profile upstream of the cylinder is affected (figure 5a bottom) as it shows a slight asymmetric shape. Finally for larger rotation speeds, the flame topology changes suddenly and it stabilises upstream of the cylinder, a situation which is not discussed in this study. The effect of rotation is visible on the flow streamlines, as shown in figure 6 which displays the flow and temperature fields for α = 1.16 and 4.10. For the moderate rotation rate (figure 6a), the streamlines are deflected in the direction of the cylinder rotation. When the rotation rate is increased (figure 6b), flow streamlines wrap around the cylinder and also shift the temperature iso-contours in the direction of rotation. Moreover, the induced cylinder shear flow may be sufficiently high to continuously carry burned gases to the leading edge of the cylinder.

The five flames of figure 4 exhibit non-trivial stabilisation mechanisms. At moderate rotation rates, the upper branch of the flame is quenched over a long distance so that it is anchored far downstream of the cylinder. In contrast, the lower branch exhibits higher reaction rates in the vicinity of the cylinder (figure 4b,c).

Figure 7 depicts the different mechanisms involved in this configuration. When the cylinder is rotating, burned gases are carried towards the upper branch while they transfer their energy to heat up the cylinder (Φ bg-c ). In parallel, the warm cylinder also heats up the incoming fresh gases (Φ c-ug ). As a consequence, the flame roots location and the flame structures are conditioned by heat transfer between the cylinder and the surrounding flow. Moreover, the flow induced by the cylinder rotation mixes fresh and burned gases, and that is dependent on the level of rotation. Figure 8 shows how this mixing affects chemistry. It compares scatter plots of CH 4 , CO 2 , CH 3 and CO mass fractions versus the progress variable c in two regions of the upper branch: zone A is far away from the cylinder (black dots) whereas region C is located in the quenched zone of the upper flame (grey dots). Flame structures in region A are all similar and collapse on a single curve, showing that, far from the cylinder, there is a unique flame structure (which could be tabulated for example). However, the scatter plots in region C reveal multiple flame structures. Burned gases recirculation and subsequent dilution lead to higher CO 2 and lower CH 4 mass fractions. In turn, formation of radicals and intermediates (CH 3 and CO) is inhibited. In this zone, lowtemperature chemistry plays a major role, which is captured by the 'Lu 19' scheme and would not be correctly predicted using tabulation of a global mechanism. The aim of the next sections is to analyse the two branches of figure 7 in terms of flame structure and flame/wall interaction when the wall (here the cylinder) moves.

Analysis of the flame structures

Excess or defect of flow enthalpy

As shown in figure 7, the flow is not adiabatic in the vicinity of the cylinder because of heat transfer with the flame holder. The proper quantity to measure the departure from adiabadicity is the reduced enthalpy defect/excess L given by:

L = h t -h in t QY in F , (4.1)
where h t is the total chemical enthalpy taking into account sensible, chemical and kinetic energies. h in t is a reference enthalpy, taken in the fresh gases. The normalizing term QY in F in (4.1) is the heat of the reaction. Positive values of L correspond to sensible or chemical energy given to the fresh gases: the flame is therefore superadiabatic. Conversely, negative values of L correspond to an enthalpy loss and a subadiabatic flame. Figure 9 presents reduced enthalpy fields for different rotation rates α. Far away from the cylinder and the flame fronts, the flow is adiabatic and L = 0. When the cylinder is not rotating (α = 0, figure 9a), L is negative in the wake of the cylinder because hot gases loose energy to the cylinder through conductive heat transfer. Conversely, L is positive upstream of the cylinder because incoming gases are heated by the cylinder. The latter trend is still valid for moderate values of α and even prominent as L goes up to 0.2. However, this region of reduced enthalpy gain narrows and eventually vanishes for larger rotation rates (figure 9e). The distribution of L around the cylinder is shown in figure 10 (θ is defined in figure 7). When the cylinder is at rest or has a low rotation rate (i.e. α = 0.00 or 1.16), L is negative in the wake of the cylinder where burned gases are located, and it is positive on the upstream part where fresh gases impact the hot cylinder. Increasing the rotation rate reduces L all along the cylinder so that it is almost negative at all angles for the higher rotation rate (α = 4.10). Therefore, it is speculated that the flow induced by the rotation (figure 6) may explain this trend.

Dilution effect at high rotation rates

Figure 11 displays the field of CO 2 mass fraction at α = 1.16 (a) and α = 4.10 (b). For α = 4.10, CO 2 is found all around the cylinder. This recirculation of burned gases, also observed in figure 6(b), changes the local composition by diluting incoming fresh gases and diminish L to negative values (figure 10).

To separate preheating and dilution effects, L can be written as:

L = 1 QY in F T T 0 C p dT L s + 1 QY in F k h 0 f ,k (Y k -Y in ) L c + 1 2QY in F [(u 2 i -u in 2 i )] L k , (4.2)
where L s measures gas preheating by the cylinder while L c measures the effect of dilution and L k accounts for the kinetic effect. Because the Mach number is 0.003, L k can be neglected. Figure 12 plots the reduced sensible L s (black line) and reduced chemical L c enthalpies (dashed line) contributing to L, for two rotation rates α = 1.16 (figure 12a) and 4.10 (figure 12b). The shape or the amplitude of L s is similar in both cases because gases near the wall receive energy from the hot cylinder. However, the contribution of L c differs from the two cases. At a moderate rotation rate (figure 12a), the reduced chemical enthalpy parameter L c is null on the upstream part of the cylinder, meaning that only fresh unburned gases are found at this location.

In contrast, when the rotation rate is increased (figure 12b), L c is negative on the leading edge of the cylinder. Therefore, the dilution observed in figure 11(b) is the parameter that decreases L at high rotation rates. In summary, the excess or defect of the reduced enthalpy L around the cylinder evolves due to the competition between two effects. First, the warmer cylinder heats up near-wall gases that results in positive values of L s for all rotation rates investigated. Second, the induced flow created by the cylinder carries burned gases that diminish the reduced chemical enthalpy L c . This latter effect is prominent when the rotation rate exceeds a certain threshold (α c ≈ 3). 4.3. Super-and subadiabatic flame structure This section focuses on the differences between the two flame branches observed in figure 4 and their link with the reduced enthalpy field (figure 9). Figure 13 shows two flame profiles extracted from the DNS (α = 1.16 and z/d = 1.25). The comparison between the overdriven lower branch (dashed lines with symbols) and an adiabatic unstrained flame (lines) in figure 13(a,b) shows a good agreement for the species profiles but the maximum heat-release rate is much higher for the flame of the DNS. This observation is also verified with the reduced enthalpy parameter as this branch is located in a superadiabatic flow (L > 0, figure 13b). In contrast, the upper branch (figure 13c,d) presents smoother species profiles and a heat-release rate profile quite low. Similarly, as the upper branch stabilises in a subadiabatic flow, no self-sustained or weak flame can exist at this location (L < 0, figure 13d). Figure 14 presents the evolutions of the maximum heat release rate along the lower flame branch for different rotation rates. The location of maximum heat-release rate is detected with a crest sensor, for both flames and these maximum values are plotted along the curvilinear abscissa s normalised by the flame thickness of a planar unstretched adiabatic flame δ 0 L = (T b -T u )/ max(dT/dx) (0.680 mm for the φ = 0.7 equivalence ratio studied here). When the cylinder is not rotating (figure 14, solid line), the heat-release rate monotonically increases to reach the value obtained in an adiabatic flame. Close to the cylinder (s/δ 0 L < 10) the reduced heat-release rate drops, due to the heat losses to the cylinder and the recirculation zones. The zone over which the cylinder inhibits chemical reactions is about 10δ 0 l , being larger in comparison to canonical head-on quenching scenarios [START_REF] Baum | Progress in knowledge of flamelet structure and extinction[END_REF][START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]. When the cylinder is rotating at a moderate level (α = 1.16, 2.30), the quenching zone narrows to about δ 0 l , and the heat-release rate increases near the wall, being locally twice as intense as in an adiabatic flame. When the rate of rotation is much higher (α = 3.07, 4.10), the flame intensity starts to decrease, and it is even lower than in the non-rotating case when α = 4.10. This behaviour is an indicator of the dilution effect mentioned previously.

Figure 15 presents the reduced heat-release rate along the upper flame. In contrast with the lower superadiabatic flame branch, all cases exhibit heat-release rates lower than adiabatic conditions. The zone over which the flame is quenched becomes very large. For the case α = 2.30, the flame has not recovered its adiabatic structure after 40δ 0 l . However, beyond a critical rotation rate (α c ≈ 3), as dilution of burned gases shield the cylinder with burned gases, the flame starts to retrieve its adiabatic behaviour.

At this point, L may be an adequate parameter explaining the flame topologies near the cylinder. In other words, could the two flame structures be explained by considering only local values of L? It is therefore relevant to compute non-adiabatic stretched planar flame with a one-dimensional solver and to compare the local flame structures obtained with DNS flame profiles.

Two different counterflow flame prototypes are considered to mimic the two flame branches only with the reduced enthalpy parameter L. The lower branch corresponds to a superadiabatic flame where fresh gases have been preheated (figure 16a) whereas the upper branch corresponds to a flame where fresh gases interact with burned gases that have lost energy to the cylinder (figure 16b). However, the computation of a one-dimensional counterflow flame also introduces a strain effect. The range of strain rates to be used for the Cantera computations is evaluated by post-processing the DNS. Figure 17 shows the probability density functions of strain rate κ for the two branches for α = 0 to 4.1. For all rotation rates, strain rate values are quite low and the possible mechanism of flame extinction due to an important strain rate is not supported here [START_REF] Vagelopoulos | Laminar flame speeds and extinction strain rates of mixtures of carbon monoxide with hydrogen, methane, and air[END_REF][START_REF] Williams | Progress in knowledge of flamelet structure and extinction[END_REF][START_REF] Shanbhogue | Flame macrostructures, combustion instability and extinction strain scaling in swirl-stabilized premixed CH 4 /H 2 combustion[END_REF]. A small value κ = 300 s -1 is therefore used in the Cantera calculations integrating the heat release rate on profiles normal to s:

S c = 1 ρY 0 F +∞ -∞ ωF dn, (4.3)
where ωF , ρ and Y 0 F are the fuel source term, density and fuel mass fraction in the fresh gases, respectively. Data from the DNS have been obtained from profiles taken at s/δ 0 L = 0.7, 4.7 and 40 for the overdriven branch, and s/δ 0 L = 0.1, 5.0 and 40 for the cooled branch, and for different rotation rates. The corresponding Cantera results for the two counterflow flame prototypes are shown with the black line. The agreement between the one-dimensional and DNS flame properties is very good, showing that the flame structures observed locally in the DNS can be described to first order by the reduced enthalpy loss parameter L. However, some discrepancies are observed for the lower, superadiabatic branch where curvature effects could modify the flame structure (label 3). This comparison with one-dimensional flames proves that L is a robust and adequate parameter to describe the behaviour of the observed flames. Another approach to study this peculiar flame stabilisation is to look at the details of chemical kinetics, as for example done by [START_REF] Michaels | Impact of the bluff-body material on the flame leading edge structure and flame-flow interaction of premixed CH 4 /air flames[END_REF]. This analysis for the rotating cylinder is left to a further study.

Flame/moving wall interaction

Flame-wall interaction (FWI) is a classical topic in combustion (Ezekoye, Greif & Lee 1992;[START_REF] Poinsot | Direct simulation and modeling flame-wall interaction for premixed turbulent combustion[END_REF][START_REF] Wichman | Head-on quenching of a premixed flame by a cold wall[END_REF][START_REF] Bruneaux | Premixed flame wall interaction in a turbulent channel flow: budget for the flame surface density evolution equation and modelling[END_REF]. In all these studies, the wall is fixed and the main objective is to determine the maximum heat flux it must sustain during FW interaction. Reproducing numerically flame/wall interactions is dependent on the chemical scheme, that needs to mimic the exothermic zero-activation-energy radical recombination, and the chemical reactivity of the wall [START_REF] Popp | Heterogeneous/homogeneous reactions and trasport coupling during flame-wall interaction[END_REF][START_REF] Popp | An analysis of wall heat fluxes, reaction mechanisms and unburnt hydrocarbons during the head-on quenching of a laminar flame methane flame[END_REF][START_REF] Gruber | Turbulent flame-wall interaction: a direct numercial simulation study[END_REF][START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]. The present experiment allows to investigate an interaction between a flame and a moving wall (FMW), which is not the case for classical FWI. For instance, in side wall quenching (SWQ), a flame propagates along a fixed wall (figure 20b) and the maximum flux to the wall (scaled by the flame power ρs 0 l QY F ) is of the order of 0.15. In the case of head-on quenching (HOQ) where a flame impinges normal to the wall, It can go up to 0.35 (figure 20a) on fixed walls.

Figure 21 presents the reduced wall heat flux F along the cylinder (θ is the angular position), for different rotation rates. All fluxes are scaled by the flame power ρs 0 l QY F , and the azimuthal locations of the lower and upper flame roots are marked with a circle and a square, respectively. When the cylinder is at rest (α = 0), the wall heat flux is symmetric. As the cylinder heats up gases, negative reduced fluxes are found on the upstream part of the cylinder. In contrast, the positive bump in the centre denotes the flux from burned gases to the downstream face of the cylinder. The maximum reduced wall heat flux is found to be approximately 0.20, which agrees with a typical SWQ flame wall/interaction [START_REF] Lu | Unsteady heat transfer during side wall quenching of a laminar flame[END_REF][START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]. When the cylinder is rotating with moderate rotation rates (α = 1.16 or 2.30), the profile of the wall heat flux becomes asymmetric. A maximum is visible in the vicinity of the lower flame root (circles). The maximum flux obtained at this location reaches a value of 1.1 for α = 2.30, which exceeds values for SWQ (0.15) or HOQ (0.35) cases [START_REF] Wichman | Head-on quenching of a premixed flame by a cold wall[END_REF][START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]. This shows that the flow induced by the rotating cylinder leads to a more intense flame/wall interaction than found in usual flames. It also indicates that the maximum flux which can be reached in a combustion chamber with moving parts [START_REF] Labarrere | Experimental and numerical study of cyclic variations in a constant volume combustion chamber[END_REF]) can exceed SWQ or HOQ values by almost a factor of three. After the maximum reached at the lower branch location, the local wall heat flux decays in the direction of the cylinder rotation, and the upper flame root does not create a new local flame maximum as this flame is subadiabatic and stabilises far downstream of the cylinder. For higher rotation rates (α > 3), the maximum heat flux decreases, even if it remains higher than for classical FW interaction. This trend can be explained by the dilution effect seen in figures 10 and 12. The flow induced by the cylinder creates a layer of burned gases that insulate the cylinder from the flame roots. This leads to a thermal shielding and explains that the cylinder temperature T c does not increase any more (table 1).

Conclusion

The stabilisation of a lean premixed laminar methane/air flame on a rotating cylindrical bluff body was investigated with experiments and DNS to obtain more insights on flame stabilisation. All flames studied were steady both in the experiment and in the DNS. Computations were carried out on a two-dimensional grid with an analytically chemical scheme that could accurately capture flame/wall interactions scenarios and low-temperature chemistry in the wake of the cylinder. In addition, a lumped model was implemented to adapt the cylinder temperature with an energy balance coupled to the DNS. An excellent agreement is found between DNS and experiments over a wide range of rotation rates (ratio of the cylinder's surface speed to the incoming cold flow) and results show that the bluff-body rotation breaks the symmetry of the flow and creates two distinct flame branches. Surprisingly, a bifurcation takes place at high rotation rates. Beyond a certain rotation rate, the flow induced by the rotation of the cylinder circulates burned gases all around the cylinder, thus creating a layer of gases that changes locally the gas composition. As a result, it reduces the wall heat flux and the cylinder temperature does not increase any more.

DNS results show that the reduced enthalpy defect/gain L is an appropriate control parameter to describe the structure of the two branches. Strain remains small and does not contribute to the flame stabilisation. The superadiabatic branch (positive L) burns fresh gases that have been preheated while flowing along the hot cylinder. In contrast, the subadiabatic branch (negative L) interacts with combustion products that have lost their enthalpy while being convected by the flow induced by the cylinder rotation. The local structures of the two branches can be predicted in terms of maximum heat-release rate and consumption speed using planar unstretched non-adiabatic flames (computed with Cantera) and considering only the local value of the reduced enthalpy parameter L.

This study also unveiled a configuration to investigate FMW (flame/moving wall), something that was not studied before. The wall heat fluxes generated by the interaction of the subadiabatic flame with the moving walls of the cylinder are lower than what is observed when a premixed flame hits a wall in head-on (HOQ) or side-wall quenching (SWQ) configurations. However, the superadiabatic flame creates a wall heat flux to the cylinder that exceeds HOQ by a factor of almost three, showing that the interaction between a flame and a moving wall can be very different and more intense than it is in classical HOQ or SWQ situations.

  FIGURE 1. (Colour online) Iso-contours of vorticity, coloured by vorticity in the case of a non-reacting flow around (a) a non-rotating cylinder and (b) a rotating cylinder (α = 4.10, α is the rotation rate, see (2.3)).

FIGURE 2 .

 2 FIGURE 2. (a) Schematic of the rotating cylinder in the present configuration. (b) Experimental field of view (EFOV) and definition of the controlling parameters. The shaded area corresponds to the computed domain.

  FIGURE 3. (a) Experimental inlet velocity profile (hot-wire anemometry) measured at z/d = -4.6 and fit used in the simulation (figure 2a). (b) Laminar burning velocities for air-methane mixtures. Comparison between different solvers, chemical schemes and the experiment of Varea et al. (2012).

  FIGURE 4. (Colour online) Comparison between experimental and DNS flame fronts for increasing rotation rates α. CH* traces from experiments are shown in grey scale and red dashed lines are an iso-contour of the heat release rate from the DNS (20 % of the maximum). (a) α = 0.00; (b) α = 1.16; (c) α = 2.30; (d) α = 3.07; (e) α = 4.10. The temperatures of the cylinder are given in table 1.

FIGURE 5 .

 5 FIGURE 5. Experimental (circles) and DNS (line) mean axial velocity profiles (u/u b ) at three axial locations (columns) and for three different α (rows). The dashed line corresponds to u/u b = 0.0. (a) z/d = -0.9; (b) z/d = 0.0; (c) z/d = 0.9.

  FIGURE 6. (Colour online) Temperature iso-contours (black lines) and flow streamlines (grey lines) for cases: (a) α = 1.16 and (b) α = 4.10. The flame location is marked by the heat release crest (red dashed line).

FIGURE 8 .

 8 FIGURE 7. Schematic representation of flame stabilisation mechanisms on a rotating cylinder.

FIGURE 10 .

 10 FIGURE 9. (Colour online) Field of enthalpy loss L for different rotation rates α. Isocontour lines of L are shown in black. The red thick lines are the heat-release rate crest of each branch. (a) α = 0.00; (b) α = 1.16; (c) α = 2.30; (d) α = 3.07; (e) α = 4.10.

FIGURE 12 .

 12 FIGURE 11. (Colour online) Reduced carbon dioxide mass fraction c co 2 field for (a) α = 1.16 and (b) α = 4.10. Black lines are iso-contour of c co 2 .

  FIGURE 13. (Colour online) Flame profiles extracted from the DNS at z/d = 1.25 for the case α = 1.16. n is an axis normal to the flame front. Solid lines correspond to an adiabatic unstrained flame. Dashed lines with markers represent the overdriven and cooled branches. (a) Species profiles on the overdriven branch. (b) Reduced heat-release rate profiles on the overdriven branch. (c) Species profiles on the cooled branch. (d) Reduced heat-release rate profiles on the cooled branch.

FIGURE 14 .

 14 FIGURE 14. Normalised maximum heat-release rate along the lower flame front (overdriven) branch, for rotation rates α from 0.00 to 4.10.

FIGURE 15 .FIGURE 16 .FIGURE 17 .FIGURE 18 .

 15161718 FIGURE 15. Normalised maximum heat-release rate along the curvilinear abscissa of the upper flame (subadiabatic) branch, for various rotation rates α.

FIGURE 19 .

 19 FIGURE 19. Evolution of the reduced consumption speed S c /S adia c versus enthalpy loss L, for various rates of rotation. Planar unstretched flames computed with Cantera are plotted for comparison (solid black line).

FIGURE 20 .FFIGURE 21 .

 2021 FIGURE 20. Two different types of flame/wall interaction. (a) Head-on quenching; (b) side wall quenching.

TABLE 1 .

 1 Operating conditions investigated in the present study. The last column gives the temperature of the cylinder (T c ) obtained with (2.2).S/(ρC p V) = 0.003 m 2 KJ -1 and the emissivity of the cylinder is set to = 0.9, corresponding to an unpolished used surface. Equation (2.2) is solved at run time and provides the value of the equilibrium wall temperature T c used as boundary condition for the DNS. The value of S/(ρC p V) is only relevant for the resolution of the transient evolution of T c . Since we are only interested in the steady state, in practice this value is increased in order to speed up the transient phase. Depending on the operating conditions, the radiative heat flux accounts from 12 % to 25 % of the total

		Bulk	Equivalence	Fresh gas	Pressure
		velocity	ratio	temperature
	0.00	1.07	0.70	288	0.99	508
	1.16	-	-	-	-	707
	2.30	-	-	-	-	783
	3.07	-	-	-	-	787
	4.10	-	-	-	-	773

heat flux. Neglecting its contribution would result in incorrect cylinder temperature and wrong flame stabilisation pattern

[START_REF] Mejia | Joint experimental and numerical study of the influence of flame holder temperature on the stabilization of a laminar methane flame on a cylinder[END_REF]

.

† Email address for correspondence: pradip.xavier@imft.fr

Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement ERC-AdG 319067-INTECOCIS. This work was granted access to the high-performance computing resources of IDRIS under the allocation x20162b7036 made by Grand Equipement National de Calcul Intensif. Access to the PIV/CH* acquisition system was obtained through the Toulouse Fermat Federation. The CFD team at Cerfacs is acknowledged for its support on the AVBP code.

The authors wish also to thank M. Marchal and S. Cazin from IMFT for their precious help with the experimental diagnostics, and G. Daviller (Cerfacs), C. Lapeyre (IMFT), C. Kraus (IMFT) and B. Bedat (IMFT) for useful discussions.


Yohann Uguen
email: yohann.uguen@insa-lyon.fr

Florent De

Steven Derrien
email: steven.derrien@univ-rennes1.fr

High-Level Synthesis Using Application-Specific Arithmetic: A Case Study

On the one hand, a strength of FPGAs is their ability to perform non-standard computations not supported by classical microprocessors. Many libraries of highly customizable application-specific IPs have been developed to exploit this strength.

On the other hand, HLS tools, which allow to program an FPGA using a dialect of the C language, are gaining traction. However, the ease of use of the C language becomes a hindrance when one wants to express non-standard computations. Indeed, the C language was designed for programming microprocessors and carries with it many restrictions of the microprocessor paradigm. This is especially true when computing with floatingpoint, whose data-types and evaluation semantics are defined by the IEEE-754 and C11 standards. If the high-level specification was a computation on the reals, HLS imposes a very restricted implementation space.

This work attempts to bridge FPGA application-specific efficiency and HLS ease of use. This is illustrated on the ubiquitous floating-point summation-reduction pattern. A source-to-source compiler rewrites, inside critical loop nests of the input C code, selected floating-point additions into sequences of simpler operators using non-standard arithmetic formats.

Evaluation of this method demonstrates that the benefits of application-specific operators (better performance and better accuracy) can be brought to HLS workflows while keeping their ease of use.

1 Catapult C Synthesis, Mentor Graphics, 2011, http://calypto.com/en/ products/catapult/overview/ 2 Arbitrary-size floating-point should follow some day, it is well supported by mature libraries and tools

I. INTRODUCTION

Many case studies have demonstrated the potential of Field-Programmable Gate Arrays (FPGAs) as accelerators for a wide range of applications, from scientific or financial computing to signal processing and cryptography. FPGAs offer massive parallelism and programmability at the bit level. This enables programmers to exploit a range of techniques that avoid many bottlenecks of classical von Neumann computing: dataflow operation without the need of instruction decoding; massive register and memory bandwidth, without contention on a register file and single memory bus; operators and storage elements tailored to the application in nature, number and size.

However, to unleash this potential, development costs for FPGAs are orders of magnitude higher than classical programming. High performance and high design costs are the two faces of the same coin.

Hardware design flow and High-level synthesis: To address this, languages such as C or Java are increasingly being considered as hardware description languages. This has many advantages. The language itself is more widely known than any HDL. The sequential execution model makes designing and debugging much easier. One can even use software execution on a processor for simulation. All this drastically reduces development time.

The process of compiling a software program into hardware is called High-Level Synthesis (HLS), with tools such as Vivado HLS [START_REF]Vivado Design Suite User Guide: High-Level Synthesis[END_REF] or Catapult C 1 among others [START_REF] Nane | A survey and evaluation of fpga high-level synthesis tools[END_REF]. These tools are in charge of turning a C description into a circuit. This task requires to extract parallelism from sequential programs constructs (e.g. loops) and expose this parallelism in the target design. Today's HLS tools are reasonably efficient at this task, and can automatically synthesize highly efficient pipelined dataflow architectures.

They however miss one important feature: they are not able to tailor operators to the application in size, and even less in nature. This comes from the C language itself: its high-level datatypes and operators are limited to a small number (more or less matching the hardware operators present in mainstream processors). Indeed, such a high-level language have been designed to be compiled and run on hardware and not to describe hardware. Therefore, HLS tools performs better for general purpose computing whereas FPGAs performs better for application-specific computing.

The broader objective of this work is to show that HLS tools can produce application-specific hardware and therefore, unleash FPGAs potential again.

Arithmetic in HLS: To better exploit the freedom offered by hardware and FPGAs, HLS vendors have enriched the C language with integer and fixed-point types of arbitrary size 2 . However the operations on these types remain limited to the basic arithmetic and logic ones. Exotic or complex operators (for instance for finite-field or floating-point arithmetic) may The case study in this work is a program transformation that applies to floating-point additions on a loop's critical path. It decomposes them into elementary steps, resizes the corresponding sub-components to guarantee some userspecified accuracy, and merges and reorders these components to improve performance. The result of this complex sequence of optimizations could not be obtained from an operator generator, since it involves global loop information.

For this purpose, we envision a compilation flow involving one or several source-to-source transformations, as illustrated by Figure 1. Before detailing it, we must digress a little on the subtleties of the management of floating-point arithmetic by compilers.

HLS faithful to the floats: Most recent compilers, including the HLS ones [START_REF] Hrica | Floating-point design with vivado HLS[END_REF], attempt to follow established standards, in particular C11 and, for floating-point arithmetic, IEEE-754. This brings the huge advantage of almost bit-exact reproducibility -the hardware will compute exactly the same results as the software. However, it also greatly reduces the freedom of optimization by the compiler. For instance, as floating point addition is not associative, C11 mandates that code written a+b+c+d should be executed as ((a+b)+c)+d, although (a+b)+(c+d) would have shorter latency. This also prevents the parallelization of loops implementing reductions. A reduction is an associative computation which reduces a set of input values into a reduction location. Listing 1 provides the simplest example of reduction, where acc is the reduction location.

The first column of Table I shows how Vivado HLS synthesizes Listing 1 on Kintex7. The floating-point addition takes 7 cycles, and the adder is only active one cycle out of 7 due to the loop-carried dependency. Listing 2 shows a different version of Listing 1 that we coded such that Vivado HLS expresses more parallelism. Vivado HLS will not transform Listing 1 into Listing 2, because they are not semantically equivalent (the floating-point additions are reordered as if they were associative). However as Listing 2 expresses more parallelism, Vivado HLS is able to exploit (second column of Table I). The main adder is now active at each cycle on a different sub-sum. Note that a parallel execution with the sequential semantics is also possible, but very expensive [START_REF] Kapre | Optimistic parallelization of floating-point accumulation[END_REF].

Note that Listing 2 is only here as an example and might need more logic if N was not a multiple of 10. Towards HLS faithful to the reals: Another point of view, chosen in this work, is to assume that the floatingpoint C program is intended to describe a computation on real numbers when the user specifies it. In other words, the floats are interpreted as real numbers in the initial C, thus recovering the freedom of associativity (among other). Indeed, most programmers will perform the kind of nonbit-exact optimizations illustrated by Listing 2 (sometimes assisted by source-to-source compilers or "unsafe" compiler optimizations). In a hardware context, we may also assume they wish they could tailor the precision (hence the cost) to the accuracy requirements of the application -a classical concern in HLS [START_REF] Gort | Range and bitmask analysis for hardware optimization in high-level synthesis[END_REF], [START_REF] Caffarena | High-level synthesis of multiple word-length DSP algorithms using heterogeneousresource FPGAs[END_REF]. In this case, a pragma should specify the accuracy of the computation with respect to the exact result. A high-level compiler is then in charge of determining the best way to ensure the prescribed accuracy.

The proposed approach uses number formats that are larger or smaller than the standard ones. These, and the corresponding operators, are presented in Section II. The contribution of this paper, which are compiler transformations to generate C description of these operators in a HLS workflow, is presented in Section III. Section IV evaluates our approach on the FPMark benchmark suite.

II. THE ARITHMETIC SIDE: AN APPLICATION-SPECIFIC ACCUMULATOR IN VIVADO HLS

The accumulator that we used for this paper is based on a more general idea developed by Kulisch. He advocated a very large floating-point accumulator [START_REF] Kulisch | The exact dot product as basic tool for long interval arithmetic[END_REF] whose 4288 bits would cover the entire range of double precision floating-point. Such an accumulator would remove rounding errors from all the possible floating-point additions and sums of products, with the added bonus that addition would become associative. So far, Kulisch's full accumulator has proven too costly to appear in mainstream processors. However, in the context of application acceleration with FPGAs, it can be tailored to the accuracy requirements of applications. Its cost then becomes comparable to classical floating point operators, although it vastly improves accuracy [START_REF] De Dinechin | An FPGA-specific approach to floating-point accumulation and sum-of-products[END_REF]. This operator can be found in the FloPoCo [START_REF] De Dinechin | High-Performance Computing Using FPGAs, chapter Reconfigurable Arithmetic for High-Performance Computing[END_REF] generator in Altera DSP Builder Advanced. Its core idea, illustrated on Figure 2, is to use a large fixed-point register into which the mantissas of incoming floating-point summands are shifted (top) then accumulated (middle). A third component (bottom) converts the content of the accumulator back to the floating-point format. The sub-blocks visible on this figure (shifter, adder, and leading zero counter) are essentially the building blocks of a classical floating-point adder.

The accumulator described is the one offered in FloPoCo [START_REF] De Dinechin | An FPGA-specific approach to floating-point accumulation and sum-of-products[END_REF], although it is not the contribution of this paper, we included two improvements to it:

• In FloPoCo, Float-to-Fix and Accumulator form a single component, which restricts its application to simple accumulations similar to Listing 1. The two components of Figure 2 enable a generalization to arbitrary summations within a loop, as Section III will show. • Our implementation supports subnormal numbers (special floating-point numbers with leading zeros to their significand to fill the underflow gap around 0). Note that we could have implemented any other nonstandard operator performing a reduction such as [START_REF] Luo | Accelerating pipelined integer and floatingpoint accumulations in configurable hardware with delayed addition techniques[END_REF], [START_REF] Kadric | Accurate parallel floating-point accumulation[END_REF].

A. The parameters of a large accumulator

The main feature of this approach is that the internal fixed-point representation is configurable in order to control accuracy. It has two parameters:

• MSB A is the weight of the most significant bit of the accumulator. For example, if MSB A = 20, the accumulator can accommodate values up to a magnitude of 2 20 ≈ 10 6 . • LSB A is the weight of the least significant bit of the accumulator. For example, if LSB A = -50, the accumulator can hold data accurate to 2 -50 ≈ 10 -15 . Such a fixed-point format is illustrated in Figure 3. The accumulator width w a is then computed as MSB A -LSB A +1, 71 bits in the previous example. 71 bits represents a wide range and high accuracy, and still additions on this format

- Shifter Negate Exponent Mantissa Sign MaxMSB X Shift value w f w e MaxMSB X -LSB A + 1 FloatToFix + w A Accumulator Registers Negate LZC + Shifter Exponent Mantissa Sign FixToFloat w A w e w f
-8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 2 MSB A 2 MSB A -1 2 0 2 LSB A Figure 3:
The bits of a fixed-point format, here (MSB A , LSB A ) = (7, -8).

will have one-cycle latency for practical frequencies on recent FPGAs. If this is not enough the frequency can be improved thanks to partial carry save [START_REF] De Dinechin | An FPGA-specific approach to floating-point accumulation and sum-of-products[END_REF] but this was not useful in the present work. For comparison, for the same frequency, a floating-point adder has a latency of 7 to 10 cycles, depending on the target.

In the following of this paper, we refer as latency of a circuit the number of cycles needed for the entire application to complete.

B. Implementation within a HLS tool

This accumulator has been implemented in C, using arbitrary-precision fixed point types (ap_int). The leading zero count, bit range selections and other operations are implemented using Vivado HLS built-in functions. For modularity purposes, the FloatToFix and FixToFloat are wrapped into C functions (28 LoC for the FloatToFix, 22 LoC for FixToFloat) whose calls inlined to enable HLS optimizations.

Because the internal accumulation is performed on a fixedpoint integer representation, the combinational delay between two accumulations is lower compared to a direct floating point addition. We expect HLS tools to benefit from this delay reduction by taking advantage of more agressive loop pipelining (with shorter Initiation Interval pipeline) resulting in a design with a shorter overall latency.

C. Validation

To evaluate and refine this implementation, we used Listing 3, which we compared to Listings 1 and 2. In the latter, the loop was unrolled by a factor 7, as it is the latency of a floating-point adder on our target FPGA (Kintex-7).

For test data, we use as in Muller et al. [START_REF] Muller | Handbook of Floating-Point Arithmetic[END_REF] the input values c[i]=(float)cos (i), where i is the input array's index. Therefore the accumulation computes

i c[i].
The parameters chosen for the accumulator are:

• MSBA = 17. Indeed, as we are adding cos(i) 100K times, an upper bound is 100K, which can be encoded in 17 bits. • MAXMSBx = 1 as the maximum input value is 1.

• LSBA = -50: the accumulator itself will be accurate to the 50th fractional bit. Note that a float input will see its mantissa rounded by FloatToFix only if its exponent is smaller than 2 -25 , which is very rare. In other words, this accumulator is much more accurate than the data that is thrown to it. The results are reported in Table I for simple and double precision. The Accuracy line of the table reports the number of correct bits of each implementation, after the result has been rounded to a float. All the data in this table was obtained by generating VHDL from C synthesis using Vivado HLS followed by place and route from Vivado v2015.4, build 1412921. This table also reports synthesis results for the corresponding FloPoCo-generated VHDL, which doesn't include the array management.

Vivado HLS uses DSPs to implement the shifts in its floating-point adders. Even if the shifts were implemented in LUTs, the first column would remain well below 500 LUTs: it has the best resource usage. However the latency of one iteration is 7 cycles, hence 100K iterations takes 700K cycles. When unrolling the loop, Vivado HLS is using almost 4 times more LUTs for floats, and 3 times more for doubles. The unrolled versions improves latency over naive versions. Nevertheless, our approach gets even betters latencies for a reasonable LUT usage. Also, we achieve maximum accuracy for the float format which caps at 24 bits (the internal representations of the double, unrolled double and our approach have a higher accuracy than 24 bits, but are then casted to the 24 bits of the float format). Finally, our results are very close to FloPoCo ones, both in terms of LUTs usage, DPSs and latency. Using this implementation method, we also created an exact floating-point multiplier with the final rounding removed as in [START_REF] De Dinechin | An FPGA-specific approach to floating-point accumulation and sum-of-products[END_REF]. This function is called ExactProductFloatToFix. Due to lack of space we do not present it in detail. As the output of this multiplier is not standard, we also created an adapted Float-to-fix block called ExactProductFloatToFix. These functions represent 44 lines of code for ExactProduct and 21 lines of code for ExactProductFloatToFix.

III. THE COMPILER SIDE: GECOS SOURCE-TO-SOURCE

TRANSFORMATIONS

We have shown in previous Section that Vivado HLS can be used to synthesize very efficient specialized floating point operators which rival in quality with those generated by FloPoCo. Our goal is now to study how such optimization can benefit from automation. More precisely, we aim at being able to automatically optimize Listing 1 into Listing 3, and generalizes this transformation to many more situations.

For convenience, we chose to develop our optimization as a source-to-source transformation implemented within the open source GeCoS compiler framework [START_REF] Floc'h | GeCoS: A framework for prototyping custom hardware design flows[END_REF], and aim at making our tool publicly available. Source-to-source compiler are very convenient in a HLS context since they can be used as optimization front-ends on top of closed-source commercial tools.

This work focuses on two computational patterns, namely the accumulation and the sum of product. Both are specific instances of the reduction pattern, which can be optimized by many compilers or parallel run-time environments. Reduction pattern are exposed to the compiler/runtime either though user directives (e.g #pragma reduce in openMP), or automatically inferred using static analysis techniques [START_REF] Redon | Detection of scans in the polytope model[END_REF], [START_REF] Doerfert | Polly's polyhedral scheduling in the presence of reductions[END_REF].

As the problem of detecting reductions is not the main focus on this work, our tool uses a straightforward solution to the problem using a combination of user directive and (simple) program analysis. More specifically, the user must identify a target accumulation variable through a pragma, and provide additional information such as the dynamic range of the accumulated data along with the target accuracy (in the future, we expect to automate our flow such that the two later parameter could systematically be omitted).

We found this approach easier, more general and less invasive than those attempting to convert a whole floatingpoint program into a fixed-point implementation [START_REF] Sentieys | Automatic Fixed-Point Conversion: a Gateway to High-Level Power Optimization[END_REF].

A. Proposed compiler directive

In imperative languages such as C, reductions are implemented using for or while loop constructs. Our compiler directive must therefore right outside such a construct. Listing 4 illustrates its usage on the code of Listing 1.

The pragma must contain the following information:

• The keyword FPacc, which triggers the transformations.

• The name of the variable in which the accumulation is performed, preceded with the keyword VAR. In the example, the accumulation variable is acc.

• The maximum value that can be reached by the accumulator through the use of the MaxAcc keyword. This value is used to determine the weight MSB A . • The desired accuracy of the accumulator using the epsilon keyword. This value is used to determine the weight LSB A . • Optional: The maximum value among all inputs of the accumulator in the MaxInput field. This value is used to determine the weight MaxMSB X . If this information is not provided, then MaxMSB X is set to MSB A . In the case when no size parameters are given, a full Kulisch accumulator is produced. Also note that the user can quietly overestimate the maximum value of it's accumulator without major impact on area. For instance, overestimating MaxAcc by a factor 10 only adds 3 bits to the accumulator width.

B. Proposed code transformation

The proposed transformation operates on the compiler program intermediate representation (IR), and rely on the ability to identify loops constructs and expose def/use relations between instructions of a same basic block in the form of an operation dataflow graph (DFG).

To illustrate our transformation, consider the sample program shown in Listing 5. This program performs a reduction into the variable sum, involving both sums and sums of product operations. The operation dataflow graph associated to the basic block forming the loop body in this program is depicted in Figure 4a. In this Figure, dotted arrows represent loop carried dependencies between operations belonging to distinct loop iterations. Such loop carried dependencies have a very negative impact on the kernel latency as they prevent loop pipelining. For example, when using a pipelined floating-point adder with a seven cycle latency, the HLS tool will schedule a new iteration of the loop at best every seven cycle.

As illustrated in Figure 5a, our proposed optimization hoists the floating-point normalization step out of the loop, and performs the accumulation using fixed point arithmetic. Since integer add operations can generally be implemented with a 1-cycle delay, the HLS tool may now be able to initiate a new iteration every cycle, improving the overall latency by a factor of 7.

; i++){ sum+=in1[i] * in2[i-1]; sum+=in1[i]; sum+=in2[i+1]; } return sum; }
The code transformation first identifies all relevant basic blocks (i.e those associated to the pragma directive). It then performs a backward traversal of the dataflow graph, starting from a Float Add node that writes to the accumulation variable identified by the #pragma.

During this traversal, the following actions are performed depending on the visited nodes:

• A node with the summation variable is ignored.

• A Float Add node is transformed to an accurate fixedpoint adder. The analysis is then recursively launched on that node. • A Float Mul node is replaced with a call to the ExactProduct function followed by a call to ExactProdFloatToFix. • Any other node has a call to FloatToFix inserted. This algorithm rewrites the DFG from Figure 4a into the new DFG shown on Figure 5a. In addition, a new basic block containing a call to FixToFloat is inserted immediately after the transformed loop, in order to expose the floating point representation of the results to the remainder of the program.

From there, it is then possible to regenerate the corresponding C code. As an illustration of the whole process, the synthesized codes from before and after the transformations result in the architectures from Figure 4b and Figure 5b respectively.

C. Evaluation of the toy example of Listing 5

The proposed transformations work on non-trivial examples such as Listing 5. Table II shows how resource consumption depends on epsilon, all the other parameters being those given in the pragma of Listing 5. All these versions where synthesised for 100 MHz. Indeed, some circuits might need more work than others to increase their frequency and would not give a fair comparison.

in2[i-1] in1[i] Float Mul Float Add in2[i+1] Float Add Float Add sum sum (a) Loop body dataflow graph in2[i-1] in1[i] in2[i+1]
Our transformed code makes Vivado HLS use more LUTs for less DSPs compared to the classical IEEE-754 implementation. This is due to the smaller shifter and multiplier of our implementation. In all cases, on this example, the transformed code has its latency reduced by a factor 20. This is due to the fact that Vivado HLS is able to perform the Float Mul and the Float Adder in a single block with a latency of 10 cycles for the naive version. The two following Float Adder are accelerated from 14 cycles to 10 cycles due to a very short pipeline that follows the IEEE-754 standard. Comparatively, the transformed code has 1 cycle latency operators that can be pipelined.

IV. EVALUATION ON FPMARK BENCHMARKS

In order to evaluate the relevance of the proposed transformations on real-life programs, we used the EEMBC FPMark Table II: Comparison between the naive code from Listing 5 and its transformed equivalent. All these versions run at 100MHz.

in2[i-1] in1[i] Exact Product FloatToFix ExactProduct FloatToFix Fixed Add in2[i+1] FloatToFix Fixed Add Fixed Add long accumulator long accumulator (a) Loop body dataflow graph in2[i-1] in1[i] in2[i+1] Exact Product FloatToFix ExactProduct FloatToFix Fixed Add FloatToFix Fixed Add
benchmark suite [START_REF]EEMBC -the embedded microprocessor benchmark consortium[END_REF]. This suite consists of 10 programs. A first result is that half of these programs contain visible accumulations:

• Enhanced Livermore Loops (1/16 kernels contains one accumulation) • LU Decomposition (multiple accumulations) • Neural Net (multiple accumulations)

• Fourier Coefficients (one accumulation)

• Black Scholes (one accumulation)

The following focuses on these, and ignores the other half (Fast Fourier Transform, Horner's method, Linpack, ArcTan, Ray Tracer).

Most benchmarks come in single-precision and doubleprecision versions, and we focus here on the single-precision ones.

A. Benchmarks and accuracy: methodology

Each benchmark comes with a golden reference against which the computed results are compared. As the proposed transformations are controlled by the accuracy, it may happen that the transformed benchmark is less accurate than the original. In this case, it will not pass the benchmark verification test, and rightly so.

A problem is that the transformed code will also fail the test if it is more accurate than the original. Indeed, the golden reference is the result of a certain combination of rounding errors using the standard FP formats, which we do not attempt to replicate.

To work around this problem, each benchmark was first transformed into a high-precision version where the accumulation variable is a 10,000-bit floating-point numbers using the MPFR library. We used the result of this highly-accurate version as a "platinum" reference, against which we could measure the accuracy of the benchmark's golden reference. This allowed us to choose our epsilon parameter such that the transformed code would be at least as accurate as the golden reference. This way, the epsilon of the following results is obtained through profiling. The accuracy of the obtained results are computed as the number of correct bits of the result.

We first present the benchmarks that are improved by our approach before discussing the reasons why we can't prove that the others are.

B. Benchmarks improved by the proposed transformation

Enhanced Livermore Loops: This program contains 16 kernels of loops that compute numerical equations. Among these kernels, there is one that performs a sum-of-product (banded linear equations). This kernel computes 20000 sumsof-products. The values accumulated are pre-computed. This is a perfect candidate for the proposed transformations.

For this benchmark, the optimal accumulation parameters were found as: MaxAcc=50000.0 epsilon=1e-5 MaxInput=22000.0

Synthesis results of both codes (before and after transformation) are given in Table III Vivado HLS is unable to predict a latency for these implemented designs due to their non-constant loop counts, therefore we do not present complete results for these two benchmarks. Still, in order to show that the approach works on these examples, the LU inner loops were transformed and synthesized. Table III shows the results obtained for the smallest (8 terms) and the largest (45 terms) sums-of-products in lines LU-8 and LU-45 respectively. The latency is vastly improved even for the smallest one. The accuracy results of the original code here varies from 8 to 23 bits between different instances of the loops. To have a fair comparison, we generated a conservative design that performs 23 bits accuracy on all loops, using a sub-optimal amount of resources.

C. Benchmarks that exposed the limitations of HLS tools

Black Scholes: This program contains an accumulation that sums 200 terms. The result of this computation is divided by a constant (that could be optimized by using transformations based on [START_REF] De Dinechin | Table-Based Division by Small Integer Constants[END_REF]). This process is performed 5000 times.

Here the optimal accumulator parameters are the following: MaxAcc=245000.0 epsilon=1e-4 MaxInput=278.0 This gives us an accumulator that uses 19 bits for the integer part and 10 bits for the fractional part. The result of the synthesis are provided in Table III.

For comparable area, accuracy is vastly improved but latency could not be obtained from Vivado HLS. Indeed, the Black Scholes algorithm uses the mathematical function power. Such a function is not implemented in Vivado HLS, therefore we coded it using a loop with a non-constant count. As the latency of this operator depends on the input data, the latency of the all circuit cannot be statically recovered.

Fourier Coefficients: The Fourier coefficients program, which computes the coefficients of a Fourier series, contains an accumulation which is performed in single precision. This program comes in three different configurations: small, medium and big. Each of them computes the same algorithm but with a different amount of iterations. The "big" version is supposed to compute the most accurate answer. We obtain similar results for the three versions of this program, as a consequence we only present the "big" version here. In this version, there are multiple instances of 2K terms accumulations. The accumulator is reset at every call.

The parameters determined for this benchmark were the following: MaxAcc=6000.0 epsilon=1e-7 MaxInput=10.0

This results in an accumulator using 14 bits for the integer part and 24 bits for the fractional part. The synthesis results obtained for the original and transformed codes are given in Table III.

Here, area is again comparable, accuracy is improved by 5 bits (which represents one order of magnitude), but latency could not be obtained. Vivado HLS faces the same problem as in Black Scholes, it cannot compute the latency of the power function.

Note that our operators have a shorter latency by nature. Therefore, even if Vivado HLS is not able to provide latency results, the circuits indeed have a shorter latency.

The pow operator could be implemented such as in [START_REF] De Dinechin | Floating-Point Exponentiation Units for Reconfigurable Computing[END_REF] in the near future. This would allow us to obtain the latency of the two last benchmarks.

V. CONCLUSION

The main result of this work is that HLS tools have the potential to generate efficient designs for handling floatingpoint computations in a completely non-standard way. The use of application-specific intermediate formats can provide both performance and accuracy at a competitive cost. For this, we have to sacrifice the strict respect of the IEEE-754 and C11 standards. It is replaced by the strict respect of a high-level accuracy specification.

Classically, designers have to face a trade-off between performance and cost. This approach adds computation accuracy to this trade-off. Some designers may not like this. To convince them, consider that established performance benchmarks compute results which are accurate only to a few bits. If only a few bits are important, do we really need to instantiate 32bit or 64-bit floating-point operators to compute them ? Isn't this accuracy information worth investigating and exploiting?

This work also provides a practical tool that improves a given C program. The input to the tool is application-specific information representing high-level domain knowledge such as the range and desired accuracy of a variable. The resulting code is compatible with Vivado HLS.

The proposed transformation already works very well on all the FPMarks that contains a reduction where it improves both latency and accuracy by an order of magnitude for comparable area.

In the longer term, we believe there is much more to come. The arithmetic optimizations that a classical compiler can do are very limited by the fixed hardware of classical processors. With compilers of high-level software to hardware, there is much more freedom, hence many more opportunities to build application-specific arithmetic operators.

Future work will attempt to explore this new realm, starting with operator specialization; operator fusion such as in [START_REF] Ye | MixFX-SCORE: Heterogeneous fixed-point compilation of dataflow computations[END_REF] but at a more coarse grain allowing more aggressive fusion; compile-time generation of application-specific cores; error analysis such as in [START_REF] Langhammer | FPGA floating point datapath compiler[END_REF] benefiting from compilers static analysis and more generaly building upon compiler progresses in program analysis.

Figure 1 :

 1 Figure 1: The proposed compilation flow

Listing 1 :

 1 Naive reduction #define N 100000 float acc = 0; for(int i=0; i<N; i++){ acc+=in[i]; } Listing 2: Parallel reduction #define N 100000 float acc = 0, tmp1=0, ... , tmp10=0; for(int i=0; i<N; i+=10){ tmp1+=in[i]; ... tmp10+=in[i+9]; } acc=tmp1+...+tmp10;

Figure 2 :

 2 Figure 2: The conversion from float to fixed-point (top), the fixed-point accumulation (middle) and the conversion from the fixed-point format to a float (bottom).

Listing 3 :

 3 Sum of floats using the large fixed-point accumulator #define N 100000 float acc = 0; ap_int<68> long_accumulator = 0; for(int i = 0; i < N; i++) { long_accumulator += FloatToFix(in[i]); } acc = FixToFloat(long_accumulator);

Listing 4 :

 4 Illustration of the use of a pragma for the naive accumulation #define N 100000 float accumulation(float in[N]){ float acc = 0; #pragma FPacc VAR=acc MaxAcc=100000.0 epsilon=1E-15 MaxInput=1.0 for(int i=0; i<N; i++){ acc+=in[i]; } return acc; }

Listing 5 :

 5 Simple reduction with multiple accumulation statements #define N 100000 float computeSum(float in1[N], float in2[N]){ float sum = 0; #pragma FPacc VAR=sum MaxAcc=300000.0 epsilon=1e-15 MaxInput=3.0 for (int i=1; i<N-1

Figure 4 :

 4 Figure 4: DFG of the loop body from Listing 5 (top) and it's corresponding architecture (bottom). Keywords float mul and float add correspond to floating-point multipliers and adders respectively

Figure 5 :

 5 Figure 5: DFG of the loop body from Listing 5 (top) and it's corresponding architecture (bottom) after transformations

Table I :

 I Different ways of implementing a simple accumulation.

		Listing 1	Listing 2	Listing 1	Listing	Listing 3	FloPoCo VHDL
		(float)	(float)	(double)	(double)	(71 bits)	(71 bits)
	LUTs	266	907	801	2193	736	719
	DSPs	2	4	3	6	0	0
	Latency	700K	142K	700K	142K	100K	100K
	Accuracy	17 bits	17 bits	24 bits	24 bits	24 bits	24 bits

Table III :

 III . As in the previous toy examples, latency and accuracy are vastly improved for comparable area.LU Decomposition and Neural Net: Both the LU decomposition and the neural net programs contain multiple nested small accumulations. In the LU decomposition program, an inner loop accumulates between 8 and 45 values. Such accumulations are performed more than 7M times. In the neural net program, inner loops accumulate between 8 and 35 values, and such accumulations are performed more than 5K times. Synthesis results of benchmarks before and after transformations Both of these programs accumulate values from registers or memory that are already computed. It makes these programs good candidates for the proposed transformations.

	Benchmark	Type	LUTs	DSPs Latency Accuracy
	Livermore	Original	384	5	80K	11 bits
		Transformed	576	2	20K	13 bits
	LU-8	Original	809	5	82	8-23 bits
		Transformed	1007	2	17	23 bits
	LU-45	Original	819	5	452	8-23 bits
		Transformed	1034	2	54	23 bits
	Scholes	Original	15640	175	N/A	19 bits
		Transformed 15923	175	N/A	23 bits
	Fourier	Original	34596	64	N/A	6 bits
		Transformed 34681	59	N/A	11 bits