zxkjack123's blog
0%

使用PyNE生成MAGIC权窗的例子

Posted on Edited on In Views: Valine:
None

说明

本篇文章使用例子来说明使用PyNE来进行MAGIC的权窗文件的流程。

本流程需要提前安装PyNE (PyNE的安装参考这里)

例子描述

这里直接使用用于R2S计算停机剂量的例子。例题几何与材料描述参考这里。

DAGMC input文件准备

input文件由于仅用于迭代权窗，而不是为了计算精确的通量数据。因此与普通输入文件有几点不同：

  1. 使用的mesh tally是为了生成权窗，网格voxel的尺寸不要超过材料的平均自由程(mean free path)；
  2. mesh tally的能群一般设置为两群就够了(如 0.1 20)；
  3. 迭代过程中一般使用较高的能量截断(energy cut-off)，随着迭代可以逐步降低截断能量，逐步获得更精确的数据。

对于此例子的一个简单input文件如下：

1
2
3
4
5
6
7
8
9
10
11
12
C sample input for MAGIC weight window
C neutron source: point source, 14 MeV
sdef pos=0 0 0  erg=14  
C 
nps 1e4
cut:n j 0.01 -0.5 -0.25 1
C Tally card
fmesh4:n geom=xyz origin= -5 -5 -5
         imesh=15 iints=2
         jmesh=5 jints=1
         kmesh=15 kints=2
         emesh 0.1 20

第一次中子输运计算

使用DAGMC进行输运计算时，需要通过i=**来指定输入文件，通过g=**来指定加载了材料的几何文件，详细运行指令可以参考官方文档。

对于本例题，可以使用下面的方法运行本例子。
单线程运行：

1
2
# 需要先保障系统能够找到正确的mcnp5.mpi，找不到的话需要给出完整路径
mcnp5.mpi i=input g=geom.h5m

多线程运行：

1
mpirun -np n mcnp5.mpi i=input g=geom.h5m

备注：input文件和几何文件的名称长度都不能超过8个字符，可以使用ln -sf original_long_name.h5m geom.h5m把原文件链接到一个较短名字的文件后再运行。

计算完成之后，我们便得到了一个meshtal文件。

生成权窗

MAGIC生成权窗的说明可以参考这里。实际使用时，我们通过PyNE中的工具读取meshtal文件中的内容，生成MAGIC权窗，然后写入到文件中。下面是一个生成权窗的python脚本例子：

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
import sys
from pyne import variancereduction as vr
from pyne.mcnp import Meshtal, Wwinp

if __name__ == "__main__":

    f = sys.argv[1] #must be a meshtal file
    M = Meshtal(f)
    # User can also define the tally as 4, 14, ...
    tally_number = M.tally.keys()[0]
    meshtally = M.tally[tally_number]
    tags = meshtally.tag_names

    res_tag = tags[0]
    err_tag = tags[1]

    vr.magic(meshtally, res_tag, err_tag)

    meshtally.write_hdf5("mesh-ww-tags.h5m")

    wwinp = Wwinp()
    wwinp.read_mesh(meshtally.mesh)

    # write out the wwinp file
    wwinp.write_wwinp("wwinp-magic")

bash执行下面命令：

1
python magic.py meshtal

之后便生成了我们所需要的权窗文件wwinp-magic。

使用权窗

使用权窗时，使用权窗时，在输入文件中添加合适的WWP参数读取生成的权窗文件中的权窗数据即可。下面是第二轮计算的输入文件例子：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
C sample input for MAGIC 2nd run
C neutron source: point source, 14 MeV
sdef pos=0 0 0  erg=14
C 
nps 1e4
C weight window parameters
wwp:n 5 j 5 0 -0.5
cut:n j 1e-3 -0.5 -0.25 1
C Tally card
fmesh4:n geom=xyz origin= -5 -5 -5
         imesh=15 iints=2
         jmesh=5 jints=1
         kmesh=15 kints=2
         emesh 0.1 20

备注： wwp和cut有变化。

使用权窗计算后，会有新一轮的结果。可以通过对比meshtal的变化发现，使用权窗后的结果误差更小(这里的例子规模太小，提升效果并不明显)。

最终计算

当迭代出满意的权窗后，就可以使用那个权窗进行真正的通量计算了。那个时候再根据实际情况使用恰当的cut卡，使用真正的tally卡，计算得到想要的结果。

备注：在计算得到结果之后，还可以把几何图形及计算结果的网格数据进行可视化绘图。绘图方法可以参考DAGMC数据绘图。