基于PETSc的有限元高性能求解方法

高怀玉，张 峰，秦忠国

(1.河海大学力学与材料学院，江苏 南京 210098;2.南京市水利规划设计院有限责任公司工程分院，江苏南京 210006)

随着工程技术领域的不断开拓和发展，对大型复杂结构的分析要借助高阶数值分析模型［1］和大规模数学计算，以保证数值解的精度［2］。在普通的计算设备上采用传统的串行计算面临着因计算资源不足而无法计算以及计算时间过长、工作效率低的问题。使用高性能计算机进行并行有限元分析是解决上述难题的一个方案［3-4］。自从美国国家航空航天局NASA的Noor等［4］在1975年发表第1篇有限元并行计算的论文以来，有限元并行算法的研究越来越具体和广泛，适用于不同类型计算机的有限元并行算法也相继产生。20世纪90年代以来，单处理机的计算速度不断提高，并行计算机的体系结构趋于成熟，数据传输网络的标准化和传输速率的大幅度提升，使得并行计算机的研制周期从几年缩短至几个月，为研制并行计算机系统创造了有利条件［5-7］。高性能计算机的发展势头异常迅猛，拥有并行计算功能的有限元软件［8］也相继出现。李根国等［9］将NASTRAN在“神威I”超级计算机系统上进行移植和并行计算功能的二次开发，获得了成功。一些优秀的并行数值算法包也日益引起关注。借助并行算法包开发有限元并行程序无需考虑复杂的数据分布和通信，大大降低了有限元并行计算的难度和成本，缩短了开发周期［10］。

可移植可扩展科学计算工具箱PETSc(portable，extensible toolkit for scientific computation)由一套数据结构和程序组成，由美国Argonne国家实验室承担开发并获得成功，主要用于分布式存储环境下高效求解偏微分方程组及相关问题的求解［3］。PETSc在软件实现上采用ANSI C和消息传递标准接口MPI［11-12］。针对单处理机多核计算机，笔者利用PETSc提供的高层应用程序开发平台开发出一种线性并行求解方法，并应用编写的相应求解程序对求解器进行了测试。

1 PETSc的数据结构

PETSc的方程组求解器组件定义在各种数据结构组件之上，包括线性方程组求解器KSP、非线性方程组求解器SNES、时间依赖方程求解器TS。PETSc基于基础线性代数子程序库BLAS、线性代数包LAPACK、MPI等实现计算功能。PETSc为用户提供了丰富的Krylov子空间迭代方法和预条件子［13-14］，并具有错误检测、性能统计和图形打印等功能。基于PETSc提供的大量对象和求解库，用户可以灵活地开发自己的应用程序，还可随意添加和完善某些功能。

PETSc的数据结构包括向量Vec、矩阵Mat、网格管理对象、索引与排序、结构网格DA以及无结构网格的IS对象。用户在PETSc环境下编写串行或并行应用程序时，要按照固定的框架结构来编写，即对于一般的应用对象均需要声明、创建、设置、使用和销毁。

2 基于PETSc的有限元求解器概况

基于PETSc的有限元并行求解器(以下简称求解器)以C++类［15］的形式存在，求解器的类封装了一般线性有限元分析所必需的数据及操作，用户按照给定的格式提供数据，依次调用成员函数，即可完成一个有限元的并行分析，得到位移解。使用求解器时，可以通过改变变量nfreedom值来改变有限元分析的类别，如nfreedom值设为1时，对应的是一维有限元问题，nfreedom值为2和3则分别对应平面和空间有限元问题。自定义的单元结构体StructElement，对外提供一个单元刚度矩阵计算程序接口，用户可以自由地将自定义单元接入程序。求解器的类能方便地扩展单元类，具有一定的通用性和可扩展性，它继承了PETSc工具箱所提供的运行时选项(runtime option)、性能分析等功能，即在求解时可以改变求解程序的迭代方法、预条件方法等设置，这些功能为测试程序的效率提供了方便。

2.1 求解器结构

求解器在耗时较多的分析部分进行并行处理，具体流程如图1所示。当程序运行到插入节点力、形成荷载列阵和引入边界条件、修改方程组模块式时，进程0要向后层其他各进程广播这些信息。当程序运行到计算矩阵预分配参数模块时，各进程之间要保持相互通信。最后，各进程把各自计算出的结果返回到进程0，并保存结果，结束计算。并行分析部分包括预分配、形成总刚度矩阵及求解方程组，这些部分消耗的时间一般情况下占整个过程的90%以上，而且这一比例会随着计算规模的增大而增大，用如上并行方式求解线弹性问题，理论上能取得很好的效果。

2.2 求解器调用

求解器提供的可选设置选项说明如下:

a.SetElement(intType，intmaxnode)与 SetElement(int Type，proc－ESElementStiffer)是提供接口给用户接入自己的单元类型。其中参数Type是单元类型号，求解器默认类型为八节点六面体空间等参元;参数maxnode为该类型单元的节点数;自定义指针函数类型 proc－ElementStiffer则为指向该单元类型的单刚求解程序。通常这2条设置应该设置在初始化函数SolverInitialize()之后。

b.SetKSPType(char*ktype)与 SetPCType(char*ptype)语句用于设置子空间迭代法的类型，和PETSc工具箱的命名一样。子空间迭代法的类型 ktype可以是KSPCG，KSPGMRES，KSPTCQMR，KSPBCGS，KSPTFQMR，KSPCR，KSPLSOR，KSPBICG中的任何一种。预条件法的类型 ktype包括 PCJACOBI，PCBJACOBI，PCSOR，PCEISENSTAT。

图1 求解器流程Fig.1 Flow chart of solver

c.SetTolerances(double abstol，double rtol，double dtol，int maxits)用于设置迭代收敛参数，前3个参数是基于－范数 ( (∑)1/2)的剩余范数相对下降量rtol，剩余范数的绝对值atol及剩余中的相对增加量dtol。Maxits则为最大迭代次数。

d.SetTimeInfo()函数用于保存求解时间明细，包括各主要部分的时间消耗，以此来分析求解效率。SaveBinary()与SaveInfo()是将刚度矩阵和等效节点力列阵分别保存为文本形式和二进制形式，但只在单进程下有效。其中二进制保存的格式为PETSc规定的格式。使用求解器的类，即为依次调用类中的成员函数来实现有限元分析的过程。求解器的一个完整的使用示例如下:

SOLVER S;声明一个求解器对象

S.Solver Initialize(argc，args);求解器初始化

if(S.rank=0)

S.Read Information(Char*path);由根进程读入数据

Set Element(int Type，int maxnode);设置单元节点数，Type为单元类型编号

Set Element(int Type，proc－ESElement Stiffer);设置单刚求解程序

S.Preallocation();创建矩阵、向量结构、矩阵预分配

S.Add KERE();计算单刚、组集总刚和节点力列阵

Set KSP Type(char*ktype);设置求解迭代方法

Set PC Type(char*ptype);设置预条件方法

S.Set Bound Condition();根据边界条件修正支配方程

Set Tolerances(double abstol，double rtol，double dtol，int maxits);设置收敛误差

Save Binary();保存系数矩阵和右端顶，仅单进程运行时有效

肾功能指标水平：取上述血清3 ml，用肌氨酸氧化酶检测血肌酐(Scr)法酶偶联速率法检测血尿素氮(BUN)；采用核素动态显像法测定肾小球过滤率(GFR)。

Save ASCII();文本保存系数矩阵和右端项，仅单进程运行时有效

S.SOLVE();求解方程组

S.Get Answer();保存结果

Set Time Info();保存时间耗用明细

S.Finalize();结束

3 求解器的输入输出

求解器的输入目前提供了2种形式:读入文件和传入网格信息结构体指针。用户在调用Solver Initialize()函数进行初始化以后，即可在进程0里接受该结构体指针或执行读入文件函数。如以读入文件的方式输入，需要读入节点信息、单元信息、材料信息、约束信息、节点力信息。

求解器的输出包括3个部分:(a)方程构成。求解器能够在一个进程运行时输出方程的系数矩阵和右端项，默认的状态是不输出，用户需要在求解方程之前将其设置成输出。输出文件的格式有文本文档和二进制格式。文本文档可供查阅，保存时间较长，存储较大;二进制格式存储小，速度快，可由PETSc程序读入。(b)求解信息。求解网格信息的简要说明，求解设置，以及可选的时间信息，供性能分析使用。(c)方程的解。包括方程的解和节点的位移矢量大小，在文件的开始则为位移的最大值及其节点号。

4 算例分析

笔者以计算某软土地基上水闸在自重荷载下所产生的位移为算例，对本文所开发的并行求解器进行测试。水闸结构由3部分组成:闸身、桩基、软土地基，材料参数见表1。本文仅作线弹性分析，故将地基和桩基密度置为零，闸身密度为2500kg/m3，以计算由水闸自身重力产生的位移。因水闸结构对称分布，因此取其一半结构进行建模［16］计算(图2)。从求解器求解规模、并行效率、对于病态方程组的迭代收敛情况来测试本文开发的并行求解器。求解环境为一台四核I5、主频为2.8 GHz的32位台式机，可使用内存为3 G。

表1 水闸材料参数Table 1 Material parameters of water lock

4.1 求解规模

按网格疏密将水闸结构离散为八节点六面体单元网格，表2为6组有限元模型的网格信息，图2为网格单元数为51112的水闸网格图。对每一组有限元模型都进行四进程并行求解，计算的时间如表2所示。

第6组模型在计算时由于每个进程的内存使用都超过600 M，系统需要使用虚拟内存，导致计算时间达到2 409.33 s。在同样的计算环境下，单进程能计算节点数不超过40万的网格，Ansys软件能求解节点数不超过20万的线弹性问题，可见本文求解器能在普通的多核计算机下处理更大规模的网格。

取节点数为58 622的网格模型比较本文求解器和Ansys软件计算得到的节点位移，两者的x，y，z方向最大位移及总位移间的误差分别为0，0，0.007%，0.001%，最大位移值的误差不超过1/10000(表3)。将节点位移导入软件Tecplot中，绘制位移云图。本文求解器得到的位移云图与Ansys结果分布规律基本一致。求解器在提高收敛精度下能得到更精确的解，但因迭代次数增加而消耗更多的时间。

图2 水闸结构及整体网格Fig.2 Structure of water locks and overall grid

表2 不同规模的网格信息Table 2 Grid information with different scales

表3 本文求解器和Ansys软件计算的部分节点的结果对比Table 3 Comparison of node calculation results obtained by solver developed in this study and ANSYS software

不同规模网格计算得出的最大节点位移值如图3所示，可看出随着单元总数增加，节点最大位移值也逐渐增大，此现象符合有限单元法中位移近似解是从下方逼近精确解的理论，因而水闸结构的最大位移应该在0.067 m附近。

4.2 并行加速效果

为验证并行加速效果，选取节点总数为178 007的网格数据，分别进行1～4个进程的并行计算。迭代方法为共轭梯度法(CG)［17］，预条件子为块雅可比方法(bjacobi)。图4为程序在单机多进程并行求解时的并行加速比、并行效率与进程数相关关系曲线。

程序的并行分析部分主要包括矩阵预分配参数计算、形成总刚度矩阵和迭代求解方程组。表4列出了程序主要并行部分的并行加速比和并行效率。由表4可知，求解器用于计算求解矩阵分配参数的时间较少，对于一个节点数为17万的网格，计算时间不到1 s。计算单元刚度矩阵、组集总刚度矩阵时，由于并行程度较高，其加速效果比较理想。另外，求解的迭代次数随着进程数目的增加而变大。

图3 不同网格规模下最大节点位移Fig.3 Maximum displacement of node for different scales of grid

求解器由于其总刚度矩阵的计算采取了精确预分配空间的方法，只存非零元素，故程序运行时占用内存较小。多个进程运行时，每个进程占用内存会明显减小。

4.3 对于病态方程组的迭代收敛

在有限元中，若结构材料参数相差过大会导致支配方程呈病态，即出现迭代次数过多或者迭代不收敛的问题。PETSc提供的迭代法的预条件子，能有效地解决或缓解这一问题。不考虑实际情况，将闸身的弹性模量作表5所示改变，即改变闸身弹性模量使之分别相差10倍、100倍至100万倍。将这7组模型分别进行4个进程的并行迭代求解，得出的位移与Ansys软件计算出的位移进行对比(表6)。由表6可得，采用预条件子技术后，求解程序的迭代方法对因材料参数相差过大(最大相差10万倍)而产生的病态方程组具有良好的适应能力。

图4 1～4个进程计算时的并行加速比与并行效率Fig.4 Parallel speedup ratio and efficiency for one to four processes

表4 并行分析各部分的并行加速比和并行效率Table 4 Parallel speedup ratio and efficiency for various parts

表5 闸身弹性模量变化Table 5 Lock body elastic modulus change Pa

表6 修改材料参数后求解器的计算结果Table 6 Calculated results of solver after material parameters were modified

5 结 语

基于PETSc的数据结构(矩阵Mat和向量Vec)和线性求解器(KSP)组件开发了有限元求解器，并以水闸在自重荷载下所产生的位移为算例，对该求解器进行了测试。测试结果表明，该求解器能在普通的内存为3G的四核计算机上，计算高达80万节点的线弹性有限元问题，并行加速效果也较为明显(4个进程计算时加速比能达到3.24)。另外由于使用了PETSc提供的较为成熟的预条件子技术，求解器对病态方程组也有良好的适应能力，为求解由于材料参数相差10万倍而产生的病态方程组提供了新方法。由于现阶段各方面的限制，本文开发的求解器仅仅使用了PETSc工具箱的很小一部分功能，程序也有很多地方需要改进。基于PETSc的有限元高性能求解方法的研究和开发还可以进行以下几方面的深入研究:

a.对程序进一步完善扩展，能够对线性有限元问题作进一步的求解，如求出应力、应变。

b.将PETSc中的非线性求解器SNES和时间步进器TS对应到非线性有限元、时间依赖的动力学问题的求解上。

c.本文所开发的求解器是一个没有界面的程序，使用起来不方便，若能将其与图像处理软件结合起来或自己形成界面将会使该求解器更接近实际应用。

［1］范颖，王磊，章青.多尺度有限元法及其应用研究进展［J］.水利水电科技进展，2012，32(3):90-94.(FAN Ying，WANG Lei，ZHANG Qing.Research progress and application of multiscale finite element method［J］.Advances in Science and Technology of Water Resources，2012，32(3):90-94.(in Chinese))

［2］张汝清.概说并行计算结构力学［J］.计算结构力学及其应用，1995，12(4):447-484.(ZHANG Ruqing.Introduction to parallel computational structural mechanics［J］.Computational Structural Mechanics and Applications，1995，12(4):447-484.(in Chinese))

［3］张健飞.采用PETSc的有限元并行计算实现与优化［J］.计算机工程与应用，2010，46(10):57-59.(ZHANG Jianfei.Parallel implementation and tuning of finite element computing with PETSc［J］.Computer Engineering and Applications，2010，46(10):57-59.(in Chinese))

［4］姜弘道.水利高性能计算的进展［J］.水利水电科技进展，2006，26(2):70-76.(JIANG Hongdao.Progress in highperformance computing in water science［J］.Advances in Sciences and Technology of Water Resources，2006，26(2):70-76.(in Chinese))

［5］NOOR A K，VOIGT S J.Hypermatrix scheme for finite element system on CDC STAR-100 computer［J］.Computers and Structures，1975，5(5):287-296.

［6］ DIMITRIS Z.Parallel processing techniques for FE analysis:stiffness，loads and stress evaluation［J］.Computers and Structures，1988，28(2):247-260.

［7］ BITZARAKIS S，PAPADRAKAKIS M.Parallel solution techniques in computational structural mechanics［J］.Computer Methods in Applied Mechanics and Engineering，1997，148:75-104.

［8］CHIEN L S，SUNCT.Parallel processing techniques for finite element analysis of nonliear large truss structures［J］.Computers and Structures，1989，31(6):1023-1029.

［9］李丽君，金先龙，李渊印，等.一种新型并行化有限元结构模态分析集成系统［J］.计算力学学报，2004，21(5):546-550.(LI Lijun，JIN Xianlong，LI Yuanyin，et al.A new parallel integration system of finite element model analysis［J］.Chinese Journal of Computational Mechanics，2004，21(5):546-550.(in Chinese))

［10］张林波，迟学斌，莫则尧，等.并行计算导论［M］.北京:清华大学出版社，2006:1-487.

［11］Michael JQ.MPI与OpenMP并行程序设计［M］.北京:清华大学出版社，2004:34-40.

［12］都志辉.高性能计算并行编程技术MPI并行程序设计［M］.北京:清华大学出版社，2001:1-336.

［13］KRYLOV A N.On the numerical solution of the equation by which the frequency of small oscillations is determined in technical problems［J］.Izv Akad NaukAzerb Ser Fiz-Tekh Mat Nauk，1931，4:491-539.

［14］李晓梅，吴建平.Krylov子空间方法及其并行计算［J］.计算机科学，2005，32(1):19-20.(LI Xiaomei，WU Jianping.Krylovsubspace methods and parallel computation［J］.Computer Science，2005，32(1):19-20.(in Chinese))

［15］付成华，周洪波.面向对象的铁路路基加固设计系统集成与实现［J］.水利水电科技进展，2012，32(4):70-73.(FU Chenghua，ZHOU Hongbo.Development and implementation of object-oriented railway subgrade strengthening design system［J］.Advances in Science and Technology of Water Resources，2012，32(4):70-73.(in Chinese))

［16］林江，母德伟，余葵，等.万州长江公路大桥防撞装置结构稳定性及通航净空尺度［J］.水利水电科技进展，2013，33(2):59-62.(LIN Jiang，MU Dewei，YU Kui，et al.Structural stability and navigable clearance of anti-collision device of Wanzhou Yangtze River Highway Bridge［J］.Advances in Science and Technology of Water Resources，2013，33(2):59-62.(in Chinese))

［17］HESTENSEM R，STIEFEL E.Methods of conjugate gradients for solving linear systerms［J］.Tournal of Reaserch of the Natural Bureau of Standards，1952，49(6):409-436.