diff --git a/.devcontainer/base.Dockerfile b/.devcontainer/base.Dockerfile index 2c30cc67..040b6207 100644 --- a/.devcontainer/base.Dockerfile +++ b/.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 && \ diff --git a/tasks/task_09_CSG_dose_tallies/1_surface_dose_from_gamma_source.ipynb b/tasks/task_09_CSG_dose_tallies/1_surface_dose_from_gamma_source.ipynb index a90006f0..2353d2bc 100644 --- a/tasks/task_09_CSG_dose_tallies/1_surface_dose_from_gamma_source.ipynb +++ b/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", diff --git a/tasks/task_09_CSG_dose_tallies/2_surface_dose_from_gamma_source_study.ipynb b/tasks/task_09_CSG_dose_tallies/2_surface_dose_from_gamma_source_study.ipynb index 5210e9be..2a89b467 100644 --- a/tasks/task_09_CSG_dose_tallies/2_surface_dose_from_gamma_source_study.ipynb +++ b/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", diff --git a/tasks/task_09_CSG_dose_tallies/3_cell_dose_from_neutron.py b/tasks/task_09_CSG_dose_tallies/3_cell_dose_from_neutron.py index 935e68ff..bb5553f3 100644 --- a/tasks/task_09_CSG_dose_tallies/3_cell_dose_from_neutron.py +++ b/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) diff --git a/tasks/task_09_CSG_dose_tallies/4_cell_dose_from_photon.py b/tasks/task_09_CSG_dose_tallies/4_cell_dose_from_photon.py index 83652e96..af178750 100644 --- a/tasks/task_09_CSG_dose_tallies/4_cell_dose_from_photon.py +++ b/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) diff --git a/tasks/task_09_CSG_dose_tallies/5_mesh_dose_from_neutrons.py b/tasks/task_09_CSG_dose_tallies/5_mesh_dose_from_neutrons.py new file mode 100644 index 00000000..aad0d5c1 --- /dev/null +++ b/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() diff --git a/tasks/task_09_CSG_dose_tallies/compare_dose_simulation_with_back_of_envelope.py b/tasks/task_09_CSG_dose_tallies/compare_dose_simulation_with_back_of_envelope.py index 86365d0c..705701d7 100644 --- a/tasks/task_09_CSG_dose_tallies/compare_dose_simulation_with_back_of_envelope.py +++ b/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) diff --git a/tasks/task_14_activation_transmutation_depletion/1_example_transmutation_isotope_build_up.ipynb b/tasks/task_14_activation_transmutation_depletion/1_example_transmutation_isotope_build_up.ipynb index 5d7e2d94..6baf39aa 100644 --- a/tasks/task_14_activation_transmutation_depletion/1_example_transmutation_isotope_build_up.ipynb +++ b/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" } } }, diff --git a/tasks/task_14_activation_transmutation_depletion/2_example_tally_change_with_burnup.ipynb b/tasks/task_14_activation_transmutation_depletion/2_example_tally_change_with_burnup.ipynb index eb766c91..327fd4ff 100644 --- a/tasks/task_14_activation_transmutation_depletion/2_example_tally_change_with_burnup.ipynb +++ b/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" } } },