本發(fā)明涉及一種基于粒子密度不均勻性的迭代蒙卡全局權(quán)窗參數(shù)生成方法，屬于核物理、核技術(shù)應(yīng)用等領(lǐng)域中的輻射屏蔽計(jì)算分析方向。
背景技術(shù)：
越來越多的蒙特卡羅輻射輸運(yùn)計(jì)算需要得到全局收斂的計(jì)算結(jié)果，比如裂變堆的全堆功率分布計(jì)算和聚變堆的停堆劑量計(jì)算。但對于厚屏蔽和幾何復(fù)雜的問題，粒子難以穿透屏蔽區(qū)域，若使用直接模擬法進(jìn)行計(jì)算，不僅屏蔽區(qū)域的統(tǒng)計(jì)誤差會遠(yuǎn)高于其他區(qū)域，計(jì)算耗時也難以接受。
使用網(wǎng)格權(quán)窗是一種常用的蒙特卡羅減方差方法，權(quán)窗是相空間(空間-能量)的分裂和俄羅斯輪盤賭技巧。根據(jù)粒子權(quán)重和粒子當(dāng)前所在網(wǎng)格的權(quán)窗參數(shù)，低權(quán)重粒子通過俄羅斯輪盤賭淘汰，高權(quán)重粒子分裂。這項(xiàng)技術(shù)可以根據(jù)用戶定義的相空間重要性調(diào)整蒙特卡羅粒子分布，同時能夠控制粒子權(quán)重波動，保證粒子權(quán)重分布在一個合理的界限內(nèi)。
手動進(jìn)行權(quán)窗參數(shù)設(shè)置需要大量的使用經(jīng)驗(yàn)，同時難以保證權(quán)窗的加速效果。而傳統(tǒng)的權(quán)窗產(chǎn)生方法雖然在一定程度上輔助生成權(quán)窗參數(shù)，但存在若干缺陷：1)只能優(yōu)化相空間單一位置的計(jì)算效率；2)需要用戶手動迭代優(yōu)化；3)生成的權(quán)窗偏向于將粒子沿著屏蔽最小的路徑輸運(yùn)；4)難以優(yōu)化厚屏蔽區(qū)域的計(jì)算效率。
低質(zhì)量的權(quán)窗參數(shù)不僅不能優(yōu)化計(jì)算效率，還可能導(dǎo)致某些區(qū)域采樣不足，而得到錯誤的計(jì)算結(jié)果。因而如何自適應(yīng)地生成高效的權(quán)窗參數(shù)，是基于蒙特卡羅的復(fù)雜模型輸運(yùn)計(jì)算的關(guān)鍵技術(shù)之一。
目前國際上主流的自動權(quán)窗生成方法主要有兩種，一種為FW-CADIS(Forward Weighted-Consistent Adjoint Driven Importance Sampling)方法，通過確定論程序進(jìn)行伴隨計(jì)算得出每個柵元的重要性，再根據(jù)柵元的重要性生成權(quán)窗來加速計(jì)算。FW-CADIS方法的主要不足為蒙特卡羅輸運(yùn)計(jì)算程序的用戶擁有需要高效的確定論輸運(yùn)計(jì)算程序，并且熟悉確定論程序的使用方法，才能根據(jù)FW-CADIS方法進(jìn)行輸運(yùn)加速。另外一種方法為MAGIC(Method of Automatic Generation of Importance by Calculation)方法，通過迭代計(jì)算正向通量，用正向通量的倒數(shù)作為權(quán)窗參數(shù)加速輸運(yùn)計(jì)算。其主要不足在于權(quán)窗參數(shù)生成時間相對較長，且適用性沒有FW-CADIS方法廣泛。
技術(shù)實(shí)現(xiàn)要素：
本發(fā)明技術(shù)解決問題：克服傳統(tǒng)的需依靠用戶經(jīng)驗(yàn)進(jìn)行優(yōu)化，同時無法對全局的計(jì)算效率做優(yōu)化的權(quán)窗產(chǎn)生器的不足，提供一種基于粒子密度不均勻性的迭代蒙卡全局權(quán)窗參數(shù)生成方法，用戶不需要依賴經(jīng)驗(yàn)手動調(diào)試權(quán)窗，提高了復(fù)雜模型全局蒙卡計(jì)算的效率。
本發(fā)明的技術(shù)方案如下：一種基于粒子密度不均勻性的迭代蒙卡全局權(quán)窗參數(shù)生成方法，包括以下步驟：
(1)生成初始權(quán)窗網(wǎng)格：
a)建立覆蓋整個模型空間的網(wǎng)格；
b)將所有網(wǎng)格的權(quán)窗參數(shù)都初始化為0；
(2)在粒子輸運(yùn)計(jì)算過程中，統(tǒng)計(jì)進(jìn)入每個網(wǎng)格的粒子數(shù)量(Ni，i指第i個網(wǎng)格)，粒子權(quán)重(Wi)，和每個網(wǎng)格對粒子密度均勻性的貢獻(xiàn)，以及所有網(wǎng)格的平均粒子數(shù)量(Na)，在每個粒子輸運(yùn)過程中需要：
a)記錄粒子的徑跡信息，即粒子所穿過的所有網(wǎng)格的編號；
b)當(dāng)粒子第一次進(jìn)入網(wǎng)格i后，如果粒子權(quán)重w大于當(dāng)前所有網(wǎng)格權(quán)窗參數(shù)的最小值，將進(jìn)入網(wǎng)格i的粒子數(shù)量Ni增加1，進(jìn)入網(wǎng)格i的總粒子權(quán)重Wi增加當(dāng)前粒子的權(quán)重w；
c)將粒子每一步輸運(yùn)中，終點(diǎn)所在的網(wǎng)格的粒子數(shù)量Ni和所有網(wǎng)格的平均粒子數(shù)量Na做比較，如果Ni小于Na，則該粒子徑跡上的所有網(wǎng)格的貢獻(xiàn)增加，增加的貢獻(xiàn)即為粒子進(jìn)入一個網(wǎng)格對粒子密度均勻性的貢獻(xiàn)；
d)步驟c)中的對粒子密度均勻性的貢獻(xiàn)的初始值為1.0，并根據(jù)網(wǎng)格粒子密度權(quán)重和粒子輸運(yùn)時間權(quán)重進(jìn)行加權(quán)，粒子密度權(quán)重Fdensity定義為：
公式中Na為所有網(wǎng)格平均粒子數(shù)量，Ni為進(jìn)入網(wǎng)格i的粒子數(shù)量。
粒子輸運(yùn)時間則通過粒子當(dāng)前由于權(quán)窗而分裂的次數(shù)進(jìn)行近似，粒子輸運(yùn)時間權(quán)重Ftime定義為：
公式中Ts為粒子分裂的次數(shù)。
從而對粒子密度的均勻性的貢獻(xiàn)定義為：
ΔCi＝Fdensity×Ftime
(3)計(jì)算權(quán)窗參數(shù)
a)根據(jù)公式
計(jì)算每個網(wǎng)格的重要性Impi。公式中Ci為所有計(jì)算的粒子進(jìn)入網(wǎng)格i對粒子密度均勻性的總貢獻(xiàn)，Wi為進(jìn)入網(wǎng)格i的粒子總權(quán)重；
b)將每個網(wǎng)格的重要性根據(jù)源所在網(wǎng)格的平均重要性進(jìn)行歸一化；
c)將每個網(wǎng)格的權(quán)窗參數(shù)設(shè)為歸一化后的網(wǎng)格重要性的倒數(shù)的一半；
(4)迭代優(yōu)化
a)在權(quán)窗生成過程中，將輸運(yùn)分為若干個周期，根據(jù)前一個周期的網(wǎng)格重要性生成新的權(quán)窗參數(shù)，并應(yīng)用于當(dāng)前周期的輸運(yùn)，同時每個周期生成的權(quán)窗參數(shù)都通過文本進(jìn)行輸出；
b)在每個周期輸運(yùn)開始時，將每個網(wǎng)格的貢獻(xiàn)和進(jìn)入權(quán)重都乘以一個折扣因子f(取0～1之間的值)；
c)迭代自動進(jìn)行，在達(dá)到了用戶設(shè)定的時間上限后停止計(jì)算，當(dāng)前的權(quán)窗參數(shù)即為最終生成的全局權(quán)窗參數(shù)；
本發(fā)明與現(xiàn)有技術(shù)具有如下優(yōu)點(diǎn)：
(1)本發(fā)明自適應(yīng)地生成全局權(quán)窗參數(shù)，用戶只需提供一個覆蓋整個模型的網(wǎng)格，簡化了用戶的操作，而且避免了用戶由于經(jīng)驗(yàn)缺乏而生成的低質(zhì)量權(quán)窗參數(shù)可能導(dǎo)致的錯誤計(jì)算結(jié)果；
(2)本發(fā)明生成的權(quán)窗能極大的提升蒙特卡羅粒子輸運(yùn)在模型全空間的收斂速度，實(shí)現(xiàn)對厚屏蔽和大空間等模型的準(zhǔn)確快速粒子輸運(yùn)。
附圖說明
圖1為本發(fā)明的實(shí)現(xiàn)流程圖。
具體實(shí)施方式
國際熱核聚變實(shí)驗(yàn)堆ITER是目前全球規(guī)模最大的國際科研合作項(xiàng)目之一，ITER裝置是一個能產(chǎn)生大規(guī)模核聚變反應(yīng)的超導(dǎo)托卡馬克。本專利選用ITER官方發(fā)布的中子學(xué)模型Alite.004作為基準(zhǔn)例題。測試用的模型是一個40度的Alite模型，半徑約20米，高約40米。其放射源為完整的聚變等離子體源。
針對該例題，如圖1所示，本發(fā)明的實(shí)施方式為：
(1)生成初始權(quán)窗網(wǎng)格，網(wǎng)格在R、Z、Theta三個維度上分別劃分了100，100，20份，覆蓋了整個模型空間，每個網(wǎng)格的權(quán)窗參數(shù)都初始化為0；
(2)在粒子輸運(yùn)計(jì)算過程中，對每個網(wǎng)格，統(tǒng)計(jì)進(jìn)入該網(wǎng)格的粒子數(shù)量和粒子權(quán)重，以及所有網(wǎng)格的平均粒子數(shù)量，在每個粒子輸運(yùn)過程中需要：
a)記錄該粒子的徑跡信息，即該粒子所穿過的所有網(wǎng)格的編號；
b)當(dāng)該粒子第一次進(jìn)入一個網(wǎng)格后，如果該粒子權(quán)重大于當(dāng)前權(quán)窗參數(shù)的最小值，該網(wǎng)格的粒子數(shù)量增加1，該網(wǎng)格的進(jìn)入權(quán)重增加粒子的權(quán)重；
c)將該粒子每一步輸運(yùn)中，最后進(jìn)入的網(wǎng)格的粒子數(shù)量和所有網(wǎng)格平均粒子數(shù)量做比較，如果網(wǎng)格粒子數(shù)平均粒子數(shù)，則該粒子徑跡上的所有網(wǎng)格的貢獻(xiàn)增加，該貢獻(xiàn)即為粒子進(jìn)入該網(wǎng)格對粒子密度均勻性的貢獻(xiàn)；
d)步驟c)中的對粒子密度均勻性的貢獻(xiàn)的初始值為1.0，并根據(jù)網(wǎng)格粒子密度權(quán)重和粒子輸運(yùn)時間權(quán)重進(jìn)行加權(quán)，其中粒子密度權(quán)重Fdensity定義為：
公式中Na為所有網(wǎng)格平均粒子數(shù)量，Ni為進(jìn)入網(wǎng)格i的粒子數(shù)量。
粒子輸運(yùn)時間則通過粒子當(dāng)前由于權(quán)窗而分裂的次數(shù)進(jìn)行近似，粒子輸運(yùn)時間權(quán)重Ftime定義為：
公式中Ts為粒子分裂的次數(shù)。
從而對粒子密度的均勻性的貢獻(xiàn)ΔCi定義為：
ΔCi＝Fdensity×Ftime
(3)計(jì)算權(quán)窗參數(shù)
a)根據(jù)公式
計(jì)算每個網(wǎng)格的重要性Impi。公式中Ci為所有計(jì)算的粒子進(jìn)入網(wǎng)格i對粒子密度均勻性的總貢獻(xiàn)，Wi為進(jìn)入網(wǎng)格i的粒子總權(quán)重；
b)將每個網(wǎng)格的重要性根據(jù)源所在網(wǎng)格的平均重要性進(jìn)行歸一化；
c)將每個網(wǎng)格的權(quán)窗參數(shù)設(shè)為歸一化后的網(wǎng)格重要性的倒數(shù)的一半；
(4)迭代優(yōu)化
a)在權(quán)窗生成過程中，共使用了3臺計(jì)算機(jī)(共60核)進(jìn)行輸運(yùn)模擬，每經(jīng)過5分鐘，則根據(jù)當(dāng)前的網(wǎng)格重要性生成權(quán)窗，并應(yīng)用于下一個周期的輸運(yùn)，同時每個周期生成的權(quán)窗都通過文本進(jìn)行輸出；
b)在每個周期輸運(yùn)開始時，將每個網(wǎng)格的貢獻(xiàn)和進(jìn)入權(quán)重都乘以一個折扣因子；
c)迭代自動進(jìn)行，在計(jì)算時間達(dá)到100分鐘后停止計(jì)算，并輸出最終的權(quán)窗。
(5)基于最終生成的權(quán)窗進(jìn)行輸運(yùn)計(jì)算，與直接模擬法相比，全局的計(jì)算效率提高了200余倍。
本發(fā)明根據(jù)一個粒子進(jìn)入網(wǎng)格后對粒子密度分布均勻性所做的貢獻(xiàn)計(jì)算得到每個網(wǎng)格的重要性，再依據(jù)網(wǎng)格重要性計(jì)算出權(quán)窗參數(shù)。同時，通過自適應(yīng)迭代計(jì)算逐步優(yōu)化權(quán)窗參數(shù)，將粒子均勻地輸運(yùn)到整個模型空間，從而降低全局計(jì)算結(jié)果的方差。本發(fā)明避免了傳統(tǒng)權(quán)窗生成方法僅能優(yōu)化模型空間單一位置方差的局限性，用戶不需要依賴經(jīng)驗(yàn)手動調(diào)試權(quán)窗，提高了復(fù)雜模型全局蒙卡計(jì)算的效率。
提供以上實(shí)施例僅僅是為了描述本發(fā)明的目的，而并非要限制本發(fā)明的范圍。本發(fā)明的范圍由所附權(quán)利要求限定。不脫離本發(fā)明的精神和原理而做出的各種等同替換和修改，均應(yīng)涵蓋在本發(fā)明的范圍之內(nèi)。