zxkjack123's blog
0%

使用PyNE进行基于网格的计算停机计量的例子

Posted on Edited on In Views: Valine:
None

说明

本篇文章使用例子来说明使用PyNE来进行基于网格的计算停机剂量率的流程。

本流程需要提前安装Trelis, DAGMC, ALARA和PyNE。
DAGMC的安装可以参考这里， PyNE的安装可以参考这里。ALARA的安装可以参考这里。

例子的重要性

使用的工具功能越多，那么使用起来一般就越复杂，新手就越难上手，新手总会在各种环节出现各种开发者意想不到的困难。即便软件官网上介绍了软件的各个工具的使用方法，但如果没有一个生动的例子来把这些应用工具的使用方法串起来的话，还是非常难的。

我自己作为一个PyNE系统的使用者，现在也参与了其中的部分开发工作，对此深有体会。对于新手的教学，告诉他们这个软件的各种功能是远远不够的，必须辅以实例，手把手的教会他们怎么用这个功能。用例子来说明这些问题是比较容易让新手接受的。我作为一个PyNE的使用者，也是一个新手。

这篇文章介绍了我如何从Trelis开始，一步一步实现计算一个简单模型的停机计量的过程的。

例题描述

现有如下形状的几何体：



图中每个小快代表着$ 10 \times 10 \times 10 \, \mathrm{cm}^3$的立方体。简单起见，我们把栅元1和3的材料设为密度为$1.0 \, \mathrm{g/cm^3}$的水(Water)，栅元2的材料设为密度为$7.8 \, \mathrm{g/cm^3}$的钢(Steel)。除了这几个栅元外其它区域都是真空。栅元1的正中心有一个点中子源，各项同性的发射能量为14 MeV的中子。

让点中子源以$10^{10}$个中子每秒的中子发射率辐照这几个栅元1年，然后停止辐照，考察停止辐照后1小时的栅元1, 2, 3的停机剂量率。

输入文件准备

要进行蒙特卡罗输运计算，我们需要提供以下信息：

  1. 几何信息，包括栅元和曲面。
  2. 材料信息。每个区域（栅元）是什么材料（材料的组成份额）
  3. 中子源信息。
  4. Tally信息。想统计的东西
  5. 其他信息。如截断卡，权重及权窗，计算终止条件等。

与普通MCNP计算不同，使用DAGMC计算时这些信息并不是全部都由一个输入文件提供，而是分别通过不同的文件提供。一个加载了几何和材料信息的h5m文件提供几何和材料信息。而其他信息则通过另一个input文本文件提供。

Trelis建模

几何体的建立

使用Trelis建立3个立方体，并在立方体外围设置一个graveyard区域，所有到达graveyard区域的中子将被杀死。graveyard的设置可以参考DAGMC关于graveyard的官方说明。



上图所示的三个小立方体块而就是我们建立的几个栅元，外面一个稍大的立方体薄壳就是我们想要的graveyard。 几何体分别标记为： 1. 蓝色小立方体，cell 1， trelis vol id: 1, 材料为密度为1的水，材料号设定为1 2. 黄色小立方体，cell 2, trelis vol id: 2, 材料为密度为7.8的铁，材料号设定为2 3. 红色小立方体块，cell 3, trelis vol id: 3， 材料为密度为1的水，材料号设定为1 4. 黄色大立方体薄壳，graveyard, trelis vol id: 8

分组并赋材料

几何体设置好了之后，对这几个栅元进行分组并赋材料。
在trelis界面输入下面命令：

1
2
trelis> group "mat:Water/rho:1.0" add vol 1 3
trelis> group "mat:Steel/rho:7.8" add vol 2

备注1：’trelis>’代表该命令是在trelis窗口中执行。

设置问题边界

边界包括终止边界和反射边界等。终止边界指的是我们不关心的区域，中子进入到终止边界后就会被杀死，这和MCNP中imp=0的区域相似。反射边界就是反射面，中子到达反射面后会被反射回去，在对称几何中很常用。

终止边界：graveyard

终止边界的设置实际上是通过对这个区域设置特殊材料实现的。比如此例题中需要将8号volume (可以简写为vol)设置为graveyard,可以使用下面命令进行设置：

1
trelis> group "mat:Graveyard" add volume 8
反射边界

反射边界是通过指定反射面的方式来设置的。与volume类似，我们可以把一些面加入到一个group，然后将这个group设置为反射面即可。反射面属性的设置也是通过把group名字设置为特殊字符串的方式实现。
比如将10, 11号曲面设置为反射面，可以通过下面的命令：

1
trelis> group "boundary:Reflecting" add surf 10 11

本例题中不需要将任何曲面设置为反射面。

imprint & merge

Imprint和merge是trelis中处理曲面的两个操作，他们的作用包括合并重合曲面和将无穷大的曲面根据几何体分割为总多的小的facet。功能说明可以参考trelis imprint and merge说明。

我们需要做的是，在trelis界面执行下面的操作：

1
2
trelis> imprint body all
trelis> merge all

导出几何文件

安装了DAGMC插件的trelis可以使用内置命令到处几何文件，使用方法可以通过下面命令查看帮助信息：

1
trelis> help dagmc

也可以参考官方文档。

对于本例题，我们可以在trelis中执行下面命令将几何信息存放到”geometry1.h5m”中：

1
trelis> export dagmc geometry1.h5m

至此，我们完成了几何文件的制作。

注意：在没有设置trelis默认导出文件路径的情况下，此时生成的”geometry1.h5m”文件并不一定在当前文件夹，可能存放于$HOME/Documents中，也可能在其他常用默认文件夹中

制作材料

我们需要首先制作本例题需要使用的材料库，包括Water和Steel两种材料。可以参考PyNE材料库制作。

下面是制作本例子中的”Water”和”Steel”用到的python脚本文件：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#!/usr/bin/python
#
from pyne.material import Material,MaterialLibrary
print "Welcome!"
mat_lib=MaterialLibrary()
#
mat2 = Material({'Fe':0.655,'Cr':0.170,'Ni':0.120,'Mo': 0.025,'Mn': 0.02, 'Si':.01},density=7.92)
mat2=mat2.expand_elements()
#
# define a simple water since O-18 not in mcnp xs libs
watervec={10010:2,80160:1} # simple water
water = Material()
water.density = 1.0
water.from_atom_frac(watervec)
#
mat_lib["Steel"]=mat2
mat_lib["Water"]=water
#
mat_lib.write_hdf5("test_material_lib.h5")
#
print "All done!"

备注：制作材料的时候会调用PyNE，需要保证此时系统LD_LIBRARY_PATH中的MOAB是4.9.1版本的。

整合几何与材料

按照UWUW流程制作材料，使用UWUW_preprocessor：

1
uwuw_preproc <dagmc h5m filename> -v -l <path to nuclear data library>

假设我们制作了材料库并存放在”material_lib1.h5”中，我们接下来就可以执行：

1
uwuw_preproc geometry1.h5m -v -l material_lib1.h5

DAGMC input文件准备

DAGMC的input文件中需要给定普通MCNP中除了栅元卡和曲面卡以及材料卡的其他部分，包括但不限于：

  1. 中子源定义，必要项
  2. 计算终止条件，必要项
  3. Tally卡，可选项，但一般都会有
  4. 其他功能性辅助卡，可选项，根据实际情况添加

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

1
2
3
4
5
6
7
8
9
10
11
C test1 input
C neutron source: point source, 14 MeV
sdef pos=0 0 0  erg=14  
C 
nps 1e5
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文件。

活化计算

中子输运完成后，我们开始正式调用PyNE的R2S计算模块的功能进行停机剂量率计算了。

进行停机剂量率计算，我们需要的主要信息包括：

  1. 每个网格单元的材料成分。必须项。由同时包含了几何和材料的geom.h5m文件提供。（此文件就是进行DAGMC计算的那个）
  2. 每个网格单元的中子通量能谱。必须项。由meshtal文件提供。
  3. 活化计算能力。必须项。由Analytic and Laplacian Adaptive Radioactivity Analysis软件:ALARA提供。
  4. 辐照及冷却方案。必须项。需要在R2S过程中指定。

在我们的工作文件夹下新建一个子文件夹用来存放R2S计算相关的文件。

1
2
mkdir r2s_run
cd r2s_run

初始化配置:setup

初始化配置r2s。执行：

1
r2s.py setup

执行之后我们应该能够看到下面的信息：

1
2
3
4
5
# 一些warning
File "config.ini" has been written
File "alara_params.txt" has been written
Fill out the fields in these filse then run ">> r2s.py step1"

修改配置文件

setup这一步会生成两个文件：config.ini和alara_params.txt，我们需要修改这两个文件中的一些设置，包括辐照方案的设置和一些计算细节的设置。

修改config.ini设置

一般情况下，我们需要注意的几个设置是：

  1. structured: 是否是结构化网格。结构化网格就是指MCNP中用的xyz或cyl网格。我们应该选True
  2. sub_voxel: 是否采用sub_voxel模式，sub_voxel模式是我正在添加的代码，还没有完成，暂时不能选用，因此是False
  3. reverse: 反转中子通量顺序。DAGMC输出的中子通量是按照从低能到高能排序，而用ALARA进行活化计算时需要按照从高能到低能读取，因此需要把中子通量逆序排列。这里填True。
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
[general]
# Specify whether this problem uses structured or unstructured mesh
structured: True
# Specify whether this problem uses sub-voxel r2s
sub_voxel: False

[step1]
# Path to MCNP MESHTAL file containing neutron fluxes or a DAG-MCNP5
# unstructured mesh tally .h5m file.
meshtal: meshtal
# Tally number within the meshtal file containing the fluxes for activation.
tally_num: 4
# The name of the tag used to store flux data on the mesh. For unstructured
# mesh this tag must already exist within the file specified in <meshtal>.
flux_tag: n_flux
# Path to the DAGMC material-laden geometry file (.h5m).
geom: geom.h5m
# If True the fluxes in the fluxin file will be printed in the reverse
# order of how they appear within the flux vector tag. Since MCNP and
# the Meshtal class order fluxes from low energy to high energy, this
# option should be true if the transmutation data being used is
# ordered from high-energy to low-energy.
reverse: True
# Number of rays to fire down each mesh row in each direction to calculate
# cell volume fractions.
num_rays: 10
# If true, rays will be fired down mesh rows in evenly spaced intervals.
# In this case <num_rays> must be a perfect square. If false, rays are fired
# down mesh rows in random intervals.
grid: False

[step2]
# List of decays times, seperated by commas. These strings much match exactly
# with their counterparts in the phtn_src file produced in step1. No spaces
# should appear in this line except the space between the time and the time unit
# for each entry.
decay_times:1 h
# The prefix of the .h5m files containing the source density distributations for
# each decay time.
output: source
# The name of the output files containing the total photon source intensities for
# each decay time
tot_phtn_src_intensities : total_photon_source_intensites.txt

修改alara_params.txt

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
material_lib alara_matlib
element_lib data/nuclib
data_library alaralib data/truncated_fendl2bin

cooling
    1 h
end

output zone
       integrate_energy
       # Energy group upper bounds. The lower bound is always zero.
       photon_source  data/truncated_fendl2bin  phtn_src 24 1.00E4 2.00E4 5.00E4 1.00E5
       2.00E5 3.00E5 4.00E5 6.00E5 8.00E5 1.00E6 1.22E6 1.44E6 1.66E6
       2.00E6 2.50E6 3.00E6 4.00E6 5.00E6 6.50E6 8.00E6 1.00E7 1.20E7
       1.40E7 2.00E7
end

#     flux name    fluxin file   norm   shift   unused
flux  my_flux     alara_fluxin  1e10     0      default

# Specify the irradiation schedule below.
# Syntax is found in the ALARA user manual
# This example is for a single 1 y pulse
schedule    my_schedule
    1 y my_flux my_pulse_history 0  s
end
pulsehistory  my_pulse_history
    1    0.0    s
end

#other parameters
truncation 1e-12
impurity 5e-6 1e-3
dump_file dump.file

r2s.py step1

r2s.py step1是用来根据r2s.py setup生成的配置文件生成ALARA输入文件，中子通量文件和ALARA材料库文件的。

当我们修改好上面两个文件后，就可以使用r2s.py生成ALARA计算需要的文件了：

1
r2s.py step1

材料数据库准备

ALARA输入文件中有三行是材料数据库相关的信息。

  1. material_lib alara_matlib 这个alara_matlib是alara需要读取的材料成功定义的文件名。这个文件是会由程序从test1_geometry.h5m (geom.h5m)中读取的，不需要做修改。
  2. element_lib data/nuclib 这个文件是定义元素的丰度用的。默认是自然丰度，也可以自己定义修改丰度。自己定义一个元素丰度数据表显然是非常麻烦的，在ALARA的安装文件中$ALARA/data下有一个elelib.std，可以把这个文件拷贝到工作文件夹下并重命名为nucib即可（也可以设置element_lib elelib.std)。但是这个elelib.std只能用于材料定义是按照元素定义的，如果有针对核素的定义，则这个就不能用了。需要将元素扩展成核素，需要后续介的方法单独处理。
  3. data_library alaralib truncated_fendl2bin （ALARA应该配备fendl3bin)但是并没有，ALARA安装文件中的sample/datam文件夹下只有truncated_fendl2bin数据库。我们可以把这个数据库相关的文件拷贝到工作文件夹。

当材料定义是针对核素定义了的化，普通的elelib.std就没办法使用了，会出现错误：

1
2
3
4
5
6
7
8
310: Could not find element <string> in element library.

The element string was not found in the element library. This could be due to an error in the material library, incorrect user input, or an omission in the element library.
``` 
这时我们需要使用工具**elelib_to_nuclib.py**来将元素扩展成核素数据。这个工具是一个python脚本，调用pyne数据模块扩展元素数据。工具的位置在ALARA源文件的tools文件夹下，使用方法：
```bash
cd r2s_run/data
python <path_to_elelib_to_nuclib.py> elelib_std -o nuclib

alara计算

当step1完成后，程序会给出运行alara的提示，按照提示运行：

1
alara alara_inp > output.txt

等待alara运行完成，会生成几个文件：

  1. output.txt: ALARA输出文件，里面有活化后材料的核素组成以及光子发射率数据
  2. phtn_src: ALARA输出文件，里面有详细的各个各个材料各个核素缠身的各个能群的光子数据

r2s step2

ALARA计算完成后，运行：

1
r2s.py step2

执行完成后，我们可以看到工作文件夹中多了几个文件：

  1. phtn_src.h5: HDF5格式的光子源数据，无法直接打开读取，只能程序读。
  2. source_i.h5m: 第i个冷却时间点的网格光子源数据。
  3. total_photon_source_intensites.txt: 里面是各个冷却时间节点的总光子发射率。
    你

光子输运

获得光子源后，我们就可以进行光子输运计算了。计算前必须确认DAGMC是已经包含了源子程序的版本，如果没有，需要先参考编译带用于R2S的源子程序的DAGMC安装DAGMC。

文件准备

进行光子输运计算，我们需要提供几何文件，材料，光子源，以及需要统计的信息。因此需要准备的文件有：

  1. geom.h5m: 就是中子输运时用的那个文件，即本例子中的test1_geometry.h5m，它包含了几何和材料信息
  2. source.h5m:光子源信息文件，源子程序需要从这里读取信息。由于MCNP无法从输入文件中传递字符串，所以这个文件名称不能改变，必须是source.h5m。
  3. input: 提供tally信息，另外需要在这个文件中定义idum，用于控制源子程序过程。

准备好这些文件后，我们还需要对这些文件进行一些小修改。

  1. input文件修改：首先，需要加上mode p， 然后修改tally部分。
  2. input文件修改：注释掉原本的SDEF卡

我们还需要修改栅元的重要性，因为中子输运时的文件里面的重要性是默认为中子的重要性，没有光子的重要性。但是这个重要性信息是写在h5m文件里面的，没法直接修改。

我们需要先运行一次这个文件：

1
mcnp5.mpi i=input g=geom.h5m

现在这个运行肯定会出现错误，说粒子产生子重要性为0的区域。但是同时这个过程会生成一个可编辑的lcad文件，我们可以修改这个lcad文件，在重要性部分添加上imp:p=1/imp:p=0，将这个文件重命名为lcad_m，删掉过程文件fcad, outp, runtpe后再次运行

1
mcnp5.mpi i=input g=geom.h5m l=lcad_m

之后就得到结果，至此停机剂量的计算完成。

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