Three-dimensional Forward Modeling Of Time-domain Electromagnetic Response On GPUs

Time domain electromagnetic method has been widely used in resources and environment areas, such as exploration works for oil and gas, mineral resources, coal mines, geothermal and ground water resources, and monitoring and evaluation of engineering geology, natural and environmental disasters for its features of high resolution, strong ability of anti-noise and cost effectiveness comparing to other electromagnetic methods. It is also the preferred choice to use for predicting and evaluation of oil-bearing traps in exploration stage and for dynamic monitoring of developing reservoirs in petroleum industry. With the deepening of petroleum exploration, higher density sites for data acquizition and finer grids of large scale three dimensional models for data inversion are required to further improve the precision of data imaging and interpretation, resulting in a huge demand for computer memory capacity and long computation time. It is difficult to meet the computing needs by a single computer. Time domain finite difference (FDTD) scheme is adopted for forward modeling of time domain electromagnetic responses of three dimensional complex model, and parallel algorithm running on CPU/GPU platform are used to speedup large scale matrix calculations. This configuration of three dimensional forward modelling is fast enough as a useful tool for theoretical studies and for three dimensional inversions of large scale data sets and grids in time domain electromagnetic methods.Most commonly used methods of numerical simulation in electromagnetic responses are integral equation method, boundary element method, finite difference method and finite element method, and finite difference and finite element method are the main schemes for time domain EM response modelling. FDTD scheme based on staggered grids and leapfrog time advancing is widely used in computational electromagnetics for time domain electromagnetic response modelling since the scheme of staggered grid discretization was proposed by Yee in1966. In spatial domain, the Yee’s scheme divides the simulation region into rectangular meshs, the electric field is defined at the midpoint of the cell edges, and the magnetic field at the center of the cell surface, thus the spatial derivatives of electric and magnetic fields governed by Maxwell’s equations can be easily expressed by central difference. In time domain, electric and magnetic fields are calculated at interlaced moments and advanced gradually in leapfrog style, and finally all field components over the whole model space at various times are obtained by iteratively solving the Maxwell’s equations. This algorithm is simple and concise, but some constrains on grid size, time step and initial time must be meet to ensure the convergence and stability of solution. In order to simulate the responses of opening region, absorbing boundary conditions must be set for calculation of electromagnetic fields at truncated boundaries. There are various absorbing boundary conditions, such as the interpolation boundary, Mur absorbing boundary conditions, and perfectly matched layer (PML) which has beeb widely used. The effect of absorption of fields at boundaries is getting better with the increase in methods studied.The earth’s medium interested by geophysical exploration is a kind of lossy medium to electromagnetic fields. The electromagnetic signals for deep penetrating in earth are long period or low frequency variationsand can be treated as quasi-stationary fields. Based on the three-dimensional finite difference time domain electromagnetic simulation algorithm proposed by Tsili Wang etc, the electromagnetic response is decomposed into the primary field and the secondary field, and difference equations of the decomposed fields are derived for three dimensional model and the secondary fields governed by Maxwell’s equations are solved by FDTD scheme. Therefore the analytic solutions of field responses of homogeneous half space model for dipole source stimulating above ground, on surface and under ground have been derived, and then the primary fields in time domain are calculated at different point for different components by G-S transformation and Anderson digital filtering scheme for solving the Bessel integrals. For models with same background conductivity, primary fields for different source location can be obtained by simply interpolation of existing results, thus computer time can be greatly saved. Only recalculations of primary fields are needed for different type of source, such as steping source, Gauss pulse source, and the difference format of the secondary fields are the same and need no rederived for forward modeling and inversion. The adaptability to different type of sources is one of the advantages of this decomposed fields treatment.In dealing with the boundary conditions, different approaches are adopted for different boundaries. On the ground-air boundary, one grid is extended into the air, and the time domain magnetic fields on the top boundary are obtained by inverse Fourier transform the spectrum of field components on surface in wavenumber domain. Dirichlet boundary condition is imposed on the bottom boundary and the side boundaries since the source and anomalous bodies are assumed to be far enough away to these boundaries. The initial value of the secondary field can be set to zero for early enough initial time for this decomposed field approach. The initial time and time steps must meet the Courant stability condition. The secondary electric and magnetic fields with spatial and temporal distribution can be calculated by using the finite diffrernce equations of the secondary fields after setting the primary field distribution, initial condition and boundary conditions. The total field of the station for output is the summation of the primary field and the secondary field.Discretization of all computation area is required in FDTD algorithm, and thus large computer capacities in memories and CPU speed are expected. One of the key issues to practicability of the FDTD algorithm is to further improve the efficiency of computing. Graphics processor unit (GPU) is a new hardware acceleration technology arising recently and in a rapid growth trend, computing capability is different on different hardware configuration. The main advantages of GPU hardware acceleration include cost effective, surport Fortarn language, easy for development for cross platform. The FDTD algorithm based on CPU/GPU framework has been developed in this work through optimizing the algorithm and parallel programming. The feasibility and effectiveness of the GPU acceleration is illustrated by some sample model computations. The software environment for parallel algorithm development is Visual Stadio with PGI Fortran compiler, which not only support GPU and CUDA, but also OpenACC and OpenMP parallel instruction sets. The speedup ability of FDTD on GPU is tightly related with the hardware environment, and will be readily improved as the performance improved as hardware upgraded. The ability of upgrading and expansion can ensure this algorithm meet the need for high precision three dimensional inversion of large scale data set and model grids, and have broad application prospects.Various types of testing and validation have been carried out after implementation of this algoritm. The correctness of this algorithm is checked with the analytical solutions of homogeneous half space and numerical solutions of layer earth model, validated with calculation results by integral equation method for three dimensional model with conductive bodies, resistive bodies and body with covered layers. By comparison, it is shown that this algorithm is correct, efficient, and feasible for practical use in three dimensional forward modelling and inversion.