基于MIC平臺的三維瑞雷面波有限差分并行模擬

曾瑞庚， 熊章強,2， 張大洲,2

(1.中南大學 地球科學與信息物理學院，長沙 410083;2.中南大學 有色金屬成礦預測教育部重點實驗室, 長沙 410083)

基于MIC平臺的三維瑞雷面波有限差分并行模擬

曾瑞庚1， 熊章強1,2， 張大洲1,2

(1.中南大學 地球科學與信息物理學院，長沙 410083;2.中南大學 有色金屬成礦預測教育部重點實驗室, 長沙 410083)

三維地震波動方程有限差分數值模擬，受到計算量大且內存消耗大地制約而計算效率較低，嚴重影響其發展和應用。目前，一種全新的異構眾核協處理器，為三維并行計算提供了一條解決計算效率和數據存儲問題的重要途徑。基于三維有限差分計算的特點，介紹了面向MIC(Many Integrated Core)平臺的三維瑞雷波有限差分并行算法移植和性能優化，采用ofload模式將計算核函數加載到MIC上，在MIC協處理上使用openMP多線程并行方式，并通過循環合并、nocopy數據傳輸、SIMD向量化和CPU+MIC協同計算方式進行優化。通過模型計算可知：三維數值計算程序在MIC上具有近線性加速比， MIC+CPU協同并行計算的性能可以達到純CPU節點的4.5倍，大大提高了計算效率。通過對模擬結果與解析解的比較，以及對頻散曲線特征的分析，驗證了數值模擬的正確性。

MIC平臺； 瑞雷面波； 有限差分； offload模式

0 引言

瑞雷面波具有衰減小、信噪比高、分辨率高以及在層狀介質中的頻散特性等優點，在近地表工程勘察以及石油勘探靜校正等領域有著廣泛地應用。目前，有關瑞雷面波勘探的研究在空間上從二維轉向三維，在反演方面的研究也從傳統的頻散曲線反演逐漸轉向更為精確的波形反演。快速、高效地獲取不同地質模型下的波場信息,是當前瑞雷面波研究領域亟待解決的問題。

有關瑞雷面波的有限差分數值模擬，國內、外學者已經進行了大量的研究工作[1-2]。根據當前的計算機硬件及軟件的發展趨勢，采用并行計算技術增大數據存儲量，減小正演模擬時間是一條非常有效的途徑。當前，一種基于MPI和OpenMP的并行算法得到了眾多研究者的青睞。方金偉等[3]研究了一種基于OpenMP的一階聲波方程波場并行算法，該方法在微機商用多核CPU上獲得了2倍多的性能提升；何兵壽等[4]引入消息傳遞接口技術(MPI)，實現了二維彈性波動方有限差分并行求解，大大提高了彈性波動方程的計算效率；張明財等[5]基于MPI并行算法對三維瑞雷面波進行正演模擬，闡述了三維模型區域劃分方法和進程間的通信方式，達到了擴大模型規模、加快模擬速度的目的；宋鵬等[6]針對三維聲波方程模擬的大計算量和大內存消耗問題，研究了基于MPI+OpenMP的三維聲波方程數值模擬并行算法，大大節省了內存消耗。

由于當前計算機以及計算技術的迅速發展，高性能計算開始進入異構、眾核時代。隨著CUDA的日趨成熟，以GPU為加速器的GPGPU已逐漸成為研究的熱點。Komatish D等[7]研究了三維有限差分波場模擬在NVIDIA GPU上的實現，取得了理想的加速效果；龍桂華等[8]采用了MPI+CUDA的方式在GPU集群上實現了大尺度三維地震波數值模擬；Intel于2012年推出的基于MIC(Many Integrated Core)眾合架構的處理器Intel Xeon Phi協處理器,因其成本低、性能高等特點得到了重大關注。與GPGPU不同的是，MIC是基于傳統CPU架構和X86指令集開發的協處理器，在兼容性、編程模式高度一致性和編程簡易性等方面具有優勢[9]。

筆者從三維交錯網格有限差分出發，基于CPU+MIC平臺，使用MIC的offload編程模式，對有限差分算法進行了移植和優化。

1 三維瑞雷面波波動方程

在三維各向同性介質中，一階速度-應力彈性波波動方程可表示為[10]：

(1)

其中:i、j、k∈{x,y,z}；ρ為介質的密度；fi表示體力i的分量；vi表示速度i的分量；σij表示應力分量；δij為克羅內克函數；λ、μ為拉梅系數。

筆者所采用的是時間2階空間4階交錯網格有限差分算法，其各分量定義方式和具體公式參考文獻[5]。自由邊界條采用聲學-彈性介質邊界代替自由邊界法[11]，吸收邊界條件采用非分列式完全匹配層(NPML)邊界條件[12]。

2 MIC平臺移植和優化

MIC是基于與CPU一致的x86指令架構，即可以看成一個協處理器，也可以當做一個獨立的節點，具有非常靈活的編程模式(圖1)。與常規CPU相比MIC擁有50個以上核心，其中每個核心支持4個硬件線程，其性能優于超線程。MIC協處理器與CPU結合緊密，可以在不改動源代碼或者改動很少的情況下使用該協處理器。本文采用CPU為主，MCI為輔的編程模式，即MIC offload模式[13]。這種模式應用最為普遍，符合傳統的協處理器應用方式。MIC支持現有的編程標準工具和語言(C、C++、FORTRAN等)和常用的并行技術(如Pthread、MPI、OpenMP等)，筆者使用了MPI+ OpenMP混合并行方法以實現MIC和CPU間的協同計算。

圖1 MIC架構圖Fig.1 Architecture of MIC

2.1 MIC平臺實現

offload模式適合串行程序中含有計算量大且并行度高的函數或者代碼段。有限差分正演模擬中速度和應力的計算是程序中最耗時的部分，這部分的差分計算具有良好的并行度和空間局部性，適合在MIC上并行處理。如圖2所示，程序主函數在CPU上運行，運行到需要并行計算部分時，通過在速度和應力計算函數上添加offload編譯指導語句，將該函數加載到MIC上并行處理。offload引語本身不會指示代碼并行執行，需要借助于OpenMP、MPI、Pthreads等并行編程方式實現。

圖2 offload模式流程圖Fig.2 Flow chart of offload mode

三維有限差分算法總體上由一個四重循環組成，最外層為時間循環迭代，其次分別為空間上的x、y、z循環，這部分循環計算采用OpenMP多線程并行方式可以很好地并行化。OpenMP是通過引語的方式對程序的循環語句并行化，其要求是循環語句不存在相互依賴關系。由于在時間循環上當前的結果依賴于上一次的循環結果，所以不適合OpenMP的并行。對于空間循環任一x、y、z方向上的循環都是獨立的，不存在相關性，因此可以在任一方向上對循環指定并行。使用offload指導語句結合openMP，即可實現在MIC上的并行計算。

2.2 性能優化

針對MIC架構的特點，筆者從循環合并、nocopy數據傳輸、SIMD向量化、CPU+MIC協同計算四個角度對程序進行優化。

1) 循環合并。為了提高OpenMP在MIC上多線程并行計算性能，將空間方向上的3層循環合并為2層循環，提高每個線程上循環計算的長度以減少線程并行的開銷。圖3為對150×150×100的網格模型在開啟不同進程數量情況下，循環合并和未進行循環合并對運算速度的影響。當線程數量小于循環計算長度(150)的時候，循環合并后比未合并需要更多的計算時間。然而線程數超過循環計算的長度時(>150)，循環合并后的計算時間比循環未合并時少，最高獲得1.7倍的計算速度，循環合并具有更好的計算性能。

圖3 循環合并對計算時間的影響Fig.3 Effect of loop confution with different threads number

2) nocopy數據傳輸優化。CPU和MIC之間在硬件上是通過PCI-E進行數據傳輸，這直接導致了CPU和MIC之間的傳輸速度較慢，因此要盡量減少CPU和MIC之間的數據傳輸。在默認情況下offload的in、out等語句每次開始時在MIC上申請空間，結束后釋放空間。但計算過程中大部分變量不需要每次都在MIC上重新申請和釋放空間，這些空間和數據在時間循環過程中可以重復利用。可以利用nocopy語句來減少CPU和MIC的數據傳遞。如圖4所示，首先在時間循環開始之前為MIC上所需的變量分配空間，并把變量的值傳入到MIC中。在時間迭代過程中，CPU和MIC之間只傳輸必要的數據，而中間變量則使用nocopy來表示以減少數據傳遞，時間循環迭代完畢后釋放MIC上的空間。這種方法可大大減少MIC和CPU之間的數據通信，提高并行效率。

圖4 數據傳輸優化Fig.4 Optimization of data transfer

3) SIMD向量化。向量化是使用向量處理單元進行批量計算的方法，MIC的處理器支持512bit寬的指令，利用這些指令級別的并行可以很好地提高程序性能。有限差分算法的循環體是大量而重復的循環計算，且涉及到連續的存儲空間訪問。通過在循環語句前添加!DECMYM IVDEP語句使循環自動向量化，同時編譯器將忽略假定的數據依賴關系強制循環向量化。

4) CPU+MIC協同計算。MIC的offload模式僅把程序中高度并行計算任務提交給MIC，而CPU只負責控制管理，這使得CPU的計算資源沒有得到充分利用。如果把計算任務同時交給CPU和MIC，就可以充分利用兩者計算資源，提升計算效率。由于MIC卡和CPU之間計算能力不同，為了達到CPU和MIC的之間的負載均衡，需要對MIC和CPU的計算任務進行合理的分配。

筆者采用MPI多進程方式進行任務分配，主進程把模型數據廣播給所有進程，其中一個進程負責調用CPU并啟用OpenMP多線程進行計算，另一個進程負責調用MIC并把數據傳輸到MIC中進行計算。CPU和MIC之間的負載均衡采用靜態任務劃分方法[14]，即讓CPU和MIC分別對同樣網格大小的模型做一次計算得到程序在CPU和MIC的計算時間比值，根據這個時間比值來對模型任務進行劃分。MIC+CPU協同運算示意圖如圖5所示。

圖5 CPU+MIC協同運算示意圖Fig.5 Cooperative computation with CPU and MIC

項目配置CPU2個Intel(R)Xeon(R)E5-2609CPU,主頻2.50GHz,4核,10ML3cacheMIC1個IntelXeonPhi7110p,主頻1.238GHz,61核,244線程,內存CPU:32GB MIC:16GB軟件CPU:CentOS6.3x64編譯器:Intel○RFortranComposerXEMIC驅動:MPSS3.4.3

3 模型測試與計算

3.1 模型測試

為了測試程序在MIC中的并行效率，設計一個大小為150×150×100的模型，網格間距dx=dy=dz=1，總網格數為2.25×106。采樣間隔為0.1 ms，采樣時長為0.3 s。使用了一臺MIC服務器進行計算，計算環境配置如表1所示。

加速比(speedup)，是同一個任務在單處理器系統和并行處理器系統中運行消耗的時間的比率，用來衡量并行系統或程序并行化的性能和效果，其計算公式為：

Sp=T1/TP

(2)

其中：Sp是加速比；T1是單處理器下的運行時間；TP是在有P個處理器并行系統中的運行時間。

圖6為程序在單個MIC上運行的加速比和計算時間，以測試程序在MIC中的可擴展性和穩定性。從圖6中可以看出，計算加速比隨著線程數的增加幾乎呈線性增長，每個核的線程均開啟的情況下，最高加速比達到34.4，多核計算加速比幾乎達到線性。

圖6 程序在MIC卡上的運行時間和加速比Fig.6 The runtime and speed-up ratio on MIC

圖7 異構并行的計算效率對比Fig.7 Performance comparision of heterogeneous parallel computation

圖7為程序分別在CPU上串行計算(單核)、CPU上openMP多線程并行計算(8個線程)、MIC上多線程并行計算(240個線程)、MIC+CPU協同計算的計算時間和加速比。對比圖6和圖7可以看出，MIC的單核計算能力比CPU單核弱很多，在MIC線程數為240的情況下的計算能力是單核Xeon CPU的4.7倍，是CPU節點在開啟8個OpenMP線程情況下的3.7倍。采用CPU+MIC協同運算相對于純CPU節點加速比達到了4.5，異構并行獲得了更好的加速效果。

3.2 均勻半空間瑞雷面波正演計算

為了驗證本文有限差分并行計算算法的正確性，設計一均勻各向同性介質模型。模型大小為300 m×300 m×200 m ，網格間距dx=dy=dz=1。縱波速度vp=1 000 m/s，橫波速度vs=577 m/s，密度ρ=2 100 kg/m3。震源為25 Hz的高斯一階導函數，位于模型的(150，150，0)處。采樣間隔為0.1 ms，采樣時長為300 ms，延時為40 ms。

圖8為檢波器距離震源100 m處水平速度分量vx和垂直速度分量vz的數值解和解析解的對比，其中解析解由Lamb問題解析法得到，其中格林函數的計算采用Cagniard-de Hoop 解析法[15-16]。從圖8可以看出，vx和vz的數值解與對應的解析解吻合較好，其均方差為分別為4.1%和4.0%。

圖9為在y=150 m的平面上速度分量vx、vy、vz分別在110 ms、160 ms、240 ms時刻的波場快照。 從圖9中可以很清楚地看到瑞雷面波、縱波和橫波。沿深度方向，瑞雷面波的能量集中于大約25 m內，深度超過這個距離后，P波開始變得清晰，這說明瑞雷面波的穿透深度約為一個波長。此外，在吸收邊界處沒有產生強烈的反射，說明本文中使用的邊界條件是有效的。

圖10為各分量波場記錄，其中瑞雷面波的同向軸皆為直線，這說明瑞雷面波在均勻各向同性介質中沒有發生頻散。圖11為用垂直速度分量vz提取的頻散曲線，從圖中可知，均勻介質中的頻散曲線為一條直線，其相速度約為530 m/s，且未發生頻散現象，這符合面波相關理論。以上實例計算，說明了數值模擬的正確性。

圖8 vx和vz分量數值解與解析解對比(距離震源100 m)Fig.8 Compersion of numerical solution and analytical solution with vx and vz component(a)水平速度分量vx; (b)垂直速度分量vz

圖9 各方向速度分量在x=150 m平面上不同時刻的波場快照Fig.9 The snapshot at plane of x=150 m with different velocity component(a)vx分量在110 ms時刻;(b) vx分量在160 ms時刻;(c)vx分量在240 ms時刻;(d)vy分量在110 ms時刻;(e) vy分量在160 ms時刻;(f)vy分量在240 ms時刻;(g)vz分量在110 ms時刻;(h) vz分量在160 ms時刻;(i)vz分量在240 ms時刻

3.3 水平層狀介質模型瑞雷面波正演計算

設計一個3層水平層狀速度遞增模型，模型參數如表2所示。震源采用25 Hz的高斯一階導數，模型大小為300 m×300 m×200 m，網格間距取dx=dy=dz=1。

表2 水平層狀模型參數

圖12為在y=150 m的平面上速度分量vx、vy、vz分別在120 ms、180 ms、260 ms時刻的波場快照。由于層狀介質中各層的參數不同，從圖12中可以很明顯看出頻散，且大部分能量聚集在地表附近。由于反射界面的存在，由震源激發的波與在各界面發生多次反射的波相互疊合形成了多模態的瑞雷面波。

圖13(a)為x=150 m平面上vz分量的單炮波場記錄。為了分析在瑞雷面波在層狀介質中的頻散特性，在vz分量波場記錄上利用相移法提取頻散曲線，如圖12(b)所示。圖中白色方塊為該模型參數下計算得到的解析解[17]，從圖中可見基階模式和高階模式的頻散曲線吻合得較好。

4 結論

計算效率的高低是影響三維數值計算在實際應用中的重要因素。筆者利用Intel最新推出的MIC平臺協處理器Xeon Phi具有較高的加速比和良好可擴展性的特點，將三維瑞雷面波數值模擬在MIC平臺的offload模式程序進行移植并優化，包括循環合并、nocopy數據傳輸優化、編譯自動向量化以及CPU+MIC協同計算。通過對所設計的均勻介質模型進行波場模擬測試可知：基于OpenMP的并行程序在MIC上具有近線性加速比，MIC+CPU協同并行計算的性能為單個CPU節點的4.5倍。因此，利用MIC+CPU協同并行計算方式，可極大提高計算效率，為后續的全波形反演、波動方程疊前深度偏移等提高計算效率。此外，針對大規模計算問題，利用MPI多進程方式可以將程序擴展到單節點多MIC卡或者MIC集群上計算。

致謝

感謝中南大學高性能計算中心提供計算平臺。

圖12 各方向速度分量在x=150 m平面上不同時刻的波場快照Fig.12 The snapshot at plane of x=150 m with different velocity component(a)vx分量在120 ms時刻;(b) vx分量在180 ms時刻;(c)vx分量在260 ms時刻;(d)vy分量在120 ms時刻;(e) vy分量在180 ms時刻;(f)vy分量在260 ms時刻;(g)vz分量在120 ms時刻;(h) vz分量在180 ms時刻;(i)vz分量在260 ms時刻

圖13 單炮vz波場記錄及其頻散曲線Fig.13 Wave record of vz conponant and its dispersion curves(a)x=150 m平面上的波場記錄;(b)頻散曲線

[1] 董良國, 馬在田, 曹景忠, 等. 一階彈性波動方程交錯網格高階差分法[J]. 地球物理學報, 2000,43(3):411-419.

DONG L G, MA Z T, CAO J Z, et al, Astaggered-grid high-order difference method of one-order el-astic wave equation[J]. Chinese Journal Geophysi-cs, 2000,43(3):411-419.(In Chinese)

[2] PETER M, JOZEF K, LADISLAV H. 3D fourth-order staggered-grid finite-difference schemes: stabil-ity and grdi dispersion[J]. Bulletin of the Seismol-ogical Society of Americal, 2000,90(3):587-603.

[3] 方金偉, 白洪濤, 孫惠敏. 基于OpenMP的一階聲波方程波場并行算法[J]. 物探化探計算技術, 2015, 37(5)： 622-627.

FANG J W, BAI H T, SUN H M, A parallel algori-thm of wave field using first-order acoustic wave equation based on OpenMP[J]. Computing Tech-niques for Geophysical and Geochemical Explorat-ion, 2015,37(5)：622-627.(In Chinese)

[4] 何兵壽,張會星,韓令賀.彈性波方程正演的粗粒度并行算法[J]. 地球物理學進展，2010,25(2)：650-656.

HE B S, ZHANG H X, HAN L H. Forward modeling of elastic wave equation with coarse-grained parallel algorithm[J]. Progress in Geophysics, 2010,25(2): 650-656.(In Chinese)

[5] 張明財,熊章強,張大洲. 基于MPI的三維瑞雷面波有限差分并行模擬[J]. 石油物探，2013,52(4)：354-362.

ZHANG M C, XIONG Z Q, ZHANG D Z. 3D finite difference parallel simulation of Rayleigh wave based on Message Passing Interface[J]. Geophysical Prospecting for petroleum, 2013,52(4):354-362.(In Chinese)

[6] 宋鵬, 解闖, 李金山, 等. 基于MPI+OpenMP的三維聲波方程正演模擬[J]. 中國海洋大學學報(自然科學版), 2015(09): 097-102.

SONG P, XIE C, LI J S, et al, The 3D modeling of acoustic wave equation based on MPI+OpenMP[J]. Periodical of Ocean University of China, 2015(09):097-102.(In Chinese)

[7] Komatish D, Tomatitsch D. Accelerating a three-dimensional finite-difference wave propagation code using GPU graphics cards. Geophys[J], 2010, 182(1): 389-402.

[8] 龍桂華,趙宇波, 李小凡, 等. 三維交錯網格有限差分地震波模擬的GPU集群實現[J]. 地球物理學進展, 2011(06): 1938-1949.

LONG G H, ZHAO Y B, LI X F, et al, Accelerati-ng 3D staggered-grid finite-difference seismic w-ave modeling on GPU cluster[J]. Progress in Ge-ophysicals, 2011(06): 1938-1949.(In Chinese)

[9] INTEL. Intel Xeon Ph Coprocessor DEVELOP-ER’S QUICK START GUIDE[OL].https://software.intel.com/sites/default/files/managed/ee/4e/intel-xeon-phi-coprocessor-quick-start-developers-guide.pdf,2013.

[10] HELSTHOLM S O. Finate-difference seismic m-odeling including surface topograrhy Texas U-niverzones [J].Physics of the Earth and Plantetary Interiors,2002,132(1):219-234.

[11] XU Y X, XIA J H, MILLER RD, Numerical invertigation of implemention of air-earth boundry by acoustic-elastic boundry approach[J]. Geophysics,2007,72(5):147-153.

[12] QIN ZHEN, LU MINGHUI, ZHENG XIAODONG,et al. The implemention of an improved NPML absorbing boundry condition in elastic wave modeling[J]. Applied Geophysics, 2009:113-121.

[13] 沈鉑, 張廣勇, 吳韶華, 等. 基于MIC平臺的 offload并行方法研究[J]. 計算機科學, 2014,41(6A): 477-480.

SHEN B, ZHANG G Y, WU S H, et al. Research of offload paralel method based on MIC platform[J]. Computer Science, 2014,41(6A):477-480.(In Chinese)

[14] 王恩東, 張清, 沈鉑, 等.MIC高性能計算編程指南[M].北京：中國水利水電出版社，2012.

WANG E D, ZHANG Q, SHENG B, et al. Many Integrated Core High Performance Programming Guid[M].Beijing: China Water&Power Press, 2012.(In Chinese)

[15] HOOP, A.T.DE. A modification of cagniard's method for solving seismic pulse problems[J]. Applied Scientific Research, 1960, 8(1):349-356.

[16] 巴特.地震學的數學問題[M]. 鄭治真，譯.北京: 科學出版社, 1976.

BART M.Mathematical problem of seismology[M].Zhen Z Z, translator.Beijing: Science Press, 1976.(In Chinese)

[17] HISADA,Y.An efficient method for computing green's functions for a layered half-space with sources and receivers at close depths[J].Bulletin of the Seismological Society of America,1994:1456-1472.

3DfinitedifferenceparallelsimulationofRayleighwavebasedonMICplatform

ZENG Ruigeng1, XIONG Zhangqiang1,2, ZHANG Dazhou1,2

(1.Central South University School of Geosciences and Info-Physics, Changsha 410083, China;2.Central South University Key Laboratory of Metallogenic Prediction of Nonferrous Metals, Ministry of Education, Changsha 410083, China)

The development and application of Numerical simulation of three-dimensional seismic wave equation is hindered for its huge computation and large computer memory. In recent years, MIC (Many Integrated Core), as a new kind of heterogeneous core coprocessor, provides a better choice for three-dimensional parallel modeling. Based on the characteristic of 3D finite difference computation, we introduce the code migration and performance optimization on MIC platform. We offload the kernel function to MIC with openMP parallel method, and optimize the performance by loop confusion, data transfer optimization, SIMD vectorization and CPU+MIC cooperative computation. Through the model computation, the speed-up ratio on MIC is nearly linear. The performance of the MIC+CPU cooperative computation was 4.5 times faster than pure CPU compute node. The method is verified by the comparison of numerical simulation result and analytical solution, as well as by the analysis of dispersion curves from wave record.

Many integrated core; Rayleigh wave; finite difference; offload model

P 631.4

A

10.3969/j.issn.1001-1749.2017.05.11

2016-09-14 改回日期： 2016-11-06

國家自然科學基金(41274123)；博士點基金(20130162110066)

曾瑞庚(1990- )，男，碩士，主要從事瑞雷波正演研究，E-mail：rgzeng01@163.com。

1001-1749(2017)05-0649-08