弹性波动方程数值解的有限元并行算法
2006年 第30卷 中国石油大学学报(自然科学版) V ol. 30 No. 5 第5期 Jour nal of China U niv ersity of Petroleum Oct. 2006
文章编号:1673-5005(2006) 05-0027-04
弹性波动方程数值解的有限元并行算法
王月英, 孙成禹
(中国石油大学地球资源与信息学院, 山东东营257061)
摘要:在求解弹性波动方程中, 有限元法的高内存量和巨大运算量的需求在基于单CP U 串行算法中一直难于满足, 制约其优势的发挥。根据有限元法的 化整为零、集零为整 的基本思想与并行处理技术的 分而治之 的原则基本一致, 采用基于多CP U 的并行算法, 从有限元参数矩阵计算和线性方程组求解两个方面入手, 把求解区域分到多个CPU 上并行计算参数矩阵, 对线性方程组采用循环块三对角线方程组进行并行求解。对比了不同大小空间和不同CPU 个数下的加速比, 证实了多CPU 的并行算法能够克服基于单CP U 串行算法的物理限制, 满足了有限元法的巨大空间量和运算量的需求。此算法具有理论上的正确性和实践上的可行性。关键词:有限元; 并行算法; 弹性波动方程; 数值模拟; 块三对角矩阵中图分类号:P 315. 31 文献标识码:A
Fin ite element parallel algorith m for nu merical solution of elastic w ave equ ation
WANG Yue -ying , SUN Cheng -yu
(Faculty of Geo -Resource and I nf or mation in China University of Petroleum ,
Dongying 257061, Shandong Pr ovince, China)
Abstract :T he numerical solution of the elast ic w av e equatio n r equires too much computer memory and calculatio n t ime to achiev e by finite element metho d based on single CPU (central processing unit). T he basic idea of finite element method is consistent w ith the basic principle of parallel algorithm. Based on mult-i CPU par allel algorithm and car ried out from matrix calculation of finite element and solution of linear equat ions, t he so lution r eg ion was located on many CPU s, and system lin -ear equations were parallel solved by using circle block -tr idiago nal equations. T hrough contr asting t he speedup -ratios of dif -ferent size of regions and different CPU numbers, it shows that t he parallel algorit hm on mult-i CPU can overcome the limits o f the single CPU o n memo ry and calculation rate, and be satisfied all requirements of the finite element method in the ca-l culation of the practical data processing. T hi s algorithm has theoretical cor rectness and feasibility in practice. Key words :finite element; par allel algorit hm; elastic w av e equation; numerical simulation; block -tridiagonal matr ix
自1976年哥伦比亚大学研究人员首次将有限元法引入到地震波研究领域以来, 许多学者就该方法做了大量研究, 先后对比了模拟弹性波传播的有限差分法和有限元法的精度[1], 研究了地震波在不同介质中的传播性质, 以及处理三维弹性波动方程的思路; 研究了三角形有限单元、正方形有限单元、不规则四边形有限单元以及矩形和三角形的混合有限单元
[2]
模拟任意震源; 能适应变速不均匀介质情况等优点。
然而有限元法也有其巨大的局限性, 需要巨大的空间量和计算量, 而计算机的空间和运算能力是有限的, 外加地震波方程较其他领域的控制方程复杂, 所以, 有限元不能满足在实际生产中的要求。随着计算机硬件的发展, 高性价比微机群的出现和基于消息传递的并行平台的不断完善, 为有限单元法在地震波场正演模拟中的发展提供了新的空间。有限元法的 化整为零、集零为整 的基本思想[4]与并行处理技术的基本原则 分而治之 基本一致[5]。笔者采用多CPU 的并行处理方法, 以有效满足有限元法
下地震波数值模拟, 以及地震波模拟的低
[3]
阶和高阶有限元法。有限元法有其独特的优势:
能够模拟任意不规则的地质体, 不受边界几何条件
的限制; 能够简单地处理自由边界条件; 可以方便地
收稿日期:2006-03-09
() , () , , ,
28
的高内存量和巨大运算量的需求。
中国石油大学学报(自然科学版) 2006年10月
各研究所和高等院校。在新的计算环境下继续研究有限元法, 增强了有限元法在实际资料处理中的可行性。采用并行求解, 把求解区域分到多个CPU 上, 既克服了单CPU 的物理限制, 又提高了运算速度。并行求解的思路可以从两个方面并行, 一方面是有限元参数矩阵的并行计算; 另一方面是线性方程组的并行求解。
2. 1 参数矩阵并行计算
把求解空间 划分为一系列正方形有限单元, 假设在x 和z 方向上有限单元个数分别为nx 和nz , 共有nx nz 个有限单元。各个有限单元的参数矩阵相互独立, 可以把空间 分成p 份, 在p 个CPU 上并行计算, 则每个CPU 上有限单元个数为nx nz /p , 先计算局部单元的参数矩阵, 最后集合成该
i 2i i
区域的整体参数矩阵M i MN , K 1 MN , L 1MN , K MN 和
2i L MN ,
p
1 弹性波动方程数值解有限元算法
以二维问题为例, 根据变分原理和虚功原理, 在无外力作用的情况下处于运动平衡状态的弹性体, 每一有限单元应满足文献[6]中式(18) , 化简式
(18) , 得
21e 22e e 21e
M e MN u e M +(v P K MN +v S K MN ) u M +(v P L MN -1e 22e e 2v 2S L MN +v S L MN ) w M =0;
22e 21e e 22e M e MN w e M +(v P K MN +v S K MN ) w M +(v P L MN -
2v S L MN +v S L MN ) u M =0.
(1)
其中M e MN =
e K 2MN e L 1MN
22e 21e e
= d x d z ,
= d x d z , L = d x d z .
z x x z
1e m n d x d z , K MN =m
n
A
A
A
m n
d x d z ,
R =
m n
2e MN
m n
i =1
R MN .
i
(4)
1
A A
式中, 数。
e u M
和
e w M
分别为x 和z 方向的位移函数; u M
e
式中, i 为CPU 序号(0 i
K 2, L 1和L 2。
求解空间的分配情况是将原串行算法中一个CPU 上的空间量和计算量分担到p 个CPU 上, 如图1所示。
为u e M 对时间的二阶导数; i 为局部有限元插值函
得到整体有限元运动方程为
M U i +K 1U i + L 1W i =0; 22
M W i +K W i +L U i =0. 其中
1e 22e
K 1=v 2P K MN +v S K MN ,
2e 21e K 2=v 2P K MN +v S K MN ,
(2)
L =v P L MN -2v S L MN +v S L MN ,
2e 22e 21e 2=v 2L P L MN -2v S L MN +v S L MN .
121e 21e 22e
U 和W 分别为x 和z 方向上的位移矢量。通常将式(2) 统一为一个运动方程式为
M V i +K V i =0.
节点位移矢量; V 为V 对时间的二阶导数。采用有限差分法来求解运动方程(3) , 将最终化为求解线性方程问题。
(3)
图1 整体有限元分到p 个C PU 上的分配情况
式中, M 为系统质量矩阵; K 为系统刚度矩阵; V 为
2. 2 线性方程组并行求解
线性方程求解是正演模拟中重要部分, 也是耗时最多的部分。在时间上需要递推求解运动方程, 采用传统的串行算法, 巨大的内存量使得巨型方程组的求解速度极慢。以159 159的剖分为例, 节点数为160 160个, 每一节点有两个位移分量u 和w , 每一个总体参数矩阵行数为51200, 带宽为323, 采用带状存储, 共需要内存约80MB 。此外, 若根据需要细分网格, 减小时间上的递推步长, 则需要更多的内存量和运算时间, 这是单CPU 串行算法的局限。
2 有限元并行求解
基于单CPU 的串行算法的物理限制, 难以满足有限元法高内存量和大运算量的要求, 制约了有限元法的进一步发展。计算机硬软件技术的发展, 高性
第30卷 第5期 王月英, 等:弹性波动方程数值解的有限元并行算法
29
为了节约空间, 在此保持式(1) 中的两式独立, 对其分别求解。从总体上看, 线性方程组系数矩阵为块状三对角矩阵, 如方程式A 1D 1C 2
A 2C
3
D 2A 3
D 3 C n-1
A n-1C n
D n-A n
T
1
部系数矩阵做一系列转化, 各CPU 之间的数据联合, 得到局部X i 值, 进而求得全部未知量数值, 详细求解思路参考文献[5]。对于不规则的地质构造, 多采用三角形有限单元来逼近, 此时生成的线性方程组系数矩阵不再是规则的块三对角矩阵, 则可以采用迭代法并行求解。
U 1U 2U 3 U n-1U n
(5) =
3 实现方案
采用多CPU 的有限元并行算法求解弹性波动方程, 有效地克服了传统的基于单一CUP 的串行算法的不足, 具体实现思路如下:
(1) 把求解空间离散为一系列有限单元, 并对节点标号。
(2) 根据求解空间确定CPU 的个数p (p =2q ) , 将求解空间等分到各个CPU 上, 以确保生成的局部整体参数矩阵大小相同。
(3) 各个CPU 先计算各个有限单元的参数矩阵, 再集合成该区域的局部整体参数矩阵M MN , K MN ,
2i 1i 2i K MN , L MN 和L MN , 最后生成局部线性方程组, 方程组
i
1i
(B 1 B 2 B 3 B n-1 B n ) .
式中, A i , C i 和D i 为160 160的三对角矩阵; U i 和
B i 为矢量, 0 i
线性方程组并行求解方法分为直接法和迭代法两类。线性方程组并行直接解法的优点是算法稳定性好、精度高。在此选用直接解法中的LU 分解法。采用循环块三对角线方程组的并行求解
[5]
的思路,
基于 分而治之 的思想, 各CPU 分别对所分配的局
中i 为CPU 序号(0 i
X im
C im+1A im +1C im+2
D im+1A im+2
D im+2 C (i+1) m-1
A (i +1) m-1C (i +
1) m
X i m+1X i m+2
D (i+1) m-1A (i+1) m
D (i +
1) m
B im+B im+
=
12
.
1
(6)
X (i+1) m-1X (i+1) m X (i+1) m+1
B (i+1) m-B (i+1) m
(4) 各个CPU 对局部线性方程组进行转化, 形成缩减方程组, 其明显的特征是从im +2到(i +1) m -1行对角线子矩阵为单位矩阵。
(5) 将前q 个CPU 中方程组中首、尾两行方程传给第1个CPU , 即方程式im +1和(i +1) m ; 将后q 个CPU 中方程组中首、尾两行方程传给第p 个CPU , 在首、尾两个CPU 中形成新的方程组。
(6) 首、尾CPU 对新方程组做转化, 得到系数矩阵的对角线子矩阵为单位矩阵的方程组, 将此方程组中首、尾两个方程相联合, 得到未知向量X 0和X (q+
1) m+1
(8) 重复步骤(4) ~(7) , 直至完成时间上的递推。
4 数值并行算例
以弹性波动方程的并行正演模拟为例, v P =6000m/s 和v S =3600m/s, 空间网格步长为 x = z =13m, 节点数为160 160, 采用等能量纵、横波源, 在有限元的节点上施加不同方向的力, 可以分
解为圆环切向和法向的分力, 其综合作用的效果既有法向上的爆炸源, 又有圆环切向上的纯剪切力源, 等效于激发能量相当的纵波和横波震源。子波选主频为20Hz 的雷克子波, 时间步长 t 为1ms 。图2(a) 和(b) 分别是在x 和z 方向上的位移分量u 和w 在t =200ms 时刻的波场快照。
的值, 代入此方程组, 求得全部解。
(7) 首、尾CPU 将X i m+1和X (i+1) m+1再传回对应的CPU, 各个CPU 根据X im+1和X (i+1) m+1, 求解局部线性方程组。
30
中国石油大学学报(自然科学版) 2006年10月
图2 位移分量u 和w 在200ms 时刻的波场快照
分别取2, 4, 8, 10个CPU 来运算, 在时间上采用10个递推步长, 得到的加速比如表1所示。表中m 为每个子矩阵的行列数, n 为总体矩阵所包含的矩阵行数, p 为CPU 个数。
表1 不同m , n 和p 下的加速比
m , n m =160, n =160m =160, n =640m =160, n =1280
p
21. 61. 71. 7
43. 23. 43. 5
85. 17. 36. 3
105. 38. 28. 4
参考文献:
[1] M AR FU RT Kurt J. Accuracy of finite -difference and f-i
nite -element modeling of the scalar and elastic w av e equa -tions[J]. Geophysics, 1984, 49(5) :533-549.
[2] K E Ben -x i, ZHAO Bo. 2-D finite element acoust ic wave
modeling including rugged topog raphy[R ].Presented at SEG Int l Exposition and A nnual M eeting, San A ntonio, T ex as, 2001.
[3] PADO VAN I E, PR IOL O E, SER IAN I G. L ow and
hig-t o rder finite element method, ex perience in sei smic modeling[J]. J Comp Acoust, 1994(2) :371-422. [4] 王瑁成. 有限单元法[M ].北京:清华大学出版社,
2003.
[5] 骆志刚. 典型结构大型线性方程组的分布式并行算法
研究[D]. 长沙:国防科学技术大学, 2000.
[6] 杜世通. 变速不均匀介质中波动方程的有限元法数值
解[J]. 华东石油学院学报, 1982, 6(2) :1-20.
DU Sh-i tong. Finite element numerical solution of wave propagation in non -homogenous medium with variable ve -locities[J]. Journal of East China Petroleum Institute, 1982, 6(2) :1-20.
[7] 杜世通. 固体连续介质中地震波微分方程式及其有限
单元法数值解[J].华东石油学院学报, 1985, 9(2) :1-9. DU Sh -i tong. Di fferential equations seismic w aves in soli d medium with continuously varying velocities [J].Journal of East China Petroleum Institute, 1985, 9(2) :1-9.
[8] DI M IT RI Komatitsch, CHRIST OPHE Barnes, JEROEN
T romp. Simulatio n o f anistropic w ave pr opagation based upon a spectral element method[J].Geophysics, 2000, 65(4) :1251-1260.
从表1可以看出, 采用多CPU 的并行算法, 提高了运算速度, 缩减了运算时间; 随着CPU 个数的增加, 加速比逐渐增大。CPU 之间的通信需要一定时间, 当求解空间一定时, 即m 和n 一定, 存在最佳
并行效率的CPU 个数; m 和n 不同, 相同个数CPU 下其加速比也不同, 因此合理的编号有助于算法的经济化。通常选有限单元个数少的一方为编号起始方, 即此时m 值小。
5 结束语
随着计算机技术的发展, 高性价比的PC Clus -ter 开始走进各个研究所和高等院校, 给有限元法提供了新的空间。采用基于多CPU 的并行算法做地震波有限元法正演模拟, 克服了基于单CPU 的串行算法物理限制, 扩大求解区域, 提高了运算速度, 增强了有限元法实际资料处理的可行性, 使得有限元法的优势得到了很好的发挥。
致谢 在本研究过程中, 得到了杜世通教授的悉心指导和帮助, 在此表示诚挚的感谢。
(