Merge pull request #137 from fusion-energy/develop

General updates for new openmc API

Dockerfile

@@ -102,7 +102,8 @@ RUN pip install neutronics_material_maker[density] \
                 openmc-tally-unit-converter \
                 regular_mesh_plotter \
                 spectrum_plotter \
-                openmc_source_plotter
+                openmc_source_plotter \
+                openmc_depletion_plotter
 
 # Python libraries used in the workshop
 RUN pip install cmake\

tasks/task_04_make_sources/4_neutron_tracks.ipynb

@@ -146,7 +146,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "tracks.write_tracks_to_vtk('tracks.vtp')"
+    "tracks.write_to_vtk('tracks.vtp')"
    ]
   },
   {
@@ -270,7 +270,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.12"
+   "version": "3.9.13"
   }
  },
  "nbformat": 4,

tasks/task_09_CSG_surface_tally_dose/1_surface_dose_from_gamma_source.ipynb

@@ -226,9 +226,13 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Dose coeffients can then be used in a neutronics tally with the openmc.EnergyFunctionFilter.\n",
+    "Dose coefficients can then be used in a neutronics tally with the openmc.EnergyFunctionFilter.\n",
     "\n",
-    "This will effectivly multilpy the particle energy spectra with the effictive dose coefficient to produce a single number for effective dose."
+    "This will multiply the particle energy spectra with the effective dose coefficient to produce a single number for effective dose.\n",
+    "\n",
+    "ICRP recommend the use of 'cubic' interpolation.\n",
+    "\"For interpolations of absorbed dose and effective dose per fluence, a four-point (cubic) Lagrangian interpolation formula is recommended, and alog–log graph scale is more appropriate. Interpolations of absorbed dose and effective dose per air kerma of photons should be carried out using a four-point(cubic) Lagrangian interpolation formula, and a linear–log graph scale is more appropriate.\" \n",
+    "https://journals.sagepub.com/doi/pdf/10.1177/ANIB_40_2-5"
    ]
   },
   {
@@ -238,7 +242,10 @@
    "outputs": [],
    "source": [
     "energy_function_filter_n = openmc.EnergyFunctionFilter(energy_bins_n, dose_coeffs_n)\n",
+    "energy_function_filter_n.interpolation == 'cubic'\n",
+    "\n",
     "energy_function_filter_p = openmc.EnergyFunctionFilter(energy_bins_p, dose_coeffs_p)\n",
+    "energy_function_filter_p.interpolation == 'cubic'\n",
     "\n",
     "photon_particle_filter = openmc.ParticleFilter([\"photon\"])\n",
     "surface_filter = openmc.SurfaceFilter(sphere_1)\n",

tasks/task_14_activation_transmutation_depletion/1_example_transmutation_isotope_build_up.ipynb

@@ -13,7 +13,7 @@
    "id": "d4ccbdf3-2acc-45e8-8ece-62e5efb08dd0",
    "metadata": {},
    "source": [
-    "This first cell imports the packages needed, not the extra import openmc.deplete import"
+    "This first cell imports the packages needed, note the extra import openmc.deplete import"
    ]
   },
   {
@@ -24,9 +24,7 @@
    "outputs": [],
    "source": [
     "import openmc\n",
-    "import openmc.deplete\n",
-    "import matplotlib.pyplot as plt\n",
-    "import math"
+    "import openmc.deplete"
    ]
   },
   {
@@ -47,18 +45,18 @@
    "outputs": [],
    "source": [
     "\n",
-    "iron_sphere_radius = 250\n",
+    "sphere_radius = 250\n",
     "\n",
     "# MATERIALS\n",
     "\n",
     "mats = openmc.Materials()\n",
     "\n",
-    "# makes a simple material from Iron\n",
+    "# makes a simple material from Silver\n",
     "my_material = openmc.Material(name=\"my_material\") \n",
-    "my_material.add_nuclide('Co59', 1, percent_type='ao')\n",
-    "my_material.set_density('g/cm3', 7.7)\n",
-    "my_material.volume = (4/3) * math.pi * iron_sphere_radius**3\n",
-    "my_material.depletable = True\n",
+    "my_material.add_element('Ag', 1, percent_type='ao')\n",
+    "my_material.set_density('g/cm3', 10.49)\n",
+    "my_material.volume = 13502000  # a volume is needed so openmc can find the number of atoms in the cell/material\n",
+    "my_material.depletable = True  # depletable = True is needed to tell openmc to update the material with each time step\n",
     "\n",
     "materials = openmc.Materials([my_material])\n",
     "materials.export_to_xml()\n",
@@ -67,12 +65,11 @@
     "# GEOMETRY\n",
     "\n",
     "# surfaces\n",
-    "sph1 = openmc.Sphere(r=iron_sphere_radius, boundary_type='vacuum')\n",
+    "sph1 = openmc.Sphere(r=100, boundary_type='vacuum')\n",
     "\n",
     "# cells, makes a simple sphere cell\n",
     "shield_cell = openmc.Cell(region=-sph1)\n",
     "shield_cell.fill = my_material\n",
-    "shield_cell.volume = (4/3) * math.pi * sph1.r**3\n",
     "\n",
     "# sets the geometry to the universe that contains just the one cell\n",
     "universe = openmc.Universe(cells=[shield_cell])\n",
@@ -108,17 +105,11 @@
     "settings = openmc.Settings()\n",
     "settings.batches = 2\n",
     "settings.inactive = 0\n",
-    "settings.particles = 5000\n",
+    "settings.particles = 10000\n",
     "settings.source = source\n",
     "settings.run_mode = 'fixed source'\n",
     "\n",
-    "\n",
-    "tallies = openmc.Tallies()\n",
-    "geometry.export_to_xml()\n",
-    "settings.export_to_xml()\n",
-    "# tallies.export_to_xml()  # running in depletion mode doesn't need the tallies to be written out\n",
-    "materials.export_to_xml()\n",
-    "model = openmc.model.Model(geometry, materials, settings, tallies)"
+    "model = openmc.model.Model(geometry, materials, settings)"
    ]
   },
   {
@@ -127,11 +118,8 @@
    "metadata": {},
    "source": [
     "This is the depletion specific part of the model setup.\n",
-    "Here we:\n",
     "\n",
-    "    specify the chain file, this tells openmc the decay paths between isotopes including proabilities of different routes and half lives\n",
-    "\n",
-    "    set the time steps and corresponding source rates "
+    "This section specifies the chain file, this tells openmc the decay paths between isotopes including probabilities of different routes and half lives"
    ]
   },
   {
@@ -143,38 +131,110 @@
    "source": [
     "# This chain_endfb71 file was downloaded using the download_nndc_chain script that is included in the openmc_data package https://github.com/openmc-data-storage/openmc_data\\n\",\n",
     "# this file tells openmc the decay paths between isotopes including probabilities of different routes and half lives\n",
-    "chain_filename = '/tasks/task_14_activation_transmutation_depletion/chain-nndc-b7.1.xml'\n",
+    "openmc.config['chain_file'] = 'chain-nndc-b7.1.xml'\n",
+    "\n",
     "operator = openmc.deplete.Operator(\n",
     "    model=model,\n",
-    "    chain_file=chain_filename,\n",
     "    normalization_mode=\"source-rate\",  # set for fixed source simulation, otherwise defaults to fission simulation\n",
     "    dilute_initial=0,  # set to zero to avoid adding small amounts of isotopes, defaults to adding small amounts of fissionable isotopes\n",
-    "    reduce_chain=True # reduced to only the isotopes present in depletable materials and their possible progeny\n",
+    "    reduce_chain=True,  # reduced to only the isotopes present in depletable materials and their possible progeny\n",
+    "    reduce_chain_level=5,\n",
     ")\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "601d7b58",
+   "metadata": {},
+   "source": [
+    "This next stage sets the time steps and corresponding source rates for the irradiation schedule.\n",
+    "\n",
+    "An output file will be produced with showing the material composition at every time step.\n",
     "\n",
+    "We are irradiating the Silver for multiple half lives to show build up and saturation\n",
     "\n",
-    "operator = openmc.deplete.Operator(model, chain_filename)\n",
+    "Saturation happens when decay is = to creation of the particular isotope\n",
+    "\n",
+    "Ag110 half life is 24 seconds so it will start to become saturated after 120 seconds\n",
+    "\n",
+    "Ag108 half life is 145 seconds so it will not be saturated"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c35d5f9b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# We define timesteps together with the source rate to make it clearer\n",
+    "timesteps_and_source_rates = [\n",
+    "    (24, 1e20),\n",
+    "    (24, 1e20),\n",
+    "    (24, 1e20),\n",
+    "    (24, 1e20),\n",
+    "    (24, 1e20),  # should saturate Ag110 here as it has been irradiated for over 5 halflives\n",
+    "    (24, 1e20),\n",
+    "    (24, 1e20),\n",
+    "    (24, 1e20),\n",
+    "    (24, 1e20),\n",
+    "    (24, 0),\n",
+    "    (24, 0),\n",
+    "    (24, 0),\n",
+    "    (24, 0),\n",
+    "    (24, 0),\n",
+    "    (24, 0),\n",
+    "    (24, 0),\n",
+    "    (24, 0),\n",
+    "    (24, 0),\n",
+    "    (24, 0),\n",
+    "    (24, 0),\n",
+    "]\n",
+    "\n",
+    "# Uses list Python comprehension to get the timesteps and source_rates separately\n",
+    "timesteps = [item[0] for item in timesteps_and_source_rates]\n",
+    "source_rates = [item[1] for item in timesteps_and_source_rates]\n",
     "\n",
-    "# 1e9 neutrons per second for 5 years then 5 years of no neutrons (shut down cooling time)\n",
-    "time_steps = [365*24*60*60] * 5 + [365*24*60*60] * 5\n",
-    "source_rates = [1e9]*5 + [0] * 5\n",
     "\n",
     "# PredictorIntegrator has been selected as the depletion operator for this example as it is a fast first order Integrator\n",
     "# OpenMC offers several time-integration algorithms https://docs.openmc.org/en/stable/pythonapi/deplete.html#primary-api\\n\",\n",
     "# CF4Integrator should normally be selected as it appears to be the most accurate https://dspace.mit.edu/handle/1721.1/113721\\n\",\n",
     "integrator = openmc.deplete.PredictorIntegrator(\n",
-    "    operator=operator, timesteps=time_steps,source_rates=source_rates\n",
-    ")\n",
+    "    operator=operator,\n",
+    "    timesteps=timesteps,\n",
+    "    source_rates=source_rates\n",
+    ")\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ee0f0913",
+   "metadata": {},
+   "source": [
+    "This next section starts the depletion simulation and produces the output files"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f6a75285",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\n",
+    "integrator.integrate()\n",
     "\n",
-    "integrator.integrate()\n"
+    "# bash command to show the output files produce\n",
+    "!ls -lh *.h5"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "2cb387d1-d9bb-45bb-89ff-dea703fa3c3e",
    "metadata": {},
    "source": [
-    "This section extracts the results of the depletion simulation from the h5 file and gets the amount of Cobolt 60 (Co60) in the material at each of the time steps"
+    "This section extracts the results of the depletion simulation from the h5 file and gets the amount of Ag110 in the material at each of the time steps"
    ]
   },
   {
@@ -186,15 +246,18 @@
    "source": [
     "results = openmc.deplete.ResultsList.from_hdf5(\"depletion_results.h5\")\n",
     "\n",
-    "times, number_of_co60_atoms = results.get_atoms(my_material, 'Co60')"
+    "times, number_of_Ag110_atoms = results.get_atoms(my_material, 'Ag110')\n",
+    "\n",
+    "for time, num in zip(times, number_of_Ag110_atoms):\n",
+    "    print(f\" Time {time}s. Number of Ag110 atoms {num}\")"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "072b7db4-e37c-47bd-9c59-40f2a5a85e83",
    "metadata": {},
    "source": [
-    "This section plots the results"
+    "In addition to Ag110 other atoms get created. This section plots the number of nuclides in the material excluding the original nuclides in the unirradiated material"
    ]
   },
   {
@@ -204,22 +267,48 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import matplotlib.pyplot as plt\n",
+    "import openmc_depletion_plotter\n",
+    "# this package provides convenient plotting methods for depletion simulations like this one\n",
+    "# more details here https://github.com/fusion-energy/openmc_depletion_plotter\n",
     "\n",
-    "fig, ax = plt.subplots()\n",
-    "ax.plot(times, number_of_co60_atoms)\n",
+    "results.plot_atoms_vs_time(excluded_material=my_material)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "819a86c6",
+   "metadata": {},
+   "source": [
+    "Not all nuclide are unstable and the unstable ones have a different half life. This next plot shows the specific activity (activity per unit mass) as a function of time.\n",
     "\n",
-    "ax.set(xlabel='time (s)', ylabel='Number of atoms',\n",
-    "       title='Build up of atoms saturates when decay is equal to activation this occurs at circa 5 half lives')\n",
-    "ax.grid()\n",
-    "plt.savefig('atoms.png')\n",
-    "plt.show()"
+    "This is useful for identifying a suitable waste repository for activated waste."
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1aecca70",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "results.plot_activity_vs_time()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "966f5b78",
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
+  "interpreter": {
+   "hash": "5fc6d323269991ea8221605d53bd3c8fe6e348e3e76883e363c0d5ab9e3d0fe2"
+  },
   "kernelspec": {
-   "display_name": "Python 3.8.12 ('paramak_dev')",
+   "display_name": "Python 3.10.4 ('openmc_pauls_fork')",
    "language": "python",
    "name": "python3"
   },
@@ -233,7 +322,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.12"
+   "version": "3.10.4"
   },
   "vscode": {
    "interpreter": {

tasks/task_14_activation_transmutation_depletion/2_example_tally_change_with_burnup.ipynb

@@ -17,7 +17,6 @@
    "source": [
     "import openmc\n",
     "import openmc.deplete\n",
-    "import matplotlib.pyplot as plt\n",
     "import math\n",
     "\n",
     "lithium_orthosilicate_radius = 250\n",
@@ -122,11 +121,10 @@
    "source": [
     "# This chain_endfb71 file was downloaded using the download_nndc_chain script that is included in the openmc_data package https://github.com/openmc-data-storage/openmc_data\\n\",\n",
     "# this file tells openmc the decay paths between isotopes including probabilities of different routes and half lives\n",
-    "chain_filename = '/tasks/task_14_activation_transmutation_depletion/chain-nndc-b7.1.xml'\n",
+    "openmc.config['chain_file'] = 'chain-nndc-b7.1.xml'\n",
     "\n",
     "operator = openmc.deplete.Operator(\n",
     "    model=model,\n",
-    "    chain_file=chain_filename,\n",
     "    normalization_mode=\"source-rate\",  # set for fixed source simulation, otherwise defaults to fission simulation\n",
     "    dilute_initial=0,  # set to zero to avoid adding small amounts of isotopes, defaults to adding small amounts of fissionable isotopes\n",
     "    reduce_chain=True # reduced to only the isotopes present in depletable materials and their possible progeny\n",
@@ -186,8 +184,11 @@
   }
  ],
  "metadata": {
+  "interpreter": {
+   "hash": "5fc6d323269991ea8221605d53bd3c8fe6e348e3e76883e363c0d5ab9e3d0fe2"
+  },
   "kernelspec": {
-   "display_name": "Python 3.8.12 ('paramak_dev')",
+   "display_name": "Python 3.10.4 ('openmc_pauls_fork')",
    "language": "python",
    "name": "python3"
   },
@@ -201,7 +202,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.12"
+   "version": "3.10.4"
   },
   "vscode": {
    "interpreter": {