激光推進數值模擬預處理并行算法研究

摘 要:激光推進數值模擬預處理程序為數值模擬主程序提供網格信息和設置邊界條件，實現并行計算可以縮短大規模網格生成時間，從而提升高分辨率數值模擬的效率。分析并利用公用數據的特點改進原有串行算法，進而實現并行計算。算法測試結果表明，該并行算法有效地縮短了網格生成時間。

關鍵詞:并行計算; 激光推進; 數值模擬; 網格生成

中圖分類號:TP301.6文獻標志碼:A

文章編號:1001-3695(2010)06-2054-03

doi:10.3969/j.issn.1001-3695.2010.06.016

Parallel computing research on preprocessing of numerical simulation for laser propulsion

DU Yang， ZHAO Jun， ZHAO Wen-tao

(School of Computer Science， National University of Defense Technology， Changsha 410073， China)

Abstract:Preprocessing of the numerical simulation for laser propulsion provides grid information and boundary condition of flow field for main program. Parallel computing abbreviates the time of large scale grid generating， then advanced the simulation efficiency in high resolution. The paper made use of the characteristic of public data to improve the algorithm， and rea-lized a parallel algorithm. The algorithm test result indicates that parallel algorithm abbreviates the time of grid generating effectively.

Key words:parallel computing; laser propulsion; numerical simulation; grid generating

激光推進是一種新興的推進技術，由美國學者Kantrowitz[1]于1972年首次提出。與傳統推進技術相比，激光推進技術有比沖大、成本低、污染小以及有效載荷高等許多優勢，因此有著十分廣闊的應用前景，受到各航天大國的廣泛重視。目前對激光推進研究的主要手段是數值模擬，這不僅可以節約研究經費，更能精細地描繪出實驗中難以測量的流場時空演化過程。經過多年的研究發展，激光推進領域涌現出了許多優秀的數值模擬程序[2~7]，極大地促進了激光推進技術的發展。

激光推進數值模擬遵循計算流體力學(CFD)[8]的過程，流場離散為網格，針對定義在網格點的特征量進行反復迭代計算。本文所涉及的激光推進數值模擬預處理程序對Gambit生成的網格數據進行處理，包括網格點的重排和邊界條件的設置等，從而為主程序提供可識別的計算數據，其生成的網格規模直接影響模擬的精細度。隨著激光推進技術研究的深入，對流場演化描繪的精細度要求也逐漸提高，這就需要減小網格間距，從而導致網格點數劇增，需要生成更大規模的網格。然而原有程序的時間復雜度為O(n2)(其中n為網格點規模)，顯然網格點數增大導致的時間代價增加過于迅速。例如，在CPU主頻為3.0 GHz的系統上生成30萬網格點規模耗時766 s，而生成120萬網格點規模耗時達到了14 015 s，這顯然嚴重影響了研究進展，有必要實現并行計算。但是原串行算法程序中頻繁地修改公用數組信息，導致并行計算實現困難。本文對輸入的公用數組信息進行分析，利用信息中各網格點編號的惟一性來替代原程序中頻繁對數組修改而惟一確定后繼點的算法，由于減少比較次數，提高了程序的運行速度;另外無須對公用數組進行修改，使得并行計算易于實現，從而極大地縮短了大規模網格生成時間。

1 串行算法分析

激光推進數值模擬采用Gambit軟件構造網格，但是其所構造的網格不能完全滿足主程序讀取網格信息的輸入順序要求，因此需掃描網格文件，按照要求對網格點進行重排。

輸入的網格文件包含兩部分數據，即網格點信息(網格點的坐標)和網格單元信息(每個單元包含的四個網格點編號)。程序采用三個數組保存這些信息，即X[PN]、Y[PN]、POINT[CN][4]。其中PN表示網格點數(point number);CN表示網格單元數(cell number)。另外使用兩個大小為PN的數組PX、PY保存順序化之后的網格點坐標。網格點的位置有三種情況，即頂點、邊界點和內點，如圖1所示的點A、B、C。

顯然此三類網格點在網格單元信息中的情況也有三種:分別出現1、2、4次，也即三類網格點分別包含在1、2、4個網格單元中。原串行算法步驟設計如下:

a)從文件讀入數據到相應數組。

b)程序采用findminx函數找到第一個點，在二維條件下即為i最小、j最大的點。

c)對于i=1，imax-1:(a)列的第一個網格點保存至臨時變量;(b)查找下一個網格點;(c)按順序將網格點坐標保存至數組PX、PY;(d)該列是否結束，是轉(e)，否轉(b);(e)一列掃描完畢，使用findprp函數找到位于下一列列首的點;(f)if(i

d)對于i=imax，根據首個網格點findLast查找下一個，直到完成整個網格的掃描。

e)輸出順序網格信息(數組PX和PY)到文件。

其中:步驟c)中的(b)查找下一個點的算法為:掃描整個網格單元數組，如果網格單元第j個網格點是前網格點(j=1，2，3，4)，則下一網格點為網格單元的第1(j=4)or第j+1(1，2，3)個網格點，并設置網格點編號取相反數作為訪問標志。查找網格點算法的時間復雜度為O(n)。由于整個流場的每個點都要進行搜索，網格掃描算法的時間復雜度為O(n2)。

算法中需要對point數組進行兩種修改:a)找到一個點后將網格單元信息對應位置取相反數;b)在一列網格點掃描完成后，將該列網格單元的信息置0。以第一列為列，先掃描到包含A的網格單元，則D可惟一確定，按第一種方式修改該網格單元信息，則包含D的網格單元惟一。掃描到本列結束，按第二種方式修改本列網格單元。通過這樣的修改，將網格單元設置為三種狀態，即已掃描、未掃描和正在掃描。已掃描的網格單元無網格點信息，而正在掃描的網格單元，由于前一列網格單元已掃描，則與其公有的網格點信息只存在于本列網格單元之中，這樣處于掃描過程中的每一列都可以看成第一列，由上至下進行本列掃描，只有一個網格單元包含該點，從而保證搜索的正確性。然而這樣的數據修改不僅耗時長，也使并行計算難以實現。因為無論數據如何劃分，每個進程對point數組進行的修改都應該適時地通知其他進程，以維護各進程中數組拷貝的一致性。如果通知不及時，在一個時刻就可能有多個網格單元包含目標網格點，從而導致錯誤的產生，由此可見各數據塊有很強的數據依賴關系。如果使用集合通信對point數組進行更新，由于通信非常頻繁而導致巨大的通信開銷讓并行效率異常低下，甚至不如串行算法，該算法不適合進行并行計算，需要加以改進。

2 串行算法改進

實現并行計算的關鍵在于保證公用數據point數組的穩定性，即不能進行寫操作，這樣才可以使各個進程并行掃描point數組尋找屬于本區域的網格點。

由于每一個點在point中最多出現四次，且在四個網格單元中的編號不同，網格點在各網格單元中按逆時針順序排列的編號有兩種情況:a)網格點的編號在各網格單元中完全一致，如圖2所示。橫向掃描根據如下條件確定下一個網格單元:前一個網格單元右上的網格點編號為k(圖2中第2個網格點)是下一個網格單元左上的網格點，編號為k%4+1;縱向掃描則根據如下條件確定:前一個網格單元左下網格點編號為k是下一個網格單元左上網格點，編號為j-1(k=1是為4)。b)編號順序有兩種，有公共邊的不同，只有公共點的網格單元相同，如圖3 所示。這種情況的掃描與情況a)正好相反。由于網格點的編號在各網格單元中是惟一的，從而根據編號搜索到的網格單元也是惟一的。

在新算法中，新增一個seq數組記錄網格點順序，如seq[i]=j表示原有網格中第j點在新的網格中的排序為i，可以設定如下掃描算法:

a)從文件讀入數據;

b)搜索第一個點，并確定網格單元的網格點編號方式;

c)根據兩種不同編號方式選擇不同的橫向掃描算法，得到第一行的各網格單元;

d)由得到的各列首個網格單元依次進行縱向掃描，找到一個網格點順序號賦予seq當前指標指向位置，并將seq的指標加1;

e)最后一列網格點的掃描與該編號方式下的橫向掃描類似，按照上述相同的方式修改seq數組;

f)根據seq數組中網格點的順序，將x，y中的數據保存至文件， 也即save x[seq[i]]，y[seq[i]] to file。

算法比較:原算法在整個網格掃描中，需對各網格單元編號為1、2、3、4的點依次進行判斷，改進后的算法只需對網格單元的一個編號進行判斷，減少了冗余判斷次數。另外，改進后的算法無須對point數組進行修改即可準確搜索目標網格點，減少了修改數組的耗時。在第一行掃描完畢后，各列之間無數據依賴關系，可獨立執行縱向掃描，因此易實現并行計算。

原程序使用五個數組，即四個一維雙精度浮點數組X[PN]、Y[PN]、PX[PN]、PY[PN]，一個二維整形數組point[CN][4]。由于雙精度占8 Byte，整形數占2 Byte，原算法占用的內存可以近似表示為(2+8×4)×PN=34×PN。改進的算法中需增加一個整形數組seq[PN]，刪除了兩個浮點數組PX、PY，內存占用可以表示為(2+2+8×2)×PN= 20×PN。由此可見，改進后的算法優化內存的使用。

3 并行算法設計

3.1 數據劃分

本文采用輸出數據劃分方法，對保存掃描結果的seq數組按列進行一維劃分，使用靜態任務映射方案[9]將任務分配到各進程。設進程數為n，列數為j，則將j/n分配到各個進程而余下的平均分配到前j%n個進程，這樣各個進程任務大小相差最多不超過一列，因此在行數i較小的情況下負載均衡性較好。

各進程所占數據的長度l:

if(p≤j%n) l=i×(j/n+1);

else l=i×(j/n);

其中:j為列數;p為進程號;n為總進程數。

各進程都須對整個網格掃描，因此point數組在每個進程上都應有拷貝，減少通信開銷。各處理單元分別掃描point數組，搜索屬于本區域的網格點，保存至seq數組對應的位置。在整個掃描過程中無須與其他進程交互，由此可以預計并行程序有良好的可擴展性，并可以達到較高的加速比。

3.2 數據通信

程序運行過程中，只需進行兩次通信:a)第一行橫向掃描完成后，主進程將掃描結果廣播至各進程;b)各進程在本區域掃描結束后將結果發送至主進程。因此通信開銷對系統影響較小。具體設計如下:

a)由于point數組需要在每個進程中都有拷貝，同時讀取文件可能造成沖突，本文設計為由主進程讀取之后進行廣播，主進程在進行列首網格單元的掃描之后，將結果廣播。通信偽代碼可以表示如下:

if(myid==0)

readfile to point，x，y;

first row scan ;

broadcast array point;

broadcast scan result first;

b)主進程在各從進程掃描完成后，收集從進程局部seq數組，獲得整個網格的網格點順序，至此并行計算結束。偽代碼表示如下:

if(myid==0)

for(i=1;i≤myid;i++)

receive seq from process i;

else send seq to process 0;

4 算法測試結果與分析

測試環境:Linux集群系統，每個節點兩個Intel 4核Xeon E5450微處理器，主頻3.0 GHz，CPU cache容量為6 MB，內存總容量16 GB，節點之間使用Infiniband互連。

測試方案如下:

a)串行程序測試。在網格規模為600×500 和1200×1000兩種情況下進行測試，記錄運行時間如表1所示，并對各種輸出結果與原有串行程序的輸出進行比較，結果完全一致。由表1可以看出，修改后的程序運行時間明顯地縮短了。

表1 串行程序優化效果測試

網格規模原串行程序/s優化程序/s優化加速比

600×5007667210.64

1 200×1 00014 0154 4633.14

b) 并行程序測試。在同樣的兩種網格規模下，對并行程序進行測試，并以優化后的串行程序作為標準計算加速比和并行效率，結果如表2所示。

表2 并行程序加速效果測試

進程數

網絡規模

600×500

運行時間/s加速比并行效率/%

1 200×1 000

運行時間/s加速比并行效率/%

17211446311

2381.89594.7522761.9698

4233.1378.2511603.8596.25

8145.1464.256357.0287.75

168956.2532713.6585.31

3261237.526017.1753.7

將加速比和并行效率的變化更直觀地表示如圖4、5所示。

并行程序性能受三方面因素影響，即通信開銷、負載均衡和并行代碼在程序中的比重。程序中的通信僅需兩次，因此開銷較小;與串行代碼開銷一樣比較固定，受進程數和網格規模變化的影響較小。因此加大網格規模和減少進程數，使進程的計算量增加，可以提升并行效率。導致負載不均衡的因素可能來自兩方面:首先，各網格點在數組中的位置不同，這就導致搜索各點時間的差別，而且網格點具有一定的局部性，即空間相鄰的點在數組中的位置也很可能相鄰。因此，即使兩個進程擁有的網格點數量相同，搜索時間的差別也可能很大。其次，采用數據劃分方法雖然使進程之間任務量只相差一列，但是隨著行數增大，也可能產生不小影響。

從圖4中可以看出，在進程數較少的情況下，并行程序可以達到近似線性加速比。另外也可以看出網格規模越大，在進程數相等的情況下，加速效果越明顯。這都是因為單個進程計算量大，而導致通信開銷和串行開銷比重減小。

圖5表明相同網格規模情況下，進程數增多導致并行效率的降低。原因之一是進程數增多使進程的計算量減少，通信開銷與串行開銷比重上升;另外，由于進程任務量減少，進程之間計算量差距導致的影響就會更加明顯;運行進程數相等的情況下，網格規模增大會提升并行效率，因為規模增大使得進程計算量增加，不僅提高計算通信比，同時也削弱了負載不均衡的影響，還減少了串行執行部分所占比重。由此可以預計，網格規模進一步增大時可以獲得更好的并行性能。

5 結束語

本文通過對激光推進預處理串行算法進行深入分析，找出其不易實現并行計算的原因是由于頻繁地對公共數組進行寫操作;然后利用存在于輸入數據中各網格點在網格單元中編號的特點進行算法改進，減少了冗余的比較操作，優化了內存使用，更重要的是避免了對公共數組的寫操作，從而縮短了網格生成時間。實現的并行算法使得大規模網格的構造時間進一步縮減，從而對高分辨率下激光推進數值模擬研究提供了很好的支持作用。然而如第4章所分析，程序存在一些影響負載平衡的因素，因此在后續工作中，考慮使用動態負載均衡的方法并實行二維劃分，以提高并行程序的加速比和并行效率。

參考文獻:

[1]KANTROWITZ A. Propulsion to orbit ground-based lasers[J]. Astronautics Aeronautics，1972，10(5):74-75.

[2]THOMPSON S L， LAUSON H S. Improvements in the chart-D radiation-hydrodynamic code III:revised analytical equation of state[R].1972.

[3]CHRISTIANSEN J P， WINSOR N K. Castor 2:a 2-D laser target code[J]. Computer Physics Communications，1979，17(4):397.

[4]GONG Ping， TANG Zhi-ping. Numerical simulation for laser propulsion of air breathing mode considering moving boundaries and multi-pulses[C]//Proc of the 4th International Symposium on Beamed Energy Propulsion.2006:87-94.

[5]李倩，洪延姬，曹正蕊.吸氣式激光推進推力產生機理的數值模擬[J].爆炸與沖擊， 2006，26(6):550-555.

[6]許仁萍，唐志平.大氣模式下多脈沖激光推進的數值模擬[J].強激光與粒子束， 2007，19(3):369-372.

[7]童慧峰，唐志平，張凌.燒蝕模式激光推進的數值模擬[J].爆炸與沖擊， 2007，27(2):165-170.

[8]王承堯，王正華，楊曉輝.計算流體力學及其并行算法[M].長沙:國防科技大學出版社，2000:2-10.

[9]GRAMA A， GUPTA A， KARYPIS G， et al. Introduction to parallel computing [M].2版.張武，毛國勇，程海英，等譯.北京:機械工業出版社，2005:85-91.