面向对象有限元并行计算框架PANDA的并行机制

摘 要：为实现面向对象有限元并行计算框架PANDA对高性能计算的支持，分别从并行计算流程、区域分解、分区信息和通信封装等部分设计PANDA框架在并行计算方面的数据结构.在计算流程中建立区域分解和并行求解器的配合协作方式，进而描述进行区域分割的3种网格剖分方法；对分区边界单元和节点信息的组织以及对并行通信操作的封装使复杂的并行通信调用简单、易行. 用PANDA框架开发有限元程序时，只需以串行方式编写代码就可得到并行计算服务.对某夹具的算例表明PANDA框架并行计算实现的正确性和可靠性.

关键词：并行计算框架； PANDA； 有限元； 软件开发； 并行机制

中图分类号：TP319； O246 文献标志码：A

Parallel strategy of object-oriented finite element parallel computation framework PANDA

LI Yufeng, ZHANG Yalin

(Institute of Computer Application, China Academy of Engineering Physics, Mianyang 621900， Sichuan, China)

Abstract： To support high performance computation in object-oriented finite element parallel computation framework PANDA, the data structures which are relative to parallel computation are designed in the parts such as computation flow, domain decomposition, partition information, communication encapsulation and so on. The collaborative approach between domain decomposition and parallel solver is established in computation flow, and then three meshing methods are described for domain partition. Due to the good organization of boundary elements and node information of partition with parallel communication encapsulation, the complicated communication invoking becomes simple and feasible. While developing finite element programs, parallel computation service can be obtained by programming in serial codes. An example about clamp indicates that the implementation on parallel computation of PANDA is correct and reliable.

Key words： parallel computation framework; PANDA; finite element; software development; parallel strategy

0 引 言

有限元法是当今主流的数值计算方法，被广泛应用于工程领域的各个方面.大量高性能计算机的出现，特别是集群系统的广泛应用，使大规模复杂数值模拟得以实现，但必须采用并行技术进行程序开发.有限元并行计算框架开发的目的是使在复杂网格上进行的大型数值模拟程序开发变得更加快捷、高效.近年来，面向对象有限元技术［1-4］发展迅速，而开发框架却为数不多，且功能模块有较强耦合［5］.为此，中国工程物理研究院总体工程研究所和计算机应用研究所联合开发有限元并行计算框架PANDA［6］.框架中各组件松散耦合，可适应各种线性和非线性有限元法，不仅提供单元公式和本构模型的添加接口，也提供各种复杂边界条件的实现方法.

基于消息传递在分布存储模型的并行机上开发并行程序，能显著控制数据分配和通信，这种控制又促进在大规模分布存储并行机上的高性能编程.［7］ 针对目前高性能集群的广泛应用，采用MPI(Message Passing Interface)以分布存储方式实现有限元并行是正确的选择.

面向对象技术对并行编程过程有重要影响(如对数据和函数的封装、继承和多态性)，许多基于类等的编程思想都来源于串行编程领域，但也开始在并行编程中使用，尤其是恰当的抽象可隐藏那些使编程复杂的并行性.［7］ 以面向对象技术为基础设计的PANDA框架正体现该优点.

1 框架的数据结构

框架的主要数据结构有网格、场变量、单元组和通信数据管理等.

网格是整个模拟计算区域的离散，支持线段、三角形、四边形、四面体和六面体等几何形状.一个计算区域可包含不同形状的网格子集.所有网格单元被划分为若干单元集，一个单元属于且仅属于一个单元集.

场变量即自由度及其微分，框架以场的形式组织变量.任意一个自由度变量都归属于某个场，场通常定义在网格节点或单元上.根据时间积分算法，变量可能存储的数量不一，即场可能含有d阶微分，均要开设存储.所有的场均由节点管理器管理.

单元组是PANDA框架的算法核心类，实现一类有限元公式，在特定单元形状、形函数（含单元空间积分方式）、计算参考系（Euler或Lagrange或ALE）和物理固有非线性（如接触）等情况下形成一种特殊的计算模式，即如何形成系统方程的单个单元贡献.每个单元组都会关联一个场，同时关联若干个单元集，负责遍历所施加的单元，进行单元矩阵和右端项的计算，并提供组装方程组的接口.

通信数据管理基于网格剖分，所有单元集作为一个整体进行剖分.在采用镜像节点后，通信只发生在节点信息上.

2 并行的实现策略

PANDA框架基于MPI实现并行，具体策略为：首先，基于几何网格的区域分解，将子域映射到各进程；然后，每个进程在其子域上计算单元矩阵和单元右端项并提供组装接口，借助并行线性代数求解器实现大规模矩阵的分布式存储、组装接口以及方程求解，求解后，交换场数据边界信息；最后，进程在局部进行时间积分更新.所有的有限元相关计算都是分布式并行的，而大规模方程求解的主要并行任务交由线性求解器库负责.目前，已在PANDA框架中集成Aztec，SPOOLES，SuperLU，HYPRE［8］，PETSc，Trilinos和PHG（由中国科学院科学与工程计算国家重点实验室开发）等著名求解器库.

3 框架的并行计算流程

框架并行设计的关键是数据变量的分布式存储和场数据时间更新算法的并行化.框架默认对所有的场变量基于网格区域分解进行并行通信.框架并行计算流程见图1.

图 1 框架并行计算流程

Fig.1 Parallel computation process of framework

由图1可知，并行计算首先需进行区域分解，然后将局部模型信息读入框架，再进行场变量的编号(局部编号在进程内部确定，全局编号则通过全局通信确定). 基于编号的局部和全局映射关系，每个子域将遍历其所存单元，进行单元矩阵和右端残力的计算和组装，并组装到在初始阶段就设置好的并行求解器中.

有限元计算对每个自由度都形成一个代数方程，方程的编号很重要，框架采用的编号方式为：局部编号按节点进行，节点内按场顺序进行，场内按自由度顺序进行；全局编号按进程号排序，并保持局部编号的连续性.框架的并行是基于节点的剖分，所以按节点编号的方式对于并行处理十分方便.得到自由度编号和活动方程的编号(包括全局编号和局部编号对应关系)，就可将从单元算法得到的单元矩阵和残量组装得到整体方程的矩阵和残量.通过并行求解器得解的本地部分.

并行求解器主要完成A(u)=f的求解.作为一个独立的外部接口模块，其输入包含向量f（或者还有初始值u）的分布式存储、算子A的操作对u的影响方式.其中，可能涉及稀疏矩阵的分布式存储、A(u)的分布式计算和规约；输出为解向量u的分布式最终结果，其在误差允许范围内近似满足方程或报告不收敛状态.

求得分布式存储后，下一步对场的时间积分算法更新.为计算内力，需从相邻分区获得外部影像节点场信息，需用到分区管理和通信管理.PANDA框架将分区和通信进行封装，使上层只需1~2个函数的调用即可完成外部信息的交换，设计灵活、巧妙.

在进行场信息的更新以及内力的计算后，就可推进时间步进行下一次的循环迭代：根据算法选择进入局部单元遍历计算和组装，或只进行右端残量的遍历计算和组装，再到并行求解器求解，最后到场的更新和内力计算，循环往复.

3.1 区域分解

框架采用网格区域剖分的方法实现计算任务的分解.目前，框架网格分解以串行实现，且以节点分解为基础，进而对单元和通信作出安排.

PANDA框架提供的分解方法有3种：图分解方式、索引分解方式和空间分解方式.图分解方式先将网格组织成图的形式，然后按照图的分解方法进行分解，最后将结果转变为网格单元和节点的分解；索引分解方式是按照模型读入的节点顺序依次分配给p个进程；空间分解方式是按照空间坐标进行箱子分解，每个箱子的尺寸都相同，先形成很多小箱子，然后按照箱子的邻接关系，借助图剖分将箱子分成p个组.由于有限元计算大多是以单元为单位进行的计算密集型作业，对于复杂几何构型的离散化，图剖分方式似乎是最好的选择.

将网格组成图的方式由单元连接关系和节点场关联信息共同确定.场关联信息存在于一些复杂相依的约束条件(如边界条件或接触条件)中.目前，PANDA框架调用METIS［9］开源软件包进行多水平图剖分，也可采用其他方式分解.由图剖分得到网格节点的划分情况，然后需确定如何对有限元计算形成分区信息，用这种信息指导数据通信.

3.2 分区信息

框架采用类PartitionT管理所有的分区信息.首先明确的是哪些节点或单元信息需进行通信，哪些在本地局部范围不需要与其他分区通信.

从区域分解得到节点-分区号映射后，结合图的表示进行节点的类别划分，再结合单元连接关系进行单元的类别划分.根据图和单元连接关系可将本地存储的节点和单元进行分类，形成3类节点和2类单元，见图2.图中网格被虚线分为2个子域，节点和单元的分类基于右边的子域.内部节点指没有与其他分区节点有连接关系的本分区所拥有的节点；边界节点指至少与某个其他分区节点有连接关系的本分区所拥有的节点；外部节点指与本分区边界节点有连接关系的其他分区的节点.内部单元特征是其所有节点均为内部节点，边界单元特征是至少一个节点为边界节点或外部节点.

图 2 并行计算中分区的单元和节点分类

Fig.2 [WB]Node and element classification in a partition

of parallel computation

在节点和单元分类后，就可以明确消息的来源和目的.为方便管理，采用分区对象处理这些信息.调用分区函数，得到分区信息的初始化，其中主要函数有4个：（1）Set()设置分区总数、ID号和本地3类节点，接收信息的分区号和节点号等；（2）SetOutgoing()设置发送信息的节点号；（3）InitElement Blocks()设置2类单元数组，单元映射和逆映射数组的主维数为单元块数目；（4）SetElements()按照单元连接关系将单元划分为2类，设置各单元数组及映射数组.

通过这些函数可对分区对象进行必要的设置，另一种方法是读取分区文件.在设置完成后就可使用PartitionT对象，得到的重要信息包含本分区所有节点号、所有单元号、编号局部到全局的映射及其逆映射等.得到分区信息后，可根据需要写几何输出文件.设总分区数目为M，本分区号为m(0≤m

在框架中有限元管理器读入分区文件，构造本地分区PartitionT对象.通信管理器设置其成员fPartition.在并行模拟中，可通过PartitionT实现对通信消息的掌握，通过层次封装，上层只能得到通信管理器和有限元管理器掌握的节点的处理器位置信息.

3.3 通信封装

边界节点信息的传递是整个通信的核心任务.在图剖分的区域分解中，通过分区类PartitionT的NodesOut和NodesIn函数可获得本分区需互换信息的节点号.

节点数据的分布在框架中影响以下重要的计算步骤：大型稀疏矩阵的分布式存储和场信息的更新.求解器的初始化函数设置方程的全局方程总数、局部方程总数和局部方程起始号；在单元贡献集成时需判断方程号是否在本地范围内，不在则不集成；场信息的更新由节点管理器负责，在其成员函数的实现中，调用通信管理器的AllGather函数（已经过封装，在图剖分情况下属于点到点的非阻塞通信）.

图 3 通信的封装层次

Fig.3 Hierarchy of communication encapsulation分区信息通过通信的封装影响系统方程的并行求解和场信息的更新.通信在框架中的重点类有CommunicatorT，MessageT和CommManagerT，封装层次见图3，上层的类使用下层的类提供服务，应用程序可使用以下各层，但最方便的是使用通信管理类，其与分区之间的协作见图4.图 4 通信管理类与分区之间的协作

Fig.4 Relationship between communication classes and

partition

通信器类CommunicatorT是对标准MPI函数的封装.涉及的全局通信有全规约、收集、全收集、广播和同步等；局部通信有非阻塞和阻塞的发送与接收.开发者在需要其他数据通信时，可使用这种封装的通信函数进行.该方法体现模块性，提高错误处理的能力，但用户使用时需理解具体的通信细节.

MessageT是对各种并行消息通信模式的封装，有PointToPointT，AllGatherT和CartesianShiftT等3个派生类，分别完成点对点、全收集和笛卡儿网格移位通信，对应于图剖分、索引分解和空间分解等3种区域分解方式的通信.

PointToPointT提供点对点的阻塞、非阻塞通信.fPartition成员提供剖分信息，其他成员包括非阻塞的接收缓冲区、发送缓冲区、接收和发送请求.成员函数Initialize()按分区的进出节点信息检查或分配接收和发送缓冲区，接收发送请求.AllGatherT()依据分区进行非阻塞的通信，先投递接收信息，然后对要发送信息的节点信息进行发送投递，接着等待处理接收到的消息并组装到欲通信的数组里，最后等待发送的完成.

AllGatherT类提供全收集的通信功能，经过通信，每个进程都得到相同的全局数据场，主要服务于基于索引的区域分解方式.成员函数Initialize()设置本进程发送数据的长度以及确定本进程数据在目标数组中的位移.AllGatherT()执行全收集操作.应用于索引分解的全收集通信只是计算功能的并行，在数据空间上并没有分布存储，每个进程都存储全局数据.

CartesianShiftT提供笛卡儿网格上的移位通信，主要成员为构造函数和AllGatherT().针对空间分解方式，在空间的方向和位移上交换数据.该类的主要成员有周围进程号，通信移位模式以及发送、接收信息的节点号.在实际应用中，可利用框架提供的3种消息模式进行通信，只需创建对象，初始化数据分布，然后调用AllGatherT()就可实现.

CommManagerT是框架中对并行操作的顶层封装.它维护分区信息、管理消息列表，作为框架的重要组成部分负责节点场信息的通信.使用CommManagerT进行数值模拟物理量的通信非常方便：其将所有的通信消息管理起来，开发者只需提供欲通信的数组，在进行注册后就可通过简单的调用完成通信.

CommManagerT的主要接口函数有Init_AllGatherT()和AllGatherT().前者根据区域分解方式初始化相应的消息通信方式，将分配好的消息追加存储到成员fCommunications数组；后者调用相应消息的通信函数（对于图剖分则调用点对点消息的AllGatherT()函数）.

典型的CommManagerT使用方式为

nArray2DTvalues(num_nodes,num_dofs);

int id = fCommManagerT.Init_AllGatherT(values);

… /* values assingment locally */

fCommManagerT.AllGatherT(id, values);

fCommManagerT.Clear_AllGatherT(id);

其中：nArray2DT为框架工具箱中二维数值数组模板类.使用CommManagerT很方便，隐藏具体的分区通信方式，易于上层应用，并行通信在框架中逐层抽象封装，应用程序可选择使用不同封装层次的类.

4 算 例

用基于PANDA框架开发的有限元模拟程序对某夹具进行计算，采用开源并行预条件库HYPRE求解代数方程组.模型自由度约为10万，得到的位移云图见图5.

图 5 某夹具位移云图

Fig.5 Displacement clound picture of a clamp

由图5可知，并行计算的结果与串行计算结果一致，并行计算中单计算节点的时间和内存消耗见表1，测试平台为某Linux集群，单计算节点处理器为Intel(R) Xeon(TM) CPU 3.40 GHz，内存4 GB.可知，基本呈现出线性加速比.

5 结 论

PANDA框架可实现区域分解、分区管理、通信封装、场的通信等的无缝集成、相互协作，完成工程计算的并行通信任务，屏蔽具体通信底层实现的复杂性，与外部线性代数求解器的接口设计巧妙，可配置多种求解器.基于PANDA框架开发有限元程序，开发者只需以串行的方式编写与物理相关的代码，就可得到并行的计算服务.在PANDA框架中利用封装好的CommManagerT进行自定义变量的通信也十分方便，甚至只需1~2个函数的调用就能实现多计算节点间复杂拓扑关系的通信.

参考文献：

［1］ BOSE A, CAREY G F. A class of data structures and object-oriented implementation for finite element methods on distributed memory systems［J］. Comput Methods Appl Mech Eng, 1999, 171(1-2): 109-121.

［2］ PATZAK B, BITTNAR Z. Design of object-oriented finite element code［J］. Adv Eng Software, 2001, 32(10-11): 759-767.

［3］ PANTALI O, CAPERAA S, RAKOTOMALALA R. Development of an object-oriented finite element program: application to metal-forming and impact simulations［J］. J Comput Appl Math, 2004, 168(1-2): 341-351.

［4］ COMMEND S, ZIMMERMANN T. Object-oriented nonlinear finite element programming: a Primer［J］. Adv Eng Software, 2001, 32(8): 611-628.

［5］ YU Lichao, KUMMAR A V. An object-oriented modular framework for implementing the finite element method［J］. Computers & Structures, 2001, 79(9): 919-928.

［6］ 史光梅, 何颖波, 吴瑞安, 等. 面向对象有限元并行计算框架PANDA［J］. 计算机辅助工程, 2010, 19(4): 8-14.

[KG*2]SHI Guangmei, HE Yinbo, WU Ruian, et al. Object-oriented finite element parallel computation framework PANDA［J］. Comput Aided Eng, 2010, 19(4): 8-14.

［7］ DONGARRA J, FOSTER I, FOX G, et al. 并行计算综论［M］. 莫则尧，陈军，曹小林, 等, 译. 北京： 电子工业出版社, 2005: 198-199.

［8］ FALGOUT R D, YANG U M. Hypre: a library of high performance preconditioners［C］//SLOOT P M A, TAN C J K, DONGARRA J J, et al. Proc ICCS 2002. Berlin: Spinger-Verlag, 2002: 632-641.

［9］ KARYPIS G, KUMAR V. METIS 4.0: unstructured graph partitioning and sparse matrix ordering system［R］. Minneapolis: Univ Minnesota, 1998.