Introduction

The problem considered in this paper is finding several interfaces between layers in a multilayered medium using known gravitational data [1,2]. This problem is described by a nonlinear integral Fredholm equation of the first kind; so, it is ill-posed. The real gravity measurements are carried out over a large area producing the large-scale grids. Processing the gravitational data is a time consuming process and requires a lot of memory. So, it is necessary to develop parallel algorithms for parallel computing systems.

In works [3,4], for solving the structural gravimetry problem, efficient parallel algorithms based on the modified conjugate gradient method were proposed and implemented multicore CPUs. The modification is based on approximation of the Jacobian matrix by calculating only significant elements. This method reduces the computation time in comparison with the full calculation approach.

Here, we construct a parallel algorithm on the basis of this method and implement it for GPUs using the CUDA technology.

We compare the CPU and GPU implementations in terms of computation time in solving the model problem with synthetic gravitational data.

Statement of the Inverse Gravity Problem for the Model of Multilayered Medium

We assume that the lower half-space is composed of several layers with constant densities, which are separated by the sought surfaces S l , l = 1..L, where L is the number of boundary surfaces (see, Fig. 1).

Fig. 1. Model of multilayer medium

The gravitational field due to this half-space is equal to the sum of the gravitational fields due to each surface. Let the boundary surfaces be specified by the functions z = ζ l (x, y), let the density contrasts on them be ∆σ l , and let the surfaces have the horizontal asymptotic planes z = H l . The field ∆g produced by the superposition of the boundary surfaces and measured on the Earth's surface z = 0 is found by the following equation (with accuracy up to a constant term of summation) [1]:

f l=1..L ∆σ l ∞ −∞ ∞ −∞ 1 (x − x ) 2 + (y − y ) 2 + ζ 2 (x , y ) − − 1 (x − x ) 2 + (y − y ) 2 + ζ 2 (x , y ) dx dy = ∆g(x, y, 0), (1)

where f is the gravitational constant.

This equation is the Fredholm nonlinear equation of the first kind of functions ζ l and is ill-posed: it has a nonunique solution which unstably depends on the initial data.

After the discretization of equation ( 1) on the rectangular grid n = M × N , where the right-hand side ∆g(x, y, 0) is given, and the approximation of the left-side integral operator A(ζ 1 , ..., ζ l ) by the quadrature rules, we obtain the vector of the right-hand side F with the length n, the combined surfaces vector

z = [ζ 1 (x 1 , y 1 ), ..., ζ 1 (x M , y N ), ..., ζ 2 (x 1 , y 1 ), ..., ζ 2 (x M , y N ), ..., ζ L (x M , y N )]

with the length Ln, and the system of nonlinear equations having the form

A(z) = F. (2)

3 Algorithm for Solving the Inverse Problem

In this work, to solve problem (2), we use the approach based on the modified linearized conjugate gradient method with weighting factors [5]. This method has the following form:

z k+1 = z k − ψ p k , S(z k ) A (z k ) − p k p k , p k = v k + β k p k−1 , p 0 = v 0 , β k = max v k , v k − v k−1 v k−1 , 0 v = γ • S(z), S(z) = A (z) * (A(z) − F ),(3)

where z k is the solution estimate at kth iteration, ψ is the damping factor, γ is the vector of weighting factors, A (z) is the Jacobean matrix of the discretized integral operator A(z), • is operation of componentwise vector multiplication.

In this work, we propose the following rules for selection of the weighting factors:

[F 1 , ..., F L ] → [f 1 , ..., f Ln ] → [γ 1 , ..., γ Ln ], F l → [γ n(l−1) , γ n(l−1)+1 , ..., γ nl ], γ l = f 2 i + µ max i=1..n { f 2 i + µ} , 0 < µ < 1, (4)

where µ is the smoothing parameter.

The fields F l are extracted from the total field F using heightwise transformation technique from [6].

The condition A(z) − F / F < ε for sufficiently small ε is used as the termination criterion for the iterative process.

Numerical Implementation

In works [3,4], the efficient technique for reducing the computing time was proposed and implemented. The main idea is to drop out small elements and replace matrix A with the block-band matrix A as shown in Fig. 2. The elements inside each square block depend on terms (x − x ) and (y − y ). Thus, the farther the matrix element is from the main diagonal of block, the less its value is. For each block, we get the bandwidth parameter 0 < β l ≤ 1 by automatic adjustment procedure. First, we find the maximal element a max l of the block. To find it, we just need to look over the main diagonal of the block. Then, we set the parameter β l in such way that the resulting band matrix elements would be larger then some threshold, i.e. a ij > θa max l for the threshold parameter θ. Most expensive part is calculation of the Jacobian matrix A and vector A(z k ) at each iteration. The matrix and vector are divided into a number of fragments, and each fragment is processed by its own thread. The elements of the Jacobian matrix are calculated on-the-fly, which means that the value of an element is computed when calling this element, without storing it in memory.

The adjustment of the kernel execution parameters for the grid size is an important problem. In previous works [7,8], we implemented the original method for automatic adjustment of parameters. For the reference 128 × 128 grid and M2090 GPUs, the optimal parameters were found manually. For the grid sizes divisible by 128, the reference parameters are multiplied by the coefficient. When using multiple GPUs, the x dimension is divided by the number of GPUs; i.e., the number of threads in the block is reduced while number of blocks in the grid remains constant.

This imposes some constraints on the input data and GPUs configuration: grid size should be divisible by 128 (128, 256, 512, 1024 ...) and GPUs number should be a power of 2 (1, 2, 4, 8, ...).

Numerical Experiments

To test the constructed modified algorithm and to compare it with the unmodified one in terms of execution time, we use the model problem of reconstructing three boundary surfaces on a large grid (512 × 512 nodes) using the quasi-real data.

The model gravitational field shown in Fig. 3 was obtained by solving the forward problem using three surfaces S 1 , S 2 , S 3 (see, Fig. 4a) with asymprotic planes H 1 = 10km, H 2 = 20km, H 3 = 30km from work [3] and density contrasts ∆σ 1 = ∆σ 2 = ∆σ 3 = 0.2g/cm 3 . This surfaces were obtained from gravity maps [9] for an area of 600x600 km near Ekaterinburg, Russia. The implementation for multiple GPU demonstrates an excellent scaling; the efficiency is more than 90% for eight GPUs.

Conclusion

For solving the inverse gravimetry problem of finding several boundary surfaces in multilayered medium, the parallel algorithm was constructed and implemented for multiple GPUs using the CUDA technology. The algorithm is based on the nonlinear conjugate gradient method and on the approach with weighting factors previously proposed by authors. The numerical implementation uses the modification based on approximation the Jacobian matrix by dropping out the less significant elements and to replacing the matrix by a block-band one.

The model problem of reconstructing three surfaces using the quasi-real gravitational data was solved on a large grid. The GPU implementation reduces the computation time by two orders of magnitude in comparison with CPU. The multi-GPU implementation demonstrates an excellent scaling and nearly 90% efficiency.