Fabien Casenave

Laurent Métivier

Gwendoline Pajot-Métivier

Isabelle Panet

Fast computation of general forward gravitation problems

Keywords: Gravitation anomaly computation, general 3D geometry, Fast Multipole Method

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The problem of the forward computation of the gravity potential, its gradient, and its second derivatives, generated by a mass density distribution of general 3D geometry is a requisite for many geophysical applications. These applications range from global topographic correction or local terrain correction (the correction of gravitation measurements from the effect of the topography), to the gravity disturbance generated by any 3D object. Geophysicists and oil, gas and mining companies have been investigating the gravity field (and other techniques) for a long time to explore the subsurface for energy and mineral resources, [START_REF] Blakely | Potential theory in gravity and magnetic applications[END_REF]. Terrain corrections and inverse gravity problems are to be solved. In both cases, one needs a fast and reliable method for the forward computation of gravitation anomalies associated to a given mass distribution. For a complete a recent review on these methods, see [START_REF] Jekeli | Comparison of methods to model the gravitational gradients from topographic data bases[END_REF][START_REF] Nabighian | The historical development of the magnetic method in exploration[END_REF][START_REF] Nabighian | The historical development of the magnetic method in exploration[END_REF].

The Newton integral for computing the gravity potential at a point x for a mass distribution 1 with support included in a domain V ⊂ R 3 writes φ(x) := G y∈V ρ(y)

x -y dV (y), [START_REF] Balmino | Spherical harmonic modelling to ultrahigh degree of bouguer and isostatic anomalies[END_REF] where G denotes the gravitational constant, ρ(y) is the density distribution at y and dV (y) is the considered volume element. The first and second derivatives of φ are respectively

g i (x) := ∂φ ∂x i (x) = G y∈V ρ(y) y i -x i x -y 3 dV (y), 1 ≤ i ≤ 3, (2)
and

T i,j (x) := ∂ 2 φ ∂x i ∂x j (x) = G y∈V ρ(y) 3(x i -y i )(x j -y j )
x -y 5 -δ i,j x -y 3 dV (y), 1 ≤ i, j ≤ 3, (3) where δ i,j denotes the Kronecker delta. Note that the computation of magnetic anomalies lead to the same kind of integrals. Many methods have been proposed to compute [START_REF] Barnett | Theoretical modeling of the magnetic and gravitational fields of an arbitrarily shaped three-dimensional body[END_REF] and [START_REF] Bhattacharyya | Magnetic anomalies due to prism-shaped bodies with arbitrary polarization[END_REF]. Closed-form equations were first developed for simple sources: magnetic anomaly of a vertical right cylinder [START_REF] Singh | Magnetic anomaly due to a vertical right circular cylinder with arbitrary polarization[END_REF] and a prism [START_REF] Bhattacharyya | Magnetic anomalies due to prism-shaped bodies with arbitrary polarization[END_REF]; gravity anomaly of a prism [START_REF] Nagy | The gravitational attraction of a right rectangular prism[END_REF]. Terrain correction by computing the gravity anomaly of a finite-thickness horizontal plate was proposed in [START_REF] Plouff | Gravity and magnetic fields of polygonal prisms and application to magnetic terrain corrections[END_REF]. Computation where sources are modeled by stacked laminae were proposed in [START_REF] Talwani | Rapid computation of gravitational attraction of threedimensional bodies of arbitrary shape[END_REF] (gravity) and [START_REF] Talwani | Computation with the help of a digital computer of magnetic anomalies caused by bodies of arbitrary shape[END_REF] (magnetic). Polyhedral elements have been used to model general shaped 3D sources, see [START_REF] Furness | A physical approach to computing magnetic fields1[END_REF][START_REF] Holstein | Gravimagnetic similarity in anomaly formulas for uniform polyhedra[END_REF][START_REF] Holstein | Invariance in gravimagnetic anomaly formulas for uniform polyhedra[END_REF][START_REF] Okabe | Analytical expressions for gravity anomalies due to homogeneous polyhedral bodies and translations into magnetic anomalies[END_REF][START_REF] Pajot | Characterization, analysis and interpretation of gravity gradiometry data[END_REF][START_REF] Petrovic | Determination of the potential of homogeneous polyhedral bodies using line integrals[END_REF][START_REF] Pohánka | Optimum expression for computation of the gravity field of a homogeneous polyhedral body1[END_REF][START_REF] Singh | New method for fast computation of gravity and magnetic anomalies from arbitrary polyhedra[END_REF][START_REF] Tsoulis | On the singularities of the gravity field of a homogeneous polyhedral body[END_REF] and [START_REF] Barnett | Theoretical modeling of the magnetic and gravitational fields of an arbitrarily shaped three-dimensional body[END_REF][START_REF] Hansen | Simplified frequency-domain expressions for potential fields at arbitrary three-dimensional bodies[END_REF] for magnetic data of uniform magnetization. Finite length geological structures have also been approximated as 2D bodies with strike extension perpendicular to a polygonal cross section [START_REF] Cady | Calculation of gravity and magnetic anomalies of finite-length right polygonal prisms[END_REF][START_REF] Rassmussen | End corrections in potential field modeling*[END_REF][START_REF] Shuey | End corrections in magnetic profile interpretation[END_REF][START_REF] Talwani | Rapid gravity computations for twodimensional bodies with application to the mendocino submarine fracture zone[END_REF]. Fourier transforms have been widely used [START_REF] Bhattacharyya | Continuous spectrum of the total-magnetic-field anomaly due to a rectangular prismatic body[END_REF][START_REF] Bhattacharyya | Some general properties of potential fields in space and frequency domain: A review[END_REF][START_REF] Cooley | An algorithm for the machine calculation of complex fourier series[END_REF][START_REF] Forsberg | Gravity field terrain effect computations by FFT[END_REF][START_REF] Parker | The rapid calculation of potential anomalies[END_REF][START_REF] Pedersen | The gravity and magnetic fields from ellipsoidal bodies in the wavenumber domain*[END_REF][START_REF] Sideris | Rigorous gravimetric terrain modeling using Molodensky operator[END_REF][START_REF] Wieczorek | SHTOOLS -Tools for working with spherical harmonics[END_REF]. Prism-based methods lead to very time-consuming computations when the size of the sources and observation points data become large. Fast Fourier Transforms (FFT) enable efficient computations; however, they are limited to periodic sources and observation points have to lie on a plane. Tesseroid-based computations have also been implemented in [START_REF] Grombein | Optimized formulas for the gravitational field of a tesseroid[END_REF][START_REF] Heck | A comparison of the tesseroid, prism and point-mass approaches for mass reductions in gravity field modelling[END_REF][START_REF] Uieda | Optimal forward calculation method of the marussi tensor due to a geologic structure at goce height[END_REF] for spherical or ellipsoid geometries. Recent developments aim at increasing the speed of prism-based methods [START_REF] Davis | Rapid gravity and gravity gradiometry terrain corrections via an adaptive quadtree mesh discretization[END_REF]. The mesh has an adapted resolution using a quad-tree. The geometry is limited to quasi 2D objects for local terrain corrections and the mesh changes depending on the location of the observation point.

In this work, we propose an efficient way to compute (2) and (3) for general 3D mass distributions and general sets of observation points, using an accelerated prism-based method. For the sake on conciseness, we restrict the presentation to g. In Section 2, we consider two methods to compute (2), analytic and numerical, and investigate their accuracy, with a tetrahedral mesh to model the mass density distribution. In Section 3, we present the Fast Multipole Method (FMM) and explain how it can be used to accelerate the computation of long-range interactions. We present our algorithm in Section 4: depending on the close-or long-range nature of the involved interactions, the algorithm automatically switches between analytic integration formulae and numerical quadratic formulae, and relies on the Fast Multipole Method to drastically increase the computation of long-range interactions while guarantying a chosen accuracy. In Section 5, we give details on how the parameters of this method are chosen. Section 6 proposes various topographic correction applications: over the Himalaya in Section 6.2, over France in Section 6.3, and Earth's global correction is presented in Section 6.4. A computation of the gravity field generated by the comet Churyumov-Gerasimenko is proposed in Section 6.5.

Direct computation

For general mass distributions, there exists no analytic formula to compute [START_REF] Barnett | Theoretical modeling of the magnetic and gravitational fields of an arbitrarily shaped three-dimensional body[END_REF]. Usually, the geometry is first modeled as a mesh. In this work, we use the tetrahedron as volume element. If we denote by N the number of tetrahedrons of the mesh, and {V k } 1≤k≤N , the set of these tetrahedrons, the formula (2) writes

g i (x) = G N k=1 y∈V k ρ(y) y i -x i x -y 3 dV (y), 1 ≤ i ≤ 3. (4)
In what follows, we discuss the integration over one volume element:

g i (x) V k = G y∈V k ρ(y) y i -x i x -y 3 dV (y), 1 ≤ i ≤ 3. ( 5
)

Analytical integration

We consider in this work the analytical formulae derived in [START_REF] Okabe | Analytical expressions for gravity anomalies due to homogeneous polyhedral bodies and translations into magnetic anomalies[END_REF]. They enable to compute the gravity potential, and its first and second derivatives, generated by homogeneous polyhedral bodies. These formulae are derived by applying the divergence theorem to the Newton integrals ( 2) to obtain surface integrals on the faces of the polyhedrons. Then, the divergence theorem is applied a second time to transform the surface integrals into 1D integrals over the edges of each face. To avoid introducing new notations, we refer to [START_REF] Okabe | Analytical expressions for gravity anomalies due to homogeneous polyhedral bodies and translations into magnetic anomalies[END_REF][Equation ( 31) and (51)] for the analytic formulae of respectively the first and second derivatives of the gravity potential. We denote the formula for computing g the "Okabe formula".

In [START_REF] Li | Three-dimensional gravity modeling in all space[END_REF][Chapter 5], it is specified that such analytic formulae are sensitive to round-off errors (errors occurring when using computer finite-precision arithmetics) when the distance from the target increases. Actually, since the signal deceases with the distance, there exists a distance beyond which the signal is dominated by rounding errors, as we will see in Section 2.3. Hence, we can evaluate the integral (5) accurately when x is not too far from the tetrahedron V k and with a uniform density in V k .

Numerical integration

Numerical integration schemes can be used to evaluate the integral (5):

g i (x) V k ≈ G S s=1 ω V k ,s ρ(y V k ,s ) (y V k ,s ) i -x i x -y V k ,s 3 v(V k ), 1 ≤ i ≤ 3, ( 6
)
where S is the order of quadrature, and {ω V k ,s } 1≤s≤S , {y V k ,s } 1≤s≤S are the sets of quadrature weights and points associated to the tetrahedron V k , and v(V k ) is the volume of V k . We consider quadrature formulae of order 1, 4, 5, 11, and 15 on tetrahedrons, which integrate exactly monomials respectively up to order 1, 2, 3, 4, and 5. We denote these formulae to compute g the quadrature formulae.

In general, error estimates using numerical integration scheme exist for regular enough functions. For instance, from [START_REF] Stoer | Introduction to Numerical Analysis[END_REF][Theorem 3.6.24], the error of order S 1D Gaussian quadrature is controlled by the 2S-th derivative of the considered function. Therefore, the singularity of (5) for y = x prevents the use of numerical integration where the observation points are too close to the tetrahedron V k . Moreover, as we will see in Section 2.3, numerical integration are more accurate as the distance between the source and the observation point increases.

Accuracy of the integration

It is well-known that Okabe formula are accurate for close-range interactions and inaccurate for long-range interactions, and that the opposite behavior holds for quadrature formula. In this section, we numerically investigate these assertions. Quantifying these behaviors is required to know when to switch between these two formulae in the proposed algorithm.

Consider a trirectangular tetrahedron, which legs (the three edge that meet at the right angle) have a length of 0.1, and consider the four meshes of this tetrahedron represented in Figure 1. Define {x i,r } 1≤i≤144 , a set of 144 points located on a sphere of radius r centered at the barycenter of the tetrahedron. We evaluate the first derivative of the gravity potential g using the Okabe and quadrature formulae on the sets {x i,r } 1≤i≤144 for some values of r, and denote them respectively g O (x i,r ) and g q (x i,r ). We then define the relative difference between the Okabe and quadrature formulae by

E r := i g q (x i,r ) -g O (x i,r ) 2 2 i g O (x i,r ) 2 2 , (7)
where • 2 denotes the Euclidian norm in R 3 .

Figure 2 shows the relative difference E r for various values of r and using a 1-point and a 15-points quadrature formula. The quadrature values are computed on Meshes 1 to 4, whereas the Okabe values are always computed in Mesh 1. For the sake of the presentation, we do not provide the results for the 5-points and the 11-points quadrature formula: they are simply the expected intermediate scenarios. The relative difference is globally much smaller for the 15-point quadrature formula than for the 1-point one. We notice that for small values of r, the difference is important, then decreases down to a minimum, and finally increases for large values of r. For small values of r, the Okabe formula is known to be accurate, but the quadrature formula fails due to the singularity of (5) at x = y. Notice that as the number of the mesh increases, the ratio r/c decreases, where c denotes the mean edge of the mesh, and the relative error decreases. It is due to the fact that since the integrated function is homogeneous, the error made by the quadrature does not depend on the absolute distance r to the tetrahedron, but on the distance ratio r/c. For large values of r, the contrary happens: the integrated function is regular enough for quadrature formulae to be accurate, whereas the Okabe formula suffers from important round-off errors.

In this work, we propose to separate, for each observation point, the tetrahedrons in two groups: the close-range and the long-range ones, and apply an Okabe formula for the first group and a quadrature formula for the second group. This way, we are able to ensure a given accuracy for the computed result. In Section 3, we present the Fast Multipole Method, an

The Fast Multipole Method (FMM)

The FMM is an algorithm that drastically speeds-up the computation of particular summations, to the price of a controlled approximation. In this section, we give some elements on the functioning of the FMM and we refer to the founding article [START_REF] Greengard | A fast algorithm for particle simulations[END_REF] for a complete description. In this funding article, there were already motivations for the fast computation of the N -body gravitation problem, where point masses where considered. In this work, the FMM will be applied to quadrature formulae that model volume integrals for long-range interactions only, and not on mass points models.

Let us denote by {x i } 1≤i≤M a set of observation points and by {y j } 1≤j≤N a set of source points, and consider the following summation:

f (x i ) := M j=1 σ j K(x i , y j ), 1 ≤ i ≤ N, ( 8
)
where K is called the kernel. In an ideal case, we suppose that there exist d functions g k (x) and

h k (y), 1 ≤ k ≤ d, such that K(x i , y j ) = d k=1 g k (x i )h k (y j ), (9)
then the following factorization holds

f (x i ) = d k=1 g k (x i )   M j=1 σ j h k (y j )   , 1 ≤ i ≤ N. (10)
The quantities Q k := M j=1 σ j h k (y j ) being independent from x i , they can be precomputed, so that the summation (8) can be computed in complexity O(max(N, M )), leading to important execution time gains.

The FMM deals with cases where the kernel K is not separable and has a singularity at x = y. Under certain assumptions, the asymptotic complexity O(max(N, M )) can be reached as well. Suppose that the domain is decomposed in an octree of given height. An octree is a recursive subdivision of the space of the following structure: first, the level zero, called the root, consists in the smallest cube containing the whole domain. Then, the first level is obtained by subdividing this cube in 8 cubes called boxes: 2 per direction. The following levels are constructed recursively in the same fashion, up to the last level, where the boxes are called the leaves, see Figure 5. The points (source and observation) have been added in the octree, meaning that each points lies in one box per level, and that the FMM algorithm knows the list of points lying in each box at each level. Figure 4 shows an example of tree in 1D. Definition 3.1 (Two kinds of "close"-and "long"-range interactions) At this points, it is mandatory to distinguish between the interactions considered "close-range" in the "approximation sense" and in the "octree sense". An interaction is called close-range in the approximation sense if is it well-approximated by an Okabe formula, see Section 2.3. On the contrary, an interaction is called close-range in the octree sense is the involved source and observation points lie in the same leaf (box at the last level) of the tree or in adjacent leaves, see Figure 4. This distinction is of paramount importance: for an octree of low height, the leaves are very large, and an approximation that is close-range in the octree sense has no guaranty to be close-range in the approximation sense.

Then, we need to obtain the separation of variables [START_REF] Casenave | A volumic mesh for computing global earth topographic corrections[END_REF] approximately. Due to the singularity at x = y, only the long-range interactions can be efficiently approximated by low-rank decompositions, i.e. such that d in ( 9) is not too large. Factorization of approximate decompositions for long-range interactions and precomputations are at the foundation of FMM algorithms. Some key quantities at each levels are computed by recursion formulae in the octree, see Figure 5. The short-range interactions are computed using the exact formula (8)

4 2 3 3 4 g I4 m g I3 m g I2 m=l I2 m (σ j ,y j ) W I4 m W I3 m W I2 m l I3 m l I4m (σ j ,y j ) f(x i )
Figure 5: Schematic representation of the computation steps in the FMM algorithm. The first line represents a 1D tree of height 4, and is the starting points of the algorithm. The following lines represent the following steps in the algorithm. First, the vectors W (interim calculation object) are computed using local σ j and y j for each box (P2M step). Then, in the downward stage, the local W vectors are assembled at each level (the size and location of the boxes change from one level to another), using W vectors from the previous level (M2M). At each level, the g vectors combine the W vectors from the interaction list (in green) of the considered box (in red) (M2L step). The upward stage consists in computing vectors l, using the vectors g of the considered level, and the vectors W of the previous level (L2L step). Finally, at the last level, the vectors l correspond to a component of the long-range interactions, to which a local quadratic step is added to compute the searched quantity f (x i ) (P2P step). in a local quadratic complexity. The height of the octree is increased as N and M increase. With the right choice of the octree height with respect to N and M , the summation ( 8) is computable in linear complexity and with a controllable error.

Remark 1 (The octree height determines the nature of the interaction) In Figure 5, we notice that the height of the octree determine the size of the leaves (the boxes at the last level). This size determines the limit at which the interaction between an observation point and a source point is close-or long-range in the octree sense.

Among the available FMMs, some are kernel specific and other are kernel independent. In our case, a kernel independent algorithm is desirable, since the first and second derivatives of the potential require the use of two different kernels. We use the Scalfmm C++ library [START_REF] Coulaud | Scalfmm, C++ Fast Multipole Method Library for HPC[END_REF], and among the available algorithms, we choose the BlackBox FMM which was proposed in [START_REF] Fong | The Black-box fast multipole method[END_REF]. In the BlackBox FMM, the long-range interactions in the octree sense are approximated using a tensorization of 1D Chebyshev polynomials. For a given octree height, the accuracy of the computation is controlled by the number of Chebyshev node per dimension, and by a tolerance parameter used to accelerate a certain step in the FMM algorithm. Usually, the number of Chebyshev nodes and the tolerance parameter are changed together, in a way that minimizes the execution time while guarantying a given accuracy for the result. In what follows, when the tolerance parameter is not specified, it means that it is chosen accordingly to the number of Chebyshev nodes.

Coupling FMM accelerated numerical quadrature and analytic Okabe formulae

Since we are going to use the FMM, the sets of tetrahedrons and observation points are organized in an octree. We take for the location of a tetrahedron the coordinates of its barycenter. A tetrahedron is close to x if it lies in the same box of the octree or in an adjacent box, and away from x otherwise, see Figure 4.

The quadrature formula (6) for long-range interactions (in the approximation sense) can be rewritten, for 1 ≤ ≤ 3,

g (x) V k ≈ G S s=1 [ω V k ,s ρ(y V k ,s )(y V k ,s ) v(V k )] 1 x -y V k ,s 3 -Gx S s=1 [ω V k ,s ρ(y V k ,s )v(V k )] 1 x -y V k ,s 3 , (11)
which consists of two summations of the form [START_REF] Cady | Calculation of gravity and magnetic anomalies of finite-length right polygonal prisms[END_REF], where

σ j is ω V k ,s ρ(y V k ,s )(y V k ,s ) v(V k ) and ω V k ,s ρ(y V k ,s )v(V k
) respectively in the first and second summation, and K(x, y) = 1 x-y 3 in both summations. Hence, the quadrature formula can be accelerated using the FMM as explained in Section 3. The difference is that the close-range interaction (in the approximation sense) are not computed using the quadrature formulae [START_REF] Cooley | An algorithm for the machine calculation of complex fourier series[END_REF], but using the analytic Okabe formulae.

Remark 2 (Issue with high-order quadrature formulae) The box I containing a given tetrahedron V k of the mesh is determined solely by the location of the barycenter of V k . This means that when computing the interaction between V k and a far-away observation point x using a quadrature formula, all the quadrature points y V k ,s , 1 ≤ s ≤ S, in [START_REF] Cooley | An algorithm for the machine calculation of complex fourier series[END_REF] are seen by the FMM as being located in I. However, this may not be the case, see Figure 6. The consequence is that, in the FMM algorithm, some Chebyshev polynomials will be evaluated outside their domain in which the approximations are trustworthy. Notice however that in the case of a one-point quadrature formula, the unique quadrature point is the barycenter of V k , and the phenomenon described here never happens.

Constraints on the octree height

In an FMM algorithm, there is a competition between the short-range interactions (in the octree sense), which are computed faster as the octree height increases, and the long-range interactions (in the octree sense), which are computed slower as the octree height increases. There is a trade-off value for the octree height that minimizes the overall execution time. In the FMM summation as described in [START_REF] Cady | Calculation of gravity and magnetic anomalies of finite-length right polygonal prisms[END_REF], with homogeneous kernels, the accuracy does not depend on the octree height. As a consequence, the octree height can be freely chosen to optimize execution time.

Recall the definition 3.1 on close-and long-range interactions in the octree and approximation senses. The ability to determine, for all possible interactions in the considered problem, which one in close-or long-range in the approximation sense is crucial since it determines which formula (Okabe or quadrature) should be used. As said at the beginning of this section, only the quadrature formulae can be FMM accelerated, due to their particular form. In the present algorithm, we want to be able to use the Okabe formulae for close-range interactions in the octree sense and quadrature formulae for long-range interactions in the octree sense, that will be accelerated by the FMM. The stake is hence to find a regime where the closeand long-range interactions in the approximation and octree sense are the same. The octree height determines the size of the leafs in the octree, and therefore the limiting distance between short-range and long-range interactions in the octree sense, see Remark 1. Equating close-and long-range interactions in the approximation and octree sense actually imply that we add a constraint in the determination of the octree height. In what follows, we no longer distinguish the definitions for close-and long-range interaction since we now enforce, through a constraint on the octree height, the two regime to be equal.

We define the "aspect ratio" of the mesh by this limiting distance divided by the mean edge of the mesh. From Figure 2, we see that for a chosen accuracy, there is an interval of values for admissible (i.e. leading to accurate enough results) aspect ratio, and corresponding admissible octree height. It may happen that the optimal octree height (in terms of execution time) is not among the admissible octree heights. In that case, we choose the admissible octree height that is the closest to the optimal one. Then, we might not reach in all cases the asymptotic linear complexity of the FMM.

Influence of the other parameters

The execution time is increased when the mesh size, the order of quadrature and the number of Chebyshev nodes increase, by increasing the size of the considered summation. These parameter improve the accuracy of the result in the following ways: (i) the mesh size improves the quality of the model, (ii) the order of quadrature improves the approximation of the volume integration of the long-range interactions, (iii) the number of Chebyshev nodes improves the approximation of the kernel decomposition used in the FMM.

In Section 6, we present a campaign of numerical experiments that aims to check the effects of these parameters and find the set of parameters leading to the fastest computation guarantying a given relative error.

Empirical selection of the parameters

In what follows, we present a numerical campaign aiming to establish empirical landmarks to select the values of the parameters: octree height, order of quadrature and number of Chebyshev nodes. The goal is to minimize the execution time for the computation of g, while guarantying a given accuracy. Then, the scalability of the algorithm is checked by measuring execution time on large test cases.

Numerical campaign

The mesh of this test case is represented on Figure 7. It consists of a cube denoted C, centered at the origin, of edge length 0.2 (C = (-0.1, 0.1) 3 ) and containing 13,712 tetrahedrons. We consider 1,000 observation points randomly taken in (-0.2, 0.2) 3 \C. The reference value is computed using the Okabe formula on the whole cube directly (the observation are close enough for the Okabe formula to be very accurate) and the relative accuracy is computed using [START_REF] Bouman | Goce gravity gradients: Combination with grace and satellite altimetry[END_REF]. In Table 1 are reported relative accuracy of the computation with respect to the octree height and the order of quadrature. The BlackBox FMM is used with 8 Chebyshev nodes per dimension and a tolerance of 10 -10 on the compression of the M2L operators, which enables an relative accuracy of at least 10 -8 on the FMM approximation. Since the relative accuracy measured in Table 1 are higher than 10 -8 , they are free of FMM approximation errors. We observe that the lower the octree height is, the more accurate the computation is. This indicates that, among the two competing effects (i) inaccuracy of quadrature formula for short-range interactions and (ii) inaccuracy of Okabe formula for long-range interactions, the Relative accuracy of the FMM computation Order of quad. \ Octree height Table 1: Influence of the octree height and the order of quadrature in the accuracy of the result. The reference value is computed using Okabe directly on the whole cube. The BlackBox FMM is used with 8 Chebyshev nodes per direction, and a tolerance parameter of 10 -10 first one is predominant. Indeed, in small height octrees, the leafs are large, and the limiting distance between short-range and long-range interactions is large. This means that, in the trade-off, it is better to compute some long-range using Okabe than some short-range using quadratures. We also observe that the higher the order of quadrature is, the more accurate the result is, which was expected. We observe that for a relative accuracy of 10 -4 , we can choose either an order of 1 with an octree of height 5 and below, or higher orders of quadrature with an octree of height 6 and below. Figures 8 and9 show the execution time of the computation of the short-range interactions (Okabe) and the long-range interactions (far FMM) with respect of the octree height, for various values for order of quadrature and number of Chebyshev nodes. As expected, the execution time for Okabe is not sensible to the order of quadrature of the number of Chebyshev nodes, since these parameters only affect the computation of long-range interactions. The smaller the octree height, the larger the execution time: since the leaves are larger, the ratio of short-range interactions divided by the total of interactions is larger. For the longrange interactions, the last effect is inverted. Finally, as expected, the larger the order of quadrature and the number of Chebyshev points, the larger the execution time for the longrange interactions. As a conclusion, for fast computation, one should choose some optimal octree height, a low order of quadrature and few Chebyshev points.

In Table 1, the FMM approximation (via the number of Chebyshev nodes) is chosen to be very accurate in order to evaluate the effect of the octree height in the accuracy of the algorithm. In Figure 10 are reported the relative accuracy of the algorithm with respect of the octree height, for various values of the order of quadrature and the number of Chebyshev nodes. This time, the effect of the number of Chebyshev nodes is visible: taking 4, 5, and 6 Chebyshev points introduce lower bounds on the relative accuracy of respectively roughly 10 -3 , 10 -4 , and 10 -5 . Consider we want the fastest set of parameters that allows a relative accuracy of approximately 10 -4 . From Figure 8, the fastest parameters are order of quadrature 1, 4 Chebyshev nodes and octree height 5. From Figure 9, the relative accuracy of this combination is close to 10 -3 , which is too large. The second fastest set of parameters is order of quadrature 1, 5 Chebyshev nodes and octree height 5. From Figure 9, the relative accuracy of this combination is less than 2 × 10 -4 , which we decide to choose. With an octree height of 5, the aspect ratio (the edge of a leaf in the tree divided by the mean edge of the mesh) is approximately 2.

Remark 3 A natural trick to reduce the computation time consists in saying that due to the rapid decreasing of the kernel, on can neglect the physical long-range interactions. However, with an aspect ratio of 2, a tetrahedron located at more than twice the mean edge of the mesh of an observation has a long-range interaction, but cannot be neglected. Hence, the physical longrange limit is much larger than the one in the octree sense, for an octree of recommended height. Computing physical long-range interactions with no FMM acceleration results in computation of still quadratic asymptotic complexity (with a better constant than when computing all the interactions). As a conclusion, the empirically selected set of parameters for fast computations and relative accuracy of approximately 10 -4 is: order of quadrature 1, 5 Chebyshev nodes and an octree height such that the aspect ratio is 2. We recall that with on order of quadrature of 1, the issues described in Remark 2 do nor occur. In what follows, we denote by "recommended height", the octree height that corresponds to an aspect ratio of 2.

The modeling error

When deriving the mesh of a 3D geometry, an modeling error is introduced, in particular for objects with curved boundaries. In this section, we consider a 3D ball, and a set of four meshes represented in Figure 11. A set of 1,000 observation points are taken uniformly along a radius of the ball. Figure 12 shows the known exact solution for the x-component of the vector g, and the difference between the exact and computed solution using our algorithm, on Meshball 4, with order of quadrature 15, an octree of size 6 (corresponding to an aspect ratio of 2), and 5 Chebyshev nodes.

Figure 12: Analytic solution and difference between the computed and analytic solutions for the x-component of g along a radius of the ball on MeshBall 4 for an octree of height 6 and order of quadrature 15 Figure 13 presents the relative accuracy of the x-component of the vector g with respect to one parameter, while the other two are kept constant, for MeshBall 1 to 4. We observe that in this case, the variation of the parameter have a little influence on the relative accuracy, whereas the choice of the mesh has an important effect. This indicates that while we have derived an empirical criterion for ensuring a relative accuracy of 10 -4 for a computation on any mesh, the modeling error was not considered, and can sometimes be significant. Hence, it does not appear useful to consider an accuracy better than 10 -4 .

6 Numerical experiments

Performance on large scale computations

In what follows, we assess the scalability of the algorithm by running large test-cases on a 32-procs cluster AMD Opteron Processor 6276@2.3 GHz. The test case is the same as in Section 5.1, and we checked that the computed relative error was approximately 10 -4 using our empirically selected set of parameters.

Figure 14 shows execution times and acceleration factors for test-cases ranging from 10 3 to 10 7 tetrahedrons. The number of observation points is the same as the number of tetrahedrons. The unaccelerated curve is extrapolated in the following way: the computation is done on one processor using the Okabe formula on a 6 tetrahedron mesh of the cube depicted in Figure 7. Then, the obtained execution time is scaled to the number of tetrahedrons of the test case, then divided by 32, the number of processors used to apply our algorithm. For 10 7 tetrahedrons and observation points, our algorithm takes 47 min, while the unaccelerated one would take more than 3 years. We observe that the execution time increases with respect to the number of tetrahedrons and observation points is linear starting from roughly 10 5 (and slower before). The 10 7 limit is due to the memory size of our cluster. A test-case with 10 8 tetrahedrons and 2 × 10 6 observation points (with corresponds to one observation points per pixel on a full HD screen) has been computed on our cluster.

Topographic correction over the Himalaya

In this section, we present the topographic correction over the Himalaya. The topographic correction is the computation of the gravitation effects created by the difference between the real Earth and the reference WGS84 ellipsoid. Each volume element above the altitude zero of the reference ellipsoid is supposed to have a density of 2700 kg/m 3 , and each element of volume below the altitude zero is supposed to have a density of -1900 kg/m 3 (which corresponds to the density of the oceanic crush subtracted to the density of water). The algorithm used to construct the mesh from the ETOPO1 database is detailed in [START_REF] Casenave | A volumic mesh for computing global earth topographic corrections[END_REF], see Figure 15 for an example of a large region around the Himalaya at an input resolution of 10 min. In this algorithm, a particular attention is paid to the modeling of the coasts: each time the altitude changes sign between two points, the intersection with the reference ellipsoid is computed, see Figure 16.

Figure 17 presents the result of the topographic correction over the Himalaya. The mesh has an input resolution of 1 min, and consists of 1.2 × 10 7 tetrahedrons. There are 10 6 observation points, located 1 cm above the topographic points. The octree has a height of 11, and has 1,327,186 leaves. The computation takes < 20 minutes on a 16-procs cluster AMD Opteron Processor 6276@2.3 GHz. The zoom on mount Everest is computed using 3.7 × 10 3 topographic points (from the same dataset ETOPO1) and 4.5 × 10 5 observation points.

Figure 18 shows separately the close-and long-range interactions for the topographic correction over the Himalaya. First, we notice that both components contain a physical signal. Then, the circular artifacts illustrates than when two observation points are close to each other, but in different leaves of the octree, the sets of tetrahedrons in close-range interactions associated to the two observation points are translated from one another by at least twice the mean edge of the mesh.

Computation of the topographic correction over France

Figure 19 show the mesh used for the computation. Notice that no geometry approximation is done, in particular the Earth curvature is visible on the mesh.

Figure 20 shows the radial component of g, the observation points being located 1 cm above Figure 20: The effect of the topography's gravity field over France, radial component of g, in mGal (1 mGal = 10 -5 m.s -2 ). The observation points are taken 1 cm above the topographic points on the ground, and 1 cm above the sea level in maritime areas

The data is composed of 2,804,433 tetrahedrons and 422,521 observation points. The octree height is 10, the order of quadrature is 1 and the number of Chebyshev nodes 5. The computation takes 2 min 32 s on our 32-procs cluster for the long-range interactions, and 4 min 32 s for the short-range interactions. The postreatment to convert the vector g from cartesian to latitude-longitude coordinates takes 58 s on one processor.

Computation of the global topographic correction at the altitude of GOCE

The satellite GOCE (Gravity field and steady-state Ocean Circulation Explorer), whose mission was to map in detail the Earth's gravity field [START_REF] Bouman | Goce gravity gradients: Combination with grace and satellite altimetry[END_REF][START_REF] Panet | Mapping the mass distribution of Earth's mantle using satellite-derived gravity gradients[END_REF], flew at an altitude of approximately 250 km. We compute here Earth's global topographic correction at such an altitude. Highresolution global topography and topographic corrections have been considered in [START_REF] Balmino | Spherical harmonic modelling to ultrahigh degree of bouguer and isostatic anomalies[END_REF][START_REF] Claessens | Ellipsoidal topographic potential new solutions for spectral forward gravity modelling of topography with respect to a reference ellipsoid[END_REF][START_REF] Hirt | Earth2014: 1' shape, topography, bedrock and ice-sheet modelsavailable as gridded data and degree 10,800 spherical harmonics[END_REF].

The computation is carried out on a 16-procs cluster AMD Opteron Processor 6276@2.3 GHz. Some number on the computation are given in Figure 21 shows the gravity fields generated by the difference between the real Earth and the WGS84 ellipsoid at 250 km. In this figure, the red color refers to high values, whereas the blue color refers to low values. Figure 21 shows that at 250 km, the signal is mainly composed of large wavelengths. To demonstrate the capabilities of our methodology and the resolution of our mesh, we include the global topographic correction at 20 km of g r and T rr in respectively Figure 22 and 23.

67P/Churyumov-Gerasimenko (67P) gravity field

Figure 24 shows the mesh used to compute the gravity fields created by the comet 67P, containing 989,950 tetrahedrons. The largest dimension of the comet is 5 km, and we compute the vector g on 1,994,004 observation points located at the surface of the sphere centered in the geometric center of the comet; see Figure 25. The computation is carried out on a 4-core Intel Core i5@2.67 GHz, using a height-8 octree. BlackBox FMM is taken with 5 Chebyshev nodes per dimension. The FMM summation is computed in 251 s, while the postreatment to plot the components of g takes 68 s. Figure 26 shows a 3D representation of 1,000 vectors g, and the tangential components of 2,048 vectors g is represented in Figure 27. The asymmetry of the latter field is due to the asymmetry of the two lobes of the comet.

Conclusion

In this work, we have presented a new algorithm to compute gravitation anomalies for a general 3D source. Depending on the close-or long-range nature of the involved interactions, the algorithm automatically switches between analytic integration formulae and numerical quadratic formulae, and relies on the Fast Multipole Method to drastically increase the speed of the computation of long-range interactions while guarantying a chosen accuracy.

Figure 1 :

 1 Figure 1: Representations of Mesh 1 to 4

Figure 2 :

 2 Figure 2: Relative difference (7) between the Okabe formula on Mesh 1 and quadrature formulae on Mesh 1 to 4 for various distance r to the center of the tetrahedron; left: 1-point quadrature, right: 15-point quadrature

Figure 3 :

 3 Figure 3: Example of an octree of height 1: in black is the level 0 (the root) containing the whole domain ; in red are represented the boxes of the level 1 (8 boxes: 2 per direction).

Figure 4 :

 4 Figure 4: 1D tree; consider a box I, N (I) denotes the set of close boxes and F (I) denotes the boxes far-away from I. Let x ∈ I, the interaction between x and a point y is considered close-range if y ∈ N (I) and long-range if y ∈ F (I), in the approximation sense.

Figure 6 :

 6 Figure 6: 2D illustration of Remark 2: the square represents the box containing the triangle, since its barycenter (the cross) lies in it. Any quadrature points in the red part of the triangle, outside the box, will be considered inside the box by the FMM algorithm.

Figure 7 :

 7 Figure 7: Meshcube, centered at the origin, edge length: 0.2, 13,712 tetrahedrons

3 × 3 . 5 ×

 335 10 -6 3.9 × 10 -5 1.3 × 10 -4 5.5 × 10 -3 3.7 × 10 -2 5 10 -7 4.4 × 10 -7 1.5 × 10 -6 3.3 × 10 -5 2.9 × 10 -4 11 3.4 × 10 -7 5.5 × 10 -7 2.2 × 10 -6 9.5 × 10 -5 0.23[START_REF] Forsberg | Gravity field terrain effect computations by FFT[END_REF] 2.1 × 10 -6 2.0 × 10 -6 2.1 × 10 -6 3.1 × 10 -5 1.4 × 10 -2

Figure 8 :

 8 Figure 8: Short-range interactions (Okabe) and the long-range interactions (far FMM) execution time with respect to the octree height for orders of quadrature 1 and 5, and number of Chebyshev nodes 4, 5 and 6.

Figure 9 :

 9 Figure 9: Short-range interactions (Okabe) and the long-range interactions (far FMM) execution time with respect to the octree height for orders of quadrature 11 and 15, and number of Chebyshev nodes 4, 5 and 6.

Figure 10 :

 10 Figure 10: Relative accuracy with respect to Octree height for orders of quadrature 1, 5 and 11, and number of Chebyshev nodes 4, 5 and 6.

Figure 11 :

 11 Figure 11: Representations of MeshBall 1 to 4

Figure 13 :

 13 Figure 13: Difference between the computed and analytic solutions for the x-component of g along a radius of the ball on MeshBall 1 to 4. Left: 5 Chebyshev interpolation points par dimension and recommended octree height, middle: 5 Chebyshev interpolation points par dimension and order of quadrature 5; right: recommended octree height, order of quadrature 5, and 5 Chebyshev nodes

Figure 14 :

 14 Figure 14: Execution time and acceleration factors for the unaccelerated computation and the FMM accelerated one; computation on a 32-processors cluster

Figure 15 :

 15 Figure 15: Triangular mesh (shaded representation) of the topography for 0 • ≤ θ ≤ 50 • and 60 • ≤ φ ≤ 120 • , input resolution of 10 min

Figure 16 :

 16 Figure 16: Top: left: the triangle element on the topographic mesh does not intersect the reference ellipsoid, right: the triangle element on the topographic mesh intersects the reference ellipsoid; bottom: coarse mesh of the coasts in South-East Asia: the locally refined coasts illustrate the cases with intersection.

Figure 17 :

 17 Figure 17: Topographic correction over the Himalaya, with a zoom over mount Everest

Figure 18 :Figure 19 :

 1819 Figure 18: Topographic correction over the Himalaya; left: close-range interactions, right: long-range interactions

Figure 21 :

 21 Figure 21: Visualization of the gravity fields generated by the difference between the real Earth and the WGS84 ellipsoid (i.e. topographic correction) at 250 km; l.1: φ (in m 2 s -2 ), l.2: g r , g θ , g φ (in ms -2 ), l.3: T rr , T rθ , T rφ , l.4: T θθ , T θφ , l.5: T φφ (in s -2 ). The numbers written under each picture correspond to the min and max reached values.

Figure 22 :

 22 Figure 22: Topographic correction at 20 km: g r (in ms -2 ).

Figure 23 :

 23 Figure 23: Topographic correction at 20 km: T rr (in s -2 ).

25

 25

Figure 24 :

 24 Figure 24: Mesh for the comet 67P

Figure 25 :

 25 Figure 25: g r , g t , and g p at an altitude of one diameter (in ms -2 ).

Figure 26 :

 26 Figure 26: Representation of the vector field g; the color indicates the magnitude of g: red for the larger magnitudes and blue for the smaller ones.

Figure 27 :

 27 Figure 27: Representation of the tangential component of g; the color indicates the magnitude of g: red for the larger magnitudes and blue for the smaller ones.

Table 2 .

 2 Observation points 2 × 10 6 points, 250 km altitude Mesh 1.6 × 10 7 vertices, 3.3 × 10 7 triangles, 10 8 tetrahedrons FMM 10-level octree, 16 processors, 6m 30s computation time per component

Table 2 :

 2 Some numbers on the Earth topographic correction application computation