基于非結(jié)構(gòu)網(wǎng)格求解三維D′Alembert介質(zhì)中聲波方程的并行加權(quán)Runge-Kutta間斷有限元方法

賀茜君, 楊頂輝, 仇楚鈞, 周艷杰, 常蕓凡

1 北京工商大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院, 北京 100048 2 清華大學(xué)數(shù)學(xué)科學(xué)系, 北京 100084

0 引言

基于地震波動(dòng)方程的正演是計(jì)算地球物理學(xué)的重要研究?jī)?nèi)容，它不僅能為我們提供精確的波場(chǎng)模擬結(jié)果，而且也是基于波動(dòng)方程的地震勘探和地震學(xué)的重要基礎(chǔ)(牟永光和裴正林, 2005).在最近幾十年里，隨著計(jì)算機(jī)能力的快速提升，涌現(xiàn)了許許多多優(yōu)秀的數(shù)值算法，例如有限差分法(Finite difference method，簡(jiǎn)稱FDM)、有限元法(Finite element method，簡(jiǎn)稱FEM)、偽譜法(Pseudo-spectral method，簡(jiǎn)稱PSM)、譜元法(Spectral element method，簡(jiǎn)稱SEM)、間斷有限元法(Discontinuous Galerkin method，簡(jiǎn)稱DGM)等.這些算法均能滿足我們對(duì)于數(shù)值算法的部分要求，它們也有其獨(dú)特的優(yōu)勢(shì).對(duì)于數(shù)值求解3D波動(dòng)方程來說，這些算法都得到了廣泛應(yīng)用.FDM編程簡(jiǎn)單、計(jì)算速度快，且存儲(chǔ)量小，是求解3D波動(dòng)方程最為常用的一種數(shù)值算法(e.g.,董良國(guó)等, 2000; Moczo et al., 2000, 2002; Zhang and Liu, 2002; Kristek and Moczo, 2003; 牟永光和裴正林, 2005; Yang et al., 2007; 張金海等, 2007; Zhang and Gao, 2009, Liu and Sen, 2009; Yang and Wang, 2010; Zhang et al., 2012; Igel, 2017; Shragge and Tapley, 2017; Wang et al., 2019a, 2019b).FEM最早由Feng(1965)和西方學(xué)者獨(dú)立提出，它能處理復(fù)雜區(qū)域和邊界條件，在求解3D波動(dòng)方程也得到了一定的應(yīng)用(Galis et al., 2008; Moczo et al., 2011; Igel, 2017)，但是由于其計(jì)算量大且并行性差，作為單一方法用于求解波傳播問題已不多見，它常與有限差分方法結(jié)合形成的混合方法(e.g., Ma et al., 2004; Galis et al., 2008).PSM利用快速傅里葉變換(Fast Fourier transform，簡(jiǎn)稱FFT)來處理全部波場(chǎng)中的空間導(dǎo)數(shù), 精度很高，且數(shù)值頻散小，在求解波動(dòng)方程領(lǐng)域也得到了廣泛應(yīng)用(Furumura et al., 1998; Igel, 1999; Klin et al., 2010; Carcione et al., 2018)，但是3D情形的PSM需要利用全局?jǐn)?shù)據(jù)進(jìn)行FFT，因此并行性較差(Pelz, 1991).SEM是計(jì)算地球物理領(lǐng)域近十年來發(fā)展最為迅速的數(shù)值算法，它具有FEM的基本特征，在單元內(nèi)部的基函數(shù)基于特定積分節(jié)點(diǎn)——Gauss-Lobatto-Legendre(GLL)點(diǎn)，因此在使用GLL積分準(zhǔn)則的情況下質(zhì)量矩陣是對(duì)角陣，SEM的這種處理不僅使得方法的精度提高，而且并行性好，在3D石油勘探領(lǐng)域及大尺度模擬方面都有很好的表現(xiàn)(e.g., Komatitsch and Vilotte, 1998; Komatitsch and Tromp, 1999; Komatitsch et al., 1999; Tromp et al., 2008; Liu et al., 2017).

DGM也是近十年來迅速發(fā)展起來的一種數(shù)值算法，它最早是由Reed 和 Hill(1973)在求解中子輸運(yùn)方程時(shí)提出，其后，經(jīng)過Cockburn and Shu(1989, 2001)，Rivière 和 Wheeler(2003)等人的不斷發(fā)展，以基于數(shù)值通量的龍格-庫塔間斷有限元方法(Runge-Kutta discontinuous Galerkin method, 簡(jiǎn)稱RKDGM)及基于內(nèi)部罰函數(shù)的間斷有限元方法(Interior penalty discontinuous Galerkin method，簡(jiǎn)稱IP DGM)為代表的一系列DGM方法被提出，在計(jì)算流體力學(xué)、計(jì)算熱力學(xué)、計(jì)算電磁學(xué)、計(jì)算地球物理等領(lǐng)域得到了廣泛應(yīng)用(e.g., Li, 2006; Hesthaven and Warburton, 2008; Rivière, 2008; 張鐵, 2015; 孟雄等, 2015).間斷有限元方法是傳統(tǒng)有限元方法的一種推廣，它最主要的特征是未知量以及基函數(shù)是在每個(gè)網(wǎng)格單元上獨(dú)立定義的，在其它網(wǎng)格單元上取值為0，這種定義使得它具有許多良好的性質(zhì).相對(duì)于FDM，它的網(wǎng)格剖分靈活使得它可以處理任意復(fù)雜邊界；相對(duì)于FEM，它可以使用任意階基函數(shù)從而可以達(dá)到任意階精度，避免了FEM在高階格式上構(gòu)造的困難，且具有良好的并行性；相對(duì)于SEM，它可以采用更靈活的網(wǎng)格剖分——非相容(non-conforming)網(wǎng)格.特別地，DGM允許變量在單元之間存在間斷，因而特別適合模擬強(qiáng)地面運(yùn)動(dòng)問題，它不僅能適應(yīng)非平面斷層及速度反差較大的介質(zhì)，而且能有效避免滑移率時(shí)間序列中存在的虛假高頻振蕩的影響(De La Puente et al., 2009; Pelties et al., 2012).另外，DGM引入的數(shù)值頻散較小，且可以采用非均勻網(wǎng)格，因而特別適用于大尺度非均勻介質(zhì)中的波場(chǎng)模擬.DGM的質(zhì)量和剛度矩陣可以提前計(jì)算并存儲(chǔ)，因此在實(shí)際計(jì)算過程中不需要計(jì)算質(zhì)量矩陣和剛度矩陣，這種無積分(quadrature-free)的技巧使得計(jì)算量大為減少(Atkins and Shu, 1998).然而，盡管DGM可以采用無積分技巧，但是其計(jì)算量和存儲(chǔ)量相對(duì)于FEM及SEM還是大很多，這是由于DGM引入了較多的自由度從而導(dǎo)致存儲(chǔ)量增大，而且相比于FEM及SEM來說，DGM需要更小的時(shí)間步長(zhǎng)才能保持格式穩(wěn)定(de Basabe et al., 2008)，進(jìn)一步導(dǎo)致了計(jì)算量的增大.

DGMs在計(jì)算地球物理領(lǐng)域得到了廣泛應(yīng)用.其中應(yīng)用最為廣泛的是由K?ser and Dumbser(2006)提出的任意階導(dǎo)數(shù)的時(shí)間推進(jìn)步間斷有限元方法(arbitrary high-order derivatives time stepping discontinuous Galerkin method, 簡(jiǎn)稱ADER-DGM),這種方法基于迎風(fēng)數(shù)值通量，采用具有任意階精度的Lax-Wendroff時(shí)間推進(jìn)方式，從而使得時(shí)間精度和空間精度均能達(dá)到任意階精度.ADER-DGM及其衍生方法已被成功用于數(shù)值求解3D彈性、粘彈性、孔隙介質(zhì)、流固耦合等模型的波傳播問題及地震斷層破裂數(shù)值模擬中(Dumbser and K?ser, 2006; Dumbser et al., 2007; K?ser and Dumbser, 2008; De La Puente et al., 2007, 2008).此外，也涌現(xiàn)了許多基于其它形式的3D DGMs用于求解非均勻介質(zhì)、粘彈性介質(zhì)、聲波-彈性波界面等復(fù)雜介質(zhì)的波傳播模擬及偏移成像問題(Wilcox et al., 2010; L?hivaara and Huttunen, 2010; Etienne et al., 2010; Pelties et al., 2012; Minisini et al., 2013; Mu et al., 2013a, 2013b; Mulder et al., 2014; Ferroni et al., 2017; Ye et al., 2016; Lambrecht et al., 2017; Geevers et al., 2018).由于基于通量函數(shù)的DGM主要針對(duì)于求解雙曲方程，因此在求解地震波方程領(lǐng)域使用更為廣泛，本研究也主要基于數(shù)值通量形式的DGM.另外，本文中提到的DGM方法均指不帶限制器的DGM，因?yàn)榈卣鸩▌?dòng)方程大部分是線性方程，Cockburn 和 Shu(2001)指出當(dāng)雙曲系統(tǒng)為線性系統(tǒng)時(shí)，可以不加限制器，但是，如果考慮的是非線性方程或者解存在間斷時(shí)，必須使用限制器或者特殊的數(shù)值通量才能保證數(shù)值格式的精度(Chabot et al., 2018)，這已超出本文的研究范圍，在此也不作討論.

在國(guó)內(nèi)研究方面，汪文帥等(2013)、廉西猛等(2013)、薛昭等(2014)、賀茜君等(2014)最早將DGM應(yīng)用到求解波動(dòng)方程，隨后He 等(2015)、Yang 等(2016)、Meng 等(2018)、張金波等(2018)、He 等(2019a, 2019b, 2020)、Zhang 等(2019)等將其應(yīng)用到2D各種復(fù)雜情形的模擬和分析中.對(duì)于3D情形，He 等(2020)已經(jīng)開始著手研究和分析工作，他們將一種加權(quán)的DGM推廣至3D各向異性介質(zhì)中彈性波的傳播，采用了MPI并行策略，但是使用的是3D規(guī)則的六面體單元.本研究針對(duì)3D非結(jié)構(gòu)網(wǎng)格，發(fā)展了求解3D D′Alembert介質(zhì)中的并行WRKDG方法.D′Alembert介質(zhì)是一種粘彈性介質(zhì)，它直接對(duì)運(yùn)動(dòng)方程加入粘性項(xiàng)來刻畫粘性，因而具有簡(jiǎn)潔的表達(dá)式.Li 等(2015)、Cai等(2017)、L?hivaara 和 Huttunen(2010)等人對(duì)這種介質(zhì)進(jìn)行了數(shù)值模擬研究，本文的研究也是基于此.另外，我們對(duì)并行算法的設(shè)計(jì)及編程也給出了較為詳細(xì)的分析.特別地，我們根據(jù)常微分方程理論推導(dǎo)了3D D′Alembert介質(zhì)中滿足數(shù)值穩(wěn)定性條件的一般經(jīng)驗(yàn)公式.同時(shí)，由于針對(duì)基于數(shù)值通量的DGM，目前還沒有相關(guān)的數(shù)值穩(wěn)定性和數(shù)值頻散、耗散的分析結(jié)果，因此本文首次對(duì)該方法的數(shù)值頻散和耗散進(jìn)行了詳細(xì)的研究，且考慮了耗散參數(shù)對(duì)結(jié)果的影響.最后，我們給出了包含均勻模型、非規(guī)則幾何模型以及非均勻Marmousi模型在內(nèi)的數(shù)值模擬算例以說明方法的有效性.

1 WRKDG方法

1.1 空間離散

我們考慮3D D′ Alembert介質(zhì)中聲波方程(牛濱華和孫春巖, 2007; L?hivaara and Huttunen, 2010)：

(1)

其中u是位移場(chǎng)，c是聲波速度，f(t)是源項(xiàng)，r>0是D′ Alembert介質(zhì)中引入的耗散參數(shù).參數(shù)r與頻率ω有如下關(guān)系式(牛濱華和孫春巖, 2007)：

其中，Q為D′Alembert 介質(zhì)縱波的品質(zhì)因子.當(dāng)r/ω≤1時(shí)，我們有r≈ωQ-1.設(shè)3D求解區(qū)域?yàn)棣?我們將區(qū)域Ω劃分為不重疊的子區(qū)域{Ωi}.由于本文主要研究非結(jié)構(gòu)網(wǎng)格，因此主要采取四面體網(wǎng)格剖分.需要指出的是，為方便起見，本文僅考慮相容網(wǎng)格的情況，也就是不允許有“懸點(diǎn)”的存在.

DGM基于Galerkin近似，所以我們首先需選取試驗(yàn)函數(shù)空間.我們使用的試驗(yàn)函數(shù)空間是多項(xiàng)式空間Vh={v∈L1(Ω):v|Ωi∈Pκ(Ωi)}，其中Pκ(Ωi)表示區(qū)域Ωi上次數(shù)不超過κ次的多項(xiàng)式.從定義中可以看出，測(cè)試函數(shù)v在Ωi以外的區(qū)域上不定義或者令其為0，因此，它不必在整個(gè)區(qū)域上連續(xù)，可以在子區(qū)域{Ωi}之間有間斷，這也是間斷有限元與傳統(tǒng)有限元最大的區(qū)別.

為了將方程(1)改寫成一階雙曲系統(tǒng)的形式，我們引入三個(gè)變量p、q和s，且p、q、s滿足條件pt=ux,qt=uy,st=uz，其中ux、uy和uz分別表示位移u對(duì)空間變量x、y和z的偏導(dǎo)數(shù).則方程(1)經(jīng)過一次時(shí)間積分后，可改寫成如下形式：

(2)

即：

(3)

(4)

則可以將3D D′ Alembert介質(zhì)中聲波方程化成適于用DGM求解的如下形式：

(5)

其中F(u)表示通量，B表示耗散項(xiàng)，當(dāng)B=0時(shí)表示無耗散，方程(5)退化為普通的聲波方程.需要指出的是，本文發(fā)展的方法也適用于一階速度-應(yīng)力方程.在(5)式兩邊同乘以試驗(yàn)函數(shù)v，并在子區(qū)域Ωi上積分，利用格林公式，我們得到

+??ΩivF(u)·ndS=?ΩivfdV,

(6)

(7)

(8)

在上式中，由于我們考慮的是線性問題，不妨將通量F(u)寫成F(u)=(A1u,A2u,A3u)的形式(Dumbser and K?ser, 2006)，其中

(9)

(10)

其中，

(12)

I表示單位矩陣.令u具有如下的基函數(shù)展開形式：

(12)

將(12)和(10)式代入(8)式中，我們可以得到如下半離散化的方程：

l′=1,…,N.

(13)

本研究中采用的是無積分(quadrature-free)的DGM(Atkins and Shu, 1998)，因此只需提前計(jì)算好參考單元上的所有質(zhì)量矩陣和剛度矩陣.我們選取如圖1中所示的參考單元(Dumbser and K?ser, 2006)，且假設(shè)參考單元內(nèi)的坐標(biāo)軸三分量為：ξ，η和ζ.對(duì)于任意四面體單元，均將之通過坐標(biāo)變換到如圖1所示的參考單元，圖1a中的頂點(diǎn)與圖1b中的頂點(diǎn)嚴(yán)格對(duì)應(yīng)，具體的變換過程可參考附錄A.我們仿照Dumbser 和 K?ser(2006)的流程計(jì)算所有的體質(zhì)量矩陣和剛度矩陣，以及面與面之間關(guān)聯(lián)的質(zhì)量矩陣.具體的過程可參考附錄A.

圖1 一般四面體單元變換到參考單元示意圖(Dumbser and K?ser, 2006)，其中參考單元內(nèi)1、2、3和4這四個(gè)頂點(diǎn)的坐標(biāo)分別是(0, 0, 0)、(1, 0, 0)、(0, 1, 0)和(0, 0, 1)Fig.1 Transformation from the general tetrahedron to the reference tetrahedron(Dumbser and K?ser, 2006)with four vertices (0, 0, 0), (1, 0, 0), (0, 1, 0), and (0, 0, 1)

1.2 時(shí)間離散

(14)

(15)

(16)

(17)

圖2 Von Neumann分析中用到的網(wǎng)格構(gòu)型.在(a)中所示的規(guī)則六面體剖分的基礎(chǔ)上，每個(gè)六面體中有圖(b)中所示六個(gè)四面體網(wǎng)格(Mulder et al., 2014)Fig.2 Grid configuration used in Von Neumann analysis. Based on the (a) regular hexahedral division, each hexahedron has (b) six tetrahedrons (Mulder et al., 2014)

2 Von Neumann分析

3D情形下的數(shù)值穩(wěn)定性、數(shù)值頻散和耗散主要基于Von Neumann分析.為此，假設(shè)一個(gè)簡(jiǎn)諧波傳播經(jīng)過3D區(qū)域，并假定此區(qū)域是均勻、無界區(qū)域，且方程右端為無源項(xiàng).我們主要考察四面體網(wǎng)格剖分，在如圖2a中所示的規(guī)則六面體剖分的基礎(chǔ)上(Mulder et al., 2014; Ferroni et al., 2017; He et al., 2020)，在一個(gè)立方體單元內(nèi)有6個(gè)四面體單元，如圖2b所示，其中每個(gè)四面體及其面的順序參考Mulder 等(2014)編號(hào).我們假設(shè)如下平面簡(jiǎn)諧波在所考慮的區(qū)域內(nèi)傳播：

Cn,m,l(t)=C(t)exp[i(ksinθcosφnh+ksinθsinφmh

+kcosθlh-iωt)],

(18)

其中k是波數(shù)，ω是頻率，h是立方體的邊長(zhǎng)，θ和φ決定了平面波的傳播方向，θ∈[0,π]表示波傳播方向與z軸的夾角，而φ∈[0,2π]則表示波傳播方向在xy平面內(nèi)的投影與x軸的夾角.將(18)式代入數(shù)值格式中，即可用于分析數(shù)值穩(wěn)定性和數(shù)值頻散及耗散.

2.1 數(shù)值穩(wěn)定性分析

為了保持?jǐn)?shù)值算法的穩(wěn)定性，我們引入庫朗數(shù)條件，也稱Courant-Friedrichs-Lewy(簡(jiǎn)稱CFL)條件：α=cΔt/h≤αmax，其中α是庫朗數(shù)，αmax即是最大庫朗數(shù)(也稱CFL條件數(shù)).不妨令Λ表示代入簡(jiǎn)諧波(18)之后的數(shù)值格式增長(zhǎng)矩陣的特征值，則Λ與波數(shù)k、庫朗數(shù)α、傳播方向θ和φ有關(guān).要使得數(shù)值格式穩(wěn)定，必須滿足|Λ|≤1對(duì)所有的kh∈[0,2π]、θ∈[0,π]和φ∈[0,2π]都成立.這在數(shù)學(xué)上可以用下面一個(gè)優(yōu)化問題來表示：

maxα

s.t.|Λ(kh,α,θ,φ)|≤1 for ?kh∈[0,2π],

θ∈[0,π] andφ∈[0,2π],

α≥0.

(19)

上述求α的最大值問題是一個(gè)非線性優(yōu)化問題，通常情況下不容易求出αmax的解析表達(dá)式.一般我們通過數(shù)值遍歷算法來求其近似解，使α從0開始以小步長(zhǎng)增長(zhǎng)，直至|Λ(kh,α,θ,φ)|>1則跳出循環(huán)，獲得αmax的值.

在本文中，我們僅考慮基函數(shù)為1、2次多項(xiàng)式的情形，對(duì)應(yīng)空間精度分別為2、3階.首先，我們考慮不帶耗散的情形，即在方程(1)中r= 0的情形.通過數(shù)值求解上述非線性優(yōu)化問題，我們得到當(dāng)η=0.5時(shí)的WRKDG方法的最大庫朗數(shù)αmax：對(duì)于P1階WRKDG方法，cΔt/h≤αmax≈0.244；對(duì)于P2階WRKDG方法，cΔt/h≤αmax≈0.144.

以上分析給出的最大庫朗數(shù)αmax中所采用的網(wǎng)格步長(zhǎng)h是圖2a中立方體邊長(zhǎng)，若是采用圖2b中四面體內(nèi)切球體直徑d，則此時(shí)的時(shí)間步長(zhǎng)所需滿足的條件是

(20)

其中因子0.3597是對(duì)應(yīng)圖2中，當(dāng)立方體邊長(zhǎng)為h時(shí)，此時(shí)最小的四面體內(nèi)切球體直徑約為d≈0.3597h.

對(duì)于帶耗散情形，即方程(1)中r>0的情形，也可以根據(jù)(19)式中同樣的流程進(jìn)行分析，求出3D D′Alembert介質(zhì)中的庫朗數(shù)條件.然而此時(shí)由于加入了耗散參數(shù)r，我們希望能得到一個(gè)包含r的顯式的αmax的表達(dá)式.為此，在已知聲波方程庫朗數(shù)條件的基礎(chǔ)上，我們進(jìn)行如下分析.

我們首先注意到，經(jīng)過空間離散之后形成的半離散系統(tǒng)(13)可以分解為兩部分，一部分是無耗散的雙曲型系統(tǒng)，另一部分則與耗散有關(guān).因此，相應(yīng)地，我們可以將常微分方程組(14)分成兩部分，忽略震源項(xiàng)后，寫成如下形式(Carcione and Quiroga-Goode，1995)：

(21)

其中，算子L1對(duì)應(yīng)無耗散情形的聲波傳播，算子L2對(duì)應(yīng)系統(tǒng)中的耗散項(xiàng)，實(shí)際上，L2只與耗散參數(shù)r有關(guān).我們將算子L、L1和L2對(duì)應(yīng)的譜半徑分別記作λ、λ1和λ2.則由(21)，我們可得λ≤λ1+λ2.對(duì)于雙曲型系統(tǒng)：

(22)

我們?cè)谏衔闹幸呀o出關(guān)于Δt≤αmaxh/c的結(jié)果.利用求解常微分方程的數(shù)值分析理論(李慶楊等, 2008)，我們知道對(duì)于加權(quán)的Runge-Kutta時(shí)間離散格式(15)，當(dāng)η=0.5時(shí)，要使得格式穩(wěn)定，必須滿足如下條件：

(23)

從而可得

(24)

對(duì)于帶耗散的系統(tǒng)：

(25)

由于算子L2只與矩陣耗散參數(shù)r有關(guān)，所以算符的譜半徑等于r，即

λ2≈r.

(26)

由于λ≤λ1+λ2，所以由方程(24)及(26)可得

(27)

因此，我們得到利用WRKDG方法求解3D D′ Alembert介質(zhì)中的數(shù)值穩(wěn)定性條件的經(jīng)驗(yàn)公式：

(28)

方程(28)中給出的時(shí)間步長(zhǎng)限制是保持計(jì)算穩(wěn)定的充分條件，但不是必要條件.我們?cè)诒?列出了當(dāng)波速c=1 km·s-1, 隨hr變化時(shí)，P1階和P2階WRKDG方法的真實(shí)最大庫朗數(shù)(αmax)D′Ale與從經(jīng)驗(yàn)公式(28)獲得的值(αmax)′D′Ale進(jìn)行對(duì)比，從表中可以看出兩者之間相差不大.因此，在實(shí)際應(yīng)用中，(28)式是一個(gè)合理的估計(jì).若是采用四面體內(nèi)切球體直徑d，則(28)式可寫為

(29)

表1 3D D′ Alembert介質(zhì)中P1階和P2階WRKDG方法的真實(shí)最大庫朗數(shù)(αmax)D′Ale與從經(jīng)驗(yàn)公式(28)獲得的值(αmax)′D′Ale的對(duì)比

另外，從方程(28)我們可以看出，隨著耗散參數(shù)r的增大，時(shí)間步長(zhǎng)減??；尤其當(dāng)r的值很大時(shí)，此時(shí)系統(tǒng)(21)是一個(gè)剛性(stiff)系統(tǒng)，一般的時(shí)間推進(jìn)方法將導(dǎo)致極小的時(shí)間步長(zhǎng)，因而不再適用，此時(shí)應(yīng)尋找特殊的時(shí)間推進(jìn)方法.

2.2 數(shù)值頻散、耗散分析

在本節(jié)中，我們將討論求解D′Alembert介質(zhì)中聲波方程的3D WRKDG方法的數(shù)值頻散和數(shù)值耗散情況.對(duì)于方程(18)中的簡(jiǎn)諧波表達(dá)式，當(dāng)代入數(shù)值格式中，如果我們假定波數(shù)k是實(shí)數(shù)，則一般說來角頻率ω是復(fù)數(shù)，不妨假設(shè)ω=ωr-iωi，其中實(shí)部ωr表示ω中與頻散有關(guān)的量，非負(fù)的虛部ωi表示ω中與耗散有關(guān)的量.數(shù)值波速可由cnum=ωr/k估計(jì).我們引入采樣率δ=kh/(2π)，并將數(shù)值頻散R定義為數(shù)值速度與真實(shí)物理速度之比，則R的表達(dá)式為

(30)

其中α=cΔt/h是庫朗數(shù).我們將振幅在一個(gè)時(shí)間步內(nèi)的衰減記作S=e-ωiΔt，S可以反映D′ Alembert介質(zhì)中的衰減情況.

圖3給出了用WRKDG方法計(jì)算的數(shù)值頻散情況.我們?cè)O(shè)置參數(shù)c=2 km·s-1，h=0.05 km，α=cΔt/h=0.1.在3D情形，數(shù)值頻散R的大小與傳播方向θ和φ有關(guān)，在圖3中，我們僅給出了θ=π/2、φ=0,π/12,π/4情況下，WRKDG方法的數(shù)值頻散隨采樣率δ的變化圖.注意此時(shí)由于h固定，δ的大小與波數(shù)成正比.圖3中(a—b)、(c—d)分別表示P1階、P2階WRKDG方法，(a)、(c)表示無耗散r=0的情形，而(b)、(d)表示耗散參數(shù)r=10的情形.從圖3中我們可以明顯觀察到數(shù)值頻散的各向異性特征.引起數(shù)值頻散各向異性的因素較多，數(shù)值算法本身、算法精度、網(wǎng)格剖分方式、網(wǎng)格步長(zhǎng)大小等都會(huì)影響數(shù)值頻散各向異性的產(chǎn)生及幅度.例如，在進(jìn)行數(shù)值頻散分析時(shí)，引入了方位角θ和φ，當(dāng)空間網(wǎng)格分布存在不對(duì)稱性時(shí)，會(huì)導(dǎo)致數(shù)值頻散出現(xiàn)各向異性特征.其次，不同算法精度產(chǎn)生的數(shù)值頻散各向異性也不相同.一般來說，減小網(wǎng)格的各向異性、提高算法精度、減小網(wǎng)格步長(zhǎng)可以有效降低數(shù)值頻散的各向異性.

圖3 在θ=π/2時(shí)，P1階和P2階WRKDG方法取φ=0,π/12,π/4時(shí)數(shù)值頻散R隨采樣率δ的變化情況.(a—b)P1階WRKDG方法，(c—d)P2階WRKDG方法，(a)、(c)對(duì)應(yīng)耗散參數(shù)r=0，而(b)、(d)對(duì)應(yīng)耗散參數(shù)r=10Fig.3 Numerical dispersion R as a function of sampling rate δ for P1 and P2 WRKDG methods when θ=π/2 and φ=0,π/12,π/4. (a—b) are for P1 WRKDG method, while (c—d) are for P2 WRKDG method. Panels (a) and (c) show the cases for r=0; panels (b) and (d) show the cases for r=10

另外，我們從圖3b與3d中發(fā)現(xiàn)，在δ較小時(shí)D′ Alembert介質(zhì)存在比較大的頻散，隨著δ的增加，R的值接近1，也就是說，隨著δ的增加，頻散變小.對(duì)比r=0與r=10這兩種情形，我們可以得出下述結(jié)論：在δ(波數(shù))較小時(shí)，D′ Alembert介質(zhì)中存在主要由耗散參數(shù)r引起的頻散；隨著δ(波數(shù))的增加，耗散參數(shù)引起的頻散變小，而數(shù)值算法引起的數(shù)值頻散占主導(dǎo).

圖4顯示了振幅在一個(gè)時(shí)間步Δt內(nèi)的耗散S=e-ωiΔt隨采樣率的變化圖.圖中參數(shù)的選取和上文中數(shù)值頻散分析中相同.當(dāng)r=0時(shí)(圖4a與圖4c)，S表示聲波介質(zhì)中的衰減情況，此時(shí)S完全由數(shù)值離散方法引起，而當(dāng)r=10時(shí)(圖4b與圖4d)，S的大小由數(shù)值離散方法引起的數(shù)值耗散與D′ Alembert介質(zhì)中的耗散參數(shù)r共同決定.圖4也體現(xiàn)了數(shù)值耗散的各向異性特征.在圖4中，對(duì)比耗散因子r=0與r=10這兩種情形，我們可以觀察到r=10時(shí)的耗散曲線有一個(gè)整體的下降，其值約為0.9876，約等于D′ Alembert介質(zhì)中理論耗散因子e-rΔt/2，體現(xiàn)了D′Alembert介質(zhì)的衰減特性.

圖4 在θ=π/2時(shí)，P1階和P2階WRKDG方法取φ=0,π/12,π/4時(shí)數(shù)值耗散S隨采樣率δ的變化情況.(a—b)P1階WRKDG方法，(c—d)P2階WRKDG方法，(a)、(c)對(duì)應(yīng)耗散參數(shù)r=0，而(b)、(d)對(duì)應(yīng)耗散參數(shù)r=10Fig.4 Numerical dissipation S as a function of sampling rate δ for P1 and P2 WRKDG methods when θ=π/2 and φ=0,π/12,π/4. (a—b) are for P1 WRKDG method, while (c—d) are for P2 WRKDG method. Panels (a) and (c) show the cases for r=0; panels (b) and (d) show the cases for r=10

3 并行計(jì)算

3.1 網(wǎng)格準(zhǔn)備

結(jié)構(gòu)網(wǎng)格(六面體)和非結(jié)構(gòu)網(wǎng)格(四面體)均可用于實(shí)際計(jì)算，一般說來，結(jié)構(gòu)網(wǎng)格精度較高，但是對(duì)于復(fù)雜幾何模型來說，生成結(jié)構(gòu)網(wǎng)格會(huì)花費(fèi)較大代價(jià)；非結(jié)構(gòu)網(wǎng)格計(jì)算精度不如結(jié)構(gòu)網(wǎng)格，但是其網(wǎng)格生成較簡(jiǎn)單.本研究中選取非結(jié)構(gòu)網(wǎng)格——四面體網(wǎng)格.網(wǎng)格剖分可以使用開源軟件或商業(yè)軟件.當(dāng)模型比較簡(jiǎn)單或網(wǎng)格量比較小的時(shí)候選擇開源軟件，當(dāng)模型比較復(fù)雜或者問題規(guī)模較大時(shí)，優(yōu)先選擇商業(yè)軟件.我們使用開源軟件或者商業(yè)軟件生成了四面體網(wǎng)格，獲得了所有網(wǎng)格節(jié)點(diǎn)的信息以及單元連接關(guān)系矩陣，為了計(jì)算的便利性，我們還需要計(jì)算面與面之間的關(guān)聯(lián)矩陣、面與單元之間的關(guān)聯(lián)矩陣、每個(gè)面中第一個(gè)頂點(diǎn)對(duì)應(yīng)的于相鄰面中的局部編號(hào)矩陣.這些矩陣的計(jì)算過程可參考Hesthaven 和 Warburton(2008)的研究.

3.2 網(wǎng)格分劃

由于3D情形計(jì)算量和存儲(chǔ)量均較大，所以要進(jìn)行并行處理.在并行之前，我們需要對(duì)整體的3D網(wǎng)格進(jìn)行分劃，將其分給不同的進(jìn)程.我們使用開源網(wǎng)格分劃程序包Metis(Karypis and Kumar, 1998).在使用時(shí)只需要輸入單元數(shù)、頂點(diǎn)數(shù)、單元連接矩陣、進(jìn)程數(shù)等數(shù)據(jù)，即可用一行命令對(duì)網(wǎng)格進(jìn)行快速分劃.例如，圖5顯示了Metis分劃的結(jié)果，圖5a代表將3D區(qū)域[0,2]×[0,2]×[0,2]離散成6000個(gè)四面體網(wǎng)格，用Metis劃分給5個(gè)不同的進(jìn)程，圖5b給出了運(yùn)行結(jié)果，圖5b中不同顏色代表分給不同進(jìn)程.

圖5 (a) 6000個(gè)四面體網(wǎng)格；(b)利用Metis將(a)中所有網(wǎng)格分給5個(gè)進(jìn)程Fig.5 (a) 6000 tetrahedrons; (b) Metis divides all tetrahedrons in (a) into 5 processors

3.3 并行流程

在本研究中，我們使用Message Passing Interface(MPI)編程策略對(duì)程序進(jìn)行并行化處理，整個(gè)流程可參考圖6.在程序開始階段，由于我們采用無積分(quadrature-free)的DGM，因此可將方程(13)中所需的參考單元內(nèi)的質(zhì)量、剛度矩陣以及四面體面積分矩陣提前計(jì)算出來，然后導(dǎo)入到程序中.同時(shí)，我們也將生成的網(wǎng)格文件導(dǎo)入，并利用Metis進(jìn)行網(wǎng)格劃分，同時(shí)需要計(jì)算出網(wǎng)格單元連接矩陣，以備后用.在用Metis進(jìn)行網(wǎng)格劃分完畢后，對(duì)于每一個(gè)進(jìn)程我們需要記錄三類網(wǎng)格及其信息.以圖7來進(jìn)行說明，第一類網(wǎng)格是每個(gè)進(jìn)程的內(nèi)網(wǎng)格，例如，以第一個(gè)進(jìn)程為例，圖7中藍(lán)色區(qū)域內(nèi)的網(wǎng)格即是第一個(gè)進(jìn)程的內(nèi)網(wǎng)格；第二類網(wǎng)格是該進(jìn)程的輔助網(wǎng)格，這類網(wǎng)格屬于其他進(jìn)程的內(nèi)網(wǎng)格，位于進(jìn)程1所有網(wǎng)格的邊界處，如圖7a中的綠色網(wǎng)格部分(圖7b是圖7a的一個(gè)更清晰的展示)，作用是接收來自于其他進(jìn)程更新后的變量信息；第三類網(wǎng)格是屬于該進(jìn)程內(nèi)，但是作為其他進(jìn)程的輔助網(wǎng)格，這類網(wǎng)格在消息傳遞時(shí)需要由本進(jìn)程傳遞給其他進(jìn)程使用，如圖7c中紅色網(wǎng)格顯示.

圖6 并行算法流程圖Fig.6 Flow chart of the parallel algorithm

開始階段結(jié)束后，如流程圖6中的說明，我們將變量賦初值，然后進(jìn)入時(shí)間迭代.在每個(gè)時(shí)間迭代步中，首先更新每個(gè)進(jìn)程內(nèi)網(wǎng)格(即第一類網(wǎng)格)的變量信息.這一步結(jié)束后，我們需要進(jìn)行兩步消息傳遞：將所在進(jìn)程內(nèi)屬于其他進(jìn)程的輔助網(wǎng)格(即第三類網(wǎng)格)上的Cn+1發(fā)送給相應(yīng)進(jìn)程，并接收來自其他進(jìn)程的Cn+1并存放于輔助網(wǎng)格(即第二類網(wǎng)格)中.這樣我們就完成了一步完整的時(shí)間迭代.時(shí)間迭代結(jié)束后，輸出數(shù)據(jù)結(jié)果，并利用畫圖軟件等進(jìn)行畫圖.

圖7 (a) 所考慮進(jìn)程的輔助網(wǎng)格示意圖，即圖中綠色網(wǎng)格部分，這類網(wǎng)格屬于其他進(jìn)程的內(nèi)網(wǎng)格，位于該進(jìn)程所有網(wǎng)格的邊界處；(b) 圖a的側(cè)面圖；(c) 屬于該進(jìn)程內(nèi)、作為其他進(jìn)程的輔助網(wǎng)格，即圖中紅色網(wǎng)格部分Fig.7 (a) Illustration of the auxiliary grids of this processor (the green part in the figure). This type of grids belongs to the internal grids of other processors, and is located at the boundary of this process; (b) the side view of figure a; (c) the auxiliary grids of other processors which belongs to this processor (the red part in the figure)

4 精度測(cè)試及并行效率分析

我們考慮一個(gè)帶有解析解的模型來驗(yàn)證該方法的正確性和收斂性.對(duì)于無源的方程(1)，考慮如下解析解(Cai et al., 2017)

u(t,x,y,z)=e-rt/2cos(K1x+K2y+K3z-Wt),

(31)

u(t,x,y,z)=e-rt/2cos(K1x+K2y+K3z-Wt),

p(t,x,y,z)=

q(t,x,y,z)=

s(t,x,y,z)=

(32)

在本例中，選取計(jì)算區(qū)域?yàn)?≤x,y,z≤2 km，速度c=2 km·s-1，K1=K2=K3=π,，以t=0 s時(shí)刻的值作為初始條件，采用周期邊界條件.由于我們主要考察空間離散的誤差和收斂精度，因此我們?cè)O(shè)置時(shí)間步長(zhǎng)為0.1 ms，此時(shí)，由時(shí)間離散引起的誤差可以忽略，數(shù)值誤差主要來自于空間離散.我們采用如圖2所示的網(wǎng)格剖分及四面體單元，圖中每個(gè)小立方體的邊長(zhǎng)為h，且每個(gè)立方體含有6個(gè)四面體單元.我們令N表示在x，y及z三個(gè)方向劃分的立方體個(gè)數(shù)，則四面體總數(shù)目為6N3.定義L2與L1誤差為

(33)

其中uex是方程(31)給出的精確解.我們令N=4,8,10,16,20，在表2中列出了T=0.1 s時(shí)刻在耗散參數(shù)r=1以及r=10兩種情形下的數(shù)值誤差和相應(yīng)的收斂階.從表2可以看出，數(shù)值誤差隨著空間步長(zhǎng)的減小而減小，說明WRKDG方法是收斂的.由于我們使用了二次基函數(shù)，因此得到了預(yù)期的三階空間精度.同時(shí)，從表2中可以看出，隨著耗散參數(shù)r的增大，誤差也隨之減小，體現(xiàn)了D′Alembert介質(zhì)的衰減特性.

表2 3D D′ Alembert介質(zhì)中P2階WRKDG的誤差及收斂階Table 2 Convergence rates of u for the acoustic wave in 3D D′Alembert medium

圖8 3D WRKDG算法的并行加速比(圖中線型“-o”表示).其中橫軸表示進(jìn)程數(shù)，縱軸表示并行加速比.圖中線型“-*”表示理想情形下的并行加速比Fig.8 Parallel speed-ups of the 3D WRKDG algorithm (see the line type “-o”). The horizontal axis is the number of processors, and the vertical axis is the parallel speed-ups. The line type “-*” in the figure represents the parallel speed-up in the ideal case

接下來，為了考察3D WRKDG算法的并行效率，我們將上述精度測(cè)試中N=20的數(shù)值模擬程序進(jìn)行并行化處理，此時(shí)計(jì)算區(qū)域被離散成48000個(gè)四面體，利用Metis將網(wǎng)格分別分劃給2、4、8、16、32、64個(gè)進(jìn)程，記錄運(yùn)行的CPU時(shí)間數(shù).假設(shè)每個(gè)處理器的計(jì)算性能相同，在此條件下，并行程序的加速比(speed-up)可定義為：SP=TS/TP,其中TS是串行程序運(yùn)行的時(shí)間，TP是并行程序在P個(gè)進(jìn)程下運(yùn)行的時(shí)間，SP在通常情況下總是小于P，因?yàn)樵诓⑿谐绦蛟O(shè)計(jì)時(shí)往往會(huì)引入一些通信時(shí)間及其他管理花費(fèi).圖8給出了3D WRKDG方法的并行程序加速比及理想情形下的加速比，從圖中可以看出，當(dāng)進(jìn)程數(shù)比較小的時(shí)候，加速比接近理論情形，但是隨著進(jìn)程數(shù)越增大，實(shí)際加速比越偏離理論情形，這是由于進(jìn)程數(shù)增加引起進(jìn)程與進(jìn)程之間的通信時(shí)間開銷大大增大，同時(shí)進(jìn)程與進(jìn)程之間等待的時(shí)間開銷也增大，從而降低了并行效率.本文所使用的計(jì)算平臺(tái)是國(guó)家超級(jí)計(jì)算天津中心的天河TH-1A高性能機(jī)群系統(tǒng)，每個(gè)節(jié)點(diǎn)上有12個(gè)主頻為 2.93 GHz 的核，每個(gè)節(jié)點(diǎn)內(nèi)存為 24 GB.

5 波場(chǎng)模擬

在本節(jié)中，我們通過幾個(gè)數(shù)值例子來說明本文所給出的WRKDG方法在求解3D D′Alembert介質(zhì)中聲波傳播問題的有效性.考慮到更高階格式的庫朗數(shù)條件更為嚴(yán)格，且存儲(chǔ)量和計(jì)算量也越大，因此，如果沒有特殊說明，在本節(jié)中我們僅使用空間精度為3階的P2階WRKDG方法進(jìn)行波場(chǎng)模擬.

5.1 均勻介質(zhì)模型

在這個(gè)例子中，我們考查D′Alembert介質(zhì)中聲波在3D均勻介質(zhì)中的傳播.選取計(jì)算區(qū)域?yàn)?≤x,y,z≤2 km的立方體區(qū)域，我們將其離散為3072000個(gè)四面體，每個(gè)四面體邊長(zhǎng)平均約為25 m.聲波速度c=3 km·s-1, 時(shí)間步長(zhǎng)取Δt≈1.29 ms.震源函數(shù)是：

f(t)=-9.6f0(0.6f0t-1)exp[-8(0.6f0t-1)2].

(34)

其中f0=45 Hz，主頻約為20 Hz.震源位于(0.981251,1.00625,1.00625)km處，在坐標(biāo)(1.35,1.35,1.35)km處設(shè)置一個(gè)接收點(diǎn)用于記錄波形信息.我們首先考查無耗散情形，即r=0.圖9給出了T=0.5 s時(shí)刻接收點(diǎn)接收到的歸一化的波形記錄，圖中紅色實(shí)線是用Cagnidard-de-Hoop方法(Aki and Richards，2002)計(jì)算得到的解析解，而藍(lán)色虛線及黑色實(shí)線分別表示利用P2及P1方法得到的數(shù)值解.從圖中可以看出，P1方法出現(xiàn)少許數(shù)值頻散，而P2方法與解析解吻合較好，這說明了提高算法精度有助于降低數(shù)值頻散.圖10給出了P2方法的波場(chǎng)快照?qǐng)D，此時(shí)波已經(jīng)傳至邊界.波場(chǎng)快照?qǐng)D中無明顯可見數(shù)值頻散.

圖9 在耗散參數(shù)r=0的3D均勻介質(zhì)模型中，T=0.5 s 時(shí)刻接收點(diǎn)處的歸一化的波形記錄圖，圖中紅色實(shí)線代表解析解，藍(lán)色虛線及黑色實(shí)線分別表示利用P2及P1方法得到的數(shù)值解Fig.9 Comparisons of normalized waveforms at time T=0.5 s for the 3D homogeneous model with dissipation parameter r=0, in which the red solid line in the figure represents the analytical solution, and the blue dotted line and black solid line represent the numerical solution computed by the P2 and P1 methods, respectively

圖10 在耗散參數(shù)r=0的3D均勻介質(zhì)模型中，使用P2方法計(jì)算得到的T=0.5 s 時(shí)刻的波場(chǎng)快照?qǐng)DFig.10 Snapshots of the acoustic wave fields computed by the P2 method at time T=0.5 s for the 3D homogeneous model with dissipation parameter r=0

圖11 在均勻介質(zhì)模型中，不同耗散參數(shù)r=0、2、4、8和16對(duì)應(yīng)的接收器處的聲波波形記錄Fig.11 Waveform records at the receiver with different dissipation parameters r=0,2,4,8 and 16 for the homogeneous model

表3 3D D′ Alembert介質(zhì)中波形記錄的波谷處的振幅和衰減系數(shù)Table 3 The amplitudes and attenuation ratios at the trough at the receiver for the acoustic wave in D′Alembert medium

5.2 非規(guī)則幾何模型

在這個(gè)例子中，我們主要利用3D WRKDG方法模擬波在非規(guī)則幾何模型中的傳播，模型如圖12所示，在計(jì)算區(qū)域?yàn)?≤x,y,z≤2 km的立方體中，有一個(gè)球狀區(qū)域，球中心坐標(biāo)(1, 1, 0.5)km，半徑0.2 km.球內(nèi)介質(zhì)波速1.5 km·s-1，球外介質(zhì)波速3 km·s-1.我們采用2179529個(gè)四面體的非均勻網(wǎng)格離散，其中球外的四面體最大邊長(zhǎng)不超過0.05 km，在球面上四面體最大邊長(zhǎng)不超過0.02 km.圖13a給出了球體部分四面體網(wǎng)格的3D顯示，圖13b給出了清晰的2D剖面y=0處的網(wǎng)格剖分示意圖，從圖中可以看出，四面體網(wǎng)格可以貼合內(nèi)邊界——球面生成，且在包裹球體的一個(gè)立方體內(nèi)的網(wǎng)格密度大，而在遠(yuǎn)離此立方體的地方網(wǎng)格密度小.時(shí)間步長(zhǎng)取Δt≈0.52 ms.震源函數(shù)與方程(34)中相同，其中f0=60 Hz， 主頻約為25 Hz.圖14給出了在T=0.3 s時(shí)刻下的波場(chǎng)快照?qǐng)D，圖14a對(duì)應(yīng)于無耗散情形， 而圖14b對(duì)應(yīng)于耗散參數(shù)r=4.從圖中可以看出，圖14b較圖14a暗一些，證明了 D′Alembert介質(zhì)中的衰減效應(yīng).

圖12 非規(guī)則幾何模型示意圖，在計(jì)算區(qū)域0≤x,y,z≤2中，有一個(gè)球狀區(qū)域，球中心坐標(biāo)(1,1,0.5)km，半徑0.2 kmFig.12 Illustration of the irregular geometric model. In the computational domain 0≤x,y,z≤2 km, there is a spherical area with spherical center coordinates (1,1,0.5) km and a radius of 0.2 km

圖13 (a) 球體部分四面體網(wǎng)格的3D示意圖； (b) 二維剖面y=0處的網(wǎng)格剖分示意圖Fig.13 Illustration of (a) tetrahedrons in the ball and (b) the grid division at the cross section y=0

圖14 T=0.3 s 時(shí)刻的波場(chǎng)快照?qǐng)D，其中(a)對(duì)應(yīng)于無耗散情形r=0， 而(b)對(duì)應(yīng)于耗散參數(shù)r=4Fig.14 Snapshots of the seismic waves at T=0.3 s with dissipation parameters (a) r =0 and (b) r=4

5.3 Marmousi模型

在這個(gè)例子中，我們選取Marmousi速度模型(Versteeg and Grau, 1991)以測(cè)試WRKDG方法在非均勻復(fù)雜介質(zhì)情況下的計(jì)算效果.為了簡(jiǎn)化3D模型，我們采取2D Marmousi模型在z方向進(jìn)行平移得到3D模型，模型尺寸是9.216×2.928×2.928 km，其速度結(jié)構(gòu)如圖15所示.Marmousi模型有很強(qiáng)的非均勻性，其速度變化范圍是1.5～5.5 km·s-1.本實(shí)驗(yàn)采用2250000個(gè)四面體，四面體平均邊長(zhǎng)為24 m，震源函數(shù)如方程(34)所示，其中f0=30 Hz， 主頻約為13 Hz，震源位于(4.577,0.015,1.449)km處.模擬中時(shí)間步長(zhǎng)取Δt≈1.69 ms，D′Alembert介質(zhì)中耗散參數(shù)r=2.圖16給出了T=1.0 s時(shí)刻的波場(chǎng)快照，從圖中可以看出并無明顯可見的數(shù)值頻散.這說明3D WRKDG方法可以有效模擬復(fù)雜非均勻D′Alembert介質(zhì)中的聲波波場(chǎng).

圖15 3D Marmousi模型Fig.15 3D Marmousi model

圖16 對(duì)于3D Marmousi模型，T=1.0 s 時(shí)刻的波場(chǎng)快照?qǐng)D，其中耗散參數(shù)r=2Fig.16 Snapshots at T=1.0 s for the 3D Marmousi model with dissipation parameter r=2

6 結(jié)論

本文將加權(quán)Runge-Kutta間斷有限元(WRKDG)方法應(yīng)用于求解3D D′Alembert介質(zhì)中的聲波方程，空間離散采用了基于數(shù)值通量的間斷有限元公式，時(shí)間離散基于對(duì)角隱式Runge-Kutta方法，我們通過兩步迭代過程將其轉(zhuǎn)化為顯式方法，并在時(shí)間離散化過程中引入加權(quán)因子，最終獲得了求解3D D′Alembert介質(zhì)中聲波方程的WRKDG新方法.進(jìn)一步，我們?cè)敿?xì)研究了該方法的數(shù)值穩(wěn)定性條件，給出了四面體情形下的最大庫朗數(shù).由于D′Alembert介質(zhì)中耗散參數(shù)的引入，我們也推導(dǎo)了一種帶有耗散參數(shù)的數(shù)值穩(wěn)定性條件經(jīng)驗(yàn)公式.數(shù)值試驗(yàn)表明，該經(jīng)驗(yàn)公式是一種正確的估計(jì).此外，我們也深入研究了WRKDG方法在四面體情形下的數(shù)值頻散及數(shù)值耗散，研究表明D′ Alembert介質(zhì)中的數(shù)值頻散和耗散由耗散參數(shù)r及數(shù)值算法共同決定，存在一個(gè)理論耗散因子e-rΔt/2.同時(shí)，我們也觀察到數(shù)值頻散和數(shù)值耗散具有明顯的各向異性特征，這主要是由于所用網(wǎng)格的各向異性特征導(dǎo)致的.

我們通過數(shù)值模擬實(shí)驗(yàn)證明了WRKDG方法的收斂性，給出了3D并行WRKDG算法基于MPI并行策略下的加速比曲線，從中可以看出WRKDG算法具有良好的并行性能.為了進(jìn)一步驗(yàn)證3D WRKDG方法的正確性和有效性，我們模擬了聲波在D′Alembert介質(zhì)中均勻、非均勻介質(zhì)及非規(guī)則幾何模型中的傳播，且針對(duì)均勻介質(zhì)給出了理論耗散因子與觀測(cè)衰減因子，二者較為吻合.這些結(jié)果均表明3D WRKDG方法能夠正確且有效地模擬衰減介質(zhì)中的聲波傳播，能充分體現(xiàn)D′ Alembert介質(zhì)中波的衰減特征.

最后需要指出的是，盡管3D WRKDG方法應(yīng)用了并行策略能夠有效節(jié)省計(jì)算時(shí)間，但是其計(jì)算量和存儲(chǔ)量相對(duì)于其它數(shù)值方法還是很大，因此，今后我們應(yīng)重點(diǎn)考慮如何從多種途徑聯(lián)合提高其計(jì)算效率，以真正實(shí)現(xiàn)復(fù)雜介質(zhì)中大規(guī)模地震模擬、逆時(shí)偏移和基于波動(dòng)方程反演成像的快速計(jì)算和實(shí)際應(yīng)用，這些都是值得我們繼續(xù)深入研究的方向.

附錄A 方程(13)所需矩陣的具體表達(dá)式

根據(jù)方程(13)，我們需要計(jì)算如下四面體上的體積分矩陣：

(A1)

其中，角標(biāo)l和m表示矩陣的l行m列元素，上角標(biāo)i表示所考慮的單元為Ωi.在計(jì)算之前，我們首先將一般的四面體單元Ωi變換到如圖1所示的參考單元E內(nèi).如圖1所示，假設(shè)原單元四個(gè)頂點(diǎn)1、2、3及4的坐標(biāo)分別為 (x1,y1,z1)、(x2,y2,z2)、(x3,y3,z3)和(x4,y4,z4)，變換到參考單元E內(nèi)的四個(gè)頂點(diǎn)坐標(biāo)分別是(0,0,0)、(1,0,0)、(0,1,0)和(0,0,1).原坐標(biāo)三分量是x,y和z，且假設(shè)參考單元內(nèi)的坐標(biāo)軸三分量為：ξ，η和ζ，則任意四面體均將通過如下坐標(biāo)變換成如圖1所示的參考單元：

(A2)

易知：dxdydz=|J|dξdηdζ，其中|J|是四面體Ωi的體積，且容易得到如下偏導(dǎo)數(shù)的值(Dumbser and K?ser, 2006)：

(A3)

在(A2)所示的坐標(biāo)變換下，要計(jì)算方程(A1)中的矩陣，只需要在參考單元中計(jì)算如下矩陣即可：

(A4)

例如，要計(jì)算原質(zhì)量矩陣M1，可以利用關(guān)系式，M1=|J|M′1.

(A5)

(A6)

至此，我們已將方程(13)中所有積分矩陣的表達(dá)式闡述完畢.對(duì)于四面體的四個(gè)面F1、F2、F3和F4的外法向量n1,n2,n3和n4的計(jì)算，有如下公式，

n1=L13×L12,n2=L12×L14,

n3=L14×L13,n4=L23×L24，

(A7)

其中Lij表示以頂點(diǎn)i為起點(diǎn)，j為終點(diǎn)的向量.

表A1 四面體單元所屬四個(gè)面的定義順序(Dumbser and K?ser, 2006)Table A1 Face Definition on tetrahedrons (Dumbser and K?ser, 2006)

表A2 (a)三維坐標(biāo)軸ξ，η，和ζ與面積分用到的參變量和τ之間的對(duì)應(yīng)關(guān)系； (b)對(duì)于不同的h，在Ωi中的參變量和τ與相鄰單元Ωj中參變量′和τ′的對(duì)應(yīng)關(guān)系(Dumbser and K?ser, 2006)

Table A2 (a) Relationship between the three-dimensional coordinate axes ξ, η, and ζ and the face parameters and τ used in the area integrals; (b)Relationship between the face parameters and τ in the tetrahedron Ωi and the face parameters ′and τ′ in the adjacent tetrahedron Ωj (Dumbser and K?ser, 2006)

表A2 (a)三維坐標(biāo)軸ξ，η，和ζ與面積分用到的參變量和τ之間的對(duì)應(yīng)關(guān)系； (b)對(duì)于不同的h，在Ωi中的參變量和τ與相鄰單元Ωj中參變量′和τ′的對(duì)應(yīng)關(guān)系(Dumbser and K?ser, 2006)

Face1234ξτ01--τη0τζ0ττ(a)h123′τ1--ττ′τ1--τ(b)