多物理耦合計算中動態輸運問題高效蒙特卡羅模擬方法*

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

(北京應用物理與計算數學研究所,北京 100094)

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

1 引言

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

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

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

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

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

表1 常見算法優缺點比較Table 1.Relative merits of common algorithms.

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

2 計數規約優化算法

2.1 輸運燃耗耦合計算模式

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

2.2 計數規約新算法

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

3 樣本數自適應算法

3.1 研究基礎

對于動態輸運問題蒙特卡羅模擬進行樣本數自適應研究,起源于三個方面的研究基礎,即香農熵概念引入定態臨界計算以指示裂變源分布的收斂、臨界計算全局計數算法的研究和臨界計算全局計數問題的樣本數自適應算法研究.

3.1.1 香農熵

對于定態輸運臨界計算問題,理論上只有當裂變源分布已經收斂時才能進行各種計數,否則結果是包含系統性誤差的.當迭代計算進行時,可以利用信息論中香農熵的概念對裂變源分布是否收斂進行后驗檢驗.假設系統包含N個柵元,在某一步迭代計算后產生了M個裂變點,每個柵元包含Mi(i=1,2,···,N)個裂變點,則裂變源分布對應的香農熵H為

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

3.1.2 臨界計算全局計數算法

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

3.1.3 臨界計算全局計數問題的樣本數自適應算法

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

3.2 樣本數自適應算法

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

3.2.1 基本思路

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

3.2.2 基本工具

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

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

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

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

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

則可認為對應熵值序列已經收斂(p,m,ε固定為20,50,0.1).

4 數值結果及討論

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

4.1 計數規約優化算法結果

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

如圖1 所示,計數規約時間和燃耗計算時間之和減少了大約30%.表2 是該模型一個完整的計算結果,其中把原程序對總釋放能量的計算結果取為標準1.從表2 可以看出,總釋放能量只有約1%的變化,而總計算時間減少了約15%.可見,計數規約優化算法可以在基本保持原結果不變的前提下實現計算時間的減少,從而提高了效率.

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

表2 對于一個包含141 萬非結構六面體網格的模型,全過程計算結果與計算時間比較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 樣本數自適應算法結果

首先考慮一個包含550 萬非結構六面體網格的模型,含7 種不同的輕重介質,其中包含裂變材料.對其進行定態迭代計算求解中子時間增殖常數:

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

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

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

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

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

5 結論

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