多物理耦合計算中動態(tài)輸運(yùn)問題高效蒙特卡羅模擬方法*

上官丹驊 閆威華 魏軍俠 高志明 陳藝冰 姬志成

(北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100094)

多物理耦合計算在眾多領(lǐng)域都有重要應(yīng)用.如果其包含粒子輸運(yùn)過程,用蒙特卡羅方法模擬粒子輸運(yùn)常占據(jù)大部分的計算時間,因此多物理耦合計算中動態(tài)輸運(yùn)問題的高效蒙特卡羅模擬方法意義重大,其不可避免地依賴于大規(guī)模并行.基于動態(tài)輸運(yùn)問題的特點(diǎn),本文提出了兩種新方法:一是針對輸運(yùn)燃耗耦合計算的新型計數(shù)規(guī)約算法;二是動態(tài)輸運(yùn)計算樣本數(shù)自適應(yīng)算法.兩種算法都能在保持計算結(jié)果基本不變的前提下使計算時間大幅減少,從而提高了效率.

1 引言

在定量討論各種核系統(tǒng)(如反應(yīng)堆等)的動態(tài)行為時,中子的時-空及能量分布是很關(guān)鍵的.中子與原子核反應(yīng)導(dǎo)致核能釋放,從而導(dǎo)致介質(zhì)溫度、密度和運(yùn)動速度等物理量發(fā)生變化;反過來,因為中子與核的反應(yīng)依賴于二者之間的相對運(yùn)動能量,從而介質(zhì)中原子核熱運(yùn)動和流體運(yùn)動將影響中子的輸運(yùn)過程.基于物理上的考慮,可以用Maxwell分布來描述核的熱運(yùn)動.在對核的密度分布函數(shù)作為權(quán)重求平均以后,中子輸運(yùn)方程中的系數(shù)將只依賴于中子與流體運(yùn)動之間的相對速度.基于相空間中子數(shù)守恒原理,可得考慮流體宏觀運(yùn)動的中子輸運(yùn)方程[1]:

其中V=v?u是相對速度,v是中子速度,u是流體速度,N(r,V,t)drdV是t時刻r處體積元 dr內(nèi)相對速度V處 dV范圍內(nèi)的中子可幾數(shù),Σ′f′V ′N′=Σ(r,V′)f(r,V′→V)V′N(r,V′,t),Σ(r,V)是宏觀總截面,是一個中子在單位距離內(nèi)與i類核經(jīng)受某一由x標(biāo)志的特定反應(yīng)的概率,表示r處速度為V′的中子與i類核發(fā)生x類反應(yīng)產(chǎn)生速度為V的中子的概率,q(r,V,t)drdV表示r處 dr體積元內(nèi)產(chǎn)生的速度為V處 dV范圍內(nèi)的中子數(shù).嚴(yán)格來說,由于宏觀截面依賴于中子通量,所以上述方程中宏觀截面是隨時間變化的.所以,輸運(yùn)方程其實(shí)是一個非線性的微分-積分方程.由于流體運(yùn)動的時間尺度較中子運(yùn)動的時間尺度長,所以可選一個時間間隔,其相對于中子輸運(yùn)相當(dāng)長,而對于流體運(yùn)動來說相當(dāng)短,從而在該間隔內(nèi)可以忽略核數(shù)密度改變的影響.在時間間隔末通過解線性輸運(yùn)方程得到中子角通量以后,就可以利用所得中子角通量、流體力學(xué)方程組和燃耗方程組求解新的物質(zhì)密度和核數(shù)密度.這樣就把非線性輸運(yùn)方程轉(zhuǎn)化為多時間步的線性輸運(yùn)問題,形成了多物理耦合計算(本文僅考慮流體力學(xué)、輸運(yùn)和燃耗的多物理耦合計算).

隨著科學(xué)技術(shù)的發(fā)展,特別是并行計算能力的飛躍,粒子輸運(yùn)問題的建模(包括幾何建模、物理建模等)越來越逼近真實(shí)情況,以前由于計算能力不足而不得不采取的各種近似逐漸被替換,這就給蒙特卡羅粒子輸運(yùn)模擬帶來了巨大的挑戰(zhàn).以反應(yīng)堆pin-by-pin 模型蒙特卡羅粒子輸運(yùn)模擬領(lǐng)域著名的Kord-Smith 挑戰(zhàn)[2]為例,由于其高達(dá)百萬以上柵元,幾十億數(shù)量級的計數(shù)規(guī)模的設(shè)置,使得蒙特卡羅模擬面臨存不下、算不快和算不準(zhǔn)的難題.對于多物理耦合計算中的動態(tài)蒙特卡羅粒子輸運(yùn)模擬而言,其面臨挑戰(zhàn)的難度相較Kord-Smith 挑戰(zhàn)而言只高不低.這主要是由動態(tài)輸運(yùn)問題高達(dá)千萬以上網(wǎng)格的幾何建模、由于極端反應(yīng)條件造成的數(shù)十萬以上的計算步、由于輸運(yùn)與燃耗耦合而必須求解的大量網(wǎng)格計數(shù)所造成的,而且動態(tài)輸運(yùn)問題由于時間累積效應(yīng),定態(tài)輸運(yùn)蒙特卡羅模擬所常用的許多技巧并不能直接應(yīng)用,需要根據(jù)動態(tài)問題的特點(diǎn)設(shè)計新的算法.

在公開發(fā)表的文獻(xiàn)中,反應(yīng)堆大型pin-bypin 模型的高效蒙特卡羅模擬方法呈現(xiàn)出蓬勃發(fā)展的趨勢.區(qū)域分解[3,4]和數(shù)據(jù)分解[5,6]方法的提出主要是為了解決存不下的問題;為了解決算不快的問題,提出了若干高效樣本數(shù)并行算法;為了解決全局計數(shù)算不準(zhǔn)的問題,除了加大樣本并高效并行以外,還提出了UFS (uniform fission site method)[7,8]和UTD (uniform tally density method)[9,10]等算法(這些算法都屬于偏倚算法).表1 列出了常見的大型反應(yīng)堆pin-by-pin 模型的蒙特卡羅模擬方法對于存不下、算不快和算不準(zhǔn)這3 個難點(diǎn)各自的優(yōu)缺點(diǎn)分析.其中“√”表示左邊的算法對上面難點(diǎn)的解決有促進(jìn)作用,“×”表示表示左邊的算法對上面難點(diǎn)的解決有阻礙作用,“—”表示沒有影響.框內(nèi)的語句表示的是促進(jìn)或阻礙作用最本質(zhì)的原因.這些研究對于動態(tài)輸運(yùn)問題的高效蒙特卡羅模擬也具有重要的借鑒作用.

從表1 可以看出,各種算法在解決一個問題的同時常引入新的問題.對于本文關(guān)注的動態(tài)輸運(yùn)問題,由于其不同的特點(diǎn),上述算法可以借鑒但不能照搬.鑒于動態(tài)輸運(yùn)蒙特卡羅模擬極大的計算量,設(shè)計高效算法的一個基本思路就是在基本不改變大樣本直接模擬結(jié)果的基礎(chǔ)上盡可能減少計算時間.為此,本文開展了針對輸運(yùn)燃耗耦合計算的計數(shù)規(guī)約優(yōu)化算法和針對樣本并行的樣本數(shù)自適應(yīng)算法研究.這些研究是在北京應(yīng)用物理與計算數(shù)學(xué)研究所自主開發(fā)的極端條件下三維多物理多介質(zhì)問題數(shù)值模擬軟件JUPITER (jointed numerical simulation software for multi-physics and multimaterials problems under extreme conditions)上進(jìn)行的,此平臺對模型采取三維非結(jié)構(gòu)網(wǎng)格建模,在此基礎(chǔ)上進(jìn)行多物理耦合計算.其中的蒙特卡羅粒子輸運(yùn)模擬模塊(以下簡稱為JUPITERMC 程序)以區(qū)域分解加樣本并行為基本架構(gòu),可支持?jǐn)?shù)千萬三維非結(jié)構(gòu)多面體網(wǎng)格、數(shù)百億樣本的中子輸運(yùn)及燃耗耦合計算,具備數(shù)萬處理器核的高效并行計算能力.

表1 常見算法優(yōu)缺點(diǎn)比較Table 1.Relative merits of common algorithms.

本文第2 節(jié)介紹了一種基于輸運(yùn)燃耗耦合計算模式的計數(shù)規(guī)約新算法;第3 節(jié)基于蒙特卡羅模擬分批計算提出了一種樣本數(shù)自適應(yīng)算法;第4 節(jié)的數(shù)值結(jié)果表明兩種算法都可以在基本不改變計算結(jié)果的前提下大幅減少計算時間,從而提高了效率;最后是總結(jié).

2 計數(shù)規(guī)約優(yōu)化算法

2.1 輸運(yùn)燃耗耦合計算模式

對于大規(guī)模動態(tài)輸運(yùn)蒙特卡羅模擬問題,為解決存不下的難題,JUPITER-MC 程序采取了區(qū)域分解加樣本并行的二級并行模式.整個模型進(jìn)行網(wǎng)格分片,每個進(jìn)程只存儲單個網(wǎng)格片(包含若干非結(jié)構(gòu)網(wǎng)格)并分配若干樣本粒子,一旦粒子逃逸出該網(wǎng)格片,就將粒子屬性信息發(fā)送給相應(yīng)的其他進(jìn)程(其存儲的是逃逸粒子所應(yīng)到達(dá)的網(wǎng)格片).在輸運(yùn)燃耗耦合計算時,燃耗方程求解的是網(wǎng)格片上所有網(wǎng)格內(nèi)所有核素新的核數(shù)密度.由于燃耗方程的輸入是網(wǎng)格片上所有網(wǎng)格的多種微觀反應(yīng)率等計數(shù),存儲同一網(wǎng)格片的不同進(jìn)程需要將這些計數(shù)進(jìn)行規(guī)約平均并全局廣播,以此保證求解同一網(wǎng)格片的不同進(jìn)程上的燃耗計算輸入是一致的,由此可以保證求解同一網(wǎng)格片的不同進(jìn)程的燃耗計算的輸出(即新的核數(shù)密度)是一致的(并不需要規(guī)約平均).

2.2 計數(shù)規(guī)約新算法

由2.1 節(jié)所述,同一網(wǎng)格片上所有網(wǎng)格的微觀反應(yīng)率需要規(guī)約平均并全局廣播,而燃耗計算所得核數(shù)密度自然保持一致,并不需要費(fèi)時的規(guī)約平均操作.但是,若網(wǎng)格片所含網(wǎng)格數(shù)目較多,或者存儲同一網(wǎng)格片的進(jìn)程數(shù)較多,需要規(guī)約平均并全局廣播的微觀反應(yīng)率數(shù)據(jù)本身就是網(wǎng)格數(shù)目的多倍,由此導(dǎo)致并行通信的時間是不可忽略的(同一進(jìn)程上,網(wǎng)格片所含網(wǎng)格數(shù)相對于樣本數(shù)越多,計數(shù)規(guī)約所占的時間比例越高).仔細(xì)分析輸運(yùn)燃耗耦合模式,可以發(fā)現(xiàn)燃耗計算的輸出數(shù)據(jù)幾乎必然遠(yuǎn)少于輸入數(shù)據(jù).因此,本文提出了一種計數(shù)規(guī)約的新算法,即燃耗計算只利用本地的微觀反應(yīng)率數(shù)據(jù)進(jìn)行計算,所以其輸入端不需要費(fèi)時的通信,而輸出端(即新的核數(shù)密度)由于不同進(jìn)程不再能保持一致,所以將存儲同一網(wǎng)格片的不同進(jìn)程的燃耗輸出端按相應(yīng)進(jìn)程數(shù)進(jìn)行規(guī)約平均并全局廣播.雖然此種算法只是將費(fèi)時的通信從燃耗輸入端轉(zhuǎn)移到輸出端,但由于通信的數(shù)據(jù)量可以大幅減少,預(yù)期可以從總體上減少通信時間.值得注意的是,這種算法從原理上就不能保證串并行嚴(yán)格一致,但在樣本數(shù)趨于無窮時,串并行結(jié)果是趨于一致且與原算法的結(jié)果是趨于一致的.

3 樣本數(shù)自適應(yīng)算法

3.1 研究基礎(chǔ)

對于動態(tài)輸運(yùn)問題蒙特卡羅模擬進(jìn)行樣本數(shù)自適應(yīng)研究,起源于三個方面的研究基礎(chǔ),即香農(nóng)熵概念引入定態(tài)臨界計算以指示裂變源分布的收斂、臨界計算全局計數(shù)算法的研究和臨界計算全局計數(shù)問題的樣本數(shù)自適應(yīng)算法研究.

3.1.1 香農(nóng)熵

對于定態(tài)輸運(yùn)臨界計算問題,理論上只有當(dāng)裂變源分布已經(jīng)收斂時才能進(jìn)行各種計數(shù),否則結(jié)果是包含系統(tǒng)性誤差的.當(dāng)?shù)嬎氵M(jìn)行時,可以利用信息論中香農(nóng)熵的概念對裂變源分布是否收斂進(jìn)行后驗檢驗.假設(shè)系統(tǒng)包含N個柵元,在某一步迭代計算后產(chǎn)生了M個裂變點(diǎn),每個柵元包含Mi(i=1,2,···,N)個裂變點(diǎn),則裂變源分布對應(yīng)的香農(nóng)熵H為

其中pi=Mi/M.每一迭代步將產(chǎn)生一個熵值,當(dāng)計算完成時,若判斷出K步以后熵序列已經(jīng)收斂而開啟計數(shù)的迭代步號大于等于k,則認(rèn)為所得計數(shù)沒有系統(tǒng)性誤差.當(dāng)然,也存在一些基于熵概念對裂變源分布是否收斂的實(shí)時診斷方法研究[11?13].

3.1.2 臨界計算全局計數(shù)算法

在定態(tài)臨界計算時常常伴隨著全局計數(shù)的求解,由于全局計數(shù)問題目標(biāo)眾多且核系統(tǒng)各部分之間反應(yīng)性的不均勻,全局計數(shù)問題單靠增加樣本數(shù)是低效的.為此,基于對裂變次級中子數(shù)進(jìn)行不同標(biāo)準(zhǔn)的偏倚,產(chǎn)生了UFS 和UTD 等算法,這些算法以若干全局計數(shù)整體效率指標(biāo)進(jìn)行衡量,都大幅提高了計算效率.

3.1.3 臨界計算全局計數(shù)問題的樣本數(shù)自適應(yīng)算法

由于定態(tài)臨界計算時何時啟動全局計數(shù)是憑經(jīng)驗設(shè)置的,且每步樣本數(shù)是固定不變的,所以可能存在過早啟動全局計數(shù)(即裂變源分布還未收斂就啟動全局計數(shù))和冗余計算(即全局計數(shù)結(jié)果已經(jīng)滿足要求還在計算)的問題.文獻(xiàn)[14]提出了一種臨界計算全局計數(shù)問題的樣本數(shù)自適應(yīng)算法,其基于對香農(nóng)熵序列的一種實(shí)時收斂性診斷方法和對一種全局計數(shù)精度指標(biāo)的實(shí)時監(jiān)測,可以保證在更合理的地方啟動全局計數(shù)并且節(jié)省不必要的計算.

3.2 樣本數(shù)自適應(yīng)算法

在國際著名的定態(tài)蒙特卡羅粒子輸運(yùn)模擬軟件MCNP 中,存在根據(jù)某一計數(shù)的相對誤差變化進(jìn)行樣本數(shù)自適應(yīng)調(diào)整的方法.但這一方法并不適合動態(tài)輸運(yùn)蒙特卡羅模擬.因為動態(tài)計算中反饋給其他物理過程的全局反饋量很難確定哪一個重要,無法根據(jù)某個計數(shù)結(jié)果的相對誤差進(jìn)行樣本數(shù)自適應(yīng)調(diào)整.基于定態(tài)輸運(yùn)臨界計算和動態(tài)多步輸運(yùn)計算的類比,即定態(tài)臨界計算的裂變源中子對應(yīng)于動態(tài)輸運(yùn)計算上一步傳遞給本步的未死亡粒子(稱之為時間源粒子),定態(tài)臨界計算的全局計數(shù)對應(yīng)于多物理耦合計算中動態(tài)輸運(yùn)反饋給其他物理過程的全局反饋量,可以設(shè)計動態(tài)輸運(yùn)蒙特卡羅模擬的一種樣本數(shù)自適應(yīng)算法,詳述如下.

3.2.1 基本思路

由于JUPITER-MC 程序是樣本數(shù)分批計算,每批計算都會形成時間源中子和全局反饋量計數(shù).時間源中子形成的時間源分布和全局反饋量計數(shù)對于整個多物理耦合計算來說是最重要的兩個因素.所以,一旦某批次計算完成,根據(jù)某種判斷標(biāo)準(zhǔn)表明時間源分布已經(jīng)抽樣充分且全局反饋量計數(shù)已經(jīng)整體收斂,則可以提前結(jié)束本步的計算,不再進(jìn)行剩余批次的計算.這是多物理耦合計算中動態(tài)輸運(yùn)蒙特卡羅模擬樣本數(shù)自適應(yīng)算法的基本思想.

3.2.2 基本工具

上述基本思路的實(shí)現(xiàn)需要3 個基本工具.一是根據(jù)什么定量指標(biāo)表明時間源分布已經(jīng)抽樣充分;二是根據(jù)什么定量指標(biāo)表明全局反饋量計數(shù)已經(jīng)充分收斂;三是如何根據(jù)批次的增加實(shí)時判斷上述兩個指標(biāo)已經(jīng)收斂.前兩個工具的實(shí)現(xiàn)離不開香農(nóng)熵概念的推廣.第3 個工具則是一種熵值序列的實(shí)時收斂性診斷方法.

由于時間源中子是時間源分布的抽樣,對于某一固定的相空間網(wǎng)格(這里僅考慮空間和能量組成的相空間),假設(shè)共包含N個網(wǎng)格(Si,i=1,2,···,N)、G個能群(gi,i=1,2,···,G)和M個時間源中子,令pij(i=1,2,···,N;j=1,2,···,G)是包含在相空間網(wǎng)格Si×gj中的時間源中子數(shù)除以總時間源中子數(shù)再除以Vi(Vi是相應(yīng)網(wǎng)格的體積),則

Hs是一個指示時間源中子是否已經(jīng)足夠代表時間源分布的類香農(nóng)熵指標(biāo).基于同樣的考慮和全局反饋量都是通量與某響應(yīng)函數(shù)的積分的事實(shí),可設(shè)計指示全局反饋量計數(shù)是否整體收斂充分的類香農(nóng)熵指標(biāo)如下:

其中fi是第i個網(wǎng)格上的體通量除以總通量.注意到上述兩個指標(biāo)是第一次提出用以指示時間源中子是否充分和全局反饋量計數(shù)是否整體收斂,雖然指示全局反饋量計數(shù)是否整體收斂的指標(biāo)不止一個.

對于分批計算,每完成一批計算就會得到新的熵值Hs和Hf.可以設(shè)計熵值序列的實(shí)時收斂性診斷方法如下:對每一種熵值序列,在第n批計算完成后可構(gòu)造隨機(jī)振蕩指標(biāo)Kn,其定義為[15]

則可認(rèn)為對應(yīng)熵值序列已經(jīng)收斂(p,m,ε固定為20,50,0.1).

4 數(shù)值結(jié)果及討論

如引言中所述,本文只關(guān)注力學(xué)、輸運(yùn)和燃耗的多物理耦合計算問題,因此本節(jié)將給出若干此類模型的計算結(jié)果.必須指出的是,同樣的方法是可以應(yīng)用到其他多物理耦合計算過程的.

4.1 計數(shù)規(guī)約優(yōu)化算法結(jié)果

考慮一個包含141 萬非結(jié)構(gòu)六面體網(wǎng)格的模型,該模型包含5 種介質(zhì),其中有含裂變物質(zhì)的材料.對其進(jìn)行多介質(zhì)ALE (arbitrary Lagrange-Euler)力學(xué)、中子輸運(yùn)和燃耗的耦合計算.其中蒙特卡羅粒子輸運(yùn)模擬采取每步1600 萬樣本,160 核并行,不進(jìn)行區(qū)域分解(經(jīng)物理分析,此樣本數(shù)對于本模型是足夠的).由于計數(shù)規(guī)約優(yōu)化算法在減少燃耗輸入端規(guī)約數(shù)據(jù)的同時增加了燃耗輸出端的規(guī)約數(shù)據(jù),所以比較計數(shù)規(guī)約時間和燃耗計算時間之和(包含燃耗輸出端的規(guī)約時間)是合理的.

如圖1 所示,計數(shù)規(guī)約時間和燃耗計算時間之和減少了大約30%.表2 是該模型一個完整的計算結(jié)果,其中把原程序?qū)傖尫拍芰康挠嬎憬Y(jié)果取為標(biāo)準(zhǔn)1.從表2 可以看出,總釋放能量只有約1%的變化,而總計算時間減少了約15%.可見,計數(shù)規(guī)約優(yōu)化算法可以在基本保持原結(jié)果不變的前提下實(shí)現(xiàn)計算時間的減少,從而提高了效率.

圖1 新舊方法計數(shù)與燃耗計算時間之和的對比Fig.1.Comparison of the total time of tally reduce and burnup calculation by using old and new methods.

表2 對于一個包含141 萬非結(jié)構(gòu)六面體網(wǎng)格的模型,全過程計算結(jié)果與計算時間比較Table 2. Comparison of the results and calculation time for the whole simulation for the model including 1.41 million unstructured hexahedral meshes.

4.2 樣本數(shù)自適應(yīng)算法結(jié)果

首先考慮一個包含550 萬非結(jié)構(gòu)六面體網(wǎng)格的模型,含7 種不同的輕重介質(zhì),其中包含裂變材料.對其進(jìn)行定態(tài)迭代計算求解中子時間增殖常數(shù):

其中N(t)是t時刻系統(tǒng)中的中子數(shù).蒙特卡羅模擬迭代100 步,每步3200 萬樣本,64 核并行.表3是其計算結(jié)果及時間.可見,最終的λ只變化了0.5%,而計算時間減少了約57%.

表3 對于一個包含550 萬非結(jié)構(gòu)六面體網(wǎng)格的模型,定態(tài)迭代計算結(jié)果與計算時間比較Table 3. Comparison of the λ and calculation time for the iteration calculation for the model including 5.50 million unstructured hexahedral meshes.

其次考慮一個包含2233 萬非結(jié)構(gòu)六面體網(wǎng)格的模型,同樣含7 種不同的輕重介質(zhì),其中包含裂變材料.對其進(jìn)行多介質(zhì)ALE 力學(xué)、中子輸運(yùn)和燃耗的耦合計算.蒙特卡羅模擬每步1.6 億樣本,8000 核并行,區(qū)域分解為250 個網(wǎng)格片(注意,由于此模型足夠精細(xì),若不進(jìn)行區(qū)域分解,大型機(jī)的內(nèi)存仍然不夠).表4 是其計算結(jié)果及計算時間.其中,原程序的總釋放能量取為標(biāo)準(zhǔn)1.可見,最終的總釋放能量只變化了0.4%,而計算時間減少了約39%.

表4 對于一個包含2233 萬非結(jié)構(gòu)六面體網(wǎng)格的模型,全過程計算結(jié)果與計算時間比較Table 4. Comparison of the results and calculation time for the whole simulation for the model including 22.33 million unstructured hexahedral meshes.

從上面兩個模型的計算結(jié)果可以明顯看出,樣本數(shù)自適應(yīng)算法可以在基本保持原結(jié)果不變的前提下實(shí)現(xiàn)計算時間的大幅減少,從而提高了整體計算效率.

5 結(jié)論

基于多物理耦合計算中動態(tài)輸運(yùn)問題的特點(diǎn),本文在香農(nóng)熵概念的啟發(fā)下通過設(shè)計新的指標(biāo)和利用一種實(shí)時熵值序列診斷方法,實(shí)現(xiàn)了對動態(tài)輸運(yùn)蒙特卡羅模擬的樣本數(shù)進(jìn)行自適應(yīng)調(diào)節(jié)的方法,避免了過多的冗余計算.同時,基于對多物理耦合計算中輸運(yùn)燃耗耦合模式的分析,提出了一種新的MPI 環(huán)境下的計數(shù)規(guī)約算法,可以在保持結(jié)果統(tǒng)計等價的前提下減少必須規(guī)約的數(shù)據(jù)量,從而減少時間.經(jīng)若干大型模型的檢驗,兩種算法都能在保持計算結(jié)果基本不變的前提下使計算時間大幅減少,從而提高了效率.