mcnp权窗生成器
简介
Mcnp自带有降低方差技术，很多其他的蒙卡代码因为没有这项技术，使得计算问题无法足够快速的完成。人们趋向于学习与他们手边工作相关的降低方差技术。MCNP使用说明书中叙述了降低方差技术在使用过程中的细节，但是使用者经常会犯愁如何将这些降低方差技巧用于他们的实际问题中。对于MCNP的使用者，《A Sample Problem for Variance Reduction in MCNP》是一本值得一读的报告。事实上，当前的报告成功的提供了一些推动，他们通过例子涉及了降低方差技术的教学，至少是部分。
人们希望这些报告中的例子能够对MCNP的使用者有帮助，但是例子中只有很少部分设计方差降低技术。对于那些研究数学的的使用者，标准的书籍通过积分方程讨论的降低方差理论很用用处。读这些书的人会警告大家，这些书的作者也发现这些书中的降低方差的理论同使用大多数的MCNP计算过程没有直接的关系。尽管如此，这些书还是能够介绍一些基础的理论。
一些例子对用户在课程或MCNP论坛的问题给出了直接的回应。其他的例子有更多风格的。
权窗发生器的输出全是0
最常见的一个挫折是使用者把权窗发生器仅仅当做一个黑箱来使用而不去理解他怎么工作。发生器是一个相当简单的对例子在不同相空间的平均分数（重要性）统计估计。对于那些感兴趣的，在第43和44页给出了细节。就当前目的而言，使用者应当明白，发生器输出的是粒子进入特定相空间的平均数的倒数。输出0表示在给定相空间中没有粒子（或者次级粒子）计数。有两个原因使得在区域内没有粒子计数：
1、没有粒子进入这一区域，所以没有计数；
2、粒子计入这一区域，但没有一个被计数。
例如：
C 200 cm concrete slab problem
1 1 -2.03 -
2 1
2 1 -2.0
3 -3 2
3 1 -2.03 -
4 3
4 1 -2.03 -
5 4
5 1 -2.03 -
6 5
6 1 -2.03 -
7 6
7 1 -2.03 -8 7
8 1 -2.03 -9 8
9 1 -2.03 -10 9
10 1 -2.03 -11 10
11 1 -2.03 -12 11
12 1 -2.03 -13 12
13 1 -2.03 -14 13
14 1 -2.03 -15 14
15 1 -2.03 -16 15
16 1 -2.03 -17 16
17 1 -2.03 -18 17
18 1 -2.03 -19 18
19 1 -2.03 -20 19
20 1 -2.03 -21 20
21 0 -1
22 0 21
1 py 0
2 py 10
3 py 20
4 py 30
5 py 40
6 py 50
7 py 60
8 py 70
9 py 80
10 py 90
11 py 100
12 py 110
13 py 120
14 py 130
15 py 140
16 py 150
17 py 160
18 py 170
19 py 180
20 py 190
21 py 200
4
mode n
imp:n 1 19r 0 0
c the following is schaeffer portlan
d concret
e m1
1001.50c -.010
8016.50c -.529
11023.51c -.016
12000.51c -.002
13027.50c -.034
14000.51c -.337
19000.51c -.013
20000.51c -.044
26000.55c -.014
6012.50c -.001
sdef x=0 y=1.e-6 z=0 cel=1 wgt=1 erg=14
c generate spatial windows
wwg 1 1 .5
f1:n 21
cut:n 1.e20 .01 0 0
ctme 5
nps 100000
prdmp j -60 j 2
使用这个例子运行100000个粒子，生成如下的权窗值：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
-1.0000E+00 -1.0000E+00
这是一个正常的结果，特别是对那些初学者而言。粒子径迹输入126表格表明了问题所在：
从表格中可以看出，没有粒子进入18栅元，因此没有粒子穿过计数器的21面，所以没有计数。
这里有若干种方法使用者可以尝试解决这一问题。显而易见的一种是增加运行粒子数目。要增加多少呢？这个有点难说，但是通过126表格中的径迹衰减数据而言，10倍的粒子数应当足够了。但是这一判断在使用过程中并不对。
除了暴力的增加运行粒子总数外，还有别的方法来解决权窗值为0的情况。一些方法是：
1、分配问题中的源的特性（Distribute the source），至少可以使得部分栅元获得较好的重要性评价。再使用这个重要性参数（权窗）来运行另一个问题以便获得更好的评价值。
这种方法事实上使得重要性方程不依赖于源，而取决于计数器。
2、保持源不变，降低密度使得一些例子可以被计数。再使用权重窗去运行实际密度的问题。有时人们不得不使用成对的密度阶来一步步运行从低密度到高密度的结果。
3、保持源不变，使用指数传输方法来减小有效截面。
4、根据观察模拟计算过程，对权窗或者重要性参数做一个初始的假设。这经常是最好的方法来解决权窗值或者栅元重要性参数为一个很小的数。如果在一个问题中有1000个栅元，这可能变得冗长而费时。
5、使用另一个同最终计数器高度相关，但更易于观察的计数器。用计数器来用于其他计数器的窗口。那些用于其他计数器的窗口，尽管没有对目标窗口进行优化，但总比什么信息都没有好很多。用于其他计数器的窗口，它通常可以用于对最终窗口的窗口生成。
第一种方法生成的窗口：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 7.2866E-04 3.6386E-04 1.8151E-04
1.2053E-04 4.8785E-05
2.8172E-05 1.0236E-05 4.3976E-06
2.0247E-06 9.3715E-07 4.5544E-07 2.1243E-07 1.0000E-07
-1.0000E+00 -1.0000E+00
从slab0的运行结果可得到，粒子轻易的通过了超过半数的平板。因此，从栅元1到栅元8获得粒子计数应该不存在问题，因此权窗信息得以开始。唯一的一个值得注意的问题是所有进入栅元8的粒子拥有一个期望大约为0.001的权重，因为窗口的权重下限为
7.2866E-4。因此，在slab3运行中，slab2中获得的权窗被使用，且源的权值为0.001.
Slab3对每一个栅元都建立了权窗值
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 5.0000E-01 1.1527E-01 6.2555E-02 3.2895E-02 1.7015E-02
8.4250E-03 4.2100E-03 2.0900E-03 1.0050E-03 5.0000E-04
2.4250E-04 1.1773E-04 5.2577E-05 2.1667E-05 9.5579E-06
4.7931E-06 1.6319E-06 7.8125E-07 3.6458E-07 1.4757E-07
-1.0000E+00 -1.0000E+00
把这个窗值插入到输入文件中并进行一次迭代会对结果有所帮助，因为表格160表明这个窗口只基于6个计数，但这一示范也可使用，因为他对于现有问题足够了。
对于方案2，从slab0得知，粒子至少通过了一半的平板。这表明如果密度下降一个因子，2，一些粒子将被计数，同事权窗将生成。对于初始密度，生成的权窗并不正确，但对于实际密度而言，相对于全是0的权窗，至少是一个较好的降低方差的权窗。Slab4生成的权窗：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 5.0000E-01 1.5584E-01 1.1426E-01 8.4810E-02 6.1680E-02
4.4995E-02 3.2525E-02 2.3585E-02 1.6785E-02 1.1930E-02
8.4200E-03 5.8800E-03 4.0200E-03 2.7900E-03 1.9600E-03
1.3950E-03 1.0200E-03 7.1000E-04 4.3000E-04
2.8000E-04
-1.0000E+00 -1.0000E+00
将Slab4生成的权窗带回初始密度（Slab5），Slab5运行获得如下窗值：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 5.0000E-01 1.1528E-01 6.2738E-02 3.3215E-02 1.7151E-02
8.5562E-03 4.2231E-03 2.0497E-03 1.0004E-03 4.8758E-04
2.3900E-04 1.1243E-04 5.5200E-05 2.5313E-05 1.1992E-05
5.5469E-06 2.4414E-06 1.0547E-06 4.4922E-07 2.2461E-07
-1.0000E+00 -1.0000E+00
通过迭代可以改进窗值的结果，但他表明如何使用方法2。
方法3使用指数传递法来减小前进方向上的有效截面。它使得粒子在前进方向上优先跳过一段距离。明白这一点后，这一方法同方法2是不同的，因为方案2使用一个虚构的密度，而方案3使用真实的密度以及有区分的抽样方案。在无区分的问题（Slab0）中，粒子至少通过了一般的路程，因此，在前进方向上通过一个因子2，使得指数传递被设置为减弱一半有效截面（Slab6）。Slab6生成了一个警告信息“警告，指数传递通常需要权窗值（warning.exponential transform usually needs weight window）,”但这一信息将不予理会，只到获得可靠的结果。例如，如果修正值表明权窗值在使用因子2时不可靠，这并不是一个严重的问题。但如果最终评估值在使用因子2时不可靠，这是一个严重的问题。在建立权窗时，警告信息可以被忽略。Slab6生成权窗如下：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 5.0000E-01 1.1458E-01 6.2718E-02 3.3055E-02 1.6950E-02
8.5871E-03 4.1946E-03 2.0540E-03 1.0088E-03 4.9506E-04
2.3645E-04 1.1029E-04 4.7777E-05 2.4353E-05 9.8767E-06
5.1597E-06 2.3141E-06 8.7391E-07 4.3701E-07 2.9811E-07
-1.0000E+00 -1.0000E+00
权窗值可通过迭代方法得到改进。
方案4，通过观察slab0的径迹输入可以发现粒子径迹从一个栅元到另一个栅元时大约按照因子2递减。这点说明径迹可以一个相对的常数2设置窗口来解决这一问题，因此，没一个权重值可以依据栅元逐个以因子2递减。使用这一假设获得权窗（slab7），slab7生成的窗值如下：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 5.0000E-01 1.1592E-01 6.3000E-02 3.3363E-02 1.7097E-02
8.5240E-03 4.2772E-03 2.0897E-03 1.0277E-03 4.9595E-04
2.4282E-04 1.1738E-04 5.6387E-05 2.7182E-05 1.2855E-05
5.9984E-06 2.8597E-06 1.3302E-06
6.2955E-07 2.9390E-07
-1.0000E+00 -1.0000E+00
应当注意，生成的权窗值对假设的要求并不高，例如，slab8中栅元间的因子设定为0.707会生成如下窗值：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 5.0000E-01 1.1593E-01 6.2942E-02 3.3075E-02 1.7017E-02
8.3292E-03 4.1146E-03 2.0437E-03 9.8542E-04 4.8021E-04
2.2500E-04 1.0990E-04 5.5729E-05 2.5781E-05 1.1198E-05
5.8594E-06 2.5391E-06 1.6276E-06 8.1380E-07 4.8828E-07
-1.0000E+00 -1.0000E+00
Slab9将因子设定为0.333，生成窗值如下：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 5.0000E-01 1.1800E-01 6.2000E-02 3.1500E-02 1.5000E-02
8.4722E-03 3.9969E-03 1.9918E-03 1.0394E-03 5.2995E-04
2.6478E-04 1.3001E-04 6.4132E-05
3.1505E-05 1.5440E-05
7.4715E-06 3.5575E-06 1.6693E-06 7.8905E-07 3.6547E-07
-1.0000E+00 -1.0000E+00
重点在于使用者不需要一个很好的假设来开始进程。
方案5使用了这样一个事实，粒子从源到计数器首先必须经过从源到计数器的部分路程。在slab10中，生成器被设置为基于面计数器（f11：n，一个在100cm处的中间计数器）生成窗口，11号面后的区域的重要性暂时被设为0，slab10生了了如下窗口：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 5.0000E-01 1.1528E-01 6.2555E-02 3.2895E-02 1.7015E-02
8.4250E-03 4.2100E-03 2.0900E-03 1.0050E-03 5.0000E-04
-1.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
-1.0000E+00 -1.0000E+00
现在讲11栅元的-1.0的输入改为0，因为在接下来的过程中并不希望粒子在11栅元被终结。输入卡片为：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 5.0000E-01 1.1528E-01 6.2555E-02 3.2895E-02 1.7015E-02
8.4250E-03 4.2100E-03 2.0900E-03 1.0050E-03 5.0000E-04
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
-1.0000E+00 -1.0000E+00
这一窗口由slab10得到，然后再带入到slab11中，使用初始的探测器（f1：n 21）。slab11将生成如下窗值：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 5.0000E-01 1.1621E-01 6.3145E-02 3.3654E-02 1.7062E-02
8.5618E-03 4.3042E-03 2.1257E-03 1.0423E-03 5.0584E-04
2.4569E-04 1.1841E-04 5.5736E-05 2.5203E-05 1.2232E-05
5.7225E-06 3.0158E-06 1.3830E-06 5.6641E-07 2.3438E-07
-1.0000E+00 -1.0000E+00
值得注意的是由slab10和slab11生成的栅元10的窗值非常相似。因为他们相当高的相关度，在slab10中的中间的数器（f11：n 11）事实上扮演了所需探测器（f1：n 21）的替身。
Slab12回到了这个问题上，问什么slab10对栅元11的重要性设置为0。特别是slab12同slab10除了栅元11的重要性不同外，其余的都相同。slab12生成的窗值如下：
wwp:n 5 3 5 0 0 0
wwe:n 1.0000E+02
wwn1:n 5.0000E-01 1.1528E-01 6.2565E-02 3.2895E-02 1.7015E-02
8.4250E-03 4.2100E-03 2.0900E-03 1.0100E-03 5.0000E-04
8.4632E-04 1.5410E-03 1.4517E-03 0.0000E+00 0.0000E+00
0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00 0.0000E+00
-1.0000E+00 -1.0000E+00
值得注意的是，超过栅元10的窗值比栅元10的值更高。这是因为，如果这一窗值被用做（替代slab11窗值的）输入卡片，那么粒子在从面11（中间面）发出到面21（最终面）时，实际上会发生轮盘赌。这是可以理解的，因为为中间计数器生成的结果将会自然的拥有一个窗值的增长，伴随着粒子相更远的地方走，从面11向两个方向传播。