并行有限元技术研究进展

并行有限元技术研究进展

张巍

上海市隧道工程轨道交通设计研究院上海200235

摘要：对于工程中越来越多的大型复杂结构，并行化处理已成为实现高性能计算的重要方式。本文首先简述了高性能计算机的发展历程，然后总结了针对系统方程组的并行化求解的技术机理与特点，并对结构动力分析和非线性分析的并行算法进行了讨论，指出了该领域存在的问题和研究发展方向，对相关研究提供参考。

关键词：有限元；高性能计算；并行计算；集群

1引言

高性能计算(HPC,HighPerformanceComputing)能力最原始的获取方法是以高成本为代价的，即从单方面提高元件处理速度和时钟频率等方法入手。对于现在越来越多的超大化、复杂化的工程结构，进行整体动力分析时要保证数值解的精度，不可避免会遇到巨大计算量的问题，时常还需考虑非线性本构关系、流固耦合、大变形与复杂边界条件等复杂情况，普通计算机的处理能力成为分析问题的瓶颈。并行技术作为更便捷、更有效的提升计算速度的方法，已成为高性能计算的重要实现方式[1,2]。

本文首先介绍了高性能计算机的发展历程，然后分析了利用并行数值求解方法针对系统方程组求解过程进行并行化和分割系统方程组进行并行求解的方法，对动力分析直接积分策略的隐式、显式方法和非线性分析方法进行了讨论。最后对未来并行有限元技术的发展方向进行了展望。

2高性能计算机的发展历程

20世纪70年代出现的向量计算机可以看作是第一代HPC。其中较著名的有CDC系列、CRAY系列、NEC的SX系列和中国的银河一号及中科院计算所的757计算机。它们都属于单指令多数据流(SIMD,SingleInstructionMultipleData)类型，基本上都获得了很好的向量计算效果。

20世纪80年代出现的MIMD(MultipleInstructionMultipleData)并行计算机的性能有了较大提高。早期的MIMD并行机属于共享存储式并行架构，无论是对称多处理机(SMP,SymmetricMultipleProcessor)结构还是分布共享存储(DSM,DistributedSharedMemory)MIMD并行计算机，都可以保持较好的负载平衡。但共享式存储MIMD并行机的最大缺陷在于内存访问瓶颈造成的扩展性问题，之后出现的分布式存储MIMD并行机很好地克服了这一点。分布式并行计算机又分为两类：一类是用专门网络紧密耦合多个计算节点而构成的大规模并行处理机(MPP,MassivelyParallelProcessor)，例如TMC的CM-5、IntelParagon等；另一类是由多台高性能计算节点甚至是普通微机经网络连接组成的集群并行机(COW,ClusterofWorkstation)，例如IBMSP2、国内的曙光4000、5000、联想深腾6800、天河一号等。这些典型的分布式存储计算机系统已经逐渐成为目前高性能计算机领域的重要应用模式。尤其是集群系统，由于其投资风险小、性价比高、结构灵活、可扩展性强、通用性能好、异构能力强等较多优点，得到迅速普及应用[3,4]。在目前的世界高性能计算机TOP500中，绝大多数都属于集群和MPP的架构。

3并行有限元计算的研究现状

有限元化整为零、积零为整的基本思想，与并行处理技术的基本原则“分而治之”(pideandconquer)是高度一致的，从NASA的A.K.Noor于1975年发表第一篇关于有限元并行化计算的文章以来，有限元高性能计算的研究已经历了30余年的发展[5-8]，获得了许多重要的研究成果。

3.1针对系统方程组的并行化方法

并行有限元计算最直接简便的实现方法是针对系统方程组进行并行计算。即单元和系统方程的形成是串行过程，而系统方程组的求解是并行过程，其原理如图1所示。针对线性方程组的并行求解策略有直接法和迭代法两种。直接解法包括Gauss消去法、LU分解法、Cholesky分解法和Givens分解法等。直接解法的主要优点是算法稳定性好，精度高，同时便于对现有串行程度进行修改，但主要缺点是在求解方程过程中，由于多次进行通讯，会花费不必要的进程等待信息和传递信息的时间，从而降低并行效率。迭代解法包括预处理共轭梯度法(PCG)、最小残余法、逐次超松弛(SOR)法和Gauss-Seidel法等。与直接法相比，迭代解法占用内存较小，在理论上可以达到更高的收敛精度，但迭代法也存在有时收敛慢甚至无法收敛的问题。

目前已有一些成熟的软件包可供使用，使得此种并行策略更为简便。表1列出了较为典型的一些软件包[9]。夏江[10]基于PETSc函数库开发了PSAP程序，在集群系统上实现了与通用有限元程序ANSYS的无缝连接，对场地地震反应过程中方程组的求解过程进行了并行化计算。李根国等[11]将NASTRAN在“神威I”超级计算机系统进行了移植和并行计算功能的二次开发，他们的思路是将有限元的求解过程(线性方程组的求解和矩阵的特征值问题)采用PETSc软件包进行并行化，前后的处理工作由NASTRAN软件完成，模块间通过文件进行数据交换。

另一种做法是利用区域分解方法对任务进行“分割”，各处理机存储与每块范围相应的有关数据(单元、节点、荷载等)，这样单元分析和系统方程组装过程就能以少量的重复计算为代价并行进行[12]，原结构划分成若干子区域，划分面为“共享单元”。每个处理机保存各自子区域的相关信息(包括“共享单元”的信息)，子区域内的单元分析和总刚的集成是完全并行的，而只有“共享单元”涉及重复计算，而系统方程组的求解一般采用迭代并行求解器。计算过程如图2所示。

3.4结构动力分析和非线性分析的并行算法

动力和非线性问题的计算成本更高，所以对它们进行并行分析是十分重要且有意义的。动力分析可分为结构固有特性分析和结构动力响应分析。模态分析主要涉及矩阵特征值问题，计算方法主要有Lanczos方法和子结构模态综合法等，目前已有一些成熟的算法和数值软件包供使用。

(1)

对于结构动力学方程(式1)的直接积分求解策略，可分显式积分方法(中心差分法、线性加速度法等)、隐式积分方法(Newmark法、Wilson-?法等)和混合积分方法几种。隐式算法的扩展性较差，一般在8-16cpu以上计算效率就很难获得进一步的提升。相对于隐式算法，显式算法仅仅使用前一时间步的数据来更新所求的解而不需要求解整个方程系统，比隐式方法更简单和占用更小的内存、通信带宽，故扩展性较好，一般在64-128cpu以内都能获得较好的并行性能。然而显式算法是有条件稳定的，一般需要很小的时间步。当采用集中质量矩阵时，显式积分法可针对对角阵简化方程组求解过程，从而提高并行效率。程建钢等[13]在分布式环境下进行了结构动力分析两种积分方法的并行化讨论，对空间板弯曲问题进行了瞬态并行计算。

非线性有限元问题的常用分析方法有增量法、迭代法和混合法等，其过程均相当于等效为一系列静力问题的求解。对于弹塑性、粘塑性、蠕变损伤、接触、岩土渗透固结等具有复杂本构关系的非线性问题，经典变分原理不能适应，钟万勰和张柔雷[14]、曾攀[15]等相继在这方面作过系统的研究工作，提出了这些问题相应的参变量变分原理。

4结语与展望

并行有限元技术的实现和高性能计算机与相应并行算法的发展密不可分。基于硬件技术和网络通信的快速发展，基于高速互联网的网格计算和云计算是未来高性能计算技术的重要应用方向，它们可以解决仅靠本地资源无法解决的问题，但数据的安全性和计算的容错能力仍然面临巨大挑战[16,17]。此外正在探索和研究中的完全属于另一类并行的有：以分子作为处理单元来开发并行性的生物元计算机、使用叠加原理开发并行性的量子计算机，包括利用超导元件制造的CPU等技术，它们在未来也会对高性能计算产生巨大的影响。

并行算法方面，如何充分考虑存储、计算和通信之间的平衡，实现最优的计算能力是目前需要注意的一个重要问题，目前最简便的实现方式是针对系统方程组的并行求解。自动并行化也是研究热点之一，目前已有一些相关研究[18,19]。但其计算效率受HPC架构、任务划分、通信速度、计算规模等影响较大，如何在有限元问题自动并行化求解时实现对高性能计算机的最有效利用，仍需进一步研究。另外，目前衡量高性能计算能力的主要指标仅是峰值速度和并行效率，并行有限元求解系统性能的标准尚不明确，缺乏权衡考虑通信性能、问题规模、分析效率、计算成本等整体性价比因素的一个评价体系。

参考文献

[1]KUMARV,etal.IntroductiontoParallelComputing:AlgorithmDesignandAnalysis[M].BenjaminCummings/AddisonWesley,RedwoodCity,CA,1994.

[2]李晓梅,程建刚,李明瑞.并行计算环境与有限元分析并行算法[J].中国科学基金,1995,9(1):15-21.

[3]RajkumarBuyya.高性能集群计算[M].北京:电子工业出版社,2001.

[4]MARKB.ClusterComputingWhitePaper[M].UniversityofPortsmouth,UK.2000.

[5]阮红河,袁勇,柳献.分布存储环境并行有限元研究进展[J].同济大学学报.2005,33(1):21-27.

[6]付朝江,张武.结构有限元分析并行处理的研究进展[J].力学进展.2006,36(3):354-362.

[7]金先龙,李渊印.结构动力学并行计算方法及应用[M].北京:国防工业出版社.2008.

[8]SBitzarakis,MPapadrakakis,AKotsoulos.Parallelsolutiontechniquesincomputationalstructuralmechanics[J].Comput.MethodsAppl.Mech.Engrg.1997(148):75-104.

[9]谷同祥.并行数值分析软件综述[R].北京:中国科学院软件研究所.2003.

[10]夏江.基于PC集群系统的场地地震反应并行计算研究[D].上海:同济大学.2007.

[11]李根国,蔡贵新,李渊印等.商业性有限元分析软件在神威超级计算机上的并行开发[J].结构工程师,2003,66：318-323.

[12]PapadrakakisM,BitzarakisS.DomaindecompositionPCGmethodsforserialandparallelprocessing[J].AdvancesinEngineeringSoftware,1996(25):291-307.

[13]程建钢,姚振汉,李明瑞等.结构动力分析显式积分并行算法与实现[J].清华大学学报(自然科学版),1996,36(10):80-85.

[14]ZhongWX,ZhangRL.Aparametricvariationalprinciplesforidealelastoplasticproblems.ActaMechanicaSinica,1988,4(2):125-131.

[15]曾攀,钱令希.岩土中弹塑性渗透固结问题的参变量变分原理.力学学报,1991,23(4):484-490.

[16]SEGALB.GridComputing:TheEuropeanDataProject[C]IEEENuclearScienceSymposiumandMedicalImagingConference.Lyon:IEEESocietyPress,2000:15-20.

[17]网络并行：刘勇,朱国强.基于TCP/IP的网络并行计算的实现.计算机工程.2000,26(8):187-189.

[18]WANShui,HuHong,CHENJianpin.FiniteElementProgramGeneratorandItsApplicationinEngineering.船舶力学,2004,8(3):107-115.

[19]王立涛,陈厚群,马怀发.人工黏弹性边界的接触非线性问题在FEPG中的实现.水力发电学报,2009,28(5):179-181.

作者简介

张巍(1982-)，男，山西太原人，博士，研究方向：工程结构系统抗震与地下结构设计。