CUDA Parallel Algorithms for Forward and Inverse Structural Gravity Problems

Abstract

This paper describes usage of CUDA parallelization scheme for forward and inverse gravity problems for structural boundaries. Forward problem is calculated using the nite elements approach. This means that the whole calculation volume is split into parallelepipeds and then the gravity e ect of each is calculated using known formula. Inverse problem solution is found using iteration local corrections method. This method requires only forward problem calculation on each iteration and does not use the operator inversion. Obtained results show that even cheap consumer video cards are highly e ective for algorithm parallelization.

Keywords

• geophysics
• parallel algorithm
• forward gravity problem inverse gravity problem
• CUDA

-

Alexander Tsidaev

Introduction

Gravity eld measurement and interpretation is one of the simplest geophysical methods. It does not require additional preparations (unlike, e.g., the method of arti cal magnetization in magnetometry) or large amount of sensors (unlike seismometry). But despite this apparent simplicity, the interpretation of big amounts of data is very hard. For a long time only 2D density models along pro les were constructed, but now we have a need in 3D models construction as they are much more informative.

3D geophysical problems solution usually requires high-performance calculations, especially if input data is speci ed on a big grid. In many cases the performance of regular personal computers is not enough, this leads to the need of supercomputers usage. But currently the supercomputers are not widely distributed and accessing one of them can be problematic. Even taking into account that the data can be sent over the Internet to any distant supercomputer, the big amount of data to send makes the process very slow.

But now we can use a di erent parallelization approach. With big popularity of videogames, a standard video card of personal computer became much more than just a display adapter. It is a mini-computer with higher performance than central processing unit (CPU) itself. Performance of top-class graphical processing units (GPU) is similar to the performance of 10-years old supercomputers. But even cheap consumer GPU can provide a high speed up of calculations if program was parallelized correctly.

Compute Uni ed Device Architecture (CUDA) is a parallelization technology by one of the leading GPU producers, nVidia. Current paper is dedicated to the usage of CUDA for the creation of e cient programs for gravity geophysical problems. Section 2 contains description of gravity structural forward problem and algorithm for its parallelization using CUDA. Section 3 describes the inverse gravity problem for structural boundaries and the idea of local corrections method. In the 4 section the comparison of di erent calculation schemes is provided. 2

Forward gravity problem for structural boundary In general, gravity potential that is produced by some object is calculated as (for details one may refer to, e.g., [1])

W (x; y; z) =

( ; ; )dV p( x)2 + ( y)2 + ( z)2 ( 1 ) Z V

Z

V

Here x; y; z - point of observation; ( ; ; ) is the density value at a point ( ; ; ) under Earth surface; = 6:67408 10 11m2 kg 1 s 2 is a gravitational constant. z axis direction is downwards.

If we neglect the sphericity of Earth, the rst derivative by z of ( 1 ) can be taken as anomaly gravity eld g on Earth surface. This is the main formula for gravity forward problem: g = Wz0 (x; y; z) = ( ; ; )(

z)dV [( x)2 + ( y)2 + ( z)2]3=2 ( 2 )

It is clearly seen from ( 2 ) that gravity eld on Earth level depends only on density distribution. While density values may vary between neighbour points, in regional studies of big squares we can ignore this variance. In such scales Earth gradiental density distribution could be replaced with the structural model of homogeneous layers with constant densities (Fig. 1(a)). In [2] it is shown that the topography of structural boundaries between these layers is the only element, which produces non-constant gravity eld in this model (Fig. 1(b)). This is 3D analogue of Strakhov boundaries parameterization, which was presented in [3].

As the result, we do not need to know absolute density values for layers and the gravity e ect of structural boundary could be calculated up to a constant as 2 1 z(x; y) z(x; y)

H = 0 H = 0 Hasympt here (x0; y0) - point of observation on Earth level (where z0 = 0); z(x; y) structural boundary surface; H - asymptote of z(x; y), - density jump on a surface z(x; y), it is equal to the di erence ( 1 2) between density values of layers above ( 1 ) and below ( 2 ) the structural boundary.

Direct calculation using formula ( 3 ) is complicated by the weak singularity in the point (x = x0; y = y0; z(x; y) [ H = 0), so the formula is inapplicable for the shallow surfaces. To avoid this limitation in [2] we o ered an idea of more stable and e cient structural gravity problem calculation using the nite elements approach. We split the observation plane into a set of rectangular elements xi = i; yi = i (0 i < N; 0 j < M , the same discretization is used for the boundary too) and calculate a eld g(x0; y0) as a sum X ij g(x0; y0) = gij (x0; y0) ( 4 ) of gravity e ects of parallelepiped elements (whose formula could be obtained analytically) gij (x0; y0) = ( y0) ln(j

x0j + R) + ( arctg ( x0) (

R y0)

x0) ln(j i+1 j+1 z( i; j) i j

H y0j + R) ( 5 ) Here R = p( x)2 + ( y)2 + 2.

In addition to the possibility to calculate eld of near-surface boundaries, this method has another advantage over direct calculation using ( 3 ): there is no need for the \real" asymptote for a boundary. Since the boundary gravity e ect calculation is actually replaced with calculation of gravity e ect for the limited body, any horizontally located plane could be used instead of asymptote H. But it is better to take the average value of the boundary: for all other H the eld will contain additional constant component, which is equal to the gravity e ect of a at layer located between H and boundary average depth.

It is obvious that the algorithmic implementation of ( 4 ) should contain two nested loops: one iterates through points of observation (x0; y0), and another iterates through the surface discretes ( ; ). This is very close to the CUDA internal structure, which has blocks and threads [4]. Programmer can run any number of calculation blocks, and each block will execute the same program code in parallel in prede ned (and permanent for all blocks) number of threads. So if we run calculation of ( 5 ) on N M blocks in N M threads and then summarize results for all threads inside each block, the combination of the block results will be exactly the gravity e ect of the structural boundary (up to ). The described scheme is presented on Fig. 2.

This method also provides a good scalability across multiple GPUs. Each graphic card can calculate gravity eld in a separate set of grid rows, and host then can easily combine results from di erent GPUs together. Number of graphic cards in single PC is limited by the number of physical PCIe slots and PCIe channels, but it is almost always possible to insert at least two GPUs in one motherboard. Another option is to use several CUDA-enabled PCs connected with MPI (Message Passing Interface) over a local network. But one should control the performance of the process due to a big overhead of network operations. 3

Inverse gravity problem for structural boundary Inverse gravity problem occurs when it is needed to restore density distribution by known (observed) gravity eld. This problem has big practical necessity because it allows one to explore Earth structure mathematically, without such costly operations as borehole drilling. But the problem is ill-posed so the solution is not unique (there are continuum of solutions that produce the same gravity e ect) and the solution process of the problem is unstable.

Direct inversion of ( 3 ) gives us the Fredholm equation of rst kind. It can be reduced to a system of linear equations, which then could be solved using Tikhonov regularization scheme and approximate methods. This approach has a good description in [5], one of possible parallelization schemes is provided there as well. But for the big grids this system of linear equations will be very large (i.e. N 2 M 2) and so some additional algorithms would be needed to operate such a big data. This is the reason to use iterative method of local corrections.

The method of local corrections [6, 7] was developed by I. Prutkin and its idea is that the main contribution to the gravity eld at some point (x; y) is introduced by the masses located directly under that point. Since the boundaries between density layers rarely have large gradients, this assumption can be considered as true. Main formula of the method is simple: zxny+1 =

zxny 1 + c zxny Uxy

Uxny ( 6 ) here zxny is the position of a boundary z in a point (x; y) on nth iteration; Uxy is the observed eld in the point (x; y); Uxny is the gravity eld of nth iteration boundary zxny; c is cubature formula coe cient; is regularization parameter.

Parameter is very important. This is not the regularization parameter in Tikhonov's sense. It acts more like stabilizing factor, which does not allow boundary to deviate too far from the previous approximation. It should be selected relatively small, 0 < 1.

I. Prutkin selected the estimated asymptote level for a z(x; y) as the initial approximation z0. In [2] we o ered modi ed technique. We construct initial approximation using a priori data obtained using another geophysical methods (e.g. seismic surveys). This reduces the number of iterations and also forces the method to produce more geologically meaningful solutions.

It is clearly seen from ( 6 ) that the only calculation that is needed on each iteration of local corrections method is the forward gravity problem solution for a density boundary of previous approximation. So the parallelization idea is quite simple and relies on parallelization algorithm described in section 2: 1. Calculate eld of the boundary of nth approximation using CUDA parallel scheme. 2. Calculate di erence U U n between observed and obtained elds. 3. If mean value of the di erence is less than some preselected ", then exit.

Current approximation of z(x; y) is the solution. 4. Else calculate (n + 1)th approximation of the boundary by ( 6 ) and repeat algorithm from the rst step. 4

Calculations and results

Algorithms were implemented for di erent hardware to evaluate e ciency. Program for CPU was written in C++ and compiled with gcc (g++ -o3) compiler. GPU implementation is also in C++ and CUDA part was compiled with proprietary nvcc compiler while gcc has been used for the host code. The grid with of 256x256 nodes was taken as input data. Tesla M2050 GPUs were accessed on \Uran" supercomputer, which is located in Krasovskii Institute of Mathematics and Mechanics (Yekaterinburg). Results of testing are presented in Table 1.

So for a big grids calculations on CPU are almost unacceptable. Though CPU still could be used for a single forward problem calculation, its usage in iterative process of inverse problem solution entails long waiting time. Local corrections method requires parameter preselection which could be performed only empirically, so the number of required iterations (usually 20-30) is multiplied by the amount of experiments. This makes CPU usage for inverse problem calculation highly problematic.

Conversely, even relatively cheap outdated consumer GPU nVidia GTX580 has acceptable time of calculation. And modern consumer devices like nVidia GTX780 Ti can compete even with previous generation of specialized CUDA calculators (Tesla). The one who has no access to supercomputer can easily build high-performance computer himself using these widely available devices.

As it is seen from Table 1, the dependence between number of GPUs and calculation time is almost linear for 1, 2 and 4 Tesla cards. Regularity break on 8 GPUs is related to the increased overhead of host operations. For a bigger grids the linear nature of dependence returns. For example, for a grid of size 1335x1404 the calculation times on 4 and 8 GPUs are 20m12s and 9m47s respectively. This linear dependence con rms very good scalability of CUDA calculations: with twofold GPU count increase the duration of calculation is reduced in 2 times too. 5

Conclusion

Gravity modeling is very important geophysical task. In this paper the algorithms for forward and inverse problems parallelization were o ered. Obtained results show that GPU calculations are very e ective way to lower the calculation time. Even widespread graphic cards, which initially were intended for a videogame use, can be utilized for a high-performance calculations. Simpleness of CUDA make this technology available for a mass use while scalability allows one to increase performance easily.

The research was supported by the Russian Science Foundation (project 1427-00059).

References

1. Blokh , Yu.I. : Interpretation of gravity and magnetic anomalies (In Russian) , http://sigma3d.com/pdf/books/blokh-interp. pdf ( 2009 )

2. Martyshko , P.S. , Ladovskiy , I.V. , Tsidaev , A.G. : Constructional of Regional Geophysical Models Based on the Joint Interpretation of Gravity and Seismic Data . Izvestiya, Physics of the Solid Earth , 2010 , Vol. 46 , No. 11 , 931 { 942 ( 2010 )

3. Strakhov , V.N. : On the parameterization problem in gravity inversion (In Russian) . Fizika Zemli, No. 6 , 39 { 50 ( 1978 )

4. Wilt , N.: CUDA handbook: A comprehensive guide to GPU programming . AddisonWesley Professional , Boston ( 2013 )

5. Akimova , E.N. , Vasin , V.V. , Misilov , V.E. : Algorithms for solving inverse gravimetry problems of nding the interface between media on multiprocessing computer systems (In Russian) . Vestnik UGATU , vol. 18 , no. 2 ( 63 ), 208 { 217 ( 2014 )

6. Prutkin , I.L. : The solution of three-dimensional inverse gravimetric problem in the class of contact surfaces by the method of local corrections . Izvestiya. Phys. Solid Earth 22 ( 1 ), 49 - 55 ( 1986 )

7. Prutkin , I. , Saleh , A. : Gravity and magnetic data inversion for 3D topography of the Moho discontinuity in the northern Red Sea area, Egypt . Journal of Geodynamics 47 ( 5 ), 237 - 245 ( 2009 )