{
 "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": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With UO2 finished, let's now create materials for the coolant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see if `inside_sphere` actually contains points inside the sphere:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run OpenMC in plotting mode by calling the `plot_geometry()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can simply re-run our model with this additional tally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": []
  },
  {
   "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": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": []
  },
  {
   "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": []
  },
  {
   "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": []
  },
  {
   "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": []
  },
  {
   "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": []
  },
  {
   "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": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Guide Tube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": []
  },
  {
   "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": []
  },
  {
   "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": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 17x17 Fuel Assembly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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_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": []
  },
  {
   "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": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}
