crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-22 06:55:35 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5
OpenMC/docs/source/pythonapi/examples/pandas-dataframes.ipynb

1119 lines
37 KiB
Text
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This notebook demonstrates how systematic analysis of tally scores is possible using Pandas dataframes. A dataframe can be automatically generated using the `Tally.get_pandas_dataframe(...)` method. Furthermore, by linking the tally data in a statepoint file with geometry and material information from a summary file, the dataframe can be shown with user-supplied labels.\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": [
"import glob\n",
"from IPython.display import Image\n",
"import matplotlib.pylab as pylab\n",
"import scipy.stats\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": 2,
"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": 3,
"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": 4,
"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. This problem will be a square array of fuel pins, which we can use OpenMC's lattice/universe feature for. The basic universe will have three regions for the fuel, the clad, and the surrounding coolant. The first step is to create the bounding surfaces for fuel and clad, as well as the outer bounding surfaces of the problem."
]
},
{
"cell_type": "code",
"execution_count": 5,
"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=-10.71, boundary_type='reflective')\n",
"max_x = openmc.XPlane(x0=+10.71, boundary_type='vacuum')\n",
"min_y = openmc.YPlane(y0=-10.71, boundary_type='vacuum')\n",
"max_y = openmc.YPlane(y0=+10.71, boundary_type='reflective')\n",
"min_z = openmc.ZPlane(z0=-10.71, boundary_type='reflective')\n",
"max_z = openmc.ZPlane(z0=+10.71, 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": 6,
"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": [
"Using the pin cell universe, we can construct a 17x17 rectangular lattice with a 1.26cm pitch."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create fuel assembly Lattice\n",
"assembly = openmc.RectLattice(name='1.6% Fuel - 0BA')\n",
"assembly.dimension = (17, 17)\n",
"assembly.pitch = (1.26, 1.26)\n",
"assembly.lower_left = [-1.26 * 17. / 2.0] * 2\n",
"assembly.universes = [[pin_cell_universe] * 17] * 17"
]
},
{
"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 = assembly\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": true
},
"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 minimum active batches each with 2500 particles. We also tell OpenMC to turn tally triggers on, which means it will keep running until some criterion on the uncertainty of tallies is reached."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# OpenMC simulation parameters\n",
"min_batches = 20\n",
"max_batches = 200\n",
"inactive = 5\n",
"particles = 2500\n",
"\n",
"# Instantiate a SettingsFile\n",
"settings_file = openmc.SettingsFile()\n",
"settings_file.batches = min_batches\n",
"settings_file.inactive = inactive\n",
"settings_file.particles = particles\n",
"settings_file.output = {'tallies': False, 'summary': True}\n",
"settings_file.trigger_active = True\n",
"settings_file.trigger_max_batches = max_batches\n",
"source_bounds = [-10.71, -10.71, -10, 10.71, 10.71, 10.]\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 = [21.5, 21.5]\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+wAAAAFiS0dEAIgFHUgAAAAJcEhZcwAAAEgAAABIAEbJaz4AAAPZSURB\nVGje7Zs7buMwEIZ9iey50gyNjQpXKTYudIScgkdQYTfut1idwkdQkQNsYQO2Qj0sPiVK+mlQDmwg\nwIcgg8Cc4fCTSK5W4OeFkM8rHv+2I/rgxPZEPZgR7XtQxKdXYuUXJSUnBQ/9WCgo4vOSJ+WFUvF7\nE08mlia+rn7VcKXP8sRszFX8b2MdX2y6v1Tw6MZUw4H4ojfIjD8mvn/qRL5p4+vvlMqvp2EhR8WB\nzfiz20hXORmP9fi/bM9EeUFvV5H/0yRkeSbiGRfFJErxD9ENdz7Mbhig/h89fvtFdMiI/ePUIXV4\nlXju8K3DKv9NThOZ3q2KmUy6grxFES8rjeyic+FFQav+ncg3fXjH+Ts+/iibztFqOiZuZP/Z3Oaf\nPX40NGgST2r+uvQkXXp6cKvmr+r0e1Eef5um3+JHP3IFF1D/seNZJgaDmvY0Gav1s+2f1fqpIcub\nlfKGt6apotG/NVx3SInWtLX+7Vg/Pv1YqOsnun6JSVdOXT/X7vk75f938QP+8OmSBs0fXtymMhJb\nf8qlPynYmpKCh7OB1fzNalOj1sl0ZAruHLiA+RM73pDe/VjMVP89+aTXwjyc/x5n+u991895/utr\nJTy8/06TXh0r/5JOa2JmYmqi4r/vUm/H4wLmT+z4anhr05X+q6KUXhtzr/9qSff5L5uMT//V/NdU\n4YuBTPa/8P67l/6r44ds+hYuoP5jx9ciy6XTWlibBrmx8V/TdMfjkP+6pOsu/lvM9N90sf7r+f6m\n/65n+S8p/itN15v0UkW3/+48+PRfJX6S9Joo4g+G/1qYG9KroqP/WypcuvyXPf13wH89/hHef7MB\n6R3Cqn55U4rv4kfH3zaSgQuYP7HjVf89tXrbO+hfLdr+Ozv/SP1dgtQ/Ov8C+i/3+q/Zf2D/HWi6\nbjT6rym9I/v/03/b+LHS4cTg/utTsV7/net/Afzz4f0XGX84/2j9xZ4/sePR/of2X7D/o+vPo/sv\n6h9B/Bfxr9j1Hz2eN/hO8/wfff4A848+f/1A/530/I0+/8PvH9D3H9HnT+R49P0b+v4PfP/4E/wX\nfP8Mvf9G37/D/ovuP8SeP7Hj0f0vdP8tqP9O339cyv7p3P1fdP8Z3v9G999j13/seMax8x/o+ZN7\n+O+E8zdP/8XOf8Hnz9Dzb7HnT+x49PxlCp7/BM+fOv13wvnXBfivt2lMvD8TyH/Hnb+Gz3+j589j\nz5/Y8ej9h4D+W7qQmf57efqv239n3T+C7z+h969i13/seMax+3/o/cMcu/8Y2H9n3p+J6r98pv8m\n4fwXuH+M3n+OO3++AX9clR+4PhbRAAAAJXRFWHRkYXRlOmNyZWF0ZQAyMDE1LTEwLTAzVDAxOjE0\nOjM0LTA0OjAwS1fWPgAAACV0RVh0ZGF0ZTptb2RpZnkAMjAxNS0xMC0wM1QwMToxNDozNC0wNDow\nMDoKboIAAAAASUVORK5CYII=\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 array of pin cells 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": true
},
"outputs": [],
"source": [
"# Instantiate an empty TalliesFile\n",
"tallies_file = openmc.TalliesFile()\n",
"tallies_file._tallies = []"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Instantiate a fission rate mesh Tally"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Instantiate a tally Mesh\n",
"mesh = openmc.Mesh(mesh_id=1)\n",
"mesh.type = 'regular'\n",
"mesh.dimension = [17, 17]\n",
"mesh.lower_left = [-10.71, -10.71]\n",
"mesh.width = [1.26, 1.26]\n",
"\n",
"# Instantiate tally Filter\n",
"mesh_filter = openmc.Filter()\n",
"mesh_filter.mesh = mesh\n",
"\n",
"# Instantiate energy Filter\n",
"energy_filter = openmc.Filter()\n",
"energy_filter.type = 'energy'\n",
"energy_filter.bins = np.array([0, 0.625e-6, 20.])\n",
"\n",
"# Instantiate the Tally\n",
"tally = openmc.Tally(name='mesh tally')\n",
"tally.add_filter(mesh_filter)\n",
"tally.add_filter(energy_filter)\n",
"tally.add_score('fission')\n",
"tally.add_score('nu-fission')\n",
"\n",
"# Add mesh and Tally to TalliesFile\n",
"tallies_file.add_mesh(mesh)\n",
"tallies_file.add_tally(tally)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Instantiate a cell Tally with nuclides"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Instantiate tally Filter\n",
"cell_filter = openmc.Filter(type='cell', bins=[fuel_cell.id])\n",
"\n",
"# Instantiate the tally\n",
"tally = openmc.Tally(name='cell tally')\n",
"tally.add_filter(cell_filter)\n",
"tally.add_score('scatter-y2')\n",
"tally.add_nuclide(u235)\n",
"tally.add_nuclide(u238)\n",
"\n",
"# Add mesh and tally to TalliesFile\n",
"tallies_file.add_tally(tally)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create a \"distribcell\" Tally. The distribcell filter allows us to tally multiple repeated instances of the same cell throughout the geometry."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Instantiate tally Filter\n",
"distribcell_filter = openmc.Filter(type='distribcell', bins=[moderator_cell.id])\n",
"\n",
"# Instantiate tally Trigger for kicks\n",
"trigger = openmc.Trigger(trigger_type='std_dev', threshold=5e-5)\n",
"trigger.add_score('absorption')\n",
"\n",
"# Instantiate the Tally\n",
"tally = openmc.Tally(name='distribcell tally')\n",
"tally.add_filter(distribcell_filter)\n",
"tally.add_score('absorption')\n",
"tally.add_score('scatter')\n",
"tally.add_trigger(trigger)\n",
"\n",
"# Add mesh and tally to TalliesFile\n",
"tallies_file.add_tally(tally)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"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": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"rm: cannot remove ‘statepoint.*’: No such file or directory\n",
"\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:14:34\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 0.59998 \n",
" 2/1 0.65473 \n",
" 3/1 0.67452 \n",
" 4/1 0.66458 \n",
" 5/1 0.70093 \n",
" 6/1 0.70726 \n",
" 7/1 0.65977 0.68351 +/- 0.02375\n",
" 8/1 0.68457 0.68387 +/- 0.01372\n",
" 9/1 0.70024 0.68796 +/- 0.01053\n",
" 10/1 0.64895 0.68016 +/- 0.01128\n",
" 11/1 0.68744 0.68137 +/- 0.00929\n",
" 12/1 0.68037 0.68123 +/- 0.00786\n",
" 13/1 0.64865 0.67715 +/- 0.00793\n",
" 14/1 0.71415 0.68127 +/- 0.00811\n",
" 15/1 0.65717 0.67886 +/- 0.00764\n",
" 16/1 0.71598 0.68223 +/- 0.00769\n",
" 17/1 0.67285 0.68145 +/- 0.00707\n",
" 18/1 0.69329 0.68236 +/- 0.00656\n",
" 19/1 0.65696 0.68055 +/- 0.00634\n",
" 20/1 0.65500 0.67884 +/- 0.00615\n",
" Triggers unsatisfied, max unc./thresh. is 1.21110 for absorption in tally 10002\n",
" The estimated number of batches is 28\n",
" Creating state point statepoint.020.h5...\n",
" 21/1 0.67090 0.67835 +/- 0.00577\n",
" 22/1 0.69025 0.67905 +/- 0.00546\n",
" 23/1 0.66113 0.67805 +/- 0.00525\n",
" 24/1 0.67934 0.67812 +/- 0.00496\n",
" 25/1 0.67203 0.67781 +/- 0.00472\n",
" 26/1 0.66928 0.67741 +/- 0.00451\n",
" 27/1 0.70271 0.67856 +/- 0.00445\n",
" 28/1 0.70233 0.67959 +/- 0.00437\n",
" Triggers satisfied for batch 28\n",
" Creating state point statepoint.028.h5...\n",
"\n",
" ===========================================================================\n",
" ======================> SIMULATION FINISHED <======================\n",
" ===========================================================================\n",
"\n",
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 1.2000E+00 seconds\n",
" Reading cross sections = 2.5000E-01 seconds\n",
" Total time in simulation = 1.8967E+01 seconds\n",
" Time in transport only = 1.8921E+01 seconds\n",
" Time in inactive batches = 2.8760E+00 seconds\n",
" Time in active batches = 1.6091E+01 seconds\n",
" Time synchronizing fission bank = 3.0000E-03 seconds\n",
" Sampling source sites = 3.0000E-03 seconds\n",
" SEND/RECV source sites = 0.0000E+00 seconds\n",
" Time accumulating tallies = 1.0000E-03 seconds\n",
" Total time for finalization = 0.0000E+00 seconds\n",
" Total time elapsed = 2.0192E+01 seconds\n",
" Calculation Rate (inactive) = 4346.31 neutrons/second\n",
" Calculation Rate (active) = 2330.50 neutrons/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
" k-effective (Collision) = 0.68196 +/- 0.00427\n",
" k-effective (Track-length) = 0.67959 +/- 0.00437\n",
" k-effective (Absorption) = 0.67957 +/- 0.00402\n",
" Combined k-effective = 0.67943 +/- 0.00295\n",
" Leakage Fraction = 0.34370 +/- 0.00201\n",
"\n"
]
},
{
"data": {
"text/plain": [
"0"
]
},
"execution_count": 20,
"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": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# We do not know how many batches were needed to satisfy the \n",
"# tally trigger(s), so find the statepoint file(s)\n",
"statepoints = glob.glob('statepoint.*.h5')\n",
"\n",
"# Load the last statepoint file\n",
"sp = StatePoint(statepoints[-1])"
]
},
{
"cell_type": "code",
"execution_count": 22,
"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-22-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": [
"**Analyze the mesh fission rate tally**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Find the mesh tally with the StatePoint API\n",
"tally = sp.get_tally(name='mesh tally')\n",
"\n",
"# Print a little info about the mesh tally to the screen\n",
"print(tally)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the new Tally data retrieval API with pure NumPy"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Get the relative error for the thermal fission reaction \n",
"# rates in the four corner pins \n",
"data = tally.get_values(scores=['fission'], filters=['mesh', 'energy'], \\\n",
" filter_bins=[((1,1),(1,17), (17,1), (17,17)), \\\n",
" ((0., 0.625e-6),)], value='rel_err')\n",
"print(data)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Get a pandas dataframe for the mesh tally data\n",
"df = tally.get_pandas_dataframe(nuclides=False)\n",
"\n",
"# Print the first twenty rows in the dataframe\n",
"df.head(20)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create a boxplot to view the distribution of\n",
"# fission and nu-fission rates in the pins\n",
"bp = df.boxplot(column='mean', by='score')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Extract thermal nu-fission rates from pandas\n",
"fiss = df[df['score'] == 'nu-fission']\n",
"fiss = fiss[fiss['energy [MeV]'] == '(0.0e+00 - 6.3e-07)']\n",
"\n",
"# Extract mean and reshape as 2D NumPy arrays\n",
"mean = fiss['mean'].reshape((17,17))\n",
"\n",
"pylab.imshow(mean, interpolation='nearest')\n",
"pylab.title('fission rate')\n",
"pylab.xlabel('x')\n",
"pylab.ylabel('y')\n",
"pylab.colorbar()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Analyze the cell+nuclides scatter-y2 rate tally**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Find the cell Tally with the StatePoint API\n",
"tally = sp.get_tally(name='cell tally')\n",
"\n",
"# Print a little info about the cell tally to the screen\n",
"print(tally)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Get a pandas dataframe for the cell tally data\n",
"df = tally.get_pandas_dataframe()\n",
"\n",
"# Print the first twenty rows in the dataframe\n",
"df.head(100)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the new Tally data retrieval API with pure NumPy"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Get the standard deviations for two of the spherical harmonic\n",
"# scattering reaction rates \n",
"data = tally.get_values(scores=['scatter-Y2,2', 'scatter-Y0,0'], \n",
" nuclides=['U-238', 'U-235'], value='std_dev')\n",
"print(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Analyze the distribcell tally**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Find the distribcell Tally with the StatePoint API\n",
"tally = sp.get_tally(name='distribcell tally')\n",
"\n",
"# Print a little info about the distribcell tally to the screen\n",
"print(tally)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Use the new Tally data retrieval API with pure NumPy"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Get the relative error for the scattering reaction rates in\n",
"# the first 30 distribcell instances \n",
"data = tally.get_values(scores=['scatter'], filters=['distribcell'],\n",
" filter_bins=[(i,) for i in range(10)], value='rel_err')\n",
"print(data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Print the distribcell tally dataframe **without** OpenCG info"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Get a pandas dataframe for the distribcell tally data\n",
"df = tally.get_pandas_dataframe(nuclides=False)\n",
"\n",
"# Print the last twenty rows in the dataframe\n",
"df.tail(20)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Print the distribcell tally dataframe **with** OpenCG info"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Get a pandas dataframe for the distribcell tally data\n",
"df = tally.get_pandas_dataframe(summary=su, nuclides=False)\n",
"\n",
"# Print the last twenty rows in the dataframe\n",
"df.head(20)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Show summary statistics for absorption distribcell tally data\n",
"absorption = df[df['score'] == 'absorption']\n",
"absorption[['mean', 'std. dev.']].dropna().describe()\n",
"\n",
"# Note that the maximum standard deviation does indeed\n",
"# meet the 5e-4 threshold set by the tally trigger"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Perform a statistical test comparing the tally sample distributions for two categories of fuel pins."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Extract tally data from pins in the pins divided along y=x diagonal \n",
"multi_index = ('level 2', 'lat',)\n",
"lower = df[df[multi_index + ('x',)] + df[multi_index + ('y',)] < 16]\n",
"upper = df[df[multi_index + ('x',)] + df[multi_index + ('y',)] > 16]\n",
"lower = lower[lower['score'] == 'absorption']\n",
"upper = upper[upper['score'] == 'absorption']\n",
"\n",
"# Perform non-parametric Mann-Whitney U Test to see if the \n",
"# absorption rates (may) come from same sampling distribution\n",
"u, p = scipy.stats.mannwhitneyu(lower['mean'], upper['mean'])\n",
"print('Mann-Whitney Test p-value: {0}'.format(p))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that the symmetry implied by the y=x diagonal ensures that the two sampling distributions are identical. Indeed, as illustrated by the test above, for any reasonable significance level (*e.g.*, $\\alpha$=0.05) one would **not reject** the null hypothesis that the two sampling distributions are identical.\n",
"\n",
"Next, perform the same test but with two groupings of pins which are not symmetrically identical to one another."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Extract tally data from pins in the pins divided along y=-x diagonal\n",
"multi_index = ('level 2', 'lat',)\n",
"lower = df[df[multi_index + ('x',)] > df[multi_index + ('y',)]]\n",
"upper = df[df[multi_index + ('x',)] < df[multi_index + ('y',)]]\n",
"lower = lower[lower['score'] == 'absorption']\n",
"upper = upper[upper['score'] == 'absorption']\n",
"\n",
"# Perform non-parametric Mann-Whitney U Test to see if the \n",
"# absorption rates (may) come from same sampling distribution\n",
"u, p = scipy.stats.mannwhitneyu(lower['mean'], upper['mean'])\n",
"print('Mann-Whitney Test p-value: {0}'.format(p))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that the asymmetry implied by the y=-x diagonal ensures that the two sampling distributions are *not* identical. Indeed, as illustrated by the test above, for any reasonable significance level (*e.g.*, $\\alpha$=0.05) one would **reject** the null hypothesis that the two sampling distributions are identical."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Extract the scatter tally data from pandas\n",
"scatter = df[df['score'] == 'scatter']\n",
"\n",
"scatter['rel. err.'] = scatter['std. dev.'] / scatter['mean']\n",
"\n",
"# Show a scatter plot of the mean vs. the std. dev.\n",
"scatter.plot(kind='scatter', x='mean', y='rel. err.', title='Scattering Rates')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Plot a histogram and kernel density estimate for the scattering rates\n",
"scatter['mean'].plot(kind='hist', bins=25)\n",
"scatter['mean'].plot(kind='kde')\n",
"pylab.title('Scattering Rates')\n",
"pylab.xlabel('Mean')\n",
"pylab.legend(['KDE', 'Histogram'])"
]
}
],
"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
}
Reference in a new issue View git blame Copy permalink