一種求解地震波方程的高效并行譜元格式

林　燈　崔　濤　冷　偉　張林波

(中國科學(xué)院數(shù)學(xué)與系統(tǒng)科學(xué)研究院　北京　100190)(科學(xué)與工程計(jì)算國家重點(diǎn)實(shí)驗(yàn)室(中國科學(xué)院數(shù)學(xué)與系統(tǒng)科學(xué)研究院)　北京　100190)(lind@lsec.cc.ac.cn)

一種求解地震波方程的高效并行譜元格式

林燈崔濤冷偉張林波

(中國科學(xué)院數(shù)學(xué)與系統(tǒng)科學(xué)研究院北京100190)(科學(xué)與工程計(jì)算國家重點(diǎn)實(shí)驗(yàn)室(中國科學(xué)院數(shù)學(xué)與系統(tǒng)科學(xué)研究院)北京100190)(lind@lsec.cc.ac.cn)

摘要地震波數(shù)值模擬在地震學(xué)和地震勘探中扮演著非常重要角色.在已有工作的基礎(chǔ)上，提出1種高效并行的地震波PML方程譜元格式.PML被引入地震波方程以吸收外向波進(jìn)而模擬無界區(qū)域.進(jìn)一步，為了適應(yīng)復(fù)雜地形同時(shí)允許時(shí)間顯式推進(jìn)，譜元方法被用來離散地震波PML方程.由此得到地震波PML方程譜元格式.在此基礎(chǔ)上，闡述了單元?jiǎng)偠染仃嚪纸庑再|(zhì)，并說明了利用單元?jiǎng)偠染仃嚪纸饪梢源蠓鶞p少剛度矩陣存儲(chǔ)量同時(shí)顯著加速剛度矩陣與向量乘積，進(jìn)而顯著減少格式的計(jì)算量和存儲(chǔ)量.此外，算法復(fù)雜性分析表明格式無論在計(jì)算量上還是在存儲(chǔ)量上都優(yōu)于幾種已知的1階地震波PML方程譜元格式.結(jié)合并行技術(shù)，給出了高效并行的地震波PML方程譜元格式.數(shù)值實(shí)驗(yàn)驗(yàn)證了格式的正確性、良好的強(qiáng)弱并行可擴(kuò)展性以及對(duì)復(fù)雜地形的適應(yīng)性.

關(guān)鍵詞地震波方程；數(shù)值模擬；完美匹配層；譜元方法；單元?jiǎng)偠染仃嚪纸猓徊⑿杏?jì)算

譜元方法(spectral element method, SEM)最初應(yīng)用于計(jì)算流體力學(xué)[1-2]，隨后被成功引入地震波數(shù)值模擬領(lǐng)域.其中，Legendre譜元方法離散產(chǎn)生對(duì)角質(zhì)量矩陣，使得時(shí)域地震波數(shù)值模擬得以顯式推進(jìn)，顯著減少了計(jì)算開銷；同時(shí)，譜元方法的高精度、低頻散、優(yōu)良并行性能、適應(yīng)復(fù)雜地形等優(yōu)勢也得以充分展現(xiàn)[3-7].近年來，譜元方法被廣泛應(yīng)用于全球尺度的地震模擬[8-14]及地區(qū)尺度的地震模擬[15-18].

速度-應(yīng)力(應(yīng)變)格式是地震波數(shù)值模擬的常用格式.使用速度-應(yīng)力(應(yīng)變)格式可以便捷地結(jié)合完美匹配層(perfectly matched layer, PML)[19-20].PML是1種吸收邊界層，用以在地震波數(shù)值模擬中吸收外向波，模擬無界區(qū)域.PML具有2個(gè)重要性質(zhì)：1)波場振幅在PML中指數(shù)衰減；2)在離散前對(duì)任意入射角及任意頻率都有反射系數(shù)為零[21].憑借著良好的吸收效果，PML得到了極大的發(fā)展和廣泛的應(yīng)用[22-25]，也應(yīng)運(yùn)而生了許多變種及改進(jìn)[19-20,22,26-30].

本文基于地震波方程速度-應(yīng)變-應(yīng)力形式，應(yīng)用PML變換及譜元離散，得到地震波PML方程譜元格式.進(jìn)一步，提出利用單元?jiǎng)偠染仃嚪纸獠呗詢?yōu)化剛度矩陣(轉(zhuǎn)置)與向量乘積，有效減少了格式的計(jì)算量與存儲(chǔ)量.

1地震波方程與PML方法

1.1地震波方程

地震波方程刻畫了地震波在地層中傳播的基本規(guī)律，是地震波數(shù)值模擬的基礎(chǔ)，在天然地震學(xué)及地震勘探等領(lǐng)域扮演著重要的角色.

遵循文獻(xiàn)[31]，在線性、無粘及理想彈性假設(shè)下，地震波方程可寫成一階(時(shí)間)系統(tǒng)：

(1)

其中,v是速度，σ是應(yīng)力張量(2階對(duì)稱張量)，f為震源(外力)，ρ為密度，C是彈性張量(4階張量).在各向同性假設(shè)下，彈性張量可以表寫為

Cijkl=λδijδkl+μ(δikδjl+δilδjk)，

(2)

其中λ,μ是Lame常數(shù)，

借鑒文獻(xiàn)[19]，將一階(時(shí)間)系統(tǒng)改寫成分量形式：

(3)

其中，x∈D，

為了確保方程適定性，為1階(時(shí)間)系統(tǒng)附加初值條件：

(4)

其中，x∈D，同時(shí)為其附加自由邊界條件(也稱零載荷邊界條件)：

σ(x,t)n(x)=0, x∈?D.

(5)

1.2PML方法

PML方法是目前最為有效的吸收邊界層方法.在地震波數(shù)值模擬中，PML方法常用于吸收外向波，模擬無窮遠(yuǎn)邊界條件，進(jìn)而將無界問題轉(zhuǎn)化為有界問題求解.

PML可以解釋為復(fù)坐標(biāo)拉伸變換[22]:

(6)

等價(jià)微分形式為

(7)

亦即:

(8)

這里i為虛數(shù)單位，ω為頻率，τk在非PML區(qū)域恒為0，在PML區(qū)域定義為[32]

(9)

其中,L為PML厚度，c為最快波速，r為理論反射系數(shù)，xk0為沿xk正負(fù)方向PML起始坐標(biāo).計(jì)算區(qū)域D由PML區(qū)域Dm和非PML區(qū)域Dn m構(gòu)成，如圖1所示：

Fig. 1　Components of computational domain.圖1　計(jì)算區(qū)域組成

1.3地震波PML方程

我們把PML方法應(yīng)用于地震波方程，導(dǎo)出地震波PML方程.

遵循文獻(xiàn)[19]，對(duì)式(3)進(jìn)行Fourier變換及復(fù)坐標(biāo)拉伸變換，有:

(10)

其中，x∈Dm.引入輔助變量：

(11)

其中:

T≡diag{τ1,τ2,τ3},

從而有:

(12)

這里e=(1,1,1)T，進(jìn)而有:

(13)

(14)

其中，x∈Dm.值得注意的是，式(14)在非PML區(qū)域退化為式(3).

2地震波PML方程譜元格式

2.1譜元離散

(15)

Fig. 2　Invertible map Fi.圖2　可逆映射Fi

定義Ji為Fi的Jacobian矩陣，Ji為Ji的行列式,速度場有限元空間Vr定義為

(16)

應(yīng)力場有限元空間Σr定義為[20]

(17)

(18)

進(jìn)一步，通過張量積定義:

(19)

其中,p=(p1,p2,p3),pi=0,1,…,r；ed為第d個(gè)分量為1的單位向量.Vr的基函數(shù)滿足:

其中,q(k,i)表示q依賴于k及i.在此基礎(chǔ)上，對(duì)變分形式進(jìn)行有限元離散，有:

(20)

注意到一方面自由度定義在Gauss-Lobatto點(diǎn)，另一方面采用Gauss-Lobatto數(shù)值積分公式近似積分，我們有結(jié)論：

1) G是對(duì)角矩陣.

2) B是對(duì)角矩陣，且滿足:

其中,X=Aij,T,T′.

4) Rσi與RTVj可以寫為

Rσi=(Reσi|e)e→g，

綜上所述，我們得到地震波PML方程譜元格式(自由度形式)：

(21)

2.2格式優(yōu)化

式(21)的計(jì)算量高度集中于矩陣(轉(zhuǎn)置)乘向量Rσi與RTVj，存儲(chǔ)量主要集中于矩陣R.也就是說，如果格式該部分的計(jì)算量和存儲(chǔ)量大幅減少，格式整體的計(jì)算量和存儲(chǔ)量就會(huì)顯著減少.通過單元?jiǎng)偠染仃嚪纸鈨?yōu)化式(21)，注意到單元?jiǎng)偠染仃嚲哂芯仃嚪纸猓?/p>

(22)

其中:

KT∈3n×n，

利用單元?jiǎng)偠染仃嚪纸庥?jì)算剛度矩陣(轉(zhuǎn)置)與向量乘積可以簡潔表達(dá)為

(23)

以下我們說明，利用單元?jiǎng)偠染仃嚪纸庥?jì)算剛度矩陣(轉(zhuǎn)置)與向量乘積，無論在存儲(chǔ)量上還是在計(jì)算量上，都顯著優(yōu)于直接計(jì)算剛度矩陣(轉(zhuǎn)置)與向量乘積.

基于地震波PML方程譜元格式，結(jié)合單元?jiǎng)偠染仃嚪纸鈨?yōu)化策略，我們得到:

(24)

需要指出的是，我們可以通過變量替換進(jìn)一步消去質(zhì)量矩陣B，使得計(jì)算量和存儲(chǔ)量進(jìn)一步略微減少.

以下是3種地震波PML方程譜元格式在3維非結(jié)構(gòu)網(wǎng)格下內(nèi)存開銷和單個(gè)時(shí)間步計(jì)算量比較.可以看到，無論在計(jì)算量還是存儲(chǔ)量上，式(24)都顯著優(yōu)于其他2種格式.

Table 1　Storage and Calculation of Spectral Element Schemes for 3D Seismic Wave Equation with PML

Notes：nais the number of PML cells;nogives the amount of non-PML cells;ncis the size of Vi, which satisfiesnc<(r+1)3(na+no).

2.3格式并行化

基于區(qū)域分解及消息傳遞技術(shù)[33]，我們給出式(24)并行算法.

算法1. 式(24)并行算法.

將網(wǎng)格劃分成若干子網(wǎng)格分發(fā)到各進(jìn)程上;將震源及接收點(diǎn)分發(fā)給所在進(jìn)程;計(jì)算K，遍歷單元計(jì)算Me、組裝計(jì)算G、計(jì)算Ce,T,T′等(僅限PML單元).

t=0；

Whilet

For All Elements

If (PML elements)

Else

End If

End For All Elements

將非本地自由度數(shù)據(jù)發(fā)送給鄰居進(jìn)程；

計(jì)算Vj；

從鄰居進(jìn)程獲取非本地自由度數(shù)據(jù)；

t=t+Δt；

End While

由于該并行算法的計(jì)算量與自由度個(gè)數(shù)成線性關(guān)系且MPI通信僅涉及鄰居進(jìn)程，因此該并行算法具有良好的并行可擴(kuò)展性.

3數(shù)值實(shí)驗(yàn)

本節(jié)算例均在中國科學(xué)院數(shù)學(xué)與系統(tǒng)科學(xué)研究院科學(xué)和工程計(jì)算國家重點(diǎn)實(shí)驗(yàn)室LSSC-Ⅲ機(jī)群上測試.LSSC-Ⅲ機(jī)群擁有282個(gè)計(jì)算刀片，每個(gè)刀片包含2顆Intel X5550處理器和24 GB內(nèi)存，單核雙精度浮點(diǎn)峰值性能為10.68GFLOPS，282個(gè)計(jì)算結(jié)點(diǎn)的總浮點(diǎn)峰值性能為24.09GFLOPS.所有結(jié)點(diǎn)通過DDR Infiniband網(wǎng)絡(luò)互聯(lián).程序部分依賴于3維并行自適應(yīng)有限元平臺(tái)PHG[34].

在地震波數(shù)值模擬中，必須約束時(shí)空采樣以抑制時(shí)空頻散.遵循文獻(xiàn)[35]，程序采用8階譜元作空間離散并讓空間采樣不低于每個(gè)最短波長5個(gè)采樣點(diǎn)，同時(shí)采用2階leapfrog格式作時(shí)間離散并讓時(shí)間采樣不低于每個(gè)最短周期12個(gè)采樣點(diǎn).注意到式(24)是顯格式，空間步長與時(shí)間步長還受CFL條件約束[36].

3.1驗(yàn)證格式正確性

本節(jié)通過解析解與數(shù)值解的比較驗(yàn)證式(24)的正確性.考慮2維算例如圖3所示.每個(gè)方向64個(gè)單元，共4 096個(gè)單元.時(shí)間步長為0.4×10-3，推進(jìn)步數(shù)為5 000.采用P波震源如下：

其中,δ是delta函數(shù)，I為單位矩陣，

Fig. 3　Specification of the 2D benchmark.圖3　2維算例設(shè)定

震源主頻f0=10 Hz，震源時(shí)移t0=0.12 s，震源位于(0,0.375)，比例因子m=1.接收點(diǎn)位于(0.61,-0.61).PML厚度為0.25，理論反射系數(shù)為0.000 1.計(jì)算結(jié)果如圖4、圖5所示.圖4給出了0.6 s時(shí)刻速度場垂直分量快照.

其中，透射P波(a)、透射P-to-S波(b)、反射P-to-S波(c)、反射P波(d)及直達(dá)P波(e)清晰可見.圖5給出了速度場垂直分量在接收點(diǎn)(0.61,-0.61)處的解析解及數(shù)值解.解析解由Gar6more2D[37]計(jì)算得到,可以看到，解析解與數(shù)值解吻合得很好.

Fig. 4　Snapshot of vertical component of velocity field at 0.6 s.圖4　0.6 s時(shí)刻速度場垂直分量快照

Fig. 5　Comparison between the analytic and numerical solutions of vertical component of velocity field at the receiver.圖5　速度場豎直分量在接收點(diǎn)處解析解與數(shù)值解比較

3.2驗(yàn)證格式性能優(yōu)勢

本節(jié)通過比較在同一數(shù)值算例下3種地震波PML方程譜元格式的內(nèi)存開銷和計(jì)算時(shí)間驗(yàn)證式(24)的性能優(yōu)勢.考慮模型問題如下：計(jì)算區(qū)域?yàn)閇-0.5,0.5]3，每個(gè)方向16個(gè)單元，共4 096個(gè)單元.時(shí)間步長為0.75×10-3，推進(jìn)步數(shù)為1 000.采用速度震源如下：

其中,δ及S的定義見3.1節(jié).震源主頻f0=10 Hz，震源時(shí)移t0=0.12 s，震源位于區(qū)域中心，比例因子m為1.區(qū)域?yàn)閱我痪鶆蚪橘|(zhì)，介質(zhì)密度為1，P波速度為1.732，S波速度為1.PML厚度為0.25，理論反射系數(shù)為0.000 1.表2比較了3種地震波PML方程譜元格式的內(nèi)存開銷和計(jì)算時(shí)間.從表2可以看出，式(24)確實(shí)在計(jì)算量和存儲(chǔ)量上顯著優(yōu)于其他2種格式.

Table 2　　　Performance of Spectral Element Schemes for

3.3格式弱可擴(kuò)展性測試

本節(jié)對(duì)式(24)進(jìn)行弱可擴(kuò)展性測試.考慮模型問題如下：網(wǎng)格單元均為邊長為0.062 5的立方體，時(shí)間步長為0.000 75，推進(jìn)步數(shù)為1 000.震源形式與3.2節(jié)相同，震源主頻f0=10 Hz，震源時(shí)移t0=0.12 s，震源位于區(qū)域中心，比例因子為1.區(qū)域?yàn)閱我痪鶆蚪橘|(zhì)，介質(zhì)密度為1，P波速度為1.732，S波速度為1,不設(shè)置PML,網(wǎng)格規(guī)模隨進(jìn)程數(shù)變化而變化，確保每個(gè)進(jìn)程所擁有的網(wǎng)格單元個(gè)數(shù)近似相等，從而，每個(gè)進(jìn)程所擁有的自由度個(gè)數(shù)大體相同.計(jì)算結(jié)果如表3所示.從表3可以看出，式(24)具有良好的弱可擴(kuò)展性.

Table 3　Weak Scalability Test of Formula(24)

3.4格式強(qiáng)可擴(kuò)展性測試

本節(jié)對(duì)式(24)進(jìn)行強(qiáng)可擴(kuò)展性測試.考慮模型問題如下：計(jì)算區(qū)域?yàn)閇-2,2]3，每個(gè)方向加密6次，共有262 144個(gè)單元，譜元階數(shù)為8，時(shí)間步長為0.000 75，推進(jìn)步數(shù)為1 000.震源設(shè)定及介質(zhì)設(shè)定均與3.2節(jié)相同.速度分量自由度個(gè)數(shù)為135 005 697，應(yīng)力分量自由度個(gè)數(shù)為573 308 928.表4給出了采用不同進(jìn)程數(shù)進(jìn)行計(jì)算的推進(jìn)用時(shí)及并行效率.這里，并行效率定義為

其中,s是參與計(jì)算的進(jìn)程數(shù).之所以與16個(gè)進(jìn)程的計(jì)算結(jié)果比較是由于算例內(nèi)存開銷較大.從表4可以看出，式(24)具有良好的強(qiáng)可擴(kuò)展性.

Table 4　Strong Scalability Test of Formula(24)

Fig. 6　Computational domain and its unstructured hexahedral mesh圖6　計(jì)算區(qū)域及其非結(jié)構(gòu)六面體網(wǎng)格

3.5起伏地表

式(24)能夠有效模擬地震波在實(shí)際地形中傳播.考慮計(jì)算問題如下：計(jì)算區(qū)域?yàn)?0.91°N～31.11°N，103.32°E～103.52°E，地下深度12.5 km.Digital Elevation Model數(shù)據(jù)來源于Shuttle Radar Topographic Mission，如圖6所示.時(shí)間步長Δt=0.002 s，時(shí)間步數(shù)為4 000.震源形式與3.2節(jié)相同.震源主頻f0=2 Hz，震源時(shí)移t0=0.6 s，震源位于31.01°N，103.42°E，震源深度為8 km，比例因子m=108.區(qū)域?yàn)閱我痪鶆蚪橘|(zhì)，密度ρ=3.0 gcm3，P波速度為3.5 kms，S波速度為2.02 kms.PML厚度L=2.5 km，理論反射系數(shù)r=0.000 1.速度分量自由度個(gè)數(shù)為26 651 025，應(yīng)力分量自由度個(gè)數(shù)為37 558 080.圖7顯示了4.8 s時(shí)刻速度豎直分量波場圖.

Fig. 7　Snapshot of vertical component of velocity field at 4.8 s.圖7　4.8 s時(shí)刻速度場豎直分量快照

4結(jié)束語

本文從地震波方程速度-應(yīng)變-應(yīng)力形式出發(fā)，應(yīng)用PML方法導(dǎo)出地震波PML方程.進(jìn)一步，應(yīng)用譜元方法得到地震波PML方程譜元格式.在此基礎(chǔ)上，提出利用單元?jiǎng)偠染仃嚪纸鈨?yōu)化剛度矩陣(轉(zhuǎn)置)與向量乘積，避免了剛度矩陣存儲(chǔ)，有效減少了格式的計(jì)算量與存儲(chǔ)量.數(shù)值實(shí)驗(yàn)驗(yàn)證了該格式的正確性、性能優(yōu)勢及良好的并行強(qiáng)弱可擴(kuò)展性，并把該格式應(yīng)用于實(shí)際地形中.下一步我們將結(jié)合Local Time Stepping技術(shù)把該格式應(yīng)用于真實(shí)的大規(guī)模的地震波數(shù)值模擬.

參考文獻(xiàn)

[1]Patera A T. A spectral element method for fluid dynamics: Laminar flow in a channel expansion[J]. Journal of Computational Physics, 1984, 54(3): 468-488

[2]Maday Y, Patera A T. Spectral element methods for the incompressible Navier-Stokes equations[C]Proc of State-of-the-Art Surveys on Computational Mechanics. New York: ASME, 1989: 71-143

[3]Cohen G, Joly P, Tordjman N. Construction and analysis of higher-order finite elements with mass lumping for the wave equation[C]Proc of the 2nd Int Conf on Mathematical and Numerical Aspects of Wave Propagation. Philadelphia, PA: SIAM, 1993: 152-160

[4]Faccioli E, Maggio F, Paolucci R, et al. 2D and 3D elastic wave propagation by a pseudo- spectral domain decomposition method[J]. Journal of Seismology, 1997, 1(3): 237-251

[5]Komatitsch D. Spectral and spectral element methods for the 2D and 3D elastodynamics equations in heterogeneous media[D]. Paris: The Institute of Earth Physics of Paris, 1997

[6]Komatitsch D, Vilotte J P. The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures[J]. Bulletin of the Seismological Society of America, 1998, 88(2): 368-392

[7]Komatitsch D, Tromp J. Introduction to the spectral-element method for 3-D seismic wave propagation[J]. Geophysical Journal International, 1999, 139(3): 806-822

[8]Komatitsch D, Tromp J. Spectral-element simulations of global seismic wave propagation-I: Validation[J]. Geophysical Journal International, 2002, 149(2): 390-412

[9]Komatitsch D, Tromp J. Spectral-element simulations of global seismic wave propagation-Ⅱ: Three-dimensional models, oceans, rotation and self-gravitation[J]. Geophysical Journal International, 2002, 150(1): 303-318

[10]Komatitsch D, Ritsema J, Tromp J. The spectral-element method, Beowulf computing, and global seismology[J]. Science, 2002, 298(5599): 1737-1742

[11]Komatitsch D, Tsuboi S, Ji C, et al. A 14.6 billion degrees of freedom, 5 teraflops, 2.5 terabyte earthquake simulation on the Earth Simulator[C]Proc of the 2003 ACMIEEE Conf on Supercomputing. New York: ACM, 2003

[12]Chaljub E. Numerical modeling of the propagation of seismic waves in spherical geometry: Applications to global seismology[D]. Paris: University of Paris Ⅶ Denis Diderot, 2000

[13]Chaljub E, Capdeville Y, Vilotte J P. Solving elastodynamics in a fluid-solid heterogeneous sphere: A parallel spectral element approximation on non-conforming grids[J]. Journal of Computational Physics, 2003, 187(2): 457-491

[14]Chaljub E, Valette B. Spectral element modelling of three-dimensional wave propagation in a self-gravitating earth with an arbitrarily stratified outer core[J]. Geophysical Journal International, 2004, 158(1): 131-141

[15]Komatitsch D, Liu Q, Tromp J, et al. Simulations of ground motion in the Los Angeles basin based upon the spectral-element method[J]. Bulletin of the Seismological Society of America, 2004, 94(1): 187-206

[16]Liu Q, Polet J, Komatitsch D, et al. Spectral-element moment tensor inversions for earthquakes in southern California[J]. Bulletin of the Seismological Society of America, 2004, 94(5): 1748-1761

[17]Lee S J, Chen H W, Liu Q, et al. Three-dimensional simulations of seismic-wave propagation in the Taipei basin with realistic topography based upon the spectral-element method[J]. Bulletin of the Seismological Society of America, 2008, 98(1): 253-264

[18]Chaljub E, Moczo P, Tsuno S, et al. Quantitative comparison of four numerical predictions of 3D ground motion in the Grenoble Valley, France[J]. Bulletin of the Seismological Society of America, 2010, 100(4): 1427-1455

[19]Cohen G, Fauqueux S. Mixed spectral finite elements for the linear elasticity system in unbounded domains[J]. SIAM Journal on Scientific Computing, 2005, 26(3): 864-884

[20]Festa G, Vilotte J P. The Newmark scheme as velocity-stress time-staggering: An efficient PML implementation for spectral element simulations of elastodynamics[J]. Geophysical Journal International, 2005, 161(3): 789-812

[21]Berenger J P. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of Computational Physics, 1994, 114(2): 185-200

[22]Chew W C, Weedon W H. A 3D perfectly matched medium from modified Maxwell’s equations with stretched coordinates[J]. Microwave and Optical Technology Letters, 1994, 7(13): 599-604

[23]Chew W C, Liu Q H. Perfectly matched layers for elastodynamics: A new absorbing boundary condition[J]. Journal of Computational Acoustics, 1996, 4(4): 341-359

[24]Liu Q H, Tao J. The perfectly matched layer for acoustic waves in absorptive media[J]. The Journal of the Acoustical Society of America, 1997, 102(4): 2072-2082

[25]Zeng Y Q, He J Q, Liu Q H. The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media[J]. Geophysics, 2001, 66(4): 1258-1266

[26]Komatitsch D, Tromp J. A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation[J]. Geophysical Journal International, 2003, 154(1): 146-153

[27]Komatitsch D, Martin R. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[J]. Geophysics, 2007, 72(5): 155-167

[28]Meza-Fajardo K C, Papageorgiou A S. A nonconvolutional, split-field, perfectly matched layer for wave propagation in isotropic and anisotropic elastic media: Stability analysis[J]. Bulletin of the Seismological Society of America, 2008, 98(4): 1811-1836

[29]Martin R, Komatitsch D, Gedney S D, et al. A high-order time and space formulation of the unsplit perfectly matched layer for the seismic wave equation using auxiliary differential equations (ADE-PML)[J]. Computer Modeling in Engineering and Sciences, 2010, 56(1): 17-42

[30]Chen Z, Cui T, Zhang L. An adaptive anisotropic perfectly matched layer method for 3-D time harmonic electromagnetic scattering problems[J]. Numerische Mathematik, 2013, 125(4): 639-677

[31]Fichtner A. Full Seismic Waveform Modelling and Inversion[M]. Berlin: Spinger, 2010

[32]Collino F, Tsogka C. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[J]. Geophysics, 2001, 66(1): 294-307

[33]Zhang Linbo, Chi Xuebin, Mo Zeyao, et al. An Introduction to Parallel Computing[M]. Beijing: Tsinghua University Press, 2006 (in Chinese)(張林波, 遲學(xué)斌, 莫?jiǎng)t堯, 等. 并行計(jì)算導(dǎo)論[M]. 北京: 清華大學(xué)出版社, 2006)

[34]Zhang Linbo. Parallel hierarchical grid (PHG)[CPOL]. [2014-07-03]. http:lsec.cc.ac.cnphg

[35]Seriani G, Oliveira S P. Dispersion analysis of spectral element methods for elastic wave propagation[J]. Wave Motion, 2008, 45(6): 729-744

[36]Courant R, Friedrichs K, Lewy H. On the partial difference equations of mathematical physics[J]. IBM Journal of Research and Development, 1967, 11(2): 215-234

[37]Diaz J, Ezziani A. Gar6more2D[CPOL]. [2014-05-20]. http:gar6more2d.gforge.inria.fr

Lin Deng, born in 1989. PhD candidate in Academy of Mathematics and Systems Science, Chinese Academy of Sciences. His main research interests include parallel computing and seismic exploration.

Cui Tao, born in 1979. PhD and associate professor in Academy of Mathematics and Systems Science, Chinese Academy of Sciences. His main research interests include parallel computing and electromagnetic computing.

Leng Wei, born in 1984. PhD and assistant professor in Academy of Mathematics and Systems Science, Chinese Academy of Sciences. His main research interests include parallel computing and seismic exploration.

Zhang Linbo, born in 1962. Professor and PhD supervisor in Academy of Mathematics and Systems Science, Chinese Academy of Sciences. His main research interests include parallel computing especially 3D parallel (adaptive) finite element computing.

An Efficient Parallel Spectral Element Scheme for Solving Seismic Wave Equation

Lin Deng, Cui Tao, Leng Wei, and Zhang Linbo

(AcademyofMathematicsandSystemSciences,ChineseAcademyofSciences,Beijing100190)(StateKeyLaboratoryofScientificandEngineeringComputing(AcademyofMathematicsandSystemSciences,ChineseAcademyofSciences),Beijing100190)

AbstractNumerical simulation of seismic waves plays an essential role in seismology and seismic exploration. We propose here an efficient parallel spectral element scheme for seismic wave equation with perfectly matched layer (PML). PML is integrated into the seismic wave equation to absorb out-going waves and mimic unbounded domain. Ulteriorly, to enable adapting complex topography and explicit time stepping, the spectral element method (SEM) is used to discretize seismic wave equation with PML, which results in a spectral element scheme. In addition, we demonstrate that element stiffness matrices can be decomposed, which can be used to greatly reduce the storage of stiffness matrix and accelerate stiffness matrix-vector multiplication and thus remarkably speed up the scheme and cut down memory cost. Furthermore, we study several spectral element schemes known and show that our scheme is superior to others in both calculation and storage. Combined with parallel technique, an efficient parallel spectral element solver for seismic wave equation with PML is present. Numerical experiments show that our scheme is correct, well strongly?weakly scalable and of good adaptation to complex topography.

Key wordsseismic wave equation; numerical simulation; perfectly matched layer (PML); spectral element method (SEM); element stiffness matrices decomposition; parallel computing

收稿日期：2014-12-30；修回日期：2015-04-03

基金項(xiàng)目：國家“九七三”重點(diǎn)基礎(chǔ)研究發(fā)展計(jì)劃基金項(xiàng)目(2011CB309703);國家“八六三”高技術(shù)研究發(fā)展計(jì)劃基金項(xiàng)目(2012AA01A309);國家自然科學(xué)基金項(xiàng)目(11171334,11321061,11101417);中國科學(xué)院國家數(shù)學(xué)與交叉科學(xué)中心資助課題

中圖法分類號(hào)TP391

This work was supported by the National Basic Research Program of China (973 Program) (2011CB309703), the National High Technology Research and Development Program of China (863 Program) (2012AA01A309), the National Natural Science Foundation of China (11171334, 11321061, 11101417), and the Foundation of the National Center for Mathematics and Interdisciplinary Sciences, Chinese Academy of Sciences.