zxkjack123's blog
0%

On-the-flux GVR使用方法

Posted on Edited on In Views: Valine:
On-the-flux GVR使用方法

在蒙卡输运计算中，减方差方法对于屏蔽问题时非常重要的。这篇博客介绍了郑俞等开发的一种on-the-flux全局减方差技术[1]的使用方法。

编译

该方法修改了部分MCNP源代码，添加了一些源文件，需要重新编译。
为了不混淆可执行代码名称，下面我们用mcnp.mpi指代不带GVR的可执行文件，用mcnp_gvr.mpi指代GVR版代码。

使用

使用该方法，主要有以下几步：

  • 产生初始权窗
  • 使用GVR优化权窗
  • 使用优化后的权窗进行计算

下面以例题的形式来展示具体过程。我们用一个非常简单的几何，来针对某一个网格区域产生一个初始的(两群）权窗，然后用GVR优化这个权窗，最后使用优化后的全长进行计算。例题的输入文件如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
C example for on-the-flux GVR
C cell cards
1 1 -1.0 -1 imp:n=1
2 0      1 -2 imp:n=1
3 0      2 imp:n=0

C surface cards
1 so 10
2 so 20

C data crads
m1
      1001 -0.11
      8016 -0.89
sdef
nps 1000
C tally
f4:n 1
wwg   4  0  0  j j j j 0
wwge:n 0.1 20
mesh origin=-1 -1 -1 geom=xyz ref=0.0 0.0 0.0
     imesh=1 iints=10
     jmesh=1 jints=10
     kmesh=1 kints=10
C wwp:n 5.0 J 5 0 -0.5
print

产生初始权窗

首先用不带GVR的代码进行一次初始计算，产生一个初始化权窗。这个权窗是否精细准确不重要，只需要格式正确即可。

1
mcnp.mpi i=input

这个过程生成的wwout就是我们需要的初始权窗文件，将它重命名为wwinp1备用：

1
mv wwout wwinp1

使用GVR代码优化权窗

现在，我们使用带GVR的代码对初始权窗进行优化。

  • 将输入文件中的产生权窗的部分注释掉，使用wwp
  • 创建一个和权窗使用的网格对应的fmesh，注意网格和能群需要与权窗一致
  • 暂时只支持NPS设定history数目，不支持CTME为终止计算条件
  • wwp参数支持两种模式：
    • 权窗迭代计算模式：wwp第5个参数设置为-0.5，计算时会对权窗进行迭代优化，同时开启长历史控制；
    • 长历史控制模式：wwp第5个参数设置为非-0.5的值，如-1，计算将使用wwinp读取的权窗，不会进行权窗迭代。同时开启长历史控制。

改造后的输入文件如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
C example for on-the-flux GVR
C cell cards
1 1 -1.0 -1 imp:n=1
2 0      1 -2 imp:n=1
3 0      2 imp:n=0

C surface cards
1 so 10
2 so 20

C data cards
m1
      1001 -0.11
      8016 -0.89
sdef
nps 1000
C meshtally for GVR WW
fmesh4:n origin=-1 -1 -1 geom=xyz
         imesh=1 iints=10
         jmesh=1 jints=10
         kmesh=1 kints=10
         emesh 0.1 20
C meshtally for flux
fmesh14:n origin=-1 -1 -1 geom=xyz
         imesh=1 iints=20
         jmesh=1 jints=20
         kmesh=1 kints=20
         emesh 0.1 20
C wwg   4  0  0  j j j j 0
C wwge:n 0.1 20
C mesh origin=-1 -1 -1 geom=xyz ref=0.0 0.0 0.0
C      imesh=1 iints=10
C      jmesh=1 jints=10
C      kmesh=1 kints=10
C 
C 5th parameter use -0.5 if you want to update the WW, otherwise, use a negtive value other than -0.5
C 3rd reverse parameter: number of scored mesh element / total mesh element number
C 2nd reverse parameter: number of mesh element with rel_err < 0.1 / number of scored mesh element
C 1st reverse parameter: mesh coverage / valid geometry
wwp:n 5.0 J 5 0 -0.5 j j j 0.95 0.95 0.3
C use explicite capture when GVR is on
cut:n j 1e-11 0
print

运行代GVR可执行文件：

1
mcnp_gvr.mpi i=input wwinp=wwinp1

运行结束，将会生成weightwindow文件。

这里的weightwindow里面包含的数据，就是我们需要的权窗数据（但不包含权窗头文件）。我们将wwinp1的权窗头文件和weightwindow里的权窗数据组合起来，就是优化后的权窗文件了，命名它为wwinp2。

输出信息

权窗参数复现

在程序进行粒子输运前，会输出权窗信息，如下所示为一个使用np耦合输运的权窗信息：

1
2
3
4
5
newtronWWP 5.00000E+00 3.00000E+00 5.00000E+00 1.00000E+00-5.00000E-01 0.00000E+00
newtronWWP 0.00000E+00-1.00000E+00 0.00000E+00 1.00000E+00 1.00000E+33
PHOTONWWP 5.00000E+00 3.00000E+00 5.00000E+00 1.00000E+00-5.00000E-01 0.00000E+00
PHOTONWWP 0.00000E+00-1.00000E+00 0.00000E+00 1.00000E+00 1.00000E+33
stop WW updating parameters 3.00000E-01 1.00000E-01 3.00000E-01

这6行信息分别是中子、光子的权窗信息，及终止权窗迭代的信息。

权窗迭代现状

使用GVR时，在nps输出信息后面，还会跟上一些参数，表明当前权窗的迭代现状信息，以下列2(n_group=2, p_group=2)群的权窗信息为例：

1
2
3
 nps =       4000          1.84 minutes       445410        7185        5461        4803        3003           0           0
 nps =       4000          1.84 minutes       445410        9581       12967        3006        6156           0           0
stop WW updating parameters 1.33623E+05 3.00000E-01 1.00000E-01 5.37707E-02 0.00000E+00

其中第1行信息的含义为：

- 第一个参数：总网格数目，数据量1。 此处为445410，即网格共445410个网格单元。
- 第二组参数：有计数的网格数目, 数据量n_group。此处为7185 5461.即当nps=4000时，中子计数的第1群有7185个网格单元有计数，第2群有5461个网格有计数。
- 第三组参数：统计误差小于等于0.7的网格数目。此处为4803 3003，即中子计数的第1群有4803个网格统计误差小于0.7，中子计数的第2群有3003个网格统计误差小于0.7。
- 第四组参数：统计误差小于0.1的网格数目。此处为0 0，即中子计数的第1、2群都没有网格计数误差小于0.1。

第二行信息的含义与第一行一样，只是粒子为光子。

第三行信息的含义为：

1
2
nertron split      168025      922530     1090555                 79979   1.00000E+00 7.57956E+00 1.24783E+08 1.06568E+10 0.00000E+00
 photon split      639081       26141      665222                 79979   1.00000E+00 7.57956E+00

信息提示

1
=======///First clean for removing the badflux\\\\\====================

出现以下信息，表明权窗已经达到设定的条件，终止迭代。后续计算将使用当前权窗信息。

1
=======///From Now On weightwindow Is not update anymore\\\\\====================

常见错误

Recieved 16 but expected 32

在合并wwinp1和weightwindow中信息时，如果表头与内容不匹配，那么可能出现以下错误：

1
2
3
4
5
6
7
8
9
Fatal error in PMPI_Bcast: Other MPI error, error stack:
PMPI_Bcast(448)...............: MPI_Bcast(buf=0x55908bd959d0, count=4, dtype=0x4c000829, root=0, comm=MPI_COMM_WORLD) failed
PMPI_Bcast(434)...............:
MPIR_Bcast_impl(310)..........:
MPIR_Bcast_intra_auto(189)....:
MPIR_Bcast_intra_smp(105).....:
MPIR_Bcast_impl(310)..........:
MPIR_Bcast_intra_auto(223)....:
MPIR_Bcast_intra_binomial(123): message sizes do not match across processes in the collective routine: Received 16 but expected 32

出现这个错误的原因，可能是使用mcnp初始生成权窗生成时，运行模式为n, p。此时虽然没有要求生成光子权窗，但GVR会自动生成默认光子权窗(1群，能量上限100MeV,网格与中子网格相同)，此时wwout表头第一行第三个变量为2,第二行第二个变量为1。
但是在迭代优化时，如果没有明确要求迭代光子权窗，GVR不会自动生成迭代的光子权窗，此时生成的weightwindow只包含中子权窗信息，不会有光子权窗信息。
因此，在将weightwindow中的内容拷贝进入wwinp2时，表头信息和内容信息不匹配，表头中说有光子信息，实际内容没有，导致程序错误的想传递没有的信息，出现错误。
解决方法：修改表头使其和内容一致。

其他

本博客仅用于以例题说明该方法的使用，对该方法的原理，实现过程及其中涉及的其他文件的作用不做任何说明。使用方法最终解释权归代码原作者所有。本博客如有错误，敬请指正。

[1] Zheng, Yu, Yuefeng Qiu, Peng Lu, Yixue Chen, Ulrich Fischer, and Songlin Liu. “An Improved On-the-Fly Global Variance Reduction Technique by Automatically Updating Weight Window Values for Monte Carlo Shielding Calculation.” Fusion Engineering and Design 147 (October 1, 2019): 111238. https://doi.org/10.1016/j.fusengdes.2019.06.011.