肖永浩，莫則堯

(1.中國工程物理研究院 計算機應(yīng)用研究所，四川 綿陽 621900;2北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100088)

接觸問題廣泛存在于汽車碰撞，金屬成型等工程應(yīng)用和現(xiàn)實生活中，對接觸問題的數(shù)值模擬對于安全問題、武器模擬等非常重要。作為非線性大變形問題，接觸問題的并行計算受到廣泛關(guān)注。國外接觸問題及其并行算法已有多年研究，PRONTO3D［1］，ALE3D［2］等程序都可進行接觸問題的大規(guī)模并行計算，其并行算法值得借鑒。國內(nèi)主要依賴DYNA、ABAQUS等商業(yè)軟件，寇哲君［3］，亓文果［4］等在集群環(huán)境下進行了汽車碰撞的小規(guī)模并行模擬。

并行程序運行于并行機上，并行機的體系結(jié)構(gòu)直接影響并行程序的性能。當(dāng)前，節(jié)點內(nèi)部多處理器共享內(nèi)存而節(jié)點間分布式內(nèi)存的多處理器集群(CLUsters of MultiProcessors，CLUMPs)成為主流。在多處理器集群存在兩種編程模式，稱之為單一編程模式和混合編程模式［5］。單一編程模式采用單一的API描述多處理器節(jié)點間和節(jié)點內(nèi)部的通信，由于并行計算主要通過處理器間的消息傳遞(如 Message Passing Interface，MPI)來完成，因此這種模式也稱為消息傳遞模式?；旌暇幊棠Ｊ綄⒐?jié)點內(nèi)的共享存儲模式與節(jié)點間的消息傳遞結(jié)合起來。單一編程模式可采用MPI實現(xiàn)，混合編程模式可采用MPI+OpenMP實現(xiàn)。

雖然在大部分情形下，基于MPI實現(xiàn)的單一模式的并行程序的性能要優(yōu)于混合并行模式的程序［5］，但是在諸如全局通信量很大，內(nèi)存受限等條件下，混合并行模式的程序在CLUMPs系統(tǒng)上的性能要好于單一模式的程序。本文在接觸計算程序IPDC的基礎(chǔ)上，根據(jù)接觸計算的特點，采用混合編程模式實現(xiàn)了接觸問題的并行計算。

1 接觸問題

接觸計算一般采用Lagrangian有限元方法，其基本方程通過有限元離散后寫成矩陣形式為:

時間方向采用顯式時間積分，即對式(1)，若已知第n步的節(jié)點力，第n+1步的節(jié)點位移可以用下式計算:

在一個時間步內(nèi)，計算過程分為內(nèi)力計算和接觸力計算。其中內(nèi)力計算通過計算單元應(yīng)變率，應(yīng)用本構(gòu)方程計算單元應(yīng)力，最終通過對單元節(jié)點的力裝配得到節(jié)點內(nèi)力。接觸力計算用于更新接觸點的接觸力，其計算過程分為全局搜索，局部搜索和接觸力更新，全局搜索確定離接觸面片最近的接觸點，局部搜索確定接觸點在接觸面片上具體的穿透點，最后更新接觸力。接觸問題的計算流程如圖1所示。

圖1 接觸計算流程Fig.1 Outline of contact problem solving

接觸力計算中全局搜索是最耗時部分，point-in-box(捕捉盒)［1］、bucket-sorting(桶排序)［6］等算法用于加速接觸點之間臨近關(guān)系的搜索，其原理都是通過對搜索空間的劃分，減少全局搜索量，本文程序采用桶排序算法。而接觸力更新算法有罰函數(shù)法，動態(tài)約束法，分布參數(shù)法等［7］，本文程序采用罰函數(shù)法，即根據(jù)接觸點穿透接觸面片的穿透深度確定接觸點所受到的接觸力。

2 雙重區(qū)域剖分算法

根據(jù)接觸問題算法，時間步之間是嚴(yán)格順序依賴關(guān)系，并行計算只能在時間步內(nèi)進行。在一個時間步內(nèi)，內(nèi)力計算和接觸力計算的計算對象是不同的，內(nèi)力計算針對的是模型中所有的有限元網(wǎng)格，包括單元和節(jié)點，而接觸力計算只針對模型外表面的面片和構(gòu)成面片的節(jié)點。并且內(nèi)力計算和接觸力計算的對象間數(shù)據(jù)依賴性質(zhì)不同，內(nèi)力計算中數(shù)據(jù)依賴基于有限元網(wǎng)格拓撲連接，而接觸力計算基于接觸點的空間鄰近關(guān)系。為了達到兩部分計算的負載平衡，兩部分的并行計算采取兩種區(qū)域剖分算法，稱之為雙重區(qū)域剖分算法。

在碰撞計算過程中單元拓撲連接關(guān)系不變，因此內(nèi)力計算區(qū)域只需要剖分一次，達到靜態(tài)負載平衡即可。但是接觸力計算中接觸點之間的空間臨近關(guān)系是時變的，易引起負載的動態(tài)不平衡。雙重區(qū)域剖分算法不僅能較好的處理兩種計算區(qū)域不同性質(zhì)的剖分問題，而且可以通過對接觸力計算區(qū)域的動態(tài)重剖分實現(xiàn)動態(tài)負載平衡。

接觸問題的雙重區(qū)域剖分算法步驟如下:

(1)內(nèi)力計算區(qū)域剖分。

(2)內(nèi)力 fint并行計算。

(3)接觸力計算區(qū)域剖分。

(4)接觸力并行計算。

(5)節(jié)點力更新，即計算f=fext－ fint。

(6)節(jié)點加速度、速度、位移更新，轉(zhuǎn)至第(2)步。

2.1 內(nèi)力并行算法

內(nèi)力并行計算是直接的，在一個顯式時間步中，每個單元只和其有限元網(wǎng)格拓撲關(guān)系上相鄰單元有關(guān)系。每個內(nèi)力計算子區(qū)域只需和這個子區(qū)域邊界單元相鄰單元所在的子域進行通信。

具體的，在每個時間步，單元應(yīng)力、應(yīng)變等的計算只需組成本單元節(jié)點的速度，位移等信息，可以并行更新。節(jié)點內(nèi)力是由所有與這個節(jié)點相連接的單元應(yīng)力裝配得到的，因此，子區(qū)域間共享節(jié)點的節(jié)點內(nèi)力裝配需要擁有這個節(jié)點所有子區(qū)域得到的節(jié)點內(nèi)力相加得到。

2.2 接觸力并行算法

接觸力計算中，局部搜索和接觸力更新是局部計算，并行算法的核心是并行全局搜索算法，本文采用基于并行桶排序算法的全局搜索算法。假設(shè)已有分布于多個處理器上的內(nèi)力計算子區(qū)域，而接觸點是分布在各個內(nèi)力計算子區(qū)域中的，接觸力并行計算流程如下:

(1)根據(jù)內(nèi)力計算子區(qū)域中接觸點坐標(biāo)構(gòu)造全局桶。

(2)基于全局桶進行區(qū)域剖分，得到接觸力計算子區(qū)域。

(3)內(nèi)力計算子區(qū)域確定接觸點所屬接觸力計算子區(qū)域，得到內(nèi)力計算子區(qū)域和接觸力計算子區(qū)域間通信關(guān)系。

(4)內(nèi)力計算子區(qū)域發(fā)送接觸點信息到接觸力計算子區(qū)域。

(5)局部搜索和接觸力更新。

(6)接觸力計算子區(qū)域發(fā)送接觸力至內(nèi)力計算子區(qū)域。

3 并行實現(xiàn)

3.1 MPI實現(xiàn)

內(nèi)力計算的區(qū)域剖分采用圖剖分(Metis)實現(xiàn)，即將網(wǎng)格抽象為圖，圖的頂點代表單元，圖的邊代表單元間存在公共面，基于不同材料而導(dǎo)致單元計算量的不同，可通過施加圖中頂點的權(quán)值來達到負載平衡。

接觸力計算時，接觸點空間坐標(biāo)的時變導(dǎo)致接觸點間空間臨近關(guān)系未知，因此接觸力計算的區(qū)域剖分面向并行桶排序。得到全局桶后，將每個桶抽象為圖的頂點，若兩個桶相鄰，則相應(yīng)的頂點間存在一條邊，由此得到的圖是接觸點間臨近關(guān)系圖的一個超圖。每個桶中接觸點的數(shù)目作為相應(yīng)頂點的權(quán)重，用于表示計算負載，最后通過圖剖分(Metis)實現(xiàn)接觸力計算區(qū)域的剖分。

接觸點空間坐標(biāo)的時變性質(zhì)會導(dǎo)致接觸力并行計算時負載的動態(tài)不平衡現(xiàn)象，為了實現(xiàn)負載平衡，需要進行接觸力計算區(qū)域的重剖分。但是從接觸力并行算法可知，其區(qū)域剖分需要耗時的全局操作，為了減少這種區(qū)域重剖分的次數(shù)，在確定內(nèi)力計算區(qū)域和接觸力計算區(qū)域通信關(guān)系時考慮發(fā)送冗余信息，內(nèi)力計算子區(qū)域?qū)⒔佑|力計算子區(qū)域中每個面片相鄰的多個接觸點發(fā)送給接觸力計算子區(qū)域，使得在若干個時間步內(nèi)，接觸力計算子區(qū)域中的面片能在本區(qū)域找到最近的接觸點。

并行程序中消息傳遞采用MPI實現(xiàn)，基于MPI的并行程序可運行在分布式內(nèi)存和共享內(nèi)存并行機上。

3.2 OpenMP實現(xiàn)

OpenMP是科學(xué)計算領(lǐng)域中共享內(nèi)存并行機的主要并行化范式，運行在共享內(nèi)存節(jié)點上的OpenMP程序在內(nèi)存需求、編程效率等各方面都有它的優(yōu)點［8］。OpenMP并行化最簡單的方式是針對程序中的循環(huán)采用指導(dǎo)語句$OMP PARALLEL DO并行化。由于OpenMP指令本身帶來的同步等系統(tǒng)開銷，并行化所有的循環(huán)并不是最好的選擇，可采用基于性能檢測的增量并行模式，其流程見圖2。這種方法通過測試典型算例來確定哪些循環(huán)需要被并行化，能保證所有并行化都能減少總的執(zhí)行時間。

除了采用增量并行化方法外，增加OpenMP并行的粒度，是提高程序并行效率的另外一個途徑。本文采用分塊數(shù)據(jù)結(jié)構(gòu)對單元進行分塊操作。通過將循環(huán)中間計算數(shù)據(jù)的成塊緩存，不僅能充分利用計算機的Cache，而且能增加 OpenMP 并行粒度［9］。

圖2 OpenMP中增量并行模式Fig.2 Incremental parallelization model with OpenMP

以具體程序為例，假設(shè)要對NUM_ELE個單元進行計算，原始程序表示為:

假設(shè)單元塊大小為NBLOCK，假設(shè)單元分為NUM_ELE_BLOCK塊，有 NUM_ELE_BLOCK=NUM_ELE/NBLOCK?；诜謮K數(shù)據(jù)結(jié)構(gòu)的單元塊計算為BLOCK_OPERATION，經(jīng)過分塊和OpenMP并行化后的程序表示為:

采用順序分塊時，即根據(jù)通過對單元順序遍歷，將每NBLOCK個單元為一塊。BLOCK_OPERATION(J)有如下形式:

分塊數(shù)據(jù)結(jié)構(gòu)提供了數(shù)據(jù)遍歷的另一種方法，可以根據(jù)OPERATION中操作的數(shù)據(jù)對象確定分塊方法，以增加每次OPERATION中從內(nèi)存中讀取和存儲的數(shù)據(jù)的局部性，從而提高程序效率。本文程序采用順序分塊，今后將對其他分塊方式進行研究測試。分塊大小NBLOCK是可變的，其大小會影響程序運行效率，具體計算時可根據(jù)實際計算機環(huán)境和具體算例測試進行調(diào)整。

3.3 MPI+OpenMP 實現(xiàn)

基于雙重區(qū)域剖分的接觸并行算法在接觸力計算和內(nèi)力計算子區(qū)域之間的通信是全局通信，隨著處理器數(shù)的增多，全局通信的開銷增加的速度會遠遠大于通過并行計算而減小的計算開銷，通過MPI+OpenMP混合并行模式的計算可有效的減少MPI中的全局通信開銷。

MPI+OpenMP混合并行模式在每個MPI進程計算時增加OpenMP線程。MPI可以在CLUMPs節(jié)點間和節(jié)點內(nèi)部實現(xiàn)信息傳遞，而OpenMP在節(jié)點內(nèi)部利用共享內(nèi)存實現(xiàn)并行計算。OpenMP可使用最大線程數(shù)為節(jié)點內(nèi)部共享內(nèi)存的處理器(核)數(shù)。通過程序輸入?yún)?shù)調(diào)整，程序可以為串行程序，MPI并行程序，OpenMP并行程序，和MPI+OpenMP并行程序，使用何種并行模式取決于可用硬件資源，問題規(guī)模和并行程序特性。圖3給出了接觸并行計算中MPI+OpenMP混合并行計算模式，圖中是在節(jié)點內(nèi)部采用2個OpenMP線程對內(nèi)力計算部分并行實現(xiàn)。

圖3 接觸并行計算中MPI+OpenMP混合并行計算模式Fig.3 MPI+OpenMP hybrid model in contact parallel computing

4 數(shù)值模擬實驗

在某國產(chǎn)大型機上對接觸并行程序進行了數(shù)值模擬測試實驗。實驗平臺計算節(jié)點由8個處理器核共享32G內(nèi)存。模擬錐形桶撞擊平板過程，模型采用8節(jié)點實體單元，模擬自由度為1千萬，計算5000時間步。圖4給出了計算模型。測試時每個計算節(jié)點內(nèi)部可使用1～8個處理器核共享內(nèi)存并行，即OpenMP并行最多可使用8個處理器核，程序中分塊大小取128。圖5給出了從總處理器核數(shù)為P=16～512時，分別采用OMP=1，2，4，8個核共享內(nèi)存時的程序并行效率，若記MPI進程數(shù)為MPIP，有P=OMP*MPIP。

從圖5可知，在總處理器核數(shù)較少時，采用單一模式的MPI程序(OMP=1)的性能好于混合模式的性能。隨著總處理器核數(shù)的增多，采用適當(dāng)?shù)腛penMP共享內(nèi)存并行將提升并行效率?？偺幚砥骱嗽蕉?，為了達到較好的并行效率，每個節(jié)點采用共享內(nèi)存的核數(shù)需要越多，至總處理器核數(shù)512時，每個節(jié)點采用8個核共享內(nèi)存時并行效率最好，此時有64個MPI進程參與并行計算。

為了深入分析產(chǎn)生以上現(xiàn)象的原因，給出總處理器核數(shù)為32，128，512，采用 OpenMP 線程數(shù)為 1，2，4，8時，總計算時間，內(nèi)力計算時間和全局通信時間的趨勢圖，見圖6。從圖中可以看出:① 總處理器核數(shù)一定的情形下，全局通信時間隨著OpenMP線程數(shù)的增加而下降。原因是MPI進程數(shù)隨著OpenMP線程數(shù)的增加而減少。② 總處理器核數(shù)一定的情形下，內(nèi)力計算部分開銷隨著線程數(shù)增加而增加。原因是內(nèi)力并行計算部分通信簡單，基于MPI的并行效率要高于基于OpenMP的并行效率。在總處理器核數(shù)為512時，采用4個線程的計算時間要小于采用MPI并行的時間，其原因可能是由于此時并行粒度較小，OpenMP并行效率在粒度較小時占優(yōu)，但是隨著線程數(shù)的繼續(xù)增加，OpenMP系統(tǒng)開銷增加較多，其運行時間又多于MPI并行程序。③ 隨著總處理器核數(shù)的增多，全局通信在總計算時間中占有比例增加。如32處理器核時內(nèi)力計算時間大于全局通信時間，而128處理器核時兩條線有交叉，至512處理器核時全局通信時間大于內(nèi)力計算時間。④總計算時間隨著OpenMP線程增加的趨勢從左至右依次為單調(diào)增加，先下降再增加，單調(diào)下降。

圖6 不同處理器規(guī)模下的總計算時間(﹡線)，內(nèi)力計算時間(●線)和全局通信時間(▲線)Fig.6 Total time(line with﹡)，Computing time of internal force(line with●)，and communication time(line with▲)with different CPUs

接觸問題性質(zhì)是產(chǎn)生以上現(xiàn)象的主要原因。內(nèi)力計算子區(qū)域和接觸力計算子區(qū)域間需要較多的通信，而且接觸力計算區(qū)域剖分過程中需要全局通信，全局通信在分布式并行機上會嚴(yán)重影響并行效率，隨著參與全局通信的處理器核數(shù)的增多，其影響越大。也就是說接觸力并行計算中采用的處理器核數(shù)越多，全局通信的負面影響將越大。

位于每個節(jié)點內(nèi)部的OpenMP并行模式下的并行程序可擴展性有限。原因之一是OpenMP并行只能針對部分程序并行化，其并行效率受限于串行運行部分;另一個原因是采用較多處理器核時OpenMP并行程序的系統(tǒng)開銷，如隱式同步，Cache一致性等開銷較大，降低了并行程序的效率。在本文算例中，通過增加OpenMP線程得到的內(nèi)力計算部分并行效率在8核時只有50%～55%。在總處理器核數(shù)較少的情形下，OpenMP的并行效率低于MPI并行程序效率。

隨著接觸力并行計算采用的處理器核越多，由接觸力計算部分產(chǎn)生的全局通信的開銷會增加，最終通信開銷的增加會抵消由處理器核數(shù)增加帶來的計算時間好處。而在總處理器核較多時，采用OpenMP并行程序能通過內(nèi)力計算部分的并行計算減少內(nèi)力計算部分的運行時間，但接觸力計算中全局通信量并沒有增加。

5 結(jié)論

針對接觸計算特點，在設(shè)計雙重區(qū)域剖分并行算法的基礎(chǔ)上，采用MPI+OpenMP的混合編程模式實現(xiàn)了接觸問題的并行計算。采用基于分塊的數(shù)據(jù)結(jié)構(gòu)增加了OpenMP的并行粒度，OpenMP并行實現(xiàn)可有效的減少MPI并行程序在處理器核數(shù)較多時接觸力并行計算中的全局通信開銷。數(shù)值實驗表明該并行程序能在數(shù)百處理器上計算上千萬自由度的接觸問題，而且具有較好的并行效率。今后將針對算法特點進行進一步的優(yōu)化。

［1］Attaway S W，Hendrickson B A，Plimpton S J，et al.A parallel contact detection algorithm for transient solid dynamics simulations using PRONTO3D［J］.Computational Mechanics，1998，22(4):143 －159.

［2］Pierce T，Rodrigue G.A parallel two-sided contact algorithm in ALE3D［J］.Computer Methods in Applied Mechanics and Engineering，2005，194(27 －29):3127 －3146.

［3］寇哲君，程建剛，姚振漢.機群環(huán)境下汽車碰撞的并行模擬［J］.清華大學(xué)學(xué)報(自然科學(xué)版)，2003，43(5):666－668.

［4］亓文果，金先龍，張曉云.沖擊－接觸問題有限元仿真的并行計算［J］.振動與沖擊，2006，25(4):68－72.

［5］Cappello F，Etiemble D.MPI versus MPI+OpenMP on IBM SP for the NAS Benchmarks.Supercomputing，ACM/IEEE，2000 Conference［C］.Dallas，04 －10 Nov，2000.

［6］Har J，F(xiàn)ulton R E.A parallel finite element procedure for contact-impact problems［J］.Engineering with Computers，2003，19(2):67－84.

［7］白金澤.LS－DYNA3D理論基礎(chǔ)與實例分析［M］.北京:科學(xué)出版社，2005.

［8］Bader D A，都志輝譯，面向千萬億次計算的算法與應(yīng)用［M］.北京:清華大學(xué)出版社，2008.

［9］肖永浩，黃清南.基于分塊數(shù)據(jù)結(jié)構(gòu)的沖擊問題并行計算［J］.計算機輔助工程，2011，20(1):33 －36.