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

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1 研究基礎(chǔ)

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

3.1.1 香農(nóng)熵

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

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

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

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

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

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

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

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

3.2.1 基本思路

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

3.2.2 基本工具

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

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

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

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

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

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

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

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

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

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

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

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

表2 對(duì)于一個(gè)包含141 萬非結(jié)構(gòu)六面體網(wǎng)格的模型,全過程計(jì)算結(jié)果與計(jì)算時(shí)間比較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é)果

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

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

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

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

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

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

5 結(jié)論

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