MSRE/msre_docker/MSRE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a28f447d-a699-4495-bf71-bace1a4c0833",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Notebook for running MSRE calculations directly from a CAD drawing\n",
    "This notebook show an example of running computations of a model of the Moten Salt Reactor Experiment of MSRE in short.\n",
    "The model itself has been generated from tehe original drawings from Oak Ridge National lab \n",
    "\n",
    "The CAD-model is available on github at https://github.com/openmsr/msre, which in turn is generated from a long list of documents which have been compiled at https://github.com/openmsr/msr-archive\n",
    "\n",
    "The simulation backend is run using the Open Source Monte Carlo particle transport code OpenMC (https://openmc.org), through its' python interface.\n",
    "\n",
    "**Important: If you want your work to be available after you shutdown the docker, you must copy your notebooks to a location mounted on your local machine.**\n",
    "\n",
    "If you started the docker using the supplied ```run_docker.sh```-script, the ```notebooks```-directory has been mounted like this."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb624881-ed0a-4962-b528-fb1ee5141f35",
   "metadata": {},
   "source": [
    "## The (obvious) 1st step is to import the OpenMC python interface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85306562-ad94-4072-9cd4-1734fc7b0ab0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import openmc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51730852-f7a5-4491-8b1d-de7609dce96f",
   "metadata": {},
   "source": [
    "Next we define a set of materials objects that form the core of the MSRE, graphite,  hastelloy N / inor-8, inconel, the fuel salt, and helium. Lastly these are exported to am OpenMC-xml control file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6df036db-7b00-4a73-a4e0-dd9a407d7395",
   "metadata": {},
   "outputs": [],
   "source": [
    "graphite=openmc.Material(name='graphite')\n",
    "graphite.add_element('C',1.0,'ao')\n",
    "graphite.set_density('g/cc',2.26)\n",
    "\n",
    "#Hastelloy N / INOR-8 nominal material composition from\n",
    "#ORNL-TM-4189\n",
    "inor=openmc.Material(name='inor')\n",
    "inor.add_element('Ni',0.72)\n",
    "inor.add_element('Mo',0.16)\n",
    "inor.add_element('Cr',0.07)\n",
    "inor.add_element('Fe',0.05)\n",
    "inor.set_density('g/cc',9)\n",
    "\n",
    "# LiF,BeF2,UF4,ZrF4 [0.67,0.23,0.05,0.0079] mol % @ 33% enrichment \n",
    "molar_comp={'LiF':0.67,'BeF2':0.23, 'ZrF4':0.05, 'UF4':0.0079}\n",
    "enrichment=0.3333\n",
    "salt=openmc.Material(name='salt')\n",
    "salt.add_element('F',molar_comp['LiF']*1/2+molar_comp['BeF2']*2/3+molar_comp['ZrF4']*4/5+molar_comp['UF4']*4/5,'ao')\n",
    "salt.add_nuclide('Li7',molar_comp['LiF']*1/2,'ao')\n",
    "salt.add_element('Be',molar_comp['BeF2']*1/3,'ao')\n",
    "salt.add_element('Zr',molar_comp['ZrF4']*1/5,'ao')\n",
    "salt.add_nuclide('U235',enrichment*molar_comp['UF4']*1/5,'ao')\n",
    "salt.add_nuclide('U238',(1-enrichment)*molar_comp['UF4']*1/5,'ao')\n",
    "salt.set_density('g/cc',2.2)\n",
    "\n",
    "# The natural isotopes have been used for this alloy\n",
    "# The density is set to that of Ni.\n",
    "inconel=openmc.Material(name='inconel')\n",
    "inconel.add_element('Ni',0.72,'ao')\n",
    "inconel.add_element('Cr',0.20,'ao')\n",
    "inconel.add_element('Fe',0.08,'ao')\n",
    "inconel.set_density('g/cc',8.9)\n",
    "\n",
    "helium=openmc.Material(name='helium')\n",
    "helium.add_nuclide('He4',1.0,'ao')\n",
    "helium.set_density('g/cc',1.0e-4)\n",
    "\n",
    "materials=openmc.Materials([helium,salt,graphite,inconel,inor])\n",
    "materials.export_to_xml()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e21759c-4d3f-44f4-8b61-eb1a519d7d68",
   "metadata": {},
   "source": [
    "As a control we can inspect the materials object. Notice how OpenMC has in the revant cases expanded our material definition to consist the naturally occurring  isotope concentrations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82b7ea1e-28d6-450d-9bad-e5f5262c5272",
   "metadata": {},
   "outputs": [],
   "source": [
    "materials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a5f8f08-ec15-4ab4-b091-c3e0b2147490",
   "metadata": {},
   "outputs": [],
   "source": [
    "#geometry\n",
    "h5m_filepath=\"../msre/msre_simple.h5m\"\n",
    "dag_univ = openmc.DAGMCUniverse(h5m_filepath)\n",
    "geom = openmc.Geometry(root=dag_univ)\n",
    "geom.export_to_xml()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be73fe16-ea81-4efa-939f-3934bc5b77bc",
   "metadata": {},
   "source": [
    "We can now plot our geometry to verify that this is in fact the geometry we want. We plot two slices (xz and xy) through the centre of the MSRE core, and color the geometry by constituent material."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dde16460-008e-4e6f-ac81-32e82dd066d4",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xwidth = 350\n",
    "yheight = 350\n",
    "material_colors={salt:'red', inor:'lightblue', inconel:'blue',helium:'white',graphite:'gray'}\n",
    "#xz plot\n",
    "p1 = openmc.Plot()\n",
    "p1.background='white'\n",
    "p1.basis = 'xz'\n",
    "p1.width = (xwidth,yheight)\n",
    "p1.origin=(0,0,125)\n",
    "p1.pixels = (800, 800)\n",
    "p1.color_by = 'material'\n",
    "p1.colors=material_colors\n",
    "#xy plot\n",
    "p2 = openmc.Plot()\n",
    "p2.background='white'\n",
    "p2.basis='xy'\n",
    "p2.width=(xwidth, yheight)\n",
    "p2.pixels = (800,800)\n",
    "p2.origin=(0,0,100)\n",
    "p2.color_by='material'\n",
    "p2.colors=material_colors\n",
    "\n",
    "plots=openmc.Plots([p1,p2])\n",
    "openmc.plot_inline(plots)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d02c4c6-2f47-44ac-989c-1ecaf9441450",
   "metadata": {},
   "source": [
    "Now we need to define som settings for our calculations.\n",
    "\n",
    "First of all - we need to some neutrons to kick-start or chain reaction. In OpenMC this is doen by defining a source region. Here this is simply defined as being a region that encloses the MSRE core.\n",
    "\n",
    "Next we define some settings for the Monte Carlo-computation, such as how many particles we would initially run with. \n",
    "\n",
    "After the members of the settings python object have been filled to our desires, we export this to a settings xml-file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a17a765-10dd-4666-9f66-2d7f7931bb3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a neutron source for kick-starting\n",
    "source_volume=openmc.stats.Box([-125,-125,0],[125,125,500], only_fissionable=True)\n",
    "source = openmc.Source(space=source_volume)\n",
    "source.angle=openmc.stats.Isotropic()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee524ba8-d65f-48ab-8024-e05f3409a6ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Finally we build a settings object for OpenMC where we define parameters for the run.\n",
    "settings = openmc.Settings()\n",
    "settings.source = source\n",
    "settings.batches = 20\n",
    "settings.inactive = 5\n",
    "settings.particles = 20000\n",
    "settings.export_to_xml()\n",
    "\n",
    "openmc.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfadbeb7-077b-42c6-8714-c92829536998",
   "metadata": {},
   "source": [
    "Now that we have a running model let's try to do some more useful work with and extract some data from the model. To do this we need to specify what information we want to extract before \n",
    "starting simulations. In many Monte Carlo particle transport codes, we add objects known as tallies to our models. In this respect OpenMC is no different.\n",
    "\n",
    "We will add tallies to monitor the neutron flux, and the fission sites in volumes along the geomtrical slices through our reactor that we plotted earlier.\n",
    "\n",
    "A tally needs to know what to measure and where to measure that. In OpenMC the \"where\" is known as a filter and the \"what\" is known as a score.\n",
    "In or case we'd like to spatially resolve the flux so we first generate mesh object as filters and then add that to tally objects. In addtion we assign a list of scores to the score-member of the tallies. Lastly (as always) we export this to an xml-file which will be read by OpenMC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3383b0ba-4cec-4733-a271-a032079d0b00",
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh1=openmc.RegularMesh()\n",
    "mesh1.dimension =  [400,400,1]\n",
    "mesh1.lower_left = [-125, -125, 90]\n",
    "mesh1.upper_right = [125, 125, 110]\n",
    "mesh1_filter = openmc.MeshFilter(mesh1)\n",
    "\n",
    "mesh2=openmc.RegularMesh()\n",
    "mesh2.dimension =  [400,1,400]\n",
    "mesh2.lower_left = [-175, -10, -50]\n",
    "mesh2.upper_right = [175,  10, 400]\n",
    "mesh2_filter = openmc.MeshFilter(mesh2)\n",
    "\n",
    "t1 = openmc.Tally(name='flux1')\n",
    "t1.filters =[mesh1_filter]\n",
    "t1.scores = ['flux','fission']\n",
    "\n",
    "t2 = openmc.Tally(name='flux2')\n",
    "t2.filters =[mesh2_filter]\n",
    "t2.scores = ['flux','fission']\n",
    "\n",
    "tallies=openmc.Tallies([t1,t2])\n",
    "tallies.export_to_xml()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e62c1ebc-75bf-44eb-8dde-90fef06e0030",
   "metadata": {},
   "source": [
    "We have to re-run our simulation after generating the ```tallies.xml``` file.\n",
    "\n",
    "Note that we forcibly remove the _old_ datafiles first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e624e13-f734-4e58-9307-f3b4d56151a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.system(\"rm -f summary.h5 statepoint.20.h5\")\n",
    "openmc.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "430aff88-5238-4a68-8d7a-2d0de99eb77d",
   "metadata": {},
   "source": [
    "After the run has finished the data we are after resides in the \"statepoint\" file that OpenMC saves.\n",
    "In the below code, we will open that and extract the mean values for neutron flux."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7bebead-b04a-4a38-8802-0d01239e2c57",
   "metadata": {},
   "outputs": [],
   "source": [
    "sp=openmc.StatePoint('statepoint.20.h5')\n",
    "\n",
    "tl1=sp.get_tally(name='flux1')\n",
    "tl2=sp.get_tally(name='flux2')\n",
    "\n",
    "flux1=tl1.get_slice(scores=['flux'])\n",
    "flux2=tl2.get_slice(scores=['flux'])\n",
    "\n",
    "flux1.mean.shape=(400,400)\n",
    "flux2.mean.shape=(400,400)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8205bf93-0351-458e-9d44-f8b43b0509ac",
   "metadata": {},
   "source": [
    "The last 2 lines are necessary to reshape the flux maps into a 400x400 grid.\n",
    "\n",
    "In the end we plot the maps using matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84b1be81-188c-467a-9da5-74e7054792d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "fig,(ax1,ax2)=plt.subplots(ncols=2,figsize=(20,16), constrained_layout=True)\n",
    "\n",
    "ax1.set_xticks(np.arange(0,401,399/4))\n",
    "ax1.set_xticklabels(np.arange(-175,176,350./4))\n",
    "ax2.set_xticks(np.arange(0,400,399/4))\n",
    "ax2.set_xticklabels(np.arange(-175,176,350./4))\n",
    "ax1.set_yticks(np.arange(0,400,399/4))\n",
    "ax1.set_yticklabels(np.arange(-175,176,350./4))\n",
    "ax2.set_yticks(np.arange(0,400,399/4))\n",
    "ax2.set_yticklabels(np.arange(-175,176,350./4))\n",
    "ax1.set_xlabel('X / cm')\n",
    "ax1.set_ylabel('Y / cm')\n",
    "\n",
    "ax2.set_xlabel('X / cm')\n",
    "ax2.set_ylabel('Z / cm')\n",
    "im1=ax1.imshow(flux1.mean)\n",
    "fig.colorbar(im1,ax=ax1,shrink=0.4)\n",
    "im2=ax2.imshow(flux2.mean)\n",
    "fig.colorbar(im2,ax=ax2,shrink=0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1099e99a-e4a8-46c8-b3e4-00259de64c82",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Similarly we can plot fission reactions:\n",
    "\n",
    "fission1=tl1.get_slice(scores=['fission'])\n",
    "fission2=tl2.get_slice(scores=['fission'])\n",
    "\n",
    "fission1.mean.shape=(400,400)\n",
    "fission2.mean.shape=(400,400)\n",
    "\n",
    "fig,(ax1,ax2)=plt.subplots(ncols=2,figsize=(20,16), constrained_layout=True)\n",
    "\n",
    "ax1.set_xticks(np.arange(0,401,399/4))\n",
    "ax1.set_xticklabels(np.arange(-175,176,350./4))\n",
    "ax2.set_xticks(np.arange(0,400,399/4))\n",
    "ax2.set_xticklabels(np.arange(-175,176,350./4))\n",
    "ax1.set_yticks(np.arange(0,400,399/4))\n",
    "ax1.set_yticklabels(np.arange(-175,176,350./4))\n",
    "ax2.set_yticks(np.arange(0,400,399/4))\n",
    "ax2.set_yticklabels(np.arange(-175,176,350./4))\n",
    "ax1.set_xlabel('X / cm')\n",
    "ax1.set_ylabel('Y / cm')\n",
    "\n",
    "ax2.set_xlabel('X / cm')\n",
    "ax2.set_ylabel('Z / cm')\n",
    "im1=ax1.imshow(fission1.mean)\n",
    "fig.colorbar(im1,ax=ax1,shrink=0.4)\n",
    "im2=ax2.imshow(fission2.mean)\n",
    "fig.colorbar(im2,ax=ax2,shrink=0.4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd6f56f4-6ade-4e4f-b751-e67df40b3d8c",
   "metadata": {},
   "source": [
    "Suppose we now would like to see what the energy spectrum of the neutrons generated in our reactor is. To explore this we will add another talliy to our simulation. This time however, instead of a spatial regular mesh, the tally will have an energy filter. Furthermore, we restrict the tally to neutron flux and fission events within the fuel salt, by means of a material filter. Unfortunately to fill the new tally we have to re-run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e440a7a0-3f4d-4067-ba47-09245d82400e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#define a lograrithmic binning from 1keV to 10 MeV\n",
    "ef=energyrange=openmc.EnergyFilter(np.logspace(3,7,200))\n",
    "te = openmc.Tally(name='energy')\n",
    "\n",
    "sf = openmc.MaterialFilter(salt)\n",
    "te.filters =[ef,sf]\n",
    "te.scores = ['flux','fission']\n",
    "\n",
    "tallies.append(te)\n",
    "tallies.export_to_xml()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa856b84-144c-49fc-8184-1c35eedc0449",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.system(\"rm -f summary.h5 statepoint.20.h5\")\n",
    "openmc.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f431317-75ca-4960-ba3a-fae1c13e643b",
   "metadata": {},
   "outputs": [],
   "source": [
    "sp=openmc.StatePoint('statepoint.20.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4920fa40-9468-4fdc-a7e6-4cfe1a88612d",
   "metadata": {},
   "outputs": [],
   "source": [
    "tl=sp.get_tally(name='energy')\n",
    "\n",
    "e1=tl.get_slice(scores=['flux'])\n",
    "e2=tl.get_slice(scores=['fission'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22a6b009-f2d9-441b-ac98-4cb71c6c99a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax=plt.subplots(figsize=(7,3.5),squeeze=True)\n",
    "ax.plot(ef.values[:199],e1.mean[:,0,0])\n",
    "ax.plot(ef.values[:199],e2.mean[:,0,0])\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9aefede-8fdf-4494-b0c0-d19e41fb6c4f",
   "metadata": {},
   "source": [
    "A reactor with a k_{eff} significantly higher than 1 is likely not what you want. Therefore we need to modify our initial model to also include a set of control rods. In the MSRE these rods consisted of three sets of cylindrical elements made from a Al_2O_3/Gd_2O_3-mixture. These absorb neutrons to \"dampen\" the nuclear process - something also known as poisoning.\n",
    "The cylindrical elements were stacked to form the control rods (~=80''), which can be inserted into the reactor core in 3 of the 4 voids visible in the centre of the XZ-geometry of the reactor.\n",
    "\n",
    "To run our reactor model with control-rods we will now simply point openmc at a different geometry-file, and append the missing materials to the materials list. To visualize we also add tallies to track the absorption."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57d11f47-2c4a-4549-8944-ef1aacf81c7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#geometry\n",
    "os.system('rm -f materials.xml geometry.xml tallies.xml')\n",
    "h5m_filepath=\"../msre/msre_control_in.h5m\"\n",
    "dag_univ = openmc.DAGMCUniverse(h5m_filepath)\n",
    "geom = openmc.Geometry(root=dag_univ)\n",
    "geom.export_to_xml()\n",
    "\n",
    "b_w_conc={'Al2O3':0.3,'Gd2O3':0.7} # Wt. (D. Shen et.al. Nucl. Sc. & Eng., v. 195, pp. 825, 2021)\n",
    "A_w={'Al':26.9815385,'Gd':157.25, 'O':15.99} #g/mol, (Webelements.com)\n",
    "rho={'Al2O3':3.987,'Gd2O3':7.07} # g /cc, (Wikipedia: Aluminium_oxide & Gadolinium(III)_oxide)\n",
    "b_mol_w={'Al2O3':2*A_w['Al']+ 3*A_w['O'],'Gd2O3':2*A_w['Gd']+3*A_w['O']}\n",
    "bush_mol_comp={'Al2O3':(b_w_conc['Al2O3']/b_mol_w['Al2O3'])/( b_w_conc['Al2O3']/b_mol_w['Al2O3'] + b_w_conc['Gd2O3']/b_mol_w['Gd2O3'] ),\n",
    "                'Gd2O3':(b_w_conc['Gd2O3']/b_mol_w['Gd2O3'])/( b_w_conc['Al2O3']/b_mol_w['Al2O3'] + b_w_conc['Gd2O3']/b_mol_w['Gd2O3'] )}\n",
    "\n",
    "bush=openmc.Material(name='bush')\n",
    "bush.add_element('Al',2/5*bush_mol_comp['Al2O3'],'ao')\n",
    "bush.add_element('Gd',2/5*bush_mol_comp['Gd2O3'],'ao')\n",
    "bush.add_element('O',3/5*bush_mol_comp['Al2O3']+3/5*bush_mol_comp['Gd2O3'],'ao')\n",
    "bush.set_density('g/cc',b_w_conc['Al2O3']*rho['Al2O3'] +b_w_conc['Gd2O3']*rho['Gd2O3'] )\n",
    "\n",
    "t1 = openmc.Tally(name='flux1')\n",
    "t1.filters =[mesh1_filter]\n",
    "t1.scores = ['flux','fission','absorption']\n",
    "\n",
    "t2 = openmc.Tally(name='flux2')\n",
    "t2.filters =[mesh2_filter]\n",
    "t2.scores = ['flux','fission','absorption']\n",
    "\n",
    "tallies=openmc.Tallies([t1,t2])\n",
    "tallies.export_to_xml()\n",
    "\n",
    "\n",
    "materials=openmc.Materials([helium,salt,graphite,inconel,inor,bush])\n",
    "materials.export_to_xml()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea8d4900-f24e-4f0e-b5af-39ad3342bebc",
   "metadata": {},
   "outputs": [],
   "source": [
    "xwidth = 350\n",
    "yheight = 350\n",
    "material_colors={salt:'red', inor:'lightblue', inconel:'blue',helium:'white',graphite:'gray',bush:'purple'}\n",
    "#xz plot\n",
    "p1 = openmc.Plot()\n",
    "p1.background='white'\n",
    "p1.basis = 'xz'\n",
    "p1.width = (xwidth,yheight)\n",
    "p1.origin=(0,0,125)\n",
    "p1.pixels = (800, 800)\n",
    "p1.color_by = 'material'\n",
    "p1.colors=material_colors\n",
    "#xy plot\n",
    "p2 = openmc.Plot()\n",
    "p2.background='white'\n",
    "p2.basis='xy'\n",
    "p2.width=(xwidth, yheight)\n",
    "p2.pixels = (800,800)\n",
    "p2.origin=(0,0,100)\n",
    "p2.color_by='material'\n",
    "p2.colors=material_colors\n",
    "\n",
    "p3 = openmc.Plot()\n",
    "p3.background='white'\n",
    "p3.basis='xy'\n",
    "p3.width=(xwidth/10, yheight/10)\n",
    "p3.pixels = (800,800)\n",
    "p3.origin=(0,0,100)\n",
    "p3.color_by='material'\n",
    "p3.colors=material_colors\n",
    "\n",
    "plots=openmc.Plots([p1,p2,p3])\n",
    "openmc.plot_inline(plots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c024f9cd-3ca8-4f8d-a5d0-60eca95cac3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.system(\"rm -f summary.h5 statepoint.20.h5\")\n",
    "openmc.run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "046eb354-e11f-4dfa-90e8-59768cc4c063",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot maps of the absorption\n",
    "sp=openmc.StatePoint('statepoint.20.h5')\n",
    "\n",
    "tl1=sp.get_tally(name='flux1')\n",
    "tl2=sp.get_tally(name='flux2')\n",
    "\n",
    "abs1=tl1.get_slice(scores=['absorption'])\n",
    "abs2=tl2.get_slice(scores=['absorption'])\n",
    "\n",
    "abs1.mean.shape=(400,400)\n",
    "abs2.mean.shape=(400,400)\n",
    "\n",
    "fig,(ax1,ax2)=plt.subplots(ncols=2,figsize=(20,16), constrained_layout=True)\n",
    "\n",
    "ax1.set_xticks(np.arange(0,401,399/4))\n",
    "ax1.set_xticklabels(np.arange(-175,176,350./4))\n",
    "ax2.set_xticks(np.arange(0,400,399/4))\n",
    "ax2.set_xticklabels(np.arange(-175,176,350./4))\n",
    "ax1.set_yticks(np.arange(0,400,399/4))\n",
    "ax1.set_yticklabels(np.arange(-175,176,350./4))\n",
    "ax2.set_yticks(np.arange(0,400,399/4))\n",
    "ax2.set_yticklabels(np.arange(-175,176,350./4))\n",
    "ax1.set_xlabel('X / cm')\n",
    "ax1.set_ylabel('Y / cm')\n",
    "\n",
    "ax2.set_xlabel('X / cm')\n",
    "ax2.set_ylabel('Z / cm')\n",
    "im1=ax1.imshow(abs1.mean)\n",
    "fig.colorbar(im1,ax=ax1,shrink=0.4)\n",
    "im2=ax2.imshow(abs2.mean)\n",
    "fig.colorbar(im2,ax=ax2,shrink=0.4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}