基于非结构网格求解三维D′Alembert介质中声波方程的并行加权Runge-Kutta间断有限元方法

贺茜君, 杨顶辉, 仇楚钧, 周艳杰, 常芸凡

1 北京工商大学数学与统计学院, 北京 100048 2 清华大学数学科学系, 北京 100084

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ä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={v∈L1(Ω):v|Ωi∈Pκ(Ωi)}，其中Pκ(Ωi)表示区域Ωi上次数不超过κ次的多项式.从定义中可以看出，测试函数v在Ωi以外的区域上不定义或者令其为0，因此，它不必在整个区域上连续，可以在子区域{Ωi}之间有间断，这也是间断有限元与传统有限元最大的区别.

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

(2)

即：

(3)

(4)

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

(5)

其中F(u)表示通量，B表示耗散项，当B=0时表示无耗散，方程(5)退化为普通的声波方程.需要指出的是，本文发展的方法也适用于一阶速度-应力方程.在(5)式两边同乘以试验函数v，并在子区域Ω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具有如下的基函数展开形式：

(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)的流程计算所有的体质量矩阵和刚度矩阵，以及面与面之间关联的质量矩阵.具体的过程可参考附录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)中所示的规则六面体剖分的基础上，每个六面体中有图(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情形下的数值稳定性、数值频散和耗散主要基于Von Neumann分析.为此，假设一个简谐波传播经过3D区域，并假定此区域是均匀、无界区域，且方程右端为无源项.我们主要考察四面体网格剖分，在如图2a中所示的规则六面体剖分的基础上(Mulder et al., 2014; Ferroni et al., 2017; He et al., 2020)，在一个立方体单元内有6个四面体单元，如图2b所示，其中每个四面体及其面的顺序参考Mulder 等(2014)编号.我们假设如下平面简谐波在所考虑的区域内传播：

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

+kcosθlh-iωt)],

(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π]都成立.这在数学上可以用下面一个优化问题来表示：

maxα

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

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

α≥0.

(19)

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

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

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

(22)

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

(23)

从而可得

(24)

对于带耗散的系统：

(25)

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

λ2≈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的对比

另外，从方程(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)

其中α=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)分别表示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=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中发现，在δ较小时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介质的衰减特性.

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

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

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

图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)

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

(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

图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

接下来，为了考察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.

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.震源函数是：

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方法得到的数值解.从图中可以看出，P1方法出现少许数值频散，而P2方法与解析解吻合较好，这说明了提高算法精度有助于降低数值频散.图10给出了P2方法的波场快照图，此时波已经传至边界.波场快照图中无明显可见数值频散.

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

图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 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)对应于耗散参数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，震源函数如方程(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=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介质中的声波方程，空间离散采用了基于数值通量的间断有限元公式，时间离散基于对角隐式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的坐标分别为 (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的体积，且容易得到如下偏导数的值(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)