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
区域 的整体 参数 矩阵
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 上, 如图
自 1976 年哥伦比亚大学研究人员首次将有限 元法引入到地震波研究领域以来, 许多学者就该方 法做了大量研究, 先后对比了模拟弹性波传播的有 限差分法和有限元法的精度[ 1] , 研究了地震波在不 同介质中的传播性质, 以及处理三维弹性波动方程 的思路; 研究了三角形有限单元、正方形有限单元、 不规则四边形有限单元以及矩形和三角形的混合有 限单元[ 2] 下地震波数值模拟, 以及地震波模拟的低 阶和高阶有限元法[ 3] 。有限元法有其独特的优势: 能够模拟任意不规则的地质体, 不受边界几何条件 的限制; 能够简单地处理自由边界条件; 可以方便地
大空间量和运算量的需求。此算法具有理论上的正确性和实践上的可行性。
关键词: 有限元; 并行算法; 弹性波动方程; 数值模拟; 块 三对角矩阵
中图分类号: P 315. 31
文献标识码: A
Finite element parallel algorithm for numerical solution of elastic wave equation
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)
模拟任意震源; 能适应变速不均匀介质情况等优点。 然而有限元法也有其巨大的局限性, 需要巨大的空 间量和计算量, 而计算机的空间和运算能力是有限 的, 外加地震波方程较其他领域的控制方程复杂, 所 以, 有限元不能满足在实际生产中的要求。随着计 算机硬件的发展, 高性价比微机群的出现和基于消 息传递的并行平台的不断完善, 为有限单元法在地 震波场正演模拟中的发展提供了新的空间。有限元 法的 化整为零、集零为整 的基本思想[ 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
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
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
,
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.
( 3)
式中, M 为系统质量矩阵; K 为系统刚度矩阵; V 为
节点位移矢量; V 为 V 对时间的二阶导数。
采用有限差分法来求解运动方程( 3) , 将最终化
为求解线性方程问题。
2 有限元并行求解
基于单 CPU 的串行算法的物理限制, 难以满足 有限元法高内存量和大运算量的要求, 制约了有限 元法的进一步发展。计算机硬软件技术的发展, 高性 价比微机群逐步代替了昂贵的工作站、巨型机进入
3 实现方案
采用多 CPU 的有限元并行算法求解弹性波动
方程, 有效地克服了传统的基于单一 CUP 的串行算 法的不足, 具体实现思路如下:
( 1) 把求解空 间离散为一系列有限单元, 并对
节点标号。
( 2) 根据 求解空 间确 定 CPU 的个 数 p ( p = 2q ) , 将求解空间等分到各个 CPU 上, 以确保生成的
王月英, 孙成禹
( 中国石油大学 地球资源与信息学院, 山东 东营 257061)
摘要: 在求解弹 性波动方程中, 有限元法的高 内存量和巨大运算 量的需求在基于单 CP U 串行算法中 一直难于满足,
制约其 优势的发挥。根据有限元法的 化整为零、集零 为整 的 基本思 想与并行 处理技 术的 分 而治之 的 原则基 本
收稿日期: 2006- 03- 09 作者简介: 王月英( 1977- ) , 女( 汉族) , 山东沾化人, 博士研究生, 主要从事地震波理论与方法研究。
28
中国石油大学学报( 自然科学版)
2006 年 10 月
的高内存量和巨大运算量的需求。
1 弹性波动方程数值解有限元算法
以二维问题为例, 根据变分原理和虚功原理, 在 无外力作用的情况下处于运动平衡状态的弹性体,
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
弹性波动方程数值解的有限元并行算法
一致, 采用基于多 CP U 的并行算法, 从有限 元参数矩阵计算和线 性方程组求解两个 方面入手, 把求解区域分到多 个
CPU 上并行计算参数矩阵, 对线性方程组采用循环块三 对角线方 程组进 行并行 求解。对比 了不同 大小空 间和不 同
CPU 个数下的加速比, 证实了多 CPU 的并行算法 能够克服基于单 CP U 串行算法的 物理限制, 满足了有 限元法的 巨
线性方程组并行求解方法分为直接法和迭代法
两类。线性方程组并行直接解法的优点是算法稳定