Alexandra Lelong

Pierre Guiffant

Jacques André Astolfi
email: jacques-andre.astolfi@ecole-navale.fr

S Y Mpo

O Tatin

An Experimental Analysis of the Structural Response of Flexible Lightweight Hydrofoils in Various Flow Conditions

Keywords:

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

In hydraulic rotating machines, the perspective of using lightweight structures with a certain degree of flexibility is a challenging way to reduce the weight, to enhance the performance and opens possibility to control the machine. For instance an improvement can come from the use of flexible lifting structures that could adapt passively or actively to variable operating conditions. However using flexible lightweight structural in water gives rise to basic questions related to Fluid Structure Interaction in a high density fluid that is quite different from aeroelasticity for which the fluid density is much smaller than the mean structure density in many applications as aerospace or wind engineering. For hydraulic rotative machine operating at relative high velocity, structures are subjected to high loadings due to the high density of water and relative large velocity. Using nonmetallic and lighter materials can lead to high stresses and significant deformation inducing strong coupling with the flow. Moreover from a dynamical point of view, several points need to be well understood as the dynamic structural response related to added-mass effects, damping, stiffness together with possible flow-structure instabilities in turbulent high Reynolds and potentially cavitating flows. The latter is of primary importance for hydraulic machine design . The control of cavitation occurring as the minimum pressure in the fluid is lower than the va-por pressure ( [START_REF] Brennen | Cavitation and bubble dynamics[END_REF], [START_REF] Franc | Physics and control of cavitation[END_REF])is a very challenging task and using active or passive flexible structures is a way that should be of great interest. In hydraulic machine, operating flow conditions induce strong low pressure on the suction side of blades leading to various kinds of cavitation according to the angle of attack ( [START_REF] Astolfi | An experimental investigation of cavitation inception and development on a two-dimensional hydrofoil[END_REF]). For instance, it is known that for small angles of incidence, travelling bubble occur on the suction side and that leading edge cavitation (namely partial or sheet cavitation) occurs for larger angles of incidence. Moreover under specific conditions, the vaporized area becomes unstable, with cavity break-off and periodical shedding of large bubble clusters. This configuration, named "cloud cavitation", generates pressure fluctuations downstream of the cavity that was studied largely on rigid structures ( [START_REF] Leroux | A joint experimental and numerical study of mechanisms associated to instability of partial cavitation on two-dimensional hydrofoil[END_REF]) but is unknown when coupling with flexible structures. Recent efforts were made in order to approach through simulations both the fluid and the structure dynamics. This comes from the recent development of numerical methods and models which allow to compute complex flows together with coupled simulation with reasonable CPU time, (see [START_REF] Young | Fluid-structure interaction analysis of flexible composite marine propellers[END_REF]). Some numerical developments were developed to analyze fluid-structure interaction (FSI) of flexible composite structure in hydraulic applications mainly for marine propellers ( [START_REF] Sontvedt | Propeller blade stresses, application of finite element methods[END_REF], [START_REF] Lin | Strength evaluation of a composite marine propeller blade[END_REF], [START_REF] Young | Hydroelastic behavior of flexible composite propellers in wake inflow[END_REF], [START_REF] Young | Fluid-structure interaction analysis of flexible composite marine propellers[END_REF], [START_REF] Young | Hydroelastic response of wind or tidal turbines[END_REF], [START_REF] Mulcahy | Flexible composite hydrofoils and propeller blades[END_REF], [START_REF] Motley | Performance-based design and analysis of flexible composite propulsors[END_REF]). However the numerical developments suffer from a lack of experimental data. For instance, although a large number of experimental studies deals with the analysis of cavitation over rigid hydrofoils ([4], [START_REF] Coutier-Delgosha | Analysis of cavitating flow structure by experimental and numerical investigations[END_REF]) little referred to experimental cavitation on flexible structures. Ausoni at al. ( [START_REF] Ausoni | Cavitation influence on von karman vortex shedding and induced hydrofoil vibrations[END_REF]) studied experimentally the vibration response of a metallic hydrofoil but only recently, cavitation over a flexible hydrofoil has been studied numerically and experimentally ( [START_REF] Ducoin | Numerical and experimental investigation of hydrodynamic characteristics of deformable hydrofoils[END_REF], [START_REF] Reese | Vibration and damping of hydrofoil in uniform flow[END_REF], [START_REF] Ducoin | An experimental analysis of fluid structure interaction on a flexible hydrofoil in various flow regimes including cavitating flow[END_REF], [START_REF] Gaugain | Numerical and experimental study of hydroelastic behaviour of a hydrofoil[END_REF], [START_REF] Akcabay | Influence of cavitation on the hydroelastic stability of hydrofoils[END_REF] and [START_REF] Akcabay | Cavity induced vibration of flexible hydrofoils[END_REF]) or theoretically ( [START_REF] Benaouicha | A numerical study of cavitation induced vibration[END_REF], [START_REF] Benaouicha | Analysis of added mass in cavitating flow[END_REF]). Theses studies deal with the vibration induced by the cavitation on a 3D hydrofoil or a two degrees of freedom model of a rigid hydrofoil in pitching and heaving motion. They indicated that cavitation can change the modal response of the structure and that a lock-in of the foil's frequencies in cavitating flows can occur in specific conditions but that experimental data and physical observations are still required for a better understanding of a complex FSI phenomena to be beneficial for simulation of FSI in hydraulic applications.

The present paper deals with a collaborative research program between the French Naval Academy and the Department of Naval Architecture and Marine Engineering at the University of Michigan in order to analyze the hydroelastic response of a flexible homogeneous hydrofoil in various flow conditions including cavitating flow. The objective of the program is to perform physical analysis as well as to improve simulations of FSI for heavy fluid applications.

EXPERIMENTAL SETUP

The experiments were carried out in the cavitation tunnel of Naval Academy Research Institute (IRENav) fitted with a 192 mm square test section of 1m long (Figure 1). The flow velocity ranges between 3 and 12m.s-1 and the pressure in the tunnel test section is between 3 bars to 100 mbar mainly to control cavitation inception. The inlet turbulence intensity measured about two chords from the hydrofoil leading edge by LDV is close to 2%. Two NACA 0015 flexible foils of c = 100mm chord and b = 191mm span were used. The first one was used for strain and vibration measurements. It was constituted of a cylindrical base that insure a quasi-perfect clamped condition at the foil root forming a cantilevered flexible hydrofoil. It was set into a cylindrical rigid aluminium base (Figure 2). The foil was obtained by milling a Polyoxymethylene plastic (POM, Table 1). The other foil was used to measure static deformation. In that case, the flexible part was fabricated using stereolithography 3D printing using POLY1500 resin. In both case the rotation axis was located at x/c = 0.5 from the leading edge. The characteristics of the material are given in the Table 1 and compared with the steel's ones for information.Note the uncertainties on the POLY1500 resin. The incidence of the foil in the tunnel is controlled by a Baldor drive system that allows us to control with high accuracy the angle step ∆α, the acceleration and the rotational speed. The theoretical accuracy of the adjustment is 6.10 -4 degree, however the uncertainties of the angle of attack was founded to be 0.1 • , considering the mechanical mounting system. A high speed camera was fixed over the test section to record movies when cavitation develops over the foil. The model used is a Fastcam SA3 120K from Photron. It is equipped with a CMOS image sensor, which sensor size is 17.4mm*17.4mm. The frame rate is spread from 60fps to 2000fps for the full sensor resolution (1024*1024 pixels), when the pixel size is 17µm.

Static deformation measurements

The static deformation was measured using a Laser distance measurement system mounted on a 2D translation system on the uper side of the test section. The system allowed us to scan the hydrofoil surface for a given flow condition along various sections selected along the span.

In that case five sections from the root to the tip was selected. Then the deformation was obtained by comparing the scan without flow and the scan obtained at a given velocity. With this method it was possible to extract the bending and twisting of the hydrofoil for a given angle of incidence and velocity. However, a specific geometrical system was used to calibrate the laser distance measurement in order to take into account of the deflection of the beam laser crossing the Plexiglas and the water layers.

The calibration consisted to measure the position of a well-known geometrical system in air then to measure the same system position with water and Plexiglas layers.

A correction factor was then found in such a way the two mesaurements give the same results. The correction factors were aplied to the measurements.The results are presented on Figure 3 for the static bending at mid-chord and on Figure 4 for twisting. As shown the hydrofoil experienced mostly bending ranging from 1mm to 10mm than twisting ranging within 0.05 • and 0.1 • according the angle of incidence. Bending (mm) Span (mm) Twisting (degrees)

0 • 2 • 4 • 6 • 8 • 10 •
0 • 2 • 4 • 6 • 8 • 10 •

Strain measurements

To determine components of plane strain, the hydrofoil is equipped with three identical strain gauges build in a rosette way. They are glued in a hole close to the root of the foil. The cavity is filled with a resin whose mechanical characteristics are close to the foil's material. The cavity is then polished with sand paper to obtain a smooth continuous surface (Figure 2).

The strain gauges are L2A 13-125WW-120 from Vishay Micro Measurements. They are made in constantan and have a K-factor equal to 2.11. Their accuracy is ±6%. The strain gauges are assembled in a Quarter Bridge . The strain gauges measure the strains (ε 1,2,3 in the direction of their own axis. The gauges' axes are separated by a 45 • angle and the first one is given by the span direction (ε 1 ), the third by the chord direction ((ε 3 ). The data are first amplified (Quantum mx16, HBM) and then recorded by the acquisition software CATMAN (HBM).

The principal strains ε I,II can be calculated by:

ε I,II = ε 1 + ε 3 2 ± 1 2 (ε 3 -ε 1 ) 2 + 4 * (ε 2 - 1 2 (ε 1 + ε 3 )) 2 (1)
as well as the principal stresses σ I,II :

σ I,II = λ(ε I + ε II ) + 2µε I,II (2)
with λ = Eν (1+ν)(1-2ν) and µ = E 2(1+ν) the Lamé parameters.

The Von Mises stress can be computed as well by:

σ V M = σ 2 I + σ 2 II -σ I σ II (3)

Vibration measurements

The measurements were performed with two vibrometers from Polytec: the first one is fixed and is taken as a reference whereas the second one is a PSV-400 scanning vibrometer. This model can detect vibrations up to 10m.s -1 with a HeNe laser (λ = 633nm). It is equipped with two analog velocity decoder VD-04 and VD-06. The scanner is a high precision scan unit, with an angular resolution lower than 0.002 • and an angular stability lower than 0.01 • per hour. Because of laser light diffusion in the platic material, reflecting tapes were glued on the foil's surface to enhance the signal to noise ratio. It allows us to measure the vibration level on a user-defined grid over the structure surface. The cross-spectrum between the reference laser point and a scanned measured point is computed to preserve the phase at a given frequency. From this phase, the phase between two scanned points measured at two different time could be computed. This allowed us to get the vibration shape of the structure at a given frequency. It is very convenient to infer the modal shape associated to a modal frequency.

An electrodynamic shaker was used to measure the response of the foil to an impulse and identify the natural frequencies of the foil in air and in still water. Eight measurements were performed on each point, with one impulse per measurement, and the mean spectrum was calculated. The frequency resolution of the following measurements is ∆f = 0.625Hz. The repeatability of the experiment was tested and the accuracy was then lower than 2%.

The modes have been identified thanks to the phase preservation between the reference and the scanning vibrometer. The spectra are shown on figure 5 and the values of the three first modes' frequencies are given in Table 2. As shown, the frequencies of the bending (f 1 ) and twisting (f 2 ) modes in water decrease strongly compared to in the air as a result of added mass effects. The modal shapes are given on Figures 6(a) and 6(b) for bending and twisting. The mean bending and twisting along the chord for various sections is compared to a beam theoretical bending and twsting modes showing that the 3D hydrofoil has a beam like behavior for the two first modes of vibration (Figures 7 and8)

f 1 f 2 f 3 f 1 f 2 f 3 Frequency (Hz) Magnitude (m.s -1 )
air water

STRAINS AND STRESSES ANALYSIS

Non cavitating flow

Concerning strains' measurements, the experiment was performed from -10 • to 10 • and with a step ∆α = 0.5 The results are shown on Figure 9. Again it is observed that the hydrofoil behaves like a built-in beam: ε 1 and ε 2 increase with the angle of attack. The higher strain is ε 1 , which corresponds to the span direction (blue line on Figure 9). It is due to the hydrodynamic lift force. The strain ε 3 , in the chord direction (green line), is opposed to the others because of the twist and the shearing stress. It is observed that the strains are linear up to 3 • . Then they discard from the linear trend. This is particularly true for ε 1 that is highly related to the lift. Figure 10 shows the principal stresses (2) and the Von Mises stresses (3). As shown, an inversion of the principal stresses' curves appears passing through 0 • . This is due to a change of the principal directions. The Von Mises stress is symmetrical and it is linear up to 4 • (respectively -4 • for negative angle of incidence). Beyond 4 • or -4 • it discards from the linear evolution.This is a consequence of the lift force evolution that was found to evolve in the same way as a result of boundary layer transition due to a Laminar Separation Bubble at a moderate chord length Reynolds number [START_REF] Delafin | Effects of the laminar separation bubble induced transition on the hydrodynamic performance of a hydrofoil[END_REF]. It was observed that triggering the boundary layer at the leading edge suppress this peculiar nonlinear behavior.

σ I σ II σ V M

Cavitating flows

Measurements have been carried out with a decreasing number of cavitation in order analyze the behavior of the strains and stresses in cavitating flows. It was observed that depending on the strain direction the mean values of the strains ε i increases with σ (see ε 1 ), and decreases for the lowest cavitation number as it is illustrated on Figure 11. In the same time, the strain fluctuations increase as shown by the vertical bars, determined by the standard deviation. This is particularly true for ε 1 directly related to the lift forces. The same observation is clearly shown on the Von Mises stress evolution as function of the cavitation number (Figure 12). As an interesting result, it is shown that strain gauges' response can be a very interesting way to get information on forces acting on the structures. It could be a good way for control lifting surfaces in operating condition. In fact as the cavitation number decreases, the strain signals experience large fluctuations at a frequency related to the oscillation of the cavity. The frequency, named f c , decreases as the cavity increases. This is clearly shown on Figure 13, which represents the time series of the Von Mises stresses depending on the cavitation number. It is particularly clear on the Von Mises stresses' signal for σ = 1.70 (Figure 13). The frequencies of the strain gauge's signals were found to be similar to the vibration's frequencies of the foil reported in the next session.

FOIL'S VIBRATION

Non cavitating flows

Measurements were carried out at a constant pressure in tunnel section close to the atmospheric pressure, for which cavitation does not develop. They were performed for angles of incidence from 0 • to 8 • with a step ∆α = 2 • , the flow velocity was 3, 4, 5 and 6 m.s -1 , corresponding to Reynolds numbers based on the chord length ranging from 3.10 5 to 6.10 5 . The frequencies of the first two modes, which correspond to bending and twisting modes, are reported on Table 3 for each incidence and each flow velocity. An example of vibration spectra is given on Figure 14, which represents the spectra with α = 8 • depending on the flow velocity.

The first mode's frequency (bending) is constant with the flow speed and the angle of attack, whereas the second (twisting) and the third one's increase with the flow speed (Figure 14). For the lowest incidences (0 to 4 • ) and low flow velocity (3m.s -1 ), a high peak appears near the twisting mode's frequency (Figure 15). As mentioned previously it is a consequence of the unsteadiness of the Laminar Separation Bubble (LSB) inducing transition at low angle of attack and moderate flow velocities developping on the rear part of the suction side ( [START_REF] Ducoin | An experimental study of boundary layer transition induced vibrations on a hydrofoil[END_REF]). This peak disappear when the incidence or the velocity increases.

It is explained by the sudden displacement of the LSB and transition towards the leading edge [START_REF] Delafin | Effects of the laminar separation bubble induced transition on the hydrodynamic performance of a hydrofoil[END_REF].

Cavitating flows

To study the behavior of the flexible foil with cavitation, measurements have been carried out by decreasing the pressure. The cavitation number is defined by

σ = (P 0 -P v ) 0.5ρU 2 (4)
where P 0 is the pressure at middle of the test section, P v the vapor pressure at the water temperature, U the flow velocity and ρ the water density. The cavitation was controlled by decreasing or increasing the pressure P 0 at a fixed velocity.

Tests have been performed with a fixed flow velocity (U = 6 m.s -1 ) and two angles attack (α 8, 10 • ). For α equal to 8 • , eight measurements have been performed, that corresponds to cavitation numbers ranging between σ = 1.52 and σ = 5.68 (wetted flow). For α equal to 10 • , six measurements have been carried out with σ decreasing from 2.51 to 1.52. The spectra are shown on Figure 16 and Figure 17. As shown on Figure 17, the frequencies of the twisting and the third modes increase as σ decreases. On the contrary, even it is not very significant, the bending mode's frequency tends to decrease slightly as the cavitation develops. For the lowest cavitation numbers (σ lower than 1.41), a low frequency peak and its harmonics appear. For α = 8 • , the first peak corresponds to the bending mode's frequency and is constant until σ = 1.6, then it falls to 8Hz. The same trend is visible on the curves with α = 10 • , but the fall seems to occur later.

To study more deeply this phenomenon, high speed movies was recorded to observe the cavitating flow developing on the suction side. It was observed that for α equal to 8 • and σ = 2.42, attached bubbles appear on the foil mainly due to surface irregularities (Figure 18). As the pressure decrease, a sheet cavitation attached to the leading edge takes places and oscillate between 20 and 50% of the chord for σ = 2.08. In that case, the frequency of cavity oscillation is close to the bending mode's frequency of the foil that implies a lock-in between both frequencies and a strong increase of the magnitude of the peak. A second peak corresponding to the harmonic of the cavity oscillation is observed as cavitation develops more.

For σ = 1.81 and 1.63 (Figure 19), the cavity oscillation amplitude becomes very large. The cavity oscillates from 0% to 70% of the chord length before being carried away downstream. The frequency of the cavity oscillation, noted f o , can be deduced from the high speed camera records between two periods of the cavity evolution. The frequencies correspond to the higher peaks of the frequency responses, that explains the decrease of the first peak's frequency from the bending mode frequency (around 34.4Hz) to 30Hz. Harmonics of this frequency are observed for σ 1.63.

As the cavitation number becomes lower than σ = 1.41, a peak appears at f = 8Hz (Figure 20) and subsequent harmonics are observed. This low frequency peak corresponds to the oscillation frequency of the cavity. In that case the vapor cavity oscillates from fully wetted to nearly 100% of the chord length. Its higher order harmonics are relatively low but can interact with the bending 30Hz-peak. As the pressure decreases again (σ = 1.37 and 1.25), this interaction disappears and the first harmonic becomes stronger.

The Table 4 gathers the oscillation frequencies f o obtained from high speed visualization and the main frequencies of the foil vibration spectrum, called f m , according the cavitation number. The same phenomenon is observed with α = 10 • . For σ = 2.51 and σ = 2.07, a sheet cavitation oscillates between 20% and 40%, and between 0 and 50% respectively, but no frequency appears distinctly on the films. Then we can observe a growth of the cavity until the shedding with a measurable frequency value of approximately 30.3Hz for σ = 1.75 and 33Hz for σ = 1.59, very close to bending frequency. The experiments were performed in the cavitation tunnel of the French Naval Academy. The hydrofoil was equipped with strain gauges embedded close to the root of the foil. Vibrations measurements were carried out using vibrometers. Many new observations can be reported. It was observed that the cantilevered rectangular hydrofoil has a beam-like behavior with bending and twisting deformation mainly. The strains evolve linearly with the angle of incidence as the hydrodynamic loading. In partial cavitating flow, the average value of the strains and stresses tend to increase in the first stage of cavitation together with an increase of fluctuations. As the vapor cavity increased again, the strains and stresses felt down together with a strong increase of fluctuations. This was clearly observable on the Von-Mises stress measurement. Concerning vibrations, the bending and mode's frequencies were clearly identified thanks to impulse tests on the foil's surface in the air and in still water. The effect of added-mass in heavy fluid is clearly determined by The measurements carried out in unsteady cavitating flow allowed us to notice that the bending mode's frequency changes slightly depending on the oscillation frequency of the cavity, resulting from a frequency lock-in with the cavity oscillation's frequency and harmonics or subharmonics . Moreover, the twisting mode's frequency tends to increase with the cavity length that could be related to a decrease of the added mass effects in presence of a vapor cavity on the foil surface. This work brings new informations about the physics of FSI on lightweight flexible structures in a heavy fluid that should be very interesting for physical analysis and simulation of relative complex flow in FSI. This work has to be followed with different enhancements. In particular, some observations need to be studied more deeply. The hydrodynamic force measurements on flexible foil should be also examined. Moreover, tests should be performed with other foils made of other flexible materials.

Figure 1 .

 1 Figure 1. Cavitation tunnel of IRENav.

Figure 2 .

 2 Figure 2. Cantilevered flexible hydrofoil equipped with strain gauges at the root

Figure 3 .

 3 Figure 3. Deformation of the POLY1500 hydrofoil with U = 5m.s -1 -bending.

Figure 4 .

 4 Figure 4. Deformation of the POLY1500 hydrofoil with U = 5m.s -1 -twisting.

Figure 5 .

 5 Figure 5. Response of the foil in POM to an impulse in m.s -1 .

Figure 6 .

 6 Figure 6. Modal shapes.

Figure 8 .

 8 Figure 8. Comparison of the theoretical and experimental modal shapes -twisting.

Figure 9 .

 9 Figure 9. Strains (µm/m) depending on the AoA, U = 5m.s -1 .

Figure 10 .

 10 Figure 10. Stresses (MPa) depending on the AoA, U = 5m.s -1 .

Figure 13 .

 13 Figure 13. Time-history of Von Mises stresses for different cavitation numbers.

Figure 14 .

 14 Figure 14. Spectra with α = 8 • , U = 3 to 6 m.s -1 at atmospheric pressure (non cavitating flow).

- 1 Figure 15 .

 115 Figure 15. Spectra with α = 2 • , U = 3 to 6 m.s -1 at atmospheric pressure (non cavitating flow).

Figure 16 .

 16 Figure 16. Response of a foil in cavitating flow, α = 10 • , U = 6m.s -1 .

Figure 17 .

 17 Figure 17. Vibration spectra in cavitating flow, α = 8 • , U = 6m.s -1 (Spectra are displayed on two subfigures (a) and (b) depending on σ to improve visibility, the green curve (σ = 1.63) is represented on both figures)

Figure 18 .

 18 Figure 18. Vibration spectra in cavitating flow, α = 8 , U = 6m.s -1 (Spectra are displayed on two subfigures (a) and (b) depending on σ to improve visibility, the green curve (σ = 1.63) is represented on both figures)

Figure 19 .

 19 Figure 19. Vibration spectra in cavitating flow, α = 8 • , U = 6m.s -1 (Spectra are displayed on two subfigures (a) and (b) depending on σ to improve visibility, the green curve (σ = 1.63) is represented on both figures)

Figure 20 .

 20 Figure 20. Vibration spectra in cavitating flow, α = 8 • , U = -1 (Spectra are displayed two subfigures (a) and (b) depending on σ to improve visibility, the green curve (σ = 1.63) is represented on both figures)

Table 1 .

 1 Characteristics of the flexible materials. Steel is given for comparison.

		Steel POM	POLY1500
	Young mod. (E, GPa)	203	2.9	1.227-1.462
	Poisson coef. (ν)	0.30	0.35	--
	Density (ρ, kg.m -3 )	8010	1420 1180 -12000

Table 2 .

 2 First modes' frequencies of the foil in POM as response to series of eight impulses (Hz).

		Air Still water
	f 1	80.6	34.4
	f 2	390	183.5
	f 3 556.5	292
	•10 -4
	8
	6
	4
	2
	0
	50 100 150 200 250 300 350 400 450 500 550 600

Table 3 .

 3 First two modes frequencies (Hz) depending on the flow speed and the AoA

				3m.s -1	4m.s -1		5m.s -1	6m.s -1
		8 •		36.5 Hz 182 Hz 188.5 Hz 35.5 Hz		37 Hz 191 Hz	37 Hz 194 Hz
		6 •		37 Hz 185 Hz 189.5 Hz 190.5 Hz 190.5 Hz 36.5 Hz 37 Hz 36.5 Hz
		4 •		37.5 Hz 187 Hz 187.5 Hz 35.5 Hz		36.5 Hz 187 Hz	36 Hz 189 Hz
		2 •		36.5 Hz 186.5 Hz	37.5 Hz 187 Hz 186.5 Hz 189.5 Hz 37 Hz 37 Hz
		0 •		35 Hz 186 Hz 186.5 Hz 37 Hz		36 Hz 190 Hz 184.5 Hz 37 Hz
	Magnitude (dB)	-100 -80								3m.s -1 4m.s -1 5m.s -1 6m.s -1
		-120	50	100	150	200	250	300	350	400
					Frequency (Hz)
						(a) Spectrum

Table 4 .

 4 Oscillation frequencies of the cavity sheet and main frequencies of the foil depending on the cavitation number, α = 8 • , U = 6m.s -1 .

	σ	2.42 2.08 1.81 1.63 1.41 1.37 1.25
	f o	-	-	33 31.4 8.51 8.54 7.84
	f m	35	35	35	31	8	8	8

ACKNOWLEDGMENTS

The authors wish to acknowledge the Office of Naval Research for the support under the grant number N62909-12-1-7076, as well as the technical staff of the Naval Academy Research Institute, Alain Boulch is gratefully acknowledged.


Julie Coloigner

Jean-Marie Batail

Isabelle Corouge

Dominique Drapier

published or not. The documents may come

White matter connectivity analysis in patients suffering

from depression.

Introduction

• Major depressive disorder (MDD) is characterized by a profound dysregulation of affect and mood: cognitive dysfunction, insomnia, fatigue and appetite disturbance.

• Up to 80% of patients will suffer from a relapse.

• Understanding the neural correlates underlying the depression is critical for improving the efficacy of diagnostic and treatment strategies.

• Dysfunction of the circuits connecting frontal and subcortical brain regions, leads to a "disconnection syndrome" [1].

• Using graph theory-based analysis, we examined white matter changes in the organization of networks in patients suffering from depression.

Figure 1. Overview of the method

Silent strokes SCD versus CTL Data:

• LONGIDEP is a routine care cohort of patients suffering from MDD who underwent clinical, neuropsychological testing and imaging study.

• DTI scans: 30 directions using an EPI sequence with a b-value of 1000s/mm 2 .

• 20 MDD patients and 20 healthy controls (CTL) subjects.

Network node definition:

• Parcellation of the cerebral cortex into 66 cortical and 14 subcortical regions using Freesurfer image analysis [2].

• Each regions represents a node of the cortical network

Diffusion tractography:

• Fiber Assignment by Continuous Tracking (FACT) algorithm using Anima (https://github.com/Inria-Visages/Anima-Public/wiki).

Network analysis:

• Computation of various metrics that characterize the global organization of anatomical networks

• Fiber density (number of fibers normalized by ROI volume) characterizes the strength of inter-regional connections.

• Regional nodal properties: Strength (reflects how the node is interacting with many other nodes in the network), Betweenness (in integrating dispersed networks) and Clustering (in local networks).

• Two-sample Student t-test between the two groups with age, gender, duration of disease and medication load [3] as nuisance covariates, combined with a permutation test.

Material and Methods

White matter connectivity analysis in patients suffering from depression

Results

Freesurfer processing stream

Diffusion images Structural images

Skullstripping Distortion correction Denoising

Tractography

Graph analysis Fiber density

• Widespread white matter abnormalities in patients with MDD: reduced connection between the frontal and parietal lobes, and between the frontal pole and limbic regions

• In the frontal pole: reduced strength and reduced clustering

• Higher betweenness in right and left thalamus and in right putamen for MDD

• Decreased fiber density in circuits connecting subcortical brain regions with the frontal and parietal cortex, supporting the theory of limbic-frontal circuit dysfunction in MDD.

• Less segregated network organization in the frontal lobe, implicated in abnormalities of emotion regulation and cognition in MDD.

• Thalamus and putamen: highly interactive regions that likely participate in more functional interactions.

Conclusion

Fiber density

. Regions exhibited significant between-group differences in regional nodal parameters between MDD and control groups. The blue color represented the higher values of regional nodal metrics in MDD group, and the red color represented the higher values of regional nodal metrics in CTL group (p<0.05, with permutation test).

[1] Liao, Y., Huang, X., Wu, Q., Yang, C., Kuang, W., Du, M., .

J. Coloigner 1 , J.-M. Batail 1,2 , I. Corouge 1 , D. Drapier 2 and C. Barillot 1 1 Univ Rennes, Inria, CNRS, Inserm, IRISA UMR 6074, VISAGES ERL U-1228, F-35000 Rennes, France 2 Academic Psychiatry Department, Centre Hospitalier Guillaume Régnier, Rennes, France, EA 4712 Behavior and Basal Ganglia, CHU Rennes, Rennes 1 University, Rennes, France


Ksander Ejjaaouani
email: ksander.ejjaaouani@inria.fr

Olivier Aumage
email: olivier.aumage@inria.fr

Julien Bigot
email: julien.bigot@cea.fr

Michel Mehrenberger

Hitoshi Murai
email: h-murai@riken.jp

Masahiro Nakao
email: masahiro.nakao@riken.jp

Mitsuhisa Sato
email: msato@riken.jp

InKS, a Programming Model to Decouple Performance from Algorithm in HPC Codes

Keywords: programming model, separation of concerns, HPC, DSL

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

It is more and more common to identify simulation as the third pillar of science together with theory and experimentation. Parallel computers provide the computing power required by the more demanding of these simulations. The complexity and heterogeneity of these architectures do however force scientists to write complex code (using vectorization, parallelization, accelerator specific languages, etc.) These optimizations heavily depend on the target machine and the whole code has to be adapted whenever it is ported to a new architecture.

As a result, scientists have to become experts in the art of computer optimizations in addition to their own domain of expertise. It is very difficult in practice to maintain a code targeting multiple distinct architectures. One fundamental cause for this situation is the tight interleaving of two distinct concerns imposed by most programming models. On the one hand, the algorithm comes from the expertise of the domain scientist and does not depend on the target architecture. On the other hand, optimization is the expertise of optimization specialists and has to be adapted for every new architecture.

Many approaches have been proposed to improve this situation in the form of libraries or languages. Approaches based on automated optimization processes typically isolate the algorithmic aspects very well but restrict their domain of applicability and/or the range of supported optimizations. Approaches based on optimization tools and libraries enable optimization specialists to express common optimizations very efficiently but leave others mixed with the algorithm.

In this paper, we propose the Independent Kernel Scheduling (InKS) programming model to separate algorithm from optimization choices in HPC simulation codes. We define the InKS pia language used to express the algorithm of an application independently of its optimization. Such a program can be optimized with the InKS O family of domain-specific languages (DSLs) for common optimizations or C++ for less common optimizations.

This paper makes the following contributions: 1) it defines the InKS programming model and its platform-independent algorithmic language InKS pia ; 2) it proposes an implementation of InKS with two optimization DSLs, InKS O/Loop and InKS O/XMP ; and 3) it evaluates the approach on the synthetic NAS parallel benchmarks [START_REF] Bailey | The nas parallel benchmarks[END_REF] and on the 6D Vlasov-Poisson solving with a semi-Lagrangian method.

The remaining of the paper is organized as follows. Section 2 presents and analyzes related work. Section 3 describes the InKS programming model and its implementation. Section 4 presents the 6D Vlasov-Poisson problem and its implementation using InKS while Section 5 evaluates the approach. Finally, Section 6 concludes the paper and identifies some perspectives.

Related works

Many approaches are used to implement optimized scientific simulations. A first widely used approach is based on imperative languages such as Fortran, C or C++. Libraries like MPI extend this to distributed memory with message passing. Abstractions very close to the actual execution machine make fine-tuning possible to achieve good performance on any specific architecture. It does however require encoding complex optimizations directly in the code. As there is no language support to separate the algorithm and architecture-specific optimizations, tedious efforts have to be applied [START_REF] Höhnerbach | The vectorization of the tersoff multibody potential: An exercise in performance portability[END_REF] to support performance portability. Algorithm and optimizations are instead often tightly bound together in codes.

A second approach is thus offered by tools (libraries, frameworks or language extensions) that encode classical optimizations. OpenMP [START_REF] Chandra | Parallel Programming in OpenMP[END_REF] or Kokkos [4] support common shared-memory parallelization patterns. E.g, Kokkos offers multidimensional arrays and iterators for which efficient memory mappings and iteration orders are selected independently. E.g, UPC [START_REF] El-Ghazawi | UPC: Distributed Shared Memory Programming[END_REF] or XMP [START_REF] Lee | Implementation and performance evaluation of xcalablemp: A parallel programming language for distributed memory systems[END_REF] support the partitioned global address space paradigm. In XMP, directives describe array distribution and communications between nodes. These tools offer gains of pro-ductivity when the optimization patterns they offer fit the requirements. The separation of optimizations from the main code base also eases porting between architectures. Even if expressed more compactly optimizations do however remain mixed with the algorithm and only cover part of the requirements.

A third approach pushes this further with tools that automate the optimization process. For example, PaRSEC [START_REF] Hoque | Dynamic task discovery in parsec: A data-flow task-based runtime[END_REF] or StarPU [START_REF] Augonnet | Starpu: A unified platform for task scheduling on heterogeneous multicore architectures[END_REF] support the many-tasks paradigm. In StarPU, the user expresses its code as a DAG of tasks with data dependencies that is automatically scheduled at runtime depending on the available resources. Another examples are SkeTo [START_REF] Tanno | Parallel skeletons for variable-length lists in sketo skeleton library[END_REF] or Lift [START_REF] Steuwer | Lift: A functional data-parallel ir for high-performance gpu code generation[END_REF] that offer algorithmic skeletons. Lift offers a limited set of parallel patterns whose combinations are automatically transformed by an optimizing compiler. Automating optimization improves productivity and clearly separate these optimizations which improves portability. The tools do however not cover the whole range of potential optimizations such as the choice of work granularity inside tasks. The algorithm remains largely interleaved with optimization choices even with this approach.

A last approach is based on DSLs that restrict optimizations, such as Pochoir [START_REF] Tang | The pochoir stencil compiler[END_REF] or PATUS [START_REF] Christen | PATUS: A Code Generation and Autotuning Framework for Parallel Iterative Stencil Computations on Modern Microarchitectures[END_REF], DSLs for stencil problems. In Pochoir, the user only specifies a stencil (computation kernel and access pattern), boundary conditions and a space-time domain while all optimizations are handled by a compiler. DSLs restrict the developer to the expression of the algorithm only, while optimizations are handled independently. This ensures a very good separation of these aspects. The narrower the target domain is, the more efficient domain and architecturespecific optimizations are possible. However, it makes it less likely for the tool to cover the needs of a whole application. Real-world applications can fall at the frontier between the domains of distinct DSLs or not be covered by a single one. Performance then depends on the choice of DSLs to use and the best choice depends on the target architecture leading to new portability issues.

To summarize, one can consider a continuum of approaches from very general approaches where the optimization process is manual to more and more domain specific where the optimization process can be automated. The more general approaches support a large range of optimizations and application domains but yield high implementation costs and low separation of concerns and portability. The more automated approaches reduce implementation costs and offer good separation of concerns and portability but restrain the range of supported domains and optimizations. Ideally, one would like to combine all these advantages: 1) the domain generality of imperative languages, 2) the ease of optimization offered by dedicated tools and 3) the separation of concerns and performance portability offered by DSLs. The following section describes the InKS programming model that aims to combine these approaches to offer such a solution.

The InKS programming model

This section describes the core of our contribution, the InKS programming model. This approach is based on the use of distinct languages to express algo-rithm and optimization choices; thus enforcing their separation. The algorithm of the simulation consists in the set of values computed, the formula used to produce each of them as well as the simulation inputs and outputs. We include in optimization choices all that is not the algorithm, such as the choice of a computing unit for each computation, their ordering, the choice of a location in memory for each value, etc. Multiple optimization choices can differ in performance but simulation results depend on the algorithm only.

Fig. 1: The InKS model

The InKS approach is summarized in Figure 1. The InKS pia language is used to express the algorithm with no concern for optimization choices. A compiler generates non-optimized, generic choices automatically from this specification for test purposes. The InKS O family of DSLs is used to define common optimizations efficiently while C++ is used to describe arbitrarily complex optimizations. Many versions of the optimization choices can be devised to optimize for multiple targets.

The remaining of the section describes InKS pia and proposes two preliminary InKS O DSLs. InKS O/XMP handles domain decomposition and inter-node communications while InKS O/Loop focuses on efficient loops.

The InKS pia language. In InKS pia [START_REF] Aumage | Inks, a programming model to decouple performance from semantics in simulation codes[END_REF] (illustrated in Listing 1), values are stored in infinite multidimensional arrays based on dynamic single assignment (DSA, each coordinate can only be written once). Memory placement of each coordinate is left unspecified. Computations are specified by kernel procedures that 1) take as parameters data arrays and integer coordinates; 2) specify the coordinate they might read and will write in each array; and 3) define either a C++ or InKS implementation. An InKS implementation defines kernels validity domains: coordinates where C++ kernels can generate values in arrays. Kernel execution order is left unspecified. The simulation entry point is a kernel marked public. This specifies a parameterized task graph (PTG) [START_REF] Cosnard | Compact dag representation and its dynamic scheduling[END_REF]. This representation covers a large range of problems but imposes a few limitations. Mostly, all the problem parameters must be known at launch time, it does for example not support adaptive mesh or time-steps. It is still possible to express these concerns outside InKS pia and call the InKS implementation multiple times with different parameters.

The InKS pia compiler [START_REF] Aumage | Inks, a programming model to decouple performance from semantics in simulation codes[END_REF] extracts C++ functions from InKS pia kernels and generates a function with correct but non-optimized loops and memory allocations to execute them. Arbitrarily complex versions of these can also be written

[n][m][l][k][j][i] with t6d[n][m][l][k][j][i] 7 #pragma xmp template_fix t6d[N][M][L][K][J][I] 8 f6d = (double(*)[M][L][K][J][I]) xmp_malloc(xmp_desc_of(f6d), N, M, L, K, J, I);
Listing 3: XMP code generated for the f6D decomposition of Listing 2 manually in plain C++ and rely on existing optimization tools. These functions can be called from any existing code whose language supports the C calling convention. However, that approach requires information present in InKS pia to be repeated. The InKS O DSLs thus interface optimization tools with InKS pia .

The InKS O/XMP optimization language. InKS O/XMP (illustrated in Listing 2) handles distributed memory domain decomposition by combining C and directives based on XMP and adapted for InKS. The compiler replaces these directives by C and XMP code (Listing 3). The inks decompose directive supports static or dynamic allocation of arrays from the algorithm. The domain size is extracted from InKS pia source and the user only has to specify its mapping onto memory. As in XMP, InKS O/XMP supports domain decomposition mapped onto an XMP topology. In addition, it supports dimension reordering and folding which consists in reusing the same memory address for subsequent indices in a given dimension. This feature is important to reuse memory due to the DSA nature of InKS pia arrays. The exchange directive supports halo exchanges. The required halo size is automatically extracted from the algorithm and the user only has to specify when to execute the communications and in what dimension. While XMP requires halo values to be stored contiguously with the domain, InKS O/XMP support a dynamic halo extension where halo values are stored in dedicated, dynamically allocated buffers to reduce memory footprint.

The InKS O/Loop optimization language. InKS O/Loop (illustrated in Listing 4) offers to specify manually loop nests for which the compiler generates plain

The 6D Vlasov/Poisson problem

The 6D Vlasov-Poisson equation, presented in (1), describes the movement of particles in a plasma and the resulting electric field. We study its resolution for a single species on a 6D Cartesian mesh with periodic boundary conditions. We solve the Poisson part using a fast Fourier transform (FFT) and rely on a Strang splitting (order 2 in time) for the Vlasov part. This leads to 6 1D advections: 3 in the space dimensions (x 1 , x 2 , x 3 ) and 3 in the velocity dimensions (v 1 , v 2 , v 3 ). Each 1D advection relies on a Lagrange interpolation of degree 4. In the space dimensions, we use a semi-Lagrangian approach where the stencil is not applied around the destination point but at the foot of characteristics, a coordinate known at runtime only. This is described in more details in [START_REF] Mehrenberger | Vlasov on gpu (vog project)******[END_REF].

                 ∂f (t, x, v) ∂t + v.∇xf (t, x, v) -E(t, x).∇vf (t, x, v) = 0 -∆φ(t, x) = 1 -ρ(t, x) E(t, x) = -∇φ(t, x) ρ(t, x) = f (t, x, v)dv (1)
The main unknown is f (f6D in the code), the distribution function of particles in 6D phase space. Due to the Strang splitting, a first half time-step of advections is required after f6D initialization but before the main time-loop. These advections need the electric field E as input. E is obtained through the FFT-based Poisson solver that in turn needs the charge density ρ as input. ρ is computed by a reduction of f6D. The main time-loop is composed of 3 steps: advections in space dimensions, computation of the charge density (reduction) and electric field (Poisson solver) and advections in velocity dimensions.

The 6D nature of f6D requires a lot of memory, but the regularity of the problem means it can be distributed in blocks with good load-balancing. Halos are required to hold values of neighbors for the advections. Connected halo zones would increase the number of points in all dimensions and consume too much memory. Split advections mean that halos are required in a single dimension at a time though. We therefore use dynamic halos composed of two buffers, one for each boundary of the advected dimension (denoted "right" and "left") as shown on Figure 2. Listing 2 corresponds to the InKS O/XMP implementation of this strategy.

Advections are the main computational cost of the problem, accounting for 90% of the sequential execution time. Six loops surround the stencil computation of each advection and in a naive version, the use of a modulo to handle periodicity and application along non-contiguous dimensions slow down the computation. To enable vectorization and improve cache use, we copy f6D elements into contiguous buffers along with the left and right halos. Advections are applied on these buffers before copying them back into f6D. Blocking further improves performance by copying 16 elements at a time. Listing 5 corresponds to the InKS O/Loop implementation of these optimizations.

Synthetic benchmarks

We have implemented 4 out of the 5 sequential NAS benchmark kernels (IS, FT, EP and MG) in InKS pia and optimized them with plain C++. The C++ version [START_REF] Griebler | Efficient nas benchmark kernels with c++ parallel programming[END_REF], is used as reference. We have also implemented a 3D heat equation resolution by finite differences (7-point stencil) with two distinct C++ optimization strategies from a single InKS pia source. Both strategies comes from the reference [START_REF] Kamil | Stencilprobe: A microbenchmark for stencil applications[END_REF]. One uses double buffering (Heat/Buf) and the other implements a cache oblivious strategy (Heat/Obl).

The NAS CG kernel relies on indirections not expressible in the PTG model of InKS pia . Its implementation would thus have to rely on a large C++ kernel whose optimization would be mixed with the algorithm. InKS pia can however be used to express all other NAS kernels as well as the 3D heat equation solver. Even if not as expressive as C or Fortran, InKS pia covers the needs of a wide range of simulation domains and offers abstractions close to the execution machine rather than from a specific simulation domain. Among others, it can express computations such as FFTs or stencils with input coordinates unknown at compile-time.

InKS separates the specification of algorithm and optimization in distinct files. Multiple optimization strategies can be implemented for a single algorithm, Finding the right metric to evaluate the easiness of writing a code is a difficult question. As illustrated in Listing 4 however, algorithm expression in InKS pia is close to the most naive C implementation where loops are replaced by InKS validity domains with no worry for optimization. The specification of optimization choices is close to their expression in C++. Table 1 compares the GNU complexity score of InKS optimizations to the reference code. InKS scores are slightly better because kernels extracted from the algorithm hide computations and thus, part of the complexity. In addition, the use of C++ to write optimizations let optimization specialists reuse their preexisting knowledge of this language. These considerations should not hide the fact that some information has to be specified both in the InKS pia and C++ files with this approach leading to more code overall.

Regarding performance, the InKS approach makes it possible to express optimizations that do not change the algorithm. Optimizations of the four NAS parallel benchmarks and 3D heat equation solver in InKS were trivial to implement and their performance match or improve upon the reference as presented in Table 1. Investigation have shown that Intel ICC 18 does not vectorize properly the reference versions of Heat/Obl and NAS/MG. The use of the Intel ivdep directive as done on the InKS versions leads to slightly better performance.

Vlasov-Poisson

We evaluate InKS O/XMP and InKS O/Loop on Vlasov-Poisson separately as they target different optimizations and are not usable together currently. A first experiment focuses on the sequential aspects with the intra-node optimization of the v 1 advection using either InKS O/Loop or plain C++. A second experiment focuses on the parallel aspects with the charge density computation, the Poisson solver and a halo exchange optimized either with C/XMP or with InKS O/XMP . The reference is the Fortran/MPI implementation from the Selalib library [START_REF] Kormann | Massively parallel semilagrangian solution of the 6d vlasov-poisson problem[END_REF]. For the v 1 advection, both the C++ and InKS O/Loop optimizations of the InKS code achieve performance similar to the reference as shown in Table 2. For the parallel aspects, the InKS O/XMP optimization offers performance similar to XMP as shown on Figure 3. The performance is comparable to MPI on the reduction operation but MPI is faster on the Poisson solver and the halo exchanges. At the moment, it seems that XMP does not optimize local copies which slows down the Poisson solver. Besides, XMP directives used for the halo exchanges are based on MPI RMA which make the comparison with MPI Send/Receive complex. Still, MPI is much harder to program: more than 350 lines of MPI and Fortran are required to handle domain decomposition, remapping for FFT and halo exchange in Selalib vs. 50 lines in XMP and 15 in InKS O/XMP .

The InKS O family of DSLs enables the developer to specify optimization choices only while algorithmic information is extracted from InKS pia code. This is illustrated by Listings 2, presenting the InKS O/XMP 6D domain decomposition, and 3, presenting the XMP version. Both are equivalent, but the InKS O/XMP expects only optimization choices parameters. Hence, one can test another memory layout, such as a different dimension ordering, by changing only a few parameters, while multiple directives must be modified in XMP. Similarly, with InKS O/Loop (Listing 4), developers can easily test different optimization choices that would be tedious in plain C++. Since InKS O/XMP and InKS O/Loop are respectively usable with C and C++, InKS does not restrict the expressible optimization choices: one can still implement optimizations not handled by InKS O in C/C++. Moreover, operations such as halo size computation or vectorization capabilities detection are automatized using the algorithm. In summary, the approach enables optimization specialists to focus on their specialty which make the development easier.

Conclusion and future works

In this paper, we have presented the InKS programming model to separate algorithm (InKS pia ) and optimization choices (InKS O & C++) and its implementation supporting two DSLs : InKS O/Loop for loop optimizations and InKS O/XMP for domain decomposition. We have evaluated InKS on synthetic benchmarks and on the Vlasov-Poisson solving. We have demonstrated its generality and its advantages in terms of separation of concerns to improve maintainability and portability while offering performance on par with existing approaches.

While this paper demonstrates the interest of the InKS approach, it still requires some work to further develop it. We plan to apply the InKS model on a range of different problems. We will improve the optimization DSLs; base InKS O/Loop on existing loop optimization tools and ensure good interactions with InKS O/XMP . We also want to target different architectures to demonstrate the portability gains of the InKS approach.

  (a) Buffers used as halo regions in the first dimension (b) Buffers used as halo regions in the second dimension

Fig. 2 :

 2 Fig. 2: Dynamic halo exchange representation on a 2D domain

  Dynamic halo exchange on the 4th dimension, halo sizes are computed automatically 6 #pragma inks exchange periodic f6d(4, advection4) to R and L 7 /*R and L are now allocated buffers and contain the halo values: */ foo(R, L);

	1 kernel stencil3(x, t): ( double A(2) {in: (x-1:x+1, t) | out: (x, t+1)} )
	2	#CODE(C) A(x, t+1) = 0.5*A(x, t-1) + 0.25*(A(x-1, t-1)+A(x+1, t-1)); #END
	3 public kernel inks_stencil(DIM_X, N_ITER):
	4	( double A_decl(2) {in: (0:DIM_X, 0) | out: (1:DIM_X-1, N_ITER-1)} )
	5	#CODE(INKS) stencil3(1:DIM_X-1, 1:N_ITER-1):(A_decl) #END
		Listing 1: 1D stencil computation on a 2D domain in InKS pia
	1 #pragma xmp nodes p6d[pV3][pV2][pV1][pZ][pY][pX]
	2 #pragma inks decompose dynamic % f6d // dynamic allocation of f6d algorithmic array
	3	(8:, 7:, 6, 5, 4, 3, 2, 1) // dimension reordering, dim 7 & 8 are folded
	4	with t6d onto p6d // block decomposition mapped on the XMP topology
	5 // Listing 2: 6D decomposition and dynamic halo exchange in InKS O/XMP
	1 double (f6d*)[][][][][]; // need to declare f6d global, valid in xmp
	2 #pragma xmp nodes p6d[pV3][pV2][pV1][pZ][pY][pX] // xmp 6d cartesian node topology
	3 #pragma xmp template t6d[:][:][:][:][:][:] // xmp 6d logical array
	4 #pragma xmp distribute t6d [block][block][block][block][block][block] onto p6d
	5 // map element of f6d to element of t6d
	6 #pragma xmp align f6d

  1 /*** advec.iks InKSpia algorithmic file ***/ 2 kernel advection3 (i, j, k, l, m, n, t, K, step) : ( It takes as parameters the size of block for the loops starting from the innermost one. If there are less block sizes than loops, the remaining loops are not blocked. The Buffer keyword supports copying data in a local buffer before computation and back after to ensure data continuity and improve vectorization. The compiler uses data dependencies from the algorithm to check the validity of the loops order and generate vectorization directives where possible.

	3	double f6d { in: (i, j, 0:K, l, m, n, t, step-1) | out: (i, j, k, l, m, n, t, step) },
	4	int disp {in: (n)}, double coef {in: (n, 0:4)} )
	5 #CODE (C) /* ... */ #END
	6
	7 public kernel main_code(t, I, J, K, L, M, N, Niter): ( double coef(2), int disp(1)
	8	double f6d(6) { in: (i, j, 0:K, l, m, n, t, step-1) | out: (i, j, k, l, m, n, t,
		step) } )
	9 #CODE(INKS)
	10	/* ... */
	11	advection3 (0:I, 0:J, 0:K, 0:L, 0:M, 0:N, 1:Niter, K, 2) : (f6d, disp, coef)
	12	/* ... */
	13 #END
	14
	15 /*** advec.iloop InKSo/Loop optimization choices file ***/
	16 loop advection3_loops(t, I, J, K, L, M, N, Niter) : advection3 { // set "t" value
	17	// "Set" not specified -> loop bounds are computed, with a fixed "t"
	18	Order: n, m, l, j, i, k; // order of the loop
	19	Block: 16; // blocking on the inner dimension k
	20	Buffer: f6d(3); } // copy the third dimension of f6d to a 1d buffer

Listing 4: A loop nest in InKS pia optimized in InKS O/Loop 1 /* for all N, M, L, J, do */ for (int i=0; i<I; i+=blockSize){ 2 for(int ci=-halo_size; ci<0; ++ci) /*Copy to buffer*/ 3 for(int ii=0; ii<blockSize; ++ii) buff.buff_in[...] = left[...]; 4 for(int ci=0; ci<K; ++ci) 5 for(int ii=0; ii<blockSize; ++ii) buff.buff_in[...] = f6d[...]; 6 for(int ci=sizeK; ci<sizeK+halo_size; ++ci) 7 for(int ii=0; ii<blockSize; ++ii) buff.buff_in[...] = right[...]; 8 for(int idb=0; idb<size_block; ++idb) /*Computation*/ 9 for(int k=0; k<K; ++k) advection3(buff, idb+i, j, k, l, m, n); 10 for(int ci=0; ci<K; ++ci) /*Copy to f6d*/ 11 for(int ii=0; ii<block_size; ++ii) f6d[...] = buff.buff_out[...];} Listing 5: C++ code generated for the loop nest of Listing 4 C++ loops (Listing 5). Plain C++ is usable with InKS O/Loop . The loop keyword introduces a nest optimization with a name, the list of parameters from the algorithm on which the loop bounds depend and a reference to the optimized kernel. Loop bounds can be automatically extracted from InKS pia , but the Set keyword makes it possible to restrict these bounds. The Order keyword specifies the iteration order on the dimensions and the Block keyword enables the user to implement blocking.

Table 1 :

 1 Execution time of the C++ and InKS implementations of the sequential NAS benchmark, class B -Time/iteration of the 3D heat equation (7point stencil), size (1024 3 ) -Median and standard deviation of 10 executions -GNU Complexity score of the implementation. as shown for the 3D heat equation where each relies on a specific memory layout and scheduling. It thus offers a clear separation of algorithm and optimization.

		Execution time (second)	Complexity
	Benchmark	Reference	InKS	Rel. dev.	Ref. InKS
	NAS/FT 62.43 (±0.57) 62.86 (±0.71)	0.68%	6	5
	NAS/IS	3.39 (±0.00)	3.44 (±0.00)	1.47%	55	52
	NAS/MG	5.10 (±0.02)	4.73 (±0.06)	-7.25%	20	12
	NAS/EP 76.43 (±0.21) 76.47 (±0.22)	0.05%	19	19
	Heat/Buf	3.05 (±0.01)	2.97 (±0.06)	-2.58%	5	3
	Heat/Obl	2.43 (±0.01)	2.05 (±0.02)	-15.59%	22	13

  Version time/advection GFLOP/s % peak core perf.

	Selalib (Fortran)	4.81 s	1.12	5.36%
	InKS (C++)	3.76 s	1.43	6.87%
	InKS (InKS O/Loop )	3.61 s	1.49	7.16%

Table 2 :

 2 Comparison between Fortran, C++ and InKS O/Loop version of the v 1 advection on a 32 6 grid (double precision) on a single E5-2670 core with vectorization. Median of 12 executions Weak and strong scaling for 3 parts of the Vlasov-Poisson solver up to 64 nodes (1 process/node) on a 32 6 grid divided among processes (strong scaling) or 166 grid per process (weak scaling). Median of 10 executions.Distinct files for both concerns in InKS makes possible to write a unique InKS pia algorithm and multiple versions of optimization choices. Four optimization choices are implemented based on one InKS pia source : (1) InKS O/XMP , (2) C with XMP, (3) InKS O/Loop and (4) plain C++. This proves, to some extent, that the separation of concerns is respected. As of now, InKS O/XMP and InKS O/Loop are not usable together since InKS O/Loop relies on C++ that XMP does not support. We plan to address this limitation in the future.

	1.3							1					1
	1.2
								0.3
	1												0.1
								0.1
				XMP weak
	0.8		InKSO/XMP weak MPI weak								0.01
				XMP strong				0.03
			InKSO/XMP strong MPI strong										Local copies		Remote copies
	0.6	1	2	4	8	16	32	64	8	16	32	64	0.001	1	2	4	8	16	32	64
		(a) f6D reduction			(b) Poisson solver		(c) Dynamic halo exchanges
	Fig. 3: