地球物理学报  2019, Vol. 62 Issue (7): 2660-2670   PDF    
基于多重网格预条件求解平均导数法离散的Helmholtz方程
袁雨欣1,2,3,4, 李阿满1,2,3,4, 胡婷1,2,3,4, 郭鹏5, 刘洪1,2,3,4     
1. 中国科学院地质与地球物理研究所, 北京 100029;
2. 中国科学院地球科学研究院, 北京 100029;
3. 中国科学院油气资源研究院重点实验室, 北京 100029;
4. 中国科学院大学, 北京 100049;
5. 天津大学水利工程仿真与安全国家重点实验室, 天津 300072
摘要:有限差分法求解Helmholtz方程,依赖于两点:1差分格式的构造;2高效的求解算法.本文采用平均导数法离散Helmholtz方程.该差分格式有三点好处:1能适用于横纵不等间距采样;2在完全匹配层区域(PML),差分方程与微分方程逐点相容;3能将一个波长内的采样点数减少至少于4.求解离散的Helmholtz方程的算法一般分为直接法和迭代算法.直接法由于内存需求太大而无法适用于大规模问题;基于Krylov子空间的迭代方法结合多重网格预条件算法是一种快速高效求解方法,然而对于横纵不等间距采样(在多重网格中称为各向异性问题),经典的多重网格方法失效.本文分析了经典多重网格的三个重要组成部分:完全加权限制算子,点松弛技术以及双线性延拓算子,进而采用了半粗化技术代替全粗化技术,线松弛技术代替点松弛技术以及依赖差分算子的延拓算子代替双线性延拓算子,使得各向异性问题变得收敛;而且对于非均匀介质中-低频率的迭代问题,我们获得了较为满意的收敛速度.
关键词: 多重网格      各向异性问题      Helmholtz方程      平均导数法     
A multigrid-based preconditioner for solving the Helmholtz equation with average-derivative optimal scheme
YUAN YuXin1,2,3,4, LI AMan1,2,3,4, HU Ting1,2,3,4, GUO Peng5, LIU Hong1,2,3,4     
1. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
2. Institutes of Earth Science, Chinese Academy of Sciences, Beijing 100029, China;
3. Key Laboratory of Petroleum Resources Research, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
4. University of Chinese Academy of Sciences, Beijing 100049, China;
5. State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University, Tianjin 300072, China
Abstract: An efficient finite-difference method for solving Helmholtz equation depends on two points:one is discrete scheme, the other is efficient algorithm. In this paper, we adopt the average-derivative scheme, which owns three advantages:Firstly, it can be applied to unequal directional sampling intervals for Helmholtz equation. Secondly, the scheme is pointwise consistent with Helmholtz equation in a perfect matched layer. And thirdly, it requires less than 4 grid points sampling per wavelength. To solve the discrete Helmholtz equation, which is extremely large and indefinite, direct methods cannot resolve well, and the Krylov subspace iterative methods, such as Bi-CGSTAB and GMRES combining a multigrid-based preconditioner, are good choices. However, the standard multigrid algorithm fails to converge when it encounters unequal directional sampling intervals, which is called anisotropy in multigrid. We analyze the most important three parts of standard multigrid:full weighting restriction operator, point relaxation methods and bilinear interpolation operator, and then we replace them with semi-coarsening, line relaxation and operator-dependent interpolation to make it convergent in anisotropic problems. Consequently, we obtain a satisfactory convergence speed for low and moderate frequency iterative problems in heterogeneous media.
Keywords: Multigrid    Anisotropic problem    Helmholtz equation    Average-derivative optimal scheme    
0 引言

Helmholtz方程是一类非常重要的椭圆型偏微分方程,它可以用来描述波传播和散射现象,因此在电磁辐射、声学、光学以及地震学方面具有重要的研究价值.

为了在有界区域数值求解Helmholtz方程,我们引入PML(Bérenger, 1994; Tsynkov and Turkel, 2001;Chen et al., 2011, 2013).可以采用有限差分方法离散该方程:Jo等(1996)提出旋转九点差分格式,相比经典的五点差分格式,它在保持同样计算精度的前提下,能减少一个波长内的采样点数,从而减少计算量;由于Jo的方法不适用于横纵不等间距网格,Chen(2012)提出了平均导数法,该方法能适应于任意横纵网格间距; 与此同时,Chen等(2013)证明了Jo的方法不适用于PML,并且指出若采用平均导数法,差分方程在PML区域与微分方程逐点相容.

采用平均导数法离散该方程之后,我们得到了一个稀疏的线性系统,由于大波数的影响,这个线性系统是不定的,并且不具有厄密性.求解线性方程组的方法,一般可以分为两类:直接法(LU分解)和迭代法.二者各有优缺点:LU分解得到的结果可以适用于多炮并行计算,然而受到计算机内存和计算量的限制,LU分解不适用于求解大规模问题;基于Krylov子空间的迭代方法,例如Bi-CGSTAB方法和GMRES方法可以适用于大规模问题,但是如果线性系统没有好的谱分布,迭代方法的收敛速度很慢甚至停滞.因此我们需要一个好的预条件算子来改善线性系统的谱分布,许多学者为此做了很多工作:Bayliss等(1983)提出拉普拉斯预条件算子,结合迭代法CGNR高效地求解了Helmholtz方程;Laird等(2002)在拉普拉斯预条件算子的基础上增加了一个实数项,构造出来了实移位拉普拉斯预条件算子;Erlangga等(2004)总结了以上两种预条件算子,提出复移位拉普拉斯预条件算子,Van Gijzen等(2007)分析了此类预条件后的算子的谱分布,给出了如何选择最优的“移位”值,并且结合Krylov子空间迭代方法加速了Helmholtz方程的收敛;Erlangga等(2006b)比较了不完全LU分解预条件方法和复位移拉普拉斯算子预条件方法,并得到结论:用多重网格方法得到复移位拉普拉斯算子的近似逆后,结合Bi-CGSTAB迭代法,是一种快速且稳健的方法.

对于预条件的Krylov子空间方法,至关重要的一步就是如何求取预条件算子的逆.获得逆算子的显式形式是没有必要的,我们在尽量减少计算量的前提下,只需求取一个近似逆.因此多重网格方法是一个很好的选择.对于横纵网格等间距的情形,多重网格中又称为各向同性问题,反之称为各向异性问题(Briggs et al., 2000).

对于各向同性问题,Erlangga等(2006a, 2008)提出V(1, 1)循环的多重网格法去近似求解预条件算子的逆,这种多重网格方法改进了De Zeeuw(1990)提出的依赖差分算子的延拓算子,提出了一种适用于Helmholtz方程的延拓算子;Riyanti等(2007)提出将Erlangga的方法与Bi-CGSTAB方法结合,取得了很好的收敛速度;Riyanti等(2007)Plessix(2007)将上述方法应用于三维Helmholtz方程求解,证明了迭代法对于大规模问题的有效性.Chen等(2011)针对Helmholtz方程给出了一种不同于De Zeeuw(1990)Erlangga等(2006a)的全新的依赖于差分算子的延拓算子,取得了比经典多重网格方法收敛速度更快的效果.

对于各向异性问题,曹健等(2015)Cao等(2018)指出采用经典的多重网格方法结合Bi-CGSTAB方法,并不能获得很好的收敛速度,甚至完全不收敛;本文分析了经典多重网格的三个重要组成成分:延拓算子、光滑算子、限制算子.进一步提出了适用于各向异性问题的多重网格算法:半粗化技术、线迭代技术以及依赖差分算子的延拓算子.本文的两个数值算例说明:经典多重网格方法不适用于各向异性问题;本文的算法对于中低频率的各向异性问题,收敛速度很快.

1 Helmholtz方程的离散以及谱分析 1.1 Helmholtz问题的数学定义

求解在二维声介质中波的传播问题,相应的Helmholtz方程为:

(1)

其中,u代表单频压力场,是拉普拉斯算子,是虚数单位.α是一个实数,表征介质的阻尼性质,通常0≤α < < 1,本文设定为0.k是波数,定义为k:=2πf/v(x),f是频率,v(x)是介质速度.g(x, ω)是源项.

为了在有界区域数值求解Helmholtz问题,我们引入PML:

(2)

其中,是角频率,σx是仅与x坐标有关的函数,定义为:

(3)

σy也可以类似定义.其中f0是雷克子波主频.如图 1所示:LPML是PML区域的宽度.lx是当前点距离∂Ω0的距离,a0是经验常数,一般取为a0=1.79(Chen et al., 2011).

图 1 完全匹配层模型,Ω0是有界计算区域,Ω1Ω2Ω3是PML区域 Fig. 1 The perfectly matched layer model.Ω0 is the bounded physical domain, Ω1Ω2Ω3 is perfectly matched domain

采用平均导数法离散Helmholtz方程(Chen, 2012),从而得到线性系统:

(4)

上述线性系统,可以通过矩阵A表示;另外一种较为简洁的表示是差分模板:在(xm, yn)点,差分模板记为:

(5)

写出差分模板的各个元素:

(6)

其中,hx, hy分别是xy方向的采样步长;α, β, c, d, e是通过最小化相速度误差得到的优化权重系数(Chen, 2012);分别定义为:

(7)

其中:

(8)

从差分模板可以得到阻抗矩阵的一些性质,例如,若要求矩阵对称,需满足

(9)

若忽略PML,阻抗矩阵是对称的实矩阵,根据线性代数的理论,它的特征值全部是实数;由于PML的存在,阻抗矩阵是非对称复矩阵.

1.2 预条件

阻抗矩阵A具有A=L-z1D的形式,z1=1-αi(本文α取为0),LD分别是离散的Laplacian算子和零阶项exeyk2.由于PML的层数相比内部区域的层数小得多,所以A是一个复数域的、带状分布的、稀疏、近对称矩阵.

为了求解上述线性系统,我们采用Krylov子空间方法,特别是Bi-CGSTAB方法.但是该方法要求线性系统有较好的谱分布,从而保证问题快速收敛.对于Helmholtz方程而言,预条件算子可以是(Erlangga et al., 2006a;Chen et al., 2011):

(10)

至此,我们得到了预条件的Helmholtz微分方程:

(11)

同样采用平均导数法在Ω上离散,从而得到预条件矩阵M以及预条件线性系统

(12)

其中.若(β1, β2)=(0, 0),则它是拉普拉斯预条件算子(Bayliss et al., 1983);Laird和Erlangga提到的预条件算子对应于.如果β1>0, β2>0,则它是复移位拉普拉斯算子.本文选取β1=1, β2=0.5.

1.3 谱分析

首先,忽略PML,考虑零值Dirichlet边界条件,形如的差分模板离散得到的系数矩阵的特征值为:

(13)

Nx, Ny分别为x, y方向的网格点数,l=1, 2,…,Ny, m=1, 2,…,Nx.对于本文的线性系统而言:

(14)

实际数值模拟中,我们需要考虑PML,此时,阻抗矩阵是近对称复矩阵,特征值没有解析表达式.

我们采用如下的试验计算两种情况下矩阵的谱分布:设定参数:Nx=Ny=70, 均匀介质v=3000 m·s-1f=30Hz(或者设定波数k=0.0628),定义横纵采样间距比 .图 2所示为由(13)式计算得到的矩阵的谱分布(未考虑PML,取γ=2),图 3所示为数值计算(考虑PML,每边层数设定为5,取γ=4)的矩阵的谱分布.

图 2 设定Nx=Ny=71, k=0.0628,hy=5 m,γ=2,下列矩阵的谱分布 (a) A; (b) AM-1; (c) D-1L; (d)是(c)中的特征值作分式线性变换得到. Fig. 2 Spectral pictures of the following matrices for Nx=Ny=71, k=0.0628, hy=5 m, γ=2 (a) A; (b) AM-1; (c) D-1L; (d) Fractional linear transform of characteristic value in (c).
图 3 PML情况下,图 2中各个矩阵的谱分布 Fig. 3 Same as Fig. 2 but for PML case

图 2a可以看到,A是实对称矩阵,特征值全部是实数,而且特征值在(0, 0)多有聚集,条件数较大.从图 2b可以看到,用复移位拉普拉斯算子预条件后的阻抗矩阵AM-1,聚集在(0, 0)处的特征值被“推”至复平面,从而改善了系统的特征值分布.可以证明(Chen et al., 2011),如果μD-1L的特征值,当且仅当AM-1的特征值,因此图 2b图 2d等价;并且这个分式线性变换(共形变换)能将实轴变换成一个以为圆心,R= 为半径的圆,虚部小于0的点被映射至此圆的内部,反之被映射至此圆的外部.可以从图 2c图 2d看出,AM-1的特征值全部位于圆心C=(0.5, 0), 半径R=0.5的圆上.

图 3a可以看到,由于PML的影响,A的特征值不再只有实数,而且矩阵的条件数仍然较大.图 3bAM-1的谱.图 3cD-1L的特征值分布,图 3d是将图 3c中的特征值作分式线性变换得到,它与图 3b等价.可以看出,由于D-1L的所有特征值虚部全部小于0,因此经过分式线性变换之后,这些特征值全部位于圆心C=(0.5, 0), 半径R=0.5的圆的内部,如图 3d所示.换言之,我们证明了AM-1的特征值全部位于圆心C=(0.5, 0), 半径R=0.5的圆的内部.可以看到,复移位拉普拉斯预条件算子对于改善Helmholtz方程的谱分布是很有效果的.

2 多重网格求解M-1

为了行文方便,我们给出一些记号,如表 1所示.

表 1 多重网格的参数记号 Table 1 Parameter notations of multigrid
2.1 多重网格的思想

我们可以通过研究最简单的两网格法来理解多重网格的思想,假设待求解的线性方程组为:

(15)

表示某近似解,ul表示精确解.则的修正值为

(16)

残差为

(17)

残差与修正值之间满足残量方程:

(18)

现在需要用Ml的一个近似来求解el.

经典迭代方法,例如Jacobi迭代,Gauss-Seidel迭代是通过寻找下述方程的一个近似解来完成的:

(19)

是一个"较简单"的算子. 例如, 在Jacobi迭代中, Ml的对角线部分; 在Gauss-Seidel迭代中, Ml的下三角部分.第二个近似是根据(20)式得到的:

(20)

与经典迭代方法不同,多重网格是寻找“粗糙化”的算子Ml-1而不是“简单化”的算子.在网格大小为hl-1, x, hl-1, y的一个更粗糙的网格上(常取hl-1, x=2hl, x, hl-1, y=2hl, y,即经典的全粗化方法;也有特殊的取法,如本文采用的半粗化:hl-1, x=hl, x, hl-1, y=2hl, y),构造Ml的某个适当近似算子Ml-1.则残量方程近似为:

(21)

由于Ml-1的维数较小,该方程比原方程更易求解.为了定义粗网格上的亏损rl-1,需要一个限制算子Rl-1,将细网格上的亏损rl限制在粗网格上

(22)

一旦得到(22)式的一个解后,还需要一个延拓算子Pl将修正值内插到细网格Ωl

(23)

最后ul近似解可以校正为

(24)

总而言之,多重网格的逻辑是:

(1) 误差可划分为高频摆动分量和低频光滑分量;

(2) 通过简单的定常迭代去消除高频摆动分量;

(3) 在粗网格上消除那些顽固的低频光滑分量.

2.2 各向异性问题中的光滑、限制以及延拓算子

最流行,而且也是最应该试验的光滑方法是Gauss-Seidel方法,因为它能产生很好的收敛速度.然而Gauss-Seidel方法的精确形式依赖于网格点所选取的排序.对于典型的等间距差分格式,最好使用红-黑排序;当沿某一维比沿另一维的耦合程度更强时,例如本文γ=2, 4的情形,沿y维比沿x维的耦合程度更强,我们可以考虑沿着y维同时将整条线松弛:

(25)

我们称之为y-line Gauss-Seidel.最邻近耦合的线松弛涉及一个三对角系统的求解,因而仍然具有很高的效率;如果选择红黑斑马线迭代,则根据m的奇偶性,对网格点涂色分类:当m为奇数,则为红色;反之为黑色.迭代可以分为两步,首先对红网格点更新:

(26)

再对黑网格点更新:

(27)

红与黑两部分迭代是分开进行的,因此适合并行.

延拓和限制算子的一个简洁的表示方法是,给出它们的符号.选择粗网格Ωk-1上的脉冲函数uk-1=δk-1Uk-1, 它在某个(x, y)∈Ωk-1上取1,其他点取0.再定义k-1,从而得到符号P.最流行的拓展算子是简单的双线性插值算子,这种算子适用于各向同性问题,它给出Ωk上九个点(x, y), (x+hk, x, y), …, (xhk, x, yhk, y)上的非零值,分别为1,1/2,…,1/4.因此符号为:

对于本文的各向异性问题,我们选择的半粗化延拓算子的符号为:

关于限制算子的选择,我们通常可以选为延拓算子的伴随算子.根据算子的符号记法,限制算子和延拓算子的符号元素之间的关系为:

(28)

因此,完全加权限制算子,就是双线性插值算子的伴随.对于非均匀不连续介质问题,我们可以选择依赖差分算子的延拓算子(De Zeeuw, 1990Erlangga et al., 2006a).本文采用Chen(2012)提出的适用于Helmholtz的延拓算子.

2.3 本文的算法框架

针对γ=2的情形,我们给出多重网格V循环求解Mlul=fl的递归算法, 如表 2所示.

表 2 针对各向异性问题的多重网格算法 Table 2 Multigrid algorithm for anisotropic problem

这里,粗网格Ωk-1上的阻抗矩阵Mk-1,通常可以通过两种方式获得:

(ⅰ)类似于Ml,可以通过平均导数法,在Ωk-1上离散算子得到.

(ⅱ) Galerkin方法:Mk-1=Rk-1MkPk.

本文选用第二种方法来获取Mk-1,它的好处在于,我们无须考虑边界条件.

3 数值算例

本节选取了两个模型,来说明本文算法的有效性.我们采用V(1, 1)多重网格,层数设定为9层.PML层数设定为每一个维度的网格点数的1%,震源设置在PML吸收层以下第二层中间的位置.收敛的准则是相对残差的二范数达到10-6,如果迭代次数超过10000次仍未达到此准则,我们认定它不收敛.所有测试均在Intel 8核,内存16 GB的机器上,由Fortran 77程序完成.

3.1 均匀模型

首先,考虑一个速度大小为3000 m·s-1的均匀模型,横纵采样点数均为401;为了测试本文算法对于各向异性问题的有效性,我们选定γ=1, 2, 4,从低到高选取了6个频率,分别是f=10, 20, 30, 40, 50, 60 Hz.本文的方法我们记为Ⅰ,经典的多重网格算法我们记为Ⅱ,测试结果如表 3图 4所示,表中数据为Bi-CGSTAB的迭代次数以及CPU耗时.

表 3 γ=1, 2, 4,均匀模型的Bi-CGSTAB的迭代次数以及CPU耗时(圆括号内) Table 3 Number of Bi-CGSTAB iteration and CPU time (in parentheses) for γ=1, 2, 4
图 4 γ=1, 2, 4,均匀模型的Bi-CGSTAB的迭代次数(a)以及CPU耗时(b)与频率的关系 Fig. 4 Relationship between the number of Bi-CGSTAB iteration (a) and CPU time (b) with frequency for γ=1, 2, 4

从测试结果可以看出,对于各向同性问题,本文算法和经典多重网格算法一致,收敛速度较快;对于弱各向异性问题(γ=2),经典多重网格仅在低频段缓慢收敛,中频段不收敛,而本文算法仍然能快速收敛;对于强各向异性问题(γ=4),经典多重网格已经失效,而本文算法在中低频段仍然能有令人满意的收敛速度.

3.2 部分Marmousi模型

本节选用部分Marmousi模型,来说明本文算法对于强不均匀模型的有效性;我们考虑的模型横纵采样点数为301×301,速度模型如图 5所示,选定γ=1, 2, 4,测试了两个频率f=10, 30 Hz,结果如图 6所示.

图 5 部分Marmousi纵波速度模型 Fig. 5 Part of the Marmousi velocity model
图 6 部分Marmousi模型在10 Hz和30 Hz的收敛情况 (a) γ=1, 2, 4, f=10 Hz, Δz=4 m相对残差随着迭代次数的变化情况;(b)γ=1, 2, 4, f=30 Hz, Δz=4 m相对残差随着迭代次数的变化情况;(c),(e),(g)是f=10 Hz下,γ=1, 2, 4, Helmholtz方程解的实部;(d),(f),(h)是f=30 Hz下,γ=1, 2, 4, Helmholtz方程解的实部. Fig. 6 Convergence situations for the 10 Hz and 30 Hz with part of Marmousi model (a) γ=1, 2, 4, f=10 Hz, Δz=4 m; (b) γ=1, 2, 4, f=30 Hz, Δz=4 m; (c), (e), (g) are real part of the mono-frequency wavefield at 10 Hz for γ=1, 2, 4 respectively; (d), (f), (h) are real part of the mono-frequency wavefield at 30 Hz for γ=1, 2, 4 respectively.

经典的多重网格方法对于本例已经失效,从测试结果可以看到,本文的算法对于强不均匀介质的中低频段,收敛速度仍然很快;随着频率的提高,由于零特征值的影响,收敛速度变慢.

3.3 内存消耗和耗时对比

本小节测试了本文方法的内存需求以及计算速度,采用的模型是3.1节的均匀模型和3.2节的部分Marmousi模型,并将计算结果与直接分解方法作了对比,结果如表 4所示.

表 4 直接法和本文迭代法内存消耗情况以及CPU耗时对比 Table 4 Comparison of the memory occupation and CPU time between direct solver and iterative method

表中,MLUTLU分别代表LU分解法的内存消耗以及CPU耗时,MiterTiter分别代表本文的迭代法内存消耗以及耗时.迭代法内存消耗远小于直接法,并且在弱各向异性的情况下,中低频率的迭代问题,迭代法相比直接法也明显更省时.

4 结论与建议

Helmholtz方程的高效数值求解,不仅依赖于离散方法,更与算法的选择密切相关.本文的研究重点是:横纵不等间距采样情况下,Helmholtz方程的高效迭代解法.本文采用平均导数法离散Helmholtz方程;进一步,我们分析了原有线性系统的谱,发现如果不做预条件处理,条件数太大,特征值在0处多有聚集,导致经典的迭代方法失效.复移位拉普拉斯算子是一种高效的针对Helmholtz方程的预条件算子,根据谱分析结果,我们知道,预条件后的线性系统的特征值将全部落在圆心C=(0.5, 0),半径R=0.5的圆内,极大地改善了系统的谱分布,从而使得Krylov子空间方法奏效.

对于预条件算子的近似逆,我们通过多重网格的方法获得.对于横纵不等间距采样的情况,多重网格中称为各向异性问题,从我们的数值试验可以看出,经典的多重网格方法对于此类问题已经失效.本文提出的方法对于中低频率的各向异性问题的求解是高效的;并且与直接法的内存消耗与耗时对比也证明了本文方法相比直接法的优势.

对于Helmholtz方程的高频迭代问题,仍然是一个难题,可能结合deflation技术(Nicolaides, 1987;Morgan, 1995; Frank and Vuik, 2001),消除小特征值对于迭代的影响,能加速高频Helmholtz问题的收敛.

致谢  感谢审稿专家和编辑提出的宝贵意见
References
Bayliss A, Goldstein C I, Turkel E. 1983. An iterative method for the Helmholtz equation. Journal of Computational Physics, 49(3): 443-457. DOI:10.1016/0021-9991(83)90139-0
Bérenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics, 114(2): 185-200. DOI:10.1006/jcph.1994.1159
Briggs W L, Van Emden H, McCormick S F. 2000. A Multigrid Tutorial. 2nd ed. Philadelphia, PA: Society for Industrial and Applied Mathematics.
Cao J, Chen J B, Cao S H. 2015. Studies on iterative algorithms for modeling of frequency-domain wave equation based on multi-grid precondition. Chinese Journal of Geophysics (in Chinese), 58(3): 1002-1012. DOI:10.6038/cjg20150325
Cao J, Chen J B, Dai M X. 2018. Modeling of frequency-domain scalar wave equation with the average-derivative optimal scheme based on a multigrid-preconditioned iterative solver. Journal of Applied Geophysics, 148: 70-82. DOI:10.1016/j.jappgeo.2017.10.006
Chen J B. 2012. An average-derivative optimal scheme for frequency-domain scalar wave equation. Geophysics, 77(6): T201-T210. DOI:10.1190/geo2011-0389.1
Chen Z Y, Cheng D S, Feng W, et al. 2011. A multigrid-based preconditioned Krylov subspace method for the Helmholtz equation with PML. Journal of Mathematical Analysis and Applications, 383(2): 522-540. DOI:10.1016/j.jmaa.2011.05.054
Chen Z Y, Cheng D S, Feng W, et al. 2013. An optimal 9-point finite difference scheme for the helmholtz equation with PML. International Journal of Numerical Analysis & Modeling, 10(2): 389-410.
De Zeeuw P M. 1990. Matrix-dependent prolongations and restrictions in a blackbox multigrid solver. Journal of Computational and Applied Mathematics, 33(1): 1-27.
Erlangga Y A, Vuik C, Oosterlee C W. 2004. On a class of preconditioners for solving the Helmholtz equation. Applied Numerical Mathematics, 50(3-4): 409-425. DOI:10.1016/j.apnum.2004.01.009
Erlangga Y A, Oosterlee C W, Vuik C. 2006a. A novel multigrid based preconditioner for heterogeneous Helmholtz problems. SIAM Journal on Scientific Computing, 27(4): 1471-1492. DOI:10.1137/040615195
Erlangga Y A, Vuik C, Oosterlee C W. 2006b. Comparison of multigrid and incomplete LU shifted-Laplace preconditioners for the inhomogeneous Helmholtz equation. Applied Numerical Mathematics, 56(5): 648-666. DOI:10.1016/j.apnum.2005.04.039
Erlangga Y A. 2008. Advances in iterative methods and preconditioners for the Helmholtz equation. Archives of Computational Methods in Engineering, 15(1): 37-66. DOI:10.1007/s11831-007-9013-7
Frank J, Vuik C. 2001. On the construction of deflation-based preconditioners. SIAM Journal on Scientific Computing, 23(2): 442-462.
Jo C H, Shin C, Suh J H. 1996. An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator. Geophysics, 61(2): 529-537. DOI:10.1190/1.1443979
Laird A L, Giles M. 2002. Preconditioned iterative solution of the 2D Helmholtz equation. Technical Report, NA 02-12, Computing Lab, Oxford University.
Morgan R B. 1995. A restarted GMRES method augmented with eigenvectors. SIAM Journal on Matrix Analysis and Applications, 16(4): 1154-1171. DOI:10.1137/S0895479893253975
Nicolaides R A. 1987. Deflation of conjugate gradients with applications to boundary value problems. SIAM Journal on Numerical Analysis, 24(2): 355-365.
Plessix R E. 2007. A Helmholtz iterative solver for 3D seismic-imaging problems. Geophysics, 72(5): SM185-SM194. DOI:10.1190/1.2738849
Riyanti C D, Kononov A, Erlangga Y A, et al. 2007. A parallel multigrid-based preconditioner for the 3D heterogeneous high-frequency Helmholtz equation. Journal of Computational Physics, 224(1): 431-448. DOI:10.1016/j.jcp.2007.03.033
Tsynkov S, Turkel E. 2001. A cartesian perfectly matched layer for the helmholtz equation.//Absorbing Boundaries and Layers, Domain Decomposition Methods: Applications to Large Scale Computers. Huntington, NY: Nova Science, 279-309.
Van Gijzen M B, Erlangga Y A, Vuik C. 2007. Spectral analysis of the discrete Helmholtz operator preconditioned with a shifted Laplacian. SIAM Journal on Scientific Computing, 29(5): 1942-1958. DOI:10.1137/060661491
曹健, 陈景波, 曹书红. 2015. 频率域波动方程正演基于多重网格预条件的迭代算法研究. 地球物理学报, 58(3): 1002-1012. DOI:10.6038/cjg20150325