mg-mode-part-i.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multigroup Mode Part I: Introduction\n",
    "This Notebook illustrates the usage of OpenMC's multi-group calculational mode with the Python API. This example notebook creates and executes the 2-D [C5G7](https://www.oecd-nea.org/jcms/pl_17882) benchmark model using the `openmc.MGXSLibrary` class to create the supporting data library on the fly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate MGXS Library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "import numpy as np\n",
    "\n",
    "import openmc\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a model to tie together geometry, materials, settings, and tallies\n",
    "model = openmc.Model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now create the multi-group library using data directly from Appendix A of the [C5G7](https://www.oecd-nea.org/jcms/pl_17882) benchmark documentation. All of the data below will be created at 294K, consistent with the benchmark.\n",
    "\n",
    "This notebook will first begin by setting the group structure and building the groupwise data for UO2.  As you can see, the cross sections are input in the order of increasing groups (or decreasing energy).\n",
    "\n",
    "*Note*: The C5G7 benchmark uses transport-corrected cross sections. So the total cross section we input here will technically be the transport cross section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a 7-group structure with arbitrary boundaries (the specific boundaries are unimportant)\n",
    "groups = openmc.mgxs.EnergyGroups(np.logspace(-5, 7, 8))\n",
    "\n",
    "uo2_xsdata = openmc.XSdata('uo2', groups)\n",
    "uo2_xsdata.order = 0\n",
    "\n",
    "# When setting the data let the object know you are setting the data for a temperature of 294K.\n",
    "uo2_xsdata.set_total([1.77949E-1, 3.29805E-1, 4.80388E-1, 5.54367E-1,\n",
    "                      3.11801E-1, 3.95168E-1, 5.64406E-1], temperature=294.)\n",
    "\n",
    "uo2_xsdata.set_absorption([8.0248E-03, 3.7174E-3, 2.6769E-2, 9.6236E-2,\n",
    "                           3.0020E-02, 1.1126E-1, 2.8278E-1], temperature=294.)\n",
    "uo2_xsdata.set_fission([7.21206E-3, 8.19301E-4, 6.45320E-3, 1.85648E-2,\n",
    "                        1.78084E-2, 8.30348E-2, 2.16004E-1], temperature=294.)\n",
    "\n",
    "uo2_xsdata.set_nu_fission([2.005998E-2, 2.027303E-3, 1.570599E-2, 4.518301E-2,\n",
    "                           4.334208E-2, 2.020901E-1, 5.257105E-1], temperature=294.)\n",
    "\n",
    "uo2_xsdata.set_chi([5.87910E-1, 4.11760E-1, 3.39060E-4, 1.17610E-7,\n",
    "                    0.00000E-0, 0.00000E-0, 0.00000E-0], temperature=294.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now add the scattering matrix data. \n",
    "\n",
    "*Note*: Most users familiar with deterministic transport libraries are already familiar with the idea of entering one scattering matrix for every order (i.e. scattering order as the outer dimension). However, the shape of OpenMC's scattering matrix entry is instead [Incoming groups, Outgoing Groups, Scattering Order] to best enable other scattering representations. We will follow the more familiar approach in this notebook, and then use numpy's `numpy.rollaxis` function to change the ordering to what we need (scattering order on the inner dimension)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The scattering matrix is ordered with incoming groups as rows and outgoing groups as columns\n",
    "# (i.e., below the diagonal is up-scattering).\n",
    "scatter_matrix = \\\n",
    "    [[[1.27537E-1, 4.23780E-2, 9.43740E-6, 5.51630E-9, 0.00000E-0, 0.00000E-0, 0.00000E-0],\n",
    "      [0.00000E-0, 3.24456E-1, 1.63140E-3, 3.14270E-9, 0.00000E-0, 0.00000E-0, 0.00000E-0],\n",
    "      [0.00000E-0, 0.00000E-0, 4.50940E-1, 2.67920E-3, 0.00000E-0, 0.00000E-0, 0.00000E-0],\n",
    "      [0.00000E-0, 0.00000E-0, 0.00000E-0, 4.52565E-1, 5.56640E-3, 0.00000E-0, 0.00000E-0],\n",
    "      [0.00000E-0, 0.00000E-0, 0.00000E-0, 1.25250E-4, 2.71401E-1, 1.02550E-2, 1.00210E-8],\n",
    "      [0.00000E-0, 0.00000E-0, 0.00000E-0, 0.00000E-0, 1.29680E-3, 2.65802E-1, 1.68090E-2],\n",
    "      [0.00000E-0, 0.00000E-0, 0.00000E-0, 0.00000E-0, 0.00000E-0, 8.54580E-3, 2.73080E-1]]]\n",
    "scatter_matrix = np.array(scatter_matrix)\n",
    "scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)\n",
    "uo2_xsdata.set_scatter_matrix(scatter_matrix, temperature=294.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the UO2 data has been created, we can move on to the remaining materials using the same process.\n",
    "\n",
    "However, we will actually skip repeating the above for now.  Our simulation will instead use the `c5g7.h5` file that has already been created using exactly the same logic as above, but for the remaining materials in the benchmark problem.\n",
    "\n",
    "For now we will show how you would use the `uo2_xsdata` information to create an `openmc.MGXSLibrary` object and write to disk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the library\n",
    "mg_cross_sections_file = openmc.MGXSLibrary(groups)\n",
    "\n",
    "# Add the UO2 data to it\n",
    "mg_cross_sections_file.add_xsdata(uo2_xsdata)\n",
    "\n",
    "# And write to disk\n",
    "mg_cross_sections_file.export_to_hdf5('mgxs.h5')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate 2-D C5G7 Problem Input Files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To build the actual 2-D model, we will first begin by creating the `materials.xml` file.\n",
    "\n",
    "First we need to define materials that will be used in the problem. In other notebooks, either nuclides or elements were added to materials at the equivalent stage. We can do that in multi-group mode as well.  However, multi-group cross-sections are sometimes provided as macroscopic cross-sections; the C5G7 benchmark data are macroscopic. In this case, we can instead use the `Material.add_macroscopic` method to specify a macroscopic object. Unlike for nuclides and elements, we do not need provide information on atom/weight percents as no number densities are needed.\n",
    "\n",
    "When assigning macroscopic objects to a material, the density can still be scaled by setting the density to a value that is not 1.0. This would be useful, for example, when slightly perturbing the density of water due to a small change in temperature (while of course ignoring any resultant spectral shift). The density of a macroscopic dataset is set to 1.0 in the `openmc.Material` object by default when a macroscopic dataset is used; so we will show its use the first time and then afterwards it will not be required.\n",
    "\n",
    "Aside from these differences, the following code is very similar to similar code in other OpenMC example Notebooks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For every cross section data set in the library, assign an openmc.Macroscopic object to a material\n",
    "material_dict = {}\n",
    "for xs in ['uo2', 'mox43', 'mox7', 'mox87', 'fiss_chamber', 'guide_tube', 'water']:\n",
    "    material_dict[xs] = openmc.Material(name=xs)\n",
    "    material_dict[xs].set_density('macro', 1.)\n",
    "    material_dict[xs].add_macroscopic(xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can go ahead and produce a `materials.xml` file for use by OpenMC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate a Materials collection, register all Materials, and export to XML\n",
    "materials = openmc.Materials(material_dict.values())\n",
    "\n",
    "# Set the location of the cross sections file to our pre-written set\n",
    "materials.cross_sections = 'c5g7.h5'\n",
    "\n",
    "model.materials = materials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our next step will be to create the geometry information needed for our assembly and to write that to the `geometry.xml` file.\n",
    "\n",
    "We will begin by defining the surfaces, cells, and universes needed for each of the individual fuel pins, guide tubes, and fission chambers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the surface used for each pin\n",
    "pin_surf = openmc.ZCylinder(x0=0, y0=0, r=0.54, name='pin_surf')\n",
    "\n",
    "# Create the cells which will be used to represent each pin type.\n",
    "cells = {}\n",
    "universes = {}\n",
    "for material in material_dict.values():\n",
    "    # Create the cell for the material inside the cladding\n",
    "    cells[material.name] = openmc.Cell(name=material.name)\n",
    "    # Assign the half-spaces to the cell\n",
    "    cells[material.name].region = -pin_surf\n",
    "    # Register the material with this cell\n",
    "    cells[material.name].fill = material\n",
    "    \n",
    "    # Repeat the above for the material outside the cladding (i.e., the moderator)\n",
    "    cell_name = material.name + '_moderator'\n",
    "    cells[cell_name] = openmc.Cell(name=cell_name)\n",
    "    cells[cell_name].region = +pin_surf\n",
    "    cells[cell_name].fill = material_dict['water']\n",
    "    \n",
    "    # Finally add the two cells we just made to a Universe object\n",
    "    universes[material.name] = openmc.Universe(name=material.name)\n",
    "    universes[material.name].add_cells([cells[material.name], cells[cell_name]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to take our universes (representing the different pin types) and lay them out in a lattice to represent the assembly types"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "lattices = {}\n",
    "\n",
    "# Instantiate the UO2 Lattice\n",
    "lattices['UO2 Assembly'] = openmc.RectLattice(name='UO2 Assembly')\n",
    "lattices['UO2 Assembly'].dimension = [17, 17]\n",
    "lattices['UO2 Assembly'].lower_left = [-10.71, -10.71]\n",
    "lattices['UO2 Assembly'].pitch = [1.26, 1.26]\n",
    "u = universes['uo2']\n",
    "g = universes['guide_tube']\n",
    "f = universes['fiss_chamber']\n",
    "lattices['UO2 Assembly'].universes = \\\n",
    "    [[u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],\n",
    "     [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],\n",
    "     [u, u, u, u, u, g, u, u, g, u, u, g, u, u, u, u, u],\n",
    "     [u, u, u, g, u, u, u, u, u, u, u, u, u, g, u, u, u],\n",
    "     [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],\n",
    "     [u, u, g, u, u, g, u, u, g, u, u, g, u, u, g, u, u],\n",
    "     [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],\n",
    "     [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],\n",
    "     [u, u, g, u, u, g, u, u, f, u, u, g, u, u, g, u, u],\n",
    "     [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],\n",
    "     [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],\n",
    "     [u, u, g, u, u, g, u, u, g, u, u, g, u, u, g, u, u],\n",
    "     [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],\n",
    "     [u, u, u, g, u, u, u, u, u, u, u, u, u, g, u, u, u],\n",
    "     [u, u, u, u, u, g, u, u, g, u, u, g, u, u, u, u, u],\n",
    "     [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u],\n",
    "     [u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u, u]]\n",
    "    \n",
    "# Create a containing cell and universe\n",
    "cells['UO2 Assembly'] = openmc.Cell(name='UO2 Assembly')\n",
    "cells['UO2 Assembly'].fill = lattices['UO2 Assembly']\n",
    "universes['UO2 Assembly'] = openmc.Universe(name='UO2 Assembly')\n",
    "universes['UO2 Assembly'].add_cell(cells['UO2 Assembly'])\n",
    "\n",
    "# Instantiate the MOX Lattice\n",
    "lattices['MOX Assembly'] = openmc.RectLattice(name='MOX Assembly')\n",
    "lattices['MOX Assembly'].dimension = [17, 17]\n",
    "lattices['MOX Assembly'].lower_left = [-10.71, -10.71]\n",
    "lattices['MOX Assembly'].pitch = [1.26, 1.26]\n",
    "m = universes['mox43']\n",
    "n = universes['mox7']\n",
    "o = universes['mox87']\n",
    "g = universes['guide_tube']\n",
    "f = universes['fiss_chamber']\n",
    "lattices['MOX Assembly'].universes = \\\n",
    "    [[m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m],\n",
    "     [m, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, m],\n",
    "     [m, n, n, n, n, g, n, n, g, n, n, g, n, n, n, n, m],\n",
    "     [m, n, n, g, n, o, o, o, o, o, o, o, n, g, n, n, m],\n",
    "     [m, n, n, n, o, o, o, o, o, o, o, o, o, n, n, n, m],\n",
    "     [m, n, g, o, o, g, o, o, g, o, o, g, o, o, g, n, m],\n",
    "     [m, n, n, o, o, o, o, o, o, o, o, o, o, o, n, n, m],\n",
    "     [m, n, n, o, o, o, o, o, o, o, o, o, o, o, n, n, m],\n",
    "     [m, n, g, o, o, g, o, o, f, o, o, g, o, o, g, n, m],\n",
    "     [m, n, n, o, o, o, o, o, o, o, o, o, o, o, n, n, m],\n",
    "     [m, n, n, o, o, o, o, o, o, o, o, o, o, o, n, n, m],\n",
    "     [m, n, g, o, o, g, o, o, g, o, o, g, o, o, g, n, m],\n",
    "     [m, n, n, n, o, o, o, o, o, o, o, o, o, n, n, n, m],\n",
    "     [m, n, n, g, n, o, o, o, o, o, o, o, n, g, n, n, m],\n",
    "     [m, n, n, n, n, g, n, n, g, n, n, g, n, n, n, n, m],\n",
    "     [m, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, m],\n",
    "     [m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m, m]]\n",
    "    \n",
    "# Create a containing cell and universe\n",
    "cells['MOX Assembly'] = openmc.Cell(name='MOX Assembly')\n",
    "cells['MOX Assembly'].fill = lattices['MOX Assembly']\n",
    "universes['MOX Assembly'] = openmc.Universe(name='MOX Assembly')\n",
    "universes['MOX Assembly'].add_cell(cells['MOX Assembly'])\n",
    "    \n",
    "# Instantiate the reflector Lattice\n",
    "lattices['Reflector Assembly'] = openmc.RectLattice(name='Reflector Assembly')\n",
    "lattices['Reflector Assembly'].dimension = [1,1]\n",
    "lattices['Reflector Assembly'].lower_left = [-10.71, -10.71]\n",
    "lattices['Reflector Assembly'].pitch = [21.42, 21.42]\n",
    "lattices['Reflector Assembly'].universes = [[universes['water']]]\n",
    "\n",
    "# Create a containing cell and universe\n",
    "cells['Reflector Assembly'] = openmc.Cell(name='Reflector Assembly')\n",
    "cells['Reflector Assembly'].fill = lattices['Reflector Assembly']\n",
    "universes['Reflector Assembly'] = openmc.Universe(name='Reflector Assembly')\n",
    "universes['Reflector Assembly'].add_cell(cells['Reflector Assembly'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now create the core layout in a 3x3 lattice where each lattice position is one of the assemblies we just defined.\n",
    "\n",
    "After that we can create the final cell to contain the entire core."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "lattices['Core'] = openmc.RectLattice(name='3x3 core lattice')\n",
    "lattices['Core'].dimension= [3, 3]\n",
    "lattices['Core'].lower_left = [-32.13, -32.13]\n",
    "lattices['Core'].pitch = [21.42, 21.42]\n",
    "r = universes['Reflector Assembly']\n",
    "u = universes['UO2 Assembly']\n",
    "m = universes['MOX Assembly']\n",
    "lattices['Core'].universes = [[u, m, r],\n",
    "                              [m, u, r],\n",
    "                              [r, r, r]]\n",
    "\n",
    "# Create boundary planes to surround the geometry\n",
    "min_x = openmc.XPlane(x0=-32.13, boundary_type='reflective')\n",
    "max_x = openmc.XPlane(x0=+32.13, boundary_type='vacuum')\n",
    "min_y = openmc.YPlane(y0=-32.13, boundary_type='vacuum')\n",
    "max_y = openmc.YPlane(y0=+32.13, boundary_type='reflective')\n",
    "\n",
    "# Create root Cell\n",
    "root_cell = openmc.Cell(name='root cell')\n",
    "root_cell.fill = lattices['Core']\n",
    "\n",
    "# Add boundary planes\n",
    "root_cell.region = +min_x & -max_x & +min_y & -max_y\n",
    "\n",
    "# Create root Universe\n",
    "root_universe = openmc.Universe(name='root universe', universe_id=0)\n",
    "root_universe.add_cell(root_cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we commit to the geometry, we should view it using the Python API's plotting capability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f22b8fc6e50>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh4AAAITCAYAAABIXRWUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABEuklEQVR4nO3df6w113Xe92eXopUgFmEKFmSJZGoZZmkygqJEBCGjRZTaLkwHRui0jSmDrZXEAGFARmMghStFpJ3UUuBAQNOmidMSiREFICwKVhITrl3rB4yo/UNmmIRxJNFvyMQ1JEGUmpiyZbhQSHn3j3tGGs47M2fNzJq91+zz/QCCXs5Zd/ZeZ+ae2ffcM89NOWcBAACU8B/VngAAALgcLDwAAEAxLDwAAEAxLDwAAEAxLDwAAEAxLDwAAEAxL6s9gb5X3HxzftUtt9SeBgAA2OA3PvnJf5dzftXYY6EWHq+65Rb9tQ/+nG567uHrHvudb/rJ67ZFros8N2td5LlZ67zHfOpNN1637Y3/7IXidTXG9K6LPDdrXaS5AZH8wLfd+ZtTj4X7VcvYBWBse+S6yHOz1kWe25I6z32NXQDGtpeoqzGmd13kuVnroswNOJJwC485lovKkrq991tyrCP30NXttd+99tHXXQisFwSPur3GoodldbV7AI7mUAuPpc5dNDwuwnvz6qFmr15j1+zB6yJQ82JCD/77qTk2CxMcVdMLD3zN2OcVAAAoremFx7mL7REuxvSwvG4P5z7MZ/2wX80PBVp78Op1D5fQg7WOD5jiqMItPOYuLntdoDz223+MHpbVee63dg8eF5X+Y9Y6jzkt5d2DR69LXdLxAqJIOefac/iqb3n96/Nf++DPSbr+9/lTFwBL3dpbKreM6V23962o9GCrq3F7pncdPbTZAxDJD3zbnf8s53z32GNhFx4AAOCY5hYe4X7VAgAA2hUqubQzddtkhGTNLXOz1tFD3B6sb3s/8Myd12179Panz9bUqhvWWOvoYf8e+FULWhPuHY+tuRQ1kjCHda33ECHhdMvcrHVrk0unLjxT2+fq5r7GWrdmbvQQoweSS9GizQuPlNIfSCk9kVL6lymlT6aU/upp++tSSr+aUno2pfRYSunrto7lnYS5V2Kmxz5qp356uKQeugvB1gvPHpbOiR72sXROLC7QKo93PL4s6Ttyzn9U0hsl3ZtSerOkvy7pb+Scv1XS85J+yGGsKkj9jDF2C8mlJS90e41NDzHGZmGCo9q88MhXfvf0nzee/pclfYek7haV90n6vq1j1XKEkK5z6CGGsd/h47JwDuDSuXzGI6V0Q0rpKUlfkPRhSf9G0hdzzi+eSj4j6ZaJr30wpfRkSunJLz3/vMd0vsqamBn5gkYPMVhTJCNfVM7NrXvcWldDCz2c082N5FK0ymXhkXP+Ss75jZJulXSPpG9b8LWP5Jzvzjnf/Yqbby6aSmmxZkx6OP/YUi300Dd34es/5l3niR7W1wGXzD1ALKX045L+P0n/g6Rvyjm/mFL6dkl/Jef83XNfu1dy6bBu7mJiqYuQcEoP9Xro/2597qfO/u/ypy46w9/316ibuyDSQ7m6uR6s5xwQxa7JpSmlV0l6Ief8xZTSH5T0IV19sPRtkj6Yc35/Sul/k/RrOeefntsXyaUAABzf3MLDI0DsNZLel1K6QVe/uvlAzvkXUkqfkvT+lNK7Jf0LSX9v6Y67n0C3/mS8Z925t+Kv3f+QJOmOx95dbG7dmNZxz/VgqVs65rn9eR8HzzrPdzz6defeivessx4v7x4s3w/W/VnntrTXkseBdzxwiUL+rZYaqZTeY/Zf7PqGL3yePZQY03vcGomkW+qmbmHsXwwiJWGO1VmOV40xvcetMaZ3neV8AyI61N9q2Zo2ad3fnkmYUy94w8dK9WAd07vXqa+xbPeYm/e5tDU3wTsJc01ipvV4Wca01lm/H0qFflmfgwjHawoZHjiycAsPD3sll5a0tIetF5Q9dHNq4Tice6Ffmlxa0tJFSQupn2sWF1FYk0tZfOComlx4eIt4QVwq4sJkqRaOQ80Lndc50ELqZ83vh4iLHaAkFh4GkYOtrM59mO8IWjgONbMcvM6ByAFiVjW/H8jzwKVrcuFhTcKMnJi5tIeIC4tuTi0chyMnl1rn1B2viD1YdXO3fj9E7JXkUrQu3MLDI23Ssr81X2+d29yLXv+xUj1Yx/Tudepr1qiVXNq39YV+zyTMpYuLtdbMzfr94N2rZcw5EY7XFBYdOLKQt9N2yPFYV0eOx7o665iRczz2yrbwHrfGmN7jkuMBTNs1udQTyaUAABzf3smlu6jxtz8k209Hkf9+ibUu+t9qqX0c5upq/K2W4V0YU8+JdX9PvPCWr/77nhv/yWiNta5fU6tuaw+1jgPveOAShXzHwzNtcu4WzH5drdRPzx4ip362cBzmchP2TC71TOAcXtT7+hflrXXDC7xnnXcPkdNXreccEM1FJ5dutWfqp7Vuq1qpn549tHAc1iRceqR+7pUbMXdh7z/mXeep1nEgywOXLNzCw0MLqZ+dI6d+dlo4DkdOLu2cu3h3j1vramihh3NILkXrmlx4eIt8QbSKvDCxauE4RF6YoAzOAVw6Fh4GEcO5lqoZzuU1ds3j4PX79JqBVXMfwKyxn5pj1+zB6xzgMx44qkMtPJYGbnmnfpa8eO+VvlqjB2td5ONwztLk0pILEOtY3cXYelEuefFeOidrfcTjYE0uBY4q3MLDO83TWjd10VtzMdyamLm2hxoprS0fh6GpC8Fw+9QFZs2ixCP1s2/qgjzcbq2zjDP3Nda6NXPb2sOex8F6jljPOeBIQt5Ou9RNzz1sunhY667d/1Dxt/W9e7DWeap1HDx7te7rqTfdaHrxf+CZO00XI2udhXVfT7zwFtNF2LvOotbcahwHa531nAMiILkUAAAUczHJpVN1Y3d0WOoipX7Sw75zs9aN3cI49lPo2J0LYz/Vrkm4nKuzpH6O3UrqWVciubRWD57Hy3qOWM854ChCvuNRIs3Tu65GIql3HT3M121NLpVeemHxTLgslUhqrTtScql3wqmlznqOkFyKo2oyuXTtPvbar3cSJj2sq6vdgzWjwSPh0jrWXmFZLSSXWvdb8niR84HWhVt4eLImZkZGD8vr9uCVIlnzotJC6mcLPVjPAZJL0aqmFx74miMsUAAA7Wt64WEN34rMq4cWkktr9kByqf9+ao5NcilQz6EWHntdZEtevOlhuu7IPeyVXDpXtzU0bKrOO7nUGiDmOZb3fj2OA8mlwJVwC49aaZ6edZHnZq2LPDdrnfeY3smlWxJOhzXeaZ4l6iLPbWz7muMytZ3kUlyykLfTAgCA4zrU7bQAAKBdIZNL1yZcTtVdu/+h67aN/Q0Qz3FrjFli3CjP75Zxt4zpnVxq7cEzzXPtmN7j1hjTe9wtY5JciksV7h2Pqds+h9utdWMvAmPbPcetMWapcSM8v1vHXTvmVG7CcPtUTsNwu7WHsYvicNtULsVw+5YxvcetMab3uGvHtJ4j1nMOOJJwC485kRMzraZemLaO5b1fj9TPpXPaa79LH1tT110IrOFQ1h480jytdd2c9trv0sfWjNXV7bXfLfuwzqk7h1hcoFWHWngsde4b3fpCUJNX6mfNXr3CyyInl1rVTC71OgciJ5da1fx+8DoHWJjgqJpeeOBrxn7HDABAaU0vPM5dbI9wMSZ9dXndHs59mM/6Yb+ayaVe3w9HSC6N/L1vPQe8zjkgmnALj7mLS+TUz/5jcy9q1he8pT147NcSvnXusTVz8txv/7Ea55L1omLtoWTq51779fh+iNxr/zGPXmsuToESwgaIWW6BtNbVuMV0bNypFx7P2z09bzFdUmfpde/babc8v9a6vW+nnerBchunta7GLaZj427pde/babc8v2PjcjstLs1cgFjYhQcAADgmkksBAEAIh0kula5/69v7rXtL3Za5WevoIW4Pnr9qmbqtclg3dQvp8G1+79RPz7pL6sF6XPlVCy5VuHc8tgY7rU2ltNZZ5tN6D5a6FnvwTi6d0q+zBld5p3561l1SD2uO8dx2kkvRonALjzlLLz7W8K2SyaXeY9XoweqSeliaXFoySGxpOqd3cqmHFnpYem6wuECrDrXwqKVmYmYLqZ8t9NBCcqnXRbaF5NKaPZBcikvHwsPgCCFd59ADACCCphce1sTMyBc0eojBmiJ5LvypZjjUubCs7nFrXQ0t9GA9R0guRavCLTxKJmsunY/33OhhnRZ6mLv49B+zJmZ6JJwu5T23FnqwHlfgkoUNEPNO1uzXee6rVl3kuVnr5i7qkXvt/2597qfO/u/ypy46w9/3T9V5J2v26+Yuqpa6PedmrYvWg/W4Ws4RyX7OAVGQXAoAAIqZW3iEDBDrdD+Bnnu721Ln+ZO299ysdfSw79ysdZ7vePTrzr0V3/20fe7XDpY6z3cLvOdmrYveg+W48o4HLlHIdzxqpFKWGNO7LvLcrHWR5zZWZ0mR3JJcOVa3JfXTWue5r+h1e4/pffxJLsURHepvtWxNm7Tur0YS5pK6qa+xbN97bpfaw9YUyf6FZi7LwTO5dO5rLNs96kr1YK3bswfrcbVmeZBcihaFW3h42Cu5tCR6iKGb07kX+qXJpSW1kPrZQg9W1uRSFh84qiYXHt4iXhCXoocYSC6NMXYLyaXAUbHwMIgcbGVFDzFEDhArvZ+aY0cOEANa1+TCw5qEGTkxkx5i6ObklVxag/Uia039jGxpDxF7JbkUrQu38PBIm7Tsb83Xl0zMnPqaNehhfV3f1hf6/qKkZHLp1NessWciqXevln1Z67yTS7cuUFl04MhC3k7bqZkLYa2LnG3hXRd5btY6j33tleNxrq5mtoV3XeS5WevmaryPPTkeOBqSSwEAQDGHTC71/tsf1+5/6Kv/vuOxd5vG3fo3Pfpjzo3r/dP9p//E737137d97OvP7mtuf2t63fr8WuvWPL9z+6vxjoe1B++/JWI9Xp5/56R/XkrT56b332qx9Lrn8zs3Lu944BKFfMfDM21y+ALQ138x8E64nBp3+AJk2d/cbaT9uuELe2f4Al+jV2sPNZ5fa91cbkL/YmBNpLT24Jmsaf1+8E7ztJ6bLfRqPa6W88R6zgHRXHRy6VZrkjDnXvD6j1n3t1X/Rd879dPa61Z7Pr+ljkP/QmPtwSP10zOzYs2YU4uO4WPeCacW3s+v9biS5YFLFm7h4aG7WJy78HWPRw6msqZ+zr2413bk5NLOkZNLO9bvh3MX9sgBYt3j1l4jIrkUrWty4eEt8gXRKvLCxKqF4xB5YYIyOAdw6Vh4GEQMtlpq6oN8JXg9fzWPg9fv02sGjM19iHSJFpJLvZ6LNbzOAT7jgaNqcuHRXaDOvbh0j5cKwVpiaeqndWFRsgfri/vS5NKSF42l4WLWi0rJHqwX625OEVM/l87J+vyW7ME6J2tyKXBU4RYeW9Mmh9unvtmH2z0TTudeYPqPWce01k0tPtYsSqzHYeni4tz2rcfB8vwOH9uaXDrcPrX4WLMo8U7ztH4/TO1vzaJk7vzrP7Z1wTPcbu3VMo71OCxdXJzbbj3ngCMJeTtt56bnHjZdFKx11+5/yPTC4DlujTGlq890WBYclv1Ffn6t4y5JmrWM+dSbbjS9+D/wzJ2mBYelB0tiZldnuYhbj5d1f5a67rNG585NzzEl316tx8E6pvUcsZ5zQAQklwIAgGIuJrl0qm7sbghL3dY0T3qYr7HWReph7BbGsZ9Ch3cuTP1Ea63zTNYcuyXVUrc1zbNGXa1ePY+/9ZwDjiLkOx6eaZ6l6tYkYdKDf92ePWxNLpVeemGxJpx6JmvO5WB41q1N/WyhB8txtZ4jJJfiqHZNLk0p3ZZS+pWU0qdSSp9MKf3F0/ZXppQ+nFJ65vT/N1v255EiuTTvwWO/3kmY9LCurnYPVnMXnv5j0ZJLrftdM7cWerAeV+CSedzV8qKkv5RzvkvSmyW9PaV0l6R3SPpozvl2SR89/XdR1sTMyOhhed0evJJLa158rKmfLSSXRu7Beo6QXIpWbV545Jw/l3P+56d/f0nS05JukXSfpPedyt4n6fu2joX1jrBAAQC0zzXHI6X0zZL+mKRflfTqnPPnTg89J+nVE1/zYErpyZTSk196/nnP6ZjDtyJrIfWzhR5aSC71CstqIbm0Zg8kl+LSuS08UkpfL+mDkn405/w7/cfy1SdYRz/FmnN+JOd8d8757lfcPP8xkL0SRksml3qPtTT1kx6W7WOv5FJrnTW4au0+xur22u9e+xiri9jD0nODhQVa5bLwSCndqKtFx6M553942vz5lNJrTo+/RtIXLPvyTrisURd5bta6yHOz1nmP6Z1cuiXhdHixs6Z51qqziN6DpSfvY09yKVq0+XbalFLS1Wc4fivn/KO97e+V9O9zzj+VUnqHpFfmnH9sbl8EiAEAcHx7B4j9p5L+W0n/KqX01GnbX5b0U5I+kFL6IUm/Ken7HcYCAAAHtnnhkXP+vyWliYe/c80+1yZcRq+LPDdrXeS5Weu27Gttcqk0/va6tW7PlM5adZHnZq0bq/E+9iSXojXh/jrt1G2fw+1HrIs8N2td5LlZ69buayo3Ybh9KqdhuN1aN3ZRHG6byqWIXBd5bta64TbvY28954AjCbfwmBM5MdNaRw/b6yP2YA0Q61jrPNI8l9bttd+99jFWF7GHpecGiwu06lALj6VI/Vxet4cWevC6CEROLi29n5pjR04utWJhgqNqeuGBrzlCWBoAoH1NLzwuIbn0knqInFxq/bBf5ORS79CtPXjN8QjJpV7nHBBNuIXH3MUlcnKpJbjKa6y99ksP6+s6R04u3Wu/NZJLrSInlwKt2hwg5qkfIHak2zinLk6X0oO1roUe9r6dduqiY7mN01pX4xZT77roPViOK7fTomVzAWJhFx4AAOCY5hYe4X7VAgAA2uURme5u6rbJ4Vvf1rfkr93/0HXb7njs3aZxLWOO1Y2NOTau968fLL1ueX7H6jyfX2vdlud3bH+1kkstz93UrZ/Dt/mtv1awHi/PX2dYj5f3r1A8e7UeB+uY/KoFlyrcOx5bMx2G26de8IbbrfuzzGdqTOu+1yZrWnu1zMfat/fzu/U49Me1BojVSi61PHfW4Cpr+qb1eHmmfs6df5Ze1yaNevZqPQ7WMUkuxSULt/CYs3QRcO6C2z1eMsBq66Jkqs57vx77WDon7/162Cu5tGQP1rCsbk7eqZ8els7J+vyW7ME6J5JL0bpDLTxqiZz6aVXyQjfk1QPJpdt4nQMtJJfW/H4guRSXjoWHwRFCus4Z+x3z0bRwHADg0jW98Dh3sT3Cxdia+hm5F2sPkRcW1hTJc+FPNcOhrN8PkVM/rXOL/L1vPUdILkWrwi08SqZS7jXm3Ita/zHPHqxq9erJ+/mtcRysPVgTMz2SNZeyjundqyfv57fG9wNwNGEDxLxTKfu/0516AdhzTO9xa4zpPW7k53eurv+79amfOq2JpNYevJM1Ld8Pw/15jjk37ppe5xYGll73fH7nxrWeJ5ZzDoiE5FIAAFDM3MIjZIBYp/sJ9Nzb3ZY6y0/a3mN619FDjDrrT5/dT7PnfqfvXdf9VD73LoD13YKlded+JeJZV6uHGseVdzzQkpDveHgna3omkkaqizw3a13kuY3VWVIkp26XHF5YvOu8E0490zxL1O3dQ63jSnIpjuhQf6tla9qkdX9rvt5jbvRg27733Nb0sDVF0prf0K+b+5r+Y1tDv9amg1rG2Zq+ap3Dnj1Yj8OaYzyH5FK0KNzCw4M1CXNpYmZJ9BDDXsmlJYPEll5oW0gujdjD0nODxQVa1eTCw1vEC+JS9LBNC8mlXhfZFpJLa/ZAcikuHQsPg8jBVlb0EEPNADEvkQPEjqCFcwDYosmFhzUJM3JiJj3E0M3JK7m0Jmvq55Ev7tYeIvdIcilaF27h4ZE2ad3f0q/3TsKkB7sIPWx9oe8vSuYWKGvqPBI4LdakeZass1gzpvfx2rpAZdGBIwt5O620PO8hYl3kua2pizw3a53HvpbmeEjzF5qldXM1e2VbeNdFntu5Ostx8D725HjgaEguBQAAxRwyudTz75csqavxN11qvxvjXdf6cajxjseav/0R7Z0R77/Bstc7GVK548A7HrhEId/x8EybnLsFs183fEHpDF9YaiRmWnuInPrZwnGYu32xfzHwTq60PieeCaeREknX1nknknofB0ud9ZwDorno5NKt+i82tRIzt6qV+unZX63j4NnDmkTSqYvd8DHvhNMp3omkeyacWvZV6zjUzHMBagu38PCwNDFz7kWllhZSPztHPg6dc2FNS5NLazh3wbamfkYOEDtCD+dYk0sJEMNRNbnw8Bb5gmgVeWFi1cJxiLwwQRmcA7h0LDwMpj5AdiQ1w7m8xq55HLx+n14zYMwrNKuF5NKaPXidA3zGA0fV5MJjaWKm9YJW8uK9V+pnjR6sdZGPwzlLk0tLLkCsYy1NLi158V46J2t9xONgTS4FjircwmNr2uRwu7Vu6qK35mJYKvVzba9b57ZnD5GPw9SFYLh96gKzZlEy13f/saWLi3PbrXWWcbwTSUv1sOdxsJ4j1nMOOJKQt9N2LGmTXZ3l4mGtu3b/Q2cvdEvmZq3z7MFSF72HyMfhqTfdaHrxf+CZO00XI0td9xmXc8+JdcwnXniL6SLsWWdJLq01N6nOcbDWWc85IAKSSwEAQDEXk1w6VTd2R4elbmsSpjXl0LMHzzGX1FnGrXUcPOvGbmEc+yl0eOfC1E+0teosx6vGmN7j1hjTu856zgFHEfIdjxJpnt51w4uUJeXQe25zt5taxl2T5rl1zOH+Iqe0bk0ulV56YfFOOLXUWY+Xdw+eqZ/WuW3tdc/jYO2B5FIcVZPJpWv3sdd++49ZUw495uS5X+80z6W5G97HIVpyaZ9HwuWaOk+1Uj89tXAcgKMJt/DwZE3MjIweltftwSu5tObF59xCsHs8cg/WuVl7rcHaA8mlaFXTCw8vR7i4n9NC6icA4PiaXnhYw7e27mdPLaR+evVQ8zi0kFzqdQ7U7MFr7JrfDySX4tIdauGx10Kh5ALEO52zRuqndR9L51TyOHiPVSO51Hus7njttd+lj60Zq6vba79bkFwKXAm38PBO6axRZ03frJW+atlfjTG967xTVb2TSy113kmY1uPl2YOV95jevU6xHK+1PZBcihaFvJ0WAAAc16FupwUAAO06dHLp2iTMWnWR52atizw3a92Wfa1NLpXG315fW+e5r1p1kedmrSsxJsmlaE24dzxueu7h0YvAcNvULa6R6yLPzVoXeW7WurX7mspNGG6fymkYbt9S57mvWnWR52at23tM6zkHHEm4hcecyMml1jp62F4fsQdrgFjHo26vsehhWd1eY7G4QKsOtfBYitTP5XV7aKEHr4tA5NTP0vupOXYLPbAwwVE1vfAAAACxNL3w8EourekS0leP0MO5D/NZP+wXOfWzZJjWWl5zPEIPXuccEE24hcfcxSVycqk1uIoettd79+oxVsf74j1Xt9dC4cg91Oi1dg/A0YQNEDvSbZxTFyd6KFu3pQdr3d63005ddDzratxi6l13ST1wOy2OaC5ALOzCAwAAHBPJpQAAIISQyaVTt00O3/qu8Rb/lrlZ6+ghRq9bftXiWTd126Slbu+5WevoYf8egKMI947H1kyHvZMwLfOhhxjJpd49WFMkPevmshosdXvOzVpHD/v3ABxJuIXHnL2SLUsGWNHDdN2Re+guBNYLQskLx9I50cM+9uoBOJpDLTyW8kq6jJz6WXo/NcduIbm05sWEHvz3U3NsFiY4qqYXHl6OEDR2Dj0AACJoeuFxScmlkXtpoQdrimTktEl6sO9nTy30AGwRbuFRMllzab333Dx7sKKH9XWe5i4a/ce86zzRw/o64JKFDRDbM5WyRMLl3nWR52atizy3ubr+79anLibD378fsS7y3Kx1kee2to4FDI6A5FIAAFDM3MIjZIBYp1vln1vhW+qsPzF4juld10IP1rrIc5O+9u7I3K9iLDXR6yLPzVoXeW7WOuu+gCNweccjpfQzkr5X0hdyzq8/bXulpMckfbOk/0fS9+ecn5/bT/eOx5aEQGud575q1UWem7Uu8tzG6qKkqnrX0cO+c7PWbUnyBSIp8bda/r6kewfb3iHpoznn2yV99PTfZ21NCLTub83Xe8yNHmzb957bmh6iJchGSJrdOjd6mN++tQcgIpeFR875Y5J+a7D5PknvO/37fZK+z2MsC2vyX+SEQHqIYemipGT6qlXkBFkregDasefttK/OOX/u9O/nJL16x7F2FfGCuBQ9bNNC+io9+O+n5tgsTHBURXI88tUHSUY/TJJSejCl9GRK6ckvPT/7EZBqWrh9jR5iaOF38zV74PkDjm/PhcfnU0qvkaTT/39hrCjn/EjO+e6c892vuPlml4GXJv9FvKDRQwzdnFpIX22hh3OsPUTusYUegDl7Ljwel/S207/fJunnLV/kkRDoac/0QmsPW3ulh/V1fVtf6PdMVfWos6CH/essWHTgyLxup/1ZSX9S0jdK+rykn5D0jyV9QNIflvSburqddvgB1JfoB4gtzayIWBd5bta6yHOz1nmPaUlf9a6rMeaedZHn5lFXa25AFCSXAgCAYg6ZXGr9GwYPPHPnV//96O1PT+7Ps65f411Xqgdr3SX1cJR3PKx1kd8FsNbRw7I64AhCvuNhTZEcXlA6/QvLVI21bniRqlFHD/v3sCW5VHrpxSBSEubaushzs9ZFnpu1zrovIJoSyaVurGmTcxcyT/1x5sb0rvNUqwfP/vacm/Wc25qb4J2EuWeyZrS50cP0voCjCbfw8NBdVM5d+Kx1NbTQQ6eFHs690B8hbdKrh8jhW0fo4ZwWegDmNLnw8Bb5gmhFDzFwsQDnAC4dCw+DuQ83HkXNHrzGrtmD1+/UW0j9pIcYY/M5DxxVkwuP7gJ17kJlrathrx5K9rp0ThF7WBouZr0YlLxoLJ0TPexjrx6Aowm38LCmTW69SA23W+rmxvSuWzq3JXVb59ZaD1uTS4fbrXWWcbyTMOnhsnoAIgp5O22nu6Pg3IXhgWfuNF2YPOu6zxvUqKOH/eZmPedueu5h04u/pc6S5rlXHT3sOzdrnVcPQBQklwIAgGIuJrlUGv+JduxuCEvd1iTMa/c/9JL/vuOxd6/en7UHzzGX1O3Za6kEWes5Z0mvXJNwWbLOcrxqjOk9bo0xveus+wKOIuQ7HpYUya1JmN51w4vZ8AWv03/h857b1JjWcdemftbotUZy6dxtkJa0SWvdnkmY1nPEuwfLOWLdn3VuW3utkUi6tg6Ipsnk0jlL8x6s9dYkzLkXvLnH1szJc79rUj89el2qVnKpp1qJmZ6sY1rPkcg9RD4OwNGEW3h4siZmRnZJPUTupYXEzHMLwe7xyD1Y52bttYYWziVgi6YXHl4iXxCtar7QemnhOADApWt64WEN39q6nz15jT31oboSSC71388aXudAC6mfNb8fWjiXgC0OtfCwBj0tvUCVXIBYX/CWzsl7vx77WDqnksfBuo/IyaXeY3XHq4XUT+u5F/l4sbBAq8ItPKZe6IfbvZMwPeumXvSG273n5jnu1sXF3r169mA957zTJi113kmY1uPl2cPcIqD/mPfz693rFMvx2vMcAY4m5O20AADguA51Oy0AAGjXoZNL1yaSTtWN3fkx9natNQnziRfe8pL/vufGf3K2plbdWI21zjrm3sdhan+W4zWW1zF2zo3dwmhJm9xatyZptNTc1tZFnttUXeTjABxFyHc8xi4Cw21Tt1YOt1vrpm43HW4f298Dz9x53faxi/Fw21hNrbonXniLuW7NmCWOw9j2qePVNxUSNtw+lZsw3O5dN/acDLfVmtuWushzG9se+TgARxJu4VEjudQ7TXTqYnzusTV1e+23Rg/ex8GaXDrHes5ZQ5+G9ed4JMMunZN3Dx5pnrV7iHwcgKMJt/Dw1ELq57mL916LAE9exyFy2mTp/dQcmx5ijM3CBEfV9MIDAADE0vTCwyu5tKapD31aH19atwev4xA5bTJi6NZeY9NDmbG9zjkgmnALj7mkyL2SSz1SP/uPzV3k91ooeOy3/1iNHvZKdd2yD+s519krTXTpY1vm5N3DXN1eY7VwHIBWhQ0Q43basnUt3E47dRyOdDvt1EUn8m2c1h486yLfTju2P26nxaWZCxALu/AAAADHRHIpAAAIIWRy6VR2wvCtb+9ftVjqpm5LHf5qYUtdlF/J1OrB83hN3ao7rIvyq5axuqnbJi11UX5NQQ/8qgXohHvHY2tI2NrkUmvdlP5FdWu2RoSE01o9eB6vrSFhEZJLreFbkVM/6YHkUqAv3MJjjndyaVfnkUhqtVfCaMkgsb168EgktTpKcqmH2qmfHugBaMehFh5LeSWT1kw49UocrZlc6jV2zePQQtokPfjvp+bYLExwVE0vPLwcIWjsnJoBYl5aOA4AcOmaXniQXBqDtYfIx8OaIhk5bZIe7PvZUws9AFuEW3jUSC611nunfnqkgy6tj9yDlTWR1Frncc4tZU3z9K7zRA/r64BLFjZAzJJcKq1Lr9yScCn5pn6urasxZqk6z0TSJXXWc25NKmWJNE/vushzs9ZFnpu1zrovIBKSSwEAQDGHTS596k03mm5nfOCZO013PIz9DZChJ154i+kuDGvdtfsfch+3xpilx/Ue03qOWM+5m557+OxdBZYayd6DdX+e49YY03vcyM+vdX/WMYEjCPmOx5bkUun6t9Itf3TMO5F06gXHc9waY3qPuzbh1DrmluRS6fpzzjNt0trDlsRM67hr/+jcljHHxvVOB7X0WuL5tY5Lcilacah3POZ+2uw/Zs10mHoRsPwkItnTPPuPze27/5h3OujWMb17nfqaue2evc6dI/3HloaEndu+NW2y38PWxMy5fc9t39pDv876/bA19Gu4fev3/prk0jX7nttOcilaFG7h4aG7qJx7EegerxmuNcWaSNo9bn3BK8n6/C5NXy3Z617JpSV7WDqnFlI/ty4A9rB0Tiwu0KomFx7eIi5Mloq4MFmq5nFoIW3S6xxoIfWz5vdDC+cSsAULD4MjhHSdM/Y7/KNp4Ti08Lv5yOFbR9BCD8AWTS48ug8OnrvYdo9HvKB1c7KmfkZcWFifX2uvNXmlTdZk/X6I3MM53dytvUZEcilaF27hYU2RLBWvvSalc+5Frf+YR9qohXVM716nvmaNNc+vd3Lp1hf6/tdbe/BOzNx6wfUe07tXi1rPr3VcCxYdOLKQt9NKL/1Q39yFwZJIKr30d7pzLw6WdFBrXY0xvcetMab3uNZzxHrO9X+3PncBsNZZeliThLl1zCX7s9QNP1dh6bXU82vdn/U41Hh+gUhILgUAAMXMLTxeVnoyVt7veFjrPH8i9/wbLLXqovfg+Td4ar3j4fmTdo25WevoYf864AhCvuPhmVw6FyLVr/NO8/Ssm7uN1FIXoQdrnXVf1kTSEsml0ksvBt5JmDXqIs/NWhd5btY6676AaC46uXSrNWmeHnWeIvSwtdf+11sTST2SS/uPbc1NWJOEWbLOgh72r7MgwwNHFm7h4aG7qJxbnHSPRwwIW5pcSg/7svwRL0tdTS30cI61h8g9ttADMKfJhYe3yBdEK3qIoYWLRQvJpTW10AOwBQsPg8jBVlb0sI3X79RbSP2khxhj8zkPHFWTC4/ug4PnQsa6xyNelJcml5YKI1tirx5KhcdJ83e39C1NLi150Vg6J3rYx149AEcTbuHhnVw6VWf9+jVpnh51U18zt91at3VukXvoH1eP5NK+qQvBcLu1zjKOd5onPVxWD0BEIW+n7XR3FJy7MHQfEj23mHjgmTvP1nSfIzh3oatZF3lu1jqvfS059pY66zl303MPn33x736XH7nOcgGL3msLPVjrWHDgKEguBQAAxRw+uVSa/gnUmkppqVuT0hmtLvLcrHXWfXkee8l+zlnSK9ckXEarizw3a13kuVnrrPsCjiLkOx6WFElrIunW5FKpfjro1rlZ647UQ43k0rnbIC1pk9a6GkmY3nX0sO/cxsYFImkyuXTO0lRTj/TNPdNBt87NWle7B6tayaWeaiVmeqKH9XXAJQu38PBkTS6NzJr6GRnJpfb97Ike7PvZUws9AFs0vfDwEvmCaEUPAIAIml54eAVN1QwY8xqbHrZpIW2SHvz3U3NsPuOBozrUwsMa9LQ0XKxk6qf3WDWSS2v34LGgtO4jcnKp91iXmPoZ+XixsECrwi08pl7oh9utiaTWOs9kzVppnpfUg+W4bl1cDLd7p01a6molYZbqwVIXoYc5EXoAjiTk7bQAAOC4DnU7LQAAaNfuyaUppXsl/S+SbpD0d3POP2X5OkuK5NjtsGNvr1+7/6Hrtt3x2Luv2zZ218SWZM3huJ5jTtV9+k/87nXbbvvY15/d35ak0b2fX+u4Y2NKtuTSsbyOsXNu7BZGS9rkVJ31uVuTcOk95tZxh+fm2HkZ5fm1jhvp+QWOYtd3PFJKN0j625K+R9Jdkn4gpXTXua8buwgMt01lcAy3j70IjG2fulVzuH2sbmzb2LieY45tH1t0jG2f6mFNryWeX+u41+5/6LrtY+fJcNtUSNhw+1RuwnC7tc763I3tr8aYW8YdOzeH2yI9v9Zxozy/wJHs/auWeyQ9m3P+tznn/yDp/ZLum/uCFpJLp158zj22Zk6dqUXH2v16JJJae91rv9bk0jlLk0utFwRrDx5JmEvn5L3fuXPTet4unZPH82utK3nMgRbsvfC4RdKne//9mdO2Is59oy+9MNbglVxqfYHfg1cPNQPEvFIka15UvM73yKmfVjW/973OFRYoOKrqHy5NKT2YUnoypfTkl55/vvZ0AADAjvZeeHxW0m29/771tO2rcs6P5Jzvzjnf/Yqbb3YdfOrDhtbHIziXp2HN2xj7IF8pXj1ETi6NGLo15HW+t5D6WfN73+tc4QOmOKq9Fx7/VNLtKaXXpZS+TtJbJT0+9wVzSZHeyaUdjwtf/7G5FzXrC97Si6x1YeHd65ylL+7e+507D7yTSzvWi4G1B2v4luecvPc7d25az9ulc/J4fq11JY850IJdFx455xcl/YikX5b0tKQP5Jw/ee7rxl7s1yaXTn2zD7dvSdYc2zY2rueYY9unXsSH26d6WNNriefXOu4dj737uu1j58lwW63kUutzN7a/GmNuGXfs3Bxui/T8WseN8vwCR0JyKQAAcEVyKQAACGH35NI1rCmS1uRSa51nsubUrZ+Wui1pnpfUg+W4TuV1DOtqJZda6qZum7TU7T03ax097N8DcBTh3vGwpkhak0utdZ7JmltzKdameV5SD5bjujUkbO/kUkvd1syHPedmrbOGb0XuYU6EHoAjCbfwmOOdXNrVlQyw8h6rq7ukHpYm027Zh/Wc6y4EJYPEvMeq0YNV7R5qHC8WF2jVoRYeS3lcoKS6iZleY9PDNl4XgRZSP+khxtgsTHBUTS88vNQMrvJCDwCACJpeeJwLiVoaNFaDV+pnTdYeIvdiTZGMnDZJD/b97KmFHoAtwi08Wkgu9a7bOjdrXe0erKyJpNY6j3NuKWsSpnedJ3pYXwdcsrABYsMP9U1dAIaf45i62Fjqhp8hmLpwRq6LPDdrnXVfnsdesp9zw9+tW26BnLroRK6LPDdrXeS5Weus+wIimQsQC7vwAAAAx3TY5NKn3nSj6XbGB56503QHi6XmiRfeYrp7omZd5LmV7mHJsbfUWc85yx0FNz33cPg6C3qIUwe0IOQ7HlMv/MO3vq2plJaEyy0pnbXqaiSSetet7WHLsR+rs55zUZIwvevoYd+5WetILkUrDvWOx9xPm/3HloaEWbcPWdM8veumvmZuu7Vu69wi92BNLu0/tjQk7Nz2EkmYS+rWzI0e2ugBiCjcwsNDd1E5t7hYmlxakjXNs0ZyqdVePXgFw1lETi61qp366YEegHY0ufDwFnFhshQ9bNNC2iQ9+O+n5tgsTHBULDwMIgdbWdFDDC38bj5y+NYRtNADsEWTC4/ug4PW5NKIFzRrmmfk1M8Weuh4pU3W1EIP57SQ+tlCD8CccAsPa4pkqbjzPdM8S11oI/Swtdf+19dKLt36Qr9nEqZHnQU97F9nwaIDRxbydlrppR/qm7sw9D9oOHehsdb1P0cwd7G01K1J89w6pndd9B4sx3VNcuncOdf/3frcBcCzbk0SZqm5WevoYf86IAqSSwEAQDFzC4+XlZ6Mlfc7Htfuf+ir/77jsXdP1nn+RF5jTO9xa4zpPa71HKn1joelB++ftK3Hy7PX/phz49Z4fq37sx6HGs8vcBQh3/HwTC4dvtj19V8QvFM6p8Ydvgh5jltjTO9xvccskVwqvfRiYE2utPbgnZhpGbfGmN7jRn5+rfuzjglEc9HJpVutSfOcW+z0HysV+mUd07vXqa9ZY83z65Fc2n9sa25C/+utPXgkZlrHtfAe07tXi1rPr3VcCzI8cGThFh4euovKuRfZ7vGI4VpLUz+3XlD2YH1+lyaX1nDuhf4IaZPW74fIPZzTzd3aa0TWc+nIxwmXrcmFh7fIF0SryC+0Vi0chxYuFi2kftbUQg/AFiw8DCIHW1nNfcDtKGoeB6/fqdf83bzXOdBCcmnN74cWziVgiyYXHt0HB8+9uHSPR1xYLE39jLiwsD6/S5NLS/Y6d3dL39LUz5I9LJ2Ttb7khW/pnKzPb8kels6JhQVaFW7h4Z1cOvXNbn0RWJPmObfv/mNbUz/XXKi9k0vXPI9z2z179Ugu7Zu6EAy3W+um9HvwTsK0fj9s7aFfZ/1+2LrgGW7f+r2/JpF064LH+1wCIgp5O22nu6Pg3IWh+zDpucXItfsfOvvC0H2O4NyFzlrXfbbCc9waY1rqPMf1HtN6jljPue739JbshXMXCWsP1v15jltjTO9xIz+/1v1ZxwSiILkUAAAUc/jkUmn6J1Dr3+Hw/Jsea/5+iXddjTFL1a05ph511nPOkl65Jmk0Wl3kuVnrIs/NWmfdF3AUId/xsKRIzoVD9S8s3nXeqZ9b0jy96yL0sPV4DRcVnsmlc7dBWtImrXVrEzO3JGF619HDvnMbGxeIpMnk0jlLU02t9d6pn9a6rXOz1tXuwcqaSOqdXOrJOzHTWueJHtbXAZcs3MLD07kFRanY9S2sqZ+RWXuIfDy8kksjh2/RQxkt9ABs0fTCw0vkC6LVERYo57RwHADg0jW98LBmfZTazxpe4WY1Q9K8xq55HFpIm6QH//3UHJvPeOCoDrXwsAY9WS9QXZ213uMCujU0bKrOe78e+1hat/S4bWHdx17JpZFTP+lhH3v1ABxNuIXH1sXFcLt33RRr6ufU18xtr1FXqwfP47V1cTHc7p02aanbmly659ysdfSwfw/AkYS8nRYAABzXoW6nBQAA7Tp0cunYXQ5jb69b67q/rdA39ncWPBNOx+42qVG3JWnUOmat42CpG8vrGDvnxm5htKRNWuum3kIfPidjz8fec7PW1UjzLFVnOQ5j+ysxN+AoQv6qpUZy6djFrtN/cbEmYUZKJLXWRUgk3XochvuzHK+5kLAayaXDuqnnpP981JqbZzpotLkN6yIfByCaQ/2qpUZy6dzFzrrf/mMeaZ5HSS71mFPH4zhY66z7WJpcag11stbNPSfW52vpnLx78EjzrN1DjeMAtCrcwsMTyaXL6/bgdRysL/B78EqRjJyYWXo/Ncc+Qg8kl6JVTS88AABALE0vPM7lONRMwrQ6l6dRMjhsLa/jMPVBvhJaSJukB//91Bybz3jgqMItPOYCxLyTSzvWC5pHwuleC4W90kRLjNXxPg5zdUdJLp17TqzPV+3UT2v4ludYl3QcgKMJt/CQxl/sh9u8E0mnXjyG28f29+jtT1+3feyiGyGRdG5u1ro1Y5Y4DmPbp45XX+TkUmn8ORlui5z6OVUXeW5j2yMfB+BIQt5OCwAAjutQt9MCAIB2hUwutaZIeieXetatTd/cOjfPca1habV69ewhSnLpWJ0lfGvJmNbj5dmDJXzLe8ypcdf26n0cSC7FpQr3jsdUYNNw+9QFZbi9Rt3Ui+xwu/fcPMfdmq2xd6+ePVjPuakLz3C7Z93WzIfhduvx8uzBGr7l/fx69zrFcrz2PEeAowm38JizV3KpRxKm1V7pnN779djH0jmVPA7eyaV7JWZ67GNpOmfEMLSlc4qYJrpX+ipwNIdaeCzllZhZM+HUa+yaqZ9ePdQ8Di0kZnqdA0dI/Twncgpu6f0ApTW98PByhKCxc2qGb3lp4TgAwKVreuHRQnLpJfUQuZdzH+azhj7V/FDgucVn93jkHqxzs/ZaQwvnErBFuIVHjeRS7yRMj5TDGumr/cdK9rqUdw8e59xS1jRP7zpP1jGt50jkHiIfB+BowgaIDT/UN3UBGP7e3/P2zKkLl7Vu+HvkqRdgzx48x1xSt2evW+dmrbOec8PfrVtugZy66NSqsxyvGmN6j1tjTO86676ASOYCxMIuPAAAwDEdNrn0qTfdaLqd0fvuFEvdA8/cWa3Ogh7W1VnPOc/bMG967uFqdRb0sH+dBXexoBUh3/GYeuEfvvVtTaWskYRZoo4e/OZmPeeiJJJ619HDvnOz1pFcilYc6h2PuZ82+49t/anZOwnTu27p3JbUbZ1baz0sDQk7t71EEuaSujVzo4c2egAiCrfw8GBNwlyamFnSXj2U7HXpnCL2EDm51GqvxEx6WIbkUuBKkwsPbxEXJku1kL5KcmmMsekhxtgsTHBULDwMIgdbWdFDDPxuHpwDuHRNLjysSZiREzNb6KHTQg9eaZM1tZCY2UIP57TQAzAn3MLDmiJZ6iK1JgnTo85TrR48+9tzbtZzbusL/Z5JmB51FvSwf50Fiw4cWcjbaaV1yaVzFxrPuj2TNUv1YK27pB7659zcYqT/u/W5C4Bn3ZokzFJzs9bRw/51QBQklwIAgGLmFh4v27LjlNKflfRXJN0p6Z6c85O9x94p6YckfUXSf5dz/uUl+7b+9Bm5LvLcrHWR52at8x6z9k/a0d4FWFMXeW4edbzjAUzb9I5HSulOSb8v6X+X9N93C4+U0l2SflbSPZJeK+kjkv6TnPNX5va3NLnUUjeXw2CpWzNmrbrIc7PWRZ6bZEuljJSEubYu8tysdZHnZq2z7guIZrfk0pzz0znnayMP3Sfp/TnnL+ecf0PSs7pahJxlTS61BjttZR3Tu27qa9agh/V1fVtzE7yTMPdM1ow2N3qY3hdwNHvd1XKLpE/3/vszp21FdBeScxcUa10N9BBDN6dzL/RHSJtsoYdzrD1E7rGFHoA5ZxceKaWPpJQ+MfK/+zwmkFJ6MKX0ZErpyS89/7zHLt1FvCAuRQ8xtHCxaCH1s6YWegC2OLvwyDl/V8759SP/+/mZL/uspNt6/33radvY/h/JOd+dc777FTffvGz2hcx90PAo6GEbr9+pRw7fKr2fmmPTA1DPXr9qeVzSW1NKL08pvU7S7ZKe2Gms63QXqHMXKmtdDfQQg3VOS5NLS140ls6JHvaxVw/A0WxaeKSU/kxK6TOSvl3S/5FS+mVJyjl/UtIHJH1K0v8p6e3n7mjpWFMkrReEqbo1X+8xN3qwbd97bmt6mLoQDLdb6yzjeCdh0sNl9QBEFDpArPud/rkLg6VuaZaDx5jedS30YK2LPDdpWR7FuYtE5LrIc7PWRZ6btc66LyAKkksBAEAxu+V4AAAALMHCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFLNp4ZFSem9K6ddTSr+WUvpHKaVv6D32zpTSsymlayml7948UwAAcHhb3/H4sKTX55zfIOlfS3qnJKWU7pL0Vkl/RNK9kn46pXTDxrEAAMDBbVp45Jw/lHN+8fSfH5d06+nf90l6f875yznn35D0rKR7towFAACOz/MzHn9B0i+d/n2LpE/3HvvMadt1UkoPppSeTCk9+aXnn3ecDgAAiOZl5wpSSh+R9E0jD70r5/zzp5p3SXpR0qNLJ5BzfkTSI5L0La9/fV769QAA4DjOLjxyzt8193hK6c9J+l5J35lz7hYOn5V0W6/s1tM2AABwwbbe1XKvpB+T9Kdzzr/Xe+hxSW9NKb08pfQ6SbdLemLLWAAA4PjOvuNxxt+S9HJJH04pSdLHc84/nHP+ZErpA5I+patfwbw95/yVjWMBAICD27TwyDl/68xj75H0ni37BwAAbSG5FAAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFMPCAwAAFLNp4ZFS+smU0q+llJ5KKX0opfTa0/aUUvqbKaVnT4//cZ/pAgCAI9v6jsd7c85vyDm/UdIvSPrx0/bvkXT76X8PSvo7G8cBAAAN2LTwyDn/Tu8//5CkfPr3fZL+Qb7ycUnfkFJ6zZaxAADA8b1s6w5SSu+R9IOSflvSf37afIukT/fKPnPa9rmRr39QV++K6Btf+9qt0wEAAIGdfccjpfSRlNInRv53nyTlnN+Vc75N0qOSfmTpBHLOj+Sc78453/2Km29e3gEAADiMs+945Jy/y7ivRyX9oqSfkPRZSbf1Hrv1tA0AAFywrXe13N77z/sk/frp349L+sHT3S1vlvTbOefrfs0CAAAuy9bPePxUSukOSb8v6Tcl/fBp+y9K+lOSnpX0e5L+/MZxAABAAzYtPHLO/9XE9izp7Vv2DQAA2kNyKQAAKIaFBwAAKIaFBwAAKIaFBwAAKIaFBwAAKIaFBwAAKIaFBwAAKIaFBwAAKIaFBwAAKIaFBwAAKIaFBwAAKCZd/VmVGFJK/6+u/tjcFt8o6d85TOfS8Lwtx3O2Ds/bOjxv6/C8rbP1efuPc86vGnsg1MLDQ0rpyZzz3bXncTQ8b8vxnK3D87YOz9s6PG/r7Pm88asWAABQDAsPAABQTIsLj0dqT+CgeN6W4zlbh+dtHZ63dXje1tnteWvuMx4AACCuFt/xAAAAQTWx8Egp/WRK6ddSSk+llD6UUnrtaXtKKf3NlNKzp8f/eO25RpJSem9K6ddPz80/Sil9Q++xd56et2sppe+uOM1wUkp/NqX0yZTS76eU7h48xvM2I6V07+m5eTal9I7a84kqpfQzKaUvpJQ+0dv2ypTSh1NKz5z+/+aac4wmpXRbSulXUkqfOn1//sXTdp63GSmlP5BSeiKl9C9Pz9tfPW1/XUrpV0/fq4+llL7Oa8wmFh6S3ptzfkPO+Y2SfkHSj5+2f4+k20//e1DS36kzvbA+LOn1Oec3SPrXkt4pSSmluyS9VdIfkXSvpJ9OKd1QbZbxfELSfynpY/2NPG/zTs/F39bV9+Vdkn7g9Jzhen9fV+dQ3zskfTTnfLukj57+G1/zoqS/lHO+S9KbJb39dH7xvM37sqTvyDn/UUlvlHRvSunNkv66pL+Rc/5WSc9L+iGvAZtYeOScf6f3n39IUvfBlfsk/YN85eOSviGl9JriEwwq5/yhnPOLp//8uKRbT/++T9L7c85fzjn/hqRnJd1TY44R5ZyfzjlfG3mI523ePZKezTn/25zzf5D0fl09ZxjIOX9M0m8NNt8n6X2nf79P0veVnFN0OefP5Zz/+enfX5L0tKRbxPM263R9/N3Tf954+l+W9B2Sfu603fV5a2LhIUkppfeklD4t6QF97R2PWyR9ulf2mdM2XO8vSPql07953tbheZvH87PNq3POnzv9+zlJr645mchSSt8s6Y9J+lXxvJ2VUrohpfSUpC/o6p3wfyPpi70fTF2/Vw+z8EgpfSSl9ImR/90nSTnnd+Wcb5P0qKQfqTvbOM49b6ead+nqbcpH6800FsvzBtSSr25H5JbEESmlr5f0QUk/Ong3nOdtQs75K6ePKtyqq3cmv23P8V6258495Zy/y1j6qKRflPQTkj4r6bbeY7eetl2Mc89bSunPSfpeSd+Zv3ZvNc+b/Xzru/jn7Qyen20+n1J6Tc75c6dfGX+h9oSiSSndqKtFx6M553942szzZpRz/mJK6VckfbuuPprwstO7Hq7fq4d5x2NOSun23n/eJ+nXT/9+XNIPnu5uebOk3+695XbxUkr3SvoxSX865/x7vYcel/TWlNLLU0qv09WHc5+oMceD4Xmb908l3X76tPzX6eqDuI9XntORPC7pbad/v03Sz1ecSzgppSTp70l6Ouf8P/Ue4nmbkVJ6VXdHY0rpD0r6L3T1+ZhfkfRfn8pcn7cmAsRSSh+UdIek39fVX7f94ZzzZ08n4t/S1afDf0/Sn885P1lvprGklJ6V9HJJ//606eM55x8+PfYuXX3u40VdvWX5S+N7uTwppT8j6X+V9CpJX5T0VM75u0+P8bzNSCn9KUn/s6QbJP1Mzvk9dWcUU0rpZyX9SV39hdDP6+od3H8s6QOS/rCuXue+P+c8/ADqxUop/WeS/i9J/0pX1wJJ+su6+pwHz9uElNIbdPXh0Rt09WbEB3LO/2NK6Vt09QHwV0r6F5L+m5zzl13GbGHhAQAAjqGJX7UAAIBjYOEBAACKYeEBAACKYeEBAACKYeEBAACKYeEBAACKYeEBAACKYeEBAACK+f8BllMA+Ag0K90AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 645.161x662.252 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# set the mgxs environment variable for plotting, use openmc.config in the future\n",
    "os.environ['OPENMC_MG_CROSS_SECTIONS'] = os.getcwd() + '/c5g7.h5'\n",
    "root_universe.plot(origin=(0., 0., 0.), width=(3 * 21.42, 3 * 21.42), pixels=(500, 500),\n",
    "                   color_by='material')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, it looks pretty good, let's go ahead and write the file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Geometry and set root Universe\n",
    "geometry = openmc.Geometry(root_universe)\n",
    "\n",
    "model.geometry = geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now create the tally file information. The tallies will be set up to  give us the pin powers in this notebook. We will do this with a mesh filter, with one mesh cell per pin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "tallies = openmc.Tallies()\n",
    "\n",
    "# Instantiate a tally Mesh\n",
    "mesh = openmc.RegularMesh()\n",
    "mesh.dimension = [17 * 2, 17 * 2]\n",
    "mesh.lower_left = [-32.13, -10.71]\n",
    "mesh.upper_right = [+10.71, +32.13]\n",
    "\n",
    "# Instantiate tally Filter\n",
    "mesh_filter = openmc.MeshFilter(mesh)\n",
    "\n",
    "# Instantiate the Tally\n",
    "tally = openmc.Tally(name='mesh tally')\n",
    "tally.filters = [mesh_filter]\n",
    "tally.scores = ['fission']\n",
    "\n",
    "# Add tally to collection\n",
    "tallies.append(tally)\n",
    "\n",
    "model.tallies = tallies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the geometry and materials finished, we now just need to define simulation parameters for the `settings.xml` file. Note the use of the `energy_mode` attribute of our `settings_file` object. This is used to tell OpenMC that we intend to run in multi-group mode instead of the default continuous-energy mode. If we didn't specify this but our cross sections file was not a continuous-energy data set, then OpenMC would complain.\n",
    "\n",
    "This will be a relatively coarse calculation with only 500,000 active histories. A benchmark-fidelity run would of course require many more!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OpenMC simulation parameters\n",
    "batches = 150\n",
    "inactive = 50\n",
    "particles = 5000\n",
    "\n",
    "# Instantiate a Settings object\n",
    "settings = openmc.Settings()\n",
    "settings.batches = batches\n",
    "settings.inactive = inactive\n",
    "settings.particles = particles\n",
    "\n",
    "# Tell OpenMC this is a multi-group problem\n",
    "settings.energy_mode = 'multi-group'\n",
    "\n",
    "# Set the verbosity to 6 so we dont see output for every batch\n",
    "settings.verbosity = 6\n",
    "\n",
    "# Create an initial uniform spatial source distribution over fissionable zones\n",
    "bounds = [-32.13, -10.71, -1e50, 10.71, 32.13, 1e50]\n",
    "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n",
    "settings.source = openmc.IndependentSource(space=uniform_dist)\n",
    "\n",
    "# Tell OpenMC we want to run in eigenvalue mode\n",
    "settings.run_mode = 'eigenvalue'\n",
    "\n",
    "model.settings = settings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's go ahead and execute the simulation! You'll notice that the output for multi-group mode is exactly the same as for continuous-energy. The differences are all under the hood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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-2022 MIT, UChicago Argonne LLC, and contributors\n",
      "         License | https://docs.openmc.org/en/latest/license.html\n",
      "         Version | 0.13.1\n",
      "        Git SHA1 | 33bc948f4b855c037975f16d16091fe4ecd12de3\n",
      "       Date/Time | 2022-10-05 22:19:28\n",
      "   MPI Processes | 1\n",
      "  OpenMP Threads | 2\n",
      "\n",
      " Reading settings XML file...\n",
      " Reading cross sections HDF5 file...\n",
      " Reading materials XML file...\n",
      " Reading geometry XML file...\n",
      " Loading cross section data...\n",
      " Loading uo2 data...\n",
      " Loading mox43 data...\n",
      " Loading mox7 data...\n",
      " Loading mox87 data...\n",
      " Loading fiss_chamber data...\n",
      " Loading guide_tube data...\n",
      " Loading water data...\n",
      " Reading tallies XML file...\n",
      " Preparing distributed cell instances...\n",
      " Reading plot XML file...\n",
      " Writing summary.h5 file...\n",
      " Initializing source particles...\n",
      "\n",
      " ====================>     K EIGENVALUE SIMULATION     <====================\n",
      "\n",
      " Creating state point statepoint.150.h5...\n",
      "\n",
      " =======================>     TIMING STATISTICS     <=======================\n",
      "\n",
      " Total time for initialization     = 1.6491e-02 seconds\n",
      "   Reading cross sections          = 8.1206e-03 seconds\n",
      " Total time in simulation          = 1.0349e+01 seconds\n",
      "   Time in transport only          = 1.0309e+01 seconds\n",
      "   Time in inactive batches        = 2.4253e+00 seconds\n",
      "   Time in active batches          = 7.9239e+00 seconds\n",
      "   Time synchronizing fission bank = 2.9981e-02 seconds\n",
      "     Sampling source sites         = 2.5394e-02 seconds\n",
      "     SEND/RECV source sites        = 4.4660e-03 seconds\n",
      "   Time accumulating tallies       = 3.2036e-04 seconds\n",
      "   Time writing statepoints        = 1.7235e-03 seconds\n",
      " Total time for finalization       = 7.8247e-04 seconds\n",
      " Total time elapsed                = 1.0370e+01 seconds\n",
      " Calculation Rate (inactive)       = 103082 particles/second\n",
      " Calculation Rate (active)         = 63100 particles/second\n",
      "\n",
      " ============================>     RESULTS     <============================\n",
      "\n",
      " k-effective (Collision)     = 1.18598 +/- 0.00158\n",
      " k-effective (Track-length)  = 1.18629 +/- 0.00194\n",
      " k-effective (Absorption)    = 1.18569 +/- 0.00111\n",
      " Combined k-effective        = 1.18574 +/- 0.00111\n",
      " Leakage Fraction            = 0.00184 +/- 0.00006\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Run OpenMC\n",
    "statepoint_filename = model.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Results Visualization\n",
    "\n",
    "Now that we have run the simulation, let's look at the fission rate and flux tallies that we tallied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAEICAYAAADhtRloAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAl+ElEQVR4nO2df5RdVZXnPzs/CAmJFCQSAgTCj7SKqAEj4ICC0LSAP5DVaoM9iI7LODY42ktnibYjoR271W6xccEgUWjQhYItopGm20YGENY0SMAA4ZckWBhCfhBiQULIj6rs+ePegqqXs0+99+rVq3cr389ab9V7+969z7nvx65zzz57H3N3hBCiSowb7Q4IIUSjyHEJISqHHJcQonLIcQkhKocclxCicshxCSEqhxxXhTGzTWZ2yGj3Q4h2I8fV4ZhZt5m9VDqptWZ2tZlNBXD3qe7+ZBM255iZlzY3lW1c0PreCzEyyHFVg/e4+1TgKGA+8KUW2e0q7Z4NfNnMTm2R3boxswntblNUHzmuCuHuq4B/A44AKEdNh5XPrzazy8zsX81so5ndY2aH1mn3P4GHgSPMbJyZfcnMnjKzdWb2fTPbs2zjGjP7bPl8/7L988rXh5rZBjMbV75+t5ktNbMeM/t/ZvbG/vbKEd7nzexB4EUzm1C+XlX2/XEzO7llb5wYc8hxVQgzmw2cDvw2OOUs4CJgL2A58NU6bJqZHQe8vrT7kfLxDuAQYCpwaXn6HcCJ5fMTgCeBtw94fae77zCzI4GrgE8A04ErgMVmNmlA02cD7wK6gEOB84G3uPs04J1A91B9F7suclzV4Gdm1gPcReE8/i4470Z3/4279wLXAvOGsLse2AB8D7jA3W8F/hK42N2fdPdNwBeAs8pbujuA48tR1duBbwDHlbZOKI8DLACucPd73L3P3a8BtgLHDmj72+6+0t1fAvqAScDhZjbR3bvdfUUd74vYRZHjqgbvc/cudz/I3f+q/LGnWDPg+WaK0VKOGe6+l7u/zt2/Xcr2A54acM5TwARgZulMXqRwiG8DbgKeMbPXMNhxHQR8trxN7Cmd7uzSdj8r+5+4+3LgM8BCYJ2ZXWdmA88VYhByXKKWZygcTz8HAr3A2vL1HcD7gd3KObc7gHMpbk+XluesBL5aOtv+xxR3/9EAu4PKkrj7D939+LJtB77e2ssSYwk5LlHLj4C/NrODy2UXfwdcX95+QuGozgd+Xb6+vXx9l7v3lbLvAv/dzI4p59D2MLN3mdm0VINm9hozO6mcA9sCvATsGJGrE2MCOS5Ry1XADygc0+8pHMmnBhy/A5jGK47rLmDKgNe4+xLg4xST+n+kCBR8JNPmJOBrFHNua4B9KObWhEhiKiQohKgaGnEJISqHHJcQonLIcQkhKocclxCicrQ1wXXGePM5qRYnNmFse+ZYdFXjA3mu/a2BfFIgz7FH5ljU56Cd3j0sNLWd3ZLyvvANgBdIrlRgAn1J+XNMD21t9XSn+7Znvm7bg+tpZlHElkAeLduF+F94TifqWyTflrEVxcjijxn2TMg2deNb1ue0huQwM99c57mr4Zfu3vbk/LY6rjkTYMn+iQP7NGFsXeZY9Jt6VSCfmbEVFY3JVcGK/MOxgRzSX0KAuWnxhmNiz7mS2Un5xsA5AdzKnybl03kuKf/nzOqGFVsPS8qfXxM7O9YE/z02xSohywP50oxOlGOwLKMT9S1ynE9nbEU6u2d0Uu7iF/MzCvXxEnBened+CWYMu8EmUEkRIcQgjOZugtqJHJcQYhBG5zuGTu+fEKLNaMQlhKgc44DJo92JIWiv45oEHJyQ5yYg9w7kz2R0gonOFx5Ky1+V6lM/0ST83RmdaBI+pxNFHF9Mi/feM5rNhXWvTYevnsuEQjczJSnvoicp/wuuD23dOentSfnmg+Kfw7aD0n3roSspX7YiMwk9Iwiq5YIjvwrk+2Z0ogn96HPOTWNHE/1dGZ0Rogq3ilrHJYQYRP+tYj2PIW2ZzTaz28zsETN72Mw+nTjnRDN7viz1vdTMvjyU3U53rEKINtPiEVcv8Fl3v78sa3Sfmd3i7o/UnHenu7+7XqNyXEKIQbRyct7dVwOry+cbzexRYH+g1nE1hG4VhRCDaPBWcYaZLRnwWBDaNZsDHAnckzj8VjN7wMz+zcxeP1QfNeISQgzCaCiquN7dh1yuX1bTvQH4jLu/UHP4fuAgd99kZqcDPyPMGSnQiEsIMYj+Oa56HnXZM5tI4bSudfef1h539xfKHaVw95uBiWaWTSVq74hrK0Ux4AZ6sfaJtLzWZQ9k72CpxPTedCbrExbnpE54LC0/OFc59sjAXpSPCHB7YO+jga0o7xKYM/uppHzjHnGu4uxXNt0ZRC6ZOmI665PySZlcySiP8k5OScrHTw3WiQA7NqXXlvibQhUsym/siXX8Y4GtZpZDRORyNVPLO5pJ/q+hlXNcZmbAlcCj7n5xcM6+wFp3dzM7mmJAlU6SLdGtohBiEC2OKh4HnAM8ZGZLS9kXKXaPwt2/Q7Fr1CfNrJcix/ssH6KmvByXEGIQLY4q3kW+OA/ufimv7JZeF3JcQohBKOVHCFE5qpDy0+n9E0K0GVWHqGUc6WTiTJJ1FFo4PDN3d1cQJZwYyI9JFwwtiPr2usxte5SYfUKmnY8H9qLS0WfH17/7yrStKXvEBXmXMi8pjyqgHhaWGYX1QSTyaj4Z6vwJDyTlJ3NToPGO0FbE5Oc3xAd7gmz+TFTPfh4cCL4zfnnG1tuCA6+NdZJJ3nF17rqpguMach2Xme1uZr8pV7U+bGYXlfKDzeweM1tuZtebWbrQuRCicrRyHddIUM8C1K3ASe7+JmAecKqZHQt8HfiWux9Gsc16sKpFCFElDJg4ob7HaDGk4/KC/gFzf4qSAycBPynl1wDvG4kOCiHay7hxMHlSfY/Roi6faWbjgfuAw4DLgBVAj7v3lqc8TZHxndJdACwAOLDTb5yFEJjBhA4P29XVPXfvA+aZWRdwI/kpw1rdRcAigPmTLbsaVggx+vTfKnYyDXXP3XvM7DbgrUCXmU0oR10HAKtGooNCiDZjtCQ6OZIM6bjM7NXA9tJpTQZOoZiYv40ix+g64FwgCg6/wmTgdQl5phZ7Ol0YfpxJjP7LQN4dKeR2mI6Ic3zjDz0uE9/whrDcmFmOEcw9bJwdJzkfyoqkPKr5vj6TMdzDXkn5p/hGqDObw5Pyl4Ja+FmCb/WW5dEGBrFOlqWNnW7pfPGCrkDek9FJbTCb2y27XiqwArWe7s0CrinnucYBP3b3m8zsEeA6M/vfwG8pMsCFEFVnLDgud3+QomphrfxJ4OiR6JQQYhQxWlIeZyTpcL8qhGg7Y2HEJYTYxZDjEkJUkqpHFVvKS0BqN+lMkvVpQTL1YZmoYsQpga3tXbGtifsEB/6QWZL2zsDeuliFSwN7lwS2chHKM9O2ZhFdDPQF39T9gi3D92FtaOsLfCsp/zwXhTozgzfnh/y3pNy2bA9t0ZsW+5tjlbB0c7oKdWHvK4GtLwQKx2dsXRjY+misk0wA35E5v1404hJCVA5NzgshKodGXEKIyiHHJYSoHHJcQohKoqjiAMaRLjebyxV8QzqqNrfu+hSvEEUPJx6cUYqOHZuJar4hkO+XaefCwF4UCYvKQwMsS9taecRRocpW0gVsNwad7mZOaOuXpGsUrw7KQwM8zmuS8kN5OK3Q8/rQFmvS4mwcOooeZiLedllwIJVDCNlIsEXfmdz3PBU97cucXy8acQkhKsc4FFUUQlQMjbiEEJWkwz1Dh3dPCNF2xkIhQSHELoZuFYUQlUOOq4bdaWCbjZInAnluKBuEsMNlD5mKvmE4PLcc4T2B/LGMTpSAnVo+MgQeLLuYEe4LDtPYmJSvDGoKd2VqCufKOjfK+ChjemomyX1CsPAh922PPufc+x8toYjaCS4FiEs3Z3bSTvat8doDaRuKKgohKoVGXEKIylEBxzXkTtZCiF2Q8XU+hsDMZpvZbWb2iJk9bGafTpxjZvZtM1tuZg+aWZziUdLhflUI0XZaO+LqBT7r7veb2TTgPjO7xd0fGXDOaRQb8c0FjgEuL/+GyHEJIQYzjmyOZiO4+2pgdfl8o5k9CuwPDHRcZwDfd3cH7jazLjObVeomaa/j6gU2JOS5CMaqIHo0MxM+id703wa2Ppqx9fu02O4Iau0CPi8oUZwbWv9T0LfLgr49H5uyI9K2Huedoc7vgiTnKHo4KbPz6A/4eFL+CS4JdTYHG7/eSXoXVcuVro6CmkcsjHV+FhyLSjoD/vdpuZ0TKHRlbN3YoC1IJ223onQzNLIAdYaZLRnwepG7L0qdaGZzKLY6vKfm0P7AygGvny5lzTsuM5sNfB+YCXjZsUvMbCHwceDZ8tQvuvvNQ9kTQnQ4jd0qrnf3+UOaNJsK3AB8xt1faL5zBfV0L3mPWh77lrv/43A7IYToIFocVTSziRRO61p3/2nilFXA7AGvDyhlIUNGFd19tbvfXz7fCPTfowohxiqtiyoacCXwqLtfHJy2GPhwGV08Fng+N78FDfrVmnvU44DzzezDwBKKUdkfEzoLgAUAB7Zowk8IMYK0dsR1HHAO8JCZLS1lXwQOBHD37wA3A6dTzChuBnKbskEj3au9RzWzy4GvUMx7fQX4Juy8CV45UbcIYP6elsnTEEJ0BC0sJOjudzFEIlIZTTyvEbt1Oa7UPaq7rx1w/LvATY00LIToUCqwcr6eqGLyHrVmncWZwLIhW5sMHJGQvy6j897AWZ+Q0Ymu6sjA1skZW/PSYj8z3pWZYwN5elPogquDvh2XFm/P1Mlfx/Sg+Sj7GzYyLSmfHiRmL+fQ0NbJwf+wKZmi+1H7b+DetEJXVKQdeCyYk7hpYawT1YnvilXsU8GBKJk6Zyu6nIwO+2aODZeqOy7ie9SzzWwexa1iN/CJEeifEKLdjIURV+YeVWu2hBiLqAKqEKJyjIURlxBiF0OFBIUQlUMjrhqmkI7S1aZcDmRmIM/9R9gayKNCGbmdtKPo3YsZnehYLnoaJQ0HgbiNe8areaOE5Vy55clsTspXBNHD3TJJ1nsF7fRmJk56gvDZxr50tJH1mdXMUYQwV1E6OpYrndxoVK8ncyy6nNwv9ICEbGLdvYmR4xJCVBJNzgshKoVGXEKIytHCQoIjhRyXEGJndKsohKgUulUUQlQOOa4a+oBU0dY9MzpRLfZcnfhow+bFga3PZWxFSxvOzlTouTmw91CswgcCe8vStjbuEywTAOYOKt/9Cj/k86HOykEFKF9hAn1J+Ry6Q1tXsNMOVAB8iKtCnfFBO8+NT9estK7toS260msC/F2xil0fHMjUtve/CWxFydfRzteAB7nk9oFYJ7m8Iv02No5uFYUQlUIjLiFE5VBUUQhRRVy3ikKIKuEGfR3uGTq8e0KItiPHVcNEYJ+EfG5G58IgQhclUgNE0aOPN1EG+shA/n+biERGJZ2BLS+m7U0Kkqz7MmGfe0nXAe7iv4Q624Ks9SlBZvCUICkb4D38S1K+ksNDndlBJHR6X1daYVNmh7wgEmj/GquEidmRHLD/FRyIErMzEUo7LTgQlYGGdFnn7LYU9eEGveOH3LmwpFVbZzdGh/tVIUS7cTP6JtTrGuIqISOJHJcQYhCOsW38bnWeLcclhOgAHMvWTusE5LiEEDvR1+GuobN7J4RoO45lgz+dQHsd11TwRATv5d0aU0R7iF6UyRW8IgitpDajhXxU81VpsaeioyX2ZHAgyqEEdg++J2v2SydyHsLqpBxgBW9LyteGdbBhM5OT8ija90gmQhiVjl6aCat2bV2TlO82KZhDyX1zozLMPRmdVBlkyEYCQ4J2/LZYJdwQNtd+V0LWglzFMeG4zGw28H2K6u8OLHL3S8xsb+B6YA7FhrAfdPc/jlxXhRDtwDG2Uu/k/OhQz2KNXuCz7n44xUqk88zscOAC4FZ3nwvcWr4WQlScYsQ1oa7HaFHPTtarobgvcfeNZvYosD9wBnBiedo1wO2QqZsihKgMlb9VHIiZzaFYS34PMLN0agBrCDYSM7MFwAKAA6N5BCFExzAm5rj6MbOpwA3AZ9z9BbNXJsDd3c0sOVvu7ouARQDz56XPEUJ0Dk5+D8xOoC7HZWYTKZzWte7+01K81sxmuftqM5sFrBupTgoh2olVfx2XFUOrK4FH3f3iAYcWA+cCXyv//nwoW5sn7M5v9955a+gpJ8UJu6+d/VRS7hvibFKbFRx4b3rA9ySRQlxS+KCMn35y77S9/V5Mh/wBdt8j3bdngmUH9/DO0NYjzE/Ku5kT6jwTrDuJvsCTMlnu0XKIffhDqDN+Qvo//DoOTMotl3wc4H8eH7MrggM9GXv/ENgKSjfbmbGtMDE79wtNFftrRZI1xrYWRRXN7Crg3cA6d99pQZKZnUjhO35fin7q7n87lN163OpxwDnAQ2Yvr7j6IoXD+rGZfQx4CvhgHbaEEB1Oi+e4rgYupVhSFXGnu7+7EaP1RBXvIvbjJzfSmBCi82llrqK7/7oM6rWUzr6RFUKMCg3Mcc0wsyUDXi8qA3KN8FYzewB4Bvicuz88lIIclxBiEA3eKq539/Skan3cDxzk7pvM7HTgZ+ST8AA5LiFEDe1M+XH3FwY8v9nM/o+ZzXD3zC6UbXZcW9g9mZybi1Atn3toUj4pU8Bs+nvT2cxdQfTwHo4Jbc0MooermRfqPMSfpg/sEaowndOT8ud4c1K+JIgcQj4BOqLn+a6kfMae6fdyMi+Ftv4QbC6bK/f83JbpSbmtDzZ+XZPe9BUII3RNlW7OlVu+qEGd3K8tOrZvRid1nS2opOxtXA5hZvsCa8u1oEdTpCFmyhEUaMQlhBhEK6OKZvYjitTAGWb2NHAhxe4TuPt3gPcDnzSzXuAl4Cx3H3KhuhyXEGInWuW43P3sIY5fSrFcoiHkuIQQg1DpZiFE5WjnHFezdHbvhBBtp5UpPyOFHJcQYhBjqqxNK9jC7jzOn+wkX5ethZ5O2B2f2eJ3QpAYPTkIx+eWVmxkWlK+W2YJRzc7J5ID2f9iW4OdpCOd3Bdred9hSfmGNeklBwBsSbdPuuQ9dzx8amwrShjOrcyJvok9DbYBRSHxRmzl7OWSuZc1aCtH1E53Ruf4hOzuJtpOoDkuIUSl0ByXEKJy6FZRCFFJ5LiEEJViB+PCOddOQY5LCLETGnENk5msTcpzb2wUEYkSpqOyxRCXbt6W+Y805+UqtINZnWkninhGUc1cJHbm+PR7ttv+cSR0RpDXuuyGt6QVUmWD++kO5HMyOj2BPEpYzkXuIlu5qGZ0PVMzOpG9qM+5CGUzpZuXJmRxHnvdaI5LCFFJ5LiEEJVCuYpCiMpRpPxocl4IUSE0xyWEqBy6VaxhI9O4nXfsJP9wZsu1x3lNUj49U911ThDWiiKUD/GG0NbaIHoX2YI4qpnLiYyilyuDMshRhBRI5oNCHKEEWLk13U5YOjj3zYkicblIYBTVm9GgHOK+5dqPSjfnoorplNC4nZ6Mra5Anoveptq/L3N+A3R6ys+4oU4ws6vMbJ2ZLRsgW2hmq8xsaflIF0wXQlSO/lvFeh6jxZCOi2In2lQpgG+5+7zycXNruyWEGC2q4Ljq2cl6RHaiFUJ0Ju3cnqxZ6hlxRZxvZg+Wt5J7RSeZ2QIzW2JmS7Y/+/wwmhNCtIP+sjb1PEaLZh3X5cChwDxgNfDN6ER3X+Tu8919/sRXB1XphBAdReVvFVO4+8shNTP7LnBTy3okhBhVxuw6LjOb5e6ry5dnEhexHcSmbdO466kTd5JPOSjODI2WCRzOI6HO0mCX6X/kS0n5r4i2JIYVpHfS/joXhjp/xcUN2QL4JWck5W/jlqQ8Vzp66dojk/IpU+P3edPdr04fiBKDDwhNwb7B7tO9ma/bFkuK/fXp0+2XmfZfmxb7ybGKXREciJZJAH5jYOucQKEnthUmZueWQ6SSvHOJ3HUyJtZxBTvRnmhm8wCnqAXwiZHrohCi3XT6Oq56ooqpnWivHIG+CCE6gB2M0/ZkQojqUflbRSHEroV2+RFCVI4xG1VsNbnC/F38MSmPIoQA7+TnSfkX+XJSviJISoY4yTlqA2ArhyflPWEmbWzvmSCTdwVBuA04fOb9Sfn6TGbypiiq1ZMWTz3+2dhWTzqZ2w+aGOrYU+lIpN0a6OTKMPcEbeRi34EOj8UqdmZwIPhV+b0ZWzvXHijoiXVGKqoIqoAqhKgYOyqQ8iPHJYSoofPnuIaTqyiEGIO0sjpEqixWzXEzs2+b2fIy9/moevooxyWE2IkW5ipeTbosVj+nAXPLxwKKPOgh6ezxoBCi7bQy5aeOslhnAN93dwfuNrOumpTCJHJcQohBtHkd1/7AygGvny5lHeS4thos3zm83XNQV6gSLYc4kX/PNJReXnEfb07Ko7r2AM+9OD0pn7XHM5n206zbGu8+3Tcp/R9uM5OT8kN5OLTV05duZ2OwTAGIk3m70uJN3UFSNsDUYGnDA7EKW4JlD0H74W7ZENd8zyUs9wTyXG37BpeQ2HsytqI+R21A+npaMPlTbE9Wd1RxhpktGfB6kbsvGn4v8mjEJYQYhGP07aj7VnG9u88fRnOrYNBiyQNKWRY5LiHEYBx6e9u2AHUxRTXl64BjgOeHmt8COS4hRA3uRl+udloDBGWxJhbt+HeAm4HTgeXAZuCj9diV4xJCDKJwXC2LKqbKYg087sB5jdqV4xJCDMZpmeMaKdrruPpIRk+iyBnAz0g77FyS8++Zk5T/jjcl5bYqDt3sHpQ7jmwBTO9Lzy1OnxTvvv0Qb2nI1oTx6ZLWAM+N3z8p340XQp0wShdF4jK7Qvvr0xFCm5AJkS1PN+QnpE+3TEnlCP9YfMwuCQ50Z+wFSyXD6GFmV+woAdvSVbgLRiqq6OPYtiUufNAJaMQlhBiMAxpxCSEqhZsclxCiYjjQm951qVOQ4xJC7EyLChKOFHJcQojB7CCfatQB1LOv4lXAu4F17n5EKdsbuB6YQxF3+aC7p5MKBzKOfL5Yqv0VnpTvPSfOMti8KR2ltDWRRtyp6PPbZ88/hDobnp6Vlveko30As9/0RFpn2dy0PLwWsCh61fuqWKkn0olVwvajN21L5sMPSjGHG792ZzoQvDf2hYxOdJ2ZEtH2geBAdyDP/Nos/THn6UnIWjFSciDY07dTqCd4ejU719O5ALjV3ecCt5avhRBjAadYulTPY5QY0nG5+6+BDTXiM4BryufXAO9rbbeEEKNKb52PUaLZOa6ZAxIh1wBxvRYhRLVwxv7kvLu7maUnogAzW0BRkhX2OXC4zQkhRpoKOK5mEwTWmtksgPLvuuhEd1/k7vPdfT57ZorPCSE6g/6oYj2PUaJZx7UYOLd8fi5kEgeFENWj6nNcQT2drwE/NrOPAU8BH6yrtW0UFaVreGLFG2OdYAXvhuXx0oIoHM4BwR3tlswq4SDZ9Nmn0kseAHg6KEM8J/4Xtfa5YJpweaDQEzfPAYE892lH9qIvZ85WVyCPrgXizyxqP5OwHI4EcktxomUPPU3oRH3O2WqmzyNFBZZDDOm4MvV0Tm5xX4QQnUD/cogORivnhRCDqcDkvByXEGIwYyHlRwixC6IRlxCiUuhWsYaXgGUJ+fFxVM+DvVp3ey4uQ7yddDKxH5pux24ITcFhQd/mfTXWuXthWr4mDhFte3P6mEURwq64eQ/CJnZtrBNG6aJoX65087uaaD+4NfHPB7ZyCdM9gfzqhbHO8cGxzC/E/zMtt7cGCpmEbX82sJVb+tiTkLViUl2OSwhROcbCcgghxC6GlkMIISqHo6iiEKJiaI5LCFE5NMclhKgcmuOqIao5f1esYlFi7vpM/fTg/tyWBudnQvtEOrcvjHWiXZYz8wYWvQfRcoQZGVvdwYHcdUZ9eyyQR8s0APufwYEmEqPtU8H5mZr74bf6vy6MdboDeWbHbEtvPp7vW2RrWnAgl2Tdk5C1yuHoVlEIUSk0xyWEqBzKVRRCVA6NuIQQlUSOSwhRKbQcooagdDNHZHSWNNFOFHFrpgxxdK9/e929qa+dXMQvRU/mWNTnnM7SQB71a9+MrVaWLo4Sk3PvVzOfc2SvJ6MTEfW5mXmj3HWmIo47mmijlhYvhzCzU4FLgPHA99z9azXHPwL8A7CqFF3q7t/L2dSISwgxmBZOzpvZeOAy4BSKYcu9ZrbY3R+pOfV6dz+/XrtyXEKIwbT2VvFoYLm7PwlgZtcBZwC1jqshmt2eTAgxlumr81Hs/rVkwGNBjaX9gZUDXj9dymr5czN70Mx+Ymazh+qeRlxCiME0thxivbvPH2aLvwB+5O5bzewTwDXASTkFjbiEEIPpd1yt2RB2FTBwBHUAr0zCF825P+fuW8uX3wPePJTRYY24zKwb2EgxaOxtgecVQow2rZ3juheYa2YHUziss4APDTzBzGa5++ry5XuBR4cy2opbxXe4e6aa9gC2kU5mTclK/MK03K7MtBMk8/pfBLYuydgKllZE/QKwvw4OZJKM/e8DW58MFDLJt/6VwNY5sU74TYjkPZn2/zloP3pfIP7MLg9sRXXdIVyq4TfGKva64MCcWMfvDWztFShkInUebLJuFhwAWJPYD6EVC0cd2DrkWfWZcu81s/OBX1Ish7jK3R82s78Flrj7YuB/mNl7KXq/AfjIUHY1xyWEGEyLU37c/Wbg5hrZlwc8/wKQ2/5kJ4Y7x+XAf5jZfYloghCiivTfKtbzGCWGO+I63t1Xmdk+wC1m9pi7/3rgCaVDK5zapAOH2ZwQYsSpQCHBYY243H1V+XcdcCPFYrPacxa5+3x3n8/E3CZxQoiOoLVRxRGhacdlZnuYFXUbzWwP4M9Ib/cqhKgaHe64zKNwxlCKZodQjLKguOX8obtntncGmzLfeU0iazq3iKIrkDeafApxSd1cXtZrA3nuQ4vaicpQQ3ydUTJzM6WLc32OShRH72WmdHR4LJeYHbUf9Tn3XkbfjWaS6XPfjSiWHr7PmzPGJjYoJx2J3Twf71sSbw1fBzZhvvOqOqsb/NHuG41lUE3PcZW5R29qYV+EEJ2ACgkKISqHHJcQonKokKAQopJ0+HIIOS4hxM40F7NrG+11XL2kIzGZTTeje+0oHw6GyMlL0ZM5Fm2IOqfBNiC/uWdEd1rst8Uq9rYm2on61nDkLD7mP4hV7D3BgShCmIv2Rcea6HNTJaID3KeEx8yiiOPa2OCmzK68YxyVtRFCVA7dKgohauj82Xk5LiFEDZ2/HkKOSwhRg0ZcQojKsQN4abQ7kUWOSwhRg0Zc9ZFLGA6O2TsyOlEZ4F8Etg7O2DossPWtWMWODA5k3u1oeYOdEsijUsMQL23ITVs0upN2Ewnbdlrj9sL3JfeZBbxc1TxlL0pLzpXIDsstR/LM0obwy7F3Rie1hKIVW1mD5riEEBVDIy4hROVQVFEIUTk04hJCVA5FFYUQlUO3ioMZT7pEca6kby6ZNiKIBNmZwflReWYII5ShLYg/88x3IbQX6XQ10X5Op9Fk6lzCeBSh7M7oBN9Em9tgGxBep83K6DRBvFnruiasRW90d0Yn9eYMq2pziW4VhRCVQyMuIUTl0IhLCFE5HE3OCyEqhkZcQojK0flzXMOqgGpmp5rZ42a23MwuaFWnhBCjSf+Iq57H6ND0iMvMxgOXAadQpNXea2aL3f2RUGkv4P0JeVTXHeJ68HMyOlFoP1p2kXsXon88URsAxzahEy0vmBfIc4npEbl/oqnPBeJofG6ZSjO19aPPIHrPuppoI0fU52A5THEsWHqwfmZanv2eRUsr9skopdpv1XKIzh5xDedW8WhgebmjNWZ2HXAGEDsuIUQFGNtzXPsDKwe8fho4ZnjdEUKMPkr5wcwWAAsA2PPAkW5OCDFsxvat4ipg9oDXB5SyQbj7ImARgO03v8O3mRRCjPVbxXuBuWZ2MIXDOgv4UEt6JYQYRTp/xGUe1Z+tR9nsdOCfKNKnr3L3rw5x/rPAU+XLGeTjbCON2lf7Y7H9g9z91cMxYGb/TtG/eljv7qcOp71mGJbjGlbDZkvcff6oNK721f4u3n7VGdYCVCGEGA3kuIQQlWM0HdeiUWxb7av9Xb39SjNqc1xCCNEsulUUQlQOOS4hROUYFcc12uVwzKzbzB4ys6VmtqQN7V1lZuvMbNkA2d5mdouZPVH+3avN7S80s1Xle7C0XJM3Em3PNrPbzOwRM3vYzD5dytty/Zn223X9u5vZb8zsgbL9i0r5wWZ2T/kbuN7MdhuJ9scs7t7WB8Vi1RXAIcBuwAPA4W3uQzcwo43tvR04Clg2QPYN4ILy+QXA19vc/kLgc2249lnAUeXzacDvgMPbdf2Z9tt1/QZMLZ9PBO6hKHz0Y+CsUv4d4JPt+j6OhcdojLheLofj7tuA/nI4YxZ3/zWwoUZ8BnBN+fwa4H1tbr8tuPtqd7+/fL4ReJSiskhbrj/Tflvwgv7N1CaWDwdOAn5Sykf08x+LjIbjSpXDadsXqcSB/zCz+8rqFaPBTHdfXT5fAwTV50aU883swfJWcsRuVfsxsznAkRSjjrZff0370KbrN7PxZraUYsPFWyjuOHrcvT8hcDR+A5VmV52cP97djwJOA84zs7ePZme8uF9o97qUy4FDKWqsrga+OZKNmdlU4AbgM+7+wsBj7bj+RPttu35373P3eRQVVI4mvwWxqIPRcFx1lcMZSdx9Vfl3HXAjxZep3aw1K/ZWLv82s/1x07j72vIHtQP4LiP4HpjZRAqnca27/7QUt+36U+238/r7cfce4DbgrUCXmfVXZ2n7b6DqjIbjerkcThlJOQtY3K7GzWwPM5vW/xz4M2BZXmtEWAycWz4/F/h5OxvvdxolZzJC74GZGXAl8Ki7XzzgUFuuP2q/jdf/ajPrKp9Pptij4VEKB9Zf6b/tn3/lGY2IAHA6RXRnBfA3bW77EIpI5gPAw+1oH/gRxe3Idor5jI8B04FbgSeAXwF7t7n9HwAPAQ9SOJFZI9T28RS3gQ8CS8vH6e26/kz77br+NwK/LdtZBnx5wPfwN8By4F+ASSP9PRxLD6X8CCEqx646OS+EqDByXEKIyiHHJYSoHHJcQojKIcclhKgcclxCiMohxyWEqBz/H6T0F7pe9v53AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Load the last statepoint file and keff value\n",
    "with openmc.StatePoint(statepoint_filename) as sp:\n",
    "    # Get the OpenMC pin power tally data\n",
    "    mesh_tally = sp.get_tally(name='mesh tally')\n",
    "    fission_rates = mesh_tally.get_values(scores=['fission'])\n",
    "\n",
    "# Reshape array to 2D for plotting\n",
    "fission_rates.shape = mesh.dimension\n",
    "\n",
    "# Normalize to the average pin power\n",
    "fission_rates /= np.mean(fission_rates[fission_rates > 0.])\n",
    "\n",
    "# Force zeros to be NaNs so their values are not included when matplotlib calculates\n",
    "# the color scale\n",
    "fission_rates[fission_rates == 0.] = np.nan\n",
    "\n",
    "# Plot the pin powers and the fluxes\n",
    "plt.figure()\n",
    "plt.imshow(fission_rates, interpolation='none', cmap='jet', origin='lower')\n",
    "plt.colorbar()\n",
    "plt.title('Pin Powers')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There we have it! We have just successfully run the C5G7 benchmark model!"
   ]
  }
 ],
 "metadata": {
  "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}