基于新離散公式和OpenMP優(yōu)化的有限差分聲波數(shù)值模擬

鄭如秋, 王波濤，馮永照, 余衛(wèi)江

中海油田服務(wù)股份有限公司 特普公司, 廣東 湛江 542051

0 引言

有限差分算法是偏微分方程的主要數(shù)值解法之一，是最早的數(shù)值模擬方法。該方法的基本原理是采用離散后的差分算子取代相應(yīng)的連續(xù)微分算子，具有計(jì)算效率高、占用內(nèi)存相對較小的優(yōu)勢，并且對于近遠(yuǎn)場及復(fù)雜邊界都有廣泛的適用性，能夠準(zhǔn)確地模擬波在各種介質(zhì)及復(fù)雜結(jié)構(gòu)地層中的傳播規(guī)律，是勘探地震學(xué)中應(yīng)用最廣泛的數(shù)值計(jì)算方法[1]。

目前，國內(nèi)外學(xué)者對基于波動(dòng)方程的地震波正演模擬技術(shù)進(jìn)行了大量的研究工作，地震波正演模擬技術(shù)十分成熟。但在模擬計(jì)算過程中由于傳統(tǒng)的彈性波動(dòng)方程數(shù)值模擬存在參數(shù)繁多、計(jì)算量大、耗時(shí)長和占用內(nèi)存空間大等問題，因而聲波方程數(shù)值模擬得到廣泛應(yīng)用[2]。1971年，Claerbout et al.將有限差分法在波場模擬中的應(yīng)用做了很大改進(jìn)，提高了模擬效率[3]；1974年，Alford et al.研究聲波方程數(shù)值模擬，給出精度與差分階數(shù)和網(wǎng)格步長關(guān)系[4]；1976年Madariaga首次提出交錯(cuò)網(wǎng)格差分算法[5]；1984年，Virieux et al.推導(dǎo)了一階應(yīng)力--速度方程有限差分算法，提高了計(jì)算精度[6]；1986年，Dablain et al.給出聲波方程交錯(cuò)網(wǎng)格算法差分格式[7]； 1997年， Clayton將吸收邊界應(yīng)用到聲波方程數(shù)值模擬當(dāng)中[8]； 2001年，Collino et al.將PML邊界應(yīng)用到應(yīng)力--速度方程的數(shù)值模擬當(dāng)中[9]；2004年，Saenger將交錯(cuò)網(wǎng)格進(jìn)一步發(fā)展，提出了旋轉(zhuǎn)交錯(cuò)網(wǎng)格法，這一方法對復(fù)雜介質(zhì)體的模擬精度有顯著提高[10]；2000年，董良國等詳細(xì)給出了基于交錯(cuò)網(wǎng)格的一階應(yīng)力--速度的彈性波方程有限差分法[11]；2004年，裴正林等進(jìn)一步研究了空間任意偶數(shù)階精度的交錯(cuò)網(wǎng)格差分法，使其更加適應(yīng)不同復(fù)雜介質(zhì)的數(shù)值模擬[12]。

通過對有限差分聲波波場數(shù)值模擬研究的分析發(fā)現(xiàn)，基于一階應(yīng)力--速度方程的有限差分聲波模擬是業(yè)界應(yīng)用最為廣泛的算法，很對學(xué)者對其進(jìn)行優(yōu)化，例如優(yōu)化差分系數(shù)的求解方法，使用旋轉(zhuǎn)網(wǎng)格代替常規(guī)網(wǎng)格。雖然這些方法可以提高模擬精度，不容易出現(xiàn)頻散現(xiàn)象，但是計(jì)算效率會降低，而在工業(yè)生產(chǎn)中計(jì)算效率更為重要。通常情況下可以通過使用高階有限差分算法來提高計(jì)算精度，克服頻散現(xiàn)象，因此，常規(guī)網(wǎng)格下有限差分算法是保證較高的模擬精度條件下，計(jì)算效率最高的算法。隨著計(jì)算機(jī)水平的發(fā)展，利用OpenMP并行計(jì)算可以大大提高數(shù)值模擬速度，但是常規(guī)網(wǎng)格下的有限差分算法的PML邊界公式復(fù)雜，包含三階偏導(dǎo)數(shù)，需要將計(jì)算區(qū)域分為有效計(jì)算區(qū)域和邊界區(qū)域，不利于實(shí)現(xiàn)并行計(jì)算，容易產(chǎn)生內(nèi)存錯(cuò)誤。對此，模仿一階應(yīng)力--速度方程PML邊界公式格式，改進(jìn)不含有三階偏導(dǎo)數(shù)倒數(shù)的PML邊界。保證了計(jì)算區(qū)域和邊界區(qū)域計(jì)算代碼統(tǒng)一，并行計(jì)算編程易于實(shí)現(xiàn)，對OpenMP的制導(dǎo)語句參數(shù)進(jìn)行測試優(yōu)化，進(jìn)一步提高了計(jì)算效率。

1 聲波方程有限差分算法原理

利用有限差分法進(jìn)行數(shù)值模擬，對空間采用高階差分離散，而在時(shí)間上只采用二階中心差分，推導(dǎo)出常規(guī)網(wǎng)格條件下計(jì)算區(qū)域與邊界區(qū)域一致的離散公式，使得OpenMP并行計(jì)算過程中效率提高，所以對于方程離散原理，差分系數(shù)計(jì)算，頻散問題不做詳細(xì)論述。

標(biāo)準(zhǔn)二維聲波方程可以表示為：

(1)

式中：P=P(x,z;t) 表示聲波波場；V=V(x,z)表示模型速度；ρ=ρ(x,z)表示模型密度。相應(yīng)的PML邊界條件下邊界區(qū)域計(jì)算方程是包含三階偏導(dǎo)數(shù)的復(fù)雜方程組，以右邊界為例[13]：

(2)

其中：

(3)

式中：L為PML區(qū)域?qū)挾?；R為反射系數(shù)；lx為該點(diǎn)到邊界的距離。通過公式(2)可以得出，對于模擬的邊界區(qū)域要分別判斷所屬區(qū)域，左、右邊界使用水平方向的衰減系數(shù)d(x)，而對于上、下邊界需要使用垂直方法的衰減系數(shù)d(z)，對于四個(gè)交點(diǎn)有的學(xué)者提出對以對角線進(jìn)行劃分成兩部分，上部分使用d(z)，下部分使用d(x)，但是這些都將導(dǎo)致需要對計(jì)算區(qū)域進(jìn)行判定，導(dǎo)致編程實(shí)現(xiàn)過程繁瑣，不利于并行計(jì)算。

通常情況下實(shí)現(xiàn)編程計(jì)算是通過對方程(1)進(jìn)行離散可以得到空間2N階精度的計(jì)算區(qū)域的波動(dòng)方程有限差分格式。

(4)

Collino提出將二階聲波方程轉(zhuǎn)換為一階應(yīng)力--速度方程后，再加入PML邊界。這個(gè)方法大大降低了直接對二階聲波方程加載PML邊界而增加的復(fù)雜度。一階應(yīng)力--速度方程加載PML邊界條件的控制方程具有如下形式：

(5)

其中vx和vz是質(zhì)點(diǎn)震動(dòng)的速度分量，離散公式如下：

(6)

一階應(yīng)力速度方程有限差分算法其相應(yīng)的PML邊界區(qū)域離散公式與計(jì)算區(qū)域離散公式相同，只有阻尼系數(shù)在邊界區(qū)域d(x)≠0,d(z)≠0 ，而計(jì)算區(qū)域阻尼系數(shù)為0。編程過程中不需要對模型進(jìn)行區(qū)域劃分，有利于OpenMP并行計(jì)算。但由于引入中間變量vx和vz，相比于常規(guī)網(wǎng)線條件下的有限差分算法，增加了大約兩倍的計(jì)算量。

目前關(guān)于數(shù)值模擬中PML邊界加載方式十分成熟，但是絕大部分都是基于一階應(yīng)力--速度方程，交錯(cuò)網(wǎng)格框架下的研究。筆者研究了一階應(yīng)力--速度方程PML邊界離散過程，SPML和NPML的推導(dǎo)原理，結(jié)合參考文獻(xiàn)[14--15]，引入二階阻尼系數(shù)，給出直接對常規(guī)網(wǎng)格下二階聲波方程加載PML邊界的方法：

(7)

(8)

(9)

Px項(xiàng)經(jīng)過數(shù)值模擬測試，不需要展開為 0.5×[(Px)k+1+(Px)k-1]，直接保留為原型，此時(shí)來自邊界反射能量最小。

最終獲得邊界反射能量最小的常規(guī)網(wǎng)格二階聲波方程離散公式為：

(10)

計(jì)算機(jī)實(shí)現(xiàn)代碼：

/*波動(dòng)方程計(jì)算函數(shù) 定義*/

void wavequ( int i, int j, int k, double *cm, double (*u1t0)[XN], double (*u2t0)[XN], double (*u1tp1)[XN],double (*u2tp1)[XN], double (*u1tn1)[XN], double (*u2tn1)[XN], double (*u1tmp)[XN], double (*u2tmp)[XN],double (*ut0)[XN], double (*utp1)[XN], double (*dx)[XN], double (*dz)[XN], double (*v)[XN])

{

int s;

double sumequ1,sumequ2;

for(s=0,sumequ1=0.0,sumequ2=0.0;s

{

sumequ1=sumequ1+cm[s]*(v[i][j]*v[i][j]*ut0[i][j+s+1]+v[i][j]*v[i][j]*ut0[i][j-s-1]-v[i][j]*v[i][j]*2*ut0[i][j]);

sumequ2=sumequ2+cm[s]*(v[i][j]*v[i][j]*ut0[i+s+1][j]+v[i][j]*v[i][j]*ut0[i-s-1][j]-v[i][j]*v[i][j]*2*ut0[i][j]);

}

u1tmp[i][j]=DT*DT/DH/DH*sumequ1;

u1tp1[i][j]=1/(1+2*dx[i][j]*DT)*((2-dx[i][j]*dx[i][j]*DT*DT-2*dx[i][j]*DT)*u1t0[i][j]- u1tn1[i][j]+u1tmp[i][j]);

u2tmp[i][j]=DT*DT/DH/DH*sumequ2;

u2tp1[i][j]=1/(1+2*dz[i][j]*DT)*((2-dz[i][j]*dz[i][j]*DT*DT-2*dz[i][j]*DT)*u2t0[i][j]-u2tn1[i][j]+u2tmp[i][j]);

utp1[i][j]=u1tp1[i][j]+u2tp1[i][j];

}

推導(dǎo)出來的離散有限差分公式(10)，在計(jì)算過程中計(jì)算區(qū)域與邊界區(qū)域計(jì)算機(jī)實(shí)現(xiàn)代碼統(tǒng)一，不需要對邊界進(jìn)行判斷而改變計(jì)算公式，更加有利于并行計(jì)算。推導(dǎo)的公式(10)與不含邊界的有限差分公式(4)相比，只是增加了阻尼系數(shù)的計(jì)算，這個(gè)阻尼系數(shù)在速度模型確定時(shí)候即可計(jì)算，不需要每次都特殊計(jì)算。

2 OpenMP多核并行算法

本文提出基于OpenMP及改進(jìn)的常規(guī)網(wǎng)格有限差分聲波數(shù)值模擬優(yōu)化，在第1章中主要通過改進(jìn)邊界條件，從算法原理上是代碼簡潔為并行計(jì)算奠定基礎(chǔ)，提高計(jì)算效率。本章主要對于OpenMP制導(dǎo)語句、子句等細(xì)節(jié)進(jìn)行調(diào)整，可以進(jìn)一步提高計(jì)算效率。

OpenMP標(biāo)準(zhǔn)通過定義編譯制導(dǎo)、庫例程和環(huán)境變量規(guī)范的方法，基于fork-jorn的并行執(zhí)行模型，將程序主體劃分為并行區(qū)和串行區(qū)，不同線程之間可以通過共享變量實(shí)現(xiàn)數(shù)據(jù)交換，多核CPU進(jìn)行并行運(yùn)算時(shí)，多個(gè)CPU可同時(shí)對程序各個(gè)不相關(guān)進(jìn)程進(jìn)行運(yùn)算，從而提高程序的運(yùn)行效率。OpenMP可直接在串行代碼上通過編寫并行制導(dǎo)語句實(shí)現(xiàn)程序的并行，以良好的簡潔性和可移植性成為并行編程的工業(yè)標(biāo)準(zhǔn)。需要并行化的程序如果比較復(fù)雜，需要注意程序內(nèi)部結(jié)構(gòu)和變量之間的邏輯關(guān)系，以防止并行后內(nèi)存讀寫沖突。

本章主要分析討論利用OpenMP技術(shù)進(jìn)行并行設(shè)計(jì)的兩個(gè)關(guān)鍵問題：一是如何將原來的串行程序改為并行程序，避免內(nèi)存讀寫沖突；二是如何更為有效地使用制導(dǎo)語句，提高并行效率。聲波波場數(shù)值模擬過程當(dāng)中，任意一點(diǎn)下一時(shí)刻的波場值的計(jì)算是通過當(dāng)前時(shí)刻的波場值、前一時(shí)刻的波場值和當(dāng)前時(shí)刻波場在空間方向的導(dǎo)數(shù)獲得的，這三個(gè)變量之間互不相關(guān)，計(jì)算過程中不會產(chǎn)生影響，可以共享至全部線路，具有可并行性。但是控制波場點(diǎn)位的變量i,j只能對應(yīng)于對應(yīng)線路計(jì)算的波場值，需要進(jìn)行私有變量處理。

OpenMP使用C編輯器的#pragma擴(kuò)展機(jī)制來定義制導(dǎo)，基本格式為：

#pragma omp directive-name[clause[[,]clause]…] new-line

其中，每個(gè)制導(dǎo)均以#pragma omp開始，directive-name是制導(dǎo)名，[]可以加入相應(yīng)功能的子句。OpenMP最基本的單元是parallel結(jié)構(gòu)，由parallel制導(dǎo)。使用OpenMP中private數(shù)據(jù)子句解決簡單變量的內(nèi)存讀寫沖突問題。private子句可將簡單變量聲明為本線程的私有變量，每個(gè)線程都存有變量的副本，其他線程無法訪問。即使在并行區(qū)域外有同名的共享變量，此共享變量也不會對并行區(qū)域造成影響，且并行區(qū)域內(nèi)部計(jì)算也不會改變外部共享的變量值。因此在并行起始語句中加入private(i,j)即可以解決簡單的內(nèi)存讀寫沖突問題。

OpenMP采用fork-join(分叉--合并)的并行執(zhí)行模式，在運(yùn)行過程中，遇到并行開始代碼(由parallel制導(dǎo))，計(jì)算機(jī)的主線程則調(diào)用多個(gè)線程來執(zhí)行并行任務(wù)。如果沒有nowait子句，所有線程在共享計(jì)算結(jié)構(gòu)的結(jié)束處隱式同步。為了提高并行計(jì)算效率，子句需要使用nowait，使參與計(jì)算的線程無需同步而繼續(xù)執(zhí)行隨后的程序段，減少多余的同步以提高計(jì)算效率。OpenMP的調(diào)度模式使用靜態(tài)調(diào)度schedule(static)，同時(shí)仍然需要使用num_threads()子句保證并行區(qū)域的占用線程不會影響主線程的計(jì)算效率，減少線程調(diào)度的時(shí)間消耗。通過測試以計(jì)算機(jī)總線程數(shù)為8為例，在并行時(shí)盡量不要使用全部線程數(shù)量，在并行制導(dǎo)語句中加入num_threads(7)，保證并行過程最多使用7個(gè)線程，以保證主線程的計(jì)算效率。

本文中的有限差分算法，計(jì)算每一時(shí)刻的波場快照需要通過for循環(huán)結(jié)構(gòu)完成，需要使用parallel for制導(dǎo)，OpenMP對程序的循環(huán)部分進(jìn)行并行化，要求用于并行的循環(huán)部分每次循環(huán)互不相關(guān)。二維模擬的計(jì)算任務(wù)分為三層循環(huán)(最外層是時(shí)間，內(nèi)層依次為x，z方向上的空間循環(huán))，由于每個(gè)時(shí)間切片上的當(dāng)前循環(huán)，依賴于上兩個(gè)時(shí)間切片上的計(jì)算結(jié)果，即時(shí)間循環(huán)中各個(gè)時(shí)間切片存在相關(guān)性，因此對于時(shí)間不適合使用OpenMP并行而對于x，z方向上的空間循環(huán)，具備應(yīng)用OpenMP并行計(jì)算條件，因此每一次空間循環(huán)都可以被指定進(jìn)行并行計(jì)算。

使用有限差分算法進(jìn)行地震波場聲波數(shù)值模擬區(qū)域時(shí)，需要在邊界加入特殊的邊界條件，以消除來自邊界的虛假反射問題。使用上一節(jié)推導(dǎo)出的有限差分公式，在空間循環(huán)外層加上并行制導(dǎo)指令即可。并行計(jì)算核心代碼結(jié)構(gòu)如下所示：

for(k=0;k

{

#pragma omp parallel num_threads(7)

{

#pragma omp for private(i,j) schedule(static) nowait

for(i=M;i

for(j=M;j

逐點(diǎn)計(jì)算的波動(dòng)方程計(jì)算函數(shù)代碼

}

}

波場轉(zhuǎn)存儲代碼

}

}

3 數(shù)值模擬計(jì)算

本章主要為了驗(yàn)證兩個(gè)問題：①改進(jìn)的離散的波場模擬公式在任何模型中可以準(zhǔn)確模擬波場；②結(jié)合OpenMP可以進(jìn)一步提高計(jì)算效率。目前一階應(yīng)力--速度方程時(shí)聲波模擬的標(biāo)準(zhǔn)方法，為了證明改進(jìn)的離散公式(10)可以準(zhǔn)確模擬波場，以下的試驗(yàn)將從波場模擬的準(zhǔn)確程度和計(jì)算效率與其對比。

試驗(yàn)1：

使用一個(gè)500×500，dx=dz=10 m，dt=800 ms，速度為4 000 m/s的均勻速度模型，震源放置在中心(250,250)。分別使用本文方法進(jìn)行模擬并與一階應(yīng)力--速度方程的模擬結(jié)果進(jìn)行對比。如圖1所示，在t=0.56 s，波場未到達(dá)邊界，模型為均勻模型，因此波場快照表現(xiàn)為圓形，兩種方法波形一致。在t=0.8 s時(shí)，波場觸碰到人工邊界，均為產(chǎn)生反射現(xiàn)象，說明改進(jìn)的離散計(jì)算公式可以正確模擬波場，并且在邊界處無明顯反射。值得注意的是，由于波場傳播的控制方程和模擬過程的差分格式差異，波場在能量振幅存在差異。

a. 本文方法(t=0.56 s);b. 本文方法(t=0.8 s);c. 一階應(yīng)力--速度方程(t=0.56 s); d. 一階應(yīng)力--速度方程(t=0.8 s)。圖1 簡單模型的波場快照Fig.1 Snapshots of wavefield of a simple model

為進(jìn)一步證明推導(dǎo)出的公式(10)在復(fù)雜模型中也可以準(zhǔn)確模擬波場，使用Sigsbee模型進(jìn)一步驗(yàn)證。模型大小為1 001×501，dx=dz=10 m，dt=0.001 s。速度模型和波場快照如圖2所示，sigsbee速度模型中包含一個(gè)高速鹽體，鹽體上表面起伏，容易產(chǎn)生繞射波場，將震源放置于模型中間(5 000,2 600)位置。理論上，震源上方的崎嶇界面會產(chǎn)生一系列繞射波，而鹽體的下界面由于界面平緩，界面處速度差異明顯，會產(chǎn)生強(qiáng)反射波和沿著界面滑行的折射波，折射波的波前表現(xiàn)為直線。選擇t=0.8 s時(shí)的波場快照進(jìn)行展示，兩種算法的波場快照中都可以看到震源激發(fā)的透射波場波前(綠色箭頭所指)、鹽體界面產(chǎn)生的反射波(藍(lán)色箭頭所指)、折射波(黑色箭頭所指)和繞射波(紅色箭頭所指)。通過對于sigsbee模型進(jìn)行模擬可以充分說明所推導(dǎo)的離散公式(10)可以準(zhǔn)確模擬波場。

圖2 Sigsbee速度模型Fig.2 Sigsbee velocity model

a.本文方法； b.一階應(yīng)力--速度方程。圖3 t=0.8 s時(shí)的波場快照Fig.3 Snapshots at t=0.8 s

試驗(yàn)2：

該實(shí)驗(yàn)中的程序均為并行程序，使用Linux系統(tǒng)機(jī)群的一個(gè)節(jié)點(diǎn)進(jìn)行并行計(jì)算效率測試對比。節(jié)點(diǎn)主頻2.4 GHz，十六線程。分別使用基于一階應(yīng)力--速度方程，本文方法公式(10)和并行代碼優(yōu)化后的程序進(jìn)行計(jì)算效率測試。測試使用上一個(gè)試驗(yàn)Sigsbee模型進(jìn)行模擬，總采樣點(diǎn)數(shù)設(shè)置為1 000，總炮數(shù)設(shè)置為1炮、20炮和100炮。測試結(jié)果如表1所示。

表1 計(jì)算效率對比表Table 1 Comparison of computation efficiency /s

從試驗(yàn)結(jié)果可以看出，本文改進(jìn)的有限差分公式(10)相比于一階應(yīng)力--速度方程的計(jì)算效率已經(jīng)有所提高。而對于并行代碼優(yōu)化后的就算效率進(jìn)一步提高，計(jì)算時(shí)間只需要一階應(yīng)力--速度方程的60%。

4 結(jié)論

(1)規(guī)則網(wǎng)格下的有限差分算法，適當(dāng)提高差分階數(shù)，既可以滿足模擬精度要求，相比于一階應(yīng)力--速度方程的計(jì)算效率更高。通過改進(jìn)邊界條件公式，推導(dǎo)出新的離散公式，相比于傳統(tǒng)二階聲波方程的PML邊界控制方程，邊界計(jì)算區(qū)域中不再包含波場的三階偏導(dǎo)數(shù)和阻尼系數(shù)的一階導(dǎo)數(shù)。并且，新的離散公式在計(jì)算區(qū)域與邊界區(qū)域公式一致，計(jì)算機(jī)實(shí)現(xiàn)過程中不需要加入條件語句判斷。相比于一階應(yīng)力--速度方程，減少中間變量引入，模擬的控制方程只需要3個(gè)，可以提高計(jì)算效率。

(2)通過結(jié)合OpenMP并行計(jì)算，可以進(jìn)一步提高計(jì)算效率。并行過程容易產(chǎn)生變量的內(nèi)存問題，串行代碼轉(zhuǎn)換為并行代碼后對程序變量優(yōu)化，減少內(nèi)存訪問。通過OpenMP制導(dǎo)語句中使用適當(dāng)?shù)膮?shù)，可以進(jìn)一步加快計(jì)算效率。