OpenMC/docs/source/pythonapi/examples/tally-arithmetic.ipynb

1025 lines
35 KiB
Text

{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook shows the how tallies can be combined (added, subtracted, multiplied, etc.) using the Python API in order to create derived tallies. Since no covariance information is obtained, it is assumed that tallies are completely independent of one another when propagating uncertainties. The target problem is a simple pin cell.\n",
"\n",
"**Note:** that this Notebook was created using the latest Pandas v0.16.1. Everything in the Notebook will wun with older versions of Pandas, but the multi-indexing option in >v0.15.0 makes the tables look prettier."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import glob\n",
"from IPython.display import Image\n",
"import numpy as np\n",
"\n",
"import openmc\n",
"from openmc.statepoint import StatePoint\n",
"from openmc.summary import Summary\n",
"\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generate Input Files"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First we need to define materials that will be used in the problem. Before defining a material, we must create nuclides that are used in the material."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Instantiate some Nuclides\n",
"h1 = openmc.Nuclide('H-1')\n",
"b10 = openmc.Nuclide('B-10')\n",
"o16 = openmc.Nuclide('O-16')\n",
"u235 = openmc.Nuclide('U-235')\n",
"u238 = openmc.Nuclide('U-238')\n",
"zr90 = openmc.Nuclide('Zr-90')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the nuclides we defined, we will now create three materials for the fuel, water, and cladding of the fuel pin."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# 1.6 enriched fuel\n",
"fuel = openmc.Material(name='1.6% Fuel')\n",
"fuel.set_density('g/cm3', 10.31341)\n",
"fuel.add_nuclide(u235, 3.7503e-4)\n",
"fuel.add_nuclide(u238, 2.2625e-2)\n",
"fuel.add_nuclide(o16, 4.6007e-2)\n",
"\n",
"# borated water\n",
"water = openmc.Material(name='Borated Water')\n",
"water.set_density('g/cm3', 0.740582)\n",
"water.add_nuclide(h1, 4.9457e-2)\n",
"water.add_nuclide(o16, 2.4732e-2)\n",
"water.add_nuclide(b10, 8.0042e-6)\n",
"\n",
"# zircaloy\n",
"zircaloy = openmc.Material(name='Zircaloy')\n",
"zircaloy.set_density('g/cm3', 6.55)\n",
"zircaloy.add_nuclide(zr90, 7.2758e-3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With our three materials, we can now create a materials file object that can be exported to an actual XML file."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Instantiate a MaterialsFile, add Materials\n",
"materials_file = openmc.MaterialsFile()\n",
"materials_file.add_material(fuel)\n",
"materials_file.add_material(water)\n",
"materials_file.add_material(zircaloy)\n",
"materials_file.default_xs = '71c'\n",
"\n",
"# Export to \"materials.xml\"\n",
"materials_file.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's move on to the geometry. Our problem will have three regions for the fuel, the clad, and the surrounding coolant. The first step is to create the bounding surfaces -- in this case two cylinders and six reflective planes."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create cylinders for the fuel and clad\n",
"fuel_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.39218)\n",
"clad_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, R=0.45720)\n",
"\n",
"# Create boundary planes to surround the geometry\n",
"# Use both reflective and vacuum boundaries to make life interesting\n",
"min_x = openmc.XPlane(x0=-0.63, boundary_type='reflective')\n",
"max_x = openmc.XPlane(x0=+0.63, boundary_type='reflective')\n",
"min_y = openmc.YPlane(y0=-0.63, boundary_type='reflective')\n",
"max_y = openmc.YPlane(y0=+0.63, boundary_type='reflective')\n",
"min_z = openmc.ZPlane(z0=-0.63, boundary_type='reflective')\n",
"max_z = openmc.ZPlane(z0=+0.63, boundary_type='reflective')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the surfaces defined, we can now create cells that are defined by intersections of half-spaces created by the surfaces."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create a Universe to encapsulate a fuel pin\n",
"pin_cell_universe = openmc.Universe(name='1.6% Fuel Pin')\n",
"\n",
"# Create fuel Cell\n",
"fuel_cell = openmc.Cell(name='1.6% Fuel')\n",
"fuel_cell.fill = fuel\n",
"fuel_cell.add_surface(fuel_outer_radius, halfspace=-1)\n",
"pin_cell_universe.add_cell(fuel_cell)\n",
"\n",
"# Create a clad Cell\n",
"clad_cell = openmc.Cell(name='1.6% Clad')\n",
"clad_cell.fill = zircaloy\n",
"clad_cell.add_surface(fuel_outer_radius, halfspace=+1)\n",
"clad_cell.add_surface(clad_outer_radius, halfspace=-1)\n",
"pin_cell_universe.add_cell(clad_cell)\n",
"\n",
"# Create a moderator Cell\n",
"moderator_cell = openmc.Cell(name='1.6% Moderator')\n",
"moderator_cell.fill = water\n",
"moderator_cell.add_surface(clad_outer_radius, halfspace=+1)\n",
"pin_cell_universe.add_cell(moderator_cell)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"OpenMC requires that there is a \"root\" universe. Let us create a root cell that is filled by the pin cell universe and then assign it to the root universe."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create root Cell\n",
"root_cell = openmc.Cell(name='root cell')\n",
"root_cell.fill = pin_cell_universe\n",
"\n",
"# Add boundary planes\n",
"root_cell.add_surface(min_x, halfspace=+1)\n",
"root_cell.add_surface(max_x, halfspace=-1)\n",
"root_cell.add_surface(min_y, halfspace=+1)\n",
"root_cell.add_surface(max_y, halfspace=-1)\n",
"root_cell.add_surface(min_z, halfspace=+1)\n",
"root_cell.add_surface(max_z, halfspace=-1)\n",
"\n",
"# Create root Universe\n",
"root_universe = openmc.Universe(universe_id=0, name='root universe')\n",
"root_universe.add_cell(root_cell)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now must create a geometry that is assigned a root universe, put the geometry into a geometry file, and export it to XML."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create Geometry and set root Universe\n",
"geometry = openmc.Geometry()\n",
"geometry.root_universe = root_universe"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Instantiate a GeometryFile\n",
"geometry_file = openmc.GeometryFile()\n",
"geometry_file.geometry = geometry\n",
"\n",
"# Export to \"geometry.xml\"\n",
"geometry_file.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the geometry and materials finished, we now just need to define simulation parameters. In this case, we will use 5 inactive batches and 15 active batches each with 2500 particles."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# OpenMC simulation parameters\n",
"batches = 20\n",
"inactive = 5\n",
"particles = 2500\n",
"\n",
"# Instantiate a SettingsFile\n",
"settings_file = openmc.SettingsFile()\n",
"settings_file.batches = batches\n",
"settings_file.inactive = inactive\n",
"settings_file.particles = particles\n",
"settings_file.output = {'tallies': True, 'summary': True}\n",
"source_bounds = [-0.63, -0.63, -0.63, 0.63, 0.63, 0.63]\n",
"settings_file.set_source_space('box', source_bounds)\n",
"\n",
"# Export to \"settings.xml\"\n",
"settings_file.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us also create a plot file that we can use to verify that our pin cell geometry was created successfully."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Instantiate a Plot\n",
"plot = openmc.Plot(plot_id=1)\n",
"plot.filename = 'materials-xy'\n",
"plot.origin = [0, 0, 0]\n",
"plot.width = [1.26, 1.26]\n",
"plot.pixels = [250, 250]\n",
"plot.color = 'mat'\n",
"\n",
"# Instantiate a PlotsFile, add Plot, and export to \"plots.xml\"\n",
"plot_file = openmc.PlotsFile()\n",
"plot_file.add_plot(plot)\n",
"plot_file.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the plots.xml file, we can now generate and view the plot. OpenMC outputs plots in .ppm format, which can be converted into a compressed format like .png with the convert utility."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Run openmc in plotting mode\n",
"executor = openmc.Executor()\n",
"executor.plot_geometry(output=False)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6AgMAAAD1grKuAAAABGdBTUEAALGPC/xhBQAAAAFzUkdC\nAK7OHOkAAAAgY0hSTQAAeiYAAICEAAD6AAAAgOgAAHUwAADqYAAAOpgAABdwnLpRPAAAAAxQTFRF\n////chIS6YCRTb/E6kGE+wAAAAFiS0dEAIgFHUgAAAAJcEhZcwAAAEgAAABIAEbJaz4AAALKSURB\nVGje7dpLcqQwDAbgHHE2YeEj+D4cwQucBUfo+3CEXoSp8OhuhF70T4qpKXmdr21LogK2Pj7A8QmN\nP+HDhw8fPnz48Kf6VH9G+66vy+je8k19jnf8C5dXIPv86ms56lPdjvaYbyodx3ze+XLE76cXFiD4\nzPji99z0/AJ4n1lfvJ6fnl0A6x+578efMSg1wPr172/jPO5yFXM+Ef78gdblM+WPHyguP//t1/g6\npA0wfln+ho/fwgYYn19C/xwDvwHGc9OvC+hs37DTrwuwfWanXxdQTC9Mvyygs3wjTL8uwPJpn/tN\nDbSGz7T0SBEWw4vLXzbQ6b6RoveIoO6TvPxlA63qs7z8ZQPF9F+SH22vbX8OQKf5Rtv+EgDNJ3X5\n8wZaxWd1+fMGiuFvir8bvjp8J/tGy/6jAmRvhW8fwL3vVT+o3grfPoB7r/IpALI3tz8FoJN84/NV\n873hB8UnM3xzANtf8nb4dwmg3grfFEDJO8JPE0i9Ff4pAYL3pI8mkHor/HMCeO9JH00g9SafEsh7\nT/ppARBvp48UwJnelT5SACd7O31TAlnvKx9SQCd7B58KgPO+8iMFuPWe9E8F8BveWX7bAjzX9y4/\n/Jve+fhsH6Ctv7n8PTzjvY/v9gEOHz58+PBX+6v/f/wPvnd54f3j6venE/yl769Xv7+j3x/o98/V\n32/o9+fl389Xnx+g5x/o+Qt6/oOeP6HnX+j5G3z+h54/ouefV5/foufP6Pk3ev4On/+j9w/o/Qd6\n/4Le/6D3T/D9V67Y/ZsVQBq+s+8f0ftP+P41axXguP9NWgDuu/Cdfv+N3r/D9/9TAID+A7T/Ae2/\ngPs/0P4TtP8F7r9J3AIO9P+g/Udw/9Oygbf7r9D+L7j/DO1/Q/vv4P4/tP8Q7n9E+y/h/k+0/xTu\nf4X7b+H+X7T/+BPuf3aM8OHDhw8fPnz4w/4vzcvgeY10sY0AAAAldEVYdGRhdGU6Y3JlYXRlADIw\nMTUtMTAtMDNUMDE6MDM6MjktMDQ6MDDFeHPZAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDE1LTEwLTAz\nVDAxOjAzOjI5LTA0OjAwtCXLZQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Convert OpenMC's funky ppm to png\n",
"!convert materials-xy.ppm materials-xy.png\n",
"\n",
"# Display the materials plot inline\n",
"Image(filename='materials-xy.png')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we can see from the plot, we have a nice pin cell with fuel, cladding, and water! Before we run our simulation, we need to tell the code what we want to tally. The following code shows how to create a variety of tallies."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Instantiate an empty TalliesFile\n",
"tallies_file = openmc.TalliesFile()"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create Tallies to compute microscopic multi-group cross-sections\n",
"\n",
"# Instantiate energy filter for multi-group cross-section Tallies\n",
"energy_filter = openmc.Filter(type='energy', bins=[0., 0.625e-6, 20.])\n",
"\n",
"# Instantiate flux Tally in moderator and fuel\n",
"tally = openmc.Tally(name='flux')\n",
"tally.add_filter(openmc.Filter(type='cell', bins=[fuel_cell.id, moderator_cell.id]))\n",
"tally.add_filter(energy_filter)\n",
"tally.add_score('flux')\n",
"tallies_file.add_tally(tally)\n",
"\n",
"# Instantiate reaction rate Tally in fuel\n",
"tally = openmc.Tally(name='fuel rxn rates')\n",
"tally.add_filter(openmc.Filter(type='cell', bins=[fuel_cell.id]))\n",
"tally.add_filter(energy_filter)\n",
"tally.add_score('nu-fission')\n",
"tally.add_score('scatter')\n",
"tally.add_nuclide(u238)\n",
"tally.add_nuclide(u235)\n",
"tallies_file.add_tally(tally)\n",
"\n",
"# Instantiate reaction rate Tally in moderator\n",
"tally = openmc.Tally(name='moderator rxn rates')\n",
"tally.add_filter(openmc.Filter(type='cell', bins=[moderator_cell.id]))\n",
"tally.add_filter(energy_filter)\n",
"tally.add_score('absorption')\n",
"tally.add_score('total')\n",
"tally.add_nuclide(o16)\n",
"tally.add_nuclide(h1)\n",
"tallies_file.add_tally(tally)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# K-Eigenvalue (infinity) tallies\n",
"fiss_rate = openmc.Tally(name='fiss. rate')\n",
"abs_rate = openmc.Tally(name='abs. rate')\n",
"fiss_rate.add_score('nu-fission')\n",
"abs_rate.add_score('absorption')\n",
"tallies_file.add_tally(fiss_rate)\n",
"tallies_file.add_tally(abs_rate)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Resonance Escape Probability tallies\n",
"therm_abs_rate = openmc.Tally(name='therm. abs. rate')\n",
"therm_abs_rate.add_score('absorption')\n",
"therm_abs_rate.add_filter(openmc.Filter(type='energy', bins=[0., 0.625]))\n",
"tallies_file.add_tally(therm_abs_rate)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Thermal Flux Utilization tallies\n",
"fuel_therm_abs_rate = openmc.Tally(name='fuel therm. abs. rate')\n",
"fuel_therm_abs_rate.add_score('absorption')\n",
"fuel_therm_abs_rate.add_filter(openmc.Filter(type='energy', bins=[0., 0.625]))\n",
"fuel_therm_abs_rate.add_filter(openmc.Filter(type='cell', bins=[fuel_cell.id]))\n",
"tallies_file.add_tally(fuel_therm_abs_rate)"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Fast Fission Factor tallies\n",
"therm_fiss_rate = openmc.Tally(name='therm. fiss. rate')\n",
"therm_fiss_rate.add_score('nu-fission')\n",
"therm_fiss_rate.add_filter(openmc.Filter(type='energy', bins=[0., 0.625]))\n",
"tallies_file.add_tally(therm_fiss_rate)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Instantiate energy filter to illustrate Tally slicing\n",
"energy_filter = openmc.Filter(type='energy', bins=np.logspace(np.log10(1e-8), np.log10(20), 10))\n",
"\n",
"# Instantiate flux Tally in moderator and fuel\n",
"tally = openmc.Tally(name='need-to-slice')\n",
"tally.add_filter(openmc.Filter(type='cell', bins=[fuel_cell.id, moderator_cell.id]))\n",
"tally.add_filter(energy_filter)\n",
"tally.add_score('nu-fission')\n",
"tally.add_score('scatter')\n",
"tally.add_nuclide(h1)\n",
"tally.add_nuclide(u238)\n",
"tallies_file.add_tally(tally)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Export to \"tallies.xml\"\n",
"tallies_file.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we a have a complete set of inputs, so we can go ahead and run our simulation."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" .d88888b. 888b d888 .d8888b.\n",
" d88P\" \"Y88b 8888b d8888 d88P Y88b\n",
" 888 888 88888b.d88888 888 888\n",
" 888 888 88888b. .d88b. 88888b. 888Y88888P888 888 \n",
" 888 888 888 \"88b d8P Y8b 888 \"88b 888 Y888P 888 888 \n",
" 888 888 888 888 88888888 888 888 888 Y8P 888 888 888\n",
" Y88b. .d88P 888 d88P Y8b. 888 888 888 \" 888 Y88b d88P\n",
" \"Y88888P\" 88888P\" \"Y8888 888 888 888 888 \"Y8888P\"\n",
"__________________888______________________________________________________\n",
" 888\n",
" 888\n",
"\n",
" Copyright: 2011-2015 Massachusetts Institute of Technology\n",
" License: http://mit-crpg.github.io/openmc/license.html\n",
" Version: 0.7.0\n",
" Git SHA1: e0c2aace2e73367536fa03e153b67a2d038cd2b3\n",
" Date/Time: 2015-10-03 01:03:29\n",
" MPI Processes: 1\n",
"\n",
" ===========================================================================\n",
" ========================> INITIALIZATION <=========================\n",
" ===========================================================================\n",
"\n",
" Reading settings XML file...\n",
" Reading cross sections XML file...\n",
" Reading geometry XML file...\n",
" Reading materials XML file...\n",
" Reading tallies XML file...\n",
" Building neighboring cells lists for each surface...\n",
" Loading ACE cross section table: 92238.71c\n",
" Loading ACE cross section table: 8016.71c\n",
" Loading ACE cross section table: 92235.71c\n",
" Loading ACE cross section table: 5010.71c\n",
" Loading ACE cross section table: 1001.71c\n",
" Loading ACE cross section table: 40090.71c\n",
" Initializing source particles...\n",
"\n",
" ===========================================================================\n",
" ====================> K EIGENVALUE SIMULATION <====================\n",
" ===========================================================================\n",
"\n",
" Bat./Gen. k Average k \n",
" ========= ======== ==================== \n",
" 1/1 1.00279 \n",
" 2/1 1.03320 \n",
" 3/1 1.04467 \n",
" 4/1 1.09693 \n",
" 5/1 1.05008 \n",
" 6/1 1.08426 \n",
" 7/1 1.05363 1.06894 +/- 0.01531\n",
" 8/1 0.97961 1.03917 +/- 0.03106\n",
" 9/1 1.06444 1.04549 +/- 0.02285\n",
" 10/1 1.08345 1.05308 +/- 0.01926\n",
" 11/1 1.06871 1.05568 +/- 0.01594\n",
" 12/1 1.03183 1.05228 +/- 0.01390\n",
" 13/1 1.04486 1.05135 +/- 0.01207\n",
" 14/1 1.06468 1.05283 +/- 0.01075\n",
" 15/1 1.04185 1.05173 +/- 0.00968\n",
" 16/1 1.01268 1.04818 +/- 0.00944\n",
" 17/1 1.04129 1.04761 +/- 0.00864\n",
" 18/1 1.01127 1.04481 +/- 0.00843\n",
" 19/1 1.03738 1.04428 +/- 0.00782\n",
" 20/1 1.04410 1.04427 +/- 0.00728\n",
" Creating state point statepoint.20.h5...\n",
"\n",
" ===========================================================================\n",
" ======================> SIMULATION FINISHED <======================\n",
" ===========================================================================\n",
"\n",
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 4.1300E-01 seconds\n",
" Reading cross sections = 9.0000E-02 seconds\n",
" Total time in simulation = 2.1398E+01 seconds\n",
" Time in transport only = 2.1378E+01 seconds\n",
" Time in inactive batches = 2.0260E+00 seconds\n",
" Time in active batches = 1.9372E+01 seconds\n",
" Time synchronizing fission bank = 2.0000E-03 seconds\n",
" Sampling source sites = 0.0000E+00 seconds\n",
" SEND/RECV source sites = 1.0000E-03 seconds\n",
" Time accumulating tallies = 1.0000E-03 seconds\n",
" Total time for finalization = 3.0000E-03 seconds\n",
" Total time elapsed = 2.1823E+01 seconds\n",
" Calculation Rate (inactive) = 6169.79 neutrons/second\n",
" Calculation Rate (active) = 1935.78 neutrons/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
" k-effective (Collision) = 1.04044 +/- 0.00527\n",
" k-effective (Track-length) = 1.04427 +/- 0.00728\n",
" k-effective (Absorption) = 1.04794 +/- 0.00535\n",
" Combined k-effective = 1.04628 +/- 0.00467\n",
" Leakage Fraction = 0.00000 +/- 0.00000\n",
"\n"
]
},
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Remove old HDF5 (summary, statepoint) files\n",
"!rm statepoint.*\n",
"\n",
"# Run OpenMC with MPI!\n",
"executor.run_simulation()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Tally Data Processing"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our simulation ran successfully and created a statepoint file with all the tally data in it. We begin our analysis here loading the statepoint file and 'reading' the results. By default, the tally results are not read into memory because they might be large, even large enough to exceed the available memory on a computer."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [],
"source": [
"# Load the statepoint file\n",
"sp = StatePoint('statepoint.20.h5')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You may have also noticed we instructed OpenMC to create a summary file with lots of geometry information in it. This can help to produce more sensible output from the Python API, so we will use the summary file to link against."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [
{
"ename": "KeyError",
"evalue": "10003",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mKeyError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-25-3cd87b8121ef>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[1;31m# Load the summary file and link with statepoint\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 2\u001b[0m \u001b[0msu\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mSummary\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'summary.h5'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 3\u001b[1;33m \u001b[0msp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlink_with_summary\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msu\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[1;32m/usr/local/lib/python2.7/dist-packages/openmc-0.7.0-py2.7.egg/openmc/statepoint.pyc\u001b[0m in \u001b[0;36mlink_with_summary\u001b[1;34m(self, summary)\u001b[0m\n\u001b[0;32m 610\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mtally_id\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtally\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtallies\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 611\u001b[0m \u001b[1;31m# Get the Tally name from the summary file\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 612\u001b[1;33m \u001b[0mtally\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mname\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0msummary\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtallies\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mtally_id\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 613\u001b[0m \u001b[0mtally\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwith_summary\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mTrue\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 614\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mKeyError\u001b[0m: 10003"
]
}
],
"source": [
"# Load the summary file and link with statepoint\n",
"su = Summary('summary.h5')\n",
"sp.link_with_summary(su)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We have a tally of the total fission rate and the total absorption rate, so we can calculate k-infinity as:\n",
"$$k_\\infty = \\frac{\\langle \\nu \\Sigma_f \\phi \\rangle}{\\langle \\Sigma_a \\phi \\rangle}$$\n",
"In this notation, $\\langle \\cdot \\rangle^a_b$ represents an OpenMC that is integrated over region $a$ and energy range $b$. If $a$ or $b$ is not reported, it means the value represents an integral over all space or all energy, respectively."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Compute k-infinity using tally arithmetic\n",
"fiss_rate = sp.get_tally(name='fiss. rate')\n",
"abs_rate = sp.get_tally(name='abs. rate')\n",
"keff = fiss_rate / abs_rate\n",
"keff.get_pandas_dataframe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that even though the neutron production rate and absorption rate are separate tallies, we still get a first-order estimate of the uncertainty on the quotient of them automatically!\n",
"\n",
"Often in textbooks you'll see k-infinity represented using the four-factor formula $$k_\\infty = p \\epsilon f \\eta.$$ Let's analyze each of these factors, starting with the resonance escape probability which is defined as $$p=\\frac{\\langle\\Sigma_a\\phi\\rangle_T}{\\langle\\Sigma_a\\phi\\rangle}$$ where the subscript $T$ means thermal energies."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Compute resonance escape probability using tally arithmetic\n",
"therm_abs_rate = sp.get_tally(name='therm. abs. rate')\n",
"res_esc = therm_abs_rate / abs_rate\n",
"res_esc.get_pandas_dataframe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The fast fission factor can be calculated as\n",
"$$\\epsilon=\\frac{\\langle\\nu\\Sigma_f\\phi\\rangle}{\\langle\\nu\\Sigma_f\\phi\\rangle_T}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Compute fast fission factor factor using tally arithmetic\n",
"therm_fiss_rate = sp.get_tally(name='therm. fiss. rate')\n",
"fast_fiss = fiss_rate / therm_fiss_rate\n",
"fast_fiss.get_pandas_dataframe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The thermal flux utilization is calculated as\n",
"$$f=\\frac{\\langle\\Sigma_a\\phi\\rangle^F_T}{\\langle\\Sigma_a\\phi\\rangle_T}$$\n",
"where the superscript $F$ denotes fuel."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Compute thermal flux utilization factor using tally arithmetic\n",
"fuel_therm_abs_rate = sp.get_tally(name='fuel therm. abs. rate')\n",
"therm_util = fuel_therm_abs_rate / therm_abs_rate\n",
"therm_util.get_pandas_dataframe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The final factor is the number of fission neutrons produced per absorption in fuel, calculated as $$\\eta = \\frac{\\langle \\nu\\Sigma_f\\phi \\rangle_T}{\\langle \\Sigma_a \\phi \\rangle^F_T}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Compute neutrons produced per absorption (eta) using tally arithmetic\n",
"eta = therm_fiss_rate / fuel_therm_abs_rate\n",
"eta.get_pandas_dataframe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we can calculate $k_\\infty$ using the product of the factors form the four-factor formula."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"keff = res_esc * fast_fiss * therm_util * eta\n",
"keff.get_pandas_dataframe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that the value we've obtained here has exactly the same mean as before. However, because of the way it was calculated, the standard deviation appears to be larger.\n",
"\n",
"Let's move on to a more complicated example now. Before we set up tallies to get reaction rates in the fuel and moderator in two energy groups for two different nuclides. We can use tally arithmetic to divide each of these reaction rates by the flux to get microscopic multi-group cross sections."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [],
"source": [
"# Compute microscopic multi-group cross-sections\n",
"flux = sp.get_tally(name='flux')\n",
"flux = flux.get_slice(filters=['cell'], filter_bins=[(fuel_cell.id,)])\n",
"fuel_rxn_rates = sp.get_tally(name='fuel rxn rates')\n",
"mod_rxn_rates = sp.get_tally(name='moderator rxn rates')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fuel_xs = fuel_rxn_rates / flux\n",
"fuel_xs.get_pandas_dataframe()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that when the two tallies with multiple bins were divided, the derived tally contains the outer product of the combinations. If the filters/scores are the same, no outer product is needed. The `get_values(...)` method allows us to obtain a subset of tally scores. In the following example, we obtain just the neutron production microscopic cross sections."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Show how to use Tally.get_values(...) with a CrossScore\n",
"nu_fiss_xs = fuel_xs.get_values(scores=['(nu-fission / flux)'])\n",
"print(nu_fiss_xs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The same idea can be used not only for scores but also for filters and nuclides."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Show how to use Tally.get_values(...) with a CrossScore and CrossNuclide\n",
"u235_scatter_xs = fuel_xs.get_values(nuclides=['(U-235 / total)'], \n",
" scores=['(scatter / flux)'])\n",
"print(u235_scatter_xs)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Show how to use Tally.get_values(...) with a CrossFilter and CrossScore\n",
"fast_scatter_xs = fuel_xs.get_values(filters=['energy'], \n",
" filter_bins=[((0.625e-6, 20.),)], \n",
" scores=['(scatter / flux)'])\n",
"print(fast_scatter_xs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A more advanced method is to use `get_slice(...)` to create a new derived tally that is a subset of an existing tally. This has the benefit that we can use `get_pandas_dataframe()` to see the tallies in a more human-readable format."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# \"Slice\" the nu-fission data into a new derived Tally\n",
"nu_fission_rates = fuel_rxn_rates.get_slice(scores=['nu-fission'])\n",
"nu_fission_rates.get_pandas_dataframe()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# \"Slice\" the H-1 scatter data in the moderator Cell into a new derived Tally\n",
"need_to_slice = sp.get_tally(name='need-to-slice')\n",
"slice_test = need_to_slice.get_slice(scores=['scatter'], nuclides=['H-1'],\n",
" filters=['cell'], filter_bins=[(moderator_cell.id,)])\n",
"slice_test.get_pandas_dataframe()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 0
}