CN105467443A - A three-dimensional anisotropy elastic wave numerical simulation method and system - Google Patents

Abstract

The invention relates to a three-dimensional anisotropy elastic wave numerical simulation method and system. The method comprises the following steps: 1, a medium model is established to carry out grid discretization on the medium model to obtain a plurality of grid points; 2, a seismic source function is calculated, and a pressure value at each grid point is calculated according to the seismic source function; 3, a three-dimensional anisotropy elastic wave equation is converted into a propagation equation, and the pressure value at each grid point is taken into the propagation equation as substitutes to conduct calculating to obtain a wave field value at each moment; and 4, according to the wave field value, a calculating area of each grid point is determined; area division is carried out; and data exchange is carried out on divided area boundary data to complete elastic wave numerical simulation. According to the invention, acceleration of three-dimensional elastic wave numerical simulation under a complex medium is realized through utilization of a GPU; a realization scheme of accelerating data transmission through utilization of GPU Direct technology is realized; a large amount of data copying from a CPU to the GPU and from the GPU to the CPU is avoided; and communication bottleneck optimization is realized.

Description

A kind of three dimensional anisotropic elastic-wave numerical modeling method and system

Technical field

The present invention relates to a kind of three dimensional anisotropic elastic-wave numerical modeling method and system, belong to field of geophysical exploration.

Background technology

Seismic wave numerical modeling, based on the communication theory of seismic event in underground medium, is widely used in prospecting seismology and earthquake seismology.Three-dimensional acoustic wave equation, the elastic wave isotropy numerical simulation of current routine have been widely used in each geophysics fields such as numerical simulation, imaging and inverting.But very large challenge is still existed for the three dimensional elasticity wave field numerical simulation study how effectively realizing large scale acquisition mode (such as wide-azimuth collection) and complex anisotropic medium (such as horizontal cross isotropy HTI or orthotropy), is not widely used in actual applications.In addition conventional extensive three-dimensional acoustic wave, elastic-wave numerical modeling based on CPU, need a large amount of special PC cluster resource usually.From hardware cost or its cost is all very expensive calculating energy consumption.Therefore improving calculated performance, significantly reduction fast assesses the cost significant for the practical application realizing the simulation of three dimensional elasticity ripple anisotropy values.

In nearly ten years, the acceleration utilizing Graphics Processing Unit (GPU) to carry out compute-intensive applications realizes the development having obtained formula of advancing by leaps and bounds.Graphics Processing Unit (GPU) has memory bandwidth at a high speed due to it, the computing core of two orders of magnitude is at least exceeded compared to CPU, be more suitable for single instruction multiple data (SIMD) computation schema of parallel computation, and lower energy consumption cost, be applied to the association area of computational science just widely.For exploration geophysics field, for using the interest of Graphics Processing Unit (GPU) also in remarkable enhancing, GPU is used for the core algorithm accelerated in seismic processing by increasing research, such as earthquake numerical simulation, seismic imaging, the inverting of earthquake high precision etc.

Summary of the invention

Find through research, prior art and there is following problem:

Time finite element method method is applied to the research that GPU equipment realizes on complex dielectrics wave propagation algorithm little, especially for how to process there is a large amount of multinode exchanges data extensive three-dimensional problem research just less (Heinrichetal, 2014).

Some GPU seismic wave propagations of current existence simulation implementation, for three-dimensional GPU cluster implementation we be described for the scheme of (Long Guihua etc.).The program utilizes Region Decomposition technology that the geologic body that single GPU can not calculate is carried out coarseness decomposition along Z-direction, adopt message passing interface (MPI) to exchange data boundary, thus use the mode of MPI+CUDA to achieve the numerical simulation of large scale 3-D seismics wave field.But there is a very large problem in the method, the speed-up ratio adopting GPU cluster to calculate significantly declines compared to single GPU and CPU, causes the reason of this result to be in the major part consumption consuming time of the GPU calculating data copy between internodal GPU to CPU and CPU to the GPU equipment of GPU cluster.

Technical matters to be solved by this invention is, for the deficiencies in the prior art, there is provided a kind of three dimensional anisotropic elastic-wave numerical modeling method and system based on GPUDirect optimize communicate GPUDirect technology and finite difference numerical simulation methods combining obtained, for numerical simulation provides strong guarantee.

The technical scheme that the present invention solves the problems of the technologies described above is as follows: a kind of three dimensional anisotropic elastic-wave numerical modeling method, specifically comprises the following steps:

Step 1: set up dielectric model, carries out to dielectric model that grid is discrete obtains multiple net point;

Step 2: calculate source function, calculate the force value on each net point according to source function;

Step 3: three dimensional anisotropic equations for elastic waves is converted to propagation equation, brings the force value on each net point into propagation equation and calculates, and obtains the wave field value at per a moment;

Step 4: the zoning determining each net point according to wave field value, carries out subregion and carries out exchanges data to partition boundaries data, completing elastic-wave numerical modeling.

The invention has the beneficial effects as follows: from ess-strain equation, achieve the three dimensional elasticity ripple numerical simulation under utilizing Graphics Processing Unit (GPU) to accelerate complex dielectrics, and for three-dimensional problem introduce GPU equipment faced by the multinode caused due to Region Decomposition between, intra-node communication bottleneck problem conducts in-depth research and analyzes, propose the implementation utilizing GPUDirect technology expedited data to transmit, avoid a large amount of CPU to GPU, the data copy of GPU to CPU, achieves the problem of optimize communicate bottleneck.GPUDirect communication optimization strategy is introduced by GPU figure acceleration equipment, overall computational performance can be promoted significantly, the three dimensional anisotropic elasticity numerical simulation that can realize with lower hardware cost and less time, for the various algorithm depending on Wave equation forward modeling is as reverse-time migration, the application of full waveform inversion provides strong guarantee.

On the basis of technique scheme, the present invention can also do following improvement.

Further, described step 4 specifically comprises the following steps:

Step 4.1: the zoning determining each net point according to all wave field value, and adopt the mode of absorbing boundary to determine border, zoning, in zoning, simulate wave propagation in underground medium;

Step 4.2: carry out subregion to all zonings, carries out exchanges data to the data boundary of adjacent sectors, obtains simulation elasticity wave datum, completes elastic-wave numerical modeling.

Further, in described step 4.2, the employing GPU-Direct technology of data boundary carries out exchanges data.

Adopt the beneficial effect of above-mentioned further scheme to be utilize GPU to accelerate three dimensional elasticity ripple numerical simulation, significantly improve computing velocity, realize compared to conventional CPU cluster, 20-30 speed doubly can be reached and promote.

Further, the exchanges data in described step 4.2 specifically comprises:

Carry out Region Decomposition along the axis calculating zone boundary data variation in dielectric model the slowest, obtain multiple subregion, be assigned to by all subregions in multiple GPU, each subregion is independent performs calculating on a GPU; The data boundary of every two adjacent sectors exchanges again.

The beneficial effect of above-mentioned further scheme is adopted to be, GPUDirect technology is utilized to accelerate the implementation of data transmission between node, in node, our scheme because which obviating a large amount of CPU to GPU, the data copy of GPU to CPU, thus the problem achieving optimize communicate bottleneck.

Further, described step 1 specifically comprises the following steps:

Step 1.1: set up dielectric model according to geology background condition, the actual physical test of rock data that record and Well logging Data;

Step 1.2: adopt the 3D grid of regular shape to carry out grid to dielectric model discrete, obtain net point.

The beneficial effect adopting above-mentioned further scheme is that dielectric model is separated into net point, and net point, as the carrier of wave field value, forms wave field by multiple net points of three-dimensional structure, as the starting condition that wave field is propagated.

Further, described source function spatially adopts Gaussian function, and adopt Ricker wavelet in time, the concrete formula of described source function is:

S (x, y, z, t)=g (x, y, z) f (t) formula (1)

Wherein, f (t)=(1-2 (π f ₀t) ²) exp (-(π f ₀t) ²) formula (2)

$g(x,y,z)=\frac{\exp \left\{\right[-{(x-x_0)}^2+{(z-z_o)}^2+{(y-y_o)}^2\left]\right\}}{{\mathrm{\β}}^2}$ Formula (3)

In formula: t represents the time, f ₀represent the centre frequency of Ricker wavelet, f during model calculates ₀=15Hz, β are constant; (x ₀, y ₀, z ₀) be the locus of focus, x, y and z are respectively the position on x-axis, y-axis and z-axis direction.

Further, described step 3 specifically comprises the following steps:

Step 3.1: the differential difference approximation in three dimensional anisotropic equations for elastic waves is substituted, obtains the propagation equation of corresponding finite difference scheme, in described propagation equation, space sampling step length and time-sampling step-length meet the stability condition of this numeric format;

Step 3.2: adopt the mode of Region Decomposition to bring the force value on each net point into propagation equation and calculate, obtain the wave field value at per a moment.

The technical scheme that the present invention solves the problems of the technologies described above is as follows: a kind of three dimensional anisotropic elastic-wave numerical modeling system, comprises MBM, focus module, propagation module and data exchange module;

Described MBM is used for setting up dielectric model, dielectric model is carried out to grid is discrete obtains multiple net point;

Described focus module, for calculating source function, calculates the force value on each net point according to source function;

Described propagation module is used for three dimensional anisotropic equations for elastic waves to be converted to propagation equation, brings the force value on each net point into propagation equation and calculates, obtain the wave field value at per a moment;

Described data exchange module is used for the zoning determining each net point according to wave field value, carries out subregion and carries out exchanges data to partition boundaries data, completing elastic-wave numerical modeling.

The invention has the beneficial effects as follows: from ess-strain equation, achieve the three dimensional elasticity ripple numerical simulation under utilizing Graphics Processing Unit (GPU) to accelerate complex dielectrics, and for three-dimensional problem introduce GPU equipment faced by the multinode caused due to Region Decomposition between, intra-node communication bottleneck problem conducts in-depth research and analyzes, propose the implementation utilizing GPUDirect technology expedited data to transmit, avoid a large amount of CPU to GPU, the data copy of GPU to CPU, achieves the problem of optimize communicate bottleneck.GPUDirect communication optimization strategy is introduced by GPU figure acceleration equipment, overall computational performance can be promoted significantly, the three dimensional anisotropic elasticity numerical simulation that can realize with lower hardware cost and less time, for the various algorithm depending on Wave equation forward modeling is as reverse-time migration, the application of full waveform inversion provides strong guarantee.

On the basis of technique scheme, the present invention can also do following improvement.

Further, described data exchange module comprises area determination module and subregion Switching Module;

Described area determination module is used for the zoning determining each net point according to all wave field value, and adopts the mode of absorbing boundary to determine border, zoning, simulates wave propagation in underground medium in zoning;

Described subregion Switching Module is used for carrying out subregion to all zonings, carries out exchanges data to the data boundary of adjacent sectors, obtains simulation elasticity wave datum, completes elastic-wave numerical modeling.

Further, in described subregion Switching Module, the employing GPU-Direct technology of data boundary carries out exchanges data.

The beneficial effect of above-mentioned further scheme is adopted to be, GPU is utilized to accelerate three dimensional elasticity ripple numerical simulation, significantly improve computing velocity, realize compared to conventional CPU cluster, the speed that can reach 20 ~ 30 times promotes, and utilize GPUDirect technology to accelerate the implementation of data transmission between node, in node, our scheme is because which obviating a large amount of CPU to GPU, the data copy of GPU to CPU, thus the problem achieving optimize communicate bottleneck.

Accompanying drawing explanation

Fig. 1 is a kind of three dimensional anisotropic elastic-wave numerical modeling method flow diagram described in the embodiment of the present invention 1;

Fig. 2 utilizes the first order derivative of eight rank limited precision difference approximation central points to need 25 grid point diagrams in the embodiment of the present invention;

Fig. 3 is the implementing procedure figure of three dimensional anisotropic elastic-wave numerical modeling on single GPU equipment in the embodiment of the present invention;

Fig. 4 is data volume Region Decomposition and signal intelligence general introduction figure in the embodiment of the present invention;

Fig. 5 a adopts MPI data transfer mode schematic diagram in prior art interior joint;

Fig. 5 b adopts GPUDirectP2P data transfer mode schematic diagram in the embodiment of the present invention;

Fig. 6 a adopts MPI data to transmit schematic diagram between prior art interior joint;

Fig. 6 b adopts GPUDirectRDAM mode to carry out data transmission schematic diagram in the embodiment of the present invention;

Fig. 7 is a kind of three dimensional anisotropic elastic-wave numerical modeling system architecture diagram described in the embodiment of the present invention 1;

Fig. 8 a-8c has different TI symmetric elastic medium stiffness coefficient matrix schematic diagram for three kinds of testing in the embodiment of the present invention;

Fig. 9 a-9c is the three-dimensional impulse response schematic diagram that in the embodiment of the present invention, three kinds have a symmetric elastic medium of different TI;

Figure 10 is the Measurement results comparison diagram of conventional MPI communication mode and GPUDirect communication mode in embodiment of the present invention interior joint;

Figure 11 adopts conventional MPI communication mode and the test result comparison diagram adopting GPUDirectRDMA communication mode between embodiment of the present invention interior joint.

In accompanying drawing, the list of parts representated by each label is as follows:

1, MBM, 2, focus module, 3, propagation module, 4, data exchange module.

Embodiment

Be described principle of the present invention and feature below in conjunction with accompanying drawing, example, only for explaining the present invention, is not intended to limit scope of the present invention.

As shown in Figure 1, a kind of three dimensional anisotropic elastic-wave numerical modeling method described in the embodiment of the present invention 1, specifically comprises the following steps:

Step 1: set up dielectric model, carries out to dielectric model that grid is discrete obtains multiple net point;

Step 2: calculate source function, calculate the force value on each net point according to source function;

Step 3: three dimensional anisotropic equations for elastic waves is converted to propagation equation, brings the force value on each net point into propagation equation and calculates, and obtains the wave field value at per a moment;

Step 4: the zoning determining each net point according to wave field value, carries out subregion and carries out exchanges data to partition boundaries data, completing elastic-wave numerical modeling.

In embodiment 2, on the basis of embodiment 1, described step 4 specifically comprises the following steps:

Step 4.1: the zoning determining each net point according to all wave field value, and adopt the mode of absorbing boundary to determine border, zoning, in zoning, simulate wave propagation in underground medium;

Step 4.2: carry out subregion to all zonings, carries out exchanges data to the data boundary of adjacent sectors, obtains simulation elasticity wave datum, completes elastic-wave numerical modeling.

In embodiment 3, on the basis of embodiment 1 or 2, in described step 4.2, the employing GPU-Direct technology of data boundary carries out exchanges data.

In embodiment 4, on the basis of embodiment 1-3 any embodiment, the exchanges data in described step 4.2 specifically comprises:

Carry out Region Decomposition along the axis calculating zone boundary data variation in dielectric model the slowest, obtain multiple subregion, be assigned to by all subregions in multiple GPU, each subregion is independent performs calculating on a GPU; The data boundary of every two adjacent sectors exchanges again.

In this embodiment, further, carry out exchanges data to the computation bound data acquisition GPU-Direct technology of each computing node in each moment to be implemented as follows:

When process solves extensive three-dimensional model, limited global memory makes single GPU equipment cannot store whole model meshes.For this reason, our algorithm by CUDA code extensions to operating in many GPU equipment, in the isomeric architecture of multinode, must carry out high-performance calculation to make full use of GPU computing power.Based on this reason, we adopt a Region Decomposition scheme, our Region Decomposition scheme as shown in Figure 3, Region Decomposition is carried out along the slowest Y direction of computing grid change, thus whole zoning is decomposed into N/M part, N is the lattice number in Y-direction, and M is the GPU number of current computing cluster.Each subregion is independent performs calculation procedure 1-3 at a GPU equipment.

Because finite difference method (as shown in Figure 2, for 8 rank precision) need on y direction forward and backward four points, therefore the GPU equipment carrying out finite difference operation in every sub regions needs to obtain the value with its zoning consecutive point on other equipment, therefore the data boundary in adjacent subarea territory must exchange on each time step, its data transfer path as shown in Figure 3, namely be divided in node and two types between node, as shown in Figure 4, data volume Region Decomposition and signal intelligence general introduction, term explanation in figure, PCI-E, it is the abbreviation of PCI-Express, refer to computer bus and interface standard.InfiniBand, it is a kind of multi-concurrency network framework being specifically designed to server-side network and connecting, GPUDirectP2P, implementation (being called the direct-connected P-2-P technology of GPU) in the node of GPUDirect, GPUDirectRDMA, RDMA technology full name is remote direct data access, implementation (becoming the long-range direct-connecting technology of GPU) between the node that GPUDirectRDMA refers to GPUDirect technology.Utilize GPUDirectP2P to communicate in node, between node, utilize GPUDirectRDAM to communicate.Usually, the communication of this different GPU equipment room needs the data copy of carrying out GPU to CPU and CPU to GPU due to data, and whole communication process is very consuming time, and this is the Main Bottleneck that conventional GPU cluster calculates.

To this, the present invention utilizes GPUDirect to optimize these data communication, and it is specifically implemented as follows:

A) ensure that CUDA version is on 6.0, and video card computing power is on 2.5, to ensure the support to GPUDirect.

B) CUDA-AwareMPI is adopted, CUDA-AwareMPI supports that the MPI for GPU calculation optimization of GPUDirect characteristic realizes, present stage supports that the MPI of CUDA-AwareMPI realizes having, and more than MAVPCH21.9 version or more than OpenMPI1.7 version, we are for MAVPICH2.

C) cudaSetDevice function needs to call before MPI_Init, and utilizes the appointment of MPI environmental variance realization to GPU, MVAPICH2:MV2_COMM_WORLD_LOCAL_RANK.

D) MV2_USE_CUDA needs to add when performing mpirun.

E) adopt GPUDirectP2P mode (GPU is direct-connected point-to-point) for the exchanges data in node, can directly utilize cudaMemcoycopy to carry out the data transmission of GPU to GPU after completing above-mentioned configuration, it specifically describes see Fig. 5.

F) exchanges data for cross-node adopts (the RDMA technology full name remote direct data access of GPUDirectRDMA mode, produce to solve the delay of servers' data process in Internet Transmission, GPU is direct-connected is referred to as GPUDirectRDMA to utilize Infiniband network directly to carry out), MPI_Isend and MPI_Irecv can be directly adopted to send and accepting device internal memory after completing above-mentioned configuration, but not need first device memory to be explicit cudaMemcpy to host memory again by MPI_Isend and MPI_Irecv transmission and reception as in the past, again by explicit cudaMemcpy to equipment end, communicate, two kinds of visible Fig. 6 of communication difference.

In embodiment 5, on the basis of embodiment 1-4 any embodiment, described step 1 specifically comprises the following steps:

Step 1.1: set up dielectric model according to geology background condition, the actual physical test of rock data that record and Well logging Data;

Step 1.2: adopt the 3D grid of regular shape to carry out grid to dielectric model discrete, obtain net point.

In embodiment 6, on the basis of embodiment 1-5 any embodiment, described source function spatially adopts Gaussian function, and adopt Ricker wavelet in time, the concrete formula of described source function is:

S (x, y, z, t)=g (x, y, z) f (t) formula (1)

Wherein, f (t)=(1-2 (π f ₀t) ²) exp (-(π f ₀t) ²) formula (2)

$g(x,y,z)=\frac{\exp \left\{\right[-{(x-x_0)}^2+{(z-z_o)}^2+{(y-y_o)}^2\left]\right\}}{{\mathrm{\β}}^2}$ Formula (3)

In formula: t represents the time, f ₀represent the centre frequency of Ricker wavelet, f during model calculates ₀=15Hz, β are constant; (x ₀, y ₀, z ₀) be the locus of focus, x, y and z are respectively the position on x-axis, y-axis and z-axis direction.

In embodiment 6, on the basis of embodiment 1-5 any embodiment, described step 3 specifically comprises the following steps:

Step 3.1: the differential difference approximation in three dimensional anisotropic equations for elastic waves is substituted, obtains the propagation equation of corresponding finite difference scheme, in described propagation equation, space sampling step length and time-sampling step-length meet the stability condition of this numeric format;

Step 3.2: adopt the mode of Region Decomposition to bring the force value on each net point into propagation equation and calculate, obtain the wave field value at per a moment.

In this embodiment, further, adopt the mode of Region Decomposition to utilize GPU to enter finite difference formulations at each computing node to three dimensional anisotropic equations for elastic waves in described step 3.2, the concrete grammar obtaining the wave field value at per a moment is as follows:

First, provide the mathematical description of the elastic wave wave equation under 3-D elastic anisotropic media, its form of three dimensional anisotropic Time Migration of Elastic Wave Equation is:

$\mathrm{\ρ}{\∂}_{\mathrm{tt}}^2{u}_i={\∂}_j{\mathrm{\σ}}_{\mathrm{ij}}+F_i,$ Formula (4)

In formula, u is wave field function; σ is cauchy stress tensor; F _ibe the physical item of unit volume, it can be equivalent to stress riser, and ρ is medium material density, represent second-order time local derviation.Linear elasticity theory is utilized to obtain for stress tensor:

σ _ij=c _ijklε _kl, formula (5)

ε in formula _klfor linear strain tensor, c _ijklfor the elastic stiffness matrix of the fourth-order tenstor, this parameter needs to provide according to the actual conditions of medium; And for ε _klwe have

${\mathrm{\ϵ}}_{\mathrm{kl}}=\frac{1}{2}[{\∂}_k{u}_l+{\∂}_l{u}_k],k,l=\mathrm{1,2,3},$ Formula (6)

Wherein ε _klrepresent linear strain tensor, represent and space derivative is asked, u to kth dimension _lrepresent the displacement of lth wave field.According to this formula, we just can the numerical simulation of row Time Migration of Elastic Wave Equation, and the concrete calculation procedure of the elastic-wave numerical modeling utilizing GPU to carry out is as follows:

3.2.1 we define a N _x× N _y× N _zcomputing grid, and specify N _tindividual time step, just can utilize a four-tuple to represent at each net point of room and time: namely like this, [x, y, z|t]=[p Δ x, q Δ y, r Δ z|n Δ t].Therefore discrete net point just can be utilized to represent at the continuous print wave field displacement record of specified point:

U _i| _{x, y, z|t}≈ u _i ^{p, r, q|n}formula (7)

We are similar to single order local derviation by a following M rank precision center difference operator D _x[], can obtain:

${\∂}_x{u}_j\≈D_x\left[u_j^{p,q,r|n}\right]=\frac{1}{\mathrm{\Δx}}\underset{\mathrm{\α}=1}{\overset{M}{\mathrm{\Σ}}}{W}_{\mathrm{\α}}(u_j^{p+\mathrm{\α},q,r|n}-u_j^{p-\mathrm{\α},q,r|n})$ Formula (8)

W in formula _αbe a weight coefficient polynomial expression, i.e. finite difference coefficient, M is finite difference formulations exponent number;

Space difference operator D _y[] and D _z[] and D _x[] similar partial derivative representing y and z direction respectively.This three's difference operators all are substituted into equation (6), in order to solve the partial derivative in constitutive equation.

3.2.2 for time discrete, we use the space second order accuracy of a standard to be similar in accounting equation (6)

${\∂}_{\mathrm{tt}}^2{u}_j\≈D_{\mathrm{tt}}{u}_j^{p,q,r|n}=\frac{1}{\mathrm{\Δ}{t}^2}[u_j^{p,q,r|n+1}-2u_j^{p,q,r|n}+u_j^{p,q,r|n-1}]$ Formula (9)

Above-mentioned difference operator is substituted into equation (4) and by rearranging these discrete items, the mode that we just can be increased by time step is carried out unknown wave field and solved.Such as we need its wave field value current and before given and x-, y-, z-direction consecutive point required for masterplate are calculated according to finite difference formulations space.Fig. 3 describe limited space Difference Calculation masterplate at current time step about point required computing grid point.

As shown in Figure 7, be a kind of three dimensional anisotropic elastic-wave numerical modeling system described in the embodiment of the present invention 1, comprise MBM 1, focus module 2, propagation module 3 and data exchange module 4;

Described MBM 1, for setting up dielectric model, carries out to dielectric model that grid is discrete obtains multiple net point;

Described focus module 2, for calculating source function, calculates the force value on each net point according to source function;

Described propagation module 3, for three dimensional anisotropic equations for elastic waves is converted to propagation equation, is brought the force value on each net point into propagation equation and is calculated, and obtains the wave field value at per a moment;

Described data exchange module 4, for determining the zoning of each net point according to wave field value, carries out subregion and carries out exchanges data to partition boundaries data, completing elastic-wave numerical modeling.

In embodiment 2, on the basis of embodiment 1, described data exchange module comprises area determination module and subregion Switching Module;

Described area determination module is used for the zoning determining each net point according to all wave field value, and adopts the mode of absorbing boundary to determine border, zoning, simulates wave propagation in underground medium in zoning;

Described subregion Switching Module is used for carrying out subregion to all zonings, carries out exchanges data to the data boundary of adjacent sectors, obtains simulation elasticity wave datum, completes elastic-wave numerical modeling.

In embodiment 3, on the basis of embodiment 2, in described subregion Switching Module, the employing GPU-Direct technology of data boundary carries out exchanges data.

Basis of the present invention is communication theory and the high-performance GPU method for numerical simulation thereof of seismic event, the present invention is based on the GPUDirect communication technology and achieve high performance 3-D elastic anisotropic media elastic-wave numerical modeling method, in concrete example, implementation step is respectively:

1. dielectric model is set up:

Using the three dimensional elasticity medium of different TI symmetry (isotropy, VTI, HTI) as test model, its elastic parameter is as follows: identical isotropy parameter, v _p=2.0km/s, and ρ=2000kg/m ³, but use different Thomsen anisotropic parameterses.Make VTI medium, [ε ₁, δ ₁, γ ₁]=[0.2 ,-0.1,0.2], HTI medium, [ε ₂, δ ₂, γ ₂]=[0.2 ,-0.1,0.2].These parameters are all passed through corresponding transformation for mula (Thomsen, 1986) and are changed in stiffness matrix.Fig. 8 a-c represents the stiffness matrix C of ISO, VTI, HTI medium 6 × 6 for impulse response respectively _ij, represent different values by different color block.

2. model discretize:

Carry out grid to the dielectric model of design discrete, test geological data dimension size is N _x× N _y× N _z=200 ³, mesh spacing is Δ x=Δ y=Δ z=0.005km.

3. source function:

Source function spatially adopts Gaussian function, and the time adopts Ricker wavelet, and form is as follows:

S (x, y, z, t)=g (x, y, z) f (t) formula (1)

Wherein, f (t)=(1-2 (π f ₀t) ²) exp (-(π f ₀t) ²) formula (2)

$g(x,y,z)=\frac{\exp \left\{\right[-{(x-x_0)}^2+{(z-z_o)}^2+{(y-y_o)}^2\left]\right\}}{{\mathrm{\β}}^2}$ Formula (3)

In formula: t represents the time, f ₀represent the centre frequency of Ricker wavelet, f during model calculates ₀=15Hz, β are constant; (x ₀, y ₀, z ₀) be the locus of focus, x, y and z are respectively the position on x-axis, y-axis and z-axis direction.

4. the numeric format of Time Migration of Elastic Wave Equation:

Adopt 8 rank precision to carry out finite difference numerical simulation, then corresponding the wave field that we utilize formula (4)-(6) to carry out each step calculates.

5. Region Decomposition and GPUDirect communication:

Adopt aforesaid Region Decomposition scheme, carry out Region Decomposition in the Y direction, adopt GPUDirectP2P to communicate in node, between node, adopt GPUDirectRDMA to communicate.Test environment is for using the machine to be that three two-way GPU service nodes are as test platform.Each Joint Enterprise has two-way 10 core XeonE5-2680V22.80GHz processor, 128GBDDR3 internal memory.GPU equipment is NVIDIATeslaK40c (Kepler), and Infiniband adapter is MellanoxConnectX-3IBQDRMT26428, and all nodes operate in RHEL6.3, and the driving version of Infiniband is OFED-2.3-1.0.1.All tests are based on the MPI of MVAPICH2-2.1a release version.

6. numerical simulation result analysis:

Fig. 9 is the wave field snapshot that three test models obtain, and can be clearly seen that the impact that anisotropy causes wave field.Figure 10 is the communication test in the GPUDirectP2P node that carries out, ordinate in figure is speed-up ratio, obtained computing time by computing time of GPU and single CPU, can be clearly seen that, in the node that have employed GPUDirectP2P, data transmission has significant counting yield to promote compared to traditional MPI method.Figure 11 is the internodal communication test of GPUDirectRDMA carried out, the meaning of what the MPI-GDR wherein shown in legend represented is GPUDirectRDMA, carried out the test of two platforms, each node of platform A 1 piece of K40C totally three nodes, each node of platform B 3 pieces of K40C are three phases altogether.Test result shows that the communication of GPUDirectRDMA to cross-node has and promotes effect significantly, and along with its effect of increase of amount of communications obvious all the more.

The foregoing is only preferred embodiment of the present invention, not in order to limit the present invention, within the spirit and principles in the present invention all, any amendment done, equivalent replacement, improvement etc., all should be included within protection scope of the present invention.

Claims (10)

1. a three dimensional anisotropic elastic-wave numerical modeling method, is characterized in that, specifically comprises the following steps:

Step 1: set up dielectric model, carries out to dielectric model that grid is discrete obtains multiple net point;

Step 2: calculate source function, calculate the force value on each net point according to source function;

Step 3: three dimensional anisotropic equations for elastic waves is converted to propagation equation, brings the force value on each net point into propagation equation and calculates, and obtains the wave field value at per a moment;

Step 4: the zoning determining each net point according to wave field value, carries out subregion and carries out exchanges data to partition boundaries data, completing elastic-wave numerical modeling.

2. a kind of three dimensional anisotropic elastic-wave numerical modeling method according to claim 1, it is characterized in that, described step 4 specifically comprises the following steps:

Step 4.1: the zoning determining each net point according to all wave field value, and adopt the mode of absorbing boundary to determine border, zoning, in zoning, simulate wave propagation in underground medium;

Step 4.2: carry out subregion to all zonings, carries out exchanges data to the data boundary of adjacent sectors, obtains simulation elasticity wave datum, completes elastic-wave numerical modeling.

3. a kind of three dimensional anisotropic elastic-wave numerical modeling method according to claim 2, is characterized in that, in described step 4.2, the employing GPU-Direct technology of data boundary carries out exchanges data.

4. a kind of three dimensional anisotropic elastic-wave numerical modeling method according to claim 3, it is characterized in that, the exchanges data in described step 4.2 specifically comprises:

Carry out Region Decomposition along the axis calculating zone boundary data variation in dielectric model the slowest, obtain multiple subregion, be assigned to by all subregions in multiple GPU, each subregion is independent performs calculating on a GPU; The data boundary of every two adjacent sectors exchanges again.

5. a kind of three dimensional anisotropic elastic-wave numerical modeling method according to any one of claim 1-4, it is characterized in that, described step 1 specifically comprises the following steps:

Step 1.1: set up dielectric model according to geology background condition, the actual physical test of rock data that record and Well logging Data;

Step 1.2: adopt the 3D grid of regular shape to carry out grid to dielectric model discrete, obtain net point.

6. a kind of three dimensional anisotropic elastic-wave numerical modeling method according to claim 5, it is characterized in that, described source function spatially adopts Gaussian function, and adopt Ricker wavelet in time, the concrete formula of described source function is:

S (x, y, z, t)=g (x, y, z) f (t) formula (1)

Wherein, f (t)=(1-2 (π f ₀t) ²) exp (-(π f ₀t) ²) formula (2)

$g(x,y,z)=\frac{\exp \{\[-{(x-x_0)}^2+{(z-z_o)}^2+{(y-y_o)}^2\]\}}{{\mathrm{\β}}^2}$ Formula (3)

In formula: t represents the time, f ₀represent the centre frequency of Ricker wavelet, f during model calculates ₀=15Hz, β are constant; (x ₀, y ₀, z ₀) be the locus of focus, x, y and z are respectively the position on x-axis, y-axis and z-axis direction.

7. a kind of three dimensional anisotropic elastic-wave numerical modeling method according to claim 5, it is characterized in that, described step 3 specifically comprises the following steps:

Step 3.1: the differential difference approximation in three dimensional anisotropic equations for elastic waves is substituted, obtains the propagation equation of corresponding finite difference scheme, in described propagation equation, space sampling step length and time-sampling step-length meet the stability condition of this numeric format;

Step 3.2: adopt the mode of Region Decomposition to bring the force value on each net point into propagation equation and calculate, obtain the wave field value at per a moment.

8. a three dimensional anisotropic elastic-wave numerical modeling system, comprises MBM, focus module, propagation module and data exchange module;

Described MBM is used for setting up dielectric model, dielectric model is carried out to grid is discrete obtains multiple net point;

Described focus module, for calculating source function, calculates the force value on each net point according to source function;

Described propagation module is used for three dimensional anisotropic equations for elastic waves to be converted to propagation equation, brings the force value on each net point into propagation equation and calculates, obtain the wave field value at per a moment;

Described data exchange module is used for the zoning determining each net point according to wave field value, carries out subregion and carries out exchanges data to partition boundaries data, completing elastic-wave numerical modeling.

9. a kind of three dimensional anisotropic elastic-wave numerical modeling system according to claim 8, it is characterized in that, described data exchange module comprises area determination module and subregion Switching Module;

Described area determination module is used for the zoning determining each net point according to all wave field value, and adopts the mode of absorbing boundary to determine border, zoning, simulates wave propagation in underground medium in zoning;

Described subregion Switching Module is used for carrying out subregion to all zonings, carries out exchanges data to the data boundary of adjacent sectors, obtains simulation elasticity wave datum, completes elastic-wave numerical modeling.

10. a kind of three dimensional anisotropic elastic-wave numerical modeling system according to claim 8 or claim 9, it is characterized in that, in described subregion Switching Module, the employing GPU-Direct technology of data boundary carries out exchanges data.