A Menasria

P Brenner
email: pierre.brenner@ariane.group

P Cinnella
email: paola.cinnella@ensam.eu

G Pont

Ph D Student

Toward an improved wall treatment for multiple-correction k-exact schemes

Improved wall boundary treatments are investigated for a family of high-order Godunovtype finite volume schemes based on k-exact polynomial reconstructions in each cell of the primitive variables, via a successive corrections procedure. We focus more particularly on the 1-exact and 2-exact schemes which offer a good trade-off between accuracy and computational efficiency. In both cases, the reconstruction stencil needs to be extended to the boundaries. Additionally, information about wall curvature has to be taken into account, which is done by using a surface model based on bicubic Bézier patches for the walls. The performance of the proposed models is presented for two compressible cases, namely the inviscid flow past a Gaussian bump and the viscous axisymmetric Couette flow.

I. Introduction

H -order finite volume methods for Computational Fluid Dynamics have received considerable interest in recent years, because of their potential for providing high-fidelity solutions of challenging configurations involving, e.g., flow unsteadiness, boundary layer separation or steep gradients.

A widely studied family of high-order schemes is represented by the so-called k-exact schemes, initially introduced by Barth and Frederickson [START_REF] Barth | Higher Order Solution of the Euler Equations on Unstructured Grid Using Quadratic Reconstruction[END_REF] [START_REF] Barth | Recent Developments in High Order k-exact Reconstruction on Unstructured Meshes[END_REF] and Vandersbilck and Deconinck [START_REF] Vankeirsbilck | Higher order upwind finite volume schemes with ENO-properties for general unstructured meshes[END_REF], which make use of high-order piecewise polynomial reconstructions of the solution and of its derivatives on mesh cells, along with high-order formulae for the computation of flux integrals (see [START_REF] Ollivier-Gooch | A High-Order Accurate Unstructured Mesh Finite-Volume Scheme for the Advection-Diffusion Equation[END_REF] [START_REF] Caraeni | Unstructured-Grid Third-Order Finite Volume Discretization Using a Multistep Quadratic Data-Reconstruction Method[END_REF] [6] [START_REF] Pont | Multiple-correction hybrid k-exact schemes for High-Order compressible RANS-LES simulations on fully unstructured grids[END_REF] for instance).

In past work [START_REF] Haider | Efficient Implementation of High Order Reconstruction in Finite Volume Methods[END_REF] [8] [START_REF] Pont | Multiple-correction hybrid k-exact schemes for High-Order compressible RANS-LES simulations on fully unstructured grids[END_REF], a family of k-exact schemes ensuring high-accuracy on fully unstructured grids was constructed thanks to an efficient successive corrections procedure. The schemes were applied to the computation of several steady and unsteady flow configurations [START_REF] Pont | Multiple-correction hybrid k-exact schemes for High-Order compressible RANS-LES simulations on fully unstructured grids[END_REF], showing their capability of ensuring high-order accuracy on various grid topologies.

In this context, making sure that these schemes model the flow close to the wall as accurately as possible is of critical importance: indeed, near-wall regions are characterized by several complex phenomena.

Additionally, these exhibit low local values of the Reynolds and Mach numbers. The accuracy of low-Mach flow simulations is directly affected by the ability of compressible solvers to reproduce flows in the incompressible limit [START_REF] Qu | Investigation into the influences of the low-speed flows' accuracy on RANS simulations[END_REF].

Another important issue for high-order simulations of complex geometries is the treatment of curved boundaries. Approximating a smooth surface (i.e. without ridge or fault) with flat faces introduces a spatial error of order O(h 2 ), where h is the mesh size [START_REF] Giesselmann | Geometric Error of Finite Volume Schemes for Conservation Laws on Evolving Surfaces[END_REF]. The real shape of the boundary has to be approximated to at least the order of accuracy of the discretization scheme [START_REF] Ollivier-Gooch | Obtaining and Verifying High-Order Unstructured Finite Volume Solutions to the Euler Equations[END_REF].

Finally, an accurate discretization of the boundary conditions is a crucial ingredient for the overall accuracy, although it is not straightforward to achieve.

In the k-exact schemes proposed in [START_REF] Pont | Multiple-correction hybrid k-exact schemes for High-Order compressible RANS-LES simulations on fully unstructured grids[END_REF] [START_REF] Haider | Efficient Implementation of High Order Reconstruction in Finite Volume Methods[END_REF] [START_REF] Brenner | Unsteady flows about bodies in relative motion[END_REF], a simplified reconstruction stencil has been used close to solid walls, leading to a local loss of accuracy that can propagate to the entire domain, thus affecting the global order of accuracy of the simulation. Moroever, wall curvature is not taken into account (only polygonal faces have been considered up to now).

In the aim of progressing toward a more accurate representation of flow configurations with curved boundaries, two major improvements of the wall boundary treatment are investigated in the present paper. On the one hand, a high-order representation of wall faces is introduced by means of bicubic Bézier patches, which are then employed to calculate the geometrical moments, that are notably necessary in the fluxes integration in our methodology. On the other hand, solutions to enrich the reconstruction stencil in the vicinity of solid walls in order to improve the local order of accuracy are proposed and analyzed.

All of the present developments are implemented within the industrial code FLUSEPA * , developed by the European aerospace company Ariane Group to model all the phases of flight of a space launcher such as takeoff, stage separation, and reentry [START_REF] Pont | CFD modeling of the Ariane 6 space launch from ground to space[END_REF]. The solver can simulate tridimensional, unsteady, compressible, viscous and reactive flows loaded with particles over bodies in relative motion. A CHIMERA-like strategy based on 3D intersections coupled with an Adaptive Mesh Refinement module [START_REF] Limare | Toward a high-order mesh adaptation strategy for unsteady flows on overlapping grids in a finite-volume framework[END_REF] is used to deal with complex geometries.

The present paper is organized as follow. In Section II, the off-wall spatial discretization implemented in FLUSEPA is recalled. The surface model chosen to account for boundary curvature and its implementation within the solver are described in Section III. Strategies for near-wall discretization are presented in Section IV. Numerical validations for selected inviscid and viscous flows are shown and discussed in Section V. Finally, conclusions and perspectives for future work are made in Section VI.

II. Numerical scheme

In this Section, the spatial discretization scheme of FLUSEPA [START_REF] Pont | Multiple-correction hybrid k-exact schemes for High-Order compressible RANS-LES simulations on fully unstructured grids[END_REF] [6] is described. Source terms are not considered in the present paper.

Let q(x, t) be a variable verifying a set of conservation equations of the form

∂q ∂t + ì ∇ • F = 0 (1)
F is the flux function of q, composed of a convective part that depends only on the solution variable q, and of a diffusive part that depends on the solutions variable q and on its gradient ì ∇q.

F = F c (q) + F d (q, ì ∇q) (2)
The modification of the calculation of diffusive fluxes has not been considered, and will not be detailed in the present article (F d is thus assumed null in the following). Note however that improving the accuracy of the solution reconstruction close to walls is also expected to have a beneficial effect on the viscous flux approximation.

In the cell-centered finite volume framework, the domain is divided into N non-overlapping cells

(Ω J ) J ∈[|1, N |] of volume |Ω J | and cell center x j = 1 |Ω J | ∬ Ω J
xdV. The computational variables are the averaged values

q J = 1 |Ω J | ∭ Ω J qdV. (3)
By integrating Eq. 1 on Ω J and applying Gauss' theorem, we get

|Ω J | dq J dt + P K=1 ∬ A J K F c (q) • ndS (4)
with A J K the interface between cell J and its neighbor K, P the total number of neighbors of J and n the local outward-pointing normal.

The solver consists in three steps, following the pattern of k-exact schemes [START_REF] Barth | Higher Order Solution of the Euler Equations on Unstructured Grid Using Quadratic Reconstruction[END_REF] [2]:

• Reconstruction of the solution in each cell;

• Flux integration along the faces;

• Advance in time. The time stepping scheme is a second-order Heun scheme with consistent local time stepping (see [START_REF] Brenner | Unsteady flows about bodies in relative motion[END_REF] for details). In the following, time variable t will be omitted to simplify notations.

A. Solution reconstruction

To reconstruct the solution in a given cell Ω J to the k th order, a Taylor expansion of order k at cell centroid x j is performed. For a third-order reconstruction q(x) = q j + D (1)

q j • (x -x j ) + 1 2 D (2) q j • (x -x j ) ⊗ (x -x j ) + O(h 3 ) (5)
D (1) q and D (2) q are respectively the gradient vector and the second derivatives tensor of q. Index j refers to pointwise value at x j . • is the contracted product and ⊗ is the tensor product.

A classical method to determine a reconstruction is to solve the system of Eq. 5 such that the derivatives verify conservativity relations in a least-squares sense [START_REF] Jalali | Accuracy Analysis of Unstructured Finite Volume Discretization Schemes for Diffusive Fluxes[END_REF] [START_REF] Wang | Compact high order finite volume method on unstructured grids I: Basic formulations and one-dimensional schemes[END_REF]. Such method has the drawback to be sensitive to mesh irregularities [START_REF] Mavriplis | Revisiting the Least-Squares Procedure for Gradient Reconstruction on Unstructured Meshes[END_REF].

Let us call s(J) the direct neighborhood of cell Ω J , comprising the cells Ω K that share a common edge with Ω J , and Ω J itself. The reconstruction procedure used in this work is described below.

Linear reconstruction of the solution

For a second-order accurate reconstruction, we need to calculate q j at second order and D (1) q j at first order. For this purpose we construct a 1-exact reconstruction operator for the gradient D (1) 1 [q] J . In this notation: • the italic font means that the function is a linear approximation operator;

• the index letter refers to the reconstruction stencil. Here it is J, hence the operator is defined on s(J); • q = (q K ∈s(J) ) are the averaged variables of the conservative variables on the stencil s(J). The left and right brackets in [q] means that the gradient is constructed as a function of these variables; • the exponent (1) refers to the derivation order. Here it is first order, which means the operator approximates the first derivatives (solution gradient); • the index number refers to the accuracy order. Here we are looking for a 1-exact [START_REF] Pont | Multiple-correction hybrid k-exact schemes for High-Order compressible RANS-LES simulations on fully unstructured grids[END_REF] operator on s(J), i.e. such that ∀x ∈ s(J), D (1) 1

[q] J (x) = D (1) q(x) + O(h). (6)
In the solver, the reconstruction is applied to the primitive variables since this represents a more robust choice than the conservative ones. For each primitive variable q, an average q is defined, which can be expressed as a function of the averaged conservative variables. Reconstruction of q from the conservative variables in the k-exact framework is described in [START_REF] Pont | Multiple-correction hybrid k-exact schemes for High-Order compressible RANS-LES simulations on fully unstructured grids[END_REF]. For a second-order accurate reconstruction, we want

q = q + O(h 2 ). (7)
To compute the gradient, a technique based on a modification of the Green-Gauss method is used [7] [8]. Gauss' theorem is applied to the exact gradient to bring back the problem to the calculation of surface integrals, which are then approximated by using a linear interpolation between the left and right averaged cells values.

1 |Ω J | ∭ Ω J D (1) q(x)dV = 1 |Ω J | P K=1 ∬ A J K qndS = D (1) q(x j ) + O(h) (8)
∬ A J K qndV = β K q K + (1 -β K )q J A J K , with β K = ||x j -x Γ || ||x j -x Γ || + ||x k -x Γ || ∈ [0, 1] (9)
with

A J K = O(h 2 )n J K the surface vector of interface A J K , |Ω J | = O(h 3
) the volume of cell J, and x Γ a chosen point related to the interface called integration point † . Note that this formula is not used for wall faces: the original solver uses β K = 0 for these faces such that the contribution of the wall is totally excluded from the recontruction stencil. This will be modified in Section IV. Additional conditions for the 1-exact gradient operator are given by the k-exact approach

K ∈s(J) W J K = 0 (0-exactness) K ∈s(J) (α • (x k -x j ))W J K = α (1-exactness).
(10) † The choice of such a point will be discussed in Section III.

3 where (W J K ) K ∈s(J) are the weight vectors such that D (1) 1 [q] J = K ∈s(J) W J K q K . By combining equations 8 and 9, we obtain a first discrete approximation for the gradient given by D (1) 0

[q] J = 1 |Ω J | K ∈s(J) (β K q K + (1 -β K )q J )A J K (11)
However, if approximation 9 is of order O(h 2 )A J K when x j , x Γ and x k are aligned, it is only of order O(h)A J K otherwise. Thus the gradient operator 11 is no more consistent on general meshes. To recover consistency on unstructured grids, a corrective operator is built, with its associated matrix M 1 (J) called simple correction matrix [START_REF] Brenner | Unsteady flows about bodies in relative motion[END_REF].

The matrix M 1 (J) of the linear operator 11 is established by applying the conditions 10 to the canonical basis of the polynomials of degree 1 xx j , yy j , zz j . We get

M 1 (J)D (1) 1 [q] J = D (1) 0 [q] J (12)
and finally obtain what we call a Quasi-Green 1-exact gradient operator

D (1) 1 [q] J = M -1 1 (J) 1 |Ω J | K ∈s(J) (β K q K + (1 -β K )q J )A J K , with M 1 (J) = K ∈s(J) β K A J K ⊗ (x k -x j ) (13)
and M -1 1 (J) the inverse of matrix M 1 (J). Note that for cartesian meshes, M 1 (J) is the identity matrix. This procedure leads a second-order accurate scheme.

Parabolic reconstruction of the solution

For a third-order accurate reconstruction we need to calculate q j at third order, D (1) q j at second order and D (2) q j at first order.

First, we must compute the primitive averaged values q with third-order accuracy. A correction term ∆q can be established using what has been done in part 1: details of this procedure can be found in [START_REF] Pont | Multiple-correction hybrid k-exact schemes for High-Order compressible RANS-LES simulations on fully unstructured grids[END_REF].

Afterwards, successive corrections of the derivatives are carried out. The 1-exact gradient operator described above is applied twice, leading to an inconsistent approximation of the second derivatives in x j D (2) 0 [q] J = D (1) 1

D (1) 1 [q] J K (14)
The consistency of the second-order accurate gradient operator is restored on general meshes using the linearity properties of the numerical gradient operator and 2-exactness constraints (see [START_REF] Haider | Applications of Efficient Parallel k-exact Finite Volume Reconstruction on Unstructured Grids[END_REF]). This gives the third-order tensor M 2 (J), which depends on geometrical quantities of two levels of neighbors of cell J (neighbors of J and their neighbors).

D (2) 1 [q] J = M 2 (J) -1 D (2) 0 [q] J (15)
We then correct the 1-exact approximation gradient to increase its accuracy to the second order, by calculating the truncation error between the exact gradient and its 1-exact reconstruction operator. For each component m ∈ {1, 3} of the gradient, a matrix H 1 (J) m is calculated, which depends only on geometrical quantities of the cells of the direct neighborhood s(J):

(D (1) q j ) m = (D (1) 2 [q] J ) m (D (1) 1 [q] J ) m -D (2) q j • H 1 (J) m +O(h 2 ) (16)
In 16, D (2) q j is established using the approximation 15. This procedure leads to a third-order scheme.

B. Flux integration

Flux integration along face A J K is achieved by using a one-point integration formula, based on a Taylor expansion to the third order:

∬ A J K F c • ndS = F c |Γ • S (0) A J K + D (1) F c |Γ • S (1) A J K + 1 2 D (2) F c |Γ • S (2) A J K + O(h 3 )A J K (17)
where x Γ is the integration point and where S (m) A J K are the surface moments of order m of the interface

A J K S (m) A J K = ∬ A J K (x -x Γ ) ⊗m ndS, m ∈ {0, 1, 2}; (18)
Surface moments are fundamental for the solver. They are tensors of order m which are sized in O(h m ). As will be seen in Section III, particular attention must be paid to the calculation of surface moments on the one hand, and on the choice of the integration point on the other hand when taking into account curvature. Eq. 17 is different from the classical choice of Gauss points to integrate fluxes. For a third-order reconstruction process, we need F c |Γ at order 3, D (1) F c |Γ at order 2 and D (2) F c |Γ at order 1. These approximations are obtained by solving a Riemann problem at the interface A J K . Left and right values of the variable are given by the solution reconstruction and derivative approximates established previously, following a MUSCL-type approach [START_REF] Van Leer | Towards the ultimate conservative difference scheme. V. A second-order Sequel to Godunov method[END_REF].

Note that the overall k-exact scheme presented in this Section can formally be derived to any order of accuracy. For a scheme of order m, one would have to perform the fluxes' expansion to the m th order, and to derive the corresponding adapted approximations of the derivatives using the successive corrections procedure generalized from what has been done for the 2-exact reconstruction scheme.

• x k • x j • x Γ • x m Ω J Ω K A J K • x j Ω J • x k Ω K • x l Ω L • x Γ • x m

III. Wall curvature

Modelling the wall curvature is a prerequisite to high-order boundary treatment. In the baseline implementation of FLUSEPA, boundary faces are represented by flat faces ‡ . This introduces a second-order spatial error compared to the actual surface [START_REF] Giesselmann | Geometric Error of Finite Volume Schemes for Conservation Laws on Evolving Surfaces[END_REF]. However, the approximation of the physical surface has to be at least of the order of the numerical scheme to ensure the nominal order of accuracy close to the wall [START_REF] Ollivier-Gooch | Obtaining and Verifying High-Order Unstructured Finite Volume Solutions to the Euler Equations[END_REF]. In the case of parametric surfaces using a polynomial basis, the polynomials have to be at least of degree two.

For this purpose, we use the geometric module BSHAPE developed by INRIA [START_REF] Borouchaki | BSHAPE, interpolant quasi G1 d'une surface discrète : descriptif technique[END_REF], which generates from a surface mesh a composite surface made of bicubic Bézier patches.

A bicubic Bézier patch is a parametric polynomial surface of degree three originated by De Casteljau and Bézier [START_REF] Bézier | The mathematical basis of the UNISURF CAD system[END_REF] around the same time. Given a set of 16 control points (P i j ) i ∈ {0,1,2,3}, j ∈ {0,1,2,3} , the equation of the Bézier patch is given by [START_REF] Farin | Curves and Surfaces for Computer-Aided Geometric Design -A practical guide[END_REF]: ‡ More precisely, triangular faces are used. Quadrangular faces are treated by being split into two triangles.

∀(u, v) ∈ [0, 1] × [0, 1], f (u, v) = 3 i=0 3 j=0 B 3 i (u)B 3 j (v)P i j ( 19
)
where ∀u

∈ [0, 1], ∀i ∈ {0, 1, 2, 3}, B 3 i (u) = 3! i!(3-i)!
(1u) 3-i u i are the Bernstein polynomials of degree three. BSHAPE module is based on a point normal interpolation meshing strategy [START_REF] Walton | A triangular G1 patch from boundary curves[END_REF]: the patches are generated with the points coordinates and associated normals of the mesh. The modelization ensures G1 continuity, i.e. tangent plane continuity, between the patches. The module has been incorporated in FLUSEPA. The geometry routine selects the wall faces and the side faces sharing an edge with the wall faces, and generates a set of control points corresponding to a patch for each wall face (see Figure 2 ). The composite surface is not directly seen by the solver, that is, a new mesh is not created, but the information about the patches, namely the control points, is used to integrate the boundary geometric quantities required for the scheme introduced in Section II. The surface moments of the wall faces are derived directly by integrating analytically the constitutive equation of the patch. Then volume moments and centroids of the boundary cells are deduced by the application of Gauss' theorem.

This strategy was tested on simple configurations for which the analytical equation of the surface is known (parabolic channel, sphere). Geometrical quantities (such as the surface and volume moments required by the k-exact scheme) were calculated using 19 and we observed that the error decreased from 1 up to 2 orders of magnitude when the module was used.

In addition to computing more accurately the geometrical quantities, flux integration on curved wall faces needs further modifications. In the original solver, the integration point x Γ is constructed to minimize the first-order error in the flux integration: x Γ is chosen so that the associated surface moments of order 1 are orthogonal to A J K . Consequently, the first-order error in the normal face direction n J K =

A J K | | A J K | | cancels out by construction. x Γ = 1 | | A J K | | ∬ A J K xndS • n J K 1 | | A J K | | ∬ A J K yndS • n J K 1 | | A J K | | ∬ A J K zndS • n J K (20)
In case of a flat face, n is constant equal to n J K thus x Γ is the center of gravity of the face and lies on the face. However, when the face is curved, n is not constant: x Γ is not equal a priori to the center of gravity of the face. More importantly, it does not generally lie on the face.

In case of a curved wall, flux and its derivatives needed in 17 are evaluated at a particular point of the Bézier patch. In a first approach, the middle point of the patch x m is chosen (point of the patch with parameters u = 0.5 and v = 0.5). This is not a problem for the error generated during the flux integration, as long as the flux and its derivatives are approximated to the required order of accuracy. Concerning the resolution of the Riemann problem at the wall, the local normal at point x m is used instead of the integrated normal n J K .

The radius of curvature is also computed for each patch and will be used in Section IV in the boundary conditions of Euler cases. In two-dimensional space, any twice differentiable curve can be locally approximated by a circular arc whose radius is called the radius of curvature. In three-dimensional space, several types of curvatures can be defined. In the present work, we consider the normal curvature associated to a given direction. The tangential velocity (generally non-zero for Euler cases) is calculated at the middle point of the patch, and the radius of curvature R associated to this vector and to the normal at this point is estimated (see Chapter 19 of [START_REF] Farin | Curves and Surfaces for Computer-Aided Geometric Design -A practical guide[END_REF]).

IV. Near-wall discretization

In this Section, a near-wall discretization consistent with the second and third-order accurate schemes described in Section II is derived. FLUSEPA originally does not take into account the wall in the polynomial reconstruction of the solution variable: for wall faces, the parameter β introduced in Section II is set to 0, which corresponds to a complete decentering to the inside of the domain.

The original stencils used respectively for the second and third-order schemes corresponding to a boundary cell Ω J are shown in Figure 1.

In the following, we want to account for the wall contribution in the calculation of the derivatives approximation operators. Denote by Ω K 0 the wall boundary face of a given cell Ω J , A J K 0 . This would correspond to a volumeless cell and is actually just a notation that will permit to extend formally the previous formulae. Let us denote s * (J) the extended neighborhood of Ω J as

s * (J) = s(J) ∪ Ω K 0 ( 21
)
with the convention that Ω K 0 = ∅ when Ω J is not a boundary cell, so that the extended stencil s * (J) is equal to s(J) when Ω J is inside the domain, and symmetrizes the stencil when Ω J is a boundary cell. The extension of the numerical scheme is explained hereafter for the 1-exact and for the 2-exact reconstructions of the solution variables. Treatment of the wall boundary conditions for the Euler and Navier-Stokes equations is also discussed.

A. Modification of the linear reconstruction for boundary cells

In this part, a simple approach for achieving a truly second-order accurate reconstruction of the solution in boundary cells is presented.

To build a second-order polynomial representation of the variable q that uses the wall information, the Quasi-Green gradient is modified consistently with the methodology introduced in the first part of Section II. Assume that a Dirichlet condition function q D (x ∈ ∂Ω) is prescribed at the wall. We need to find a value at the wall q k 0 such that, similarly to Eq. 8

∬ A J K 0 q(x)n(x)dS ≈ (β K 0 q k 0 + (1 -β K 0 )q J )A J K 0 (22)
A Taylor expansion of q D at any point of the wall face x W ∈ A J K 0 in the integral of 22 gives:

∬ A J K 0 q(x)n(x)dS = ∬ A J K 0 (q D (x W ) + O(h)) ndS, with x W ∈ A J K 0 . (23)
thus:

∬ A J K 0 qndS = q D,W A J K 0 + O(h)A J K 0 . ( 24
)
and finally:

∬ A J K 0 qndS = (β K 0 q D,W + (1 -β K 0 )q J )A J K 0 + O(h)A J K 0 with β K 0 = 1 (25)
This wall contribution approximation is only first-order accurate on general grids, however this is not a problem since by construction, the consistency of the overall gradient operator will be restored thanks to the simple correction matrix M 1 (J). More precisely, M 1 (J) must be modified accordingly to 25 in order to satisfy the 1-exactness consistency condition of Eq. 10. By following the same reasoning as in Eq. 13, we obtain the simple correction matrix with a reconstruction stencil extended to the wall, noted M * 1 (J):

M * 1 (J) = 1 |Ω J | K ∈s * (J) β K A J K ⊗ (x k -x j ), with β K 0 = 1 and x K 0 = x W (26)
Finally, a new 1-exact gradient operator D * (1)

1
based on a more symmetrical stencil is proposed

D * (1) 1 [q] J = M * -1 1 (J) 1 |Ω J | K ∈s * (J) (β K q K + (1 -β K )q J )A J K . ( 27
)
with the notation q K 0 = q D,W . Let's now consider the case for which a Neumann condition ∂q ∂n (x ∈ ∂Ω) = D 1 q(x) • n(x) = q N (x) is prescribed at the boundary. To avoid multiple definitions of the simple correction matrix M 1 , a second-order accurate equivalent fictive value at the wall q f ict W is calculated in order to fall back into the Dirichlet case. However, one has to be careful to evaluate this fictive value at a unique point x W , independently on the type of boundary condition. For Dirichlet-type conditions, we have seen before that any point could be chosen for x W , as long as it lies on the wall face. For Neumann-type conditions, we perform a Taylor expansion in the normal direction of the chosen wall point to obtain this fictive value. Hence, our choice is to choose the wall point as the point x W = x N whose corresponding normal is colinear to the vector x Nx j . In these conditions, q f ict W is given by

q f ict W = q J - ∂q ∂n |W ((x j -x W ) • n W ) + O(h 2 ) (28)

B. Modification of the parabolic reconstruction for boundary cells

As can be seen in Figure 1 (b),the Quasi-Green operator D * (1)

1

[q] J established in part A can not be applied twice since gradients in the neighborhood of Ω J are missing (assuming a Dirichlet-type boundary condition).

We choose to apply the enriched gradient operator and its associated correction matrix whenever possible, and to keep the baseline approximation otherwise. With this choice, the initial (inconsistent) approximation of the second derivatives is obtained by applying the modified 1-exact operator D * (1) 1 onto the original operator D (1) 1 :

D * (2) 0 [q] J = D (1) 1 D * (1) 1 [q] J K (29)
The correction matrix M 2 (J) is modified in a similar way and a matrix M * 2 (J) is derived. Then, the geometrical matrices (H 1 (J)) i ∈ {1,3} appearing in Eq. 16 which give the truncation error between the exact gradient and its 1-exact representation are changed using what has been done in paragraph A: (H 1 (J)) i ∈ {1,3} originally depends on the parameters β K and on the second-order volume moments of the cells of s(J). We thus take in the wall contribution β K 0 = 1 and the volume moments for cell Ω K 0 null to calculate (H * 1 (J)) i ∈ {1,3} . The delicate point of this method is the accuracy of the boundary value q(x W ). Indeed, thanks to the correction matrices, the derivatives will reach their expected order of accuracy in boundary cells (second-order for the gradient, first order for the second derivatives) if and only if q W is at least third-order accurate. Boundary conditions for q W are studied in the next paragraph.

In part A, we chose to directly make the surface integral of Eq. 25 (which corresponds to the wall contribution of the operator 12) consistent, which led us to β W = 1, and then we used the extended 1-exact gradient in the successive corrections procedure. Other approaches are possible. For example, a reasoning on the global method in one dimension shows that β W = 3 enables a direct obtention of a second-order accurate gradient at the boundary cell centroid (even if the surface integral is no more consistent). In three dimensions this method is a bit more complicated: choosing a β W such that the truncation error of order 1 cancels out is not straightforward, since the system to solve is non-linear, due to the inverse of the M * 1 (J) matrix which depends also on β W . An idea is to linearize the system and to project it in a privileged direction (normal direction for instance) in which the first-order error on the approximation gradient will cancel out. In all cases, the prerequisite for a truly k-exact reconstruction when β W 0 is to approximate the pointwise value at the wall point at least to O(h n+1 ).

C. Boundary conditions

In the previous parts, the reconstruction stencil of the k-exact schemes of order 2 and 3 have been extended for boundary cells, using the information available at the wall through boundary (Dirichlet and Neumann) conditions. An implicit assumption was made that these boundary conditions were known exactly at the wall. In practice, this is not generally the case. In the following, proper approximations of the physical wall boundary conditions are investigated for the Euler and Navier-Stokes equations.

For the scheme with 1-exact reconstruction of the solution, pointwise boundary values of the variables at x W must be known to the second order of accuracy.

For the Euler equations, a slip condition at the wall is prescribed for the velocity:

∀x ∈ A J K 0 , u(x) • n(x) = u mesh ( 30
)
with u mesh the wall velocity.

Hence, if a Dirichlet-type information is available on the normal velocity, no information is a priori known about the tangential velocity at the wall. This problem is likely to deteriorate the accuracy of the gradient approximation. An apparent solution to calculate u W would be to remove its normal part to the extrapolated value of the velocity in the boundary cell: u W = u Ju J • n J K . However, this approximation for the velocity at the wall is at most of order one.

A way of dealing with this issue is to use the computed velocity and velocity gradients of the boundary cell at the previous time to correct the wall velocity value:

u W 0 = u J + D * (1) 1 [u] J • (x W -x j ) + O(h 2 ) ( 31
)
where the right-hand side of this equation is taken at the previous timestep; and then to remove the normal part of this value

u W = u W 0 -(u W 0 • n W )n W (32)
where n W is the normal at point x W .

For the pressure, two methods are considered. The first approach consists of extrapolating the pressure to the second order using a Taylor expansion similarly to what has been done for velocity:

P W = P J + D * (1) 1 P J • (x W -x j ) + O(h 2 ) ( 33
)
When curvature is taken into account, a second approach more physical can be used, which consists in using the normal momentum equation verified for the Euler equations at the wall [START_REF] Dadone | Some Solutions of the Equations of Flow of a Viscous Compressible Fluid[END_REF]:

∂P ∂n |W = - ρu 2 tan R (34)
with R the local radius of curvature of the wall face (which depends on the tangential velocity at the wall point and on the corresponding normal as explained in Section III), u tan the tangential velocity magnitude and ρ the density. Eq.34 is discretized using the previous calculation of u W (Eqs. 31 and 32), leading to the second order formula

P W = P 1 - ρ J u 2 W R ||x W -x j || + O(h 2 ). (35)
Finally, ρ W is calculated using the gas law.

For the Navier-Stokes equations, the no-slip boundary condition is prescribed for the velocity:

∀x ∈ A J K 0 , u(x) = u mesh . ( 36
)
This can be directly applied in the discretized boundary condition for velocity. Generally, a Dirichlet-type condition is also given for temperature. Pressure, however, can only be defined at the wall through its gradient. This Neumann boundary condition is changed into an equivalent fictive Dirichlet boundary condition by performing a Taylor expansion of order 2 in the normal direction of point x W , as explained in part A.

For the scheme with 2-exact reconstruction of the solution, pointwise boundary values of the variables at x W must be known to the third order of accuracy. In a first approach, the same iterative method as previously is used.

For the Euler equations, the slip condition on velocity leads to

u W 0 = u J + D * (1) 2 [u] J • (x W -x j ) + 1 2 D * (2) 1 [u] J • (x W -x j ) ⊗2 + O(h 3 ) (37)
with the right-hand side taken at the previous timestep. Then the normal part of this value is removed:

u W = u W 0 -(u W 0 • n W )n W . (38)
For pressure

P W = P J + D * (1) 2 P J • (x W -x j ) + 1 2 D * (2) 1 P J • (x W -x j ) ⊗2 + O(h 3 ) (39)
For the Navier-Stokes equations, the velocity is again prescribed directly (no-slip condition), and so is the temperature. For the Neumann condition on the pressure, a fictive value can be evaluated by a normal Taylor expansion to the third order, using the normal gradient given by the Neumann boundary condition on the one hand, and the normal component of D * (2) 1 P J at the previous timestep on the other hand. Note that another possible approach would be to build several differential operators, one for each type of boundary condition, and a corresponding set of correction matrices. This method is however cumbersome and costly (this is even more true when the order k of the k-exact scheme increases) and is thus not considered here.

V. Numerical results

A. Inviscid flow past a smooth Gaussian bump

The first case is an inviscid flow through a channel with a Gaussian bump. The inflow Mach number is M = 0.5 with 0 degree angle of attack. Left and right boundaries are respectively a subsonic inflow, and a subsonic outflow (imposed pressure); upper and lower walls are slip walls. The flow is isentropic, thus entropy is used as an indicator of the accuracy of the simulation, since the analytical solution is unknown.

A mesh convergence study is carried out for the entropy error, using grids of 400, 1600, 6400, 25600 and 102400 cells, corresponding respectively to 40, 80, 160 320, and 640 wall faces along the bump. Results are provided for both the 1-exact and 2-exact reconstruction schemes, using the baseline solver, the solver with curved boundaries, and the solver with curved boundaries and an enriched discretization stencil at the wall.

Fig. 3 Inviscid Gaussian bump: entropy error generated downstream of the bump with the 2-exact scheme (a) with no modification; (b) with the wall curvature model

The grid convergence study for the linear reconstruction scheme is shown in Figure 4 using the L 2 norm. Firstly, we see that the convergence slope when adding the curvature model to the scheme increases compared to using the unmodified scheme. However, the entropy error reduction is low. This is not surprising: as it was stressed before, the spatial error introduced when approximating curved surfaces with flat faces is of order 2. Yet, the approximation of the boundary should at least represent the real boundary with the same order of accuracy as the discretization scheme [START_REF] Ollivier-Gooch | Obtaining and Verifying High-Order Unstructured Finite Volume Solutions to the Euler Equations[END_REF]. Here the original scheme with 1-exact reconstruction exhibits a convergence of order 2 at most, as can be seen in Figure 4. Even if the Bézier patches approximate the wall with third-order accuracy, the gain is expected to be small when this only improvement is used in the solver. When using the model with extended representation of the solution in boundary cells along with curved walls, much larger improvements are obtained: the entropy error drops and the convergence slope is again increased. The information at the wall is accurate enough thanks to the inclusion of curvature and of a proper handling of boundary conditions (even if the information available at the wall for the Euler equations is limited), and is used consistently to calculate a 1-exact approximation of the gradient.

A grid convergence study is also conducted for the scheme with 2-exact reconstruction (in Figure 5). Here the error is given in L 2 norm and in L ∞ norm. The L 2 norm measures the average error over the domain, hence giving a global error. On the contrary, the L ∞ norm gives the maximum error reached in the domain. In our case, this norm will emphasize the local error in the near-wall region, which is the focal point of this work.

Figure 3 compares the entropy error for the 2-exact reconstruction scheme on the mesh with 80 wall faces along the bump. A generation of numerical entropy is observed downstream of the bump, which is significantly reduced when adding the curvature model. In the rest of the domain where curvature is smaller, entropy is already correctly computed by the baseline method.

Figure 5 shows that very satisfying improvements are obtained with the curvature model, which approximates the physical surface with a spatial error of order 3. The entropy error is reduced and more importantly, the convergence slope is significantly increased with this model (getting closer to third order in L 2 norm), especially when considering the L ∞ norm: this confirms that the improvement of the curved model enhances the scheme near walls. Note that because the L ∞ norm is more constraining than the L 2 norm, the convergence slope obtained with this norm is smaller than the one obtained with the L 2 norm.

When considering the curvature model coupled with the parabolic reconstruction with a more symmetrical stencil in boundary cells, results are not deteriorated but no significant gain is observed. This may be explained by the fact that the boundary information available at the wall is too limited, and is therefore not sufficient to improve the overall parabolic reconstruction of the solution. To verify this, we test our models on a viscous test case, for which the boundary conditions are better posed than for inviscid flows. the flow is constant, and a Reynolds number of Re = 10 based on the inner cylinder radius is chosen as in [START_REF] Ollivier-Gooch | Unstructured High-Order Accurate Finite-Volume Solutions of the Navier-Stokes Equations[END_REF], in order to be sure that the flow does not exhibit toroidal vortices characteristic of the Taylor-Couette instability.

In this stable regime, the flow converges to a steady state for which the velocity is purely azimuthal and can be derived analytically [START_REF] Illingworth | Some Solutions of the Equations of Flow of a Viscous Compressible Fluid[END_REF]:

u θ,re f (r) = ω i R i R o r -r R o R o R i -R i R o
with ω i the angular velocity and r the radial coordinate in the polar coordinate system.

A mesh convergence study of the error on the azimuthal velocity Err u θ = |u θ -u θ, r e f | |u θ, r e f | is performed for a set of five hexahedral grids whose characteristics are presented in the table of Fig. 7. Once again, results are given for the L 2 and L ∞ norms, using the 1-exact (Fig. 8) and 2-exact (Fig. 9) reconstruction schemes. For both schemes, the results are compared for the three models (no modification, with curvature taken into account, with curvature plus enriched stencil close to the wall).

Let us first analyze the results for the scheme with linear polynomial reconstruction. The error on the azimuthal velocity on Grid 2 is presented in Figure 6, for the original solver with 1-exact reconstruction (left) and for the solver with 1-exact reconstruction coupled with the curvature model and the scheme with enriched stencil at the wall (right). It is clear in this Figure that an important error is generated on the solution close to the wall. This error is locally reduced when the wall curvature is taken into account and the modified reconstruction at boundary cells is introduced.

A grid convergence study is made for the 1-exact scheme is displayed in Figure 8. First-order convergence is observed for all of the compared models, which is due to the fact that diffusive fluxes are in general one order less accurate than convective fluxes.

Figure 8 (b) shows the L ∞ error. We see that the scheme coupled with the curvature model produces a lower error than the original scheme. However, the convergence slope appears to be slightly smaller than the one obtained by using the unmodified solver. This might be explained by the fact that the effect of wall curvature corrections is reduced when refining the mesh. Indeed a very good improvement is achieved for the coarsest grid for which the approximation of the wall by flat faces is more detrimental to the scheme.

When the scheme is enriched with this enhanced wall information, the grid convergence shows very good properties. The error is subsequently reduced, and the convergence slope is restored to a value that is actually slightly higher than with the unchanged solver. Notice that for the coarsest grid however, extending the stencil at the wall in addition to the curvature model decreases the performance of the scheme compared to the curvature model alone. This is due to the fact that asymptotical convergence has not been reached yet.

Similar features are revealed by the L 2 error curves, with less difference observed since the averaged error over all the domain is here displayed.

Results of the the grid convergence study performed for the scheme with quadratic polynomial reconstruction are reported in Fig. 9. Due to the sensitivity to the viscous fluxes near wall, the expected second-order of accuracy is not retrieved. This reveals the need to provide a special treatment for viscous fluxes, which have not been modified in the present work. However, this does not affect the relative effect of the three models, that will be analyzed in the following. Indeed, the same trends as for the 1-exact reconstruction scheme are observed again for the 2-exact reconstruction scheme asymptotically: the curvature model decreases substantially the error on the azimuthal velocity, but has a slightly smaller order of accuracy; however the correct order of accuracy is restored when the curvature model is used along with the extended stencil.

VI. Conclusion and Perspectives

High-order boundary treatment is investigated for a family of k-exact Godunov-type schemes to maintain the order of accuracy in the near-wall region, specifically for configurations with curved boundaries. Two main modifications are proposed to handle the 1-exact and 2-exact reconstruction schemes.

Firstly, a high-order composite representation of the wall based on bicubic Bézier patches is incorporated into the solver. The computation of boundary geometrical quantities is modified to conform to the patches, notably the volume and surface moments, fundamental for the reconstruction and for the fluxes' integration.

Secondly, a numerical strategy is presented for modifying the existing schemes so that the stencil is enriched with the boundary information. For the 1-exact reconstruction scheme, the Quasi-Green operator is extended by adding a wall contribution to the approximation, with a wall pointwise value and a β parameter set to 1, and by modifying the simple correction matrix accordingly. The wall point x W whose normal is colinear to the direction x jx W (x j being the boundary cell centroid) is chosen to evaluate the solution variable. This way, the same methodology can be applied for Dirichlet-type boundary conditions and for Neumann-type boundary conditions, where a fictive Dirichlet value is derived at the wall point by performing a normal Taylor expansion. For the 2-exact scheme, the approximate for the gradient with stencil incorporating the wall is applied only one time in the composition of the operators of the successive corrections procedure. The matrices used to recover consistency are modified in a similar way.

In the proposed procedure, the boundary conditions have to be prescribed at the wall point with at least the same order of accuracy as the reconstruction.

The proposed numerical treatments were first assessed for a compressible inviscid flow past a smooth Gaussian bump. For the 1-exact reconstruction scheme, the improvement due to the curvature model alone is small, but when this is used in conjunction with the extended solver, the entropy is reduced substantially and exhibits a higher convergence slope. This is due to the fact that the error commited when the wall surface is approximated with polygonal faces is already of the same order as the numerical scheme itself. On the contrary, when the 2-exact scheme is used, significant increase of the order of accuracy of the method is observed with the curvature model. In order to test the 2-exact scheme with enriched stencil close to the wall, the laminar Couette flow is then chosen, notably because a true Dirichlet condition is given on the velocity through the no-slip condition. For both schemes, the curvature model decreases the error on the azimuthal velocity. Moroever, the order of accuracy of the method is slightly increased when the curvature model is coupled with the extended stencil versions of the schemes.

These numerical results are encouraging: the modifications are consistent with the existing schemes, and the capability of the method to improve the accuracy of the flow in near-wall regions has been shown. Diffusive fluxes have not been investigated in the present work. However, the integration of these fluxes is directly affected by the reconstruction of the solution and of its derivatives on the one hand, and by the exactness of the geometrical moments on the other hand. Finally, the proposed strategy also leads to improvements for a Couette flow case. Although the diffusive fluxes close to the wall have to be computed more accurately, it does not alter the relative results obtained by using curvature corrections and the enriched stencil. Future work will comprise the modification of the calculation of diffusive fluxes in boundary cells, but also a complete study of low-Mach and low-Reynolds recentring of the scheme in near-wall regions.

Fig. 1

 1 Fig. 1 Discretization stencil at the wall and curved wall face for the: (a) 1-exact reconstruction scheme; (b) 2-exact reconstruction scheme

Fig. 2 (

 2 Fig. 2 (a) Representation of a Bézier patch; (b) Control points (P i, j ) i ∈[|1,3 |], j ∈[|1,3 |] and middle point x m generated by the solver for the wall faces of the bump test case

Fig. 4 Fig. 5

 45 Fig. 4 Mesh convergence on the L 2 norm of entropy error on the bump test case with the different models, for the 1-exact reconstruction scheme

Fig. 6

 6 Fig. 6 Plot of the error Err u θ on azimuthal velocity for Grid 2, using the 1-exact reconstruction scheme: (a) unmodified solver (b) solver augmented with wall curvature and enriched stencil for bounday cells.

Fig. 7 CircularFig. 8

 78 Fig. 7 Circular Couette flow case: (a) description of the meshes used for the grid convergence study; (b) azimuthal velocity u θ simulated for Grid 2 with the third-order accurate scheme.

Fig. 9

 9 Fig. 9 Grid convergence of the error on the azimuthal velocity for the viscous circular Couette flow: (a) 2-exact reconstruction of the variables, L 2 norm; (b) 2-exact reconstruction of the variables, L ∞ norm.


Tiago Carneiro
email: tiago.carneiro-pessoa@inria.fr

Nouredine Melab
email: nouredine.melab@univ-lille.fr

An Incremental Parallel PGAS-based Tree Search Algorithm

Keywords: High-productivity Language, PGAS, Chapel, MPI+OpenMP, Tree Search Algorithms

In this work, we show that the Chapel highproductivity language is suitable for the design and implementation of all aspects involved in the conception of parallel tree search algorithms for solving combinatorial problems. Initially, it is possible to hand-optimize the data structures involved in the search process in a way equivalent to C. As a consequence, the single-threaded search in Chapel is on average only 7% slower than its counterpart written in C. Whereas programming a multicore tree search in Chapel is equivalent to C-OpenMP in terms of performance and programmability, its productivityaware features for distributed programming stand out. It is possible to incrementally conceive a distributed tree search algorithm starting from its multicore counterpart by adding few lines of code. The distributed implementation performs load balancing among different computer nodes and also exploits all CPU cores of the system. Chapel presents an interesting tradeoff between programmability and performance despite the high level of its features. The distributed tree search in Chapel is on average 16% slower and reaches up to 80% of the scalability achieved by its C-MPI+OpenMP counterpart.

I. INTRODUCTION

Tree-based search algorithms are strategies that implicitly enumerate a solution space, dynamically building a tree. This class of algorithms is often used for the exact resolution of permutation combinatorial optimization problems (COP), and it is present in many areas, such as operations research, artificial intelligence, bioinformatics, and machine learning [START_REF] Zhang | Branch-and-bound search algorithms and their computational complexity[END_REF], [START_REF] Grama | A survey of parallel search algorithms for discrete optimization problems[END_REF]. Algorithms that belong to this class are compute-intensive and highly irregular, which demands hand-optimized data structures for efficient single-core utilization and load balancing between processes [START_REF] San Segundo | Recent developments in bit-parallel algorithms[END_REF]- [START_REF] Gmys | Ivm-based parallel branch-and-bound using hierarchical work stealing on multi-gpu systems[END_REF].

High-productivity languages historically suffer from severe performance penalties, do not provide low-level features, and are not suited to parallelism [START_REF] Costa | Exascale machines require new programming paradigms and runtimes[END_REF], [START_REF] Lusk | Languages for high-productivity computing: the darpa hpcs language project[END_REF]. Therefore, they are not often employed within the scope of parallel tree search. Instead, this kind of algorithm is frequently coded in either C or C++, and different libraries and programming models are combined for exploiting parallelism [START_REF] Carneiro | Um levantamento na literatura sobre a resoluc ¸ão de problemas de otimizac ¸ão combinatória através do uso de aceleradores gráficos[END_REF], [START_REF] Crainic | Parallel branch-and-bound algorithms[END_REF].

Among the high-productivity languages, Chapel is one that stands out. It was designed for high-performance computing, and it is competitive to both C-OpenMP and C-MPI+OpenMP in terms of performance, considering different benchmarks [START_REF] Chamberlain | Chapel comes of age: Making scalable programming productive[END_REF]. The objective of the present research is to investigate whether Chapel is suitable for the design and implementation of all aspects involved in the conception of a parallel tree search algorithm for solving combinatorial problems. To the best of our knowledge, the present research is the first one that investigates the use of a high-productivity language for this purpose.

The experimental results show that Chapel is a suitable language for the design and implementation of parallel tree search algorithms. It is possible to hand-optimize the data structures involved in the search process. As a consequence, the single-threaded search in Chapel is on average only 7% slower than its counterpart written in C. Whereas programming a tree search in Chapel is equivalent to C-OpenMP in terms of performance and programmability, its productivity-aware features for distributed programming stand out.

Thanks to Chapel's global view of the control flow and data structures, it is possible to conceive a distributed tree search starting from its multicore counterpart by incrementally adding few lines of code. The distributed implementation performs load balancing among different processes and also uses all CPU cores that a computer node has. Despite the high level of its features, the distributed tree search in Chapel is on average 16% slower and reaches up to reaches 80% of the scalability reached by its C-MPI+OpenMP counterpart. Finally, the distributed load balancing strategies provided are effective: the dynamic load balancing version is up to 1.5× faster than its static counterpart.

The remainder of this paper is structured as follows. Section II brings background information and the related works. Section III presents the incremental and PGAS-based distributed tree search algorithm. In turn, Section IV and Section V present the multicore and distributed evaluations, respectively. Next, Section VI brings a discussion of the results obtained in Sections IV and VI. Finally, conclusions are outlined in Section VII.

II. BACKGROUND AND RELATED WORKS A. The Chapel Programming Language

Chapel is an open-source parallel programming language designed to improve the programmability for highperformance computing. It incorporates features from compiled languages such as C, C++, and Fortran, as well as highlevel elements related to Python and Matlab. The parallelism is expressed in terms of lightweight tasks, which can run on several locales or a single one. In this work, the term locale refers to a symmetric multiprocessing computer in a parallel system [START_REF] Inc | Chapel language specification v.986[END_REF].

In Chapel, both global view of control flow and global view of data structures are present [START_REF] Chamberlain | Chapel comes of age: Making scalable programming productive[END_REF]. Concerning the first one, the program is started with a single task and parallelism is added through data or task parallel features. Moreover, a task can refer to any variable lexically visible, whether this variable is placed in the same locale on which task is running, or in the memory of another one. Concerning the second one, indexes of data structures are globally expressed, even in case the implementation of such data structures distributes them across several locales. Thus, Chapel is a language that realizes the Partitioned Global Address Space (PGAS) programming model [START_REF] Almasi | Pgas (partitioned global address space) languages[END_REF].

Finally, indexes of data structures are mapped to different locales using distributions. Contrasting to other PGAS-based languages, such as UPC and Fortran, Chapel also supports user-defined distributions [START_REF] Chamberlain | Userdefined parallel zippered iterators in chapel[END_REF].

B. Tree-based Search Algorithms

Tree-based search algorithms are strategies that implicitly enumerate a solution space, dynamically building a tree [START_REF] Grama | A survey of parallel search algorithms for discrete optimization problems[END_REF]. The internal nodes of the tree are incomplete solutions, whereas the leaves are solutions. Algorithms that belong to this class start with an initial node, which represents the root of the tree, i.e., the initial state of the problem to be solved. Nodes are branched during the search process, which generates children nodes more restricted than their parent node. As shown in Fig. 1, generated nodes are evaluated, and then, the valid and feasible ones are stored in a data structure called Active Set.

At each iteration, a node is removed from the active set according to the employed search strategy [START_REF] Zhang | Branch-and-bound search algorithms and their computational complexity[END_REF]. The search generates and evaluates nodes until the data structure is empty or another termination criterion is reached. If an undesirable state is reached, the algorithm discards this node and then chooses an unexplored (frontier) node in the active set. This action prunes some regions of the solution space, keeping the algorithm from unnecessary computation. The degree of parallelism of tree-based search algorithms is potentially very high, as the solution space can be partitioned into a large number of disjoint portions, which can be explored in parallel. As these algorithms are compute intensive, diverse strategies have been used for improving performance, such as instruction-level parallelism, architecture-specific code optimizations and problem-specific data structures [START_REF] San Segundo | Recent developments in bit-parallel algorithms[END_REF]- [START_REF] Gmys | Ivm-based parallel branch-and-bound using hierarchical work stealing on multi-gpu systems[END_REF], [START_REF] Feinbube | Nqueens on cuda: Optimization issues[END_REF]. Thus, parallel tree-based search algorithms are frequently written in C/C++, due to their low-level features and supported parallel computing libraries [START_REF] Carneiro | Um levantamento na literatura sobre a resoluc ¸ão de problemas de otimizac ¸ão combinatória através do uso de aceleradores gráficos[END_REF]. In the context of distributed algorithms, the performance-aware strategies above mentioned are combined with distributed programming libraries for implementing load balancing and explicit communication between processes [START_REF] Crainic | Parallel branch-and-bound algorithms[END_REF], [START_REF] Tschoke | Solving the traveling salesman problem with a distributed branch-and-bound algorithm on a 1024 processor network[END_REF], [START_REF] Mezmaz | A grid-enabled branch and bound algorithm for solving challenging combinatorial optimization problems[END_REF]. As a consequence, programming distributed tree search algorithms can be challenging and time-consuming.

III. PARALLEL TREE-BASED SEARCH ALGORITHMS IN CHAPEL

A major objective of Chapel concerning productivity is allowing distributed programming using concepts close to the ones of shared-memory programming [START_REF] Chamberlain | Chapel comes of age: Making scalable programming productive[END_REF]. In this section, a multicore and single-locale tree search algorithm is initially proposed. Then, it is incrementally extended using Chapel's productivity-aware features for distributed programming.

A. Algorithm Overview

This work focuses on permutation combinatorial problems, for which an N -sized permutation represents a valid and complete solution. Permutation combinatorial problems are used to model diverse real-world situations, and their decision versions are often NP-Complete [START_REF] Zhang | Branch-and-bound search algorithms and their computational complexity[END_REF], [START_REF] Mezmaz | A grid-enabled branch and bound algorithm for solving challenging combinatorial optimization problems[END_REF].

This section presents two backtracking algorithms for enumerating all complete and feasible solutions of the N-Queens. Backtracking is a fundamental problem-solving paradigm that consists in dynamically enumerating a solution space in a depth-first fashion. Due to its low memory requirements and its ability to quickly find new solutions, depth-first search (DFS) is often preferred [START_REF] Zhang | Branch-and-bound search algorithms and their computational complexity[END_REF].

The N-Queens problem consists in placing N non-attacking queens on a N × N chessboard, and it is often used as a benchmark for novel tree-based search algorithms [START_REF] Feinbube | Nqueens on cuda: Optimization issues[END_REF], [START_REF] Bell | A survey of known results and research areas for n-queens[END_REF]. The N-Queens is easily modeled as a permutation problem: position r of a permutation of size N designates the column in which a queen is placed in row r. Furthermore, the concepts herein presented are similar to any permutation combinatorial problem and can be adapted for solving other problems of this class with straightforward modifications [START_REF] Carneiro Pessoa | GPU-accelerated backtracking using CUDA dynamic parallelism[END_REF], [START_REF] Gmys | Ivm-based parallel branch-and-bound using hierarchical work stealing on multi-gpu systems[END_REF].

B. The Single-locale Multicore Implementation

Algorithm 1 presents a pseudocode for the single-locale backtracking in Chapel. The algorithm starts receiving the problem to be solved (line 1) and the cutoff depth (line 2). Then, it is required to generate an initial load for the parallel search. For this purpose, task 0 performs backtracking from depth 1 (initial problem configuration) until the cutoff depth cutof f , storing all feasible, valid, and incomplete solutions at depth cutof f in the active set A (line 4). After generating the initial load, the parallel search strategy begins through a forall statement (line 5).

As one can see in Fig. 2, nodes in the centralized active set A are assigned to tasks in chunks. Each task has its active set and executes a backtracking search strategy. In turn, nodes are used to initialize the backtracking, which enumerates the solution space rooted by a node. The load balancing is done through the iterator (DynamicIters) used to assign indexes of A to tasks, like in OpenMP.

Metrics are reduced through Reduce Intents. In Chapel, it is possible to use the Tuple data type (equivalent to Cstructs) and reduce all metrics at once (line 6). Differently from OpenMP, it is not required to define a tuple reduction. Finally, the parallel search finishes when the active set A is empty.

Algorithm 1: The multicore tree search algorithm.

1 I ← get problem() 2 cutof f ← get cutof f depth() 3 A ← ∅ 4 A ← generate initial active set(cutof f, I) 5 forall node in A with(+ reduce metrics) do 6 metrics+ = tree search(node, cutof f, I) 7 end

C. The Multi-locale Implementation

One can see in Algorithm 2 a pseudocode for the distributed tree-based search algorithm in Chapel. Thanks to Chapel's global view of control flow, the search also starts serially, with task 0 generating the initial load to populate the active set A (line 4). To make it possible to distributed the nodes of A across several locales, it is required to define a domain (line 5) and to indicate how the indexes of this domain are mapped across different locales (line 6). In this work, only standard distributions are used 1 . Finally, the distributed active set A d of type Node is defined over the mapped domain D (line 8).

After the initial load generation, the nodes of A are distributed by using a parallel forall (line 9), which generates the distributed active set A d . Thanks to Chapel's global view of A d , the indexes of both active sets are directly accessed in line 10. Moreover, as shown in Fig. 3, A d is an abstraction. The distributed active set A d consists of several sets A i d , i ∈ {0, ..., l -1}, where l is the number of locales on which the application is going to run.

The parallel search takes place in line 12. As one can see in Algorithm 2, its forall is similar to the one of Algorithm 1. However, distributed iterators are used instead (DistributedIters). Additionally, the distributed search exploits two levels of parallelism, and the compiler is also responsible for generating the code that exploits all CPU cores a locale has. Finally, the metrics are reduced in the same way as in the single-locale algorithm.

Algorithm 2: The multi-locale tree search algorithm.

1 I ← get problem() 2 cutof f ← get cutof f depth() 3 A ← ∅ 4 A ←

D. Search Procedure and Data Structures

The kernel of both parallel algorithms previously presented is based on a serial and hand optimized backtracking for solving permutation combinatorial problems, originally written in C [START_REF] Carneiro Pessoa | GPU-accelerated backtracking using CUDA dynamic parallelism[END_REF]. The serial backtracking was then adapted to Chapel, obeying the handmade optimizations, instruction-level parallelism, data structures, and C-types.

The data structure Node is similar to any permutation combinatorial problem. It contains an unsigned 8-bit integer vector of size cutof f , identified by board, and an unsigned integer variable. The vector board stores the feasible and valid incomplete solution. In turn, the integer variable, identified by bitset, keeps track of board lines by setting its bit n to 1 each time a queen is placed in the n-th line. The search performed by the kernel (Algorithm 1, line 6 and Algorithm 2, line 13) is a non-recursive backtracking that does not use dynamic data structures, such as stacks. Initially depth receives the value of cutof f . Next, board and bitset are initialized with the incomplete solution that Node[i] contains.

The semantics of a stack is obtained by using a variable depth and by trying to increment the value of the vector board at position depth. If this increment results in a feasible and valid incomplete solution, the depth variable is then incremented, and the search proceeds to the next depth. After trying all configurations for a given depth, the search backtracks to the previous one.

IV. A SINGLE-LOCALE PERFORMANCE EVALUATION OF CHAPEL

The primary objective of this section is to investigate the single-locale programming features and performance of Chapel.

A. Protocol

For this evaluation, the following programs were conceived for enumerating all valid and complete solutions of the N-Queens problem.

• Multicore: Chapel and C-OpenMP implementations of the backtracking search algorithm described in Section III-B. • Serial: Chapel and C implementations corresponding to the kernel of the multicore programs above listed (refer to Section III-D). All implementations apply the data structures and search procedure detailed in Section III-D.

B. Parameters Settings

In the experiments, both multicore and serial implementations enumerate all complete and valid solutions of the N-Queens problem, for which sizes (N ) range from 10 to 19. The experiments take from few milliseconds to several hours of parallel processing.

The testbed Operates under SMP Debian 4.9.65 64 bits, and it is composed of two AMD EPYC 7301, with 32 cores @2.7 GHz, 64 threads, and 128 GB RAM. All C programs were compiled with gcc 6.3.0. The Chapel version used was 1.8.0.

Chapel provides three task layer implementations: qthreads (default), Tokio's University Massive Threads, and POSIX Threads (Pthreads). A preliminary experiment was performed to verify which thread implementation is the most advantageous in the context of this work. It is important to point out that the task layer is chosen in terms of environment variables and this action means no coding efforts. Fig. 4 shows that both massive threads and Pthreads are much heavier than qthreads for the smaller tree sizes. All task layer implementations perform similarly as the size of the solution space grows. As OpenMP, Chapel makes available various load balancing strategies, which are implemented as built-in iterators used in forall statements. They are close to OpenMP's scheduling policies, such as guided and dynamic. A preliminary experiment was carried on to figure it out the best Chapel's built-in load balance strategy for solving the N-Queens. It is shown in Fig. 5 the average execution time required by the multicore Chapel backtracking to solve the N-Queens problem, taking into account different built-in load balance strategies. According to the results, the dynamic approach is the fastest one.

Experiments were also carried out to choose a suitable cutoff depth for both Chapel and OpenMP implementations.

One can see in Table I the best parameters experimentally found for the Chapel and C-OpenMP implementations.

C. Results

One can see in Fig. 6 a comparison between Chapel and C serial implementations. As previously pointed out, it is possible to write in Chapel a hand-optimized code similar to C. For sizes ranging from 10 to 11, the serial implementation in Chapel is on average 1.4× and 1.21× slower than its counterpart written in C, respectively. However, as the size of the problem grows, this difference becomes much smaller. In turn, taking into account problem sizes ranging from 12 to 18, the Chapel serial implementation is on average 7% slower than its C counterpart.

It is shown in Fig. 4 the average execution time of spent by Chapel for solving the N-Queens compared to its C-OpenMP counterpart. The results also consider all task implementations of Chapel and use the best parameters found, summarized in Table I. The version running over the qthreads task layer is comparable to C-OpenMP even for the smallest sizes (10 to 13). For sizes ranging from 14 to 18, the version over qthreads is on average 8% slower than the search in C-OpenMP.

Both massive threads and Pthreads task layer implementations contrast to qthreads. For these two task layers, the overhead of managing threads amounts negatively when enumerating small solution spaces, and they perform poorly for sizes ranging from 10 to 13. The massive threads version is from 29× (N = 10) to 2.52× (N = 13) slower than its C-OpenMP counterpart. In turn, the implementation over Pthreads is from 16× (N = 10) to 1.8× (N = 13) slower than its counterpart written in C-OpenMP. As the size of the solution space grows, both Pthreads and massive threads version stand out. From sizes ranging from 15 to 19, both implementations are on average on average 5% faster than C-OpenMP and 13% faster than its counterpart over qthreads.

V. A MULTI-LOCALE PERFORMANCE EVALUATION OF CHAPEL

In this section, the incrementally conceived distributed algorithm presented in Section III-C is evaluated. The primary goal of this section is to show that it is possible to use a highproductivity language for programming distributed tree search algorithms and achieve metrics similar to MPI+X.

A. Protocol

The following applications were programmed for enumerating all valid and complete configurations of the N-Queens problem.

• Chapel: implementation of the multi-locale backtracking search algorithm described in Algorithm 2, written in Chapel.

• MPI+X: single program -multiple data (SPMD) counterpart written in C of the program above introduced. In this case, MPI is applied for communication, and X means the use of OpenMP for exploiting all CPU cores/threads a locale has. Both applications implement the data structures and search procedure detailed in Section III-D.

In this evaluation, it is investigated how the applications scale according to the number of locales. Furthermore, the influence of the PGAS data structure distribution on the application execution time is also studied. Moreover, the impact of the distributed load balancing strategies on the overall performance of the application is also investigated. Finally, all metrics collected for the implementation in Chapel are compared to the ones achieved by its MPI+X counterpart.

B. Parameters Settings

Problems of size (N ) ranging from 15 to 20 are considered. The experiments take from few seconds to several hours of parallel processing. The number of locales ranges from 1 to 32, and the application is the same for either one or more than one computer node(s). The number of locales is passed to the application using Chapel's built-in command line parameter -nl l (-np l for MPI), where l is the number of locales on which the application is executed.

All computer nodes are symmetric and operate under Debian 4.9.130 -2, 64 bits. They are equipped with two Intel Xeon X5670 @ 2.93 GHz (a total of 12 cores/24 threads), and 96 GB RAM. Thus, up to 384 cores/768 threads are used in the experiments. All locales are interconnected through an Infiniband network: Mellanox Technologies MT26428 (Con-nectX VPI PCIe 2.0 5GT/s -IB QDR / 10GigE).

The Chapel implementation was programmed in its current version (1.18), and the default task layer (qthreads) is the one employed. Chapel's multi-locale code runs on top of GASNet, and several environment variables should be set with the characteristics of the system the multi-locale code is supposed to run. Concerning the MPI+X implementation, OpenRTE 2.0.2 along with gcc 6.3.0 and OpenMP 4.5 were used for compilation and execution. One can see in Table II a summarization of the runtime configurations for multi-locale execution. The Infiniband GASNet implementation is the one used for communication (CHPL_COMM_SUBSTRATE) along with MPI, which is responsible for getting the executables running on different locales (GASNET_IBV_SPAWNER).

Chapel provides several standard distributions to map data structures onto locales. Different tests were also carried out to identify the best option in the context of this work. The one chosen was the one-dimension BlockDist, which horizontally maps elements across locales. For instance, in case l = 3 and |A d | = 8, elements 0, ..., 2 are on locale l 0 , 3, .., 5 on locale l 1 , and 6, 7 on locale l 2 . In the scope of the present research, choosing a different standard distribution does not lead to performance improvements.

Experiments were carried out to choose a suitable cutoff depth (Algorithm 2, line 2). This parameter directly influences the size of A d , and therefore the time spent in distributing the active set across locales. As observed in Fig. 7, the fastest data structure distribution is observed for cutof f = 3. However, such a cutoff value limits parallelism, resulting in a slow distributed search. In contrast, when the cutoff is set to 6, the distribution of A d becomes 10× slower than the search procedure itself. This behavior happens due to the combinatorial nature of N-Queens: a cutoff depth twice deeper results in an active set 725× bigger. When choosing cutof f = 5, the search takes the same time as for cutof f = 4. Despite that, the distribution of A d is on average 9× slower for cutof f = 5. Thus, the cutoff depth chosen is 4. Preliminary experiments also show that cutof f = 4 is the best value for the MPI+X implementation.

Chapel also provides two different distributed load balancing iterators: guided and dynamic, which are also similar to OpenMP's schedules of the same name. Experiments were carried out to identify the best chunk for both load balancing strategies. They present the best performance when using the default chunk size.

C. Results

First of all, the benefits of using distributed load balancing are not observed for the smallest solution space, i.e., for the problem of size N = 15. In such a situation, the static search performs slightly better because there is no communication among locales during the search. As shown in Fig. 8, the overhead of data structure initialization and distribution becomes less detrimental as the solution space grows, and the benefits of using distributed load balancing can be observed. For sizes bigger than 15, using the dynamic iterator is from 1.17× to 1.51× times faster than using no load balancing (static version). Moreover, the guided iterator does not seem a suitable load balancing in the scope of this work: it shows benefits compared to the static version only for sizes ranging from 18 to 20. For these problem sizes, using the guided iterator makes the search up to 1.21× faster than its static counterpart. In turn, using the dynamic distributed iterator results in a search from 1.21× to 1.25× faster than using the guided one.

It is shown in Fig. 9 how the distributed searches in Chapel and MPI+X scales according to the number of locales. The worst scalability is observed for the smallest size (N = 15). In such a situation, the initialization and distribution of A d amount for almost the whole execution time (see Fig. 8). For problem sizes ranging from 17 to 20, the dynamic version scales up to 20.5× (N = 19), whereas guided and static scale up to 16.8× and 16.9×, respectively (also N = 19). The MPI+X version scales up to 25.4× (N = 18). Therefore, the distributed search in Chapel achieves up to 80% of the scalability observed for its MPI+X counterpart. It is worth to mention that the time spent on distributing A d does not grow linearly according to the number of locales, as shown in Fig. 10. The time required to distribute A d grows up to size N = 16, then it becomes almost constant. This behavior comes from the fact that the size of A d is the same for one or more locale(s). Thus, as the number of locales grows, the number of messages sent grows as well, but their size decreases. Moreover, the A d distribution is performed in parallel (Algorithm 2, line 9), and the Infiniband GASNet implementation supports one-sided communication.

In terms of wall-clock time, Chapel is equivalent to MPI+OpenMP when running on one locale. For the smaller solution space (N = 15), Chapel stands out, and it is up to 25% faster than MPI+X. In such a situation, A d is not distributed, and the program behaves like a single-locale and multicore one. Moreover, MPI implements the SPMD programming model. This way, MPI is started, and its functions are called even for one locale. Additionally, it is worth to mention that Chapel is a compiled language and it is possible to program in Chapel both search strategy and data structures equivalent to the ones present in its counterpart written in C. In contrast, for multiple locales and bigger problem sizes, the Chapel distributed search is on average 16% slower than its MPI+X counterpart.

VI. DISCUSSION

In this work, all aspects of the search process were programmed in Chapel, even though C code can be incorporated into a Chapel program. It was possible to hand-optimize the search kernel in a way similar to C. Both codes are equivalent in terms of types, data structures and code size, which resulted in a single-threaded performance competitive to C. This fact is essential in the context of this work, otherwise using the parallel features of Chapel would result in low-performance.

Programming a multicore search in Chapel involves almost the same effort as using C-OpenMP. Both provide built-in load balancing features and reduction of variables. However, Chapel presents some advantages, as it provides more load balancing strategies, and it is possible to reduce all metrics at once. Moreover, there are several task layer implementations, which may be advantageous for some users. Concerning performance, the multicore Chapel implementation using the default task layer is competitive to C-OpenMP even when solving the smallest problems. Thanks to Chapel's global view of the control flow and data structures, the main difference between the multi-and single-locale versions lies mainly in the use of the PGAS data structures and distributed iterators for load balancing. There is no need for explicitly dealing with communication, metrics reduction, or distributed load balancing. Furthermore, the compiler generates code for exploiting all CPU cores a locale has. Differently from the classic MPI+X, there is no need for an additional library to exploit each level of parallelism.

Concerning the program size, Chapel's multi-locale implementation is only 8 lines longer than its single-locale counterpart, which results in a code 33% bigger. Consequently, the two communication behaviors presented in Fig. 11 are achieved by the same program, but different parameters. In contrast, it is required to add 24 lines to the backtracking written in C-OpenMP to use MPI, which almost doubles the program size and also incorporates the SPMD programming model. Therefore, Chapel presents an interesting trade-off between programmability and performance.

The most significant limitations found concern neither programmability nor performance. Instead, they are related to technical issues. For instance, it took much more time to configure the GASNet library for running on a cluster than programming the multi-locale backtracking itself. In our case, a modification in the GASNet source code was necessary to run the Chapel distributed search on an MXM network with a non-default partition key. This problem would keep a not so enthusiastic user from Chapel. The bright side is that it was not a Chapel-only effort, as other PGAS libraries, such as UPC, Fortran, SHMEM use GASNet as communication layer.

VII. CONCLUSION

This work has investigated the use of Chapel highproductivity language for the design and implementation of all aspects involved in the conception of parallel tree search algorithms. This research covered from instruction-level parallelism used to improve the single-threaded search to the distributed and multi-level parallelism. According to the results, Chapel is a suitable language for programming such a complex and compute-intensive application. It is possible to handoptimize the data structures involved in the search process in a way equivalent to C. Moreover, Chapel's multicore features are similar to OpenMP. Additionally, programmers familiarized with shared memory programming can incrementally conceive a multi-level and distributed tree search.

Chapel presented an interesting trade-off between performance and programmability, despite the high level of its features for distributed programming. One would argue that it could be possible to program an MPI+X version faster than the one used; however, that is also the case for Chapel. For instance, the code for exploiting all CPU cores a locale has could be programmed by hand, as well as the communication and load balancing among locales. However, the latter does not seem necessary, as the use of high-productivity features resulted in performance competitive to MPI+OpenMP.

It is worth to point out that the parallel optimization community already possesses legacy code mainly written in C/C++. Therefore, programmers may be resistant to learn another language and translate programs to Chapel [START_REF] Lusk | Languages for high-productivity computing: the darpa hpcs language project[END_REF]. The capacity of Chapel to include C code can be a partial solution for this situation. One could use C could along with Chapel's high-productivity features for distributed programming. Finally, graphics processing units are crucial for solving big and challenging combinatorial optimization problems [START_REF] Gmys | Ivm-based parallel branch-and-bound using hierarchical work stealing on multi-gpu systems[END_REF]. The adoption of Chapel by the parallel optimization community, besides performance and productivity, also may also depend on the support of GPUs.

Fig. 1 .

 1 Fig. 1. Visual representation of a tree-based search algorithm (Own representation based on [9]).

  1 https://chapel-lang.org/docs/modules/layoutdist.html

Fig. 2 .Fig. 3 .

 23 Fig.2. Task 0 is responsible for managing the centralized active set A and performing load balancing. The searches are independent, and metrics are reduced using the Reduce Intents of Chapel (Own representation adapted from[START_REF] Crainic | Parallel branch-and-bound algorithms[END_REF]).

  generate initial active set(cutof f, I) 5 Space ← {0..(|A| -1)} 6 D ← Space mapped according to a standard distribution 7 A d ← ∅ 8 A d ← [D] : N ode 9 forall s in Space do 10 A d [s] ← A[s] 11 end 12 forall node in A d following the iterator with(+ reduce metrics) do 13 metrics+ = tree search(node, cutof f, I) 14 end

Fig. 4 .Fig. 5 .

 45 Fig. 4. Average Chapel execution time of compared to its C-OpenMP counterpart for silving the N-Queens problem. Both implementations use the parameters of Table I Results are shown for all three Chapel task layer implementations: qthreads, massive threads and Pthreads. Problem sizes (N) range from 10 to 19.

Fig. 6 .

 6 Fig. 6. Average execution time of the serial backtracking in Chapel compared to its counterpart written in C. Problem sizes (N) range from 10 to 18.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Influence of the cutof f choice on the execution time. Values are for the N-Queens of size N = 17 and 32 locales. On the X axis: the cutoff depth and the distributed active set size in parentheses (in 10 6 nodes).

Fig. 9 .

 9 Fig. 9. The highest speedup achieved by Chapel and MPI+X implementations when executed on 2 to 32 locales.

Fig. 10 .

 10 Fig. 10. Normalized time required to initialize and distribute the PGAS-based active set (A d ). Results are for 2 to 32 locales.

Fig. 11 .

 11 Fig. 11. Output of Chplvis [18] for (a) static load distribution (b) Dynamic distributed iterator. Results are for size N = 18 running on 16 locales. In the figure, when the color of an edge {x, y} is red closer to x means that much more communication happen from x to y rather than y to x.

TABLE I LIST

 I OF BEST PARAMETERS FOUND EXPERIMENTALLY FOR THE CHAPEL AND C-OPENMP IMPLEMENTATIONS.

			Parameters Settings
	Implementation	Load Balancing	Chunk	Optimization cutof f
	Chapel	Dynamic	def ault	--fast	4
	C -OpenM P	Dynamic	def ault	-O3	4

TABLE II SUMMARIZATION

 II OF THE ENVIRONMENT CONFIGURATION FOR MULTI-LOCALE EXECUTION AND COMPILATION.

	Variable	Value
	CHPL_RT_NUM_THREADS_PER_LOCALE	24
	CHPL_TARGET_ARCH	native
	CHPL_COMM	gasnet
	CHPL_COMM_SUBSTRATE	ibv
	GASNET_IBV_SPAWNER	mpi

http://www.grid5000.fr we also thank Paul Hargrove for the modifications in GAS-Net infiniband implementation necessary to run GASNet on GRID'5000 MXM infiniband networks.

ACKNOWLEDGMENTS

The experiments presented in this paper were carried out on the Grid'5000 testbed [19], hosted by INRIA and including several unother organizations 2 . We thank Bradford Chamberlain, Elliot Ronaghan (Cray inc.) and Paul Hargrove (Berkeley lab.) for helping us to run GASNet on GRID5000. Moreover,