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近似，所以我们首先需选取试验函数空间.我们使用的试验函数空间是多项式空间V_h={v∈L¹(Ω): v|_Ωi∈P^κ(Ω_i)}，其中P^κ(Ω_i)表示区域Ω_i上次数不超过κ次的多项式.从定义中可以看出，测试函数v在Ω_i以外的区域上不定义或者令其为0，因此，它不必在整个区域上连续，可以在子区域Ω_i之间有间断，这也是间断有限元与传统有限元最大的区别.

为了将方程(1)改写成一阶双曲系统的形式，我们引入三个变量p、q和s，且p、q、s满足条件p_t=u_x, q_t=u_y, s_t=u_z，其中u_x、u_y和u_z分别表示位移u对空间变量x、y和z的偏导数.则方程(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=(n₁, n₂, n₃)^T表示边界∂Ω_i处的单位外法向量.由于u在∂Ω_i处有可能间断，所以(6)式中的F(u)·n没有定义，我们将F(u)·n替换为数值通量：F(u)·n|_∂Ωi=, n)，其中u^int表示u在区域Ω_i内部逼近∂Ω_i处的值，u^ext表示u在区域Ω_i外部逼近∂Ω_i处的值.这里我们采用一种常用通量，即Local Lax-Friedrichs(LLF)通量(Cockburn and Shu, 1989, 2001)，其表达式为

(7)

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

(8)

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

(9)

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

(10)

其中，

(12)

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

(12)

其中w_lⁱ(x)表示Ω_i内的基函数，C_lⁱ(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.2 时间离散

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

(14)

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

(15)

其中η∈[0, 1]是加权系数.在下文中若无特殊说明，我们均取η=0.5进行计算.K₁⁽ⁿ⁾和K₂⁽ⁿ⁾可以用如下公式计算：

(16)

其中.得到K⁽ⁿ⁾的近似值之后，令T⁽ⁿ⁾=C⁽ⁿ⁾+(1－2r)△tK⁽ⁿ⁾，按如下公式求得K₁⁽ⁿ⁾和K₂⁽ⁿ⁾：

(17)

将如上所得的K₁⁽ⁿ⁾、K₂⁽ⁿ⁾、K₁⁽ⁿ⁾和K₂⁽ⁿ⁾代入方程(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)

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

2.1 数值稳定性分析

为了保持数值算法的稳定性，我们引入库朗数条件，也称Courant-Friedrichs-Lewy(简称CFL)条件：α=c△t/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：对于P¹阶WRKDG方法，c△t/h≤α_max≈0.244；对于P²阶WRKDG方法，c△t/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)

其中，算子L₁对应无耗散情形的声波传播，算子L₂对应系统中的耗散项，实际上，L₂只与耗散参数r有关.我们将算子L、L₁和L₂对应的谱半径分别记作λ、λ₁和λ₂.则由(21)，我们可得λ≤λ₁+λ₂.对于双曲型系统：

(22)

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

(23)

从而可得

(24)

对于带耗散的系统：

(25)

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

(26)

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

(27)

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

(28)

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

(29)

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

2.2 数值频散、耗散分析

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

(30)

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

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

另外，我们从图 3b与3d中发现，在δ较小时D′ Alembert介质存在比较大的频散，随着δ的增加，R的值接近1，也就是说，随着δ的增加，频散变小.对比r=0与r=10这两种情形，我们可以得出下述结论：在δ(波数)较小时，D′ Alembert介质中存在主要由耗散参数r引起的频散；随着δ(波数)的增加，耗散参数引起的频散变小，而数值算法引起的数值频散占主导.

图 4显示了振幅在一个时间步Δt内的耗散S=e^－ω_i△t随采样率的变化图.图中参数的选取和上文中数值频散分析中相同.当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介质中理论耗散因子e^－r△t/2，体现了D′Alembert介质的衰减特性.

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中不同颜色代表分给不同进程.

3.3 并行流程

在本研究中，我们使用Message Passing Interface(MPI)编程策略对程序进行并行化处理，整个流程可参考图 6.在程序开始阶段，由于我们采用无积分(quadrature-free)的DGM，因此可将方程(13)中所需的参考单元内的质量、刚度矩阵以及四面体面积分矩阵提前计算出来，然后导入到程序中.同时，我们也将生成的网格文件导入，并利用Metis进行网格划分，同时需要计算出网格单元连接矩阵，以备后用.在用Metis进行网格划分完毕后，对于每一个进程我们需要记录三类网格及其信息.以图 7来进行说明，第一类网格是每个进程的内网格，例如，以第一个进程为例，图 7中蓝色区域内的网格即是第一个进程的内网格；第二类网格是该进程的辅助网格，这类网格属于其他进程的内网格，位于进程1所有网格的边界处，如图 7a中的绿色网格部分(图 7b是图 7a的一个更清晰的展示)，作用是接收来自于其他进程更新后的变量信息；第三类网格是属于该进程内，但是作为其他进程的辅助网格，这类网格在消息传递时需要由本进程传递给其他进程使用，如图 7c中红色网格显示.

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

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

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

(31)

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

(32)

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

(33)

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

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

5 波场模拟

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

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

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

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)中相同，其中f₀=60 Hz，主频约为25 Hz.图 14给出了在T=0.3 s时刻下的波场快照图，图 14a对应于无耗散情形，而图 14b对应于耗散参数r=4.从图中可以看出，图 14b较图 14a暗一些，证明了D′Alembert介质中的衰减效应.

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)所示，其中f₀=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介质中的声波波场.

6 结论

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

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

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

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

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

(A1)

其中，角标l和m表示矩阵的l行m列元素，上角标i表示所考虑的单元为Ω_i.在计算之前，我们首先将一般的四面体单元Ω_i变换到如图 1所示的参考单元E内.如图 1所示，假设原单元四个顶点1、2、3及4的坐标分别为(x₁, y₁, z₁)、(x₂, y₂, z₂)、(x₃, y₃, z₃)和(x₄, y₄, z₄)，变换到参考单元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的体积，且容易得到如下偏导数的值(Dumbser and KäserKäser, 2006)：

(A3)

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

(A4)

例如，要计算原质量矩阵M₁，可以利用关系式，M₁=|J|M′₁.

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

(A5)

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

(A6)

至此，我们已将方程(13)中所有积分矩阵的表达式阐述完毕.对于四面体的四个面F₁、F₂、F₃和F₄的外法向量n₁, n₂, n₃和n₄的计算，有如下公式，

(A7)

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