光滑粒子流體動力學的一種并行數值計算方案

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

（國防科學技術大學 理學院 工程物理研究所，長沙 410073）

0 引言

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

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

1 SPH方法并行計算方案

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

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

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

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

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

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

為了簡便，下面以np=5為例展開討論。此時，所有鄰域粒子對被分為如下 15類：(1, 1)、(1,2)、……、(1, 5)，(2, 2)、……、(2, 5)，(3, 3)、……、(3, 5)、……、(5, 5)。每個進程承擔3類，各個子任務在每個進程上面的分配如表1所示。

表1 在5個進程上的任務劃分Table 1 Task partition of five processors

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

在搜索的同時，每個進程上如果發現粒子i、j是近鄰粒子，那么立刻計算有關i、j粒子的核函數及其導數，并且在array2(i)和array2(j)上面各加上一項。這樣，進程1就計算了第1、2、3組粒子隨體導數的部分信息，進程2計算了第2、3、4組粒子隨體導數的部分信息，……，進程5計算了第5、1、2組粒子隨體導數的部分信息，如此處理就不需要一個超大數組來存儲所有粒子對信息。

1.2 每一次循環的計算步驟

在循環之前先建立2m+1個通信子，每個通信子包含的進程對應于表2中每行的所有進程，如通信子1包含進程1、m+2、m+3、……、2m+1，通信子 2包含進程 2、m+3、m+4、……、2m+1、1等等。計算中每一次循環包括以下四個步驟。

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

表2 第一步中的信息流動方向Table 2 Information flow in step one

可以看出，這一步中只需要每個進程在循環之前創建的通信子中調用若干次廣播函數。

第二步：粒子搜索與隨體導數更新。每個進程用自己維護的那部分粒子造一棵樹，然后找到所有分配給自己的粒子對，并更新相應粒子的隨體導數。對應于表1，進程1用第1組粒子造一棵樹，找出粒子對(1, 1)、(1, 2)、(1, 3)并更新第1、2、3組粒子的隨體導數；……；進程5用第5組粒子造一棵樹，找出粒子對(5, 5)、(5, 1)、(5, 2)并更新第5、1、2組粒子的隨體導數。每一次搜索完成之后要編一個遞歸子程序銷毀之前造的樹，否則計算一段時間之后會因為虛擬內存不足而導致計算停止。這一步是 SPH并行計算的主要步驟，占了整個計算時間的絕大部分。

第三步：交換粒子隨體導數信息。對應于表1，進程4和5有部分關于第1組粒子隨體導數的信息需要傳遞給進程1，進程5和1有部分關于第2組粒子隨體導數的信息需要傳遞給進程2等等。這一步的信息流動方向與第一步剛好相反，所以只需要每個進程在循環之前創建的通信子中調用若干次求和歸約函數。這一步結束后，每個進程就在array2中相應位置保存了自己維護的那部分粒子的隨體導數。

第四步：每個進程利用粒子隨體導數信息和時間步長更新自己維護的那部分粒子的信息。

為了測試以上4步的時間消耗情況，我們設置了1 825 600個粒子，分別申請不同數量的進程，得到的結果如圖2所示。

圖2 不同進程數的時間消耗對比Fig. 2 Running time with different numbers of processors

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

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

2 計算實例

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

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

由圖3可以看到，計算結果和試驗結果基本吻合。試驗測量的碎片云長度為55.4 mm，而計算結果為55.1 mm。

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

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

由圖4可以看到，碎片云的計算結果和試驗結果基本吻合。

3 結束語

SPH方法計算三維超高速碰撞很有優勢，但運算量很大。針對這個情況，我們設計了一種并行計算方案并且實現了數值模擬。這種方案簡單直接，易于編程，整個并行三維SPH程序只有1 000多行，而且和并行相關的只有主程序中的200多行，從幾萬到幾百萬個粒子規模以內能夠提供大約10的加速比，即并行計算時間為串行計算的1/10，顯著減少了SPH方法的計算時間。

（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]李寶寶. 超高速碰撞下相變效應的數值模擬研究[D].長沙: 國防科學技術大學碩士學位論文, 2010

[4]陳國良. 并行計算——結構、算法、編程[M]. 北京: 高等教育出版社, 2003

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

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

[7]徐金中. 空間碎片超高速碰撞特性及其防護結構優化設計的SPH研究[D]. 長沙: 國防科學技術大學博士學位論文, 2008

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