蒙特卡洛模擬的并行實現及并行效率研究

柳 建 孫勝軍 毛國平 石秀安

1(成都理工大學工程技術學院 四川 樂山 614007) 2(中科華核電技術研究院 廣東 深圳 518026)

0 引 言

為了縮短計算周期，同時保證結果的可信度，有必要利用日益成熟的并行計算技術[2]。通過分析MC算法流程，選用OpenMP多線程并行技術對原串行程序進行了并行化[3-4]，并在大型SMP服務器Oracle M9000進行了實驗。結果顯示，并行計算結果除與串行結果完全一致，還充分發揮了M9000的強大性能，獲得滿意的加速性能，極大縮短了計算周期。

1 MC模擬算法流程及其并行化

MC模擬的算法流程如圖1所示(以中子為例)。通過跟蹤源粒子庫中的所有中子在系統中的運動歷史，用統計方法得出中子通量密度、劑量率、沉積能或各種反應率等所需要的結果。

圖1 中子跟蹤流程

具體而言，MC計算程序遵循如圖2所示的流程，實現對每個中子的歷史跟蹤過程，其中值得特別注意的有兩個主要問題：次級粒子的處理、判斷中子歷史是否結束。當中子發生(n,2n),(n,3n),(n,f)等反應時，將有次級中子產生，例如(n,2n)反應會產生2個次級中子，連同原來的中子，共跟蹤3個中子。但每次計算只能跟蹤一個中子，因此需要把其余兩個中子先存起來。同樣，在(n,γ)反應和(n,f)反應中產生的次級γ-光子，也需要先存起來。當一個中子歷史結束后，再跟蹤源中的其他中子，直到全部中子跟蹤完，完成一個循環步。然后再同理跟蹤次級粒子。

圖2 MC模擬代碼流程框圖

通過以上的程序流程可知，MC算法中的中子的運動歷史之間沒有數據相關性，因而非常適合并行計算[5]以縮短計算所需桌面時間。

當前的高效并行技術主要有OpenMP、MPI和PVM[5-8]。其中OpenMP是支持共享存儲并行方式的庫，適合在多核CPU上進行并行，采用共享存儲系統，編程和執行效率都很高，因而得到了廣泛應用。而在蒙特卡羅模擬中，在算法上，源粒子經抽樣后，它們在計算中是數據無關的，因而可以采用任意并行方式。由于改進的程序將主要運行在共享存儲的SMP結構硬件平臺上，同時為減小線程間的通信數據量、獲得好的負載平衡，此處采用來了OpenMP并行庫為基礎的數據并行方式[3-4]。

對圖2中的串行代碼進行多線程并行化時，為提高并行效率，一般需要進行粗粒度的并行，以減少線程喚醒、掛起和等待的時間。為此將并行區域設置為流程圖中最大的循環區域，即以每個源粒子的歷史過程作為劃分線程粒度的標準。程序中通過將該區域使用OpenMP的parallel關鍵字來指定為并行區域。為了各線程總體計算量負載平衡，通常每個線程每次抽取1個源粒子。由于每個線程在累計已抽樣粒子數、調用截面數據、匯總計數結果、累加運行錯誤等情形下會使用到公共數據，這將導致線程之間會出現數據競爭，為避免這個問題，在這些地方通過critical關鍵字或atomic關鍵字加以防止。

最終的并行程序結構中具體采用了如圖3的fork & join并行方式[8]，思路和數據流程簡潔清晰。

圖3 采用OpenMP方式并行的程序結構圖

2 測試算例及并行結果

大型SMP服務器Oracle M9000擁有64個獨立處理器，每個處理器有4個獨立核，能夠同時處理512個線程，非常適合運行多線程并行程序。將改進的并行程序在M9000上運行，對一個在工程中遇到的輻射防護問題進行了計算，其中源中子數給定為24 373 075個，并行線程數從1開始，直到256個線程結束。該問題的幾何模型如圖4所示。

圖4 反應堆外殼幾何模型

為定量考察并行代碼的計算結果精度、加速性能，對幾個主要闡述：不同線程數時的總碰撞次數、線程數n、桌面耗時Ct(分鐘)和CPU使用率Usg進行記錄和分析。結果顯示，在各線程數下的總碰撞次數和串行程序的計數結果完全一致，并行結果的準確性得到驗證，下面僅對并行后的加上效率進行討論。

測試結果如表1所示。其中，Ct是通過直接記錄程序運行完成時間得到，Usg是通過使用M9000上的SOLARIS操作系統的“prstat”命令統計得到。通過表1，對時耗和系統CPU占用率變化趨勢作圖得到圖5和圖6。

表1 不同線程數量下的計算時長和CUP使用率

圖5 不同線程數下的計算時長

圖6 不同線程數下的總CPU資源使用率

3 測試結果分析

首先對結果進行定性分析。由圖5可見，隨著并行線程數的增加，桌面計算時間快速下降，并以接近負指數趨勢收斂到定值，其中在1至16個線程階段，桌面計算時間的下降基本呈線性趨勢，并行加速效果最為明顯。由圖6可見，系統CPU資源使用效率隨著線程數增長呈接近線性增長的趨勢，當線程數小于64時，CPU占用率維持較為穩定的狀態，且線程數越少，CPU占用率越穩定。

從定量上對加速倍數NS和加速比S進行分析。一般可以直接以式(1)得出實際加速倍數：

NS=Ct1/Ct2

(1)

式中：n表示并行線程的個數。由式(1)可得，在256個并行線程時，加速倍數為100.7倍。

為考察大型機上的極限加速比，可采用阿姆達赫(Amdahl)法則對并行加速性能和程序并行、串行耗時比例進行分析，以式(2)得出加速比：

(2)

式中：f是程序中可以并行執行部分所占用的計算時間在整個程序執行時間中的比例；n是并行線程數量；S是理論上可以獲得的最大加速倍數。

然而，采用式(2)對表1的數據進行分析時，首先需要得到f值。此處通過趨勢近似的方法得到，即設1個線程時，程序將耗時1，則串行部分執行比例s和并行部分執行比例f滿足：

s+f=1

(3)

而當n個線程時，可先設s部分耗時不變，而并行部分耗時將減小為原來的1/n，則s和f滿足：

(4)

由式(3)、式(4)可以得到：

(5)

將表1中的數據代入式(5)，得出一個f值序列，收斂時f≈0.993，對f序列作圖，得到圖7的可并行部分比例圖。

圖7 MC算法中并行部份占比隨線程數的變化

由圖7可以見，在源粒子數較多時，MC模擬的f部分在線程數較大時基本維持在99.3%～99.4%之間。取對比f值分別為99.3%和99.4%，代入式(2)，得到如圖8所示的理想與實際加速倍數(SUT)對比。

圖8 實際加速比與理想加速比的比較

4 結 語

采用OpenMP并行技術，對MC模擬程序進行了并行化， 并在ORACLE M9000大型SMP服務器上進行了實驗計算，結果顯示并行計算結果與串行計算結果完全一致，并且并行程序能夠獲得極好的加速性能，極大縮短桌面計算時間。同時通過對測試結果的分析，可以得知MC模擬中，當抽樣粒子數較大時，其可并行執行部分的比例為99.3%左右。以上結論為工程計算中提高MC模擬的時間效率給出了可行的并行方案，并為工期的制定提供了準確的參考。

[1] 許淑艷.蒙特卡洛方法在實驗核物理中的應用[M].原子能出版社,2010:7-8.

[2] Balakrishnan S,Rajwar R,Upton M,et al.The Impact of Performance Asymmetry in Emerging Multicore Architectures[C]//32st International Symposium on Computer Architecture (ISCA 2005),4-8 June 2005,Madison,Wisconsin,USA.2005:506-517.

[3] Oracle Inc.OpenMP Application Program Interface(Version 3.0)[M].2008.

[4] Sun Microsystems Inc.OpenMP API用戶指南[M].2005.

[5] 王嫻,王秋旺,陶文銓,等.直接模擬蒙特卡羅方法的并行方案設計[J].西安交通大學學報,2003,37(1):105-107.

[6] 陳國良.并行計算[M].高等教育出版社,2003.

[7] 唐兵,Laurent BOBELIN,賀海武.基于MPI和OpenMP混合編程的非負矩陣分解并行算法[J].計算機科學,2017,44(3):51-54.

[8] 杜根遠,張火林,苗放.一種基于MPI和OpenMP的剖分遙感影像并行分割方法[J].計算機應用與軟件,2016,33(9):180-183.