Three-dimension holographic numerical simulation of gravity anomalies

In gravity-anomaly-based prospecting, it is difficult to achieve large-scale and high-precision inversion imaging of complex geological models. To address this issue, this paper proposes a three-dimensional holographic numerical simulation method for gravity anomalies. This method transforms the 3D partial differential equation of gravity potential into many independent 1D differential equations with different wavenumbers by performing a 2D Fourier transform along the horizontal direction. It decomposes a large-scale 3D numerical modeling problem into many 1D numerical modeling problems, which greatly reduces the computation and memory requirements of modeling, and each 1D differential equation is independent, so it has high parallelism. The vertical direction is reserved as the spatial domain, which has strict upper and lower boundary conditions; The horizontal 2D Fourier transform uses a holographic Fourier transform (Holo-FT) with a full interval of integration, consistent with the physical boundary in the horizontal direction, thus the physical information of the gravity potential described by the three-dimensional partial differential equation is fully and accurately simulated. Using the high parallelism of the algorithm, the CPU is used for parallel solving ordinary differential equations, while the GPU is used for parallel computing of Holo-FT, which implements the CPU-GPU parallel acceleration scheme and further improves the efficiency of this algorithm.


Introduction
Gravity-based prospecting is widely used in the fields of oil and gas exploration, metal mineral exploration, deep Earth structure research, regional geological structure delineation and environmental geophysics.The study of efficient and high-precision three-dimensional numerical simulation of gravity anomalies under large-scale and complex conditions plays an important role in gravity anomaly inversion and interaction interpretation.
Numerical simulation methods of gravity anomalies can be divided into spatial and frequency domain methods.For the spatial domain, there are mainly analytical and numerical methods.Numerical simulation methods are mainly categorized into finite difference method (Farquharson and  Mosher, 2009)  [1] , finite volume method (Jahandari and Farquharson, 2013) [2] and finite element method (Jahandari and Farquharson, 2013) [2] .The spatial domain methods work well for models with simple anomalies, and the frequency domain approach has greater advantages for 3D numerical simulation of gravity fields under complex conditions (Chai, 1988  [3] ; Parker, 2010 [4] ; Wu & Tian,

A 3-D holographic numerical simulation method of gravity anomalies
The gravity potential satisfying the Poisson equation can be described as 2 ( , , ) 4 ( , , ) ) where ( , , ) U x y z represents the gravity potential in spatial, G is the Newtonian gravitational constant, ( , , ) x y z


is the residual density.Equation( 1) was converted into the spatial-wavenumber mixed domain by a 2D Fourier transform in the x-and y-directions.From this, we obtained where ( ) boundary conditions(3) and ( 4) are consistent with the gravity physical boundary conditions, which accurate represent the upper and lower boundary of gravity anomalies in the mixed domain.The finite element method is used to solve the equations to obtain the gravity potential in the mixed domain.The gravity field can be obtained from the relationship between the gravity potential and the field in the mixed domain.The final step is to calculate the gravity fields in the spatial domain, which is performed using Holo-FFT.

Holo-FT
The two-dimensional Fourier transform is defined as follows: () ( , ) ( , ) where ( , ) f x y is the field in the spatial domain, ( , ) xy F k k represents the two-dimensional spectrum in space-wavenumber domain.
Let's take the two-dimensional Fourier transform(5) as an example to explain the computation process of Holo-FT based on quadratic interpolation shape function.We transform Equation( 5) into two one-dimensional integrals, the one-dimensional Fourier transform integration interval is ( , ) − + , thus all information in the horizontal direction is preserved.We discretize the calculation domain using either uniform or non-uniform grids along the x and y directions, dividing it into Mx×My elements.Within each element, the function change in the space domain is approximated using quadratic interpolation shape function.The analytical expression for the element integral is derived to obtain the integration coefficients, and the final result is expressed as follows: , , When the grid divisions along the x-and y-directions are fixed, the shape function remain the same in the element.The integration coefficients along both directions can be pre-calculated and stored, allowing for reuse and reducing redundant calculations, thus enhancing computational efficiency.Moreover, the integration coefficients are independent for different wavenumbers.The derivation process for the two-dimensional inverse Fourier transform( 6) is similar to the forward transform and would not be elaborated here.
The vertical direction is retained as spatial domain with strict upper and lower boundary conditions, which is consistent with the physical upper and lower boundaries; the two-dimensional Fourier forward and inverse transforms are adopted as Holo-FT, and all the information is retained, which is consistent with the physical boundaries in the horizontal direction, and the physical information of the gravity potential is completely and accurately simulated, thus this method is called the threedimensional holographic numerical simulation of gravity anomalies.
In the algorithm of three-dimensional holographic numerical simulation of gravity anomalies, the computation of different wavenumbers involves independent integration coefficients and independent solving of five-diagonal equation systems.Both cores exhibit good independence.Therefore, we can use OpenMP to parallelize the solving of the five-diagonal equation systems and GPU to accelerate Holo-FT, achieving the acceleration scheme of CPU-GPU.

Numerical example
In order to verify the correctness and high efficiency of the proposed algorithm, different models were designed to carry out forward calculation.The following operations are implemented using Fortran95 programming language.All examples are carried out on a computer with CPU Intel(R) Core(TM) i9-9980XE, 128 GB RAM, NVIDIA TITAN RTX.
The model is shown in Figure 1, and the simulation area is 20 km (x) × 20 km (y) × 7.5 km (z).The computational region has one prism anomaly located in the middle and four prism anomalies centrally symmetric about the origin, all of which have the size of 2 km × 2 km × 2 km and are buried at the top surface at a depth of 0.5 km, four anomalies symmetric about the origin are 200 m from the boundary, and the residual density of the anomaly at the middle of the area is , and that of the other four anomalies is . The number of nodes is 251(x)×251(y)×151(z), and the grid is uniformly divided, with the grid spacing of 80 m in the horizontal direction and 50 m in the vertical direction.According to the sampling theorem, the range of wavenumber of x k and y k are both 22 3.9 10 ~3.9 10 , the number of wavenumber are selected as follows: 217 ( x k )×217( y k ).Comparing of numerical simulation results with analytical solutions, Figure 2 shows the contours of the numerical and analytical solutions for z = 0 m on the ground.The contours of the numerical and analytical solution overlap and the field match well, and the relative root-mean-square error (RMES) of three components in the z = 0 m plane are 0.057%, 0.057%, and 0.058%, respectively, indicating that the 3D holographic numerical simulation is correct.

Figure 2. The contours of numerical and analytical solutions of gravity field
In order to compare the effect of CPU-GPU parallel acceleration, the number of nodes in the zdirection is kept as 151, and the number of computational nodes in the horizontal direction is changed, and the number of waves is kept consistent with the number of nodes in the spatial domain.The calculation results are shown in the following table 1, where the acceleration ratio refers to the serial time consumed/parallel computation time consumed.From the table, with the increase of computational nodes, the acceleration ratio of CPU-GPU is getting bigger and bigger, and it reaches 191.63 at the highest, which is because the GPU has many computational units, and the computational scale has not yet reached the upper limit of GPU computation.

Conclusion
In this paper, a three-dimensional holographic numerical simulation algorithm of gravity anomalies is realized, using the holographic Fourier transform, all the physical information in the horizontal direction is preserved, the boundary is consistent with the physical field boundary, the threedimensional problem is transformed into a one-dimensional ordinary differential problem solving, the computation and storage requirements are reduced, the vertical is preserved as spatial domain, the upper and lower boundary conditions are consistent with the physical upper and lower boundaries, and the physical information of the gravity potential is completely and accurately simulated.We realize the CPU-GPU parallel scheme.The algorithm is applied to the model, and the results show that the algorithm has high computational accuracy and efficiency, and is suitable for large-scale model numerical simulation.The algorithm is also applied to the numerical simulation of direct current and electromagnetic fields.

Table 1 .
The efficiency of different schemes