Merge pull request #191 from fusion-energy/develop

fixing interpolation for dose tally

.devcontainer/base.Dockerfile

@@ -232,7 +232,13 @@ RUN cd /opt && \
 RUN pip install openmc_data && \
     mkdir -p /nuclear_data && \
     download_nndc_chain -d nuclear_data -r b8.0 && \
-    download_nndc -d nuclear_data -r b8.0
+    download_nndc -d nuclear_data -r b8.0 --cleanup && \
+    rm -rf nndc-b8.0-download
+# The last line here deletes the downloaded compressed nuclear datafile
+# not sure why but cleanup is not working, trackign on issue
+# https://github.com/openmc-data-storage/openmc_data/issues/29
+
+
 
 # install WMP nuclear data
 RUN wget https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/WMP_Library_v1.1.tar.gz && \

tasks/task_09_CSG_dose_tallies/1_surface_dose_from_gamma_source.ipynb

@@ -242,10 +242,10 @@
    "outputs": [],
    "source": [
     "energy_function_filter_n = openmc.EnergyFunctionFilter(energy_bins_n, dose_coeffs_n)\n",
-    "energy_function_filter_n.interpolation == 'cubic'\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",
+    "energy_function_filter_p.interpolation = 'cubic'  # cubic interpolation is recommended by ICRP\n",
     "\n",
     "photon_particle_filter = openmc.ParticleFilter([\"photon\"])\n",
     "surface_filter = openmc.SurfaceFilter(sphere_1)\n",

tasks/task_09_CSG_dose_tallies/2_surface_dose_from_gamma_source_study.ipynb

@@ -185,6 +185,7 @@
     "    energy_bins_p,\n",
     "    dose_coeffs_p\n",
     ")\n",
+    "energy_function_filter_p.interpolation = 'cubic'  # cubic interpolation is recommended by ICRP\n",
     "\n",
     "photon_particle_filter = openmc.ParticleFilter([\"photon\"])\n",
     "\n",

tasks/task_09_CSG_dose_tallies/3_cell_dose_from_neutron.py

@@ -66,14 +66,14 @@
     # openmc native units for length are cm so volume is in cm3
     phantom_volume = math.pi * math.pow(10.782, 2) * 169.75
 
-    # these are the dose coeffients coded into openmc
-    # originall from ICRP https://journals.sagepub.com/doi/10.1016/j.icrp.2011.10.001
+    # these are the dose coefficients coded into openmc
+    # originally from ICRP https://journals.sagepub.com/doi/10.1016/j.icrp.2011.10.001
 
     energy_bins_n, dose_coeffs_n = openmc.data.dose_coefficients(
         particle="neutron", geometry="AP"
     )
     energy_function_filter_n = openmc.EnergyFunctionFilter(energy_bins_n, dose_coeffs_n)
-    energy_function_filter_n.interpolation == "cubic"
+    energy_function_filter_n.interpolation = "cubic"  # cubic interpolation is recommended by ICRP
 
     neutron_particle_filter = openmc.ParticleFilter("neutron")
     cell_filter = openmc.CellFilter(phantom_cell)

tasks/task_09_CSG_dose_tallies/4_cell_dose_from_photon.py

@@ -86,7 +86,7 @@
         particle="photon", geometry="AP"
     )
     energy_function_filter_p = openmc.EnergyFunctionFilter(energy_bins_p, dose_coeffs_p)
-    energy_function_filter_p.interpolation == "cubic"
+    energy_function_filter_p.interpolation = "cubic"  # cubic interpolation is recommended by ICRP
 
     photon_particle_filter = openmc.ParticleFilter("photon")
     cell_filter = openmc.CellFilter(phantom_cell)

tasks/task_09_CSG_dose_tallies/5_mesh_dose_from_neutrons.py

@@ -0,0 +1,145 @@
+# This simulation a dose map and exports it as a vtk file and an image
+
+import openmc
+import matplotlib.pyplot as plt
+
+
+mat = openmc.Material()
+mat.add_element("Al", 1)
+mat.set_density("g/cm3", 2.7)
+my_materials = openmc.Materials([mat])
+
+cylinder_surface = openmc.ZCylinder(r=10)
+cylinder_upper_surface = openmc.ZPlane(z0=100)
+cylinder_lower_surface = openmc.ZPlane(z0=0)
+
+outer_surface = openmc.Sphere(r=200, boundary_type="vacuum")
+
+cylinder_region = -cylinder_surface & -cylinder_upper_surface & +cylinder_lower_surface
+
+# void region is below the outer surface and not the cylinder region
+void_region = -outer_surface & ~cylinder_region
+
+void_cell = openmc.Cell(region=void_region)
+cylinder_cell = openmc.Cell(region=cylinder_region)
+cylinder_cell.fill = mat
+
+my_geometry = openmc.Geometry([cylinder_cell, void_cell])
+
+# 14MeV point source
+source = openmc.Source()
+source.angle = openmc.stats.Isotropic()
+source.energy = openmc.stats.Discrete([14e6], [1])
+source.space = openmc.stats.Point((0.0, 50.0, 50.0))
+
+# Instantiate a Settings object
+my_settings = openmc.Settings()
+# when running a mesh tally simulation you might want to tell openmc not to save
+# the tallies.out file which is a ASCII file containing the tally results.
+# for mesh tallies this can get very large and take a long time to write.
+# the statepoint.h5 is smaller and quicker as it is a binary file
+my_settings.output = {"tallies": False}
+my_settings.batches = 2
+my_settings.particles = 500000
+my_settings.run_mode = "fixed source"
+my_settings.source = source
+
+# these are the dose coefficients coded into openmc
+# originally from ICRP https://journals.sagepub.com/doi/10.1016/j.icrp.2011.10.001
+energy_bins_n, dose_coeffs_n = openmc.data.dose_coefficients(
+    particle="neutron",
+    geometry="ISO",  # we are using the ISO direction as this is a dose field with dose
+)
+energy_function_filter_n = openmc.EnergyFunctionFilter(energy_bins_n, dose_coeffs_n)
+energy_function_filter_n.interpolation = "cubic"  # cubic interpolation is recommended by ICRP
+
+# just getting the dose for neutrons, not photons or other particles
+neutron_particle_filter = openmc.ParticleFilter("neutron")
+
+mesh = openmc.RegularMesh().from_domain(my_geometry, dimension=(30, 30, 30))
+mesh_filter = openmc.MeshFilter(mesh)
+
+# Create tally to score dose
+dose_cell_tally = openmc.Tally(name="neutron_dose_on_mesh")
+# note that the EnergyFunctionFilter is included as a filter
+dose_cell_tally.filters = [
+    mesh_filter,
+    neutron_particle_filter,
+    energy_function_filter_n,
+]
+dose_cell_tally.scores = ["flux"]
+my_tallies = openmc.Tallies([dose_cell_tally])
+
+model = openmc.Model(my_geometry, my_materials, my_settings, my_tallies)
+
+statepoint_filename = model.run()
+
+# makes use of a context manager "with" to automatically close the statepoint file
+with openmc.StatePoint(statepoint_filename) as statepoint:
+    my_mesh_tally_result = statepoint.get_tally(name="neutron_dose_on_mesh")
+
+# tally.mean is in units of pSv-cm3/source neutron
+# multiplication by neutrons_per_second changes units to neutron to pSv-cm3/second
+neutrons_per_second = 1e8  # units of neutrons per second
+
+# multiplication by pico_to_milli converts from (pico) pSv/second to (milli) mSv/second
+pico_to_milli = 1e-9
+
+# exports the mesh tally result to a vtk file with unit conversion
+mesh.write_data_to_vtk(
+    datasets={
+        "Dose [milli Sv per second]": my_mesh_tally_result.mean.flatten()
+        * neutrons_per_second
+        * pico_to_milli,
+    },
+    volume_normalization=True,  # this converts from dose-cm3/second to dose/second
+    filename="dose_on_mesh.vtk",
+)
+
+# this part of the script plots the images, so these imports are only needed to plot
+import openmc_geometry_plot  # extends openmc.Geometry class with plotting functions
+import regular_mesh_plotter  # extends openmc.Mesh class with plotting functions
+import matplotlib.pyplot as plt
+from matplotlib.colors import LogNorm
+import numpy as np
+
+# gets a 2d slice of data to later plot
+data_slice = mesh.slice_of_data(dataset=my_mesh_tally_result.mean, view_direction="x")
+
+data_slice = data_slice * neutrons_per_second * pico_to_milli
+
+
+plot_1 = plt.imshow(
+    data_slice,
+    extent=mesh.get_mpl_plot_extent(view_direction="x"),
+    interpolation=None,
+    norm=LogNorm(
+        vmin=1e-12,  # trims out the lower section of the colors
+        vmax=max(data_slice.flatten()),
+    ),
+)
+cbar = plt.colorbar(plot_1)
+cbar.set_label(f"Dose [milli Sv per second]")
+
+
+# gets unique levels for outlines contour plot and for the color scale
+material_ids = my_geometry.get_slice_of_material_ids(view_direction="x")
+# gets unique levels for outlines contour plot and for the color scale
+levels = np.unique([item for sublist in material_ids for item in sublist])
+
+plt.contour(
+    material_ids,
+    origin="upper",
+    colors="k",
+    linestyles="solid",
+    levels=levels,
+    linewidths=2.0,
+    extent=my_geometry.get_mpl_plot_extent(view_direction="x"),
+)
+xlabel, ylabel = my_geometry.get_axis_labels(view_direction="x")
+plt.xlabel(xlabel)
+plt.ylabel(ylabel)
+plt.title('Dose map showing some shielding of the source')
+
+
+plt.show()

tasks/task_09_CSG_dose_tallies/compare_dose_simulation_with_back_of_envelope.py

@@ -11,29 +11,44 @@
 import numpy as np
 
 
-def manual_dose_calc(particles_per_shot, distance_from_source, particle, particle_energy):
-    # as the model is so simple we can calculate the dose manually by finding
-    # neutrons across the surface area of the sphere.
-    # Conversion factor from fluence to dose at 14.1MeV = 495pSv cm2 per neutron (AP)
+def manual_dose_calc(
+    particles_per_shot:int,
+    distance_from_source:float,
+    particle:str,
+    particle_energy:float
+):
+    """Finds the dose in Sv a given distance from a point source 
+
+    Args:
+        particles_per_shot (int): the number of particles emitted by the source
+        distance_from_source (float): the distance between the person and source in cm
+        particle (str): the particle type "photon" or "neutron"
+        particle_energy (float): the particle energy in eV
+
+    Returns:
+        float: the dose in Sv
+    """
+    # As the model is so simple we can calculate the dose manually by finding
+    # neutrons across the surface area of the sphere and assuming no shielding
 
     sphere_surface_area = 4 * math.pi * math.pow(distance_from_source, 2)
 
     particles_per_unit_area = particles_per_shot / sphere_surface_area
 
-    # this obtains the does coefficients for the particle
+    # this obtains the does coefficients for the particle, AP is front facing (worst case)
     energy, dose_coeffs = openmc.data.dose_coefficients(particle, geometry='AP')
 
     # this gets the index of the particle energy
     closest_index = np.argmin(np.abs(np.array(energy)-particle_energy))
 
+    #dose coefficient from ICRP
+    # example conversion factor from fluence to dose at 14.1MeV = 495pSv cm2 per neutron (AP)
     dose_coeff = dose_coeffs[closest_index]
 
-    print('dose_coeff', dose_coeff)
 
     return particles_per_unit_area * dose_coeff * 1e-12
 
 
-
 mat_tissue = openmc.Material()
 mat_tissue.add_element("O", 76.2)
 mat_tissue.add_element("C", 11.1)

tasks/task_14_activation_transmutation_depletion/1_example_transmutation_isotope_build_up.ipynb

@@ -26,6 +26,12 @@
     "import openmc\n",
     "import openmc.deplete\n",
     "import os\n",
+    "\n",
+    "# this command makes use of a terminal command from the openmc_data package\n",
+    "# more info on openmc_data package can be found here https://github.com/openmc-data-storage/openmc_data\n",
+    "# this command downloads a chain file \n",
+    "# the chain file has information on half-lives, branching rations and gamma spec of isotopes\n",
+    "\n",
     "os.system('download_nndc_chain')\n"
    ]
   },
@@ -47,7 +53,7 @@
    "outputs": [],
    "source": [
     "\n",
-    "sphere_radius = 250\n",
+    "import math\n",
     "\n",
     "# MATERIALS\n",
     "\n",
@@ -57,7 +63,11 @@
     "my_material = openmc.Material() \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",
+    "\n",
+    "\n",
+    "sphere_radius = 100\n",
+    "volume_of_sphere = (4/3) * math.pi * math.pow(sphere_radius, 3)\n",
+    "my_material.volume = volume_of_sphere  # 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",
@@ -67,7 +77,7 @@
     "# GEOMETRY\n",
     "\n",
     "# surfaces\n",
-    "sph1 = openmc.Sphere(r=100, boundary_type='vacuum')\n",
+    "sph1 = openmc.Sphere(r=sphere_radius, boundary_type='vacuum')\n",
     "\n",
     "# cells, makes a simple sphere cell\n",
     "shield_cell = openmc.Cell(region=-sph1)\n",
@@ -306,11 +316,8 @@
   }
  ],
  "metadata": {
-  "interpreter": {
-   "hash": "5fc6d323269991ea8221605d53bd3c8fe6e348e3e76883e363c0d5ab9e3d0fe2"
-  },
   "kernelspec": {
-   "display_name": "Python 3.10.4 ('openmc_pauls_fork')",
+   "display_name": "Python 3.10.6 ('openmc_plot_dev')",
    "language": "python",
    "name": "python3"
   },
@@ -324,11 +331,11 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.4"
+   "version": "3.10.6"
   },
   "vscode": {
    "interpreter": {
-    "hash": "c47dabf1cc2568b64caa04441372f68f02228a3c450a6b1274acdfdc3a93f19e"
+    "hash": "c4091e6bedd918551d58365c181e219409ee0bad4691fd7a9c5816fcde549631"
    }
   }
  },

tasks/task_14_activation_transmutation_depletion/2_example_tally_change_with_burnup.ipynb

@@ -19,7 +19,7 @@
     "import openmc.deplete\n",
     "import math\n",
     "\n",
-    "lithium_orthosilicate_radius = 250\n",
+    "\n",
     "\n",
     "# MATERIALS\n",
     "\n",
@@ -29,8 +29,10 @@
     "breeding_material = openmc.Material() \n",
     "breeding_material.add_elements_from_formula('Li4SiO4')\n",
     "breeding_material.set_density('g/cm3', 2.5)\n",
-    "breeding_material.volume = (4/3) * math.pi * lithium_orthosilicate_radius**3\n",
-    "breeding_material.depletable = True\n",
+    "\n",
+    "lithium_orthosilicate_radius = 250\n",
+    "breeding_material.volume = (4/3) * math.pi * lithium_orthosilicate_radius**3 # a volume is needed so openmc can find the number of atoms in the cell/material\n",
+    "breeding_material.depletable = True  # depletable = True is needed to tell openmc to update the material with each time step\n",
     "\n",
     "materials = openmc.Materials([breeding_material])\n",
     "materials.export_to_xml()\n",
@@ -184,11 +186,8 @@
   }
  ],
  "metadata": {
-  "interpreter": {
-   "hash": "5fc6d323269991ea8221605d53bd3c8fe6e348e3e76883e363c0d5ab9e3d0fe2"
-  },
   "kernelspec": {
-   "display_name": "Python 3.10.4 ('openmc_pauls_fork')",
+   "display_name": "Python 3.10.6 ('openmc_plot_dev')",
    "language": "python",
    "name": "python3"
   },
@@ -202,11 +201,11 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.4"
+   "version": "3.10.6"
   },
   "vscode": {
    "interpreter": {
-    "hash": "c47dabf1cc2568b64caa04441372f68f02228a3c450a6b1274acdfdc3a93f19e"
+    "hash": "c4091e6bedd918551d58365c181e219409ee0bad4691fd7a9c5816fcde549631"
    }
   }
  },