36
内蒙古石油化工 2011年第24期
多重网格法在求解泊松方程中的应用进展
杨金凤, 邓居智, 陈 辉
(东华理工大学放射性地质与勘探技术国防重点学科实验室, 江西抚州 344000)
摘 要:多重网格(Multigrid, 简称M G) 法是求解复杂科学运算中由偏微分方程离散以后所得的大型代数线性方程组的最高效、最快速的迭代求解方法。本文详述了多重网格法在求解二维、三维泊松方程求解问题研究中的应用, 与传统迭代方法相比在收敛速度、精度和稳定性上多重网格方法都有其固有的高效性。并初步的讨论了多重网格并行算法求解泊松方程正演问题, 针对多重网格法应用于地球物理的各种数值模拟研究, 提出了四点具体的研究发展方向, 希望对今后研究有所帮助。
关键词:多重网格法; 泊松方程; 正演; 多重网格并行算法 中图分类号:O343 文献标识码:A 文章编号:1006—7981(2011) 24—0036—03 多重网格法最初是由苏联计算数学家Fedo renk 在20世纪60年代提出的, 但当时并未引起人们足够的注意。其真正兴起是在近30年, Hackbusch 等在数学上证明了对于线性椭圆型偏微分方程, 多重网格算法是最优的, 其计算量仅与网格节点数的一次方成正比, 把解方程所需计算量降到最低量级, 在计算效率上远远超过了其它方法。国内学者许进超在用多重网格算法求解线性偏微分方程时, 引进泛函分析工具研究这一问题, 取得了较为满意的效果, 自此这方面的研究得到了快速的发展。该方法是现今数值计算领域最活跃的分支之一, 在计算结构力学、流体力学、波动方程等诸多领域也得到广泛的应用,
[1][2]
许多成果总结在Hackbusch 、刘超群等著作中。在计算地球物理领域, 该方法的应用也逐渐得到关注, 不过还多限于二维问题。1 多重网格法的基本思想
多重网格算法的最初动力来源于对网格方程迭代求解[3], 误差的各个Fourier 分量的不同衰减程度, 在分析线性方程组的收敛性时, 不同的频率分量具有不同的误差和不同的收敛速度。频率高、变化快的分量称为高频分量, 其振荡误差是局部行为, 来源于附近几个网格点之间的相互耦合, 与边界或距离较远的网格点信息无关, 因而在传统的迭代过程中衰减很快; 而把频率低、变化慢、较光滑的分量称光滑分量, 其光滑误差是全局行为, 主要来源于边界信息, 因而在传统的迭代过程中衰减很慢。传统的点或块松弛都是局部性较强的方法, 因此它们能迅速抹平局部性的高频振荡误差, 但对全局性的低频光滑误差却衰减缓慢。实际上, 经过初始几次迭代后, 误差将呈现光滑性。这就是迭代法在求解线性方程组开始时迭代速度相当快, 而在随后的迭代过程中收
敛速度逐渐变慢的原因。多重网格法整个迭代过程可以概括如下:①在细网格层上用传统迭代法迅速衰减其中的那些高频误差分量; ②接着通过插值算子转到粗网格层上, 这时的低频误差分量变为高频误差分量, 同样利用迭代法进行衰减; ③再利用延拓算子返回到细网格层对初始解进行修正得到近似解, 再进行下一个多重网格循环以达到一定的收敛精度为止。2 泊松方程的正演问题
众所周知, 与地球物理有关的数学问题实际上就是偏微分方程的正问题与反问题。正问题就是求解这个微分方程的定解问题, 从而得到过程与状态的数学描述。泊松方程是一类应用范围极广的偏微分方程, 它是流体力学、传热传质学、弹性力学、水动力学及电磁学等基础理论研究的基本控制方程之一。在直流电法勘探和直流激化法中所建立的地下电场均可作为稳定电流场处理, 若地下充满各向同性均匀的介质, 则从稳定电场强度与电势梯度的关
系E =-gradU 和高斯定理的微分形 ·E =, 即
t
可得出稳定电场的泊松方程(1):
2222
++= U = x y z t
(1)
·(σ U)=-I δ(x-x 0) (y-y 0) (z-z 0) (2)
式中U 表示电位, q 表示自由电荷的密度, t 表示时间。(2) 式为非均匀介质的泊松方程, 也是直流电法的基本微分方程, 二维的泊松方程类似。
正演计算主要就是利用有限差分法、有限元法、有限体积法等将泊松方程离散为一个线性代数方程
收稿日期:2011-09-12
基金项目:东华理工大学研究生创新基金(DY CA10015) 资助。:() 女, 。
2011年第24期 杨金凤等 多重网格法在求解泊松方程中的应用进展组, 再采用经典迭代算法求出泊松方程的数值解。由于只有少量的的泊松方程有解析解, 因此研究该方程的数值求解方法具有非常重要的意义。3 泊松方程正演中的多重网格法
多重网格方法在数学界的理论研究, 通过松弛过程、限制算子、插值算子等相应的改进, 使得在各种复杂情况下, 其有效性并没有丧失, 为其应用于实
[4]
际复杂问题奠定了坚实的理论基础。葛永斌等利用六阶紧致差分格式结合多重网格法求解了二维泊松方程的Dirichlet 边值问题, 并采用不同的松弛算子将其与四阶精度格式的多重网格方法进行了比较。得出了在不增加计算量的前提下四阶精度格式的多重网格方法具有更高的精确度和收敛阶, 而且ZLGS 迭代对四阶精度还和六阶精度格式的多重网格算法, 都是一种较其他松弛算子更加有效的一种算子选择, 对误差残量的“光滑”效果都较好。
[5]
Robert 研究了多重网格法应用于直流电阻率二维正演问题, 利用多重网格方法求解稳定电流场的泊松方程, 并计算了几个实际的电磁问题, 指出此方法虽能以较高的速度和精度进行数值模拟, 但是其最终的收敛性受均匀导体与围岩的电性差异的影响; 此外, 网格剖分(多重网格所需的) 时应该考虑导体的各向异性的问题, 在典型的多重网格法中并未得以解决, 而这一步对其随后的反演问题非常重要。针对这些问题, Robert 用泊松方程九点有限差分算子代替传统的五点有限差分算子来提高多重网格算法的稳定性和精确性, 并且研究了新的模型粗化方法, 以保持导体的各向异性, 同时也能正确反映粗网格
[6]
层上的电性差异。为电阻率二维正演研究提供了新的方向, 具有实际的意义, 使得其快速反演成为可能。多重网格一般是建立在有限差分的基础上, 王青平等人[7]将多重网格引入到有限元方法中, 并以二维泊松方程为例, 验证了算法的精确性与高效性, 得出了多重网格方法与有限元可以结合的结论, 并且与单重网格中高斯-赛德尔迭代(GS ) 、共轭梯度法(CG ) 在迭代次数、程序运行时间以及误差等方面进行比较, 随着网格加密多重网格法有较明显的优势。但是也伴随着区域不宜太重杂, 无法进行多网格剖分等不足, 以及从数学角度上讲对其收敛性的分析比较困难。
对于二维泊松方程, 已存在不少高精度紧致格式多重网格算法。但是对于三维泊松方程的高精度紧致差分法, 由于格式的构造比二维具有更高的代
[8]
数复杂性要求, 所以对此研究不是很多。葛永斌等利用对称网格点泰勒展开式中各阶导数项明显对称的性质, 得到了数值求解三维泊松方程的四阶和六阶精度的紧致差分格式, 同时采用多重网格加速迭代技术, 设计了相应的多重网格算法, 求解三维泊松, 37
解高维问题时计算量大、收敛速度慢的缺陷。相应的数值试验表明多重网格的加速效果非常明显, 高精度紧致格式也达到了所期望的精度, 为精确求解复杂问题数值模拟中的泊松方程提供了理论依据。
随着三维勘探的展开, 地球物理三维正反演越来越重要, 而首先要解决的就是三维正演问题。因此, 快速、精确地求解三维泊松方程, 在计算地球物理中具有重要意义。乔中林将多重网格方法应用于地球物理模型不等距剖分中, 推导出三维模型不等距剖分情况下的限制算子和插值算子的公式形式, 同时编写出了不等距剖分多重网格方法求泊松方程类偏微分方程的C 程序代码。多重网格方法在电法模型正演上的应用迈出了一步, 使得正演迭代速度得到了提高[9]。对于现阶段地球物理三维勘探面临越来越复杂的问题, 其要求更精细的三维网格剖分, 为克服这一问题, 2009年鲁晶津, 吴小平[10]等人将多重网格方法应用于三维泊松方程的数值模拟中, 网格尺度逐渐增加至257×257×257, 并与不完全Chelesky 共轭梯度法(ICCG) 和Gauss 直接解法进行比较, 结果表明无论是在所耗内存方面还是在收敛速度与网格关系上, 多重网格方法的效果都相当明显, 是精确高效的三维数值模拟方法。随后鲁晶津[11]等人又用代数多重网格(AM G ) 方法求解三维直流电阻率法正演模拟形成的有限差分线性方程组, 研究了求解二次场的方法, 通过求解二次场的方法解决了一般的直接求解实际场所带来的总场中由点电源导致的奇异性的问题, 从而获得快速、精确的三维电阻率数值模拟, 也克服了常用的几何多重网格法在处理带跃变系数的偏微分方程时存的缺陷, 并对两个存在很大电性差异的模型进行模拟计算, 与不完全Chelesky 共轭梯度法相比具有更高的计算效率和稳定性[12]。汤井田等[13]将多重网格方法应用到直流电阻率三维自适应有限元法(AFEM ) 中弥补了人为剖分大型三维地球物理模型不能满足的高精度和高效率需求的不足, 模型验证了算法的有效性。
近年来, 随着高性能计算机系统的日益普及和性能的大幅提升, 受实际应用领域中大规模科学计算问题的驱动, 多重网格算法并行化引起了人们的
[14]
普遍重视, 也将成为求解大规模地球物理正反演问题的一种强大趋势。
为保持原有算法的数值高效性, 应在尽可能不影响算法收敛性前提下考虑并行, 多重网格法并行化工作以串行程序的并行化为主。并行化方法主要包括:①对规则区域采用网格划分, 将原始网格划分成几块, 分配给各台处理机, 并且相互包含相邻网格块的拟边界信息; ②对复杂区域采用块结构方法, 根据几何外形将区域分成几块到几十块, 每块适当变,
38
算。
[15]
内蒙古石油化工 2011年第24期
[2] 刘超群. 多重网格法及其在计算流体力学中的
应用[M ].北京:清华大学出版社, 1995.
[3] 李晓梅, 莫则尧. 多重网格算法综述[J ].中国
科学基金, 1996(1):4～11.
[4] 葛永斌, 吴文权, 卢曦. 基于二维泊松方程六阶
紧致格式的多重网格方法[J ].上海理工大学学报, 2002, 24(4):337～344.
[5] M oucha R . Multig rid methods fo r forw ard
and inv erse resistivity problem s in geoph ysics [D ].Ca nada :Univ ersity of To ronto , 2003.
[6] M oucha R B R C. An Accurate And Robust
M ultig rid Algo rithm Fo r 2D Forw ard Resistivity M odelling [J ].Geophy sical Prospecting , 2004, 52(3):197～212.
[7] 王青平, 白武明, 王洪亮. 多重网格在二维泊松
方程有限元分析中的应[J].地球物理学进展, 2010, 25(4):1467～1471.
[8] 葛永斌, 田振夫, 马红磊. 三维泊松方程的高精
度多重网格解法[J].应用数学, 2006, 19(2):313～318.
[9] 乔中林. 直流电法三维正演多重网格算法研究
[D].北京:中国地质大学, 2007.
[10] 鲁晶津, 吴小平, spitzerk. 三维泊松方程数值
模拟的多重网格方法[J].地球物理学进展, 2009, 24(1):154～158.
[11] 鲁晶津, 吴小平, spitzerk. 直流电阻率三维正
演的代数多重网格方法[J].地球物理学报, 2010, 53(3):700～707.
[12] 吴小平, 徐果明, 李时灿. 利用不完全
Cholesky 共轭梯度法求解点源三维地电场[J].地球物理学报, 1998, 41(6):848～855.
[13] TAN G Jing -tian , W AN G Fei -yan , etc . 3
-D direct current resistivity forw ard m odeling by adaptiv e multig rid finite element method [J ].Journa l o f Central So uth University of Technolog y , 2010, 17(3):587～592.
[14] Xu Z Q . Distributing multig rid computation
:on m essage passing system [D ].Wuhan
Hua zho ng Univ ersity of Science a nd Technolog y , 1993.
[15] 罗铁祥. 多重网格区域分裂分布式计算[J ].
地球科学, 2001, 26(3):323～327.
[16] 郭庆平, 王高, 卫高. 多重网格方法的一种并
行模式-交错多重网格方法剖析[J ].武汉理工大学学报(交通科学与工程版) , 2002, 26(4):429～431.
2001年罗铁祥对多重网格分布式计算进行
了研究, 并以二维泊松方程为例, 基于区域分列式实现了算法的并行化, 其顺序算法描述的是非递归形式。网状拓扑结构组织在多处理机上, 并行算法映射到多进程上, 在一定程度上显著提高并行化速度和并行化效率。算法思想简单, 具有广泛的实用性。郭庆平等[16]剖析了多重网格方法的一种并行模式--交错多重网格方法, 该方法以二维泊松方程的边值问题为例, 将区域划分成网格后, 在网格中交错取点分布到不同处理器上去形成较粗网格, 而在各个粗网格上用多重网格方法计算。通过交错取点的划分方式来加快边界信息的传递。对这种并行处理模式的性状作了相应的分析和比较, 对于交错多重网格方法中的计算复杂度和通信复杂度, 以及相应的处理机效率, 也进行了定量的分析。4 结论与展望
综上所述, 利用多重网格方法求解泊松方程可以得到较高精度的数值解, 而将其应用于电阻率正演数值模拟中, 算例结果表明与传统的几种迭代法相比, 多重网格法可获得较精确的三维数值模拟, 在收敛速度、稳定性、所耗内存与网格关系上, 也具有其固有的高效性。现如今, 随着找矿的程度的不断深入, 电磁场勘探中的各种方法(M T 、CSAM T) 的正演问题逐渐成为热点, 也是难点问题, 而将多重网格方法应用到其中的却很少。应进一步从以下几个方面展开多重网格法应用于地球物理(如地电) 的各种数值模拟研究中:①要不断的关注数学界关于多重网格的最新研究成果, 如对限制、插值算子、松弛方法等的改进, 以便能够使多重网格在泊松方程求解中得到更广泛的应用, 得到更高精度的数值解。如:多重网格梯度法(M Gg M ) 、多重网格共轭梯度法(M Gcg M ) 。②多重网格的研究大多都是针对有限差分法进行的, 还应考虑与有限元、有限体积、积分方程等正演方法结合。还可以将非结构网络、残差校正与双离散技术、自适应技术与多重网格法结合, 寻求其中最适合高效求解复杂介质中的泊松方程的方法。③随着计算机多核CPU 的普及和盛行, 迫切需要研究多重网格并行算法, 将多重网格并行算法结合并行计算机的特点, 进行大规模的并行计算, 实现对复杂介质中的泊松方程离散后所得的大型稀疏线性方程组更快速、更高效的求解。④多重网格算法在泊松方程中求解有着快速, 高效的特点, 因此可以将多重网格算法引入到求解复杂的M AXW ELL 方程组中, 以便实现快速、高精度、适应性强电磁正演模拟。
[参考文献]
[1] Hackbusch W . M ultig rid methods and
application [M].Berlin :Spring er, 1985.