基于OpenMP的一階聲波方程波場并行算法

方金偉， 白洪濤， 孫惠敏

(1. 吉林大學 a.地球探測科學與技術學院, b.公共計算機教學與研究中心, 長春 130012)

?

基于OpenMP的一階聲波方程波場并行算法

方金偉1a， 白洪濤1b*， 孫惠敏1a

(1. 吉林大學 a.地球探測科學與技術學院, b.公共計算機教學與研究中心, 長春 130012)

針對高階交錯網(wǎng)格技術的計算效率瓶頸，研究了一種基于OpenMP的一階聲波方程波場并行算法。通過對多組不同規(guī)模的模型測試，分析了并行效率在不同并行粒度下模型和CPU核心數(shù)目的關系。實驗結果表明，該方法達到了與串行方法相同的數(shù)據(jù)精度，在微機商用多核CPU上獲得了2倍多的性能提升。

一階聲波方程; 高階交錯網(wǎng)格; OpenMP; 并行效率

0 引言

地震波場數(shù)值模擬是一種研究復雜地區(qū)地震資料的有用手段。常見的波場模擬的方法：有限差分法、有限元法和偽譜法等，其中有限元法精度高，但計算量大，對計算機的內存要求高；有限差分法計算速度較快，占用內存小，編程相對簡單，而且通過加密網(wǎng)格可以很好地壓制頻散效應，是目前應用最廣泛的方法。憑借著一階速度-應力彈性波方程無需對彈性常數(shù)進行空間微分的優(yōu)點，所以此類波動方程在彈性波波場模擬中常被使用。交錯網(wǎng)格技術最早由Madariaga R[1]提出，Virieux在同性介質正演中用到的差分精度為o(Δt2+Δx2)的交錯網(wǎng)格，其相比于常規(guī)網(wǎng)格，在計算精度和效率均有所提高。Levande[4]把這種差分網(wǎng)格的精度進一步提高到了o(Δt2+Δx4)。Dablain[5]提出了高階差分方法，解決了聲波方程求解的精度和效率的矛盾。隨后，Crase[6]又在求解二階彈性波方程中用到了該方法。高階差分可以通過選取合適的時間和空間階數(shù)，采用空間的階數(shù)高于時間的階數(shù)的差分形式，既可以保證較高的計算精度，而且計算效率可以得到一定的提高，另外，可以選取盡可能大的空間步長，在保證精度的前提下，獲得效率的提升。其后很多學者又對其進一步發(fā)展，董良國等[7]給出了一階交錯網(wǎng)格高階差分解法，并詳細討論了其穩(wěn)定性條件；劉洋等[8]從Taylor級數(shù)展開式出發(fā)，推導出任意階導數(shù)的任意偶數(shù)階精度差分格式，并給出相應差分系數(shù)的公式。高階交錯網(wǎng)格有限差分理論有很高的差分精度，但其計算效率還有一定的提升空間，這里將結合現(xiàn)代微機CPU多核技術，來提高該方法的實際計算效率。

在現(xiàn)有的工藝條件下，摩爾定律限制了通過提升指令的吞吐量和時鐘速度來改善CPU性能的方法，故超線程、多核技術和計算機的緩存等技術孕育而生，尤其在微型商用機芯片中使用較多，多核技術更為矚目。AMD、Intel等主要微機CPU廠商通過從提高主頻轉向整合多個核心引擎，來改變處理器的性能。目前常見的多核多線程開發(fā)工具有pThreads庫、Win32 線程庫以及OpenMP等。pThreads庫是Linux下最常用的多線程支持庫，Windows操作系統(tǒng)也有其移植版pthreads-win32，它具有很好的可移植性，但使用難度比較大；Win32 線程庫擁有完善而復雜的函數(shù)庫，目前比較成熟，但對編程人員有較高的要求；OpenMP用于共享內存并行系統(tǒng)的多線程程序設計的一套指導性的編譯處理方案，降低了并行編程的難度和復雜度，靈活性強易使用。目前支持多線程開發(fā)工具OpenMP的有：MicrosoftVisualStudio、IntelC++編譯器。OpenMP在并行在細粒度與粗粒度技術[9]上均具有很高的計算效率。這里采用一階聲波方程進行交錯網(wǎng)格高階差分方法模擬計算均勻層狀介質波場，通過OpenMP并行化，測試不同模型下的不同粒度的并行效率。

1 一階聲波方程交錯網(wǎng)格高階差分方法

1.1 一階聲波方程

在不考慮體應力的情況下，二維一階聲波方程具體形式為式(1)：

(1)

其中：K=pv2為體積模量；v為縱波速度；vx為x方向上的速度分量；vz為z方向上的速度分量；p為聲壓；ρ為密度。

1.2 高階有限差分原理

空間上2N階差分，運用導數(shù)的高階展式，

(2)

(3)

時間上2M階差分原理如公式(4)所示：

(4)

其中:Δt為時間步長。

1.3 邊界吸收

有限差分數(shù)值模擬中常采用完全匹配層(PML)邊界吸收條件[10]，即在邊界處加一匹配層，采用衰減效應的波動方程，通過改變阻尼因子來控制邊界反射吸收的效果。常規(guī)的PML吸收條件針對一階波動方程，將波場分裂成沿著垂直和平行于傳播方向的兩組，在人工邊界處使得平行界面法向傳播的平面波快速衰減，對其他方向波場不衰減，從而邊界吸收的目的。對于一階聲波方程，只對聲壓分量進行分裂，因為速度分量是聲壓分量對相應的坐標軸的一階空間導數(shù)。聲壓可以分解為px、pz兩個部分，如公式(5)所示。

p=px+pz

(5)

對應的邊界吸收方程為公式(6)。

(6)

速度邊界吸收方程為公式(7)：

(7)

其中： d(x)=-1.5vmaxlog(R)x2/δ3； d(z)=-1.5vmaxlog(R)z2/δ3；vmax為最大縱波速度；δ為匹配層厚度；R為理想狀態(tài)下介質的反射系數(shù)。

1.4 穩(wěn)定性條件

二維交錯網(wǎng)格高階有限差分格式的穩(wěn)定性條件[11]，可由傅里葉分析方法可得到，其公式為式(8)。

(8)

式(8)中：vmax為最大縱波速度；Δt為時間采樣間隔；Δx、Δz分別為x、z兩個方向的網(wǎng)格間距。

2 基于OpenMP的并行化

2.1 并行算法流程

我們通過運用一階速度-應力彈性波方程為波動方程，使用時間上二階，空間八階的交錯網(wǎng)格，采用阻尼衰減對網(wǎng)格邊界處理，對水平層狀介質進行波場計算，并行程序的算法流程如圖1所示。

2.2 程序算法正確性對比

通過不同模型的串行和并行算法分別得到波場記錄(圖2、圖3)。圖2是模型(800×1 200×800)的串行和細粒度并行結果對比；圖3是模型(800×1 200×2 000)的串行和粗粒度并行結果對比。

由圖2、圖3可以看出，無論在細粒度還是粗粒度，并行算法獲得與串行算法完全相同的計算結果。

圖1 波場并行算法流程圖

圖2 細粒度串、并行結果

圖3 粗粒度串、并行結果

2.3 并行效率分析

通過進行多組模型的測試，從加速比和并行效率兩方面來討論并行算法的效率。測試平臺為CPU 4核心, 主頻1.5 GHz，集成開發(fā)環(huán)境Visual Studio 2010。

從表1可以看出，在細粒度并行情況下(當核心數(shù)目為1時，為串行結果)，隨著核心數(shù)目的增加，加速比有較大的變化，當核心數(shù)目為3時，加速比最大，加速效率在兩倍以上；隨著模型的改變，加速比總體呈現(xiàn)增加的趨勢。

從表2可以看出，在粗粒度并行情況下，隨著核心數(shù)目的增加，加速比有大的變化，但總體上當核心數(shù)目為3時，加速比最大，且比細粒度下的加速比大；隨著模型的改變，加速比并沒有體現(xiàn)遞增趨勢。

表1 細粒度并行不同模型的測試結果

表2 粗粒度并行不同模型的測試結果

3 結論

采用了基于多核多線程技術，實現(xiàn)了基于OpenMP的一階聲波方程高階有限差分數(shù)值并行模擬計算。在保證計算結果正確可靠的前提下，獲得了較高的并行計算效率。作為一種簡單又實用的技術，OpenMP具有可以跨平臺、跨語言使用的特點，在波場模擬及一系列計算比較大量數(shù)據(jù)預算中使用，可以提高效率，從而縮短計算時間。

這里基于均勻層狀介質的波場模擬計算，有望于在復雜介質中應用OpenMP技術探究其加速效果；另外，在地震記錄中，往往多炮記錄間有并行化處理，但在單炮記錄中少見并行算法，所以在地震單炮記錄的數(shù)據(jù)處理中有一定的應用空間。

[1] MADARIAGA R. Dynamics of an expanding circular fault[J]. Bulletin of the Seismological Society of America, 1976, 66(3): 639-666.

[2] VIRIEUX J. SH-wave propagation in heterogeneous media: velocity-stress finite-difference method[J]. Geophysics, 1984, 49(11): 1933-1942.

[3] VIRIEUX J. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method[J]. Geophysics, 1986, 51(4): 889-901.

[4] LEVANDER A R. Fourth-order finite-difference P-SV seismograms[J]. Geophysics, 1988, 53(11): 1425-1436.

[5] DABLAIN M A. The application of high-order differencing to the scalar wave equation[J]. Geophysics, 1986, 51 (1): 54-66.

[6] CRASE E.High-order (space and time) finite-difference modeling of the elastic wave equation. 60th SEG meeting, San Francisco, USA[J]. Expanded Abstracts, 1990,987-991.

[7] 董良國, 馬在田, 曹景忠, 等. 一階彈性波方程交錯網(wǎng)格高階差分解法[J]. 地球物理學報, 2000, 43(3): 411-419. DONG L G,MA Z T,CAO J ZH,et al. The staggered-grid high-order difference method of one-order elastic equation. Geophys, 2000,43(3):411-419.(In Chinese)

[8] 劉洋, 李承楚, 牟永光. 任意偶數(shù)階精度有限差分法數(shù)值模擬[J]. 石油地球物理勘探, 1998,33(1):1-11. LIU Y, LI CH CH, MOU Y G. Finite-difference numerical modeling of any even-order accuracy.OGP, 1998,33(1):1-11.(In Chinese)

[9] 白洪濤. 基于GPU的高性能并行算法研究[D].長春:吉林大學, 2010. BAI H T. Research on high performance parallel algorithms based on GPU. Changchun: Doctoral Dissertation of Jilin University, 2010.(In Chinese)

[10]王永剛, 邢文軍, 謝萬學, 等. 完全匹配層吸收邊界條件的研究[J]. 中國石油大學學報: 自然科學版, 2007, 31(1): 19-24. WANG Y G, XING W J, XIE W X, et al. Study of absorbing boundary condition by perfectly matched layer. Journal of China university of petroleum:natural science, 2007, 31(1): 19-24.(In Chinese)

[11]董良國, 馬在田, 曹景忠. 一階彈性波方程交錯網(wǎng)格高階差分解法穩(wěn)定性研究[J]. 地球物理學報, 2000, 43(6): 856-864. DONG L G, MA Z T, CAO J ZH. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation. Geophysics, 2000, 43(6): 856-864.(In Chinese)

[12]黃易, 師學明, 范建柯, 等. 并行計算技術及其在勘探地球物理學中的現(xiàn)狀與展望[J]. 地球物理學進展, 2010,25(2): 642-649. HUANG Y, SHI X M, FAN J K, et al. Review on parallel computing and its application in exploration geophysics. Progress in Geophys, 2010, 25(2): 642-649.(In Chinese)

[13]杜增利, 李亞林, 尹成, 等. 虛譜法一階應力-速度方程地震數(shù)值模擬[J]. 石油地球物理勘探, 2009, 44(5): 637-641. DU Z L,LI Y L,YIN CH, et al. Pseudo-spectral method seismic numerical modeling of first-order stress-velocity equation. OGP, 2009, 44(5): 637-641.(In Chinese)

[14]牟永光, 裴正林. 三維復雜介質地震數(shù)值模擬[M].北京:石油工業(yè)出版社, 2005. MOU Y G, PEI ZH L. Seismic numerical modeling for 3-D complex media. Beijing: Petroleum Industry Press, 2005.(In Chinese)

[15]潘海濱. 交錯網(wǎng)格地震波場模擬及頻散校正策略[J]. 物探化探計算技術, 2009, 31(4): 369-373. PAN H B. Staggered-grid seism wave IC wave simulation and the correction of numerical disperse. Computing Techniques for Geophysical and Geochemical Exploration, 2009, 31(4): 369-373.(In Chinese)

[16]李振春, 張慧, 張華. 一階彈性波方程的變網(wǎng)格高階有限差分數(shù)值模擬[J]. 石油地球物理勘探,2008.43(6):711-718 LI ZH CH, ZHANG H, ZHANG H. Variable grid high order finite difference numeric simulation of first order elastic wave equation. OGP, 2008,43(6):711-718.(In Chinese)

[17]寧剛, 熊章強, 陳持遜. 波動方程有限差分正演模擬誤差來源分析[J]. 物探與化探, 2008, 32(2): 203-206. NING G, XIONG ZH Q, CHEN CH X. An analysis of the error source in the wave propagation forward numerical simulation. Geophysical and Geochemical Exploration, 2008, 32(2): 203-206.(In Chinese)

A parallel algorithm of wave field using first-order acousticwave equation based on OpenMP

FANG Jin-wei1a， BAI Hong-tao1b*, SUN Hui-min1a

(1. Jilin University, a.College of Geo-Exploration Science and Technology,b.Center for Computer Fundamental Education, Changchun 130012, China)

The implementation of wave field forward modeling parallel computing algorithm is achieved using first-order acoustic wave equation based on OpenMP in this paper, aiming to getting a higher efficiency in comparison with that was in technique of the high-order staggered grid. The relationship between the number of cores and parallel efficiency is analyzed in changing models. The experiments demonstrate that our parallel algorithm can not only acquire the effective data accuracy, but also gain at least two times than the serial version in the multi-core CPU of commercial computer.

first-order acoustic wave equation; high-order staggered grid; OpenMP; parallel efficiency

2014-09-15改回日期：2014-11-25

國家自然科學基金(61272208)；中央高校基本科研業(yè)務費專項資金(JCKY-QKJC47,JCKY-QKJC49)；吉林省科技發(fā)展計劃項目基金(201101039)

方金偉(1991-)，男，本科，主要研究波場正演模擬和參數(shù)反演，E-mail：fangjw2311@mails.jlu.edu.cn。

*通信作者：白洪濤(1975-)，男，博士后，主要從事高性能計算研究，E-mail：baiht@jlu.edu.cn。

1001-1749(2015)05-0622-06

P 631.4

A

10.3969/j.issn.1001-1749.2015.05.13