弹性波动方程数值解的有限元并行算法

K
1e MN
)
w
e M
+
(
v
2 P
L
2e MN
-
2v
2 S
L
2e MN
+
v
2 S
L
1e MN
)
u
e M
=
0.
( 1) 其中
MeMN =
m
ndx dz ,
K
1e MN
=
m
x
n
x
d
x
d
z
,
A
A
K
2e MN
=
A
m
z
zndxdz ,
L
1e MN
=
A
m
z
n
x
d
x
dz
,
L
2e MN
=
A
m
x
n
z
d
x
d
z
.
式中,
u
大空间量和运算量的需求。此算法具有理论上的正确性和实践上的可行性。
关键词: 有限元; 并行算法; 弹性波动方程; 数值模拟; 块 三对角矩阵
中图分类号: P 315. 31
文献标识码: A
Finite element parallel algorithm for numerical solution of elastic wave equation
e M
和
w
e M
分别为 x
和z
方向的位移函数;
u
e M
为 ueM 对时间的二阶导数; i 为局部有限元插值函
数。
得到整体有限元运动方程为
M Ui + K 1 U i + L1 Wi = 0;
MW i + K 2 Wi + L 2 Ui = 0.
( 2)
其中
K1 =
v
2 P
K
1e MN
+
v
2 S
KM2eN
3 实现方案
采用多 CPU 的有限元并行算法求解弹性波动
方程, 有效地克服了传统的基于单一 CUP 的串行算 法的不足, 具体实现思路如下:
( 1) 把求解空 间离散为一系列有限单元, 并对
节点标号。
( 2) 根据 求解空 间确 定 CPU 的个 数 p ( p = 2q ) , 将求解空间等分到各个 CPU 上, 以确保生成的
收稿日期: 2006- 03- 09 作者简介: 王月英( 1977- ) , 女( 汉族) , 山东沾化人, 博士研究生, 主要从事地震波理论与方法研究。
28
中国石油大学学报( 自然科学版)
2006 年 10 月
的高内存量和巨大运算量的需求。
1 弹性波动方程数值解有限元算法
以二维问题为例, 根据变分原理和虚功原理, 在 无外力作用的情况下处于运动平衡状态的弹性体,
模拟任意震源; 能适应变速不均匀介质情况等优点。 然而有限元法也有其巨大的局限性, 需要巨大的空 间量和计算量, 而计算机的空间和运算能力是有限 的, 外加地震波方程较其他领域的控制方程复杂, 所 以, 有限元不能满足在实际生产中的要求。随着计 算机硬件的发展, 高性价比微机群的出现和基于消 息传递的并行平台的不断完善, 为有限单元法在地 震波场正演模拟中的发展提供了新的空间。有限元 法的 化整为零、集零为整 的基本思想[ 4] 与并行处 理技术的基本原则 分而治之 基本 一致[ 5] 。笔者 采用多 CPU 的并行处理方法, 以有效满足有限元法
各研究所和高等院校。在新的计算环境下继续研究
有限元法, 增强了有限元法在实际资料处理中的可
行性。采用并行求解, 把求解区域分到多个 CPU 上, 既克服了单 CPU 的物理限制, 又提高了运算速度。
并行求解的思路可以从两个方面并行, 一方面是有
限元参数矩阵的并行计算; 另一方面是线性方程组
的并行求解。
每一 有 限单 元应 满足 文献 [ 6] 中式 ( 18) , 化 简 式
( 18) , 得
MeMN
u
e M
+
(
v
2 P
K
1e MN
+
v
2 S
K
2e MN
)
u
Baidu Nhomakorabea
e M
+
(
v
2 P
L
1e MN
-
2v
2 S
L
1e MN
+
v
2 S
L
2e MN
)
w
e M
=
0;
MeMN
w
e M
+
(
v
2 P
KM2eN
+
v
2 S
一致, 采用基于多 CP U 的并行算法, 从有限 元参数矩阵计算和线 性方程组求解两个 方面入手, 把求解区域分到多 个
CPU 上并行计算参数矩阵, 对线性方程组采用循环块三 对角线方 程组进 行并行 求解。对比 了不同 大小空 间和不 同
CPU 个数下的加速比, 证实了多 CPU 的并行算法 能够克服基于单 CP U 串行算法的 物理限制, 满足了有 限元法的 巨
( 3)
式中, M 为系统质量矩阵; K 为系统刚度矩阵; V 为
节点位移矢量; V 为 V 对时间的二阶导数。
采用有限差分法来求解运动方程( 3) , 将最终化
为求解线性方程问题。
2 有限元并行求解
基于单 CPU 的串行算法的物理限制, 难以满足 有限元法高内存量和大运算量的要求, 制约了有限 元法的进一步发展。计算机硬软件技术的发展, 高性 价比微机群逐步代替了昂贵的工作站、巨型机进入
,
K2 =
v
2 P
K
2e MN
+
v
2 S
KM1eN
,
L1 =
v
2 P
LM1eN
-
2v
2 S
L M1 eN
+
v
2 S
LM2eN
,
L2 =
v
2 P
LM2eN
-
2v
2 S
L M2 eN
+
v
2 S
LM1eN
.
U 和 W 分别为 x 和 z 方向上的位移矢量。通常将式
( 2) 统一为一个运动方程式为
MVi + K Vi = 0.
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 linear equations were parallel solved by using circle block-tr idiago nal equations. T hrough contr asting t he speedup- ratios of different 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 his 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
2006 年 第 30 卷 第 5期
中国石油大学学报( 自然科学版) Jour nal of China U niv ersity of Petroleum
文章编号: 1673- 5005( 2006) 05- 0027- 04
V ol. 30 No. 5 Oct. 2006
弹性波动方程数值解的有限元并行算法
线性方程组并行求解方法分为直接法和迭代法
两类。线性方程组并行直接解法的优点是算法稳定
性好、精度高。在此选用直接解法中的 LU 分解法。 采用循环块三对角线方程组的并行求解[ 5] 的思路,
基于 分而治之 的思想, 各 CPU 分别对所分配的局
部系数矩阵做一系列 转化, 各 CPU 之 间的数据联 合, 得到局部 X i 值, 进而求得全部未知量数值, 详细 求解思路参考文献[ 5] 。对于不规则的地质构造, 多 采用三角形有限单元来逼近, 此时生成的线性方程 组系数矩阵不再是规则的块三对角矩阵, 则可以采 用迭代法并行求解。
( B1 B 2 B3
Cn A n
Un
Bn- 1 Bn ) T .
( 5)
式中, Ai , Ci 和Di 为 160 160 的三对角矩阵; Ui 和 Bi 为矢量, 0 i < 160。总内存空间约为 60 MB, 分
担到 p 个 CPU 上, 每一台并行机分担约 60/ p M B,
减小了串行算法中单个 CPU 上的内存压力。
自 1976 年哥伦比亚大学研究人员首次将有限 元法引入到地震波研究领域以来, 许多学者就该方 法做了大量研究, 先后对比了模拟弹性波传播的有 限差分法和有限元法的精度[ 1] , 研究了地震波在不 同介质中的传播性质, 以及处理三维弹性波动方程 的思路; 研究了三角形有限单元、正方形有限单元、 不规则四边形有限单元以及矩形和三角形的混合有 限单元[ 2] 下地震波数值模拟, 以及地震波模拟的低 阶和高阶有限元法[ 3] 。有限元法有其独特的优势: 能够模拟任意不规则的地质体, 不受边界几何条件 的限制; 能够简单地处理自由边界条件; 可以方便地
区域 的整体 参数 矩阵
MiM N
,
K
1i MN
,
K
2i MN
,
L
1i MN
和
L
2i MN
,
p
R=
R
i MN
.
i= 1
式中, i 为 CPU 序号( 0
K 2, L 1 和 L 2。
( 4) i < p ); R 代表 M , K1,
求解空间的分配 情况是将原串行算法中一个
CPU 上的空间量和计算量分担到 p 个 CPU 上, 如图
第 30 卷 第 5 期
王月英, 等: 弹性波动方程数值解的有限元并行算法
29
为了节约空间, 在此保持式( 1) 中的两式独立, 对其分别求解。从总体上看, 线性方程组系数矩阵为 块状三对角矩阵, 如方程式
A1 D1
U1
C2 A2 D2
U2
C3 A 3 D3
U3 =
Cn- 1 A n- 1 Dn- 1 Un- 1
王月英, 孙成禹
( 中国石油大学 地球资源与信息学院, 山东 东营 257061)
摘要: 在求解弹 性波动方程中, 有限元法的高 内存量和巨大运算 量的需求在基于单 CP U 串行算法中 一直难于满足,
制约其 优势的发挥。根据有限元法的 化整为零、集零 为整 的 基本思 想与并行 处理技 术的 分 而治之 的 原则基 本
2. 1 参数矩阵并行计算
把求解空间 划分为一系列正方形有限单元,
假设在 x 和z 方向上有限单元个数分别为 nx 和nz , 共有 nx nz 个有限单元。各个有限单元的参数矩
阵相互独立, 可以把空间 分成 p 份, 在 p 个 CPU
上并行计算, 则每个 CPU 上有限单元个数为 nx
nz / p , 先计算局部单元的参数矩阵, 最后集合成该
1 所示。
图 1 整体有限元分到 p 个 CPU 上的分配情况
2. 2 线性方程组并行求解 线性方程求解是正演模拟中重要部分, 也是耗
时最多的部分。在时间上需要递推求解运动方程, 采 用传统的串行算法, 巨大的内存量使得巨型方程组 的求解速度极慢。以 159 159 的剖分为例, 节点数 为 160 160 个, 每一节点有两个位移分量 u 和 w , 每一个总体参数矩阵行数为 51 200, 带宽为 323, 采 用带状存储, 共需要内存约 80 MB。此外, 若根据需 要细分网格, 减小时间上的递推步长, 则需要更多的 内存量和运算时间, 这是单 CPU 串行算法的局限。
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)