基于GPU的桿系離散元并行算法在大型工程結構中的應用

葉繼紅，王 佳

(1. 中國礦業(yè)大學江蘇省土木工程環(huán)境災變與結構可靠性重點實驗室，徐州 221116；2. 中國礦業(yè)大學徐州市工程結構火災安全重點實驗室，徐州 221116)

20世紀70年代，F(xiàn)eng等[1]率先提出離散單元法(discrete element method，DEM)，隨后在巖土工程領域獲得了廣泛應用。DEM屬于顯式求解算法，DEM允許單元間發(fā)生相對運動，不需要滿足位移連續(xù)條件和變形協(xié)調條件，計算中不需要集成整體剛度矩陣，相比于傳統(tǒng)有限單元法，求解計算不存在收斂性問題，非常適合求解強非線性問題。盡管離散單元法己被證明是求解強非線性問題的有效方法，并且廣泛運用在各種研究領域，更是擁有傳統(tǒng)有限單元法難以企及的獨特優(yōu)勢，但是DEM計算效率普遍偏低一直是一個亟待解決的問題。

目前，在并行計算領域基于多處理器并行技術[2]與CUDA(compute unified device architecture)計算架構是解決大規(guī)模離散元數(shù)值計算的有效途徑。對于多處理器并行技術而言，要實現(xiàn)多機協(xié)同工作，針對解決DEM這種全局問題，編程難度偏大；而基于CUDA的GPU(graphic processing unit)高性能計算，近幾年應用領域則越來越廣，例如在DEM相關領域，分子動力學、地質工程、原子間作用等方面均得到了廣泛應用。

在國際上，He等[3]開發(fā)出基于GPU的離散單元法，用于大規(guī)模粉末壓實模擬，大幅縮短了仿真計算時間。Chase等[4]將GPU并行技術運用到Yade開放式DEM軟件中，運用GPU多線程技術加速矩陣因子分解，將多孔彈性滲流數(shù)值分析效率提高了170倍。 Spellings等[5]運用GPU加速DEM算法，用于模擬重力作用下各向異性顆粒系統(tǒng)，獲得了比較顯著的加速效果。Liu等[6]采用GPU加速Blaze-DEM算法，并將其運用到大規(guī)模顆粒材料破損分析中。在國內，付帥旗等[7]充分利用GPU眾核多線程的計算優(yōu)勢，實現(xiàn)了大規(guī)模顆粒離散元接觸模擬。在基于大規(guī)模連續(xù)介質離散元(continuum-based discrete element method，CDEM)計算方面，中科院力學研究所的研究表明[8]，在CDEM中使用GPU并行技術可以使計算效率提高100倍以上，GPU并行架構的優(yōu)越性不言而喻。

在結構工程領域，桿系DEM是求解結構強非線性問題的有效方法，相比于巖土工程領域廣泛應用的CDEM計算方法，桿系DEM在單元變形計算與單元內力求解上更為復雜，且數(shù)據計算復雜度更高。因此相比于CDEM并行算法，桿系DEM并行算法的設計難度更大。

為了設計高效的桿系DEM并行算法，本研究提出單元級并行、節(jié)點級并行的計算方法，并對桿系DEM的數(shù)據存儲方式、GPU線程計算模式、節(jié)點物理量集成方式以及數(shù)據傳輸模式進行了詳細設計?；贑PU-GPU異構平臺建立了桿系DEM并行計算框架并編制了相應的幾何非線性計算程序，并將其運用到大型工程結構數(shù)值計算中，獲得了顯著的加速效果。

1 桿系離散元幾何非線性計算理論——運動方程的建立與求解

運用桿系離散元求解結構力學問題[9]，首要步驟是建立結構幾何模型，主要包括確定關鍵點坐標和桿件連接，然后將模型離散為一系列連續(xù)的虛擬球體。如圖1所示，球心為節(jié)點位置，兩相鄰球體球心所包含的區(qū)域為單元所在位置；以球心為中心，離散球體所包含的單元總質量等于節(jié)點的集中質量。

圖 1 空間桿系結構離散模型Fig. 1 Discrete model of space framed structure

三維桿系離散元接觸本構方程基于接觸力學基本理論[10]，單元與單元發(fā)生接觸后，將產生接觸變形和接觸力，根據式(2)可將接觸力和接觸力矩分解為沿接觸平面方向的Fτ、Mτ和垂直于接觸平面方向的Fn、Mn。

2 基于異構平臺桿系離散元幾何非線性并行算法設計

2.1 桿系離散元計算理論并行性分析

桿系離散元計算理論以單元和節(jié)點計算為基礎，在單元和節(jié)點之上建立以時間步為基準的迭代計算，單元與節(jié)點計算幾乎完全獨立，并行計算特征十分鮮明。根據桿系DEM計算理論潛在的并行性，可將桿系離散元數(shù)據計算分為3大類：單元計算、節(jié)點計算、單元-節(jié)點計算，如圖2所示。

圖 2 桿系離散元并行性分析Fig. 2 Parallelism analysis of discrete elements

對于單元與節(jié)點計算，如單元內力計算、節(jié)點位移計算等，所有參與運算的單元或者節(jié)點相互獨立，具備天然并行性，完全符合GPU多線程并行計算模式。對于單元-節(jié)點計算，如節(jié)點物理量集成，所有參與計算的單元和節(jié)點之間需要進行相互轉換，此時的數(shù)據計算并不獨立于單元或者節(jié)點，若直接采用GPU多線程計算，則會出現(xiàn)線程沖突，導致數(shù)據計算錯誤，此時需要采用輔助算法才能正確計算。

2.2 桿系DEM數(shù)據存儲方式

桿系DEM幾何非線性計算程序涉及大量多維數(shù)組，本研究對算法中所有多維數(shù)組采用降維存儲，即一維數(shù)組存儲多維數(shù)據。圖3為桿系DEM數(shù)據存儲形式，“一維數(shù)組存儲多維數(shù)據”按照各維度順序依次存儲在GPU全局內存中，每一維度數(shù)據存儲完成再進入下一維度數(shù)據存儲，直至所有維度數(shù)據存儲完成為止。

2.3 桿系DEM線程計算模式

桿系DEM采用一維數(shù)組存儲多維數(shù)據，而在多維數(shù)據的每一維度上，數(shù)據量往往是隨機的，并不是32的整數(shù)倍，自然也不滿足GPU線程對數(shù)據的合并訪問要求[12]，這將導致GPU線程對數(shù)據的訪問效率大打折扣。為了滿足GPU線程對數(shù)據的合并訪問要求，本研究對一維數(shù)組的多維數(shù)據做對齊處理。如圖4所示，采用CUDA中的cudaMallocPitch( )函數(shù)，在多維數(shù)據每一維度的最后添加一個或多個單位數(shù)據內存空間，使每一維度的數(shù)據量為32的整數(shù)倍，對齊后的一維數(shù)組可滿足線程對內存的合并訪問要求。

2.4 節(jié)點物理量集成方式

在桿系DEM計算理論中，每一時間步都需要將單元內力、單元轉動慣性矩由單元集成到節(jié)點上，此過程涉及單元-節(jié)點轉換計算，運算形式將不再獨立于單元或者節(jié)點，采用GPU多線計算模式將會導致計算結果錯誤。

圖 3 桿系DEM數(shù)據存儲Fig. 3 Framed system DEM data storage

圖 4 桿系DEM數(shù)據對齊Fig. 4 Framed system DEM data alignment

為了使GPU線程能夠正確計算節(jié)點物理量，本研究提出兩種解決方案。方案一是按單元集成節(jié)點物理量，即以單元計算為中心，通過單元對應節(jié)點的方式，在同一節(jié)點上累加單元所對應的物理數(shù)據，使用原子加操作[13](對應CUDA中atomicAdd( )函數(shù))完成運算，如圖5、圖6(a)所示。方案二是按節(jié)點集成節(jié)點物理量，即以節(jié)點計算為中心，將不同節(jié)點對應的單元完全分開，同一節(jié)點所對應的不同單元物理量直接累加，從而使各節(jié)點之間的計算相互獨立，在GPU中采用多線程循環(huán)計算，如圖5、圖6(b)所示。

圖 5 桿系DEM單元連接模型Fig. 5 The connection model of the DEM unit of the framed system

圖 6 節(jié)點物理量集成方式Fig. 6 Node physical quantity integration method

對于上述2種計算方案，若按單元集成節(jié)點物理量，使用原子加操作將導致所有線程排列執(zhí)行，因而需要耗費很長時間；若按節(jié)點集成節(jié)點物理量，需要事先建立節(jié)點-單元索引數(shù)組，當節(jié)點對應的單元數(shù)目差異較大時，GPU計算將會產生大量線程分支以及嚴重的負載不平衡現(xiàn)象，這將引起較大的計算性能損耗。

兩種方案實際的GPU程序運行結果如表1所示，可以看出兩者在計算時間上幾乎相同?？紤]到按節(jié)點集成節(jié)點物理量需要采用輔助索引數(shù)組，而創(chuàng)建索引數(shù)組不僅增加了內存需求，還需要額外的計算時間，所以本研究選用按單元集成節(jié)點物理量。

表 1 兩種方案集成節(jié)點物理量所需的計算時間Table 1 The calculation time required to integrate the physical quantities of the two schemes

2.5 多線程數(shù)據傳輸優(yōu)化

為了獲得桿系DEM計算過程中各時間子步產生的結構響應，需要在執(zhí)行數(shù)個計算步后從GPU全局內存中將節(jié)點位移計算結果存儲至硬盤上。傳統(tǒng)的異構平臺數(shù)據傳輸模式如圖7(a)所示，GPU計算完成后，必須等待數(shù)據在硬盤上全部存儲完成后，才能進行下一步計算。由于數(shù)據傳輸過程始終為串行執(zhí)行，根據Amdahl定律[14]，數(shù)據傳輸所消耗的時間將嚴重限制計算性能的提升。

本研究為了突破數(shù)據存儲與數(shù)據計算同步執(zhí)行帶來的性能瓶頸，運用C++11中thread類函數(shù)接口與CUDA中cudaMemcpyAsync( )函數(shù)接口，實現(xiàn)了異構平臺數(shù)據計算與數(shù)據存儲異步執(zhí)行。異步執(zhí)行模式如圖7(b)所示，CPU主線程控制GPU執(zhí)行數(shù)據計算，CPU輔助線程執(zhí)行數(shù)據存儲，主線程與輔助線程異步執(zhí)行。

2.6 桿系離散元幾何非線性整體并行程序設計

圖 7 GPU數(shù)據計算與數(shù)據存儲模式Fig. 7 GPU data calculation and data storage mode

本研究基于CPU-GPU異構平臺和CUDA計算架構開發(fā)出桿系DEM并行算法并編制了相應的計算程序。桿系DEM整體并行計算框架如圖8所示。

圖 8 桿系離散元并行計算架構Fig. 8 Framed system discrete element parallel computing architecture

3 桿系DEM并行算法在大型工程結構中的應用

3.1 CPU-GPU異構平臺參數(shù)

本文采用的CPU與GPU硬件參數(shù)如表2、表3所示；軟件開發(fā)環(huán)境為：Windows 7 64位操作系統(tǒng)、Visual Studio 2012軟件開發(fā)環(huán)境、CUDA8.0程序包。

表 2 CPU硬件參數(shù)Table 2 CPU hardware parameters

表 3 GPU硬件參數(shù)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所示的動荷載，持續(xù)時間為10 s，為模擬結構大變形將動荷載峰值加速度調至3.2g。阻尼比取ξ=5%，選用Rayleigh阻尼，阻尼計算公式可參考文獻[11]。

圖 9 某6層空間框架結構 /mFig. 9 A 6-story space frame structure

圖11為ANSYS有限元軟件計算結果與桿系DEM串、并行算法計算結果對比，可以看出三者計算得到的位移時程曲線在波形和幅值上完全吻合，表明了桿系DEM并行算法具有較高的計算精度。

為了測試桿系DEM并行算法的計算性能，將框架計算模型單元數(shù)目逐步增加，計算結果見表4，同一工況下且單元離散模型相同時，相比于桿系DEM串行算法，桿系DEM并行算法加速比最高達到了12.3倍，表明桿系DEM并行算法具有良好的加速效果。

3.3 大型球殼結構幾何非線性動力分析

圖 10 動荷載Fig. 10 Dynamic load

圖 11 動荷載作用下A點時程-位移曲線(結構模型單元數(shù)為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個節(jié)點組成，在球殼底部設有40個分段連續(xù)的固定支座。桿件采用Q235圓形鋼管，材料彈性模量E=210 GPa，泊松比ν=0.30，材料密度7850 kg/m3，桿件規(guī)格為：最外圈支座環(huán)梁Φ114 mm×4 mm；其他桿件Φ23 mm×2 mm。所有桿件節(jié)點均采用剛性連接，每個節(jié)點附加30 kg質量塊。

在ANSYS有限元軟件中，采用PIPE20單元模擬桿件力學行為；采用MASS21單元模擬節(jié)點質量塊。在結構基座底部沿x軸方向上施加如圖13所示的正弦波，并將峰值加速度上調至0.9g，阻尼比取ξ=2%，選用Rayleigh阻尼。

圖 13 正弦波Fig. 13 Sine wave

計算結果如圖14所示，桿系DEM串、并行算法與ANSYS有限元軟件計算結果在波形和幅值上完全吻合，再次驗證了桿系DEM并行算法的計算精度。

為了測試桿系DEM并行算法的計算性能，將球殼計算模型單元數(shù)目逐步調多，測試結果見表5，可以看出桿系DEM并行算法加速比最高達到12.7倍。

圖 14 正弦波作用下A點時程位移曲線(結構模型單元數(shù)為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)

表 5 球殼結構模型計算加速比(取計算時長為0.02 s)Table 5 Calculated acceleration ratio of spherical shell structure model (the calculation time is 0.02 s)

4 結論

本文提出桿系DEM的GPU并行算法，并給出了并行算法的設計過程，最后通過算例驗證了并行算法的計算精度和計算效率，得到如下結論：

(1) 桿系DEM計算理論以單元和節(jié)點計算為基礎，單元與節(jié)點計算幾乎完全獨立，具備天然并行性。本研究針對桿系DEM計算理論提出單元級并行、節(jié)點級并行的計算方法，將單元、節(jié)點與GPU多線程建立映射關系，實現(xiàn)了桿系DEM的GPU多線程并行計算。

(2) 基于CPU-GPU異構平臺建立了桿系DEM并行計算框架，并編制相應的計算程序。對桿系DEM并行算法的設計主要包括數(shù)據存儲方式、GPU線程計算模式、節(jié)點物理量集成方式以及數(shù)據傳輸優(yōu)化。

(3) 采用大型三維框架、球殼結構模型分別驗證了桿系DEM并行算法的計算精度，并對桿系DEM并行算法進行了計算性能測試，測試結果表明桿系DEM并行算法加速比最高可達到12.7倍，并具有優(yōu)良的計算精度。