etDevice函數(shù)需要在MPI_Init之前調(diào)用，并利用MPI環(huán)境變量實現(xiàn)對GPU 的指定，MVAPICH2 :MV2_C0MM_W0RLD_L0CAL_RANK。
[0081 ] d)MV2_USE_CUDA 需要在執(zhí)行 mpirun 時添加。
[0082] e)對于節(jié)點內(nèi)的數(shù)據(jù)交換采用GPU Direct P2P方式(GPU直連點對點），在完成上 述配置后可以直接利用cudaMemcoycopy進(jìn)行GPU至ijGPU的數(shù)據(jù)傳輸，其具體描述參見圖5。 [0083] f)對于跨節(jié)點的數(shù)據(jù)交換采用GPU Direct RDMA方式(RDMA技術(shù)全稱遠(yuǎn)程直接數(shù) 據(jù)存取，就是為了解決網(wǎng)絡(luò)傳輸中服務(wù)器端數(shù)據(jù)處理的延遲而產(chǎn)生的，利用Infiniband網(wǎng) 絡(luò)直接進(jìn)行GPU直連稱之為GPU Direct RDMA)，在完成上述配置后可以直接采用MPI_Isend 與MPI_IreCV發(fā)送和接受設(shè)備內(nèi)存，而非像以往一樣需要先將設(shè)備內(nèi)存做顯式的 cudaMemcpy到主機內(nèi)存再通過MPI_I send與MPI_Irecv發(fā)送和接收，再通過顯式的 cudaMemcpy到設(shè)備端，進(jìn)行通信，兩種通信區(qū)別可見圖6。
[0084] 實施例5中，在實施例1-4任一實施例的基礎(chǔ)上，所述步驟1具體包括以下步驟：
[0085] 步驟1.1:根據(jù)地質(zhì)背景條件、實際測得的巖石物理測試數(shù)據(jù)和測井資料數(shù)據(jù)建立 介質(zhì)模型；
[0086]步驟1.2:采用形狀規(guī)則的三維網(wǎng)格對介質(zhì)模型進(jìn)行網(wǎng)格離散，得到網(wǎng)格點。
[0087] 實施例6中，在實施例1-5任一實施例的基礎(chǔ)上，所述震源函數(shù)在空間上采用高斯 函數(shù)，在時間上采用Ricker子波，所述震源函數(shù)的具體公式為：
[0088] s(x，y，z，t)=g(x，y，z) · f(t)公式（1)
[0089] 其中，以七）=(1-2(對成）2)6叉〇(-(對呲）2) 公式（2)
[0090]
[0091 ] 式中：t表示時間，fo表示Ricker子波的中心頻率，模型計算中fo = 15Ηζ，β為常數(shù)； (義0,7〇,2())為震源的空間位置，1、7和2分別為1軸、7軸和2軸方向上的位置。
[0092]實施例6中，在實施例1-5任一實施例的基礎(chǔ)上，所述步驟3具體包括以下步驟： [0093]步驟3.1:將三維各向異性彈性波方程中的微分用差分近似替代，得到相應(yīng)的有限 差分格式的傳播方程，所述傳播方程中空間采樣步長和時間采樣步長滿足該數(shù)值格式的穩(wěn) 定性條件；
[0094] 步驟3.2:采用區(qū)域分解的方式將每個網(wǎng)格點上的壓力值帶入傳播方程進(jìn)行計算， 得到每一刻的波場值。
[0095] 此實施例中，進(jìn)一步，所述步驟3.2中采用區(qū)域分解的方式利用GPU在各個計算節(jié) 點對三維各向異性彈性波方程進(jìn)有限差分計算，得到每一刻的波場值的具體方法如下： [0096]首先，給出三維各向異性介質(zhì)下的彈性波波動方程的數(shù)學(xué)描述，三維各向異性彈 性波動方程其形式為：
[0097] =4%+5，公式⑷
[0098]式中，u為波場函數(shù);〇為柯西應(yīng)力張量;Fi是單位體積的體力項，其可以被等效為 應(yīng)力源，P是介質(zhì)物質(zhì)密度，巧表示二階時間偏導(dǎo)。對于應(yīng)力張量利用線性彈性理論可以得 到：
[0099] 〇ij = cijkieki,公式（5)
[0100] 式中ekl為線性應(yīng)變張量，cljkl為四階張量的彈性剛度矩陣，該參數(shù)需要根據(jù)介質(zhì) 的實際情況給出；而對于eid我們有
[0101]
[0102] 其中ekl代表線性應(yīng)變張量，代表對kth維求空間導(dǎo)數(shù)，ui表示1 th波場位移。根 據(jù)此式我們就可以得行彈性波動方程的數(shù)值模擬，利用GPU進(jìn)行的彈性波數(shù)值模擬的具體 計算步驟如下：
[0103] 3.2.1我們定義一個仏\心\隊的計算網(wǎng)格，并且指定隊個時間步長，這樣在空間 和時間的每一個網(wǎng)格點就可以利用一個四元組進(jìn)行表示：即，[x，y，z|t] = [pAx，qAy，rA z | η △ t]。因此在指定點的連續(xù)的波場位移記錄就可以利用離散的網(wǎng)格點進(jìn)行表示：
[0104] ιη|χ,υ,ζ|ριηρ'γ'+ 公式（7)
[0105] 我們近似一階偏導(dǎo)4通過一個如下的Μ階精度的的中心差商算子Dx[ ·]，可以得 到：
[0106]
[0107] 式中Wa是一個權(quán)重系數(shù)多項式，即有限差分系數(shù)，Μ為有限差分計算階數(shù)；
[0108] 空間差分算子Dy[ ·]和Dz[ ·]與Dx[ ·]相似分別表示y和ζ方向的偏導(dǎo)數(shù)。將所有 這三者差分算子代入方程(6)，用以求解本構(gòu)方程中的偏導(dǎo)數(shù)。
[0109] 3.2.2對于時間離散，我們使用一個標(biāo)準(zhǔn)的空間二階精度近似來計算方程(6)中的
[0110]
[0111] 將上述的差分算子代入方程(4)并且通過重新排列這些離散項，我們就可以通過 時間步增的方式進(jìn)行未知波場求解。例如對于< ΑΦ+1，我們需要給定其當(dāng)前和之前的波場 值彳^劣^1以及根據(jù)有限差分計算空間計算模版所需要 x-，y-，z-方向相鄰點。圖 3表述了空間有限差分計算模版的在當(dāng)前時間步長關(guān)于點<#"所需要的計算網(wǎng)格點。
[0112] 如圖7所示，為本發(fā)明實施例1所述的一種三維各向異性彈性波數(shù)值模擬系統(tǒng)，包 括建模模塊1、震源模塊2、傳播模塊3和數(shù)據(jù)交換模塊4;
[0113] 所述建模模塊1用于建立介質(zhì)模型，對介質(zhì)模型進(jìn)行網(wǎng)格離散得到多個網(wǎng)格點； [0114]所述震源模塊2用于計算震源函數(shù)，根據(jù)震源函數(shù)計算每個網(wǎng)格點上的壓力值； [0115]所述傳播模塊3用于將三維各向異性彈性波方程轉(zhuǎn)換為傳播方程，將每個網(wǎng)格點 上的壓力值帶入傳播方程進(jìn)行計算，得到每一刻的波場值；
[0116] 所述數(shù)據(jù)交換模塊4用于根據(jù)波場值確定每個網(wǎng)格點的計算區(qū)域，進(jìn)行分區(qū)并對 分區(qū)邊界數(shù)據(jù)進(jìn)行數(shù)據(jù)交換，完成彈性波數(shù)值模擬。
[0117] 實施例2中，在實施例1的基礎(chǔ)上，所述數(shù)據(jù)交換模塊包括區(qū)域確定模塊和分區(qū)交 換模塊；
[0118] 所述區(qū)域確定模塊用于根據(jù)所有波場值確定每個網(wǎng)格點的計算區(qū)域，并采用吸收 邊界的方式確定計算區(qū)域邊界，在計算區(qū)域內(nèi)模擬地下介質(zhì)中波的傳播；
[0119] 所述分區(qū)交換模塊用于對所有計算區(qū)域進(jìn)行分區(qū)，對相鄰分區(qū)的邊界數(shù)據(jù)進(jìn)行數(shù) 據(jù)交換，得到模擬彈性波數(shù)據(jù)，完成彈性波數(shù)值模擬。
[0120] 實施例3中，在實施例2的基礎(chǔ)上，所述分區(qū)交換模塊中邊界數(shù)據(jù)的采用GPU-Direct 技術(shù)進(jìn)行數(shù)據(jù)交換。
[0121] 本發(fā)明的基礎(chǔ)是地震波的傳播理論及其高性能GPU數(shù)值模擬方法，本發(fā)明基于GPU Direct通信技術(shù)實現(xiàn)了高性能的三維各向異性介質(zhì)彈性波數(shù)值模擬方法，在具體示例中實 施步驟分別為：
[0122] 1.介質(zhì)模型建立：
[0123] 以不同TI對稱性(各向同性、VTI、HTI)的三維彈性介質(zhì)作為測試模型，其彈性參數(shù) 如下：相同的各向同性參數(shù)，VpzS.Okm/sA =νρ /λ/?，以及p = 2000kg/m3,但使用不同的 丁11〇1118611各向異性參數(shù)。令￥1'1介質(zhì)，[￡1，6 1，丫1] = [0.2，-0.1，0.2]，!11'1介質(zhì)，[￡2,62，丫2]= [0.2，-0.1，0.2]。這些參數(shù)都通過相應(yīng)的變換公式(Thomsen，1986)轉(zhuǎn)換進(jìn)剛度矩陣中。圖 8a-c分別表示用于脈沖響應(yīng)的IS0、VTI、HTI介質(zhì)6X6的剛度矩陣Cij，用不同的顏色塊來表 示不同的值。
[0124] 2.模型離散化：
[0125] 對設(shè)計的介質(zhì)模型進(jìn)行網(wǎng)格離散，測試地震數(shù)據(jù)維度大小為NxXNyXNz = 2003,網(wǎng) 格間距為 Δ χ= Δ y= Δ z = 0.005km。
[0126] 3.