基于残余平滑-预处理共轭梯度算法的有限元并行计算

基于残余平滑-预处理共轭梯度算法的有限元并行计算
付朝江;陈洪均
【摘要】针对弹塑性问题的有限元分析非常耗时,基于消息传递接口(MPI)集群环境,提出了残余平滑的子结构预处理共轭梯度并行算法.采取区域分解,将子结构通过界面条件处理为独立的有限元模型.整体分析时,每个处理器仅存储与其相关的子结构信息并生成局部刚度矩阵.采用对角存储方式和最小残余平滑法,设计出了结合残余平滑(MR)的并行子结构预处理共轭梯度(PCG)算法.并行算法中对负载平衡进行了探讨,对处理器间的通信进行了优化.利用于步法对弹塑性应力应变进行积分,根据预定的容许值自动调整每个子步的大小来控制积分过程的误差.在工作站集群上实现了数值算例,分析了算法的性能,计算性能与传统的PCG算法进行了比较.算例显示:所提算法具有良好的加速比和效率,优于传统的PCG算法,对弹塑性问题的有限元分析,是一种有效的并行求解算法.
【期刊名称】《计算机应用》
【年(卷),期】2015(035)012
【总页数】5页(P3387-3391)
【关键词】预处理共轭梯度法;消息传递接口;并行计算;区域分解;有限元
【作者】付朝江;陈洪均
【作者单位】福建工程学院土木工程学院,福州350108;福建省土木工程新技术与信息化重点实验室,福州350108;福建工程学院土木工程学院,福州350108
【正文语种】中文
【中图分类】TP301.6;TP31
网络并行计算具有良好的性价比和可扩展性，以及灵活的体系结构等优点，已成为高性能计算的主要趋势[1]。

在弹塑性分析中，由于材料的非线性，使得问题求解变得非常复杂和耗时。

有限元法在非线性分析中的应用取得的进展和并行计算的出现为高效快速求解这类耗时的弹塑性问题开辟了一条新途径。

有限元结构分析的高性能计算成为一个重要的研究方向[2-3]。

对弹塑性问题的有限元分析、探讨和开发高效的并行求解算法具有重要的意义。

非线性弹塑性问题的分析是以增量方式进行，因为求解过程不仅取决于结构的当前位移而且与加载历史有关。

在大型三维弹塑性有限元分析中，应变和应力关系的积分和方程组的求解成为两个关键的阶段。

文献[4-5]在并行计算机系统上实现了弹塑性分析，但积分方式是基于传统的有限元方法，没有采取任何措施来控制积分过程中的误差。

在非线性弹塑性分析中，计算应力的误差是很关键的，因为太大的误差可能会导致长时间的迭代才能收敛甚至导致解发散。

文献[6-7]介绍了准确积分弹塑性本构关系的算法，是利用修正的欧拉方法和Runge-Kutta-England方法[8]，该方法在积分过程中利用不同阶的截断误差来确定子步的大小。

为此，本文采取修正的子步法对弹塑性应力进行积分，该方法适用于任何类型的本构模型，它能通过自动调整每个子步的大小控制积分过程的误差。

位移计算成为并行实现的瓶颈，因为在这一阶段需要进行整体矩阵和向量的组集。

从并行处理的角度希望可以避免整体方程的组集，这样可以在单元级进行运算。

在有限元分析中迭代法的优点是它不会改变刚度矩阵结构并保持其稀疏性。

因此，通过使用适当的存储方式，可以大大减小计算成本和涉及零填充的存储空间[9]。

文献[10]采取多重网格法实现弹塑性问题的有限元计算。

文献[11]采取多级区域分解方法对弹塑性问题的有限元分析进行计算，利用共轭梯度(Conjugate Gradient，CG)法求解线性方程组，Newton-Raphson法求解非线性问题。

为此，本文将预
处理共轭梯度(Preconditioned Conjugate Gradient, PCG)法和最小残余平滑(Minimal Residual, MR)法相结合, 采取修正的子步法对弹塑性问题的并行求解算法进行研究。

本文算法不需要进行整体系统矩阵的组集，每个子结构的位移直接进行计算。

为减少计算成本和存储量，使用对角存储方式来存储每个子结构的刚度矩阵。

编程中使用单元编号和处理器间通信优化。

对弹塑性问题，假定弹塑性模型是基于弹性和各向同性硬化的塑性模型，其应变增量可表示为：
Δε=Δεe+Δεp
式中：Δε为总应变增量；Δεe和Δεp分别为弹性和塑性应变增量。

当应力满足屈服准则时屈服才会出现：
F(σ,k)=0
式中k为硬化参数，其取决于塑性加载历程。

当塑性屈服出现时，应力满足式(2)的屈服面。

进行微分可得到：
式(2)和式(3)称为连续条件，当塑性屈服出现时必须满足。

利用式(1)～(3)进行变换，弹塑性增量应力-应变关系可表示为：
Δσ=De:Δεe-ΔλDe:aT
式中
将弹塑性应力-应变积分过程划分为多个子步，计算每个子步的应力-应变反应。

当然，随着子步数的增加，其精度增加。

但每个子步的计算非常耗时，因为向量和弹塑性矩阵必须在每个子步计算。

传统上，子步数由经验确定并且每个子步假定为相同大小。

这种积分一般会导致应力变化偏离屈服面，采用一定比例的应力缩放使其恢复到屈服面。

本文采取的子步法是在进行弹塑性应力-应变关系积分时，通过自动调整每个子步的大小来控制误差。

这样，可以避免出现为满足适当的精度而子步数盲目地增加，
积分过程可得以控制。

为了在积分时实现更好的精度，构成弹塑性反应的应变增量通常被划分成很多子步数。

由于确定子步数基于试错法，在整个积分过程中没有误差控制。

通常利用应力校正将应力回归到屈服面。

本文采用Runge-Kutta-Dormand-Prince (RKDP)方
法[12]。

应力-应变积分需要初值问题的解，如：
式中：σ(T=0)为满足屈服准则的相应应力；σ(T=1)表示在荷载增量结束时的应力
状态。

为估算在求解中的局部截断误差，采取嵌入技术。

局部误差定义为采用不同阶的两种方法获得的解之差。

因此，获得的局部截断误差估算为：
在Tk到Tk+1=Tk+ΔTk子步的局部误差估算，定义该子步的相对误差为：
Rk+1=‖Ek+1‖/‖‖
将Rk+1与给定的容许值TOL比较，如果Rk+1≤TOL，该步可接受;否则不采用。

并利用Rk+1值进行渐近优化子步大小的估算在不采用时，用ΔTk+1替代ΔTk;接受时，使用ΔTk+1继续积分。

通过控制每个子步的局部相对误差，可实现求解过程的整体误差控制。

并行有限元分析最常用的方法是根据可用处理器的数目将整个结构划分为多个子结构。

涉及有限元分析的绝大部分计算是在单元级进行。

因此，这些计算可以并行进行，不需要任何对子结构级的同步[13]。

然而，在整体系统中涉及节点变量的组集和求解过程是串行进行。

这样，方程组的求解就可以在子结构级进行。

将单元接单元的PCG算法[14]扩展为子结构接子结构的算法，即子结构接子结构的PCG算法，并采用对角预处理。

这样，求解过程中不需要形成整体方程系统。

为进一步加快PCG算法的收敛，采用并行最小残余(MR)平滑法[15-16]。

这样得到结合MR平
滑的并行子结构预处理共轭梯度(Parallel Preconditioned Conjugate Gradient, PPCG)算法，描述如下：
其中： g代表组集向量；C为条件矩阵；上标s(s=1,2,…,n)表示子结构编号；表示局部通信求和；∑表示全局通信求和。

4.1 并行网格生成
并行有限元方法的第一步是根据可用的处理器数将有限元区域划分为子结构，将子结构单元的计算工作分配给各个处理器。

为了避免这一阶段的通信，每个处理器执行相同的指令，不同的数据。

对每个子结构建立单独的输入和输出文件。

这些文件读取和写入本地副本是并行进行。

读完相应的输入文件，每个处理器自动生成子结构，这个过程无需任何通信。

此方法的优点是，每个处理器只需存储各自使用的数据，这可以减少存储量。

并行有限元程序框图如图1所示。

4.2 单元编号
单元编号采用一定的编号顺序，即单元编号从一端开始，在整个子结构中按顺序进行，这样，子结构的终端是另一相邻子结构的起始端。

这种编号的优点是：编程相对容易；容易实现两个相邻子结构间的通信。

4.3 矩阵存储方式
通常存储一个n×n的稠密矩阵使用一个n×n数组。

对稀疏矩阵，因为在矩阵中大多数的元素为零而不需要显式存储，只需存储非零元素值和其在矩阵中的位置。

这样不仅节省了存储空间，还减小了计算。

因为矩阵中非零元素的位置是已知的，可以避免不必要的与零元素的运算。

有多种存储方案可用于存储和处理稀疏矩阵。

本文采用对角存储方案。

以图2为例：考虑一个n×n的稀疏矩阵，如图2(a)，有非零元素的d个对角线；这些非零的对角线都存储在n×d数组DARRAY，如图
2(b)；一个d×1的OFFSET数组用于存储相对主对角线的每个对角线的偏移值，如图2(c)。

由于主对角线以外的所有对角线都有少于n个元素，在DARRAY数组中会有未使用的位置。

显然，稀疏度很高的矩阵，d远远小于n，这样可节省零填充的很多时间。

4.4 负载平衡
处理器间的负载平衡是实现高性能的关键。

好的负载平衡取决于如何将整个结构分配给各个处理器而且不会涉及太多的通信。

以图3所示的悬臂梁为例，有三种划
分方式，即垂直划分、水平划分和块划分，分别如图4～6。

若悬臂梁左端为固定端，力作用于悬臂梁右端。

当增加作用力时，塑性域将首先在最左端开始，然后沿水平向右发展。

可以看出，第一种划分方式可以减少界面节点，但在初始几个加载增量步，当结构的左端出现塑性而右端依然是弹性时，可能会导致负载不平衡。

第二个划分方式能更好地实现负载平衡，因为它将塑性区均匀地分配给每个处理器。

但在分析中涉及较多的界面节点，这样将导致更多的通信。

第三种划分方式，为图6的分块划分，不仅会导致负载不平衡，而且使得编程更复杂。

因而，本文只采用第一和第二种划分方式进行算例计算。

4.5 处理器间通信
对第一种和第二种两种划分方式，每个子结构仅有左端(底部)和右端(顶部)与其他
处理器相邻，所以通信可通过以下步骤实现：
1)每个处理器发送右端(顶部)界面的数据给它右边(顶部)相邻的处理器，然后接收
其左端(底部)界面相应的数据。

2)每个处理器按相反方向重复相同的过程。

为有效地通过消息传递接口 (Message Passing Interface,MPI) 实现此通信，通
常采取图7所示的发送和接收方式。

在开始时，只有一个处理器(“顶部”处理器
P7)不发送给任何处理器，它可以接收它下面的处理器信息。

本文采取按发送和接收顺序方式实现这种通信，以便如果一个进程正在发送给另一个，目的进程将接收相匹配的数据，发送之前完成自身数据的接收。

图8显示了使用这种方法的通信
模式。

可以看到，偶数处理器先发送，奇数处理器先接收。

之后，奇数处理器发送，偶数处理器接收。

这将减少造成的发送-接收匹配的闲置时间。

图7～8中阴影区
域表示处理时间的闲置。

数值计算的并行计算环境是在DELL工作站机群上进行。

该机群由4台双CPU的DELL工作站通过100 Mb/s以太网连接而成，共有8个CPU(2.4 GHz Xeon chips,512 KB cache),每个节点内存1 GB。

每台工作站都有真实的互联网协议(Internet Protocol, IP)地址，使用消息传递接口(MPI)作为并行开发环境[16]。

将算法应用于悬臂梁的非线性弹塑性有限元分析。

矩形截面悬臂梁如图3所示，长1 600 mm、宽200 mm、高400 mm。

材料常数为E=2.0×104 N/mm2、各向同性硬化，硬化后弹模E1=0.36×104 N/mm2、ν=0.3、σy=24.0 N/mm2,右端自由端作用力P=100 kN。

采用冯·米塞斯屈服准则，有效应力和有效塑性应变关系采用Swift硬化规则。

对该模型使用两种划分方式，即垂直划分方式和水平划分方式。

将整个结构划分为与处理器数相等的子结构，每个子结构有相等的单元数。

每个子结构映射到各相应的处理器。

各处理器得到子结构的单元数据后，相应子结构中的计算在各处理器中进行。

采用20节点的三维块单元。

网格划分如图4、5所示。

为评价并行性能，分别用1、2、4、6、8个处理器对所求问题进行求解，计算结果如表1和表2所示。

从表1和表2所示的结果可以看出，虽然垂直划分方式由于塑性区出现可能导致处理器的工作负载不平衡，但仍然实现良好的加速比和效率。

当随着问题规模增大(单元数增加)，加速比和效率增加，算法的性能提高。

这良好的性能归因于采取的子步法和结合残余平滑的预处理共轭梯度法及使用优化处理器间的通信和对角存储方案。

相反，水平划分方式使每个处理器参与计算塑性区域，但结果变差。

从表2可以看出，性能不能随着问题规模的增加而提高。

这是因为在大规模问题的总体分析中涉及更多的通信。

对大规模问题，虽然相当多的时间花在计算上，但很多时间也花费在处理器间的通信。

求解方程是有限元分析最耗时的部分，需要相当多的通
信。

随共享节点数量变大，在这一阶段要耗费大量的时间，这就是增加问题规模而加速比和效率下降的原因。

这也说明：负载不平衡引起的闲置时间会降低性能，如果并行分析涉及太多的通信，良好的负载平衡造成的这种划分优势会因涉及太多的通信而抵消，会降低其性能。

本文的残余平滑并行PCG算法是采取子步法和结合残余平滑的预处理共轭梯度法，并使用对角存储方式和优化处理器间通信技术。

为验证本文算法的良好并行性能，采用传统的PCG算法(即没有采取残余平滑法和对角存储方式和优化处理器间通信技术的PCG算法)进行比较，结果如表3和表4所示。

可以看出，对相同规模的
问题，本文算法比传统PCG算法的CPU耗时少，相应的加速比和效率高，表明
本文算法较传统的PCG算法的性能更优越，显示出本文算法具有良好的性能。

在MPI集群环境下，针对弹塑性有限元分析提出了子步法和结合残余平滑的并行
预处理共轭梯度算法。

在积分过程中采用子步法，该方法依据弹塑性本构关系的特性改变每个子步的大小，从而控制积分过程中的误差。

研究了采用MR平滑法的
子结构接子结构的并行PCG算法，该算法不需要组集整体系统方程组。

每个处理器分配到子结构后，存储子结构的相关信息，编程中使用单元编号和处理器间通信优化。

将所提算法应用于弹塑性有限元分析，实验结果表明，该算法可得到良好的加速和效率，是一种有效的切实可行的并行求解算法。

【相关文献】
[1] MO Z, YUAN G. Message passing parallel programming environment MPI [M]. Beijing: Science Press, 2001:1-11. (莫则尧，袁国兴. 消息传递并行编程环境MPI[M]. 北京：科学出版社，2001:1-11).
[2] FERIAI A, FRANCHI A, GENNA F. An incremental elastic-plastic finite element solver in a workstation cluster environment Part 1. Formulations and parallel processing [J]. Computer Methods in Applied Mechanics and Engineering, 1996,130(3/4): 289-298.
[3] GUO Y, JIN X, DING J. Parallel numerical simulation with domain decomposition for
seismic response analysis of shield tunnel [J]. Advances in Engineering Software, 2006,
37(7): 450-456.
[4] YU T, JIANG H. Elasto-plastic substructure parallel finite element method with multi-front parallel processing [J].Chinese Journal of Computational Mechanics, 1999, 16(4):
493-496. (余天堂,姜弘道. 多波前并行处理的弹塑性子结构并行有限元[J].计算力学学报,1999,16(4): 493-496.)
[5] MIYAMURA T, NOGUCHI H, SHIOYA R, et al. Elastic-plastic analysis of nuclear structures with million of DOFs using the hierarchical domain decomposition method [J]. Nuclear Engineering and Design, 2012, 212(1/2/3): 335-355.
[6] PEZESHK S, CAMP C V. An explicit time integration technique for dynamic analysis [J]. International Journal for Numerical Methods in Engineering, 1995, 38(13):2265-2281. [7] YAMANAKA N, OGITA T, RUMP S M. A parallel algorithm for accurate dot product [J]. Parallel Computing, 2008, 34(6/7/8): 392-410.
[8] SLOAN S W. Sub-stepping schemes for the numerical integration of elastoplastic stress-strain relations [J]. International Journal for Numerical Methods in Engineering, 1987, 24(5): 893-911.
[9] FARHAT C, CRIVELLI L. A general approach to nonlinear FE computations on share-memory multiprocessors [J]. Computer Methods in Applied Mechanics and Engineering, 1989, 72(2): 153-171.
[10] KACOU S, PARSONS I D. A parallel multigrid method for history-dependent elastoplasticity computations [J]. Computer Methods in Applied Mechanics and Engineering, 1993, 108(1/2): 1-21.
[11] WISSMANN J W, HAUCK C. Efficient elastic-plastic finite element analysis with higher order stress-point algorithms [J]. Computers & Structures, 1988, 17(1): 89-95.
[12] COELHO P G, CARDOSO J B, FERNANDES P R, et al. Parallel computing techniques applied to the simultaneous design of structure and material [J]. Advances in Engineering Software, 2011, 42(5): 219-227.
[13] PAVARINO L F. BDDC and FETI-DP preconditioners for spectral element discretizations [J]. Computer Methods in Applied Mechanics and Engineering, 2007,
196(8): 1380-1388.
[14] DING Z, KALYANASUNDARAM S, GROSZY L, et al. Development of a new integration algorithm for parallel implementation of the finite element elasto-plastic analysis [J]. Australian & New Zealand Industrial and Applied Mathematics Journal, 2006, 42(E): C561-C585.
[15] KWAK J Y, CHUN T Y, SHIN S J, et al. Domain decomposition approach to flexible multibody dynamics simulation [J]. Computational Mechanics, 2014, 53(1):147-158. [16] FU C. The research on parallel computation of finite element structural analysis based
on MPI cluster [D]. Shanghai: Shanghai University, 2006: 18-35. (付朝江.集群MPI环境下有限元结构分析并行计算研究[D].上海：上海大学,2006:18-35).。