结构力学有限元分层并行计算方法

规模和条件数也随之急剧增加，从而给其求解带来 困难[9]。另外在迭代求解界面方程的过程中，由于 所有进程都要参与全局通信，进程间通信和同步开 销的增加也会极大地降低并行效率[10]。 在硬件方面，具有分布式存储体系结构的大规 模并行计算机（Massively Parallel Processors, MPP） 是目前进行科学和工程计算最广泛使用的并行计算 机之一。 对于 MPP 来说， 最为重要的部分就是数据 的存储以及处理器之间的相互通信和协作，这也是 影响并行效率的很重要因素。 因此， 利用 MPP 提高 并行效率的关键在于处理好大规模数据的存储以及 处理器间的相互通信和协作问题。 为解决采用区域分解法利用大量 CPU 求解大 规模问题并行效率低的问题， 本文针对 MPP 分布式 存储体系结构提出了一种结构力学有限元分层并行 计算方法。该方法在两级分区两次缩聚的基础上不 仅实现了大量数据的分布式存储，提高了数据的内 存访问速率；而且实现了计算过程的三层并行，将 节点内通信与节点间通信分离从而有效提高了通信 效率；此外，它还进一步降低了界面方程的规模， 加快了其迭代收敛速度。鉴于有限元分析的类型十 分广泛，为有针对性地分析和解决问题本文将研究 重点集中在结构线性静力问题的分析上。
3
数值算例
图 3 三层并行机制 Fig.3 Scheme of three-layer paraBiblioteka Baidulelization
第二层并行，各节点机内分别同时进行相应 1 级子区域的组集、并行缩聚和回代求解，通信只存 在于同一节点机内的不同核之间。各节点机内可分 别指定一个局部主进程，统一负责相应 1 级子区域 的组集，数据分发，并行缩聚结果的汇总以及内部 自由度的回代求解。各节点机内相应 1 级子区域的 并行缩聚过程则由同一节点内所有进程同时参与。
则有限元网格首先被剖分为 2 个 1 级子区域，然后 每个 1 级子区域又进一步被剖分为 4 个 2 级子区域。
（1）
式中， UI, UB, PI, PB 分别为内部节点和边界节点 对应的位移和外部载荷，K 为刚度矩阵。 消去内部节点自由度后，得到只含边界节点自 由度未知量的界面方程
{U I } [ K II ]1 ({PI } [ K IB ]{U B })
（5 ）
为节省内存和减少计算量，实际计算时一般不 直接采用矩阵求逆缩聚，而是通过对稀疏存储的刚 度矩阵进行分解实现缩聚。本文各子区域局部总体 刚度矩阵全部采用 skyline 格式存储，缩聚根据 Han[16]等提出的改进乔列斯基分解法实现。
1
传统区域分解法
—
—
收到第 1 稿，
—
—
收到修改稿.
1）国家 863 高技术研究发展计划（2012AA01A307） ，国家自然科学基金（11072150,61073088）. 2） 金先龙， 教授， 博士生导师， 主要研究方向： 高性能计算、 数值仿真； Professor Xianlong Jin, Ph.D. His research interest is high performance computing and numerical simulation. E-mail ： jxlong@sjtu.edu.cn
上海超级计算中心的 “魔方”超级计算机是基于 国际主流的 MPP 概念设计的大型计算机系统。 本文 所有算例都是在“魔方”C 区测试运行通过。C 区每 个节点共计 16 核，64GB 内存，平均每个处理器核 心 4GB 内存。节点机间采用 Infiniband 光纤网络互 联，理论带宽 20Gbits 每秒。并行计算程序采用 FORTRAN 编写，进程间通信根据 MPI 标准实现。 3.1 单工况问题分析 并行算法评估的有限元模型分别如图 4 和图 5 所示。图 4 为平面应力问题的有限元模型，平板左 端固定，右端受 q=10MPa 均布载荷作用。平板长 30.72m，宽 5.12m，材料的弹性模量为 2.1e5MPa， 泊松比为 0.3。采用四边形等参单元进行网格剖分 后，该模型具有 11,796,480 单元，11,827,585 节点， 23,655,170 自由度。图 5 为某高层建筑框架结构有 限元模型，需分析该模型在楼板压力下的变形。模 型总高度为 48m，16 层，层高为 3m。采用六面体 单元进行网格剖分后，该模型具有 3,009,600 单元，
**
（上海交通大学 机械与动力工程学院，上海 200240)
（上海超级计算中心，上海 201203)
摘要：根据 MPP 分布式存储的体系结构特点，提出了一种结构力学有限元分层并行计算方法。该方法在两级分 区两次缩聚策略的基础上不仅实现了大量数据的分布式存储，提高了数据的内存访问速率；而且实现了计算过程 的三层并行，有效提高了通信效率；此外，它还进一步降低了界面方程的规模，大幅度减少了界面方程的求解时 间。因此,它能够充分利用 MPP 分布式存储的体系结构特点提升大规模并行计算效率。最后通过典型数值算例验 证了该方法的正确性和有效性。 关键词：大规模并行计算机，分层并行计算，区域分解，缩聚 中图分类号：TP311，O246 文献标识码：A DOI：10.6052/0459-1879-13-335
每个 1 级子区域局部总体刚度矩阵也采用 skyline 格式存储，但缩聚时采用并行改进乔列斯基 分解实现。具体过程为：首先由局部主进程通过组 集同一节点内所有 2 级子区域的界面方程形成相应 1 级子区域的系统方程，并以循环列的方式将刚度 矩阵轮廓数据平均分配给同一节点内的所有进程； 同一节点内的各进程根据分配到的分布式轮廓数据 采用活动列求解技术[17-18]进行并行改进乔列斯基分 解，然后将分解结果发送给相应的局部主进程；局 部主进程将分解结果收集汇总，完成相应 1 级子区 域有效刚度矩阵和有效载荷向量的计算。 第三层并行， 利用并行 PCG 算法求解 1 级子区 域界面方程，每个节点机只有一个进程——局部主 进程参与求解和通讯。该方法不需要组装形成总体 界面方程。求解时各 1 级子区域的有效刚度矩阵和 有效载荷向量仍分布式存储在相应的节点机上，中 间计算结果也以矩阵向量积的形式分布存储。局部 通信仅存在于相邻 1 级子区域间，只有少量的点积 操作和整体迭代误差的计算需要全局通信。因此， 它能够有效减少通信量，提高并行计算效率。 综上所述，分层并行计算方法能够将大量局部 通信限制在各节点机内部，并确保各节点机只有一 个进程参与全局通信。这就实现了节点内通信与节 点间通信的分离， 减少了进程间通信和同步的开销， 有效提高了通信效率。
引
言
随着科学研究和工程技术的不断发展，出现了 诸如跨江隧道，高速列车，摩天大厦等大型和超大 型的复杂结构。这些结构不但自由度高，而且还含 有非线性和复杂的边界条件等多种因素，因此对它 们的分析需要借助大规模有限元模型进行大量计 算。由于计算规模巨大，求解复杂，这些问题在传 统串行机上无法得到满意的解答。利用并行计算机 研究和开发相应的并行算法则为这类问题的解决提 供了切实可行的方法。 直接法[1-2]和迭代法[3-4]是有限元并行求解的两 种基本算法。直接法通过对系数矩阵进行三角分解 能够在预期的操作内得到方程的解。然而随着有限 元规模的增加，直接法所需的内存空间和计算量也 会迅速增加。 迭代法所需内存小， 容易实现并行化。 但迭代法并不能保证在合理的时间内收敛，对于条 件数很大的病态问题也可能是不收敛的。将直接法 和迭代法结合在一起的区域分解法 [5-8]能够很好地 利用两者的优点提高并行效率。它首先通过三角分 解进行缩聚消去各子区域内部自由度以降低界面方 程的阶数，然后利用迭代法求解系统界面方程。界 面方程的规模比系统整体方程的规模小了很多，因 而更容易求解。然而对于大规模问题采用传统区域 分解法求解时，随着子区域数目的增多界面方程的
2
分层并行计算方法
分层并行计算方法通过两级分区两次缩聚在分 布式数据存储的基础上既将节点内通信与节点间通 信分离有效提高了通信效率，又进一步降低了界面 方程的规模，加快了其迭代收敛速度。 2.1 两级分区 在两级分区里， 有限元网格首先被剖分为 M 个 1 级子区域，然后每个 1 级子区域再进一步被剖分 为 N 个 2 级子区域。为与 MPP 分布式存储体系结 构相适应，M 应为并行计算每次启动节点机总数， N 为单个节点机内处理器总核数。并行计算时，每 个 1 级子区域被分配到 1 个节点机上，同时派生于 该 1 级子区域的所有 2 级子域分别被分配到同一节 点机的不同核上。如图 1 所示，若 M 为 2，N 为 4，
图 1 两级分区 Fig.1 Two-level partitioning
将所有子区域的界面方程联立求解，即可得到 本文采用迭代法中广泛 各子区域边界节点位移 UB。 使用的并行 PCG 算法[14-15]求解界面方程。 各子区域 内部位移由下式回代求解
第 卷
第 期
年 月
力 学 学 报 Chinese Journal of Theoretical and Applied Mechanics
Vol. No. ，
结构力学有限元分层并行计算方法 1）
苗新强*,+，金先龙*,+,2)，丁峻宏**
* +
（上海交通大学 机械系统与振动国家重点实验室，上海 200240）
在传统区域分解法中，有限元网格仅被剖分 1 次，然后在生成的若干个子区域上应用 Schur 补过 程完成求解[11-13]。按照先内部自由度后边界自由度 的编号原则，每个子区域的系统方程可表示为:
K II K BI K IB U I PI K BB U B PB
图 2 两次缩聚 Fig.2 Twice condensation
各子区域的数据信息均可通过多文件流存储在 相应的节点机上实现存储局部化，这就有效提高数
据的内存访问速率。另外相对于传统区域分解法， 分层并行计算方法通过两次缩聚进一步降低了界面 方程的规模，加快了其迭代收敛速度。 2.3 三层并行计算的实现 由于 MPP 节点内通信速率远高于节点间通信 速率，故要提高通信效率就要实现节点内通信与节 点间通信的分离以减少进程间通信和同步的开销。 即要将大量局部通信限制在各节点内部，并最大限 度减少不同节点间的全局通信。而分层并行计算方 法刚好满足这些条件，如图 3 所示它在两级分区两 次缩聚的基础上实现了计算过程的三层并行。 第一层并行，每个进程独立负责 1 个 2 级子区 域的处理过程，各进程同时运行，进程间没有任何 数据交互。这些处理过程包括：读取子区域模型数 据，形成子区域系统方程，子区域缩聚，回代内部 自由度以及计算应变、应力。由于这些过程都是基 于单元的操作，各进程可根据相应 2 级子区域模型 信息单独运算，所以进程间不需要任何数据交互。 为节省内存空间和减少计算量，各 2 级子区域 局部总体刚度矩阵均采用 skyline 格式存储。缩聚 时，通过在每个进程上单独执行串行改进乔列斯基 分解过程实现。计算时只对列高之内的元素进行操 作，从而大幅度减少了计算量。
2.2 两次缩聚 两次缩聚通过先后在 2 级子区域和 1 级子区域 上应用 Schur 补方法实现缩聚过程。如图 2 所示， 首先形成每个 2 级子区域的系统方程，并经过缩聚 消去其内部自由度；然后通过组集同一节点机内所 有 2 级子区域界面方程，形成相应 1 级子区域系统 方程； 再通过缩聚消去每个 1 级子区域内部自由度， 得到其仅与边界自由度相关的界面方程。
]{U } {P } [K B
（2 ） （3） （4 ）
式中有效刚度矩阵
] [ K ] [ K ][ K ]1[ K ] [K BB BI II IB
有效载荷向量
} {P } [ K ][ K ]1{P } {P B BI II I