Skip to content

Instantly share code, notes, and snippets.

@paulromano

paulromano/Distributed Cell Tally Visualization.ipynb

Last active June 15, 2021 09:39
Show Gist options
  • Save paulromano/f2fbf3d4731e324b6f5ab31ef3fcaa26 to your computer and use it in GitHub Desktop.
Save paulromano/f2fbf3d4731e324b6f5ab31ef3fcaa26 to your computer and use it in GitHub Desktop.
Download ZIP
Example of visualizing tallies with distributed cell filters
Display the source blob
Display the rendered blob
Raw
Distributed Cell Tally Visualization.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Distributed Cell Tally Visualization\n",
"\n",
"This example demonstrates how a tally with a `DistribcellFilter` can be plotted using the `openmc.lib` module to determine geometry information. First, we'll begin by creating a simple model with a hexagonal lattice."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from pprint import pprint\n",
"\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import openmc\n",
"import openmc.lib"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our model will have two materials, fuel and water:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"fuel = openmc.Material(name='fuel')\n",
"fuel.add_nuclide('U235', 1.0)\n",
"fuel.set_density('g/cm3', 10.0)\n",
"\n",
"water = openmc.Material(name='water')\n",
"water.add_nuclide('H1', 2.0)\n",
"water.add_nuclide('O16', 1.0)\n",
"water.set_density('g/cm3', 1.0)\n",
"\n",
"mats = openmc.Materials((fuel, water))\n",
"mats.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We'll create two universes for use in our hexaongal lattice, one which has a pin consisting of fuel and one which is entirely water."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"r_pin = openmc.ZCylinder(r=0.25)\n",
"fuel_cell = openmc.Cell(fill=fuel, region=-r_pin)\n",
"water_cell = openmc.Cell(fill=water, region=+r_pin)\n",
"pin_universe = openmc.Universe(cells=(fuel_cell, water_cell))\n",
"\n",
"all_water_cell = openmc.Cell(fill=water)\n",
"water_universe = openmc.Universe(cells=(all_water_cell,))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's create a hexagonal lattice with three rings that is otherwise surrounded by water."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"lat = openmc.HexLattice()\n",
"lat.center = (0., 0.)\n",
"lat.pitch = [1.25]\n",
"lat.outer = water_universe\n",
"outer_ring = [pin_universe]*12\n",
"middle_ring = [pin_universe]*6\n",
"inner_ring = [pin_universe]\n",
"lat.universes = [outer_ring, middle_ring, inner_ring]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's put our lattice inside a cylindrical cell so that is properly contained. This completes the geometry."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"outer_radius = 4.0\n",
"outer_surface = openmc.ZCylinder(r=outer_radius, boundary_type='vacuum')\n",
"main_cell = openmc.Cell(fill=lat, region=-outer_surface)\n",
"geom = openmc.Geometry([main_cell])\n",
"geom.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To see what our geometry looks like, we can plot it:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAGQAgMAAAD90d5fAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAACVBMVEX///8AAP///wCN6GYzAAAAAWJLR0QAiAUdSAAAAAd0SU1FB+UGDxccGcYLhQ8AAAb8SURBVHja7Z3NmaMwDIYzB5fgfijBB8iBEugnJeSwVLmbZPd5MoMlffpByTLjq03eSJ+Rwdjy6aQqZbqXprtKVaankoDYCfMxbcoQzShTp7RYRp26ZUxghFLKRJYwjzGMMMrHxJYhBDIJJYJRJUiA+EViRMgiM/wOE50V4LAPhOHtYRjD5zBAdb/2KMNjCmyIxxScYTdFYYjdFA3DaorKEKspOobNFKUhNlOqFmKIYGDUei7D/oYYTDEYojdFLbtFegtD24tN3tL6yyD7reiktzF0/ipWiEZ6o7dU/jLKfiu49GZvafxl9pbCXw5v4f5yeAv3l8NbsL9c3kL95fIW6i+Xt1B/+RhY/HJKgoni9BbmLy8D8ZfbW4i/ih8id2K3JIgofoYsSoAksiglAiKJUiMgkigRDEmUEEkkUUoMhBelxkB4UWIYvChBkvCiFPny87pe5FacKFW+fFnXX3IrThTAEeufAjRzSXK+QQB/0aIU+eL5BrnK7WhRqnzxcoO4RAF8jUImD2RdncojtyIMoZQvkRBK+RoJoZQHLsUhhChQdMQhw+sgJRbSV77GQvrKI1cqIJMdAoeVPgQbehWQnvIFuhIN9VNf+QpdiQ5aU195DIIOvwQEuxJ8kLgXq+7oI9G9DGaIomwhJR6y7V41HjK+BhLP2HavHXTfKp8CKUzbmQ5WMxvHvnavSjc9M4GED2QjDpnpuDvzIfkrhDF6oUeQRRhccAgz4EpjMdy5zhLkQl876CAXZVUPUuiGswS50td+7sNVglyVVY8y5kPodo+HIA4C9+EMCBceXZAhG1KYdh7hP/XhFEjdCzKCEE9Y+QRhmvkgEwjxhPpnCP8UwXTTRYQMIGThh1/+GX+AevDD8xd11d/SQMjE/NtFeiNqUA++/1/mD/CGPPVhAeIpI9aDneV4kF2e6P+VIRNStlXMfbwwVZ2+3kjITI8S2ioawoSkha3amtKoe/FMB1duAqdbNVKQmf4lZirq3K0aqdtkoX9pFqq2rmQhfc/rqyjIaodspaTuRWbKVF81vBBy9kAufUjZF9LUEJkfADFbUveFjMRjRGjvSoQQzaPu+ImAhMYuFnLt/pI+ClOQyPHkAfkgmk/0L/X5E1E1EA9EzBjPvPYs/SoKEvi08oCU7v+lH66YN5KlW9XkNyB/ORak7g0ZvxnE8/YLQzzv8X8hgCGOGYlbQSCuuRUU4pol0kDM810gxDdz94CI79d+yCBD7gOudTb1+0EWCSLeKIM8nPgh7TtB/MK/F+TSqxK/M8IQf1iBIf06FFKlJv5QPwKQZTV/j8ch82pdWaCA2NdIKCD21R4aiLf8QN4Rsj/jz4PX6yDMEkFLVR+yWL869Kv6ECa4Wqq6EGaWxlTVhUiTWtqqLoSZVGNm7may6qT8JVNVF8I8szGzqQtZdVL+kqmqB4md4f6B/EB0v6Sset3NmBJWUgJkSqhPGbRSht+cB4mUR6LgkgOp+zOO9OrwrSAZr9gpkwUp0x4pEzgBU1FFahIxqQZCLsoqJWSWINcwyFVZ9Y6QRYL8N3P1EZCUjzTv8bnpOJCIj5kpn2Xf4wNzzqfyZU346B+xfKEC/nIvxJAhKUtKnOVYkLI35FiLyLqQuA0VD0jKEkXlYstFuaGChMQvG01ZAPuqpbyx04N5y6trH/IfrkbfQPZYvF/2hRxtQ0XKJpfj7AlK2UL1os1goVE4c4Ne7Taf6F/q8/mthqXXPGpDxYu3fwZuqPgHOc6+3+Psxc6B1P0Y8E7/mHQChW8YkxiBhwSleEhJVpGSdiMlgUhOKpS6FwTNHHMfcCPS05S9IC0bcpwMS8fJenWgTGRFglx6VcrEbSkp6FKS6aWkBcxJcFjplvNqXVmgyAcZl3SyMEaHpc88TrbRlOSsB8plW+IhbQNJSWKcko45JbF0TopsDOJM9l2gK2dhwH0qrQNJScCekko+Jyl+jYX00/uXWEjrQo5zrkPKMRg5B3qUSEgjICmHrKQcF5Nz8E2Ng9BH+BT5YjTUNxKScqxSygFROUddFchfgLcaA0k5fizlILWcI+FKDKSxkJRj+lIOHMw5OrFEQJoASTnOMuVgzpwjRosfIkmSdOxrygG2OUfxphwqnHI8cs5Bz8UHkTtwgCiQJDnHiOcciJ5ytHvKIfUef6HecvkL9pbDX7i3HP7CvXWyxy8Nw+ovjbfM0itkt/tLx7BJr5Ld6i+lt0zS62Q3mqI2xGCK3hCD9FrZ70ULsTC0ppgM0ZpiY+hMMRqiM8XK0JhiNkRjip2Bm+IwBDfFw0Aj2OCCYBHMErXUDvMyEO2bGyI7zO0sxGERDKmHDSEQXpYWw2ApYQxG/BDRBUoog/BYoK8epdPHhmjGaXu/7ID4itkJcSvFJMZvJ6jEfJPW6pAAAAAldEVYdGRhdGU6Y3JlYXRlADIwMjEtMDYtMTVUMTY6Mjg6MjUrMDc6MDCADMiJAAAAJXRFWHRkYXRlOm1vZGlmeQAyMDIxLTA2LTE1VDE2OjI4OjI1KzA3OjAw8VFwNQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<IPython.core.display.Image object>"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = openmc.Plot.from_geometry(geom)\n",
"p.color_by = 'material'\n",
"p.colors = {\n",
" fuel: 'yellow',\n",
" water: 'blue'\n",
"}\n",
"p.to_ipython_image()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Great! Now to finish things, we just need to create some simulation settings and setup a tally with a distribcell filter. Beginning with the simulation settings:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"settings = openmc.Settings()\n",
"settings.batches = 25\n",
"settings.inactive = 5\n",
"settings.particles = 10000\n",
"settings.source = openmc.Source(space=openmc.stats.Point((0., 0., 0.)))\n",
"settings.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our tally will use a distribcell filter over the `fuel_cell`, and we'll tally the flux as the only score."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"tally = openmc.Tally()\n",
"tally.filters = [openmc.DistribcellFilter(fuel_cell)]\n",
"tally.scores = ['flux']\n",
"\n",
"tallies = openmc.Tallies([tally])\n",
"tallies.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We're ready to run the transport simulation and get tally results to plot."
]
},
{
"cell_type": "code",
"execution_count": 9,
"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-2021 MIT and OpenMC contributors\n",
" License | https://docs.openmc.org/en/latest/license.html\n",
" Version | 0.12.2\n",
" Git SHA1 | cbfcf908f8abdc1ef6603f67872dcf64c5c657b1\n",
" Date/Time | 2021-06-15 16:28:25\n",
" OpenMP Threads | 12\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 /opt/data/hdf5/nndc_hdf5_v15/U235.h5\n",
" Reading H1 from /opt/data/hdf5/nndc_hdf5_v15/H1.h5\n",
" Reading O16 from /opt/data/hdf5/nndc_hdf5_v15/O16.h5\n",
" Minimum neutron data temperature: 294.0 K\n",
" Maximum neutron data temperature: 294.0 K\n",
" Reading tallies XML file...\n",
" Preparing distributed cell instances...\n",
" Writing summary.h5 file...\n",
" Maximum neutron transport energy: 20000000.0 eV for U235\n",
" Initializing source particles...\n",
"\n",
" ====================> K EIGENVALUE SIMULATION <====================\n",
"\n",
" Bat./Gen. k Average k\n",
" ========= ======== ====================\n",
" 1/1 0.24098\n",
" 2/1 0.21492\n",
" 3/1 0.21007\n",
" 4/1 0.21650\n",
" 5/1 0.20297\n",
" 6/1 0.20552\n",
" 7/1 0.20198 0.20375 +/- 0.00177\n",
" 8/1 0.21231 0.20660 +/- 0.00303\n",
" 9/1 0.19424 0.20351 +/- 0.00376\n",
" 10/1 0.20340 0.20349 +/- 0.00291\n",
" 11/1 0.20235 0.20330 +/- 0.00239\n",
" 12/1 0.20565 0.20364 +/- 0.00204\n",
" 13/1 0.21515 0.20508 +/- 0.00228\n",
" 14/1 0.21415 0.20608 +/- 0.00225\n",
" 15/1 0.19936 0.20541 +/- 0.00212\n",
" 16/1 0.21481 0.20627 +/- 0.00210\n",
" 17/1 0.22121 0.20751 +/- 0.00229\n",
" 18/1 0.21317 0.20795 +/- 0.00215\n",
" 19/1 0.19631 0.20711 +/- 0.00216\n",
" 20/1 0.20903 0.20724 +/- 0.00201\n",
" 21/1 0.20843 0.20732 +/- 0.00188\n",
" 22/1 0.20758 0.20733 +/- 0.00177\n",
" 23/1 0.20702 0.20732 +/- 0.00167\n",
" 24/1 0.20964 0.20744 +/- 0.00158\n",
" 25/1 0.20801 0.20747 +/- 0.00150\n",
" Creating state point statepoint.25.h5...\n",
"\n",
" =======================> TIMING STATISTICS <=======================\n",
"\n",
" Total time for initialization = 4.7180e-01 seconds\n",
" Reading cross sections = 4.5829e-01 seconds\n",
" Total time in simulation = 5.6116e+00 seconds\n",
" Time in transport only = 5.3339e+00 seconds\n",
" Time in inactive batches = 8.6052e-01 seconds\n",
" Time in active batches = 4.7511e+00 seconds\n",
" Time synchronizing fission bank = 4.9694e-02 seconds\n",
" Sampling source sites = 4.7703e-02 seconds\n",
" SEND/RECV source sites = 1.9507e-03 seconds\n",
" Time accumulating tallies = 1.9557e-01 seconds\n",
" Time writing statepoints = 3.9933e-03 seconds\n",
" Total time for finalization = 5.9375e-04 seconds\n",
" Total time elapsed = 6.1034e+00 seconds\n",
" Calculation Rate (inactive) = 58104.4 particles/second\n",
" Calculation Rate (active) = 42095.5 particles/second\n",
"\n",
" ============================> RESULTS <============================\n",
"\n",
" k-effective (Collision) = 0.20881 +/- 0.00110\n",
" k-effective (Track-length) = 0.20747 +/- 0.00150\n",
" k-effective (Absorption) = 0.20813 +/- 0.00159\n",
" Combined k-effective = 0.20898 +/- 0.00115\n",
" Leakage Fraction = 0.89048 +/- 0.00078\n",
"\n"
]
}
],
"source": [
"openmc.run()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now that the simulation is finished, we need to pull the resulting flux data out of the statepoint file that was produced. Since we only have a single score and a single filter, this is especially easy because we can just take the `.mean` attribute of a `Tally` object and flatten it using the `ravel()` method for numpy arrays. While we're at it, we'll get a Pandas dataframe to see what the results look like."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" level 1 level 2 level 3 distribcell nuclide score mean \\\n",
" univ cell lat univ cell \n",
" id id id x y id id \n",
"0 4 4 3 0 -2 1 1 0 total flux 3.71e-02 \n",
"1 4 4 3 1 -2 1 1 1 total flux 4.00e-02 \n",
"2 4 4 3 2 -2 1 1 2 total flux 3.46e-02 \n",
"3 4 4 3 -1 -1 1 1 3 total flux 3.97e-02 \n",
"4 4 4 3 0 -1 1 1 4 total flux 4.73e-02 \n",
"5 4 4 3 1 -1 1 1 5 total flux 4.79e-02 \n",
"6 4 4 3 2 -1 1 1 6 total flux 3.99e-02 \n",
"7 4 4 3 -2 0 1 1 7 total flux 3.62e-02 \n",
"8 4 4 3 -1 0 1 1 8 total flux 4.83e-02 \n",
"9 4 4 3 0 0 1 1 9 total flux 5.13e-02 \n",
"10 4 4 3 1 0 1 1 10 total flux 4.79e-02 \n",
"11 4 4 3 2 0 1 1 11 total flux 3.49e-02 \n",
"12 4 4 3 -2 1 1 1 12 total flux 4.12e-02 \n",
"13 4 4 3 -1 1 1 1 13 total flux 4.74e-02 \n",
"14 4 4 3 0 1 1 1 14 total flux 4.81e-02 \n",
"15 4 4 3 1 1 1 1 15 total flux 4.09e-02 \n",
"16 4 4 3 -2 2 1 1 16 total flux 3.77e-02 \n",
"17 4 4 3 -1 2 1 1 17 total flux 3.96e-02 \n",
"18 4 4 3 0 2 1 1 18 total flux 3.69e-02 \n",
"\n",
" std. dev. \n",
" \n",
" \n",
"0 4.11e-04 \n",
"1 6.21e-04 \n",
"2 5.68e-04 \n",
"3 5.90e-04 \n",
"4 5.93e-04 \n",
"5 4.61e-04 \n",
"6 3.57e-04 \n",
"7 3.91e-04 \n",
"8 4.51e-04 \n",
"9 4.54e-04 \n",
"10 4.58e-04 \n",
"11 6.19e-04 \n",
"12 5.54e-04 \n",
"13 5.99e-04 \n",
"14 4.38e-04 \n",
"15 4.50e-04 \n",
"16 3.89e-04 \n",
"17 4.67e-04 \n",
"18 3.21e-04 \n"
]
}
],
"source": [
"with openmc.StatePoint('statepoint.25.h5') as sp:\n",
" # Get the Tally object\n",
" t = sp.tallies[1]\n",
" \n",
" # Get the mean value of the flux for each instance of the fuel cell as a flattened (1D) numpy array\n",
" flux = t.mean.ravel()\n",
" \n",
" # Show a Pandas dataframe\n",
" df = t.get_pandas_dataframe() \n",
" print(df)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now it's time to actually try to plot the data. There's no \"easy\" way to do this, but the technique we'll use is to loop over a bunch of (x, y, z) points, determine what cell/instance corresponds to each point, and create a 2D array with the flux values where the cell/instance matches one of the bins of our tally. First, let's set the resolution of our desired image and create an array of the appropriate size filled in with \"null\" values indicated by `np.nan` (matplotlib won't plot data anywhere it sees nan values). When an (x, y, z) point corresponds to a cell/instance that we tallied, we'll set the corresponding pixel in the image."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"resolution = (600, 600)\n",
"img = np.full(resolution, np.nan)\n",
"xmin, xmax = -3., 3.\n",
"ymin, ymax = -3., 3."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In order to determine the cell/instance for a given (x, y, z) point, we need to use the `find_cell` function from `openmc.lib`. Calls to `find_cell` return a tuple of two values, a `Cell` object and distribcell index that can be used to determine the appropriate bin from the tally. Note that we actually have to initialize the OpenMC library (which reads the XML files) in order to use this function. We'll use the `run_in_memory()` context manager, which will call `openmc.lib.init()` and `openmc.lib.finalize()` before and after the `with` block, respectively."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"with openmc.lib.run_in_memory():\n",
" for row, y in enumerate(np.linspace(ymin, ymax, resolution[0])):\n",
" for col, x in enumerate(np.linspace(xmin, xmax, resolution[1])):\n",
" try:\n",
" # For each (x, y, z) point, determine the cell and distribcell index\n",
" cell, distribcell_index = openmc.lib.find_cell((x, y, 0.))\n",
" except openmc.exceptions.GeometryError:\n",
" # If a point appears outside the geometry, you'll get a GeometryError exception.\n",
" # These lines catch the exception and continue on\n",
" continue\n",
"\n",
" if cell.id == fuel_cell.id:\n",
" # When the cell ID matches, we set the corresponding pixel in the image using the\n",
" # distribcell index. Note that we're taking advantage of the fact that the i-th element\n",
" # in the flux array corresponds to the i-th distribcell instance.\n",
" img[row, col] = flux[distribcell_index]\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now plot the image. Because we looped over y values from the bottom to top, we need to specify that the origin corresponds to the lower-left corner rather than the upper-left corner. We'll also set a few optional arguments to make the plot nicer."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.colorbar.Colorbar at 0x152845135460>"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAEKCAYAAAB0cRxpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhNUlEQVR4nO3df7BcZZ3n8ffn3vDDBQQMrLIESFyiiAyiXkEddFTEDbuWUQvKMJbgDk6KclBnXUdRZhCZ0RKtldGBcs0Ca0AdYEHWrIZBFEe0hJgbDYQQ0RBUQqEkgEgGwST3u3/06aRvp2/3093n9ulz+vOqOkX3OU+f+73RfPI853nOaUUEZmY2s7GiCzAzG3YOSjOzDhyUZmYdOCjNzDpwUJqZdeCgNDProLCglLSvpB9LukvSekmfKKoWMxtekhZJuk/SRknntzi+j6TrsuOrJM1vOHa8pDuyjFknad9s/8uz9xslfUGS2tVQZI/yGeANEfES4ARgkaRXFliPmQ0ZSePA5cBpwLHAmZKObWp2DvB4RBwNXApckn12DvAV4NyIeDHwOmB79pkvAn8JLMy2Re3qKCwoo2Zb9navbPPqdzNrdCKwMSI2RcQfgWuBxU1tFgPLs9c3AKdkPcQ3AXdHxF0AEfFoROyUdBjw7Ii4M2p33FwNvLVdEXNy+3V6kP1rsQY4Grg8Ila1aLMUWAqw3377vfyYY44ZbJFmI2TNmjVbI+LQfs7xEu0XT7Izqe0DPLMeeLph17KIWNbw/nDgwYb3m4GTmk6zq01E7JD0BDAXeAEQkm4BDgWujYjPZO03N53z8HZ1FhqUEbETOEHSQcBNko6LiHua2iwDlgFMTEzE5OTk4As1GxGSftXvOZ5kJ//AUUlt38nPn46IiX5/5gzmACcDrwCeAr4raQ3wRLcnGopZ74j4HfA9OlwnMLNyGBtL2xI8BBzR8H5etq9lm+y65IHAo9R6irdHxNaIeApYCbwsaz+vwzmn/z5Jpc4CSYdmPUkkPQs4FfhZUfWYWT4kmDMnbUuwGlgoaYGkvYElwIqmNiuAs7PXpwO3ZdcebwH+RNK/ywL0z4B7I+Jh4PeSXpldyzwL+Ea7Iooceh8GLM+uU44B10fENwusx8xyIJJ7ix1l1xzPoxZ648BVEbFe0sXAZESsAK4ErpG0EXiMWpgSEY9L+hy1sA1gZUR8Kzv1e4EvA88Cbs62GRUWlBFxN/DSon6+mc2esbarErsTESupDZsb913Y8Ppp4IwZPvsVakuEmvdPAsel1lDoZI6ZVZDy61EOCwelmeUqz6H3sHBQmlm+3KM0M2tPJM9ol0bFfh0bFp+69b5p7z926gsLqsQGzj1Ks/aaA7J5vwOz+qp4jbJiv44VaaaQ7LaNlZ+kpK0sHJSWi24C0GFZccr1FsahUKJSbVj1EnwOy+rK+RbGoeCgtL70E3gOy2qqX6OsUo+yRJluZmVRphBM4aA0s3x5eZCZWXtVXB7koDSzfLlHaWbWnoA540VXka+K5b4NWj932vgunYqq4DpK9yjNLFdVvEZZsV/HitBLz9C9ySoTY2NpW1k4KC0X3QSfQ7LiBBpX0lYWDkrLTUoAOiSrT4KxOWNJW1n4GqXlykFoQKl6iykclGaWLwmV6PpjCgelmeXOPUozszYkSjWjncJBaWa589DbzKwdibG9qnUPo4PSzHIl+RplbiQdAVwNPBcIYFlEfL6oesrmW5u27rHvvzz/kAIqKZ/J3z65x76J5x5QQCXVlefQW9Ii4PPAOHBFRHy66fg+1LLk5cCjwDsi4peS5gMbgPqj9O+MiHOzz/wrcBjwh+zYmyLikZlqKLJHuQP47xHxE0kHAGsk3RoR9xZY09BrFZDNxxyYrbUKyOZjDswcKL+glDQOXA6cCmwGVkta0ZQT5wCPR8TRkpYAlwDvyI7dHxEnzHD6d0bEZEodhS2Nj4iHI+In2esnqSX/4UXVUwbtQrKXdqOkXUj20s7aSbt9MXF4fiKwMSI2RcQfgWuBxU1tFgPLs9c3AKco5+/CHYp7iLIu8kuBVQWXMrS6DT+H5W7dhp/Dsj/15UE5PRTjcODBhveb2bNDtatNROwAngDmZscWSPqppO9Lek3T5/63pLWS/q5TsBYelJL2B24E/joift/i+FJJk5Imt2zZMvgCS8xh6dArhGBsr7GkDTik/vc725bmWMnDwJER8VLgg8DXJD07O/bOiPgT4DXZ9q52Jyo0KCXtRS0kvxoRX2/VJiKWRcREREwceuihgy1wSDjwBs8B2x+NKWkDttb/fmfbsqZTPQQc0fB+XravZRtJc4ADgUcj4pmIeBQgItYA9wMvyN4/lP33SeBr1Ib4MyosKLOu7pXAhoj4XFF1VN0oh6zDriDK9RrlamChpAWS9gaWACua2qwAzs5enw7cFhEh6dBsMghJzwcWApskzZF0SLZ/L+DNwD3tiihy1vtPqXV310lam+37WESsLK6k4TPKQVe0yd8+6VnwHghQTo84j4gdks4DbqG2POiqiFgv6WJgMiJWUOtwXSNpI/AYtTAFeC1wsaTtwBRwbkQ8Jmk/4JYsJMeB7wD/q10dhQVlRPyQ2p+pmVVJjsuDALLO08qmfRc2vH4aOKPF526kdmmvef+/UVtzmazwyRxrz2sii+PeZG8kMb7XWNJWFuWp1HoyykHroCuOvwrCBm6Uw64oDtk+qKtZ71JwUFaYA9aBV5hxpW0l4aAsiW5DzyG5W7dh6XDtk3uUVqTU8HNI7ik1/BySeRCMj6VtJeHnUZZMPQRnWl/pkJzZxHMPmHERugMyRwKVaEY7hYOypByIvXEgDoCAEg2rUzgozSxn5ZqoSeGgNLNcKec7c4aBg9LM8leiiZoUDkozy5d7lGZmHUiwd7V6lNX6bWbJP/3ogaJLKKWpBy5l6oFLiy6jlN6z7I6iS+hL1Racu0c5g+ZwbHz/vlcvGHQ5pdEqGBv3jS34b4Msp1Saw7Hx/RVLXzXocnonPOtddSm9x3obB+Z0Kb3HqQcudVg2Sek91tuUIzAFOT24d1hU67fpk4fYvetmiO3heO9KMSSXH7NmDRysNb0En8OyphTB1y0Be42lbSVRnkpnWa+hN+ph2U/gjXpY9hqSwx+uaRM5nswxs9ElvODczKytCs56Vyv2e9Tv8HnUh9/WvX6Hz8M8/K59XW21ht4Oyhx4mZB1qxzLfHqkxK+BKFGv00GJg87KZ6iDNntwb8pWFr5GaWb584Lzauq1VznqvdF+7rIZ9Tt0eu0VDnVvEth1Z07KVhLlqdSGVi+BN+ohWWnCQVll3fYOR7032aib4HNI7tZt73D4e5OZMaVtJVFoUEq6StIjku4pso5G73v1go4BmNJmFKUEoENyT1csfVXHAExpMzSU79Bb0iJJ90naKOn8Fsf3kXRddnyVpPnZ/vmS/iBpbbb9z4bPvFzSuuwzX5DUNrWLnsz5MnAZcHXBdeyhHoR+vFp3GoOwfouiwzFNPQhL+3i1XQRz8okWSePA5cCpwGZgtaQVEXFvQ7NzgMcj4mhJS4BLgHdkx+6PiBNanPqLwF8Cq4CVwCLg5pnqKDQoI+L2evoPK4dj7xyQvSlnODaoX6PMx4nAxojYBCDpWmAx0BiUi4GLstc3AJe16yFKOgx4dkTcmb2/GngrbYJy6K9RSloqaVLS5JYtW4oux8xSpA+9D6n//c62pU1nOhx4sOH95mxfyzYRsQN4ApibHVsg6aeSvi/pNQ3tN3c45zRFD707iohlwDKAiYmJKLgcM+ukfo0yzdaImJilSh4GjoyIRyW9HPi/kl7cy4mGPijNrITym9F+CDii4f28bF+rNpslzQEOBB6NiACeAYiINZLuB16QtZ/X4ZzTDP3Q28xKRtlkTsrW2WpgoaQFkvYGlgArmtqsAM7OXp8O3BYRIenQbDIISc8HFgKbIuJh4PeSXpldyzwL+Ea7IopeHvTPwB3ACyVtlnROkfWYWT6ksaStk+ya43nALcAG4PqIWC/pYklvyZpdCcyVtBH4IFBfQvRa4G5Ja6lN8pwbEY9lx94LXAFsBO6nzUQOFD/rfWaRP7/s4uFlu17rsOZr4NZO/PoLu17ryPcXWEkV5fvlYhGxktoSnsZ9Fza8fho4o8XnbgRunOGck8BxqTX4GmUJNQZkq30OzdYaw7HVfgdmTvJdHjQUHJQl0iogZ2rnsJxuppBs1caBmYMS3Z6YolqxX2GpIdlr+ypLCcl+2luTnG9hHAblqdS65rB06BUmv1nvoeCgLAEH3uA5YPvgHqWVzSiHrMOuQA5KG6RRDrqiOWh7JCr3PMryXCQws5LIdx3lMKjWb1NBXuZTHC8T6lG+tzAOhfJUamblkXB7YplU67exPYxyj9Q9wgJpLG0rifJUOsJGOeyK4pDtQwWXB3noXRI6bGlXM+AO15p64HUzg+2QzEH77+oqnfJEujn8BsAhmZOKDb3doyyZTj1Lh2lrKT1Lh2ROJBivVrRU67cZEQ7D3jkMB0Gl6i2mmDEoJb0s4fPbI2JdjvWYWRWMSlAC36f2fRXtrsouAObnWZCZVUDFJnPaBeXqiHhDuw9Lui3nesys9EZo6N0pJFPbmNmIGdWvgpB0PLUh9q72EfH1WarJzEpNaKxa88QdfxtJVwHHA+uBqWx3AA5KM9uTGJ2hd4NXRsSxs17JkLl767YZjx1/yP4DrKR8dl57dsv940uWD7iScrn1V4/NeOzUo54zwEr6pcpN5qTE/h2SRioo24VkyvFRNlNIdjo26tqFZMrxoTOC93pfTS0sfwM8Q61jHRFx/KxWVpDUELx76zb3LBukhmC9nXuXu6WG4K2/eqw8PcsRHHpfCbwLWMfua5SV1G1P0WFp/eq2p1iKsJSgYpM5KbG/JSJWRMQDEfGr+jbrlQ1Yr8NpD8N7G1J7GN77cHr4h+HVe8xaSqU/lfQ1SWdKent9y+OHS1ok6T5JGyWdn8c5bbD6CTyHZYXl+PSgTjkhaR9J12XHV0ma33T8SEnbJH2oYd8vJa2TtFbSZKcaUvrHz6J2bfJNDfv6Xh4kaRy4HDgV2AyslrQiIu7t57xmVjCR26x3Yk6cAzweEUdLWgJcAryj4fjngJtbnP71EbE1pY6OQRkR/zXlRD04EdgYEZsAJF0LLAYGHpT9Dp99rdK61e/webivVeZ6C2NKTiwGLspe3wBcJkkREZLeCjwA/Fs/RXT8bSQtl3RQw/uDs0Xo/ToceLDh/eZsX/PPXyppUtLkli1bcvixZjbr0ofeh9T/fmdb8zMEU3JiV5uI2AE8AcyVtD/wEeATLSoM4NuS1rT4mXtIGXofHxG/23X2iMclvTThc7mIiGXAMoCJiYmYjZ9x/CH799WrdG/SunXqUc/pq1c5vL1JADGV/qjbrRExMUuFXARcGhHbtOelgJMj4iFJ/x64VdLPIuL2mU6U8tuMSTo4Ih4HkPScxM918hBwRMP7edk+Myu5YDyvU6XkRL3NZklzgAOBR4GTgNMlfQY4CJiS9HREXBYRDwFExCOSbqI2xJ8xKFMuJPwPagvO/17S3wM/Aj6T8LlOVgMLJS2QtDewBFiRw3ltgPpZOO5F59UUiIixpC1BSk6sAOpLKE4Hboua10TE/IiYD/wj8KmIuEzSfpIOAJC0H7WJ6nvaFdGx0oi4Gng78Ntse3tEXJPyG3Y47w7gPOAWYANwfUSs7/e8vep1+Oxhd2+B55Dsffg83MPummAsaet4nhlyQtLFkt6SNbuS2jXJjcAHgU5LDZ8L/FDSXcCPgW9FxL+0+0DSEDqbis99NjoiVgIr8z5vr7q9VumQtH51e62yDCEJYiqtt5ikVU5ExIUNr58GzuhwjosaXm8CXtJNDTP+NpJ+0unDKW3KJjX8HJLTjS9ZntRLTG03SlLDrxwhWROMJ21l0a5H+SJJd7c5LmoXTSunU8/SITmz8SXL/Zi1HnTqWZYqJENMRbXu9W732xyT8PmdeRUybByGvXMg9qZMYdjJVNI8cXm0+86cyj34wswGQakz2qVRrf6xmRUuyHUd5VBwUJpZzqrXo0y51/t9kg4eRDFmVn6B2BlzkraySIn951J7tNH12XPhqvWtQWaWu6otD0q5M+dvgYXUVr+/G/iFpE9J+o+zXJuZlVKutzAOhdQ7cyL7crHfADuAg4EbJN0aER+ezQJtT1O/+OyMx8YW/s0AKymf+PUXZjymI98/wEqqLeX2xDLpGJSSPgCcBWwFrgD+JiK2SxoDfgE4KAeoXUg2HndgTtcuIBvbOCxzEJSqt5gipUf5HGoPwpi2rjIipiS9eXbKslY6haT1z2HZv0BMlej6Y4qUr4L4eJtjG/Itx2bSbUi6Z1mT0pNs9RmHZT+qdwtjtfrHFeWe5OD1ErC2W9Umc8pTqfVklEPWYVeM+tA7ZSuLavWPzWwolKm3mKJav00FjXKPsGjukfZKuT3hfFi4R2lmuQpgKmblC1MLU55IH1GjPmtdJM989yiCHVNpW1k4KCtulIPWQVeMWo8ybSsLB2UJjHLYFcUh25+piKStLByUFeaAdeAVIcJBaQXpNvQckrt1G5YO1/7tjLStLDzrXSL18Ou0ZMghuScd+f6Oy30ckPkIKNVETQoHZQm1CkyHY2eNQdgYmg7I/JVpWJ3CQVliDsfeORxnT5RsRjuFg9LMclauiZoUhUzmSDpD0npJU5ImiqjBzGZH/c6cvGa9s+/quk/SRknntzi+j6TrsuOrJM1vOn6kpG2SPpR6zmZFzXrfA7wduL2gn29msyivBeeSxoHLgdOAY4EzJR3b1Owc4PGIOBq4FLik6fjngJu7POc0hQRlRGyIiPuK+NlmNrsiyPMWxhOBjRGxKSL+CFwLLG5qsxhYnr2+ATil/m2xkt4KPACs7/Kc0wz9OkpJSyVNSprcsmVL0eWYWYIuht6H1P9+Z9vSplMdDjzY8H5ztq9lm4jYATwBzJW0P/AR4BM9nHOaWZvMkfQd4HktDl0QEd9IPU9ELAOWAUxMTAz8CvFHbrp71+tL3nb8oH98qcVvrgBAz3tPwZWUy7s+/4Ndr6/5wGsKrKQ30d1kztaImK15iouASyNiW9bB7NmsBWVEvHG2zj0IjQHZvM+BObN6OM60z6E5s8aAbN5XtsDMcXnQQ8ARDe/nZftatdksaQ5wIPAocBJwuqTPAAcBU5KeBtYknHMaLw9q0iogZ2rjwJyuVUi2auOwnK5VQM7UpgyBGcDO/JYHrQYWSlpALcyWAH/e1GYFcDZwB3A6cFtEBLDrD0vSRcC2iLgsC9NO55ymqOVBb5O0GXgV8C1JtxRRR7OUkOynfZWlhGQvbasuJST7aV+IHCdzsmuO5wG3ABuA6yNivaSLJb0la3YltWuSG4EPAm2X+8x0znafKaRHGRE3ATcV8bMtf70En3uW1ZX3E84jYiWwsmnfhQ2vnwbO6HCOizqds52hn/UelF57h6Peq+yndzjqPctee4dl6FVW7cG9vkZpZrnyd+ZU1Kj3Cq18hrpXGTA1lbaVhYMyBw5a69ZQB13fgqmptK0sHJT0v8zHy4SsW/0u8xnmZUIRsH3nVNJWFr5GaWa5CihVbzGFg9LMcufJnIrqdfg86sPuftZCjvo6yl6Hz8M87AaI8DVKazDqIVnXS+CNekjWDXvo9cpBWWEOvt51E3wOyd6VIVgjqheUvkbZpB6W7Zb8OFBb0/Pe0/FuG4fknurh127JUBkCsi6gVDPaKRyUM2gMw4/cdLfDMVFjEPp5lN1pDMN3ff4HpQrHaaJcvcUUDsoEDsneOCB7V9qQzDgozczaqOK93g5KM8tXuEdpZtZWEJ7MMTNrJ6JcTwZK4aA0s9x56G1D4eYHHm25/7QFcwdcSbl899ePt9x/ypEHD7iSanNQWqFmCsjG4w7L1mYKycZjDsz+RVRv1tu3MJZIp5Dstt0oaReSvbSz9qp2C6ODsiS6DT+H5W7dhp/Dsj8RwfYdU0lbWTgoK8xh6dArQn3BecpWFg7KEnDgDZ4Dtj8eelupjHLIOuwKUsHHrDkoh9woB13RHLS9iQp+C6OXB5lZvgKmdpYnBFMU0qOU9FlJP5N0t6SbJB1URB1l4DWRxfGayt5EwM4dO5O2FJIWSbpP0kZJ57c4vo+k67LjqyTNz/afKGlttt0l6W0Nn/mlpHXZsclONRQ19L4VOC4ijgd+Dny0oDoqb5SD1kFXlPyG3pLGgcuB04BjgTMlHdvU7Bzg8Yg4GrgUuCTbfw8wEREnAIuAL0lqHEW/PiJOiIiJTnUUEpQR8e2I2JG9vROYV0QdZTHKYVcUh2zvcv7OnBOBjRGxKSL+CFwLLG5qsxhYnr2+AThFkiLiqYac2ZfayqWeDMNkzl8ANxddRBU5YB14RZnaGUlbgsOBBxveb872tWyTBeMTwFwASSdJWg+sA85tCM4Avi1pjaSlnYqYtaCU9B1J97TYFje0uQDYAXy1zXmWSpqUNLlly5bZKnfodRt6Dsndug1Lh2ufuutRHlL/+51tHUOrq1IiVkXEi4FXAB+VtG926OSIeBm1If1fSXptu/PM2qx3RLyx3XFJ7wbeDJwSMfMS/YhYBiwDmJiYqNZUWpdOWzA3abmQQ3JPpxx5cNJyH4dk/4JInqgBtna4RvgQcETD+3nZvlZtNmfXIA8Epv1FiYgNkrYBxwGTEfFQtv8RSTdRG+LfPlMRhSwPkrQI+DDwZxHxVBE1lFU9BFsFpgOyvXoItgpMB2SO8v0qiNXAQkkLqAXiEuDPm9qsAM4G7gBOB26LiMg+82BE7JB0FHAM8EtJ+wFjEfFk9vpNwMXtiihqHeVlwD7ArZIA7oyIcwuqpZQcir1zKM6uIL91lFnInQfcAowDV0XEekkXU+sZrgCuBK6RtBF4jFqYApwMnC9pOzAFvDcitkp6PnBTlj1zgK9FxL+0q6OQoMym8c2sinL+crGIWAmsbNp3YcPrp4EzWnzuGuCaFvs3AS/ppgbfmWNmuarfwlglDkozy1dAOCjNzGYWATu2l+ehvCkclGaWs2CqYt9X66A0s1wF/hZGs7Y+/d2ftz1+/ikvGFAlVpicZ72HwTDc620V0SkkU9tY+VXtwb0OSstFNwHosKy2CIidkbSVhYPS+tZL8DksKyyCHTumkray8DVKM8tVAFGxWW/3KK0v/fQM3ausqAoOvd2jNLOche/MMTNrK+tRVomD0szyFRDbkx/cWwoOSjPLVXjobWbWQQWH3p71tr70c0uib2esrpiKpK0sHJTWt14CzyFZcVNTaVtJOCgtF90En0Oy4iJtDWWZhucOSstNSgA6JKsvAqa2TyVtZeHJHMtVPQib77pxQI6WMl1/TOGgtFnhYBxhEbCzPL3FFA5KM8uXv1zMzKyzMk3UpHBQmlm+3KM0M2svInyvt5lZJx56m5m1U8GhdyELziX9vaS7Ja2V9G1J/6GIOsxsNuR7Z46kRZLuk7RR0vktju8j6brs+CpJ87P9J2YZs1bSXZLelnrOZkXdmfPZiDg+Ik4AvglcWFAdZpa3ILd7vSWNA5cDpwHHAmdKOrap2TnA4xFxNHApcEm2/x5gIsuZRcCXJM1JPOc0hQRlRPy+4e1+1P5ozawCcr6F8URgY0Rsiog/AtcCi5vaLAaWZ69vAE6RpIh4KiJ2ZPv3ZXfOpJxzmsKuUUr6JHAW8ATw+jbtlgJLs7fPSLpnAOWlOgTYWnQRDYatHhi+mlxPey/s9wR/eOrXt/z0x391SGLzfSVNNrxfFhHLGt4fDjzY8H4zcFLTOXa1iYgdkp4A5gJbJZ0EXAUcBbwrO55yzmlmLSglfQd4XotDF0TENyLiAuACSR8FzgM+3uo82R/asuyckxExMVs1d8v1dDZsNbme9ppCqycRsSiPWvIQEauAF0t6EbBc0s29nGfWgjIi3pjY9KvASmYISjMbaQ8BRzS8n5fta9Vms6Q5wIHAo40NImKDpG3AcYnnnKaoWe+FDW8XAz8rog4zG3qrgYWSFkjaG1gCrGhqswI4O3t9OnBbRET2mTkAko4CjgF+mXjOaYq6RvlpSS8EpoBfAecmfm5Z5yYD5Xo6G7aaXE97Q1VPdk3xPOAWYBy4KiLWS7oYmIyIFcCVwDWSNgKPUQs+gJOB8yVtp5Y1742IrQCtztmuDkV4wtnMrB0/4dzMrAMHpZlZB6ULymG7/VHSZyX9LKvpJkkHFVzPGZLWS5qSVNiyk25vERtAPVdJemRY1uFKOkLS9yTdm/3v9YGC69lX0o+zW/3WS/pEkfUMm9Jdo5T07PqdPZLeDxwbEamTQbNRz5uozbLtkHQJQER8pMB6XkTtwvWXgA9FRN/r4nqoYRz4OXAqtcW8q4EzI+LeQdfSUNNrgW3A1RFxXFF1NNRzGHBYRPxE0gHAGuCtRf0ZSRKwX0Rsk7QX8EPgAxFxZxH1DJvS9SiH7fbHiPh2w21Sd1Jbk1VkPRsi4r4ia6CHW8RmW0TcTm1GdChExMMR8ZPs9ZPABmp3mBRVT0TEtuztXtlWrl7ULCpdUELt9kdJDwLvZLgeqPEXQE8r/yum1S1ihYXAsMuedvNSYFXBdYxLWgs8Atya3dViDGlQSvqOpHtabIsBIuKCiDiC2l095xVdT9bmAmBHVlPh9Vg5SNofuBH466bR0sBFxM7sSTvzgBMlFX6JYlgM5YN7h+32x071SHo38GbglBjARd8u/nyK0vUtYqMouxZ4I/DViPh60fXURcTvJH2P2qPJhmLyq2hD2aNsZ9huf5S0CPgw8JaIeKrIWoZI17eIjZps8uRKYENEfG4I6jm0vmJD0rOoTcT51uJMGWe9b6T2KKhdtz9GRGG9ley2qX3YfRP+nQXPwr8N+CfgUOB3wNqI+E8F1PGfgX9k9y1inxx0DU31/DPwOmqPNfst8PGIuLLAek4GfgCso/b/ZYCPRcTKguo5ntozHcepdaCuj4iLi6hlGJUuKM3MBq10Q28zs0FzUJqZdeCgNDPrwEFpZtaBg9LMrAMHpZlZBw5K65uk+ZL+kN0nnMf5vidpW5GPiTNr5KC0vNyf3Sfct4h4PTDwx8OZzcRBaW1JekX2UOJ9Je2XPdS148MSJJ2Vfe4uSddk+74s6YuS7pS0SdLrsgfqbpD05Vn/Zcx6NJQPxbDhERGrJa0A/gF4FvCViGj7oARJLwb+Fnh1RGyV9JyGwwcDrwLeQu3+7z8F3gOslnRCRKydhV/DrC8OSktxMbUHXTwNvD+h/RuA/1P/atCIaHxg7v/LvnN5HfDbiFgHIGk9MB9Ym2PdZrnw0NtSzAX2Bw4A9u3zXM9k/51qeF1/73+4bSg5KC3Fl4C/o/b8z0sS2t8GnCFpLkDT0NusdPwvuLUl6Sxge0R8LfvSsB9JekNE3DbTZyJivaRPAt+XtBP4KfDuwVRslj8/Zs36ln3nyzfz/HZDSf9KQd8iadbMQ2/Lw07gwDwXnAPPB7bncT6zfrlHaWbWgXuUZmYdOCjNzDpwUJqZdeCgNDPr4P8Dx6SuSes6KRMAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"options = {\n",
" 'origin': 'lower',\n",
" 'extent': (xmin, xmax, ymin, ymax),\n",
" 'vmin': 0.03,\n",
" 'vmax': 0.06,\n",
" 'cmap': 'RdYlBu_r',\n",
"}\n",
"plt.imshow(img, **options)\n",
"plt.xlabel('x [cm]')\n",
"plt.ylabel('y [cm]')\n",
"plt.colorbar()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Success! As expected, the highest flux is observed in the center of the model."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Cell Instance Filters\n",
"\n",
"This same technique can be used with a `CellInstanceFilter` with a slight modification. With a `CellInstanceFilter`, we cannot assume the index in the flux array will corresponding directly to the distribcell index. As such, we'll need to create mapping of the distribcell index to the flux array index. Let's go through such an example. First, let's create a new tally with a `CellInstanceFilter` that requests every other instance of the fuel cell rather than all of them (as is the default for a `DistribcellFilter`):"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Create a cell instance\n",
"cell_instances = [(fuel_cell, i) for i in range(0, 18, 2)]\n",
"cellinst_filter = openmc.CellInstanceFilter(cell_instances)\n",
"\n",
"instance_tally = openmc.Tally()\n",
"instance_tally.filters = [cellinst_filter]\n",
"instance_tally.scores = ['flux']\n",
"\n",
"# Add to existing Tallies object and re-export\n",
"tallies.append(instance_tally)\n",
"tallies.export_to_xml()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's re-run OpenMC to generate a new statepoint and get the results from the statepoint as before:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" cellinstance nuclide score mean std. dev.\n",
" cell instance \n",
"0 1 0 total flux 3.71e-02 4.11e-04\n",
"1 1 2 total flux 3.46e-02 5.68e-04\n",
"2 1 4 total flux 4.73e-02 5.93e-04\n",
"3 1 6 total flux 3.99e-02 3.57e-04\n",
"4 1 8 total flux 4.83e-02 4.51e-04\n",
"5 1 10 total flux 4.79e-02 4.58e-04\n",
"6 1 12 total flux 4.12e-02 5.54e-04\n",
"7 1 14 total flux 4.81e-02 4.38e-04\n",
"8 1 16 total flux 3.77e-02 3.89e-04\n"
]
}
],
"source": [
"openmc.run(output=False)\n",
"\n",
"with openmc.StatePoint('statepoint.25.h5') as sp:\n",
" t = sp.tallies[2]\n",
" flux_inst = t.mean.ravel()\n",
" \n",
" df = t.get_pandas_dataframe()\n",
" print(df)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see that in this case, the i-th value in the array of fluxes is **not** the i-th distribcell instance, so we need to create a dictionary mapping the instance to the index to be used later when plotting."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"{0: 0,\n",
" 2: 1,\n",
" 4: 2,\n",
" 6: 3,\n",
" 8: 4,\n",
" 10: 5,\n",
" 12: 6,\n",
" 14: 7,\n",
" 16: 8}\n"
]
}
],
"source": [
"instance = df['cellinstance']['instance']\n",
"instance_to_index = dict(zip(instance.values, instance.index))\n",
"pprint(instance_to_index, width=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The code for plotting will look exactly the same, except that when we set `img[row, col]` we first need to use the `instance_to_index` mapping to determine the appropriate index:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.colorbar.Colorbar at 0x152844ffb4f0>"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAEKCAYAAAB0cRxpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdvUlEQVR4nO3df5BlZX3n8fdnevjhgoICRhcYZlxGEVjE0IIx6KqIO25ZjFpQQizBDWaKclGzrlESEkQ2poLWwsZAGWeFFVADBCTO6hDEYERLGWdGhx/DgDagYaZQGUBkguDM9Gf/uKfxzqX73tN9T/e55/bnVXUq957z3NPfO9GPz3Oe85wr20RExNQW1F1ARMSgS1BGRPSQoIyI6CFBGRHRQ4IyIqKHBGVERA+1BaWkPSV9X9LtkjZK+nhdtUTE4JK0TNK9ksYknTPJ8T0kXVMcXyNpcduxoyR9r8iYOyXtWew/png/JunTktSthjp7lE8Db7T9CuBoYJmkV9dYT0QMGEkjwKXAW4DDgdMkHd7R7EzgMduHAhcDFxafXQh8ATjL9hHA64HtxWc+A/wRsLTYlnWro7agdMu24u1uxZa73yOi3bHAmO37bf8GuBpY3tFmOXBF8fo64ISih/hm4A7btwPYfsT2TkkvBp5n+za3VtxcCbytWxELK/s6M1D8r8V64FDgUttrJmmzAlgBsNdeex1z2GGHzW2REfPI+vXrt9o+oJ9zvEJ7+Ql2lmr7AE9vBJ5q27XS9sq29wcCD7a93wwc13GaZ9rY3iHpcWA/4KWAJd0EHABcbfuTRfvNHec8sFudtQal7Z3A0ZL2BW6QdKTtuzrarARWAoyOjnrdunVzX2jEPCHpp/2e4wl28pccUqrtu/jRU7ZH+/2bU1gIHA+8CngS+GdJ64HHp3uigZj1tv1L4Jv0uE4QEc2wYEG5rYQtwMFt7w8q9k3aprguuQ/wCK2e4q22t9p+ElgN/G7R/qAe59z1+5QqdRZIOqDoSSLpOcCJwD111RMR1ZBg4cJyWwlrgaWSlkjaHTgVWNXRZhVwRvH6ZOCW4trjTcB/lPTvigD9T8Ddth8CfiXp1cW1zNOBr3Qros6h94uBK4rrlAuAa21/tcZ6IqIConRvsafimuPZtEJvBLjc9kZJFwDrbK8CLgOukjQGPEorTLH9mKSLaIWtgdW2v1ac+n3A54HnADcW25RqC0rbdwCvrOvvR8TsWdD1rsTpsb2a1rC5fd95ba+fAk6Z4rNfoHWLUOf+dcCRZWuodTInIoaQqutRDooEZURUqsqh96BIUEZEtdKjjIjoTpSe0W6MIfs6EVG79CgjIrrLNcqIiBJ6PLWscRKUEVGtDL0jIrqbWMI4TIbs60RE3XKNMiKihARlREQ3uUYZEdFdht4REb2kRxkR0Z2AhSN1V1GtBGVEVCs9yoiI7nKNMiKiJ7GgykecD4AEZURUS6CRBGVExJQkWLBwuMbeCcqIqFx6lBER3Ugo1ygjIrpLjzIioguJzHpHRPSSoXdERDcSC3YbrjWMCcqIqJSG8D7K2m52knSwpG9KulvSRkkfrKuWiKiWFqjUVupc0jJJ90oak3TOJMf3kHRNcXyNpMXF/sWSfi1pQ7H9Xdtn/qU458SxF3aroc4e5Q7gf9j+gaTnAusl3Wz77hprioh+qbprlJJGgEuBE4HNwFpJqzpy4kzgMduHSjoVuBB4Z3HsPttHT3H6d9leV6aO2nqUth+y/YPi9RPAJuDAuuqJiKoIjZTbSjgWGLN9v+3fAFcDyzvaLAeuKF5fB5ygin8vdyDWGRVd5VcCa2ouJSL6NHF7UJmthAOBB9veb+bZHapn2tjeATwO7FccWyLph5K+Jem1HZ/7v8Ww+y96BWvtkzmS9gauB/7Y9q8mOb4CWAGwaNGiOa4uIqZNsGC30n2w/SW1D39X2l5ZUSUPAYtsPyLpGOAfJR1R5My7bG8pLvtdD7wbuHKqE9UalJJ2o1XkF21/ebI2xT/aSoDR0VHPYXkRMUPTuEa51fZol+NbgIPb3h9U7JuszWZJC4F9gEdsG3gawPZ6SfcBLwXW2d5S7H9C0pdoDfGnDMo6Z70FXAZssn1RXXVERMVU6TXKtcBSSUsk7Q6cCqzqaLMKOKN4fTJwi21LOqCYDELSS4ClwP2SFkrav9i/G/BW4K5uRdTZo/x9Wt3dOyVtKPb9me3V9ZUUEf0SoIoecW57h6SzgZuAEeBy2xslXUCrZ7iKVofrKkljwKO0whTgdcAFkrYD48BZth+VtBdwUxGSI8A3gP/TrY7agtL2d2j9m0bEMKnw9iCAovO0umPfeW2vnwJOmeRz19O6tNe5/9+AY6ZTQ+2TORExXCQxUn4ypxESlBFRuWFbwpigjIhqVTz0HgQJyoioXnqUERFdpEcZEdGLYCSTORERUxMos94REV0IyNA7IqIbZTInIqIbZTInIqKETOZERHSRHmVERA8S7D5cPcrh+jaz5G+/+0DdJTTS+AMXM/7AxXWXETWo8lcYB0F6lFPoDMf29+9/zZK5LqcxJgvG9n0Llvz3uSwn6iAy6z3syvQeJ9okMHdVpvc4/sDFCcuhJ6jowb2DYri+TZ8yxJ656QyxMxwfcqLKn4IYCAnKPiRYW2YSfAnLISZgtwXltoZoTqWzbKahN9/Dsp/AS1gOq3ITOZnMiYj5S+SG84iIroZw1nu4Yn+G+h0+z/fhd0S71s/VZugdHXKbUEQbDd/Tg9KjJEEXUaniwb1ltqZIjzIiqpcbzofTTHuV87032s8qm6zQGVbFypwyW0M0p9IYWDMJvITkEBMJymE23d7hfO9NtptO8CUk54EFKrc1RK1BKelySb+QdFeddbR7/2uW9AzAMm3mozIBmJCcB1Tt0FvSMkn3ShqTdM4kx/eQdE1xfI2kxcX+xZJ+LWlDsf1d22eOkXRn8ZlPS+qa2nVP5nweuAS4suY6nmUiCPN4telpD8KJJYoJx/lGsLCaaJE0AlwKnAhsBtZKWmX77rZmZwKP2T5U0qnAhcA7i2P32T56klN/BvgjYA2wGlgG3DhVHbUGpe1bJ9J/UCUcZy4BOU9NXKOsxrHAmO37ASRdDSwH2oNyOXB+8fo64JJuPURJLwaeZ/u24v2VwNvoEpQDf41S0gpJ6ySte/jhh+suJyLKKD/03n/iv9/FtqLjTAcCD7a931zsm7SN7R3A48B+xbElkn4o6VuSXtvWfnOPc+6i7qF3T7ZXAisBRkdHXXM5EdHLxDXKcrbaHp2lSh4CFtl+RNIxwD9KOmImJxr4oIyIBqpuRnsLcHDb+4OKfZO12SxpIbAP8IhtA08D2F4v6T7gpUX7g3qccxcDP/SOiIZRMZlTZuttLbBU0hJJuwOnAqs62qwCzihenwzcYtuSDigmg5D0EmApcL/th4BfSXp1cS3zdOAr3Yqo+/agvwe+B7xM0mZJZ9ZZT0RUQ1pQauuluOZ4NnATsAm41vZGSRdIOqlodhmwn6Qx4EPAxC1ErwPukLSB1iTPWbYfLY69D/gcMAbcR5eJHKh/1vu0Ov9+zF/+108/81qLPlBjJcOo2h8Xs72a1i087fvOa3v9FHDKJJ+7Hrh+inOuA44sW0OuUca80R6Ok+1PYFak2tuDBkKCMuaFqUJysjYJzAo0aHliGcMV+xGTKBOS/bSPDhUvYRwE6VHGUEvo1aSiJYyDojmRHjGHErB9SI8yojkSdjVqUAiWMVzfJqJCCdoZEkP3PMr0KCOiYtXeRzkIEpQRU8htQjOk6p5HOSiG69tExGAosTyxSYbr20S0SY+wRlpQbmuI9CgjJpGQ7cP0nkfZCAnKGGoTgTedGeyEZAW6/1ZX4yQoI9okJCvSoGF1GQnKmBfK9CwTkhWRYGS4omW4vk1EDwnDuaD506OU9LslPr/d9p0V1hMRw2C+BCXwLVq/V9HtquwSYHGVBUXEEJhHkzlrbb+x24cl3VJxPRHRePNo6N0rJMu2iYh5Zr7+FISko2gNsZ9pb/vLs1RTRDSa0ILhmifu+W0kXQ4cBWwExovdBhKUEfFsYv4Mvdu82vbhs17JgLlj67Ypjx21/95zWEnMFzf/9NEpj514yAvmsJJ+aegmc8rE/vckzaug7BaSZY5HTFe3kCxzfOAM2U9BlKn0Slphea+kOyTdKemO2S6sLmVDMGEZVSkbgo0Ky3n49KDLgHcDd/Lba5RDabrhd8fWbRmGR1+mG343//TRwR+GSzBkkzllIv1h26tsP2D7pxPbrFc2x2baQ0zPMmZqpj3Ewe9ZDt+vMJap9IeSviTpNEnvmNiq+OOSlhVD+jFJ51RxzogYABUOvXvlhKQ9JF1THF8jaXHH8UWStkn6cNu+nxSXETdIWterhjL94+cATwNvbtvX9+1BkkaAS4ETgc3AWkmrbN/dz3kjomaislnvkjlxJvCY7UMlnQpcCLyz7fhFwI2TnP4NtreWqaNnUNr+r2VONAPHAmO27weQdDWwHJjzoOx3+JxrlTFd/Q6fB/taZaVLGMvkxHLg/OL1dcAlkmTbkt4GPAD8Wz9F9Pw2kq6QtG/b++cXN6H360Dgwbb3m4t9nX9/haR1ktY9/PDDFfzZiJh15Yfe+0/897vYVnScqUxOPNPG9g7gcWA/SXsDHwU+PkmFBr4uaf0kf/NZygy9j7L9y2fObj8m6ZUlPlcJ2yuBlQCjo6Oejb9x1P5799WrTG8ypuvEQ17QV69ycHuTAGK8/KNut9oenaVCzgcutr1Nz74UcLztLZJeCNws6R7bt051ojLfZoGk59t+DEDSC0p+rpctwMFt7w8q9kVEw5mRqk5VJicm2myWtBDYB3gEOA44WdIngX2BcUlP2b7E9hYA27+QdAOtIX5fQfm/aN1w/g/F+1OAT5T4XC9rgaWSltD6oqcCf1DBeSOiRkbYlV2jLJMTq4AzgO8BJwO32Dbw2okGks4Httm+RNJewALbTxSv3wxc0K2Int/G9pXAO4CfF9s7bF9V6it2P+8O4GzgJmATcK3tjf2ed6ZmOnzOsDtmaqbD58EedreYBaW2nueZIickXSDppKLZZbSuSY4BHwJ63Wr4O8B3JN0OfB/4mu1/6vYBtYK3GUZHR71uXc9bnvoynWuVCcmownSuVc52SEpa3+81w2NGD/d311xZqu2eC1/V99+bC1NGuqQf9PpwmTZNUzb8EpJRlbLh14Se5AQzUmprim7XKF/e4+EXonXRdOj0mgVPSEbVes2CNyokLcY9XGu9u32bw0p8fmdVhQyahGHMtSaFYS/jpVZHN0e338wZugdfRMRcqHTWeyAMV/84ImpnKr2PciAkKCOiYsPXoyyz1vv9kp4/F8VERPMZsdMLS21NUSb2f4fWo42uLZ4LN1y/GhQRlRu224PKrMz5c2Aprbvf3wP8WNJfSfoPs1xbRDRSa+hdZmuKUn3f4rluPwN+BuwAng9cJ+lm2x+ZzQLj2cZ//Kkpjy1Y+idzWEnE5MosT2ySnkEp6YPA6cBW4HPAn9jeLmkB8GMgQTmHuoVk+/EEZtTGNKq3WEaZHuULaD0IY5f7Km2PS3rr7JQVk+kVkhGDwIjxBl1/LKPMT0F8rMuxTdWWE1OZbkimZxn1mV9LGGNApCcZTTNsQ+/h+jbxLAnZmGsTQ+8yW1OkRxkRlRu2HmWCcsClRxjNo/l3e1BExHQYGG/QLyeUMVyxP4Qyax2NY7NjvNzWFOlRDrkEbcy1Vo+y7iqqlR5lAyTsomnG7VJbU6RHOcQSsFEHO9cooybTDb2EZNRpp8ttTZEeZYNMhF+vW4YSklEnQ6MmaspIUDbQZIGZcIxBMmxD7wRlgyUcYxC1rlHWXUW1EpQRUbFmzWiXUctkjqRTJG2UNC5ptI4aImJ2TKzMqer2oOK3uu6VNCbpnEmO7yHpmuL4GkmLO44vkrRN0ofLnrNTXbPedwHvAG6t6e9HxCwad7mtF0kjwKXAW4DDgdMkHd7R7EzgMduHAhcDF3Ycvwi4cZrn3EUtQWl7k+176/jbETG7bKpcwngsMGb7ftu/Aa4Glne0WQ5cUby+Djhh4tdiJb0NeADYOM1z7mLg76OUtELSOknrHn744brLiYgSpjH03n/iv9/FtqLjVAcCD7a931zsm7SN7R3A48B+kvYGPgp8fAbn3MWsTeZI+gbwokkOnWv7K2XPY3slsBJgdHR0uK4QDzn/7HMA6EXvrbmSZnn333z7mddXffC1NVYyM57eZM5W27M1T3E+cLHtbUUHc8ZmLShtv2m2zh2DayIcp9qX0Jxae0B27mtaYFZ4e9AW4OC29wcV+yZrs1nSQmAf4BHgOOBkSZ8E9gXGJT0FrC9xzl3k9qCozGQhOVmbhOWuJgvIqdo0ITAN7Kzu9qC1wFJJS2iF2anAH3S0WQWcAXwPOBm4xbaBZ/6xJJ0PbLN9SRGmvc65i7puD3q7pM3A7wFfk3RTHXVEdcqE5EzaDrsyIdlP+1pUOJlTXHM8G7gJ2ARca3ujpAsknVQ0u4zWNckx4ENA19t9pjpnt8/U0qO0fQNwQx1/O6o3k+BLz3J4Vf2Ec9urgdUd+85re/0UcEqPc5zf65zdDPysdwy2fnqH871nOdPeYRN6lVXdRzkoco0yIiqV38yJiIEw0L1Kw/h4ua0p0qOMqMFAB13fzHiTxtUlpEcZUYN+b/MZ5NuEbNi+c7zU1hTpUUZEpQxD16NMUEZE5TKZE9Gmn3sh5/t9lDMdPg/ysBvAbl2jLLM1RYIy+jaTwJvvITlh0ENvphKUEZOYTvAlJGeuCcFqD19Q5hplVEYvem/P1TYJyWebCL9utww1ISAnGBo1o11GgjIq1R6EeR7l9LSH4bv/5tuNCsdduFm9xTISlDFrEpAz19iQLCQoIyK6GMa13gnKiKiW06OMiOjKOJM5ERHd2M16MlAZCcqIqFyG3hERPSQoIyK6sDPrHRHRU3qUERFd2Gb7juGazUlQRkSlcsN5REQJGXpHRHSTlTkREd15CH+FMUEZEdUyjO8crqCs5Qnnkj4l6R5Jd0i6QdK+ddQREdWzYeeOnaW2MiQtk3SvpDFJ50xyfA9J1xTH10haXOw/VtKGYrtd0tvbPvMTSXcWx9b1qqGun4K4GTjS9lHAj4A/ramOiKhcdT8uJmkEuBR4C3A4cJqkwzuanQk8ZvtQ4GLgwmL/XcCo7aOBZcBnJbWPot9g+2jbo73qqCUobX/d9o7i7W3AQXXUERHVq/g3c44Fxmzfb/s3wNXA8o42y4EritfXASdIku0n23JmT1p3Ls3IIPy42B8CN9ZdRERUZ3ynS20lHAg82PZ+c7Fv0jZFMD4O7Acg6ThJG4E7gbPagtPA1yWtl7SiVxGzNpkj6RvAiyY5dK7trxRtzgV2AF/scp4VwAqARYsWzUKlEVGp6d0etH/HNcKVtldWVoq9BjhC0suBKyTdaPsp4HjbWyS9ELhZ0j22b53qPLMWlLbf1O24pPcAbwVOsKe+jb/4R1sJMDo6OlxTaRFDyLj0RA2wtcc1wi3AwW3vDyr2TdZmc3ENch/gkV1qsjdJ2gYcCayzvaXY/wtJN9Aa4k8ZlHXNei8DPgKcZPvJOmqIiFlS7TXKtcBSSUsk7Q6cCqzqaLMKOKN4fTJwi20Xn1kIIOkQ4DDgJ5L2kvTcYv9ewJtpTfxMqa77KC8B9qDV5QW4zfZZNdUSERUy1d1HaXuHpLOBm4AR4HLbGyVdQKtnuAq4DLhK0hjwKK0wBTgeOEfSdmAceJ/trZJeAtxQZM9C4Eu2/6lbHbUEZTGNHxHDqOIljLZXA6s79p3X9vop4JRJPncVcNUk++8HXjGdGrIyJyIqlSWMERG9GJygjIiYmg07tufBvRERXZjxIfu92gRlRFTK5HmUEV399T//qOvxc0546RxVErUZwgf3DsJa7xgSvUKybJtovgpvOB8ICcqoxHQCMGE53GzwTpfamiJBGX2bSfAlLIeYzY4d46W2psg1yoiolAEP2ax3epTRl356hulVDqkhHHqnRxkRFXNW5kREdFX0KIdJgjIiqmXw9tIP7m2EBGVEVMoZekdE9DCEQ+/Mekdf+lmSmOWMw8vjLrU1RYIy+jaTwEtIDrnx8XJbQyQooxLTCb6E5JBzuXsomzQ8T1BGZcoEYEJy+Nkwvn281NYUmcyJSk0EYeeqmwTk/NKk649lJChjViQY5zEbdjant1hGgjIiqpUfF4uI6K1JEzVlJCgjolrpUUZEdGc7a70jInrJ0DsiopshHHrXcsO5pP8p6Q5JGyR9XdK/r6OOiJgN1a7MkbRM0r2SxiSdM8nxPSRdUxxfI2lxsf/YImM2SLpd0tvLnrNTXStzPmX7KNtHA18FzqupjoiomqlsrbekEeBS4C3A4cBpkg7vaHYm8JjtQ4GLgQuL/XcBo0XOLAM+K2lhyXPuopagtP2rtrd70fqnjYghUPESxmOBMdv32/4NcDWwvKPNcuCK4vV1wAmSZPtJ2zuK/Xvy25wpc85d1HaNUtIngNOBx4E3dGm3AlhRvH1a0l1zUF5Z+wNb6y6izaDVA4NXU+rp7mX9nuDXT/7rTT/8/n/bv2TzPSWta3u/0vbKtvcHAg+2vd8MHNdxjmfa2N4h6XFgP2CrpOOAy4FDgHcXx8uccxezFpSSvgG8aJJD59r+iu1zgXMl/SlwNvCxyc5T/KOtLM65zvbobNU8Xamnt0GrKfV01xFaM2J7WRW1VMH2GuAISS8HrpB040zOM2tBaftNJZt+EVjNFEEZEfPaFuDgtvcHFfsma7NZ0kJgH+CR9ga2N0naBhxZ8py7qGvWe2nb2+XAPXXUEREDby2wVNISSbsDpwKrOtqsAs4oXp8M3GLbxWcWAkg6BDgM+EnJc+6irmuUfy3pZcA48FPgrJKfW9m7yZxKPb0NWk2pp7uBqqe4png2cBMwAlxue6OkC4B1tlcBlwFXSRoDHqUVfADHA+dI2k4ra95neyvAZOfsVofsTDhHRHSTJ5xHRPSQoIyI6KFxQTloyx8lfUrSPUVNN0jat+Z6TpG0UdK4pNpuO5nuErE5qOdySb8YlPtwJR0s6ZuS7i7+//XBmuvZU9L3i6V+GyV9vM56Bk3jrlFKet7Eyh5JHwAOt112Mmg26nkzrVm2HZIuBLD90RrreTmtC9efBT5su+/74mZQwwjwI+BEWjfzrgVOs333XNfSVtPrgG3AlbaPrKuOtnpeDLzY9g8kPRdYD7ytrn8jSQL2sr1N0m7Ad4AP2r6tjnoGTeN6lIO2/NH219uWSd1G656sOuvZZPveOmtgBkvEZpvtW2nNiA4E2w/Z/kHx+glgE60VJnXVY9vbire7FVuzelGzqHFBCa3lj5IeBN7FYD1Q4w+BGd35P2QmWyJWWwgMuuJpN68E1tRcx4ikDcAvgJuLVS3BgAalpG9IumuSbTmA7XNtH0xrVc/ZdddTtDkX2FHUVHs90QyS9gauB/64Y7Q052zvLJ60cxBwrKTaL1EMioF8cO+gLX/sVY+k9wBvBU7wHFz0nca/T12mvURsPiquBV4PfNH2l+uuZ4LtX0r6Jq1Hkw3E5FfdBrJH2c2gLX+UtAz4CHCS7SfrrGWATHuJ2HxTTJ5cBmyyfdEA1HPAxB0bkp5DayIuS4sLTZz1vp7Wo6CeWf5ou7beSrFsag9+uwj/tppn4d8O/C1wAPBLYIPt/1xDHf8F+N/8donYJ+a6ho56/h54Pa3Hmv0c+Jjty2qs53jg28CdtP6zDPBntlfXVM9RtJ7pOEKrA3Wt7QvqqGUQNS4oIyLmWuOG3hERcy1BGRHRQ4IyIqKHBGVERA8JyoiIHhKUERE9JCijb5IWS/p1sU64ivN9U9K2Oh8TF9EuQRlVua9YJ9w3228A5vzxcBFTSVBGV5JeVTyUeE9JexUPde35sARJpxefu13SVcW+z0v6jKTbJN0v6fXFA3U3Sfr8rH+ZiBkayIdixOCwvVbSKuAvgecAX7Dd9UEJko4A/hx4je2tkl7Qdvj5wO8BJ9Fa//37wHuBtZKOtr1hFr5GRF8SlFHGBbQedPEU8IES7d8I/MPET4Pabn9g7v8rfnP5TuDntu8EkLQRWAxsqLDuiEpk6B1l7AfsDTwX2LPPcz1d/N/xttcT7/M/3DGQEpRRxmeBv6D1/M8LS7S/BThF0n4AHUPviMbJ/4JHV5JOB7bb/lLxo2HflfRG27dM9RnbGyV9AviWpJ3AD4H3zE3FEdXLY9aib8Vvvny1yl83lPQv1PQrkhGdMvSOKuwE9qnyhnPgJcD2Ks4X0a/0KCMiekiPMiKihwRlREQPCcqIiB4SlBERPfx/UZfxW7KuEjkAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"img[:] = np.nan\n",
"with openmc.lib.run_in_memory():\n",
" for row, y in enumerate(np.linspace(ymin, ymax, resolution[0])):\n",
" for col, x in enumerate(np.linspace(xmin, xmax, resolution[1])):\n",
" try:\n",
" # For each (x, y, z) point, determine the cell and distribcell index\n",
" cell, distribcell_index = openmc.lib.find_cell((x, y, 0.))\n",
" except openmc.exceptions.GeometryError:\n",
" # If a point appears outside the geometry, you'll get a GeometryError exception.\n",
" # These lines catch the exception and assign a \"null\" value for the array\n",
" continue\n",
"\n",
" if cell.id == fuel_cell.id:\n",
" # When the cell ID matches, we set the corresponding pixel in the image using the\n",
" # distribcell index. NOTE: In this case, we need to use the dictionary to find the\n",
" # proper index in the flux array\n",
" index = instance_to_index.get(distribcell_index)\n",
" if index is not None:\n",
" # If the distribcell instance was specified in our filter, get the flux value\n",
" img[row, col] = flux_inst[index]\n",
"\n",
"# Plot the image using the same options as before. Notably, this ensures that the colors are consistent.\n",
"plt.imshow(img, **options)\n",
"plt.xlabel('x [cm]')\n",
"plt.ylabel('y [cm]')\n",
"plt.colorbar()"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python 3",
"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.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment