PyNE unstructured R2S

这篇博客描述了使用非结构网格的PyNE R2S流程。首先确保Trelis，并安装DAGMC，安装PyNE, 安装ALARA。

基本流程

- 处理CAD文件，生成DAGMC使用的集成材料的几何文件
- 生成并导出tally用的非结构网格
- 中子输运计算
- PyNE R2S计算，中间使用ALARA进行活化计算
- 光子输运计算

使用Trelis处理DAGMC用的几何

Build up geometry, make several groups graveyard) (including and assign material for each group.

Imprint and merge

Use the following command:

1 2	trelis> imprint body all trelis> merge all

Export the geometry

The following command will export the geometry file into geom.h5m.

1	trelis> export dagmc geom.h5m

Build materials

Write a python script to build materials using PyNE, reference: define PyNE materials。

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

#!/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!"

整合几何与材料

按照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

生成并导出Tally用非结构网格

在”Command Panel”选择”Mode-Mesh” -> “Entity-Volumes” -> “Action-Mesh” -> “Tetmesh”，
然后选择需要划分网格的体积，可以使用这些方式进行选择：
- 输入体积的id号
- 将需要画网格的体积选中并加入picked, 然后在”Select Volumes”中输入”in picked”
- 将需要画网格的体积选中并设为一个block,然后在”Select Volumes”中输入”in block [id]”
(optional, suggested) 勾选”Number of Tets in Proximity”，根据需要的网格精细程度输入数值，如500
(optional, suggested) 取消勾选”Use Geomtric Sizing”。Geometric Sizing将会自动根据网格几何在边界区域划分更加精细的网格。对蒙卡计算来说并不需要，建议取消勾选。
(optional, suggested) 勾选”Advanced”, 并根据需要调整”Tet Optimization Level”,如选择”Standard”或”Medium”。蒙卡网格不需要太高的Optimization Level。
Apply Scheme
(optional, suggested) 如果检查出有重叠几何，建议对集合进行修正后再重新进行网格划分
(optional, suggested) 设置最大网格大小，如”vol xxx size 10”
Mesh
导出文件为”exodus”格式，如”tetmesh.exo”
使用mbconvert将网格文件转为h5m格式：mbconvert tetmesh.exo tetmesh.h5m

DAGMC input文件准备

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

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

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

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

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

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

C test1 input
C neutron source: point source, 14 MeV
sdef pos=0 0 0  erg=14  
C 
nps 1e5
C normal fmesh 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
C dagmc tetmesh tally card
fmesh14:n geom=dag emesh=0.1 20
fc14 dagmc type=unstr_track inp=tetmesh.h5m out=tetmesh_out.h5m

中子输运计算

使用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文件，和一个tetmesh_out.h5m文件。

PyNE R2S

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

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

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

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

1 2	mkdir r2s_run cd r2s_run

初始化配置:setup

初始化配置r2s。执行：

1	r2s.py setup

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

# 一些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设置

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

structured: 是否是结构化网格。结构化网格就是指MCNP中用的xyz或cyl网格。我们应该选False
reverse: 反转中子通量顺序。DAGMC输出的中子通量是按照从低能到高能排序，而用ALARA进行活化计算时需要按照从高能到低能读取，因此需要把中子通量逆序排列。这里填True。

[general]
# Specify whether this problem uses structured or unstructured mesh
structured: False
# 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: tetmesh.h5m
# 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

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计算

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

1	alara alara_inp > output.txt

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

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

r2s step2

ALARA计算完成后，运行：

1	r2s.py step2

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

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

光子输运

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

文件准备

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

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

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

input文件修改：首先，需要加上mode p，然后修改tally部分。

我们还需要修改栅元的重要性，因为中子输运时的文件里面的重要性是默认为中子的重要性，没有光子的重要性。但是这个重要性信息是写在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

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