武泗海, 趙 虎, 尹 成， 賈 鵬

(1.西南石油大學(xué) a.地球科學(xué)與技術(shù)學(xué)院, b.天然氣地質(zhì)四川省重點(diǎn)實驗室,成都 610500;2.川慶鉆探地球物理勘探公司 成都 610213)

基于C++AMP并行加速的三維彈性波波動方程數(shù)值模擬

武泗海1a,b, 趙 虎1a,b, 尹 成1a,b， 賈 鵬2

(1.西南石油大學(xué) a.地球科學(xué)與技術(shù)學(xué)院, b.天然氣地質(zhì)四川省重點(diǎn)實驗室,成都 610500;2.川慶鉆探地球物理勘探公司 成都 610213)

在地震勘探方面,地震數(shù)值模擬技術(shù)在地球物理反演及觀測系統(tǒng)設(shè)計領(lǐng)域有著廣泛地應(yīng)用。這里主要提出并實現(xiàn)了利用C++AMP并行加速內(nèi)存/計算密集型的三維彈性波的正演模擬,使用C++AMP編寫的代碼可以有效地在圖形處理器(GPU)等支持?jǐn)?shù)據(jù)并行的硬件上運(yùn)行。在重點(diǎn)研究C++AMP三維彈性波數(shù)值模擬的并行算法實現(xiàn)的基礎(chǔ)上,結(jié)合CPU串行及CPU多核并行進(jìn)行加速對比，三維鹽丘實驗的結(jié)果說明，利用GPU的強(qiáng)大并行計算能力的C++AMP，可以極大地提高數(shù)值模擬的效率。此外，C++AMP在Windows平臺具有良好的可移植性，同時兼容NVIDA和AMD的設(shè)備，極大地提升了開發(fā)效率。

C++AMP; GPU并行計算; 三維彈性波; 正演

0 引言

波動方程數(shù)值模擬技術(shù)是研究地震波傳播特征的重要工具之一,同時也是地震反演的基礎(chǔ),在研究復(fù)雜地區(qū)地震資料采集、處理和解釋等領(lǐng)域應(yīng)用廣泛。Alterman[1]將有限差分法首次應(yīng)用于地震波模擬; Vireux[2]通過將速度和應(yīng)力分量置于空間格點(diǎn)和半格點(diǎn)的方式提高了空間微分的精度,并完成了對彈性介質(zhì)P-SV波和SH波的正演工作。實踐證明,對于三維復(fù)雜構(gòu)造,地質(zhì)體和巖性的地震模擬而言,交錯網(wǎng)格差分法以其綜合性能最好,是目前實用性最好的方法。

盡管基于交錯網(wǎng)格差分法的數(shù)值模擬方法具有效率高,實現(xiàn)方便等特點(diǎn),但對于高精度數(shù)值模擬,疊前波動方程逆時偏移,波形反演等嚴(yán)重依賴于正演效率的技術(shù),效率問題依然是制約上述技術(shù)的瓶頸。隨著GPU通用計算技術(shù)的飛速發(fā)展,上述局面得以改觀。張兵等[3]在CUDA平臺實現(xiàn)裂步法疊前深度偏移；龍桂華等[4]在GPU集群上實現(xiàn)了三維交錯網(wǎng)格有限差分?jǐn)?shù)值模擬。GPU并行技術(shù)的應(yīng)用極大地推進(jìn)了地震勘探新技術(shù)的應(yīng)用與推廣。

C++加速大規(guī)模并行計算技術(shù)(C++ Accelerated Massive Parallelism,C++AMP)是微軟公司發(fā)布的一套C++應(yīng)用程序GPU加速技術(shù),C++AMP是眾多的GPU并行編程技術(shù)之一,它是微軟基于 DirectX 11 技術(shù)實現(xiàn)的并行計算庫。其建立在C++語言規(guī)范之上,大幅降低實現(xiàn)異構(gòu)計算的難度。與常規(guī)的CPU串行和多核并行不同的是: C++AMP的主體在CPU端執(zhí)行,在執(zhí)行數(shù)據(jù)密集型的計算時,在GPU設(shè)備端創(chuàng)建眾多的輕量級線程,并交由GPU的計算單元加以執(zhí)行。在計算結(jié)束后,拷貝數(shù)據(jù)并返回CPU端。與CUDA和OpenCL等主流GPU并行技術(shù)相比, C++AMP是基于CPU/GPU共享內(nèi)存的設(shè)計,在實際工作中,普通計算機(jī)的GPU設(shè)備內(nèi)存受限,且不易于實現(xiàn)大規(guī)模的GPU集群并行。而這一共享內(nèi)存的設(shè)計,使得C++AMP不再受限于GPU內(nèi)存。另外,C++AMP的良好硬件兼容性，使得程序可以運(yùn)行于AMD和NVIDIA設(shè)備,這是前兩種技術(shù)所不具備的優(yōu)勢。

筆者重點(diǎn)研究C++AMP在加速三維彈性波數(shù)值模擬技術(shù)中并行算法的實現(xiàn)，及數(shù)據(jù)傳輸?shù)膬?yōu)化等問題。通過對比CPU串行及多核并行技術(shù)的加速結(jié)果,驗證了C++AMP優(yōu)異的加速性能。同時,在不同計算設(shè)備的對比結(jié)果也更進(jìn)一步說明了其良好的移植性。

1 C++ AMP簡介

GPU計算由來已久,目前主流的GPU并行技術(shù)包括OpenCL和CUDA。C++ AMP是微軟2012年推出的基于 DirectX 11 技術(shù)實現(xiàn)的并行計算庫,它是Visual Studio和C++編程語言的新擴(kuò)展,可真正意義上實現(xiàn)支持跨Windows平臺及GPU設(shè)備的編譯運(yùn)行。OpenCL作為第一個異構(gòu)計算標(biāo)準(zhǔn),更注重各平臺的移植性；CUDA是獲得最佳性能提升的首選,但其僅支持NVIDIA設(shè)備。C++AMP則更注重易編程性,具有實現(xiàn)簡單,主機(jī)端與設(shè)備端數(shù)據(jù)的隱式拷貝,GPU端線程間的自動負(fù)載均衡等特點(diǎn),可以快速地用于幫助開發(fā)人員充分適應(yīng)高度并行和異構(gòu)計算環(huán)境。因此,C++AMP無疑是在Windows平臺上進(jìn)行異構(gòu)計算的最佳選擇。

以向量加法的C++AMP與CUDA并行實現(xiàn)為例(圖2、圖3)，說明C++AMP的基本數(shù)據(jù)類型,核函數(shù)的執(zhí)行與線程的調(diào)度及其與常規(guī)GPU代碼的差異。

圖1 GPU并行編程接口示意圖Fig.1 The interfaces of GPU parallel programming

1.1 array與array_view

在C++AMP中,array與array_view類均可用于表示N維、T類型的多維數(shù)據(jù)。其區(qū)別在于：array用于申請設(shè)備端內(nèi)存,array中的數(shù)據(jù)位于加速器上(通常是GPU或虛擬加速器設(shè)備等)。array_view是數(shù)據(jù)的包裝器,其數(shù)據(jù)源既可位于主機(jī)端,也可以位于設(shè)備端。利用array_view包裝的主機(jī)端數(shù)據(jù),通過array_view在設(shè)備端訪問時,可自動實現(xiàn)隱式數(shù)據(jù)傳輸,這極大地方便了主機(jī)與設(shè)備的數(shù)據(jù)遷移。

1.2 extent與index

在GPU并行計算中,GPU線程規(guī)模和索引與數(shù)據(jù)一一對應(yīng)，才能保證結(jié)果的準(zhǔn)確性。 extent表示array與array_view數(shù)據(jù)的尺寸, index表示array與array_view中某一位置元素的索引。在圖2中示例程序中,線程的布局與數(shù)據(jù)的布局是一致的。另外,三維線程通常采用計算線程偏移量索引,C++AMP中則簡便地使用index<3>idx(x,y,z)對象進(jìn)行索引。

1.3 核函數(shù)parallel_for_each

parallel_for_each函數(shù)是C++AMP的核心,其函數(shù)參數(shù)為Lambda表達(dá)式,它定義了核函數(shù)的基本計算功能,還決定了線程的分配與布局。在串行程序中,通過循環(huán)進(jìn)行遍歷計算, C++AMP將數(shù)組分配給眾多輕量級線程,通過GPU大量的計算單元并行處理。在圖2中,parallel_for_each代替了for循環(huán),參數(shù)dev_c.exetent表示創(chuàng)建與dev_c維度一致的線程布局,表示變量按值傳入,(index<1> idx)表示生成1D的線程索引對象idx,resrict(amp)用于指定加速器類型。

1.4 C++AMP與常規(guī)GPU代碼的差異分析

為了對比C++AMP與常規(guī)GPU代碼的差異,筆者以向量加法的C++AMP與CUDA并行實現(xiàn)為例(分別如圖2，圖3代碼所示)說明C++AMP的基本數(shù)據(jù)類型,核函數(shù)的執(zhí)行與線程的調(diào)度及其與常規(guī)GPU(CUDA)代碼的差異。

通過對比兩種代碼可以看出,array_view的設(shè)計極大地簡化了數(shù)據(jù)傳輸?shù)拇a,因此C++AMP的代碼更加簡潔。在線程配置上,C++AMP利用Lambda表達(dá)式可以方便利用extent布局,idx索引的方式,特別是多維線程,索引更加方便。

2 基本理論與并行方法實現(xiàn)

2.1 三維彈性波波動方程交錯網(wǎng)格差分格式

對于非均勻各向同性介質(zhì)中的三維彈性波波動方程一階應(yīng)力-速度方程組,采用差分交錯網(wǎng)格,給出2L階空間差分精度,二階時間差分精度高階有限差分格式[6]：

(1)

(2)

式中：T表示應(yīng)力；V表示速度；ρ表示密度；λ、μ表示彈性參數(shù)。其余Txy、Txz、Tyy、Tyz、Tzz及vy、vz的差分格式由同理可得。

2.2 并行算法思路及實現(xiàn)

基于C++AMP實現(xiàn)三維彈性波數(shù)值模擬的時域GPU并行算法,通過對串行程序的分析可知,波場迭代部分為主要耗時計算,是并行加速的關(guān)鍵。并行實現(xiàn)的主要問題在于以下幾個方面。

1)顯存需求及分配。隨著模型網(wǎng)格的增加,三維彈性波數(shù)值模擬對內(nèi)存需求急劇增加。在常規(guī)GPU并行方法中,通過拷貝模型數(shù)據(jù)至GPU顯存中,實現(xiàn)計算單元快速訪問。然而，對于大尺寸模型而言，顯存受限使得這一做法不再適用。C++AMP通過array_view實現(xiàn)GPU與CPU內(nèi)存的共享,在一定程度上解決了顯存不足的問題。另外,對于訪問頻率較高的變量,可根據(jù)其大小分配在GPU端寄存器,共享存儲器等高速緩存中,以獲得性能的提升。

2)線程的布局與索引。C++AMP通過多線程并行執(zhí)行的方式展開循環(huán)方式加快計算速度,因此,線程配置與尋址非常重要,否則將導(dǎo)致錯誤的計算結(jié)果。線程布局是由模型數(shù)據(jù)的extent屬性確定,這基于模型與線程的維度及尺寸保持一致的情況。但實際情況并非如此,當(dāng)對數(shù)據(jù)的部分進(jìn)行操作時,如高階差分的邊界區(qū)域處理等。針對這一問題,如果在GPU采用if語句進(jìn)行分支運(yùn)算,分支結(jié)構(gòu)無法并行執(zhí)行而降低加速效果。如圖4代碼所示,C++AMP提供利用extent類實例化對象的方式,創(chuàng)建特定的線程布局,為上述操作提供了極大地方便。

3)GPU命令的排隊與異步執(zhí)行。對于三維波場模擬,數(shù)據(jù)的復(fù)制操作與波場的計算核函數(shù)的命令以排隊提交GPU,并排隊執(zhí)行。C++AMP提供立即隊列和自動隊列兩種方式,前者是將命令立即提交,而后者則批量提交。通常情況下,前者適用存在時延的需求操作。后者以時延為代價換取吞吐率,批量模式更能有效地使用硬件。任務(wù)在GPU上排序,其形式是以包含命令(核)的DMA緩沖區(qū)及該命令所使用的內(nèi)存資源的引用。在該緩沖執(zhí)行前,所需的全部資源均分配在GPU內(nèi)存。因此,對于較大模型的計算,建議采用立即隊列,以減小DMA緩沖區(qū)對顯存的需求。另外,C++AMP提供copy_async等異步執(zhí)行函數(shù)簡化了異步執(zhí)行的實現(xiàn)。通過異步執(zhí)行的方式在不同的GPU引擎上實現(xiàn)復(fù)制操作與核函數(shù)計算的重疊,可以進(jìn)一步提升性能。

2.3 三維并行程序的實現(xiàn)

C++AMP的并行算法主要實現(xiàn)為：在CPU端執(zhí)行程序主體,當(dāng)場值迭代計算時,通過parallel_for_each定義并提交GPU任務(wù),在GPU端執(zhí)行多線程計算,場值計算結(jié)束后保存數(shù)據(jù)。為了在代碼中實現(xiàn)在GPU端的并行計算,我們需要包含頭文件,并使用 concurrency名稱空間簡化代碼；再利用CPU 端的array_view加載場值變量及模型數(shù)據(jù)在GPU端實現(xiàn)并行執(zhí)行,其代碼如圖4所示。

3 性能測試

在小尺寸模型情況下,GPU顯存足夠存儲模型數(shù)據(jù)及場值迭代變量,利用array數(shù)組分配GPU內(nèi)存,可顯著降低訪問延遲。為了測試這種情況下的C++AMP的加速效果,采用CPU串行及OpenMP多核并行與C++AMP并行對比。另外,不同的GPU設(shè)備加速測試說明C++AMP具有良好的移植性,實驗硬件參數(shù)見表1。

表1 性能實驗中設(shè)備的參數(shù)

通過對不同尺寸模型的計算,統(tǒng)計運(yùn)行時結(jié)果如圖5所示,隨著模型尺寸的增加,串行計算耗時急劇增加；OpenMP并行利用多核處理器提升了近2倍效率。C++AMP利用GPU強(qiáng)大的并行處理能力，大大縮短了計算用時,與串行運(yùn)算相比，獲得近40倍的加速效果。

由于C++AMP支持不同的GPU設(shè)備,這使得其在提高計算效率的同時,也可以增加開發(fā)效率。因此,我們分別在三個不同的GPU設(shè)備進(jìn)行加速測試,測試結(jié)果如圖6所示。

圖5 串行,OpenMP及C++AMP計算效率對比示意圖Fig.5 Comparison of computational efficiency between CPU serial computation, OpenMP and C++AMP

圖6 不同GPU設(shè)備并行計算效率對比Fig.6 Comparison of computational efficiency between the difference GPU device

通過圖6中所反映的不同GPU設(shè)備的加速表現(xiàn),GPU可以顯著地提高計算效率,流處理單元越多,效率提升越明顯。尤其是GTX970,隨著模型尺寸的增加,依然表現(xiàn)出良好的加速效果。同時,HD6700顯卡在150*150*150模型尺寸時出現(xiàn)異常偏低效率，而且當(dāng)模型大于180*180*180后,其效率卻低于GT540,這也不符合其之前的高效表現(xiàn)。

4 三維鹽丘模型實例

在數(shù)值模擬實例中,為了測試C++AMP針對大尺寸模型的加速效果,利用三維鹽丘模型(圖7)進(jìn)行測試,其模型網(wǎng)格尺寸為676*201*676。 由于顯存受限,array_view可實現(xiàn)主機(jī)與設(shè)備內(nèi)存共享,但主機(jī)與設(shè)備間存在頻繁的數(shù)據(jù)傳輸。

圖7 三維鹽丘模型示意圖Fig.7 The 3-D Salt model

本次模擬的正演參數(shù)分別為：采樣間隔為0.5 ms；空間采樣間隔為8*9*8 m；震源坐標(biāo)為(2 704 m,2 704 m,18 m)；震源函數(shù)為30 Hz雷克子波；PML吸收層為20層；差分精度采用2階時間精度和16階空間精度。

圖8 三維鹽丘的波場示意圖Fig.8 The wave field snapshot of 3-D Salt model

圖8為C++AMP并行實現(xiàn)3D鹽丘模型的波場快照。為了驗證C++AMP的加速效果,由于三維正演串行計算耗時過長,我們將C++AMP與OpenMP并行計算進(jìn)行對比正演單炮的計算時間,采集時間均為3 s的計算時間統(tǒng)計如表2所示，與OpenMP并行相比，C++AMP獲得近3倍的性能提升。實驗設(shè)備：CPU為Intel(R) Xeon E5-2609 V3 @1.9GHz(12*核,32G),GPU為NVIDIA GTX 970(4G)。

本例中的數(shù)據(jù)傳輸采用常規(guī)模式：CPU-緩沖區(qū)-GPU,其中存在雙重數(shù)據(jù)拷貝,因此加速效果一般。由此可見,主機(jī)和設(shè)備間的數(shù)據(jù)傳輸是提升效率的瓶頸問題。通過array_view包裝主機(jī)端臨時數(shù)組,可消除緩沖區(qū)和拷貝操作,實現(xiàn)CPU與GPU內(nèi)存直接訪問，降低數(shù)據(jù)傳輸?shù)挠绊慬7]。

表2 OpenMP并行與C++AMP計算用時

5 結(jié)論及建議

1)利用C++AMP加速并行,在保證足夠精度的基礎(chǔ)上,可以極大地提升三維數(shù)值模擬的效率,在顯存足夠條件下，實驗中效率提升了近40倍。

2)針對大尺寸模型,GPU顯存受限,但基于共享內(nèi)存的設(shè)計依然可以實現(xiàn)加速并行,在3D鹽丘正演實例中,獲得了一定地性能提升。其性能受限于頻繁的數(shù)據(jù)傳輸。CPU臨時數(shù)組可消除數(shù)據(jù)拷貝，從而提升計算效率。

3)GPU并行效率的提升主要依賴于硬件設(shè)備及算法的優(yōu)化,主要體現(xiàn)在硬件設(shè)備架構(gòu)及流式處理器的數(shù)量方面及數(shù)據(jù)存儲與拷貝優(yōu)化等方面。

4)GPU集群技術(shù)適用于疊前偏移或全波形反演等大規(guī)模計算問題。而對于PC單機(jī), C++AMP具有良好的設(shè)備兼容性和并行實現(xiàn)的便捷性等優(yōu)勢。因此,C++AMP更適于加速小規(guī)模的正演,非線性反演的優(yōu)化和觀測系統(tǒng)評價計算等PC機(jī)能承受的常規(guī)問題的加速。

[1] ALTERMAN Z, KARAL F C. Propagation of elastic wave sin layer media by finite difference methods[J]. Bull. Seism. Soc. AM. , 1968,58(1):367-398.

[2] VIRIEUX, JEAN. P-SVwave propagation in heterogeneous media: Velocity‐stress finite‐difference method[J]. GEOPHYSICS , 1986,51(4):889-901.

[3] 張兵,趙改善，黃駿，等.地震疊前深度偏移在CUDA平臺上的實現(xiàn)[J].勘探地球物理,2008(6):427-432+403.

ZHANG B, ZHAO G S，HUANG J, et al. Seismic prestack depth migration on CUDA platform[J]. Progress in Geophysics,2008(6):427-432+403.(In Chinese)

[4] 龍桂華,趙宇波,李小凡,等. 三維交錯網(wǎng)格有限差分地震波模擬的GPU集群實現(xiàn)[J]. 地球物理學(xué)進(jìn)展,2011(6):1938-1949.

LONG G H, ZHAO Y B, LI X F, et al. Three dimensional staggered grid finite difference forward modeling on GPU cluster[J]. Progress in Geophysics,2011(6):1938-1949. (In Chinese)

[5] KATE GREGORY, ADE MILLER. C++AMP 用Visual C++加速大規(guī)模并行計算[M]. 車皓陽,譯.北京：人民郵電出版社,2014.

KATE GREGORY, ADE MILLER. C++ AMP Accelerated massive parallelism with Microsoft Visual C++[M].CHE H.Y., T rans,Beijing：The People's Posts and Telecommunications Press,2014. (In Chinese)

[6] 牟永光,裴正林.三維復(fù)雜介質(zhì)地震波場數(shù)值模擬[M].北京:石油工業(yè)出版社,2005.

MOU Y G,PEI Z L. Seismic numerical modeling for 3-D complex media [M]. Beijing：Petroleum Industry Press, 2005. (In Chinese)

[7] AMIT K AGARWAL. Shared memory support in C++ AMP - array_view[EB/OL].https://blogs.msdn.microsoft.com/nativeconcurrency. 2013-7-11.

[8] 李振春.地震疊前成像理論與方法[M]. 東營：中國石油大學(xué)出版社，2007.

LI Z C. Theory and methods of pre-stack migration and imaging[M]. Dongying: China University of Petroleum Press,2007.(In Chinese)

[9] 孫成禹.地震波理論與方法[M].東營：中國石油大學(xué)出版社，2007.

SUN C Y. Theory and methods of seismic waves[M]. Dongying: China University of Petroleum Press,2007.(In Chinese)

Acceleratedparallelisminnumericalsimulationfor3DelasticwaveequationwithC++AMP

WU Sihai1a,b, ZHAO Hu1a,b, YIN Cheng1a,b, JIA Peng2

(1.School of Geoscience and Technology, a.South-west Petroleum University，b.Sichuan Key Laboratory of Natural Gas Geology, Southwest Petroleum University, Chengdu 610500, China;2.Chuanqing Drilling & Exploration Corporation, Sichuan Geophysical Company, Chengdu 610213, China)

Seismic forward numerical simulations have extensive applications in the inversions and the geometry design. This paper introduces the parallel simulation of the memory/computing-intensive and large-scale 3D elastic wave equation with C++ AMP. C++AMP is Microsoft's technology for accelerating C++ applications by allowing code run on data-parallel hardware like graphic-processing units(GPUs). On the basis of research on parallel algorithm with C++AMP, we compared it with CPU serial algorithms and OpenMP. The results of numerical experiment demonstrate the effectiveness and portability of this method. Furthermore, C++ AMP has a good portability on windows platform because the same executable can run on NVIDIA and AMD devices. It can significantly improve develop efficiency.

C++AMP; GPU parallel computing; three-dimensional elastic wave; forward

P 631.4

A

10.3969/j.issn.1001-1749.2017.05.10

2016-10-08 改回日期： 2016-11-22

國家重大專項(2016ZX05024-001-003)；中國石油天然氣集團(tuán)公司科學(xué)研究與技術(shù)開發(fā)項目(2016A-33)

武泗海(1991-)，男，碩士，主要研究地震采集與高性能計算，E-mail：wsh_18980796471@sina.com。

1001-1749(2017)05-0643-06