面向对象有限元并行计算框架PANDA求解器的服务构件化设计与集成

摘 要：针对很多结构有限元程序需要使用数值方法进行系统方程组求解的共性特点，提出对数值求解部分进行构件化设计，形成求解器服务，用于面向对象有限元并行计算框架PANDA为不同应用类型的有限元程序开发提供数值算法.将美国Livermore国家实验室开发的数值求解器库HYPRE作为数值求解构件集成到PANDA框架中，利用其提供的数值解法进行结构静力有限元并行求解，获得较理想的结果.目前，PANDA框架已经能够为应用程序的开发提供较丰富的数值解法方面的选择.

关键词：并行计算框架； PANDA； 有限元； 软件开发； 构件化软件设计； 求解器服务； 构件集成

中图分类号：TP311； O242 文献标志码：A

Component-based design and integration of solver service of object-oriented finite element parallel computation

framework PANDA JI Xiaoyu, HAO Zhiming, MO Jun, WANG Keying

(Institute of Systems Engineering, China Academy of Engineering Physics, Mianyang 621900, Sichuan, China)

Abstract： According to the common feature that many structural finite element programs need to use numerical methods to solve system formula, component-based design for numerical solving is proposed to build solver service, which is used for object-oriented finite element parallel computation framework PANDA to provide numerical algorithms for different types of finite element application development. The numerical solver library HYPRE as a numerical solving component is integrated into PANDA framework, by which the parallel finite element solving on structure static is performed and disired results are achieved. At present, PANDA framework can provide abundant numerical solving methods for finite element applications.

Key words： parallel computation framework; PANDA; finite element; software development; component-based software design; solver service; component integration

0 引 言

重大工程结构与装备在服役和极端环境下高保真、高精度模拟的需求越来越强烈，越来越需要开发大型集成化CAE分析软件.要在统一的软件平台上解决几何、材料和状态非线性、不同物理场耦合、不同学科、跨越时间和空间多尺度系统分析等多种工程计算领域的前沿问题是个庞大而复杂的系统工程，构建这样的CAE分析软件系统需充分利用已有的成熟理论、算法和代码，并进行有效的管理.传统的结构化程序设计方法已很难适应这种涉及全物理场、高保真、高精度模拟所需的复杂软件系统的开发，引入面向对象的设计思想，基于构件化软件开发（Component-Based Software Development, CBSD）方法开发大型软件已成为不可逆转的趋势.^［1-4］

在过去的10多年中，利用面向对象的程序设计方法在不同领域的工程分析中已取得丰富的成果^［5-10］.经过近60年的发展，有限元法已成为工程结构分析领域的主要方法，本文基于CBSD开发思想，搭建并行自适应非线性分析框架(以下称PANDA框架)^［11］，为工程结构有限元程序开发提供基础支撑.有限元计算可简单地分为前处理、计算和后处理等3部分.力学工作者关心的计算部分可分为单元刚度矩阵与载荷向量、组集方程组、求解方程组、计算应变与应力等几部分.实践证明：线性系统求解是整个计算过程最耗时的部分，解决该问题应首先关注数值求解算法.^［11］与此同时，不同应用类型的有限元程序往往需要使用不同的求解器库进行系统方程的求解工作，如并行稀疏直接求解需使用SuperLU和SPOOLES等；并行迭代求解需使用Aztec，PETSc和HYPRE等；特征值问题求解需使用SLEPc和PARPACK等；多重网格求解需使用UG和HYPRE等.这些开放源代码的数值求解器库都针对不同应用类型进行设计，并经实践验证其有效性和并行效率.如果每编制一种类型的应用程序都研究一种数值方法或研究一种数值求解器库并完成对接工作，会花费较多的精力.将数值求解部分进行构件化设计，形成求解器服务沉淀到PANDA框架中，在其中为众多的求解器库定义统一的接口，在应用程序开发过程中只需根据框架的要求调用求解函数即可实现并行求解.这样可使应用程序开发人员专注于物理力学算法的研究，不必太过关注数值算法，从而缩短程序的开发周期.

本文介绍PANDA框架求解器服务构件化设计与实现，为大规模并行有限元程序的数值求解提供丰富的数值求解算法.通过构件化设计并定义统一接口，可方便集成经实践检验、广受欢迎的开源数值求解器库作为求解器构件.在并行有限元程序开发过程中，可借用这些成熟的并行数值分析软件进行有限元方程组的数值求解以缩短开发周期.作为求解器服务构件化设计的集成示例，介绍集成美国Livermore国家实验室开发的著名并行求解器库HYPRE的基本思路和方法，为基于PANDA框架开发的应用程序提供数值支持.

1 构件化设计

PANDA框架的层次结构和分层架构参见文献［12］，其充分借鉴美国Sandia国家实验室开发的SIERRA框架和开源代码Tahoe，遵循面向对象、层次化和模块化的设计思想，通过高度抽象与提炼，形成一系列基础服务构件，从而为结构有限元并行程序开发提供支持.求解器服务位于PANDA框架的数值共性层.该层包括数据I/O服务、用户输入解析、成熟的计算方法和数值算法等，在科学与工程计算中具有较强的普适性，可直接支持并行有限元程序的研制.

大型稀疏线性方程组并行数值求解是实现大规模并行有限元分析的关键技术之一，目前尚未有一种数值方法或一个求解器库能解决所有类型的数值求解问题，不同类型的有限元分析需在不同的求解器库和数值算法中进行选择.对数值求解部分进行构件化设计，为不同数值算法或求解器库定义统一的接口，将数值方法沉淀到框架的数值共性层，便于应用程序选择、组集适合其应用类型的数值算法.基于此便利性，应用程序开发人员可将主要精力用于物理力学问题本身的研究，只需根据应用类型的数值求解需求选择相应的求解器库和数值求解算法.PANDA框架求解器服务构件化设计架构见图1.

图 1 求解器服务构件化设计架构

Fig.1 Component-based design architechture of

solver services

基于PANDA框架开发的应用程序通过有限元管理器进行有限元运行过程的管理：在求解部分通过求解器管理器针对不同问题类型实施具体的求解策略，对有限元计算流程（如问题类型、求解驱动、结果与异常处理等）进行抽象的定义，具体实现由其派生类LinearSolver（线性求解）和NLSolver（非线性求解）等进行定义（见图2），主要负责向应用程序提供数值求解方面的支持，并通过句柄操作GlobalMatrix对象，实现PANDA框架数据对象到GlobalMatrix数据对象的传递.纯虚基类GlobalMatrix是面向第三方求解器库的抽象接口，是数据对象到第三方求解器库转换的抽象实现，针对每个求解器库需要一个具体的实现.目前，已经在PANDA框架中集成Aztec，SPOOLES，PHG，HYPRE，SuperLU，Epetra和PSPASES等著名的开源求解器库，并将一些数值算法作为其构件沉淀到求解器服务中(如CCS和CCNS等).

图 2 求解器服务的构成

Fig.2 Structure of solver services

2 构件集成

PANDA框架为并行有限元程序开发提供基础支撑，开发人员在框架的基础上定制、组装应用程序.构件化领域框架的可定制性和可扩展性体现在框架的核心结构保持稳定的前提下，在框架中提供足够多的反映应用领域易变性的热点来满足不同应用类型的特殊要求.^［3］求解器服务是PANDA框架的重要扩展点之一，向结构有限元应用程序提供选择、组集不同数值解法方面的支持.

由于基于继承思想和多态性机制实现求解器服务的构件化设计，尤其是使用多态性机制的虚基类GlobalMatrix为程序开发人员和用户提供方便的抽象接口函数，只需根据该类的规范编写PANDA框架与数值求解器库的接口代码就可完成在框架中集成第三方求解器库的工作.本文给出集成HYPRE的方法，作为在PANDA框架中集成求解器库成为其求解器服务构件的示例.HYPRE为应用人员提供功能强大的预条件子，同时提供一组迭代求解法(如PCG等).将其集成到框架中作为求解器构件的步骤如下：

(1)从虚基类GlobalMatrix中派生出一个名为HYPREMatrix的类，它继承类GlobalMatrix的属性和方法.针对不同的求解器库需声明不同的私有成员以适应求解器库的要求，如针对HYPRE需声明

HYPRE_IJMatrix A;/* HYPRE Matrix */

HYPRE_ParCSRMatrix parcsr_A;

HYPRE_IJVector b; /* RHS */

HYPRE_ParVector par_b;

HYPRE_IJVector x; /* Solution vector */

HYPRE_ParVector par_x;

HYPRE_Solver solver, precond;

其中：parcsr_A，par_b和par_x是HYPRE求解器可直接调用的数据对象，通过HYPRE提供的函数分别从A，b和x转化而来.要使Solver能识别新的求解器库，需在函数Solver::NewSub中添加参数HYPRE_Matrix，在函数Solver::SetGlobal Matrix中设置当求解器参数选定为HYPRE_Matrix时需执行的操作.

(2)为派生类HYPREMatrix添加属性和操作.虚基类GlobalMatrix抽象求解器的共性，而不同求解器在数据存储、数值算法等方面往往不同，需重定义GlobalMatrix中的虚函数以添加具体的属性和操作，这部分代码的编写工作是集成HYPRE的重点.HYPRE Matrix.cpp给出成员函数的具体实现，其中最重要的是组装函数和求解驱动函数.为实现有限元程序分析中总体刚度矩阵的并行组装和方程求解，首先需使用初始化函数Initialize传递3个形式参数：总的方程数tot_num_eq，本地管理的方程数loc_num_eq和起始方程号star_eq.Initialize函数以ilower=start_eq-1和iupper=ilower+loc_num_eq-1这2个参数记录当前管理的本地起止方程号，以参数ilower和iupper为基础创建A，b和x.组装函数完成PANDA生成的单元刚度矩阵向HYPRE的传递，并按照HYPRE软件要求的并行稀疏行压缩格式进行存储.它需要从框架中获取elMat（单刚）和eqnos（全局方程编号）2个数组.当矩阵对称时，调用函数elMat.CopySymmetric先将矩阵补全，然后统计矩阵的行和列信息，并创建临时数组记录矩阵元素，以方便调用HYPRE求解器库中的函数HYPRE_IJMatrixAddToValues，将单刚组装成局部总体刚度矩阵.求解驱动函数向HYPRE传递线性系统的载荷向量，并完成线性系统求解、返回求解结果等工作.该函数需调用1个双精度数组result，此时存放的是线性系统的右端项.使用函数HYPRE_IJVectorSetValues将result中的值赋给右端项b，求解线性系统获得解向量x.将x的值再赋给result，释放临时数组指针x_values.这样返回到框架和应用程序的result存放的是x.

(3)设置与求解器相关的编译制导文件.PANDA框架采用编译制导文件makefile进行分级管理.需在GlobalMatrixT所在的目录/…/globalmatrix下创建目录HYPRE，将派生类HYPREMatrix的头文件和实现文件放入该目录；在/…/globalmatrix/makefile文件中指定编译子目录HYPRE；另外，还要在PANDA框架主目录下的编译制导文件makefile中加入HYPRE的预编译指令、头文件目录、链接参数和库文件等信息.

按照上述步骤即可完成集成HYPRE的工作，可按照Linux操作系统的要求编译和安装PANDA框架，形成可应用的PANDA框架.数值求解的2个重要步骤（组集系统总体方程组和方程组并行求解）都已在步骤(2)中实现，集成其他求解器库只需根据各自库函数实施同样的步骤即可.

3 数值算例

基于框架的应用程序开发工作将变得快捷、高效，只需编写程序驱动代码组装一个针对不同应用的程序即可.以结构非线性静力有限元分析为例，为Solver类定义一个派生类(如NLSolver)是重要的工作之一，在其中定义问题求解流程，需调用的模块及其接口函数均由框架提供；另一个重要工作是编写主驱动程序，根据框架提供的规范接口，编写驱动代码实现模型与数据的输入/输出、区域分解和通信初始化、问题求解等过程.

以某简支受拉蜂窝梁为数值研究对象(见图3)，其横截面的外轮廓为正方形，长宽比为10∶1.

图 3 某简支受拉蜂窝梁模型

Fig.3 Model of a freely supported tensioned ceu beam

使用八节点六面体单元对其进行离散，根据单元尺寸不同分为2种规模，分别为286万自由度和1 080万自由度.在图3所示蜂窝梁的左端施加简支约束，在右端表面施加10 MPa的拉力.取模型中上层圆孔的正上方节点为特征点（共计20个），使用基于PANDA的静力分析程序和ANSYS进行对比分析，其位移和第一主应力的计算结果见表1和图4.其中，基于PANDA框架开发的静力有限元分析程序采用HYPRE提供的PCG方法进行线性方程组的求解.可知，2种程序在模型中20个特征点上位移结果的平均误差约为1.6%，第一主应力的平均误差约为2.6%.

基于PANDA框架开发的应用程序在可移植性方面的测试工作参见文献［13］，在中国科学院系统科学与数学研究院、中国工程物理研究院仿真中心和中国工程物理研究院总体工程研究所的不同计算平台上均成功运行，计算效率良好.将本算例的286万自由度模型和1 080万自由度模型放在中国工程物理研究院总体工程研究所的曙光TC2600刀片服务器上进行测试，结果见表2，加速比曲线见图5.其中，1 080万自由度模型不能实现4进程以下的计算，因此，以4进程计算时间为基准进行加速比计算，可称为伪加速比.

图 5 2种规模的加速比曲线

Fig.5 Speedup curves of two different models2种规模的加速比测试结果表明：基于PANDA框架开发的静力有限元分析程序具备较好的加速比性能.在进程数少于16个时，2种规模的加速比曲线都接近斜率为1的基准线，1 080万自由度模型在32个进程的测试中仍然获得较好的加速比，286万自由度模型在32个进程的测试中加速比略有下降.

4 结束语

构件化领域框架极大地方便软件集成，利用CBSD思想已成功建立PANDA框架的雏形.求解器构件化设计为PANDA框架提供求解器服务方面的支撑，已成为框架的一个重要扩展点.通过多方集成，已能为基于PANDA框架开发的应用程序提供较为丰富的数值解法方面的选择.

PANDA框架为工程数值模拟程序开发工作提供高性能数据结构，封装并行算法和并行数值求解技术，使应用程序开发人员和用户可专注于物理力学本质的研究.在不同计算平台上进行多种规模的数值计算表明，基于PANDA框架开发的静力有限元程序具备较好的可移植性、可扩展性和并行效率.

致谢：在求解器构件化设计和数值算法的测试等方面得到中国科学院数学与系统科学研究院张林波研究员的悉心指导和热情帮助，张林波研究员还向本文作者提供PHG数值求解器库，并为相关测试提供硬件资源.

参考文献：

［1］ WHITEHEAD K. Component-based development: principles and planning for business systems［M］. London: Addison-Wesley, 2002: 2-16.

［2］ SZYPERSKI C, GRUNTZ D, MURER S. Component software: beyond object-oriented programming［M］. 2nd ed. London: Addison-Wesley, 2002: 105-178.

［3］ 肖刚, 徐俊, 张元鸣, 等. 基于软件配置模型的构件化领域框架研究［J］. 计算机应用与软件, 2008, 25(5): 90-92.

 XIAO Gang, XU Jun, ZHANG Yuanming, et al. Research on compenent-based domain-specific framework on software configration model basis［J］. Comput Applications & Software, 2008, 25(5): 90-92.

［4］ 郭秋萍, 焦允, 王全兰. 大规模系统构架建模及其开发技术［M］. 北京: 北京航空航天大学出版社, 2008: 4-20.

［5］ AMUNDSON J F, DECHOW D, MCINNES L, et al. Multiscale, multiphysics beam dynamics framework design and applications［J］. J Phys: Conf Ser, 2008, 125(1): 1-6.

［6］ PATZK B, BITTNAR Z. Design of object-oriented finite element code［J］. Adv Eng Software, 2001, 32(10-11): 759-767.

［7］ FORDE B W R, FOSCHI R O, STIEMER S F. Object-oriented finite element analysis［J］. Computers & Structures, 1990, 34(3): 355-374.

［8］ ZABARAS N, SRIKANTH A. Using objects to model finite deformation plasticity［J］. Eng Computers, 1999, 15(1): 37-60.

［9］ ZABARAS N, SRIKANTH A. An object-oriented programming approach to the Lagrangian FEM analysis of large inelastic deformations and metal-forming processes［J］. Int J Numer Meth Eng, 1999, 45(4): 399-445.

［10］ SILVA E J, MESQUITA R C, SALDAHHA R， et al. An object-oriented finite-element program for electro-magnetic field computation［J］. IEEE Trans Magn, 1994, 30(5): 3618-3621.

［11］ 张向, 许晶月, 沈启彧, 等. 面向对象的有限元程序设计［J］. 计算力学学报, 1999, 16(2): 216-225.

 ZHANG Xiang, XU Jingyue, SHEN Qiyu, et al. Object-oriented finite element programming［J］. Chin J Comput Mech, 1999, 16(2): 216-225.

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

 SHI Guangmei, HE Yingbo, WU Ruian, et al. Object-oriented finite element parallel computing framework PANDA［J］. Comput Aided Eng, 2010, 19(4): 8-14.

［13］ 范宣华, 吴瑞安, 郝志明, 等. 基于Tahoe框架的某夹具并行计算［J］. 计算物理, 2009, 26(5): 699-702.

 FAN Xuanhua, WU Ruian, HAO Zhiming, et al. Parallel computing of clamp structure in Tahoe frame［J］. Chin J Comput Phys, 2009, 26(5): 699-702.

(编辑 陈锋杰)