From dce2c74158098964175fcd5de9e2d32f4de3ecea Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 8 May 2025 06:22:18 -0500 Subject: [PATCH 01/29] Function to clone materials for mesh-material --- openmc/deplete/microxs.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index df78f8a267c..684dbb4182a 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -225,6 +225,21 @@ def get_microxs_and_flux( return fluxes, micros +def get_material_copies(model: openmc.Model, mmv: openmc.MeshMaterialVolumes): + mat_ids = mmv._materials[mmv._materials > -1] + elems, _ = np.where(mmv._materials > -1) + # TODO: Handle DAGMC case + material_dict = model.geometry.get_all_materials() + materials = [] + for elem, mat_id in zip(elems, mat_ids): + mat = material_dict[mat_id] + new_mat = mat.clone() + new_mat.depletable = True + new_mat.name = 'Element {elem}, Material {mat_id}' + materials.append(new_mat) + return materials + + class MicroXS: """Microscopic cross section data for use in transport-independent depletion. From d9c448d1c9ebdff59d942928bf3d12f7151dadca Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 24 May 2025 15:40:09 -0500 Subject: [PATCH 02/29] Start writing R2S module --- openmc/deplete/microxs.py | 15 -------- openmc/deplete/r2s.py | 81 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 81 insertions(+), 15 deletions(-) create mode 100644 openmc/deplete/r2s.py diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index 684dbb4182a..df78f8a267c 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -225,21 +225,6 @@ def get_microxs_and_flux( return fluxes, micros -def get_material_copies(model: openmc.Model, mmv: openmc.MeshMaterialVolumes): - mat_ids = mmv._materials[mmv._materials > -1] - elems, _ = np.where(mmv._materials > -1) - # TODO: Handle DAGMC case - material_dict = model.geometry.get_all_materials() - materials = [] - for elem, mat_id in zip(elems, mat_ids): - mat = material_dict[mat_id] - new_mat = mat.clone() - new_mat.depletable = True - new_mat.name = 'Element {elem}, Material {mat_id}' - materials.append(new_mat) - return materials - - class MicroXS: """Microscopic cross section data for use in transport-independent depletion. diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py new file mode 100644 index 00000000000..9a5b3361fba --- /dev/null +++ b/openmc/deplete/r2s.py @@ -0,0 +1,81 @@ +import numpy as np +import openmc +from .results import Results + + +def get_material_copies(model: openmc.Model, mmv: openmc.MeshMaterialVolumes): + mat_ids = mmv._materials[mmv._materials > -1] + volumes = mmv._volumes[mmv._materials > -1] + elems, _ = np.where(mmv._materials > -1) + # TODO: Handle DAGMC case + material_dict = model.geometry.get_all_materials() + materials = [] + for elem, mat_id, vol in zip(elems, mat_ids, volumes): + mat = material_dict[mat_id] + new_mat = mat.clone() + new_mat.depletable = True + new_mat.name = f'Element {elem}, Material {mat_id}' + new_mat.volume = vol + materials.append(new_mat) + return materials + + +# TODO: put as part of a R2S class +def get_mesh_source( + model: openmc.Model, + mesh: openmc.MeshBase, + materials: list[openmc.Material], + results: Results, + mat_vols: openmc.MeshMaterialVolumes, +): + """Create decay photon source for a mesh-based R2S calculation. + + This function creates N :class:`MeshSource` objects where N is the maximum + number of unique materials that appears in a single mesh element.""" + mat_dict = model.geometry.get_all_materials() + + # Some MeshSource objects will have empty positions; create a "null source" + # that is used for this case + null_source = openmc.IndependentSource(particle='photon', strength=0.0) + + # Lists to hold sources and rejection domains for each MeshSource. That is, + # the length of source_lists and domain_lists is N. + source_lists = [] + domain_lists = [] + + # Index in the overall list of activated materials + index_mat = 0 + + # Total number of mesh elements + n_elements = mat_vols.num_elements + + for index_elem in range(n_elements): + for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): + # Check whether a new MeshSource object is needed + if j >= len(source_lists): + source_lists.append([null_source]*n_elements) + domain_lists.append([]) + + # Get activated material composition + original_mat = materials[index_mat] + activated_mat = results[-1].get_material(str(original_mat.id)) + + # Create decay photon source source; note that domain rejection is + # not used here because the MeshSource implementation ignores + # spatial constraints at the level of an individual element. + energy = activated_mat.get_decay_photon_energy() + source_lists[j][index_elem] = openmc.IndependentSource( + energy=energy, + particle='photon', + strength=energy.integral(), + ) + domain_lists[j].append(mat_dict[mat_id]) + + # Increment index of activated material + index_mat += 1 + + # Return list of mesh sources + return [ + openmc.MeshSource(mesh, sources, constraints={'domains': domains}) + for sources, domains in zip(source_lists, domain_lists) + ] From 763f669b214f34e8bda02ca4e8343057230824f9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 8 Jun 2025 21:00:26 -0500 Subject: [PATCH 03/29] Use domains constraint on IndependentSource objects --- openmc/deplete/r2s.py | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 9a5b3361fba..e6d12859ce7 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np import openmc from .results import Results @@ -27,7 +29,7 @@ def get_mesh_source( materials: list[openmc.Material], results: Results, mat_vols: openmc.MeshMaterialVolumes, -): +) -> list[openmc.MeshSource]: """Create decay photon source for a mesh-based R2S calculation. This function creates N :class:`MeshSource` objects where N is the maximum @@ -38,10 +40,8 @@ def get_mesh_source( # that is used for this case null_source = openmc.IndependentSource(particle='photon', strength=0.0) - # Lists to hold sources and rejection domains for each MeshSource. That is, - # the length of source_lists and domain_lists is N. + # List to hold sources for each MeshSource (length = N) source_lists = [] - domain_lists = [] # Index in the overall list of activated materials index_mat = 0 @@ -51,31 +51,29 @@ def get_mesh_source( for index_elem in range(n_elements): for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): + # Skip void volume + if mat_id is None: + continue + # Check whether a new MeshSource object is needed if j >= len(source_lists): source_lists.append([null_source]*n_elements) - domain_lists.append([]) # Get activated material composition original_mat = materials[index_mat] activated_mat = results[-1].get_material(str(original_mat.id)) - # Create decay photon source source; note that domain rejection is - # not used here because the MeshSource implementation ignores - # spatial constraints at the level of an individual element. + # Create decay photon source source energy = activated_mat.get_decay_photon_energy() source_lists[j][index_elem] = openmc.IndependentSource( energy=energy, particle='photon', strength=energy.integral(), + constraints={'domains': [mat_dict[mat_id]]} ) - domain_lists[j].append(mat_dict[mat_id]) # Increment index of activated material index_mat += 1 # Return list of mesh sources - return [ - openmc.MeshSource(mesh, sources, constraints={'domains': domains}) - for sources, domains in zip(source_lists, domain_lists) - ] + return [openmc.MeshSource(mesh, sources) for sources in source_lists] From adef0eb2213f3ac913d601fa5c012b210d97d5b2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 9 Jun 2025 23:01:10 -0500 Subject: [PATCH 04/29] Allow use of r2s functions for DAGMC geometries --- openmc/deplete/r2s.py | 75 ++++++++++++++++++++++++++++++++++++++----- openmc/mesh.py | 33 +++++++++++-------- 2 files changed, 87 insertions(+), 21 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index e6d12859ce7..414bcb37c93 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -3,14 +3,46 @@ import numpy as np import openmc from .results import Results - - -def get_material_copies(model: openmc.Model, mmv: openmc.MeshMaterialVolumes): +from ..mesh import _get_all_materials + + +def get_activation_materials( + model: openmc.Model, + mmv: openmc.MeshMaterialVolumes +) -> list[openmc.Material]: + """Get a list of activation materials for each mesh element/material. + + When performing a mesh-based R2S calculation, a unique material is needed + for each activation region, which is a combination of a mesh element and a + material within that mesh element. This function generates a list of such + materials, each with a unique name and volume corresponding to the mesh + element and material. + + Parameters + ---------- + model : openmc.Model + The OpenMC model containing the materials. + mmv : openmc.MeshMaterialVolumes + The mesh material volumes object containing the materials and their + volumes for each mesh element. + + Returns + ------- + list of openmc.Material + A list of materials, each corresponding to a unique mesh element and + material combination. + + """ + # Get the material ID, volume, and element index for each element-material + # combination mat_ids = mmv._materials[mmv._materials > -1] volumes = mmv._volumes[mmv._materials > -1] elems, _ = np.where(mmv._materials > -1) - # TODO: Handle DAGMC case - material_dict = model.geometry.get_all_materials() + + # Get all materials in the model + material_dict = _get_all_materials(model) + + # Create a new activation material for each element-material combination materials = [] for elem, mat_id, vol in zip(elems, mat_ids, volumes): mat = material_dict[mat_id] @@ -19,11 +51,12 @@ def get_material_copies(model: openmc.Model, mmv: openmc.MeshMaterialVolumes): new_mat.name = f'Element {elem}, Material {mat_id}' new_mat.volume = vol materials.append(new_mat) + return materials # TODO: put as part of a R2S class -def get_mesh_source( +def get_decay_photon_source( model: openmc.Model, mesh: openmc.MeshBase, materials: list[openmc.Material], @@ -33,8 +66,34 @@ def get_mesh_source( """Create decay photon source for a mesh-based R2S calculation. This function creates N :class:`MeshSource` objects where N is the maximum - number of unique materials that appears in a single mesh element.""" - mat_dict = model.geometry.get_all_materials() + number of unique materials that appears in a single mesh element. For each + mesh element-material combination, and IndependentSource instance is created + with a spatial constraint limited the sampled decay photons to the correct + region. + + Parameters + ---------- + model : openmc.Model + The OpenMC model containing the geometry and materials. + mesh : openmc.MeshBase + The mesh object defining the spatial regions over which activation is + performed. + materials : list of openmc.Material + List of materials that are activated in the model. + results : openmc.deplete.Results + The results object containing the depletion results for the materials. + mat_vols : openmc.MeshMaterialVolumes + The mesh material volumes object containing the materials and their + volumes for each mesh element. + + Returns + ------- + list of openmc.MeshSource + A list of MeshSource objects, each containing IndependentSource + instances for the decay photons in the corresponding mesh element. + + """ + mat_dict = _get_all_materials(model) # Some MeshSource objects will have empty positions; create a "null source" # that is used for this case diff --git a/openmc/mesh.py b/openmc/mesh.py index 33970288440..4a478aa85be 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -317,20 +317,10 @@ def get_homogenized_materials( vols = self.material_volumes(model, n_samples, **kwargs) mat_volume_by_element = [vols.by_element(i) for i in range(vols.num_elements)] - # Create homogenized material for each element - materials = model.geometry.get_all_materials() - - # Account for materials in DAGMC universes - # TODO: This should really get incorporated in lower-level calls to - # get_all_materials, but right now it requires information from the - # Model object - for cell in model.geometry.get_all_cells().values(): - if isinstance(cell.fill, openmc.DAGMCUniverse): - names = cell.fill.material_names - materials.update({ - mat.id: mat for mat in model.materials if mat.name in names - }) + # Get dictionary of all materials + materials = _get_all_materials(model) + # Create homogenized material for each element homogenized_materials = [] for mat_volume_list in mat_volume_by_element: material_ids, volumes = [list(x) for x in zip(*mat_volume_list)] @@ -2869,6 +2859,23 @@ def _read_meshes(elem): return out +# TODO: This should really get incorporated in lower-level calls to +# get_all_materials, but right now it requires information from the Model object +def _get_all_materials(model: openmc.Model) -> dict[int, openmc.Material]: + # Get all materials from the Geometry object + materials = model.geometry.get_all_materials() + + # Account for materials in DAGMC universes + for cell in model.geometry.get_all_cells().values(): + if isinstance(cell.fill, openmc.DAGMCUniverse): + names = cell.fill.material_names + materials.update({ + mat.id: mat for mat in model.materials if mat.name in names + }) + + return materials + + # hexahedron element connectivity # lower-k connectivity offsets _HEX_VERTEX_CONN = ((0, 0, 0), From 7820a6a12d12f836be8efc1bb6142d3a7de3396d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 10 Jun 2025 09:28:08 -0500 Subject: [PATCH 05/29] Start R2SManager class --- openmc/deplete/r2s.py | 260 +++++++++++++++++++++--------------------- 1 file changed, 133 insertions(+), 127 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 414bcb37c93..3814aa693bf 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -6,133 +6,139 @@ from ..mesh import _get_all_materials -def get_activation_materials( - model: openmc.Model, - mmv: openmc.MeshMaterialVolumes -) -> list[openmc.Material]: - """Get a list of activation materials for each mesh element/material. - - When performing a mesh-based R2S calculation, a unique material is needed - for each activation region, which is a combination of a mesh element and a - material within that mesh element. This function generates a list of such - materials, each with a unique name and volume corresponding to the mesh - element and material. - - Parameters - ---------- - model : openmc.Model - The OpenMC model containing the materials. - mmv : openmc.MeshMaterialVolumes - The mesh material volumes object containing the materials and their - volumes for each mesh element. - - Returns - ------- - list of openmc.Material - A list of materials, each corresponding to a unique mesh element and - material combination. +class R2SManager: + """Manager for R2S calculations. - """ - # Get the material ID, volume, and element index for each element-material - # combination - mat_ids = mmv._materials[mmv._materials > -1] - volumes = mmv._volumes[mmv._materials > -1] - elems, _ = np.where(mmv._materials > -1) - - # Get all materials in the model - material_dict = _get_all_materials(model) - - # Create a new activation material for each element-material combination - materials = [] - for elem, mat_id, vol in zip(elems, mat_ids, volumes): - mat = material_dict[mat_id] - new_mat = mat.clone() - new_mat.depletable = True - new_mat.name = f'Element {elem}, Material {mat_id}' - new_mat.volume = vol - materials.append(new_mat) - - return materials - - -# TODO: put as part of a R2S class -def get_decay_photon_source( - model: openmc.Model, - mesh: openmc.MeshBase, - materials: list[openmc.Material], - results: Results, - mat_vols: openmc.MeshMaterialVolumes, -) -> list[openmc.MeshSource]: - """Create decay photon source for a mesh-based R2S calculation. - - This function creates N :class:`MeshSource` objects where N is the maximum - number of unique materials that appears in a single mesh element. For each - mesh element-material combination, and IndependentSource instance is created - with a spatial constraint limited the sampled decay photons to the correct - region. - - Parameters - ---------- - model : openmc.Model - The OpenMC model containing the geometry and materials. - mesh : openmc.MeshBase - The mesh object defining the spatial regions over which activation is - performed. - materials : list of openmc.Material - List of materials that are activated in the model. - results : openmc.deplete.Results - The results object containing the depletion results for the materials. - mat_vols : openmc.MeshMaterialVolumes - The mesh material volumes object containing the materials and their - volumes for each mesh element. - - Returns - ------- - list of openmc.MeshSource - A list of MeshSource objects, each containing IndependentSource - instances for the decay photons in the corresponding mesh element. + This class is responsible for managing the materials and sources needed + for mesh-based R2S calculations. It provides methods to get activation + materials and decay photon sources based on the mesh and materials in the + OpenMC model. """ - mat_dict = _get_all_materials(model) - - # Some MeshSource objects will have empty positions; create a "null source" - # that is used for this case - null_source = openmc.IndependentSource(particle='photon', strength=0.0) - - # List to hold sources for each MeshSource (length = N) - source_lists = [] - - # Index in the overall list of activated materials - index_mat = 0 - - # Total number of mesh elements - n_elements = mat_vols.num_elements - - for index_elem in range(n_elements): - for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): - # Skip void volume - if mat_id is None: - continue - - # Check whether a new MeshSource object is needed - if j >= len(source_lists): - source_lists.append([null_source]*n_elements) - - # Get activated material composition - original_mat = materials[index_mat] - activated_mat = results[-1].get_material(str(original_mat.id)) - - # Create decay photon source source - energy = activated_mat.get_decay_photon_energy() - source_lists[j][index_elem] = openmc.IndependentSource( - energy=energy, - particle='photon', - strength=energy.integral(), - constraints={'domains': [mat_dict[mat_id]]} - ) - - # Increment index of activated material - index_mat += 1 - - # Return list of mesh sources - return [openmc.MeshSource(mesh, sources) for sources in source_lists] + def __init__(self, model: openmc.Model, mesh: openmc.MeshBase): + self.model = model + self.mesh = mesh + + def get_activation_materials( + self, mmv: openmc.MeshMaterialVolumes + ) -> list[openmc.Material]: + """Get a list of activation materials for each mesh element/material. + + When performing a mesh-based R2S calculation, a unique material is needed + for each activation region, which is a combination of a mesh element and a + material within that mesh element. This function generates a list of such + materials, each with a unique name and volume corresponding to the mesh + element and material. + + Parameters + ---------- + mmv : openmc.MeshMaterialVolumes + The mesh material volumes object containing the materials and their + volumes for each mesh element. + + Returns + ------- + list of openmc.Material + A list of materials, each corresponding to a unique mesh element and + material combination. + + """ + # Get the material ID, volume, and element index for each element-material + # combination + mat_ids = mmv._materials[mmv._materials > -1] + volumes = mmv._volumes[mmv._materials > -1] + elems, _ = np.where(mmv._materials > -1) + + # Get all materials in the model + material_dict = _get_all_materials(self) + + # Create a new activation material for each element-material combination + materials = [] + for elem, mat_id, vol in zip(elems, mat_ids, volumes): + mat = material_dict[mat_id] + new_mat = mat.clone() + new_mat.depletable = True + new_mat.name = f'Element {elem}, Material {mat_id}' + new_mat.volume = vol + materials.append(new_mat) + + return materials + + def get_decay_photon_source( + self, + mesh: openmc.MeshBase, + materials: list[openmc.Material], + results: Results, + mat_vols: openmc.MeshMaterialVolumes, + ) -> list[openmc.MeshSource]: + """Create decay photon source for a mesh-based R2S calculation. + + This function creates N :class:`MeshSource` objects where N is the maximum + number of unique materials that appears in a single mesh element. For each + mesh element-material combination, and IndependentSource instance is created + with a spatial constraint limited the sampled decay photons to the correct + region. + + Parameters + ---------- + mesh : openmc.MeshBase + The mesh object defining the spatial regions over which activation is + performed. + materials : list of openmc.Material + List of materials that are activated in the model. + results : openmc.deplete.Results + The results object containing the depletion results for the materials. + mat_vols : openmc.MeshMaterialVolumes + The mesh material volumes object containing the materials and their + volumes for each mesh element. + + Returns + ------- + list of openmc.MeshSource + A list of MeshSource objects, each containing IndependentSource + instances for the decay photons in the corresponding mesh element. + + """ + mat_dict = _get_all_materials(self) + + # Some MeshSource objects will have empty positions; create a "null source" + # that is used for this case + null_source = openmc.IndependentSource(particle='photon', strength=0.0) + + # List to hold sources for each MeshSource (length = N) + source_lists = [] + + # Index in the overall list of activated materials + index_mat = 0 + + # Total number of mesh elements + n_elements = mat_vols.num_elements + + for index_elem in range(n_elements): + for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): + # Skip void volume + if mat_id is None: + continue + + # Check whether a new MeshSource object is needed + if j >= len(source_lists): + source_lists.append([null_source]*n_elements) + + # Get activated material composition + original_mat = materials[index_mat] + activated_mat = results[-1].get_material(str(original_mat.id)) + + # Create decay photon source source + energy = activated_mat.get_decay_photon_energy() + source_lists[j][index_elem] = openmc.IndependentSource( + energy=energy, + particle='photon', + strength=energy.integral(), + constraints={'domains': [mat_dict[mat_id]]} + ) + + # Increment index of activated material + index_mat += 1 + + # Return list of mesh sources + return [openmc.MeshSource(mesh, sources) for sources in source_lists] From 3b66be992e782a23c85599bd70298354fcdb4fc9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 11 Jun 2025 06:40:58 -0500 Subject: [PATCH 06/29] Restructure R2SManager class --- openmc/deplete/r2s.py | 295 ++++++++++++++++++++++++------------------ 1 file changed, 171 insertions(+), 124 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 3814aa693bf..10204f18c8a 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -2,12 +2,143 @@ import numpy as np import openmc +from . import get_microxs_and_flux, IndependentOperator, PredictorIntegrator from .results import Results from ..mesh import _get_all_materials +def get_activation_materials( + model: openmc.Model, mmv: openmc.MeshMaterialVolumes +) -> list[openmc.Material]: + """Get a list of activation materials for each mesh element/material. + + When performing a mesh-based R2S calculation, a unique material is needed + for each activation region, which is a combination of a mesh element and a + material within that mesh element. This function generates a list of such + materials, each with a unique name and volume corresponding to the mesh + element and material. + + Parameters + ---------- + model : openmc.Model + The full model containing the geometry and materials. + mmv : openmc.MeshMaterialVolumes + The mesh material volumes object containing the materials and their + volumes for each mesh element. + + Returns + ------- + list of openmc.Material + A list of materials, each corresponding to a unique mesh element and + material combination. + + """ + # Get the material ID, volume, and element index for each element-material + # combination + mat_ids = mmv._materials[mmv._materials > -1] + volumes = mmv._volumes[mmv._materials > -1] + elems, _ = np.where(mmv._materials > -1) + + # Get all materials in the model + material_dict = _get_all_materials(model) + + # Create a new activation material for each element-material combination + materials = [] + for elem, mat_id, vol in zip(elems, mat_ids, volumes): + mat = material_dict[mat_id] + new_mat = mat.clone() + new_mat.depletable = True + new_mat.name = f'Element {elem}, Material {mat_id}' + new_mat.volume = vol + materials.append(new_mat) + + return materials + + +def get_decay_photon_source( + model: openmc.Model, + mesh: openmc.MeshBase, + materials: list[openmc.Material], + results: Results, + mat_vols: openmc.MeshMaterialVolumes, +) -> list[openmc.MeshSource]: + """Create decay photon source for a mesh-based R2S calculation. + + This function creates N :class:`MeshSource` objects where N is the maximum + number of unique materials that appears in a single mesh element. For each + mesh element-material combination, and IndependentSource instance is created + with a spatial constraint limited the sampled decay photons to the correct + region. + + Parameters + ---------- + model : openmc.Model + The full model containing the geometry and materials. + mesh : openmc.MeshBase + The mesh object defining the spatial regions over which activation is + performed. + materials : list of openmc.Material + List of materials that are activated in the model. + results : openmc.deplete.Results + The results object containing the depletion results for the materials. + mat_vols : openmc.MeshMaterialVolumes + The mesh material volumes object containing the materials and their + volumes for each mesh element. + + Returns + ------- + list of openmc.MeshSource + A list of MeshSource objects, each containing IndependentSource + instances for the decay photons in the corresponding mesh element. + + """ + mat_dict = _get_all_materials(model) + + # Some MeshSource objects will have empty positions; create a "null source" + # that is used for this case + null_source = openmc.IndependentSource(particle='photon', strength=0.0) + + # List to hold sources for each MeshSource (length = N) + source_lists = [] + + # Index in the overall list of activated materials + index_mat = 0 + + # Total number of mesh elements + n_elements = mat_vols.num_elements + + for index_elem in range(n_elements): + for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): + # Skip void volume + if mat_id is None: + continue + + # Check whether a new MeshSource object is needed + if j >= len(source_lists): + source_lists.append([null_source]*n_elements) + + # Get activated material composition + original_mat = materials[index_mat] + activated_mat = results[-1].get_material(str(original_mat.id)) + + # Create decay photon source source + energy = activated_mat.get_decay_photon_energy() + source_lists[j][index_elem] = openmc.IndependentSource( + energy=energy, + particle='photon', + strength=energy.integral(), + constraints={'domains': [mat_dict[mat_id]]} + ) + + # Increment index of activated material + index_mat += 1 + + # Return list of mesh sources + return [openmc.MeshSource(mesh, sources) for sources in source_lists] + + class R2SManager: - """Manager for R2S calculations. + """Manager for mesh-based R2S calculations. This class is responsible for managing the materials and sources needed for mesh-based R2S calculations. It provides methods to get activation @@ -19,126 +150,42 @@ def __init__(self, model: openmc.Model, mesh: openmc.MeshBase): self.model = model self.mesh = mesh - def get_activation_materials( - self, mmv: openmc.MeshMaterialVolumes - ) -> list[openmc.Material]: - """Get a list of activation materials for each mesh element/material. - - When performing a mesh-based R2S calculation, a unique material is needed - for each activation region, which is a combination of a mesh element and a - material within that mesh element. This function generates a list of such - materials, each with a unique name and volume corresponding to the mesh - element and material. - - Parameters - ---------- - mmv : openmc.MeshMaterialVolumes - The mesh material volumes object containing the materials and their - volumes for each mesh element. - - Returns - ------- - list of openmc.Material - A list of materials, each corresponding to a unique mesh element and - material combination. - - """ - # Get the material ID, volume, and element index for each element-material - # combination - mat_ids = mmv._materials[mmv._materials > -1] - volumes = mmv._volumes[mmv._materials > -1] - elems, _ = np.where(mmv._materials > -1) - - # Get all materials in the model - material_dict = _get_all_materials(self) - - # Create a new activation material for each element-material combination - materials = [] - for elem, mat_id, vol in zip(elems, mat_ids, volumes): - mat = material_dict[mat_id] - new_mat = mat.clone() - new_mat.depletable = True - new_mat.name = f'Element {elem}, Material {mat_id}' - new_mat.volume = vol - materials.append(new_mat) - - return materials - - def get_decay_photon_source( - self, - mesh: openmc.MeshBase, - materials: list[openmc.Material], - results: Results, - mat_vols: openmc.MeshMaterialVolumes, - ) -> list[openmc.MeshSource]: - """Create decay photon source for a mesh-based R2S calculation. - - This function creates N :class:`MeshSource` objects where N is the maximum - number of unique materials that appears in a single mesh element. For each - mesh element-material combination, and IndependentSource instance is created - with a spatial constraint limited the sampled decay photons to the correct - region. - - Parameters - ---------- - mesh : openmc.MeshBase - The mesh object defining the spatial regions over which activation is - performed. - materials : list of openmc.Material - List of materials that are activated in the model. - results : openmc.deplete.Results - The results object containing the depletion results for the materials. - mat_vols : openmc.MeshMaterialVolumes - The mesh material volumes object containing the materials and their - volumes for each mesh element. - - Returns - ------- - list of openmc.MeshSource - A list of MeshSource objects, each containing IndependentSource - instances for the decay photons in the corresponding mesh element. - - """ - mat_dict = _get_all_materials(self) - - # Some MeshSource objects will have empty positions; create a "null source" - # that is used for this case - null_source = openmc.IndependentSource(particle='photon', strength=0.0) - - # List to hold sources for each MeshSource (length = N) - source_lists = [] - - # Index in the overall list of activated materials - index_mat = 0 - - # Total number of mesh elements - n_elements = mat_vols.num_elements - - for index_elem in range(n_elements): - for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): - # Skip void volume - if mat_id is None: - continue - - # Check whether a new MeshSource object is needed - if j >= len(source_lists): - source_lists.append([null_source]*n_elements) - - # Get activated material composition - original_mat = materials[index_mat] - activated_mat = results[-1].get_material(str(original_mat.id)) - - # Create decay photon source source - energy = activated_mat.get_decay_photon_energy() - source_lists[j][index_elem] = openmc.IndependentSource( - energy=energy, - particle='photon', - strength=energy.integral(), - constraints={'domains': [mat_dict[mat_id]]} - ) - - # Increment index of activated material - index_mat += 1 - - # Return list of mesh sources - return [openmc.MeshSource(mesh, sources) for sources in source_lists] + def run(self): + self.step1_neutron_transport() + self.step2_activation() + self.step3_decay_photon_source() + + def step1_neutron_transport(self, mat_vol_kwargs, micro_kwargs): + # Compute material volume fractions on the mesh + self.mmv = self.mesh.material_volumes(**mat_vol_kwargs) + + # Create mesh-material filter based on what combos were found + mm_filter = openmc.MeshMaterialFilter.from_volumes(self.mesh, self.mmv) + + # Get fluxes and microscopic cross sections on each mesh/material + self.fluxes, self.micros = get_microxs_and_flux( + self.model, mm_filter, **micro_kwargs) + + def step2_activation(self, timesteps, source_rates): + # Get unique material for each (mesh, material) combination + self.activation_materials = get_activation_materials(self.model, self.mmv) + + # Create depletion operator for the activation materials + op = IndependentOperator( + self.activation_materials, self.fluxes, self.micros, + normalization_mode='source-rate' + ) + + # Create time integrator and solve depletion equations + integrator = PredictorIntegrator(op, timesteps, source_rates=source_rates) + integrator.integrate(final_step=False) + + # Get depletion results + self.results = Results() + + def step3_decay_photon_source(self): + # Create decay photon source + self.model.settings.source = get_decay_photon_source( + self.model, self.mesh, self.activation_materials, self.results, self.mmv + ) + From a1f5e02255c7a7af37f59f5e14a2376f3f258b3c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 12 Jun 2025 23:57:27 -0500 Subject: [PATCH 07/29] Add some TODOs --- openmc/deplete/r2s.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 10204f18c8a..4bc3baa654b 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -149,13 +149,21 @@ class R2SManager: def __init__(self, model: openmc.Model, mesh: openmc.MeshBase): self.model = model self.mesh = mesh + # TODO: Think about directory structure and file naming + # TODO: Option to select cell-based or mesh-basedk + # TODO: Photon settings + # TODO: Think about MPI + # TODO: Voiding/changing materials between neutron/photon steps def run(self): self.step1_neutron_transport() self.step2_activation() self.step3_decay_photon_source() - def step1_neutron_transport(self, mat_vol_kwargs, micro_kwargs): + def step1_neutron_transport(self, output_dir="neutron_transport", mat_vol_kwargs=None, micro_kwargs=None): + # TODO: Run this in a "neutron_transport" subdirectory? + # TODO: Energy discretization + # Compute material volume fractions on the mesh self.mmv = self.mesh.material_volumes(**mat_vol_kwargs) From 15ea82db7054a59022b85475425818ceed8963c2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 22 Jun 2025 18:30:47 -0500 Subject: [PATCH 08/29] Add option for cell-based or mesh-based --- openmc/deplete/r2s.py | 148 +++++++++++++++++++++++++++++++++--------- 1 file changed, 118 insertions(+), 30 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 4bc3baa654b..55bd1b12e18 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -1,4 +1,5 @@ from __future__ import annotations +from collections.abc import Sequence import numpy as np import openmc @@ -55,12 +56,13 @@ def get_activation_materials( return materials -def get_decay_photon_source( +def get_decay_photon_source_mesh( model: openmc.Model, mesh: openmc.MeshBase, materials: list[openmc.Material], results: Results, mat_vols: openmc.MeshMaterialVolumes, + time_index: int = -1, ) -> list[openmc.MeshSource]: """Create decay photon source for a mesh-based R2S calculation. @@ -119,7 +121,7 @@ def get_decay_photon_source( # Get activated material composition original_mat = materials[index_mat] - activated_mat = results[-1].get_material(str(original_mat.id)) + activated_mat = results[time_index].get_material(str(original_mat.id)) # Create decay photon source source energy = activated_mat.get_decay_photon_energy() @@ -140,60 +142,146 @@ def get_decay_photon_source( class R2SManager: """Manager for mesh-based R2S calculations. - This class is responsible for managing the materials and sources needed - for mesh-based R2S calculations. It provides methods to get activation - materials and decay photon sources based on the mesh and materials in the - OpenMC model. + This class is responsible for managing the materials and sources needed for + mesh-based or cell-based R2S calculations. It provides methods to get + activation materials and decay photon sources based on the mesh/cells and + materials in the OpenMC model. """ - def __init__(self, model: openmc.Model, mesh: openmc.MeshBase): + def __init__( + self, + model: openmc.Model, + domains: openmc.MeshBase, + ): self.model = model - self.mesh = mesh + if isinstance(domains, openmc.Mesh): + self.method = 'mesh-based' + else: + self.method = 'cell-based' + self.domains = domains # TODO: Think about directory structure and file naming # TODO: Option to select cell-based or mesh-basedk # TODO: Photon settings # TODO: Think about MPI + # TODO: Option to use CoupledOperator? Needed for true burnup # TODO: Voiding/changing materials between neutron/photon steps - - def run(self): - self.step1_neutron_transport() - self.step2_activation() + self.results = {} + + def run( + self, + timesteps: Sequence[float] | Sequence[tuple[float, str]], + source_rates: float | Sequence[float], + timestep_units: str = 's', + micro_kwargs: dict | None = None, + mat_vol_kwargs: dict | None = None, + ): + self.step1_neutron_transport(micro_kwargs=micro_kwargs, mat_vol_kwargs=mat_vol_kwargs) + self.step2_activation(timesteps, source_rates, timestep_units) self.step3_decay_photon_source() def step1_neutron_transport(self, output_dir="neutron_transport", mat_vol_kwargs=None, micro_kwargs=None): + """Run the neutron transport step. + + This step computes the material volume fractions on the mesh, creates a + mesh-material filter, and retrieves the fluxes and microscopic cross + sections for each mesh/material combination. It populates the results + + This step will populate the 'mesh_material_volumes', 'fluxes', and + 'micros' keys in the results dictionary. + """ + # TODO: Run this in a "neutron_transport" subdirectory? # TODO: Energy discretization - # Compute material volume fractions on the mesh - self.mmv = self.mesh.material_volumes(**mat_vol_kwargs) + if self.method == 'mesh-based': + # Compute material volume fractions on the mesh + self.results['mesh_material_volumes'] = mmv = \ + self.mesh.material_volumes(**mat_vol_kwargs) - # Create mesh-material filter based on what combos were found - mm_filter = openmc.MeshMaterialFilter.from_volumes(self.mesh, self.mmv) + # Create mesh-material filter based on what combos were found + domains = openmc.MeshMaterialFilter.from_volumes(self.domains, mmv) + else: + # TODO: Run volume calculation for cells + domains = self.domains # Get fluxes and microscopic cross sections on each mesh/material - self.fluxes, self.micros = get_microxs_and_flux( - self.model, mm_filter, **micro_kwargs) - - def step2_activation(self, timesteps, source_rates): - # Get unique material for each (mesh, material) combination - self.activation_materials = get_activation_materials(self.model, self.mmv) + self.results['fluxes'], self.results['micros'] = get_microxs_and_flux( + self.model, domains, **micro_kwargs) + + def step2_activation( + self, + timesteps: Sequence[float] | Sequence[tuple[float, str]], + source_rates: float | Sequence[float], + timestep_units: str = 's' + ): + if self.method == 'mesh-based': + # Get unique material for each (mesh, material) combination + mmv = self.results['mesh_material_volumes'] + self.results['activation_materials'] = get_activation_materials(self.model, mmv) + else: + # Create unique material for each cell + activation_mats = [] + for cell in self.domains: + mat = cell.fill.clone() + mat.name = f'Cell {cell.id}' + mat.depletable = True + mat.volume = cell.volume + activation_mats.append(mat) + self.results['activation_materials'] = activation_mats # Create depletion operator for the activation materials op = IndependentOperator( - self.activation_materials, self.fluxes, self.micros, + self.results['activation_materials'], + self.results['fluxes'], + self.results['micros'], normalization_mode='source-rate' ) # Create time integrator and solve depletion equations - integrator = PredictorIntegrator(op, timesteps, source_rates=source_rates) + integrator = PredictorIntegrator( + op, timesteps, source_rates=source_rates, + timestep_units=timestep_units + ) integrator.integrate(final_step=False) # Get depletion results - self.results = Results() + self.results['depletion_results'] = Results() - def step3_decay_photon_source(self): - # Create decay photon source - self.model.settings.source = get_decay_photon_source( - self.model, self.mesh, self.activation_materials, self.results, self.mmv - ) + def step3_decay_photon_source( + self, + bounding_boxes: dict[int, openmc.BoundingBox] | None = None + ): + # TODO: Automatically determine bounding box for each cell + # Create decay photon source + # TODO: Add loop over cooling times + if self.method == 'mesh-based': + self.model.settings.source = get_decay_photon_source_mesh( + self.model, + self.mesh, + self.results['activation_materials'], + self.results['depletion_results'], + self.results['mesh_material_volumes'], + time_index=-1 + ) + else: + sources = [] + results = self.results['depletion_results'] + for cell, original_mat in zip(self.domains, self.results['activation_materials']): + bounding_box = bounding_boxes[cell.id] + + # Get activated material composition + activated_mat = results[-1].get_material(str(original_mat.id)) + + # Create decay photon source source + space = openmc.stats.Box(*bounding_box) + energy = activated_mat.get_decay_photon_energy() + source = openmc.IndependentSource( + space=space, + energy=energy, + particle='photon', + strength=energy.integral(), + domains=[cell] + ) + sources.append(source) + self.model.settings.source = sources From 6c0c1f8a56acdc0fa0bead655cc5dd5e0b723893 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 22 Jun 2025 19:07:49 -0500 Subject: [PATCH 09/29] Add specification of cooling times for photon transport --- openmc/deplete/r2s.py | 89 ++++++++++++++++++++++++++----------------- 1 file changed, 54 insertions(+), 35 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 55bd1b12e18..7ff09009ed4 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -160,7 +160,6 @@ def __init__( self.method = 'cell-based' self.domains = domains # TODO: Think about directory structure and file naming - # TODO: Option to select cell-based or mesh-basedk # TODO: Photon settings # TODO: Think about MPI # TODO: Option to use CoupledOperator? Needed for true burnup @@ -171,13 +170,17 @@ def run( self, timesteps: Sequence[float] | Sequence[tuple[float, str]], source_rates: float | Sequence[float], + cooling_times: Sequence[int], + dose_tallies: Sequence[openmc.Tally] | None = None, timestep_units: str = 's', + bounding_boxes: dict[int, openmc.BoundingBox] | None = None, micro_kwargs: dict | None = None, mat_vol_kwargs: dict | None = None, + run_kwargs: dict | None = None, ): self.step1_neutron_transport(micro_kwargs=micro_kwargs, mat_vol_kwargs=mat_vol_kwargs) self.step2_activation(timesteps, source_rates, timestep_units) - self.step3_decay_photon_source() + self.step3_photon_transport(cooling_times, dose_tallies, bounding_boxes, run_kwargs) def step1_neutron_transport(self, output_dir="neutron_transport", mat_vol_kwargs=None, micro_kwargs=None): """Run the neutron transport step. @@ -247,41 +250,57 @@ def step2_activation( # Get depletion results self.results['depletion_results'] = Results() - def step3_decay_photon_source( + def step3_photon_transport( self, - bounding_boxes: dict[int, openmc.BoundingBox] | None = None + cooling_times: Sequence[int], + dose_tallies: Sequence[openmc.Tally] | None = None, + bounding_boxes: dict[int, openmc.BoundingBox] | None = None, + run_kwargs: dict | None = None, ): # TODO: Automatically determine bounding box for each cell - # Create decay photon source - # TODO: Add loop over cooling times - if self.method == 'mesh-based': - self.model.settings.source = get_decay_photon_source_mesh( - self.model, - self.mesh, - self.results['activation_materials'], - self.results['depletion_results'], - self.results['mesh_material_volumes'], - time_index=-1 - ) - else: - sources = [] - results = self.results['depletion_results'] - for cell, original_mat in zip(self.domains, self.results['activation_materials']): - bounding_box = bounding_boxes[cell.id] - - # Get activated material composition - activated_mat = results[-1].get_material(str(original_mat.id)) - - # Create decay photon source source - space = openmc.stats.Box(*bounding_box) - energy = activated_mat.get_decay_photon_energy() - source = openmc.IndependentSource( - space=space, - energy=energy, - particle='photon', - strength=energy.integral(), - domains=[cell] + # Set default run arguments if not provided + if run_kwargs is None: + run_kwargs = {} + run_kwargs.setdefault('output', False) + + # Add dose tallies to model if provided + if dose_tallies is not None: + self.model.tallies.extend(dose_tallies) + + for time_index in cooling_times: + # Create decay photon source + if self.method == 'mesh-based': + self.model.settings.source = get_decay_photon_source_mesh( + self.model, + self.mesh, + self.results['activation_materials'], + self.results['depletion_results'], + self.results['mesh_material_volumes'], + time_index=time_index, ) - sources.append(source) - self.model.settings.source = sources + else: + sources = [] + results = self.results['depletion_results'] + for cell, original_mat in zip(self.domains, self.results['activation_materials']): + bounding_box = bounding_boxes[cell.id] + + # Get activated material composition + activated_mat = results[time_index].get_material(str(original_mat.id)) + + # Create decay photon source source + space = openmc.stats.Box(*bounding_box) + energy = activated_mat.get_decay_photon_energy() + source = openmc.IndependentSource( + space=space, + energy=energy, + particle='photon', + strength=energy.integral(), + domains=[cell] + ) + sources.append(source) + self.model.settings.source = sources + + # Run photon transport calculation + run_kwargs['cwd'] = f'photon_transport_{time_index}' + self.model.run(**run_kwargs) From d345386178c3624877b2de60eb4df1dfb4ce93ec Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 29 Jun 2025 21:42:45 -0500 Subject: [PATCH 10/29] Add output directories --- openmc/deplete/__init__.py | 1 + openmc/deplete/r2s.py | 55 +++++++++++++++++++++++++++++--------- 2 files changed, 44 insertions(+), 12 deletions(-) diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index 8a9509e9009..052e2245963 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -17,6 +17,7 @@ from .results import * from .integrators import * from .transfer_rates import * +from .r2s import * from . import abc from . import cram from . import helpers diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 7ff09009ed4..e9d220ccc37 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -1,10 +1,12 @@ from __future__ import annotations from collections.abc import Sequence +from pathlib import Path import numpy as np import openmc from . import get_microxs_and_flux, IndependentOperator, PredictorIntegrator from .results import Results +from ..checkvalue import PathLike from ..mesh import _get_all_materials @@ -154,7 +156,7 @@ def __init__( domains: openmc.MeshBase, ): self.model = model - if isinstance(domains, openmc.Mesh): + if isinstance(domains, openmc.MeshBase): self.method = 'mesh-based' else: self.method = 'cell-based' @@ -180,26 +182,36 @@ def run( ): self.step1_neutron_transport(micro_kwargs=micro_kwargs, mat_vol_kwargs=mat_vol_kwargs) self.step2_activation(timesteps, source_rates, timestep_units) - self.step3_photon_transport(cooling_times, dose_tallies, bounding_boxes, run_kwargs) + self.step3_photon_transport(cooling_times, dose_tallies, bounding_boxes, run_kwargs=run_kwargs) - def step1_neutron_transport(self, output_dir="neutron_transport", mat_vol_kwargs=None, micro_kwargs=None): + def step1_neutron_transport( + self, + output_dir: PathLike = "neutron_transport", + mat_vol_kwargs: dict | None = None, + micro_kwargs: dict | None = None + ): """Run the neutron transport step. This step computes the material volume fractions on the mesh, creates a mesh-material filter, and retrieves the fluxes and microscopic cross sections for each mesh/material combination. It populates the results - This step will populate the 'mesh_material_volumes', 'fluxes', and - 'micros' keys in the results dictionary. + This step will populate the 'fluxes' and 'micros' keys in the results + dictionary. For a mesh-based calculation, it will also populate the + 'mesh_material_volumes' key. """ - # TODO: Run this in a "neutron_transport" subdirectory? # TODO: Energy discretization + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) if self.method == 'mesh-based': # Compute material volume fractions on the mesh self.results['mesh_material_volumes'] = mmv = \ - self.mesh.material_volumes(**mat_vol_kwargs) + self.domains.material_volumes(self.model, **mat_vol_kwargs) + + # Save results to file + mmv.save(output_dir / 'mesh_material_volumes.h5') # Create mesh-material filter based on what combos were found domains = openmc.MeshMaterialFilter.from_volumes(self.domains, mmv) @@ -208,14 +220,25 @@ def step1_neutron_transport(self, output_dir="neutron_transport", mat_vol_kwargs domains = self.domains # Get fluxes and microscopic cross sections on each mesh/material + if micro_kwargs is None: + micro_kwargs = {} + micro_kwargs.setdefault('path_statepoint', output_dir / 'statepoint.h5') + + # Run neutron transport and get fluxes and micros self.results['fluxes'], self.results['micros'] = get_microxs_and_flux( self.model, domains, **micro_kwargs) + # Save flux and micros to file + np.save(output_dir / 'fluxes.npy', self.results['fluxes']) + for i, micros in enumerate(self.results['micros']): + micros.to_csv(output_dir / f'micros_{i}.csv') + def step2_activation( self, timesteps: Sequence[float] | Sequence[tuple[float, str]], source_rates: float | Sequence[float], - timestep_units: str = 's' + timestep_units: str = 's', + output_dir: PathLike = 'activation', ): if self.method == 'mesh-based': # Get unique material for each (mesh, material) combination @@ -245,16 +268,20 @@ def step2_activation( op, timesteps, source_rates=source_rates, timestep_units=timestep_units ) - integrator.integrate(final_step=False) + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + output_path = output_dir / 'depletion_results.h5' + integrator.integrate(final_step=False, path=output_path) # Get depletion results - self.results['depletion_results'] = Results() + self.results['depletion_results'] = Results(output_path) def step3_photon_transport( self, cooling_times: Sequence[int], dose_tallies: Sequence[openmc.Tally] | None = None, bounding_boxes: dict[int, openmc.BoundingBox] | None = None, + output_dir: PathLike = 'photon_transport', run_kwargs: dict | None = None, ): # TODO: Automatically determine bounding box for each cell @@ -273,7 +300,7 @@ def step3_photon_transport( if self.method == 'mesh-based': self.model.settings.source = get_decay_photon_source_mesh( self.model, - self.mesh, + self.domains, self.results['activation_materials'], self.results['depletion_results'], self.results['mesh_material_volumes'], @@ -301,6 +328,10 @@ def step3_photon_transport( sources.append(source) self.model.settings.source = sources + # Convert time_index (which may be negative) to a normal index + if time_index < 0: + time_index = len(self.results['depletion_results']) + time_index + # Run photon transport calculation - run_kwargs['cwd'] = f'photon_transport_{time_index}' + run_kwargs['cwd'] = Path(output_dir) / f'time_{time_index}' self.model.run(**run_kwargs) From 51c43308633b8fe8cb20cdf25ce068e0cc7950e1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 4 Jul 2025 00:10:34 -0500 Subject: [PATCH 11/29] Store all results under a single r2s_timestamp dir --- openmc/deplete/r2s.py | 36 +++++++++++++++++++++++++++--------- 1 file changed, 27 insertions(+), 9 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index e9d220ccc37..78494ac27d9 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -1,5 +1,6 @@ from __future__ import annotations from collections.abc import Sequence +from datetime import datetime from pathlib import Path import numpy as np @@ -142,7 +143,7 @@ def get_decay_photon_source_mesh( class R2SManager: - """Manager for mesh-based R2S calculations. + """Manager for Rigorous 2-Step (R2S) method calculations. This class is responsible for managing the materials and sources needed for mesh-based or cell-based R2S calculations. It provides methods to get @@ -161,7 +162,6 @@ def __init__( else: self.method = 'cell-based' self.domains = domains - # TODO: Think about directory structure and file naming # TODO: Photon settings # TODO: Think about MPI # TODO: Option to use CoupledOperator? Needed for true burnup @@ -175,14 +175,30 @@ def run( cooling_times: Sequence[int], dose_tallies: Sequence[openmc.Tally] | None = None, timestep_units: str = 's', + output_dir: PathLike | None = None, bounding_boxes: dict[int, openmc.BoundingBox] | None = None, micro_kwargs: dict | None = None, mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, ): - self.step1_neutron_transport(micro_kwargs=micro_kwargs, mat_vol_kwargs=mat_vol_kwargs) - self.step2_activation(timesteps, source_rates, timestep_units) - self.step3_photon_transport(cooling_times, dose_tallies, bounding_boxes, run_kwargs=run_kwargs) + """Run the R2S calculation.""" + + if output_dir is None: + stamp = datetime.now().strftime('%Y-%m-%dT%H-%M-%S') + output_dir = Path(f'r2s_{stamp}') + + self.step1_neutron_transport( + output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs + ) + self.step2_activation( + timesteps, source_rates, timestep_units, output_dir / 'activation' + ) + self.step3_photon_transport( + cooling_times, dose_tallies, bounding_boxes, + output_dir / 'photon_transport', run_kwargs=run_kwargs + ) + + return output_dir def step1_neutron_transport( self, @@ -194,11 +210,10 @@ def step1_neutron_transport( This step computes the material volume fractions on the mesh, creates a mesh-material filter, and retrieves the fluxes and microscopic cross - sections for each mesh/material combination. It populates the results + sections for each mesh/material combination. This step will populate the + 'fluxes' and 'micros' keys in the results dictionary. For a mesh-based + calculation, it will also populate the 'mesh_material_volumes' key. - This step will populate the 'fluxes' and 'micros' keys in the results - dictionary. For a mesh-based calculation, it will also populate the - 'mesh_material_volumes' key. """ # TODO: Energy discretization @@ -228,6 +243,9 @@ def step1_neutron_transport( self.results['fluxes'], self.results['micros'] = get_microxs_and_flux( self.model, domains, **micro_kwargs) + # Export model to output directory + self.model.export_to_model_xml(output_dir / 'model.xml') + # Save flux and micros to file np.save(output_dir / 'fluxes.npy', self.results['fluxes']) for i, micros in enumerate(self.results['micros']): From c16e60c1b25db5283acd73d15215ebd682512ba4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 13 Jul 2025 17:42:24 -0500 Subject: [PATCH 12/29] Add docstrings for methods in R2SManager --- openmc/deplete/r2s.py | 137 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 135 insertions(+), 2 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 78494ac27d9..4d572bb37ca 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -150,11 +150,32 @@ class R2SManager: activation materials and decay photon sources based on the mesh/cells and materials in the OpenMC model. + Parameters + ---------- + model : openmc.Model + The OpenMC model containing the geometry and materials. + domains : openmc.MeshBase or Sequence[openmc.Cell] + The mesh or a sequence of cells that represent the spatial units over + which the R2S calculation will be performed. + + Attributes + ---------- + domains : openmc.MeshBase or Sequence[openmc.Cell] + The mesh or a sequence of cells that represent the spatial units over + which the R2S calculation will be performed. + model : openmc.Model + The OpenMC model containing the geometry and materials. + method : {'mesh-based', 'cell-based'} + Indicates whether the R2S calculation uses mesh elements ('mesh-based') + as the spatial discetization or a list of a cells ('cell-based'). + results : dict + A dictionary that stores results from the R2S calculation. + """ def __init__( self, model: openmc.Model, - domains: openmc.MeshBase, + domains: openmc.MeshBase | Sequence[openmc.Cell], ): self.model = model if isinstance(domains, openmc.MeshBase): @@ -181,7 +202,54 @@ def run( mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, ): - """Run the R2S calculation.""" + """Run the R2S calculation. + + Parameters + ---------- + timesteps : Sequence[float] or Sequence[tuple[float, str]] + Sequence of timesteps. Note that values are not cumulative. The + units are specified by the `timestep_units` argument when + `timesteps` is an iterable of float. Alternatively, units can be + specified for each step by passing an iterable of (value, unit) + tuples. + source_rates : float or Sequence[float] + Source rate in [neutron/sec] for each interval in `timesteps`. + cooling_times : Sequence[int] + Sequence of cooling times represented as indices into the array of + `timesteps`. + dose_tallies : Sequence[openmc.Tally], optional + A sequence of tallies to be used for dose calculations in the photon + transport step. If None, no dose tallies are added. + timestep_units : {'s', 'min', 'h', 'd', 'a'}, optional + Units for values specified in the `timesteps` argument when passing + float values. 's' means seconds, 'min' means minutes, 'h' means + hours, 'd' means days, and 'a' means years (Julian). + output_dir : PathLike, optional + Path to directory where R2S calculation outputs will be saved. If + not provided, a timestamped directory 'r2s_YYYY-MM-DDTHH-MM-SS' is + created. Subdirectories will be created for the neutron transport, + activation, and photon transport steps. + bounding_boxes : dict[int, openmc.BoundingBox], optional + Dictionary mapping cell IDs to bounding boxes used for spatial + source sampling in cell-based R2S calculations. Required if method + is 'cell-based'. + micro_kwargs : dict, optional + Additional keyword arguments passed to + :func:`openmc.deplete.get_microxs_and_flux` during the neutron + transport step. + mat_vol_kwargs : dict, optional + Additional keyword arguments passed to + :meth:`openmc.MeshBase.material_volumes` during the neutron + transport step. + run_kwargs : dict, optional + Additional keyword arguments passed to :meth:`openmc.Model.run` + during the photon transport step. By default, output is disabled. + + Returns + ------- + Path + Path to the output directory containing all calculation results + """ if output_dir is None: stamp = datetime.now().strftime('%Y-%m-%dT%H-%M-%S') @@ -214,6 +282,17 @@ def step1_neutron_transport( 'fluxes' and 'micros' keys in the results dictionary. For a mesh-based calculation, it will also populate the 'mesh_material_volumes' key. + Parameters + ---------- + output_dir : PathLike, optional + The directory where the results will be saved. + mat_vol_kwargs : dict, optional + Additional keyword arguments based to + :meth:`openmc.MeshBase.material_volumes`. + micro_kwargs : dict, optional + Additional keyword arguments passed to + :func:`openmc.deplete.get_microxs_and_flux`. + """ # TODO: Energy discretization @@ -258,6 +337,32 @@ def step2_activation( timestep_units: str = 's', output_dir: PathLike = 'activation', ): + """Run the activation step. + + This step creates a unique copy of each activation material based on the + mesh elements or cells, then solves the depletion equations for each + material using the fluxes and microscopic cross sections obtained in the + neutron transport step. + + Parameters + ---------- + timesteps : Sequence[float] or Sequence[tuple[float, str]] + Sequence of timesteps. Note that values are not cumulative. The + units are specified by the `timestep_units` argument when + `timesteps` is an iterable of float. Alternatively, units can be + specified for each step by passing an iterable of (value, unit) + tuples. + source_rates : float | Sequence[float] + Source rate in [neutron/sec] for each interval in `timesteps`. + timestep_units : {'s', 'min', 'h', 'd', 'a'}, optional + Units for values specified in the `timesteps` argument when passing + float values. 's' means seconds, 'min' means minutes, 'h' means + hours, 'd' means days, and 'a' means years (Julian). + output_dir : PathLike, optional + Path to directory where activation calculation outputs will be + saved. + """ + if self.method == 'mesh-based': # Get unique material for each (mesh, material) combination mmv = self.results['mesh_material_volumes'] @@ -302,6 +407,34 @@ def step3_photon_transport( output_dir: PathLike = 'photon_transport', run_kwargs: dict | None = None, ): + """Run the photon transport step. + + This step performs photon transport calculations using decay photon + sources created from the activated materials. For each cooling time, + it creates appropriate photon sources and runs a transport calculation. + In mesh-based mode, the sources are created using the mesh material + volumes, while in cell-based mode, they are created using bounding + boxes for each cell. + + Parameters + ---------- + cooling_times : Sequence[int] + Sequence of cooling times represented as indices into the array of + timesteps used in the activation step. + dose_tallies : Sequence[openmc.Tally], optional + A sequence of tallies to be used for dose calculations in the photon + transport step. If None, no dose tallies are added. + bounding_boxes : dict[int, openmc.BoundingBox], optional + Dictionary mapping cell IDs to bounding boxes used for spatial + source sampling in cell-based R2S calculations. Required if method + is 'cell-based'. + output_dir : PathLike, optional + Path to directory where photon transport outputs will be saved. + run_kwargs : dict, optional + Additional keyword arguments passed to :meth:`openmc.Model.run` + during the photon transport step. By default, output is disabled. + """ + # TODO: Automatically determine bounding box for each cell # Set default run arguments if not provided From 2d3f52c4e8731af949de259408ed7a8766b81b13 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sun, 13 Jul 2025 17:54:50 -0500 Subject: [PATCH 13/29] Option to use separate photon settings --- openmc/deplete/r2s.py | 111 +++++++++++++++++++++++++----------------- 1 file changed, 65 insertions(+), 46 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 4d572bb37ca..5cab5c9058d 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -183,7 +183,6 @@ def __init__( else: self.method = 'cell-based' self.domains = domains - # TODO: Photon settings # TODO: Think about MPI # TODO: Option to use CoupledOperator? Needed for true burnup # TODO: Voiding/changing materials between neutron/photon steps @@ -201,6 +200,7 @@ def run( micro_kwargs: dict | None = None, mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, + photon_settings: openmc.Settings | None = None, ): """Run the R2S calculation. @@ -244,6 +244,10 @@ def run( run_kwargs : dict, optional Additional keyword arguments passed to :meth:`openmc.Model.run` during the photon transport step. By default, output is disabled. + photon_settings : openmc.Settings, optional + Custom settings to use for the photon transport calculation. If + provided, the original model settings will be temporarily replaced + during the photon transport calculations and restored afterward. Returns ------- @@ -263,7 +267,8 @@ def run( ) self.step3_photon_transport( cooling_times, dose_tallies, bounding_boxes, - output_dir / 'photon_transport', run_kwargs=run_kwargs + output_dir / 'photon_transport', run_kwargs=run_kwargs, + settings=photon_settings ) return output_dir @@ -406,15 +411,16 @@ def step3_photon_transport( bounding_boxes: dict[int, openmc.BoundingBox] | None = None, output_dir: PathLike = 'photon_transport', run_kwargs: dict | None = None, + settings: openmc.Settings | None = None, ): """Run the photon transport step. This step performs photon transport calculations using decay photon - sources created from the activated materials. For each cooling time, - it creates appropriate photon sources and runs a transport calculation. - In mesh-based mode, the sources are created using the mesh material - volumes, while in cell-based mode, they are created using bounding - boxes for each cell. + sources created from the activated materials. For each cooling time, it + creates appropriate photon sources and runs a transport calculation. In + mesh-based mode, the sources are created using the mesh material + volumes, while in cell-based mode, they are created using bounding boxes + for each cell. Parameters ---------- @@ -433,6 +439,10 @@ def step3_photon_transport( run_kwargs : dict, optional Additional keyword arguments passed to :meth:`openmc.Model.run` during the photon transport step. By default, output is disabled. + settings : openmc.Settings, optional + Custom settings to use for the photon transport calculation. If + provided, the original model settings will be temporarily replaced + during the photon transport calculations and restored afterward. """ # TODO: Automatically determine bounding box for each cell @@ -446,43 +456,52 @@ def step3_photon_transport( if dose_tallies is not None: self.model.tallies.extend(dose_tallies) - for time_index in cooling_times: - # Create decay photon source - if self.method == 'mesh-based': - self.model.settings.source = get_decay_photon_source_mesh( - self.model, - self.domains, - self.results['activation_materials'], - self.results['depletion_results'], - self.results['mesh_material_volumes'], - time_index=time_index, - ) - else: - sources = [] - results = self.results['depletion_results'] - for cell, original_mat in zip(self.domains, self.results['activation_materials']): - bounding_box = bounding_boxes[cell.id] - - # Get activated material composition - activated_mat = results[time_index].get_material(str(original_mat.id)) - - # Create decay photon source source - space = openmc.stats.Box(*bounding_box) - energy = activated_mat.get_decay_photon_energy() - source = openmc.IndependentSource( - space=space, - energy=energy, - particle='photon', - strength=energy.integral(), - domains=[cell] + # Save original settings to restore later + original_settings = self.model.settings + if settings is not None: + self.model.settings = settings + + try: + for time_index in cooling_times: + # Create decay photon source + if self.method == 'mesh-based': + self.model.settings.source = get_decay_photon_source_mesh( + self.model, + self.domains, + self.results['activation_materials'], + self.results['depletion_results'], + self.results['mesh_material_volumes'], + time_index=time_index, ) - sources.append(source) - self.model.settings.source = sources - - # Convert time_index (which may be negative) to a normal index - if time_index < 0: - time_index = len(self.results['depletion_results']) + time_index - - # Run photon transport calculation - run_kwargs['cwd'] = Path(output_dir) / f'time_{time_index}' - self.model.run(**run_kwargs) + else: + sources = [] + results = self.results['depletion_results'] + for cell, original_mat in zip(self.domains, self.results['activation_materials']): + bounding_box = bounding_boxes[cell.id] + + # Get activated material composition + activated_mat = results[time_index].get_material(str(original_mat.id)) + + # Create decay photon source source + space = openmc.stats.Box(*bounding_box) + energy = activated_mat.get_decay_photon_energy() + source = openmc.IndependentSource( + space=space, + energy=energy, + particle='photon', + strength=energy.integral(), + domains=[cell] + ) + sources.append(source) + self.model.settings.source = sources + + # Convert time_index (which may be negative) to a normal index + if time_index < 0: + time_index = len(self.results['depletion_results']) + time_index + + # Run photon transport calculation + run_kwargs['cwd'] = Path(output_dir) / f'time_{time_index}' + self.model.run(**run_kwargs) + finally: + # Restore original settings + self.model.settings = original_settings From 947b66ccef70daa9a0bdf61e6a2a4e3dc09c1404 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 14 Jul 2025 10:54:16 -0500 Subject: [PATCH 14/29] Add dose_tallies results, better argument handling --- openmc/deplete/microxs.py | 13 ++++++-- openmc/deplete/r2s.py | 62 +++++++++++++++++++++++++++++---------- openmc/mesh.py | 2 +- 3 files changed, 59 insertions(+), 18 deletions(-) diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index df78f8a267c..50235e58df3 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -47,6 +47,7 @@ def get_microxs_and_flux( reaction_rate_mode: str = 'direct', chain_file: PathLike | Chain | None = None, path_statepoint: PathLike | None = None, + path_input: PathLike | None = None, run_kwargs=None ) -> tuple[list[np.ndarray], list[MicroXS]]: """Generate microscopic cross sections and fluxes for multiple domains. @@ -59,7 +60,7 @@ def get_microxs_and_flux( .. versionadded:: 0.14.0 .. versionchanged:: 0.15.3 - Added `reaction_rate_mode` and `path_statepoint` arguments. + Added `reaction_rate_mode`, `path_statepoint`, `path_input` arguments. Parameters ---------- @@ -90,6 +91,10 @@ def get_microxs_and_flux( Path to write the statepoint file from the neutron transport solve to. By default, The statepoint file is written to a temporary directory and is not kept. + path_input : path-like, optional + Path to write the model XML file from the neutron transport solve to. + By default, the model XML file is written to a temporary directory and + not kept. run_kwargs : dict, optional Keyword arguments passed to :meth:`openmc.Model.run` @@ -108,7 +113,7 @@ def get_microxs_and_flux( check_value('reaction_rate_mode', reaction_rate_mode, {'direct', 'flux'}) # Save any original tallies on the model - original_tallies = model.tallies + original_tallies = list(model.tallies) # Determine what reactions and nuclides are available in chain chain = _get_chain(chain_file) @@ -178,6 +183,10 @@ def get_microxs_and_flux( shutil.move(statepoint_path, path_statepoint) statepoint_path = path_statepoint + # Export the model to path_input if provided + if path_input is not None: + model.export_to_model_xml(path_input) + with StatePoint(statepoint_path) as sp: if reaction_rate_mode == 'direct': rr_tally = sp.tallies[rr_tally.id] diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 5cab5c9058d..bbb76abf2de 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -183,9 +183,6 @@ def __init__( else: self.method = 'cell-based' self.domains = domains - # TODO: Think about MPI - # TODO: Option to use CoupledOperator? Needed for true burnup - # TODO: Voiding/changing materials between neutron/photon steps self.results = {} def run( @@ -243,7 +240,8 @@ def run( transport step. run_kwargs : dict, optional Additional keyword arguments passed to :meth:`openmc.Model.run` - during the photon transport step. By default, output is disabled. + during the neutron and photon transport step. By default, output is + disabled. photon_settings : openmc.Settings, optional Custom settings to use for the photon transport calculation. If provided, the original model settings will be temporarily replaced @@ -259,6 +257,14 @@ def run( stamp = datetime.now().strftime('%Y-%m-%dT%H-%M-%S') output_dir = Path(f'r2s_{stamp}') + # Set run_kwargs for the neutron transport step + if micro_kwargs is None: + micro_kwargs = {} + if run_kwargs is None: + run_kwargs = {} + run_kwargs.setdefault('output', False) + micro_kwargs.setdefault('run_kwargs', run_kwargs) + self.step1_neutron_transport( output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs ) @@ -300,7 +306,6 @@ def step1_neutron_transport( """ - # TODO: Energy discretization output_dir = Path(output_dir) output_dir.mkdir(parents=True, exist_ok=True) @@ -310,26 +315,41 @@ def step1_neutron_transport( self.domains.material_volumes(self.model, **mat_vol_kwargs) # Save results to file - mmv.save(output_dir / 'mesh_material_volumes.h5') + mmv.save(output_dir / 'mesh_material_volumes.npz') # Create mesh-material filter based on what combos were found domains = openmc.MeshMaterialFilter.from_volumes(self.domains, mmv) else: - # TODO: Run volume calculation for cells - domains = self.domains + domains: Sequence[openmc.Cell] = self.domains + + # Check to make sure that each cell is filled with a material and + # that the volume has been set + + # TODO: If volumes are not set, run volume calculation for cells + for cell in domains: + if cell.fill is None: + raise ValueError( + f"Cell {cell.id} is not filled with a materials. " + "Please set the fill material for each cell before " + "running the R2S calculation." + ) + if cell.volume is None: + raise ValueError( + f"Cell {cell.id} does not have a volume set. " + "Please set the volume for each cell before running " + "the R2S calculation." + ) - # Get fluxes and microscopic cross sections on each mesh/material + # Set default keyword arguments for microxs and flux calculation if micro_kwargs is None: micro_kwargs = {} micro_kwargs.setdefault('path_statepoint', output_dir / 'statepoint.h5') + micro_kwargs.setdefault('path_input', output_dir / 'model.xml') # Run neutron transport and get fluxes and micros self.results['fluxes'], self.results['micros'] = get_microxs_and_flux( self.model, domains, **micro_kwargs) - # Export model to output directory - self.model.export_to_model_xml(output_dir / 'model.xml') - # Save flux and micros to file np.save(output_dir / 'fluxes.npy', self.results['fluxes']) for i, micros in enumerate(self.results['micros']): @@ -347,7 +367,8 @@ def step2_activation( This step creates a unique copy of each activation material based on the mesh elements or cells, then solves the depletion equations for each material using the fluxes and microscopic cross sections obtained in the - neutron transport step. + neutron transport step. This step will populate the 'depletion_results' + key in the results dictionary. Parameters ---------- @@ -420,7 +441,8 @@ def step3_photon_transport( creates appropriate photon sources and runs a transport calculation. In mesh-based mode, the sources are created using the mesh material volumes, while in cell-based mode, they are created using bounding boxes - for each cell. + for each cell. This step will populate the 'dose_tallies' key in the + results dictionary. Parameters ---------- @@ -446,6 +468,7 @@ def step3_photon_transport( """ # TODO: Automatically determine bounding box for each cell + # TODO: Voiding/changing materials between neutron/photon steps # Set default run arguments if not provided if run_kwargs is None: @@ -461,6 +484,8 @@ def step3_photon_transport( if settings is not None: self.model.settings = settings + self.results['dose_tallies'] = {} + try: for time_index in cooling_times: # Create decay photon source @@ -501,7 +526,14 @@ def step3_photon_transport( # Run photon transport calculation run_kwargs['cwd'] = Path(output_dir) / f'time_{time_index}' - self.model.run(**run_kwargs) + statepoint_path = self.model.run(**run_kwargs) + + # Store tally results + if dose_tallies is not None: + with openmc.StatePoint(statepoint_path) as sp: + self.results['dose_tallies'][time_index] = [ + sp.tallies[tally.id] for tally in dose_tallies + ] finally: # Restore original settings self.model.settings = original_settings diff --git a/openmc/mesh.py b/openmc/mesh.py index 4a478aa85be..c501a68ed94 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -392,7 +392,7 @@ def material_volumes( # In order to get mesh into model, we temporarily replace the # tallies with a single mesh tally using the current mesh - original_tallies = model.tallies + original_tallies = list(model.tallies) new_tally = openmc.Tally() new_tally.filters = [openmc.MeshFilter(self)] new_tally.scores = ['flux'] From fcd5b1e033b79529b6981b2c04860d06c3966b11 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 18 Jul 2025 16:48:46 +0700 Subject: [PATCH 15/29] Add load_results method --- openmc/deplete/r2s.py | 65 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 64 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index bbb76abf2de..fd902e81ad0 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -1,11 +1,13 @@ from __future__ import annotations from collections.abc import Sequence from datetime import datetime +import json from pathlib import Path import numpy as np import openmc -from . import get_microxs_and_flux, IndependentOperator, PredictorIntegrator +from . import IndependentOperator, PredictorIntegrator +from .microxs import get_microxs_and_flux, MicroXS from .results import Results from ..checkvalue import PathLike from ..mesh import _get_all_materials @@ -484,6 +486,15 @@ def step3_photon_transport( if settings is not None: self.model.settings = settings + # Write out JSON file with dose tally IDs that can be used for loading + # results + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + if dose_tallies is not None: + dose_tally_ids = [tally.id for tally in dose_tallies] + with open(output_dir / 'dose_tally_ids.json', 'w') as f: + json.dump(dose_tally_ids, f) + self.results['dose_tallies'] = {} try: @@ -537,3 +548,55 @@ def step3_photon_transport( finally: # Restore original settings self.model.settings = original_settings + + def load_results(self, path: PathLike): + """Load results from a previous R2S calculation. + + Parameters + ---------- + path : PathLike + Path to the directory containing the R2S calculation results. + + """ + path = Path(path) + + # Load neutron transport results + if self.method == 'mesh-based': + neutron_dir = path / 'neutron_transport' + mmv_file = neutron_dir / 'mesh_material_volumes.npz' + if mmv_file.exists(): + self.results['mesh_material_volumes'] = \ + openmc.MeshMaterialVolumes.from_npz(mmv_file) + fluxes_file = neutron_dir / 'fluxes.npy' + if fluxes_file.exists(): + self.results['fluxes'] = list(np.load(fluxes_file, allow_pickle=True)) + self.results['micros'] = [ + MicroXS.from_csv(neutron_dir / f'micros_{i}.csv') + for i in range(len(self.results['fluxes'])) + ] + + # Load activation results + activation_dir = path / 'activation' + activation_results = activation_dir / 'depletion_results.h5' + if activation_results.exists(): + self.results['depletion_results'] = Results() + + # Load photon transport results + photon_dir = path / 'photon_transport' + + # Load dose tally IDs from JSON file + dose_tally_ids_path = photon_dir / 'dose_tally_ids.json' + if dose_tally_ids_path.exists(): + with dose_tally_ids_path.open('r') as f: + dose_tally_ids = json.load(f) + self.results['dose_tallies'] = {} + + # For each cooling time, load the statepoint and get the dose + # tallies based on dose_tally_ids + for time_dir in photon_dir.glob('time_*'): + time_index = int(time_dir.name.split('_')[1]) + for sp_path in time_dir.glob('statepoint.*.h5'): + with openmc.StatePoint(sp_path) as sp: + self.results['dose_tallies'][time_index] = [ + sp.tallies[tally_id] for tally_id in dose_tally_ids + ] From f1c384f7bfbcd86cb7cfbde5c7bdf304dd588a87 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 18 Jul 2025 20:18:10 +0700 Subject: [PATCH 16/29] Update description of cooling_times --- openmc/deplete/r2s.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index fd902e81ad0..c682b5f7e06 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -214,8 +214,11 @@ def run( source_rates : float or Sequence[float] Source rate in [neutron/sec] for each interval in `timesteps`. cooling_times : Sequence[int] - Sequence of cooling times represented as indices into the array of - `timesteps`. + Sequence of cooling times at which photon transport should be run; + represented as indices into the array of times formed by the + timesteps. For example, if two timesteps are specified, the array of + times would contain three entries, and [2] would indicate computing + photon results at the last time. dose_tallies : Sequence[openmc.Tally], optional A sequence of tallies to be used for dose calculations in the photon transport step. If None, no dose tallies are added. @@ -449,8 +452,11 @@ def step3_photon_transport( Parameters ---------- cooling_times : Sequence[int] - Sequence of cooling times represented as indices into the array of - timesteps used in the activation step. + Sequence of cooling times at which photon transport should be run; + represented as indices into the array of times formed by the + timesteps. For example, if two timesteps are specified, the array of + times would contain three entries, and [2] would indicate computing + photon results at the last time. dose_tallies : Sequence[openmc.Tally], optional A sequence of tallies to be used for dose calculations in the photon transport step. If None, no dose tallies are added. From d2483d2bd6c367bababee11d1788d9e927f95200 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 18 Jul 2025 22:15:50 +0700 Subject: [PATCH 17/29] Update photon-related arguments --- openmc/deplete/r2s.py | 109 +++++++++++++++++++++--------------------- 1 file changed, 55 insertions(+), 54 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index c682b5f7e06..2adf98379f5 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -191,15 +191,15 @@ def run( self, timesteps: Sequence[float] | Sequence[tuple[float, str]], source_rates: float | Sequence[float], - cooling_times: Sequence[int], - dose_tallies: Sequence[openmc.Tally] | None = None, timestep_units: str = 's', + photon_time_indices: Sequence[int] | None = None, + photon_tallies: Sequence[openmc.Tally] | None = None, + photon_settings: openmc.Settings | None = None, output_dir: PathLike | None = None, bounding_boxes: dict[int, openmc.BoundingBox] | None = None, micro_kwargs: dict | None = None, mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, - photon_settings: openmc.Settings | None = None, ): """Run the R2S calculation. @@ -213,19 +213,24 @@ def run( tuples. source_rates : float or Sequence[float] Source rate in [neutron/sec] for each interval in `timesteps`. - cooling_times : Sequence[int] - Sequence of cooling times at which photon transport should be run; - represented as indices into the array of times formed by the - timesteps. For example, if two timesteps are specified, the array of - times would contain three entries, and [2] would indicate computing - photon results at the last time. - dose_tallies : Sequence[openmc.Tally], optional - A sequence of tallies to be used for dose calculations in the photon - transport step. If None, no dose tallies are added. timestep_units : {'s', 'min', 'h', 'd', 'a'}, optional Units for values specified in the `timesteps` argument when passing float values. 's' means seconds, 'min' means minutes, 'h' means hours, 'd' means days, and 'a' means years (Julian). + photon_time_indices : Sequence[int], optional + Sequence of time indices at which photon transport should be run; + represented as indices into the array of times formed by the + timesteps. For example, if two timesteps are specified, the array of + times would contain three entries, and [2] would indicate computing + photon results at the last time. A value of None indicates to run + photon transport for each time. + photon_tallies : Sequence[openmc.Tally], optional + A sequence of tallies to be used in the photon transport step. If + None, no tallies are present. + photon_settings : openmc.Settings, optional + Custom settings to use for the photon transport calculation. If + provided, the original model settings will be temporarily replaced + during the photon transport calculations and restored afterward. output_dir : PathLike, optional Path to directory where R2S calculation outputs will be saved. If not provided, a timestamped directory 'r2s_YYYY-MM-DDTHH-MM-SS' is @@ -247,10 +252,6 @@ def run( Additional keyword arguments passed to :meth:`openmc.Model.run` during the neutron and photon transport step. By default, output is disabled. - photon_settings : openmc.Settings, optional - Custom settings to use for the photon transport calculation. If - provided, the original model settings will be temporarily replaced - during the photon transport calculations and restored afterward. Returns ------- @@ -277,7 +278,7 @@ def run( timesteps, source_rates, timestep_units, output_dir / 'activation' ) self.step3_photon_transport( - cooling_times, dose_tallies, bounding_boxes, + photon_time_indices, photon_tallies, bounding_boxes, output_dir / 'photon_transport', run_kwargs=run_kwargs, settings=photon_settings ) @@ -432,8 +433,8 @@ def step2_activation( def step3_photon_transport( self, - cooling_times: Sequence[int], - dose_tallies: Sequence[openmc.Tally] | None = None, + time_indices: Sequence[int] | None = None, + tallies: Sequence[openmc.Tally] | None = None, bounding_boxes: dict[int, openmc.BoundingBox] | None = None, output_dir: PathLike = 'photon_transport', run_kwargs: dict | None = None, @@ -442,24 +443,25 @@ def step3_photon_transport( """Run the photon transport step. This step performs photon transport calculations using decay photon - sources created from the activated materials. For each cooling time, it - creates appropriate photon sources and runs a transport calculation. In - mesh-based mode, the sources are created using the mesh material + sources created from the activated materials. For each specified time, + it creates appropriate photon sources and runs a transport calculation. + In mesh-based mode, the sources are created using the mesh material volumes, while in cell-based mode, they are created using bounding boxes - for each cell. This step will populate the 'dose_tallies' key in the + for each cell. This step will populate the 'photon_tallies' key in the results dictionary. Parameters ---------- - cooling_times : Sequence[int] - Sequence of cooling times at which photon transport should be run; + time_indices : Sequence[int], optional + Sequence of time indices at which photon transport should be run; represented as indices into the array of times formed by the timesteps. For example, if two timesteps are specified, the array of times would contain three entries, and [2] would indicate computing - photon results at the last time. - dose_tallies : Sequence[openmc.Tally], optional - A sequence of tallies to be used for dose calculations in the photon - transport step. If None, no dose tallies are added. + photon results at the last time. A value of None indicates to run + photon transport for each time. + tallies : Sequence[openmc.Tally], optional + A sequence of tallies to be used in the photon transport step. If + None, no tallies are present. bounding_boxes : dict[int, openmc.BoundingBox], optional Dictionary mapping cell IDs to bounding boxes used for spatial source sampling in cell-based R2S calculations. Required if method @@ -483,28 +485,28 @@ def step3_photon_transport( run_kwargs = {} run_kwargs.setdefault('output', False) - # Add dose tallies to model if provided - if dose_tallies is not None: - self.model.tallies.extend(dose_tallies) + # Add photon tallies to model if provided + if tallies is None: + tallies = [] + self.model.tallies = tallies # Save original settings to restore later original_settings = self.model.settings if settings is not None: self.model.settings = settings - # Write out JSON file with dose tally IDs that can be used for loading + # Write out JSON file with tally IDs that can be used for loading # results output_dir = Path(output_dir) output_dir.mkdir(parents=True, exist_ok=True) - if dose_tallies is not None: - dose_tally_ids = [tally.id for tally in dose_tallies] - with open(output_dir / 'dose_tally_ids.json', 'w') as f: - json.dump(dose_tally_ids, f) + tally_ids = [tally.id for tally in tallies] + with open(output_dir / 'tally_ids.json', 'w') as f: + json.dump(tally_ids, f) - self.results['dose_tallies'] = {} + self.results['photon_tallies'] = {} try: - for time_index in cooling_times: + for time_index in time_indices: # Create decay photon source if self.method == 'mesh-based': self.model.settings.source = get_decay_photon_source_mesh( @@ -546,11 +548,10 @@ def step3_photon_transport( statepoint_path = self.model.run(**run_kwargs) # Store tally results - if dose_tallies is not None: - with openmc.StatePoint(statepoint_path) as sp: - self.results['dose_tallies'][time_index] = [ - sp.tallies[tally.id] for tally in dose_tallies - ] + with openmc.StatePoint(statepoint_path) as sp: + self.results['photon_tallies'][time_index] = [ + sp.tallies[tally.id] for tally in tallies + ] finally: # Restore original settings self.model.settings = original_settings @@ -590,19 +591,19 @@ def load_results(self, path: PathLike): # Load photon transport results photon_dir = path / 'photon_transport' - # Load dose tally IDs from JSON file - dose_tally_ids_path = photon_dir / 'dose_tally_ids.json' - if dose_tally_ids_path.exists(): - with dose_tally_ids_path.open('r') as f: - dose_tally_ids = json.load(f) - self.results['dose_tallies'] = {} + # Load tally IDs from JSON file + tally_ids_path = photon_dir / 'tally_ids.json' + if tally_ids_path.exists(): + with tally_ids_path.open('r') as f: + tally_ids = json.load(f) + self.results['photon_tallies'] = {} - # For each cooling time, load the statepoint and get the dose - # tallies based on dose_tally_ids + # For each photon transport calc, load the statepoint and get the + # tally results based on tally_ids for time_dir in photon_dir.glob('time_*'): time_index = int(time_dir.name.split('_')[1]) for sp_path in time_dir.glob('statepoint.*.h5'): with openmc.StatePoint(sp_path) as sp: - self.results['dose_tallies'][time_index] = [ - sp.tallies[tally_id] for tally_id in dose_tally_ids + self.results['photon_tallies'][time_index] = [ + sp.tallies[tally_id] for tally_id in tally_ids ] From ac07884d960f14590c45d90ce9c5bf46b923f886 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 31 Jul 2025 17:35:58 +0700 Subject: [PATCH 18/29] Add HDF5 methods for MicroXS --- openmc/deplete/microxs.py | 100 ++++++++++++++++++++++++++++++++++++-- openmc/deplete/r2s.py | 9 ++-- openmc/utility_funcs.py | 19 ++++++++ 3 files changed, 120 insertions(+), 8 deletions(-) diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index 50235e58df3..c9962f59bc5 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -8,8 +8,9 @@ from collections.abc import Sequence import shutil from tempfile import TemporaryDirectory -from typing import Union, TypeAlias +from typing import Union, TypeAlias, Self +import h5py import pandas as pd import numpy as np @@ -20,6 +21,7 @@ import openmc from .chain import Chain, REACTIONS, _get_chain from .coupled_operator import _find_cross_sections, _get_nuclides_with_data +from ..utility_funcs import h5py_file_or_group import openmc.lib from openmc.mpi import comm @@ -398,8 +400,7 @@ def from_csv(cls, csv_file, **kwargs): MicroXS """ - if 'float_precision' not in kwargs: - kwargs['float_precision'] = 'round_trip' + kwargs.setdefault('float_precision', 'round_trip') df = pd.read_csv(csv_file, **kwargs) df.set_index(['nuclides', 'reactions', 'groups'], inplace=True) @@ -434,3 +435,96 @@ def to_csv(self, *args, **kwargs): ) df = pd.DataFrame({'xs': self.data.flatten()}, index=multi_index) df.to_csv(*args, **kwargs) + + def to_hdf5(self, group_or_filename: h5py.Group | PathLike, **kwargs): + """Export microscopic cross section data to HDF5 format + + Parameters + ---------- + group_or_filename : h5py.Group or path-like + HDF5 group or filename to write to + kwargs : dict, optional + Keyword arguments to pass to :meth:`h5py.Group.create_dataset`. + Defaults to {'compression': 'lzf'}. + + """ + kwargs.setdefault('compression', 'lzf') + + with h5py_file_or_group(group_or_filename, 'w') as group: + # Store cross section data as 3D dataset + group.create_dataset('data', data=self.data, **kwargs) + + # Store metadata as datasets using string encoding + group.create_dataset('nuclides', data=np.array(self.nuclides, dtype='S')) + group.create_dataset('reactions', data=np.array(self.reactions, dtype='S')) + + @classmethod + def from_hdf5(cls, group_or_filename: h5py.Group | PathLike) -> Self: + """Load data from an HDF5 file + + Parameters + ---------- + group_or_filename : h5py.Group or str or PathLike + HDF5 group or path to HDF5 file. If given as an h5py.Group, the + data is read from that group. If given as a string, it is assumed + to be the filename for the HDF5 file. + + Returns + ------- + MicroXS + """ + + with h5py_file_or_group(group_or_filename, 'r') as group: + # Read data from HDF5 group + data = group['data'][:] + nuclides = [nuc.decode('utf-8') for nuc in group['nuclides'][:]] + reactions = [rxn.decode('utf-8') for rxn in group['reactions'][:]] + + return cls(data, nuclides, reactions) + + +def write_microxs_hdf5( + micros: Sequence[MicroXS], + filename: PathLike, + names: Sequence[str] | None = None, + **kwargs +): + """Write multiple MicroXS objects to an HDF5 file + + Parameters + ---------- + micros : list of MicroXS + List of MicroXS objects + filename : PathLike + Output HDF5 filename + names : list of str, optional + Names for each MicroXS object. If None, uses 'domain_0', 'domain_1', + etc. + **kwargs + Additional keyword arguments passed to :meth:`h5py.Group.create_dataset` + """ + if names is None: + names = [f'domain_{i}' for i in range(len(micros))] + + # Open file once and write all domains using group interface + with h5py.File(filename, 'w') as f: + for microxs, name in zip(micros, names): + group = f.create_group(name) + microxs.to_hdf5(group, **kwargs) + + +def read_microxs_hdf5(filename: PathLike) -> dict[str, MicroXS]: + """Read multiple MicroXS objects from an HDF5 file + + Parameters + ---------- + filename : path-like + HDF5 filename + + Returns + ------- + dict + Dictionary mapping domain names to MicroXS objects + """ + with h5py.File(filename, 'r') as f: + return {name: MicroXS.from_hdf5(group) for name, group in f.items()} diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 2adf98379f5..4598a7f84dd 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -7,7 +7,7 @@ import numpy as np import openmc from . import IndependentOperator, PredictorIntegrator -from .microxs import get_microxs_and_flux, MicroXS +from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5 from .results import Results from ..checkvalue import PathLike from ..mesh import _get_all_materials @@ -358,8 +358,7 @@ def step1_neutron_transport( # Save flux and micros to file np.save(output_dir / 'fluxes.npy', self.results['fluxes']) - for i, micros in enumerate(self.results['micros']): - micros.to_csv(output_dir / f'micros_{i}.csv') + write_microxs_hdf5(self.results['micros'], output_dir / 'micros.h5') def step2_activation( self, @@ -577,9 +576,9 @@ def load_results(self, path: PathLike): fluxes_file = neutron_dir / 'fluxes.npy' if fluxes_file.exists(): self.results['fluxes'] = list(np.load(fluxes_file, allow_pickle=True)) + micros_dict = read_microxs_hdf5(neutron_dir / 'micros.h5') self.results['micros'] = [ - MicroXS.from_csv(neutron_dir / f'micros_{i}.csv') - for i in range(len(self.results['fluxes'])) + micros_dict[f'domain_{i}'] for i in range(len(micros_dict)) ] # Load activation results diff --git a/openmc/utility_funcs.py b/openmc/utility_funcs.py index da9f73b1651..935a589853d 100644 --- a/openmc/utility_funcs.py +++ b/openmc/utility_funcs.py @@ -3,6 +3,8 @@ from pathlib import Path from tempfile import TemporaryDirectory +import h5py + import openmc from .checkvalue import PathLike @@ -57,3 +59,20 @@ def input_path(filename: PathLike) -> Path: return Path(filename).resolve() else: return Path(filename) + + +@contextmanager +def h5py_file_or_group(group_or_filename: PathLike | h5py.Group, *args, **kwargs): + """Context manager for opening an HDF5 file or using an existing group + + Parameters + ---------- + group_or_filename : path-like or h5py.Group + Path to HDF5 file, or group from an existing HDF5 file + + """ + if isinstance(group_or_filename, h5py.Group): + yield group_or_filename + else: + with h5py.File(group_or_filename, *args, **kwargs) as f: + yield f From 6f95d958cf2146b7ef3de881791db4df1ad1bd05 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 Aug 2025 17:27:26 +0700 Subject: [PATCH 19/29] Introduce photon_model argument to replace separate settings/tallies --- openmc/deplete/r2s.py | 159 ++++++++++++++++++------------------------ 1 file changed, 69 insertions(+), 90 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 4598a7f84dd..25e9039f78e 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -1,5 +1,6 @@ from __future__ import annotations from collections.abc import Sequence +import copy from datetime import datetime import json from pathlib import Path @@ -154,19 +155,24 @@ class R2SManager: Parameters ---------- - model : openmc.Model - The OpenMC model containing the geometry and materials. + neutron_model : openmc.Model + The OpenMC model to use for neutron transport. domains : openmc.MeshBase or Sequence[openmc.Cell] The mesh or a sequence of cells that represent the spatial units over which the R2S calculation will be performed. + photon_model : openmc.Model, optional + The OpenMC model to use for photon transport calculations. If None, + a shallow copy of the neutron_model will be created and used. Attributes ---------- domains : openmc.MeshBase or Sequence[openmc.Cell] The mesh or a sequence of cells that represent the spatial units over which the R2S calculation will be performed. - model : openmc.Model - The OpenMC model containing the geometry and materials. + neutron_model : openmc.Model + The OpenMC model used for neutron transport. + photon_model : openmc.Model + The OpenMC model used for photon transport calculations. method : {'mesh-based', 'cell-based'} Indicates whether the R2S calculation uses mesh elements ('mesh-based') as the spatial discetization or a list of a cells ('cell-based'). @@ -176,10 +182,16 @@ class R2SManager: """ def __init__( self, - model: openmc.Model, + neutron_model: openmc.Model, domains: openmc.MeshBase | Sequence[openmc.Cell], + photon_model: openmc.Model | None = None, ): - self.model = model + self.neutron_model = neutron_model + if photon_model is None: + # Create a shallow copy of the neutron model for photon transport + self.photon_model = copy.copy(neutron_model) + else: + self.photon_model = photon_model if isinstance(domains, openmc.MeshBase): self.method = 'mesh-based' else: @@ -193,8 +205,6 @@ def run( source_rates: float | Sequence[float], timestep_units: str = 's', photon_time_indices: Sequence[int] | None = None, - photon_tallies: Sequence[openmc.Tally] | None = None, - photon_settings: openmc.Settings | None = None, output_dir: PathLike | None = None, bounding_boxes: dict[int, openmc.BoundingBox] | None = None, micro_kwargs: dict | None = None, @@ -224,13 +234,6 @@ def run( times would contain three entries, and [2] would indicate computing photon results at the last time. A value of None indicates to run photon transport for each time. - photon_tallies : Sequence[openmc.Tally], optional - A sequence of tallies to be used in the photon transport step. If - None, no tallies are present. - photon_settings : openmc.Settings, optional - Custom settings to use for the photon transport calculation. If - provided, the original model settings will be temporarily replaced - during the photon transport calculations and restored afterward. output_dir : PathLike, optional Path to directory where R2S calculation outputs will be saved. If not provided, a timestamped directory 'r2s_YYYY-MM-DDTHH-MM-SS' is @@ -278,9 +281,8 @@ def run( timesteps, source_rates, timestep_units, output_dir / 'activation' ) self.step3_photon_transport( - photon_time_indices, photon_tallies, bounding_boxes, - output_dir / 'photon_transport', run_kwargs=run_kwargs, - settings=photon_settings + photon_time_indices, bounding_boxes, + output_dir / 'photon_transport', run_kwargs=run_kwargs ) return output_dir @@ -318,7 +320,7 @@ def step1_neutron_transport( if self.method == 'mesh-based': # Compute material volume fractions on the mesh self.results['mesh_material_volumes'] = mmv = \ - self.domains.material_volumes(self.model, **mat_vol_kwargs) + self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs) # Save results to file mmv.save(output_dir / 'mesh_material_volumes.npz') @@ -354,7 +356,7 @@ def step1_neutron_transport( # Run neutron transport and get fluxes and micros self.results['fluxes'], self.results['micros'] = get_microxs_and_flux( - self.model, domains, **micro_kwargs) + self.neutron_model, domains, **micro_kwargs) # Save flux and micros to file np.save(output_dir / 'fluxes.npy', self.results['fluxes']) @@ -397,7 +399,7 @@ def step2_activation( if self.method == 'mesh-based': # Get unique material for each (mesh, material) combination mmv = self.results['mesh_material_volumes'] - self.results['activation_materials'] = get_activation_materials(self.model, mmv) + self.results['activation_materials'] = get_activation_materials(self.neutron_model, mmv) else: # Create unique material for each cell activation_mats = [] @@ -433,11 +435,9 @@ def step2_activation( def step3_photon_transport( self, time_indices: Sequence[int] | None = None, - tallies: Sequence[openmc.Tally] | None = None, bounding_boxes: dict[int, openmc.BoundingBox] | None = None, output_dir: PathLike = 'photon_transport', run_kwargs: dict | None = None, - settings: openmc.Settings | None = None, ): """Run the photon transport step. @@ -458,9 +458,6 @@ def step3_photon_transport( times would contain three entries, and [2] would indicate computing photon results at the last time. A value of None indicates to run photon transport for each time. - tallies : Sequence[openmc.Tally], optional - A sequence of tallies to be used in the photon transport step. If - None, no tallies are present. bounding_boxes : dict[int, openmc.BoundingBox], optional Dictionary mapping cell IDs to bounding boxes used for spatial source sampling in cell-based R2S calculations. Required if method @@ -470,10 +467,6 @@ def step3_photon_transport( run_kwargs : dict, optional Additional keyword arguments passed to :meth:`openmc.Model.run` during the photon transport step. By default, output is disabled. - settings : openmc.Settings, optional - Custom settings to use for the photon transport calculation. If - provided, the original model settings will be temporarily replaced - during the photon transport calculations and restored afterward. """ # TODO: Automatically determine bounding box for each cell @@ -484,76 +477,62 @@ def step3_photon_transport( run_kwargs = {} run_kwargs.setdefault('output', False) - # Add photon tallies to model if provided - if tallies is None: - tallies = [] - self.model.tallies = tallies - - # Save original settings to restore later - original_settings = self.model.settings - if settings is not None: - self.model.settings = settings - # Write out JSON file with tally IDs that can be used for loading # results output_dir = Path(output_dir) output_dir.mkdir(parents=True, exist_ok=True) - tally_ids = [tally.id for tally in tallies] + tally_ids = [tally.id for tally in self.photon_model.tallies] with open(output_dir / 'tally_ids.json', 'w') as f: json.dump(tally_ids, f) self.results['photon_tallies'] = {} - try: - for time_index in time_indices: - # Create decay photon source - if self.method == 'mesh-based': - self.model.settings.source = get_decay_photon_source_mesh( - self.model, - self.domains, - self.results['activation_materials'], - self.results['depletion_results'], - self.results['mesh_material_volumes'], - time_index=time_index, + for time_index in time_indices: + # Create decay photon source + if self.method == 'mesh-based': + self.photon_model.settings.source = get_decay_photon_source_mesh( + self.neutron_model, + self.domains, + self.results['activation_materials'], + self.results['depletion_results'], + self.results['mesh_material_volumes'], + time_index=time_index, + ) + else: + sources = [] + results = self.results['depletion_results'] + for cell, original_mat in zip(self.domains, self.results['activation_materials']): + bounding_box = bounding_boxes[cell.id] + + # Get activated material composition + activated_mat = results[time_index].get_material(str(original_mat.id)) + + # Create decay photon source source + space = openmc.stats.Box(*bounding_box) + energy = activated_mat.get_decay_photon_energy() + source = openmc.IndependentSource( + space=space, + energy=energy, + particle='photon', + strength=energy.integral(), + domains=[cell] ) - else: - sources = [] - results = self.results['depletion_results'] - for cell, original_mat in zip(self.domains, self.results['activation_materials']): - bounding_box = bounding_boxes[cell.id] - - # Get activated material composition - activated_mat = results[time_index].get_material(str(original_mat.id)) - - # Create decay photon source source - space = openmc.stats.Box(*bounding_box) - energy = activated_mat.get_decay_photon_energy() - source = openmc.IndependentSource( - space=space, - energy=energy, - particle='photon', - strength=energy.integral(), - domains=[cell] - ) - sources.append(source) - self.model.settings.source = sources - - # Convert time_index (which may be negative) to a normal index - if time_index < 0: - time_index = len(self.results['depletion_results']) + time_index - - # Run photon transport calculation - run_kwargs['cwd'] = Path(output_dir) / f'time_{time_index}' - statepoint_path = self.model.run(**run_kwargs) - - # Store tally results - with openmc.StatePoint(statepoint_path) as sp: - self.results['photon_tallies'][time_index] = [ - sp.tallies[tally.id] for tally in tallies - ] - finally: - # Restore original settings - self.model.settings = original_settings + sources.append(source) + self.photon_model.settings.source = sources + + # Convert time_index (which may be negative) to a normal index + if time_index < 0: + time_index = len(self.results['depletion_results']) + time_index + + # Run photon transport calculation + run_kwargs['cwd'] = Path(output_dir) / f'time_{time_index}' + statepoint_path = self.photon_model.run(**run_kwargs) + + # Store tally results + with openmc.StatePoint(statepoint_path) as sp: + self.results['photon_tallies'][time_index] = [ + sp.tallies[tally.id] for tally in self.photon_model.tallies + ] def load_results(self, path: PathLike): """Load results from a previous R2S calculation. From 31207ed2eef01b7dbd75af8886d6f0b9e6696153 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 1 Aug 2025 17:38:47 +0700 Subject: [PATCH 20/29] Use copy.copy on individual attributes of Model --- openmc/deplete/r2s.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 25e9039f78e..2a41078d5d8 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -189,7 +189,13 @@ def __init__( self.neutron_model = neutron_model if photon_model is None: # Create a shallow copy of the neutron model for photon transport - self.photon_model = copy.copy(neutron_model) + self.photon_model = openmc.Model( + geometry=copy.copy(neutron_model.geometry), + materials=copy.copy(neutron_model.materials), + settings=copy.copy(neutron_model.settings), + tallies=copy.copy(neutron_model.tallies), + plots=copy.copy(neutron_model.plots), + ) else: self.photon_model = photon_model if isinstance(domains, openmc.MeshBase): From 6d10e126af05f3ba5b278b3303387333b783deed Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 2 Aug 2025 11:02:55 +0700 Subject: [PATCH 21/29] Use photon mesh material volumes to determine if materials changed --- openmc/deplete/r2s.py | 44 +++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 42 insertions(+), 2 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 2a41078d5d8..d609324e261 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -68,6 +68,7 @@ def get_decay_photon_source_mesh( materials: list[openmc.Material], results: Results, mat_vols: openmc.MeshMaterialVolumes, + photon_mat_vols: openmc.MeshMaterialVolumes | None = None, time_index: int = -1, ) -> list[openmc.MeshSource]: """Create decay photon source for a mesh-based R2S calculation. @@ -91,7 +92,13 @@ def get_decay_photon_source_mesh( The results object containing the depletion results for the materials. mat_vols : openmc.MeshMaterialVolumes The mesh material volumes object containing the materials and their - volumes for each mesh element. + volumes for each mesh element from the neutron transport step. + photon_mat_vols : openmc.MeshMaterialVolumes, optional + The mesh material volumes object for the photon model. If provided, + sources will only be created for mesh element-material combinations + that exist in both the neutron and photon models. + time_index : int, optional + Time index for the decay photon source. Default is -1 (last time). Returns ------- @@ -116,11 +123,26 @@ def get_decay_photon_source_mesh( n_elements = mat_vols.num_elements for index_elem in range(n_elements): + # Determine which materials exist in the photon model for this element + if photon_mat_vols is not None: + photon_materials = { + mat_id + for mat_id, _ in photon_mat_vols.by_element(index_elem) + if mat_id is not None + } + else: + photon_materials = set() + for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): # Skip void volume if mat_id is None: continue + # Skip if this material doesn't exist in photon model + if mat_id not in photon_materials: + index_mat += 1 + continue + # Check whether a new MeshSource object is needed if j >= len(source_lists): source_lists.append([null_source]*n_elements) @@ -476,7 +498,7 @@ def step3_photon_transport( """ # TODO: Automatically determine bounding box for each cell - # TODO: Voiding/changing materials between neutron/photon steps + # TODO: For cell-based calculations, account for material changes in photon model # Set default run arguments if not provided if run_kwargs is None: @@ -487,6 +509,16 @@ def step3_photon_transport( # results output_dir = Path(output_dir) output_dir.mkdir(parents=True, exist_ok=True) + + # For mesh-based calculations, compute material volume fractions for + # the photon model to account for potential material changes + if self.method == 'mesh-based': + self.results['mesh_material_volumes_photon'] = photon_mmv = \ + self.domains.material_volumes(self.photon_model) + + # Save photon MMV results to file + photon_mmv.save(output_dir / 'mesh_material_volumes.npz') + tally_ids = [tally.id for tally in self.photon_model.tallies] with open(output_dir / 'tally_ids.json', 'w') as f: json.dump(tally_ids, f) @@ -502,6 +534,7 @@ def step3_photon_transport( self.results['activation_materials'], self.results['depletion_results'], self.results['mesh_material_volumes'], + self.results['mesh_material_volumes_photon'], time_index=time_index, ) else: @@ -575,6 +608,13 @@ def load_results(self, path: PathLike): # Load photon transport results photon_dir = path / 'photon_transport' + # Load photon mesh material volumes if they exist (for mesh-based calculations) + if self.method == 'mesh-based': + photon_mmv_file = photon_dir / 'mesh_material_volumes.npz' + if photon_mmv_file.exists(): + self.results['mesh_material_volumes_photon'] = \ + openmc.MeshMaterialVolumes.from_npz(photon_mmv_file) + # Load tally IDs from JSON file tally_ids_path = photon_dir / 'tally_ids.json' if tally_ids_path.exists(): From e79ca3b37094ee65ea6002159de6fda2bab4a369 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 2 Aug 2025 11:18:48 +0700 Subject: [PATCH 22/29] Skip photon MMV if photon and neutron model are the same --- openmc/deplete/r2s.py | 234 ++++++++++++++++++++---------------------- 1 file changed, 112 insertions(+), 122 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index d609324e261..dfcb42e6a45 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -62,109 +62,6 @@ def get_activation_materials( return materials -def get_decay_photon_source_mesh( - model: openmc.Model, - mesh: openmc.MeshBase, - materials: list[openmc.Material], - results: Results, - mat_vols: openmc.MeshMaterialVolumes, - photon_mat_vols: openmc.MeshMaterialVolumes | None = None, - time_index: int = -1, -) -> list[openmc.MeshSource]: - """Create decay photon source for a mesh-based R2S calculation. - - This function creates N :class:`MeshSource` objects where N is the maximum - number of unique materials that appears in a single mesh element. For each - mesh element-material combination, and IndependentSource instance is created - with a spatial constraint limited the sampled decay photons to the correct - region. - - Parameters - ---------- - model : openmc.Model - The full model containing the geometry and materials. - mesh : openmc.MeshBase - The mesh object defining the spatial regions over which activation is - performed. - materials : list of openmc.Material - List of materials that are activated in the model. - results : openmc.deplete.Results - The results object containing the depletion results for the materials. - mat_vols : openmc.MeshMaterialVolumes - The mesh material volumes object containing the materials and their - volumes for each mesh element from the neutron transport step. - photon_mat_vols : openmc.MeshMaterialVolumes, optional - The mesh material volumes object for the photon model. If provided, - sources will only be created for mesh element-material combinations - that exist in both the neutron and photon models. - time_index : int, optional - Time index for the decay photon source. Default is -1 (last time). - - Returns - ------- - list of openmc.MeshSource - A list of MeshSource objects, each containing IndependentSource - instances for the decay photons in the corresponding mesh element. - - """ - mat_dict = _get_all_materials(model) - - # Some MeshSource objects will have empty positions; create a "null source" - # that is used for this case - null_source = openmc.IndependentSource(particle='photon', strength=0.0) - - # List to hold sources for each MeshSource (length = N) - source_lists = [] - - # Index in the overall list of activated materials - index_mat = 0 - - # Total number of mesh elements - n_elements = mat_vols.num_elements - - for index_elem in range(n_elements): - # Determine which materials exist in the photon model for this element - if photon_mat_vols is not None: - photon_materials = { - mat_id - for mat_id, _ in photon_mat_vols.by_element(index_elem) - if mat_id is not None - } - else: - photon_materials = set() - - for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): - # Skip void volume - if mat_id is None: - continue - - # Skip if this material doesn't exist in photon model - if mat_id not in photon_materials: - index_mat += 1 - continue - - # Check whether a new MeshSource object is needed - if j >= len(source_lists): - source_lists.append([null_source]*n_elements) - - # Get activated material composition - original_mat = materials[index_mat] - activated_mat = results[time_index].get_material(str(original_mat.id)) - - # Create decay photon source source - energy = activated_mat.get_decay_photon_energy() - source_lists[j][index_elem] = openmc.IndependentSource( - energy=energy, - particle='photon', - strength=energy.integral(), - constraints={'domains': [mat_dict[mat_id]]} - ) - - # Increment index of activated material - index_mat += 1 - - # Return list of mesh sources - return [openmc.MeshSource(mesh, sources) for sources in source_lists] class R2SManager: @@ -277,8 +174,7 @@ def run( transport step. mat_vol_kwargs : dict, optional Additional keyword arguments passed to - :meth:`openmc.MeshBase.material_volumes` during the neutron - transport step. + :meth:`openmc.MeshBase.material_volumes`. run_kwargs : dict, optional Additional keyword arguments passed to :meth:`openmc.Model.run` during the neutron and photon transport step. By default, output is @@ -309,8 +205,8 @@ def run( timesteps, source_rates, timestep_units, output_dir / 'activation' ) self.step3_photon_transport( - photon_time_indices, bounding_boxes, - output_dir / 'photon_transport', run_kwargs=run_kwargs + photon_time_indices, bounding_boxes, output_dir / 'photon_transport', + mat_vol_kwargs=mat_vol_kwargs, run_kwargs=run_kwargs ) return output_dir @@ -465,6 +361,7 @@ def step3_photon_transport( time_indices: Sequence[int] | None = None, bounding_boxes: dict[int, openmc.BoundingBox] | None = None, output_dir: PathLike = 'photon_transport', + mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, ): """Run the photon transport step. @@ -492,6 +389,9 @@ def step3_photon_transport( is 'cell-based'. output_dir : PathLike, optional Path to directory where photon transport outputs will be saved. + mat_vol_kwargs : dict, optional + Additional keyword arguments passed to + :meth:`openmc.MeshBase.material_volumes`. run_kwargs : dict, optional Additional keyword arguments passed to :meth:`openmc.Model.run` during the photon transport step. By default, output is disabled. @@ -510,14 +410,18 @@ def step3_photon_transport( output_dir = Path(output_dir) output_dir.mkdir(parents=True, exist_ok=True) - # For mesh-based calculations, compute material volume fractions for - # the photon model to account for potential material changes + # For mesh-based calculations, compute material volume fractions for the + # photon model if it is different from the neutron model to account for + # potential material changes if self.method == 'mesh-based': - self.results['mesh_material_volumes_photon'] = photon_mmv = \ - self.domains.material_volumes(self.photon_model) + neutron_univ = self.neutron_model.geometry.root_universe + photon_univ = self.photon_model.geometry.root_universe + if neutron_univ != photon_univ: + self.results['mesh_material_volumes_photon'] = photon_mmv = \ + self.domains.material_volumes(self.photon_model, **mat_vol_kwargs) - # Save photon MMV results to file - photon_mmv.save(output_dir / 'mesh_material_volumes.npz') + # Save photon MMV results to file + photon_mmv.save(output_dir / 'mesh_material_volumes.npz') tally_ids = [tally.id for tally in self.photon_model.tallies] with open(output_dir / 'tally_ids.json', 'w') as f: @@ -528,15 +432,8 @@ def step3_photon_transport( for time_index in time_indices: # Create decay photon source if self.method == 'mesh-based': - self.photon_model.settings.source = get_decay_photon_source_mesh( - self.neutron_model, - self.domains, - self.results['activation_materials'], - self.results['depletion_results'], - self.results['mesh_material_volumes'], - self.results['mesh_material_volumes_photon'], - time_index=time_index, - ) + self.photon_model.settings.source = \ + self.get_decay_photon_source_mesh(time_index) else: sources = [] results = self.results['depletion_results'] @@ -573,6 +470,99 @@ def step3_photon_transport( sp.tallies[tally.id] for tally in self.photon_model.tallies ] + def get_decay_photon_source_mesh( + self, + time_index: int = -1 + ) -> list[openmc.MeshSource]: + """Create decay photon source for a mesh-based calculation. + + This function creates N :class:`MeshSource` objects where N is the + maximum number of unique materials that appears in a single mesh + element. For each mesh element-material combination, and + IndependentSource instance is created with a spatial constraint limited + the sampled decay photons to the correct region. + + When the photon transport model is different from the neutron model, the + photon MeshMaterialVolumes is used to determine whether an (element, + material) combination exists in the photon model. + + Parameters + ---------- + time_index : int, optional + Time index for the decay photon source. Default is -1 (last time). + + Returns + ------- + list of openmc.MeshSource + A list of MeshSource objects, each containing IndependentSource + instances for the decay photons in the corresponding mesh element. + + """ + mat_dict = _get_all_materials(self.neutron_model) + + # Some MeshSource objects will have empty positions; create a "null source" + # that is used for this case + null_source = openmc.IndependentSource(particle='photon', strength=0.0) + + # List to hold sources for each MeshSource (length = N) + source_lists = [] + + # Index in the overall list of activated materials + index_mat = 0 + + # Get various results from previous steps + mat_vols = self.results['mesh_material_volumes'] + materials = self.results['activation_materials'] + results = self.results['depletion_results'] + photon_mat_vols = self.results.get('mesh_material_volumes_photon') + + # Total number of mesh elements + n_elements = mat_vols.num_elements + + for index_elem in range(n_elements): + # Determine which materials exist in the photon model for this element + if photon_mat_vols is not None: + photon_materials = { + mat_id + for mat_id, _ in photon_mat_vols.by_element(index_elem) + if mat_id is not None + } + else: + photon_materials = set() + + for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): + # Skip void volume + if mat_id is None: + continue + + # Skip if this material doesn't exist in photon model + if mat_id not in photon_materials: + index_mat += 1 + continue + + # Check whether a new MeshSource object is needed + if j >= len(source_lists): + source_lists.append([null_source]*n_elements) + + # Get activated material composition + original_mat = materials[index_mat] + activated_mat = results[time_index].get_material(str(original_mat.id)) + + # Create decay photon source source + energy = activated_mat.get_decay_photon_energy() + source_lists[j][index_elem] = openmc.IndependentSource( + energy=energy, + particle='photon', + strength=energy.integral(), + constraints={'domains': [mat_dict[mat_id]]} + ) + + # Increment index of activated material + index_mat += 1 + + # Return list of mesh sources + return [openmc.MeshSource(self.domains, sources) for sources in source_lists] + def load_results(self, path: PathLike): """Load results from a previous R2S calculation. From 97a72e4889fd303fd9afb310f5af1db894ea936a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 2 Aug 2025 11:41:20 +0700 Subject: [PATCH 23/29] Check for changed material assignments for cell-based calculations --- openmc/deplete/r2s.py | 48 ++++++++++++++++++++++++++++--------------- 1 file changed, 32 insertions(+), 16 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index dfcb42e6a45..2107e033339 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -62,8 +62,6 @@ def get_activation_materials( return materials - - class R2SManager: """Manager for Rigorous 2-Step (R2S) method calculations. @@ -72,6 +70,13 @@ class R2SManager: activation materials and decay photon sources based on the mesh/cells and materials in the OpenMC model. + This class supports the use of a different models for the neutron and photon + transport calculation. However, for cell-based calculations, it assumes that + the only changes in the model are material assignments. For mesh-based + calculations, it checks material assignments in the photon model and any + element--material combinations that don't appear in the photon model are + skipped. + Parameters ---------- neutron_model : openmc.Model @@ -80,8 +85,8 @@ class R2SManager: The mesh or a sequence of cells that represent the spatial units over which the R2S calculation will be performed. photon_model : openmc.Model, optional - The OpenMC model to use for photon transport calculations. If None, - a shallow copy of the neutron_model will be created and used. + The OpenMC model to use for photon transport calculations. If None, a + shallow copy of the neutron_model will be created and used. Attributes ---------- @@ -398,7 +403,6 @@ def step3_photon_transport( """ # TODO: Automatically determine bounding box for each cell - # TODO: For cell-based calculations, account for material changes in photon model # Set default run arguments if not provided if run_kwargs is None: @@ -410,18 +414,20 @@ def step3_photon_transport( output_dir = Path(output_dir) output_dir.mkdir(parents=True, exist_ok=True) + # Check whether the photon model is different + neutron_univ = self.neutron_model.geometry.root_universe + photon_univ = self.photon_model.geometry.root_universe + different_photon_model = (neutron_univ != photon_univ) + # For mesh-based calculations, compute material volume fractions for the # photon model if it is different from the neutron model to account for # potential material changes - if self.method == 'mesh-based': - neutron_univ = self.neutron_model.geometry.root_universe - photon_univ = self.photon_model.geometry.root_universe - if neutron_univ != photon_univ: - self.results['mesh_material_volumes_photon'] = photon_mmv = \ - self.domains.material_volumes(self.photon_model, **mat_vol_kwargs) + if self.method == 'mesh-based' and different_photon_model: + self.results['mesh_material_volumes_photon'] = photon_mmv = \ + self.domains.material_volumes(self.photon_model, **mat_vol_kwargs) - # Save photon MMV results to file - photon_mmv.save(output_dir / 'mesh_material_volumes.npz') + # Save photon MMV results to file + photon_mmv.save(output_dir / 'mesh_material_volumes.npz') tally_ids = [tally.id for tally in self.photon_model.tallies] with open(output_dir / 'tally_ids.json', 'w') as f: @@ -429,6 +435,10 @@ def step3_photon_transport( self.results['photon_tallies'] = {} + # Get dictionary of cells in the photon model + if different_photon_model: + photon_cells = self.photon_model.geometry.get_all_cells() + for time_index in time_indices: # Create decay photon source if self.method == 'mesh-based': @@ -438,6 +448,14 @@ def step3_photon_transport( sources = [] results = self.results['depletion_results'] for cell, original_mat in zip(self.domains, self.results['activation_materials']): + # Skip if the cell is not in the photon model or the + # material has changed + if different_photon_model: + if cell.id not in photon_cells or \ + cell.fill.id != photon_cells[cell.id].fill.id: + continue + + # Get bounding box for the cell bounding_box = bounding_boxes[cell.id] # Get activated material composition @@ -527,8 +545,6 @@ def get_decay_photon_source_mesh( for mat_id, _ in photon_mat_vols.by_element(index_elem) if mat_id is not None } - else: - photon_materials = set() for j, (mat_id, _) in enumerate(mat_vols.by_element(index_elem)): # Skip void volume @@ -536,7 +552,7 @@ def get_decay_photon_source_mesh( continue # Skip if this material doesn't exist in photon model - if mat_id not in photon_materials: + if photon_mat_vols is not None and mat_id not in photon_materials: index_mat += 1 continue From 89a1e2760b55f51f436ba35b104ffc7f85ea0c84 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 6 Sep 2025 18:25:01 -0500 Subject: [PATCH 24/29] Add R2S tests and make various fixes in R2SManager --- openmc/deplete/r2s.py | 75 +++++++++++++---- tests/unit_tests/test_r2s.py | 152 +++++++++++++++++++++++++++++++++++ 2 files changed, 210 insertions(+), 17 deletions(-) create mode 100644 tests/unit_tests/test_r2s.py diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 2107e033339..8477fa3f961 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -16,7 +16,7 @@ def get_activation_materials( model: openmc.Model, mmv: openmc.MeshMaterialVolumes -) -> list[openmc.Material]: +) -> openmc.Materials: """Get a list of activation materials for each mesh element/material. When performing a mesh-based R2S calculation, a unique material is needed @@ -35,7 +35,7 @@ def get_activation_materials( Returns ------- - list of openmc.Material + openmc.Materials A list of materials, each corresponding to a unique mesh element and material combination. @@ -50,7 +50,7 @@ def get_activation_materials( material_dict = _get_all_materials(model) # Create a new activation material for each element-material combination - materials = [] + materials = openmc.Materials() for elem, mat_id, vol in zip(elems, mat_ids, volumes): mat = material_dict[mat_id] new_mat = mat.clone() @@ -137,9 +137,11 @@ def run( photon_time_indices: Sequence[int] | None = None, output_dir: PathLike | None = None, bounding_boxes: dict[int, openmc.BoundingBox] | None = None, + chain_file: PathLike | None = None, micro_kwargs: dict | None = None, mat_vol_kwargs: dict | None = None, run_kwargs: dict | None = None, + operator_kwargs: dict | None = None, ): """Run the R2S calculation. @@ -173,6 +175,9 @@ def run( Dictionary mapping cell IDs to bounding boxes used for spatial source sampling in cell-based R2S calculations. Required if method is 'cell-based'. + chain_file : PathLike, optional + Path to the depletion chain XML file to use during activation. If + not provided, the default configured chain file will be used. micro_kwargs : dict, optional Additional keyword arguments passed to :func:`openmc.deplete.get_microxs_and_flux` during the neutron @@ -184,6 +189,9 @@ def run( Additional keyword arguments passed to :meth:`openmc.Model.run` during the neutron and photon transport step. By default, output is disabled. + operator_kwargs : dict, optional + Additional keyword arguments passed to + :class:`openmc.deplete.IndependentOperator`. Returns ------- @@ -200,14 +208,21 @@ def run( micro_kwargs = {} if run_kwargs is None: run_kwargs = {} + if operator_kwargs is None: + operator_kwargs = {} run_kwargs.setdefault('output', False) micro_kwargs.setdefault('run_kwargs', run_kwargs) + # If a chain file is provided, prefer it for steps 1 and 2 + if chain_file is not None: + micro_kwargs.setdefault('chain_file', chain_file) + operator_kwargs.setdefault('chain_file', chain_file) self.step1_neutron_transport( output_dir / 'neutron_transport', mat_vol_kwargs, micro_kwargs ) self.step2_activation( - timesteps, source_rates, timestep_units, output_dir / 'activation' + timesteps, source_rates, timestep_units, output_dir / 'activation', + operator_kwargs=operator_kwargs ) self.step3_photon_transport( photon_time_indices, bounding_boxes, output_dir / 'photon_transport', @@ -248,6 +263,8 @@ def step1_neutron_transport( if self.method == 'mesh-based': # Compute material volume fractions on the mesh + if mat_vol_kwargs is None: + mat_vol_kwargs = {} self.results['mesh_material_volumes'] = mmv = \ self.domains.material_volumes(self.neutron_model, **mat_vol_kwargs) @@ -297,14 +314,15 @@ def step2_activation( source_rates: float | Sequence[float], timestep_units: str = 's', output_dir: PathLike = 'activation', + operator_kwargs: dict | None = None, ): """Run the activation step. This step creates a unique copy of each activation material based on the mesh elements or cells, then solves the depletion equations for each material using the fluxes and microscopic cross sections obtained in the - neutron transport step. This step will populate the 'depletion_results' - key in the results dictionary. + neutron transport step. This step will populate the 'depletion_results' + and 'activation_materials' keys in the results dictionary. Parameters ---------- @@ -323,6 +341,9 @@ def step2_activation( output_dir : PathLike, optional Path to directory where activation calculation outputs will be saved. + operator_kwargs : dict, optional + Additional keyword arguments passed to + :class:`openmc.deplete.IndependentOperator`. """ if self.method == 'mesh-based': @@ -331,7 +352,7 @@ def step2_activation( self.results['activation_materials'] = get_activation_materials(self.neutron_model, mmv) else: # Create unique material for each cell - activation_mats = [] + activation_mats = openmc.Materials() for cell in self.domains: mat = cell.fill.clone() mat.name = f'Cell {cell.id}' @@ -340,21 +361,27 @@ def step2_activation( activation_mats.append(mat) self.results['activation_materials'] = activation_mats + # Save activation materials to file + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + self.results['activation_materials'].export_to_xml( + output_dir / 'materials.xml') + # Create depletion operator for the activation materials + if operator_kwargs is None: + operator_kwargs = {} + operator_kwargs.setdefault('normalization_mode', 'source-rate') op = IndependentOperator( self.results['activation_materials'], self.results['fluxes'], self.results['micros'], - normalization_mode='source-rate' + **operator_kwargs ) # Create time integrator and solve depletion equations integrator = PredictorIntegrator( - op, timesteps, source_rates=source_rates, - timestep_units=timestep_units + op, timesteps, source_rates=source_rates, timestep_units=timestep_units ) - output_dir = Path(output_dir) - output_dir.mkdir(parents=True, exist_ok=True) output_path = output_dir / 'depletion_results.h5' integrator.integrate(final_step=False, path=output_path) @@ -403,6 +430,9 @@ def step3_photon_transport( """ # TODO: Automatically determine bounding box for each cell + if bounding_boxes is None and self.method == 'cell-based': + raise ValueError("bounding_boxes must be provided for cell-based " + "R2S calculations.") # Set default run arguments if not provided if run_kwargs is None: @@ -414,6 +444,11 @@ def step3_photon_transport( output_dir = Path(output_dir) output_dir.mkdir(parents=True, exist_ok=True) + # Get default time indices if not provided + if time_indices is None: + n_steps = len(self.results['depletion_results']) + time_indices = list(range(n_steps)) + # Check whether the photon model is different neutron_univ = self.neutron_model.geometry.root_universe photon_univ = self.photon_model.geometry.root_universe @@ -464,12 +499,13 @@ def step3_photon_transport( # Create decay photon source source space = openmc.stats.Box(*bounding_box) energy = activated_mat.get_decay_photon_energy() + strength = energy.integral() if energy is not None else 0.0 source = openmc.IndependentSource( space=space, energy=energy, particle='photon', - strength=energy.integral(), - domains=[cell] + strength=strength, + constraints={'domains': [cell]} ) sources.append(source) self.photon_model.settings.source = sources @@ -566,10 +602,11 @@ def get_decay_photon_source_mesh( # Create decay photon source source energy = activated_mat.get_decay_photon_energy() + strength = energy.integral() if energy is not None else 0.0 source_lists[j][index_elem] = openmc.IndependentSource( energy=energy, particle='photon', - strength=energy.integral(), + strength=strength, constraints={'domains': [mat_dict[mat_id]]} ) @@ -591,8 +628,8 @@ def load_results(self, path: PathLike): path = Path(path) # Load neutron transport results + neutron_dir = path / 'neutron_transport' if self.method == 'mesh-based': - neutron_dir = path / 'neutron_transport' mmv_file = neutron_dir / 'mesh_material_volumes.npz' if mmv_file.exists(): self.results['mesh_material_volumes'] = \ @@ -609,7 +646,11 @@ def load_results(self, path: PathLike): activation_dir = path / 'activation' activation_results = activation_dir / 'depletion_results.h5' if activation_results.exists(): - self.results['depletion_results'] = Results() + self.results['depletion_results'] = Results(activation_results) + activation_mats_file = activation_dir / 'materials.xml' + if activation_mats_file.exists(): + self.results['activation_materials'] = \ + openmc.Materials.from_xml(activation_mats_file) # Load photon transport results photon_dir = path / 'photon_transport' diff --git a/tests/unit_tests/test_r2s.py b/tests/unit_tests/test_r2s.py new file mode 100644 index 00000000000..a94f85c8c08 --- /dev/null +++ b/tests/unit_tests/test_r2s.py @@ -0,0 +1,152 @@ +from pathlib import Path + +import pytest +import openmc +from openmc.deplete import Chain, R2SManager + + +@pytest.fixture +def simple_model_and_mesh(tmp_path): + # Define two materials: water and Ni + h2o = openmc.Material() + h2o.add_nuclide("H1", 2.0) + h2o.add_nuclide("O16", 1.0) + h2o.set_density("g/cm3", 1.0) + nickel = openmc.Material() + nickel.add_element("Ni", 1.0) + nickel.set_density("g/cm3", 4.0) + + # Geometry: two half-spaces split by x=0 plane + left = openmc.XPlane(0.0) + x_min = openmc.XPlane(-10.0, boundary_type='vacuum') + x_max = openmc.XPlane(10.0, boundary_type='vacuum') + y_min = openmc.YPlane(-10.0, boundary_type='vacuum') + y_max = openmc.YPlane(10.0, boundary_type='vacuum') + z_min = openmc.ZPlane(-10.0, boundary_type='vacuum') + z_max = openmc.ZPlane(10.0, boundary_type='vacuum') + + c1 = openmc.Cell(fill=h2o, region=+x_min & -left & +y_min & -y_max & +z_min & -z_max) + c2 = openmc.Cell(fill=nickel, region=+left & -x_max & +y_min & -y_max & +z_min & -z_max) + c1.volume = 4000.0 + c2.volume = 4000.0 + geometry = openmc.Geometry([c1, c2]) + + # Simple settings with a point source + settings = openmc.Settings() + settings.batches = 10 + settings.particles = 1000 + settings.run_mode = 'fixed source' + settings.source = openmc.IndependentSource() + model = openmc.Model(geometry, settings=settings) + + mesh = openmc.RegularMesh() + mesh.lower_left = (-10.0, -10.0, -10.0) + mesh.upper_right = (10.0, 10.0, 10.0) + mesh.dimension = (1, 1, 1) + return model, (c1, c2), mesh + + +def test_r2s_mesh_expected_output(simple_model_and_mesh, tmp_path): + model, (c1, c2), mesh = simple_model_and_mesh + + # Use mesh-based domains + r2s = R2SManager(model, mesh) + + # Use custom reduced chain file for Ni + chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml") + + # Run R2S calculation + outdir = r2s.run( + timesteps=[(1.0, 'd')], + source_rates=[1.0], + photon_time_indices=[1], + output_dir=tmp_path, + chain_file=chain, + ) + + # Check directories and files exist + nt = Path(outdir) / 'neutron_transport' + assert (nt / 'fluxes.npy').exists() + assert (nt / 'micros.h5').exists() + assert (nt / 'mesh_material_volumes.npz').exists() + act = Path(outdir) / 'activation' + assert (act / 'depletion_results.h5').exists() + pt = Path(outdir) / 'photon_transport' + assert (pt / 'tally_ids.json').exists() + assert (pt / 'time_1' / 'statepoint.10.h5').exists() + + # Basic results structure checks + assert len(r2s.results['fluxes']) == 2 + assert len(r2s.results['micros']) == 2 + assert len(r2s.results['mesh_material_volumes']) == 2 + assert len(r2s.results['activation_materials']) == 2 + assert len(r2s.results['depletion_results']) == 2 + + # Check activation materials + amats = r2s.results['activation_materials'] + assert all(m.depletable for m in amats) + # Volumes preserved + assert {m.volume for m in amats} == {c1.volume, c2.volume} + + # Check loading results + r2s_loaded = R2SManager(model, mesh) + r2s_loaded.load_results(outdir) + assert len(r2s_loaded.results['fluxes']) == 2 + assert len(r2s_loaded.results['micros']) == 2 + assert len(r2s_loaded.results['mesh_material_volumes']) == 2 + assert len(r2s_loaded.results['activation_materials']) == 2 + assert len(r2s_loaded.results['depletion_results']) == 2 + + +def test_r2s_cell_expected_output(simple_model_and_mesh, tmp_path): + model, (c1, c2), _ = simple_model_and_mesh + + # Use cell-based domains + r2s = R2SManager(model, [c1, c2]) + + # Use custom reduced chain file for Ni + chain = Chain.from_xml(Path(__file__).parents[1] / "chain_ni.xml") + + # Run R2S calculation + bounding_boxes = {c1.id: c1.bounding_box, c2.id: c2.bounding_box} + outdir = r2s.run( + timesteps=[(1.0, 'd')], + source_rates=[1.0], + photon_time_indices=[1], + output_dir=tmp_path, + bounding_boxes=bounding_boxes, + chain_file=chain + ) + + # Check directories and files exist + nt = Path(outdir) / 'neutron_transport' + assert (nt / 'fluxes.npy').exists() + assert (nt / 'micros.h5').exists() + act = Path(outdir) / 'activation' + assert (act / 'depletion_results.h5').exists() + pt = Path(outdir) / 'photon_transport' + assert (pt / 'tally_ids.json').exists() + assert (pt / 'time_1' / 'statepoint.10.h5').exists() + + # Basic results structure checks + assert len(r2s.results['fluxes']) == 2 + assert len(r2s.results['micros']) == 2 + assert len(r2s.results['activation_materials']) == 2 + assert len(r2s.results['depletion_results']) == 2 + + # Check activation materials + amats = r2s.results['activation_materials'] + assert all(m.depletable for m in amats) + # Names include cell IDs + assert any(f"Cell {c1.id}" in m.name for m in amats) + assert any(f"Cell {c2.id}" in m.name for m in amats) + # Volumes preserved + assert {m.volume for m in amats} == {c1.volume, c2.volume} + + # Check loading results + r2s_loaded = R2SManager(model, [c1, c2]) + r2s_loaded.load_results(outdir) + assert len(r2s_loaded.results['fluxes']) == 2 + assert len(r2s_loaded.results['micros']) == 2 + assert len(r2s_loaded.results['activation_materials']) == 2 + assert len(r2s_loaded.results['depletion_results']) == 2 From c1cd9a8a298ad4cfaccdcf86c211adbd2bcbbb95 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 6 Sep 2025 19:20:24 -0500 Subject: [PATCH 25/29] Add user documentation for R2SManager --- docs/source/pythonapi/deplete.rst | 10 ++ docs/source/usersguide/decay_sources.rst | 211 +++++++++++++++++++---- 2 files changed, 189 insertions(+), 32 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index f112cf8ccff..25fcd898f49 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -287,6 +287,16 @@ the following abstract base classes: abc.SIIntegrator abc.DepSystemSolver +R2S Automation +-------------- + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + R2SManager + D1S Functions ------------- diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst index d5a078135ba..0158a08dc68 100644 --- a/docs/source/usersguide/decay_sources.rst +++ b/docs/source/usersguide/decay_sources.rst @@ -6,42 +6,189 @@ Decay Sources Through the :ref:`depletion ` capabilities in OpenMC, it is possible to simulate radiation emitted from the decay of activated materials. -For fusion energy systems, this is commonly done using what is known as the -`rigorous 2-step `_ (R2S) method. -In this method, a neutron transport calculation is used to determine the neutron -flux and reaction rates over a cell- or mesh-based spatial discretization of the -model. Then, the neutron flux in each discrete region is used to predict the -activated material composition using a depletion solver. Finally, a photon -transport calculation with a source based on the activity and energy spectrum of -the activated materials is used to determine a desired physical response (e.g., -a dose rate) at one or more locations of interest. - -Once a depletion simulation has been completed in OpenMC, the intrinsic decay -source can be determined as follows. First the activated material composition -can be determined using the :class:`openmc.deplete.Results` object. Indexing an -instance of this class with the timestep index returns a -:class:`~openmc.deplete.StepResult` object, which itself has a -:meth:`~openmc.deplete.StepResult.get_material` method. Once the activated -:class:`~openmc.Material` has been obtained, the -:meth:`~openmc.Material.get_decay_photon_energy` method will give the energy -spectrum of the decay photon source. The integral of the spectrum also indicates -the intensity of the source in units of [Bq]. Altogether, the workflow looks as -follows:: - +For fusion energy systems, this is commonly done using either the `rigorous +2-step `_ (R2S) method or the +`direct 1-step `_ (D1S) method. +In the R2S method, a neutron transport calculation is used to determine the +neutron flux and reaction rates over a cell- or mesh-based spatial +discretization of the model. Then, the neutron flux in each discrete region is +used to predict the activated material composition using a depletion solver. +Finally, a photon transport calculation with a source based on the activity and +energy spectrum of the activated materials is used to determine a desired +physical response (e.g., a dose rate) at one or more locations of interest. +OpenMC includes automation for both the R2S and D1S methods as described in the +following sections. + +Rigorous 2-Step (R2S) Calculations +================================== + +OpenMC includes an :class:`openmc.deplete.R2SManager` class that fully automates +cell- and mesh-based R2S calculations. Before we describe this class, it is +useful to understand the basic mechanics of how an R2S calculation works. +Generally, it involves the following steps: + +1. The :meth:`openmc.deplete.get_microxs_and_flux` function is called to run a + neutron transport calculation that determines fluxes and microscopic cross + sections in each activation region. +2. The :class:`openmc.deplete.IndependentOperator` and + :class:`openmc.deplete.PredictorIntegrator` classes are used to carry out a + depletion (activation) calculation in order to determine predicted material + compositions based on a set of timesteps and source rates. +3. The activated material composition is determined using the + :class:`openmc.deplete.Results` class. Indexing an instance of this class + with the timestep index returns a :class:`~openmc.deplete.StepResult` object, + which itself has a :meth:`~openmc.deplete.StepResult.get_material` method + returning an activated material. +4. The :meth:`openmc.Material.get_decay_photon_energy` method is used to obtain + the energy spectrum of the decay photon source. The integral of the spectrum + also indicates the intensity of the source in units of [Bq]. +5. A new photon source is defined using one of OpenMC's source classes with the + energy distribution set equal to the object returned by the + :meth:`openmc.Material.get_decay_photon_energy` method. The source is then + assigned to a photon :class:`~openmc.Model`. +6. A photon transport calculation is run with ``model.run()``. + +Altogether, the workflow looks as follows:: + + # Run neutron transport calculation + fluxes, micros = openmc.deplete.get_microxs_and_flux(model, domains) + + # Run activation calculation + op = openmc.deplete.IndependentOperator(mats, fluxes, micros) + timesteps = ... + source_rates = ... + integrator = openmc.deplete.Integrator(op, timesteps, source_rates) + integrator.integrate() + + # Get decay photon source at last timestep results = openmc.deplete.Results("depletion_results.h5") - - # Get results at last timestep step = results[-1] - - # Get activated material composition for ID=1 activated_mat = step.get_material('1') - - # Determine photon source photon_energy = activated_mat.get_decay_photon_energy() - -By default, the :meth:`~openmc.Material.get_decay_photon_energy` method will -eliminate spectral lines with very low intensity, but this behavior can be -configured with the ``clip_tolerance`` argument. + photon_source = openmc.IndependentSource( + space=..., + energy=energy, + particle='photon', + strength=energy.integral() + ) + + # Run photon transport calculation + model.settings.source = photon_source + model.run() + +Note that by default, the :meth:`~openmc.Material.get_decay_photon_energy` +method will eliminate spectral lines with very low intensity, but this behavior +can be configured with the ``clip_tolerance`` argument. + +Cell-based R2S +-------------- + +In practice, users do not need manually go through each of the steps in an R2S +calculation described above. The :class:`~openmc.deplete.R2SManager` fully +automates the execution of neutron transport, depletion, decay source +generation, and photon transport. For a cell-based R2S calculation, once you +have a :class:`~openmc.Model` that has been defined, simply create an instance +of :class:`~openmc.deplete.R2SManager` by passing the model and a list of cells +to activate:: + + r2s = openmc.deplete.R2SManager(model, [cell1, cell2, cell3]) + +Note that the ``volume`` attribute must be set for any cell that is to be +activated. The :class:`~openmc.deplete.R2SManager` class allows you to +optionally specify a separate photon model; if not given as an argument, it will +create a shallow copy of the original neutron model (available as the +``neutron_model`` attribute) and store it in the ``photon_model`` attribute. We +can use this to define tallies specific to the photon model:: + + dose_tally = openmc.Tally() + ... + r2s.photon_model.tallies = [dose_tally] + +Next, define the timesteps and source rates for the activation calculation:: + + timesteps = [(3.0, 'd'), (5.0, 'h')] + source_rates = [1e12, 0.0] + +In this case, the model is irradiated for 3 days with a source rate of +:math:`10^{12}` neutron/sec and then the source is turned off and the activated +materials are allowed to decay for 5 hours. These parameters should be passed to +the :meth:`~openmc.deplete.R2SManager.run` method to execute the full R2S +calculation. Before we can do that though, for a cell-based calculation, the one +other piece of information that is needed is bounding boxes of the activated +cells:: + + bounding_boxes = { + cell1.id: cell1.bounding_box, + cell2.id: cell2.bounding_box, + cell3.id: cell3.bounding_box + } + +Note that calling the ``bounding_box`` attribute may not work for all +constructive solid geometry regions (for example, a cell that uses a +non-axis-aligned plane). In these cases, the bounding box will need to be +specified manually. Once you have a set of bounding boxes, the R2S calculation +can be run:: + + r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes) + +If not specified otherwise, a photon transport calculation is run at each time +in the depletion schedule. That means in the case above, we would see three +photon transport calculations. To specify specific times at which photon +transport calculations should be run, pass the ``photon_time_indices`` argument. +For example, if we wanted to run a photon transport calculation only on the last +time (after the 5 day decay), we would run:: + + r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes, + photon_time_indices=[2]) + +After an R2S calculation has been run, the :class:`~openmc.deplete.R2SManager` +instance will have a ``results`` dictionary that allows you to directly access +results from each of the steps. It will also write out all the output files into +a directory that is named "r2s_/". The ``output_dir`` argument to the +:meth:`~openmc.deplete.R2SManager.run` method enables you to override the +default output directory name if desired. + +The :meth:`~openmc.deplete.R2SManager.run` method actually runs three +lower-level methods under the hood:: + + r2s.step1_neutron_transport(...) + r2s.step2_activation(...) + r2s.step3_photon_transport(...) + +For users looking for more control over the calculation, these lower-level +methods can be used in lieu of the :meth:`openmc.deplete.R2SManager.run` method. + +Mesh-based R2S +-------------- + +Executing a mesh-based R2S calculation looks nearly identical to the cell-based +R2S workflow described above. The only difference is that instead of passing a +list of cells to the ``domains`` argument of +:class:`~openmc.deplete.R2SManager`, you need to define a mesh object and pass +that instead. This might look like the following:: + + # Define a regular Cartesian mesh + mesh = openmc.RegularMesh() + mesh.lower_left = (-50., -50, 0.) + mesh.upper_right = (50., 50., 75.) + mesh.dimension = (10, 10, 5) + + r2s = openmc.deplete.R2SManager(model, mesh) + +Executing the R2S calculation is then performed by adding photon tallies and +calling the :meth:`~openmc.deplete.R2SManager.run` method with the appropriate +timesteps and source rates. Note that in this case we do not need to define cell +volumes or bounding boxes as is required for a cell-based R2S calculation. +Instead, during the neutron transport step, OpenMC will run a raytracing +calculation to determine material volume fractions within each mesh element +using the :meth:`openmc.MeshBase.material_volumes` method. Arguments to this +method can be customized via the ``mat_vol_kwargs`` argument to the +:meth:`~openmc.deplete.R2SManager.run` method. Most often, this would involve +customizing the number of rays traced to obtain better estimates of volumes. As +an example, if we wanted to run the raytracing calculation with 10 million rays, +we would run:: + + r2s.run(timesteps, source_rates, mat_vol_kwargs={'n_samples': 10_000_000}) Direct 1-Step (D1S) Calculations ================================ From 46ecb2e71c7fdee35efa2519785fe2b169e58e4d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 10 Oct 2025 23:37:12 -0500 Subject: [PATCH 26/29] Apply @eepeterson suggestions from code review Co-authored-by: Ethan Peterson --- docs/source/usersguide/decay_sources.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst index 0158a08dc68..bf84f53b5b5 100644 --- a/docs/source/usersguide/decay_sources.rst +++ b/docs/source/usersguide/decay_sources.rst @@ -67,9 +67,9 @@ Altogether, the workflow looks as follows:: photon_energy = activated_mat.get_decay_photon_energy() photon_source = openmc.IndependentSource( space=..., - energy=energy, + energy=photon_energy, particle='photon', - strength=energy.integral() + strength=photon_energy.integral() ) # Run photon transport calculation From ebe804c5138f5d7c64496655158d9cbe2c8a4518 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 10 Oct 2025 23:58:28 -0500 Subject: [PATCH 27/29] Move _get_all_materials to Model class --- openmc/deplete/r2s.py | 5 ++--- openmc/mesh.py | 19 +------------------ openmc/model/model.py | 23 +++++++++++++++++++++++ 3 files changed, 26 insertions(+), 21 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 8477fa3f961..4f67ff0b52f 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -11,7 +11,6 @@ from .microxs import get_microxs_and_flux, write_microxs_hdf5, read_microxs_hdf5 from .results import Results from ..checkvalue import PathLike -from ..mesh import _get_all_materials def get_activation_materials( @@ -47,7 +46,7 @@ def get_activation_materials( elems, _ = np.where(mmv._materials > -1) # Get all materials in the model - material_dict = _get_all_materials(model) + material_dict = model._get_all_materials() # Create a new activation material for each element-material combination materials = openmc.Materials() @@ -552,7 +551,7 @@ def get_decay_photon_source_mesh( instances for the decay photons in the corresponding mesh element. """ - mat_dict = _get_all_materials(self.neutron_model) + mat_dict = self.neutron_model._get_all_materials() # Some MeshSource objects will have empty positions; create a "null source" # that is used for this case diff --git a/openmc/mesh.py b/openmc/mesh.py index 0f33bcc9dc6..d34691bc7bd 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -318,7 +318,7 @@ def get_homogenized_materials( mat_volume_by_element = [vols.by_element(i) for i in range(vols.num_elements)] # Get dictionary of all materials - materials = _get_all_materials(model) + materials = model._get_all_materials() # Create homogenized material for each element homogenized_materials = [] @@ -2859,23 +2859,6 @@ def _read_meshes(elem): return out -# TODO: This should really get incorporated in lower-level calls to -# get_all_materials, but right now it requires information from the Model object -def _get_all_materials(model: openmc.Model) -> dict[int, openmc.Material]: - # Get all materials from the Geometry object - materials = model.geometry.get_all_materials() - - # Account for materials in DAGMC universes - for cell in model.geometry.get_all_cells().values(): - if isinstance(cell.fill, openmc.DAGMCUniverse): - names = cell.fill.material_names - materials.update({ - mat.id: mat for mat in model.materials if mat.name in names - }) - - return materials - - # hexahedron element connectivity # lower-k connectivity offsets _HEX_VERTEX_CONN = ((0, 0, 0), diff --git a/openmc/model/model.py b/openmc/model/model.py index 59e6fa511e5..66a0253341c 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -203,6 +203,29 @@ def _materials_by_name(self) -> dict[int, openmc.Material]: result[mat.name].add(mat) return result + # TODO: This should really get incorporated in lower-level calls to + # get_all_materials, but right now it requires information from the Model object + def _get_all_materials(self) -> dict[int, openmc.Material]: + """Get all materials including those in DAGMC universes + + Returns + ------- + dict + Dictionary mapping material ID to material instances + """ + # Get all materials from the Geometry object + materials = self.geometry.get_all_materials() + + # Account for materials in DAGMC universes + for cell in self.geometry.get_all_cells().values(): + if isinstance(cell.fill, openmc.DAGMCUniverse): + names = cell.fill.material_names + materials.update({ + mat.id: mat for mat in self.materials if mat.name in names + }) + + return materials + @classmethod def from_xml( cls, From 4badb20c9d7535553e1c88b91ac2fcfd567d0813 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 20 Oct 2025 09:38:08 -0500 Subject: [PATCH 28/29] Don't use constraints on a null photon distribution --- openmc/deplete/r2s.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/openmc/deplete/r2s.py b/openmc/deplete/r2s.py index 4f67ff0b52f..7b5deddc23c 100644 --- a/openmc/deplete/r2s.py +++ b/openmc/deplete/r2s.py @@ -601,13 +601,14 @@ def get_decay_photon_source_mesh( # Create decay photon source source energy = activated_mat.get_decay_photon_energy() - strength = energy.integral() if energy is not None else 0.0 - source_lists[j][index_elem] = openmc.IndependentSource( - energy=energy, - particle='photon', - strength=strength, - constraints={'domains': [mat_dict[mat_id]]} - ) + if energy is not None: + strength = energy.integral() + source_lists[j][index_elem] = openmc.IndependentSource( + energy=energy, + particle='photon', + strength=strength, + constraints={'domains': [mat_dict[mat_id]]} + ) # Increment index of activated material index_mat += 1 From f9160bf36771e1cd70cdd02ac43c2f517e762185 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 31 Oct 2025 16:16:55 -0500 Subject: [PATCH 29/29] Apply @shimwell suggestions from code review Co-authored-by: Jonathan Shimwell --- docs/source/usersguide/decay_sources.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/source/usersguide/decay_sources.rst b/docs/source/usersguide/decay_sources.rst index bf84f53b5b5..398680e7464 100644 --- a/docs/source/usersguide/decay_sources.rst +++ b/docs/source/usersguide/decay_sources.rst @@ -83,7 +83,7 @@ can be configured with the ``clip_tolerance`` argument. Cell-based R2S -------------- -In practice, users do not need manually go through each of the steps in an R2S +In practice, users do not need to manually go through each of the steps in an R2S calculation described above. The :class:`~openmc.deplete.R2SManager` fully automates the execution of neutron transport, depletion, decay source generation, and photon transport. For a cell-based R2S calculation, once you @@ -136,7 +136,7 @@ in the depletion schedule. That means in the case above, we would see three photon transport calculations. To specify specific times at which photon transport calculations should be run, pass the ``photon_time_indices`` argument. For example, if we wanted to run a photon transport calculation only on the last -time (after the 5 day decay), we would run:: +time (after the 5 hour decay), we would run:: r2s.run(timesteps, source_rates, bounding_boxes=bounding_boxes, photon_time_indices=[2]) @@ -169,7 +169,7 @@ that instead. This might look like the following:: # Define a regular Cartesian mesh mesh = openmc.RegularMesh() - mesh.lower_left = (-50., -50, 0.) + mesh.lower_left = (-50., -50., 0.) mesh.upper_right = (50., 50., 75.) mesh.dimension = (10, 10, 5)