蒙特卡羅臨界計(jì)算全局計(jì)數(shù)問(wèn)題新策略研究*

上官丹驊 姬志成 鄧力 李瑞 李剛 付元光

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

2)(中國(guó)工程物理研究院高性能數(shù)值模擬軟件中心,北京 100088)

1 引 言

計(jì)算核相關(guān)系統(tǒng)的臨界性在核科學(xué)與工程領(lǐng)域是常見(jiàn)而重要的問(wèn)題.作為一種成熟的計(jì)算方法,蒙特卡羅算法由于其幾何、物理建模的精確性以及高度的可并行性等特點(diǎn)在臨界計(jì)算中得到廣泛的運(yùn)用.

在最近的研究[1-4]中,由于計(jì)算機(jī)硬件與模擬計(jì)算方法、軟件的巨大進(jìn)步,使得過(guò)去由于能力不足而采取的近似手段有了可以進(jìn)一步精確化處理的能力,例如反應(yīng)堆設(shè)備的pin-by-pin計(jì)算等.這些趨勢(shì)反過(guò)來(lái)對(duì)于計(jì)算方法和軟件提出了更大的挑戰(zhàn).對(duì)于臨界問(wèn)題的蒙特卡羅模擬而言,由于核設(shè)備的精細(xì)化建模及多物理耦合計(jì)算需求,在計(jì)算臨界特性(即特征值)的同時(shí)進(jìn)行大規(guī)模的全局計(jì)數(shù)成為一個(gè)難以高效解決的問(wèn)題[5-7].困難一方面來(lái)自于計(jì)數(shù)規(guī)模的龐大;另一方面來(lái)自于計(jì)算目標(biāo)不能僅僅瞄準(zhǔn)于最大程度地減少最小誤差,而是必須提高全局計(jì)數(shù)的整體效率,這就使得單純?cè)黾訕颖緮?shù)是一種低效甚至無(wú)效的策略,因?yàn)樵诠β瘦^高、誤差較小的區(qū)域更容易得到新增樣本,而功率較低、誤差較大的區(qū)域反而得不到新增樣本,從而使得整體效率難以提高.

為提高臨界計(jì)算全局計(jì)數(shù)問(wèn)題的整體效率,出現(xiàn)了一些算法上的研究.這些研究以95/95準(zhǔn)則[5](即要求在一定時(shí)間內(nèi),至少95%的柵元內(nèi)計(jì)數(shù)的相對(duì)誤差要以至少95%的置信度位于1%以下)為指導(dǎo),力圖在總樣本數(shù)一定的前提下獲得更高的整體效率,其中均勻裂變點(diǎn)(UFS)算法[5-7]是最早提出的一種.該算法在假設(shè)裂變?cè)幢菊鞣植家呀?jīng)達(dá)到的前提下,通過(guò)統(tǒng)計(jì)已完成迭代步的裂變點(diǎn)數(shù)密度,以此來(lái)偏倚當(dāng)前迭代步的裂變次級(jí)平均粒子數(shù),然后通過(guò)糾偏的手段保持結(jié)果的無(wú)偏性,最終的效果是在統(tǒng)計(jì)意義上使得計(jì)算誤差較大的區(qū)域產(chǎn)生更多的低權(quán)重粒子而在計(jì)算誤差較小的區(qū)域產(chǎn)生更少的高權(quán)重粒子,從而使得整體效率得到較大的提高.但是如果在裂變?cè)捶植歼€未收斂的情況下啟動(dòng)UFS算法,由于統(tǒng)計(jì)得到的裂變點(diǎn)密度具有系統(tǒng)誤差,將會(huì)導(dǎo)致算法應(yīng)用失敗,解決這一問(wèn)題的方法一般是憑經(jīng)驗(yàn)選擇一個(gè)足夠大的非激活迭代步數(shù)來(lái)保證裂變?cè)捶植嫉氖諗?同時(shí),由于只在完成所有的設(shè)定迭代步計(jì)算后才統(tǒng)計(jì)整體效率,即使計(jì)算早已達(dá)到事先設(shè)定的整體效率標(biāo)準(zhǔn),也不得不進(jìn)行多余的計(jì)算.在后續(xù)的研究中還出現(xiàn)了效率更高的均勻計(jì)數(shù)密度算法[8,9].

本文提出了一種進(jìn)一步提高臨界計(jì)算全局計(jì)數(shù)問(wèn)題整體效率的新策略.這一策略在自主研發(fā)的定態(tài)蒙特卡羅粒子輸運(yùn)模擬軟件JMCT上進(jìn)行了驗(yàn)證.本文第2節(jié)將介紹包含兩個(gè)部分的新策略;第3節(jié)給出了相應(yīng)的數(shù)值結(jié)果;第4節(jié)給出了結(jié)論.

2 提高臨界計(jì)算全局計(jì)數(shù)問(wèn)題效率的新策略描述

這一新策略包含兩個(gè)部分.首先,基于一個(gè)裂變?cè)捶植紝?duì)應(yīng)香農(nóng)熵的實(shí)時(shí)收斂性診斷方法,為提高全局計(jì)數(shù)整體效率而設(shè)計(jì)的UFS算法將在首次激活迭代步和首次判斷收斂迭代步的最大值之后被激活,這就為保證裂變?cè)捶植家呀?jīng)收斂提供了雙重保障,從而保證了UFS算法所使用的偏倚數(shù)據(jù)更加合理;其次,在UFS算法被啟動(dòng)后,將定期監(jiān)測(cè)全局計(jì)數(shù)的一個(gè)整體精度指標(biāo),一旦這個(gè)指標(biāo)小于事先約定的數(shù)值,則認(rèn)為全局計(jì)數(shù)已經(jīng)很好地收斂,整個(gè)計(jì)算將被終止,而不必等到事先設(shè)定的最大迭代步數(shù)全部完成.下面將詳細(xì)介紹這兩部分.

2.1 基于香農(nóng)熵的裂變?cè)捶植紝?shí)時(shí)收斂性診斷方法

香農(nóng)熵最先引入蒙特卡羅臨界計(jì)算是作為判斷裂變?cè)捶植际欠袷諗康暮篁?yàn)指標(biāo),同時(shí),也存在另一些基于其他熵的實(shí)時(shí)收斂性診斷方法研究[10-12].基于香農(nóng)熵的實(shí)時(shí)收斂性診斷方法在減少非定常輸運(yùn)問(wèn)題的計(jì)算時(shí)間方面獲得了成功的應(yīng)用[13].將這一實(shí)時(shí)收斂性診斷方法應(yīng)用在定態(tài)蒙特卡羅臨界計(jì)算方面還屬首次.

裂變?cè)捶植紝?duì)應(yīng)香農(nóng)熵H的基本定義如下[10]：

基于香農(nóng)熵的實(shí)時(shí)收斂性判斷準(zhǔn)則依賴(lài)于多個(gè)香農(nóng)熵值組合出的隨機(jī)振子指標(biāo)Kn,定義為[11]

其中,Hn是當(dāng)前第n個(gè)迭代步的香農(nóng)熵值,和分別是當(dāng)前第n迭代步之前所有p個(gè)香農(nóng)熵中的最大值和最小值.在香農(nóng)熵值序列已經(jīng)收斂的前提下,Kn將在0.5附近做無(wú)規(guī)隨機(jī)漲落.實(shí)時(shí)判斷香農(nóng)熵是否收斂的準(zhǔn)則在于判斷當(dāng)前第n迭代步及之前所有m個(gè)隨機(jī)振子指標(biāo)是否滿(mǎn)足不等式[13]

一旦滿(mǎn)足,則認(rèn)為裂變?cè)捶植家呀?jīng)收斂,可以在n和事先約定的非激活迭代步的最大值之后啟動(dòng)計(jì)數(shù)過(guò)程及UFS算法;p,m,ε按參考文獻(xiàn)[11]的推薦取為20,50,0.1.

2.2 基于全局計(jì)數(shù)整體精度指標(biāo)定期監(jiān)測(cè)的提前結(jié)束計(jì)算方法

衡量全局計(jì)數(shù)整體效率的95/95準(zhǔn)則的本質(zhì)在于以盡可能少的計(jì)算時(shí)間獲得盡可能高的整體精度.整體精度指標(biāo)其實(shí)有很多種,都從一個(gè)方面反映了全局計(jì)數(shù)的整體收斂狀況.95/95準(zhǔn)則使用的整體精度指標(biāo)是Pre_95,即所有計(jì)數(shù)的相對(duì)誤差中至少95%的相對(duì)誤差都不大于該值.要在并行計(jì)算環(huán)境下定期監(jiān)測(cè)這一整體精度指標(biāo),就必須在激活UFS算法后每隔固定迭代步數(shù)計(jì)算一次所有計(jì)數(shù)的相對(duì)誤差,而通常來(lái)說(shuō),并行蒙特卡羅程序需要進(jìn)行規(guī)約之后才能得到計(jì)數(shù)的相對(duì)誤差,一般僅在所有迭代步完成之后進(jìn)行規(guī)約并輸出計(jì)數(shù).為解決這一問(wèn)題,JMCT采用了對(duì)象序列化與在線(xiàn)反序列化技術(shù)(由于JMCT采用面向?qū)ο蟮木幊谭绞綄?shí)現(xiàn),程序運(yùn)行時(shí)將產(chǎn)生多個(gè)對(duì)象實(shí)例,其內(nèi)存排布是無(wú)序的.但重啟動(dòng)及續(xù)算功能需要將對(duì)象保存在文件中,需要內(nèi)存對(duì)象是有序且連續(xù)的.在數(shù)據(jù)備份時(shí)對(duì)無(wú)序?qū)ο筮M(jìn)行的操作就是序列化,在恢復(fù)運(yùn)行時(shí)對(duì)數(shù)據(jù)的操作就是反序列化),以最小的I/O開(kāi)銷(xiāo)實(shí)現(xiàn)了并行環(huán)境的備份與恢復(fù),從而支持在任意迭代步中進(jìn)行計(jì)數(shù)規(guī)約操作.

3 數(shù)值結(jié)果

以C5G7模型(圖1)作為驗(yàn)證上述策略的基準(zhǔn)模型.該模型是NEA發(fā)布的程序檢驗(yàn)基準(zhǔn)例題.其中包含兩種組件,分別是鈾氧化物UOX組件和鈾钚混合氧化物MOX組件,交叉2×2布置.每個(gè)組件內(nèi)有17×17個(gè)pin-cell,其中又分為264個(gè)燃料pin、24個(gè)導(dǎo)向管pin和一個(gè)儀表管pin.UOX組件中,包含一種富集度;MOX組件中包含三種富集度.該模型大約包含2萬(wàn)個(gè)柵元,其詳細(xì)介紹可以參考文獻(xiàn)[14].

以全局體平均通量計(jì)數(shù)作為計(jì)算目標(biāo).由于沒(méi)有標(biāo)準(zhǔn)答案,所以,通過(guò)大規(guī)模1000迭代步的計(jì)算(每代5百萬(wàn)樣本,跳過(guò)前500代)獲得基準(zhǔn)解,這一計(jì)算規(guī)模對(duì)于該問(wèn)題是足夠的.圖2所示裂變?cè)捶植紝?duì)應(yīng)香農(nóng)熵的變化反映了這一點(diǎn).分四種情況進(jìn)行了計(jì)算：1)不啟動(dòng)策略和UFS算法(簡(jiǎn)稱(chēng)為no_stra_no_UFS);2)不啟動(dòng)策略、啟動(dòng)UFS算法(簡(jiǎn)稱(chēng)為no_stra_with_UFS);3)啟動(dòng)策略、不啟動(dòng)UFS算法(簡(jiǎn)稱(chēng)為with_stra_no_UFS);4)啟動(dòng)策略和UFS算法(簡(jiǎn)稱(chēng)為with_stra_with_UFS).最初的設(shè)置是1500迭代步,跳過(guò)200代,每代20萬(wàn)樣本數(shù).圖3顯示了不啟動(dòng)策略時(shí)有無(wú)UFS算法裂變?cè)捶植妓鶎?duì)應(yīng)香農(nóng)熵的變化,可以看出,200代的非激活迭代數(shù)目是不夠的,也就是UFS算法利用了不夠精確的數(shù)據(jù)進(jìn)行偏倚.為考察全局計(jì)數(shù)的整體精度,設(shè)計(jì)了整體精度指標(biāo)E_global,定義如下：

圖1 C5G7模型圖Fig.1.C5G7 model.

圖2 基準(zhǔn)解裂變?cè)捶植紝?duì)應(yīng)香農(nóng)熵的變化Fig.2.Shannon entropy of fission source distribution for benchmark result.

圖3 不啟動(dòng)策略時(shí)有無(wú)UFS算法的裂變?cè)捶植紝?duì)應(yīng)香農(nóng)熵的變化Fig.3.Shannon entropy of fission source distribution for two cases (without strategy and with or without UFS algorithm).

當(dāng)啟動(dòng)策略時(shí),原先約定的在第201步啟動(dòng)計(jì)數(shù)或UFS算法被推遲到第224步,顯示香農(nóng)熵的實(shí)時(shí)診斷方法判斷直到第224步裂變?cè)捶植疾攀諗康帽容^充分.如果要求只要整體精度指標(biāo)Pre_95不大于0.004就結(jié)束整個(gè)計(jì)算,且每隔40步判斷一次,則不啟動(dòng)UFS算法時(shí)在第944步達(dá)到要求,而啟動(dòng)UFS算法時(shí)在第824步就達(dá)到了要求.從表1可以看出,各種情況下計(jì)算得到的特征值keff基本沒(méi)有什么變化,說(shuō)明各種方法得到的結(jié)果都是無(wú)偏的.由于不啟動(dòng)策略時(shí)UFS算法啟動(dòng)得過(guò)早(在第201步啟動(dòng)),no_stra_with_UFS這種情況下的整體計(jì)算精度E_global是較差的,即使其已經(jīng)完成了全部1500個(gè)迭代步的計(jì)算.對(duì)于with_stra_no_UFS情形,雖然啟動(dòng)計(jì)數(shù)過(guò)程較晚,裂變?cè)捶植际諗康靡呀?jīng)比較充分,但由于總共只有944個(gè)迭代步,有效樣本較no_stra_no_UFS情形少很多,又沒(méi)有UFS算法使全局計(jì)數(shù)的相對(duì)誤差平均化,所以整體精度指標(biāo)E_global是最差的,即使此時(shí)指標(biāo)Pre_95已經(jīng)小于0.004.而對(duì)于with_stra_with_UFS情形,雖然結(jié)束計(jì)算的時(shí)刻更早(在第824步),從而有效樣本數(shù)更少,但由于啟動(dòng)UFS算法對(duì)全局計(jì)數(shù)誤差的平均化效應(yīng),其整體精度指標(biāo)E_global是最好的.

表1 結(jié)果比較Table 1.Comparison of results.

4 結(jié) 論

對(duì)于蒙特卡羅臨界計(jì)算全局計(jì)數(shù)問(wèn)題,基于已有的UFS算法,提出了一種進(jìn)一步提高效率的新策略.該策略通過(guò)基于裂變?cè)捶植紝?duì)應(yīng)香農(nóng)熵的實(shí)時(shí)判斷收斂準(zhǔn)則和對(duì)全局計(jì)數(shù)整體精度指標(biāo)的定期監(jiān)控,可進(jìn)一步確保UFS算法獲得質(zhì)量更高的數(shù)據(jù)并且在保證精度的前提下減少冗余計(jì)算,從而提高了整體效率.對(duì)C5G7基準(zhǔn)模型的計(jì)算表明本文提出的策略在參數(shù)合適時(shí)是高效的.當(dāng)然,由于新策略中還有經(jīng)驗(yàn)參數(shù),利用更智能化的方法消除這些參數(shù)是下一步工作的目標(biāo).