-
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并行计算框架并编制了相应的几何非线性计算程序,并将其运用到大型工程结构数值计算中,获得了显著的加速效果。
-
运用杆系离散元求解结构力学问题[9],首要步骤是建立结构几何模型,主要包括确定关键点坐标和杆件连接,然后将模型离散为一系列连续的虚拟球体。如图1所示,球心为节点位置,两相邻球体球心所包含的区域为单元所在位置;以球心为中心,离散球体所包含的单元总质量等于节点的集中质量。
每个离散球体都遵循牛顿第二定律,如式(1)所示:
$$\begin{split} & \left[\!\! {\begin{array}{*{20}{c}} {{M_{\rm{a}}}}&{}&{} \\ {}&{{M_{\rm{a}}}}&{} \\ {}&{}&{{M_{\rm{a}}}} \end{array}} \right]\left\{ \!{\begin{array}{*{20}{c}} {\mathop {{\ddot{x}_1}} } \\ {\mathop {{\ddot{x}_2}} } \\ {\mathop {{\ddot{x}_3}} } \end{array}} \!\right\}=\left\{ {\begin{array}{*{20}{c}} {F_1^{{\rm{int}} }} \\ {F_2^{{\rm{int}} }} \\ {F_3^{{\rm{int}} }} \end{array}} \right\} + \left\{ {\begin{array}{*{20}{c}} {F_1^{\rm{ext}}} \\ {F_2^{\rm{ext}}} \\ {F_3^{\rm{ext}}} \end{array}} \!\!\right\}\\&\qquad \left[ {\begin{array}{*{20}{c}} {{I_{11}}}&{{I_{12}}}&{{I_{13}}} \\ {{I_{21}}}&{{I_{22}}}&{I{}_{23}} \\ {{I_{31}}}&{{I_{32}}}&{{I_{33}}} \end{array}} \right]\left\{ {\begin{array}{*{20}{c}} {\mathop {{\dot\omega _1}} } \\ {{\dot\omega _2}} \\ {\mathop {{\dot\omega _3}} } \end{array}} \right\}=\left\{ {\begin{array}{*{20}{c}} {{M_1}} \\ {{M_2}} \\ {{M_3}} \end{array}} \right\} \end{split}$$ (1) 式中:
${\ddot x_i}$ 为节点沿坐标轴方向的加速度;${\dot \omega _i}$ 为节点角加速度;$F_i^{{\rm{int}}} $ 、$F_i^{{\rm{ext}}} $ 分别为作用在节点上沿坐标轴方向的内力与外力;Mi为合力矩(Mi=$M_i^{{\rm{int}}} $ +$M_i^{{\rm{exp}}} $ );Ma为节点集中质量;Iij为单元转动惯性矩。三维杆系离散元接触本构方程基于接触力学基本理论[10],单元与单元发生接触后,将产生接触变形和接触力,根据式(2)可将接触力和接触力矩分解为沿接触平面方向的Fτ、Mτ和垂直于接触平面方向的Fn、Mn。
根据单元法向接触刚度Kn和切向接触刚度Kτ[11],按式(3)可计算出局部坐标系下单元接触内力分量,经由坐标转换可得到全局坐标下单元接触内力增量,再根据式(4)叠加增量得到下一时间步单元接触内力全量,随后进入下一时间步运动方程求解,通过牛顿第二定律与中心差分法获得节点速度、位移、转角分量,再根据单元力-位移本构关系得到单元接触力、力矩分量,并以此循环直至总时间结束。
$$\left\{ \begin{aligned} & {{{F}}}={{{{{F}}}}^{n}}+{{{{{F}}}}^{\tau }} \\& {{{M}}}={{{{{M}}}}^{n}}+{{{{{M}}}}^{\tau }} \end{aligned} \right.\;\;\;\;$$ (2) $$\left\{ \begin{aligned} & \Delta {{{{{F}}}}^{n}}=-{{K}_{n}}\Delta {{{{{U}}}}^{n}} \\& \Delta {{{{{F}}}}^{\tau }}=-{{K}_{\tau }}\Delta {{{{{U}}}}^{\tau }} \end{aligned} \right.$$ (3) $$\left\{ \begin{aligned} & {{{{F}}}^{n}}\leftarrow {{{{{F}}}}^{n}}+\Delta {{{{{F}}}}^{n}} \\& {{{{F}}}^{\tau }}\leftarrow {{{{{F}}}}^{\tau }}+\Delta {{{{{F}}}}^{\tau }} \end{aligned} \right.$$ (4) -
杆系离散元计算理论以单元和节点计算为基础,在单元和节点之上建立以时间步为基准的迭代计算,单元与节点计算几乎完全独立,并行计算特征十分鲜明。根据杆系DEM计算理论潜在的并行性,可将杆系离散元数据计算分为3大类:单元计算、节点计算、单元-节点计算,如图2所示。
对于单元与节点计算,如单元内力计算、节点位移计算等,所有参与运算的单元或者节点相互独立,具备天然并行性,完全符合GPU多线程并行计算模式。对于单元-节点计算,如节点物理量集成,所有参与计算的单元和节点之间需要进行相互转换,此时的数据计算并不独立于单元或者节点,若直接采用GPU多线程计算,则会出现线程冲突,导致数据计算错误,此时需要采用辅助算法才能正确计算。
-
杆系DEM几何非线性计算程序涉及大量多维数组,本研究对算法中所有多维数组采用降维存储,即一维数组存储多维数据。图3为杆系DEM数据存储形式,“一维数组存储多维数据”按照各维度顺序依次存储在GPU全局内存中,每一维度数据存储完成再进入下一维度数据存储,直至所有维度数据存储完成为止。
-
杆系DEM采用一维数组存储多维数据,而在多维数据的每一维度上,数据量往往是随机的,并不是32的整数倍,自然也不满足GPU线程对数据的合并访问要求[12],这将导致GPU线程对数据的访问效率大打折扣。为了满足GPU线程对数据的合并访问要求,本研究对一维数组的多维数据做对齐处理。如图4所示,采用CUDA中的cudaMallocPitch( )函数,在多维数据每一维度的最后添加一个或多个单位数据内存空间,使每一维度的数据量为32的整数倍,对齐后的一维数组可满足线程对内存的合并访问要求。
-
在杆系DEM计算理论中,每一时间步都需要将单元内力、单元转动惯性矩由单元集成到节点上,此过程涉及单元-节点转换计算,运算形式将不再独立于单元或者节点,采用GPU多线计算模式将会导致计算结果错误。
为了使GPU线程能够正确计算节点物理量,本研究提出两种解决方案。方案一是按单元集成节点物理量,即以单元计算为中心,通过单元对应节点的方式,在同一节点上累加单元所对应的物理数据,使用原子加操作[13](对应CUDA中atomicAdd( )函数)完成运算,如图5、图6(a)所示。方案二是按节点集成节点物理量,即以节点计算为中心,将不同节点对应的单元完全分开,同一节点所对应的不同单元物理量直接累加,从而使各节点之间的计算相互独立,在GPU中采用多线程循环计算,如图5、图6(b)所示。
对于上述2种计算方案,若按单元集成节点物理量,使用原子加操作将导致所有线程排列执行,因而需要耗费很长时间;若按节点集成节点物理量,需要事先建立节点-单元索引数组,当节点对应的单元数目差异较大时,GPU计算将会产生大量线程分支以及严重的负载不平衡现象,这将引起较大的计算性能损耗。
两种方案实际的GPU程序运行结果如表1所示,可以看出两者在计算时间上几乎相同。考虑到按节点集成节点物理量需要采用辅助索引数组,而创建索引数组不仅增加了内存需求,还需要额外的计算时间,所以本研究选用按单元集成节点物理量。
表 1 两种方案集成节点物理量所需的计算时间
Table 1. The calculation time required to integrate the physical quantities of the two schemes
单元数 DEM迭代步 按单元集成计算时间/s 按节点集成计算时间/s 17604 1000 0.359 0.421 36240 1000 0.734 0.780 70272 1000 1.400 1.390 -
为了获得杆系DEM计算过程中各时间子步产生的结构响应,需要在执行数个计算步后从GPU全局内存中将节点位移计算结果存储至硬盘上。传统的异构平台数据传输模式如图7(a)所示,GPU计算完成后,必须等待数据在硬盘上全部存储完成后,才能进行下一步计算。由于数据传输过程始终为串行执行,根据Amdahl定律[14],数据传输所消耗的时间将严重限制计算性能的提升。
本研究为了突破数据存储与数据计算同步执行带来的性能瓶颈,运用C++11中thread类函数接口与CUDA中cudaMemcpyAsync( )函数接口,实现了异构平台数据计算与数据存储异步执行。异步执行模式如图7(b)所示,CPU主线程控制GPU执行数据计算,CPU辅助线程执行数据存储,主线程与辅助线程异步执行。
-
本研究基于CPU-GPU异构平台和CUDA计算架构开发出杆系DEM并行算法并编制了相应的计算程序。杆系DEM整体并行计算框架如图8所示。
-
本文采用的CPU与GPU硬件参数如表2、表3所示;软件开发环境为:Windows 7 64位操作系统、Visual Studio 2012软件开发环境、CUDA8.0程序包。
表 2 CPU硬件参数
Table 2. CPU hardware parameters
型号 频率/GHz 内存容量/GB 核心数 Intel i7-8700 8 16 12 表 3 GPU硬件参数
Table 3. GPU hardware parameters
型号 显存容量/
GB显存带宽/
(GB/s)计算
能力CUDA
核心数GeForece GTX1080 8 352 6.0 2560 -
如图9所示,该钢筋混凝土框架层高为4 m,框架梁和柱截面均为矩形截面,尺寸为0.3 m×0.3 m,弹性模量E=40 GPa,泊松比为0.2,材料密度为2500 kg/m3。在结构基底沿水平x轴方向输入如图10所示的动荷载,持续时间为10 s,为模拟结构大变形将动荷载峰值加速度调至3.2 g。阻尼比取ξ=5%,选用Rayleigh阻尼,阻尼计算公式可参考文献[11]。
图11为ANSYS有限元软件计算结果与杆系DEM串、并行算法计算结果对比,可以看出三者计算得到的位移时程曲线在波形和幅值上完全吻合,表明了杆系DEM并行算法具有较高的计算精度。
图 11 动荷载作用下A点时程-位移曲线(结构模型单元数为8865)
Figure 11. Time-history displacement curve at point A under dynamic load (the number of structural model elements is 8865)
为了测试杆系DEM并行算法的计算性能,将框架计算模型单元数目逐步增加,计算结果见表4,同一工况下且单元离散模型相同时,相比于杆系DEM串行算法,杆系DEM并行算法加速比最高达到了12.3倍,表明杆系DEM并行算法具有良好的加速效果。
表 4 框架结构模型计算加速比(取计算时长为0.02 s)
Table 4. Frame structure model calculation acceleration ratio (the calculation time is 0.02 s)
单元数 CPU时间/s GPU时间/s 加速比 8865 4.8 0.59 8.1 14 541 7.3 0.74 9.8 38 654 59.6 5.18 11.5 85 648 456.0 37.68 12.1 187 591 4895.0 398.00 12.3 -
选取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.9 g,阻尼比取ξ=2%,选用Rayleigh阻尼。
计算结果如图14所示,杆系DEM串、并行算法与ANSYS有限元软件计算结果在波形和幅值上完全吻合,再次验证了杆系DEM并行算法的计算精度。
图 14 正弦波作用下A点时程位移曲线(结构模型单元数为17 604)
Figure 14. Time-history displacement curve of point A under the action of sine wave (the number of structural model elements is 17 604)
为了测试杆系DEM并行算法的计算性能,将球壳计算模型单元数目逐步调多,测试结果见表5,可以看出杆系DEM并行算法加速比最高达到12.7倍。
表 5 球壳结构模型计算加速比(取计算时长为0.02 s)
Table 5. Calculated acceleration ratio of spherical shell structure model (the calculation time is 0.02 s)
单元数 CPU时间/s GPU时间/s 加速比 17 604 7.8 0.71 10.9 36 240 42.8 3.66 11.6 70 272 344.0 28.00 12.3 172 806 4302.0 339.00 12.7 -
本文提出杆系DEM的GPU并行算法,并给出了并行算法的设计过程,最后通过算例验证了并行算法的计算精度和计算效率,得到如下结论:
(1) 杆系DEM计算理论以单元和节点计算为基础,单元与节点计算几乎完全独立,具备天然并行性。本研究针对杆系DEM计算理论提出单元级并行、节点级并行的计算方法,将单元、节点与GPU多线程建立映射关系,实现了杆系DEM的GPU多线程并行计算。
(2) 基于CPU-GPU异构平台建立了杆系DEM并行计算框架,并编制相应的计算程序。对杆系DEM并行算法的设计主要包括数据存储方式、GPU线程计算模式、节点物理量集成方式以及数据传输优化。
(3) 采用大型三维框架、球壳结构模型分别验证了杆系DEM并行算法的计算精度,并对杆系DEM并行算法进行了计算性能测试,测试结果表明杆系DEM并行算法加速比最高可达到12.7倍,并具有优良的计算精度。
APPLICATION OF GPU-BASED PARALLEL COMPUTING METHOD FOR DEM IN LARGE ENGINEERING STRUCTURES
-
摘要: 杆系DEM(离散元,discrete element method)是求解结构强非线性问题的有效方法,但随着结构数值计算规模的扩大,杆系DEM所需要的计算时间也随之急剧膨胀。为了提高杆系DEM的计算效率,该研究提出单元级并行、节点级并行的计算方法,基于CPU-GPU异构平台,建构了杆系DEM并行计算框架,编制了相应的几何非线性计算程序,实现了杆系DEM的GPU多线程并行计算。对杆系DEM并行算法的设计主要包括数据存储方式、GPU线程计算模式、节点物理量集成方式以及数据传输优化。最后采用大型三维框架、球壳结构模型分别验证了杆系DEM并行算法的计算精度,并对杆系DEM并行算法进行了计算性能测试,测试结果表明杆系DEM并行算法加速比最高可达12.7倍。
-
关键词:
- 离散单元法 /
- 杆系结构 /
- 几何非线性 /
- GPU并行计算 /
- CPU-GPU异构平台
Abstract: The member DEM (discrete element method) is an effective way for the structural analyses concerning with strong nonlinear issues. However, with the expansion of structural numerical calculation model, the member DEM program has been time-consuming dramatically. It proposes element-level parallel and node-level parallel computing methods to improve the calculation efficiency of the member DEM. A parallel computing framework for the member DEM is first developed based on a CPU-GPU heterogeneous platform. Then, the member DEM program that models structural geometric nonlinearity is compiled and embedded into the framework, and finally the GPU-based multi-thread parallel algorithm for the member DEM is established. The design of this parallel algorithm mainly composed of a data storage mode, a GPU thread computing mode, a node physical quantity integration mode and the way of data transmission optimization. The accuracy of the member DEM parallel algorithm proposed in this study is verified by a large-scale space frame model and a spherical shell structure model, and the performance of the algorithm was further tested. The results show that the speedup of the member DEM parallel algorithm can reach 12.7 times at most. -
表 1 两种方案集成节点物理量所需的计算时间
Table 1. The calculation time required to integrate the physical quantities of the two schemes
单元数 DEM迭代步 按单元集成计算时间/s 按节点集成计算时间/s 17604 1000 0.359 0.421 36240 1000 0.734 0.780 70272 1000 1.400 1.390 表 2 CPU硬件参数
Table 2. CPU hardware parameters
型号 频率/GHz 内存容量/GB 核心数 Intel i7-8700 8 16 12 表 3 GPU硬件参数
Table 3. GPU hardware parameters
型号 显存容量/
GB显存带宽/
(GB/s)计算
能力CUDA
核心数GeForece GTX1080 8 352 6.0 2560 表 4 框架结构模型计算加速比(取计算时长为0.02 s)
Table 4. Frame structure model calculation acceleration ratio (the calculation time is 0.02 s)
单元数 CPU时间/s GPU时间/s 加速比 8865 4.8 0.59 8.1 14 541 7.3 0.74 9.8 38 654 59.6 5.18 11.5 85 648 456.0 37.68 12.1 187 591 4895.0 398.00 12.3 表 5 球壳结构模型计算加速比(取计算时长为0.02 s)
Table 5. Calculated acceleration ratio of spherical shell structure model (the calculation time is 0.02 s)
单元数 CPU时间/s GPU时间/s 加速比 17 604 7.8 0.71 10.9 36 240 42.8 3.66 11.6 70 272 344.0 28.00 12.3 172 806 4302.0 339.00 12.7 -
[1] Feng Z, Xu W, Lubbe R. Three-dimensional morphological characteristics of particles in nature and its application for DEM simulation [J]. Powder Technology, 2020, 364: 635 − 641. doi: 10.1016/j.powtec.2020.02.022 [2] Xua J, Qi H, Fang X, et al. Quasi-real-time simulation of rotating drum using discrete element method with parallel GPU computing [J]. Particuology, 2011, 9(4): 446 − 450. doi: 10.1016/j.partic.2011.01.003 [3] He Y, Evans T J, Yu A B, et al. A GPU-based DEM for modelling large scale powder compaction with wide size distributions [J]. Powder Technology, 2018, 333: 219 − 228. doi: 10.1016/j.powtec.2018.04.034 [4] Chase Cook, Hengyang Zhao, Takashi Sato, et al. GPU-based Ising computing for solving max-cut combinatorial optimization problems [J]. Integration, 2019, 69: 335 − 344. doi: 10.1016/j.vlsi.2019.07.003 [5] Spellings M, Marson R L, Anderson J A, et al. GPU accelerated discrete element method (DEM) molecular dynamics for conservative, faceted particle simulations [J]. Journal of Computational Physics, 2017, 334: 460 − 467. doi: 10.1016/j.jcp.2017.01.014 [6] Liu G, Xu W, Sun Q, et al. Study on the particle breakage of ballast based on a GPU accelerated discrete element method [J]. Geoscience Frontiers, 2020, 11(2): 461 − 471. doi: 10.1016/j.gsf.2019.06.006 [7] 付帅旗, 黄鹏, 丁逸飞. 基于DEM和GPU加速的颗粒运动仿真方法研究[J]. 合肥工业大学学报(自然科学版), 2019, 42(12): 1602 − 1607. Fu Shuaiqi, Huang Peng, Ding Yifei. Research on particle motion simulation method based on DEM and GPU acceleration [J]. Journal of Hefei University of Technology (Natural Science Edition), 2019, 42(12): 1602 − 1607. (in Chinese) [8] Ma Z S, Feng C, Liu T P, et al. A GPU accelerated continuous-based discrete element method for elastodynamics analysis [J]. Advanced Materials Research, 2011, 320: 329 − 334. doi: 10.4028/www.scientific.net/AMR.320.329 [9] Jihong Y, Nian Q. Progressive collapse simulation based on DEM for single-layer reticulated domes [J]. Journal of Constructional Steel Research, 2017, 128(Jan): 721 − 731. [10] Xu Ye. DEM algorithm for progressive collapse simulation of single-layer reticulated domes under multi-support excitation [J]. Journal of Earthquake Engineering, 2019, 23(1): 18 − 45. doi: 10.1080/13632469.2017.1309606 [11] 齐念, 叶继红. 弹性DEM方法在杆系结构中的应用研究[J]. 工程力学, 2017, 34(7): 11 − 20. Qi Nian, Ye Jihong. Application research of elastic DEM method in rod structure [J]. Engineering Mechanics, 2017, 34(7): 11 − 20. (in Chinese) [12] Honglin Huang. Research on the computer graphics rendering technology based on GPU and parallel computing system [C]. Proceedings of 2016 2nd Workshop on Advanced Research and Technology in Industry Applications (WARTIA 2016). Dalian, Liaoning, China, 2016: 910 − 914. [13] Gribanov I, Taylor R, Sarracino R. Parallel implementation of implicit finite element model with cohesive zones and collision response using CUDA [J]. International Journal for Numerical Methods in Engineering, 2018, 115(7): 771 − 790. doi: 10.1002/nme.5825 [14] Owens J D, Houston M, Luebke D, et al. GPU computing [J]. Proceedings of the IEEE, 2008, 96(5): 879 − 899. doi: 10.1109/JPROC.2008.917757 -