From 68fccaaee036a932a96a26a2cae21cf5b12c6b05 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Sat, 2 Dec 2023 22:35:12 +0000 Subject: [PATCH 01/13] std::bad_alloc --- ..._mesh_based_shut_down_dose_rate_example.py | 287 ++++++++++++++++++ 1 file changed, 287 insertions(+) create mode 100644 tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py new file mode 100644 index 00000000..d5d82b73 --- /dev/null +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -0,0 +1,287 @@ +# This script simulates R2S method of shut down dose rate +# on a simple sphere model. + +import numpy as np +import openmc +import openmc.deplete +from pathlib import Path +import math +from matplotlib.colors import LogNorm + +# users might want to change these to use specific xml files to use particular decay data or transport cross sections +# the chain file was downloaded with +# pip install openmc_data +# download_endf_chain -r b8.0 +# openmc.config['chain_file'] = '/nuclear_data/chain-endf-b8.0.xml' +# openmc.config['cross_sections'] = 'cross_sections.xml' + +# a few user settings +# Set up the folders to save all the data in +n_particles = 1_00000 +p_particles = 1_000 +statepoints_folder = Path('statepoints_folder') + +# First we make a simple geometry with three cells, (two with material) +sphere_surf_1 = openmc.Sphere(r=20, boundary_type="vacuum") +sphere_surf_2 = openmc.Sphere(r=1, y0=10) +sphere_surf_3 = openmc.Sphere(r=5, z0=10) + +sphere_region_1 = -sphere_surf_1 & +sphere_surf_2 & +sphere_surf_3 # void space +sphere_region_2 = -sphere_surf_2 +sphere_region_3 = -sphere_surf_3 + +sphere_cell_1 = openmc.Cell(region=sphere_region_1) +sphere_cell_2 = openmc.Cell(region=sphere_region_2) +sphere_cell_3 = openmc.Cell(region=sphere_region_3) + +# We make a iron material which should produce a few activation products +mat_iron = openmc.Material() +mat_iron.id = 1 +mat_iron.add_element("Fe", 1.0) +mat_iron.set_density("g/cm3", 7.7) +# must set the depletion to True to deplete the material +mat_iron.depletable = True +# volume must set the volume as well as openmc calculates number of atoms +mat_iron.volume = (4 / 3) * math.pi * math.pow(sphere_surf_2.r, 3) +sphere_cell_2.fill = mat_iron + +# We make a Al material which should produce a few different activation products +mat_aluminum = openmc.Material() +mat_aluminum.id = 2 +mat_aluminum.add_element("Al", 1.0) +mat_aluminum.set_density("g/cm3", 2.7) +# must set the depletion to True to deplete the material +mat_aluminum.depletable = True +# volume must set the volume as well as openmc calculates number of atoms +mat_aluminum.volume = (4 / 3) * math.pi * math.pow(sphere_surf_3.r, 3) +sphere_cell_3.fill = mat_aluminum + +my_geometry = openmc.Geometry([sphere_cell_1, sphere_cell_2, sphere_cell_3]) + +my_materials = openmc.Materials([mat_iron, mat_aluminum]) + +pristine_mat_iron = mat_iron.clone() +pristine_mat_aluminium = mat_aluminum.clone() + +# gets the cell ids of any depleted cell +activated_cell_ids = [c.id for c in my_geometry.get_all_material_cells().values() if c.fill.depletable] +print("activated_cell_ids", activated_cell_ids) +all_depletable_cells = [c for _, c in my_geometry.get_all_material_cells().items() if c.fill.depletable is True] +print("depletable_cell_ids", activated_cell_ids) +all_depletable_materials = [c.fill for _, c in my_geometry.get_all_material_cells().items() if c.fill.depletable is True] + +# 14MeV neutron source that activates material +my_source = openmc.IndependentSource() +my_source.space = openmc.stats.Point((0, 0, 0)) +my_source.angle = openmc.stats.Isotropic() +my_source.energy = openmc.stats.Discrete([14.06e6], [1]) +my_source.particle = "neutron" + +# settings for the neutron simulation(s) +my_neutron_settings = openmc.Settings() +my_neutron_settings.run_mode = "fixed source" +my_neutron_settings.particles = n_particles +my_neutron_settings.batches = 10 +my_neutron_settings.source = my_source +my_neutron_settings.photon_transport = False + +# Create mesh which will be used for tally +regular_mesh = openmc.RegularMesh().from_domain( + my_geometry, # the corners of the mesh are being set automatically to surround the geometry + dimension=[10, 10, 10] # 10 +) + +model_neutron = openmc.Model(my_geometry, my_materials, my_neutron_settings) + +hour_in_seconds = 60*60 + +# This section defines the neutron pulse schedule. +# If the method made use of the CoupledOperator then there would need to be a +# transport simulation for each timestep. However as the IndependentOperator is +# used then just a single transport simulation is done, thus speeding up the +# simulation considerably. +timesteps_and_source_rates = [ + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 1 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 2 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 3 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 4 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 4 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 4 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 5 hour +] + +timesteps = [item[0] for item in timesteps_and_source_rates] +source_rates = [item[1] for item in timesteps_and_source_rates] + +model_neutron.export_to_xml(directory=statepoints_folder/ "neutrons") + +# this does perform transport but just to get the flux and micro xs +flux_in_each_group, micro_xs = openmc.deplete.get_microxs_and_flux( + model=model_neutron, + domains=regular_mesh, + energies='CCFE-709', # different group structures see this file for all the groups available https://github.com/openmc-dev/openmc/blob/develop/openmc/mgxs/__init__.py +) + +# # constructing the operator, note we pass in the flux and micro xs +# operator = openmc.deplete.IndependentOperator( +# materials=openmc.Materials(all_depletable_materials), +# fluxes=flux_in_each_group, +# micros=micro_xs, +# reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny +# reduce_chain_level=5, +# normalization_mode="source-rate" +# ) + +# integrator = openmc.deplete.PredictorIntegrator( +# operator=operator, +# timesteps=timesteps, +# source_rates=source_rates, +# timestep_units='s' +# ) + +# this runs the depletion calculations for the timesteps +# this does the neutron activation simulations and produces a depletion_results.h5 file +integrator.integrate() +# TODO add output dir to integrate command so we don't have to move the file like this +# integrator.integrate(path=statepoints_folder / "neutrons" / "depletion_results.h5") +# PR on openmc is open +import os +os.system(f'mv depletion_results.h5 {statepoints_folder / "neutrons" / "depletion_results.h5"}') + +# Now we have done the neutron activation simulations we can start the work needed for the decay gamma simulations. + +my_gamma_settings = openmc.Settings() +my_gamma_settings.run_mode = "fixed source" +my_gamma_settings.batches = 100 +my_gamma_settings.particles = p_particles + + +# First we add make dose tally on a regular mesh + + +# creates a regular mesh that surrounds the geometry +mesh = openmc.RegularMesh().from_domain( + my_geometry, + dimension=[10, 10, 10], # 10 voxels in each axis direction (x, y, z) +) + +# adding a dose tally on a regular mesh +# AP, PA, LLAT, RLAT, ROT, ISO are ICRP incident dose field directions, AP is front facing +energies, pSv_cm2 = openmc.data.dose_coefficients(particle="photon", geometry="AP") +dose_filter = openmc.EnergyFunctionFilter( + energies, pSv_cm2, interpolation="cubic" # interpolation method recommended by ICRP +) +particle_filter = openmc.ParticleFilter(["photon"]) +mesh_filter = openmc.MeshFilter(mesh) +flux_tally = openmc.Tally() +flux_tally.filters = [mesh_filter, dose_filter, particle_filter] +flux_tally.scores = ["flux"] +flux_tally.name = "photon_dose_on_mesh" + +tallies = openmc.Tallies([flux_tally]) + +cells = model_neutron.geometry.get_all_cells() +activated_cells = [cells[uid] for uid in activated_cell_ids] + +# this section makes the photon sources from each active material at each +# timestep and runs the photon simulations +results = openmc.deplete.Results(statepoints_folder / "neutrons" / "depletion_results.h5") + +for i_cool in range(1, len(timesteps)): + + # range starts at 1 to skip the first step as that is an irradiation step and there is no + # decay gamma source from the stable material at that time + # also there are no decay products in this first timestep for this model + + photon_sources_for_timestep = [] + print(f"making photon source for timestep {i_cool}") + + all_activated_materials_in_timestep = [] + + for activated_cell_id in activated_cell_ids: + # gets the material id of the material filling the cell + material_id = cells[activated_cell_id].fill.id + + # gets the activated material using the material id + activated_mat = results[i_cool].get_material(str(material_id)) + # gets the energy and probabilities for the + energy = activated_mat.get_decay_photon_energy( + clip_tolerance = 1e-6, # cuts out a small fraction of the very low energy (and hence negligible dose contribution) photons + units = 'Bq', + ) + strength = energy.integral() + + if strength > 0.: # only makes sources for + space = openmc.stats.Box(*cells[activated_cell_id].bounding_box) + source = openmc.IndependentSource( + space=space, + energy=energy, + particle="photon", + strength=strength, + domains=[cells[activated_cell_id]], + ) + photon_sources_for_timestep.append(source) + + + my_gamma_settings.source = photon_sources_for_timestep + + + # one should also fill the cells with the activated material + # the activated material contains ALL the nuclides produced during activation + # sphere_cell_2.fill = results[i_cool].get_material("1") + # sphere_cell_3.fill = results[i_cool].get_material("2") + # my_geometry = openmc.Geometry([sphere_cell_1, sphere_cell_2, sphere_cell_3]) + + # however it is unlikely that they all appear in your transport cross_sections.xml + # so you could make use of openmc.deplete.Results.export_to_materials to export the modified activated material that + # just contains isotopes that appear in your cross_sections.xml + + # however in this example we just use the original pristine material my_materials that were cloned earlier + # my_geometry is also the same as the neutron simulation + pristine_mat_iron.id = 1 + pristine_mat_aluminium.id =2 + my_materials = openmc.Materials([pristine_mat_iron, pristine_mat_aluminium]) + + model_gamma = openmc.Model(my_geometry, my_materials, my_gamma_settings, tallies) + + model_gamma.run(cwd=statepoints_folder / "photons" / f"photon_at_time_{i_cool}") + + +# pico_to_micro = 1e-6 +# seconds_to_hours = 60*60 + +# # You may wish to plot the dose tally on a mesh, this package makes it easy to include the geometry with the mesh tally +# from openmc_regular_mesh_plotter import plot_mesh_tally +# for i_cool in range(1, len(timesteps)): +# with openmc.StatePoint(statepoints_folder / "photons" / f"photon_at_time_{i_cool}" / 'statepoint.100.h5') as statepoint: +# photon_tally = statepoint.get_tally(name="photon_dose_on_mesh") + +# # normalising this tally is a little different to other examples as the source strength has been using units of photons per second. +# # tally.mean is in units of pSv-cm3/source photon. +# # as source strength is in photons_per_second this changes units to pSv-/second + +# # multiplication by pico_to_micro converts from (pico) pSv/s to (micro) uSv/s +# # dividing by mesh voxel volume cancles out the cm3 units +# # could do the mesh volume scaling on the plot and vtk functions but doing it here instead +# scaling_factor = (seconds_to_hours * pico_to_micro) / mesh.volumes[0][0][0] + +# plot = plot_mesh_tally( +# tally=photon_tally, +# basis="xz", +# # score='flux', # only one tally so can make use of default here +# value="mean", +# colorbar_kwargs={ +# 'label': "Decay photon dose [µSv/h]", +# }, +# norm=LogNorm(), +# volume_normalization=False, # this is done in the scaling_factor +# scaling_factor=scaling_factor, +# ) +# plot.figure.savefig(f'shut_down_dose_map_timestep_{i_cool}') From de5f74694bbcad2cd747633489b8288ff5a14dd4 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Sat, 2 Dec 2023 23:37:34 +0000 Subject: [PATCH 02/13] more materials needed --- ..._mesh_based_shut_down_dose_rate_example.py | 107 ++++++++++-------- 1 file changed, 58 insertions(+), 49 deletions(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index d5d82b73..c0a530d9 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -49,6 +49,7 @@ mat_aluminum = openmc.Material() mat_aluminum.id = 2 mat_aluminum.add_element("Al", 1.0) +mat_aluminum.add_element("Fe", 1.0) mat_aluminum.set_density("g/cm3", 2.7) # must set the depletion to True to deplete the material mat_aluminum.depletable = True @@ -122,29 +123,37 @@ model_neutron.export_to_xml(directory=statepoints_folder/ "neutrons") +all_nuclides = [] +for material in my_geometry.get_all_materials().values(): + for nuclide in material.get_nuclides(): + if nuclide not in all_nuclides: + all_nuclides.append(nuclide) +# print(set(all_nuclides)) + # this does perform transport but just to get the flux and micro xs flux_in_each_group, micro_xs = openmc.deplete.get_microxs_and_flux( model=model_neutron, domains=regular_mesh, energies='CCFE-709', # different group structures see this file for all the groups available https://github.com/openmc-dev/openmc/blob/develop/openmc/mgxs/__init__.py + nuclides=all_nuclides ) # # constructing the operator, note we pass in the flux and micro xs -# operator = openmc.deplete.IndependentOperator( -# materials=openmc.Materials(all_depletable_materials), -# fluxes=flux_in_each_group, -# micros=micro_xs, -# reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny -# reduce_chain_level=5, -# normalization_mode="source-rate" -# ) - -# integrator = openmc.deplete.PredictorIntegrator( -# operator=operator, -# timesteps=timesteps, -# source_rates=source_rates, -# timestep_units='s' -# ) +operator = openmc.deplete.IndependentOperator( + materials=#TODO find the materials in each mesh voxel and their volume fractions, + fluxes=flux_in_each_group, + micros=micro_xs, + reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny + reduce_chain_level=5, + normalization_mode="source-rate" +) + +integrator = openmc.deplete.PredictorIntegrator( + operator=operator, + timesteps=timesteps, + source_rates=source_rates, + timestep_units='s' +) # this runs the depletion calculations for the timesteps # this does the neutron activation simulations and produces a depletion_results.h5 file @@ -194,43 +203,43 @@ # timestep and runs the photon simulations results = openmc.deplete.Results(statepoints_folder / "neutrons" / "depletion_results.h5") -for i_cool in range(1, len(timesteps)): +# for i_cool in range(1, len(timesteps)): - # range starts at 1 to skip the first step as that is an irradiation step and there is no - # decay gamma source from the stable material at that time - # also there are no decay products in this first timestep for this model +# # range starts at 1 to skip the first step as that is an irradiation step and there is no +# # decay gamma source from the stable material at that time +# # also there are no decay products in this first timestep for this model - photon_sources_for_timestep = [] - print(f"making photon source for timestep {i_cool}") +# photon_sources_for_timestep = [] +# print(f"making photon source for timestep {i_cool}") - all_activated_materials_in_timestep = [] +# all_activated_materials_in_timestep = [] - for activated_cell_id in activated_cell_ids: - # gets the material id of the material filling the cell - material_id = cells[activated_cell_id].fill.id +# for activated_cell_id in activated_cell_ids: +# # gets the material id of the material filling the cell +# material_id = cells[activated_cell_id].fill.id - # gets the activated material using the material id - activated_mat = results[i_cool].get_material(str(material_id)) - # gets the energy and probabilities for the - energy = activated_mat.get_decay_photon_energy( - clip_tolerance = 1e-6, # cuts out a small fraction of the very low energy (and hence negligible dose contribution) photons - units = 'Bq', - ) - strength = energy.integral() +# # gets the activated material using the material id +# activated_mat = results[i_cool].get_material(str(material_id)) +# # gets the energy and probabilities for the +# energy = activated_mat.get_decay_photon_energy( +# clip_tolerance = 1e-6, # cuts out a small fraction of the very low energy (and hence negligible dose contribution) photons +# units = 'Bq', +# ) +# strength = energy.integral() - if strength > 0.: # only makes sources for - space = openmc.stats.Box(*cells[activated_cell_id].bounding_box) - source = openmc.IndependentSource( - space=space, - energy=energy, - particle="photon", - strength=strength, - domains=[cells[activated_cell_id]], - ) - photon_sources_for_timestep.append(source) +# if strength > 0.: # only makes sources for +# space = openmc.stats.Box(*cells[activated_cell_id].bounding_box) +# source = openmc.IndependentSource( +# space=space, +# energy=energy, +# particle="photon", +# strength=strength, +# domains=[cells[activated_cell_id]], #this should be the mesh voxel and the material in the mesh voxel +# ) +# photon_sources_for_timestep.append(source) - my_gamma_settings.source = photon_sources_for_timestep +# my_gamma_settings.source = photon_sources_for_timestep # one should also fill the cells with the activated material @@ -245,13 +254,13 @@ # however in this example we just use the original pristine material my_materials that were cloned earlier # my_geometry is also the same as the neutron simulation - pristine_mat_iron.id = 1 - pristine_mat_aluminium.id =2 - my_materials = openmc.Materials([pristine_mat_iron, pristine_mat_aluminium]) + # pristine_mat_iron.id = 1 + # pristine_mat_aluminium.id =2 + # my_materials = openmc.Materials([pristine_mat_iron, pristine_mat_aluminium]) - model_gamma = openmc.Model(my_geometry, my_materials, my_gamma_settings, tallies) + # model_gamma = openmc.Model(my_geometry, my_materials, my_gamma_settings, tallies) - model_gamma.run(cwd=statepoints_folder / "photons" / f"photon_at_time_{i_cool}") + # model_gamma.run(cwd=statepoints_folder / "photons" / f"photon_at_time_{i_cool}") # pico_to_micro = 1e-6 From 4a8e80ddb06779efa573037415ca5aa2efb991e5 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Sat, 2 Dec 2023 23:41:16 +0000 Subject: [PATCH 03/13] perhaps from nuclide --- .../3_mesh_based_shut_down_dose_rate_example.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index c0a530d9..09c1ce24 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -139,8 +139,8 @@ ) # # constructing the operator, note we pass in the flux and micro xs -operator = openmc.deplete.IndependentOperator( - materials=#TODO find the materials in each mesh voxel and their volume fractions, +operator = openmc.deplete.IndependentOperator().from_nuclides( + materials=#TODO find the nuclides in each mesh voxel and their volume fractions, fluxes=flux_in_each_group, micros=micro_xs, reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny From f487381f9e74c0c83671e9a431665db0372f57ea Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 17 Jan 2024 18:06:10 +0000 Subject: [PATCH 04/13] added mesh vol step --- ..._mesh_based_shut_down_dose_rate_example.py | 145 ++++++++++-------- 1 file changed, 82 insertions(+), 63 deletions(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index 09c1ce24..af0d3cd7 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -14,6 +14,7 @@ # download_endf_chain -r b8.0 # openmc.config['chain_file'] = '/nuclear_data/chain-endf-b8.0.xml' # openmc.config['cross_sections'] = 'cross_sections.xml' +openmc.config['chain_file'] = '/home/jshimwell/ENDF-B-VIII.0-NNDC/chain-nndc-b8.0.xml' # a few user settings # Set up the folders to save all the data in @@ -135,73 +136,91 @@ model=model_neutron, domains=regular_mesh, energies='CCFE-709', # different group structures see this file for all the groups available https://github.com/openmc-dev/openmc/blob/develop/openmc/mgxs/__init__.py - nuclides=all_nuclides + nuclides=all_nuclides, + chain_file=openmc.config['chain_file'] ) -# # constructing the operator, note we pass in the flux and micro xs -operator = openmc.deplete.IndependentOperator().from_nuclides( - materials=#TODO find the nuclides in each mesh voxel and their volume fractions, - fluxes=flux_in_each_group, - micros=micro_xs, - reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny - reduce_chain_level=5, - normalization_mode="source-rate" +model_neutron.export_to_model_xml() +import openmc.lib +openmc.lib.init() +mesh = openmc.lib.RegularMesh() +mesh.dimension = regular_mesh.dimension +mesh.set_parameters( + lower_left=regular_mesh.lower_left, + upper_right=regular_mesh.upper_right ) -integrator = openmc.deplete.PredictorIntegrator( - operator=operator, - timesteps=timesteps, - source_rates=source_rates, - timestep_units='s' -) - -# this runs the depletion calculations for the timesteps -# this does the neutron activation simulations and produces a depletion_results.h5 file -integrator.integrate() -# TODO add output dir to integrate command so we don't have to move the file like this -# integrator.integrate(path=statepoints_folder / "neutrons" / "depletion_results.h5") -# PR on openmc is open -import os -os.system(f'mv depletion_results.h5 {statepoints_folder / "neutrons" / "depletion_results.h5"}') - -# Now we have done the neutron activation simulations we can start the work needed for the decay gamma simulations. - -my_gamma_settings = openmc.Settings() -my_gamma_settings.run_mode = "fixed source" -my_gamma_settings.batches = 100 -my_gamma_settings.particles = p_particles - - -# First we add make dose tally on a regular mesh - - -# creates a regular mesh that surrounds the geometry -mesh = openmc.RegularMesh().from_domain( - my_geometry, - dimension=[10, 10, 10], # 10 voxels in each axis direction (x, y, z) -) - -# adding a dose tally on a regular mesh -# AP, PA, LLAT, RLAT, ROT, ISO are ICRP incident dose field directions, AP is front facing -energies, pSv_cm2 = openmc.data.dose_coefficients(particle="photon", geometry="AP") -dose_filter = openmc.EnergyFunctionFilter( - energies, pSv_cm2, interpolation="cubic" # interpolation method recommended by ICRP -) -particle_filter = openmc.ParticleFilter(["photon"]) -mesh_filter = openmc.MeshFilter(mesh) -flux_tally = openmc.Tally() -flux_tally.filters = [mesh_filter, dose_filter, particle_filter] -flux_tally.scores = ["flux"] -flux_tally.name = "photon_dose_on_mesh" - -tallies = openmc.Tallies([flux_tally]) - -cells = model_neutron.geometry.get_all_cells() -activated_cells = [cells[uid] for uid in activated_cell_ids] - -# this section makes the photon sources from each active material at each -# timestep and runs the photon simulations -results = openmc.deplete.Results(statepoints_folder / "neutrons" / "depletion_results.h5") +vols = mesh.material_volumes(n_samples = 1000000) + +openmc.lib.finalize() + +# # # # constructing the operator, note we pass in the flux and micro xs +# operator = openmc.deplete.IndependentOperator().from_nuclides( +# volume= +# nuclides= + +# materials=#TODO find the nuclides in each mesh voxel and their volume fractions, +# fluxes=flux_in_each_group, +# micros=micro_xs, +# reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny +# reduce_chain_level=5, +# normalization_mode="source-rate" +# ) + +# integrator = openmc.deplete.PredictorIntegrator( +# operator=operator, +# timesteps=timesteps, +# source_rates=source_rates, +# timestep_units='s' +# ) + +# # this runs the depletion calculations for the timesteps +# # this does the neutron activation simulations and produces a depletion_results.h5 file +# integrator.integrate() +# # TODO add output dir to integrate command so we don't have to move the file like this +# # integrator.integrate(path=statepoints_folder / "neutrons" / "depletion_results.h5") +# # PR on openmc is open +# import os +# os.system(f'mv depletion_results.h5 {statepoints_folder / "neutrons" / "depletion_results.h5"}') + +# # Now we have done the neutron activation simulations we can start the work needed for the decay gamma simulations. + +# my_gamma_settings = openmc.Settings() +# my_gamma_settings.run_mode = "fixed source" +# my_gamma_settings.batches = 100 +# my_gamma_settings.particles = p_particles + + +# # First we add make dose tally on a regular mesh + + +# # creates a regular mesh that surrounds the geometry +# mesh = openmc.RegularMesh().from_domain( +# my_geometry, +# dimension=[10, 10, 10], # 10 voxels in each axis direction (x, y, z) +# ) + +# # adding a dose tally on a regular mesh +# # AP, PA, LLAT, RLAT, ROT, ISO are ICRP incident dose field directions, AP is front facing +# energies, pSv_cm2 = openmc.data.dose_coefficients(particle="photon", geometry="AP") +# dose_filter = openmc.EnergyFunctionFilter( +# energies, pSv_cm2, interpolation="cubic" # interpolation method recommended by ICRP +# ) +# particle_filter = openmc.ParticleFilter(["photon"]) +# mesh_filter = openmc.MeshFilter(mesh) +# flux_tally = openmc.Tally() +# flux_tally.filters = [mesh_filter, dose_filter, particle_filter] +# flux_tally.scores = ["flux"] +# flux_tally.name = "photon_dose_on_mesh" + +# tallies = openmc.Tallies([flux_tally]) + +# cells = model_neutron.geometry.get_all_cells() +# activated_cells = [cells[uid] for uid in activated_cell_ids] + +# # this section makes the photon sources from each active material at each +# # timestep and runs the photon simulations +# results = openmc.deplete.Results(statepoints_folder / "neutrons" / "depletion_results.h5") # for i_cool in range(1, len(timesteps)): From aefa7f34e6fbb102e632c394a00937b017a8c24d Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 1 Feb 2024 15:47:40 +0000 Subject: [PATCH 05/13] making materials --- ..._mesh_based_shut_down_dose_rate_example.py | 61 +++++++++++++------ 1 file changed, 41 insertions(+), 20 deletions(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index af0d3cd7..b0f0a35d 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -22,44 +22,52 @@ p_particles = 1_000 statepoints_folder = Path('statepoints_folder') -# First we make a simple geometry with three cells, (two with material) -sphere_surf_1 = openmc.Sphere(r=20, boundary_type="vacuum") -sphere_surf_2 = openmc.Sphere(r=1, y0=10) -sphere_surf_3 = openmc.Sphere(r=5, z0=10) - -sphere_region_1 = -sphere_surf_1 & +sphere_surf_2 & +sphere_surf_3 # void space -sphere_region_2 = -sphere_surf_2 -sphere_region_3 = -sphere_surf_3 -sphere_cell_1 = openmc.Cell(region=sphere_region_1) -sphere_cell_2 = openmc.Cell(region=sphere_region_2) -sphere_cell_3 = openmc.Cell(region=sphere_region_3) +al_sphere_radius = 7 +iron_sphere_radius = 4 # We make a iron material which should produce a few activation products mat_iron = openmc.Material() mat_iron.id = 1 -mat_iron.add_element("Fe", 1.0) +mat_iron.add_nuclide("Fe56", 1.0) +mat_iron.add_nuclide("Fe57", 1.0) mat_iron.set_density("g/cm3", 7.7) # must set the depletion to True to deplete the material mat_iron.depletable = True # volume must set the volume as well as openmc calculates number of atoms -mat_iron.volume = (4 / 3) * math.pi * math.pow(sphere_surf_2.r, 3) -sphere_cell_2.fill = mat_iron +mat_iron.volume = (4 / 3) * math.pi * math.pow(iron_sphere_radius, 3) # We make a Al material which should produce a few different activation products mat_aluminum = openmc.Material() mat_aluminum.id = 2 mat_aluminum.add_element("Al", 1.0) -mat_aluminum.add_element("Fe", 1.0) mat_aluminum.set_density("g/cm3", 2.7) # must set the depletion to True to deplete the material mat_aluminum.depletable = True # volume must set the volume as well as openmc calculates number of atoms -mat_aluminum.volume = (4 / 3) * math.pi * math.pow(sphere_surf_3.r, 3) -sphere_cell_3.fill = mat_aluminum +mat_aluminum.volume = (4 / 3) * math.pi * math.pow(al_sphere_radius, 3) + + + +# First we make a simple geometry with three cells, (two with material) +sphere_surf_1 = openmc.Sphere(r=20, boundary_type="vacuum") +sphere_surf_2 = openmc.Sphere(r=iron_sphere_radius, z0=10) +sphere_surf_3 = openmc.Sphere(r=al_sphere_radius, z0=-5) + +sphere_region_1 = -sphere_surf_1 & +sphere_surf_2 & +sphere_surf_3 # void space +sphere_region_2 = -sphere_surf_2 +sphere_region_3 = -sphere_surf_3 + +sphere_cell_1 = openmc.Cell(region=sphere_region_1) +sphere_cell_2 = openmc.Cell(region=sphere_region_2,fill = mat_iron) +sphere_cell_3 = openmc.Cell(region=sphere_region_3,fill = mat_aluminum) my_geometry = openmc.Geometry([sphere_cell_1, sphere_cell_2, sphere_cell_3]) +plot = my_geometry.plot(basis='xz') +import matplotlib.pyplot as plt +plt.show() + my_materials = openmc.Materials([mat_iron, mat_aluminum]) pristine_mat_iron = mat_iron.clone() @@ -90,7 +98,7 @@ # Create mesh which will be used for tally regular_mesh = openmc.RegularMesh().from_domain( my_geometry, # the corners of the mesh are being set automatically to surround the geometry - dimension=[10, 10, 10] # 10 + dimension=[2, 2, 2] # 10 ) model_neutron = openmc.Model(my_geometry, my_materials, my_neutron_settings) @@ -135,7 +143,7 @@ flux_in_each_group, micro_xs = openmc.deplete.get_microxs_and_flux( model=model_neutron, domains=regular_mesh, - energies='CCFE-709', # different group structures see this file for all the groups available https://github.com/openmc-dev/openmc/blob/develop/openmc/mgxs/__init__.py + energies=[0,30e6], # different group structures see this file for all the groups available https://github.com/openmc-dev/openmc/blob/develop/openmc/mgxs/__init__.py nuclides=all_nuclides, chain_file=openmc.config['chain_file'] ) @@ -152,13 +160,26 @@ vols = mesh.material_volumes(n_samples = 1000000) +mesh_voxel_material = [] + +for i, entry in enumerate(vols): + print(entry) + for material_volume_tuple in entry: + material = material_volume_tuple[0] + if material != None: + volume_in_cm3 = material_volume_tuple[1] + print(f' {material.id}, {volume_in_cm3}') + # units of material_atom_density are atom/b-cm + + openmc.lib.finalize() + + # # # # constructing the operator, note we pass in the flux and micro xs # operator = openmc.deplete.IndependentOperator().from_nuclides( # volume= # nuclides= - # materials=#TODO find the nuclides in each mesh voxel and their volume fractions, # fluxes=flux_in_each_group, # micros=micro_xs, From 78b13ad724bb9eef93eebb3498ca91e61148374f Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 1 Feb 2024 17:55:10 +0000 Subject: [PATCH 06/13] norming the mixed material --- ..._mesh_based_shut_down_dose_rate_example.py | 35 +++++++++++++++---- 1 file changed, 28 insertions(+), 7 deletions(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index b0f0a35d..8f8d4e8f 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -66,7 +66,7 @@ plot = my_geometry.plot(basis='xz') import matplotlib.pyplot as plt -plt.show() +# plt.show() my_materials = openmc.Materials([mat_iron, mat_aluminum]) @@ -162,15 +162,37 @@ mesh_voxel_material = [] +mat_number_offset = 100 +alls_mats = my_geometry.get_all_materials() +volume_of_material_in_voxel = [] +material_in_voxel = [] for i, entry in enumerate(vols): print(entry) + materials_in_voxel = [] + volumes_in_voxel = [] for material_volume_tuple in entry: material = material_volume_tuple[0] if material != None: volume_in_cm3 = material_volume_tuple[1] - print(f' {material.id}, {volume_in_cm3}') - # units of material_atom_density are atom/b-cm - + print(f' {material.id}, {volume_in_cm3}') + materials_in_voxel.append(alls_mats[material.id]) + volumes_in_voxel.append(volume_in_cm3) + + volume_of_material_in_voxel.append(sum(volumes_in_voxel)) + + print(materials_in_voxel, volumes_in_voxel) + if len(materials_in_voxel)==1: + voxel_mat = materials_in_voxel[0] + + elif len(materials_in_voxel)>1: + # todo check this volume fraction is correct + + norm_fracs = [v*(1/sum(volumes_in_voxel)) for v in volumes_in_voxel] + print('norm_fracs',norm_fracs) + voxel_mat = openmc.Material.mix_materials(materials_in_voxel, norm_fracs) + else: + voxel_mat = openmc.Material() + material_in_voxel.append(voxel_mat) openmc.lib.finalize() @@ -178,9 +200,8 @@ # # # # constructing the operator, note we pass in the flux and micro xs # operator = openmc.deplete.IndependentOperator().from_nuclides( -# volume= -# nuclides= -# materials=#TODO find the nuclides in each mesh voxel and their volume fractions, +# volume=volume_of_material_in_voxel +# materials= # fluxes=flux_in_each_group, # micros=micro_xs, # reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny From 4ede491f4d9ada7b9b4c4dc1795c3c8f47281a0a Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 1 Feb 2024 17:59:30 +0000 Subject: [PATCH 07/13] neutron activation draft looks ok --- ..._mesh_based_shut_down_dose_rate_example.py | 33 ++++++++++--------- 1 file changed, 17 insertions(+), 16 deletions(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index 8f8d4e8f..7b48c08b 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -192,6 +192,8 @@ voxel_mat = openmc.Material.mix_materials(materials_in_voxel, norm_fracs) else: voxel_mat = openmc.Material() + + voxel_mat.volume = sum(volumes_in_voxel) material_in_voxel.append(voxel_mat) openmc.lib.finalize() @@ -199,26 +201,25 @@ # # # # constructing the operator, note we pass in the flux and micro xs -# operator = openmc.deplete.IndependentOperator().from_nuclides( -# volume=volume_of_material_in_voxel -# materials= -# fluxes=flux_in_each_group, -# micros=micro_xs, -# reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny -# reduce_chain_level=5, -# normalization_mode="source-rate" -# ) +operator = openmc.deplete.IndependentOperator( + materials=openmc.Materials(material_in_voxel), + fluxes=[flux[0] for flux in flux_in_each_group], + micros=micro_xs, + reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny + reduce_chain_level=5, + normalization_mode="source-rate" +) -# integrator = openmc.deplete.PredictorIntegrator( -# operator=operator, -# timesteps=timesteps, -# source_rates=source_rates, -# timestep_units='s' -# ) +integrator = openmc.deplete.PredictorIntegrator( + operator=operator, + timesteps=timesteps, + source_rates=source_rates, + timestep_units='s' +) # # this runs the depletion calculations for the timesteps # # this does the neutron activation simulations and produces a depletion_results.h5 file -# integrator.integrate() +integrator.integrate() # # TODO add output dir to integrate command so we don't have to move the file like this # # integrator.integrate(path=statepoints_folder / "neutrons" / "depletion_results.h5") # # PR on openmc is open From e00e26c4a34832b835eacf30c17ae55e35a61890 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Tue, 13 Feb 2024 17:44:53 +0000 Subject: [PATCH 08/13] figure out why clone not working --- ..._mesh_based_shut_down_dose_rate_example.py | 160 ++++++++---------- 1 file changed, 71 insertions(+), 89 deletions(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index 7b48c08b..3197abfc 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -73,13 +73,6 @@ pristine_mat_iron = mat_iron.clone() pristine_mat_aluminium = mat_aluminum.clone() -# gets the cell ids of any depleted cell -activated_cell_ids = [c.id for c in my_geometry.get_all_material_cells().values() if c.fill.depletable] -print("activated_cell_ids", activated_cell_ids) -all_depletable_cells = [c for _, c in my_geometry.get_all_material_cells().items() if c.fill.depletable is True] -print("depletable_cell_ids", activated_cell_ids) -all_depletable_materials = [c.fill for _, c in my_geometry.get_all_material_cells().items() if c.fill.depletable is True] - # 14MeV neutron source that activates material my_source = openmc.IndependentSource() my_source.space = openmc.stats.Point((0, 0, 0)) @@ -149,8 +142,11 @@ ) model_neutron.export_to_model_xml() + import openmc.lib openmc.lib.init() + + mesh = openmc.lib.RegularMesh() mesh.dimension = regular_mesh.dimension mesh.set_parameters( @@ -166,6 +162,8 @@ alls_mats = my_geometry.get_all_materials() volume_of_material_in_voxel = [] material_in_voxel = [] + + for i, entry in enumerate(vols): print(entry) materials_in_voxel = [] @@ -182,7 +180,7 @@ print(materials_in_voxel, volumes_in_voxel) if len(materials_in_voxel)==1: - voxel_mat = materials_in_voxel[0] + voxel_mat = materials_in_voxel[0]#.clone() TODO find out why cloning crashes code elif len(materials_in_voxel)>1: # todo check this volume fraction is correct @@ -191,16 +189,18 @@ print('norm_fracs',norm_fracs) voxel_mat = openmc.Material.mix_materials(materials_in_voxel, norm_fracs) else: + print('no material in voxel') voxel_mat = openmc.Material() voxel_mat.volume = sum(volumes_in_voxel) + voxel_mat.id = i+mat_number_offset material_in_voxel.append(voxel_mat) openmc.lib.finalize() -# # # # constructing the operator, note we pass in the flux and micro xs +# constructing the operator, note we pass in the flux and micro xs operator = openmc.deplete.IndependentOperator( materials=openmc.Materials(material_in_voxel), fluxes=[flux[0] for flux in flux_in_each_group], @@ -219,110 +219,92 @@ # # this runs the depletion calculations for the timesteps # # this does the neutron activation simulations and produces a depletion_results.h5 file -integrator.integrate() -# # TODO add output dir to integrate command so we don't have to move the file like this -# # integrator.integrate(path=statepoints_folder / "neutrons" / "depletion_results.h5") -# # PR on openmc is open -# import os -# os.system(f'mv depletion_results.h5 {statepoints_folder / "neutrons" / "depletion_results.h5"}') +integrator.integrate(path=statepoints_folder / "neutrons" / "depletion_results.h5") # # Now we have done the neutron activation simulations we can start the work needed for the decay gamma simulations. -# my_gamma_settings = openmc.Settings() -# my_gamma_settings.run_mode = "fixed source" -# my_gamma_settings.batches = 100 -# my_gamma_settings.particles = p_particles - +my_gamma_settings = openmc.Settings() +my_gamma_settings.run_mode = "fixed source" +my_gamma_settings.batches = 100 +my_gamma_settings.particles = p_particles # # First we add make dose tally on a regular mesh +# creates a regular mesh that surrounds the geometry +mesh_photon = openmc.RegularMesh().from_domain( + my_geometry, + dimension=[10, 10, 10], # 10 voxels in each axis direction (x, y, z) +) -# # creates a regular mesh that surrounds the geometry -# mesh = openmc.RegularMesh().from_domain( -# my_geometry, -# dimension=[10, 10, 10], # 10 voxels in each axis direction (x, y, z) -# ) - -# # adding a dose tally on a regular mesh -# # AP, PA, LLAT, RLAT, ROT, ISO are ICRP incident dose field directions, AP is front facing -# energies, pSv_cm2 = openmc.data.dose_coefficients(particle="photon", geometry="AP") -# dose_filter = openmc.EnergyFunctionFilter( -# energies, pSv_cm2, interpolation="cubic" # interpolation method recommended by ICRP -# ) -# particle_filter = openmc.ParticleFilter(["photon"]) -# mesh_filter = openmc.MeshFilter(mesh) -# flux_tally = openmc.Tally() -# flux_tally.filters = [mesh_filter, dose_filter, particle_filter] -# flux_tally.scores = ["flux"] -# flux_tally.name = "photon_dose_on_mesh" +# adding a dose tally on a regular mesh +# AP, PA, LLAT, RLAT, ROT, ISO are ICRP incident dose field directions, AP is front facing +energies, pSv_cm2 = openmc.data.dose_coefficients(particle="photon", geometry="AP") +dose_filter = openmc.EnergyFunctionFilter( + energies, pSv_cm2, interpolation="cubic" # interpolation method recommended by ICRP +) +particle_filter = openmc.ParticleFilter(["photon"]) +mesh_filter = openmc.MeshFilter(mesh_photon) +dose_tally = openmc.Tally() +dose_tally.filters = [mesh_filter, dose_filter, particle_filter] +dose_tally.scores = ["flux"] +dose_tally.name = "photon_dose_on_mesh" -# tallies = openmc.Tallies([flux_tally]) +my_gamma_tallies = openmc.Tallies([dose_tally]) -# cells = model_neutron.geometry.get_all_cells() -# activated_cells = [cells[uid] for uid in activated_cell_ids] +cells = model_neutron.geometry.get_all_cells() # # this section makes the photon sources from each active material at each # # timestep and runs the photon simulations -# results = openmc.deplete.Results(statepoints_folder / "neutrons" / "depletion_results.h5") - -# for i_cool in range(1, len(timesteps)): +results = openmc.deplete.Results(statepoints_folder / "neutrons" / "depletion_results.h5") -# # range starts at 1 to skip the first step as that is an irradiation step and there is no -# # decay gamma source from the stable material at that time -# # also there are no decay products in this first timestep for this model +for i_cool in range(1, len(timesteps)): -# photon_sources_for_timestep = [] -# print(f"making photon source for timestep {i_cool}") + # we can loop through the materials in each step + # from the material ID we can get the mesh voxel id + # then we can make a MeshSource + # https://docs.openmc.org/en/develop/pythonapi/generated/openmc.MeshSource.html -# all_activated_materials_in_timestep = [] - -# for activated_cell_id in activated_cell_ids: -# # gets the material id of the material filling the cell -# material_id = cells[activated_cell_id].fill.id - -# # gets the activated material using the material id -# activated_mat = results[i_cool].get_material(str(material_id)) -# # gets the energy and probabilities for the -# energy = activated_mat.get_decay_photon_energy( -# clip_tolerance = 1e-6, # cuts out a small fraction of the very low energy (and hence negligible dose contribution) photons -# units = 'Bq', -# ) -# strength = energy.integral() + # range starts at 1 to skip the first step as that is an irradiation step and there is no + # decay gamma source from the stable material at that time + # also there are no decay products in this first timestep for this model -# if strength > 0.: # only makes sources for -# space = openmc.stats.Box(*cells[activated_cell_id].bounding_box) -# source = openmc.IndependentSource( -# space=space, -# energy=energy, -# particle="photon", -# strength=strength, -# domains=[cells[activated_cell_id]], #this should be the mesh voxel and the material in the mesh voxel -# ) -# photon_sources_for_timestep.append(source) + photon_sources_for_timestep = [] + strengths_for_timestep = [] + print(f"making photon source for timestep {i_cool}") -# my_gamma_settings.source = photon_sources_for_timestep + mat_number_offset + step = results[i_cool] + activated_mat_ids = step.volume.keys() + print(activated_mat_ids) + for activated_mat_id in activated_mat_ids: + # gets the energy and probabilities for the + activated_mat = step.get_material(activated_mat_id) + energy = activated_mat.get_decay_photon_energy( + clip_tolerance = 1e-6, # cuts out a small fraction of the very low energy (and hence negligible dose contribution) photons + units = 'Bq', + ) + strength = energy.integral() - # one should also fill the cells with the activated material - # the activated material contains ALL the nuclides produced during activation - # sphere_cell_2.fill = results[i_cool].get_material("1") - # sphere_cell_3.fill = results[i_cool].get_material("2") - # my_geometry = openmc.Geometry([sphere_cell_1, sphere_cell_2, sphere_cell_3]) + if strength > 0.: + source = openmc.IndependentSource( + energy=energy, + particle="photon", + ) + else: + source = openmc.IndependentSource() # how to make an empty source, source strength is set to 0 - # however it is unlikely that they all appear in your transport cross_sections.xml - # so you could make use of openmc.deplete.Results.export_to_materials to export the modified activated material that - # just contains isotopes that appear in your cross_sections.xml + photon_sources_for_timestep.append(source) + strengths_for_timestep.append(strength) + + mesh_source = openmc.MeshSource(regular_mesh, photon_sources_for_timestep) + mesh_source.strength = strengths_for_timestep - # however in this example we just use the original pristine material my_materials that were cloned earlier - # my_geometry is also the same as the neutron simulation - # pristine_mat_iron.id = 1 - # pristine_mat_aluminium.id =2 - # my_materials = openmc.Materials([pristine_mat_iron, pristine_mat_aluminium]) - # model_gamma = openmc.Model(my_geometry, my_materials, my_gamma_settings, tallies) + model_gamma = openmc.Model(my_geometry, my_materials, my_gamma_settings, my_gamma_tallies) - # model_gamma.run(cwd=statepoints_folder / "photons" / f"photon_at_time_{i_cool}") + model_gamma.run(cwd=statepoints_folder / "photons" / f"photon_at_time_{i_cool}") # pico_to_micro = 1e-6 From 3c2a463bd9d4b15123d90d711567aeb478a6e8e1 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 28 Feb 2024 12:54:47 +0000 Subject: [PATCH 09/13] runing n, p and plotting --- ..._mesh_based_shut_down_dose_rate_example.py | 93 +++++++++++-------- 1 file changed, 53 insertions(+), 40 deletions(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index 3197abfc..aa7350df 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -50,7 +50,7 @@ # First we make a simple geometry with three cells, (two with material) -sphere_surf_1 = openmc.Sphere(r=20, boundary_type="vacuum") +sphere_surf_1 = openmc.Sphere(r=20) sphere_surf_2 = openmc.Sphere(r=iron_sphere_radius, z0=10) sphere_surf_3 = openmc.Sphere(r=al_sphere_radius, z0=-5) @@ -62,10 +62,15 @@ sphere_cell_2 = openmc.Cell(region=sphere_region_2,fill = mat_iron) sphere_cell_3 = openmc.Cell(region=sphere_region_3,fill = mat_aluminum) -my_geometry = openmc.Geometry([sphere_cell_1, sphere_cell_2, sphere_cell_3]) +box = openmc.model.RectangularParallelepiped( + xmin=-20, xmax=20, ymin=-20, ymax=20, zmin=-20, zmax=20, boundary_type="vacuum" +) +sphere_cell_4 = openmc.Cell(region=-box & +sphere_surf_1) + +my_geometry = openmc.Geometry([sphere_cell_1, sphere_cell_2, sphere_cell_3, sphere_cell_4]) -plot = my_geometry.plot(basis='xz') -import matplotlib.pyplot as plt +# plot = my_geometry.plot(basis='xz') +# import matplotlib.pyplot as plt # plt.show() my_materials = openmc.Materials([mat_iron, mat_aluminum]) @@ -91,7 +96,7 @@ # Create mesh which will be used for tally regular_mesh = openmc.RegularMesh().from_domain( my_geometry, # the corners of the mesh are being set automatically to surround the geometry - dimension=[2, 2, 2] # 10 + dimension=[10,10,10] # 10 ) model_neutron = openmc.Model(my_geometry, my_materials, my_neutron_settings) @@ -124,6 +129,7 @@ source_rates = [item[1] for item in timesteps_and_source_rates] model_neutron.export_to_xml(directory=statepoints_folder/ "neutrons") +model_neutron.export_to_xml() all_nuclides = [] for material in my_geometry.get_all_materials().values(): @@ -180,7 +186,7 @@ print(materials_in_voxel, volumes_in_voxel) if len(materials_in_voxel)==1: - voxel_mat = materials_in_voxel[0]#.clone() TODO find out why cloning crashes code + voxel_mat = materials_in_voxel[0].clone() #TODO find out why cloning crashes code elif len(materials_in_voxel)>1: # todo check this volume fraction is correct @@ -194,6 +200,7 @@ voxel_mat.volume = sum(volumes_in_voxel) voxel_mat.id = i+mat_number_offset + # voxel_mat.depletable = True material_in_voxel.append(voxel_mat) openmc.lib.finalize() @@ -206,7 +213,7 @@ fluxes=[flux[0] for flux in flux_in_each_group], micros=micro_xs, reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny - reduce_chain_level=5, + reduce_chain_level=4, normalization_mode="source-rate" ) @@ -219,7 +226,9 @@ # # this runs the depletion calculations for the timesteps # # this does the neutron activation simulations and produces a depletion_results.h5 file -integrator.integrate(path=statepoints_folder / "neutrons" / "depletion_results.h5") +integrator.integrate( + path=statepoints_folder / "neutrons" / "depletion_results.h5" +) # # Now we have done the neutron activation simulations we can start the work needed for the decay gamma simulations. @@ -291,6 +300,7 @@ source = openmc.IndependentSource( energy=energy, particle="photon", + strength=strength ) else: source = openmc.IndependentSource() # how to make an empty source, source strength is set to 0 @@ -298,43 +308,46 @@ photon_sources_for_timestep.append(source) strengths_for_timestep.append(strength) - mesh_source = openmc.MeshSource(regular_mesh, photon_sources_for_timestep) - mesh_source.strength = strengths_for_timestep - + reshaped_photon_source_of_timestep = np.array(photon_sources_for_timestep).reshape(regular_mesh.dimension) + mesh_source = openmc.MeshSource( + regular_mesh, reshaped_photon_source_of_timestep + ) + mesh_source.strength = 1# strengths_for_timestep + my_gamma_settings.source = mesh_source model_gamma = openmc.Model(my_geometry, my_materials, my_gamma_settings, my_gamma_tallies) model_gamma.run(cwd=statepoints_folder / "photons" / f"photon_at_time_{i_cool}") -# pico_to_micro = 1e-6 -# seconds_to_hours = 60*60 +pico_to_micro = 1e-6 +seconds_to_hours = 60*60 # # You may wish to plot the dose tally on a mesh, this package makes it easy to include the geometry with the mesh tally -# from openmc_regular_mesh_plotter import plot_mesh_tally -# for i_cool in range(1, len(timesteps)): -# with openmc.StatePoint(statepoints_folder / "photons" / f"photon_at_time_{i_cool}" / 'statepoint.100.h5') as statepoint: -# photon_tally = statepoint.get_tally(name="photon_dose_on_mesh") - -# # normalising this tally is a little different to other examples as the source strength has been using units of photons per second. -# # tally.mean is in units of pSv-cm3/source photon. -# # as source strength is in photons_per_second this changes units to pSv-/second - -# # multiplication by pico_to_micro converts from (pico) pSv/s to (micro) uSv/s -# # dividing by mesh voxel volume cancles out the cm3 units -# # could do the mesh volume scaling on the plot and vtk functions but doing it here instead -# scaling_factor = (seconds_to_hours * pico_to_micro) / mesh.volumes[0][0][0] - -# plot = plot_mesh_tally( -# tally=photon_tally, -# basis="xz", -# # score='flux', # only one tally so can make use of default here -# value="mean", -# colorbar_kwargs={ -# 'label': "Decay photon dose [µSv/h]", -# }, -# norm=LogNorm(), -# volume_normalization=False, # this is done in the scaling_factor -# scaling_factor=scaling_factor, -# ) -# plot.figure.savefig(f'shut_down_dose_map_timestep_{i_cool}') +from openmc_regular_mesh_plotter import plot_mesh_tally +for i_cool in range(1, len(timesteps)): + with openmc.StatePoint(statepoints_folder / "photons" / f"photon_at_time_{i_cool}" / f'statepoint.{my_gamma_settings.batches}.h5') as statepoint: + photon_tally = statepoint.get_tally(name="photon_dose_on_mesh") + + # normalising this tally is a little different to other examples as the source strength has been using units of photons per second. + # tally.mean is in units of pSv-cm3/source photon. + # as source strength is in photons_per_second this changes units to pSv-/second + + # multiplication by pico_to_micro converts from (pico) pSv/s to (micro) uSv/s + # dividing by mesh voxel volume is not needed as the volume_normalization in the ploting function does this + # could do the mesh volume scaling on the plot and vtk functions but doing it here instead + scaling_factor = (seconds_to_hours * pico_to_micro) + + plot = plot_mesh_tally( + tally=photon_tally, + basis="xz", + # score='flux', # only one tally so can make use of default here + value="mean", + colorbar_kwargs={ + 'label': "Decay photon dose [µSv/h]", + }, + norm=LogNorm(), + volume_normalization=True, # this is done in the scaling_factor + scaling_factor=scaling_factor, + ) + plot.figure.savefig(f'mesh_shut_down_dose_map_timestep_{i_cool}') From 08341650a2b82fa088251a4b41dc90888fca8855 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 28 Feb 2024 13:00:39 +0000 Subject: [PATCH 10/13] deplete true --- .../3_mesh_based_shut_down_dose_rate_example.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index aa7350df..902ffd33 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -200,7 +200,7 @@ voxel_mat.volume = sum(volumes_in_voxel) voxel_mat.id = i+mat_number_offset - # voxel_mat.depletable = True + voxel_mat.depletable = True material_in_voxel.append(voxel_mat) openmc.lib.finalize() From 64070cd02a958249ca0decb22a193ec388d0279b Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 28 Feb 2024 14:02:22 +0000 Subject: [PATCH 11/13] arry size not right --- ..._mesh_based_shut_down_dose_rate_example.py | 100 ++++++++++-------- 1 file changed, 58 insertions(+), 42 deletions(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index 902ffd33..67e9f9cb 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -65,7 +65,7 @@ box = openmc.model.RectangularParallelepiped( xmin=-20, xmax=20, ymin=-20, ymax=20, zmin=-20, zmax=20, boundary_type="vacuum" ) -sphere_cell_4 = openmc.Cell(region=-box & +sphere_surf_1) +sphere_cell_4 = openmc.Cell(region=-box & +sphere_surf_1, fill=mat_aluminum) my_geometry = openmc.Geometry([sphere_cell_1, sphere_cell_2, sphere_cell_3, sphere_cell_4]) @@ -101,32 +101,9 @@ model_neutron = openmc.Model(my_geometry, my_materials, my_neutron_settings) -hour_in_seconds = 60*60 -# This section defines the neutron pulse schedule. -# If the method made use of the CoupledOperator then there would need to be a -# transport simulation for each timestep. However as the IndependentOperator is -# used then just a single transport simulation is done, thus speeding up the -# simulation considerably. -timesteps_and_source_rates = [ - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 1 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 2 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 3 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 4 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 4 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 4 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 5 hour -] -timesteps = [item[0] for item in timesteps_and_source_rates] -source_rates = [item[1] for item in timesteps_and_source_rates] + model_neutron.export_to_xml(directory=statepoints_folder/ "neutrons") model_neutron.export_to_xml() @@ -139,7 +116,7 @@ # print(set(all_nuclides)) # this does perform transport but just to get the flux and micro xs -flux_in_each_group, micro_xs = openmc.deplete.get_microxs_and_flux( +flux_in_each_group, all_micro_xs = openmc.deplete.get_microxs_and_flux( model=model_neutron, domains=regular_mesh, energies=[0,30e6], # different group structures see this file for all the groups available https://github.com/openmc-dev/openmc/blob/develop/openmc/mgxs/__init__.py @@ -153,14 +130,14 @@ openmc.lib.init() -mesh = openmc.lib.RegularMesh() -mesh.dimension = regular_mesh.dimension -mesh.set_parameters( +lib_mesh = openmc.lib.RegularMesh() +lib_mesh.dimension = regular_mesh.dimension +lib_mesh.set_parameters( lower_left=regular_mesh.lower_left, upper_right=regular_mesh.upper_right ) -vols = mesh.material_volumes(n_samples = 1000000) +mesh_material_volumes = lib_mesh.material_volumes(n_samples = 1000000) mesh_voxel_material = [] @@ -169,12 +146,14 @@ volume_of_material_in_voxel = [] material_in_voxel = [] +selective_flux_in_each_group = [] +selective_micro_xs = [] -for i, entry in enumerate(vols): - print(entry) +for i, (mesh_material_volume, micro_xs, flux_in_group) in enumerate(zip(mesh_material_volumes, all_micro_xs, flux_in_each_group) ): + print(mesh_material_volume) materials_in_voxel = [] volumes_in_voxel = [] - for material_volume_tuple in entry: + for material_volume_tuple in mesh_material_volume: material = material_volume_tuple[0] if material != None: volume_in_cm3 = material_volume_tuple[1] @@ -186,22 +165,33 @@ print(materials_in_voxel, volumes_in_voxel) if len(materials_in_voxel)==1: + print('2 materials found') voxel_mat = materials_in_voxel[0].clone() #TODO find out why cloning crashes code + selective_flux_in_each_group.append(flux_in_group) + selective_micro_xs.append(micro_xs) + + voxel_mat.volume = sum(volumes_in_voxel) + voxel_mat.id = i+mat_number_offset + voxel_mat.depletable = True + material_in_voxel.append(voxel_mat) elif len(materials_in_voxel)>1: # todo check this volume fraction is correct norm_fracs = [v*(1/sum(volumes_in_voxel)) for v in volumes_in_voxel] - print('norm_fracs',norm_fracs) + print('1 material found') voxel_mat = openmc.Material.mix_materials(materials_in_voxel, norm_fracs) + selective_flux_in_each_group.append(flux_in_group) + selective_micro_xs.append(micro_xs) + + voxel_mat.volume = sum(volumes_in_voxel) + voxel_mat.id = i+mat_number_offset + voxel_mat.depletable = True + material_in_voxel.append(voxel_mat) + else: print('no material in voxel') - voxel_mat = openmc.Material() - - voxel_mat.volume = sum(volumes_in_voxel) - voxel_mat.id = i+mat_number_offset - voxel_mat.depletable = True - material_in_voxel.append(voxel_mat) + openmc.lib.finalize() @@ -210,13 +200,39 @@ # constructing the operator, note we pass in the flux and micro xs operator = openmc.deplete.IndependentOperator( materials=openmc.Materials(material_in_voxel), - fluxes=[flux[0] for flux in flux_in_each_group], - micros=micro_xs, + fluxes=[flux[0] for flux in selective_flux_in_each_group], + micros=selective_micro_xs, reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny reduce_chain_level=4, normalization_mode="source-rate" ) +# This section defines the neutron pulse schedule. +# If the method made use of the CoupledOperator then there would need to be a +# transport simulation for each timestep. However as the IndependentOperator is +# used then just a single transport simulation is done, thus speeding up the +# simulation considerably. +hour_in_seconds = 60*60 +timesteps_and_source_rates = [ + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 1 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 2 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 3 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 4 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 4 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 4 hour + (1, 1e18), # 1 second + (hour_in_seconds, 0), # 5 hour +] + +timesteps = [item[0] for item in timesteps_and_source_rates] +source_rates = [item[1] for item in timesteps_and_source_rates] + integrator = openmc.deplete.PredictorIntegrator( operator=operator, timesteps=timesteps, From 6c62bf577f40f635982cf3a53b3c429357dc6b0d Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Thu, 9 May 2024 10:31:47 +0100 Subject: [PATCH 12/13] simplified example --- ..._mesh_based_shut_down_dose_rate_example.py | 276 +++++++----------- 1 file changed, 100 insertions(+), 176 deletions(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index 67e9f9cb..f80948d3 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -12,9 +12,8 @@ # the chain file was downloaded with # pip install openmc_data # download_endf_chain -r b8.0 -# openmc.config['chain_file'] = '/nuclear_data/chain-endf-b8.0.xml' -# openmc.config['cross_sections'] = 'cross_sections.xml' -openmc.config['chain_file'] = '/home/jshimwell/ENDF-B-VIII.0-NNDC/chain-nndc-b8.0.xml' +openmc.config['chain_file'] = '/nuclear_data/chain-endf-b8.0.xml' +openmc.config['cross_sections'] = '/nuclear_data/cross_sections.xml' # a few user settings # Set up the folders to save all the data in @@ -50,34 +49,29 @@ # First we make a simple geometry with three cells, (two with material) -sphere_surf_1 = openmc.Sphere(r=20) -sphere_surf_2 = openmc.Sphere(r=iron_sphere_radius, z0=10) -sphere_surf_3 = openmc.Sphere(r=al_sphere_radius, z0=-5) +sphere_surf_1 = openmc.Sphere(r=iron_sphere_radius, z0=10) +sphere_surf_2 = openmc.Sphere(r=al_sphere_radius, z0=-5) -sphere_region_1 = -sphere_surf_1 & +sphere_surf_2 & +sphere_surf_3 # void space +sphere_region_1 = -sphere_surf_1 sphere_region_2 = -sphere_surf_2 -sphere_region_3 = -sphere_surf_3 -sphere_cell_1 = openmc.Cell(region=sphere_region_1) +sphere_cell_1 = openmc.Cell(region=sphere_region_1,fill = mat_aluminum) sphere_cell_2 = openmc.Cell(region=sphere_region_2,fill = mat_iron) -sphere_cell_3 = openmc.Cell(region=sphere_region_3,fill = mat_aluminum) box = openmc.model.RectangularParallelepiped( xmin=-20, xmax=20, ymin=-20, ymax=20, zmin=-20, zmax=20, boundary_type="vacuum" ) -sphere_cell_4 = openmc.Cell(region=-box & +sphere_surf_1, fill=mat_aluminum) +box_cell = openmc.Cell(region=-box & +sphere_surf_1, fill=mat_aluminum) -my_geometry = openmc.Geometry([sphere_cell_1, sphere_cell_2, sphere_cell_3, sphere_cell_4]) +my_geometry = openmc.Geometry([sphere_cell_1, sphere_cell_2, box_cell]) +# to plot the geometry uncomment this line # plot = my_geometry.plot(basis='xz') # import matplotlib.pyplot as plt # plt.show() my_materials = openmc.Materials([mat_iron, mat_aluminum]) -pristine_mat_iron = mat_iron.clone() -pristine_mat_aluminium = mat_aluminum.clone() - # 14MeV neutron source that activates material my_source = openmc.IndependentSource() my_source.space = openmc.stats.Point((0, 0, 0)) @@ -93,7 +87,7 @@ my_neutron_settings.source = my_source my_neutron_settings.photon_transport = False -# Create mesh which will be used for tally +# Create mesh which will be used for material segmentation and activation and gamma sources regular_mesh = openmc.RegularMesh().from_domain( my_geometry, # the corners of the mesh are being set automatically to surround the geometry dimension=[10,10,10] # 10 @@ -101,10 +95,6 @@ model_neutron = openmc.Model(my_geometry, my_materials, my_neutron_settings) - - - - model_neutron.export_to_xml(directory=statepoints_folder/ "neutrons") model_neutron.export_to_xml() @@ -113,95 +103,25 @@ for nuclide in material.get_nuclides(): if nuclide not in all_nuclides: all_nuclides.append(nuclide) -# print(set(all_nuclides)) # this does perform transport but just to get the flux and micro xs -flux_in_each_group, all_micro_xs = openmc.deplete.get_microxs_and_flux( +print(f'running neutron transport to activate materials') +flux_in_each_mesh_voxel, all_micro_xs = openmc.deplete.get_microxs_and_flux( model=model_neutron, domains=regular_mesh, - energies=[0,30e6], # different group structures see this file for all the groups available https://github.com/openmc-dev/openmc/blob/develop/openmc/mgxs/__init__.py + energies=[0,30e6], nuclides=all_nuclides, chain_file=openmc.config['chain_file'] ) -model_neutron.export_to_model_xml() - -import openmc.lib -openmc.lib.init() - - -lib_mesh = openmc.lib.RegularMesh() -lib_mesh.dimension = regular_mesh.dimension -lib_mesh.set_parameters( - lower_left=regular_mesh.lower_left, - upper_right=regular_mesh.upper_right -) - -mesh_material_volumes = lib_mesh.material_volumes(n_samples = 1000000) - -mesh_voxel_material = [] - -mat_number_offset = 100 -alls_mats = my_geometry.get_all_materials() -volume_of_material_in_voxel = [] -material_in_voxel = [] - -selective_flux_in_each_group = [] -selective_micro_xs = [] - -for i, (mesh_material_volume, micro_xs, flux_in_group) in enumerate(zip(mesh_material_volumes, all_micro_xs, flux_in_each_group) ): - print(mesh_material_volume) - materials_in_voxel = [] - volumes_in_voxel = [] - for material_volume_tuple in mesh_material_volume: - material = material_volume_tuple[0] - if material != None: - volume_in_cm3 = material_volume_tuple[1] - print(f' {material.id}, {volume_in_cm3}') - materials_in_voxel.append(alls_mats[material.id]) - volumes_in_voxel.append(volume_in_cm3) - - volume_of_material_in_voxel.append(sum(volumes_in_voxel)) - - print(materials_in_voxel, volumes_in_voxel) - if len(materials_in_voxel)==1: - print('2 materials found') - voxel_mat = materials_in_voxel[0].clone() #TODO find out why cloning crashes code - selective_flux_in_each_group.append(flux_in_group) - selective_micro_xs.append(micro_xs) - - voxel_mat.volume = sum(volumes_in_voxel) - voxel_mat.id = i+mat_number_offset - voxel_mat.depletable = True - material_in_voxel.append(voxel_mat) - - elif len(materials_in_voxel)>1: - # todo check this volume fraction is correct - - norm_fracs = [v*(1/sum(volumes_in_voxel)) for v in volumes_in_voxel] - print('1 material found') - voxel_mat = openmc.Material.mix_materials(materials_in_voxel, norm_fracs) - selective_flux_in_each_group.append(flux_in_group) - selective_micro_xs.append(micro_xs) - - voxel_mat.volume = sum(volumes_in_voxel) - voxel_mat.id = i+mat_number_offset - voxel_mat.depletable = True - material_in_voxel.append(voxel_mat) - - else: - print('no material in voxel') - - -openmc.lib.finalize() - - +print(f'running transport to sample within the mesh and get material fractions') +mixed_materials_in_each_mesh_voxel = regular_mesh.get_homogenized_materials(model_neutron, n_samples=1_000_000) # constructing the operator, note we pass in the flux and micro xs operator = openmc.deplete.IndependentOperator( - materials=openmc.Materials(material_in_voxel), - fluxes=[flux[0] for flux in selective_flux_in_each_group], - micros=selective_micro_xs, + materials=openmc.Materials(mixed_materials_in_each_mesh_voxel), + fluxes=[flux[0] for flux in flux_in_each_mesh_voxel], + micros=all_micro_xs, reduce_chain=True, # reduced to only the isotopes present in depletable materials and their possible progeny reduce_chain_level=4, normalization_mode="source-rate" @@ -212,22 +132,22 @@ # transport simulation for each timestep. However as the IndependentOperator is # used then just a single transport simulation is done, thus speeding up the # simulation considerably. +# This pulse schedule has 1 second pulses of neutrons then a days of cooling steps in 1 hour blocks hour_in_seconds = 60*60 timesteps_and_source_rates = [ - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 1 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 2 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 3 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 4 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 4 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 4 hour - (1, 1e18), # 1 second - (hour_in_seconds, 0), # 5 hour + (1, 1e18), # 1 second of 1e18 neutrons/second + (2*hour_in_seconds, 0), # 2 hours after shut down + (2*hour_in_seconds, 0), # 4 hours after shut down + (2*hour_in_seconds, 0), # 6 hours after shut down + (2*hour_in_seconds, 0), # 8 hours after shut down + (2*hour_in_seconds, 0), # 10 hours after shut down + (2*hour_in_seconds, 0), # 12 hours after shut down + (2*hour_in_seconds, 0), # 14 hours after shut down + (2*hour_in_seconds, 0), # 16 hours after shut down + (2*hour_in_seconds, 0), # 18 hours after shut down + (2*hour_in_seconds, 0), # 20 hours after shut down + (2*hour_in_seconds, 0), # 22 hours after shut down + (2*hour_in_seconds, 0), # 24 hours after shut down ] timesteps = [item[0] for item in timesteps_and_source_rates] @@ -240,29 +160,29 @@ timestep_units='s' ) -# # this runs the depletion calculations for the timesteps -# # this does the neutron activation simulations and produces a depletion_results.h5 file +# this runs the depletion calculations for the timesteps +# this does the neutron activation simulations and produces a depletion_results.h5 file integrator.integrate( path=statepoints_folder / "neutrons" / "depletion_results.h5" ) -# # Now we have done the neutron activation simulations we can start the work needed for the decay gamma simulations. +# Now we have done the neutron activation simulations we can start the work needed for the decay gamma simulations. my_gamma_settings = openmc.Settings() my_gamma_settings.run_mode = "fixed source" my_gamma_settings.batches = 100 my_gamma_settings.particles = p_particles -# # First we add make dose tally on a regular mesh +# First we add make dose tally on a regular mesh -# creates a regular mesh that surrounds the geometry +# # creates a regular mesh that surrounds the geometry mesh_photon = openmc.RegularMesh().from_domain( my_geometry, - dimension=[10, 10, 10], # 10 voxels in each axis direction (x, y, z) + dimension=[20, 20, 20], # 20 voxels in each axis direction (x, y, z) ) -# adding a dose tally on a regular mesh -# AP, PA, LLAT, RLAT, ROT, ISO are ICRP incident dose field directions, AP is front facing +# # adding a dose tally on a regular mesh +# # AP, PA, LLAT, RLAT, ROT, ISO are ICRP incident dose field directions, AP is front facing energies, pSv_cm2 = openmc.data.dose_coefficients(particle="photon", geometry="AP") dose_filter = openmc.EnergyFunctionFilter( energies, pSv_cm2, interpolation="cubic" # interpolation method recommended by ICRP @@ -278,92 +198,96 @@ cells = model_neutron.geometry.get_all_cells() -# # this section makes the photon sources from each active material at each -# # timestep and runs the photon simulations results = openmc.deplete.Results(statepoints_folder / "neutrons" / "depletion_results.h5") +# this section makes the photon sources from each active material at each +# timestep and runs the photon simulations +# range starts at 1 to skip the first step as that is an irradiation step and there is no for i_cool in range(1, len(timesteps)): - # we can loop through the materials in each step # from the material ID we can get the mesh voxel id # then we can make a MeshSource # https://docs.openmc.org/en/develop/pythonapi/generated/openmc.MeshSource.html - - # range starts at 1 to skip the first step as that is an irradiation step and there is no # decay gamma source from the stable material at that time # also there are no decay products in this first timestep for this model - - photon_sources_for_timestep = [] - strengths_for_timestep = [] - print(f"making photon source for timestep {i_cool}") - - - mat_number_offset - - step = results[i_cool] - activated_mat_ids = step.volume.keys() - print(activated_mat_ids) - for activated_mat_id in activated_mat_ids: - # gets the energy and probabilities for the - activated_mat = step.get_material(activated_mat_id) - energy = activated_mat.get_decay_photon_energy( - clip_tolerance = 1e-6, # cuts out a small fraction of the very low energy (and hence negligible dose contribution) photons - units = 'Bq', - ) - strength = energy.integral() - - if strength > 0.: - source = openmc.IndependentSource( - energy=energy, - particle="photon", - strength=strength - ) - else: - source = openmc.IndependentSource() # how to make an empty source, source strength is set to 0 - - photon_sources_for_timestep.append(source) - strengths_for_timestep.append(strength) - - reshaped_photon_source_of_timestep = np.array(photon_sources_for_timestep).reshape(regular_mesh.dimension) - mesh_source = openmc.MeshSource( - regular_mesh, reshaped_photon_source_of_timestep + photon_sources_for_timestep = [] + strengths_for_timestep = [] + print(f"making photon source for timestep {i_cool}") + step = results[i_cool] + # activated_mat_ids = step.volume.keys() + activated_mat_ids = step.index_mat + # print(activated_mat_ids) + cumulative_strength_for_time_step = 0 # in Bq + for activated_mat_id in activated_mat_ids: + # gets the energy and probabilities for the + activated_mat = step.get_material(activated_mat_id) + energy = activated_mat.get_decay_photon_energy( + clip_tolerance = 1e-6, # cuts out a small fraction of the very low energy (and hence negligible dose contribution) photons + units = 'Bq', ) - mesh_source.strength = 1# strengths_for_timestep + strength = energy.integral() + cumulative_strength_for_time_step = cumulative_strength_for_time_step +strength + if strength > 0.: + source = openmc.IndependentSource( + energy=energy, + particle="photon", + strength=strength + ) + else: + print('source has no gammas') + source = openmc.IndependentSource() # how to make an empty source, source strength is set to 0 + + photon_sources_for_timestep.append(source) + strengths_for_timestep.append(strength) + + mesh_source = openmc.MeshSource( + regular_mesh, photon_sources_for_timestep + ) - my_gamma_settings.source = mesh_source - model_gamma = openmc.Model(my_geometry, my_materials, my_gamma_settings, my_gamma_tallies) + # you have options for the normalization of the source. + # you could set the mesh_source.strength to the total Bq of all the sources in that time step + # mesh_source.strength = cumulative_strength_for_time_step + # then use mesh_source.normalize_source_strengths() to update all element source strengths such that they sum to 1.0. + # or + # you can leave it so the individual sources have their own strength in Bq + # perhaps best to experiment here and check the answers, do let me know if you find one option better than the others - model_gamma.run(cwd=statepoints_folder / "photons" / f"photon_at_time_{i_cool}") + my_gamma_settings.source = mesh_source + model_gamma = openmc.Model(my_geometry, my_materials, my_gamma_settings, my_gamma_tallies) + print(f'running gamma transport on stimestep {i_cool}') + model_gamma.run(cwd=statepoints_folder / "photons" / f"photon_at_time_{i_cool}") +# this part post processes the results to get a dose map for each time step pico_to_micro = 1e-6 seconds_to_hours = 60*60 -# # You may wish to plot the dose tally on a mesh, this package makes it easy to include the geometry with the mesh tally +# You may wish to plot the dose tally on a mesh, this package makes it easy to include the geometry with the mesh tally from openmc_regular_mesh_plotter import plot_mesh_tally -for i_cool in range(1, len(timesteps)): +for i_cool in range(1, len(timesteps)): # skipping the first depletion step as we just want the part where the machine is off for the shut down dose rate with openmc.StatePoint(statepoints_folder / "photons" / f"photon_at_time_{i_cool}" / f'statepoint.{my_gamma_settings.batches}.h5') as statepoint: photon_tally = statepoint.get_tally(name="photon_dose_on_mesh") - - # normalising this tally is a little different to other examples as the source strength has been using units of photons per second. + # normalizing this tally is a little different to other examples as the source strength has been using units of photons per second. # tally.mean is in units of pSv-cm3/source photon. # as source strength is in photons_per_second this changes units to pSv-/second - # multiplication by pico_to_micro converts from (pico) pSv/s to (micro) uSv/s - # dividing by mesh voxel volume is not needed as the volume_normalization in the ploting function does this - # could do the mesh volume scaling on the plot and vtk functions but doing it here instead + # dividing by mesh voxel volume is not needed as the volume_normalization in the plotting function does this scaling_factor = (seconds_to_hours * pico_to_micro) - + print('max',(max(photon_tally.mean.flatten())*scaling_factor)/mesh_photon.volumes[0][0][0]) + print('min',(min(photon_tally.mean.flatten())*scaling_factor)/mesh_photon.volumes[0][0][0]) plot = plot_mesh_tally( tally=photon_tally, basis="xz", - # score='flux', # only one tally so can make use of default here + score='flux', # only one tally so could leave this unspecified value="mean", colorbar_kwargs={ 'label': "Decay photon dose [µSv/h]", }, - norm=LogNorm(), - volume_normalization=True, # this is done in the scaling_factor + norm=LogNorm(), # TODO find the bounds automatically in a loop above this section + volume_normalization=True, # this divides by voxel volume which is not done in the scaling_factor scaling_factor=scaling_factor, + outline=True, + geometry = my_geometry ) - plot.figure.savefig(f'mesh_shut_down_dose_map_timestep_{i_cool}') + plot.title.set_text(f"timestep {sum(timesteps[1:i_cool])/(60*60)} hours after shut down") + plot.figure.savefig(f'mesh_shut_down_dose_map_timestep_{str(i_cool).zfill(3)}') From bb785106659de776416745b4bf0ad56a27a69e94 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Mon, 20 May 2024 21:30:38 +0100 Subject: [PATCH 13/13] retrigger ci --- .../3_mesh_based_shut_down_dose_rate_example.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py index f80948d3..1465e380 100644 --- a/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py +++ b/tasks/task_11_CSG_shut_down_dose_tallies/3_mesh_based_shut_down_dose_rate_example.py @@ -181,8 +181,8 @@ dimension=[20, 20, 20], # 20 voxels in each axis direction (x, y, z) ) -# # adding a dose tally on a regular mesh -# # AP, PA, LLAT, RLAT, ROT, ISO are ICRP incident dose field directions, AP is front facing +# adding a dose tally on a regular mesh +# AP, PA, LLAT, RLAT, ROT, ISO are ICRP incident dose field directions, AP is front facing energies, pSv_cm2 = openmc.data.dose_coefficients(particle="photon", geometry="AP") dose_filter = openmc.EnergyFunctionFilter( energies, pSv_cm2, interpolation="cubic" # interpolation method recommended by ICRP