一种基于粒子密度不均匀性的迭代蒙卡全局权窗参数生成方法与流程

本发明涉及一种基于粒子密度不均匀性的迭代蒙卡全局权窗参数生成方法，属于核物理、核技术应用等领域中的辐射屏蔽计算分析方向。

背景技术：

越来越多的蒙特卡罗辐射输运计算需要得到全局收敛的计算结果，比如裂变堆的全堆功率分布计算和聚变堆的停堆剂量计算。但对于厚屏蔽和几何复杂的问题，粒子难以穿透屏蔽区域，若使用直接模拟法进行计算，不仅屏蔽区域的统计误差会远高于其他区域，计算耗时也难以接受。

使用网格权窗是一种常用的蒙特卡罗减方差方法，权窗是相空间(空间-能量)的分裂和俄罗斯轮盘赌技巧。根据粒子权重和粒子当前所在网格的权窗参数，低权重粒子通过俄罗斯轮盘赌淘汰，高权重粒子分裂。这项技术可以根据用户定义的相空间重要性调整蒙特卡罗粒子分布，同时能够控制粒子权重波动，保证粒子权重分布在一个合理的界限内。

手动进行权窗参数设置需要大量的使用经验，同时难以保证权窗的加速效果。而传统的权窗产生方法虽然在一定程度上辅助生成权窗参数，但存在若干缺陷：1)只能优化相空间单一位置的计算效率；2)需要用户手动迭代优化；3)生成的权窗偏向于将粒子沿着屏蔽最小的路径输运；4)难以优化厚屏蔽区域的计算效率。

低质量的权窗参数不仅不能优化计算效率，还可能导致某些区域采样不足，而得到错误的计算结果。因而如何自适应地生成高效的权窗参数，是基于蒙特卡罗的复杂模型输运计算的关键技术之一。

目前国际上主流的自动权窗生成方法主要有两种，一种为FW-CADIS(Forward Weighted-Consistent Adjoint Driven Importance Sampling)方法，通过确定论程序进行伴随计算得出每个栅元的重要性，再根据栅元的重要性生成权窗来加速计算。FW-CADIS方法的主要不足为蒙特卡罗输运计算程序的用户拥有需要高效的确定论输运计算程序，并且熟悉确定论程序的使用方法，才能根据FW-CADIS方法进行输运加速。另外一种方法为MAGIC(Method of Automatic Generation of Importance by Calculation)方法，通过迭代计算正向通量，用正向通量的倒数作为权窗参数加速输运计算。其主要不足在于权窗参数生成时间相对较长，且适用性没有FW-CADIS方法广泛。

技术实现要素：

本发明技术解决问题：克服传统的需依靠用户经验进行优化，同时无法对全局的计算效率做优化的权窗产生器的不足，提供一种基于粒子密度不均匀性的迭代蒙卡全局权窗参数生成方法，用户不需要依赖经验手动调试权窗，提高了复杂模型全局蒙卡计算的效率。

本发明的技术方案如下：一种基于粒子密度不均匀性的迭代蒙卡全局权窗参数生成方法，包括以下步骤：

(1)生成初始权窗网格：

a)建立覆盖整个模型空间的网格；

b)将所有网格的权窗参数都初始化为0；

(2)在粒子输运计算过程中，统计进入每个网格的粒子数量(N_i，i指第i个网格)，粒子权重(W_i)，和每个网格对粒子密度均匀性的贡献，以及所有网格的平均粒子数量(Na)，在每个粒子输运过程中需要：

a)记录粒子的径迹信息，即粒子所穿过的所有网格的编号；

b)当粒子第一次进入网格i后，如果粒子权重w大于当前所有网格权窗参数的最小值，将进入网格i的粒子数量N_i增加1，进入网格i的总粒子权重W_i增加当前粒子的权重w；

c)将粒子每一步输运中，终点所在的网格的粒子数量N_i和所有网格的平均粒子数量Na做比较，如果N_i小于Na，则该粒子径迹上的所有网格的贡献增加，增加的贡献即为粒子进入一个网格对粒子密度均匀性的贡献；

d)步骤c)中的对粒子密度均匀性的贡献的初始值为1.0，并根据网格粒子密度权重和粒子输运时间权重进行加权，粒子密度权重F_density定义为：

$F_{density}=\frac{N_a}{N_i}$

公式中N_a为所有网格平均粒子数量，N_i为进入网格i的粒子数量。

粒子输运时间则通过粒子当前由于权窗而分裂的次数进行近似，粒子输运时间权重F_time定义为：

$F_{time}=\frac{1.0}{1.0+T_s}$

公式中T_s为粒子分裂的次数。

从而对粒子密度的均匀性的贡献定义为：

ΔC_i＝F_density×F_time

(3)计算权窗参数

a)根据公式

${\mathrm{Imp}}_i=\frac{C_i}{W_i}$

计算每个网格的重要性Imp_i。公式中C_i为所有计算的粒子进入网格i对粒子密度均匀性的总贡献，W_i为进入网格i的粒子总权重；

b)将每个网格的重要性根据源所在网格的平均重要性进行归一化；

c)将每个网格的权窗参数设为归一化后的网格重要性的倒数的一半；

(4)迭代优化

a)在权窗生成过程中，将输运分为若干个周期，根据前一个周期的网格重要性生成新的权窗参数，并应用于当前周期的输运，同时每个周期生成的权窗参数都通过文本进行输出；

b)在每个周期输运开始时，将每个网格的贡献和进入权重都乘以一个折扣因子f(取0～1之间的值)；

c)迭代自动进行，在达到了用户设定的时间上限后停止计算，当前的权窗参数即为最终生成的全局权窗参数；

本发明与现有技术具有如下优点：

(1)本发明自适应地生成全局权窗参数，用户只需提供一个覆盖整个模型的网格，简化了用户的操作，而且避免了用户由于经验缺乏而生成的低质量权窗参数可能导致的错误计算结果；

(2)本发明生成的权窗能极大的提升蒙特卡罗粒子输运在模型全空间的收敛速度，实现对厚屏蔽和大空间等模型的准确快速粒子输运。

附图说明

图1为本发明的实现流程图。

具体实施方式

国际热核聚变实验堆ITER是目前全球规模最大的国际科研合作项目之一，ITER装置是一个能产生大规模核聚变反应的超导托卡马克。本专利选用ITER官方发布的中子学模型Alite.004作为基准例题。测试用的模型是一个40度的Alite模型，半径约20米，高约40米。其放射源为完整的聚变等离子体源。

针对该例题，如图1所示，本发明的实施方式为：

(1)生成初始权窗网格，网格在R、Z、Theta三个维度上分别划分了100，100，20份，覆盖了整个模型空间，每个网格的权窗参数都初始化为0；

(2)在粒子输运计算过程中，对每个网格，统计进入该网格的粒子数量和粒子权重，以及所有网格的平均粒子数量，在每个粒子输运过程中需要：

a)记录该粒子的径迹信息，即该粒子所穿过的所有网格的编号；

b)当该粒子第一次进入一个网格后，如果该粒子权重大于当前权窗参数的最小值，该网格的粒子数量增加1，该网格的进入权重增加粒子的权重；

c)将该粒子每一步输运中，最后进入的网格的粒子数量和所有网格平均粒子数量做比较，如果网格粒子数平均粒子数，则该粒子径迹上的所有网格的贡献增加，该贡献即为粒子进入该网格对粒子密度均匀性的贡献；

d)步骤c)中的对粒子密度均匀性的贡献的初始值为1.0，并根据网格粒子密度权重和粒子输运时间权重进行加权，其中粒子密度权重F_density定义为：

$F_{density}=\frac{N_a}{N_i}$

公式中N_a为所有网格平均粒子数量，N_i为进入网格i的粒子数量。

粒子输运时间则通过粒子当前由于权窗而分裂的次数进行近似，粒子输运时间权重F_time定义为：

$F_{time}=\frac{1.0}{1.0+T_s}$

公式中T_s为粒子分裂的次数。

从而对粒子密度的均匀性的贡献ΔC_i定义为：

ΔC_i＝F_density×F_time

(3)计算权窗参数

a)根据公式

${\mathrm{Imp}}_i=\frac{C_i}{W_i}$

计算每个网格的重要性Imp_i。公式中C_i为所有计算的粒子进入网格i对粒子密度均匀性的总贡献，W_i为进入网格i的粒子总权重；

b)将每个网格的重要性根据源所在网格的平均重要性进行归一化；

c)将每个网格的权窗参数设为归一化后的网格重要性的倒数的一半；

(4)迭代优化

a)在权窗生成过程中，共使用了3台计算机(共60核)进行输运模拟，每经过5分钟，则根据当前的网格重要性生成权窗，并应用于下一个周期的输运，同时每个周期生成的权窗都通过文本进行输出；

b)在每个周期输运开始时，将每个网格的贡献和进入权重都乘以一个折扣因子；

c)迭代自动进行，在计算时间达到100分钟后停止计算，并输出最终的权窗。

(5)基于最终生成的权窗进行输运计算，与直接模拟法相比，全局的计算效率提高了200余倍。

本发明根据一个粒子进入网格后对粒子密度分布均匀性所做的贡献计算得到每个网格的重要性，再依据网格重要性计算出权窗参数。同时，通过自适应迭代计算逐步优化权窗参数，将粒子均匀地输运到整个模型空间，从而降低全局计算结果的方差。本发明避免了传统权窗生成方法仅能优化模型空间单一位置方差的局限性，用户不需要依赖经验手动调试权窗，提高了复杂模型全局蒙卡计算的效率。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。