非定常输运问题的数值求解对于众多领域的研究都是必需的手段, 例如惯性约束聚变(ICF)^[1]。蒙特卡罗(MC)方法由于其强大的几何及物理过程建模能力、精细的截面参数和高度的并行可扩展性等特征，日益成为核领域研究者的必备工具。由于非定常输运问题的时间特性，MC方法对其的模拟不可避免地以多计算步的方式进行。显然，上一时间步的未死亡粒子将在下一时间步以源粒子的方式出现，这些粒子在本文中被称为幸存粒子。在每一计算步的结束时刻，显然存在着一个幸存粒子的多维属性分布，虽然其为未知的，需要通过离散幸存粒子的各种属性(如空间坐标、能量等)作为代表。本文将此类分布称为幸存粒子分布。

由于不可能事先知道对一个具体问题多大的样本量是足够的，通常的做法是凭经验事先设定一个充分大的每步样本数并对结果做后续检验。显然，对有些时间步而言，这个事先设定的样本数可能太大，以致于结果已经收敛而仍有不必要的计算。从逻辑上来说，如果可以基于某种一般性的准则自动调整每一计算步的样本数，则整个计算可以获得更高的效率。幸运的是，对于非定常输运问题而言，天然存在一个基于香农熵的自动调整样本数策略。

将信息论中香农熵的概念引入粒子输运问题的蒙特卡罗模拟首要的是作为临界计算中对裂变源分布的收敛性诊断^[2-4]，而临界计算属于定常输运问题。在临界计算的每一迭代步，都存在一个对应于裂变源分布的香农熵值，该值序列的收敛是裂变源分布收敛的必要条件。一般来说，香农熵值序列的收敛性诊断是一个后验诊断，用以在整个计算完成后判断所得裂变源分布是否已经收敛，虽然也有一些关于在线诊断收敛性的研究^{[2-3, 5]}。显然，在临界计算的裂变源粒子和幸存粒子之间存在某种相似性，它们都是在上一步产生并在下一步作为源粒子。毫无疑问，香农熵也可以作为幸存粒子分布的收敛性诊断工具并发展出一套基于此的自动调整样本数策略。这一策略分为三个步骤：首先，给每一计算步设定一个充分大的样本数，将该样本数分成若干批并逐批计算；其次，每批模拟完成后计算一个新的熵值，由此对每一计算步得到一个熵值序列；最后，如果基于某种香农熵的在线诊断方法表明熵值序列已经收敛，则表明幸存粒子分布已经收敛，可以提前结束本时间步的计算。

从上述策略的描述可以看出，有两个主要的技术问题需要解决：第一个问题是需要一种消息传递并行编程环境(MPI环境)下的高效香农熵计算方法。因为此策略需要在每一计算步频繁计算香农熵且MPI环境是重要的、我们选用的并行编程环境。如果频繁计算香农熵的时间代价太大，减少样本数的好处将被部分或全部抵消；第二个问题是需要发展一种熵值序列的在线诊断方法。由于不可避免的统计特性，这一方法需要得到数值试验的检验。本文描述了一种MPI环境下的高效香农熵计算方法, 基于金融工程中的随机振荡指标^[6]提出了一种熵值序列的在线诊断方法以实现自动调整样本数。对于一个空间一维非定常输运模型得到了一些积极的数值结果, 验证了该策略的有效性。

1.   MPI并行环境下香农熵的高效计算方法

裂变源分布对应的香农熵^[2]

式中：N是计算香农熵所需MESH网格的总数；M_h是裂变源粒子的总数；M_i是所有裂变源粒子中，其空间坐标位于第i个MESH网格内的所有粒子的总数；P_i是所有裂变源粒子中，其空间坐标位于第i个MESH网格内的所有粒子所占的份额。如果将这一概念延伸至幸存粒子分布所对应的香农熵，则P_i显然应该是所有幸存粒子中，其属性位于第i个相空间MESH网格内的所有粒子所占的份额。这里，考虑的粒子属性可以是空间、能量和方向等六个属性中的若干种，常见的只考虑空间和能量这四种属性。从上述定义可以看出，在MPI并行环境下，所有进程应该将它们自己的M_i进行归约求和才能正确地计算熵值。如果MESH网格的数目较大，频繁归约的时间代价将会很大，而对于许多实际问题，MESH网格的数目经常是较大的。

为了在MPI并行环境下高效计算香农熵值，我们在文献[7]中提出了一种新的方法。简言之，新方法使每个进程基于它们自己的局部数据计算自己独立的熵值, 再将所有熵值规约求和后以总进程数进行平均, 得到的平均值就是最终的熵值。由于频繁规约的只是一个数，所以新方法计算香农熵的时间将大为减少。可以证明，在进程数有限，总样本数趋于无穷的前提下，新方法计算的熵值将趋近于原方法计算所得的熵值。

2.   基于随机振荡指标的自动调整样本数策略

基于金融工程中常见的随机振荡指标，已经提出了一种评价蒙特卡罗临界计算中裂变源分布收敛性的后验诊断方法^[6]。在此方法中，第n迭代步的随机振荡指标可以表示为

式中：H_min^{n, p}和H_max^{n, p}是第n迭代步之前p个迭代步所对应熵值序列的最小值和最大值；Hⁿ是当前第n迭代步的熵值。如果裂变源分布已经收敛，Hⁿ应该在期望值附近随机振动，从而Kⁿ应该在0.5附近随机振动。所给出的评价裂变源分布收敛性的后验准则是看Kⁿ是否满足下列不等式

即，如果Kⁿ和它在之后m步内的平均值在0.5的上下ε范围内，则认为熵值序列已经收敛。p, m, ε是该后验准则的参数，按照文献[6]的建议，它们分别取值为20, 50, 0.1。如果临界计算的第一个激活迭代步的序号大于第一个满足不等式(4)的迭代步的序号，则整个临界计算的结果被认为可以接受。

经过修改，上述后验收敛性判断准则可以作为幸存粒子分布的在线收敛性判断准则。即，如果将总样本数划分为M批, 则可以得到熵值和随机振荡指标的序列值Hⁱ(i=1, …, M)和K^j(j=p, …, M)。当完成了第n批次的计算(n≥m+p－1)且下述不等式得到满足

则当前步的计算可以提前终止，因为幸存粒子的分布可以认为已经收敛。

3.   数值结果

通过在Godiva模型的外部增加一个²³⁸U金属球壳，这个一维非定常输运模型可以用来验证本文提出策略的有效性。在初始时刻，内球被分为6个同心球壳，外壳被分成10个同心球壳。这些球壳的密度经过调整，使得整个系统初始时处于超临界状态。在输运计算的同时系统的力学状况也将发生变化。由于开始时中子数将急剧增加，为保证计算的顺利进行，将利用舍选法挑选固定数目的总样本数，通过调整源粒子的权重保持计算的无偏性。

对这个模型的蒙特卡罗计算，初始样本数被设定为每步5×10⁶并利用50个CPU核进行并行计算。每批次的样本数被设定为2000。用来计算香农熵的相空间MESH网格数为160，其中一维空间被划分为16个网格，能量空间被划分为10个网格。整个系统所释放的能量和包含的中子数被用来作为考察计算结果是否变化剧烈的参考标准。按照以前的物理分析及多个程序的考核，可以认为每步5×10⁶样本已经足以使计算结果收敛，希望在自动调整样本数的情况下，计算结果与不调整的结果的偏差尽可能小。

图  1  调整前后的样本数比较

Figure  1.  History number before and after adjustment

下载: 全尺寸图片 幻灯片

如果以不调整样本时的计算时间作为1个单位，则调整后的计算时间只有原来的24%左右。计算时间的节省也可以从图 1看出。图 2和图 3为调整前后，总放能和总中子数随时间的变化情况，可以看出二者非常接近，最终的结果相差不到0.8%。从上述结果可以看出，基于香农熵的自动调整样本数策略实现了在数值结果基本不变的前提下计算时间大大节省的目的。

图  2  调整样本数前后总放能随时间的变化

Figure  2.  Total released energy before and after adjusting history number

下载: 全尺寸图片 幻灯片

图  3  调整样本数前后，总中子数的变化

Figure  3.  Neutron number in the whole system before and after adjustment

下载: 全尺寸图片 幻灯片

上述模型显然不是一个慢收敛问题。如果应用此策略到一个慢收敛模型上，则需要更多的样本才能满足不等式(5)。当原先设定的每步总样本数不是特别大时，自动调整样本数后的时间减少将不会那么可观，但是，如果一开始就设定了一个非常巨大的初始样本数，本策略还是预期可以减少相当数量的计算时间。

4.   结论

基于对信息论中香农熵概念和金融工程中随机振荡指标的借鉴引用，本文提出了一种在线香农熵序列收敛性诊断方法并在此基础上实现了非定常输运问题的自动调整样本数功能。对一个空间一维非定常输运模型的计算表明，该策略可以在保持计算结果基本不变的前提下实现计算时间的减少，从而提高了计算效率。该策略计划将通过更多的模型计算进行进一步的优化和分析。