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

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

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

1 引 言

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

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

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

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

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

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

2.1 基于香農(nóng)熵的裂變源分布實時收斂性診斷方法

香農(nóng)熵最先引入蒙特卡羅臨界計算是作為判斷裂變源分布是否收斂的后驗指標,同時,也存在另一些基于其他熵的實時收斂性診斷方法研究[10-12].基于香農(nóng)熵的實時收斂性診斷方法在減少非定常輸運問題的計算時間方面獲得了成功的應用[13].將這一實時收斂性診斷方法應用在定態(tài)蒙特卡羅臨界計算方面還屬首次.

裂變源分布對應香農(nóng)熵H的基本定義如下[10]：

基于香農(nóng)熵的實時收斂性判斷準則依賴于多個香農(nóng)熵值組合出的隨機振子指標Kn,定義為[11]

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

一旦滿足,則認為裂變源分布已經(jīng)收斂,可以在n和事先約定的非激活迭代步的最大值之后啟動計數(shù)過程及UFS算法;p,m,ε按參考文獻[11]的推薦取為20,50,0.1.

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

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

3 數(shù)值結果

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

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

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

圖2 基準解裂變源分布對應香農(nóng)熵的變化Fig.2.Shannon entropy of fission source distribution for benchmark result.

圖3 不啟動策略時有無UFS算法的裂變源分布對應香農(nóng)熵的變化Fig.3.Shannon entropy of fission source distribution for two cases (without strategy and with or without UFS algorithm).

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

表1 結果比較Table 1.Comparison of results.

4 結 論

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