蒙特卡洛模擬的并行實(shí)現(xiàn)及并行效率研究

柳 建 孫勝軍 毛國(guó)平 石秀安

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

0 引 言

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

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

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

圖1 中子跟蹤流程

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

圖2 MC模擬代碼流程框圖

通過(guò)以上的程序流程可知，MC算法中的中子的運(yùn)動(dòng)歷史之間沒(méi)有數(shù)據(jù)相關(guān)性，因而非常適合并行計(jì)算[5]以縮短計(jì)算所需桌面時(shí)間。

當(dāng)前的高效并行技術(shù)主要有OpenMP、MPI和PVM[5-8]。其中OpenMP是支持共享存儲(chǔ)并行方式的庫(kù)，適合在多核CPU上進(jìn)行并行，采用共享存儲(chǔ)系統(tǒng)，編程和執(zhí)行效率都很高，因而得到了廣泛應(yīng)用。而在蒙特卡羅模擬中，在算法上，源粒子經(jīng)抽樣后，它們?cè)谟?jì)算中是數(shù)據(jù)無(wú)關(guān)的，因而可以采用任意并行方式。由于改進(jìn)的程序?qū)⒅饕\(yùn)行在共享存儲(chǔ)的SMP結(jié)構(gòu)硬件平臺(tái)上，同時(shí)為減小線程間的通信數(shù)據(jù)量、獲得好的負(fù)載平衡，此處采用來(lái)了OpenMP并行庫(kù)為基礎(chǔ)的數(shù)據(jù)并行方式[3-4]。

對(duì)圖2中的串行代碼進(jìn)行多線程并行化時(shí)，為提高并行效率，一般需要進(jìn)行粗粒度的并行，以減少線程喚醒、掛起和等待的時(shí)間。為此將并行區(qū)域設(shè)置為流程圖中最大的循環(huán)區(qū)域，即以每個(gè)源粒子的歷史過(guò)程作為劃分線程粒度的標(biāo)準(zhǔn)。程序中通過(guò)將該區(qū)域使用OpenMP的parallel關(guān)鍵字來(lái)指定為并行區(qū)域。為了各線程總體計(jì)算量負(fù)載平衡，通常每個(gè)線程每次抽取1個(gè)源粒子。由于每個(gè)線程在累計(jì)已抽樣粒子數(shù)、調(diào)用截面數(shù)據(jù)、匯總計(jì)數(shù)結(jié)果、累加運(yùn)行錯(cuò)誤等情形下會(huì)使用到公共數(shù)據(jù)，這將導(dǎo)致線程之間會(huì)出現(xiàn)數(shù)據(jù)競(jìng)爭(zhēng)，為避免這個(gè)問(wèn)題，在這些地方通過(guò)critical關(guān)鍵字或atomic關(guān)鍵字加以防止。

最終的并行程序結(jié)構(gòu)中具體采用了如圖3的fork & join并行方式[8]，思路和數(shù)據(jù)流程簡(jiǎn)潔清晰。

圖3 采用OpenMP方式并行的程序結(jié)構(gòu)圖

2 測(cè)試算例及并行結(jié)果

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

圖4 反應(yīng)堆外殼幾何模型

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

測(cè)試結(jié)果如表1所示。其中，Ct是通過(guò)直接記錄程序運(yùn)行完成時(shí)間得到，Usg是通過(guò)使用M9000上的SOLARIS操作系統(tǒng)的“prstat”命令統(tǒng)計(jì)得到。通過(guò)表1，對(duì)時(shí)耗和系統(tǒng)CPU占用率變化趨勢(shì)作圖得到圖5和圖6。

表1 不同線程數(shù)量下的計(jì)算時(shí)長(zhǎng)和CUP使用率

圖5 不同線程數(shù)下的計(jì)算時(shí)長(zhǎng)

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

3 測(cè)試結(jié)果分析

首先對(duì)結(jié)果進(jìn)行定性分析。由圖5可見，隨著并行線程數(shù)的增加，桌面計(jì)算時(shí)間快速下降，并以接近負(fù)指數(shù)趨勢(shì)收斂到定值，其中在1至16個(gè)線程階段，桌面計(jì)算時(shí)間的下降基本呈線性趨勢(shì)，并行加速效果最為明顯。由圖6可見，系統(tǒng)CPU資源使用效率隨著線程數(shù)增長(zhǎng)呈接近線性增長(zhǎng)的趨勢(shì)，當(dāng)線程數(shù)小于64時(shí)，CPU占用率維持較為穩(wěn)定的狀態(tài)，且線程數(shù)越少，CPU占用率越穩(wěn)定。

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

NS=Ct1/Ct2

(1)

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

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

(2)

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

然而，采用式(2)對(duì)表1的數(shù)據(jù)進(jìn)行分析時(shí)，首先需要得到f值。此處通過(guò)趨勢(shì)近似的方法得到，即設(shè)1個(gè)線程時(shí)，程序?qū)⒑臅r(shí)1，則串行部分執(zhí)行比例s和并行部分執(zhí)行比例f滿足：

s+f=1

(3)

而當(dāng)n個(gè)線程時(shí)，可先設(shè)s部分耗時(shí)不變，而并行部分耗時(shí)將減小為原來(lái)的1/n，則s和f滿足：

(4)

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

(5)

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

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

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

圖8 實(shí)際加速比與理想加速比的比較

4 結(jié) 語(yǔ)

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

[1] 許淑艷.蒙特卡洛方法在實(shí)驗(yàn)核物理中的應(yīng)用[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] 王嫻,王秋旺,陶文銓,等.直接模擬蒙特卡羅方法的并行方案設(shè)計(jì)[J].西安交通大學(xué)學(xué)報(bào),2003,37(1):105-107.

[6] 陳國(guó)良.并行計(jì)算[M].高等教育出版社,2003.

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

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