葉繼紅，王 佳

(1. 中國(guó)礦業(yè)大學(xué)江蘇省土木工程環(huán)境災(zāi)變與結(jié)構(gòu)可靠性重點(diǎn)實(shí)驗(yàn)室，徐州 221116；2. 中國(guó)礦業(yè)大學(xué)徐州市工程結(jié)構(gòu)火災(zāi)安全重點(diǎn)實(shí)驗(yàn)室，徐州 221116)

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

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

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

在結(jié)構(gòu)工程領(lǐng)域，桿系DEM是求解結(jié)構(gòu)強(qiáng)非線性問(wèn)題的有效方法，相比于巖土工程領(lǐng)域廣泛應(yīng)用的CDEM計(jì)算方法，桿系DEM在單元變形計(jì)算與單元內(nèi)力求解上更為復(fù)雜，且數(shù)據(jù)計(jì)算復(fù)雜度更高。因此相比于CDEM并行算法，桿系DEM并行算法的設(shè)計(jì)難度更大。

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

1 桿系離散元幾何非線性計(jì)算理論——運(yùn)動(dòng)方程的建立與求解

運(yùn)用桿系離散元求解結(jié)構(gòu)力學(xué)問(wèn)題[9]，首要步驟是建立結(jié)構(gòu)幾何模型，主要包括確定關(guān)鍵點(diǎn)坐標(biāo)和桿件連接，然后將模型離散為一系列連續(xù)的虛擬球體。如圖1所示，球心為節(jié)點(diǎn)位置，兩相鄰球體球心所包含的區(qū)域?yàn)閱卧谖恢?；以球心為中心，離散球體所包含的單元總質(zhì)量等于節(jié)點(diǎn)的集中質(zhì)量。

圖 1 空間桿系結(jié)構(gòu)離散模型Fig. 1 Discrete model of space framed structure

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

2 基于異構(gòu)平臺(tái)桿系離散元幾何非線性并行算法設(shè)計(jì)

2.1 桿系離散元計(jì)算理論并行性分析

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

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

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

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

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

2.3 桿系DEM線程計(jì)算模式

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

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

在桿系DEM計(jì)算理論中，每一時(shí)間步都需要將單元內(nèi)力、單元轉(zhuǎn)動(dòng)慣性矩由單元集成到節(jié)點(diǎn)上，此過(guò)程涉及單元-節(jié)點(diǎn)轉(zhuǎn)換計(jì)算，運(yùn)算形式將不再獨(dú)立于單元或者節(jié)點(diǎn)，采用GPU多線計(jì)算模式將會(huì)導(dǎo)致計(jì)算結(jié)果錯(cuò)誤。

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

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

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

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

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

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

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

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

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

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

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

2.6 桿系離散元幾何非線性整體并行程序設(shè)計(jì)

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

本研究基于CPU-GPU異構(gòu)平臺(tái)和CUDA計(jì)算架構(gòu)開(kāi)發(fā)出桿系DEM并行算法并編制了相應(yīng)的計(jì)算程序。桿系DEM整體并行計(jì)算框架如圖8所示。

圖 8 桿系離散元并行計(jì)算架構(gòu)Fig. 8 Framed system discrete element parallel computing architecture

3 桿系DEM并行算法在大型工程結(jié)構(gòu)中的應(yīng)用

3.1 CPU-GPU異構(gòu)平臺(tái)參數(shù)

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

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

表 3 GPU硬件參數(shù)Table 3 GPU hardware parameters

3.2 空間框架結(jié)構(gòu)幾何非線性動(dòng)力分析

如圖9所示，該鋼筋混凝土框架層高為4 m，框架梁和柱截面均為矩形截面，尺寸為0.3 m×0.3 m，彈性模量E=40 GPa，泊松比為0.2，材料密度為2500 kg/m3。在結(jié)構(gòu)基底沿水平x軸方向輸入如圖10所示的動(dòng)荷載，持續(xù)時(shí)間為10 s，為模擬結(jié)構(gòu)大變形將動(dòng)荷載峰值加速度調(diào)至3.2g。阻尼比取ξ=5%，選用Rayleigh阻尼，阻尼計(jì)算公式可參考文獻(xiàn)[11]。

圖 9 某6層空間框架結(jié)構(gòu) /mFig. 9 A 6-story space frame structure

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

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

3.3 大型球殼結(jié)構(gòu)幾何非線性動(dòng)力分析

圖 10 動(dòng)荷載Fig. 10 Dynamic load

圖 11 動(dòng)荷載作用下A點(diǎn)時(shí)程-位移曲線(結(jié)構(gòu)模型單元數(shù)為8865)Fig. 11 Time-history displacement curve at point A under dynamic load (the number of structural model elements is 8865)

表 4 框架結(jié)構(gòu)模型計(jì)算加速比(取計(jì)算時(shí)長(zhǎng)為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型球殼結(jié)構(gòu)為對(duì)象進(jìn)行動(dòng)力響應(yīng)時(shí)程分析。球殼幾何構(gòu)形如圖12所示，球殼跨度為23.4 m，矢跨比為1/2，球殼一共由3660根桿件、1261個(gè)節(jié)點(diǎn)組成，在球殼底部設(shè)有40個(gè)分段連續(xù)的固定支座。桿件采用Q235圓形鋼管，材料彈性模量E=210 GPa，泊松比ν=0.30，材料密度7850 kg/m3，桿件規(guī)格為：最外圈支座環(huán)梁Φ114 mm×4 mm；其他桿件Φ23 mm×2 mm。所有桿件節(jié)點(diǎn)均采用剛性連接，每個(gè)節(jié)點(diǎn)附加30 kg質(zhì)量塊。

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

圖 13 正弦波Fig. 13 Sine wave

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

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

圖 14 正弦波作用下A點(diǎn)時(shí)程位移曲線(結(jié)構(gòu)模型單元數(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 球殼結(jié)構(gòu)模型計(jì)算加速比(取計(jì)算時(shí)長(zhǎng)為0.02 s)Table 5 Calculated acceleration ratio of spherical shell structure model (the calculation time is 0.02 s)

4 結(jié)論

本文提出桿系DEM的GPU并行算法，并給出了并行算法的設(shè)計(jì)過(guò)程，最后通過(guò)算例驗(yàn)證了并行算法的計(jì)算精度和計(jì)算效率，得到如下結(jié)論：

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

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

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