周 浩，湯文輝，冉憲文，陳 華

（國(guó)防科學(xué)技術(shù)大學(xué) 理學(xué)院 工程物理研究所，長(zhǎng)沙 410073）

0 引言

隨著航天科技的迅速發(fā)展以及航天發(fā)射任務(wù)的急劇增加，大量的空間碎片散布在近地軌道上，形成了復(fù)雜的空間碎片環(huán)境，對(duì)航天器的安全構(gòu)成了極大的威脅。如何確保航天器的飛行安全，成為當(dāng)前航天器設(shè)計(jì)中一個(gè)十分重要的問(wèn)題。對(duì)超高速碰撞問(wèn)題進(jìn)行理論研究非常困難，而現(xiàn)有的試驗(yàn)技術(shù)又很難達(dá)到超高速碰撞所需要的速度條件，因此數(shù)值模擬是一種重要的研究手段。有限元方法作為當(dāng)前計(jì)算力學(xué)領(lǐng)域主要的數(shù)值計(jì)算方法，經(jīng)過(guò)多年的發(fā)展，日漸成熟和完善，被廣泛應(yīng)用于解決工程問(wèn)題。但是在面對(duì)超高速碰撞等涉及到大變形的問(wèn)題時(shí)，應(yīng)用有限元方法經(jīng)常會(huì)產(chǎn)生網(wǎng)格畸變，嚴(yán)重影響計(jì)算精度，甚至導(dǎo)致計(jì)算終止。

光滑粒子流體動(dòng)力學(xué)（smoothed particle hydrodynamics，簡(jiǎn)稱SPH）方法[1-3]是近年來(lái)比較熱門的一種無(wú)網(wǎng)格流體動(dòng)力學(xué)算法，它既保留了拉氏計(jì)算描述物質(zhì)界面準(zhǔn)確的優(yōu)勢(shì)，又具備歐拉方法的長(zhǎng)處，而且邏輯簡(jiǎn)單，統(tǒng)一地處理不同維數(shù)的問(wèn)題，能夠避免有限元方法中的網(wǎng)格纏繞和扭曲等問(wèn)題，因而特別適合計(jì)算多物質(zhì)的大變形問(wèn)題。在SPH方法中，物質(zhì)被離散為一系列“粒子”，每個(gè)粒子攜帶如密度、壓力、速度、內(nèi)能等屬性，通過(guò)質(zhì)量、動(dòng)量、能量三個(gè)守恒方程，利用附近粒子相應(yīng)物理量插值得到各個(gè)粒子物理量的隨體導(dǎo)數(shù)，使拉氏流體力學(xué)偏微分方程組變成差分方程組，從而易于數(shù)值求解。SPH中的一個(gè)主要問(wèn)題是計(jì)算量比較大，耗時(shí)較長(zhǎng)，所以開(kāi)展并行計(jì)算研究十分必要。

1 SPH方法并行計(jì)算方案

并行編程比較復(fù)雜，面臨著執(zhí)行結(jié)果的不確定、如何在各個(gè)進(jìn)程間合理地劃分任務(wù)、如何減少通信，以及對(duì)于不同進(jìn)程數(shù)（即不同數(shù)量的核）具有伸縮性等問(wèn)題[4-6]。本文從實(shí)際需要出發(fā)，采用基于MPI（Message Passing Interface）的信息傳遞并行程序設(shè)計(jì)平臺(tái)和粒子分割算法。每個(gè)進(jìn)程中都有array1和array2兩個(gè)數(shù)組，分別用來(lái)存儲(chǔ)所有粒子信息和所有粒子隨體導(dǎo)數(shù)信息。在信息處理中，只需維護(hù)每個(gè)進(jìn)程中互不重疊的一部分信息。雖然這樣處理每個(gè)進(jìn)程會(huì)降低程序編寫復(fù)雜度（特別是在遞歸搜索粒子對(duì)以及計(jì)算各個(gè)粒子隨體導(dǎo)數(shù)時(shí)），但仍存在一些冗余信息。

1.1 近鄰粒子搜索的并行方案

SPH計(jì)算最耗時(shí)的部分是近鄰粒子搜索，而并行的 SPH方法對(duì)近鄰粒子采用了并行搜索算法，可有效減少計(jì)算時(shí)間。目前主要的并行搜索算法有全配對(duì)搜索、鏈表搜索和樹(shù)搜索3種，其中：全配對(duì)搜索算法最簡(jiǎn)單直接，但是計(jì)算量為 N2，只適合于極少數(shù)粒子的情況；鏈表搜索算法效率比較高，但是對(duì)于可變光滑長(zhǎng)度缺乏自適應(yīng)能力；而樹(shù)搜索算法效率相對(duì)較高，且對(duì)于可變光滑長(zhǎng)度具有良好的自適應(yīng)能力。因此，本文采用了樹(shù)搜索算法。

近鄰粒子樹(shù)搜索算法是 Hernquist和 Katz于1989年在TreeSPH方法研究中提出來(lái)的。它首先通過(guò)遞歸的方法構(gòu)造一棵樹(shù)，樹(shù)的每個(gè)節(jié)點(diǎn)代表一塊空間。樹(shù)的根節(jié)點(diǎn)代表整個(gè)計(jì)算空間，子節(jié)點(diǎn)代表父節(jié)點(diǎn)空間的子空間。如果某個(gè)子節(jié)點(diǎn)包含的粒子數(shù)大于1，則繼續(xù)劃分此節(jié)點(diǎn)，直到所有子節(jié)點(diǎn)只包含一個(gè)粒子。對(duì)于任意粒子i，以兩倍的光滑長(zhǎng)度為邊長(zhǎng)構(gòu)造一個(gè)立方體，使i粒子位于立方體的中心。檢驗(yàn)該立方體與樹(shù)的根節(jié)點(diǎn)所代表的空間是否有重合，如果有，則遍歷此樹(shù)節(jié)點(diǎn)的所有子節(jié)點(diǎn)，直到當(dāng)前樹(shù)節(jié)點(diǎn)所代表的空間處僅有一個(gè)粒子為止。然后判斷該粒子是否在粒子i的影響域內(nèi)，如果是，則該粒子為i粒子的近鄰粒子；如果不是，則停止在此節(jié)點(diǎn)及其所有子節(jié)點(diǎn)的搜索。對(duì)于二維情況，樹(shù)搜索算法如圖1所示[7]。

圖1 樹(shù)搜索算法示意圖Fig. 1 Tree search algorithm in two-dimensional space

假設(shè)參與計(jì)算的進(jìn)程總數(shù)為np且為奇數(shù)，于是令np=2m+1。首先把所有粒子大概平均分為np組，每個(gè)進(jìn)程維護(hù)一組粒子。設(shè)第i組和第j組之間的鄰域粒子對(duì)表示為(i, j)，其中 i=1, 2, ……, np，j=i,i+1, ……, np，于是所有的鄰域粒子對(duì)被分成了np(np+1)/2類，即將一個(gè)大任務(wù)分成為np(np+1)/2個(gè)小任務(wù)，而每個(gè)進(jìn)程將承擔(dān)(np+1)/2=m+1個(gè)小任務(wù)。

為了簡(jiǎn)便，下面以np=5為例展開(kāi)討論。此時(shí)，所有鄰域粒子對(duì)被分為如下 15類：(1, 1)、(1,2)、……、(1, 5)，(2, 2)、……、(2, 5)，(3, 3)、……、(3, 5)、……、(5, 5)。每個(gè)進(jìn)程承擔(dān)3類，各個(gè)子任務(wù)在每個(gè)進(jìn)程上面的分配如表1所示。

表1 在5個(gè)進(jìn)程上的任務(wù)劃分Table 1 Task partition of five processors

在進(jìn)程1上用第1組粒子造一棵樹(shù)，則第1、2、3組粒子在這棵樹(shù)上搜索，找到近鄰粒子對(duì)(1,1)、(1, 2)、(1, 3)；在進(jìn)程2上用第2組粒子造一棵樹(shù)，則第2、3、4組粒子在這棵樹(shù)上搜索，找到近鄰粒子對(duì)(2, 2)、(2, 3)、(2, 4)；……；在進(jìn)程5上用第5組粒子造一棵樹(shù)，第5、1、2組粒子在這棵樹(shù)上搜索，找到近鄰粒子對(duì)(5, 5)、(1, 5)、(2, 5)。由此可以看到：在每個(gè)進(jìn)程上用某組粒子構(gòu)造一棵樹(shù)時(shí)需要3組粒子的信息，而在一步搜索中需在這棵樹(shù)上搜索3組粒子。

在搜索的同時(shí)，每個(gè)進(jìn)程上如果發(fā)現(xiàn)粒子i、j是近鄰粒子，那么立刻計(jì)算有關(guān)i、j粒子的核函數(shù)及其導(dǎo)數(shù)，并且在array2(i)和array2(j)上面各加上一項(xiàng)。這樣，進(jìn)程1就計(jì)算了第1、2、3組粒子隨體導(dǎo)數(shù)的部分信息，進(jìn)程2計(jì)算了第2、3、4組粒子隨體導(dǎo)數(shù)的部分信息，……，進(jìn)程5計(jì)算了第5、1、2組粒子隨體導(dǎo)數(shù)的部分信息，如此處理就不需要一個(gè)超大數(shù)組來(lái)存儲(chǔ)所有粒子對(duì)信息。

1.2 每一次循環(huán)的計(jì)算步驟

在循環(huán)之前先建立2m+1個(gè)通信子，每個(gè)通信子包含的進(jìn)程對(duì)應(yīng)于表2中每行的所有進(jìn)程，如通信子1包含進(jìn)程1、m+2、m+3、……、2m+1，通信子 2包含進(jìn)程 2、m+3、m+4、……、2m+1、1等等。計(jì)算中每一次循環(huán)包括以下四個(gè)步驟。

第一步：粒子信息交換。對(duì)于np=5的情形（如表1所示），進(jìn)程1上面存儲(chǔ)了第1組信息，而進(jìn)程4、5為了維護(hù)第4、5組粒子信息，需要知道第1組粒子信息，所以進(jìn)程1需要將第1組粒子信息傳遞給進(jìn)程4和5。同理，進(jìn)程5需要將第5組粒子信息傳遞給進(jìn)程3和4。這一步信息流動(dòng)方向如表2所示。

表2 第一步中的信息流動(dòng)方向Table 2 Information flow in step one

可以看出，這一步中只需要每個(gè)進(jìn)程在循環(huán)之前創(chuàng)建的通信子中調(diào)用若干次廣播函數(shù)。

第二步：粒子搜索與隨體導(dǎo)數(shù)更新。每個(gè)進(jìn)程用自己維護(hù)的那部分粒子造一棵樹(shù)，然后找到所有分配給自己的粒子對(duì)，并更新相應(yīng)粒子的隨體導(dǎo)數(shù)。對(duì)應(yīng)于表1，進(jìn)程1用第1組粒子造一棵樹(shù)，找出粒子對(duì)(1, 1)、(1, 2)、(1, 3)并更新第1、2、3組粒子的隨體導(dǎo)數(shù)；……；進(jìn)程5用第5組粒子造一棵樹(shù)，找出粒子對(duì)(5, 5)、(5, 1)、(5, 2)并更新第5、1、2組粒子的隨體導(dǎo)數(shù)。每一次搜索完成之后要編一個(gè)遞歸子程序銷毀之前造的樹(shù)，否則計(jì)算一段時(shí)間之后會(huì)因?yàn)樘摂M內(nèi)存不足而導(dǎo)致計(jì)算停止。這一步是 SPH并行計(jì)算的主要步驟，占了整個(gè)計(jì)算時(shí)間的絕大部分。

第三步：交換粒子隨體導(dǎo)數(shù)信息。對(duì)應(yīng)于表1，進(jìn)程4和5有部分關(guān)于第1組粒子隨體導(dǎo)數(shù)的信息需要傳遞給進(jìn)程1，進(jìn)程5和1有部分關(guān)于第2組粒子隨體導(dǎo)數(shù)的信息需要傳遞給進(jìn)程2等等。這一步的信息流動(dòng)方向與第一步剛好相反，所以只需要每個(gè)進(jìn)程在循環(huán)之前創(chuàng)建的通信子中調(diào)用若干次求和歸約函數(shù)。這一步結(jié)束后，每個(gè)進(jìn)程就在array2中相應(yīng)位置保存了自己維護(hù)的那部分粒子的隨體導(dǎo)數(shù)。

第四步：每個(gè)進(jìn)程利用粒子隨體導(dǎo)數(shù)信息和時(shí)間步長(zhǎng)更新自己維護(hù)的那部分粒子的信息。

為了測(cè)試以上4步的時(shí)間消耗情況，我們?cè)O(shè)置了1 825 600個(gè)粒子，分別申請(qǐng)不同數(shù)量的進(jìn)程，得到的結(jié)果如圖2所示。

圖2 不同進(jìn)程數(shù)的時(shí)間消耗對(duì)比Fig. 2 Running time with different numbers of processors

圖2(a)中，第1、3步屬于通信，運(yùn)算時(shí)間隨著進(jìn)程數(shù)的增加而增加；而第2、4步屬于計(jì)算，運(yùn)算時(shí)間隨著進(jìn)程數(shù)的增加而減少。

圖2(b)表示對(duì)于大約182萬(wàn)個(gè)粒子，在23個(gè)進(jìn)程內(nèi)，計(jì)算時(shí)間隨著進(jìn)程數(shù)的增加而減少。在用23個(gè)核并行計(jì)算時(shí)，一次循環(huán)中4步計(jì)算加起來(lái)總共需要4.34 s，而在“銀河”機(jī)上串行計(jì)算需要41.5 s，所以加速比為9.56，并行效率為41.57%。但是如果進(jìn)程數(shù)繼續(xù)增加到一定數(shù)量，則粒度將繼續(xù)減小，通信時(shí)間的增加將超過(guò)每個(gè)進(jìn)程上因任務(wù)減少而減少的計(jì)算時(shí)間，因而計(jì)算時(shí)間將增加。于是對(duì)于任意固定的粒子數(shù)，總是存在一個(gè)最優(yōu)的進(jìn)程數(shù)。由于權(quán)限之故，本文只考慮了進(jìn)程數(shù)在 24個(gè)以內(nèi)的情況。

2 計(jì)算實(shí)例

實(shí)例為圓柱形鋅彈丸正面碰撞鋅板。彈丸直徑為3.98 mm，長(zhǎng)14.15 mm；鋅板厚0.965 mm；粒子之間的初始距離為0.1 mm，彈丸設(shè)置176 080個(gè)粒子，靶板設(shè)置1 243 770個(gè)粒子，總粒子數(shù)為1 419 850；彈丸速度為4.97 km/s。在“銀河”5萬(wàn)億次計(jì)算機(jī)上申請(qǐng)23個(gè)核并行計(jì)算。采取變時(shí)間步長(zhǎng)，模擬10.5 μs需要3 160步，耗時(shí)3.25 h；而串行計(jì)算則需要36.7 h，因此加速比大約為11.3，并行效率為49.1%。試驗(yàn)結(jié)果[8]與計(jì)算結(jié)果對(duì)比如圖3所示。

圖3 鋅桿彈丸撞擊鋅板試驗(yàn)與模擬對(duì)比Fig. 3 Comparison between the experiment and calculation for the impact of zinc rod on zinc slab (10.5 μs after impact )

由圖3可以看到，計(jì)算結(jié)果和試驗(yàn)結(jié)果基本吻合。試驗(yàn)測(cè)量的碎片云長(zhǎng)度為55.4 mm，而計(jì)算結(jié)果為55.1 mm。

第二個(gè)例子為銅盤正面碰撞鋁板。銅盤彈丸重3 g，其速度為5.55 km/s；直徑11.18 mm，厚3.45 mm，設(shè)置343 105個(gè)粒子。鋁板厚為2.87 mm，設(shè)置2 562 005個(gè)粒子。模擬6.4 μs需要1 737步，耗時(shí)3.95 h；而串行計(jì)算則需47.1 h，因此加速比約為11.9，并行效率為51.8%。試驗(yàn)結(jié)果[8]與計(jì)算結(jié)果對(duì)比如圖4所示。

圖4 銅盤撞擊鋁板試驗(yàn)與模擬對(duì)比Fig. 4 Comparison between the test and calculation for copper disk impacting on aluminum plane(6.4 μs after impact)

由圖4可以看到，碎片云的計(jì)算結(jié)果和試驗(yàn)結(jié)果基本吻合。

3 結(jié)束語(yǔ)

SPH方法計(jì)算三維超高速碰撞很有優(yōu)勢(shì)，但運(yùn)算量很大。針對(duì)這個(gè)情況，我們?cè)O(shè)計(jì)了一種并行計(jì)算方案并且實(shí)現(xiàn)了數(shù)值模擬。這種方案簡(jiǎn)單直接，易于編程，整個(gè)并行三維SPH程序只有1 000多行，而且和并行相關(guān)的只有主程序中的200多行，從幾萬(wàn)到幾百萬(wàn)個(gè)粒子規(guī)模以內(nèi)能夠提供大約10的加速比，即并行計(jì)算時(shí)間為串行計(jì)算的1/10，顯著減少了SPH方法的計(jì)算時(shí)間。

（References）

[1]Liu G R, Liu M B. Smoothed particle hydrodynamics[M].Singapore: World Scientific, 2003

[2]Monaghan J J. Particle methods for hydrodynamics[J].Computer Physics Report, 1985(3): 71-124

[3]李寶寶. 超高速碰撞下相變效應(yīng)的數(shù)值模擬研究[D].長(zhǎng)沙: 國(guó)防科學(xué)技術(shù)大學(xué)碩士學(xué)位論文, 2010

[4]陳國(guó)良. 并行計(jì)算——結(jié)構(gòu)、算法、編程[M]. 北京: 高等教育出版社, 2003

[5]都志輝. 高性能計(jì)算之并行編程技術(shù)——MPI并行程序設(shè)計(jì)[M]. 北京: 清華大學(xué)出版社, 2001

[6]張武生, 薛魏, 李建江, 等. MPI并行程序設(shè)計(jì)實(shí)例教程[M]. 北京: 清華大學(xué)出版社, 2009

[7]徐金中. 空間碎片超高速碰撞特性及其防護(hù)結(jié)構(gòu)優(yōu)化設(shè)計(jì)的SPH研究[D]. 長(zhǎng)沙: 國(guó)防科學(xué)技術(shù)大學(xué)博士學(xué)位論文, 2008

[8]Mullin S A, Littlefield D L, Anderson C E, et al. An examination of velocity scaling[J]. Ibid, 1995, 17:571-581