mgxs-part-iii.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multigroup Cross Section Generation Part III: Libraries\n",
    "This IPython Notebook illustrates the use of the **`openmc.mgxs.Library`** class. The `Library` class is designed to automate the calculation of multi-group cross sections for use cases with one or more domains, cross section types, and/or nuclides. In particular, this Notebook illustrates the following features:\n",
    "\n",
    "* Calculation of multi-group cross sections for a **fuel assembly**\n",
    "* Automated creation, manipulation and storage of `MGXS` with **`openmc.mgxs.Library`**\n",
    "* **Validation** of multi-group cross sections with **[OpenMOC](https://mit-crpg.github.io/OpenMOC/)**\n",
    "* Steady-state pin-by-pin **fission rates comparison** between OpenMC and [OpenMOC](https://mit-crpg.github.io/OpenMOC/)\n",
    "\n",
    "**Note:** This Notebook was created using [OpenMOC](https://mit-crpg.github.io/OpenMOC/) to verify the multi-group cross-sections generated by OpenMC. You must install [OpenMOC](https://mit-crpg.github.io/OpenMOC/) on your system to run this Notebook in its entirety. In addition, this Notebook illustrates the use of [Pandas](https://pandas.pydata.org/) `DataFrames` to containerize multi-group cross section data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Input Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import pickle\n",
    "\n",
    "from IPython.display import Image\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import openmc\n",
    "import openmc.mgxs\n",
    "from openmc.openmoc_compatible import get_openmoc_geometry\n",
    "import openmoc\n",
    "import openmoc.process\n",
    "from openmoc.materialize import load_openmc_mgxs_lib\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a model object to tie together geometry, materials, settings, and tallies\n",
    "model = openmc.Model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we need to define materials that will be used in the problem. We'll create three materials for the fuel, water, and cladding of the fuel pins."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1.6 enriched fuel\n",
    "fuel = openmc.Material(name='1.6% Fuel')\n",
    "fuel.set_density('g/cm3', 10.31341)\n",
    "fuel.add_nuclide('U235', 3.7503e-4)\n",
    "fuel.add_nuclide('U238', 2.2625e-2)\n",
    "fuel.add_nuclide('O16', 4.6007e-2)\n",
    "\n",
    "# borated water\n",
    "water = openmc.Material(name='Borated Water')\n",
    "water.set_density('g/cm3', 0.740582)\n",
    "water.add_nuclide('H1', 4.9457e-2)\n",
    "water.add_nuclide('O16', 2.4732e-2)\n",
    "water.add_nuclide('B10', 8.0042e-6)\n",
    "\n",
    "# zircaloy\n",
    "zircaloy = openmc.Material(name='Zircaloy')\n",
    "zircaloy.set_density('g/cm3', 6.55)\n",
    "zircaloy.add_nuclide('Zr90', 7.2758e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our three materials, we can now create a `Materials` object that can be exported to an actual XML file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate a Materials object\n",
    "model.materials = openmc.Materials([fuel, water, zircaloy])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's move on to the geometry. This problem will be a square array of fuel pins and control rod guide tubes for which we can use OpenMC's lattice/universe feature. The basic universe will have three regions for the fuel, the clad, and the surrounding coolant. The first step is to create the bounding surfaces for fuel and clad, as well as the outer bounding surfaces of the problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create cylinders for the fuel and clad\n",
    "fuel_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, r=0.39218)\n",
    "clad_outer_radius = openmc.ZCylinder(x0=0.0, y0=0.0, r=0.45720)\n",
    "\n",
    "# Create boundary planes to surround the geometry\n",
    "min_x = openmc.XPlane(x0=-10.71, boundary_type='reflective')\n",
    "max_x = openmc.XPlane(x0=+10.71, boundary_type='reflective')\n",
    "min_y = openmc.YPlane(y0=-10.71, boundary_type='reflective')\n",
    "max_y = openmc.YPlane(y0=+10.71, boundary_type='reflective')\n",
    "min_z = openmc.ZPlane(z0=-10., boundary_type='reflective')\n",
    "max_z = openmc.ZPlane(z0=+10., boundary_type='reflective')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the surfaces defined, we can now construct a fuel pin cell from cells that are defined by intersections of half-spaces created by the surfaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Universe to encapsulate a fuel pin\n",
    "fuel_pin_universe = openmc.Universe(name='1.6% Fuel Pin')\n",
    "\n",
    "# Create fuel Cell\n",
    "fuel_cell = openmc.Cell(name='1.6% Fuel')\n",
    "fuel_cell.fill = fuel\n",
    "fuel_cell.region = -fuel_outer_radius\n",
    "fuel_pin_universe.add_cell(fuel_cell)\n",
    "\n",
    "# Create a clad Cell\n",
    "clad_cell = openmc.Cell(name='1.6% Clad')\n",
    "clad_cell.fill = zircaloy\n",
    "clad_cell.region = +fuel_outer_radius & -clad_outer_radius\n",
    "fuel_pin_universe.add_cell(clad_cell)\n",
    "\n",
    "# Create a moderator Cell\n",
    "moderator_cell = openmc.Cell(name='1.6% Moderator')\n",
    "moderator_cell.fill = water\n",
    "moderator_cell.region = +clad_outer_radius\n",
    "fuel_pin_universe.add_cell(moderator_cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Likewise, we can construct a control rod guide tube with the same surfaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a Universe to encapsulate a control rod guide tube\n",
    "guide_tube_universe = openmc.Universe(name='Guide Tube')\n",
    "\n",
    "# Create guide tube Cell\n",
    "guide_tube_cell = openmc.Cell(name='Guide Tube Water')\n",
    "guide_tube_cell.fill = water\n",
    "guide_tube_cell.region = -fuel_outer_radius\n",
    "guide_tube_universe.add_cell(guide_tube_cell)\n",
    "\n",
    "# Create a clad Cell\n",
    "clad_cell = openmc.Cell(name='Guide Clad')\n",
    "clad_cell.fill = zircaloy\n",
    "clad_cell.region = +fuel_outer_radius & -clad_outer_radius\n",
    "guide_tube_universe.add_cell(clad_cell)\n",
    "\n",
    "# Create a moderator Cell\n",
    "moderator_cell = openmc.Cell(name='Guide Tube Moderator')\n",
    "moderator_cell.fill = water\n",
    "moderator_cell.region = +clad_outer_radius\n",
    "guide_tube_universe.add_cell(moderator_cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the pin cell universe, we can construct a 17x17 rectangular lattice with a 1.26 cm pitch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create fuel assembly Lattice\n",
    "assembly = openmc.RectLattice(name='1.6% Fuel Assembly')\n",
    "assembly.pitch = (1.26, 1.26)\n",
    "assembly.lower_left = [-1.26 * 17. / 2.0] * 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we create a NumPy array of fuel pin and guide tube universes for the lattice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create array indices for guide tube locations in lattice\n",
    "template_x = np.array([5, 8, 11, 3, 13, 2, 5, 8, 11, 14, 2, 5, 8,\n",
    "                       11, 14, 2, 5, 8, 11, 14, 3, 13, 5, 8, 11])\n",
    "template_y = np.array([2, 2, 2, 3, 3, 5, 5, 5, 5, 5, 8, 8, 8, 8,\n",
    "                       8, 11, 11, 11, 11, 11, 13, 13, 14, 14, 14])\n",
    "\n",
    "# Initialize an empty 17x17 array of the lattice universes\n",
    "universes = np.empty((17, 17), dtype=openmc.Universe)\n",
    "\n",
    "# Fill the array with the fuel pin and guide tube universes\n",
    "universes[:,:] = fuel_pin_universe\n",
    "universes[template_x, template_y] = guide_tube_universe\n",
    "\n",
    "# Store the array of universes in the lattice\n",
    "assembly.universes = universes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OpenMC requires that there is a \"root\" universe. Let us create a root cell that is filled by the assembly and then assign it to the root universe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create root Cell\n",
    "root_cell = openmc.Cell(name='root cell')\n",
    "root_cell.fill = assembly\n",
    "\n",
    "# Add boundary planes\n",
    "root_cell.region = +min_x & -max_x & +min_y & -max_y & +min_z & -max_z\n",
    "\n",
    "# Create root Universe\n",
    "root_universe = openmc.Universe(universe_id=0, name='root universe')\n",
    "root_universe.add_cell(root_cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now must create a geometry that is assigned a root universe and export it to XML."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Geometry and set root Universe\n",
    "model.geometry = openmc.Geometry(root_universe)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the geometry and materials finished, we now just need to define simulation parameters. In this case, we will use 10 inactive batches and 40 active batches each with 2500 particles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OpenMC simulation parameters\n",
    "batches = 50\n",
    "inactive = 10\n",
    "particles = 10000\n",
    "\n",
    "# Instantiate a Settings object\n",
    "settings = openmc.Settings()\n",
    "settings.batches = batches\n",
    "settings.inactive = inactive\n",
    "settings.particles = particles\n",
    "settings.output = {'tallies': False}\n",
    "\n",
    "# Create an initial uniform spatial source distribution over fissionable zones\n",
    "bounds = [-10.71, -10.71, -10, 10.71, 10.71, 10.]\n",
    "uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)\n",
    "settings.source = openmc.IndependentSource(space=uniform_dist)\n",
    "\n",
    "model.settings = settings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us also create a plot to verify that our fuel assembly geometry was created successfully."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPoAAAD6CAIAAAAHjs1qAAAUJklEQVR4nO2dzZEkuQ2FuR1zU4wXa4ccUMxtDuvIhGSCIsaRPeytQw7IjvViHdChQlXZJAg+EACJbOKdZjJfIb8EWT/JepX9yz//9mdJpc7Q226AVGqdvlT///n7r9f//viNe+2/mnmnyHxrBpE5PoPIHJ/hl+uHmZ+///r+9ft197e//iCrP4riZtC5mKF0ercMOAIDb47QNEOGN8ZaSnn/+r16ujy3g2ZRWZG5NMMmZSBrOgE/nJX5sSWbJmUgMZCyb4z1+YC2LukUqT09KUPPrGEgj9UzmzCIlE271plggC5V8fldHdJwQkwzkFv0RxkyTJfKpvkxzKzMLH5Jm5DVABiKaVrAfgZsYKW5ps1M98e1ArJreMXtpF3HZYQ3LYKCDJx508bTnbm+RoRj8UfBGdojihgYs4hhumk8A44kYlA2bXrv3BHJoyAPfxO5ebO0v625x9BuZ4DbslbAOANJtbhpSoYeld/AtTJvGrfuzixwkmbGWZr1ZtAsZdgCnAytOebA/dJmZp6XKcgT6GHGnXvNERhKNk1hVpYlpnsq9VmVEbHUQcrpnjpIOd1TB4kOAL9//c5fLJfLBcFTQ/OzrMiMMOBmc4Y54Kc2Mmxs2gTDVdMMr0tVMqqBxzJNzMp4qp+5NyTZNMYcsGmDDzN4LHPC3DviSgZl8HUlgwgjm0ZqJgAsEvNwUZZ1TWhJxMAgZdN6DOubtjQADD7QUJosq6HwLGs27amIAeBU6qaCpvt0jipIEHdLwnY6y7qlaUFCyHjT+Af29EpE9hzKRpiEb5XRRVwihi1NaxlMmiZiaJG2Dxz/cGJlRhQNdTKTIh/uxyBKqK5vGg62eOBAtvUDd91YB4DL/68P8C8gPMzXHGkEBsYcoWlt+DZC03YxMOZMRKYOUq7MpA5STvfUQcrpnjpIxHT/+fuv+Bd7ruZkmDAnA2PurswUux/MltmFjgkzz3DdW/3XkIE334ihhBk4q5+HDwLAvep+WVY/BrKCBwNjvhFDr/L6gTNkeH2YEcVxnDKcfgxkK3/89ud0+o1hIGXFoGwaziDqcE940/STB2F4BYAZtyajZxJqU7YY/x6bn8EaBpHWN61tkf6l57kXrKNpGshgfEvUCFnWjcdlGJinXLUrIPyu4+IzDZTxLVGDKEgS876K38C5SegeAF7/3GiPqHlPnxb+mo2/DygZGLV4WwZuXQCYyXCChxSJzNOKwrcmYG2bTBicsqxKBpMnUpCBI481ZHi9upMcTrHMHhnO0KvTYyBfXMmN6xlWNi0CA1nHb6ZxiUjRGytudip7O4bKHIGBN38yhgwApw5SRsRSBymne+og5XRPHaQv1f/xC6k15ggMQ3MytOYIDK2Z+9tMZbSeVcUyeyh42UItivEpWZyhBQbLWpldmzZtNmza9MAxlW2b9sZYCxszqszvir9UfzWLVqxJhl6wrjWDZa3MhWoaafNj8GuaZuCYyiRDawPLjm+JCtZ97gKLkKfXezHAw80iifK0rdmkadKwJ86gTCzjDMVt4JhXjbmmjS9VleFbkyzrcO8aBlya1LTJpHyWAp2fYOAQs3EAuIWYqK9XhAxtJfMsq6viD1zEAPCuHGnA/CrTtICB6vgD5xUA5n82O4Swem5EYMAlalrLYPUEiNC0ZQyIeRwAxktLnxjkdQxpJq+N9HOCZHDKsjo1TRoAdmqa08CZN21y3b01R1i+jbDuzjCQZsa5HfjzDRyRiAz1xdgtGETm+zKIzDEZMgCcOkgZEUsdpJzuqYOU0z11kLgAcJH8glBkvvtPHq0YROYIwBEYROZuALhd8SkWN6okzSbBVyWDFBhcJnNlKNk0OUO5NO3DhxkwlikK9/Vime3DGXMvGqph6FUQMZBFsmk8A1lkTdPGAWCnAJM0y9qaTTDaOiYBYGXTnBicmiYdOL8ZNWSALlWnc2p+54YzkFt6MkwCVqXwRKQfAyNN06zUPjfMB843ALxLAQPA91L8BoYIAFe7guRItyRscYYITYvA0B7XPDX9SkQydTUnb5LhVIaQcfF5vTUMoofnwA0ZCrky4xTLNMlwRmDA87ROd/fsmXPgeF0rDP7y3nA9HzdfnYWNyFZlReYbMdwOOD7D0EwnIpmkskaisn5mXMngzbAYOAPAqYOUEbHUQcrpnjpIOd1TB+lDAHh6ZQa5Xn7ufTzQ/AIfYXgaeIaiW5H4ZAwF6/DEwF2nhOFMg1ZmXKOh7fbeRvLL4d7NHpwYCvUd9UoGUeWes3Ru0biMgTF7MCDm14cZk2goaSbJfnRuDwvKlQHPYzgxiMK3OJgVg3LgnBiQpr0CwD0+5emJlmDxHCkjTYLNKoTsFABuD2rVNBFDi2TCoFmA57O3sgBwpS1ZVo22BJ6qc2cYql1bmrYrE1ZpOjUNyveWqKmUk1SJSF7TOaoIWdYieb00BMZfs/H3ASUDoxYveGqaf2BP41uiKk+bGeYKUZQj9XsiiRicmhaBgZHTwCk/xSEMb9UOxvqUNMNJnga5kazcW3UyYQDPTjR1nBj4Iy5jIDeKBs6EYW7g6ogY/sbaQuNm/rXZyRwBOBkmzIbAmYhMHaTMzKQOUk731EHK6Z46SF/aTc8P+8PFvutlAW5G1hCdzKEYROYIwBEYRObWyf1tpmGGE/yVcZsx5M3V98O4mQGWMkwDR2ialKFnjjBwpHm6aYMAcK+6XzR0pTkZ4jCsAR7fErWVSTSUNDPBOpCBlBVwBAa8aREYojVtfKmqnBOic5vea8iAh295humvxKUTwoThcw/cU763RN0VAA4SPL7KPMvqqvgDFzEAvCtFHSS9fRXTtICB6vgDFzEATG6Z0xoGq2Hmfxw9ZDi2ad4MMwFgk/FQZlmdGMhj9cwiBmY8RJPyqKaZM7xe3UUZzrb6t7/+EMUyybJO5oezBW5rugKTDPo8bQ4cXpYLADt91zU0T+dIIwAnw4R5GXAGgFMHKSNiqYOU0z11kHK6pw5SHQB2utqI8JPHEg84GSbMmrJcALjY3RK1LdvjFjGU5svkxQxSYGXTTBjK7qYpGaTAV+cbYy39mJHyRpW9wIOIgawjYiAxRAy4WRTuc2IonbO+18D1MBCGcQCYzHCSTpHa05MymCSrRFlWDbBITgxOA0ceq2feO3DQpep0Ts3q3DQM5Bb9UaSl8ESkHwOjCE1rnxvmDL4B4F0KGAC+l+I3MEQAuNoVP0e6jCF+0yIwTDcN1Hi6a7KsxSjDKWJQ5kgZs3cae+Lh2TSE4alxABgvLe0vnuFUMlgBg2VJZdP4sqTMm1avuxfJH0+7fn7iV4XLxw9buNmKQQMsZTABzqZ5MBCJyOdlCvIEepjBpxpuToYJcwTg4AwZAE4dpIyIpQ5STvfUQcrpnjpIHwLAz8/471+/8xfL5eNVMGh+GIbXHM/KySCtjAMfwlCZB7dEtYplttt7G0vz/fBiBlHlZGDMvam5kqEyDz7MmGRZSbLeE9QpT9tjwCuDtK4MIjA/BtDco13JUGkmACzS8H3n+l88R8pIAyxiYJCcAsDtQa2aJmJokUyapgzqhAgAgw901ZbAU3XuDEO1a0vTIuTSiiQ1zT+wp1yZSR0k31uixnzNWCP8NRt/H/BTzPdhk2ToVTO3RBUJ/6zGM2zPkbYMDJJfANipaSKGFsmkacqnHMIweHXXR0PJK+7SP7fhXS1tGfDKIK0rgwjMjwE092hXMlQiAsAlzNcrySCtHOornggMlTkTkamDlCszqYOU0z11kHK6pw7Sl3YT89unnhlc8MLNyTBhjgAcnCF/mm3PYAKcTfNgGASAe9WV8VQ/sz4tHAE4AsPtgJGy40RkK1EssyfS7MRgBQyWJZVN48uSMm/a+FJ1ZY7UhAEPXpMMePhWxCBSNk0q8FlkfEvUIAHgCIEnPMsapGkRGKabBsr4lqhBFCSJeV/Fb+DcJHQPAFs9NzSJSPMc6UQphsEvNX2vprUXoCECwCanR0ZDRQwmT6T2XERZVqfUtBOD08CRx+qZ9w7c69W9l+Eke+QUDRUxkHVEDCSGiAE3m+RplQylc9b3GrgeBsJQJyJFv6y5mnFnEHMyxGEQmTVlMwCcOkgZEUsdpJzuqYOU0z11kOoAMH5N4Gd2uojxMydD6xSZlwHXAeDqu1kmK1eaeGoPBS/rZ44AnAzbgd8Ya2FjRpX5/et3PFjXK+tkfjhb4LamKzDJQGL4MZSzB27mlqiiDGdP7elFYCCP1TNLs6wTuzQMn6Bp5gy+t0RlIKRalmUVMLEM06VMJsSzFOgM0jRvBuMA8ByEuSJkaCuZZ1ldFX/gIgaAd+VIA+ZX8URkBMUfOK8AMP/T3SGE1XNjGQNjFjFMzxhRZpA5Sg5cq3EAGC8t7a8owwkykLICjsCANy0CQ7SmcevuzAInaWacRXIPienlW35FVsQwDRyhaVKGnjnCwJHm6aYRicjn5YLTd13I89LJHIpBZI4AHIFBZG6dGQBOHaSMiKUOUk731EHK6Z46SFwAuLCXBe0XYLh5y08eIwAnw4TZEPjDLVFL890sE8tsv8VlzO12ciNZOQhDobq8kqFXWTpwHgykOQhDuQzchw8zolhme7CemSTDz43J0+oZ8OBr+/CenBj4Iy5j0CeWTRjmBm4mACwS/r7DM7RmDRUjkyyrsmkRGBg5DZwyqLM/ALwlWNeeNt5HQ2D8s2m1y4+BUYu3JadZAd8sAJxKOSlEALjaFSSIu+VNBn/NFi1oOCnILwSqczdPTb8SkUxdzQCI3k+ZYJ13EnqCgW/aNAP/8JbBpGkihhbJhEHzlOODlsTKjPKGnYyZPI3eKgEoVwZ8pJ0YmLPTgFkxKAfOiQFp2uAv7w3XVR5mPsBZnczjgWCGcxgNFTE8DTxD+bjagDMUrGn3YihYhycG7jolDGcaw5CJyNRBysxM6iDldE8dpJzuqYNET3enVVhRWT9zMsRhWAxsszKDmKcXGUTmGzHcDjg+w9D8IQBMfjHbC51qYpmisgcy3A74LgyvDzOiOI4yltkz41rP0FbwYxCZc+B4XSu8AsCMW3N6+sGQMmgSbAxDUSRDpTJJAubAXfc+/mGciIwQAG6PGzzLGqFpERja4+JNA+V7S9RdCpLEvK/iN1CViOQ1nUbkc2pOao8oCvdZYUxnWf0YGGmaZqX2AtR84Ma3RHU6bTJPaxK+FamtI2JwapoTg1PTpAPnN6OGDG/VDsb6lD6W2SPrmfGY7nQ0dI6BLJJN4xnIImuaVici8V/WtFczuJl/mXEyB2QQmSMAR2AQmStnBoBTBykjYqmDlNM9dZByuqcO0pd20/PD/nDd6npZ4GS+BYPIfF8GkTkmA/e3mQq7ltRmWXsoVdpzWLbaxfxu14/Bow+kGYzIbgH+fAPnHgCWlu2ds6ay/tT8zBEYSHOEgeuZtQFgPlAG1pXG3/D1VPzcRBLlaVtzhKbxaUSQQaSVA2fetPGlqjLDaZJlHe5dw4BLE741mZTPUqDzEwwcYva9JWqQHGkEmWdZXRV/4CIGgHflSAPmV/FEZATFHzivADD/s9khhNVzIwIDLlHTWgarJ0CEpi1jQMzjALCotOjcyOsY0kxeG+nnBMnglGUVhft60geAnZrmNHBMoHKuaa9Xd2WG89tff+ijoT86f6CnZyYZetHQ1gyWtTIXqmmkzY/Br2magWMqkwytDSzLBYBjfjG2nmFoTobWHIGhNWcAOHWQMiKWOkg53VMHKad76iDVAeBzfvK4haHMXn3uMn8yhskAMJ5T65XtcYsYSvNlMmMW5fhiMhSLpkVgKJKmKWfa1fnGWIvb3T17gQecoVenx9BbAI7AsLJpERjIOn4z7eocB4CdAkzSLKsmfMvIKQCsZHNiMAl+BRk48lhDBuhSdTqn5nduOEORRBsMk4D4J85qlx8DoxZvy8BVwObZW98A8C4FDADfS/EbGCIAXO2KnyNdxsBMINHKwxoFGTjz1PQrEcnU1Zy8SY5UGUIWvaczwBoGkdY3rW2RVfh2TdNAhteru1OWVWT2YyBnfG+1hJQyT2vFoGwaziDqcE940/STB2EY/OW9YnEPico8LCtiqMw8w3UvM8+UDLz5RgwlzMCBM23IQCQiH48Bnz2u5mSYMCcDY84AcOogZUQsdZByuqcOUk731EH6EAC+XtjyV9beZvyafQ0DY47QNOlChx9wBAbGPLglqlMss2fuOQt1hpGBD2Qo1Bf70YBfH2aUSUsTMyny4X4Moljb+qbhYIsHDmRbP3DXja8AcI/PKcvaSpQjNTmikmFL0zThWyuGFmn7wPEPXxcABh/oqi0RTjzLGqFpQVKuEQPAqdRNBU336UjQlhwpeee09RjTWdZs2gQD2LTxLVH9sqwVIs+wZhKIGPjwrQYjmzYhhGHw6q6PZTJmENqVATSLXguyaVJaV+Cr6Fuivn/9Dq7nV8fjzc+yIjPCgJvNGeaAn9rIsLFpEwxXTTNkIjJ1kHJlJnWQcrqnDlJO99RB+tJuel4ZIMtY+E+qJsqamyMwlGyawqwsy90S1eqn2W1cTvPrWp5hC3AytOaYAzcIAPeqK+Oprub//uPf7cawwGEZSil//8+/ljHogZGy41uittJnentmJlgHMhRqmdYKGGcgqRY3TcnQo/IbuFbmTRtfqq7MkU7v5Y9okmWVMkw3TTohTBiUTZveO3dE8ijIw41viRohy7rxuIzwpkVQkIEzb5rxLVGDaEuokBeeiIyggA2sNNc09wCw1VhOM5Bb9EcZMkyXyqb5McwEgE1Oj4yGihj0c4KclKIs6+IX5mzatc4Ew+vVvZfhJCf3t+Yv1ffMorIic9EFX3vdcQJ+OCvzY0s2TcpAYiBl6QDwtURblzTjX86Jyt6OQWSOzyAyx2fIAHDqIGVELHWQ/gcg1yPxJKXdEwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Instantiate a Plot\n",
    "model.export_to_xml()\n",
    "plot = openmc.Plot.from_geometry(model.geometry)\n",
    "plot.pixels = (250, 250)\n",
    "plot.color_by = 'material'\n",
    "plot.to_ipython_image()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see from the plot, we have a nice array of fuel and guide tube pin cells with fuel, cladding, and water!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create an MGXS Library"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to generate multi-group cross sections! First, let's define a 2-group structure using the built-in `EnergyGroups` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate a 2-group EnergyGroups object\n",
    "groups = openmc.mgxs.EnergyGroups(group_edges=[0., 0.625, 20.0e6])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will instantiate an `openmc.mgxs.Library` for the energy groups with the fuel assembly geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize a 2-group MGXS Library for OpenMOC\n",
    "mgxs_lib = openmc.mgxs.Library(model.geometry)\n",
    "mgxs_lib.energy_groups = groups"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we must specify to the `Library` which types of cross sections to compute. In particular, the following are the multi-group cross section `MGXS` subclasses that are mapped to string codes accepted by the `Library` class:\n",
    "\n",
    "* `TotalXS` (`\"total\"`)\n",
    "* `TransportXS` (`\"transport\"` or `\"nu-transport` with `nu` set to `True`)\n",
    "* `AbsorptionXS` (`\"absorption\"`)\n",
    "* `CaptureXS` (`\"capture\"`)\n",
    "* `FissionXS` (`\"fission\"` or `\"nu-fission\"` with `nu` set to `True`)\n",
    "* `KappaFissionXS` (`\"kappa-fission\"`)\n",
    "* `ScatterXS` (`\"scatter\"` or `\"nu-scatter\"` with `nu` set to `True`)\n",
    "* `ScatterMatrixXS` (`\"scatter matrix\"` or `\"nu-scatter matrix\"` with `nu` set to `True`)\n",
    "* `Chi` (`\"chi\"`)\n",
    "* `ChiPrompt` (`\"chi prompt\"`)\n",
    "* `InverseVelocity` (`\"inverse-velocity\"`)\n",
    "* `PromptNuFissionXS` (`\"prompt-nu-fission\"`)\n",
    "* `DelayedNuFissionXS` (`\"delayed-nu-fission\"`)\n",
    "* `ChiDelayed` (`\"chi-delayed\"`)\n",
    "* `Beta` (`\"beta\"`)\n",
    "\n",
    "In this case, let's create the multi-group cross sections needed to run an OpenMOC simulation to verify the accuracy of our cross sections. In particular, we will define `\"nu-transport\"`, `\"nu-fission\"`, `'\"fission\"`, `\"nu-scatter matrix\"` and `\"chi\"` cross sections for our `Library`.\n",
    "\n",
    "**Note**: A variety of different approximate transport-corrected total multi-group cross sections (and corresponding scattering matrices) can be found in the literature. At the present time, the `openmc.mgxs` module only supports the `\"P0\"` transport correction. This correction can be turned on and off through the boolean `Library.correction` property which may take values of `\"P0\"` (default) or `None`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify multi-group cross section types to compute\n",
    "mgxs_lib.mgxs_types = ['nu-transport', 'nu-fission', 'fission', 'nu-scatter matrix', 'chi']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we must specify the type of domain over which we would like the `Library` to compute multi-group cross sections. The domain type corresponds to the type of tally filter to be used in the tallies created to compute multi-group cross sections. At the present time, the `Library` supports `\"material\"`, `\"cell\"`, `\"universe\"`, and `\"mesh\"` domain types. We will use a `\"cell\"` domain type here to compute cross sections in each of the cells in the fuel assembly geometry.\n",
    "\n",
    "**Note:** By default, the `Library` class will instantiate `MGXS` objects for each and every domain (material, cell or universe) in the geometry of interest. However, one may specify a subset of these domains to the `Library.domains` property. In our case, we wish to compute multi-group cross sections in each and every cell since they will be needed in our downstream OpenMOC calculation on the identical combinatorial geometry mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify a \"cell\" domain type for the cross section tally filters\n",
    "mgxs_lib.domain_type = 'cell'\n",
    "\n",
    "# Specify the cell domains over which to compute multi-group cross sections\n",
    "mgxs_lib.domains = model.geometry.get_all_material_cells().values()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can easily instruct the `Library` to compute multi-group cross sections on a nuclide-by-nuclide basis with the boolean `Library.by_nuclide` property. By default, `by_nuclide` is set to `False`, but we will set it to `True` here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute cross sections on a nuclide-by-nuclide basis\n",
    "mgxs_lib.by_nuclide = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, we use the `Library` to construct the tallies needed to compute all of the requested multi-group cross sections in each domain and nuclide."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct all tallies needed for the multi-group cross section library\n",
    "mgxs_lib.build_library()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tallies can now be export to a \"tallies.xml\" input file for OpenMC. \n",
    "\n",
    "**NOTE**: At this point the `Library` has constructed nearly 100 distinct `Tally` objects. The overhead to tally in OpenMC scales as $O(N)$ for $N$ tallies, which can become a bottleneck for large tally datasets. To compensate for this, the Python API's `Tally`, `Filter` and `Tallies` classes allow for the smart *merging* of tallies when possible. The `Library` class supports this runtime optimization with the use of the optional `merge` paramter (`False` by default) for the `Library.add_to_tallies_file(...)` method, as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a \"tallies.xml\" file for the MGXS Library\n",
    "tallies = openmc.Tallies()\n",
    "mgxs_lib.add_to_tallies_file(tallies, merge=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition, we instantiate a fission rate mesh tally to compare with OpenMOC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate a tally Mesh\n",
    "mesh = openmc.RegularMesh(mesh_id=1)\n",
    "mesh.dimension = [17, 17]\n",
    "mesh.lower_left = [-10.71, -10.71]\n",
    "mesh.upper_right = [+10.71, +10.71]\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', 'nu-fission']\n",
    "\n",
    "# Add tally to collection\n",
    "tallies.append(tally)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.tallies = tallies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/mixin.py:70: IDWarning: Another Filter instance already exists with id=126.\n",
      "  warn(msg, IDWarning)\n",
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/mixin.py:70: IDWarning: Another Filter instance already exists with id=21.\n",
      "  warn(msg, IDWarning)\n",
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/mixin.py:70: IDWarning: Another Filter instance already exists with id=2.\n",
      "  warn(msg, IDWarning)\n",
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/mixin.py:70: IDWarning: Another Filter instance already exists with id=3.\n",
      "  warn(msg, IDWarning)\n",
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/mixin.py:70: IDWarning: Another Filter instance already exists with id=4.\n",
      "  warn(msg, IDWarning)\n",
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/mixin.py:70: IDWarning: Another Filter instance already exists with id=96.\n",
      "  warn(msg, IDWarning)\n",
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/mixin.py:70: IDWarning: Another Filter instance already exists with id=15.\n",
      "  warn(msg, IDWarning)\n",
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/mixin.py:70: IDWarning: Another Filter instance already exists with id=114.\n",
      "  warn(msg, IDWarning)\n"
     ]
    },
    {
     "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 23:45:30\n",
      "   MPI Processes | 1\n",
      "  OpenMP Threads | 2\n",
      "\n",
      " Reading settings XML file...\n",
      " Reading cross sections XML file...\n",
      " Reading materials XML file...\n",
      " Reading geometry XML file...\n",
      " Reading U235 from /home/pshriwise/data/xs/openmc/nndc_hdf5/U235.h5\n",
      " Reading U238 from /home/pshriwise/data/xs/openmc/nndc_hdf5/U238.h5\n",
      " Reading O16 from /home/pshriwise/data/xs/openmc/nndc_hdf5/O16.h5\n",
      " Reading H1 from /home/pshriwise/data/xs/openmc/nndc_hdf5/H1.h5\n",
      " Reading B10 from /home/pshriwise/data/xs/openmc/nndc_hdf5/B10.h5\n",
      " Reading Zr90 from /home/pshriwise/data/xs/openmc/nndc_hdf5/Zr90.h5\n",
      " Minimum neutron data temperature: 294 K\n",
      " Maximum neutron data temperature: 294 K\n",
      " Reading tallies XML file...\n",
      " Preparing distributed cell instances...\n",
      " Reading plot XML file...\n",
      " Writing summary.h5 file...\n",
      " Maximum neutron transport energy: 20000000 eV for U235\n",
      " Initializing source particles...\n",
      "\n",
      " ====================>     K EIGENVALUE SIMULATION     <====================\n",
      "\n",
      "  Bat./Gen.      k            Average k\n",
      "  =========   ========   ====================\n",
      "        1/1    1.04638\n",
      "        2/1    1.00498\n",
      "        3/1    1.00535\n",
      "        4/1    1.02695\n",
      "        5/1    1.00781\n",
      "        6/1    1.02035\n",
      "        7/1    1.03808\n",
      "        8/1    1.02532\n",
      "        9/1    1.02783\n",
      "       10/1    1.01371\n",
      "       11/1    1.03205\n",
      "       12/1    1.01947    1.02576 +/- 0.00629\n",
      "       13/1    1.01843    1.02332 +/- 0.00438\n",
      "       14/1    1.03269    1.02566 +/- 0.00388\n",
      "       15/1    1.03879    1.02829 +/- 0.00399\n",
      "       16/1    1.02251    1.02732 +/- 0.00340\n",
      "       17/1    1.01274    1.02524 +/- 0.00355\n",
      "       18/1    1.02454    1.02515 +/- 0.00307\n",
      "       19/1    1.01993    1.02457 +/- 0.00277\n",
      "       20/1    1.00665    1.02278 +/- 0.00306\n",
      "       21/1    1.02200    1.02271 +/- 0.00277\n",
      "       22/1    1.03748    1.02394 +/- 0.00281\n",
      "       23/1    1.02803    1.02425 +/- 0.00260\n",
      "       24/1    1.03052    1.02470 +/- 0.00245\n",
      "       25/1    1.02027    1.02441 +/- 0.00230\n",
      "       26/1    1.02597    1.02450 +/- 0.00216\n",
      "       27/1    1.01776    1.02411 +/- 0.00206\n",
      "       28/1    1.02652    1.02424 +/- 0.00195\n",
      "       29/1    1.01063    1.02353 +/- 0.00198\n",
      "       30/1    1.03300    1.02400 +/- 0.00194\n",
      "       31/1    1.02020    1.02382 +/- 0.00185\n",
      "       32/1    1.03720    1.02443 +/- 0.00187\n",
      "       33/1    1.02797    1.02458 +/- 0.00179\n",
      "       34/1    1.01994    1.02439 +/- 0.00172\n",
      "       35/1    1.02626    1.02446 +/- 0.00166\n",
      "       36/1    1.03183    1.02475 +/- 0.00162\n",
      "       37/1    1.03410    1.02509 +/- 0.00159\n",
      "       38/1    1.00919    1.02452 +/- 0.00164\n",
      "       39/1    1.04388    1.02519 +/- 0.00171\n",
      "       40/1    1.01254    1.02477 +/- 0.00171\n",
      "       41/1    1.01584    1.02448 +/- 0.00168\n",
      "       42/1    1.01002    1.02403 +/- 0.00169\n",
      "       43/1    1.00277    1.02339 +/- 0.00176\n",
      "       44/1    1.00672    1.02289 +/- 0.00177\n",
      "       45/1    1.03974    1.02338 +/- 0.00179\n",
      "       46/1    1.00460    1.02285 +/- 0.00181\n",
      "       47/1    1.03954    1.02331 +/- 0.00182\n",
      "       48/1    1.01414    1.02306 +/- 0.00179\n",
      "       49/1    1.02016    1.02299 +/- 0.00174\n",
      "       50/1    0.99791    1.02236 +/- 0.00181\n",
      " Creating state point statepoint.50.h5...\n",
      "\n",
      " =======================>     TIMING STATISTICS     <=======================\n",
      "\n",
      " Total time for initialization     = 1.2070e-01 seconds\n",
      "   Reading cross sections          = 1.1526e-01 seconds\n",
      " Total time in simulation          = 2.9052e+01 seconds\n",
      "   Time in transport only          = 2.9007e+01 seconds\n",
      "   Time in inactive batches        = 2.2160e+00 seconds\n",
      "   Time in active batches          = 2.6836e+01 seconds\n",
      "   Time synchronizing fission bank = 3.0553e-02 seconds\n",
      "     Sampling source sites         = 1.9077e-02 seconds\n",
      "     SEND/RECV source sites        = 1.1409e-02 seconds\n",
      "   Time accumulating tallies       = 1.0753e-03 seconds\n",
      "   Time writing statepoints        = 6.1566e-03 seconds\n",
      " Total time for finalization       = 1.0271e-05 seconds\n",
      " Total time elapsed                = 2.9182e+01 seconds\n",
      " Calculation Rate (inactive)       = 45126.1 particles/second\n",
      " Calculation Rate (active)         = 14905.4 particles/second\n",
      "\n",
      " ============================>     RESULTS     <============================\n",
      "\n",
      " k-effective (Collision)     = 1.02393 +/- 0.00174\n",
      " k-effective (Track-length)  = 1.02236 +/- 0.00181\n",
      " k-effective (Absorption)    = 1.02412 +/- 0.00167\n",
      " Combined k-effective        = 1.02362 +/- 0.00128\n",
      " Leakage Fraction            = 0.00000 +/- 0.00000\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Run OpenMC\n",
    "statepoint_filename = model.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tally Data Processing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our simulation ran successfully and created statepoint and summary output files. We begin our analysis by instantiating a `StatePoint` object. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the last statepoint file\n",
    "sp = openmc.StatePoint(statepoint_filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The statepoint is now ready to be analyzed by the `Library`. We simply have to load the tallies from the statepoint into the `Library` and our `MGXS` objects will compute the cross sections for us under-the-hood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize MGXS Library with OpenMC statepoint data\n",
    "mgxs_lib.load_from_statepoint(sp)\n",
    "# Retrieve OpenMC's k-effective value\n",
    "openmc_keff = sp.keff.nominal_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Voila! Our multi-group cross sections are now ready to rock 'n roll!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extracting and Storing MGXS Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Library` supports a rich API to automate a variety of tasks, including multi-group cross section data retrieval and storage. We will highlight a few of these features here. First, the `Library.get_mgxs(...)` method allows one to extract an `MGXS` object from the `Library` for a particular domain and cross section type. The following cell illustrates how one may extract the `NuFissionXS` object for the fuel cell.\n",
    "\n",
    "**Note:** The `MGXS.get_mgxs(...)` method will accept either the domain *or* the integer domain ID of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve the NuFissionXS object for the fuel cell from the library\n",
    "fuel_mgxs = mgxs_lib.get_mgxs(fuel_cell, 'nu-fission')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `NuFissionXS` object supports all of the methods described previously in the `openmc.mgxs` tutorials, such as [Pandas](https://pandas.pydata.org/) `DataFrames`:\n",
    "Note that since so few histories were simulated, we should expect a few division-by-error errors as some tallies have not yet scored any results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cell</th>\n",
       "      <th>group in</th>\n",
       "      <th>nuclide</th>\n",
       "      <th>mean</th>\n",
       "      <th>std. dev.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>U235</td>\n",
       "      <td>8.099261e-03</td>\n",
       "      <td>1.626934e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>U238</td>\n",
       "      <td>7.326723e-03</td>\n",
       "      <td>2.168273e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>O16</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>U235</td>\n",
       "      <td>3.613773e-01</td>\n",
       "      <td>1.025247e-03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>U238</td>\n",
       "      <td>6.739270e-07</td>\n",
       "      <td>1.911057e-09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>O16</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   cell  group in nuclide          mean     std. dev.\n",
       "3     1         1    U235  8.099261e-03  1.626934e-05\n",
       "4     1         1    U238  7.326723e-03  2.168273e-05\n",
       "5     1         1     O16  0.000000e+00  0.000000e+00\n",
       "0     1         2    U235  3.613773e-01  1.025247e-03\n",
       "1     1         2    U238  6.739270e-07  1.911057e-09\n",
       "2     1         2     O16  0.000000e+00  0.000000e+00"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = fuel_mgxs.get_pandas_dataframe()\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, we can use the `MGXS.print_xs(...)` method to view a string representation of the multi-group cross section data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Multi-Group XS\n",
      "\tReaction Type  =\tnu-fission\n",
      "\tDomain Type    =\tcell\n",
      "\tDomain ID      =\t1\n",
      "\tNuclide        =\tU235\n",
      "\tCross Sections [cm^-1]:\n",
      "            Group 1 [0.625      - 20000000.0eV]:\t8.10e-03 +/- 2.01e-01%\n",
      "            Group 2 [0.0        - 0.625     eV]:\t3.61e-01 +/- 2.84e-01%\n",
      "\n",
      "\tNuclide        =\tU238\n",
      "\tCross Sections [cm^-1]:\n",
      "            Group 1 [0.625      - 20000000.0eV]:\t7.33e-03 +/- 2.96e-01%\n",
      "            Group 2 [0.0        - 0.625     eV]:\t6.74e-07 +/- 2.84e-01%\n",
      "\n",
      "\tNuclide        =\tO16\n",
      "\tCross Sections [cm^-1]:\n",
      "            Group 1 [0.625      - 20000000.0eV]:\t0.00e+00 +/- 0.00e+00%\n",
      "            Group 2 [0.0        - 0.625     eV]:\t0.00e+00 +/- 0.00e+00%\n",
      "\n",
      "\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/pshriwise/.pyenv/versions/3.9.1/lib/python3.9/site-packages/openmc/tallies.py:1255: RuntimeWarning: invalid value encountered in true_divide\n",
      "  data = self.std_dev[indices] / self.mean[indices]\n"
     ]
    }
   ],
   "source": [
    "fuel_mgxs.print_xs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One can export the entire `Library` to HDF5 with the `Library.build_hdf5_store(...)` method as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Store the cross section data in an \"mgxs/mgxs.h5\" HDF5 binary file\n",
    "mgxs_lib.build_hdf5_store(filename='mgxs.h5', directory='mgxs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The HDF5 store will contain the numerical multi-group cross section data indexed by domain, nuclide and cross section type. Some data workflows may be optimized by storing and retrieving binary representations of the `MGXS` objects in the `Library`. This feature is supported through the `Library.dump_to_file(...)` and `Library.load_from_file(...)` routines which use Python's [`pickle`](https://docs.python.org/3/library/pickle.html) module. This is illustrated as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Store a Library and its MGXS objects in a pickled binary file \"mgxs/mgxs.pkl\"\n",
    "mgxs_lib.dump_to_file(filename='mgxs', directory='mgxs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate a new MGXS Library from the pickled binary file \"mgxs/mgxs.pkl\"\n",
    "mgxs_lib = openmc.mgxs.Library.load_from_file(filename='mgxs', directory='mgxs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Library` class may be used to leverage the energy condensation features supported by the `MGXS` class. In particular, one can use the `Library.get_condensed_library(...)` with a coarse group structure which is a subset of the original \"fine\" group structure as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a 1-group structure\n",
    "coarse_groups = openmc.mgxs.EnergyGroups(group_edges=[0., 20.0e6])\n",
    "\n",
    "# Create a new MGXS Library on the coarse 1-group structure\n",
    "coarse_mgxs_lib = mgxs_lib.get_condensed_library(coarse_groups)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>cell</th>\n",
       "      <th>group in</th>\n",
       "      <th>nuclide</th>\n",
       "      <th>mean</th>\n",
       "      <th>std. dev.</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>U235</td>\n",
       "      <td>0.074479</td>\n",
       "      <td>0.000151</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>U238</td>\n",
       "      <td>0.005950</td>\n",
       "      <td>0.000017</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>O16</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   cell  group in nuclide      mean  std. dev.\n",
       "0     1         1    U235  0.074479   0.000151\n",
       "1     1         1    U238  0.005950   0.000017\n",
       "2     1         1     O16  0.000000   0.000000"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Retrieve the NuFissionXS object for the fuel cell from the 1-group library\n",
    "coarse_fuel_mgxs = coarse_mgxs_lib.get_mgxs(fuel_cell, 'nu-fission')\n",
    "\n",
    "# Show the Pandas DataFrame for the 1-group MGXS\n",
    "coarse_fuel_mgxs.get_pandas_dataframe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Verification with OpenMOC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course it is always a good idea to verify that one's cross sections are accurate. We can easily do so here with the deterministic transport code [OpenMOC](https://mit-crpg.github.io/OpenMOC/). We first construct an equivalent OpenMOC geometry."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an OpenMOC Geometry from the OpenMC Geometry\n",
    "openmoc_geometry = get_openmoc_geometry(mgxs_lib.geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can inject the multi-group cross sections into the equivalent fuel assembly OpenMOC geometry. The `openmoc.materialize` module supports the loading of `Library` objects from OpenMC as illustrated below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ WARNING ]  Group cross sections by nuclides are not currently supported.\n",
      "[ WARNING ]  ...  Contributions from all nuclides will be summed.\n"
     ]
    }
   ],
   "source": [
    "# Load the library into the OpenMOC geometry\n",
    "materials = load_openmc_mgxs_lib(mgxs_lib, openmoc_geometry)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to run OpenMOC to verify our cross-sections from OpenMC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  NORMAL ]  Initializing a default angular quadrature...\n",
      "[  NORMAL ]  Initializing 2D tracks...\n",
      "[  NORMAL ]  Initializing 2D tracks reflections...\n",
      "[  NORMAL ]  Initializing 2D tracks array...\n",
      "[  NORMAL ]  Ray tracing for 2D track segmentation...\n",
      "[ WARNING ]  The Geometry was set with non-infinite z-boundaries and supplied\n",
      "[ WARNING ]  ...  to a 2D TrackGenerator. The min-z boundary was set to -10.00 \n",
      "[ WARNING ]  ... and the max-z boundary was set to 10.00. Z-boundaries are \n",
      "[ WARNING ]  ... assumed to be infinite in 2D TrackGenerators.\n",
      "[  NORMAL ]  Progress Segmenting 2D tracks: 0.02 %\n",
      "[  NORMAL ]  Progress Segmenting 2D tracks: 10.02 %\n",
      "[  NORMAL ]  Progress Segmenting 2D tracks: 20.02 %\n",
      "[  NORMAL ]  Progress Segmenting 2D tracks: 30.02 %\n",
      "[  NORMAL ]  Progress Segmenting 2D tracks: 40.02 %\n",
      "[  NORMAL ]  Progress Segmenting 2D tracks: 50.02 %\n",
      "[  NORMAL ]  Progress Segmenting 2D tracks: 60.02 %\n",
      "[  NORMAL ]  Progress Segmenting 2D tracks: 70.02 %\n",
      "[  NORMAL ]  Progress Segmenting 2D tracks: 80.02 %\n",
      "[  NORMAL ]  Progress Segmenting 2D tracks: 90.02 %\n",
      "[  NORMAL ]  Progress Segmenting 2D tracks: 100.00 %\n",
      "[  NORMAL ]  Initializing FSR lookup vectors\n",
      "[  NORMAL ]  Total number of FSRs 867\n",
      "[  NORMAL ]  Initializing MOC eigenvalue solver...\n",
      "[  NORMAL ]  Initializing solver arrays...\n",
      "[  NORMAL ]  Centering segments around FSR centroid...\n",
      "[  NORMAL ]  Max boundary angular flux storage per domain =   0.42 MB\n",
      "[  NORMAL ]  Max scalar flux storage per domain =   0.01 MB\n",
      "[  NORMAL ]  Max source storage per domain =   0.01 MB\n",
      "[  NORMAL ]  Number of azimuthal angles = 32\n",
      "[  NORMAL ]  Azimuthal ray spacing = 0.100000\n",
      "[  NORMAL ]  Number of polar angles = 6\n",
      "[  NORMAL ]  Source type = Flat\n",
      "[  NORMAL ]  MOC transport undamped\n",
      "[  NORMAL ]  CMFD acceleration: OFF\n",
      "[  NORMAL ]  Using 1 threads\n",
      "[  NORMAL ]  Computing the eigenvalue...\n",
      "[  NORMAL ]  Iteration 0:  k_eff = 0.823342   res = 9.831E-02  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  -17665 D.R. = 0.0983\n",
      "[  NORMAL ]  Iteration 1:  k_eff = 0.780160   res = 4.646E-02  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  -4318 D.R. = 0.4725\n",
      "[  NORMAL ]  Iteration 2:  k_eff = 0.739336   res = 9.631E-03  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  -4082 D.R. = 0.2073\n",
      "[  NORMAL ]  Iteration 3:  k_eff = 0.710750   res = 8.566E-03  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  -2858 D.R. = 0.8894\n",
      "[  NORMAL ]  Iteration 4:  k_eff = 0.689598   res = 5.205E-03  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  -2115 D.R. = 0.6076\n",
      "[  NORMAL ]  Iteration 5:  k_eff = 0.675031   res = 3.605E-03  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  -1456 D.R. = 0.6926\n",
      "[  NORMAL ]  Iteration 6:  k_eff = 0.665895   res = 2.538E-03  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  -913 D.R. = 0.7040\n",
      "[  NORMAL ]  Iteration 7:  k_eff = 0.661318   res = 1.889E-03  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  -457 D.R. = 0.7443\n",
      "[  NORMAL ]  Iteration 8:  k_eff = 0.660529   res = 1.493E-03  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  -78 D.R. = 0.7904\n",
      "[  NORMAL ]  Iteration 9:  k_eff = 0.662872   res = 1.268E-03  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  234 D.R. = 0.8490\n",
      "[  NORMAL ]  Iteration 10:  k_eff = 0.667780   res = 1.140E-03  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  490 D.R. = 0.8995\n",
      "[  NORMAL ]  Iteration 11:  k_eff = 0.674766   res = 1.065E-03  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  698 D.R. = 0.9337\n",
      "[  NORMAL ]  Iteration 12:  k_eff = 0.683410   res = 1.014E-03  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  864 D.R. = 0.9527\n",
      "[  NORMAL ]  Iteration 13:  k_eff = 0.693354   res = 9.759E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  994 D.R. = 0.9623\n",
      "[  NORMAL ]  Iteration 14:  k_eff = 0.704290   res = 9.438E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1093 D.R. = 0.9671\n",
      "[  NORMAL ]  Iteration 15:  k_eff = 0.715957   res = 9.154E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1166 D.R. = 0.9698\n",
      "[  NORMAL ]  Iteration 16:  k_eff = 0.728134   res = 8.892E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1217 D.R. = 0.9714\n",
      "[  NORMAL ]  Iteration 17:  k_eff = 0.740633   res = 8.644E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1249 D.R. = 0.9721\n",
      "[  NORMAL ]  Iteration 18:  k_eff = 0.753297   res = 8.404E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1266 D.R. = 0.9722\n",
      "[  NORMAL ]  Iteration 19:  k_eff = 0.765995   res = 8.167E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1269 D.R. = 0.9719\n",
      "[  NORMAL ]  Iteration 20:  k_eff = 0.778619   res = 7.930E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1262 D.R. = 0.9710\n",
      "[  NORMAL ]  Iteration 21:  k_eff = 0.791079   res = 7.691E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1246 D.R. = 0.9698\n",
      "[  NORMAL ]  Iteration 22:  k_eff = 0.803304   res = 7.447E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1222 D.R. = 0.9683\n",
      "[  NORMAL ]  Iteration 23:  k_eff = 0.815235   res = 7.199E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1193 D.R. = 0.9667\n",
      "[  NORMAL ]  Iteration 24:  k_eff = 0.826828   res = 6.947E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1159 D.R. = 0.9650\n",
      "[  NORMAL ]  Iteration 25:  k_eff = 0.838047   res = 6.693E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1121 D.R. = 0.9634\n",
      "[  NORMAL ]  Iteration 26:  k_eff = 0.848868   res = 6.436E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1082 D.R. = 0.9617\n",
      "[  NORMAL ]  Iteration 27:  k_eff = 0.859271   res = 6.179E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  1040 D.R. = 0.9600\n",
      "[  NORMAL ]  Iteration 28:  k_eff = 0.869247   res = 5.922E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  997 D.R. = 0.9585\n",
      "[  NORMAL ]  Iteration 29:  k_eff = 0.878788   res = 5.667E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  954 D.R. = 0.9570\n",
      "[  NORMAL ]  Iteration 30:  k_eff = 0.887894   res = 5.415E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  910 D.R. = 0.9556\n",
      "[  NORMAL ]  Iteration 31:  k_eff = 0.896566   res = 5.168E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  867 D.R. = 0.9543\n",
      "[  NORMAL ]  Iteration 32:  k_eff = 0.904810   res = 4.925E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  824 D.R. = 0.9530\n",
      "[  NORMAL ]  Iteration 33:  k_eff = 0.912635   res = 4.688E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  782 D.R. = 0.9519\n",
      "[  NORMAL ]  Iteration 34:  k_eff = 0.920049   res = 4.457E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  741 D.R. = 0.9508\n",
      "[  NORMAL ]  Iteration 35:  k_eff = 0.927066   res = 4.233E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  701 D.R. = 0.9498\n",
      "[  NORMAL ]  Iteration 36:  k_eff = 0.933696   res = 4.016E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  663 D.R. = 0.9488\n",
      "[  NORMAL ]  Iteration 37:  k_eff = 0.939954   res = 3.807E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  625 D.R. = 0.9479\n",
      "[  NORMAL ]  Iteration 38:  k_eff = 0.945855   res = 3.606E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  590 D.R. = 0.9471\n",
      "[  NORMAL ]  Iteration 39:  k_eff = 0.951412   res = 3.413E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  555 D.R. = 0.9464\n",
      "[  NORMAL ]  Iteration 40:  k_eff = 0.956641   res = 3.227E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  522 D.R. = 0.9456\n",
      "[  NORMAL ]  Iteration 41:  k_eff = 0.961556   res = 3.049E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  491 D.R. = 0.9448\n",
      "[  NORMAL ]  Iteration 42:  k_eff = 0.966174   res = 2.879E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  461 D.R. = 0.9443\n",
      "[  NORMAL ]  Iteration 43:  k_eff = 0.970507   res = 2.716E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  433 D.R. = 0.9435\n",
      "[  NORMAL ]  Iteration 44:  k_eff = 0.974571   res = 2.561E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  406 D.R. = 0.9430\n",
      "[  NORMAL ]  Iteration 45:  k_eff = 0.978380   res = 2.414E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  380 D.R. = 0.9423\n",
      "[  NORMAL ]  Iteration 46:  k_eff = 0.981948   res = 2.273E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  356 D.R. = 0.9417\n",
      "[  NORMAL ]  Iteration 47:  k_eff = 0.985287   res = 2.140E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  333 D.R. = 0.9413\n",
      "[  NORMAL ]  Iteration 48:  k_eff = 0.988410   res = 2.013E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  312 D.R. = 0.9409\n",
      "[  NORMAL ]  Iteration 49:  k_eff = 0.991331   res = 1.893E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  292 D.R. = 0.9402\n",
      "[  NORMAL ]  Iteration 50:  k_eff = 0.994060   res = 1.779E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  272 D.R. = 0.9399\n",
      "[  NORMAL ]  Iteration 51:  k_eff = 0.996609   res = 1.671E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  254 D.R. = 0.9393\n",
      "[  NORMAL ]  Iteration 52:  k_eff = 0.998988   res = 1.569E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  237 D.R. = 0.9390\n",
      "[  NORMAL ]  Iteration 53:  k_eff = 1.001209   res = 1.473E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  222 D.R. = 0.9386\n",
      "[  NORMAL ]  Iteration 54:  k_eff = 1.003280   res = 1.382E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  207 D.R. = 0.9380\n",
      "[  NORMAL ]  Iteration 55:  k_eff = 1.005211   res = 1.296E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  193 D.R. = 0.9378\n",
      "[  NORMAL ]  Iteration 56:  k_eff = 1.007012   res = 1.215E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  180 D.R. = 0.9374\n",
      "[  NORMAL ]  Iteration 57:  k_eff = 1.008689   res = 1.138E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  167 D.R. = 0.9369\n",
      "[  NORMAL ]  Iteration 58:  k_eff = 1.010251   res = 1.066E-04  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  156 D.R. = 0.9369\n",
      "[  NORMAL ]  Iteration 59:  k_eff = 1.011706   res = 9.981E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  145 D.R. = 0.9362\n",
      "[  NORMAL ]  Iteration 60:  k_eff = 1.013060   res = 9.344E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  135 D.R. = 0.9361\n",
      "[  NORMAL ]  Iteration 61:  k_eff = 1.014320   res = 8.743E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  125 D.R. = 0.9357\n",
      "[  NORMAL ]  Iteration 62:  k_eff = 1.015492   res = 8.177E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  117 D.R. = 0.9353\n",
      "[  NORMAL ]  Iteration 63:  k_eff = 1.016582   res = 7.648E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  109 D.R. = 0.9353\n",
      "[  NORMAL ]  Iteration 64:  k_eff = 1.017596   res = 7.147E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  101 D.R. = 0.9345\n",
      "[  NORMAL ]  Iteration 65:  k_eff = 1.018538   res = 6.680E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  94 D.R. = 0.9346\n",
      "[  NORMAL ]  Iteration 66:  k_eff = 1.019414   res = 6.237E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  87 D.R. = 0.9338\n",
      "[  NORMAL ]  Iteration 67:  k_eff = 1.020227   res = 5.828E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  81 D.R. = 0.9343\n",
      "[  NORMAL ]  Iteration 68:  k_eff = 1.020983   res = 5.442E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  75 D.R. = 0.9338\n",
      "[  NORMAL ]  Iteration 69:  k_eff = 1.021685   res = 5.080E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  70 D.R. = 0.9336\n",
      "[  NORMAL ]  Iteration 70:  k_eff = 1.022337   res = 4.739E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  65 D.R. = 0.9328\n",
      "[  NORMAL ]  Iteration 71:  k_eff = 1.022942   res = 4.422E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  60 D.R. = 0.9331\n",
      "[  NORMAL ]  Iteration 72:  k_eff = 1.023504   res = 4.124E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  56 D.R. = 0.9327\n",
      "[  NORMAL ]  Iteration 73:  k_eff = 1.024026   res = 3.843E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  52 D.R. = 0.9317\n",
      "[  NORMAL ]  Iteration 74:  k_eff = 1.024510   res = 3.586E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  48 D.R. = 0.9331\n",
      "[  NORMAL ]  Iteration 75:  k_eff = 1.024959   res = 3.341E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  44 D.R. = 0.9317\n",
      "[  NORMAL ]  Iteration 76:  k_eff = 1.025375   res = 3.115E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  41 D.R. = 0.9325\n",
      "[  NORMAL ]  Iteration 77:  k_eff = 1.025762   res = 2.898E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  38 D.R. = 0.9303\n",
      "[  NORMAL ]  Iteration 78:  k_eff = 1.026120   res = 2.703E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  35 D.R. = 0.9327\n",
      "[  NORMAL ]  Iteration 79:  k_eff = 1.026452   res = 2.519E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  33 D.R. = 0.9318\n",
      "[  NORMAL ]  Iteration 80:  k_eff = 1.026760   res = 2.341E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  30 D.R. = 0.9295\n",
      "[  NORMAL ]  Iteration 81:  k_eff = 1.027046   res = 2.180E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  28 D.R. = 0.9312\n",
      "[  NORMAL ]  Iteration 82:  k_eff = 1.027311   res = 2.028E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  26 D.R. = 0.9301\n",
      "[  NORMAL ]  Iteration 83:  k_eff = 1.027556   res = 1.889E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  24 D.R. = 0.9315\n",
      "[  NORMAL ]  Iteration 84:  k_eff = 1.027783   res = 1.757E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  22 D.R. = 0.9305\n",
      "[  NORMAL ]  Iteration 85:  k_eff = 1.027994   res = 1.632E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  21 D.R. = 0.9284\n",
      "[  NORMAL ]  Iteration 86:  k_eff = 1.028189   res = 1.521E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  19 D.R. = 0.9322\n",
      "[  NORMAL ]  Iteration 87:  k_eff = 1.028370   res = 1.412E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  18 D.R. = 0.9285\n",
      "[  NORMAL ]  Iteration 88:  k_eff = 1.028538   res = 1.311E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  16 D.R. = 0.9287\n",
      "[  NORMAL ]  Iteration 89:  k_eff = 1.028693   res = 1.222E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  15 D.R. = 0.9316\n",
      "[  NORMAL ]  Iteration 90:  k_eff = 1.028837   res = 1.134E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  14 D.R. = 0.9279\n",
      "[  NORMAL ]  Iteration 91:  k_eff = 1.028970   res = 1.056E-05  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  13 D.R. = 0.9314\n",
      "[  NORMAL ]  Iteration 92:  k_eff = 1.029093   res = 9.786E-06  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  12 D.R. = 0.9268\n",
      "[  NORMAL ]  Iteration 93:  k_eff = 1.029207   res = 9.093E-06  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  11 D.R. = 0.9292\n",
      "[  NORMAL ]  Iteration 94:  k_eff = 1.029313   res = 8.445E-06  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  10 D.R. = 0.9287\n",
      "[  NORMAL ]  Iteration 95:  k_eff = 1.029411   res = 7.848E-06  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  9 D.R. = 0.9294\n",
      "[  NORMAL ]  Iteration 96:  k_eff = 1.029502   res = 7.272E-06  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  9 D.R. = 0.9265\n",
      "[  NORMAL ]  Iteration 97:  k_eff = 1.029586   res = 6.769E-06  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  8 D.R. = 0.9309\n",
      "[  NORMAL ]  Iteration 98:  k_eff = 1.029663   res = 6.294E-06  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  7 D.R. = 0.9298\n",
      "[  NORMAL ]  Iteration 99:  k_eff = 1.029735   res = 5.827E-06  delta-k (pcm) =\n",
      "[  NORMAL ]  ...  7 D.R. = 0.9259\n",
      "[  NORMAL ]  Iteration 100:  k_eff = 1.029802   res = 5.413E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 6 D.R. = 0.9290\n",
      "[  NORMAL ]  Iteration 101:  k_eff = 1.029863   res = 5.032E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 6 D.R. = 0.9295\n",
      "[  NORMAL ]  Iteration 102:  k_eff = 1.029920   res = 4.648E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 5 D.R. = 0.9237\n",
      "[  NORMAL ]  Iteration 103:  k_eff = 1.029973   res = 4.328E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 5 D.R. = 0.9313\n",
      "[  NORMAL ]  Iteration 104:  k_eff = 1.030022   res = 4.004E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 4 D.R. = 0.9249\n",
      "[  NORMAL ]  Iteration 105:  k_eff = 1.030067   res = 3.734E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 4 D.R. = 0.9326\n",
      "[  NORMAL ]  Iteration 106:  k_eff = 1.030109   res = 3.452E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 4 D.R. = 0.9246\n",
      "[  NORMAL ]  Iteration 107:  k_eff = 1.030148   res = 3.195E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 3 D.R. = 0.9253\n",
      "[  NORMAL ]  Iteration 108:  k_eff = 1.030184   res = 2.979E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 3 D.R. = 0.9325\n",
      "[  NORMAL ]  Iteration 109:  k_eff = 1.030217   res = 2.750E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 3 D.R. = 0.9230\n",
      "[  NORMAL ]  Iteration 110:  k_eff = 1.030247   res = 2.541E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 3 D.R. = 0.9240\n",
      "[  NORMAL ]  Iteration 111:  k_eff = 1.030276   res = 2.364E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 2 D.R. = 0.9302\n",
      "[  NORMAL ]  Iteration 112:  k_eff = 1.030302   res = 2.178E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 2 D.R. = 0.9215\n",
      "[  NORMAL ]  Iteration 113:  k_eff = 1.030326   res = 2.022E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 2 D.R. = 0.9285\n",
      "[  NORMAL ]  Iteration 114:  k_eff = 1.030349   res = 1.896E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 2 D.R. = 0.9374\n",
      "[  NORMAL ]  Iteration 115:  k_eff = 1.030369   res = 1.753E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 2 D.R. = 0.9246\n",
      "[  NORMAL ]  Iteration 116:  k_eff = 1.030389   res = 1.619E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 1 D.R. = 0.9236\n",
      "[  NORMAL ]  Iteration 117:  k_eff = 1.030406   res = 1.500E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 1 D.R. = 0.9266\n",
      "[  NORMAL ]  Iteration 118:  k_eff = 1.030423   res = 1.406E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 1 D.R. = 0.9375\n",
      "[  NORMAL ]  Iteration 119:  k_eff = 1.030438   res = 1.287E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 1 D.R. = 0.9154\n",
      "[  NORMAL ]  Iteration 120:  k_eff = 1.030452   res = 1.193E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 1 D.R. = 0.9269\n",
      "[  NORMAL ]  Iteration 121:  k_eff = 1.030465   res = 1.095E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 1 D.R. = 0.9174\n",
      "[  NORMAL ]  Iteration 122:  k_eff = 1.030477   res = 1.033E-06  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 1 D.R. = 0.9437\n",
      "[  NORMAL ]  Iteration 123:  k_eff = 1.030488   res = 9.328E-07  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 1 D.R. = 0.9029\n",
      "[  NORMAL ]  Iteration 124:  k_eff = 1.030498   res = 8.716E-07  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 1 D.R. = 0.9344\n",
      "[  NORMAL ]  Iteration 125:  k_eff = 1.030508   res = 8.527E-07  delta-k (pcm)\n",
      "[  NORMAL ]  ...  = 0 D.R. = 0.9783\n"
     ]
    }
   ],
   "source": [
    "# Generate tracks for OpenMOC\n",
    "track_generator = openmoc.TrackGenerator(openmoc_geometry, num_azim=32, azim_spacing=0.1)\n",
    "track_generator.generateTracks()\n",
    "\n",
    "# Run OpenMOC\n",
    "solver = openmoc.CPUSolver(track_generator)\n",
    "solver.computeEigenvalue()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We report the eigenvalues computed by OpenMC and OpenMOC here together to summarize our results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "openmc keff = 1.023625\n",
      "openmoc keff = 1.030508\n",
      "bias [pcm]: 688.3\n"
     ]
    }
   ],
   "source": [
    "# Print report of keff and bias with OpenMC\n",
    "openmoc_keff = solver.getKeff()\n",
    "bias = (openmoc_keff - openmc_keff) * 1e5\n",
    "\n",
    "print('openmc keff = {0:1.6f}'.format(openmc_keff))\n",
    "print('openmoc keff = {0:1.6f}'.format(openmoc_keff))\n",
    "print('bias [pcm]: {0:1.1f}'.format(bias))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a non-trivial bias between the eigenvalues computed by OpenMC and OpenMOC. One can show that these biases do not converge to <100 pcm with more particle histories. For heterogeneous geometries, additional measures must be taken to address the following three sources of bias:\n",
    "\n",
    "* Appropriate transport-corrected cross sections\n",
    "* Spatial discretization of OpenMOC's mesh\n",
    "* Constant-in-angle multi-group cross sections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Flux and Pin Power Visualizations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will conclude this tutorial by illustrating how to visualize the fission rates computed by OpenMOC and OpenMC. First, we extract volume-integrated fission rates from OpenMC's mesh fission rate tally for each pin cell in the fuel assembly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the OpenMC fission rate mesh tally data\n",
    "mesh_tally = sp.get_tally(name='mesh tally')\n",
    "openmc_fission_rates = mesh_tally.get_values(scores=['nu-fission'])\n",
    "\n",
    "# Close the statepoint file now that we're done getting information from it\n",
    "sp.close()\n",
    "\n",
    "# Reshape array to 2D for plotting\n",
    "openmc_fission_rates.shape = (17,17)\n",
    "\n",
    "# Normalize to the average pin power\n",
    "openmc_fission_rates /= np.mean(openmc_fission_rates[openmc_fission_rates > 0.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we extract OpenMOC's volume-averaged fission rates into a 2D 17x17 NumPy array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create OpenMOC Mesh on which to tally fission rates\n",
    "openmoc_mesh = openmoc.process.Mesh()\n",
    "openmoc_mesh.dimension = np.array(mesh.dimension)\n",
    "openmoc_mesh.lower_left = np.array(mesh.lower_left)\n",
    "openmoc_mesh.upper_right = np.array(mesh.upper_right)\n",
    "openmoc_mesh.width = openmoc_mesh.upper_right - openmoc_mesh.lower_left\n",
    "openmoc_mesh.width /= openmoc_mesh.dimension\n",
    "\n",
    "# Tally OpenMOC fission rates on the Mesh\n",
    "openmoc_fission_rates = openmoc_mesh.tally_fission_rates(solver)\n",
    "openmoc_fission_rates = np.squeeze(openmoc_fission_rates)\n",
    "openmoc_fission_rates = np.fliplr(openmoc_fission_rates)\n",
    "\n",
    "# Normalize to the average pin fission rate\n",
    "openmoc_fission_rates /= np.mean(openmoc_fission_rates[openmoc_fission_rates > 0.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can easily use Matplotlib to visualize the fission rates from OpenMC and OpenMOC side-by-side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'OpenMOC Fission Rates')"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAADHCAYAAAAXg5iPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjtUlEQVR4nO3deZwcVdX/8c+XkI2sZF8gREhEQmICDouKElDZREMQENyCogFe8qiPIKD+fhLlJ264ISgP8GAQFUExLBohKCYKGkyAxCQIEiCYPYTsIXvO74+q0Z5Jd5/KdM/SlfN+veY1PVWn771VffpOdXXdujIzQggh5Nd+rd2AEEIIzSs6+hBCyLno6EMIIeeiow8hhJyLjj6EEHIuOvoQQsi56OhriKQvSLqtgud/UNK0arYphOYgaYikTZLaVVDGJkmHVrNdtarmO3pJF0qaJ+k1SSsk/UhSzxaqe5Gk7ZL6NFr+tCSTNLRg2bGSpkpaJ2mNpL9J+miJci+UtCtN1PqfG83sOjP7eFPba2Y/M7NTmvr8UiSNlbQ7bedGSc+V2rYSz58uqcnblQc1lMdvkfRo+jqvl/SgpBGNntdd0vck/SvNiRfSvxuUXxBvkjYX5Po6M/uXmXU1s11N3a70+S829fmlpPtrS9rWFZImS+qa8bkXSnqs2m3y1HRHL+ly4BvA54AewPHAIcAjkjq0UDNeAi4oaNMo4IBG7Xwz8CgwAxgG9AYuBU4vU+5f00St/7ms6i2vrmVm1hXoDvw3cKukw1u5TTWhxvJ4GnA/MAh4HTAXeLz+yDlt7x+AI4HTSPLhzcCrwLFl6h9dkOs9q7RNzek9ab6PAY4CPt+6zXGYWU3+kCTQJuC8Rsu7Aq8AH0v/ngT8Crgb2Ag8RZJU9fGDgHvT57wEfKpg3STgHuAn6XMXAHUF6xcB/weYVbDseuCLgAFD02WPATftxbZdCDxWZPkk4Kfp407AT0neQOuAWUD/gue/mLb5JeCDxcoF3pI+b336+y0F66YD1wKPp+VMA/qUaO9YYEmjZauAc9PHBwK/Sffx2vTxQem6rwK7gK3p63ljuvwNwCPAGuC5wtcZOAN4Jm3XUuCK1s7HfSSP/wz8sMg2/A74Sfr448BKoOte7AMDhjVaNjRdvr+T08NIDqDWA6uBu4uVS/IP9Cfp/nk53d79Ct8X6TavTcs/vUx7FwHvLPj7m8BvC/6+GnghbeszwPh0+RFpnu9KX/N16fKOad3/SvfdzUDndF0fkvfLOpL3wp/r271XedbaiV7BG+Q0YGd9IjRadwdwV0GS7wDOAdoDV6QvZHuSTzRPAl8COgCHpsl0asFzt5J0LO2ArwEzG7/gJB3REWnMEpKjMUuT9YD0hT1pL7btQvyO/mLgwbT8dsCbSDqNLsAG4PA0biBwZONygV5pUn8Y2J/kaG4t0DtdPz1N1tcDndO/v16ivWNJO/p0n74X2A0clS7rDbwvbWs34JfAfQXPnw58vODvLsBi4KNp244ieROPSNcvB96WPj4QOLq183FfzuP0dVqePv4FcMde7oOyHb2T03eR/EPaj+Tg54Ri5ZJ08ven+TcU+CdwUcH7YgfwiXTbLwWWASrR3kWkHT1wEDAP+H7B+nNJ/vHuB7wf2AwMLPXeBr4LPEDynuxG8r7+WrruayQdf/v0522l2lXup5ZP3fQBVpvZziLrlqfr6z1pZr8ysx3Ad0gS4njgGKCvmX3FzLZbcj7vVuD8guc+ZmZTLTlXeCcwukh9dwIfAd4F/IPkKLPegSQv+PK93L7j0/P59T/HN1q/g6QDHWZmu8zsSTPbkK7bDYyU1NnMlpvZgiLlvxt43szuNLOdZnYX8CzwnoKYH5vZP81sC8kR4Zgy7R0kaR2wBZgCfNbMngYws1fN7F4ze83MNpIcxZ9YpqwzgUVm9uO0bU+THK2eW7DtIyR1N7O1ZvZUmbLaulrJ416UzuPCdvYuEeN5qiDXbyiyvlRO7yD5hzTIzLaa2R7nv9MvdM8HPm9mG81sEfBtkoOcei+b2a3p/rmD5J9J/zLtvU/SRpIDklXANfUrzOyXZrbMzHab2d3A85Q4bSVJwETgv81sTfr+uI7/vHY70rYcYmY7zOzPlv4H2Bu13NGvBvpI2r/IuoHp+nqL6x+Y2W6So5VBpAlS2KECX6DhC7yi4PFrQKcidd4JfIDkv/VPGq1bS5KkA7Nt1r/NNLOeBT8zi9T5MPALScskfVNSezPbTHIUcQmwXNJvJb2hSPmDSD7CFnoZGFzwd+NtL/eF0zJLzq12B24ATq5fIekASf8j6WVJG4A/AT3LXFFxCHBco9flg8CAdP37SI5OX5Y0Iz13XKvykMeF7Xy1RIzn6IJc/1ThCienrwQE/E3SAkkfK1J2H5Kj4cJ8L5nrZvZa+rBcvp9lZt1IPs2+gYJ/yJI+ImlOwWsxkob/sAv1Jfm09GRB/EPpcoBvAQuBaZJelHR1mTaVVMsd/V+BbcDZhQvTb79PJ/lCqN7BBev3I/m4tYzkjfNSow61m5mdsTcNMbOXST5GnwH8utG619K2vm9vysxQ5w4z+7KZjSA5134mydEYZvawmb2L5A33LMnRXWPLSDqIQkNoeBTXlHZtA64CRkk6K118OXA4cJyZdQfeni5X/dMaFbMYmNHodelqZpemdcwys3FAP+A+kk8btapW8nhz2tZzizz1vIJ2/h44VVKXvak7Q9uK5rSZrTCzT5jZIJLTmT+UNKzR01fznyP/ehXnelr/DGAyyTl2JB2Stu0yktOgPYH5lM711SSfgo8seO16WPJFL+knkMvN7FCSU6KflfSOvW1nzXb0ZrYe+DLwA0mnSWqfXgZ2D8mRzp0F4W+SdHZ6BPMZkjfWTOBvwEZJV0nqLKmdpJGSjmlCky4CTk7fEI1dCVwo6XOSegNIGi3pF02oh/T5J0kalR4VbyBJ5N2S+ksal77RtpF86bO7SBFTgddL+oCk/SW9HxhB8sVPRcxsO8lH4y+li7qRJPM6Sb0o+JibWklyXrneb9K2fTh9XdtLOkbSEZI6KBkP0CM9hbGhxPbVhBrL46uBCZI+JambpAMl/T+Sq2q+nMbcSfKP515Jb5C0n6TeSsaA7NU/nnrlclrSuZIOSkPXknSkDfIhPR1zD/DVtN2HAJ8luZihGr4HvEvSaJLvE4zkS1+UXGY8siB2JXBQenVS/SezW4HvSuqXPmewpFPTx2dKGpae4llP8j3JXud7zXb0AGb2TZKPqNeTvOGfIEmyd6RHlvXuJ/noV//l49npEfEukiPhMSRHMquB20i+od/btrxgZrNLrPsLyamMk4EXJa0BbiHpbJtqAMlVGBtIzqfOIHmT7UeSxMtIvqU/keTLpcZtepVk2y8n+bh9JXCmma1uHNtEtwNDJL2H5I3QmWT/ziT5aFro+8A5ktZKuiE9T3kKyXnKZSQfq79BcnUCJK/hovQ00CUkp3VqVg3l8WPAqSSfPpaTnP44iuQL0OfTmG0kX+w+S3LV1AaSf0R90u1qinI5fQzwhKRNJF9oftqKXzv/XyRfir5IcoXNz0lytGJm9grJqa4vmdkzJAc5fyXp1EeRXLlW71GSq55WSKp/r11FcnpmZprTvyf5BAwwPP17U1rmD83sj3vbRjXhvH5NkTSJ5AvLD7V2W0JoqsjjUImaPqIPIYTgi44+hBByLvenbkIIYV8XR/QhhJBz0dGHEELOFRuN1+rUs48xcGj5oA3lVwO0G7DDjdm1098FPdqvdWPkRsBG6+bGdNYWN6Yd/p1bt9LJjenJOjdmvwx1Zdn65TsyDJZck+G4Y4Uf4o7LXLMI27Q6y0tWVVIvS8Y4lZPh6tYug/yYAX5Ip57FLpVvaNuujm6MbfffQ506+3Vt3eKPsapWOepQ7I4TDXVst82N2bouw7iwLDm7eVmGoFKDa+stwWxN0byuqKOXdBrJNdDtgNvM7OuN1nckub70TSTXar8/vc9EeQOHwuSil/L+x+/99nW9wt/D61f3dGPGDp7ixuyfoUP8/bZ3ujFjOs5xY7qx0Y15jte7Me/lwarUtRN/bojrVn7Rjdn90wxvmuv9EC531n+7zi2ieXL7IPyhExnmlRk1yY/JMFB+2LhZbszC9Ye5MVsX9fLrGu3XNX+BP75r2JEZypnrl9Nx6Bq/rh4v+HXdn2FM2tf9EGZOyhDkTdlQejxak0/dpCMybyIZpj0CuECNJiAgGWW31syGkdyh7RtNrS+ElhK5HfKmknP0xwILzezFdMj7L4BxjWLGkdwJDpJRnO9Ih/KG0JZFbodcqaSjH0zB3fRI7ssxuFSMJbdhXU9yG9M9SJooabak2ax7pYJmhVCxquV2g7zGP10QQnNoM1fdmNktZlZnZnX07Os/IYQa0CCv8c9lh9AcKunol1Jw21SSb5oa3/bz3zHpHfd6kHxxFUJbFrkdcqWSjn4WMFzS69Jbbp5Pcve4Qg8AE9LH5wCPNmV2lBBaWOR2yJUmX15pZjslXUYyy1E74HYzWyDpK8BsM3sA+F/gTkkLSU5Qnl+6xP/Yr/NODhhZ/jz9pnX+6Z1vdfycG/PS4Ne5McNY6MYcl+EOrId19C/X6s9KN+azy37kxiwd5J8muISb3ZiLM8Tcx3g3pltP/zLNjRnuy7h7ToZLMD/jjZ8o3x83X27X3z24HP8SXLb6IQPGFbtTb0MHN/gaorhHe5zsxvxqtD+nzkb8MSTjjzzbjZnCWW5Mt9E3uTHncK8bM+Hf37WXtnpc0a8cG1gx6VA3JtPr7uZO6TEYFV1Hb2ZTaXRhsJl9qeDxVorPSBNCmxa5HfKkzXwZG0IIoXlERx9CCDkXHX0IIeRcdPQhhJBz0dGHEELORUcfQgg5Fx19CCHkXJucM1aH1hnXlr8f/Wc++DW3nLfxZzcmyyQe4745zY3hygz7cax/c8MZ0491Y07MMDiLGX5dd5x4nhszgbv9uh7NsF0n+9u1nQ5uzCl3+6/pEe9/quz6l+o+xJbZz7T8xCNd64xRzjwLGQZDZRgzhJ2WYfOyzHUx3s/ruRnmPhg9+3m/rroM76HZ/nbNrRvut4d/+nVNybAPM8wBo4cybNd9fog7l9C8OmxT8R0UR/QhhJBz0dGHEELORUcfQgg5Fx19CCHkXHT0IYSQc5VMDn6wpD9KekbSAkmfLhIzVtJ6SXPSny8VKyuEtiRyO+RNJbcp3glcbmZPSeoGPCnpETN7plHcn83szArqCaGlRW6HXGnyEb2ZLTezp9LHG4F/sOcEyiHUnMjtkDcVTTxST9JQ4CgoOpLnzZLmkgzPuMLMFpQoYyIwEaDPkM7c8MGzytZ5Jd9y2/XdqV9wYzgjw2CGk/2BE7bGj3l8+tFuzFvW/81vz6l+XRseb+/XxV/cmDV0dmN6ZRg0kmmQ12J/uw55/7NuzD/mOvt5ywF+W1KV5nZhXtN3CFxdvr4sM0Mtf+IwN4bj/Lze0Mnf391H+THbZ410Y1bU9XBjBsz368pSznY6ujEck+E9NM8vpvtWfz8bfl0Dr/Fnn1txvzNT1eWlV1X8ZaykrsC9wGfMbEOj1U8Bh5jZaOAHlBn/ZWa3mFmdmdV16+uPkAyhuVUjtwvzmu7+9JchNIeKOnpJ7UneCD8zs183Xm9mG8xsU/p4KtBeUp9K6gyhJURuhzyp5KobkUyQ/A8z+06JmAFpHJKOTet7tal1htASIrdD3lRyjv6twIeBeZLmpMu+AAwBMLObgXOASyXtBLYA51tbvItaCA1FbodcaXJHb2aPQflvGczsRuDGptYRQmuI3A55EyNjQwgh56KjDyGEnIuOPoQQcq4qA6aqrR276cbGsjFzeaNf0M0ZKlufYRaZ4/0QZZkZKIP2GSbiKTcwol73p3f4MZOX+AVdlaE9j2WIuS3Dfr7SD8kyI1j7gxpf8t7QzvZ+Gc2hU8/NDBs3q2zMwSz2C8owM1SWwVBbtvnldPPHb3H4rufcmO73+PnIT/2QAR9a78YccN5rboxl2K4s+4csg87u8os5ijluzOJx5S/qWnjt5pLr4og+hBByLjr6EELIuejoQwgh56KjDyGEnIuOPoQQci46+hBCyLno6EMIIeeiow8hhJxrkwOmNtGVv/DWsjFPcJxbzrUPXOfGPMtQN2bo5pfdmE5d/BsXnjDFH1wxY/yxbsyJdRlma3pHhgFKXfwQDva3a9In/LomfShDXf38uh6lnxtzcu8/ll2/ZP/tGRpTfdt2dWTh+vKzQz3a42S/oPH+fsoyM1SWwVB6NUNdV2TItdl+CNMz3PxzbIYBSk9mGJyVYbv69fbrkjPpE5Dp9bqDbm7MkPXl+6Ftu0rPrFWNGaYWSZonaY6kPV5OJW6QtFDS3yX58+mF0Moir0OeVOuI/iQzW11i3enA8PTnOOBH6e8Q2rrI65ALLXGOfhzwE0vMBHpKGtgC9YbQnCKvQ82oRkdvwDRJT6Yz3jc2GBrcqWlJuiyEtizyOuRGNU7dnGBmSyX1Ax6R9KyZ/WlvC0nfTBMBug/pUYVmhVCRquc1Bx9U5SaGkE3FR/RmtjT9vQqYAjS+bGQpcHDB3welyxqXc4uZ1ZlZ3QF9D6i0WSFUpDnyWr17N1dzQyiroo5eUhdJ3eofA6cA8xuFPQB8JL1K4XhgvZktr6TeEJpT5HXIm0pP3fQHpkiqL+vnZvaQpEsAzOxmYCpwBrAQeA34aIV1htDcIq9DrlTU0ZvZi8DoIstvLnhswCf3ptz27KA/K8vGDMowzc7Wzf6Ahze85LdnxsgMg5gu8es65+Y73Zhfzf6w36AMAzl4rx/CkAwxGWaGmpTlosIMs2LxTb+uRVf6r4WXG6soP6CmufLatu/P1kW9ysb8avT73HLewuvdmO2zRroxmWaGyjAY6rHr/SEE+2eYGez4DIO8Zs7b42XZw07auTEnZNiujavauzHPtTvcjemQ4fX6C/7r7uUO20t353ELhBBCyLno6EMIIeeiow8hhJyLjj6EEHIuOvoQQsi56OhDCCHnoqMPIYSci44+hBByrk3OMLWeHjzMqWVj3sMDbjmdVmWozJ+wiBE84wd18kOGssiNub3uAjfmYzPucmOe+vQRbszRM/7hxvB9P6TBPRxLuSlDzBl+yJmbf+PGjOoyr+x6I8OAs2bQqfNmho2eVTZmY4aZhkbPft6NWVHn3xiw+z0ZZmLKMDNUpsFQH57rxkxqfJOJYjEZypl5pz+oKst2Zdk/B1/gJ/+A2evdmIfryvd3ACOd3FnYeXPJdXFEH0IIORcdfQgh5Fx09CGEkHPR0YcQQs5FRx9CCDnX5I5e0uGS5hT8bJD0mUYxYyWtL4j5UsUtDqGZRW6HvGny5ZVm9hwwBkBSO5Jp1KYUCf2zmZ3Z1HpCaGmR2yFvqnXq5h3AC2b2cpXKC6GtiNwONa9aA6bOB0qN4nmzpLnAMuAKM1tQLEjSRGAiQIch/VhHz7IVXvr4HX6r3mpuSJZZqKYz1o0593sPujFvl1/Xe8e4IfC0v11HX5dhUJA/7gam+3XxQIa6NmSoa7xf13/hnyHxcmchr2VozL9VlNuFec3AIcxfcEzZysYfebbfojp/Pw2Yn+E1+akfkuX1zzIzVKbBUObXNSnDe2jSHH9QFfMy5PW7/boGjPIHQ2V5vcY3mGe+uKsW/KB8wJYuJVdVfEQvqQPJxHW/LLL6KeAQMxsN/AC4r1Q5ZnaLmdWZWV37vj0rbVYIFatGbhfmNQf2bba2hlBONU7dnA48ZWZ7TPJqZhvMbFP6eCrQXlKfKtQZQkuI3A65UI2O/gJKfLSVNEBKPmtJOjat79Uq1BlCS4jcDrlQ0Tl6SV2AdwEXFyy7BMDMbgbOAS6VtBPYApxvluFEXAitLHI75ElFHb2ZbQZ6N1p2c8HjG4EbK6kjhNYQuR3yJEbGhhBCzkVHH0IIORcdfQgh5FybnGFqGx1ZyLCyMTPeeqxbzolT/AEPD473R7Cf+7g/qxE3+3UdYQf55VyxxI/5gF/X//35592Ya1d9za/rKr+uV77R1Y15hhFuzImz/br61010Y7bToez6/djtltEcOnXezLAjy88SNIWz3HKuzLCfsswwNeBDGQb7jPXrmjnPn9Epy8xQmQZDfcgNyTTD1PEZtus/X8OXtmJkhv2c4fWaUneZGzPSyZ2YYSqEEPZh0dGHEELORUcfQgg5Fx19CCHkXHT0IYSQc9HRhxBCzkVHH0IIORcdfQgh5FybHDDVj1Vcyg1lYy7LcD+pW8d/wo05lifcmJUnuCH0/6QfM/w2fzDU/def4saM+/40N+baaRkGQ33DD2GVH/LCNw5zY57hCDembtRsN2Y2dW7MnGS615LW8zO3jOawdUsX5s8tP8NUt9E3ueXMrRvuxmynoxtzwHn+TFvdn9zhxuyknRuTZRBTlpmhspSTpT0Z0ogN57V3YxZnmBlqZV0/N6YbG90YL3cqnmFK0u2SVkmaX7Csl6RHJD2f/j6wxHMnpDHPS5qQpb4QWkLkddhXZD11Mxk4rdGyq4E/mNlw4A/p3w1I6gVcAxwHHAtcU+qNE0IrmEzkddgHZOrozexPwJpGi8cB9TN03wFFb9JxKvCIma0xs7XAI+z5xgqhVUReh31FJV/G9jez5enjFUD/IjGDgcUFfy9Jl4XQVkVeh9ypylU36RRqFU2jJmmipNmSZm9+xf+SKITmVu28Zu0rVWpZCHunko5+paSBAOnvYtdnLIUGX0sflC7bg5ndYmZ1ZlbXpe8BFTQrhIo0W15zYN+qNzaELCrp6B8A6q82mADcXyTmYeAUSQemX1adki4Loa2KvA65k/XyyruAvwKHS1oi6SLg68C7JD0PvDP9G0l1km4DMLM1wLXArPTnK+myEFpd5HXYVyg5Ddm2HFnX0e6ZXX42prX0dMs5gSf9yr6fYaYZf0wV/Nzfj7bGr0v/yFDXWzO8Zldk2K6ZGep6LENd12Woa2p16sqyD9/eq/yAsjl1n2Tj7H9maHR17XfUGOs4/dGyMf/qcYhbTt8Mg2s4xt88e9EvRq9WKdf8sXAwPUNdWWaGyjAYiusz5FrvDO/XQzPUNcuv6xW6uTFD1r9cdv22sSez++k5RRsdt0AIIYSci44+hBByLjr6EELIuejoQwgh56KjDyGEnIuOPoQQci46+hBCyLno6EMIIefa5AxTyxjMNXy5bEwHtrvlnHCqP+DhO/5kTXz2RD/md/LrOv3HfjnPX1h+oBjA8Gl+XTdef5Ebc9kX/9dv0Al+Xb957GQ3pv8X/Kmqjpnt13VS3e/cmA5sK7teld2nrMk6ttvGsB4vlI2Z8O87JJc2dYq/nzbM89uzpfxuAqBfhkFDG1f5MzF1v8efqYp3ZxgMdbEfkmVmqG4ZtmtVhrHOnTf7Md0zvF4Txt/rxni5s7Bd6Rc0juhDCCHnoqMPIYSci44+hBByLjr6EELIuejoQwgh59yOXtLtklZJml+w7FuSnpX0d0lTJPUs8dxFkuZJmiMpy41KQ2gxkdthX5HliH4ye85w/wgw0szeCPwT+HyZ559kZmPMLMtdokNoSZOJ3A77ALejN7M/AWsaLZtmZjvTP2eSzJkZQk2J3A77imoMmPoYcHeJdQZMk2TA/5jZLaUKkTQRmAiw/5ABPMOIspW+nuf8ln3YD/nswxkGz0zNMBjqPX4xvOqHDJ+/xA/KMBhm4yn+jDVLv9rLjRl8tj9q5Mxp5WdNAuAUfz8/WzfUjZlxd+MD8CK8mbMWT/LLSFSc24V5Td8hzL//mLIVrh7X22/VID+k+9YMed2pOjMoPdfucDfm4AsWuzEDRq13Y1aM7OHGLG4wb3txxxw6343JNBgqy35+wt/PTzPGjVlxv/NirOtSclVFHb2kLwI7gZ+VCDnBzJZK6gc8IunZ9ChqD+kb5RaAznUj2t78hmGfUq3cLsxrDauLvA6toslX3Ui6EDgT+KCVmHjWzJamv1cBU4Bjm1pfCC0lcjvkTZM6ekmnAVcC7zWz10rEdJHUrf4xcArgf14KoRVFboc8ynJ55V3AX4HDJS2RdBFwI9CN5CPrHEk3p7GDJE1Nn9ofeEzSXOBvwG/N7KFm2YoQmiByO+wr3HP0ZnZBkcVFb3toZsuAM9LHLwKjK2pdCM0ocjvsK2JkbAgh5Fx09CGEkHPR0YcQQs6pxNVjrUqD6oyLyt8+ZNy1d7nlHIw/SOMH86/yG9TRD2G4vx9/nZziLevsx/0ZlHirX9cr+AOmpjPWjTmXB/32fN8fELLi0/5Al+Gbn3djNq3o47dn2EIn4GzM5mWYzqi61LXOGOXcFmdrhoLO8kPstAybtyxDXeP9XJvL692Y0bP915YswwwyzEI2t2643x7+6deVYWaoLIPX9FCG7brPD6GTs35eHbap+A6KI/oQQsi56OhDCCHnoqMPIYSci44+hBByLjr6EELIuejoQwgh56KjDyGEnIuOPoQQcq4aM0xV34HA+eVD7n+u2P2oGrro8BvdmF+PPN2Ned/cqW7Mg7zDjcky8sqO8Et5gFPdmHUZBmf1ZJ1f2V3+oJEbP32RG/Nfn7vNryvLjX4fmpQhyIvJMgKuGWxeBjMnOUHvzFDQCW7EwGtecGOOYo4bc0eGgXd/4X1uzMN1fs6OzzAz1JS6y9yYbmx0YwZl2K4J4+91Y7LMDMUlfghzHssQ9HtnfekRcFluU3y7pFWS5hcsmyRpaXob1zmSivYqkk6T9JykhZKu9uoKoSVFbod9RZZTN5OBYhN1ftfMxqQ/exzySmoH3AScDowALpBUfiLYEFrWZCK3wz7A7ejTeTD9GaL3dCyw0MxeNLPtwC+AcU0oJ4RmEbkd9hWVfBl7maS/px9/DyyyfjA0uKvYknRZUZImSpotaTZrX6mgWSFUrGq53SCvKTozYQjNrqkd/Y+Aw4AxwHLg25U2xMxuMbM6M6vjwL6VFhdCU1U1txvkNQdUoXkh7L0mdfRmttLMdpnZbuBWko+yjS2FBl+jH5QuC6HNitwOedSkjl7SwII/x1P8wrhZwHBJr5PUgeSCyQeaUl8ILSVyO+SRex29pLuAsUAfSUuAa4CxksYABiwCLk5jBwG3mdkZZrZT0mXAw0A74HYzW9AcGxFCU0Ruh31F25xhqmed8TZnJp7f7PAL6tPeDRn9ykw35jiecGOWZZhq5vfr/UFV5/TwB2l8ihvcmGv4shuzkn5uzOEZZuK5kU+6MQNeXenG7Jje3Y3JNPhkmLO+zEw8zUl6o4E3+C7DwLLjJ/kxGa7sHzlulhuzcP1hbszWRb38ukb7dc1fcIxfzpEZypnrl9NpqH+x1bAe/qCz+ff7dfF1P8QfSAfwcWf9GZj9PWaYCiGEfVF09CGEkHPR0YcQQs5FRx9CCDkXHX0IIeRcdPQhhJBz0dGHEELORUcfQgg51zYHTEmvAC8XLOoDrG6l5jRVtLn5NbW9h5hZi985r0hew76zz1vTvtLmknndJjv6xiTNTu7+Vzuizc2v1tpbTK1tQ621F6LNEKduQggh96KjDyGEnKuVjv6W1m5AE0Sbm1+ttbeYWtuGWmsvRJtr4xx9CCGEpquVI/oQQghN1OY7ekmnSXpO0kJJGe6y3bokLZI0T9KcZELotied9HqVpPkFy3pJekTS8+nvYpNit5oSbZ4kaWm6r+dIOqM127g3ai2vIXK7ubREbrfpjl5SO+Am4HRgBHCBpBGt26pMTjKzMW34kq7JwGmNll0N/MHMhgN/INPUFS1qMnu2GeC76b4eY2berB5tQg3nNURuN4fJNHNut+mOnmRi5oVm9qKZbQd+AYxr5TbVPDP7E9B4ip1xwB3p4zuAs1qyTZ4Sba5VkdfNJHK7uLbe0Q8GFhf8vSRd1pYZME3Sk5ImtnZj9kJ/M1uePl4B9G/NxuyFyyT9Pf3426Y+kpdRi3kNkdstrWq53dY7+lp0gpkdTfKx/JOS3t7aDdpbllyKVQuXY/0IOAwYAywHvt2qrcm/yO2WU9Xcbusd/VLg4IK/D0qXtVlmtjT9vQqYQvIxvRaslDQQIP29qpXb4zKzlWa2y8x2A7dSO/u65vIaIrdbUrVzu6139LOA4ZJeJ6kDcD7wQCu3qSRJXSR1q38MnALML/+sNuMBYEL6eAJwfyu2JZP6N29qPLWzr2sqryFyu6VVO7f3r6w5zcvMdkq6DHgYaAfcbmYLWrlZ5fQHpkiCZN/+3Mweat0m7UnSXcBYoI+kJcA1wNeBeyRdRHKHxfNar4V7KtHmsZLGkHwUXwRc3Frt2xs1mNcQud1sWiK3Y2RsCCHkXFs/dRNCCKFC0dGHEELORUcfQgg5Fx19CCHkXHT0IYSQc9HRhxBCzkVHH0IIORcdfQgh5Nz/B9mtw6FLaFC8AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Ignore zero fission rates in guide tubes with Matplotlib color scheme\n",
    "openmc_fission_rates[openmc_fission_rates == 0] = np.nan\n",
    "openmoc_fission_rates[openmoc_fission_rates == 0] = np.nan\n",
    "\n",
    "# Plot OpenMC's fission rates in the left subplot\n",
    "fig = plt.subplot(121)\n",
    "plt.imshow(openmc_fission_rates, interpolation='none', cmap='jet')\n",
    "plt.title('OpenMC Fission Rates')\n",
    "\n",
    "# Plot OpenMOC's fission rates in the right subplot\n",
    "fig2 = plt.subplot(122)\n",
    "plt.imshow(openmoc_fission_rates, interpolation='none', cmap='jet')\n",
    "plt.title('OpenMOC Fission Rates')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# close statepoint file to release HDF5 file handles\n",
    "sp.close()"
   ]
  }
 ],
 "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}