CN105467443B - A kind of three dimensional anisotropic elastic-wave numerical modeling method and system - Google Patents
A kind of three dimensional anisotropic elastic-wave numerical modeling method and system Download PDFInfo
- Publication number
- CN105467443B CN105467443B CN201510902580.8A CN201510902580A CN105467443B CN 105467443 B CN105467443 B CN 105467443B CN 201510902580 A CN201510902580 A CN 201510902580A CN 105467443 B CN105467443 B CN 105467443B
- Authority
- CN
- China
- Prior art keywords
- wave
- elastic
- gpu
- data
- mesh point
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The present invention relates to a kind of three dimensional anisotropic elastic-wave numerical modeling method and system, wherein method includes：Step 1：Dielectric model is set up, carries out that grid is discrete obtains multiple mesh points to dielectric model；Step 2：Source function is calculated, the pressure value on each mesh point is calculated according to source function；Step 3：Three dimensional anisotropic equations for elastic waves is converted into propagation equation, bringing the pressure value on each mesh point into propagation equation is calculated, and obtains the wave field value at per a moment；Step 4：The zoning of each mesh point is determined according to wave field value, subregion is carried out and data exchange is carried out to partition boundaries data, complete elastic-wave numerical modeling.The present invention realizes the three dimensional elasticity ripple numerical simulation accelerated using GPU under complex dielectrics, propose the implementation for accelerating data transfer using GPU Direct technologies, avoid substantial amounts of CPU to GPU, GPU to CPU data copy, the problem of realizing optimization communication performance bottleneck.
Description
Technical field
The present invention relates to a kind of three dimensional anisotropic elastic-wave numerical modeling method and system, belong to geophysical exploration neck
Domain.
Background technology
Seismic wave numerical modeling based on communication theory of the seismic wave in underground medium, exploration seismology with it is natural
It is widely used in seismology.Conventional three-dimensional acoustic wave equation, elastic wave isotropism numerical simulation are wide at present
General is applied in each geophysics field such as numerical simulation, imaging and inverting.But it is big for how effectively to realize
(such as horizontal cross isotropism HTI is either for yardstick acquisition mode (such as wide-azimuth collection) and complex anisotropic medium
Orthotropy) three dimensional elasticity wave field numerical simulation study still have very big challenge, in actual applications
To widely using.In addition the conventional extensive three-dimensional acoustic wave based on CPU, elastic-wave numerical modeling, it usually needs great Liang Zhuan
Use PC cluster resource.From hardware cost or for calculating in energy consumption, its cost is all sufficiently expensive.Therefore quickly carry
Height calculates performance, significant reduction and calculates cost for realizing that the practical application of three dimensional elasticity ripple anisotropy values simulation has
Significance.
In last decade, using graphics processing unit (GPU) carry out compute-intensive applications speed up to have been obtained for
The development for formula of advancing by leaps and bounds.Graphics processing unit (GPU) has the memory bandwidth of high speed due to it, is at least higher by compared to CPU
The calculating processing core of two orders of magnitude, is more suitable for single-instruction multiple-data (SIMD) computation schema of parallel computation, and lower
Energy consumption cost, be just widely applied to the association area of computational science.For exploration geophysics field, for using figure
The interest of processing unit (GPU) is also being significantly increased, and GPU is used to accelerate in seismic processing by increasing research
Core algorithm, such as earthquake numerical simulation, seismic imaging, the high-precision inverting of earthquake.
The content of the invention
Find prior art and its there is problems with by research：
Time finite element method method is applied to and realizes research on complex dielectrics wave propagation algorithm very in GPU equipment
It is few, particularly with how handling the extensive three-dimensional problem research just less (Heinrich that there are a large amount of multinode data exchanges
Et al, 2014).
Presently, there are some GPU seimic wave propagations simulation implementation, for three-dimensional GPU cluster implementation we
It is described by taking the scheme of (imperial osmanthus China etc.) as an example.The program utilizes the ground that Region Decomposition technology will can not be calculated on single GPU
Plastid carries out coarseness decomposition along Z-direction, data boundary is exchanged using message passing interface (MPI), so that with MPI+
CUDA mode realizes the numerical simulation of large scale 3-D seismics wave field.But this method, which has very big one, is,
The speed-up ratio calculated using GPU cluster for single GPU and CPU compared to being remarkably decreased, and the reason for causing this result is GPU
Data copy between GPU to CPU and CPU to GPU equipment of the time-consuming consumption of the major part of calculating between GPU cluster node
On.
The technical problems to be solved by the invention are, in view of the shortcomings of the prior art there is provided one kind by GPU Direct skills
Art is combined the obtained three dimensional anisotropic elastic wave for optimizing communication based on GPU Direct with finite difference numerical simulation method
Method for numerical simulation and system, strong guarantee is provided for numerical simulation.
The technical scheme that the present invention solves above-mentioned technical problem is as follows：A kind of three dimensional anisotropic elastic-wave numerical modeling side
Method, specifically includes following steps：
Step 1：Dielectric model is set up, carries out that grid is discrete obtains multiple mesh points to dielectric model；
Step 2：Source function is calculated, the pressure value on each mesh point is calculated according to source function；
Step 3：Three dimensional anisotropic equations for elastic waves is converted into propagation equation, by the pressure value band on each mesh point
Enter propagation equation to be calculated, obtain the wave field value at per a moment；
Step 4：The zoning of each mesh point is determined according to wave field value, subregion is carried out and partition boundaries data is carried out
Data exchange, completes elastic-wave numerical modeling.
The beneficial effects of the invention are as follows：From ess-strain equation, realize and added using graphics processing unit (GPU)
Three dimensional elasticity ripple numerical simulation under fast complex dielectrics, and for three-dimensional problem introduce that GPU equipment faced due to region point
Between multinode caused by solution, intra-node communication bottleneck problem conducts in-depth research and analyzes, it is proposed that utilize GPU
Direct technologies accelerate the implementation of data transfer, it is to avoid substantial amounts of CPU to GPU, GPU to CPU data copy, realization
The problem of optimization communication performance bottleneck.By GPU figures acceleration equipment and GPU Direct communication optimization strategies are introduced, can be notable
Lifting overall computational performance, the elastic numerical value of three dimensional anisotropic that can be realized with lower hardware cost and less time
Simulation, is the various such as reverse-time migrations of the algorithm dependent on Wave equation forward modeling, and the application of full waveform inversion is provided effectively
Ensure.
On the basis of above-mentioned technical proposal, the present invention can also do following improvement.
Further, the step 4 specifically includes following steps：
Step 4.1：The zoning of each mesh point is determined according to all wave field values, and it is true by the way of absorbing boundary
Devise a stratagem calculates zone boundary, and the propagation of underground medium medium wave is simulated in zoning；
Step 4.2：Subregion is carried out to all zonings, data exchange is carried out to the data boundary of adjacent sectors, obtained
The elastic wave number evidence of simulation, completes elastic-wave numerical modeling.
Further, the use GPU-Direct technologies of data boundary carry out data exchange in the step 4.2.
Beneficial effect using above-mentioned further scheme is to accelerate three dimensional elasticity ripple numerical simulation using GPU, significantly improve
Calculating speed, is realized compared to conventional CPU cluster, can reach 20-30 times of speed lifting.
Further, the data exchange in the step 4.2 is specifically included：
Data variation most slow axial direction in zone boundary is calculated along along dielectric model and carries out Region Decomposition, multiple subregions are obtained,
All subregions are assigned in multiple GPU, each subregion independently performs calculating on a GPU；The side of each two adjacent sectors
Boundary's data are swapped again.
Beneficial effect using above-mentioned further scheme is counted using between GPU Direct technologies acceleration node, in node
According to the implementation of transmission, our scheme is because which obviating substantial amounts of CPU to GPU, GPU to CPU data copy, so that
The problem of realizing optimization communication performance bottleneck.
Further, the step 1 specifically includes following steps：
Step 1.1：Set up according to geology background condition, actually measured physical test of rock data and Well logging Data
Dielectric model；
Step 1.2：It is discrete to dielectric model progress grid using the three-dimensional grid of regular shape, obtain mesh point.
Beneficial effect using above-mentioned further scheme is that dielectric model is separated into mesh point, and mesh point is used as wave field
The carrier of value, constitutes wave field, the primary condition propagated as wave field by multiple mesh points of three-dimensional structure.
Further, the source function spatially uses Gaussian function, described in time using Ricker wavelets
The specific formula of source function is：
S (x, y, z, t)=g (x, y, z) f (t) formula (1)
Wherein, f (t)=(1-2 (π f0t)2)exp(-(πf0t)2) formula (2)
Formula (3)
In formula：T represents time, f0Represent the centre frequency of Ricker wavelets, f in model calculating0=15Hz, β are constant；
(x0, y0, z0) be focus locus, x, y and z are respectively the position on x-axis, y-axis and z-axis direction.
Further, the step 3 specifically includes following steps：
Step 3.1：Differential in three dimensional anisotropic equations for elastic waves is substituted with difference approximation, obtains corresponding limited
The propagation equation of difference scheme, described hollow sampling step length of propagation equation and time sampling step-length meet the steady of the numeric format
Qualitative condition；
Step 3.2：The pressure value on each mesh point is brought into propagation equation by the way of Region Decomposition to be calculated,
Obtain the wave field value at per a moment.
The technical scheme that the present invention solves above-mentioned technical problem is as follows：A kind of three dimensional anisotropic elastic-wave numerical modeling system
System, including modeling module, focus module, propagation module and data exchange module；
The modeling module is used to set up dielectric model, carries out that grid is discrete obtains multiple mesh points to dielectric model；
The focus module is used to calculate source function, and the pressure value on each mesh point is calculated according to source function；
The propagation module is used to three dimensional anisotropic equations for elastic waves being converted to propagation equation, by each mesh point
Pressure value bring propagation equation into and calculated, obtain the wave field value at per a moment；
The data exchange module is used for the zoning that each mesh point is determined according to wave field value, carries out subregion and to dividing
Area's data boundary carries out data exchange, completes elastic-wave numerical modeling.
The beneficial effects of the invention are as follows：From ess-strain equation, realize and added using graphics processing unit (GPU)
Three dimensional elasticity ripple numerical simulation under fast complex dielectrics, and for three-dimensional problem introduce that GPU equipment faced due to region point
Between multinode caused by solution, intra-node communication bottleneck problem conducts in-depth research and analyzes, it is proposed that utilize GPU
Direct technologies accelerate the implementation of data transfer, it is to avoid substantial amounts of CPU to GPU, GPU to CPU data copy, realization
The problem of optimization communication performance bottleneck.By GPU figures acceleration equipment and GPU Direct communication optimization strategies are introduced, can be notable
Lifting overall computational performance, the elastic numerical value of three dimensional anisotropic that can be realized with lower hardware cost and less time
Simulation, is the various such as reverse-time migrations of the algorithm dependent on Wave equation forward modeling, and the application of full waveform inversion is provided effectively
Ensure.
On the basis of above-mentioned technical proposal, the present invention can also do following improvement.
Further, the data exchange module includes area determination module and subregion Switching Module；
The area determination module is used for the zoning that each mesh point is determined according to all wave field values, and using absorption
The mode on border determines zoning border, and the propagation of underground medium medium wave is simulated in zoning；
The subregion Switching Module is used to carry out subregion to all zonings, and line number is entered to the data boundary of adjacent sectors
According to exchange, obtain simulating elastic wave number evidence, complete elastic-wave numerical modeling.
Further, the use GPU-Direct technologies of data boundary carry out data exchange in the subregion Switching Module.
Beneficial effect using above-mentioned further scheme is to accelerate three dimensional elasticity ripple numerical simulation using GPU, significantly improve
Calculating speed, is realized compared to conventional CPU cluster, can be reached 20~30 times of speed lifting, be utilized GPU Direct technologies
The implementation of data transfer between acceleration node, in node, our scheme is arrived because which obviating substantial amounts of CPU to GPU, GPU
CPU data copy, it is achieved thereby that the problem of optimization communication performance bottleneck.
Brief description of the drawings
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 needs 25 nets for the first derivative in the embodiment of the present invention using eight rank limited precision difference approximation central points
Lattice point figure；
Fig. 3 is implementing procedure of the three dimensional anisotropic elastic-wave numerical modeling in single GPU equipment in the embodiment of the present invention
Figure；
Fig. 4 is data volume Region Decomposition and signal intelligence general introduction figure in the embodiment of the present invention；
Fig. 5 a are that MPI data transfer mode schematic diagrames are used in prior art interior joint；
Fig. 5 b be the embodiment of the present invention in use GPU Direct P2P data transfer mode schematic diagrames；
Fig. 6 a be prior art interior joint between use MPI data transfer schematic diagrames；
Fig. 6 b be the embodiment of the present invention in schematic diagram is carried out data transmission using GPU Direct RDAM modes；
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 are three kinds of elastic fluid rigidity systems with different TI symmetry in the embodiment of the present invention for testing
Matrix number schematic diagram；
Fig. 9 a-9c are that the three-dimensional impulse response of three kinds of elastic fluids with different TI symmetry in the embodiment of the present invention is shown
It is intended to；
Figure 10 is the test point of routine MPI communication modes and GPU Direct communication modes in interior joint of the embodiment of the present invention
Analyse comparative result figure；
Figure 11 is using routine MPI communication modes and logical using GPU Direct RDMA between interior joint of the embodiment of the present invention
The test result comparison diagram of letter mode.
In accompanying drawing, the list of parts representated by each label is as follows：
1st, modeling module, 2, focus module, 3, propagation module, 4, data exchange module.
Embodiment
The principle and feature of the present invention are described below in conjunction with accompanying drawing, the given examples are served only to explain the present invention, and
It is non-to be used to limit the 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
Comprise the following steps：
Step 1：Dielectric model is set up, carries out that grid is discrete obtains multiple mesh points to dielectric model；
Step 2：Source function is calculated, the pressure value on each mesh point is calculated according to source function；
Step 3：Three dimensional anisotropic equations for elastic waves is converted into propagation equation, by the pressure value band on each mesh point
Enter propagation equation to be calculated, obtain the wave field value at per a moment；
Step 4：The zoning of each mesh point is determined according to wave field value, subregion is carried out and partition boundaries data is carried out
Data exchange, completes elastic-wave numerical modeling.
In embodiment 2, on the basis of embodiment 1, the step 4 specifically includes following steps：
Step 4.1：The zoning of each mesh point is determined according to all wave field values, and it is true by the way of absorbing boundary
Devise a stratagem calculates zone boundary, and the propagation of underground medium medium wave is simulated in zoning；
Step 4.2：Subregion is carried out to all zonings, data exchange is carried out to the data boundary of adjacent sectors, obtained
The elastic wave number evidence of simulation, completes elastic-wave numerical modeling.
In embodiment 3, on the basis of embodiment 1 or 2, the use GPU-Direct of data boundary in the step 4.2
Technology carries out data exchange.
In embodiment 4, on the basis of embodiment 1-3 any embodiments, the data exchange in the step 4.2 is specifically wrapped
Include：
Data variation most slow axial direction in zone boundary is calculated along along dielectric model and carries out Region Decomposition, multiple subregions are obtained,
All subregions are assigned in multiple GPU, each subregion independently performs calculating on a GPU；The side of each two adjacent sectors
Boundary's data are swapped again.
In this embodiment, further, GPU- is used to the calculating data boundary of each calculate node at each moment
Direct technologies carry out data exchange and are implemented as follows：
When processing solves extensive threedimensional model, limited global memory causes single GPU equipment can not store entirely
Model meshes.Therefore, our algorithm must be by CUDA code extensions to operating in many GPU equipment, the isomery system of multinode
In structure, to make full use of GPU computing capabilitys to carry out high-performance calculation.Based on this reason, we are using a Region Decomposition
Scheme, our Region Decomposition scheme carries out Region Decomposition as shown in figure 3, changing most slow Y direction along calculating grid, from
And whole zoning is decomposed into N/M parts, N is the lattice number in Y-direction, and M is the GPU numbers of current computing cluster.Often
Individual subregion is independent to perform calculation procedure 1-3 in a GPU equipment.
Because finite difference method (as shown in Fig. 2 by taking 8 rank precision as an example) needs forwardly and rearwardly four on y directions
Point, thus in every sub-regions carry out finite difference operation GPU equipment need obtain other equipment on its zoning
The value of consecutive points, therefore the data boundary in adjacent subarea domain must swap on each time step, its data transfer path
As shown in figure 3, being divided into node the two types between node, as shown in figure 4, data volume Region Decomposition and signal intelligence are general
State, term illustrates in figure, PCI-E is PCI-Express abbreviation, refers to computer bus and interface standard.InfiniBand,
It is a kind of multi-concurrency network framework connected dedicated for server-side network, GPU Direct P2P, GPU Direct section
Implementation (being referred to as the direct-connected P-2-P technologies of GPU) in point, GPU Direct RDMA, RDMA technologies full name is long-range directly number
According to access, GPU Direct RDMA refer to implementation (turning into the long-range direct-connecting technologies of GPU) between the nodes of GPU Direct technologies.
Using GPU Direct P2P communications in node, GPU Direct RDAM communications are utilized between node.Generally, this different GPU
The communication of equipment room needs to carry out GPU to CPU and CPU to GPU data copy due to data, and whole communication process extremely consumes
When, this is the main bottleneck that conventional GPU cluster is calculated.
In this regard, the present invention optimizes these data communication using GPU Direct, its specific implementation is as follows：
A) CUDA versions are ensured on 6.0, and video card computing capability is on 2.5, to ensure to GPU Direct
Support.
B) being calculated for GPU for support GPU Direct characteristics using CUDA-Aware MPI, CUDA-Aware MPI
The MPI of optimization is realized, supports CUDA-Aware MPI MPI to realize have at this stage, more than MAVPCH21.9 versions or
The versions of OpenMPI more than 1.7, we are by taking MAVPICH2 as an example.
C) cudaSetDevice functions need to call before MPI_Init, and are realized using MPI environmental variances to GPU
Specify, MVAPICH2：MV2_COMM_WORLD_LOCAL_RANK.
D) MV2_USE_CUDA needs to add when performing mpirun.
E) GPU Direct P2P modes (GPU is direct-connected point-to-point) are used for the data exchange in node, in completion
State to match somebody with somebody to postpone and GPU to GPU data transfer can be directly carried out using cudaMemcoycopy, it is specifically described referring to Fig. 5.
F) for cross-node data exchange using GPU Direct RDMA modes (RDMA technologies full name remotely directly count
According to access, produced to solve the delay of servers' data processing in network transmission, utilize Infiniband nets
Network directly carries out that GPU is direct-connected to be referred to as GPU Direct RDMA), it can directly use MPI_Isend completing above-mentioned with postponing
Sent with MPI_Irecv and accepting device internal memory, rather than needed as in the past that device memory first to be done into explicit
CudaMemcpy is sent and received by MPI_Isend with MPI_Irecv again to host memory, then by explicit
CudaMemcpy is communicated to equipment end, two kinds of visible Fig. 6 of communication difference.
In embodiment 5, on the basis of embodiment 1-4 any embodiments, the step 1 specifically includes following steps：
Step 1.1：Set up according to geology background condition, actually measured physical test of rock data and Well logging Data
Dielectric model；
Step 1.2：It is discrete to dielectric model progress grid using the three-dimensional grid of regular shape, obtain mesh point.
In embodiment 6, on the basis of embodiment 1-5 any embodiments, the source function spatially uses Gauss
Function, in time using Ricker wavelets, the specific formula of the source function is：
S (x, y, z, t)=g (x, y, z) f (t) formula (1)
Wherein, f (t)=(1-2 (π f0t)2)exp(-(πf0t)2) formula (2)
Formula (3)
In formula：T represents time, f0Represent the centre frequency of Ricker wavelets, f in model calculating0=15Hz, β are constant；
(x0, y0, z0) be focus locus, 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 embodiments, the step 3 specifically includes following steps：
Step 3.1：Differential in three dimensional anisotropic equations for elastic waves is substituted with difference approximation, obtains corresponding limited
The propagation equation of difference scheme, described hollow sampling step length of propagation equation and time sampling step-length meet the steady of the numeric format
Qualitative condition；
Step 3.2：The pressure value on each mesh point is brought into propagation equation by the way of Region Decomposition to be calculated,
Obtain the wave field value at per a moment.
In this embodiment, further, saved in the step 3.2 by the way of Region Decomposition using GPU in each calculating
Point enters finite difference formulations to three dimensional anisotropic equations for elastic waves, and the specific method for obtaining the wave field value at per a moment is as follows：
First, the mathematical description of the elastic wave wave equation under 3-D elastic anisotropic media, three dimensional anisotropic bullet are provided
Its form of property wave equation is：
Formula (4)
In formula, u is wave field function；σ is cauchy stress tensor；FiIt is the physical item of unit volume, it can be equivalent to
Stress riser, ρ is medium material density,Represent second-order time local derviation.It can be obtained using linear elasticity theory for stress tensor
Arrive：
σij=cijklεkl, formula (5)
ε in formulaklFor linear strain tensor, cijklFor the elastic stiffness matrix of the tetradic, the parameter is needed according to medium
Actual conditions provide；And for εklWe have
Formula (6)
Wherein εklLinear strain tensor is represented,Represent and space derivation, u are asked to kth dimensionslRepresent l th wave field displacements.Root
Accordingly formula we can just obtain the numerical simulation of row Time Migration of Elastic Wave Equation, utilize the specific of the GPU elastic-wave numerical modelings carried out
Calculation procedure is as follows：
3.2.1 we define a Nx×Ny×NzCalculating grid, and specify NtIndividual time step, so in space
It can be just indicated with each mesh point of time using a four-tuple：I.e., [x, y, z | t]=[p Δs x, q Δ y, r
Δz|nΔt].Therefore it can be just indicated in the continuous wave field displacement record of specified point using discrete mesh point：
ui|X, y, z | t≈ui P, r, q | nFormula (7)
Our approximate single order local derviationsPass through the center difference operator D of a following M rank precisionx[], can obtain
Arrive：
Formula (8)
W in formulaαIt is a weight coefficient multinomial, i.e. finite difference coefficient, M is finite difference formulations exponent number；
Space difference operator Dy[] and Dz[] and Dx[] similar partial derivative for representing y and z directions respectively.Will be all
This three's difference operator substitutes into equation (6), to solve the partial derivative in constitutive equation.
3.2.2 for time discrete, we are approximately come in accounting equation (6) using the space second order accuracy of a standard
Formula (9)
Above-mentioned difference operator is substituted into equation (4) and by rearranging these discrete items, we can just pass through
The mode that time step increases carries out unknown wave field solution.For example forWe need to give its wave field current and before
Value AndX-, y-, z- directions consecutive points according to required for finite difference formulations space calculates masterplate.Figure
3 describe space finite difference formulations masterplate in current time step on pointRequired calculating mesh point.
As shown in fig. 7, be a kind of three dimensional anisotropic elastic-wave numerical modeling system described in the embodiment of the present invention 1, bag
Include modeling module 1, focus module 2, propagation module 3 and data exchange module 4；
The modeling module 1 is used to set up dielectric model, carries out that grid is discrete obtains multiple mesh points to dielectric model；
The focus module 2 is used to calculate source function, and the pressure value on each mesh point is calculated according to source function；
The propagation module 3 is used to three dimensional anisotropic equations for elastic waves being converted to propagation equation, by each mesh point
On pressure value bring propagation equation into and calculated, obtain the wave field value at per a moment；
The data exchange module 4 is used for the zoning that each mesh point is determined according to wave field value, carries out subregion and right
Partition boundaries data carry out data exchange, complete elastic-wave numerical modeling.
In embodiment 2, on the basis of embodiment 1, the data exchange module includes area determination module and subregion is handed over
Change the mold block；
The area determination module is used for the zoning that each mesh point is determined according to all wave field values, and using absorption
The mode on border determines zoning border, and the propagation of underground medium medium wave is simulated in zoning；
The subregion Switching Module is used to carry out subregion to all zonings, and line number is entered to the data boundary of adjacent sectors
According to exchange, obtain simulating elastic wave number evidence, complete elastic-wave numerical modeling.
In embodiment 3, on the basis of embodiment 2, the use GPU- of data boundary in the subregion Switching Module
Direct technologies carry out data exchange.
The basis of the present invention is the communication theory and its high-performance GPU method for numerical simulation of seismic wave, and the present invention is based on GPU
The Direct communication technologys realize high performance 3-D elastic anisotropic media elastic-wave numerical modeling method, real in specific example
Applying step is respectively：
1. dielectric model is set up：
Test model, its elastic parameter are used as using the three dimensional elasticity medium of different TI symmetry (isotropism, VTI, HTI)
It is as follows：Identical isotropism parameter, vp=2.0km/s,And ρ=2000kg/m3, but using different
Thomsen anisotropic parameterses.Make VTI media, [ε1, δ1, γ1]=[0.2, -0.1,0.2], HTI media, [ε2, δ2, γ2]
=[0.2, -0.1,0.2].These parameters are all changed into stiffness matrix by corresponding transformation for mula (Thomsen, 1986).
Fig. 8 a-c represent the stiffness matrix C of ISO, VTI, HTI medium 6 × 6 for impulse response respectivelyij, with different color blocks come
Represent different values.
2. model discretization：
Discrete to the dielectric model progress grid of design, test geological data dimension size is Nx×Ny×Nz=2003,
Grid spacing is Δ x=Δ y=Δs z=0.005km.
3. source function：
Source function spatially uses Gaussian function, and Ricker wavelets are used on the time, and form is as follows：
S (x, y, z, t)=g (x, y, z) f (t) formula (1)
Wherein, f (t)=(1-2 (π f0t)2)exp(-(πf0t)2) formula (2)
Formula (3)
In formula：T represents time, f0Represent the centre frequency of Ricker wavelets, f in model calculating0=15Hz, β are constant；
(x0, y0, z0) be focus locus, 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：
Carry out finite difference numerical simulation using 8 rank precision, then it is correspondingWe utilize public affairs
The wave field that formula (4)-(6) carry out each step is calculated.
5. Region Decomposition and GPU Direct communications：
Using foregoing Region Decomposition scheme, carry out in the Y direction in Region Decomposition, node using GPU Direct P2P
Communicated, communicated between node using GPU Direct RDMA.It using the machine is three two-way GPU clothes that test environment, which is,
Business node is used as test platform.Each node is configured with the core Xeon E5-2680V2 2.80GHz processors of two-way 10,128GB
DDR3 internal memories.GPU equipment is NVIDIA Tesla K40c (Kepler), and Infiniband adapters are Mellanox
ConnectX-3 IB QDR MT26428, all nodes operate in RHEL 6.3, and Infiniband driving version is OFED-
2.3-1.0.1.All MPIs of the test based on MVAPICH2-2.1a release versions.
6. numerical simulation result is analyzed：
Fig. 9 is the wave field snapshot that three test models are obtained, and can be clearly seen that the shadow that anisotropy is caused to wave field
Ring.Figure 10 is that the ordinate in the communication test in the GPU Direct P2P nodes of progress, figure is speed-up ratio, is by GPU
The calculating time with obtained from single CPU calculating time, can be clearly seen that, employing GPU Direct P2P node
Interior data transfer has significant computational efficiency to be lifted compared to traditional MPI methods.Figure 11 is the GPU Direct carried out
Communication test between RDMA nodes, what the MPI-GDR wherein shown in legend was represented is the GPU Direct RDMA meaning, is carried out
The test of two platforms, 1 piece of K40C of each nodes of platform A totally three nodes, each nodes of platform B 3 pieces of K40C totally three ranks
Section.Test result shows that communications of the GPU Direct RDMA to cross-node is obviously improved effect, and with amount of communications
Increase its effect it is obvious all the more.
The foregoing is only presently preferred embodiments of the present invention, be not intended to limit the invention, it is all the present invention spirit and
Within principle, any modification, equivalent substitution and improvements made etc. should be included in the scope of the protection.
Claims (8)
1. a kind of three dimensional anisotropic elastic-wave numerical modeling method, it is characterised in that specifically include following steps：
Step 1：Dielectric model is set up, carries out that grid is discrete obtains multiple mesh points to dielectric model；
Step 2：Source function is calculated, the pressure value on each mesh point is calculated according to source function；
Step 3：Three dimensional anisotropic equations for elastic waves is converted into propagation equation, the pressure value on each mesh point is brought into biography
Broadcast equation to be calculated, obtain the wave field value at per a moment；
Step 4：The zoning of each mesh point is determined according to wave field value, subregion is carried out and data is carried out to partition boundaries data
Exchange, complete elastic-wave numerical modeling；
The step 4 specifically includes following steps：
Step 4.1：The zoning of each mesh point is determined according to all wave field values, and determines by the way of absorbing boundary meter
Zone boundary is calculated, the propagation of underground medium medium wave is simulated in zoning；
Step 4.2：Subregion is carried out to all zonings, data exchange is carried out to the data boundary of adjacent sectors, simulated
Elastic wave data, complete elastic-wave numerical modeling.
2. a kind of three dimensional anisotropic elastic-wave numerical modeling method according to claim 1, it is characterised in that the step
The use GPU-Direct technologies of data boundary carry out data exchange in rapid 4.2.
3. a kind of three dimensional anisotropic elastic-wave numerical modeling method according to claim 2, it is characterised in that the step
Data exchange in rapid 4.2 is specifically included：
Data variation most slow axial direction in zone boundary is calculated along along dielectric model and carries out Region Decomposition, multiple subregions are obtained, by institute
There is subregion to be assigned in multiple GPU, each subregion independently performs calculating on a GPU；The number of boundary of each two adjacent sectors
According to swapping again.
4. a kind of three dimensional anisotropic elastic-wave numerical modeling method according to claim any one of 1-3, its feature exists
In the step 1 specifically includes following steps：
Step 1.1：Medium is set up according to geology background condition, actually measured physical test of rock data and Well logging Data
Model；
Step 1.2：It is discrete to dielectric model progress grid using the three-dimensional grid of regular shape, obtain mesh point.
5. a kind of three dimensional anisotropic elastic-wave numerical modeling method according to claim 4, it is characterised in that the shake
Source function spatially uses Gaussian function, and in time using Ricker wavelets, the specific formula of the source function is：
S (x, y, z, t)=g (x, y, z) f (t) formula (1)
Wherein, f (t)=(1-2 (π f0t)2)exp(-(πf0t)2) formula (2)
In formula：T represents time, f0Represent the centre frequency of Ricker wavelets, f in model calculating0=15Hz, β are constant；(x0,
y0,z0) be focus locus, x, y and z are respectively the position on x-axis, y-axis and z-axis direction.
6. a kind of three dimensional anisotropic elastic-wave numerical modeling method according to claim 4, it is characterised in that the step
Rapid 3 specifically include following steps：
Step 3.1：Differential in three dimensional anisotropic equations for elastic waves is substituted with difference approximation, corresponding finite difference is obtained
The propagation equation of form, described hollow sampling step length of propagation equation and time sampling step-length meet the stability of the numeric format
Condition；
Step 3.2：The pressure value on each mesh point is brought into propagation equation by the way of Region Decomposition to be calculated, and is obtained
Wave field value per a moment.
7. a kind of three dimensional anisotropic elastic-wave numerical modeling system, including modeling module, focus module, propagation module and data
Switching Module；
The modeling module is used to set up dielectric model, carries out that grid is discrete obtains multiple mesh points to dielectric model；
The focus module is used to calculate source function, and the pressure value on each mesh point is calculated according to source function；
The propagation module is used to three dimensional anisotropic equations for elastic waves being converted to propagation equation, by the pressure on each mesh point
Force value is brought propagation equation into and calculated, and obtains the wave field value at per a moment；
The data exchange module is used for the zoning that each mesh point is determined according to wave field value, carries out subregion and to subregion side
Boundary's data carry out data exchange, complete elastic-wave numerical modeling；
The data exchange module includes area determination module and subregion Switching Module；
The area determination module is used for the zoning that each mesh point is determined according to all wave field values, and uses absorbing boundary
Mode determine zoning border, in zoning simulate underground medium medium wave propagation；
The subregion Switching Module is used to carry out subregion to all zonings, and data friendship is carried out to the data boundary of adjacent sectors
Change, obtain simulating elastic wave number evidence, complete elastic-wave numerical modeling.
8. a kind of three dimensional anisotropic elastic-wave numerical modeling system according to claim 7, it is characterised in that described point
The use GPU-Direct technologies of data boundary carry out data exchange in area's Switching Module.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510902580.8A CN105467443B (en) | 2015-12-09 | 2015-12-09 | A kind of three dimensional anisotropic elastic-wave numerical modeling method and system |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510902580.8A CN105467443B (en) | 2015-12-09 | 2015-12-09 | A kind of three dimensional anisotropic elastic-wave numerical modeling method and system |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105467443A CN105467443A (en) | 2016-04-06 |
| CN105467443B true CN105467443B (en) | 2017-09-19 |
Family
ID=55605340
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510902580.8A Expired - Fee Related CN105467443B (en) | 2015-12-09 | 2015-12-09 | A kind of three dimensional anisotropic elastic-wave numerical modeling method and system |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105467443B (en) |
Families Citing this family (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106556868B (en) * | 2016-11-01 | 2019-05-07 | 中国石油天然气股份有限公司 | The quantitative identification method and device of groove |
| CN108072895B (en) * | 2016-11-09 | 2020-09-15 | 中国石油化工股份有限公司 | Anisotropic prestack reverse time migration optimization method based on GPU |
| CN106776455B (en) * | 2016-12-13 | 2020-08-21 | 苏州浪潮智能科技有限公司 | Single-machine multi-GPU communication method and device |
| CN107703538B (en) * | 2017-09-14 | 2019-08-09 | 上海交通大学 | Underground unfavorable geology survey data acquisition analysis system and method |
| CN108563802B (en) * | 2017-12-29 | 2021-12-17 | 中国海洋大学 | Method for improving numerical simulation precision of seismic converted waves |
| CN108802819B (en) * | 2018-06-26 | 2019-11-08 | 西安交通大学 | A kind of trapezoidal grid finite difference Simulation of Seismic Wave method of uniform depth sampling |
| CN112528456B (en) * | 2019-09-18 | 2024-05-07 | 曙光信息产业(北京)有限公司 | Heterogeneous node computing system and method |
| CN110866964A (en) * | 2019-11-08 | 2020-03-06 | 四川大学 | GPU accelerated ellipsoid clipping map terrain rendering method |
| US11281825B2 (en) * | 2020-06-30 | 2022-03-22 | China Petroleum & Chemical Corporation | Computer-implemented method for high speed multi-source loading and retrieval of wavefields employing finite difference models |
| CN112764105B (en) * | 2020-10-16 | 2022-07-12 | 中国石油大学(华东) | HTI medium quasi-longitudinal wave forward simulation method and device, storage medium and processor |
| CN112904417B (en) * | 2021-01-21 | 2022-05-03 | 中国石油大学(华东) | Finite difference simulation method for seismic wave propagation of prepressing solid medium |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102269820B (en) * | 2010-06-01 | 2016-01-13 | 潜能恒信能源技术股份有限公司 | A kind of 3-D seismics pre-Stack Reverse formation method |
| EP2761328A4 (en) * | 2011-09-28 | 2015-05-06 | Conocophillips Co | Reciprocal method two way wave equation targeted data selection for seismic acquisition of complex geologic structures |
| WO2016008100A1 (en) * | 2014-07-15 | 2016-01-21 | 杨顺伟 | Three-dimensional seismic anisotropic medium reverse time migration imaging method and device |
-
2015
- 2015-12-09 CN CN201510902580.8A patent/CN105467443B/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN105467443A (en) | 2016-04-06 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105467443B (en) | A kind of three dimensional anisotropic elastic-wave numerical modeling method and system | |
| Abdelkhalek et al. | Fast seismic modeling and reverse time migration on a GPU cluster | |
| Komatitsch et al. | High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster | |
| Komatitsch et al. | Modeling the propagation of elastic waves using spectral elements on a cluster of 192 GPUs | |
| Fabien-Ouellet et al. | Time-domain seismic modeling in viscoelastic media for full waveform inversion on heterogeneous computing platforms with OpenCL | |
| CN103969627B (en) | The extensive D integral pin-fin tube analogy method of GPR based on FDTD | |
| Shin et al. | 3D Laplace-domain full waveform inversion using a single GPU card | |
| Geevers et al. | New higher-order mass-lumped tetrahedral elements for wave propagation modelling | |
| CN105005072B (en) | The PML border dimensionally seismic wave propagating mode utilizing CUDA intends method | |
| Xue et al. | An efficient GPU implementation for locating micro-seismic sources using 3D elastic wave time-reversal imaging | |
| Schneider et al. | Micro-applications for communication data access patterns and MPI datatypes | |
| Karavaev et al. | A technology of 3D elastic wave propagation simulation using hybrid supercomputers | |
| Moradi et al. | Quantum computing in geophysics: Algorithms, computational costs, and future applications | |
| Abdelkhalek et al. | Fast seismic modeling and reverse time migration on a graphics processing unit cluster | |
| CN107273333A (en) | Three-dimensional mt inverting parallel method based on GPU+CPU heterogeneous platforms | |
| Noble et al. | High-performance 3D first-arrival traveltime tomography | |
| Xiang et al. | GPU acceleration of CFD algorithm: HSMAC and SIMPLE | |
| CN106662665B (en) | The interpolation and convolution of rearrangement for the processing of faster staggered-mesh | |
| Meister et al. | 2D adaptivity for 3D problems: Parallel SPE10 reservoir simulation on dynamically adaptive prism grids | |
| Komatitsch et al. | A simulation of seismic wave propagation at high resolution in the inner core of the Earth on 2166 processors of MareNostrum | |
| CN105974471B (en) | A kind of quick forward modelling method of the more GPU of seismic data based on asynchronous flow | |
| Esler et al. | GAMPACK (GPU accelerated algebraic multigrid package) | |
| Wang et al. | FP-AMR: A Reconfigurable Fabric Framework for Adaptive Mesh Refinement Applications | |
| Roca et al. | GPU-accelerated sparse matrix-vector product for a hybridizable discontinuous Galerkin method | |
| ZHAO et al. | Fast calculation of converted wave traveltime in 3‐D complex media |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20170919 Termination date: 20191209 |