基于GPU的杆系离散元并行算法在大型工程结构中的应用

2021-03-02叶继红

工程力学 2021年2期

叶继红，王佳

(1. 中国矿业大学江苏省土木工程环境灾变与结构可靠性重点实验室，徐州 221116；2. 中国矿业大学徐州市工程结构火灾安全重点实验室，徐州 221116)

20世纪70年代，Feng等[1]率先提出离散单元法(discrete element method，DEM)，随后在岩土工程领域获得了广泛应用。DEM属于显式求解算法，DEM允许单元间发生相对运动，不需要满足位移连续条件和变形协调条件，计算中不需要集成整体刚度矩阵，相比于传统有限单元法，求解计算不存在收敛性问题，非常适合求解强非线性问题。尽管离散单元法己被证明是求解强非线性问题的有效方法，并且广泛运用在各种研究领域，更是拥有传统有限单元法难以企及的独特优势，但是DEM计算效率普遍偏低一直是一个亟待解决的问题。

目前，在并行计算领域基于多处理器并行技术[2]与CUDA(compute unified device architecture)计算架构是解决大规模离散元数值计算的有效途径。对于多处理器并行技术而言，要实现多机协同工作，针对解决DEM这种全局问题，编程难度偏大；而基于CUDA的GPU(graphic processing unit)高性能计算，近几年应用领域则越来越广，例如在DEM相关领域，分子动力学、地质工程、原子间作用等方面均得到了广泛应用。

在国际上，He等[3]开发出基于GPU的离散单元法，用于大规模粉末压实模拟，大幅缩短了仿真计算时间。Chase等[4]将GPU并行技术运用到Yade开放式DEM软件中，运用GPU多线程技术加速矩阵因子分解，将多孔弹性渗流数值分析效率提高了170倍。 Spellings等[5]运用GPU加速DEM算法，用于模拟重力作用下各向异性颗粒系统，获得了比较显著的加速效果。Liu等[6]采用GPU加速Blaze-DEM算法，并将其运用到大规模颗粒材料破损分析中。在国内，付帅旗等[7]充分利用GPU众核多线程的计算优势，实现了大规模颗粒离散元接触模拟。在基于大规模连续介质离散元(continuum-based discrete element method，CDEM)计算方面，中科院力学研究所的研究表明[8]，在CDEM中使用GPU并行技术可以使计算效率提高100倍以上，GPU并行架构的优越性不言而喻。

在结构工程领域，杆系DEM是求解结构强非线性问题的有效方法，相比于岩土工程领域广泛应用的CDEM计算方法，杆系DEM在单元变形计算与单元内力求解上更为复杂，且数据计算复杂度更高。因此相比于CDEM并行算法，杆系DEM并行算法的设计难度更大。

为了设计高效的杆系DEM并行算法，本研究提出单元级并行、节点级并行的计算方法，并对杆系DEM的数据存储方式、GPU线程计算模式、节点物理量集成方式以及数据传输模式进行了详细设计。基于CPU-GPU异构平台建立了杆系DEM并行计算框架并编制了相应的几何非线性计算程序，并将其运用到大型工程结构数值计算中，获得了显著的加速效果。

1 杆系离散元几何非线性计算理论——运动方程的建立与求解

运用杆系离散元求解结构力学问题[9]，首要步骤是建立结构几何模型，主要包括确定关键点坐标和杆件连接，然后将模型离散为一系列连续的虚拟球体。如图1所示，球心为节点位置，两相邻球体球心所包含的区域为单元所在位置；以球心为中心，离散球体所包含的单元总质量等于节点的集中质量。

图 1 空间杆系结构离散模型Fig. 1 Discrete model of space framed structure

三维杆系离散元接触本构方程基于接触力学基本理论[10]，单元与单元发生接触后，将产生接触变形和接触力，根据式(2)可将接触力和接触力矩分解为沿接触平面方向的Fτ、Mτ和垂直于接触平面方向的Fn、Mn。

2 基于异构平台杆系离散元几何非线性并行算法设计

2.1 杆系离散元计算理论并行性分析

杆系离散元计算理论以单元和节点计算为基础，在单元和节点之上建立以时间步为基准的迭代计算，单元与节点计算几乎完全独立，并行计算特征十分鲜明。根据杆系DEM计算理论潜在的并行性，可将杆系离散元数据计算分为3大类：单元计算、节点计算、单元-节点计算，如图2所示。

图 2 杆系离散元并行性分析Fig. 2 Parallelism analysis of discrete elements

对于单元与节点计算，如单元内力计算、节点位移计算等，所有参与运算的单元或者节点相互独立，具备天然并行性，完全符合GPU多线程并行计算模式。对于单元-节点计算，如节点物理量集成，所有参与计算的单元和节点之间需要进行相互转换，此时的数据计算并不独立于单元或者节点，若直接采用GPU多线程计算，则会出现线程冲突，导致数据计算错误，此时需要采用辅助算法才能正确计算。

2.2 杆系DEM数据存储方式

杆系DEM几何非线性计算程序涉及大量多维数组，本研究对算法中所有多维数组采用降维存储，即一维数组存储多维数据。图3为杆系DEM数据存储形式，“一维数组存储多维数据”按照各维度顺序依次存储在GPU全局内存中，每一维度数据存储完成再进入下一维度数据存储，直至所有维度数据存储完成为止。

2.3 杆系DEM线程计算模式

杆系DEM采用一维数组存储多维数据，而在多维数据的每一维度上，数据量往往是随机的，并不是32的整数倍，自然也不满足GPU线程对数据的合并访问要求[12]，这将导致GPU线程对数据的访问效率大打折扣。为了满足GPU线程对数据的合并访问要求，本研究对一维数组的多维数据做对齐处理。如图4所示，采用CUDA中的cudaMallocPitch( )函数，在多维数据每一维度的最后添加一个或多个单位数据内存空间，使每一维度的数据量为32的整数倍，对齐后的一维数组可满足线程对内存的合并访问要求。

2.4 节点物理量集成方式

在杆系DEM计算理论中，每一时间步都需要将单元内力、单元转动惯性矩由单元集成到节点上，此过程涉及单元-节点转换计算，运算形式将不再独立于单元或者节点，采用GPU多线计算模式将会导致计算结果错误。

图 3 杆系DEM数据存储Fig. 3 Framed system DEM data storage

图 4 杆系DEM数据对齐Fig. 4 Framed system DEM data alignment

为了使GPU线程能够正确计算节点物理量，本研究提出两种解决方案。方案一是按单元集成节点物理量，即以单元计算为中心，通过单元对应节点的方式，在同一节点上累加单元所对应的物理数据，使用原子加操作[13](对应CUDA中atomicAdd( )函数)完成运算，如图5、图6(a)所示。方案二是按节点集成节点物理量，即以节点计算为中心，将不同节点对应的单元完全分开，同一节点所对应的不同单元物理量直接累加，从而使各节点之间的计算相互独立，在GPU中采用多线程循环计算，如图5、图6(b)所示。

图 5 杆系DEM单元连接模型Fig. 5 The connection model of the DEM unit of the framed system

图 6 节点物理量集成方式Fig. 6 Node physical quantity integration method

对于上述2种计算方案，若按单元集成节点物理量，使用原子加操作将导致所有线程排列执行，因而需要耗费很长时间；若按节点集成节点物理量，需要事先建立节点-单元索引数组，当节点对应的单元数目差异较大时，GPU计算将会产生大量线程分支以及严重的负载不平衡现象，这将引起较大的计算性能损耗。

两种方案实际的GPU程序运行结果如表1所示，可以看出两者在计算时间上几乎相同。考虑到按节点集成节点物理量需要采用辅助索引数组，而创建索引数组不仅增加了内存需求，还需要额外的计算时间，所以本研究选用按单元集成节点物理量。

表 1 两种方案集成节点物理量所需的计算时间Table 1 The calculation time required to integrate the physical quantities of the two schemes

2.5 多线程数据传输优化

为了获得杆系DEM计算过程中各时间子步产生的结构响应，需要在执行数个计算步后从GPU全局内存中将节点位移计算结果存储至硬盘上。传统的异构平台数据传输模式如图7(a)所示，GPU计算完成后，必须等待数据在硬盘上全部存储完成后，才能进行下一步计算。由于数据传输过程始终为串行执行，根据Amdahl定律[14]，数据传输所消耗的时间将严重限制计算性能的提升。

本研究为了突破数据存储与数据计算同步执行带来的性能瓶颈，运用C++11中thread类函数接口与CUDA中cudaMemcpyAsync( )函数接口，实现了异构平台数据计算与数据存储异步执行。异步执行模式如图7(b)所示，CPU主线程控制GPU执行数据计算，CPU辅助线程执行数据存储，主线程与辅助线程异步执行。

2.6 杆系离散元几何非线性整体并行程序设计

图 7 GPU数据计算与数据存储模式Fig. 7 GPU data calculation and data storage mode

本研究基于CPU-GPU异构平台和CUDA计算架构开发出杆系DEM并行算法并编制了相应的计算程序。杆系DEM整体并行计算框架如图8所示。

图 8 杆系离散元并行计算架构Fig. 8 Framed system discrete element parallel computing architecture

3 杆系DEM并行算法在大型工程结构中的应用

3.1 CPU-GPU异构平台参数

本文采用的CPU与GPU硬件参数如表2、表3所示；软件开发环境为：Windows 7 64位操作系统、Visual Studio 2012软件开发环境、CUDA8.0程序包。

表 2 CPU硬件参数Table 2 CPU hardware parameters

表 3 GPU硬件参数Table 3 GPU hardware parameters

3.2 空间框架结构几何非线性动力分析

如图9所示，该钢筋混凝土框架层高为4 m，框架梁和柱截面均为矩形截面，尺寸为0.3 m×0.3 m，弹性模量E=40 GPa，泊松比为0.2，材料密度为2500 kg/m3。在结构基底沿水平x轴方向输入如图10所示的动荷载，持续时间为10 s，为模拟结构大变形将动荷载峰值加速度调至3.2g。阻尼比取ξ=5%，选用Rayleigh阻尼，阻尼计算公式可参考文献[11]。

图 9 某6层空间框架结构 /mFig. 9 A 6-story space frame structure

图11为ANSYS有限元软件计算结果与杆系DEM串、并行算法计算结果对比，可以看出三者计算得到的位移时程曲线在波形和幅值上完全吻合，表明了杆系DEM并行算法具有较高的计算精度。

为了测试杆系DEM并行算法的计算性能，将框架计算模型单元数目逐步增加，计算结果见表4，同一工况下且单元离散模型相同时，相比于杆系DEM串行算法，杆系DEM并行算法加速比最高达到了12.3倍，表明杆系DEM并行算法具有良好的加速效果。

3.3 大型球壳结构几何非线性动力分析

图 10 动荷载Fig. 10 Dynamic load

图 11 动荷载作用下A点时程-位移曲线(结构模型单元数为8865)Fig. 11 Time-history displacement curve at point A under dynamic load (the number of structural model elements is 8865)

表 4 框架结构模型计算加速比(取计算时长为0.02 s)Table 4 Frame structure model calculation acceleration ratio(the calculation time is 0.02 s)

图 12 K6型球壳模型(失跨比：f / l=1∶2) /mFig. 12 K6 spherical shell model ( f / l=1∶2)

选取K6型球壳结构为对象进行动力响应时程分析。球壳几何构形如图12所示，球壳跨度为23.4 m，矢跨比为1/2，球壳一共由3660根杆件、1261个节点组成，在球壳底部设有40个分段连续的固定支座。杆件采用Q235圆形钢管，材料弹性模量E=210 GPa，泊松比ν=0.30，材料密度7850 kg/m3，杆件规格为：最外圈支座环梁Φ114 mm×4 mm；其他杆件Φ23 mm×2 mm。所有杆件节点均采用刚性连接，每个节点附加30 kg质量块。

在ANSYS有限元软件中，采用PIPE20单元模拟杆件力学行为；采用MASS21单元模拟节点质量块。在结构基座底部沿x轴方向上施加如图13所示的正弦波，并将峰值加速度上调至0.9g，阻尼比取ξ=2%，选用Rayleigh阻尼。

图 13 正弦波Fig. 13 Sine wave

计算结果如图14所示，杆系DEM串、并行算法与ANSYS有限元软件计算结果在波形和幅值上完全吻合，再次验证了杆系DEM并行算法的计算精度。

为了测试杆系DEM并行算法的计算性能，将球壳计算模型单元数目逐步调多，测试结果见表5，可以看出杆系DEM并行算法加速比最高达到12.7倍。

图 14 正弦波作用下A点时程位移曲线(结构模型单元数为17 604)Fig. 14 Time-history displacement curve of point A under the action of sine wave (the number of structural model elements is 17 604)