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

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

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

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

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

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

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

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

0 引言

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

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

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

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

1 C++ AMP簡介

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

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

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

1.1 array與array_view

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

1.2 extent與index

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

1.3 核函數parallel_for_each

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

1.4 C++AMP與常規GPU代碼的差異分析

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

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

2 基本理論與并行方法實現

2.1 三維彈性波波動方程交錯網格差分格式

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

(1)

(2)

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

2.2 并行算法思路及實現

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

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

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

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

2.3 三維并行程序的實現

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

3 性能測試

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

表1 性能實驗中設備的參數

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

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

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

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

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

4 三維鹽丘模型實例

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

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

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

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

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

本例中的數據傳輸采用常規模式：CPU-緩沖區-GPU,其中存在雙重數據拷貝,因此加速效果一般。由此可見,主機和設備間的數據傳輸是提升效率的瓶頸問題。通過array_view包裝主機端臨時數組,可消除緩沖區和拷貝操作,實現CPU與GPU內存直接訪問，降低數據傳輸的影響[7]。

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

5 結論及建議

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

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

3)GPU并行效率的提升主要依賴于硬件設備及算法的優化,主要體現在硬件設備架構及流式處理器的數量方面及數據存儲與拷貝優化等方面。

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

[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平臺上的實現[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] 龍桂華,趙宇波,李小凡,等. 三維交錯網格有限差分地震波模擬的GPU集群實現[J]. 地球物理學進展,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++加速大規模并行計算[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] 牟永光,裴正林.三維復雜介質地震波場數值模擬[M].北京:石油工業出版社,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]. 東營：中國石油大學出版社，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].東營：中國石油大學出版社，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)；中國石油天然氣集團公司科學研究與技術開發項目(2016A-33)

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

1001-1749(2017)05-0643-06