地球物理学报  2021, Vol. 64 Issue (3): 876-895   PDF    
基于非结构网格求解三维D'Alembert介质中声波方程的并行加权Runge-Kutta间断有限元方法
贺茜君1, 杨顶辉2, 仇楚钧2, 周艳杰1, 常芸凡2     
1. 北京工商大学数学与统计学院, 北京 100048;
2. 清华大学数学科学系, 北京 100084
摘要:间断有限元方法(Discontinuous Galerkin method,简称DGM)在求解地震波动方程时具有低数值频散、网格剖分灵活等优点,因此,为适应数值模拟对模拟精度和复杂地质结构的要求,本文提出一种新的加权Runge-Kutta间断有限元(weighted Runge-Kutta discontinuous Galerkin,简称WRKDG)方法,用于求解三维D'Alembert介质中声波方程.本文不仅详细推导了其数值格式,特别地,根据常微分方程理论给出了满足数值稳定性条件的一般经验公式,并首次对该方法的数值频散和耗散进行了深入分析,且考虑了耗散参数对结果的影响.同时,我们也对该方法进行了精度测试,并分析了3D情形下WRKDG方法的并行加速比,结果表明3D WRKDG方法具有良好的并行性.最后,我们给出了包含均匀模型、非规则几何模型以及非均匀Marmousi模型在内的数值模拟算例.结果表明,该方法不仅计算准确,能与解析解很好地吻合,且能有效模拟包含球体在内的非规则模型及非均匀Marmousi模型中的衰减声波波场.数值模拟实验进一步验证了WRKDG方法在求解三维D'Alembert介质中声波方程时的正确性和有效性,并获得了对这种强衰减介质中波传播特征的规律性新认识.
关键词: 间断有限元方法      三维      数值频散      D'Alembert介质      并行效率      强衰减     
A parallel weighted Runge-Kutta discontinuous galerkin method for solving acousitc wave equations in 3D D'Alembert media on unstructured meshes
HE XiJun1, YANG DingHui2, QIU ChuJun2, ZHOU YanJie1, CHANG YunFan2     
1. School of Mathematics and Statistics, Beijing Technology and Business University(BTBU), Beijing 100048, China;
2. Department of Mathematical Sciences, Tsinghua University, Beijing 100084, China
Abstract: Discontinuous Galerkin method (DGM) is a widely used numerical algorithm. It has the advantages of high accuracy, flexibility in dealing with boundary conditions, easy parallelism, and small numerical dispersion when solving seismic wave equations. In order to satisfy the numerical simulation for accuracy and complex geological structures, in this paper, we suggest a weighted Runge-Kutta discontinuous Galerkin (WRKDG) method for solving the acoustic wave equation in three-dimensional (3D) medium with strong attenuation—D'Alembert medium on unstructured meshes. The numerical scheme is derived in detail, and the general numerical stability conditions are presented based on the theory of ordinary differential equations. The numerical dispersion and dissipation of WRKDG method are also investigated for the first time, including the influence of dissipation parameters on the analysis results. In addition, we carry out a convergence test of this method, and analyze the parallel speedup ratio of the WRKDG method in 3D case. The results show that the 3D WRKDG method has good parallel capabilities. Finally, we present several numerical examples in complex media with strong attenuation, including an homogeneous model, an irregular geometric model, and the heterogeneous Marmousi model. The results show that the method is not only accurate and in good agreement with the analytical solution, but also can effectively simulate the acoustic wave field in irregular model including sphere and heterogeneous Marmousi model. Finally, we present several numerical examples. Numerical results further verify the correctness and effectiveness of the 3D WRKDG method in solving the scalar wave equation in D'Alembert medium, and they clearly show the wave propagation characteristics of this strong attenuation medium.
Keywords: Discontinuous Galerkin method    three-dimensional    numerical dispersion    D'Alembert medium    parallel computing    strong attenuation    
0 引言

基于地震波动方程的正演是计算地球物理学的重要研究内容,它不仅能为我们提供精确的波场模拟结果,而且也是基于波动方程的地震勘探和地震学的重要基础(牟永光和裴正林, 2005).在最近几十年里,随着计算机能力的快速提升,涌现了许许多多优秀的数值算法,例如有限差分法(Finite difference method,简称FDM)、有限元法(Finite element method,简称FEM)、伪谱法(Pseudo-spectral method,简称PSM)、谱元法(Spectral element method,简称SEM)、间断有限元法(Discontinuous Galerkin method,简称DGM)等.这些算法均能满足我们对于数值算法的部分要求,它们也有其独特的优势.对于数值求解3D波动方程来说,这些算法都得到了广泛应用.FDM编程简单、计算速度快,且存储量小,是求解3D波动方程最为常用的一种数值算法(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)和西方学者独立提出,它能处理复杂区域和边界条件,在求解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)来处理全部波场中的空间导数, 精度很高,且数值频散小,在求解波动方程领域也得到了广泛应用(Furumura et al., 1998; Igel, 1999; Klin et al., 2010; Carcione et al., 2018),但是3D情形的PSM需要利用全局数据进行FFT,因此并行性较差(Pelz, 1991).SEM是计算地球物理领域近十年来发展最为迅速的数值算法,它具有FEM的基本特征,在单元内部的基函数基于特定积分节点--Gauss-Lobatto-Legendre(GLL)点,因此在使用GLL积分准则的情况下质量矩阵是对角阵,SEM的这种处理不仅使得方法的精度提高,而且并行性好,在3D石油勘探领域及大尺度模拟方面都有很好的表现(e.g., Komatitsch and Vilotte, 1998; Komatitsch and Tromp, 1999; Komatitsch et al., 1999; Tromp et al., 2008; Liu et al., 2017).

DGM也是近十年来迅速发展起来的一种数值算法,它最早是由Reed和Hill(1973)在求解中子输运方程时提出,其后,经过Cockburn and Shu(1989, 2001),Rivière和Wheeler(2003)等人的不断发展,以基于数值通量的龙格-库塔间断有限元方法(Runge-Kutta discontinuous Galerkin method, 简称RKDGM)及基于内部罚函数的间断有限元方法(Interior penalty discontinuous Galerkin method,简称IP DGM)为代表的一系列DGM方法被提出,在计算流体力学、计算热力学、计算电磁学、计算地球物理等领域得到了广泛应用(e.g., Li, 2006; Hesthaven and Warburton, 2008; Rivière, 2008; 张铁, 2015; 孟雄等, 2015).间断有限元方法是传统有限元方法的一种推广,它最主要的特征是未知量以及基函数是在每个网格单元上独立定义的,在其它网格单元上取值为0,这种定义使得它具有许多良好的性质.相对于FDM,它的网格剖分灵活使得它可以处理任意复杂边界;相对于FEM,它可以使用任意阶基函数从而可以达到任意阶精度,避免了FEM在高阶格式上构造的困难,且具有良好的并行性;相对于SEM,它可以采用更灵活的网格剖分--非相容(non-conforming)网格.特别地,DGM允许变量在单元之间存在间断,因而特别适合模拟强地面运动问题,它不仅能适应非平面断层及速度反差较大的介质,而且能有效避免滑移率时间序列中存在的虚假高频振荡的影响(De La Puente et al., 2009; Pelties et al., 2012).另外,DGM引入的数值频散较小,且可以采用非均匀网格,因而特别适用于大尺度非均匀介质中的波场模拟.DGM的质量和刚度矩阵可以提前计算并存储,因此在实际计算过程中不需要计算质量矩阵和刚度矩阵,这种无积分(quadrature-free)的技巧使得计算量大为减少(Atkins and Shu, 1998).然而,尽管DGM可以采用无积分技巧,但是其计算量和存储量相对于FEM及SEM还是大很多,这是由于DGM引入了较多的自由度从而导致存储量增大,而且相比于FEM及SEM来说,DGM需要更小的时间步长才能保持格式稳定(de Basabe et al., 2008),进一步导致了计算量的增大.

DGMs在计算地球物理领域得到了广泛应用.其中应用最为广泛的是由Käser and Dumbser(2006)提出的任意阶导数的时间推进步间断有限元方法(arbitrary high-order derivatives time stepping discontinuous Galerkin method, 简称ADER-DGM), 这种方法基于迎风数值通量,采用具有任意阶精度的Lax-Wendroff时间推进方式,从而使得时间精度和空间精度均能达到任意阶精度.ADER-DGM及其衍生方法已被成功用于数值求解3D弹性、粘弹性、孔隙介质、流固耦合等模型的波传播问题及地震断层破裂数值模拟中(Dumbser and KäserKäser, 2006; Dumbser et al., 2007; Käser and Dumbser, 2008; De La Puente et al., 2007, 2008).此外,也涌现了许多基于其它形式的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).由于基于通量函数的DGM主要针对于求解双曲方程,因此在求解地震波方程领域使用更为广泛,本研究也主要基于数值通量形式的DGM.另外,本文中提到的DGM方法均指不带限制器的DGM,因为地震波动方程大部分是线性方程,Cockburn和Shu(2001)指出当双曲系统为线性系统时,可以不加限制器,但是,如果考虑的是非线性方程或者解存在间断时,必须使用限制器或者特殊的数值通量才能保证数值格式的精度(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规则的六面体单元.本研究针对3D非结构网格,发展了求解3D D′Alembert介质中的并行WRKDG方法.D′Alembert介质是一种粘弹性介质,它直接对运动方程加入粘性项来刻画粘性,因而具有简洁的表达式.Li等(2015)Cai等(2017)Lähivaara和Huttunen(2010)等人对这种介质进行了数值模拟研究,本文的研究也是基于此.另外,我们对并行算法的设计及编程也给出了较为详细的分析.特别地,我们根据常微分方程理论推导了3D D′Alembert介质中满足数值稳定性条件的一般经验公式.同时,由于针对基于数值通量的DGM,目前还没有相关的数值稳定性和数值频散、耗散的分析结果,因此本文首次对该方法的数值频散和耗散进行了详细的研究,且考虑了耗散参数对结果的影响.最后,我们给出了包含均匀模型、非规则几何模型以及非均匀Marmousi模型在内的数值模拟算例以说明方法的有效性.

1 WRKDG方法 1.1 空间离散

我们考虑3D D′ Alembert介质中声波方程(牛滨华和孙春岩, 2007; Lähivaara and Huttunen, 2010):

(1)

其中u是位移场,c是声波速度,f(t)是源项,r>0是D′ Alembert介质中引入的耗散参数.参数r与频率ω有如下关系式(牛滨华和孙春岩, 2007):

其中,Q为D′Alembert介质纵波的品质因子.当r/ω≤1时,我们有rωQ-1.设3D求解区域为Ω.我们将区域Ω划分为不重叠的子区域{Ωi}.由于本文主要研究非结构网格,因此主要采取四面体网格剖分.需要指出的是,为方便起见,本文仅考虑相容网格的情况,也就是不允许有“悬点”的存在.

DGM基于Galerkin近似,所以我们首先需选取试验函数空间.我们使用的试验函数空间是多项式空间Vh={vL1(Ω): v|ΩiPκ(Ωi)},其中Pκ(Ωi)表示区域Ωi上次数不超过κ次的多项式.从定义中可以看出,测试函数vΩi以外的区域上不定义或者令其为0,因此,它不必在整个区域上连续,可以在子区域Ωi之间有间断,这也是间断有限元与传统有限元最大的区别.

为了将方程(1)改写成一阶双曲系统的形式,我们引入三个变量pqs,且pqs满足条件pt=ux, qt=uy, st=uz,其中uxuyuz分别表示位移u对空间变量xyz的偏导数.则方程(1)经过一次时间积分后,可改写成如下形式:

(2)

即:

(3)

其中▽=(x, y, z).如果我们引入如下记号:

(4)

则可以将3D D′ Alembert介质中声波方程化成适于用DGM求解的如下形式:

(5)

其中F(u)表示通量,B表示耗散项,当B=0时表示无耗散,方程(5)退化为普通的声波方程.需要指出的是,本文发展的方法也适用于一阶速度-应力方程.在(5)式两边同乘以试验函数v,并在子区域Ωi上积分,利用格林公式,我们得到

(6)

其中∂Ωi表示区域Ωi的边界,n=(n1, n2, n3)T表示边界∂Ωi处的单位外法向量.由于u∂Ωi处有可能间断,所以(6)式中的F(un没有定义,我们将F(un替换为数值通量:F(un|∂Ωi=, n),其中uint表示u在区域Ωi内部逼近∂Ωi处的值,uext表示u在区域Ωi外部逼近∂Ωi处的值.这里我们采用一种常用通量,即Local Lax-Friedrichs(LLF)通量(Cockburn and Shu, 1989, 2001),其表达式为

(7)

其中C∂Ωi的绝对值最大的特征值,通常情况下可取相邻两个单元的最大速度.将数值通量(7)代入方程(6),我们得到DGM方法的弱形式:

(8)

在上式中,由于我们考虑的是线性问题,不妨将通量F(u)写成F(u)=(A1u, A2u, A3u)的形式(Dumbser and Käser, 2006),其中

(9)

由此我们可将数值通量写成下列形式:

(10)

其中,

(12)

I表示单位矩阵.令u具有如下的基函数展开形式:

(12)

其中wli(x)表示Ωi内的基函数,Cli(t)表示Ωi内的未知系数向量,N表示基函数的个数.一般来说,基函数有两大类(de Basabe et al., 2008),其中一类是modal型,另一类是nodal型.modal型的基函数只要求基函数能生成测试函数空间,不需要在网格节点上定义;nodal型基函数一般采用在一些特定节点上(一般是Gauss型积分节点)定义的Lagrange基函数.一般说来,只要数值积分的精度能得到保证,选什么类型的基函数并不影响数值结果,尤其是本文采用的是无积分(quadrature-free)的DGM(Atkins and Shu, 1998),只需要提前计算好参考单元上的所有质量矩阵和刚度矩阵,在计算过程中不需要再计算质量矩阵等其他矩阵,所以只需要保证参考单元上的积分是正确的即可.例如,我们可以采用3D情形下的参考单元内有序完备(order-complete)的多项式{xαyβzχ|0≤α+β+χκ}(Rivière, 2008)作为基函数,该种基函数比较简单,此时N=(κ+1)(κ+2)(κ+3)/6.例如,当κ=1,即空间精度为2阶时,N的值是4,我们称该方法为2阶DGM;当κ=2,即空间精度为3阶时,N的值是10,我们称该方法为3阶DGM.当然我们也可以选择正交的基函数,可参考De La Puente (2008)中的附录部分.

将(12)和(10)式代入(8)式中,我们可以得到如下半离散化的方程:

(13)

本研究中采用的是无积分(quadrature-free)的DGM(Atkins and Shu, 1998),因此只需提前计算好参考单元上的所有质量矩阵和刚度矩阵.我们选取如图 1中所示的参考单元(Dumbser and KäserKäser, 2006),且假设参考单元内的坐标轴三分量为:ξηζ.对于任意四面体单元,均将之通过坐标变换到如图 1所示的参考单元,图 1a中的顶点与图 1b中的顶点严格对应,具体的变换过程可参考附录A.我们仿照Dumbser和Käser(2006)的流程计算所有的体质量矩阵和刚度矩阵,以及面与面之间关联的质量矩阵.具体的过程可参考附录A.

图 1 一般四面体单元变换到参考单元示意图(Dumbser and KäserKä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äserKäser, 2006)with four vertices (0, 0, 0), (1, 0, 0), (0, 1, 0), and (0, 0, 1)
1.2 时间离散

经过上述空间离散之后,我们得到半离散化的方程(13).方程(13)是关于未知量Cli(t)的常微分方程,为了表述上的方便,我们令C表示所有单元上的未知量构成的向量,则方程(13)可以简写成下述形式:

(14)

其中L是方程(13)中与空间离散有关的线性算子,F是和震源相关的项.关于求解方程(14)有许多时间离散格式,在此我们选取一种加权的Runge-Kutta格式(He et al., 2015; 张金波等, 2018)对方程(14)进行时间推进,其格式为

(15)

其中η∈[0, 1]是加权系数.在下文中若无特殊说明,我们均取η=0.5进行计算.K1(n)K2(n)可以用如下公式计算:

(16)

其中.得到K(n)的近似值之后,令T(n)=C(n)+(1-2r)△tK(n),按如下公式求得K1(n)K2(n)

(17)

将如上所得的K1(n)K2(n)K1(n)K2(n)代入方程(15)中,形成完整的时间推进格式.相比于传统的Runge-Kutta时间离散格式,加权的Runge-Kutta格式能提高最大库朗数,从而使得计算的时间步长增大,迭代步数减小.

2 Von Neumann分析

3D情形下的数值稳定性、数值频散和耗散主要基于Von Neumann分析.为此,假设一个简谐波传播经过3D区域,并假定此区域是均匀、无界区域,且方程右端为无源项.我们主要考察四面体网格剖分,在如图 2a中所示的规则六面体剖分的基础上(Mulder et al., 2014; Ferroni et al., 2017; He et al., 2020),在一个立方体单元内有6个四面体单元,如图 2b所示,其中每个四面体及其面的顺序参考Mulder等(2014)编号.我们假设如下平面简谐波在所考虑的区域内传播:

(18)

图 2 Von Neumann分析中用到的网格构型.在(a)中所示的规则六面体剖分的基础上,每个六面体中有图(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)

其中k是波数,ω是频率,h是立方体的边长,θφ决定了平面波的传播方向,θ∈[0, π]表示波传播方向与z轴的夹角,而φ∈[0, 2π]则表示波传播方向在xy平面内的投影与x轴的夹角.将(18)式代入数值格式中,即可用于分析数值稳定性和数值频散及耗散.

2.1 数值稳定性分析

为了保持数值算法的稳定性,我们引入库朗数条件,也称Courant-Friedrichs-Lewy(简称CFL)条件:α=ct/hαmax,其中α是库朗数,αmax即是最大库朗数(也称CFL条件数).不妨令Λ表示代入简谐波(18)之后的数值格式增长矩阵的特征值,则Λ与波数k、库朗数α、传播方向θφ有关.要使得数值格式稳定,必须满足|Λ|≤1对所有的kh∈[0, 2π]、θ∈[0, π]和φ∈[0, 2π]都成立.这在数学上可以用下面一个优化问题来表示:

(19)

上述求α的最大值问题是一个非线性优化问题,通常情况下不容易求出αmax的解析表达式.一般我们通过数值遍历算法来求其近似解,使α从0开始以小步长增长,直至|Λ(kh, α, θ, φ)|>1则跳出循环,获得αmax的值.

在本文中,我们仅考虑基函数为1、2次多项式的情形,对应空间精度分别为2、3阶.首先,我们考虑不带耗散的情形,即在方程(1)中r = 0的情形.通过数值求解上述非线性优化问题,我们得到当η=0.5时的WRKDG方法的最大库朗数αmax:对于P1阶WRKDG方法,ct/hαmax≈0.244;对于P2阶WRKDG方法,ct/hαmax≈0.144.

以上分析给出的最大库朗数αmax中所采用的网格步长h图 2a中立方体边长,若是采用图 2b中四面体内切球体直径d,则此时的时间步长所需满足的条件是

(20)

其中因子0.3597是对应图 2中,当立方体边长为h时,此时最小的四面体内切球体直径约为d≈0.3597h.

对于带耗散情形,即方程(1)中r>0的情形,也可以根据(19)式中同样的流程进行分析,求出3D D′Alembert介质中的库朗数条件.然而此时由于加入了耗散参数r,我们希望能得到一个包含r的显式的αmax的表达式.为此,在已知声波方程库朗数条件的基础上,我们进行如下分析.

我们首先注意到,经过空间离散之后形成的半离散系统(13)可以分解为两部分,一部分是无耗散的双曲型系统,另一部分则与耗散有关.因此,相应地,我们可以将常微分方程组(14)分成两部分,忽略震源项后,写成如下形式(Carcione and Quiroga-Goode, 1995):

(21)

其中,算子L1对应无耗散情形的声波传播,算子L2对应系统中的耗散项,实际上,L2只与耗散参数r有关.我们将算子LL1L2对应的谱半径分别记作λλ1λ2.则由(21),我们可得λλ1+λ2.对于双曲型系统:

(22)

我们在上文中已给出关于△tαmaxh/c的结果.利用求解常微分方程的数值分析理论(李庆杨等, 2008),我们知道对于加权的Runge-Kutta时间离散格式(15),当η=0.5时,要使得格式稳定,必须满足如下条件:

(23)

从而可得

(24)

对于带耗散的系统:

(25)

由于算子L2只与矩阵耗散参数r有关,所以算符的谱半径等于r,即

(26)

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

(27)

因此,我们得到利用WRKDG方法求解3D D′ Alembert介质中的数值稳定性条件的经验公式:

(28)

方程(28)中给出的时间步长限制是保持计算稳定的充分条件,但不是必要条件.我们在表 1列出了当波速c=1 km·s-1, 随hr变化时,P1阶和P2阶WRKDG方法的真实最大库朗数(αmax)D′Ale与从经验公式(28)获得的值(αmax)′D′Ale进行对比,从表中可以看出两者之间相差不大.因此,在实际应用中,(28)式是一个合理的估计.若是采用四面体内切球体直径d,则(28)式可写为

(29)

表 1 3D D′ Alembert介质中P1阶和P2阶WRKDG方法的真实最大库朗数(αmax)D′Ale与从经验公式(28)获得的值(αmax)′D′Ale的对比 Table 1 Comparisons of the actual maximum Courant numbers (αmax)D′Aleand the values (αmax)′D′Ale given by empirical formula (28) for 3D P1 and P2 WRKDG methods in D′ Alembert media

另外,从方程(28)我们可以看出,随着耗散参数r的增大,时间步长减小;尤其当r的值很大时,此时系统(21)是一个刚性(stiff)系统,一般的时间推进方法将导致极小的时间步长,因而不再适用,此时应寻找特殊的时间推进方法.

2.2 数值频散、耗散分析

在本节中,我们将讨论求解D′Alembert介质中声波方程的3D WRKDG方法的数值频散和数值耗散情况.对于方程(18)中的简谐波表达式,当代入数值格式中,如果我们假定波数k是实数,则一般说来角频率ω是复数,不妨假设ω=ωr-iωi,其中实部ωr表示ω中与频散有关的量,非负的虚部ωi表示ω中与耗散有关的量.数值波速可由cnum=ωr/k估计.我们引入采样率δ=kh/(2π),并将数值频散R定义为数值速度与真实物理速度之比,则R的表达式为

(30)

其中α=ct/h是库朗数.我们将振幅在一个时间步内的衰减记作S=eωiΔtS可以反映D′ Alembert介质中的衰减情况.

图 3给出了用WRKDG方法计算的数值频散情况.我们设置参数c=2 km·s-1h=0.05 km,α=ct/h=0.1.在3D情形,数值频散R的大小与传播方向θφ有关,在图 3中,我们仅给出了θ=π/2、φ=0, π/12, π/4情况下,WRKDG方法的数值频散随采样率δ的变化图.注意此时由于h固定,δ的大小与波数成正比.图 3中(a-b)、(c-d)分别表示P1阶、P2阶WRKDG方法,(a)、(c)表示无耗散r=0的情形,而(b)、(d)表示耗散参数r=10的情形.从图 3中我们可以明显观察到数值频散的各向异性特征.引起数值频散各向异性的因素较多,数值算法本身、算法精度、网格剖分方式、网格步长大小等都会影响数值频散各向异性的产生及幅度.例如,在进行数值频散分析时,引入了方位角θφ,当空间网格分布存在不对称性时,会导致数值频散出现各向异性特征.其次,不同算法精度产生的数值频散各向异性也不相同.一般来说,减小网格的各向异性、提高算法精度、减小网格步长可以有效降低数值频散的各向异性.

图 3θ=π/2时,P1阶和P2阶WRKDG方法取φ=0, π/12, π/4时数值频散R随采样率δ的变化情况.(a-b)P1阶WRKDG方法,(c-d)P2阶WRKDG方法,(a)、(c)对应耗散参数r=0,而(b)、(d)对应耗散参数r=10 Fig. 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中发现,在δ较小时D′ Alembert介质存在比较大的频散,随着δ的增加,R的值接近1,也就是说,随着δ的增加,频散变小.对比r=0与r=10这两种情形,我们可以得出下述结论:在δ(波数)较小时,D′ Alembert介质中存在主要由耗散参数r引起的频散;随着δ(波数)的增加,耗散参数引起的频散变小,而数值算法引起的数值频散占主导.

图 4显示了振幅在一个时间步Δt内的耗散S=eωit随采样率的变化图.图中参数的选取和上文中数值频散分析中相同.当r=0时(图 4a图 4c),S表示声波介质中的衰减情况,此时S完全由数值离散方法引起,而当r=10时(图 4b图 4d),S的大小由数值离散方法引起的数值耗散与D′ Alembert介质中的耗散参数r共同决定.图 4也体现了数值耗散的各向异性特征.在图 4中,对比耗散因子r=0与r=10这两种情形,我们可以观察到r=10时的耗散曲线有一个整体的下降,其值约为0.9876,约等于D′ Alembert介质中理论耗散因子ert/2,体现了D′Alembert介质的衰减特性.

图 4θ=π/2时,P1阶和P2阶WRKDG方法取φ=0, π/12, π/4时数值耗散S随采样率δ的变化情况.(a-b)P1阶WRKDG方法,(c-d)P2阶WRKDG方法,(a)、(c)对应耗散参数r=0,而(b)、(d)对应耗散参数r=10 Fig. 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 网格准备

结构网格(六面体)和非结构网格(四面体)均可用于实际计算,一般说来,结构网格精度较高,但是对于复杂几何模型来说,生成结构网格会花费较大代价;非结构网格计算精度不如结构网格,但是其网格生成较简单.本研究中选取非结构网格--四面体网格.网格剖分可以使用开源软件或商业软件.当模型比较简单或网格量比较小的时候选择开源软件,当模型比较复杂或者问题规模较大时,优先选择商业软件.我们使用开源软件或者商业软件生成了四面体网格,获得了所有网格节点的信息以及单元连接关系矩阵,为了计算的便利性,我们还需要计算面与面之间的关联矩阵、面与单元之间的关联矩阵、每个面中第一个顶点对应的于相邻面中的局部编号矩阵.这些矩阵的计算过程可参考Hesthaven和Warburton(2008)的研究.

3.2 网格分划

由于3D情形计算量和存储量均较大,所以要进行并行处理.在并行之前,我们需要对整体的3D网格进行分划,将其分给不同的进程.我们使用开源网格分划程序包Metis(Karypis and Kumar, 1998).在使用时只需要输入单元数、顶点数、单元连接矩阵、进程数等数据,即可用一行命令对网格进行快速分划.例如,图 5显示了Metis分划的结果,图 5a代表将3D区域[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中蓝色区域内的网格即是第一个进程的内网格;第二类网格是该进程的辅助网格,这类网格属于其他进程的内网格,位于进程1所有网格的边界处,如图 7a中的绿色网格部分(图 7b图 7a的一个更清晰的展示),作用是接收来自于其他进程更新后的变量信息;第三类网格是属于该进程内,但是作为其他进程的辅助网格,这类网格在消息传递时需要由本进程传递给其他进程使用,如图 7c中红色网格显示.

图 6 并行算法流程图 Fig. 6 Flow chart of the parallel algorithm
图 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)

开始阶段结束后,如流程图 6中的说明,我们将变量赋初值,然后进入时间迭代.在每个时间迭代步中,首先更新每个进程内网格(即第一类网格)的变量信息.这一步结束后,我们需要进行两步消息传递:将所在进程内属于其他进程的辅助网格(即第三类网格)上的Cn+1发送给相应进程,并接收来自其他进程的Cn+1并存放于辅助网格(即第二类网格)中.这样我们就完成了一步完整的时间迭代.时间迭代结束后,输出数据结果,并利用画图软件等进行画图.

4 精度测试及并行效率分析

我们考虑一个带有解析解的模型来验证该方法的正确性和收敛性.对于无源的方程(1),考虑如下解析解(Cai et al., 2017)

(31)

其中.相应地,对于系统(2),我们有如下解析解:

(32)

在本例中,选取计算区域为0≤x, y, z≤2 km,速度c=2 km·s-1K1=K2=K3=π, ,以t=0 s时刻的值作为初始条件,采用周期边界条件.由于我们主要考察空间离散的误差和收敛精度,因此我们设置时间步长为0.1 ms,此时,由时间离散引起的误差可以忽略,数值误差主要来自于空间离散.我们采用如图 2所示的网格剖分及四面体单元,图中每个小立方体的边长为h,且每个立方体含有6个四面体单元.我们令N表示在xyz三个方向划分的立方体个数,则四面体总数目为6N3.定义L2L1误差为

(33)

其中uex是方程(31)给出的精确解.我们令N=4, 8, 10, 16, 20,在表 2中列出了T=0.1 s时刻在耗散参数r=1以及r=10两种情形下的数值误差和相应的收敛阶.从表 2可以看出,数值误差随着空间步长的减小而减小,说明WRKDG方法是收敛的.由于我们使用了二次基函数,因此得到了预期的三阶空间精度.同时,从表 2中可以看出,随着耗散参数r的增大,误差也随之减小,体现了D′Alembert介质的衰减特性.

表 2 3D D′ Alembert介质中P2阶WRKDG的误差及收敛阶 Table 2 Convergence rates of u for the acoustic wave in 3D D′Alembert medium

接下来,为了考察3D WRKDG算法的并行效率,我们将上述精度测试中N=20的数值模拟程序进行并行化处理,此时计算区域被离散成48000个四面体,利用Metis将网格分别分划给2、4、8、16、32、64个进程,记录运行的CPU时间数.假设每个处理器的计算性能相同,在此条件下,并行程序的加速比(speed-up)可定义为:SP=TS/TP, 其中TS是串行程序运行的时间,TP是并行程序在P个进程下运行的时间,SP在通常情况下总是小于P,因为在并行程序设计时往往会引入一些通信时间及其他管理花费.图 8给出了3D WRKDG方法的并行程序加速比及理想情形下的加速比,从图中可以看出,当进程数比较小的时候,加速比接近理论情形,但是随着进程数越增大,实际加速比越偏离理论情形,这是由于进程数增加引起进程与进程之间的通信时间开销大大增大,同时进程与进程之间等待的时间开销也增大,从而降低了并行效率.本文所使用的计算平台是国家超级计算天津中心的天河TH-1A高性能机群系统,每个节点上有12个主频为2.93 GHz的核,每个节点内存为24 GB.

图 8 3D WRKDG算法的并行加速比(图中线型“-o”表示).其中横轴表示进程数,纵轴表示并行加速比.图中线型“-*”表示理想情形下的并行加速比 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
5 波场模拟

在本节中,我们通过几个数值例子来说明本文所给出的WRKDG方法在求解3D D′Alembert介质中声波传播问题的有效性.考虑到更高阶格式的库朗数条件更为严格,且存储量和计算量也越大,因此,如果没有特殊说明,在本节中我们仅使用空间精度为3阶的P2阶WRKDG方法进行波场模拟.

5.1 均匀介质模型

在这个例子中,我们考查D′Alembert介质中声波在3D均匀介质中的传播.选取计算区域为0≤x, y, z≤2 km的立方体区域,我们将其离散为3072000个四面体,每个四面体边长平均约为25 m.声波速度c=3 km·s-1, 时间步长取Δt≈1.29 ms.震源函数是:

(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)计算得到的解析解,而蓝色虚线及黑色实线分别表示利用P2P1方法得到的数值解.从图中可以看出,P1方法出现少许数值频散,而P2方法与解析解吻合较好,这说明了提高算法精度有助于降低数值频散.图 10给出了P2方法的波场快照图,此时波已经传至边界.波场快照图中无明显可见数值频散.

图 9 在耗散参数r=0的3D均匀介质模型中,T=0.5 s时刻接收点处的归一化的波形记录图,图中红色实线代表解析解,蓝色虚线及黑色实线分别表示利用P2P1方法得到的数值解 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 在耗散参数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

下面考虑D′Alembert介质中耗散参数r=2、4、8和16时的波传播情形.图 11给出了此时接收点处的波形记录图,从图中可以看出,随着耗散参数r的增大,振幅明显变小,这证明了D′Alembert介质中的强衰减效应.为了定量地研究这种衰减效应,我们记录了r=0、2、4、8和16情形下波形记录中波谷处的振幅值,然后相应地除以r=0时波谷振幅,得到振幅的衰减系数.表 3列出了此时观测到的波谷处的振幅及衰减系数.此外,我们也给出了理论的衰减系数以供比较,理论衰减系数由erT/2计算(具体见第3.2节耗散分析),其中T是声波从震源传播到接收器所需的传播时间,容易算出本例中T= .从表 3可以看出,此时观测到的衰减系数与理论衰减系数符合得很好.

图 11 在均匀介质模型中,不同耗散参数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介质中波形记录的波谷处的振幅和衰减系数 Table 3 The amplitudes and attenuation ratios at the trough at the receiver for the acoustic wave in D′Alembert medium
5.2 非规则几何模型

在这个例子中,我们主要利用3D WRKDG方法模拟波在非规则几何模型中的传播,模型如图 12所示,在计算区域为0≤x, y, z≤2 km的立方体中,有一个球状区域,球中心坐标(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.震源函数与方程(34)中相同,其中f0=60 Hz,主频约为25 Hz.图 14给出了在T=0.3 s时刻下的波场快照图,图 14a对应于无耗散情形,而图 14b对应于耗散参数r=4.从图中可以看出,图 14b图 14a暗一些,证明了D′Alembert介质中的衰减效应.

图 12 非规则几何模型示意图,在计算区域0≤x, y, z≤2中,有一个球状区域,球中心坐标(1, 1, 0.5)km,半径0.2 km Fig. 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)对应于耗散参数r=4 Fig. 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,震源函数如方程(34)所示,其中f0=30 Hz,主频约为13 Hz,震源位于(4.577, 0.015, 1.449)km处.模拟中时间步长取Δt≈1.69 ms,D′Alembert介质中耗散参数r=2.图 16给出了T=1.0 s时刻的波场快照,从图中可以看出并无明显可见的数值频散.这说明3D WRKDG方法可以有效模拟复杂非均匀D′Alembert介质中的声波波场.

图 15 3D Marmousi模型 Fig. 15 3D Marmousi model
图 16 对于3D Marmousi模型,T=1.0 s时刻的波场快照图,其中耗散参数r=2 Fig. 16 Snapshots at T=1.0 s for the 3D Marmousi model with dissipation parameter r=2
6 结论

本文将加权Runge-Kutta间断有限元(WRKDG)方法应用于求解3D D′Alembert介质中的声波方程,空间离散采用了基于数值通量的间断有限元公式,时间离散基于对角隐式Runge-Kutta方法,我们通过两步迭代过程将其转化为显式方法,并在时间离散化过程中引入加权因子,最终获得了求解3D D′Alembert介质中声波方程的WRKDG新方法.进一步,我们详细研究了该方法的数值稳定性条件,给出了四面体情形下的最大库朗数.由于D′Alembert介质中耗散参数的引入,我们也推导了一种带有耗散参数的数值稳定性条件经验公式.数值试验表明,该经验公式是一种正确的估计.此外,我们也深入研究了WRKDG方法在四面体情形下的数值频散及数值耗散,研究表明D′ Alembert介质中的数值频散和耗散由耗散参数r及数值算法共同决定,存在一个理论耗散因子ert/2.同时,我们也观察到数值频散和数值耗散具有明显的各向异性特征,这主要是由于所用网格的各向异性特征导致的.

我们通过数值模拟实验证明了WRKDG方法的收敛性,给出了3D并行WRKDG算法基于MPI并行策略下的加速比曲线,从中可以看出WRKDG算法具有良好的并行性能.为了进一步验证3D WRKDG方法的正确性和有效性,我们模拟了声波在D′Alembert介质中均匀、非均匀介质及非规则几何模型中的传播,且针对均匀介质给出了理论耗散因子与观测衰减因子,二者较为吻合.这些结果均表明3D WRKDG方法能够正确且有效地模拟衰减介质中的声波传播,能充分体现D′ Alembert介质中波的衰减特征.

最后需要指出的是,尽管3D WRKDG方法应用了并行策略能够有效节省计算时间,但是其计算量和存储量相对于其它数值方法还是很大,因此,今后我们应重点考虑如何从多种途径联合提高其计算效率,以真正实现复杂介质中大规模地震模拟、逆时偏移和基于波动方程反演成像的快速计算和实际应用,这些都是值得我们继续深入研究的方向.

附录A 方程(13)所需矩阵的具体表达式

根据方程(13),我们需要计算如下四面体上的体积分矩阵:

(A1)

其中,角标lm表示矩阵的lm列元素,上角标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, yz,且假设参考单元内的坐标轴三分量为:ξηζ,则任意四面体均将通过如下坐标变换成如图 1所示的参考单元:

(A2)

易知:dxdydz=|J|dξdηdζ,其中|J|是四面体Ωi的体积,且容易得到如下偏导数的值(Dumbser and KäserKäser, 2006):

(A3)

在(A2)所示的坐标变换下,要计算方程(A1)中的矩阵,只需要在参考单元中计算如下矩阵即可:

(A4)

例如,要计算原质量矩阵M1,可以利用关系式,M1=|J|M1.

为了计算方程(13)中所有面上的积分,我们引入如下四面体所有面上质量矩阵.下面叙述的整体思路可参考Dumbser and KäserKäser(2006),我们下文仅摘录关乎编程的重要细节加以阐述.为此,我们原封不动引用Dumbser and KäserKäser(2006)文中的表 1表 2(对应本文中表 A1表 A2),分别代表四面体四个面的顺序以及在进行面积分的时候参变量对应的取值,四面体的四个面我们不妨记作F1F2F2F4,对应的参考单元上的面分别记作E1E2E3E4.方程(13)中出现了两个面积分项,其中第一项是,第二项是,下面我们分别来进行处理.先看第一项,此项包含四种情况,即,对照图 1表 A1,在F1F2F2F4上,有

(A5)

表 A1 四面体单元所属四个面的定义顺序(Dumbser and KäserKäser, 2006) Table A1 Face Definition on tetrahedrons (Dumbser and KäserKäser, 2006)
表 A2 (a) 三维坐标轴ξη,和ζ与面积分用到的参变量χτ之间的对应关系;(b)对于不同的h,在Ωi中的参变量χτ与相邻单元Ωj中参变量χ′和τ′的对应关系(Dumbser and KäserKä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äserKäser, 2006)

其中,|SFi|表示Fi的面积,表示三角形区域{(χ, τ)|0≤χ≤1-τ, 0≤τ≤1|}.下面我们来处理第二项上的积分,注意,这里的j代表的是与单元Ωi相邻单元的编号.令F=ΩiΩjΩiΩj相交的面有几种情况,根据表 A1规定的面的顺序,FΩi中有可能等于F1F2F3F4,在Ωj中有可能等于F1′F2′F3′F4′(F1′F2′F3′F4′表示按照表 A1规定的Ωj中四个面的顺序).现假设在Ωi中有:F=Fi, ∀1≤i≤4,在Ωj中有F=F′j, ∀1≤j≤4,则FiF′j的连接根据顶点对应的不同情形又可分为三种情况,Dumbser and KäserKäser(2006)表 A2(b)中引入参数h表示这三种情况,h表示在Fi中第一个顶点对应的F′j中顶点的局部编号,1 ≤h≤3.于是计算总共需要考虑4×4×3=48种情形.在表 A2中简明地表达了这48种情形,例如,假如在ΩiF=F1,在ΩjF=F′2,且h=3,则积分项的计算表达式是

(A6)

至此,我们已将方程(13)中所有积分矩阵的表达式阐述完毕.对于四面体的四个面F1F2F3F4的外法向量n1, n2, n3n4的计算,有如下公式,

(A7)

其中Lij表示以顶点i为起点,j为终点的向量.

References
Aki K, Richards P G. 2002. Quantitative Seismology. 2nd ed. Sausalito, California: University Science Books.
Atkins H L, Shu C W. 1998. Quadrature-free implementation of discontinuous Galerkin method for hyperbolic equations. AIAA Journal, 36(5): 775-782. DOI:10.2514/2.436
Cai W J, Zhang H, Wang Y S. 2017. Dissipation-preserving spectral element method for damped seismic wave equations. Journal of Computational Physics, 350: 260-279. DOI:10.1016/j.jcp.2017.08.048
Carcione J M, Quiroga-Goode G. 1995. Some aspects of the physics and numerical modeling of Biot compressional waves. Journal of Computational Acoustics, 3(4): 261-280. DOI:10.1142/S0218396X95000136
Carcione J M, Poletto F, Farina B, et al. 2018. 3D seismic modeling in geothermal reservoirs with a distribution of steam patch sizes, permeabilities and saturations, including ductility of the rock frame. Physics of the Earth and Planetary Interiors, 279: 67-78. DOI:10.1016/j.pepi.2018.03.004
Chabot S, Glinsky N, Mercerat E D, et al. 2018. A high-order discontinuous Galerkin method for 1D wave propagation in a nonlinear heterogeneous medium. Journal of Computational Physics, 355: 191-213. DOI:10.1016/j.jcp.2017.11.013
Cockburn B, Shu C W. 1989. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws Ⅱ: General framework. Mathematics of Computation, 52(186): 411-435.
Cockburn B, Shu C W. 2001. Runge-kutta discontinuous galerkin methods for convection- dominated problems. Journal of Scientific Computing, 16(3): 173-261. DOI:10.1023/A:1012873910884
De La Puente J, Käser M, Dumbser M, et al. 2007. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes-IV. Anisotropy. Geophysical Journal International, 169(3): 1210-1228. DOI:10.1111/j.1365-246X.2007.03381.x
De La Puente J. 2008. Seismic wave simulation for complex rheologies on unstructured meshes[Ph. D. thesis]. Ludwig-Maximilians-Universitat, doi: 10.5282/edoc.8074.
De La Puente J, Dumbser M, Käser M, et al. 2008. Discontinuous Galerkin methods for wave propagation in poroelastic media. Geophysics, 73(5): T77-T97. DOI:10.1190/1.2965027
De Basabe J D, Sen M K, Wheeler M F. 2008. The interior penalty discontinuous Galerkin method for elastic wave propagation: grid dispersion. Geophysical Journal International, 175(1): 83-93. DOI:10.1111/j.1365-246X.2008.03915.x
De La Puente J, Ampuero J P, Käser M. 2009. Dynamic rupture modeling on unstructured meshes using a discontinuous Galerkin method. Journal of Geophysical Research: Solid Earth, 114(B10): B10302. DOI:10.1029/2008JB006271
Dong L G, Ma Z T, Cao J Z. 2000. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation. Chinese Journal of Geophysics (in Chinese), 43(6): 856-864.
Dumbser M, Käser M. 2006. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes-Ⅱ: The three-dimensional isotropic case. Geophysical Journal International, 167(1): 319-336. DOI:10.1111/j.1365-246X.2006.03120.x
Dumbser M, Käser M, Toro E F. 2007. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes-V: Local time stepping and p-adaptivity. Geophysical Journal International, 171(2): 695-717. DOI:10.1111/j.1365-246X.2007.03427.x
Etienne V, Chaljub E, Virieux J, et al. 2010. An hp-adaptive discontinuous Galerkin finite-element method for 3-D elastic wave modelling. Geophysical Journal International, 183(2): 941-962. DOI:10.1111/j.1365-246X.2010.04764.x
Feng K. 1965. Difference schemes based on variational principle. Journal of Applied and Computational Mathematics (in Chinese), 2(4): 238-262.
Ferroni A, Antonietti P F, Mazzieri I, et al. 2017. Dispersion-dissipation analysis of 3-D continuous and discontinuous spectral element methods for the elastodynamics equation. Geophysical Journal International, 211(3): 1554-1574. DOI:10.1093/gji/ggx384
Furumura T, Kennett B L N, Takenaka H. 1998. Parallel 3-D pseudospectral simulation of seismic wave propagation. Geophysics, 63(1): 279-288. DOI:10.1190/1.1444322
Galis M, Moczo P, Kristek J. 2008. A 3-D hybrid finite-difference-finite-element viscoelastic modelling of seismic wave motion. Geophysical Journal International, 175(1): 153-184. DOI:10.1111/j.1365-246X.2008.03866.x
Geevers S, Mulder W A, Van Der Vegt J J W. 2018. Dispersion properties of explicit finite element methods for wave propagation modelling on tetrahedral meshes. Journal of Scientific Computing, 77(1): 372-396. DOI:10.1007/s10915-018-0709-7
He X J, Yang D H, Wu H. 2014. Numerical dispersion and wave-field simulation of the Runge-Kutta discontinuous Galerkin method. Chinese Journal of Geophysics (in Chinese), 57(3): 906-917. DOI:10.6038/cjg20140320
He X J, Yang D H, Wu H. 2015. A weighted Runge-Kutta discontinuous Galerkin method for wavefield modelling. Geophysical Journal International, 200(3): 1389-1410. DOI:10.1093/gji/ggu487
He X J, Yang D H, Ma X, et al. 2019a. Dispersion-dissipation analysis of the triangle-based discontinuous Galerkin method for scalar wave equation. Geophysical Journal International, 218(2): 1174-1198. DOI:10.1093/gji/ggz170
He X J, Yang D H, Ma X, et al. 2019b. Symplectic interior penalty discontinuous Galerkin method for solving the seismic scalar wave equation. Geophysics, 84(3): T133-T145. DOI:10.1190/geo2018-0492.1
He X J, Yang D H, Ma X. 2020. A weighted Runge-Kutta discontinuous Galerkin method for 3D acoustic and elastic wave-Field modeling. Communications in Computational Physics, 28(1): 372-400. DOI:10.4208/cicp.OA-2018-0072
Hesthaven J S, Warburton T. 2008. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. New York: Springer.
Igel H. 1999. Wave propagation in three-dimensional spherical sections by the Chebyshev spectral method. Geophysical Journal International, 136(3): 559-566. DOI:10.1046/j.1365-246x.1999.00758.x
Igel H. 2017. Computational Seismology: A Practical Introduction. Oxford: Oxford University Press.
Karypis G, Kumar V. 1998. Multilevelk-way partitioning scheme for irregular graphs. Journal of Parallel and Distributed Computing, 48(1): 96-129. DOI:10.1006/jpdc.1997.1404
Käser M, Dumbser M. 2006. An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes-I: The two-dimensional isotropic case with external source terms. Geophysical Journal International, 166(2): 855-877. DOI:10.1111/j.1365-246X.2006.03051.x
Käser M, Dumbser M. 2008. A highly accurate discontinuous Galerkin method for complex interfaces between solids and moving fluids. Geophysics, 73(3): T23-T35. DOI:10.1190/1.2870081
Klin P, Priolo E, Seriani G. 2010. Numerical simulation of seismic wave propagation in realistic 3-D geo-models with a Fourier pseudo-spectral method. Geophysical Journal International, 183(2): 905-922. DOI:10.1111/j.1365-246X.2010.04763.x
Komatitsch D, Vilotte J P. 1998. The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures. Bulletin of the Seismological Society of America, 88(2): 368-392.
Komatitsch D, Tromp J. 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation. Geophysical Journal International, 139(3): 806-822. DOI:10.1046/j.1365-246x.1999.00967.x
Komatitsch D, Vilotte J P, Vai R, et al. 1999. The spectral element method for elastic wave equations-application to 2-D and 3-D seismic problems. International Journal for Numerical Methods in Engineering, 45(9): 1139-1164. DOI:10.1002/(SICI)1097-0207(19990730)45:9<1139::AID-NME617>3.0.CO;2-T
Kristek J, Moczo P. 2003. Seismic-wave propagation in viscoelastic media with material discontinuities: A 3D fourth-order staggered-grid finite-difference modeling. Bulletin of the Seismological Society of America, 93(5): 2273-2280. DOI:10.1785/0120030023
Lähivaara T, Huttunen T. 2010. A non-uniform basis order for the discontinuous Galerkin method of the 3D dissipative wave equation with perfectly matched layer. Journal of Computational Physics, 229(13): 5144-5160. DOI:10.1016/j.jcp.2010.03.030
Lambrecht L, Lamert A, Friederich W, et al. 2017. A nodal discontinuous Galerkin approach to 3-D viscoelastic wave propagation in complex geological media. Geophysical Journal International, 212(3): 1570-1587.
Li B Q. 2006. Discontinuous Finite Elements in Fluid Dynamics and Heat Transfer. New York, NY: Springer Science & Business Media.
Li X F, Lu M W, Liu S L, et al. 2015. A symplectic method for structure-preserving modelling of damped acoustic waves. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2183): 20150105. DOI:10.1098/rspa.2015.0105
Lian X M, Zhang R X. 2013. Numerical simulation of seismic wave equation by local discontinuous Galerkin method. Chinese Journal of Geophysics (in Chinese), 56(10): 3507-3513. DOI:10.6038/cjg20131025
Liu S L, Yang D H, Dong X P, et al. 2017. Element-by-element parallel spectral-element methods for 3-D teleseismic wave modeling. Solid Earth, 8(5): 969-986. DOI:10.5194/se-8-969-2017
Liu Y, Sen M K. 2009. An implicit staggered-grid finite-difference method for seismic modelling. Geophysical Journal International, 179(1): 459-474. DOI:10.1111/j.1365-246X.2009.04305.x
Ma S, Archuleta R J, Liu P C. 2004. Hybrid modeling of elastic P-SV wave motion: A combined finite-element and staggered-grid finite-difference approach. Bulletin of the Seismological Society of America, 94(4): 1557-1563. DOI:10.1785/012003087
Meng W J, Fu L Y. 2018. Numerical dispersion analysis of discontinuous Galerkin method with different basis functions for acoustic and elastic wave equations. Geophysics, 83(3): T87-T101. DOI:10.1190/geo2017-0485.1
Meng X, Shu Q W, Yang Y. 2015. Superconvergence of discontinuous Galerkin methods for time-dependent partial differential equations. Scientia Sinica Mathematica, 45(7): 1041-1060. DOI:10.1360/N012014-00123
Minisini S, Zhebel E, Kononov A, et al. 2013. Local time stepping with the discontinuous Galerkin method for wave propagation in 3D heterogeneous media. Geophysics, 78(3): T67-T77. DOI:10.1190/geo2012-0252.1
Moczo P, Kristek J, Halada L. 2000. 3D fourth-order staggered-grid finite-difference schemes: stability and grid dispersion. Bulletin of the Seismological Society of America, 90(3): 587-603. DOI:10.1785/0119990119
Moczo P, Kristek J, Vavryčuk V, et al. 2002. 3D heterogeneous staggered-grid finite-difference modeling of seismic motion with volume harmonic and arithmetic averaging of elastic moduli and densities. Bulletin of the Seismological Society of America, 92(8): 3042-3066. DOI:10.1785/0120010167
Moczo P, Kristek J, Galis M, et al. 2011. 3-D finite-difference, finite-element, discontinuous-Galerkin and spectral-element schemes analysed for their accuracy with respect to P-wave to S-wave speed ratio. Geophysical Journal International, 187(3): 1645-1667. DOI:10.1111/j.1365-246X.2011.05221.x
Mou Y G, Pei Z L. 2005. Seismic Numerical Modeling for 3-D Complex Media (in Chinese). Beijing: Petroleum Industry Press.
Mu D W, Chen P, Wang L Q. 2013a. Accelerating the discontinuous Galerkin method for seismic wave propagation simulations using multiple GPUs with CUDA and MPI. Earthquake Science, 26(6): 377-393. DOI:10.1007/s11589-013-0047-7
Mu D W, Chen P, Wang L Q. 2013b. Accelerating the discontinuous Galerkin method for seismic wave propagation simulations using the graphic processing unit (GPU)-single-GPU implementation. Computers & Geosciences, 51: 282-292.
Mulder W A, Zhebel E, Minisini S. 2014. Time-stepping stability of continuous and discontinuous finite-element methods for 3-D wave propagation. Geophysical Journal International, 196(2): 1123-1133. DOI:10.1093/gji/ggt446
Niu B H, Sun C Y. 2007. Half-Space Homogeneous Isotropic-Viscoelastic Medium and Seismic Wave Propagation (in Chinese). Beijing: Geological Publishing House.
Pelties C, De La Puente J, Ampuero J P, et al. 2012. Three-dimensional dynamic rupture simulation with a high-order discontinuous Galerkin method on unstructured tetrahedral meshes. Journal of Geophysical Research: Solid Earth, 117(B2): B02309. DOI:10.1029/2011JB008857
Pelz R B. 1991. The parallel Fourier pseudospectral method. Journal of Computational Physics, 92(2): 296-312. DOI:10.1016/0021-9991(91)90212-4
Reed W H, Hill T R. 1973. Triangular mesh methods for the neutron transport equation. Los Alamos Scientific Laboratory Report. LA-UR, 473-479. Los Alamos Scientific Laboratory, N. Mex, USA.
Rivière B. 2008. Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation.//Frontiers in Applied Mathematics. Houston, Texas: Rice University.
Riviére B, Wheeler M F. 2003. Discontinuous finite element methods for acoustic and elastic wave problems. Contemporary Mathematics, 329: 271-282. DOI:10.1090/conm/329
Shragge J, Tapley B. 2017. Solving the tensorial 3D acoustic wave equation: a mimetic finite-difference time-domain approach. Geophysics, 82(4): T183-T196. DOI:10.1190/geo2016-0691.1
Tromp J, Komatitsch D, Liu Q Y. 2008. Spectral-element and adjoint methods in seismology. Communications in Computational Physics, 3(1): 1-32.
Versteeg R J, Grau G. 1991. The Marmousi Experience.//61st Ann. Internat Mtg., Soc. Expi. Geophys.. Expanded Abstracts.
Wang J, Yang D H, Jing H, et al. 2019a. Full waveform inversion based on the ensemble Kalman filter method using uniform sampling without replacement. Science Bulletin, 64(5): 321-330. DOI:10.1016/j.scib.2019.01.021
Wang N, Li J H, Borisov D, et al. 2019b. Modeling three-dimensional wave propagation in anelastic models with surface topography by the optimal strong stability preserving Runge-Kutta method. Journal of Geophysical Research: Solid Earth, 124(1): 890-907. DOI:10.1029/2018JB016175
Wang W S, Zhang H, Li X F. 2013. Review on application of the discontinuous Galerkin method for modeling of the seismic wavefield. Progress in Geophysics (in Chinese), 28(1): 171-179.
Wilcox L C, Stadler G, Burstedde C, et al. 2010. A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media. Journal of Computational Physics, 229(24): 9373-9396. DOI:10.1016/j.jcp.2010.09.008
Xue Z, Dong L G, Li X B, et al. 2014. Discontinuous Galerkin finite-element method for elastic wave modeling including surface topography. Chinese Journal of Geophysics (in Chinese), 57(4): 1209-1223. DOI:10.6038/cjg20140418
Yang D, Song G, Lu M. 2007. Optimally accurate nearly analytic discrete scheme for wave-field simulation in 3D anisotropic media. Bulletin of the Seismological Society of America, 97(5): 1557-1569. DOI:10.1785/0120060209
Yang D H, Wang L. 2010. A split-step algorithm for effectively suppressing the numerical dispersion for 3D seismic propagation modeling. Bulletin of the Seismological Society of America, 100(4): 1470-1484. DOI:10.1785/0120090200
Yang D H, He X J, Ma X, et al. 2016. An optimal nearly analytic discrete-weighted Runge-Kutta discontinuous Galerkin hybrid method for acoustic wavefield modeling. Geophysics, 81(5): T251-T263. DOI:10.1190/geo2015-0686.1
Ye R C, De Hoop M V, Petrovitch C L, et al. 2016. A discontinuous Galerkin method with a modified penalty flux for the propagation and scattering of acousto-elastic waves. Geophysical Journal International, 205(2): 1267-1289. DOI:10.1093/gji/ggw070
Zhang J B, Yang D H, He X J, et al. 2018. Discontinuous Galerkin method for solving wave equations in two-phase and viscoelastic media. Chinese Journal of Geophysics (in Chinese), 61(3): 926-937. DOI:10.6038/cjg2018L0095
Zhang J F, Liu T L. 2002. Elastic wave modelling in 3-D heterogeneous media: 3-D grid method. Geophysical Journal International, 150(3): 780-799. DOI:10.1046/j.1365-246X.2002.01743.x
Zhang J F, Gao H W. 2009. Elastic wave modelling in 3-D fractured media: an explicit approach. Geophysical Journal International, 177(3): 1233-1241. DOI:10.1111/j.1365-246X.2009.04151.x
Zhang J H, Wang W M, Zhao L F, et al. 2007. Modeling 3-D scalar waves using the Fourier finite-difference method. Chinese Journal of Geophysics (in Chinese), 50(6): 1854-1862.
Zhang W, Zhang Z G, Chen X F. 2012. Three-dimensional elastic wave numerical modelling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids. Geophysical Journal International, 190(1): 358-378. DOI:10.1111/j.1365-246X.2012.05472.x
Zhang Y J, Gao J H, Han W M, et al. 2019. A discontinuous Galerkin method for seismic wave propagation in coupled elastic and poroelastic media. Geophysical Prospecting, 67(5): 1392-1403. DOI:10.1111/1365-2478.12781
董良国, 马在田, 曹景忠. 2000. 一阶弹性波方程交错网格高阶差分解法稳定性研究. 地球物理学报, 43(6): 856-864. DOI:10.3321/j.issn:0001-5733.2000.06.015
贺茜君, 杨顶辉, 吴昊. 2014. 间断有限元方法的数值频散分析及其波场模拟. 地球物理学报, 57(3): 906-917. DOI:10.6038/cjg20140320
李庆扬, 王能超, 易大义. 2008. 数值分析. 5版. 北京: 清华大学出版社.
廉西猛, 张睿璇. 2013. 地震波动方程的局部间断有限元方法数值模拟. 地球物理学报, 56(10): 3507-3513. DOI:10.6038/cjg20131025
孟雄, 舒其望, 杨扬. 2015. 发展型偏微分方程间断有限元方法的超收敛性献给林群教授80华诞. 中国科学: 数学, 45(7): 1041-1060.
牟永光, 裴正林. 2005. 三维复杂介质地震数值模拟. 北京: 石油工业出版社.
牛滨华, 孙春岩. 2007. 半无限空间各向同性黏弹性介质与地震波传播. 北京: 地质出版社.
汪文帅, 张怀, 李小凡. 2013. 间断的Galerkin方法在地震波场数值模拟中的应用概述. 地球物理学进展, 28(1): 171-179.
薛昭, 董良国, 李晓波, 等. 2014. 起伏地表弹性波传播的间断Galerkin有限元数值模拟方法. 地球物理学报, 57(4): 1209-1223. DOI:10.6038/cjg20140418
张金波, 杨顶辉, 贺茜君, 等. 2018. 求解双相和黏弹性介质波传播方程的间断有限元方法及其波场模拟. 地球物理学报, 61(3): 926-937. DOI:10.6038/cjg2018L0095
张金海, 王卫民, 赵连锋, 等. 2007. 傅里叶有限差分法三维波动方程正演模拟. 地球物理学报, 50(6): 1854-1862. DOI:10.3321/j.issn:0001-5733.2007.06.028
张铁. 2015. 间断有限元理论与方法. 北京: 科学出版社.