In this paper we purpose a performance-e ective algorithm for solving forward gravity problem for elliptical objects. The algorithm is subject to parallel computations, it was implemented and tested with the CUDA technology. In general, the suggested method can be applied to density distribution models of arbitrary geometry that can be triangulated.
(p)dVp where is the gravitational constant, dVp is a volume element of integration, nq is the outer normal to S in the orthogonal projection of the point q to S, r and r0 are the radius vectors of the points p and q, respectively.
Assume the ellipsoid of rotation for our spherical model D has equatorial and polar radii equal to a and b. The position of the point on the surface of the ellipsoid of rotation is uniquely de ned by the geodesic latitude B 2 [ =2; =2] and the longitude L 2 ( ; ] (except for the \poles", where longitude is not de ned). As a third coordinate H we employ the distance between the given point and the ellipsoid surface; it has a plus sign if the point lies outside the ellipsoid, and the minus sign if the point lies inside it. Position of the given point in space is uniquely de ned by three values: (L; B; H), where H > N (1 e2), e = pa2 b2=a (the eccentricity of the ellipsoid), and N = a=p1 e2 sin2 B. 2
Assume Di;j;k is a descrete element of D and the density of a descrete element be a constant i;j;k. Then the eld g at point q generated by D is equals to where Gi;j;k(q) is the eld at the point q of Di;j;k up to a coe cient . In order to calculate Gi;j;k(q), we introduce in the space of a spherical model a rectangular Cartesian coordinate system O0x0y0z0. The center O0 coincides with the center of the ellipsoid, O0z0 axis is ellipsoid's axis of rotation and directed from the \south pole" to the \north" (i.e. when B = =2, then z0 > 0), the O0x0 axis points to (0; 0; 0) in (L; B; H), the axis O0y0 complements the system to the right-handed and axis set. Equations for transition from (L; B; H) to O0x0y0z0 are then as follows: 8x0 > <
y0 >:z0 = (N + H) cos B cos L = (N + H) cos B sin L = (N (1 e2) + H) sin B
; 0cos B cos L1 n = @cos B sin LA ;
sin B where n is a unit vector of the external normal to the ellipsoid at the point (L; B; 0).
Integral (
Gi;j;k(q)
G^i;j;k(q) = (p)dVp :
In formula (
Z
D^i;j;k rp jr 1 r0j
! dVp = nq;
I
Next, we divide the surface integral into parts along the faces of D^ i;j;k and take into account that the outer normal np at the integration point is constant for each face
G^i;j;k(q) =
X Si12S(D^i;j;k) (nq; ni1)
Z
Si1 jr dS r0j ; for the integral RSi1 jr where ni1 is the outer normal to the face Si1. Now we need to acquire a formula dS
over a triangle (which is the gravitational potential r0j of a triangle with a unit surface density without coe cient). Let ri (i = 1; 2; 3) be the radius vector of the vertices of the triangle hp1; p2; p3i; q is the eld computation point; ai = ri r0; a(j;i) = ai aj = ri rj ; N = a(i 1;i) a(i;i+1) is the normal to the plane of the triangle with a length equal to its doubled area; n = N=jNj is the unit normal to the plane of the triangle; A(j;i) = a(j;i)=ja(j;i)j; (ai n) is the distance from the point q to the plane of the triangle (with the appropriate sign). Then
Z
dS
+ X3 [ai 1; ai] arctanh
1. for each element Di;j;k using coordinates transformation (
In the presented algorithm, computation of formula (7) is taking most of the time of the running program. Indeed, in order to calculate the eld of body that consists of Nb discrete elements in Nf points of space, the expression should be calculated N = kNbNf times, where k is the amount of triangles in one element of the body. In practice, N is of order of 1015 for high-resolution Earth crust models (with surface area of order 1000 1000km and depth of 100km). Without using contemporary parallel technologies such computations would become senseless since it would take months or even years to run. Hence, our goal is to utilize most powerful computational devices available, which are GPUs.
Here, we assume that all the preprocessing steps have been done and as our input we have a set of triangles T (represented by triplets of points in space) that has been obtained as a result of triangulation of the discrete elements of the body. Additionally, we have a set of points Q, in which we want to calculate the normal component (or a vector) of gravitational eld of the body. We are not going into details of this step because its computational complexity does not depend on Nf ; thus, its becoming insigni cant for large values of Nf .
The eld computation requires minimum two nested cycles through Q and T . The iterations of the outer loop (over Q) are mutually independent and can be split between di erent computing nodes or/and GPUs. The inner loop (over T ) can be expressed using map-reduce pattern: gq = reduce(mapq(T )), where g is the resulting gravitational eld vector the at the point q. The trick is that those operations are implemented with high e ciency in the CUDA Thrust library [5] in the form of a single function: thrust::transform reduce(). All we were left to do is to implement formula (7) and pass it as a parameter.
All the tests have been performed using two regional Earth crust models that have been previously reconstructed in Bulashevich Institute of Geophysics (Ural Branch of Russian Academy of Sciences) as a result of solving inverse gravity problem [2{4]. Test models characteristics are shown in Table 1. Krasovsky Ellipsoid is taken as the reference one.
Comparison of our most resent CPU available (i7-6900K 8 physical cores @ 3.2GHz) and our oldest GPU (GeForce GTX TITAN Black) showed that GPU outperforms CPU from 15 to 20 times even for considerably small workloads (Workload 1). So, future tests were conducted using only the GPUs.
This work was partly supported by the Center of Excellence \Geoinformation technologies and geophisical data complex interpretation" of the Ural Federal University Program and by the The Russian Foundation for Fundamental Research (project 17-05-00916 ).