{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling a Pincell\n",
    "In this module, we'll demonstrate the basic features of the Python API for constructing input files and running OpenMC. In it, we will show how to create a basic reflective pincell model that is equivalent to modeling an infinite array of fuel pins in a pressurized water reactor. We highly recommend having a copy of the [Python API reference documentation](https://docs.openmc.org/en/stable/pythonapi/index.html) open in another browser tab that you can refer to.\n",
    "\n",
    "<img src=\"pincell.png\" alt=\"drawing\" width=\"250\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, we'll create a simple pincell that is composed of:\n",
    "- UO<sub>2</sub> with 3.5 weight% enriched in U-235 at 11 g/cm<sup>3</sup>\n",
    "- zirconium clad at 6.5 g/cm<sup>3</sup>\n",
    "- H<sub>2</sub>O moderator at 1.0 g/cm<sup>3</sup>\n",
    "\n",
    "The dimensions of our fuel pin will be as follows:\n",
    "- Fuel outer radius = 0.46955 cm\n",
    "- Clad inner radius = 0.47910 cm\n",
    "- Clad outer radius = 0.54640 cm\n",
    "- Fuel pin pitch = 1.44270 cm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import openmc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basics of Jupyter Notebook\n",
    "\n",
    "You are working within a Jupyter notebook. Some common commands which will be useful:\n",
    "- To execute a cell: `Shift+Enter`\n",
    "- To insert a cell above or below: `Esc+a`, `Esc+b`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Naming Conventions\n",
    "\n",
    "Before we start working with OpenMC's Python API, it's helpful to understand the naming convention of objects so that you can conceptualize what is a function, what is a class, etc. OpenMC's Python interface follows the same naming convention that is adopted by many/most Python projects:\n",
    "\n",
    "- Module names are `lowercase`\n",
    "- Classes are `CamelCase`\n",
    "- Functions and class attributes/variables are `lowercase_with_underscores`\n",
    "\n",
    "To give a few specific examples:\n",
    "\n",
    "- `openmc.deplete` is the depletion _module_\n",
    "- `openmc.run` is a _function_\n",
    "- `openmc.Material` is a _class_\n",
    "- `openmc.StatePoint` is a _class_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting Attributes\n",
    "\n",
    "When building OpenMC models, we will work with many different classes. Each class typically has _attributes_, a variable belonging to the class. When creating a class, you can often set those attributes directly when you instantiate the object, like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_cell = openmc.Cell(name=\"box\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or, you can assign values to attributes after you have already created the object. For example, the following is equivalent to the above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_other_cell = openmc.Cell()\n",
    "my_other_cell.name = \"box2\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can generally also mix-and-match, setting some attributes when you instantiate the object, and others at a later point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_cell = openmc.Cell(name=\"box3\")\n",
    "new_cell.temperature = 500"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## How to Get Help/Learn More\n",
    "\n",
    "When building OpenMC models, we recommend having a copy of the Python API documentation (linked earlier) open. You can also query information about classes and methods directly from Jupyter. For example, we can use `help` to get documentation on all of the valid attributes for functions and classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on class Cell in module openmc.cell:\n",
      "\n",
      "class Cell(openmc.mixin.IDManagerMixin)\n",
      " |  Cell(cell_id=None, name='', fill=None, region=None)\n",
      " |  \n",
      " |  A region of space defined as the intersection of half-space created by\n",
      " |  quadric surfaces.\n",
      " |  \n",
      " |  Parameters\n",
      " |  ----------\n",
      " |  cell_id : int, optional\n",
      " |      Unique identifier for the cell. If not specified, an identifier will\n",
      " |      automatically be assigned.\n",
      " |  name : str, optional\n",
      " |      Name of the cell. If not specified, the name is the empty string.\n",
      " |  fill : openmc.Material or openmc.UniverseBase or openmc.Lattice or None or iterable of openmc.Material, optional\n",
      " |      Indicates what the region of space is filled with\n",
      " |  region : openmc.Region, optional\n",
      " |      Region of space that is assigned to the cell.\n",
      " |  \n",
      " |  Attributes\n",
      " |  ----------\n",
      " |  id : int\n",
      " |      Unique identifier for the cell\n",
      " |  name : str\n",
      " |      Name of the cell\n",
      " |  fill : openmc.Material or openmc.UniverseBase or openmc.Lattice or None or iterable of openmc.Material\n",
      " |      Indicates what the region of space is filled with. If None, the cell is\n",
      " |      treated as a void. An iterable of materials is used to fill repeated\n",
      " |      instances of a cell with different materials.\n",
      " |  fill_type : {'material', 'universe', 'lattice', 'distribmat', 'void'}\n",
      " |      Indicates what the cell is filled with.\n",
      " |  region : openmc.Region or None\n",
      " |      Region of space that is assigned to the cell.\n",
      " |  rotation : Iterable of float\n",
      " |      If the cell is filled with a universe, this array specifies the angles\n",
      " |      in degrees about the x, y, and z axes that the filled universe should be\n",
      " |      rotated. The rotation applied is an intrinsic rotation with specified\n",
      " |      Tait-Bryan angles. That is to say, if the angles are :math:`(\\phi,\n",
      " |      \\theta, \\psi)`, then the rotation matrix applied is :math:`R_z(\\psi)\n",
      " |      R_y(\\theta) R_x(\\phi)` or\n",
      " |  \n",
      " |      .. math::\n",
      " |  \n",
      " |         \\left [ \\begin{array}{ccc} \\cos\\theta \\cos\\psi & -\\cos\\phi \\sin\\psi\n",
      " |         + \\sin\\phi \\sin\\theta \\cos\\psi & \\sin\\phi \\sin\\psi + \\cos\\phi\n",
      " |         \\sin\\theta \\cos\\psi \\\\ \\cos\\theta \\sin\\psi & \\cos\\phi \\cos\\psi +\n",
      " |         \\sin\\phi \\sin\\theta \\sin\\psi & -\\sin\\phi \\cos\\psi + \\cos\\phi\n",
      " |         \\sin\\theta \\sin\\psi \\\\ -\\sin\\theta & \\sin\\phi \\cos\\theta & \\cos\\phi\n",
      " |         \\cos\\theta \\end{array} \\right ]\n",
      " |  \n",
      " |      A rotation matrix can also be specified directly by setting this\n",
      " |      attribute to a nested list (or 2D numpy array) that specifies each\n",
      " |      element of the matrix.\n",
      " |  rotation_matrix : numpy.ndarray\n",
      " |      The rotation matrix defined by the angles specified in the\n",
      " |      :attr:`Cell.rotation` property.\n",
      " |  temperature : float or iterable of float\n",
      " |      Temperature of the cell in Kelvin.  Multiple temperatures can be given\n",
      " |      to give each distributed cell instance a unique temperature.\n",
      " |  translation : Iterable of float\n",
      " |      If the cell is filled with a universe, this array specifies a vector\n",
      " |      that is used to translate (shift) the universe.\n",
      " |  paths : list of str\n",
      " |      The paths traversed through the CSG tree to reach each cell\n",
      " |      instance. This property is initialized by calling the\n",
      " |      :meth:`Geometry.determine_paths` method.\n",
      " |  num_instances : int\n",
      " |      The number of instances of this cell throughout the geometry.\n",
      " |  volume : float\n",
      " |      Volume of the cell in cm^3. This can either be set manually or\n",
      " |      calculated in a stochastic volume calculation and added via the\n",
      " |      :meth:`Cell.add_volume_information` method. For 'distribmat' cells\n",
      " |      it is the total volume of all instances.\n",
      " |  atoms : dict\n",
      " |      Mapping of nuclides to the total number of atoms for each nuclide\n",
      " |      present in the cell, or in all of its instances for a 'distribmat'\n",
      " |      fill. For example, {'U235': 1.0e22, 'U238': 5.0e22, ...}.\n",
      " |  \n",
      " |      .. versionadded:: 0.12\n",
      " |  bounding_box : openmc.BoundingBox\n",
      " |      Axis-aligned bounding box of the cell\n",
      " |  \n",
      " |  Method resolution order:\n",
      " |      Cell\n",
      " |      openmc.mixin.IDManagerMixin\n",
      " |      builtins.object\n",
      " |  \n",
      " |  Methods defined here:\n",
      " |  \n",
      " |  __contains__(self, point)\n",
      " |  \n",
      " |  __init__(self, cell_id=None, name='', fill=None, region=None)\n",
      " |      Initialize self.  See help(type(self)) for accurate signature.\n",
      " |  \n",
      " |  __repr__(self)\n",
      " |      Return repr(self).\n",
      " |  \n",
      " |  add_volume_information(self, volume_calc)\n",
      " |      Add volume information to a cell.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      volume_calc : openmc.VolumeCalculation\n",
      " |          Results from a stochastic volume calculation\n",
      " |  \n",
      " |  clone(self, clone_materials=True, clone_regions=True, memo=None)\n",
      " |      Create a copy of this cell with a new unique ID, and clones\n",
      " |      the cell's region and fill.\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      clone_materials : bool\n",
      " |          Whether to create separate copies of the materials filling cells\n",
      " |          contained in this cell, or the material filling this cell.\n",
      " |      clone_regions : bool\n",
      " |          Whether to create separate copies of the regions bounding cells\n",
      " |          contained in this cell, and the region bounding this cell.\n",
      " |      memo : dict or None\n",
      " |          A nested dictionary of previously cloned objects. This parameter\n",
      " |          is used internally and should not be specified by the user.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      clone : openmc.Cell\n",
      " |          The clone of this cell\n",
      " |  \n",
      " |  create_xml_subelement(self, xml_element, memo=None)\n",
      " |      Add the cell's xml representation to an incoming xml element\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      xml_element : lxml.etree._Element\n",
      " |          XML element to be added to\n",
      " |      \n",
      " |      memo : set or None\n",
      " |          A set of object IDs representing geometry entities already\n",
      " |          written to ``xml_element``. This parameter is used internally\n",
      " |          and should not be specified by users.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      None\n",
      " |  \n",
      " |  get_all_cells(self, memo=None)\n",
      " |      Return all cells that are contained within this one if it is filled with a\n",
      " |      universe or lattice\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      cells : dict\n",
      " |          Dictionary whose keys are cell IDs and values are :class:`Cell`\n",
      " |          instances\n",
      " |  \n",
      " |  get_all_materials(self, memo=None)\n",
      " |      Return all materials that are contained within the cell\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      materials : dict\n",
      " |          Dictionary whose keys are material IDs and values are\n",
      " |          :class:`Material` instances\n",
      " |  \n",
      " |  get_all_universes(self)\n",
      " |      Return all universes that are contained within this one if any of\n",
      " |      its cells are filled with a universe or lattice.\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      universes : dict\n",
      " |          Dictionary whose keys are universe IDs and values are\n",
      " |          :class:`Universe` instances\n",
      " |  \n",
      " |  get_nuclide_densities(self)\n",
      " |      Return all nuclides contained in the cell and their densities\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      nuclides : dict\n",
      " |          Dictionary whose keys are nuclide names and values are 2-tuples of\n",
      " |          (nuclide, density)\n",
      " |  \n",
      " |  get_nuclides(self)\n",
      " |      Returns all nuclides in the cell\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      nuclides : list of str\n",
      " |          List of nuclide names\n",
      " |  \n",
      " |  plot(self, *args, **kwargs)\n",
      " |      Display a slice plot of the cell.\n",
      " |      \n",
      " |      .. versionadded:: 0.14.0\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      origin : iterable of float\n",
      " |          Coordinates at the origin of the plot. If left as None then the\n",
      " |          bounding box center will be used to attempt to ascertain the origin.\n",
      " |          Defaults to (0, 0, 0) if the bounding box is not finite\n",
      " |      width : iterable of float\n",
      " |          Width of the plot in each basis direction. If left as none then the\n",
      " |          bounding box width will be used to attempt to ascertain the plot\n",
      " |          width. Defaults to (10, 10) if the bounding box is not finite\n",
      " |      pixels : Iterable of int or int\n",
      " |          If iterable of ints provided, then this directly sets the number of\n",
      " |          pixels to use in each basis direction. If int provided, then this\n",
      " |          sets the total number of pixels in the plot and the number of pixels\n",
      " |          in each basis direction is calculated from this total and the image\n",
      " |          aspect ratio.\n",
      " |      basis : {'xy', 'xz', 'yz'}\n",
      " |          The basis directions for the plot\n",
      " |      color_by : {'cell', 'material'}\n",
      " |          Indicate whether the plot should be colored by cell or by material\n",
      " |      colors : dict\n",
      " |          Assigns colors to specific materials or cells. Keys are instances of\n",
      " |          :class:`Cell` or :class:`Material` and values are RGB 3-tuples, RGBA\n",
      " |          4-tuples, or strings indicating SVG color names. Red, green, blue,\n",
      " |          and alpha should all be floats in the range [0.0, 1.0], for example:\n",
      " |      \n",
      " |          .. code-block:: python\n",
      " |      \n",
      " |             # Make water blue\n",
      " |             water = openmc.Cell(fill=h2o)\n",
      " |             universe.plot(..., colors={water: (0., 0., 1.))\n",
      " |      seed : int\n",
      " |          Seed for the random number generator\n",
      " |      openmc_exec : str\n",
      " |          Path to OpenMC executable.\n",
      " |      axes : matplotlib.Axes\n",
      " |          Axes to draw to\n",
      " |      legend : bool\n",
      " |          Whether a legend showing material or cell names should be drawn\n",
      " |      legend_kwargs : dict\n",
      " |          Keyword arguments passed to :func:`matplotlib.pyplot.legend`.\n",
      " |      outline : bool\n",
      " |          Whether outlines between color boundaries should be drawn\n",
      " |      axis_units : {'km', 'm', 'cm', 'mm'}\n",
      " |          Units used on the plot axis\n",
      " |      **kwargs\n",
      " |          Keyword arguments passed to :func:`matplotlib.pyplot.imshow`\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      matplotlib.axes.Axes\n",
      " |          Axes containing resulting image\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Class methods defined here:\n",
      " |  \n",
      " |  from_xml_element(elem, surfaces, materials, get_universe) from builtins.type\n",
      " |      Generate cell from XML element\n",
      " |      \n",
      " |      Parameters\n",
      " |      ----------\n",
      " |      elem : lxml.etree._Element\n",
      " |          `<cell>` element\n",
      " |      surfaces : dict\n",
      " |          Dictionary mapping surface IDs to :class:`openmc.Surface` instances\n",
      " |      materials : dict\n",
      " |          Dictionary mapping material ID strings to :class:`openmc.Material`\n",
      " |          instances (defined in :meth:`openmc.Geometry.from_xml`)\n",
      " |      get_universe : function\n",
      " |          Function returning universe (defined in\n",
      " |          :meth:`openmc.Geometry.from_xml`)\n",
      " |      \n",
      " |      Returns\n",
      " |      -------\n",
      " |      openmc.Cell\n",
      " |          Cell instance\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Readonly properties defined here:\n",
      " |  \n",
      " |  atoms\n",
      " |  \n",
      " |  bounding_box\n",
      " |  \n",
      " |  fill_type\n",
      " |  \n",
      " |  num_instances\n",
      " |  \n",
      " |  paths\n",
      " |  \n",
      " |  rotation_matrix\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors defined here:\n",
      " |  \n",
      " |  fill\n",
      " |  \n",
      " |  name\n",
      " |  \n",
      " |  region\n",
      " |  \n",
      " |  rotation\n",
      " |  \n",
      " |  temperature\n",
      " |  \n",
      " |  translation\n",
      " |  \n",
      " |  volume\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data and other attributes defined here:\n",
      " |  \n",
      " |  next_id = 3\n",
      " |  \n",
      " |  used_ids = {1, 2, 3}\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors inherited from openmc.mixin.IDManagerMixin:\n",
      " |  \n",
      " |  __dict__\n",
      " |      dictionary for instance variables (if defined)\n",
      " |  \n",
      " |  __weakref__\n",
      " |      list of weak references to the object (if defined)\n",
      " |  \n",
      " |  id\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(openmc.Cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Attributes which have a default value will appear as `<attribute>=<the default value>` in the class/function definition. For example, in the above we see that the default name for a cell is the empty string, `''`.\n",
    "\n",
    "For classes, all of the member functions you can access will then be defined further down, with their associated documentation. To see documentation for just one specific function, you can use syntax like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on method clone in module openmc.cell:\n",
      "\n",
      "clone(clone_materials=True, clone_regions=True, memo=None) method of openmc.cell.Cell instance\n",
      "    Create a copy of this cell with a new unique ID, and clones\n",
      "    the cell's region and fill.\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    clone_materials : bool\n",
      "        Whether to create separate copies of the materials filling cells\n",
      "        contained in this cell, or the material filling this cell.\n",
      "    clone_regions : bool\n",
      "        Whether to create separate copies of the regions bounding cells\n",
      "        contained in this cell, and the region bounding this cell.\n",
      "    memo : dict or None\n",
      "        A nested dictionary of previously cloned objects. This parameter\n",
      "        is used internally and should not be specified by the user.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    clone : openmc.Cell\n",
      "        The clone of this cell\n",
      "\n"
     ]
    }
   ],
   "source": [
    "my_cell = openmc.Cell()\n",
    "help(my_cell.clone)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also query the type of a particular object using `type`. Note that we can also display the attributes for a particular object using `print`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "openmc.cell.Cell"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(my_cell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cell\n",
      "\tID             =\t4\n",
      "\tName           =\t\n",
      "\tFill           =\tNone\n",
      "\tRegion         =\tNone\n",
      "\tRotation       =\tNone\n",
      "\tTranslation    =\tNone\n",
      "\tVolume         =\tNone\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(my_cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to quickly see all of the member functions on a class, press `Tab` to perform tab-completion to view a drop-down list on all the options available to you. You can then do `Shift+Tab` to see the documentation for each function inline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#my_cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OpenMC Model\n",
    "\n",
    "The OpenMC `Model` class houses all of the pieces of a Monte Carlo simulation. We will assemble our pincell by progressively adding to a model, and then run that model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = openmc.Model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The `cross_sections.xml` file\n",
    "\n",
    "The `cross_sections.xml` tells OpenMC where it can find nuclide cross sections and $S(\\alpha,\\beta)$ tables. It serves the same purpose as MCNP's `xsdir` file and Serpent's `xsdata` file. As we mentioned, this can be set either by the `OPENMC_CROSS_SECTIONS` environment variable or the `Materials.cross_sections` attribute.\n",
    "\n",
    "Let's have a look at what's inside this file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<?xml version='1.0' encoding='utf-8'?>\n",
      "<cross_sections>\n",
      "  <library materials=\"H1\" path=\"neutron/H1.h5\" type=\"neutron\" />\n",
      "  <library materials=\"H2\" path=\"neutron/H2.h5\" type=\"neutron\" />\n",
      "  <library materials=\"H3\" path=\"neutron/H3.h5\" type=\"neutron\" />\n",
      "  <library materials=\"He3\" path=\"neutron/He3.h5\" type=\"neutron\" />\n",
      "  <library materials=\"He4\" path=\"neutron/He4.h5\" type=\"neutron\" />\n",
      "  <library materials=\"Li6\" path=\"neutron/Li6.h5\" type=\"neutron\" />\n",
      "  <library materials=\"Li7\" path=\"neutron/Li7.h5\" type=\"neutron\" />\n",
      "  <library materials=\"Be7\" path=\"neutron/Be7.h5\" type=\"neutron\" />\n",
      "    ...\n",
      "  <library materials=\"Cf253\" path=\"wmp/098253.h5\" type=\"wmp\" />\n",
      "  <library materials=\"Cf254\" path=\"wmp/098254.h5\" type=\"wmp\" />\n",
      "  <library materials=\"Es251\" path=\"wmp/099251.h5\" type=\"wmp\" />\n",
      "  <library materials=\"Es252\" path=\"wmp/099252.h5\" type=\"wmp\" />\n",
      "  <library materials=\"Es253\" path=\"wmp/099253.h5\" type=\"wmp\" />\n",
      "  <library materials=\"Es254\" path=\"wmp/099254.h5\" type=\"wmp\" />\n",
      "  <library materials=\"Es254_m1\" path=\"wmp/099254m1.h5\" type=\"wmp\" />\n",
      "  <library materials=\"Es255\" path=\"wmp/099255.h5\" type=\"wmp\" />\n",
      "  <library materials=\"Fm255\" path=\"wmp/100255.h5\" type=\"wmp\" />\n",
      "</cross_sections>\n"
     ]
    }
   ],
   "source": [
    "!cat $OPENMC_CROSS_SECTIONS | head -n 10\n",
    "print('    ...')\n",
    "!cat $OPENMC_CROSS_SECTIONS | tail -n 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining Materials\n",
    "\n",
    "Materials in OpenMC are defined as a set of nuclides with specified atom/weight fractions. To begin, we will create a material by making an instance of the `Material` class. In OpenMC, many objects, including materials, are identified by a \"unique ID\" (a positive integer). You can also give a material a `name` as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Material\n",
      "\tID             =\t1\n",
      "\tName           =\tzirconium\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\tNone [sum]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tNuclides       \n",
      "\n"
     ]
    }
   ],
   "source": [
    "zirconium = openmc.Material(1, \"zirconium\")\n",
    "print(zirconium)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One great feature in OpenMC's design is that the user interaction is designed to be as streamlined as possible. You don't *need* to set either an ID or a name (because often you can simply refer to an object with the actual object handle). If you were to create a material without any ID, OpenMC will assign a default for you."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Material\n",
      "\tID             =\t2\n",
      "\tName           =\t\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\tNone [sum]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tNuclides       \n",
      "\n"
     ]
    }
   ],
   "source": [
    "mat = openmc.Material()\n",
    "print(mat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that an ID of 2 was automatically assigned. Let's now move on to adding nuclides to our material. The `Material` object has a method `add_element()` whose first argument is the name of the nuclide and second argument is the atom or weight fraction. We see that by default it assumes we want an atom fraction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "#zirconium.add_element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "zirconium.add_element('Zr', 1.0)\n",
    "zirconium.set_density('g/cm3', 6.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Material\n",
      "\tID             =\t1\n",
      "\tName           =\tzirconium\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t6.5 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tNuclides       \n",
      "\tZr90           =\t0.5145       [ao]\n",
      "\tZr91           =\t0.1122       [ao]\n",
      "\tZr92           =\t0.1715       [ao]\n",
      "\tZr94           =\t0.1738       [ao]\n",
      "\tZr96           =\t0.028        [ao]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(zirconium)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Element Expansion\n",
    "We can see that OpenMC automatically inserted the natural abundances of the zirconium isotopes for us! How convenient! The way this feature works is as follows:\n",
    "\n",
    "- First, it checks whether `Materials.cross_sections` has been set, indicating the path to a `cross_sections.xml` file.\n",
    "- If `Materials.cross_sections` isn't set, it looks for the `OPENMC_CROSS_SECTIONS` environment variable.\n",
    "- If either of these are found, it scans the file to see what nuclides are actually available and will expand elements accordingly.\n",
    "\n",
    "Let's build our fuel material. For sake of illustration, let's suppose that we had O-16, but not natural oxygen in our fuel. Also note that OpenMC has a convenient feature to set the nuclide concentrations for weight percent enrichments in U-235."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add nuclides to uo2\n",
    "uo2 = openmc.Material(name=\"uo2\")\n",
    "uo2.add_element('U', 1.0, enrichment=3.5)\n",
    "uo2.add_nuclide('O16', 2.0)\n",
    "uo2.set_density('g/cm3', 10.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With UO2 finished, let's now create materials for the coolant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "water = openmc.Material(name=\"water\")\n",
    "water.add_element('H', 2.0)\n",
    "water.add_nuclide('O16', 1.0)\n",
    "water.set_density('g/cm3', 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An astute observer might now point out that this water material we just created will only use free-atom cross sections. We need to tell it to use an $S(\\alpha,\\beta)$ table so that the bound atom cross section is used at thermal energies. To do this, there's an `add_s_alpha_beta()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "water.add_s_alpha_beta('c_H_in_H2O')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/cross_sections.xml\n",
      "cross_sections.xml \u001b[1m\u001b[36mneutron\u001b[m\u001b[m            \u001b[1m\u001b[36mphoton\u001b[m\u001b[m             \u001b[1m\u001b[36mwmp\u001b[m\u001b[m\n",
      "Ac225.h5             Ge76.h5              Sb121.h5\n",
      "Ac226.h5             H1.h5                Sb123.h5\n",
      "Ac227.h5             H2.h5                Sb124.h5\n",
      "Ag107.h5             H3.h5                Sb125.h5\n",
      "Ag109.h5             He3.h5               Sb126.h5\n",
      "Ag110_m1.h5          He4.h5               Sc45.h5\n",
      "Ag111.h5             Hf174.h5             Se74.h5\n",
      "Al27.h5              Hf176.h5             Se76.h5\n",
      "Am240.h5             Hf177.h5             Se77.h5\n",
      "Am241.h5             Hf178.h5             Se78.h5\n",
      "Am242.h5             Hf179.h5             Se79.h5\n",
      "Am242_m1.h5          Hf180.h5             Se80.h5\n",
      "Am243.h5             Hg196.h5             Se82.h5\n",
      "Am244.h5             Hg198.h5             Si28.h5\n",
      "Am244_m1.h5          Hg199.h5             Si29.h5\n",
      "Ar36.h5              Hg200.h5             Si30.h5\n",
      "Ar38.h5              Hg201.h5             Sm144.h5\n",
      "Ar40.h5              Hg202.h5             Sm147.h5\n",
      "As74.h5              Hg204.h5             Sm148.h5\n",
      "As75.h5              Ho165.h5             Sm149.h5\n",
      "Au197.h5             Ho166_m1.h5          Sm150.h5\n",
      "B10.h5               I127.h5              Sm151.h5\n",
      "B11.h5               I129.h5              Sm152.h5\n",
      "Ba130.h5             I130.h5              Sm153.h5\n",
      "Ba132.h5             I131.h5              Sm154.h5\n",
      "Ba133.h5             I135.h5              Sn112.h5\n",
      "Ba134.h5             In113.h5             Sn113.h5\n",
      "Ba135.h5             In115.h5             Sn114.h5\n",
      "Ba136.h5             Ir191.h5             Sn115.h5\n",
      "Ba137.h5             Ir193.h5             Sn116.h5\n",
      "Ba138.h5             K39.h5               Sn117.h5\n",
      "Ba140.h5             K40.h5               Sn118.h5\n",
      "Be7.h5               K41.h5               Sn119.h5\n",
      "Be9.h5               Kr78.h5              Sn120.h5\n",
      "Bi209.h5             Kr80.h5              Sn122.h5\n",
      "Bk245.h5             Kr82.h5              Sn123.h5\n",
      "Bk246.h5             Kr83.h5              Sn124.h5\n",
      "Bk247.h5             Kr84.h5              Sn125.h5\n",
      "Bk248.h5             Kr85.h5              Sn126.h5\n",
      "Bk249.h5             Kr86.h5              Sr84.h5\n",
      "Bk250.h5             La138.h5             Sr86.h5\n",
      "Br79.h5              La139.h5             Sr87.h5\n",
      "Br81.h5              La140.h5             Sr88.h5\n",
      "C0.h5                Li6.h5               Sr89.h5\n",
      "Ca40.h5              Li7.h5               Sr90.h5\n",
      "Ca42.h5              Lu175.h5             Ta180.h5\n",
      "Ca43.h5              Lu176.h5             Ta181.h5\n",
      "Ca44.h5              Mg24.h5              Ta182.h5\n",
      "Ca46.h5              Mg25.h5              Tb159.h5\n",
      "Ca48.h5              Mg26.h5              Tb160.h5\n",
      "Cd106.h5             Mn55.h5              Tc99.h5\n",
      "Cd108.h5             Mo100.h5             Te120.h5\n",
      "Cd110.h5             Mo92.h5              Te122.h5\n",
      "Cd111.h5             Mo94.h5              Te123.h5\n",
      "Cd112.h5             Mo95.h5              Te124.h5\n",
      "Cd113.h5             Mo96.h5              Te125.h5\n",
      "Cd114.h5             Mo97.h5              Te126.h5\n",
      "Cd115_m1.h5          Mo98.h5              Te127_m1.h5\n",
      "Cd116.h5             Mo99.h5              Te128.h5\n",
      "Ce136.h5             N14.h5               Te129_m1.h5\n",
      "Ce138.h5             N15.h5               Te130.h5\n",
      "Ce139.h5             Na22.h5              Te132.h5\n",
      "Ce140.h5             Na23.h5              Th227.h5\n",
      "Ce141.h5             Nb93.h5              Th228.h5\n",
      "Ce142.h5             Nb94.h5              Th229.h5\n",
      "Ce143.h5             Nb95.h5              Th230.h5\n",
      "Ce144.h5             Nd142.h5             Th231.h5\n",
      "Cf246.h5             Nd143.h5             Th232.h5\n",
      "Cf248.h5             Nd144.h5             Th233.h5\n",
      "Cf249.h5             Nd145.h5             Th234.h5\n",
      "Cf250.h5             Nd146.h5             Ti46.h5\n",
      "Cf251.h5             Nd147.h5             Ti47.h5\n",
      "Cf252.h5             Nd148.h5             Ti48.h5\n",
      "Cf253.h5             Nd150.h5             Ti49.h5\n",
      "Cf254.h5             Ni58.h5              Ti50.h5\n",
      "Cl35.h5              Ni59.h5              Tl203.h5\n",
      "Cl37.h5              Ni60.h5              Tl205.h5\n",
      "Cm240.h5             Ni61.h5              Tm168.h5\n",
      "Cm241.h5             Ni62.h5              Tm169.h5\n",
      "Cm242.h5             Ni64.h5              Tm170.h5\n",
      "Cm243.h5             Np234.h5             U230.h5\n",
      "Cm244.h5             Np235.h5             U231.h5\n",
      "Cm245.h5             Np236.h5             U232.h5\n",
      "Cm246.h5             Np237.h5             U233.h5\n",
      "Cm247.h5             Np238.h5             U234.h5\n",
      "Cm248.h5             Np239.h5             U235.h5\n",
      "Cm249.h5             O16.h5               U236.h5\n",
      "Cm250.h5             O17.h5               U237.h5\n",
      "Co58.h5              P31.h5               U238.h5\n",
      "Co58_m1.h5           Pa229.h5             U239.h5\n",
      "Co59.h5              Pa230.h5             U240.h5\n",
      "Cr50.h5              Pa231.h5             U241.h5\n",
      "Cr52.h5              Pa232.h5             V50.h5\n",
      "Cr53.h5              Pa233.h5             V51.h5\n",
      "Cr54.h5              Pb204.h5             W180.h5\n",
      "Cs133.h5             Pb206.h5             W182.h5\n",
      "Cs134.h5             Pb207.h5             W183.h5\n",
      "Cs135.h5             Pb208.h5             W184.h5\n",
      "Cs136.h5             Pd102.h5             W186.h5\n",
      "Cs137.h5             Pd104.h5             Xe123.h5\n",
      "Cu63.h5              Pd105.h5             Xe124.h5\n",
      "Cu65.h5              Pd106.h5             Xe126.h5\n",
      "Dy156.h5             Pd107.h5             Xe128.h5\n",
      "Dy158.h5             Pd108.h5             Xe129.h5\n",
      "Dy160.h5             Pd110.h5             Xe130.h5\n",
      "Dy161.h5             Pm147.h5             Xe131.h5\n",
      "Dy162.h5             Pm148.h5             Xe132.h5\n",
      "Dy163.h5             Pm148_m1.h5          Xe133.h5\n",
      "Dy164.h5             Pm149.h5             Xe134.h5\n",
      "Er162.h5             Pm151.h5             Xe135.h5\n",
      "Er164.h5             Pr141.h5             Xe136.h5\n",
      "Er166.h5             Pr142.h5             Y89.h5\n",
      "Er167.h5             Pr143.h5             Y90.h5\n",
      "Er168.h5             Pu236.h5             Y91.h5\n",
      "Er170.h5             Pu237.h5             Zn64.h5\n",
      "Es251.h5             Pu238.h5             Zn65.h5\n",
      "Es252.h5             Pu239.h5             Zn66.h5\n",
      "Es253.h5             Pu240.h5             Zn67.h5\n",
      "Es254.h5             Pu241.h5             Zn68.h5\n",
      "Es254_m1.h5          Pu242.h5             Zn70.h5\n",
      "Es255.h5             Pu243.h5             Zr90.h5\n",
      "Eu151.h5             Pu244.h5             Zr91.h5\n",
      "Eu152.h5             Pu246.h5             Zr92.h5\n",
      "Eu153.h5             Ra223.h5             Zr93.h5\n",
      "Eu154.h5             Ra224.h5             Zr94.h5\n",
      "Eu155.h5             Ra225.h5             Zr95.h5\n",
      "Eu156.h5             Ra226.h5             Zr96.h5\n",
      "Eu157.h5             Rb85.h5              c_Al27.h5\n",
      "F19.h5               Rb86.h5              c_Be.h5\n",
      "Fe54.h5              Rb87.h5              c_Be_in_BeO.h5\n",
      "Fe56.h5              Re185.h5             c_C6H6.h5\n",
      "Fe57.h5              Re187.h5             c_D_in_D2O.h5\n",
      "Fe58.h5              Rh103.h5             c_Fe56.h5\n",
      "Fm255.h5             Rh105.h5             c_Graphite.h5\n",
      "Ga69.h5              Ru100.h5             c_H_in_CH2.h5\n",
      "Ga71.h5              Ru101.h5             c_H_in_CH4_liquid.h5\n",
      "Gd152.h5             Ru102.h5             c_H_in_CH4_solid.h5\n",
      "Gd153.h5             Ru103.h5             c_H_in_H2O.h5\n",
      "Gd154.h5             Ru104.h5             c_H_in_ZrH.h5\n",
      "Gd155.h5             Ru105.h5             c_O_in_BeO.h5\n",
      "Gd156.h5             Ru106.h5             c_O_in_UO2.h5\n",
      "Gd157.h5             Ru96.h5              c_SiO2_alpha.h5\n",
      "Gd158.h5             Ru98.h5              c_U_in_UO2.h5\n",
      "Gd160.h5             Ru99.h5              c_Zr_in_ZrH.h5\n",
      "Ge70.h5              S32.h5               c_ortho_D.h5\n",
      "Ge72.h5              S33.h5               c_ortho_H.h5\n",
      "Ge73.h5              S34.h5               c_para_D.h5\n",
      "Ge74.h5              S36.h5               c_para_H.h5\n"
     ]
    }
   ],
   "source": [
    "!ls $OPENMC_CROSS_SECTIONS\n",
    "!ls /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/\n",
    "!ls /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are done with our materials -- now we just need to register them in our model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Material\n",
      "\tID             =\t3\n",
      "\tName           =\tuo2\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t10.0 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tTrue\n",
      "\tS(a,b) Tables  \n",
      "\tNuclides       \n",
      "\tU234           =\t0.0003166930253944235 [ao]\n",
      "\tU235           =\t0.03543164439454172 [ao]\n",
      "\tU238           =\t0.964089368630351 [ao]\n",
      "\tU236           =\t0.00016229394971280895 [ao]\n",
      "\tO16            =\t2.0          [ao]\n",
      ", Material\n",
      "\tID             =\t1\n",
      "\tName           =\tzirconium\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t6.5 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tNuclides       \n",
      "\tZr90           =\t0.5145       [ao]\n",
      "\tZr91           =\t0.1122       [ao]\n",
      "\tZr92           =\t0.1715       [ao]\n",
      "\tZr94           =\t0.1738       [ao]\n",
      "\tZr96           =\t0.028        [ao]\n",
      ", Material\n",
      "\tID             =\t4\n",
      "\tName           =\twater\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t1.0 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tS(a,b)         =\t('c_H_in_H2O', 1.0)\n",
      "\tNuclides       \n",
      "\tH1             =\t1.99968852   [ao]\n",
      "\tH2             =\t0.00031148   [ao]\n",
      "\tO16            =\t1.0          [ao]\n",
      "]\n"
     ]
    }
   ],
   "source": [
    "model.materials = openmc.Materials([uo2, zirconium, water])\n",
    "print(model.materials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Material Modifications\n",
    "\n",
    "The power of OpenMC's Python API is the ability to easily modify and control your model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Material\n",
      "\tID             =\t4\n",
      "\tName           =\twater\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t1.0 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tS(a,b)         =\t('c_H_in_H2O', 1.0)\n",
      "\tNuclides       \n",
      "\tH1             =\t1.99968852   [ao]\n",
      "\tH2             =\t0.00031148   [ao]\n",
      "\tO16            =\t1.0          [ao]\n",
      "\n",
      "Material\n",
      "\tID             =\t4\n",
      "\tName           =\twater\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t1.0 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tS(a,b)         =\t('c_H_in_H2O', 1.0)\n",
      "\tNuclides       \n",
      "\tH1             =\t1.99968852   [ao]\n",
      "\tH2             =\t0.00031148   [ao]\n",
      "\n",
      "Material\n",
      "\tID             =\t4\n",
      "\tName           =\twater\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t1.0 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tS(a,b)         =\t('c_H_in_H2O', 1.0)\n",
      "\tNuclides       \n",
      "\tH1             =\t1.99968852   [ao]\n",
      "\tH2             =\t0.00031148   [ao]\n",
      "\tO16            =\t0.999621     [ao]\n",
      "\tO17            =\t0.000379     [ao]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(water)\n",
    "\n",
    "water.remove_nuclide('O16')\n",
    "\n",
    "print(water)\n",
    "\n",
    "water.add_element('O', 1.0)\n",
    "\n",
    "print(water)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that now O16 and O17 were automatically added. O18 is missing because our cross sections file (which is based on ENDF/B-VII.1) doesn't have O18."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Python Lists\n",
    "\n",
    "Many of OpenMC's classes are derived from Python's `list` type, and can therefore be appended, popped, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "print(isinstance(model.materials, list))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "new_fuel = openmc.Material(name=\"new_fuel\")\n",
    "new_fuel.add_nuclide('Pu239', 1.0)\n",
    "new_fuel.add_nuclide('Si28', 2.0)\n",
    "new_fuel.set_density('g/cm3', 9.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Material\n",
      "\tID             =\t3\n",
      "\tName           =\tuo2\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t10.0 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tTrue\n",
      "\tS(a,b) Tables  \n",
      "\tNuclides       \n",
      "\tU234           =\t0.0003166930253944235 [ao]\n",
      "\tU235           =\t0.03543164439454172 [ao]\n",
      "\tU238           =\t0.964089368630351 [ao]\n",
      "\tU236           =\t0.00016229394971280895 [ao]\n",
      "\tO16            =\t2.0          [ao]\n",
      ", Material\n",
      "\tID             =\t1\n",
      "\tName           =\tzirconium\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t6.5 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tNuclides       \n",
      "\tZr90           =\t0.5145       [ao]\n",
      "\tZr91           =\t0.1122       [ao]\n",
      "\tZr92           =\t0.1715       [ao]\n",
      "\tZr94           =\t0.1738       [ao]\n",
      "\tZr96           =\t0.028        [ao]\n",
      ", Material\n",
      "\tID             =\t4\n",
      "\tName           =\twater\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t1.0 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tS(a,b)         =\t('c_H_in_H2O', 1.0)\n",
      "\tNuclides       \n",
      "\tH1             =\t1.99968852   [ao]\n",
      "\tH2             =\t0.00031148   [ao]\n",
      "\tO16            =\t0.999621     [ao]\n",
      "\tO17            =\t0.000379     [ao]\n",
      ", Material\n",
      "\tID             =\t5\n",
      "\tName           =\tnew_fuel\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t9.0 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tTrue\n",
      "\tS(a,b) Tables  \n",
      "\tNuclides       \n",
      "\tPu239          =\t1.0          [ao]\n",
      "\tSi28           =\t2.0          [ao]\n",
      "]\n"
     ]
    }
   ],
   "source": [
    "model.materials.append(new_fuel)\n",
    "print(model.materials)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining Geometry\n",
    "\n",
    "We now need to define the geometry. One way to do so is to use constructive solid geometry (CSG), also known as combinatorial geometry. The object that allows us to assign a material to a region of space is called a `Cell` (same concept in MCNP, for those familiar). There are four stages in building a cell:\n",
    "\n",
    "#### Surfaces\n",
    "In order to define a region that we can assign to a cell, we must first define surfaces which bound the region. A *surface* is a locus of zeros of a function of Cartesian coordinates $x$, $y$, and $z$, e.g.\n",
    "\n",
    "- A plane perpendicular to the x axis: $x - x_0 = 0$\n",
    "- A cylinder parallel to the z axis: $(x - x_0)^2 + (y - y_0)^2 - R^2 = 0$\n",
    "- A sphere: $(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 - R^2 = 0$\n",
    "\n",
    "The full [list of available surfaces](https://docs.openmc.org/en/stable/pythonapi/base.html#building-geometry) is as follows.\n",
    "\n",
    "Planes:\n",
    "\n",
    "- `openmc.Plane` — An arbitrary plane of the form $Ax + By + Cz = D$\n",
    "- `openmc.XPlane` — A plane perpendicular to the x axis of the form $x - x_0 = 0$\n",
    "- `openmc.YPlane` — A plane perpendicular to the y axis of the form $y - y_0 = 0$\n",
    "- `openmc.ZPlane` — A plane perpendicular to the z axis of the form $z - z_0 = 0$\n",
    "\n",
    "Quadrics:\n",
    "\n",
    "- `openmc.XCylinder` — An infinite cylinder whose length is parallel to the x-axis of the form $(y - y_0)^2 + (z - z_0)^2 = r^2$\n",
    "- `openmc.YCylinder` — An infinite cylinder whose length is parallel to the x-axis of the form $(x - x_0)^2 + (z - z_0)^2 = r^2$\n",
    "- `openmc.ZCylinder` — An infinite cylinder whose length is parallel to the x-axis of the form $(x - x_0)^2 + (y - y_0)^2 = r^2$\n",
    "- `openmc.Sphere` — A sphere of the form $(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2$\n",
    "- `openmc.XCone` — A cone parallel to the x-axis of the form $(y - y_0)^2 + (z - z_0)^2 = r^2 (x - x_0)^2$\n",
    "- `openmc.YCone` — A cone parallel to the y-axis of the form $(x - x_0)^2 + (z - z_0)^2 = r^2 (y - y_0)^2$\n",
    "- `openmc.ZCone` — A cone parallel to the z-axis of the form $(x - x_0)^2 + (y - y_0)^2 = r^2 (z - z_0)^2$\n",
    "- `openmc.Quadric` — A generic quadric surface\n",
    "\n",
    "Torii:\n",
    "\n",
    "- `openmc.XTorus` — A torus of the form $(x - x_0)^2/B^2 + (\\sqrt{(y - y_0)^2 + (z - z_0)^2} - A)^2/C^2 - 1 = 0$\n",
    "- `openmc.YTorus` — A torus of the form $(y - y_0)^2/B^2 + (\\sqrt{(x - x_0)^2 + (z - z_0)^2} - A)^2/C^2 - 1 = 0$\n",
    "- `openmc.ZTorus` — A torus of the form $(z - z_0)^2/B^2 + (\\sqrt{(x - x_0)^2 + (y - y_0)^2} - A)^2/C^2 - 1 = 0$\n",
    "\n",
    "#### Half-Spaces\n",
    "\n",
    "<img src=\"csg_half.png\" alt=\"drawing\" width=\"600\"/>\n",
    "\n",
    "A surface *half-space* is the region whose points satisfy a positive or negative inequality of the surface equation. For example, for a sphere of radius one centered at the origin, the surface equation is $f(x,y,z) = x^2 + y^2 + z^2 - 1 = 0$. Thus, we say that the negative half-space of the sphere is defined as the collection of points satisfying $f(x,y,z) < 0$, which one can reason is the inside of the sphere. Conversely, the positive half-space of the sphere would correspond to all points outside of the sphere.\n",
    "\n",
    "#### Regions\n",
    "A region is then a combination of (or just one) half-spaces.\n",
    "\n",
    "#### Fills\n",
    "Finally, a cell is complete once we have defined what is _filling_ the cell, which may be one of:\n",
    "\n",
    "- material\n",
    "- nothing (`None`), or vacuum/void\n",
    "- universe\n",
    "- lattice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "sphere = openmc.Sphere(r=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that by default the sphere is centered at the origin so we didn't have to supply `x0`, `y0`, or `z0` arguments. Strictly speaking, we could have omitted `r` as well since it defaults to one. To get the negative or positive half-space, we simply need to apply the `-` or `+` unary operators, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "openmc.surface.Halfspace"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inside_sphere = -sphere\n",
    "outside_sphere = +sphere\n",
    "type(inside_sphere)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see if `inside_sphere` actually contains points inside the sphere:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True False\n",
      "False True\n"
     ]
    }
   ],
   "source": [
    "print((0,0,0) in inside_sphere, (0,0,2) in inside_sphere)\n",
    "print((0,0,0) in outside_sphere, (0,0,2) in outside_sphere)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Everything works as expected! Now that we understand how to create half-spaces, we can create more complex volumes by combining half-spaces using Boolean operators: `&` (intersection), `|` (union), and `~` (complement):\n",
    "\n",
    "- `&`: logical AND\n",
    "- `|`: logical OR\n",
    "- `~`: logical NOT\n",
    "\n",
    "For example, let's say we want to define a region that is the top part of the sphere (all points inside the sphere that have $z > 0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_plane = openmc.ZPlane(0)\n",
    "northern_hemisphere = -sphere & +z_plane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell = openmc.Cell()\n",
    "cell.region = northern_hemisphere\n",
    "\n",
    "# or...\n",
    "cell = openmc.Cell(region=northern_hemisphere)\n",
    "\n",
    "cell.fill = water"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Universes and in-line plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A collection of cells is known as a universe and can be used as a repeatable unit when creating a model. Although we don't need it yet, the benefit of creating a universe is that we can visualize our geometry while we're creating it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "universe = openmc.Universe()\n",
    "universe.add_cell(cell)\n",
    "\n",
    "# this also works\n",
    "universe = openmc.Universe(cells=[cell])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Universe` object has a `plot` method that will display our the universe as current constructed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='x [cm]', ylabel='y [cm]'>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 10276 characters]",
      "text/plain": [
       "<Figure size 258.065x259.74 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "universe.plot(width=(2.0, 2.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, the plot will appear in the $x$-$y$ plane. We can change that with the `basis` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='x [cm]', ylabel='z [cm]'>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 9812 characters]",
      "text/plain": [
       "<Figure size 258.065x259.74 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "universe.plot(width=(2.0, 2.0), basis='xz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have particular fondness for, say, fuchsia, we can tell the `plot()` method to make our cell that color."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='x [cm]', ylabel='z [cm]'>"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 9744 characters]",
      "text/plain": [
       "<Figure size 258.065x259.74 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "universe.plot(width=(2.0, 2.0), basis='xz',\n",
    "              colors={cell: 'fuchsia'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Boundary Conditions\n",
    "\n",
    "To specify boundary conditions, you simply need to set the `Surface.boundary_type` to one of:\n",
    "\n",
    "- `transmission` (default)\n",
    "- `vacuum`\n",
    "- `reflective`\n",
    "- `periodic` (either rotational or translational)\n",
    "- `white` (isotropic angular flux)\n",
    "\n",
    "<img src=\"mc_bcs.png\" alt=\"drawing\" width=\"600\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "northern_hemisphere.boundary_type = 'vacuum'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pin cell geometry\n",
    "\n",
    "We now have enough knowledge to create our pin-cell. We need three surfaces to define the fuel and clad:\n",
    "\n",
    "1. The outer surface of the fuel -- a cylinder parallel to the z axis\n",
    "2. The inner surface of the clad -- same as above\n",
    "3. The outer surface of the clad -- same as above\n",
    "\n",
    "These three surfaces will all be instances of `openmc.ZCylinder`, each with a different radius according to the specification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "fuel_outer_radius = openmc.ZCylinder(r=0.46955)\n",
    "clad_inner_radius = openmc.ZCylinder(r=0.47910)\n",
    "clad_outer_radius = openmc.ZCylinder(r=0.54640)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the surfaces created, we can now take advantage of the built-in operators on surfaces to create regions for the fuel, the gap, and the clad:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "fuel_region = -fuel_outer_radius\n",
    "gap_region = +fuel_outer_radius & -clad_inner_radius\n",
    "clad_region = +clad_inner_radius & -clad_outer_radius"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also create two z-planes in order to bound the geometry in the axial direction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "top = openmc.ZPlane(z0=150.0, boundary_type='vacuum')\n",
    "bot = openmc.ZPlane(z0=-150.0, boundary_type='vacuum')\n",
    "layer = +bot & -top"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can create corresponding cells that assign materials to these regions. As with materials, cells have unique IDs that are assigned either manually or automatically. Note that the gap cell doesn't have any material assigned (it is void, a.k.a. vacuum, by default)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "fuel = openmc.Cell()\n",
    "fuel.fill = uo2\n",
    "fuel.region = fuel_region & layer\n",
    "\n",
    "gap = openmc.Cell()\n",
    "gap.region = gap_region & layer\n",
    "\n",
    "clad = openmc.Cell()\n",
    "clad.fill = zirconium\n",
    "clad.region = clad_region & layer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we need to handle the coolant outside of our fuel pin. To do this, we create x- and y-planes that bound the geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "pitch = 1.44270\n",
    "left = openmc.XPlane(-pitch/2, boundary_type='reflective')\n",
    "right = openmc.XPlane(pitch/2, boundary_type='reflective')\n",
    "front = openmc.YPlane(-pitch/2, boundary_type='reflective')\n",
    "back = openmc.YPlane(pitch/2, boundary_type='reflective')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The water region is going to be everything outside of the clad outer radius and within the box formed as the intersection of four half-spaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "water_region = +left & -right & +front & -back & +clad_outer_radius & layer\n",
    "\n",
    "moderator = openmc.Cell()\n",
    "moderator.fill = water\n",
    "moderator.region = water_region"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final step is to assign the cells we created to a universe and tell OpenMC that this universe is the \"root\" universe in our geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "root_universe = openmc.Universe(cells=(fuel, gap, clad, moderator))\n",
    "model.geometry = openmc.Geometry(root_universe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: xlabel='x [cm]', ylabel='y [cm]'>"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 16196 characters]",
      "text/plain": [
       "<Figure size 407.742x410.39 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "root_universe.plot((0, 0, 0), width=(pitch, pitch), pixels=100000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Starting source and settings\n",
    "\n",
    "The Python API has a module `openmc.stats` with various univariate and multivariate probability distributions. We can use these distributions to create a starting source using the `openmc.Source` object. One can independently specify the spatial distribution (`space`), the angular distribution (`angle`), the energy distribution (`energy`), and the time distribution (`time`). For this example, we'll only specify the spatial distribution as a single point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a point source\n",
    "point = openmc.stats.Point((0, 0, 0))\n",
    "source = openmc.IndependentSource(space=point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.settings.source = source\n",
    "model.settings.batches = 100\n",
    "model.settings.inactive = 10\n",
    "model.settings.particles = 1000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running OpenMC\n",
    "\n",
    "Running OpenMC from Python can be done using the `model.run()` function. This function allows you to set the number of MPI processes and OpenMP threads, if need be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                %%%%%%%%%%%%%%%\n",
      "                           %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                    %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                     %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ##################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ###################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ####################     %%%%%%%%%%%%%%%%%%%%%%\n",
      "                #####################     %%%%%%%%%%%%%%%%%%%%%\n",
      "                ######################     %%%%%%%%%%%%%%%%%%%%\n",
      "                #######################     %%%%%%%%%%%%%%%%%%\n",
      "                 #######################     %%%%%%%%%%%%%%%%%\n",
      "                 ######################     %%%%%%%%%%%%%%%%%\n",
      "                  ####################     %%%%%%%%%%%%%%%%%\n",
      "                    #################     %%%%%%%%%%%%%%%%%\n",
      "                     ###############     %%%%%%%%%%%%%%%%\n",
      "                       ############     %%%%%%%%%%%%%%%\n",
      "                          ########     %%%%%%%%%%%%%%\n",
      "                                      %%%%%%%%%%%\n",
      "\n",
      "                 | The OpenMC Monte Carlo Code\n",
      "       Copyright | 2011-2024 MIT, UChicago Argonne LLC, and contributors\n",
      "         License | https://docs.openmc.org/en/latest/license.html\n",
      "         Version | 0.14.1-dev\n",
      "        Git SHA1 | 14ce3cec4b388ae8b7ec5b28db57dbcbff06f147\n",
      "       Date/Time | 2024-03-29 14:34:17\n",
      "   MPI Processes | 1\n",
      "  OpenMP Threads | 8\n",
      "\n",
      " Reading model XML file 'model.xml' ...\n",
      " WARNING: Other XML file input(s) are present. These files may be ignored in\n",
      "          favor of the model.xml file.\n",
      " Reading cross sections XML file...\n",
      " Reading Zr90 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr90.h5\n",
      " Reading Zr91 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr91.h5\n",
      " Reading Zr92 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr92.h5\n",
      " Reading Zr94 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr94.h5\n",
      " Reading Zr96 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr96.h5\n",
      " Reading U234 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U234.h5\n",
      " Reading U235 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U235.h5\n",
      " Reading U238 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U238.h5\n",
      " Reading U236 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U236.h5\n",
      " Reading O16 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O16.h5\n",
      " Reading H1 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H1.h5\n",
      " Reading H2 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H2.h5\n",
      " Reading O17 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O17.h5\n",
      " Reading Pu239 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Pu239.h5\n",
      " Reading Si28 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Si28.h5\n",
      " Reading c_H_in_H2O from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/c_H_in_H2O.h5\n",
      " Minimum neutron data temperature: 250 K\n",
      " Maximum neutron data temperature: 2500 K\n",
      " Preparing distributed cell instances...\n",
      " Reading plot XML file...\n",
      " Writing summary.h5 file...\n",
      " Maximum neutron transport energy: 20000000 eV for Zr90\n",
      " Initializing source particles...\n",
      "\n",
      " ====================>     K EIGENVALUE SIMULATION     <====================\n",
      "\n",
      "  Bat./Gen.      k            Average k\n",
      "  =========   ========   ====================\n",
      "        1/1    1.50315\n",
      "        2/1    1.44086\n",
      "        3/1    1.46540\n",
      "        4/1    1.37107\n",
      "        5/1    1.33015\n",
      "        6/1    1.48464\n",
      "        7/1    1.47998\n",
      "        8/1    1.28684\n",
      "        9/1    1.45556\n",
      "       10/1    1.40944\n",
      "       11/1    1.39122\n",
      "       12/1    1.32307    1.35714 +/- 0.03407\n",
      "       13/1    1.52588    1.41339 +/- 0.05959\n",
      "       14/1    1.46026    1.42511 +/- 0.04373\n",
      "       15/1    1.41125    1.42233 +/- 0.03399\n",
      "       16/1    1.45798    1.42827 +/- 0.02838\n",
      "       17/1    1.46100    1.43295 +/- 0.02444\n",
      "       18/1    1.42515    1.43197 +/- 0.02119\n",
      "       19/1    1.54914    1.44499 +/- 0.02277\n",
      "       20/1    1.43623    1.44412 +/- 0.02039\n",
      "       21/1    1.40115    1.44021 +/- 0.01885\n",
      "       22/1    1.43731    1.43997 +/- 0.01721\n",
      "       23/1    1.41573    1.43811 +/- 0.01594\n",
      "       24/1    1.46427    1.43997 +/- 0.01487\n",
      "       25/1    1.45526    1.44099 +/- 0.01389\n",
      "       26/1    1.33933    1.43464 +/- 0.01446\n",
      "       27/1    1.42147    1.43386 +/- 0.01360\n",
      "       28/1    1.36027    1.42978 +/- 0.01346\n",
      "       29/1    1.46546    1.43165 +/- 0.01287\n",
      "       30/1    1.45070    1.43261 +/- 0.01225\n",
      "       31/1    1.48473    1.43509 +/- 0.01191\n",
      "       32/1    1.47264    1.43680 +/- 0.01149\n",
      "       33/1    1.45237    1.43747 +/- 0.01100\n",
      "       34/1    1.41749    1.43664 +/- 0.01056\n",
      "       35/1    1.47971    1.43836 +/- 0.01027\n",
      "       36/1    1.45378    1.43896 +/- 0.00989\n",
      "       37/1    1.35765    1.43594 +/- 0.00998\n",
      "       38/1    1.40397    1.43480 +/- 0.00969\n",
      "       39/1    1.41697    1.43419 +/- 0.00937\n",
      "       40/1    1.40400    1.43318 +/- 0.00910\n",
      "       41/1    1.43270    1.43317 +/- 0.00881\n",
      "       42/1    1.36535    1.43105 +/- 0.00879\n",
      "       43/1    1.43240    1.43109 +/- 0.00851\n",
      "       44/1    1.46349    1.43204 +/- 0.00832\n",
      "       45/1    1.32866    1.42909 +/- 0.00860\n",
      "       46/1    1.32249    1.42613 +/- 0.00886\n",
      "       47/1    1.41734    1.42589 +/- 0.00862\n",
      "       48/1    1.39802    1.42516 +/- 0.00843\n",
      "       49/1    1.54696    1.42828 +/- 0.00878\n",
      "       50/1    1.44601    1.42872 +/- 0.00857\n",
      "       51/1    1.37546    1.42742 +/- 0.00846\n",
      "       52/1    1.48833    1.42887 +/- 0.00838\n",
      "       53/1    1.36864    1.42747 +/- 0.00830\n",
      "       54/1    1.45539    1.42811 +/- 0.00814\n",
      "       55/1    1.47463    1.42914 +/- 0.00802\n",
      "       56/1    1.43531    1.42927 +/- 0.00785\n",
      "       57/1    1.37650    1.42815 +/- 0.00776\n",
      "       58/1    1.43163    1.42822 +/- 0.00760\n",
      "       59/1    1.39161    1.42748 +/- 0.00748\n",
      "       60/1    1.48475    1.42862 +/- 0.00742\n",
      "       61/1    1.47918    1.42961 +/- 0.00734\n",
      "       62/1    1.47997    1.43058 +/- 0.00726\n",
      "       63/1    1.41811    1.43035 +/- 0.00712\n",
      "       64/1    1.41011    1.42997 +/- 0.00700\n",
      "       65/1    1.44239    1.43020 +/- 0.00688\n",
      "       66/1    1.42912    1.43018 +/- 0.00675\n",
      "       67/1    1.39312    1.42953 +/- 0.00666\n",
      "       68/1    1.47339    1.43028 +/- 0.00659\n",
      "       69/1    1.38559    1.42953 +/- 0.00652\n",
      "       70/1    1.42211    1.42940 +/- 0.00641\n",
      "       71/1    1.35744    1.42822 +/- 0.00642\n",
      "       72/1    1.40876    1.42791 +/- 0.00632\n",
      "       73/1    1.53582    1.42962 +/- 0.00645\n",
      "       74/1    1.44605    1.42988 +/- 0.00636\n",
      "       75/1    1.47428    1.43056 +/- 0.00629\n",
      "       76/1    1.37855    1.42977 +/- 0.00625\n",
      "       77/1    1.39439    1.42925 +/- 0.00618\n",
      "       78/1    1.49117    1.43016 +/- 0.00615\n",
      "       79/1    1.34815    1.42897 +/- 0.00618\n",
      "       80/1    1.30610    1.42721 +/- 0.00634\n",
      "       81/1    1.44634    1.42748 +/- 0.00625\n",
      "       82/1    1.41633    1.42733 +/- 0.00617\n",
      "       83/1    1.45974    1.42777 +/- 0.00610\n",
      "       84/1    1.45538    1.42814 +/- 0.00603\n",
      "       85/1    1.45563    1.42851 +/- 0.00596\n",
      "       86/1    1.37790    1.42785 +/- 0.00592\n",
      "       87/1    1.43465    1.42793 +/- 0.00584\n",
      "       88/1    1.32547    1.42662 +/- 0.00591\n",
      "       89/1    1.37892    1.42602 +/- 0.00587\n",
      "       90/1    1.47858    1.42667 +/- 0.00583\n",
      "       91/1    1.34703    1.42569 +/- 0.00584\n",
      "       92/1    1.47197    1.42625 +/- 0.00580\n",
      "       93/1    1.48311    1.42694 +/- 0.00577\n",
      "       94/1    1.39298    1.42653 +/- 0.00571\n",
      "       95/1    1.52197    1.42766 +/- 0.00576\n",
      "       96/1    1.35930    1.42686 +/- 0.00574\n",
      "       97/1    1.51529    1.42788 +/- 0.00577\n",
      "       98/1    1.45297    1.42816 +/- 0.00571\n",
      "       99/1    1.40311    1.42788 +/- 0.00565\n",
      "      100/1    1.42293    1.42783 +/- 0.00559\n",
      " Creating state point statepoint.100.h5...\n",
      "\n",
      " =======================>     TIMING STATISTICS     <=======================\n",
      "\n",
      " Total time for initialization     = 2.3567e+00 seconds\n",
      "   Reading cross sections          = 2.2972e+00 seconds\n",
      " Total time in simulation          = 1.7247e+00 seconds\n",
      "   Time in transport only          = 1.7078e+00 seconds\n",
      "   Time in inactive batches        = 1.6326e-01 seconds\n",
      "   Time in active batches          = 1.5614e+00 seconds\n",
      "   Time synchronizing fission bank = 6.8949e-03 seconds\n",
      "     Sampling source sites         = 5.3603e-03 seconds\n",
      "     SEND/RECV source sites        = 6.2301e-04 seconds\n",
      "   Time accumulating tallies       = 4.6476e-05 seconds\n",
      "   Time writing statepoints        = 5.4268e-03 seconds\n",
      " Total time for finalization       = 3.2870e-06 seconds\n",
      " Total time elapsed                = 4.0915e+00 seconds\n",
      " Calculation Rate (inactive)       = 61253.5 particles/second\n",
      " Calculation Rate (active)         = 57640.6 particles/second\n",
      "\n",
      " ============================>     RESULTS     <============================\n",
      "\n",
      " k-effective (Collision)     = 1.42881 +/- 0.00466\n",
      " k-effective (Track-length)  = 1.42783 +/- 0.00559\n",
      " k-effective (Absorption)    = 1.42667 +/- 0.00316\n",
      " Combined k-effective        = 1.42696 +/- 0.00299\n",
      " Leakage Fraction            = 0.00163 +/- 0.00017\n",
      "\n"
     ]
    }
   ],
   "source": [
    "statepoint = model.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geometry plotting\n",
    "\n",
    "We saw before that we could call the `Universe.plot()` method to show a universe while we were creating our geometry. There is also a built-in plotter in the codebase that is much faster than the Python plotter and has more options. The interface looks somewhat similar to the `Universe.plot()` method. Instead though, we create `Plot` instances, assign them to a `Plots` collection, export it to XML, and then run OpenMC in geometry plotting mode. As an example, let's specify that we want the plot to be colored by material (rather than by cell) and we assign yellow to fuel and blue to water."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot = openmc.Plot()\n",
    "plot.filename = 'pinplot'\n",
    "plot.width = (pitch, pitch)\n",
    "plot.pixels = (200, 200)\n",
    "plot.color_by = 'material'\n",
    "plot.colors = {uo2: 'yellow', water: 'blue'}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our plot created, we need to add it to a `Plots` collection which can be exported to XML."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.plots = openmc.Plots([plot])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run OpenMC in plotting mode by calling the `plot_geometry()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                %%%%%%%%%%%%%%%\n",
      "                           %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                    %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                     %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ##################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ###################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ####################     %%%%%%%%%%%%%%%%%%%%%%\n",
      "                #####################     %%%%%%%%%%%%%%%%%%%%%\n",
      "                ######################     %%%%%%%%%%%%%%%%%%%%\n",
      "                #######################     %%%%%%%%%%%%%%%%%%\n",
      "                 #######################     %%%%%%%%%%%%%%%%%\n",
      "                 ######################     %%%%%%%%%%%%%%%%%\n",
      "                  ####################     %%%%%%%%%%%%%%%%%\n",
      "                    #################     %%%%%%%%%%%%%%%%%\n",
      "                     ###############     %%%%%%%%%%%%%%%%\n",
      "                       ############     %%%%%%%%%%%%%%%\n",
      "                          ########     %%%%%%%%%%%%%%\n",
      "                                      %%%%%%%%%%%\n",
      "\n",
      "                 | The OpenMC Monte Carlo Code\n",
      "       Copyright | 2011-2024 MIT, UChicago Argonne LLC, and contributors\n",
      "         License | https://docs.openmc.org/en/latest/license.html\n",
      "         Version | 0.14.1-dev\n",
      "        Git SHA1 | 14ce3cec4b388ae8b7ec5b28db57dbcbff06f147\n",
      "       Date/Time | 2024-03-29 14:34:21\n",
      "   MPI Processes | 1\n",
      "  OpenMP Threads | 8\n",
      "\n",
      " Reading model XML file 'model.xml' ...\n",
      " WARNING: Other XML file input(s) are present. These files may be ignored in\n",
      "          favor of the model.xml file.\n",
      " Preparing distributed cell instances...\n",
      " Reading plot XML file...\n",
      "\n",
      " =======================>     PLOTTING SUMMARY     <========================\n",
      "\n",
      "Plot ID: 5\n",
      "Plot file: pinplot.png\n",
      "Universe depth: -1\n",
      "Plot Type: Slice\n",
      "Origin: 0 0 0\n",
      "Width: 1.4427 1.4427\n",
      "Coloring: Materials\n",
      "Basis: XY\n",
      "Pixels: 200 200\n",
      "\n",
      " Processing plot 5: pinplot.png...\n"
     ]
    }
   ],
   "source": [
    "model.plot_geometry()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can use functionality from IPython to display the `.png` image inline in our notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "AdvancedGeometry-solution.ipynb mc_bcs.png\r\n",
      "CMakeLists.txt                  mesh.i\r\n",
      "CODEOWNERS                      mesh_in.e\r\n",
      "CODE_OF_CONDUCT.md              model.xml\r\n",
      "CONTRIBUTING.md                 \u001b[1m\u001b[36mopenmc\u001b[m\u001b[m\r\n",
      "Dockerfile                      \u001b[1m\u001b[36mopenmc.egg-info\u001b[m\u001b[m\r\n",
      "LICENSE                         pincell.png\r\n",
      "Lattice-solution.ipynb          pinplot.png\r\n",
      "MANIFEST.in                     plots.xml\r\n",
      "Pincell-solution.ipynb          pyproject.toml\r\n",
      "Pincell.ipynb                   pytest.ini\r\n",
      "README.md                       \u001b[1m\u001b[36mscripts\u001b[m\u001b[m\r\n",
      "bcs.png                         settings.xml\r\n",
      "\u001b[1m\u001b[36mbuild\u001b[m\u001b[m                           \u001b[31msetup.py\u001b[m\u001b[m\r\n",
      "\u001b[1m\u001b[36mcmake\u001b[m\u001b[m                           \u001b[1m\u001b[36msrc\u001b[m\u001b[m\r\n",
      "csg_half.png                    statepoint.100.h5\r\n",
      "\u001b[1m\u001b[36mdocs\u001b[m\u001b[m                            summary.h5\r\n",
      "\u001b[1m\u001b[36mexamples\u001b[m\u001b[m                        tallies.out\r\n",
      "geometry.xml                    \u001b[1m\u001b[36mtests\u001b[m\u001b[m\r\n",
      "\u001b[1m\u001b[36minclude\u001b[m\u001b[m                         \u001b[1m\u001b[36mtools\u001b[m\u001b[m\r\n",
      "\u001b[1m\u001b[36mman\u001b[m\u001b[m                             \u001b[1m\u001b[36mvendor\u001b[m\u001b[m\r\n",
      "materials.xml\r\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAADICAIAAAAiOjnJAAAGjklEQVR4nO3dzXIbNxBF4WYq7+WsnbWfMutoHT+ZsoCLGpkiB5jB7R/M+dapROk5bMAkVb6ZvRsw2x/RPwDWRFiQICxIEBYkCAsShAUJwoIEYUGCsCBBWJAgLEgQFiQICxKEBQnCggRhQYKwIEFYkCAsSBAWJAgLEoQFCcKCxJ/RP0Be//73Y/ef+fuvfxx+kopu/MJqT0CjCO6iYSlieuaakV0oLM+YnrlOZOuHlaGnR8sXtmxYOXt6tGphC4ZVJamt9fJaKqyKSW2tlNciYVVPamuNvMqHtVJSW9XzKhzW9KS+fzv7LN9+Tv6R6uZVMqxZST0p6XboX/bFGGd1VjGvemGdr+qhp2MlvfZpqucLK9dWpbBOJvW5J0VMz3xM+GRhhfIqE9bhquJ6ejShsCptFQhrRlKxPT36NfOF88oe1oGqMq2o104tsORtpQ7rRFWZe3r0bsu1lTes0apqJrV1JK+0bSUNa6iqxHepUUfuXjnbShfW9RbVoxVWV66wDi2qlZLaGs4rVVuJfv2Lqj672eDHl6k+j8+ysfqHstCNqsfwrSvJ3koR1nhVV0hqa+xYzNBWoqNw11WrsgPHYrj4sDrX1YWragbaynDZCg6LqkZUaivyjjVS1cWT+s17530r8LIVtrGo6oRb/r0VExZVnZa9rfjL+zNUtae3rRABYfW8hjKPLJWeQYUsLe+wRqpiXe3q/XOif1vpjkKqGpT0vVPXsHZfN1R1SFdbzkvLL6zu/zGqOqBraJ5tJToKE+7zWlIN0CksDkG9XAdioo1FVaclGqBHWH3rKtFQKtt/19RnacVvrFQ3gzVkGKk8rL7XB+tqov1hOiyt4I2V4bW1pPDBasNiXQWJX1qRG4s7u1Lwdx+EYb1+TYTv6it4PWTp0oq9Y7GupCLHqwor/Mv86KF7TDEbi3PQTdSoA49CzkEHYUOWhMW1PZWQK3zUxmJduYkZdfxnhVjS/LA4BxPyPw1DNhbnoLOAgXMUQsI1LM7BQM7DnxxWx2nNORhiZ+zTr1kchZDwC4tzMJznI2BjQWJmWFywcnO9ZrGxIEFYkHAKi5t7Em4Pgo0FCc+wuLmH83sEbCxIEBYkpoXFr+UsYOJDZGNBgrAgQViQICxIEBYkCAsShAUJwoIEYUGCsCBBWJCYFlbgX5iOWSY+RDYWJAgLEp5hvTv+t/Alv0fAxoKEU1hvP/kaYApuD4KNBQnCgsTMsDreBeH+Hmhn+HPfiWRjQcIvLO7v4TwfARsLEpPD4pqVlesFy5w3FqdhIOfhcxRCIiQsTkNnAQOfH9br05rTMMTrsSu+S8dRCImosDgN3cSMWhIWp2Eq/ueghR6FLC0HYUOOCYul5SZq1Kqw+KWdEnSPKfZPhZyGUpHjFYbFFT5cyLW9idxYbz9/sLRk3mNfutqw+l4TtDXd/kjVl+Dgd945EEXCBysPi6XlLn5dWfjGsgSvrfVkGKlHWLuvD27x8+zf2X3eYozfWBu0dVKiATqF1be0LNVoqnm3jkPQ7RORRBsrw82gtFQD9Aur+7XC0jqga2ieH+C6biwORI1ch2CT6ChsaGtQV1X+vMPqed3QVrfeqvy/xRSwsUbawo6cVVnCo/COd033BH9/4bWYsDpfQ7T1XG9VUV/lDdtYtHVC9qrM7Bb72Dr/2vTv39qAbtIfpoKBPwPG/tpB8B1rZG/Z5VdXmaosPCyjrV6VqrIMYfW7cFtJ3wV9IfiOddd52bKP+5Zd48r16+n0V5VhXVmesGykLbvKdX54USWpylIdhUNDucCxWLgqS7WxmqG9ZWuuriM3qlRVWcKwmkPHotXPa/hGZfmSapKGZVdcXSssqru8Ydl4W1Y1r4PvJqStypKHZafaajIX9jH5xaqy/GHZobaaxHevI3epreRVWYmwmhl5WXRhp1ZUkz+ppkxYdqKtJq6wCT01VaqyWmE1J/Oy3wszTWSfpnr+Y75CSTX1wrIZbTUPhTXHOvtijLM+Ni5XlRUNq5mV192TzgZM/wJCxaSawmE10/NKom5STfmwmpXyqp5Us0hYTfW81kiqWSqspmJeKyXVLBhWUyWv9ZJqlg3rLmdhq/Z0t35YdxkKW76nuwuFteUZ2XVi2rpoWFuKyK4Z0xZhPdUTHAE9Q1iQSPTrX1gJYUGCsCBBWJAgLEgQFiQICxKEBQnCggRhQYKwIEFYkCAsSBAWJAgLEoQFCcKCBGFBgrAgQViQICxIEBYkCAsS/wMmEHBNh1IH6gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "!ls\n",
    "\n",
    "from IPython.display import Image\n",
    "Image(\"pinplot.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tallies\n",
    "\n",
    "In this section, we'll be looking at how to extract custom information from an OpenMC simulation in what is known as a \"tally.\" A tally accumulates statistical information during the simulation about particles when they eneter regions of phase space specified on the tally. The limits of these regions are set by \"filters\" applied to the tally. Scores and nuclides can also be applied to tallies to indicate what type of information is kept about the particle (e.g. reaction types, flux, heat, etc.).\n",
    "\n",
    "Any tally in OpenMC can be described with the following form:\n",
    "\n",
    "$$ \n",
    " X = \\underbrace{\\int d\\mathbf{r} \\int d\\mathbf{\\Omega} \\int\n",
    "    dE}_{\\text{filters}} \\underbrace{f(\\mathbf{r}, \\mathbf{\\Omega},\n",
    "    E)}_{\\text{scores}} \\underbrace{\\psi (\\mathbf{r}, \\mathbf{\\Omega}, E)}_{\\text{angular flux}}\n",
    "$$\n",
    "\n",
    "where filters set the limits of the integrals and the scoring function is convolved with particle information (e.g. reaction type, current material, etc.). For example, if you wanted to calculate the fission reaction rate caused by fast neutrons in cell 3, your tally becomes\n",
    "\n",
    "$$ \n",
    " X = \\int_\\text{cell 3} d\\mathbf{r} \\int_{4\\pi} d\\mathbf{\\Omega} \\int_{1 MeV}^{20 MeV}\n",
    "    dE \\ \\ \\Sigma_f(\\mathbf{r}, \\mathbf{\\Omega},\n",
    "    E) \\psi (\\mathbf{r}, \\mathbf{\\Omega}, E)\n",
    "$$\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "A full list of scores and their meanings can be found <a href=https://docs.openmc.org/en/stable/usersguide/tallies.html#scores >here</a>.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this exercise we'll be adding tallies to perform a few different tasks:\n",
    "\n",
    "\n",
    "  **1. Determine the energy and heat produced per fission** \\\n",
    "  **2. Plot the flux spectrum of the pincell** \\\n",
    "  **3. Plot reaction types based on material**\n",
    "  \n",
    "First, to determine the recoverable energy produced per fission we'll create a tally without filters to gather information on the fission reaction rate (\"`fission`\") and recoverable fission energy (\"`kappa-fission`\"). Because we want this information talllied throughout the model, a \"global\" tally, no filters need to be applied."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 1: Energy released per fission"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tally\n",
      "\tID             =\t1\n",
      "\tName           =\t\n",
      "\tFilters        =\t\n",
      "\tNuclides       =\t\n",
      "\tScores         =\t['fission', 'kappa-fission']\n",
      "\tEstimator      =\tNone\n",
      "\tMultiply dens. =\tTrue\n"
     ]
    }
   ],
   "source": [
    "tally1 = openmc.Tally()\n",
    "tally1.scores = ['fission', 'kappa-fission']\n",
    "print(tally1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.tallies = openmc.Tallies([tally1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can simply re-run our model with this additional tally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                %%%%%%%%%%%%%%%\n",
      "                           %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                    %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                     %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ##################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ###################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ####################     %%%%%%%%%%%%%%%%%%%%%%\n",
      "                #####################     %%%%%%%%%%%%%%%%%%%%%\n",
      "                ######################     %%%%%%%%%%%%%%%%%%%%\n",
      "                #######################     %%%%%%%%%%%%%%%%%%\n",
      "                 #######################     %%%%%%%%%%%%%%%%%\n",
      "                 ######################     %%%%%%%%%%%%%%%%%\n",
      "                  ####################     %%%%%%%%%%%%%%%%%\n",
      "                    #################     %%%%%%%%%%%%%%%%%\n",
      "                     ###############     %%%%%%%%%%%%%%%%\n",
      "                       ############     %%%%%%%%%%%%%%%\n",
      "                          ########     %%%%%%%%%%%%%%\n",
      "                                      %%%%%%%%%%%\n",
      "\n",
      "                 | The OpenMC Monte Carlo Code\n",
      "       Copyright | 2011-2024 MIT, UChicago Argonne LLC, and contributors\n",
      "         License | https://docs.openmc.org/en/latest/license.html\n",
      "         Version | 0.14.1-dev\n",
      "        Git SHA1 | 14ce3cec4b388ae8b7ec5b28db57dbcbff06f147\n",
      "       Date/Time | 2024-03-29 14:34:21\n",
      "   MPI Processes | 1\n",
      "  OpenMP Threads | 8\n",
      "\n",
      " Reading model XML file 'model.xml' ...\n",
      " WARNING: Other XML file input(s) are present. These files may be ignored in\n",
      "          favor of the model.xml file.\n",
      " Reading cross sections XML file...\n",
      " Reading Zr90 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr90.h5\n",
      " Reading Zr91 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr91.h5\n",
      " Reading Zr92 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr92.h5\n",
      " Reading Zr94 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr94.h5\n",
      " Reading Zr96 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr96.h5\n",
      " Reading U234 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U234.h5\n",
      " Reading U235 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U235.h5\n",
      " Reading U238 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U238.h5\n",
      " Reading U236 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U236.h5\n",
      " Reading O16 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O16.h5\n",
      " Reading H1 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H1.h5\n",
      " Reading H2 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H2.h5\n",
      " Reading O17 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O17.h5\n",
      " Reading Pu239 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Pu239.h5\n",
      " Reading Si28 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Si28.h5\n",
      " Reading c_H_in_H2O from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/c_H_in_H2O.h5\n",
      " Minimum neutron data temperature: 250 K\n",
      " Maximum neutron data temperature: 2500 K\n",
      " Preparing distributed cell instances...\n",
      " Writing summary.h5 file...\n",
      " Maximum neutron transport energy: 20000000 eV for Zr90\n",
      " Initializing source particles...\n",
      "\n",
      " ====================>     K EIGENVALUE SIMULATION     <====================\n",
      "\n",
      "  Bat./Gen.      k            Average k\n",
      "  =========   ========   ====================\n",
      "        1/1    1.50315\n",
      "        2/1    1.44086\n",
      "        3/1    1.46540\n",
      "        4/1    1.37107\n",
      "        5/1    1.33015\n",
      "        6/1    1.48464\n",
      "        7/1    1.47998\n",
      "        8/1    1.28684\n",
      "        9/1    1.45556\n",
      "       10/1    1.40944\n",
      "       11/1    1.39122\n",
      "       12/1    1.32307    1.35714 +/- 0.03407\n",
      "       13/1    1.52588    1.41339 +/- 0.05959\n",
      "       14/1    1.46026    1.42511 +/- 0.04373\n",
      "       15/1    1.41125    1.42233 +/- 0.03399\n",
      "       16/1    1.45798    1.42827 +/- 0.02838\n",
      "       17/1    1.46100    1.43295 +/- 0.02444\n",
      "       18/1    1.42515    1.43197 +/- 0.02119\n",
      "       19/1    1.54914    1.44499 +/- 0.02277\n",
      "       20/1    1.43623    1.44412 +/- 0.02039\n",
      "       21/1    1.40115    1.44021 +/- 0.01885\n",
      "       22/1    1.43731    1.43997 +/- 0.01721\n",
      "       23/1    1.41573    1.43811 +/- 0.01594\n",
      "       24/1    1.46427    1.43997 +/- 0.01487\n",
      "       25/1    1.45526    1.44099 +/- 0.01389\n",
      "       26/1    1.33933    1.43464 +/- 0.01446\n",
      "       27/1    1.42147    1.43386 +/- 0.01360\n",
      "       28/1    1.36027    1.42978 +/- 0.01346\n",
      "       29/1    1.46546    1.43165 +/- 0.01287\n",
      "       30/1    1.45070    1.43261 +/- 0.01225\n",
      "       31/1    1.48473    1.43509 +/- 0.01191\n",
      "       32/1    1.47264    1.43680 +/- 0.01149\n",
      "       33/1    1.45237    1.43747 +/- 0.01100\n",
      "       34/1    1.41749    1.43664 +/- 0.01056\n",
      "       35/1    1.47971    1.43836 +/- 0.01027\n",
      "       36/1    1.45378    1.43896 +/- 0.00989\n",
      "       37/1    1.35765    1.43594 +/- 0.00998\n",
      "       38/1    1.40397    1.43480 +/- 0.00969\n",
      "       39/1    1.41697    1.43419 +/- 0.00937\n",
      "       40/1    1.40400    1.43318 +/- 0.00910\n",
      "       41/1    1.43270    1.43317 +/- 0.00881\n",
      "       42/1    1.36535    1.43105 +/- 0.00879\n",
      "       43/1    1.43240    1.43109 +/- 0.00851\n",
      "       44/1    1.46349    1.43204 +/- 0.00832\n",
      "       45/1    1.32866    1.42909 +/- 0.00860\n",
      "       46/1    1.32249    1.42613 +/- 0.00886\n",
      "       47/1    1.41734    1.42589 +/- 0.00862\n",
      "       48/1    1.39802    1.42516 +/- 0.00843\n",
      "       49/1    1.54696    1.42828 +/- 0.00878\n",
      "       50/1    1.44601    1.42872 +/- 0.00857\n",
      "       51/1    1.37546    1.42742 +/- 0.00846\n",
      "       52/1    1.48833    1.42887 +/- 0.00838\n",
      "       53/1    1.36864    1.42747 +/- 0.00830\n",
      "       54/1    1.45539    1.42811 +/- 0.00814\n",
      "       55/1    1.47463    1.42914 +/- 0.00802\n",
      "       56/1    1.43531    1.42927 +/- 0.00785\n",
      "       57/1    1.37650    1.42815 +/- 0.00776\n",
      "       58/1    1.43163    1.42822 +/- 0.00760\n",
      "       59/1    1.39161    1.42748 +/- 0.00748\n",
      "       60/1    1.48475    1.42862 +/- 0.00742\n",
      "       61/1    1.47918    1.42961 +/- 0.00734\n",
      "       62/1    1.47997    1.43058 +/- 0.00726\n",
      "       63/1    1.41811    1.43035 +/- 0.00712\n",
      "       64/1    1.41011    1.42997 +/- 0.00700\n",
      "       65/1    1.44239    1.43020 +/- 0.00688\n",
      "       66/1    1.42912    1.43018 +/- 0.00675\n",
      "       67/1    1.39312    1.42953 +/- 0.00666\n",
      "       68/1    1.47339    1.43028 +/- 0.00659\n",
      "       69/1    1.38559    1.42953 +/- 0.00652\n",
      "       70/1    1.42211    1.42940 +/- 0.00641\n",
      "       71/1    1.35744    1.42822 +/- 0.00642\n",
      "       72/1    1.40876    1.42791 +/- 0.00632\n",
      "       73/1    1.53582    1.42962 +/- 0.00645\n",
      "       74/1    1.44605    1.42988 +/- 0.00636\n",
      "       75/1    1.47428    1.43056 +/- 0.00629\n",
      "       76/1    1.37855    1.42977 +/- 0.00625\n",
      "       77/1    1.39439    1.42925 +/- 0.00618\n",
      "       78/1    1.49117    1.43016 +/- 0.00615\n",
      "       79/1    1.34815    1.42897 +/- 0.00618\n",
      "       80/1    1.30610    1.42721 +/- 0.00634\n",
      "       81/1    1.44634    1.42748 +/- 0.00625\n",
      "       82/1    1.41633    1.42733 +/- 0.00617\n",
      "       83/1    1.45974    1.42777 +/- 0.00610\n",
      "       84/1    1.45538    1.42814 +/- 0.00603\n",
      "       85/1    1.45563    1.42851 +/- 0.00596\n",
      "       86/1    1.37790    1.42785 +/- 0.00592\n",
      "       87/1    1.43465    1.42793 +/- 0.00584\n",
      "       88/1    1.32547    1.42662 +/- 0.00591\n",
      "       89/1    1.37892    1.42602 +/- 0.00587\n",
      "       90/1    1.47858    1.42667 +/- 0.00583\n",
      "       91/1    1.34703    1.42569 +/- 0.00584\n",
      "       92/1    1.47197    1.42625 +/- 0.00580\n",
      "       93/1    1.48311    1.42694 +/- 0.00577\n",
      "       94/1    1.39298    1.42653 +/- 0.00571\n",
      "       95/1    1.52197    1.42766 +/- 0.00576\n",
      "       96/1    1.35930    1.42686 +/- 0.00574\n",
      "       97/1    1.51529    1.42788 +/- 0.00577\n",
      "       98/1    1.45297    1.42816 +/- 0.00571\n",
      "       99/1    1.40311    1.42788 +/- 0.00565\n",
      "      100/1    1.42293    1.42783 +/- 0.00559\n",
      " Creating state point statepoint.100.h5...\n",
      "\n",
      " =======================>     TIMING STATISTICS     <=======================\n",
      "\n",
      " Total time for initialization     = 2.0523e+00 seconds\n",
      "   Reading cross sections          = 1.9999e+00 seconds\n",
      " Total time in simulation          = 2.1030e+00 seconds\n",
      "   Time in transport only          = 2.0748e+00 seconds\n",
      "   Time in inactive batches        = 1.4362e-01 seconds\n",
      "   Time in active batches          = 1.9594e+00 seconds\n",
      "   Time synchronizing fission bank = 7.3517e-03 seconds\n",
      "     Sampling source sites         = 5.7251e-03 seconds\n",
      "     SEND/RECV source sites        = 6.3934e-04 seconds\n",
      "   Time accumulating tallies       = 1.0267e-02 seconds\n",
      "   Time writing statepoints        = 5.5837e-03 seconds\n",
      " Total time for finalization       = 2.6649e-04 seconds\n",
      " Total time elapsed                = 4.1654e+00 seconds\n",
      " Calculation Rate (inactive)       = 69627.5 particles/second\n",
      " Calculation Rate (active)         = 45932.7 particles/second\n",
      "\n",
      " ============================>     RESULTS     <============================\n",
      "\n",
      " k-effective (Collision)     = 1.42881 +/- 0.00466\n",
      " k-effective (Track-length)  = 1.42783 +/- 0.00559\n",
      " k-effective (Absorption)    = 1.42667 +/- 0.00316\n",
      " Combined k-effective        = 1.42696 +/- 0.00299\n",
      " Leakage Fraction            = 0.00163 +/- 0.00017\n",
      "\n"
     ]
    }
   ],
   "source": [
    "statepoint = model.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we list our current directory, we see that several new files have been created as a result of this run: `summary.h5`, `tallies.out`, and `statepoint.50.h5`. The `tallies.out` file contains a text output of all user-specified tallies for the simulation. The summary file contains information about the simulation's setup (geometry, materials, meshes, etc.) in an HDF5 format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " ============================>     TALLY 1     <============================\r\n",
      "\r\n",
      " Total Material\r\n",
      "   Fission Rate                         0.581783 +/- 0.00229294\r\n",
      "   Kappa-Fission Rate                   1.12634e+08 +/- 443455\r\n"
     ]
    }
   ],
   "source": [
    "!cat tallies.out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can be useful to quickly look at simple tally results, but isn't a great format to post-process simulation data. For that we'll look to the statepoint file. The statepoint file contains information about simulation results including tally specifications and data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/anovak/projects/openmc/statepoint.100.h5\n"
     ]
    }
   ],
   "source": [
    "print(statepoint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To extract information from the statepoint file we'll create an `openmc.StatePoint` object. The `statepoint.get_tally` function will search for tallies by scores, filters, nuclides, ids, and return the closest match. Exact matches can be speficied as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tally\n",
      "\tID             =\t1\n",
      "\tName           =\t\n",
      "\tFilters        =\t\n",
      "\tNuclides       =\ttotal\n",
      "\tScores         =\t['fission', 'kappa-fission']\n",
      "\tEstimator      =\ttracklength\n",
      "\tMultiply dens. =\tTrue\n",
      "Tally\n",
      "\tID             =\t1\n",
      "\tName           =\t\n",
      "\tFilters        =\t\n",
      "\tNuclides       =\ttotal\n",
      "\tScores         =\t['fission', 'kappa-fission']\n",
      "\tEstimator      =\ttracklength\n",
      "\tMultiply dens. =\tTrue\n"
     ]
    }
   ],
   "source": [
    "with openmc.StatePoint(statepoint) as sp:\n",
    "    my_tally = sp.get_tally(scores=['fission', 'kappa-fission'])\n",
    "    print(my_tally)\n",
    "    \n",
    "    my_tally = sp.get_tally(id=tally1.id)\n",
    "    print(my_tally)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>A quick aside on how statepoint objects interact with summary files:</b>\n",
    "\n",
    "\n",
    "The `openmc.statepoint` object will read information from the `summary.h5` file if one is present, keeping that file open in the Python interpreter. The open `summary.h5` file can interfere with the initialization of subsequent OpenMC simulations. It is recommended that information be extracted from statepoints within a [context manager](https://book.pythontips.com/en/latest/context_managers.html) as we do here. Alternatively, making sure to call the `openmc.StatePoint.close` method will work also. For more details please look to the [relevant section in the user's guide](https://docs.openmc.org/en/stable/usersguide/troubleshoot.html#runtimeerror-failed-to-open-hdf5-file-with-mode-w-summary-h5).   \n",
    "</div>\n",
    "\n",
    "To compute the energy released per fission event, we can simply take the tallied energy released per fission and divide it by the fission rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "#my_tally.get_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MeV per fission:  193.6014433428488\n"
     ]
    }
   ],
   "source": [
    "fission_rate = my_tally.get_values(scores=['fission']).flatten()[0]\n",
    "kappa_fission_rate = my_tally.get_values(scores=['kappa-fission']).flatten()[0]\n",
    "ev_per_fission = kappa_fission_rate / fission_rate\n",
    "mev_per_fission = ev_per_fission * 1e-6\n",
    "print('MeV per fission: ', mev_per_fission)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MeV per fission: 193.601443+/-1.078524\n"
     ]
    }
   ],
   "source": [
    "fission_rate_std_dev = my_tally.get_values(scores=['fission'], value='std_dev').flatten()[0]\n",
    "kappa_fission_rate_std_dev = my_tally.get_values(scores=['kappa-fission'],value='std_dev').flatten()[0]\n",
    "\n",
    "from uncertainties import ufloat\n",
    "f = ufloat(fission_rate, fission_rate_std_dev) \n",
    "kf = ufloat(kappa_fission_rate, kappa_fission_rate_std_dev) \n",
    "\n",
    "print('MeV per fission: {:.6f}'.format(kf / f * 1e-6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with most values coming out of an MC code, these values are per source-particle. In this case these units cancel out, but this will not be the case in our next example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 2: Plot the neutron flux spectrum\n",
    "\n",
    "To perform this task, we'll be applying a tally with an energy filter and a score. OpenMC's data module contains different group structures. For this problem we'll use the CASMO-70 group structure. An energy filter can easily be created from a pre-defined group structure in OpenMC as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['CASMO-2', 'CASMO-4', 'CASMO-8', 'CASMO-16', 'CASMO-25', 'CASMO-40', 'VITAMIN-J-42', 'SCALE-44', 'MPACT-51', 'MPACT-60', 'MPACT-69', 'CASMO-70', 'XMAS-172', 'VITAMIN-J-175', 'SCALE-252', 'TRIPOLI-315', 'SHEM-361', 'CCFE-709', 'UKAEA-1102', 'ECCO-1968'])\n"
     ]
    }
   ],
   "source": [
    "print(openmc.mgxs.GROUP_STRUCTURES.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "70"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "energy_filter = openmc.EnergyFilter.from_group_structure('CASMO-70')\n",
    "len(energy_filter.bins)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrum_tally = openmc.Tally()\n",
    "spectrum_tally.filters = [energy_filter]\n",
    "spectrum_tally.scores = ['flux']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                %%%%%%%%%%%%%%%\n",
      "                           %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                    %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                     %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ##################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ###################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ####################     %%%%%%%%%%%%%%%%%%%%%%\n",
      "                #####################     %%%%%%%%%%%%%%%%%%%%%\n",
      "                ######################     %%%%%%%%%%%%%%%%%%%%\n",
      "                #######################     %%%%%%%%%%%%%%%%%%\n",
      "                 #######################     %%%%%%%%%%%%%%%%%\n",
      "                 ######################     %%%%%%%%%%%%%%%%%\n",
      "                  ####################     %%%%%%%%%%%%%%%%%\n",
      "                    #################     %%%%%%%%%%%%%%%%%\n",
      "                     ###############     %%%%%%%%%%%%%%%%\n",
      "                       ############     %%%%%%%%%%%%%%%\n",
      "                          ########     %%%%%%%%%%%%%%\n",
      "                                      %%%%%%%%%%%\n",
      "\n",
      "                 | The OpenMC Monte Carlo Code\n",
      "       Copyright | 2011-2024 MIT, UChicago Argonne LLC, and contributors\n",
      "         License | https://docs.openmc.org/en/latest/license.html\n",
      "         Version | 0.14.1-dev\n",
      "        Git SHA1 | 14ce3cec4b388ae8b7ec5b28db57dbcbff06f147\n",
      "       Date/Time | 2024-03-29 14:34:26\n",
      "   MPI Processes | 1\n",
      "  OpenMP Threads | 8\n",
      "\n",
      " Reading model XML file 'model.xml' ...\n",
      " WARNING: Other XML file input(s) are present. These files may be ignored in\n",
      "          favor of the model.xml file.\n",
      " Reading cross sections XML file...\n",
      " Reading Zr90 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr90.h5\n",
      " Reading Zr91 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr91.h5\n",
      " Reading Zr92 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr92.h5\n",
      " Reading Zr94 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr94.h5\n",
      " Reading Zr96 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr96.h5\n",
      " Reading U234 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U234.h5\n",
      " Reading U235 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U235.h5\n",
      " Reading U238 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U238.h5\n",
      " Reading U236 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U236.h5\n",
      " Reading O16 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O16.h5\n",
      " Reading H1 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H1.h5\n",
      " Reading H2 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H2.h5\n",
      " Reading O17 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O17.h5\n",
      " Reading Pu239 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Pu239.h5\n",
      " Reading Si28 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Si28.h5\n",
      " Reading c_H_in_H2O from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/c_H_in_H2O.h5\n",
      " Minimum neutron data temperature: 250 K\n",
      " Maximum neutron data temperature: 2500 K\n",
      " Preparing distributed cell instances...\n",
      " Writing summary.h5 file...\n",
      " Maximum neutron transport energy: 20000000 eV for Zr90\n",
      " Initializing source particles...\n",
      "\n",
      " ====================>     K EIGENVALUE SIMULATION     <====================\n",
      "\n",
      "  Bat./Gen.      k            Average k\n",
      "  =========   ========   ====================\n",
      "        1/1    1.50315\n",
      "        2/1    1.44086\n",
      "        3/1    1.46540\n",
      "        4/1    1.37107\n",
      "        5/1    1.33015\n",
      "        6/1    1.48464\n",
      "        7/1    1.47998\n",
      "        8/1    1.28684\n",
      "        9/1    1.45556\n",
      "       10/1    1.40944\n",
      "       11/1    1.39122\n",
      "       12/1    1.32307    1.35714 +/- 0.03407\n",
      "       13/1    1.52588    1.41339 +/- 0.05959\n",
      "       14/1    1.46026    1.42511 +/- 0.04373\n",
      "       15/1    1.41125    1.42233 +/- 0.03399\n",
      "       16/1    1.45798    1.42827 +/- 0.02838\n",
      "       17/1    1.46100    1.43295 +/- 0.02444\n",
      "       18/1    1.42515    1.43197 +/- 0.02119\n",
      "       19/1    1.54914    1.44499 +/- 0.02277\n",
      "       20/1    1.43623    1.44412 +/- 0.02039\n",
      "       21/1    1.40115    1.44021 +/- 0.01885\n",
      "       22/1    1.43731    1.43997 +/- 0.01721\n",
      "       23/1    1.41573    1.43811 +/- 0.01594\n",
      "       24/1    1.46427    1.43997 +/- 0.01487\n",
      "       25/1    1.45526    1.44099 +/- 0.01389\n",
      "       26/1    1.33933    1.43464 +/- 0.01446\n",
      "       27/1    1.42147    1.43386 +/- 0.01360\n",
      "       28/1    1.36027    1.42978 +/- 0.01346\n",
      "       29/1    1.46546    1.43165 +/- 0.01287\n",
      "       30/1    1.45070    1.43261 +/- 0.01225\n",
      "       31/1    1.48473    1.43509 +/- 0.01191\n",
      "       32/1    1.47264    1.43680 +/- 0.01149\n",
      "       33/1    1.45237    1.43747 +/- 0.01100\n",
      "       34/1    1.41749    1.43664 +/- 0.01056\n",
      "       35/1    1.47971    1.43836 +/- 0.01027\n",
      "       36/1    1.45378    1.43896 +/- 0.00989\n",
      "       37/1    1.35765    1.43594 +/- 0.00998\n",
      "       38/1    1.40397    1.43480 +/- 0.00969\n",
      "       39/1    1.41697    1.43419 +/- 0.00937\n",
      "       40/1    1.40400    1.43318 +/- 0.00910\n",
      "       41/1    1.43270    1.43317 +/- 0.00881\n",
      "       42/1    1.36535    1.43105 +/- 0.00879\n",
      "       43/1    1.43240    1.43109 +/- 0.00851\n",
      "       44/1    1.46349    1.43204 +/- 0.00832\n",
      "       45/1    1.32866    1.42909 +/- 0.00860\n",
      "       46/1    1.32249    1.42613 +/- 0.00886\n",
      "       47/1    1.41734    1.42589 +/- 0.00862\n",
      "       48/1    1.39802    1.42516 +/- 0.00843\n",
      "       49/1    1.54696    1.42828 +/- 0.00878\n",
      "       50/1    1.44601    1.42872 +/- 0.00857\n",
      "       51/1    1.37546    1.42742 +/- 0.00846\n",
      "       52/1    1.48833    1.42887 +/- 0.00838\n",
      "       53/1    1.36864    1.42747 +/- 0.00830\n",
      "       54/1    1.45539    1.42811 +/- 0.00814\n",
      "       55/1    1.47463    1.42914 +/- 0.00802\n",
      "       56/1    1.43531    1.42927 +/- 0.00785\n",
      "       57/1    1.37650    1.42815 +/- 0.00776\n",
      "       58/1    1.43163    1.42822 +/- 0.00760\n",
      "       59/1    1.39161    1.42748 +/- 0.00748\n",
      "       60/1    1.48475    1.42862 +/- 0.00742\n",
      "       61/1    1.47918    1.42961 +/- 0.00734\n",
      "       62/1    1.47997    1.43058 +/- 0.00726\n",
      "       63/1    1.41811    1.43035 +/- 0.00712\n",
      "       64/1    1.41011    1.42997 +/- 0.00700\n",
      "       65/1    1.44239    1.43020 +/- 0.00688\n",
      "       66/1    1.42912    1.43018 +/- 0.00675\n",
      "       67/1    1.39312    1.42953 +/- 0.00666\n",
      "       68/1    1.47339    1.43028 +/- 0.00659\n",
      "       69/1    1.38559    1.42953 +/- 0.00652\n",
      "       70/1    1.42211    1.42940 +/- 0.00641\n",
      "       71/1    1.35744    1.42822 +/- 0.00642\n",
      "       72/1    1.40876    1.42791 +/- 0.00632\n",
      "       73/1    1.53582    1.42962 +/- 0.00645\n",
      "       74/1    1.44605    1.42988 +/- 0.00636\n",
      "       75/1    1.47428    1.43056 +/- 0.00629\n",
      "       76/1    1.37855    1.42977 +/- 0.00625\n",
      "       77/1    1.39439    1.42925 +/- 0.00618\n",
      "       78/1    1.49117    1.43016 +/- 0.00615\n",
      "       79/1    1.34815    1.42897 +/- 0.00618\n",
      "       80/1    1.30610    1.42721 +/- 0.00634\n",
      "       81/1    1.44634    1.42748 +/- 0.00625\n",
      "       82/1    1.41633    1.42733 +/- 0.00617\n",
      "       83/1    1.45974    1.42777 +/- 0.00610\n",
      "       84/1    1.45538    1.42814 +/- 0.00603\n",
      "       85/1    1.45563    1.42851 +/- 0.00596\n",
      "       86/1    1.37790    1.42785 +/- 0.00592\n",
      "       87/1    1.43465    1.42793 +/- 0.00584\n",
      "       88/1    1.32547    1.42662 +/- 0.00591\n",
      "       89/1    1.37892    1.42602 +/- 0.00587\n",
      "       90/1    1.47858    1.42667 +/- 0.00583\n",
      "       91/1    1.34703    1.42569 +/- 0.00584\n",
      "       92/1    1.47197    1.42625 +/- 0.00580\n",
      "       93/1    1.48311    1.42694 +/- 0.00577\n",
      "       94/1    1.39298    1.42653 +/- 0.00571\n",
      "       95/1    1.52197    1.42766 +/- 0.00576\n",
      "       96/1    1.35930    1.42686 +/- 0.00574\n",
      "       97/1    1.51529    1.42788 +/- 0.00577\n",
      "       98/1    1.45297    1.42816 +/- 0.00571\n",
      "       99/1    1.40311    1.42788 +/- 0.00565\n",
      "      100/1    1.42293    1.42783 +/- 0.00559\n",
      " Creating state point statepoint.100.h5...\n",
      "\n",
      " =======================>     TIMING STATISTICS     <=======================\n",
      "\n",
      " Total time for initialization     = 2.0435e+00 seconds\n",
      "   Reading cross sections          = 1.9552e+00 seconds\n",
      " Total time in simulation          = 1.9591e+00 seconds\n",
      "   Time in transport only          = 1.9310e+00 seconds\n",
      "   Time in inactive batches        = 1.5164e-01 seconds\n",
      "   Time in active batches          = 1.8075e+00 seconds\n",
      "   Time synchronizing fission bank = 6.6273e-03 seconds\n",
      "     Sampling source sites         = 5.1675e-03 seconds\n",
      "     SEND/RECV source sites        = 5.9234e-04 seconds\n",
      "   Time accumulating tallies       = 8.8949e-03 seconds\n",
      "   Time writing statepoints        = 8.1831e-03 seconds\n",
      " Total time for finalization       = 4.5152e-04 seconds\n",
      " Total time elapsed                = 4.0131e+00 seconds\n",
      " Calculation Rate (inactive)       = 65946 particles/second\n",
      " Calculation Rate (active)         = 49793.7 particles/second\n",
      "\n",
      " ============================>     RESULTS     <============================\n",
      "\n",
      " k-effective (Collision)     = 1.42881 +/- 0.00466\n",
      " k-effective (Track-length)  = 1.42783 +/- 0.00559\n",
      " k-effective (Absorption)    = 1.42667 +/- 0.00316\n",
      " Combined k-effective        = 1.42696 +/- 0.00299\n",
      " Leakage Fraction            = 0.00163 +/- 0.00017\n",
      "\n"
     ]
    }
   ],
   "source": [
    "model.tallies += [spectrum_tally]\n",
    "statepoint = model.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " ============================>     TALLY 1     <============================\r\n",
      "\r\n",
      " Total Material\r\n",
      "   Fission Rate                         0.581783 +/- 0.00229294\r\n",
      "   Kappa-Fission Rate                   1.12634e+08 +/- 443455\r\n",
      " ============================>     TALLY 2     <============================\r\n",
      "\r\n",
      " Incoming Energy [0, 0.005)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0453976 +/- 0.000736087\r\n",
      " Incoming Energy [0.005, 0.01)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.121742 +/- 0.00132703\r\n",
      " Incoming Energy [0.01, 0.015)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.177446 +/- 0.00165249\r\n",
      " Incoming Energy [0.015, 0.02)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.214019 +/- 0.00163567\r\n",
      " Incoming Energy [0.02, 0.025)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.237304 +/- 0.00214234\r\n",
      " Incoming Energy [0.025, 0.03)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.254086 +/- 0.00209678\r\n",
      " Incoming Energy [0.03, 0.035)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.255185 +/- 0.00217766\r\n",
      " Incoming Energy [0.035, 0.042)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.36118 +/- 0.00286327\r\n",
      " Incoming Energy [0.042, 0.05)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.390554 +/- 0.00311077\r\n",
      " Incoming Energy [0.05, 0.058)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.354584 +/- 0.00340184\r\n",
      " Incoming Energy [0.058, 0.067)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.358517 +/- 0.0027362\r\n",
      " Incoming Energy [0.067, 0.08)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.433362 +/- 0.00347373\r\n",
      " Incoming Energy [0.08, 0.1)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.505304 +/- 0.00393269\r\n",
      " Incoming Energy [0.1, 0.14)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.633421 +/- 0.00504119\r\n",
      " Incoming Energy [0.14, 0.18)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.351153 +/- 0.00338567\r\n",
      " Incoming Energy [0.18, 0.22)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.23489 +/- 0.00255307\r\n",
      " Incoming Energy [0.22, 0.25)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.136738 +/- 0.00173787\r\n",
      " Incoming Energy [0.25, 0.28)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.116409 +/- 0.00164422\r\n",
      " Incoming Energy [0.28, 0.3)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0691176 +/- 0.00123652\r\n",
      " Incoming Energy [0.3, 0.32)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0675157 +/- 0.00144419\r\n",
      " Incoming Energy [0.32, 0.35)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0906064 +/- 0.00164951\r\n",
      " Incoming Energy [0.35, 0.4)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.132927 +/- 0.00170577\r\n",
      " Incoming Energy [0.4, 0.5)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.226544 +/- 0.00227151\r\n",
      " Incoming Energy [0.5, 0.625)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.227798 +/- 0.00303226\r\n",
      " Incoming Energy [0.625, 0.78)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.224419 +/- 0.0025656\r\n",
      " Incoming Energy [0.78, 0.85)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0852675 +/- 0.00153001\r\n",
      " Incoming Energy [0.85, 0.91)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0663865 +/- 0.00138148\r\n",
      " Incoming Energy [0.91, 0.95)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0409562 +/- 0.000983402\r\n",
      " Incoming Energy [0.95, 0.972)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0221144 +/- 0.000688345\r\n",
      " Incoming Energy [0.972, 0.996)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0230513 +/- 0.000846911\r\n",
      " Incoming Energy [0.996, 1.02)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0226316 +/- 0.00070782\r\n",
      " Incoming Energy [1.02, 1.045)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0240545 +/- 0.000799047\r\n",
      " Incoming Energy [1.045, 1.071)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0236844 +/- 0.000864848\r\n",
      " Incoming Energy [1.071, 1.097)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0227409 +/- 0.000763818\r\n",
      " Incoming Energy [1.097, 1.123)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0212896 +/- 0.000663318\r\n",
      " Incoming Energy [1.123, 1.15)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.0218843 +/- 0.000730781\r\n",
      " Incoming Energy [1.15, 1.3)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.116056 +/- 0.00164544\r\n",
      " Incoming Energy [1.3, 1.5)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.138343 +/- 0.00208071\r\n",
      " Incoming Energy [1.5, 1.855)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.201696 +/- 0.00244502\r\n",
      " Incoming Energy [1.855, 2.1)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.11583 +/- 0.00196992\r\n",
      " Incoming Energy [2.1, 2.6)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.205274 +/- 0.00257039\r\n",
      " Incoming Energy [2.6, 3.3)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.22677 +/- 0.00245925\r\n",
      " Incoming Energy [3.3, 4)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.178524 +/- 0.00259963\r\n",
      " Incoming Energy [4, 9.877)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.781006 +/- 0.00494382\r\n",
      " Incoming Energy [9.877, 15.968)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.476281 +/- 0.00396726\r\n",
      " Incoming Energy [15.968, 27.7)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.522022 +/- 0.00358601\r\n",
      " Incoming Energy [27.7, 48.052)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.550981 +/- 0.00392878\r\n",
      " Incoming Energy [48.052, 75.501)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.47272 +/- 0.00427859\r\n",
      " Incoming Energy [75.501, 148.73)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.722731 +/- 0.00471049\r\n",
      " Incoming Energy [148.73, 367.26)\r\n",
      "   Total Material\r\n",
      "     Flux                                 1.01279 +/- 0.00545347\r\n",
      " Incoming Energy [367.26, 906.9)\r\n",
      "   Total Material\r\n",
      "     Flux                                 1.04418 +/- 0.00553079\r\n",
      " Incoming Energy [906.9, 1425.1)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.538633 +/- 0.00400922\r\n",
      " Incoming Energy [1425.1, 2239.5)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.537251 +/- 0.00446093\r\n",
      " Incoming Energy [2239.5, 3519.1)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.545527 +/- 0.00408891\r\n",
      " Incoming Energy [3519.1, 5530)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.553422 +/- 0.00432449\r\n",
      " Incoming Energy [5530, 9118)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.629065 +/- 0.00471568\r\n",
      " Incoming Energy [9118, 15030)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.649299 +/- 0.00469207\r\n",
      " Incoming Energy [15030, 24780)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.695172 +/- 0.00500171\r\n",
      " Incoming Energy [24780, 40850)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.747476 +/- 0.0061985\r\n",
      " Incoming Energy [40850, 67340)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.823449 +/- 0.00664113\r\n",
      " Incoming Energy [67340, 111000)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.952249 +/- 0.00742562\r\n",
      " Incoming Energy [111000, 183000)\r\n",
      "   Total Material\r\n",
      "     Flux                                 1.18162 +/- 0.0066215\r\n",
      " Incoming Energy [183000, 302500)\r\n",
      "   Total Material\r\n",
      "     Flux                                 1.53498 +/- 0.00971482\r\n",
      " Incoming Energy [302500, 500000)\r\n",
      "   Total Material\r\n",
      "     Flux                                 1.74329 +/- 0.0114778\r\n",
      " Incoming Energy [500000, 821000)\r\n",
      "   Total Material\r\n",
      "     Flux                                 2.54164 +/- 0.0154229\r\n",
      " Incoming Energy [821000, 1353000)\r\n",
      "   Total Material\r\n",
      "     Flux                                 2.43237 +/- 0.0181762\r\n",
      " Incoming Energy [1353000, 2231000)\r\n",
      "   Total Material\r\n",
      "     Flux                                 2.60684 +/- 0.0199169\r\n",
      " Incoming Energy [2231000, 3679000)\r\n",
      "   Total Material\r\n",
      "     Flux                                 2.47805 +/- 0.0227841\r\n",
      " Incoming Energy [3679000, 6065500)\r\n",
      "   Total Material\r\n",
      "     Flux                                 1.18705 +/- 0.0141164\r\n",
      " Incoming Energy [6065500, 20000000)\r\n",
      "   Total Material\r\n",
      "     Flux                                 0.303115 +/- 0.00857521\r\n"
     ]
    }
   ],
   "source": [
    "!cat tallies.out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "with openmc.StatePoint(statepoint) as sp:\n",
    "    tally = sp.get_tally(id=spectrum_tally.id)\n",
    "    tally = sp.get_tally(scores=['flux'])\n",
    "    tally = sp.get_tally(filters=spectrum_tally.filters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "values = tally.get_values(scores=['flux']).flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/anovak/projects/openmc/openmc/filter.py:1411: RuntimeWarning: divide by zero encountered in divide\n",
      "  return np.log10(self.bins[:, 1]/self.bins[:, 0])\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 18340 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "bin_boundaries = energy_filter.lethargy_bin_width\n",
    "\n",
    "#print(energy_filter.bins)\n",
    "plt.step(np.unique(energy_filter.bins)[:-1], values / bin_boundaries)\n",
    "\n",
    "plt.xscale('log')\n",
    "plt.xlabel(r'Energy (eV)')\n",
    "plt.ylabel(r'Flux ($\\frac{particle-cm}{source-particle}$)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalizing Tallies\n",
    "\n",
    "Note that the units of flux in the above plot are in $\\frac{particle-cm}{source-particle}$. As is the case with many values tallied by Monte Carlo codes, the value of the flux does not account for volume and is in terms of the number of source particles emitted. To generate this same plot in terms of absolute flux units ($\\frac{particle}{cm^{2}-s}$) we'll need to normalize this tally by:\n",
    "\n",
    "  - the volume of the region the tally covers\n",
    "  - the number of source particle emitted\n",
    "\n",
    "In this case, the volume of the region is the volume of the entire pincell, which we can obtain using a bounding box (or the dimensions we set earlier, but let's introduce the notion of a bounding box!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BoundingBox(lower_left=(-0.72135, -0.72135, -150.0), upper_right=(0.72135, 0.72135, 150.0))\n"
     ]
    }
   ],
   "source": [
    "print(root_universe.bounding_box)\n",
    "ur = root_universe.bounding_box.upper_right\n",
    "ll = root_universe.bounding_box.lower_left\n",
    "volume = np.prod(ur-ll)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determining the number of source particles per second is more complicated, however. This means computing the eV/source particle due to fission. To get the source rate, we'll need the following pieces of information:\n",
    "\n",
    "  1. the total power produced in the tally region (known a priori)\n",
    "  2. the heat produced by fission power, per source particle\n",
    "  \n",
    " To get this information we'll need to construct another tally to get additional information from the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "pincell_power = 200 # in Watts\n",
    "source_tally = openmc.Tally()\n",
    "source_tally.scores = ['kappa-fission']\n",
    "model.tallies += [source_tally]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                %%%%%%%%%%%%%%%\n",
      "                           %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                    %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                     %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ##################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ###################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ####################     %%%%%%%%%%%%%%%%%%%%%%\n",
      "                #####################     %%%%%%%%%%%%%%%%%%%%%\n",
      "                ######################     %%%%%%%%%%%%%%%%%%%%\n",
      "                #######################     %%%%%%%%%%%%%%%%%%\n",
      "                 #######################     %%%%%%%%%%%%%%%%%\n",
      "                 ######################     %%%%%%%%%%%%%%%%%\n",
      "                  ####################     %%%%%%%%%%%%%%%%%\n",
      "                    #################     %%%%%%%%%%%%%%%%%\n",
      "                     ###############     %%%%%%%%%%%%%%%%\n",
      "                       ############     %%%%%%%%%%%%%%%\n",
      "                          ########     %%%%%%%%%%%%%%\n",
      "                                      %%%%%%%%%%%\n",
      "\n",
      "                 | The OpenMC Monte Carlo Code\n",
      "       Copyright | 2011-2024 MIT, UChicago Argonne LLC, and contributors\n",
      "         License | https://docs.openmc.org/en/latest/license.html\n",
      "         Version | 0.14.1-dev\n",
      "        Git SHA1 | 14ce3cec4b388ae8b7ec5b28db57dbcbff06f147\n",
      "       Date/Time | 2024-03-29 14:34:30\n",
      "   MPI Processes | 1\n",
      "  OpenMP Threads | 8\n",
      "\n",
      " Reading model XML file 'model.xml' ...\n",
      " WARNING: Other XML file input(s) are present. These files may be ignored in\n",
      "          favor of the model.xml file.\n",
      " Reading cross sections XML file...\n",
      " Reading Zr90 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr90.h5\n",
      " Reading Zr91 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr91.h5\n",
      " Reading Zr92 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr92.h5\n",
      " Reading Zr94 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr94.h5\n",
      " Reading Zr96 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr96.h5\n",
      " Reading U234 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U234.h5\n",
      " Reading U235 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U235.h5\n",
      " Reading U238 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U238.h5\n",
      " Reading U236 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U236.h5\n",
      " Reading O16 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O16.h5\n",
      " Reading H1 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H1.h5\n",
      " Reading H2 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H2.h5\n",
      " Reading O17 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O17.h5\n",
      " Reading Pu239 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Pu239.h5\n",
      " Reading Si28 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Si28.h5\n",
      " Reading c_H_in_H2O from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/c_H_in_H2O.h5\n",
      " Minimum neutron data temperature: 250 K\n",
      " Maximum neutron data temperature: 2500 K\n",
      " Preparing distributed cell instances...\n",
      " Writing summary.h5 file...\n",
      " Maximum neutron transport energy: 20000000 eV for Zr90\n",
      " Initializing source particles...\n",
      "\n",
      " ====================>     K EIGENVALUE SIMULATION     <====================\n",
      "\n",
      "  Bat./Gen.      k            Average k\n",
      "  =========   ========   ====================\n",
      "        1/1    1.50315\n",
      "        2/1    1.44086\n",
      "        3/1    1.46540\n",
      "        4/1    1.37107\n",
      "        5/1    1.33015\n",
      "        6/1    1.48464\n",
      "        7/1    1.47998\n",
      "        8/1    1.28684\n",
      "        9/1    1.45556\n",
      "       10/1    1.40944\n",
      "       11/1    1.39122\n",
      "       12/1    1.32307    1.35714 +/- 0.03407\n",
      "       13/1    1.52588    1.41339 +/- 0.05959\n",
      "       14/1    1.46026    1.42511 +/- 0.04373\n",
      "       15/1    1.41125    1.42233 +/- 0.03399\n",
      "       16/1    1.45798    1.42827 +/- 0.02838\n",
      "       17/1    1.46100    1.43295 +/- 0.02444\n",
      "       18/1    1.42515    1.43197 +/- 0.02119\n",
      "       19/1    1.54914    1.44499 +/- 0.02277\n",
      "       20/1    1.43623    1.44412 +/- 0.02039\n",
      "       21/1    1.40115    1.44021 +/- 0.01885\n",
      "       22/1    1.43731    1.43997 +/- 0.01721\n",
      "       23/1    1.41573    1.43811 +/- 0.01594\n",
      "       24/1    1.46427    1.43997 +/- 0.01487\n",
      "       25/1    1.45526    1.44099 +/- 0.01389\n",
      "       26/1    1.33933    1.43464 +/- 0.01446\n",
      "       27/1    1.42147    1.43386 +/- 0.01360\n",
      "       28/1    1.36027    1.42978 +/- 0.01346\n",
      "       29/1    1.46546    1.43165 +/- 0.01287\n",
      "       30/1    1.45070    1.43261 +/- 0.01225\n",
      "       31/1    1.48473    1.43509 +/- 0.01191\n",
      "       32/1    1.47264    1.43680 +/- 0.01149\n",
      "       33/1    1.45237    1.43747 +/- 0.01100\n",
      "       34/1    1.41749    1.43664 +/- 0.01056\n",
      "       35/1    1.47971    1.43836 +/- 0.01027\n",
      "       36/1    1.45378    1.43896 +/- 0.00989\n",
      "       37/1    1.35765    1.43594 +/- 0.00998\n",
      "       38/1    1.40397    1.43480 +/- 0.00969\n",
      "       39/1    1.41697    1.43419 +/- 0.00937\n",
      "       40/1    1.40400    1.43318 +/- 0.00910\n",
      "       41/1    1.43270    1.43317 +/- 0.00881\n",
      "       42/1    1.36535    1.43105 +/- 0.00879\n",
      "       43/1    1.43240    1.43109 +/- 0.00851\n",
      "       44/1    1.46349    1.43204 +/- 0.00832\n",
      "       45/1    1.32866    1.42909 +/- 0.00860\n",
      "       46/1    1.32249    1.42613 +/- 0.00886\n",
      "       47/1    1.41734    1.42589 +/- 0.00862\n",
      "       48/1    1.39802    1.42516 +/- 0.00843\n",
      "       49/1    1.54696    1.42828 +/- 0.00878\n",
      "       50/1    1.44601    1.42872 +/- 0.00857\n",
      "       51/1    1.37546    1.42742 +/- 0.00846\n",
      "       52/1    1.48833    1.42887 +/- 0.00838\n",
      "       53/1    1.36864    1.42747 +/- 0.00830\n",
      "       54/1    1.45539    1.42811 +/- 0.00814\n",
      "       55/1    1.47463    1.42914 +/- 0.00802\n",
      "       56/1    1.43531    1.42927 +/- 0.00785\n",
      "       57/1    1.37650    1.42815 +/- 0.00776\n",
      "       58/1    1.43163    1.42822 +/- 0.00760\n",
      "       59/1    1.39161    1.42748 +/- 0.00748\n",
      "       60/1    1.48475    1.42862 +/- 0.00742\n",
      "       61/1    1.47918    1.42961 +/- 0.00734\n",
      "       62/1    1.47997    1.43058 +/- 0.00726\n",
      "       63/1    1.41811    1.43035 +/- 0.00712\n",
      "       64/1    1.41011    1.42997 +/- 0.00700\n",
      "       65/1    1.44239    1.43020 +/- 0.00688\n",
      "       66/1    1.42912    1.43018 +/- 0.00675\n",
      "       67/1    1.39312    1.42953 +/- 0.00666\n",
      "       68/1    1.47339    1.43028 +/- 0.00659\n",
      "       69/1    1.38559    1.42953 +/- 0.00652\n",
      "       70/1    1.42211    1.42940 +/- 0.00641\n",
      "       71/1    1.35744    1.42822 +/- 0.00642\n",
      "       72/1    1.40876    1.42791 +/- 0.00632\n",
      "       73/1    1.53582    1.42962 +/- 0.00645\n",
      "       74/1    1.44605    1.42988 +/- 0.00636\n",
      "       75/1    1.47428    1.43056 +/- 0.00629\n",
      "       76/1    1.37855    1.42977 +/- 0.00625\n",
      "       77/1    1.39439    1.42925 +/- 0.00618\n",
      "       78/1    1.49117    1.43016 +/- 0.00615\n",
      "       79/1    1.34815    1.42897 +/- 0.00618\n",
      "       80/1    1.30610    1.42721 +/- 0.00634\n",
      "       81/1    1.44634    1.42748 +/- 0.00625\n",
      "       82/1    1.41633    1.42733 +/- 0.00617\n",
      "       83/1    1.45974    1.42777 +/- 0.00610\n",
      "       84/1    1.45538    1.42814 +/- 0.00603\n",
      "       85/1    1.45563    1.42851 +/- 0.00596\n",
      "       86/1    1.37790    1.42785 +/- 0.00592\n",
      "       87/1    1.43465    1.42793 +/- 0.00584\n",
      "       88/1    1.32547    1.42662 +/- 0.00591\n",
      "       89/1    1.37892    1.42602 +/- 0.00587\n",
      "       90/1    1.47858    1.42667 +/- 0.00583\n",
      "       91/1    1.34703    1.42569 +/- 0.00584\n",
      "       92/1    1.47197    1.42625 +/- 0.00580\n",
      "       93/1    1.48311    1.42694 +/- 0.00577\n",
      "       94/1    1.39298    1.42653 +/- 0.00571\n",
      "       95/1    1.52197    1.42766 +/- 0.00576\n",
      "       96/1    1.35930    1.42686 +/- 0.00574\n",
      "       97/1    1.51529    1.42788 +/- 0.00577\n",
      "       98/1    1.45297    1.42816 +/- 0.00571\n",
      "       99/1    1.40311    1.42788 +/- 0.00565\n",
      "      100/1    1.42293    1.42783 +/- 0.00559\n",
      " Creating state point statepoint.100.h5...\n",
      "\n",
      " =======================>     TIMING STATISTICS     <=======================\n",
      "\n",
      " Total time for initialization     = 2.1963e+00 seconds\n",
      "   Reading cross sections          = 2.1396e+00 seconds\n",
      " Total time in simulation          = 1.9032e+00 seconds\n",
      "   Time in transport only          = 1.8521e+00 seconds\n",
      "   Time in inactive batches        = 1.4214e-01 seconds\n",
      "   Time in active batches          = 1.7610e+00 seconds\n",
      "   Time synchronizing fission bank = 6.1936e-03 seconds\n",
      "     Sampling source sites         = 4.8027e-03 seconds\n",
      "     SEND/RECV source sites        = 5.7378e-04 seconds\n",
      "   Time accumulating tallies       = 2.7497e-02 seconds\n",
      "   Time writing statepoints        = 1.2327e-02 seconds\n",
      " Total time for finalization       = 1.1213e-03 seconds\n",
      " Total time elapsed                = 4.1105e+00 seconds\n",
      " Calculation Rate (inactive)       = 70353.3 particles/second\n",
      " Calculation Rate (active)         = 51106.9 particles/second\n",
      "\n",
      " ============================>     RESULTS     <============================\n",
      "\n",
      " k-effective (Collision)     = 1.42881 +/- 0.00466\n",
      " k-effective (Track-length)  = 1.42783 +/- 0.00559\n",
      " k-effective (Absorption)    = 1.42667 +/- 0.00316\n",
      " Combined k-effective        = 1.42696 +/- 0.00299\n",
      " Leakage Fraction            = 0.00163 +/- 0.00017\n",
      "\n"
     ]
    }
   ],
   "source": [
    "statepoint = model.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "with openmc.StatePoint(statepoint) as sp:\n",
    "    out = sp.get_tally(id=source_tally.id)\n",
    "    heating = out.get_values(scores=['kappa-fission']).flatten()[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The combination of the following tally values and power provide us with the source normalization needed as follows:\n",
    "\n",
    "\n",
    "$$ \\text{neutron source} [\\frac{n}{s}] = \\text{power} [\\frac{J}{s}] \\times \\frac{1}{1.6\\times 10^{-19}} [\\frac{eV}{J}] \\times \\frac{1}{\\text{heat per fission} [\\frac{eV}{source}]} $$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Neutron source (n/s):  11097897947492.568\n"
     ]
    }
   ],
   "source": [
    "J_to_eV = 1 / 1.6e-19\n",
    "neutron_source = pincell_power * J_to_eV * (1 / heating)\n",
    "print('Neutron source (n/s): ', neutron_source)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "normalized_spectrum = neutron_source * values / volume"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 17568 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.step(np.unique(energy_filter.bins)[:-1], normalized_spectrum / bin_boundaries)\n",
    "plt.xscale('log')\n",
    "plt.xlabel(r'Energy (eV)')\n",
    "plt.ylabel(r'Flux ($\\frac{particle}{cm^{-2} s}$)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 3: Reaction Types by Material\n",
    "\n",
    "Looking at the different reaction types by material will require a material filter and the set of reaction types we want to score. For this example, we'll be scoring absorption, scattering and fission in each material.\n",
    "\n",
    "To start, we'll create a material filter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Material\n",
      "\tID             =\t3\n",
      "\tName           =\tuo2\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t10.0 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tTrue\n",
      "\tS(a,b) Tables  \n",
      "\tNuclides       \n",
      "\tU234           =\t0.0003166930253944235 [ao]\n",
      "\tU235           =\t0.03543164439454172 [ao]\n",
      "\tU238           =\t0.964089368630351 [ao]\n",
      "\tU236           =\t0.00016229394971280895 [ao]\n",
      "\tO16            =\t2.0          [ao]\n",
      ", Material\n",
      "\tID             =\t1\n",
      "\tName           =\tzirconium\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t6.5 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tNuclides       \n",
      "\tZr90           =\t0.5145       [ao]\n",
      "\tZr91           =\t0.1122       [ao]\n",
      "\tZr92           =\t0.1715       [ao]\n",
      "\tZr94           =\t0.1738       [ao]\n",
      "\tZr96           =\t0.028        [ao]\n",
      ", Material\n",
      "\tID             =\t4\n",
      "\tName           =\twater\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t1.0 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tFalse\n",
      "\tS(a,b) Tables  \n",
      "\tS(a,b)         =\t('c_H_in_H2O', 1.0)\n",
      "\tNuclides       \n",
      "\tH1             =\t1.99968852   [ao]\n",
      "\tH2             =\t0.00031148   [ao]\n",
      "\tO16            =\t0.999621     [ao]\n",
      "\tO17            =\t0.000379     [ao]\n",
      ", Material\n",
      "\tID             =\t5\n",
      "\tName           =\tnew_fuel\n",
      "\tTemperature    =\tNone\n",
      "\tDensity        =\t9.0 [g/cm3]\n",
      "\tVolume         =\tNone [cm^3]\n",
      "\tDepletable     =\tTrue\n",
      "\tS(a,b) Tables  \n",
      "\tNuclides       \n",
      "\tPu239          =\t1.0          [ao]\n",
      "\tSi28           =\t2.0          [ao]\n",
      "]\n"
     ]
    }
   ],
   "source": [
    "print(model.materials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "material_filter = openmc.MaterialFilter(model.materials)\n",
    "material_tally = openmc.Tally()\n",
    "material_tally.filters = [material_filter]\n",
    "material_tally.scores = ['absorption', 'scatter', 'fission']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.tallies += [material_tally]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                %%%%%%%%%%%%%%%\n",
      "                           %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                    %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                     %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ##################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ###################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ####################     %%%%%%%%%%%%%%%%%%%%%%\n",
      "                #####################     %%%%%%%%%%%%%%%%%%%%%\n",
      "                ######################     %%%%%%%%%%%%%%%%%%%%\n",
      "                #######################     %%%%%%%%%%%%%%%%%%\n",
      "                 #######################     %%%%%%%%%%%%%%%%%\n",
      "                 ######################     %%%%%%%%%%%%%%%%%\n",
      "                  ####################     %%%%%%%%%%%%%%%%%\n",
      "                    #################     %%%%%%%%%%%%%%%%%\n",
      "                     ###############     %%%%%%%%%%%%%%%%\n",
      "                       ############     %%%%%%%%%%%%%%%\n",
      "                          ########     %%%%%%%%%%%%%%\n",
      "                                      %%%%%%%%%%%\n",
      "\n",
      "                 | The OpenMC Monte Carlo Code\n",
      "       Copyright | 2011-2024 MIT, UChicago Argonne LLC, and contributors\n",
      "         License | https://docs.openmc.org/en/latest/license.html\n",
      "         Version | 0.14.1-dev\n",
      "        Git SHA1 | 14ce3cec4b388ae8b7ec5b28db57dbcbff06f147\n",
      "       Date/Time | 2024-03-29 14:34:35\n",
      "   MPI Processes | 1\n",
      "  OpenMP Threads | 8\n",
      "\n",
      " Reading model XML file 'model.xml' ...\n",
      " WARNING: Other XML file input(s) are present. These files may be ignored in\n",
      "          favor of the model.xml file.\n",
      " Reading cross sections XML file...\n",
      " Reading Zr90 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr90.h5\n",
      " Reading Zr91 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr91.h5\n",
      " Reading Zr92 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr92.h5\n",
      " Reading Zr94 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr94.h5\n",
      " Reading Zr96 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr96.h5\n",
      " Reading U234 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U234.h5\n",
      " Reading U235 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U235.h5\n",
      " Reading U238 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U238.h5\n",
      " Reading U236 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U236.h5\n",
      " Reading O16 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O16.h5\n",
      " Reading H1 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H1.h5\n",
      " Reading H2 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H2.h5\n",
      " Reading O17 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O17.h5\n",
      " Reading Pu239 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Pu239.h5\n",
      " Reading Si28 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Si28.h5\n",
      " Reading c_H_in_H2O from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/c_H_in_H2O.h5\n",
      " Minimum neutron data temperature: 250 K\n",
      " Maximum neutron data temperature: 2500 K\n",
      " Preparing distributed cell instances...\n",
      " Writing summary.h5 file...\n",
      " Maximum neutron transport energy: 20000000 eV for Zr90\n",
      " Initializing source particles...\n",
      "\n",
      " ====================>     K EIGENVALUE SIMULATION     <====================\n",
      "\n",
      "  Bat./Gen.      k            Average k\n",
      "  =========   ========   ====================\n",
      "        1/1    1.50315\n",
      "        2/1    1.44086\n",
      "        3/1    1.46540\n",
      "        4/1    1.37107\n",
      "        5/1    1.33015\n",
      "        6/1    1.48464\n",
      "        7/1    1.47998\n",
      "        8/1    1.28684\n",
      "        9/1    1.45556\n",
      "       10/1    1.40944\n",
      "       11/1    1.39122\n",
      "       12/1    1.32307    1.35714 +/- 0.03407\n",
      "       13/1    1.52588    1.41339 +/- 0.05959\n",
      "       14/1    1.46026    1.42511 +/- 0.04373\n",
      "       15/1    1.41125    1.42233 +/- 0.03399\n",
      "       16/1    1.45798    1.42827 +/- 0.02838\n",
      "       17/1    1.46100    1.43295 +/- 0.02444\n",
      "       18/1    1.42515    1.43197 +/- 0.02119\n",
      "       19/1    1.54914    1.44499 +/- 0.02277\n",
      "       20/1    1.43623    1.44412 +/- 0.02039\n",
      "       21/1    1.40115    1.44021 +/- 0.01885\n",
      "       22/1    1.43731    1.43997 +/- 0.01721\n",
      "       23/1    1.41573    1.43811 +/- 0.01594\n",
      "       24/1    1.46427    1.43997 +/- 0.01487\n",
      "       25/1    1.45526    1.44099 +/- 0.01389\n",
      "       26/1    1.33933    1.43464 +/- 0.01446\n",
      "       27/1    1.42147    1.43386 +/- 0.01360\n",
      "       28/1    1.36027    1.42978 +/- 0.01346\n",
      "       29/1    1.46546    1.43165 +/- 0.01287\n",
      "       30/1    1.45070    1.43261 +/- 0.01225\n",
      "       31/1    1.48473    1.43509 +/- 0.01191\n",
      "       32/1    1.47264    1.43680 +/- 0.01149\n",
      "       33/1    1.45237    1.43747 +/- 0.01100\n",
      "       34/1    1.41749    1.43664 +/- 0.01056\n",
      "       35/1    1.47971    1.43836 +/- 0.01027\n",
      "       36/1    1.45378    1.43896 +/- 0.00989\n",
      "       37/1    1.35765    1.43594 +/- 0.00998\n",
      "       38/1    1.40397    1.43480 +/- 0.00969\n",
      "       39/1    1.41697    1.43419 +/- 0.00937\n",
      "       40/1    1.40400    1.43318 +/- 0.00910\n",
      "       41/1    1.43270    1.43317 +/- 0.00881\n",
      "       42/1    1.36535    1.43105 +/- 0.00879\n",
      "       43/1    1.43240    1.43109 +/- 0.00851\n",
      "       44/1    1.46349    1.43204 +/- 0.00832\n",
      "       45/1    1.32866    1.42909 +/- 0.00860\n",
      "       46/1    1.32249    1.42613 +/- 0.00886\n",
      "       47/1    1.41734    1.42589 +/- 0.00862\n",
      "       48/1    1.39802    1.42516 +/- 0.00843\n",
      "       49/1    1.54696    1.42828 +/- 0.00878\n",
      "       50/1    1.44601    1.42872 +/- 0.00857\n",
      "       51/1    1.37546    1.42742 +/- 0.00846\n",
      "       52/1    1.48833    1.42887 +/- 0.00838\n",
      "       53/1    1.36864    1.42747 +/- 0.00830\n",
      "       54/1    1.45539    1.42811 +/- 0.00814\n",
      "       55/1    1.47463    1.42914 +/- 0.00802\n",
      "       56/1    1.43531    1.42927 +/- 0.00785\n",
      "       57/1    1.37650    1.42815 +/- 0.00776\n",
      "       58/1    1.43163    1.42822 +/- 0.00760\n",
      "       59/1    1.39161    1.42748 +/- 0.00748\n",
      "       60/1    1.48475    1.42862 +/- 0.00742\n",
      "       61/1    1.47918    1.42961 +/- 0.00734\n",
      "       62/1    1.47997    1.43058 +/- 0.00726\n",
      "       63/1    1.41811    1.43035 +/- 0.00712\n",
      "       64/1    1.41011    1.42997 +/- 0.00700\n",
      "       65/1    1.44239    1.43020 +/- 0.00688\n",
      "       66/1    1.42912    1.43018 +/- 0.00675\n",
      "       67/1    1.39312    1.42953 +/- 0.00666\n",
      "       68/1    1.47339    1.43028 +/- 0.00659\n",
      "       69/1    1.38559    1.42953 +/- 0.00652\n",
      "       70/1    1.42211    1.42940 +/- 0.00641\n",
      "       71/1    1.35744    1.42822 +/- 0.00642\n",
      "       72/1    1.40876    1.42791 +/- 0.00632\n",
      "       73/1    1.53582    1.42962 +/- 0.00645\n",
      "       74/1    1.44605    1.42988 +/- 0.00636\n",
      "       75/1    1.47428    1.43056 +/- 0.00629\n",
      "       76/1    1.37855    1.42977 +/- 0.00625\n",
      "       77/1    1.39439    1.42925 +/- 0.00618\n",
      "       78/1    1.49117    1.43016 +/- 0.00615\n",
      "       79/1    1.34815    1.42897 +/- 0.00618\n",
      "       80/1    1.30610    1.42721 +/- 0.00634\n",
      "       81/1    1.44634    1.42748 +/- 0.00625\n",
      "       82/1    1.41633    1.42733 +/- 0.00617\n",
      "       83/1    1.45974    1.42777 +/- 0.00610\n",
      "       84/1    1.45538    1.42814 +/- 0.00603\n",
      "       85/1    1.45563    1.42851 +/- 0.00596\n",
      "       86/1    1.37790    1.42785 +/- 0.00592\n",
      "       87/1    1.43465    1.42793 +/- 0.00584\n",
      "       88/1    1.32547    1.42662 +/- 0.00591\n",
      "       89/1    1.37892    1.42602 +/- 0.00587\n",
      "       90/1    1.47858    1.42667 +/- 0.00583\n",
      "       91/1    1.34703    1.42569 +/- 0.00584\n",
      "       92/1    1.47197    1.42625 +/- 0.00580\n",
      "       93/1    1.48311    1.42694 +/- 0.00577\n",
      "       94/1    1.39298    1.42653 +/- 0.00571\n",
      "       95/1    1.52197    1.42766 +/- 0.00576\n",
      "       96/1    1.35930    1.42686 +/- 0.00574\n",
      "       97/1    1.51529    1.42788 +/- 0.00577\n",
      "       98/1    1.45297    1.42816 +/- 0.00571\n",
      "       99/1    1.40311    1.42788 +/- 0.00565\n",
      "      100/1    1.42293    1.42783 +/- 0.00559\n",
      " Creating state point statepoint.100.h5...\n",
      "\n",
      " =======================>     TIMING STATISTICS     <=======================\n",
      "\n",
      " Total time for initialization     = 1.9115e+00 seconds\n",
      "   Reading cross sections          = 1.8639e+00 seconds\n",
      " Total time in simulation          = 1.9756e+00 seconds\n",
      "   Time in transport only          = 1.9494e+00 seconds\n",
      "   Time in inactive batches        = 1.4190e-01 seconds\n",
      "   Time in active batches          = 1.8337e+00 seconds\n",
      "   Time synchronizing fission bank = 5.9623e-03 seconds\n",
      "     Sampling source sites         = 4.6432e-03 seconds\n",
      "     SEND/RECV source sites        = 5.3352e-04 seconds\n",
      "   Time accumulating tallies       = 9.3554e-03 seconds\n",
      "   Time writing statepoints        = 6.9239e-03 seconds\n",
      " Total time for finalization       = 5.6553e-04 seconds\n",
      " Total time elapsed                = 3.8974e+00 seconds\n",
      " Calculation Rate (inactive)       = 70471.2 particles/second\n",
      " Calculation Rate (active)         = 49081 particles/second\n",
      "\n",
      " ============================>     RESULTS     <============================\n",
      "\n",
      " k-effective (Collision)     = 1.42881 +/- 0.00466\n",
      " k-effective (Track-length)  = 1.42783 +/- 0.00559\n",
      " k-effective (Absorption)    = 1.42667 +/- 0.00316\n",
      " Combined k-effective        = 1.42696 +/- 0.00299\n",
      " Leakage Fraction            = 0.00163 +/- 0.00017\n",
      "\n"
     ]
    }
   ],
   "source": [
    "statepoint = model.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll gather information from the statepoint file about each score we applied to the tally. With multiple scores and materials, we'll use a Pandas data frame to view the results in a more coherent manner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "with openmc.StatePoint(statepoint) as sp:\n",
    "    tally = sp.get_tally(filters=[material_filter])\n",
    "    absorption = tally.get_values(scores=['absorption']).flatten()\n",
    "    scatter = tally.get_values(scores=['scatter']).flatten()\n",
    "    fission = tally.get_values(scores=['fission']).flatten()\n",
    "    df = tally.get_pandas_dataframe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>material</th>\n",
       "      <th>nuclide</th>\n",
       "      <th>score</th>\n",
       "      <th>mean</th>\n",
       "      <th>std. dev.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.931082</td>\n",
       "      <td>0.003152</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>4.488983</td>\n",
       "      <td>0.007124</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.581783</td>\n",
       "      <td>0.002293</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.008032</td>\n",
       "      <td>0.000065</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>1.146501</td>\n",
       "      <td>0.001893</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>4</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.059669</td>\n",
       "      <td>0.000244</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>4</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>24.461243</td>\n",
       "      <td>0.049262</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>4</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>5</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>5</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>5</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    material nuclide       score     mean  std. dev.\n",
       "0          3   total  absorption 9.31e-01   3.15e-03\n",
       "1          3   total     scatter 4.49e+00   7.12e-03\n",
       "2          3   total     fission 5.82e-01   2.29e-03\n",
       "3          1   total  absorption 8.03e-03   6.48e-05\n",
       "4          1   total     scatter 1.15e+00   1.89e-03\n",
       "5          1   total     fission 0.00e+00   0.00e+00\n",
       "6          4   total  absorption 5.97e-02   2.44e-04\n",
       "7          4   total     scatter 2.45e+01   4.93e-02\n",
       "8          4   total     fission 0.00e+00   0.00e+00\n",
       "9          5   total  absorption 0.00e+00   0.00e+00\n",
       "10         5   total     scatter 0.00e+00   0.00e+00\n",
       "11         5   total     fission 0.00e+00   0.00e+00"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>material</th>\n",
       "      <th>nuclide</th>\n",
       "      <th>score</th>\n",
       "      <th>mean</th>\n",
       "      <th>std. dev.</th>\n",
       "      <th>normalized-mean (rxn/s)</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.931082</td>\n",
       "      <td>0.003152</td>\n",
       "      <td>1.033305e+13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>4.488983</td>\n",
       "      <td>0.007124</td>\n",
       "      <td>4.981827e+13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.581783</td>\n",
       "      <td>0.002293</td>\n",
       "      <td>6.456563e+12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.008032</td>\n",
       "      <td>0.000065</td>\n",
       "      <td>8.913624e+10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>1.146501</td>\n",
       "      <td>0.001893</td>\n",
       "      <td>1.272375e+13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>4</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.059669</td>\n",
       "      <td>0.000244</td>\n",
       "      <td>6.622006e+11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>4</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>24.461243</td>\n",
       "      <td>0.049262</td>\n",
       "      <td>2.714684e+14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>4</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>5</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>5</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>5</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    material nuclide       score     mean  std. dev.  normalized-mean (rxn/s)\n",
       "0          3   total  absorption 9.31e-01   3.15e-03                 1.03e+13\n",
       "1          3   total     scatter 4.49e+00   7.12e-03                 4.98e+13\n",
       "2          3   total     fission 5.82e-01   2.29e-03                 6.46e+12\n",
       "3          1   total  absorption 8.03e-03   6.48e-05                 8.91e+10\n",
       "4          1   total     scatter 1.15e+00   1.89e-03                 1.27e+13\n",
       "5          1   total     fission 0.00e+00   0.00e+00                 0.00e+00\n",
       "6          4   total  absorption 5.97e-02   2.44e-04                 6.62e+11\n",
       "7          4   total     scatter 2.45e+01   4.93e-02                 2.71e+14\n",
       "8          4   total     fission 0.00e+00   0.00e+00                 0.00e+00\n",
       "9          5   total  absorption 0.00e+00   0.00e+00                 0.00e+00\n",
       "10         5   total     scatter 0.00e+00   0.00e+00                 0.00e+00\n",
       "11         5   total     fission 0.00e+00   0.00e+00                 0.00e+00"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['normalized-mean (rxn/s)'] = neutron_source * df['mean']\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll add a new entry in the dataframe for our material names to make plotting easier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>material</th>\n",
       "      <th>nuclide</th>\n",
       "      <th>score</th>\n",
       "      <th>mean</th>\n",
       "      <th>std. dev.</th>\n",
       "      <th>normalized-mean (rxn/s)</th>\n",
       "      <th>mat_name</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>3</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.931082</td>\n",
       "      <td>0.003152</td>\n",
       "      <td>1.033305e+13</td>\n",
       "      <td>uo2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>4.488983</td>\n",
       "      <td>0.007124</td>\n",
       "      <td>4.981827e+13</td>\n",
       "      <td>uo2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.581783</td>\n",
       "      <td>0.002293</td>\n",
       "      <td>6.456563e+12</td>\n",
       "      <td>uo2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.008032</td>\n",
       "      <td>0.000065</td>\n",
       "      <td>8.913624e+10</td>\n",
       "      <td>zirconium</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>1.146501</td>\n",
       "      <td>0.001893</td>\n",
       "      <td>1.272375e+13</td>\n",
       "      <td>zirconium</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>zirconium</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>4</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.059669</td>\n",
       "      <td>0.000244</td>\n",
       "      <td>6.622006e+11</td>\n",
       "      <td>water</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>4</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>24.461243</td>\n",
       "      <td>0.049262</td>\n",
       "      <td>2.714684e+14</td>\n",
       "      <td>water</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>4</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>water</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>5</td>\n",
       "      <td>total</td>\n",
       "      <td>absorption</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>5</td>\n",
       "      <td>total</td>\n",
       "      <td>scatter</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>5</td>\n",
       "      <td>total</td>\n",
       "      <td>fission</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    material nuclide       score     mean  std. dev.  normalized-mean (rxn/s)  \\\n",
       "0          3   total  absorption 9.31e-01   3.15e-03                 1.03e+13   \n",
       "1          3   total     scatter 4.49e+00   7.12e-03                 4.98e+13   \n",
       "2          3   total     fission 5.82e-01   2.29e-03                 6.46e+12   \n",
       "3          1   total  absorption 8.03e-03   6.48e-05                 8.91e+10   \n",
       "4          1   total     scatter 1.15e+00   1.89e-03                 1.27e+13   \n",
       "5          1   total     fission 0.00e+00   0.00e+00                 0.00e+00   \n",
       "6          4   total  absorption 5.97e-02   2.44e-04                 6.62e+11   \n",
       "7          4   total     scatter 2.45e+01   4.93e-02                 2.71e+14   \n",
       "8          4   total     fission 0.00e+00   0.00e+00                 0.00e+00   \n",
       "9          5   total  absorption 0.00e+00   0.00e+00                 0.00e+00   \n",
       "10         5   total     scatter 0.00e+00   0.00e+00                 0.00e+00   \n",
       "11         5   total     fission 0.00e+00   0.00e+00                 0.00e+00   \n",
       "\n",
       "     mat_name  \n",
       "0         uo2  \n",
       "1         uo2  \n",
       "2         uo2  \n",
       "3   zirconium  \n",
       "4   zirconium  \n",
       "5   zirconium  \n",
       "6       water  \n",
       "7       water  \n",
       "8       water  \n",
       "9         NaN  \n",
       "10        NaN  \n",
       "11        NaN  "
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get all materials from the geometry\n",
    "materials = model.geometry.get_all_materials()\n",
    "\n",
    "# set names based on matching material IDs\n",
    "for mat_id, material in materials.items():\n",
    "    df.loc[df['material'] == mat_id, 'mat_name'] = material.name\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 4: Mesh Tallies\n",
    "\n",
    "In order to restrict a tally to a particular region in space, OpenMC supports (i) cell, (ii) structured mesh, and (iii) unstructured mesh tallies. For unstructured mesh tallies, you need to compile OpenMC with libMesh enabled. We don't have this dependency set up on the Collab instance, so we'll just work with structured mesh tallies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = openmc.RegularMesh()\n",
    "mesh.lower_left = (-pitch/2, -pitch/2)\n",
    "mesh.upper_right = (pitch/2, pitch/2)\n",
    "mesh.dimension = (50, 50)\n",
    "\n",
    "mesh_filter = openmc.MeshFilter(mesh)\n",
    "heat = openmc.Tally()\n",
    "heat.scores = ['kappa-fission']\n",
    "heat.filters = [mesh_filter]\n",
    "model.tallies += [heat]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                %%%%%%%%%%%%%%%\n",
      "                           %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                    %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                                     %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ##################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ###################     %%%%%%%%%%%%%%%%%%%%%%%\n",
      "                ####################     %%%%%%%%%%%%%%%%%%%%%%\n",
      "                #####################     %%%%%%%%%%%%%%%%%%%%%\n",
      "                ######################     %%%%%%%%%%%%%%%%%%%%\n",
      "                #######################     %%%%%%%%%%%%%%%%%%\n",
      "                 #######################     %%%%%%%%%%%%%%%%%\n",
      "                 ######################     %%%%%%%%%%%%%%%%%\n",
      "                  ####################     %%%%%%%%%%%%%%%%%\n",
      "                    #################     %%%%%%%%%%%%%%%%%\n",
      "                     ###############     %%%%%%%%%%%%%%%%\n",
      "                       ############     %%%%%%%%%%%%%%%\n",
      "                          ########     %%%%%%%%%%%%%%\n",
      "                                      %%%%%%%%%%%\n",
      "\n",
      "                 | The OpenMC Monte Carlo Code\n",
      "       Copyright | 2011-2024 MIT, UChicago Argonne LLC, and contributors\n",
      "         License | https://docs.openmc.org/en/latest/license.html\n",
      "         Version | 0.14.1-dev\n",
      "        Git SHA1 | 14ce3cec4b388ae8b7ec5b28db57dbcbff06f147\n",
      "       Date/Time | 2024-03-29 14:40:19\n",
      "   MPI Processes | 1\n",
      "  OpenMP Threads | 8\n",
      "\n",
      " Reading model XML file 'model.xml' ...\n",
      " WARNING: Other XML file input(s) are present. These files may be ignored in\n",
      "          favor of the model.xml file.\n",
      " Reading cross sections XML file...\n",
      " Reading Zr90 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr90.h5\n",
      " Reading Zr91 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr91.h5\n",
      " Reading Zr92 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr92.h5\n",
      " Reading Zr94 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr94.h5\n",
      " Reading Zr96 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Zr96.h5\n",
      " Reading U234 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U234.h5\n",
      " Reading U235 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U235.h5\n",
      " Reading U238 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U238.h5\n",
      " Reading U236 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/U236.h5\n",
      " Reading O16 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O16.h5\n",
      " Reading H1 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H1.h5\n",
      " Reading H2 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/H2.h5\n",
      " Reading O17 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/O17.h5\n",
      " Reading Pu239 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Pu239.h5\n",
      " Reading Si28 from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/Si28.h5\n",
      " Reading c_H_in_H2O from\n",
      " /Users/anovak/projects/cross_sections/endfb-vii.1-hdf5/neutron/c_H_in_H2O.h5\n",
      " Minimum neutron data temperature: 250 K\n",
      " Maximum neutron data temperature: 2500 K\n",
      " Preparing distributed cell instances...\n",
      " Writing summary.h5 file...\n",
      " Maximum neutron transport energy: 20000000 eV for Zr90\n",
      " Initializing source particles...\n",
      "\n",
      " ====================>     K EIGENVALUE SIMULATION     <====================\n",
      "\n",
      "  Bat./Gen.      k            Average k\n",
      "  =========   ========   ====================\n",
      "        1/1    1.50315\n",
      "        2/1    1.44086\n",
      "        3/1    1.46540\n",
      "        4/1    1.37107\n",
      "        5/1    1.33015\n",
      "        6/1    1.48464\n",
      "        7/1    1.47998\n",
      "        8/1    1.28684\n",
      "        9/1    1.45556\n",
      "       10/1    1.40944\n",
      "       11/1    1.39122\n",
      "       12/1    1.32307    1.35714 +/- 0.03407\n",
      "       13/1    1.52588    1.41339 +/- 0.05959\n",
      "       14/1    1.46026    1.42511 +/- 0.04373\n",
      "       15/1    1.41125    1.42233 +/- 0.03399\n",
      "       16/1    1.45798    1.42827 +/- 0.02838\n",
      "       17/1    1.46100    1.43295 +/- 0.02444\n",
      "       18/1    1.42515    1.43197 +/- 0.02119\n",
      "       19/1    1.54914    1.44499 +/- 0.02277\n",
      "       20/1    1.43623    1.44412 +/- 0.02039\n",
      "       21/1    1.40115    1.44021 +/- 0.01885\n",
      "       22/1    1.43731    1.43997 +/- 0.01721\n",
      "       23/1    1.41573    1.43811 +/- 0.01594\n",
      "       24/1    1.46427    1.43997 +/- 0.01487\n",
      "       25/1    1.45526    1.44099 +/- 0.01389\n",
      "       26/1    1.33933    1.43464 +/- 0.01446\n",
      "       27/1    1.42147    1.43386 +/- 0.01360\n",
      "       28/1    1.36027    1.42978 +/- 0.01346\n",
      "       29/1    1.46546    1.43165 +/- 0.01287\n",
      "       30/1    1.45070    1.43261 +/- 0.01225\n",
      "       31/1    1.48473    1.43509 +/- 0.01191\n",
      "       32/1    1.47264    1.43680 +/- 0.01149\n",
      "       33/1    1.45237    1.43747 +/- 0.01100\n",
      "       34/1    1.41749    1.43664 +/- 0.01056\n",
      "       35/1    1.47971    1.43836 +/- 0.01027\n",
      "       36/1    1.45378    1.43896 +/- 0.00989\n",
      "       37/1    1.35765    1.43594 +/- 0.00998\n",
      "       38/1    1.40397    1.43480 +/- 0.00969\n",
      "       39/1    1.41697    1.43419 +/- 0.00937\n",
      "       40/1    1.40400    1.43318 +/- 0.00910\n",
      "       41/1    1.43270    1.43317 +/- 0.00881\n",
      "       42/1    1.36535    1.43105 +/- 0.00879\n",
      "       43/1    1.43240    1.43109 +/- 0.00851\n",
      "       44/1    1.46349    1.43204 +/- 0.00832\n",
      "       45/1    1.32866    1.42909 +/- 0.00860\n",
      "       46/1    1.32249    1.42613 +/- 0.00886\n",
      "       47/1    1.41734    1.42589 +/- 0.00862\n",
      "       48/1    1.39802    1.42516 +/- 0.00843\n",
      "       49/1    1.54696    1.42828 +/- 0.00878\n",
      "       50/1    1.44601    1.42872 +/- 0.00857\n",
      "       51/1    1.37546    1.42742 +/- 0.00846\n",
      "       52/1    1.48833    1.42887 +/- 0.00838\n",
      "       53/1    1.36864    1.42747 +/- 0.00830\n",
      "       54/1    1.45539    1.42811 +/- 0.00814\n",
      "       55/1    1.47463    1.42914 +/- 0.00802\n",
      "       56/1    1.43531    1.42927 +/- 0.00785\n",
      "       57/1    1.37650    1.42815 +/- 0.00776\n",
      "       58/1    1.43163    1.42822 +/- 0.00760\n",
      "       59/1    1.39161    1.42748 +/- 0.00748\n",
      "       60/1    1.48475    1.42862 +/- 0.00742\n",
      "       61/1    1.47918    1.42961 +/- 0.00734\n",
      "       62/1    1.47997    1.43058 +/- 0.00726\n",
      "       63/1    1.41811    1.43035 +/- 0.00712\n",
      "       64/1    1.41011    1.42997 +/- 0.00700\n",
      "       65/1    1.44239    1.43020 +/- 0.00688\n",
      "       66/1    1.42912    1.43018 +/- 0.00675\n",
      "       67/1    1.39312    1.42953 +/- 0.00666\n",
      "       68/1    1.47339    1.43028 +/- 0.00659\n",
      "       69/1    1.38559    1.42953 +/- 0.00652\n",
      "       70/1    1.42211    1.42940 +/- 0.00641\n",
      "       71/1    1.35744    1.42822 +/- 0.00642\n",
      "       72/1    1.40876    1.42791 +/- 0.00632\n",
      "       73/1    1.53582    1.42962 +/- 0.00645\n",
      "       74/1    1.44605    1.42988 +/- 0.00636\n",
      "       75/1    1.47428    1.43056 +/- 0.00629\n",
      "       76/1    1.37855    1.42977 +/- 0.00625\n",
      "       77/1    1.39439    1.42925 +/- 0.00618\n",
      "       78/1    1.49117    1.43016 +/- 0.00615\n",
      "       79/1    1.34815    1.42897 +/- 0.00618\n",
      "       80/1    1.30610    1.42721 +/- 0.00634\n",
      "       81/1    1.44634    1.42748 +/- 0.00625\n",
      "       82/1    1.41633    1.42733 +/- 0.00617\n",
      "       83/1    1.45974    1.42777 +/- 0.00610\n",
      "       84/1    1.45538    1.42814 +/- 0.00603\n",
      "       85/1    1.45563    1.42851 +/- 0.00596\n",
      "       86/1    1.37790    1.42785 +/- 0.00592\n",
      "       87/1    1.43465    1.42793 +/- 0.00584\n",
      "       88/1    1.32547    1.42662 +/- 0.00591\n",
      "       89/1    1.37892    1.42602 +/- 0.00587\n",
      "       90/1    1.47858    1.42667 +/- 0.00583\n",
      "       91/1    1.34703    1.42569 +/- 0.00584\n",
      "       92/1    1.47197    1.42625 +/- 0.00580\n",
      "       93/1    1.48311    1.42694 +/- 0.00577\n",
      "       94/1    1.39298    1.42653 +/- 0.00571\n",
      "       95/1    1.52197    1.42766 +/- 0.00576\n",
      "       96/1    1.35930    1.42686 +/- 0.00574\n",
      "       97/1    1.51529    1.42788 +/- 0.00577\n",
      "       98/1    1.45297    1.42816 +/- 0.00571\n",
      "       99/1    1.40311    1.42788 +/- 0.00565\n",
      "      100/1    1.42293    1.42783 +/- 0.00559\n",
      " Creating state point statepoint.100.h5...\n",
      "\n",
      " =======================>     TIMING STATISTICS     <=======================\n",
      "\n",
      " Total time for initialization     = 2.0249e+00 seconds\n",
      "   Reading cross sections          = 1.9682e+00 seconds\n",
      " Total time in simulation          = 1.4976e+01 seconds\n",
      "   Time in transport only          = 1.4925e+01 seconds\n",
      "   Time in inactive batches        = 1.8773e-01 seconds\n",
      "   Time in active batches          = 1.4788e+01 seconds\n",
      "   Time synchronizing fission bank = 6.1776e-03 seconds\n",
      "     Sampling source sites         = 4.7949e-03 seconds\n",
      "     SEND/RECV source sites        = 5.6656e-04 seconds\n",
      "   Time accumulating tallies       = 2.4220e-02 seconds\n",
      "   Time writing statepoints        = 1.5840e-02 seconds\n",
      " Total time for finalization       = 2.7855e-02 seconds\n",
      " Total time elapsed                = 1.7039e+01 seconds\n",
      " Calculation Rate (inactive)       = 53268.4 particles/second\n",
      " Calculation Rate (active)         = 6086.03 particles/second\n",
      "\n",
      " ============================>     RESULTS     <============================\n",
      "\n",
      " k-effective (Collision)     = 1.42881 +/- 0.00466\n",
      " k-effective (Track-length)  = 1.42783 +/- 0.00559\n",
      " k-effective (Absorption)    = 1.42667 +/- 0.00316\n",
      " Combined k-effective        = 1.42696 +/- 0.00299\n",
      " Leakage Fraction            = 0.00163 +/- 0.00017\n",
      "\n"
     ]
    }
   ],
   "source": [
    "statepoint = model.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [],
   "source": [
    "with openmc.StatePoint(statepoint) as sp:\n",
    "    mesh_tally_out = sp.get_tally(id=heat.id)\n",
    "    \n",
    "mesh_flux = mesh_tally_out.get_values()\n",
    "mesh_flux = mesh_flux.reshape(mesh.dimension)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x12c530a10>"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 28512 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "img = plt.imshow(mesh_flux)\n",
    "plt.colorbar(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x12c9ad8d0>"
      ]
     },
     "execution_count": 115,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 28344 characters]",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mesh_flux = mesh_tally_out.get_values(value='std_dev')\n",
    "mesh_flux = mesh_flux.reshape(mesh.dimension)\n",
    "plt.figure()\n",
    "img = plt.imshow(mesh_flux)\n",
    "plt.colorbar(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Universes\n",
    "\n",
    "A universe is a collection of cells that can be used as a repeatable unit in the geometry. At a minimum, there must be one \"root\" universe (say, named `root`), which gets passed to `openmc.Geometry(root)`. But you can also use universes to repeat a collection of cells multiple times throughout a geometry. Here, we will explore some basic features of universes.\n",
    "\n",
    "We'll start by making a universe which looks similar to the pincell we built earlier - say, a cylinder of UO2 enclosed by an infinite region of water. First, we set up our materials and create our geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model.materials)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pin_surface = openmc.ZCylinder(r=1.0)\n",
    "inside_pin = -pin_surface\n",
    "\n",
    "pin_cell = openmc.Cell(region=inside_pin, fill=uo2)\n",
    "outside_cell = openmc.Cell(region=~inside_pin, fill=water)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "universe = openmc.Universe()\n",
    "universe.add_cells([pin_cell, outside_cell])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have created a universe containing a pin, enclosed in an infinite medium of water. Now let's suppose that I want to fill this universe into an enclosing cell, a cylinder of radius 5 cm. Let's first create this cylinder, and then we will fill it with our `universe`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "big_cylinder = openmc.ZCylinder(r=5.0)\n",
    "big_cell = openmc.Cell()\n",
    "big_cell.region = -big_cylinder\n",
    "big_cell.fill = universe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at our geometry. In order to visualize at this stage, we need to create a universe from our `big_cell`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "big_universe = openmc.Universe(cells=[big_cell])\n",
    "big_universe.plot(width=(10.0, 10.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that our `big_cell`, the large cylinder, has been filled with the `universe` we declared earlier. Let's increase the complexity a bit to understand how this filling works. What if the cylinder of UO2 in our `universe` is not located at the origin, but is instead shifted to a different position?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pin_surface.x0 = 1.0\n",
    "pin_surface.y0 = 1.5\n",
    "big_universe.plot(width=(10.0, 10.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that when we fill a universe inside of another cell, that there's (by default) no transformation of coordinates. You can shift the position of the universe filling a cell with the `Cell.translation` attribute. There are similar adjustments you can make, like rotations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "big_cell.translation = [-1.0, 0.0, 0.0]\n",
    "big_universe.plot(width=(10.0, 10.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lattices\n",
    "\n",
    "Lattices are a convenient way to (i) repeat a universe multiple times in space, while (ii) automatically translating that universe's origin to different positions in space. In order to explore this concept, let's work with a PWR fuel assembly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = {\n",
    "    water: 'lightblue',\n",
    "    zirconium: 'gray',\n",
    "    uo2: 'red'\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will make individual universes for each repeatable unit in the assembly design."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fuel Pin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pitch = 1.26\n",
    "fuel_or = openmc.ZCylinder(r=0.39)\n",
    "clad_ir = openmc.ZCylinder(r=0.40)\n",
    "clad_or = openmc.ZCylinder(r=0.46)\n",
    "\n",
    "fuel = openmc.Cell(fill=uo2, region=-fuel_or)\n",
    "gap = openmc.Cell(region=+fuel_or & -clad_ir)\n",
    "clad = openmc.Cell(fill=zirconium, region=+clad_ir & -clad_or)\n",
    "moderator = openmc.Cell(fill=water, region=+clad_or)\n",
    "\n",
    "fuel_pin = openmc.Universe(cells=[fuel, gap, clad, moderator])\n",
    "fuel_pin.plot(width=(pitch, pitch), color_by='material', colors=colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Guide Tube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "guide_clad_ir = openmc.ZCylinder(r=0.56)\n",
    "guide_clad_or = openmc.ZCylinder(r=0.60)\n",
    "\n",
    "guide_inner = openmc.Cell(fill=water, region=-guide_clad_ir)\n",
    "guide_clad = openmc.Cell(fill=zirconium, region=+guide_clad_ir & -guide_clad_or)\n",
    "guide_outer = openmc.Cell(fill=water, region=+guide_clad_or)\n",
    "\n",
    "guide_tube = openmc.Universe(cells=[guide_inner, guide_clad, guide_outer])\n",
    "guide_tube.plot(width=(pitch, pitch), color_by='material', colors=colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to build our fuel assembly, by stamping these two universes into a repeatable pattern. OpenMC has `RectLattice` and `HexLattice` objects, to create lattices on a Cartesian grid or on a hexagonal grid. For our fuel assembly, we need to use `RectLattice`. \n",
    "\n",
    "When creating a rectangular lattice, we need to define:\n",
    "\n",
    "1. The lower-left coordinates of the lattice (`.lower_left`)\n",
    "2. The size of each lattice element (`.pitch`)\n",
    "3. The 2D arrangement of universes (`.universes`)\n",
    "4. (_optionally_) A universe that is used outside of the defined region (`.outer`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lattice = openmc.RectLattice()\n",
    "\n",
    "# Define the lower-left coordinates and size of lattice elements\n",
    "lattice.lower_left = (-pitch, -pitch)\n",
    "lattice.pitch = (pitch, pitch)\n",
    "\n",
    "# Now we specify what is actually inside the lattice.  This 2D lattice will be a\n",
    "# list of lists like\n",
    "# lattice.universes = [ [a1, a2, ...], [b1, b2, ...], ...]\n",
    "# The inner lists specify columns from left to right.  The outer lists specify\n",
    "# rows from top to bottom.\n",
    "lattice.universes = [\n",
    "    [fuel_pin,   fuel_pin],\n",
    "    [guide_tube, fuel_pin]\n",
    "]\n",
    "\n",
    "# We also should specify what is outside of the lattice.  In this case, it is\n",
    "# the infinite water universe.\n",
    "all_water = openmc.Universe(cells=[openmc.Cell(fill=water)])\n",
    "lattice.outer = all_water"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To plot the lattice, we need to put it in a universe. For this, we'll create a single cell filled with the lattice, and then put that single cell inside a universe to plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = openmc.Universe(cells=[openmc.Cell(fill=lattice)])\n",
    "u.plot(width=(3*pitch, 3*pitch), color_by='material', colors=colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What exactly does `outer` mean?\n",
    "\n",
    "In the previous section, we set the lattice outer universe to a universe containing a single cell with only water in it. To get a better sense of what the outer universe does, let's change the outer universe to the guide tube universe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lattice.outer = guide_tube\n",
    "u.plot(width=(3*pitch, 3*pitch), color_by='material', colors=colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 17x17 Fuel Assembly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "guide_tube_positions = [\n",
    "    (2, 5), (2, 8), (2, 11),\n",
    "    (5, 2), (5, 5), (5, 8), (5, 11), (5, 14),\n",
    "    (8, 2), (8, 5), (8, 8), (8, 11), (8, 14),\n",
    "    (11, 2), (11, 5), (11, 8), (11, 11), (11, 14),\n",
    "    (14, 5), (14, 8), (14, 11)\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lattice = openmc.RectLattice()\n",
    "\n",
    "lattice.pitch = (pitch, pitch)\n",
    "lattice.outer = all_water\n",
    "\n",
    "# I want (x0, y0) = (0, 0) to be the center of the instrument tube so that means\n",
    "# the lower-left will be -half a pin pitch in x and y.\n",
    "assembly_pitch = 17*pitch\n",
    "lattice.lower_left = (-assembly_pitch/2, -assembly_pitch/2)\n",
    "\n",
    "# Most of the lattice positions are fuel pins so rather than type all of those\n",
    "# out, we use a list comprehension to start with a 17x17 array of fuel.\n",
    "lattice.universes = [[fuel_pin] * 17] * 17\n",
    "\n",
    "for row, col in guide_tube_positions:\n",
    "    lattice.universes[row, col] = guide_tube"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we just have to add the boundary conditions and root universe to finish the geometry. To create a box containing the lattice, we'll use the `rectangular_prism` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "box = openmc.model.rectangular_prism(assembly_pitch, assembly_pitch, boundary_type='reflective')\n",
    "main_cell = openmc.Cell(fill=lattice, region=box)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = openmc.Universe(cells=(main_cell,))\n",
    "u.plot(width=(assembly_pitch, assembly_pitch), color_by='material', colors=colors, pixels=(500,500))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TRISO Particles\n",
    "\n",
    "OpenMC includes a few convenience functions for generating locations of randomly packed spheres that can be used to model TRISO particles and/or pebbles in a reactor core. To be clear, this capability is not a stochastic geometry capability like that included in MCNP. It's also important to note that OpenMC does not use delta tracking, which would normally speed up calculations in geometries with tons of surfaces and cells. However, the computational burden can be eased by placing random spheres in a lattice.\n",
    "\n",
    "This capability relies on three functions/classes:\n",
    "- `openmc.model.pack_spheres` -- generate locations of random spheres\n",
    "- `openmc.model.TRISO` -- Cell-like object that holds a universe storing the internal structure of a pebble/TRISO\n",
    "- `openmc.model.create_triso_lattice` -- Creates a lattice containing `TRISO` objects for improved tracking performance\n",
    "\n",
    "Let's start with the `pack_spheres` function. This function takes an outer radius of the spheres, a containing region, and a packing fraction and will return an array of sphere coordinates. For our example, let's use spheres with a radius of 1 cm and a packing fraction of 30%. We'll put our spheres inside of a finite cylinder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#openmc.model.pack_spheres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_cylinder = 10.\n",
    "h_cylinder = 20.\n",
    "cylinder = openmc.model.RightCircularCylinder(\n",
    "    center_base=(0., 0., -10.), \n",
    "    height=h_cylinder,\n",
    "    radius=r_cylinder\n",
    ")\n",
    "\n",
    "r_sphere = 1.0\n",
    "centers = openmc.model.pack_spheres(radius=r_sphere, region=-cylinder, pf=0.3)\n",
    "centers.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "centers[:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need to actually creates cells for each of these spheres. To do so, we'll use the `TRISO` class. We'll need to define a universe that we want to fill each sphere. We'll create a universe with a single, infinite cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell = openmc.Cell(fill=uo2)\n",
    "sphere_univ = openmc.Universe(cells=[cell])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can create a `TRISO` object for each sphere center, in order to create one cell for each TRISO region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trisos = [openmc.model.TRISO(outer_radius=r_sphere, fill=sphere_univ, center=center) for center in centers]\n",
    "trisos[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's confirm that the packing fraction of our TRISOs is actually about 30%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "volume_trisos = len(trisos)*4/3*math.pi*r_sphere**3\n",
    "volume_cyl = math.pi * r_cylinder**2 * h_cylinder\n",
    "volume_trisos / volume_cyl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_background = openmc.Intersection([~t.region for t in trisos])\n",
    "background_cell = openmc.Cell(region=region_background & -cylinder)\n",
    "univ = openmc.Universe(cells=[background_cell] + trisos)\n",
    "univ.plot(width=(20., 20.))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While this works in principle, it will lead to **very** poor tracking performance; every time a particle reaches the background cell, it has to determine the distance to the boundary of _every single_ sphere. To improve tracking performance, we can use the `create_triso_lattice` function to overlay a lattice that limits how many distance checks need to be performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lower_left, upper_right = (-cylinder).bounding_box\n",
    "shape = (5, 5, 5)\n",
    "pitch = (upper_right - lower_left)/shape\n",
    "openmc.model.create_triso_lattice\n",
    "lattice = openmc.model.create_triso_lattice(\n",
    "    trisos=trisos,\n",
    "    lower_left=lower_left,\n",
    "    pitch=pitch,\n",
    "    shape=shape,\n",
    "    background=zirconium\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = openmc.Cell(fill=lattice, region=-cylinder)\n",
    "univ = openmc.Universe(cells=[c])\n",
    "univ.plot(width=(20., 20.))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "univ.plot(width=(20., 20.), color_by='material')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
