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

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

0 引言

基于地震波動方程的正演是計算地球物理學的重要研究內容，它不僅能為我們提供精確的波場模擬結果，而且也是基于波動方程的地震勘探和地震學的重要基礎(牟永光和裴正林, 2005).在最近幾十年里，隨著計算機能力的快速提升，涌現(xiàn)了許許多多優(yōu)秀的數(shù)值算法，例如有限差分法(Finite difference method，簡稱FDM)、有限元法(Finite element method，簡稱FEM)、偽譜法(Pseudo-spectral method，簡稱PSM)、譜元法(Spectral element method，簡稱SEM)、間斷有限元法(Discontinuous Galerkin method，簡稱DGM)等.這些算法均能滿足我們對于數(shù)值算法的部分要求，它們也有其獨特的優(yōu)勢.對于數(shù)值求解3D波動方程來說，這些算法都得到了廣泛應用.FDM編程簡單、計算速度快，且存儲量小，是求解3D波動方程最為常用的一種數(shù)值算法(e.g.,董良國等, 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)和西方學者獨立提出，它能處理復雜區(qū)域和邊界條件，在求解3D波動方程也得到了一定的應用(Galis et al., 2008; Moczo et al., 2011; Igel, 2017)，但是由于其計算量大且并行性差，作為單一方法用于求解波傳播問題已不多見，它常與有限差分方法結合形成的混合方法(e.g., Ma et al., 2004; Galis et al., 2008).PSM利用快速傅里葉變換(Fast Fourier transform，簡稱FFT)來處理全部波場中的空間導數(shù), 精度很高，且數(shù)值頻散小，在求解波動方程領域也得到了廣泛應用(Furumura et al., 1998; Igel, 1999; Klin et al., 2010; Carcione et al., 2018)，但是3D情形的PSM需要利用全局數(shù)據(jù)進行FFT，因此并行性較差(Pelz, 1991).SEM是計算地球物理領域近十年來發(fā)展最為迅速的數(shù)值算法，它具有FEM的基本特征，在單元內部的基函數(shù)基于特定積分節(jié)點——Gauss-Lobatto-Legendre(GLL)點，因此在使用GLL積分準則的情況下質量矩陣是對角陣，SEM的這種處理不僅使得方法的精度提高，而且并行性好，在3D石油勘探領域及大尺度模擬方面都有很好的表現(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)在求解中子輸運方程時提出，其后，經過Cockburn and Shu(1989, 2001)，Rivière 和 Wheeler(2003)等人的不斷發(fā)展，以基于數(shù)值通量的龍格-庫塔間斷有限元方法(Runge-Kutta discontinuous Galerkin method, 簡稱RKDGM)及基于內部罰函數(shù)的間斷有限元方法(Interior penalty discontinuous Galerkin method，簡稱IP DGM)為代表的一系列DGM方法被提出，在計算流體力學、計算熱力學、計算電磁學、計算地球物理等領域得到了廣泛應用(e.g., Li, 2006; Hesthaven and Warburton, 2008; Rivière, 2008; 張鐵, 2015; 孟雄等, 2015).間斷有限元方法是傳統(tǒng)有限元方法的一種推廣，它最主要的特征是未知量以及基函數(shù)是在每個網格單元上獨立定義的，在其它網格單元上取值為0，這種定義使得它具有許多良好的性質.相對于FDM，它的網格剖分靈活使得它可以處理任意復雜邊界；相對于FEM，它可以使用任意階基函數(shù)從而可以達到任意階精度，避免了FEM在高階格式上構造的困難，且具有良好的并行性；相對于SEM，它可以采用更靈活的網格剖分——非相容(non-conforming)網格.特別地，DGM允許變量在單元之間存在間斷，因而特別適合模擬強地面運動問題，它不僅能適應非平面斷層及速度反差較大的介質，而且能有效避免滑移率時間序列中存在的虛假高頻振蕩的影響(De La Puente et al., 2009; Pelties et al., 2012).另外，DGM引入的數(shù)值頻散較小，且可以采用非均勻網格，因而特別適用于大尺度非均勻介質中的波場模擬.DGM的質量和剛度矩陣可以提前計算并存儲，因此在實際計算過程中不需要計算質量矩陣和剛度矩陣，這種無積分(quadrature-free)的技巧使得計算量大為減少(Atkins and Shu, 1998).然而，盡管DGM可以采用無積分技巧，但是其計算量和存儲量相對于FEM及SEM還是大很多，這是由于DGM引入了較多的自由度從而導致存儲量增大，而且相比于FEM及SEM來說，DGM需要更小的時間步長才能保持格式穩(wěn)定(de Basabe et al., 2008)，進一步導致了計算量的增大.

DGMs在計算地球物理領域得到了廣泛應用.其中應用最為廣泛的是由K?ser and Dumbser(2006)提出的任意階導數(shù)的時間推進步間斷有限元方法(arbitrary high-order derivatives time stepping discontinuous Galerkin method, 簡稱ADER-DGM),這種方法基于迎風數(shù)值通量，采用具有任意階精度的Lax-Wendroff時間推進方式，從而使得時間精度和空間精度均能達到任意階精度.ADER-DGM及其衍生方法已被成功用于數(shù)值求解3D彈性、粘彈性、孔隙介質、流固耦合等模型的波傳播問題及地震斷層破裂數(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用于求解非均勻介質、粘彈性介質、聲波-彈性波界面等復雜介質的波傳播模擬及偏移成像問題(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主要針對于求解雙曲方程，因此在求解地震波方程領域使用更為廣泛，本研究也主要基于數(shù)值通量形式的DGM.另外，本文中提到的DGM方法均指不帶限制器的DGM，因為地震波動方程大部分是線性方程，Cockburn 和 Shu(2001)指出當雙曲系統(tǒng)為線性系統(tǒng)時，可以不加限制器，但是，如果考慮的是非線性方程或者解存在間斷時，必須使用限制器或者特殊的數(shù)值通量才能保證數(shù)值格式的精度(Chabot et al., 2018)，這已超出本文的研究范圍，在此也不作討論.

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

1 WRKDG方法

1.1 空間離散

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

(1)

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

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

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

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

(2)

即：

(3)

(4)

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

(5)

其中F(u)表示通量，B表示耗散項，當B=0時表示無耗散，方程(5)退化為普通的聲波方程.需要指出的是，本文發(fā)展的方法也適用于一階速度-應力方程.在(5)式兩邊同乘以試驗函數(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)，因此只需提前計算好參考單元上的所有質量矩陣和剛度矩陣.我們選取如圖1中所示的參考單元(Dumbser and K?ser, 2006)，且假設參考單元內的坐標軸三分量為：ξ，η和ζ.對于任意四面體單元，均將之通過坐標變換到如圖1所示的參考單元，圖1a中的頂點與圖1b中的頂點嚴格對應，具體的變換過程可參考附錄A.我們仿照Dumbser 和 K?ser(2006)的流程計算所有的體質量矩陣和剛度矩陣，以及面與面之間關聯(lián)的質量矩陣.具體的過程可參考附錄A.

圖1 一般四面體單元變換到參考單元示意圖(Dumbser and K?ser, 2006)，其中參考單元內1、2、3和4這四個頂點的坐標分別是(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 時間離散

(14)

(15)

(16)

(17)

圖2 Von Neumann分析中用到的網格構型.在(a)中所示的規(guī)則六面體剖分的基礎上，每個六面體中有圖(b)中所示六個四面體網格(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分析.為此，假設一個簡諧波傳播經過3D區(qū)域，并假定此區(qū)域是均勻、無界區(qū)域，且方程右端為無源項.我們主要考察四面體網格剖分，在如圖2a中所示的規(guī)則六面體剖分的基礎上(Mulder et al., 2014; Ferroni et al., 2017; He et al., 2020)，在一個立方體單元內有6個四面體單元，如圖2b所示，其中每個四面體及其面的順序參考Mulder 等(2014)編號.我們假設如下平面簡諧波在所考慮的區(qū)域內傳播：

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

+kcosθlh-iωt)],

(18)

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

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

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

maxα

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

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

α≥0.

(19)

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

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

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

(20)

其中因子0.3597是對應圖2中，當立方體邊長為h時，此時最小的四面體內切球體直徑約為d≈0.3597h.

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

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

(21)

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

(22)

我們在上文中已給出關于Δt≤αmaxh/c的結果.利用求解常微分方程的數(shù)值分析理論(李慶楊等, 2008)，我們知道對于加權的Runge-Kutta時間離散格式(15)，當η=0.5時，要使得格式穩(wěn)定，必須滿足如下條件：

(23)

從而可得

(24)

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

(25)

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

λ2≈r.

(26)

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

(27)

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

(28)

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

(29)

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

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

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

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

(30)

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

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

圖3 在θ=π/2時，P1階和P2階WRKDG方法取φ=0,π/12,π/4時數(shù)值頻散R隨采樣率δ的變化情況.(a—b)P1階WRKDG方法，(c—d)P2階WRKDG方法，(a)、(c)對應耗散參數(shù)r=0，而(b)、(d)對應耗散參數(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)，在δ較小時D′ Alembert介質存在比較大的頻散，隨著δ的增加，R的值接近1，也就是說，隨著δ的增加，頻散變小.對比r=0與r=10這兩種情形，我們可以得出下述結論：在δ(波數(shù))較小時，D′ Alembert介質中存在主要由耗散參數(shù)r引起的頻散；隨著δ(波數(shù))的增加，耗散參數(shù)引起的頻散變小，而數(shù)值算法引起的數(shù)值頻散占主導.

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

圖4 在θ=π/2時，P1階和P2階WRKDG方法取φ=0,π/12,π/4時數(shù)值耗散S隨采樣率δ的變化情況.(a—b)P1階WRKDG方法，(c—d)P2階WRKDG方法，(a)、(c)對應耗散參數(shù)r=0，而(b)、(d)對應耗散參數(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 并行計算

3.1 網格準備

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

3.2 網格分劃

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

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

3.3 并行流程

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

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

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

圖7 (a) 所考慮進程的輔助網格示意圖，即圖中綠色網格部分，這類網格屬于其他進程的內網格，位于該進程所有網格的邊界處；(b) 圖a的側面圖；(c) 屬于該進程內、作為其他進程的輔助網格，即圖中紅色網格部分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 精度測試及并行效率分析

我們考慮一個帶有解析解的模型來驗證該方法的正確性和收斂性.對于無源的方程(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)

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

(33)

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

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

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

5 波場模擬

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

5.1 均勻介質模型

在這個例子中，我們考查D′Alembert介質中聲波在3D均勻介質中的傳播.選取計算區(qū)域為0≤x,y,z≤2 km的立方體區(qū)域，我們將其離散為3072000個四面體，每個四面體邊長平均約為25 m.聲波速度c=3 km·s-1, 時間步長取Δ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處，在坐標(1.35,1.35,1.35)km處設置一個接收點用于記錄波形信息.我們首先考查無耗散情形，即r=0.圖9給出了T=0.5 s時刻接收點接收到的歸一化的波形記錄，圖中紅色實線是用Cagnidard-de-Hoop方法(Aki and Richards，2002)計算得到的解析解，而藍色虛線及黑色實線分別表示利用P2及P1方法得到的數(shù)值解.從圖中可以看出，P1方法出現(xiàn)少許數(shù)值頻散，而P2方法與解析解吻合較好，這說明了提高算法精度有助于降低數(shù)值頻散.圖10給出了P2方法的波場快照圖，此時波已經傳至邊界.波場快照圖中無明顯可見數(shù)值頻散.

圖9 在耗散參數(shù)r=0的3D均勻介質模型中，T=0.5 s 時刻接收點處的歸一化的波形記錄圖，圖中紅色實線代表解析解，藍色虛線及黑色實線分別表示利用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均勻介質模型中，使用P2方法計算得到的T=0.5 s 時刻的波場快照圖Fig.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 在均勻介質模型中，不同耗散參數(shù)r=0、2、4、8和16對應的接收器處的聲波波形記錄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介質中波形記錄的波谷處的振幅和衰減系數(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ī)則幾何模型

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

圖12 非規(guī)則幾何模型示意圖，在計算區(qū)域0≤x,y,z≤2中，有一個球狀區(qū)域，球中心坐標(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) 球體部分四面體網格的3D示意圖； (b) 二維剖面y=0處的網格剖分示意圖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 時刻的波場快照圖，其中(a)對應于無耗散情形r=0， 而(b)對應于耗散參數(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模型

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

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

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

6 結論

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

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

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

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

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

(A1)

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

(A2)

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

(A3)

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

(A4)

例如，要計算原質量矩陣M1，可以利用關系式，M1=|J|M′1.

(A5)

(A6)

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

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

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

(A7)

其中Lij表示以頂點i為起點，j為終點的向量.

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

表A2 (a)三維坐標軸ξ，η，和ζ與面積分用到的參變量和τ之間的對應關系； (b)對于不同的h，在Ωi中的參變量和τ與相鄰單元Ωj中參變量′和τ′的對應關系(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)三維坐標軸ξ，η，和ζ與面積分用到的參變量和τ之間的對應關系； (b)對于不同的h，在Ωi中的參變量和τ與相鄰單元Ωj中參變量′和τ′的對應關系(Dumbser and K?ser, 2006)

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