D De Santis
email: dante.de_santis@inria.fr

G Geraci

A Guardone
email: alberto.guardone@polimi.it

Finite Volume and Finite Element Schemes for the Euler Equation in Cylindrical and Spherical Coordinates

Keywords: Compressible flows, Cylindrical/Spherical Coordinates, Explosion problem, Finite Element/Volume Methods

A numerical scheme is presented for the solution of the compressible Euler equations over unstructured grids in cylindrical and spherical coordinates. The proposed scheme is based on a mixed finite volume / finite element approach. Numerical simulations are presented for the explosion problem in two spatial dimensions in cylindrical and spherical coordinates, and the numerical results are compared with the one-dimensional simulations for cylindrically and spherically symmetric explosions.

Introduction

Some gasdynamics problems exhibit such strong symmetries that could be more convenient to describe these problems in a curvilinear reference, i.e. cylindrical or spherical, rather than in a Cartesian one. These are, e.g., astrophysical flows, Inertial Confinement Fusion (ICF) applications, sonoluminescence phenomena and nuclear explosions [START_REF] Sedov | Similarity and dimensional methods in mechanics[END_REF].

To compute the numerical solution of the compressible flow equations for these kind of flows, an interesting possibility is provided by the use of a mixed finite volume (FV) / finite element (FE) approach [START_REF] Fezou | A class of implicit upwind schemes for Euler simulations with unstructured meshes[END_REF]. For example, in viscous flows, it is possible to use the FV and the FE to compute the advection and dissipation terms, respectively, within the same algorithm. Such a possibility is expected to be of use in the study of the effect of viscosity on e.g. the formation of stable shock fronts and on the determination of the onset and dynamics of Richtmyer-Meshkov instabilities in cylindrical and spherical implosions [START_REF] Chen | Stability of imploding spherical shock waves[END_REF].

The combined use of these two different techniques is made possible by the introduction of suitable equivalence conditions that relate the FV metrics, i.e. cell volumes and integrated normals, to the FE integrals. Equivalence conditions relating FV and FE schemes have been derived for Cartesian coordinates in two and three spatial dimensions [START_REF] Selmin | The node-centred finite volume approach: bridge between finite differences and finite elements[END_REF][START_REF] Selmin | Unified construction of finite element and finite volume discretizations for compressible flows[END_REF] and for cylindrical coordinates in axially symmetric two-dimensional problems [START_REF] Guardone | Finite element/volume solution to axisymmetric conservation laws[END_REF]. In both cited references, equivalence conditions are obtained by neglecting higher order FE contributions. Subsequently in [START_REF] Santis | Equivalence conditions for finite volume/element discretizations in cylindrical coordinates[END_REF], equivalence conditions for the cylindrical coordinates have been derived for the first time without introducing any approximation into the FE discrete expression of the divergence operator, and in [START_REF] Guardone | On the relation between finite element and finite volume schemes for compressible flows with cylindrical and spherical symmetry[END_REF] the differences between the consistent scheme and an alternative one violating the the equivalence conditions have been quantified for the case of one-dimensional problems in cylindrical and spherical coordinates. In the present paper the consistent formulation is presented in the cylindrical and spherical reference systems. Numerical results for the implosion and the explosion problem are also provided to demonstrate the correctness of the proposed approach.

The present paper is structured as follows. In section 2, the scalar equations in cylindrical and spherical coordinates are introduced. In sections 3 and 4 the spatial discretizations of the scalar equation by FE and FV are presented in cylindrical and spherical coordinates respectively. Equivalence conditions between FE and FV metrics are also shown. In section 6, numerical simulations are presented for the implosion and the explosion problem in the cylindrical coordinates on the R-θ and Z-R (axisymmetric) planes and in the spherical coordinates on the r-φ plane. Numerical results are also compared to one-dimensional simulations. In section 7 final remarks and comments are given.

Scalar equation in cylindrical and spherical coordinates

For simplicity we consider first a multidimensional scalar conservation law. The model equation in a threedimensional cylindrical reference reads

∂u ∂t + ∂ f Z ∂Z + 1 R ∂ ∂R (R f R ) + 1 R ∂ f θ ∂θ = 0, ( 1
)
where t is the time, Z, R and θ are the axial, radial and azimuthal coordinates, respectively, u = u(Z, R, θ, t) is the scalar unknown and

f ∅ (u) = ( f Z , f R , f θ )
is the so-called flux function. A more compact form of the above equation is obtained by introducing the divergence operator in three-dimensional cylindrical coordinates ∇ ∅ • (•) as follows

∂(u) ∂t + ∇ ∅ • f ∅ (u) = 0. ( 2
)
The model equation in a three-dimensional spherical reference reads

∂u ∂t + 1 r 2 ∂ ∂r (r 2 f r ) + 1 r sin θ ∂ ∂θ (sin θ f θ ) + 1 r sin θ ∂ f φ ∂φ = 0,
where r, θ and φ are the radial, polar and azimuthal coordinates, respectively, u = u(r, θ, φ, t) is the scalar unknown and f ⊙ (u) = ( f r , f θ , f φ ) is the flux function. In a compact form, using the divergence operator in three-dimensional spherical coordinates ∇ ⊙ • (•), it can be rearranged as

∂(u) ∂t + ∇ ⊙ • f ⊙ (u) = 0. ( 3
)

Finite Volume/Element in cylindrical coordinates

Node-pair finite element discretization

The scalar conservation law ( 2) is now written in a weak form by multiplying it by the radial coordinate R and by a suitable Lagrangian test function

φ i ∈ V h ⊂ H 1 (Ω). Integrating over the support Ω i of ϕ i gives Ω i Rϕ i ∂u ∂t dΩ i + Ω i Rϕ i ∇ ∅ • f ∅ (u) dΩ i = 0, ∀i ∈ N, (4)
where N is the set of all nodes of the triangulation. Note that by multiplying by R, the numerical singularity of the cylindrical reference system is formally removed [START_REF] Guardone | Finite element/volume solution to axisymmetric conservation laws[END_REF]. In the following, to simplify the notation, the infinitesimal volume dΩ = R dR dθ dZ is not indicated in the integrals. Integrating by parts immediately gives

Ω i Rϕ i ∂u ∂t = Ω i R f ∅ • ∇ ∅ ϕ i + Ω i ϕ i f ∅ • ∇ ∅ R - ∂Ω ∂ i Rϕ i f ∅ • n ∅ i ( 5
)
where ∂Ω ∂ i = ∂Ω i ∩ ∂Ω, with ∂Ω i and ∂Ω are the boundary of Ω i and of the computational domain Ω, respectively, and where n ∅ i = n Z Ẑ + n R R + n θ θ is the outward normal versor to Ω i . The scalar unknown is now interpolated as

u(Z, R, θ, t) ≃ u h (Z, R, θ, t) = k∈N u k (t)ϕ k (Z, R, θ),
to obtain the Bubnov-Galerkin approximation of the (1), namely

k∈N i du k dt M ∅ ik = Ω i R f ∅ (u h ) • ∇ ∅ ϕ i + Ω i ϕ i f ∅ (u h ) • ∇ ∅ R - ∂Ω ∂ i Rϕ i f ∅ (u h ) • n ∅ , (6)
where N i is the set of shape functions ϕ k whose support Ω k overlaps Ω i of ϕ i and where

M ∅ ik def = Ω ik Rϕ i ϕ k dΩ ∅ , with Ω ik = Ω i ∩ Ω k .
By resorting the so-called flux interpolation technique [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF], the flux function f ⊙ (u h ) is now expanded using the same shape function ϕ h ∈ V h as follows

f ∅ u h = f ∅         k∈N u k (t)ϕ k (Z, R, θ)         ≃ k∈N f ∅ k (t)ϕ k (Z, R, θ),
where

f ∅ k (t) = f ∅ (u k (t)), to obtain k∈N i M ∅ ik du k dt = k∈N i f ∅ k (t) • Ω ik Rϕ k ∇ ∅ ϕ i + Ω ik ϕ i ϕ k ∇ ∅ R - k∈N ∂ i f ∅ k (t) • ∂Ω ∂ ik Rϕ i ϕ k n ∅ i ,
where ∂Ω ∂ ik = ∂Ω i ∩ ∂Ω k ∩ ∂Ω and N ∂ i is the set of all boundary nodes of Ω i . Using the following identities

k∈N i f ∅ k • Ω ik Rϕ k ∇ ∅ ϕ i = k∈N i , f ∅ k + f ∅ i 2 • η ∅ ik - f ∅ k -f ∅ i 2 • ζ ∅ ik + k∈N i , f ∅ k -f ∅ i 2 • χ ∅ ik , k∈N i f ∅ k • ∂Ω ∂ ik Rϕ i ϕ k n ∅ = k∈N ∂ i, f ∅ k -f ∅ i • χ ∅ ik -f ∅ i • ξ ∅ i ,
demonstrated in Appendix A, the node-pair representation of ( 6) is find to be

k∈N i M ∅ ik du k dt = - k∈N i , f ∅ k + f ∅ i 2 • η ∅ ik - f ∅ k -f ∅ i 2 • ζ ∅ ik + f ∅ i • L ∅ i - k∈N ∂ i, f ∅ k -f ∅ i 2 • χ ∅ ik -f ∅ i • ξ ∅ i , (7)
where N i, = N i \{i} and N ∂ i, = N ∂ i \{i}, and the following FE metric quantities have been introduced

η ∅ ik def = Ω ik R(ϕ i ∇ ∅ ϕ k -ϕ k ∇ ∅ ϕ i ), ζ ∅ ik def = Ω ik ϕ i ϕ k R, L ∅ i def = k∈N i Ω ik ϕ i ϕ k R = Ω ik ϕ i R, χ ∅ ik def = ∂Ω ∂ ik Rϕ i ϕ k n ∅ and ξ ∅ i def = ∂Ω ∂ ik Rϕ i n ∅ .
In the next section, the corresponding FV metrics are derived.

Edge-based finite volume discretization

The spatially discrete form of the scalar conservation law ( 2) is now obtained according to the node-centred finite volume approach [START_REF] Leveque | Finite volume methods for conservation laws and hyperbolic systems[END_REF]. To this purpose, the integral form of (2) times the radial coordinate R is enforced over a finite number of non-overlapping finite volumes C i , with boundary ∂C i . Over each control volume C i the cell-averaged unknown is introduced as follows

u(Z, R, θ, t) ≃ u i (t) = 1 V i C i u(Z, R, θ, t),
where V i is the volume of the i-th cell. Splitting the boundary integral into interface and edge contributions

V ∅ i du i dt = - k∈N i, ∂C ik R f ∅ • n ∅ i - ∂C ∂ i R f ∅ • n ∅ i + C i f ∅ • R, (8)
where N i, is the set of the finite volume C k sharing a boundary with C i , excluding C i and where ∂C ik = ∂C i ∩ ∂C k ∅, k i, is the so-called cell interface. As it is standard practice, the flux vector is assumed to be constant over each cell interface. Under this assumption, the domain and boundary contributions read

∂C ik R f ∅ • n ∅ i ≃ f ∅ ik • ∂C ik Rn ∅ i = f ∅ ik • ν ∅ ik with ν ∅ ik def = ∂C ik Rn ∅ i and ∂C ∂ i R f ∅ • n ∅ i ≃ f ∅ i • ∂C ∂ i Rn ∅ i = f ∅ i • ν ∅ i with ν ∅ i def = ∂C ∂ i Rn ∅ i ,
respectively. If a second-order centred approximation of the fluxes is considered, namely,

f ∅ ik = ( f ∅ i + f ∅ k )/2
, the final form of the finite volume approximation of (2) reads,

V ∅ i du i dt = - k∈N i, f ∅ i + f ∅ k 2 • ν ∅ ik + f ∅ i • V ∅ i -f ∅ i • ν ∅ i with V ∅ i def = C i R, (9)
to be compared with the corresponding FE discretization (7).

Finite Element/Volume equivalence

The equivalence conditions relating the above FV metric quantities and the FE ones defined in the previous section are now derived. To this purpose, relevant properties of the FE and FV discretizations are briefly recalled.

Considering FE metric quantities first, from its definition the vector η ∅ ik is asymmetric, namely, η ∅ ik = -η ∅ ki which will be referred in the following as property FEM-a. Property FEM-b is obtained by noting that

k∈N i η ∅ ik -ζ ∅ ik +ξ ∅ i = 0 which gives immediately L ∅ i = k∈N i, η ∅ ik + ξ ∅ i . (10)
Property FEM-c stems from the following identity

3L ∅ i = Ω i Rϕ i ∇ ∅ • x ∅ = ∂Ω ∂ i Rϕ i x ∅ • n ∅ i - Ω i ϕ i x ∅ • R - Ω i Rx ∅ • ∇ ∅ ϕ i , (11)
where x ∅ is the position vector, x ∅ = Z Ẑ +R R(θ). By substituting the exact expansion x ∅ = k∈N i x ∅ k ϕ k , and applying node-pair representation described in the subsection 4.1, Eq. ( 11) can be written as

3L ∅ i = k∈N i, x ∅ k + x ∅ i 2 • η ∅ ik - x ∅ k -x ∅ i 2 • ζ ∅ ik -L ∅ i • x ∅ i + k∈N ∂ i, x ∅ k -x ∅ i 2 • χ ∅ ik + x ∅ i • ξ ∅ i .
By substituting property FEM-b in the above identity, one finally obtain property FEM-c as

3L ∅ i = k∈N i, x ∅ k + x ∅ i 2 • η ∅ ik - x ∅ k -x ∅ i 2 • ζ ∅ ik + k∈N ∂ i, x ∅ k -x ∅ i 2 • χ ∅ ik . (12)
Considering now FV metric quantities, from the fact that n ∅ i = -n ∅ k over ∂C ik , property FVM-a reads ν ∅ ik = -ν ∅ ki , which corresponds to the conservation property of the scheme. From the Gauss theorem, one also has

C i ∇ ∅ R = ∂C i Rn ∅ i ,
which, from the definition of FV metric quantities, gives property FVM-b as

V ∅ i = k∈N i, ν ∅ ik + ν ∅ i . (13)
Property FVM-c is obtained by noting that

3V ∅ i = C i R∇ ∅ • x ∅ = ∂C i Rx ∅ • n ∅ i - C i R • x ∅ . ( 14
)
The right hand side of ( 14) is now computed by means of the FV discretization described in subsection 3.2 as

3V ∅ i = k∈N i, x ∅ k + x ∅ i 2 • ν ∅ ik -x ∅ i • L ∅ i + x ∅ i • ν ∅ i .
which from property FVM-b becomes

3V ∅ i = k∈N i, x ∅ k + x ∅ i 2 • ν ∅ ik . (15)
Therefore, a FV approximation can be formally obtained from FE metric quantities defined over the same grid points by setting (see properties FEM/FVM-a and -b)

ν ∅ ik = η ∅ ik , ν ∅ i = ξ ∅ i , V ∅ i = L ∅ i .
Note that the mass lumping approximation,

k∈N i M ∅ ik du k dt ≃ L ∅ i du i dt
has to be introduced in the Eq. ( 7) for the equivalence conditions to be applicable. By subtracting Eq. ( 12) to Eq. ( 15), one finally has

V ∅ i = L ∅ i + k∈N i, x ∅ k -x ∅ i 6 • ζ ∅ ik - k∈N ∂ i, x ∅ k -x ∅ i 6 • χ ∅ ik . (16)
It is remarkable that, differently from the Cartesian case [START_REF] Selmin | The node-centred finite volume approach: bridge between finite differences and finite elements[END_REF][START_REF] Selmin | Unified construction of finite element and finite volume discretizations for compressible flows[END_REF], in the cylindrical reference the FV cell is not coincident with the FE lumped mass matrix. Moreover, the shape of the FV cells that guarantees equivalence with FE discretization still remains to be determined.

Finite Volume/Element method in spherical coordinates

Node-pair finite element discretization

In the present section the spatial discretization of conservation law in a spherical reference system is presented.

The scalar conservation law (3) is multiplied by the quantity r sin θ to formally remove the singularity of the spherical coordinate system along the zenith direction. The FE node-pair discretization is obtained as done in the cylindrical case [START_REF] Santis | Finite volueme and finite element schemes for the euler equations in spherical coordinates[END_REF], namely

M ⊙ ik du i dt = - k∈N i , f ⊙ k + f ⊙ i 2 • η ⊙ ik - f ⊙ k -f ⊙ i 2 • ζ ⊙ ik + f ⊙ i • L ⊙ i - k∈N ∂ i, f ⊙ k -f ⊙ i 2 • χ ⊙ ik -f ⊙ i • ξ ⊙ i , (17)
where the following metric quantities have been introduced

M ⊙ ik def = Ω ik r sin θϕ i ϕ k dΩ ⊙ , η ⊙ ik def = Ω ik r sin θ(ϕ i ∇ ⊙ ϕ k -ϕ k ∇ ⊙ ϕ i ) dΩ ⊙ , ζ ⊙ ik def = Ω ik ϕ i ϕ k ∇ ⊙ (r sin θ) dΩ ⊙ , L ⊙ i def = k∈N i Ω ik ϕ i ϕ k ∇ ⊙ (r sin θ) dΩ ⊙ , χ ⊙ ik def = ∂Ω ∂ ik r sin θϕ i ϕ k n ⊙ d∂Ω ⊙ , ξ ⊙ i def = ∂Ω ∂ ik r sin θϕ i n ⊙ d∂Ω ⊙ .

Edge-based finite volume discretization

The FV spatially discrete form of the scalar conservation law (3) multiplied by the quantity r sin θ reads

V ⊙ i du i dt = - k∈N i, f ⊙ i + f ⊙ k 2 • ν ⊙ ik + f ⊙ i • V ⊙ i -f ⊙ i • ν ⊙ i ( 18
)
where

V ⊙ i def = ∂C ∂ i ∇ ⊙ (r sin θ), ν ⊙ ik def = ∂C ik r sin θn ⊙ i and ν ⊙ i def = ∂C ∂ i r sin θn ⊙ i 4.

Finite Element/Volume equivalence

The equivalence conditions relating the FV and FE scheme are obtained like in the cylindrical case. The following properties are introduced for the FE metric quantities

(FEM -a) η ⊙ ik = -η ⊙ ki (19) (FEM -b) L ⊙ i = k∈N i, η ⊙ ik + ξ ⊙ i (FEM -c) 3L ⊙ i = ∂Ω ∂ i r sin θϕ i x ⊙ • n ⊙ i - Ω i ϕ i x ⊙ • r - Ω i r sin θx ⊙ • ∇ ⊙ ϕ i
and for the FV metric quantities

(FVM -a) n ⊙ i = -n ⊙ k (20) (FVM -b) V ⊙ i = k∈N i, ν ⊙ ik + ν ⊙ i (FVM -c) 3V ⊙ i = k∈N i, x ⊙ i + x ⊙ k 2 • ν ⊙ ik
Therefore, a FV approximation can be formally obtained from FE metric quantities defined over the same grid points by setting (see properties FEM/FVM-a and -b)

ν ⊙ ik = η ⊙ ik , ν ⊙ i = ξ ⊙ i , V ⊙ i = L ⊙ i .
By subtracting (FEM-c) to (FVM-c), one finally has

V ⊙ i = L ⊙ i + k∈N i, x ⊙ k -x ⊙ i 6 • ζ ⊙ ik - k∈N ∂ i, x ⊙ k -x ⊙ i 6 • χ ⊙ ik , (21)
where, like the cylindrical case, the mass lumping approximation has been introduced. It is clear that also in this case he FV cell is not coincident with the FE lumped mass matrix and the shape of the FV cells that guarantees equivalence with FE discretization still remains undetermined.

Fully discrete form of the Euler equations in cylindrical and spherical coordinates

The Euler equations in cylindrical and spherical coordinates for compressible inviscid flows are now briefly recalled. The differential form reads T momentum and E t total energy per unit volume. The flux function and the source term are defined as follow in the cylindrical case where Π is the pressure function in terms of the conservative variables.

∂u ∅ ∂t + ∇ ∅ • f ∅ = 1 R s ∅ and ∂u ⊙ ∂t + ∇ ⊙ • f ⊙ = 1 r sin θ s ⊙ , respectively, where u ∅ (Z, R, θ, t) = (ρ, m ∅ , E t ) T and u ⊙ (r, θ, φ, t) = (ρ, m ⊙ , E t ) T , with ρ density, m ∅ = (m Z , m R , m θ ) T , m ⊙ = (m r , m θ , m φ )
The FV spatially-discrete form of the Euler equations in the cylindrical and spherical refrences reads

V ∅ i du ∅ i dt = - k∈N i , f ∅ k + f ∅ i 2 • η ∅ ik + L i s ∅ i -f ∅ i • ξ ∅ i and (22)
V ⊙ i du ⊙ i dt = - k∈N i, f ⊙ i + f ⊙ k 2 • η ⊙ ik + f ⊙ i • L ⊙ i -f ⊙ i • ξ ⊙ i + s⊙ i Li + s⊙ i Ľi ,
respectively, where the following notation has been used

s ∅ i def = s ∅ (u ∅ i ), s⊙ i def = s⊙ (u ⊙ i ), s⊙ i def = s⊙ (u ⊙ i ).
The terms V ∅ i and V ⊙ i are computed from the equivalence condition in cylindrical and spherical case, respectively. In the computation, a TVD [START_REF] Harten | Self adjusting grid methods for one-dimensional hyperbolic conservation laws[END_REF] numerical flux is used, with the van Leer limiter [START_REF] Van Leer | Towards the ultimate conservative difference scheme II. Monotoniticy and conservation combined in a second order scheme[END_REF]. The fully discrete form of the Euler system is obtained by a two-step Backward Differencing Formulae. At each time level, a dual time-stepping technique is used to solve the time-implicit problem [START_REF] Venkatakrishnan | Implicit method for the computation of unsteady flows on unstructured grids[END_REF]. In all simulation the ideal gas model for the nitrogen (γ = c p /c v = 1.39) is used.

Numerical results

In the present section, numerical results for converging and diverging shock waves are presented for both the cylindrical and spherical cases.

The first problem consider is a one dimensional simulation of converging shocks in cylindrical and spherical symmetry. This kind of problems can be seen as a special Riemann problem in cylindrical and spherical coordinates for which Guderley [START_REF] Guderley | Powerful spherical and cylindrical compression shocks in the neighbourhood of the centre and of the cylinder axis[END_REF] has showed that the position of the converging shock has a similar solution as follows

R R c = A ± 1 - t t c α ,
where R = R(t) is the position of the shock as function of the time t, R c is the position of the discontinuity at initial time and t c is the time value at which the discontinuity reaches the origin of the domain. The rational number α is the exponent of the similar law and A + and A -are the coefficients for the converging and diverging shock, respectively. All these coefficients can be computed in the case of the polytropic ideal gas [START_REF] Guderley | Powerful spherical and cylindrical compression shocks in the neighbourhood of the centre and of the cylinder axis[END_REF]. The numerical simulation are performed on a one dimensional domain with wall boundary conditions at the two boundaries. The initial solution consists of a uniform distribution of the density, the velocity is zero everywhere, while the pressure has a discontinuity at the center of the domain with a value of the pressure on the outer part of the domain ten times higher that the value on the inner part.

In figure 1 the Guderley's law is compared against the numerical results on a grid of 501 (∆r = 2 × 10 -3 ) nodes with a time step ∆t = 10 -4 for both the cylindrical and spherical problems. The numerical and the analytical solutions agree very well, even near the origin of the domain where there are strong shocks. After the shocks reflections, there are complex interactions between the waves coming to the origin and those reflected that are not taken into account by the analytical model. The explosion and the implosion problem are now considered in the case of two spatial dimension. The computational domain is shown in figure 2, where a representative computational grid is also shown. Initial conditions for the explosion problem are as follows. The velocity is assumed to be zero everywhere; the density is uniform and equal to 1, whereas the pressure is uniform and equal to 10 in a circular region centered at the origin with radius r = 0.5. In the remaining portion of the domain, the pressure is uniform and equal to 1.

The thermodynamics variables are made dimensionless by the corresponding value at rest (outer region), the spatial dimensions are made dimensionless by the radius of the domain. Velocity and time are made dimensionless by the square root of the reference pressure divided by the reference density; the reference time is the reference length divided by the reference velocity.

In figure 3 and 4 are reported the numerical solution respectively for the cylindrical coordinates system on the plane R-θ and spherical coordinates on the plane r-φ, for the explosion problem. The density isolines at different time levels are shown together with the pressure profiles along the axis y = 0. The grid is made with 39 153 nodes 77 587 triangles and the time step is 2.5 × 10 4 . A shock wave propagates towards the outer boundary of the computational domain; the shock wave is followed by a contact discontinuity. A rarefaction wave propagates towards the origin and is then reflected outward. Note that the initial corrugation of the shock front, due to the un-even shape of the initial discontinuity caused by its discrete representation over an unstructured grid of triangles, is clearly visible also at later times. The pressure profiles are compared against reference one-dimensional results for the three different time levels. One-dimensional computations were performed over a evenly-spaced grid made of 2 001 nodes, which corresponds to an element spacing of 5 × 10 4 .

A grid dependence study is shown in figure 5. Pressure signals along the y = 0 axis are compared at time t = 0.16 for three different grid resolutions: the coarse grid is made of 9 551 nodes and 18 745 triangles, the medium one is made of 20 683 nodes and 40 841 triangles, the fine one is made of 39 153 nodes 77 587 triangles. Numerical results are found to be almost independent from the grid resolution. Time step dependence can be appreciated from figure 6, where the pressure signals for three different time steps is shown at time t = 0.16 for the medium grid. Numerical results are found to be independent from the chosen time step.

In figure 7 and figure 8, numerical results for the implosion problem are showed in cylindrical and spherical case respectively. In the implosion problem the initial conditions are as in the explosion problem where now the high pressure region is the outer one and the low pressure region is at r < 0.5. The grid is the fine grid and the time step is 2.5 × 10 -4 . A rarefaction wave propagates towards the outer boundary; a shock wave and a contact surface propagates inwards. The intensity of the shock increases as it moves towards the origin; when the shock wave is reflected at the origin, a region of high pressure/temperature is observed. In the spherical case, this effect is more evident. Due to the symmetry of the solution and of the computational domain, the spherical implosion and explosion problems can be also simulated by axially symmetric Euler equations formulated in a cylindrical coordinate system, as done in [START_REF] Santis | Equivalence conditions for finite volume/element discretizations in cylindrical coordinates[END_REF]. The spherical and axisymmetric solution correspond to the solution of the same problem on two different planes, i.e., the r-φ plane in the spherical problem and the Z-R plane in the axisymmetric problem, where Z is the coordinate along the axis of symmetry and R is the cooridate along the axis normal to the axis of symmetry. The solutions in spherical coordinates are shown together with the corresponding solutions in axisymmetric coordinates in figure 9 for the implosion problem. In order to reconstruct the complete spherical solution the spherical and axisymmetric solutions are represented on each own plane.

Conclusions

A novel unstructured-grid hybrid finite element/volume method in a cylindrical and spherical reference was presented. The proposed approach represents an extension of the node-pair scheme earlier for the cartesian coordinates and moves from suitable equivalence conditions linking finite element integrals to the corresponding finite volume metrics, such as the cell volume or the integrated normals. The equivalence conditions, as done for the first time in the case for the cylindrical coordinates [START_REF] Santis | Equivalence conditions for finite volume/element discretizations in cylindrical coordinates[END_REF], were derived here without introducing any approximation and allowed to determine all needed finite volume metric quantities from finite element ones. Numerical results are presented for one and two spatial dimension compressible flows: these consist in the numerical simulation of the explosion and im- plosion problems, in which an initial discontinuity in pressure results in the formation of a diverging and converging shock, respectively. The computed pressure and density profile agree fairly well with one-dimensional simulation in cylindrical and spherical symmetry over a very fine grid. The solutions obtained in spherical reference also agree very well with the corresponding solutions in a cylindrical reference where the axisymmetric condition has been used to simulate spherical explosions and implosions.

Appendix A. Node-pair finite element for a scalar conservation law

In this appendix, the following two identities are demonstrated, namely

k∈N i f ∅ k • Ω ik Rϕ k ∇ ∅ ϕ i = k∈N i , f ∅ k + f ∅ i 2 • η ∅ ik - f ∅ k -f ∅ i 2 • ζ ∅ ik + k∈N i , f ∅ k -f ∅ i 2 • χ ∅ ik (A.1)
and which together allow to recast the discrete Bubnov-Galekin equation ( 6) in its node-pair counterpart [START_REF] Santis | Equivalence conditions for finite volume/element discretizations in cylindrical coordinates[END_REF].

k∈N i f ∅ k • ∂Ω ∂ ik Rϕ i ϕ k n ∅ = k∈N ∂ i, f ∅ k -f ∅ i • χ ∅ ik -f ∅ i • ξ ∅ i , (A.2)
The proof of the identity (A.1) is considered first. The integral on the left hand side of the Eq. (A.1) is assembled considering the contributions coming from each element e in the mesh, exploiting the local support property of the shape functions, as follows

k∈N i f ∅ k • Ω ik Rϕ k ∇ ∅ ϕ i = e∈E i k∈N e f ∅ k • Ω e Rϕ k ∇ ∅ ϕ i , (A.3)
where E i is the set of the elements having the node i in common and N e is the set of the nodes of element e. The first summation on the right-hand side is limited to the elements contained in the support Ω i of node i , which are the only ones to give a nonzero contribution to integrals containing the function ϕ i . Note that Ω i = e∈E i Ω e . Considering now the following identity from the Gauss theorem

Ω e ∇ ∅ (Rϕ i ϕ k ) dΩ ∅ = ∂Ω e Rϕ i ϕ k n ∅ d∂Ω ∅ ,
which allows to write

Ω e Rϕ k ∇ ∅ ϕ i = - Ω e Rϕ i ∇ ∅ ϕ k - Ω e ϕ i ϕ k R + ∂Ω e Rϕ i ϕ k n ∅ . (A.4)
where R = ∇ ∅ R is the unit vector in the radial direction. Thanks to the previous relation one deduces

Ω e Rϕ k ∇ ∅ ϕ i = 1 2 Ω e Rϕ k ∇ ∅ ϕ i + 1 2 Ω e Rϕ k ∇ ∅ ϕ i = 1 2 Ω e Rϕ k ∇ ∅ ϕ i + 1 2 - Ω e Rϕ i ∇ ∅ ϕ k - Ω e ϕ i ϕ k R + ∂Ω e Rϕ i ϕ k n ∅ = - 1 2 η e,∅ ik - 1 2 Ω e ϕ i ϕ k R + 1 2 ∂Ω e Rϕ i ϕ k n ∅ , (A.5)
where in the last equality has been introduced the elemental contributions η e,∅ ik of the element e to the vector η

∅ ik η e,∅ ik def = Ω ik ∩Ω e R(ϕ i ∇ ∅ ϕ k -ϕ k ∇ ∅ ϕ i ), such that η ∅ ik = e∈(E i ∩E k ) η e,∅ ik
. By the relation (A.5), the integral (A.3) becomes

k∈N i f ∅ k • Ω ik Rϕ k ∇ ∅ ϕ i = - e∈E i k∈N e f ∅ k • 1 2 η e,∅ ik + 1 2 Ω e ϕ i ϕ k - 1 2 ∂Ω e Rϕ i ϕ k n ∅ . (A.6)
On the other hand, from the Eq. (A.5) follows also that

η e,∅ ik = -2 Ω e Rϕ k ∇ ∅ ϕ i - Ω e ϕ i ϕ k R + ∂Ω e Rϕ i ϕ k n ∅ ,
which can be recast, using the equation (A.4), as

η e,∅ ik = 2 Ω e Rϕ i ∇ ∅ ϕ k + Ω e ϕ i ϕ k R - ∂Ω e Rϕ i ϕ k n ∅ .
Summing up the last relation for all the nodes k belonging to the element Ω e and using the fact that k∈N e ∇ ∅ ϕ k (x ∅ ) = 0, ∀x ∅ ∈ Ω e , one obtains k∈N e η e,∅ ik -

Ω e ϕ i ϕ k R + ∂Ω e Rϕ i ϕ k n ∅ = 0. (A.7)
Summing up the previous relation for all the elements belonging to E i and multiplying by the vector f ∅ i , follows that

e∈E i k∈N e f ∅ i • η e,∅ ik - Ω e ϕ i ϕ k R + ∂Ω e Rϕ i ϕ k n ∅ = 0.
Multiplying this relation by 1/2 and adding it to the right hand side of (A.6) one has

k∈N i f ∅ k • Ω ik Rϕ k ∇ ∅ ϕ i = - e∈E i k∈N e f ∅ k + f ∅ i 2 • η e,∅ ik + f ∅ k -f ∅ i 2 • Ω e ϕ i ϕ k R + e∈E i k∈N e f ∅ k -f ∅ i 2 • ∂Ω e Rϕ i ϕ k n ∅ .
By recalling that η e,∅ ik = 0 for e (E i ∩ E k ), that η ∅ ik = e∈(E i ∩E k ) η e,∅ ik and that η ∅ ii = 0, the right hand side of the last equation can be written as

k∈N i f ∅ k • Ω ik Rϕ k ∇ ∅ ϕ i = - k∈N i, f ∅ k + f ∅ i 2 • η ∅ ik + f ∅ k -f ∅ i 2 • ζ ∅ ik + k∈N ∂ i, f ∅ k -f ∅ i 2 • ∂Ω ∂ ik Rϕ i ϕ k n ∅ ,
that is the relation (A.1).

Considering now the proof of the identity (A.2), in the left hand side of (A.2) the contribute of the node i is put into evidence, namely

k∈N ∂ i f ∅ k • ∂Ω ∂ ik Rϕ i ϕ k n ∅ = k∈N ∂ i, f ∅ k • ∂Ω ∂ ik Rϕ i ϕ k n ∅ + f ∅ i • ∂Ω ∂ ik Rϕ i ϕ i n ∅ .
The quantity

k∈N ∂ i, f ∅ i • ∂Ω ∂ ik Rϕ i ϕ k n ∅ ,
is now added and subtracted from the right hand side to obtain

k∈N ∂ i f ∅ k • ∂Ω ∂ ik Rϕ i ϕ k n ∅ = k∈N ∂ i, ( f ∅ k -f ∅ i ) • ∂Ω ∂ ik Rϕ i ϕ k n ∅ + f ∅ i • k∈N ∂ i ∂Ω ∂ ik Rϕ i ϕ k n ∅ .
By recalling that k∈N e ϕ k (x ∅ ) = 1, x ∅ ∈ Ω e , ∀e ∈ E, one has

k∈N ∂ i f ∅ k • ∂Ω ∂ ik Rϕ i ϕ k n ∅ = k∈N ∂ i, ( f ∅ k -f ∅ i ) • ∂Ω ∂ ik Rϕ i ϕ k n ∅ + f ∅ i • ∂Ω ∂ i Rϕ i n ∅ = k∈N ∂ i, ( f ∅ k -f ∅ i ) • χ ∅ ik + f ∅ i • ξ ∅ i ,
which is the relation (A.2). Using the relations (A.1), (A.2) is possible to write (6) as follow

k∈N i M ∅ ik du k dt = - k∈N i , f ∅ k + f ∅ i 2 • η ∅ ik + f ∅ k -f ∅ i 2 • ζ ∅ ik + k∈N i f ∅ k • Ω ik ϕ i ϕ k R - k∈N ∂ i, f ∅ k -f ∅ i 2 • χ ∅ ik -f ∅ i • ξ ∅ i , (A.8)
with the metric quantities defined in the subsection 4.1. In the previous relation is necessary to write in node-pair form also the second term on the right hand side. By recalling that

k∈N i f ∅ k • Ω ik ϕ i ϕ k R = k∈N i, f ∅ k • Ω ik ϕ i ϕ k R + f ∅ i • Ω ik ϕ i ϕ k R,
by adding and subtracting to the previous relation the following quantity

k∈N i, f ∅ i • Ω ik ϕ i ϕ k R, one has k∈N i f ∅ k • Ω ik ϕ i ϕ k R = k∈N i, ( f ∅ k -f ∅ i ) • Ω ik ϕ i ϕ k R + k∈N i f ∅ i • Ω ik ϕ i ϕ k R = k∈N i, ( f ∅ k -f ∅ i ) • Ω ik ϕ i ϕ k R + f ∅ i • Ω ik ϕ i R = k∈N i, ( f ∅ k -f ∅ i ) • Ω ik ϕ i ϕ k R + f ∅ i • L ∅ i .
(A.9) by substituting this relation into the (A.8) one has the node-pair FE discretization of the scalar equation ( 2)

L ∅ i du i dt = - k∈N i , f ∅ k + f ∅ i 2 • η ∅ ik - f ∅ k -f ∅ i 2 • ζ ∅ ik + f ∅ i • L ∅ i - k∈N ∂ i, f ∅ k -f ∅ i 2 • χ ∅ ik -f ∅ i • ξ ∅ i , (A.10)
where the mass matrix has been lumped

k∈N i M ∅ ik du i dt ≃ L ∅ i du i dt , with L ∅ i = k∈N i M ∅ ik .

Figure 1 :

 1 Figure 1: Position of the converging shock as function of the time for the one-dimensional simulation and for the analytical solution of Guderley: (a) cylindrical problem, (b) spherical problem.

Figure 2 :

 2 Figure 2: Exemplary grid for the explosion and implosion problems. The grid is the coarse grid.

Figure 3 :

 3 Figure 3: Density isoline for the cylindrical implosion problem: (a) t = 0.02; (b) t = 0.12; (c) t = 0.16. Each isoline corresponds to a density difference of ∆ρ/ρ re f = 0.03. (d) Pressure signal along the y = 0 axis at the same time levels: the solid line is the reference one-dimensional solutions, the dot-dashed line is the bidimensional solution.

Figure 4 :

 4 Figure 4: Density isoline for the spherical implosion problem: (a) t = 0.02; (b) t = 0.12; (c) t = 0.16. Each isoline corresponds to a density difference of ∆ρ/ρ re f = 0.03. (d) Pressure signal along the y = 0 axis at the same time levels: the solid line is the reference one-dimensional solutions, the dot-dashed line is the bidimensional solution.

Figure 5 : 4 1DFigure 6 :

 546 Figure 5: Comparison of pressure signals along the y = 0 axis for the explosion problem at time t = 0.16 for different grid resolutions against one dimensional simulations. Left: cylindrical case, right: spherical case.

Figure 7 :

 7 Figure 7: Density isoline for the cylindrical implosion problem: (a) t = 0.02; (b) t = 0.12; (c) t = 0.16. Each isoline corresponds to a density difference of ∆ρ/ρ re f = 0.03. (d) Pressure signal along the y = 0 axis at the same time levels: the solid line is the reference one-dimensional solutions, the dot-dashed line is the bidimensional solution.

Figure 8 :

 8 Figure 8: Density isoline for the spherical implosion problem: (a) t = 0.02; (b) t = 0.12; (c) t = 0.16. Each isoline corresponds to a density difference of ∆ρ/ρ re f = 0.03. (d) Pressure signal along the y = 0 axis at the same time levels: the solid line is the reference one-dimensional solutions, the dot-dashed line is the bidimensional solution.

Figure 9 :

 9 Figure 9: Pressure (a), (c) and density (b), (d) contours for the implosion problem at two different time levels: t = 0.1 first row, t = 0.32 second row. On the horizontal plane is represented the solution of the bidimensional spherical problem on the r-φ plane, while on the vertical plane is represented the solution of the cylindrical axisymmetric problem on the Z-R plane.


Benoit Da Mota

Michael Eickenberg

Soizic Laguitton

Vincent Frouin

Gaël Varoquaux

Jean-Baptiste Poline
email: jbpoline@gmail.comvincent.frouin@cea.fr

Bertrand Thirion
email: bertrand.thirion@inria.frsoizic.laguitton

A MapReduce Approach for Ridge Regression in Neuroimaging-Genetic Studies

Keywords: Bio-statistics, Neuroimaging, Genetics, Ridge regression, Permutation Tests

In order to understand the large between-subject variability observed in brain organization and assess factor risks of brain diseases, massive efforts have been made in the last few years to acquire highdimensional neuroimaging and genetic data on large cohorts of subjects. The statistical analysis of such high-dimensional and complex data is carried out with increasingly sophisticated techniques and represents a great computational challenge. To be fully exploited, the concurrent increase of computational power then requires designing new parallel algorithms. The MapReduce framework coupled with efficient algorithms permits to deliver a scalable analysis tool that deals with high-dimensional data and hundreds of permutations in a few hours. On a real functional MRI dataset, this tool shows promising results.

Introduction

Using genetics information in conjunction with brain imaging data is expected to significantly improve our understanding of both normal and pathological variability of brain organization. It should lead to the development of biomarkers and in the future personalized medicine. Among other important steps, this endeavor requires the development of adapted statistical methods to detect significant associations between the highly heterogeneous variables provided by genotyping and brain imaging, and the development of the software components that will permit large-scale computation to be done.

In current settings, neuroimaging-genetic datasets consist of a set of i) genotyping measurements at given genetic loci, such as Single Nucleotide Polymorphisms (SNPs) that represent a large amount of the genetic between-subject variability, on the one hand, and ii) quantitative measurements at given locations (voxels) in three-dimensional images, that represent e.g. either the amount of functional activation in response to a certain task or an anatomical feature, such as the density of grey matter in the corresponding brain region.

Most of the efforts so far have been focused on designing association models, and the computational procedures used to run these models on actual architectures have not been considered carefully. For instance, permutation tests of simple linear association models have been deemed as inefficient in some of these studies, e.g. [START_REF] Stein | and Alzheimer's Disease Neuroimaging Initiative. Voxelwise genome-wide association study (vGWAS)[END_REF]; however, they can be replaced by analytical tests only in very specific cases and under restrictive assumptions. Gains in sensitivity might be provided by multivariate models in which the joint variability of several genetic variables is considered simultaneously. Such models are thought to be more powerful [START_REF] Vounou | Discovering genetic associations with high-dimensional neuroimaging phenotypes: A sparse reduced-rank regression approach[END_REF][START_REF] Bunea | Penalized least squares regression methods and applications to neuroimaging[END_REF][START_REF] Kohannim | Boosting power to detect genetic associations in imaging using multi-locus, genome-wide scans and ridge regression[END_REF][START_REF] Meinshausen | Stability selection[END_REF], because they can express more complex relationships than simple pairwise association models. The cost of unitary fit is high due to high-dimensional, potentially non-smooth optimization problems and various cross-validation loops needed to optimize the parameters; moreover, permutation testing is necessary to assess the statistical significance of the results of such procedures in the absence of analytical tests. Multivariate statistical methods require thus many efforts to be tractable in this problem on both the algorithmic and implementation side, including the design of adapted dimension reduction schemes. In this work we will consider the simplest approach, ridge regression [START_REF] Kohannim | Boosting power to detect genetic associations in imaging using multi-locus, genome-wide scans and ridge regression[END_REF], that is powerful for detecting multivariate associations between large variable sets, but does not enforce sparsity in the solution.

Working in a distributed context is necessary to deal with the memory and computational loads, and yields specific optimization strategies. Once the unitary fit cost has been minimized, the main task when implementing such natural data parallel applications is to choose how to split the problem into smaller subproblems to minimize computation, memory consumption and communication overhead. For the first time, we propose an efficient framework that can manage ridge regression with numerous phenotypes and permutations.

In Section 2, we present our sequential algorithm, then we describe our framework to distribute efficiently the computation on large infrastructures. Experimental results on simulated and real data are presented in Section 3.

2 Methods : the computational framework Ridge regression of neuroimaging genetics data is clearly an embarrassingly parallel problem, which can be easily split into smaller tasks. Our computational framework relies on an adapted workflow summarized in Fig. 1, in which subtasks are optimized for the sake of efficiency. To simplify the presentation we first describe the core algorithm and then the workflow.

Optimizing the ridge regression algorithm

The map step, i.e. the scoring of ridge classifiers, is the most demanding in computation time (< 99.9% in our final implementation) and thus has to be optimized in priority. The computational burden mostly depends on the ridge regression step. Our algorithm performs Ridge Regression for multiple targets and multiple individual penalty values. It solves the following problem:

βij = argmin β y i -Xβ 2 2 + λ ij β 2 2 , i ∈ [1, p], j ∈ [1, J]
where X ∈ R n×p is the gene data matrix, y i ∈ R n is a variable extracted from brain images, βij ∈ R p is the estimated coefficient vector, and λ ij > 0 is the penalty term where j indexes J different penalties for the target y i . We obtain the solution using the singular value decomposition (SVD) of X, which we write X = U SV T , truncated to non-zero singular values. In the full rank case and for p > n we have

U ∈ R n×n and V T ∈ R n×p , while S is a diagonal matrix with entries s k , 1 ≤ k ≤ n. For one βij we have βij = V diag k s k s 2 k + λ ij U T y i
All βij are calculated with the same SVD, it is reused (and cached). For all i, U T y i is pre-calculated, which is conveniently and effectively done by multiplying the matrices U T and Y where the columns of Y are the y i . Since for a given j every target i potentially has a different penalty associated, the shrinkage operation s k s 2 k +λij is not writable as a matrix multiplication against U T Y . However, it is a linear operation on matrices, and by defining

Σ ∈ R n×p with Σ ki = s k s 2 k +λij
for a fixed j, it can be written as the pointwise matrix product

β = V (Σ • U T Y ).
These are the operations implemented by our algorithm for the J different sets of penalties, using J different matrices Σ. With a pre-calculated SVD and U T Y , the cost of this operation is O(npN ), where N is the number of target variables.

Care was taken of computational/hardware sources of optimization, like CPU cache issues. For instance, matrix-based operations are used instead of vectorbased operations to optimize the use of advanced vector extensions instructions set in new CPU. Our Python code uses the Numpy/Scipy/Scikit-learn scientific libraries, which rely on standard and optimized linear algebra libraries (Atlas or MKL) that are several order of magnitude faster than naive code. Next, we need to consider evaluation and parameter setting procedures:

the power of the procedure is measured by the ratio of explained variance, computed within a shuffle-split loop that leaves 20% of the data as a test set at each of the ten iterations; -to select the optimal shrinkage parameters, J = 5 values are tested first, then a grid refinement is performed where five other parameters are tested ; -each shrinkage parameter of the ridge regression is evaluated using an inner 5-folds cross validation loop.

This setting thus needs approximately 500 ridge regressions for one phenotype and one permutation.

The distributed algorithm

The MapReduce framework [START_REF] Chu | Map-reduce for machine learning on multicore[END_REF][START_REF] Dean | MapReduce: simplified data processing on large clusters[END_REF] seems the most natural approach to handle this problem and can easily harness large grids. The Map step yields explained variance for an image phenotype and for each permutation, while the reduce step consists in collecting all results to compute statistic distribution and corrected pvalues. Sub-tasks are created in a way that minimizes inputs/outputs (I/O). By essence, permutations imply computations on the same data after shuffling. The permutation procedure is thus embedded in the mapper, so that all permutations loops are run on the same node for a given dataset and the problem is split in the direction of the brain data. Figure 1 gives an overview of our framework.

Shared cache is a crucial feature since the Map step is dominated by costly SVDs. For instance, with 1,000 phenotypes and 1,000 permutations, each SVD in the inner CV loop is required 2 millions times and costs few tens of second. A shared cache on NFS, provided by the Joblib Python library [START_REF] Varoquaux | Joblib: running python function as pipeline jobs[END_REF], coupled with system cache saves many computations.

Results

We present three types of results. First, we present the performances of our distributed framework. Then, we illustrate the interest of our approach on simulated data with known ground truth. Finally, we present results on a real dataset.

Performance evaluation of the procedure

To illustrate the scalability of our Map-Reduce procedure, we execute the whole framework on a cluster of 20 nodes; each one is a 2 × Intel(R) Xeon(R) CPU X5650 (6 cores) @ 2.67GHz with 48GB of memory, connected with Gigabit Ethernet LAN; all files were written on the NFS storage file-system; our mapper runs with the Enthought Python Distribution (EPD 7.2-2-rh5 64 bits) with scikits-learn 0.11 [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF] and with the MKL as linear algebra library with OpenMP parallelization disabled; the workflow is described and submitted with the somaworkflow software [START_REF] Laguitton | Soma-workflow: a unified and simple interface to parallel computing resources[END_REF]. This framework makes it possible i) to describe a set of independent tasks that are executed following an execution graph and ii) to execute the code by submitting the graph to classical queuing systems operating on the cluster. We report in Fig. 2 the result of an execution with almost all the 240 cores available during all the run. The workflow is composed by 1,000 mappers and 1 reducer tasks. The mappers represent 99.9% of the total of serial computation time. Once the SVD are cached, the execution time of a map task is around 20 minutes. We can see in Fig. 2 that after 88 minutes, we use only few cores, but all the unused cores are available for other users. This comes from the number of tasks: on 240 cores, after 4 batches of 240 tasks, only 40 are left.

To improve the global speedup, we could split the problem into smaller pieces to decrease the time of the mappers or we could choose a more optimal splitting. We have not explored these possibilities yet.

Simulated Data

We simulate functional Magnetic Resonance Images (fMRI) from real genetic data obtained from the Imagen database [START_REF] Schumann | The IMAGEN study: reinforcement-related behaviour in normal brain function and psychopathology[END_REF]. We use the number of minor alleles for each SNP and we assume an additive genetic model. We use only the first chromosome in which ten random SNPs produce an effect in a spherical brain region, centered at a random position in the standard space, then intersected with the support of grey matter using a mask computed for the Imagen dataset (see below). We add i.i.d. Gaussian noise, smoothed spatially with a Gaussian kernel (σ = 3mm), to model other variability sources. The effect size and the Signal-to-Noise Ratio (SNR) can vary across simulations. Then 1,000 imaging phenotypes are obtained by computing the mean signal in brain parcels that were created using a Ward Agglomeration clustering. To assess our approach, ten different datasets were generated and were run on our framework with P=200 permutations to estimate the distribution of the maximum explained variance under the null hypothesis. Results are given in Table 1 and show that our method detects 8 effects among 10 simulations with a p-value p < .05. The results do not give evidence of the influence of the SNR simulation nor of the volume of the effect on the test sensitivity.

Results on a real dataset

We used data from Imagen, a large multi-centric and multi-modal neuroimaging database [START_REF] Schumann | The IMAGEN study: reinforcement-related behaviour in normal brain function and psychopathology[END_REF] containing functional magnetic resonance images (fMRI) associated with 99 different contrast images in more than 1,500 subjects. The dataset is built on the first batch of subjects of the study. Regarding the fMRI data, the protocol in [START_REF] Pollak | Early experience is associated with the development of categorical representations for facial expressions of emotion[END_REF] was used, which yields the [angry faces -neutral] functional contrast (i.e. the difference between watching angry faces or neutral faces).

Imaging phenotype. Standard preprocessing, including slice timing correction, spike and motion correction, temporal detrending (functional data), and spatial normalization (anatomical and functional data), were performed using the SPM8 software and its default parameters; functional images were resampled at 3mm resolution. Obvious outliers detected using simple rules such as large registration or segmentation errors or very large motion parameters were removed after this step. The [angry faces -neutral] contrast was obtained using a standard linear model, based on the convolution of the time course of the experimental conditions with the canonical hemodynamic response function, together with standard high-pass filtering procedure and temporally auto-regressive noise model. The estimation of the model parameters was carried out using the SPM8 software. A mask of the grey matter was built by averaging and thresholding the individual grey matter probability maps. Subjects with too many missing data (imaging or genetic) or not marked as good in the quality check were discarded. An outliers detection [START_REF] Fritsch | Detecting outlying subjects in high-dimensional neuroimaging datasets with regularized minimum covariance determinant[END_REF] was run and 10% of the most outlier subjects were eliminated.

Genotype. We keep only SNPs in the first chromosome with less than 2% missing data. All the remaining missing data were replaced by the median over the subjects for the corresponding variable. The age, the sex and the acquisition center were taken as confounding variables.

The final dataset contains 1,229 subjects, 1,000 brain parcels, 31,790 SNPs and 10 confounding variables. Our Map-Reduce framework was run with P=1,000 permutations to assess statistical significance. The workflow takes approximately 9 hours on the previously described 240 cores cluster, for a theoretical serial time around 75 days (i.e. a speed-up of approximately 200). Only one parcel is detected with a corrected p-value ≤ 0.05. A view of the location of the detected parcel is reported in Fig. 3.

Conclusion

Penalized linear models represent an important step in the detection of associations between brain image phenotypes and genetic data, which faces a dire sensitivity issue. Such approaches require cross validation loops to set the hyperparameters and for performance evaluation. Permutations have to be used to assess the statistical significance of the results, this yielding prohibitively expensive analyses. In this paper, we present an efficient and scalable framework that can deal with such a computational burden and that we used to provide a realistic assessment of the statistical power of our approach on simulations. Our results on simulated data highlight the potential of our method and we provide promising preliminary results on real data, including one multivariate association that reaches significance. To the best of our knowledge, this is the first result of that kind in a brain-wide chromosome-wide association study, although it needs to be reproduced to be considered as meaningful.

Fig. 1 .

 1 Fig. 1. Overview of the Map-reduce framework for the application of Ridge Regression in neuroimaging-genetics. Permutation ID is 0 for not permuted data.

Fig. 2 .

 2 Fig. 2. Setting and execution of the MapReduce algorithm on the cluster

Fig. 3 .

 3 Fig. 3. Location of the brain parcel with a significant explained variance ratio (exp. var. = 0.019, p ≃ 0.048, corrected for multiple comparisons) on the real dataset.

Table 1 .

 1 Results on the simulated datasets p-value the p-value corresponding to the given ratio of explained variance, obtained by 200 permutations)

	Simul.	Volume	Average	Best Parcel
	#	(mm 3 )	SNR	explained variance	p-value
	1	3375	0.19	0.022	0.005
	2	3348	0.66	0.042	0.005
	3	2457	0.30	0.056	0.005
	4	2754	0.54	0.033	0.005
	5	3213	0.22	0.007	0.35
	6	3348	1.50	0.027	0.005
	7	1431	0.55	0.031	0.005
	8	1890	0.66	0.005	0.5
	9	3375	0.19	0.036	0.005
	10	3132	0.41	0.026	0.005

Acknowledgment Support of this study was provided by the IMAGEN project, which receives research funding from the European Community's Sixth Framework Program (LSHM-CT-2007-037286) and coordinated project ADAMS (242-257) as well as the UK-NIHR-Biomedical Research Centre Mental Health, the MRC-Addiction Research Cluster Genomic Biomarkers, and the MRC program grant Developmental pathways into adolescent substance abuse (93558). This research was also supported by the German Ministry of Education and Research (BMBF grant # 01EV0711). This manuscript reflects only the author's views and the Community is not liable for any use that may be made of the information contained therein. This work was supported in part by the joint INRIA -Microsoft Research Center. The experiments presented in this paper were carried out using the computing facilities of the CATI.