From 8be65513b77c08e79e1dc5f4ef921d411ff990a9 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Fri, 18 Jul 2025 16:37:57 +0200 Subject: [PATCH] Adding material depletion function (#3420) Co-authored-by: Jon Shimwell Co-authored-by: Micah Gale Co-authored-by: Paul Romano --- openmc/deplete/independent_operator.py | 2 +- openmc/material.py | 174 +++++++++++++++++++++++++ tests/unit_tests/test_material.py | 43 ++++++ tests/unit_tests/test_materials.py | 63 +++++++++ 4 files changed, 281 insertions(+), 1 deletion(-) create mode 100644 tests/unit_tests/test_materials.py diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index 3eb50b05ec..c192907cf2 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -273,7 +273,7 @@ class IndependentOperator(OpenMCOperator): Returns ------- nuclides : set of str - Set of nuclide names that have cross secton data + Set of nuclide names that have cross section data """ return set(cross_sections[0].nuclides) diff --git a/openmc/material.py b/openmc/material.py index c6765178e9..cec74af627 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -6,6 +6,8 @@ from numbers import Real from pathlib import Path import re import sys +import tempfile +from typing import Sequence, Dict import warnings import lxml.etree as ET @@ -1717,6 +1719,67 @@ class Material(IDManagerMixin): return mat + def deplete( + self, + multigroup_flux: Sequence[float], + energy_group_structure: Sequence[float] | str, + timesteps: Sequence[float] | Sequence[tuple[float, str]], + source_rates: float | Sequence[float], + timestep_units: str = 's', + chain_file: cv.PathLike | "openmc.deplete.Chain" | None = None, + reactions: Sequence[str] | None = None, + ) -> list[openmc.Material]: + """Depletes that material, evolving the nuclide densities + + .. versionadded:: 0.15.3 + + Parameters + ---------- + multigroup_flux: Sequence[float] + Energy-dependent multigroup flux values, where each sublist corresponds + to a specific material. Will be normalized so that it sums to 1. + energy_group_structure : Sequence[float] | str + Energy group boundaries in [eV] or the name of the group structure. + timesteps : iterable of float or iterable of tuple + Array 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 iterable of float, optional + Source rate in [neutron/sec] or neutron flux in [neutron/s-cm^2] for + each interval in :attr:`timesteps` + timestep_units : {'s', 'min', 'h', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'min' means minutes, 'h' means hours, 'a' means Julian years + and 'MWd/kg' indicates that the values are given in burnup (MW-d of + energy deposited per kilogram of initial heavy metal). + chain_file : PathLike or Chain + Path to the depletion chain XML file or instance of openmc.deplete.Chain. + Defaults to ``openmc.config['chain_file']``. + reactions : list of str, optional + Reactions to get cross sections for. If not specified, all neutron + reactions listed in the depletion chain file are used. + + Returns + ------- + list of openmc.Material, one for each timestep + + """ + + materials = openmc.Materials([self]) + + depleted_materials_dict = materials.deplete( + multigroup_fluxes=[multigroup_flux], + energy_group_structures=[energy_group_structure], + timesteps=timesteps, + source_rates=source_rates, + timestep_units=timestep_units, + chain_file=chain_file, + reactions=reactions, + ) + + return depleted_materials_dict[self.id] + def mean_free_path(self, energy: float) -> float: """Calculate the mean free path of neutrons in the material at a given @@ -1947,3 +2010,114 @@ class Materials(cv.CheckedList): root = tree.getroot() return cls.from_xml_element(root) + + + def deplete( + self, + multigroup_fluxes: Sequence[Sequence[float]], + energy_group_structures: Sequence[Sequence[float] | str], + timesteps: Sequence[float] | Sequence[tuple[float, str]], + source_rates: float | Sequence[float], + timestep_units: str = 's', + chain_file: cv.PathLike | "openmc.deplete.Chain" | None = None, + reactions: Sequence[str] | None = None, + ) -> Dict[int, list[openmc.Material]]: + """Depletes that material, evolving the nuclide densities + + .. versionadded:: 0.15.3 + + Parameters + ---------- + multigroup_fluxes: Sequence[Sequence[float]] + Energy-dependent multigroup flux values, where each sublist corresponds + to a specific material. Will be normalized so that it sums to 1. + energy_group_structures': Sequence[Sequence[float] | str] + Energy group boundaries in [eV] or the name of the group structure. + timesteps : iterable of float or iterable of tuple + Array 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 iterable of float, optional + Source rate in [neutron/sec] or neutron flux in [neutron/s-cm^2] for + each interval in :attr:`timesteps` + timestep_units : {'s', 'min', 'h', 'd', 'a', 'MWd/kg'} + Units for values specified in the `timesteps` argument. 's' means + seconds, 'min' means minutes, 'h' means hours, 'a' means Julian years + and 'MWd/kg' indicates that the values are given in burnup (MW-d of + energy deposited per kilogram of initial heavy metal). + chain_file : PathLike or Chain + Path to the depletion chain XML file or instance of openmc.deplete.Chain. + Defaults to ``openmc.config['chain_file']``. + reactions : list of str, optional + Reactions to get cross sections for. If not specified, all neutron + reactions listed in the depletion chain file are used. + + Returns + ------- + list of openmc.Material, one for each timestep + + """ + + import openmc.deplete + from .deplete.chain import _get_chain + + # setting all materials to be depletable + for mat in self: + mat.depletable = True + + chain = _get_chain(chain_file) + + # Create MicroXS objects for all materials + micros = [] + fluxes = [] + + with openmc.lib.TemporarySession(): + for material, flux, energy in zip( + self, multigroup_fluxes, energy_group_structures + ): + temperature = material.temperature or 293.6 + micro_xs = openmc.deplete.MicroXS.from_multigroup_flux( + energies=energy, + multigroup_flux=flux, + chain_file=chain, + temperature=temperature, + reactions=reactions, + ) + micros.append(micro_xs) + fluxes.append(material.volume) + + # Create a single operator for all materials + operator = openmc.deplete.IndependentOperator( + materials=self, + fluxes=fluxes, + micros=micros, + normalization_mode="source-rate", + chain_file=chain, + ) + + integrator = openmc.deplete.PredictorIntegrator( + operator=operator, + timesteps=timesteps, + source_rates=source_rates, + timestep_units=timestep_units, + ) + + with tempfile.TemporaryDirectory() as tmpdir: + # Run integrator + results_path = Path(tmpdir) / "depletion_results.h5" + integrator.integrate(path=results_path) + + # Load depletion results + results = openmc.deplete.Results(results_path) + + # For each material, get activated composition at each timestep + all_depleted_materials = { + material.id: [ + result.get_material(str(material.id)) + for result in results + ] + for material in self + } + + return all_depleted_materials diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index db5f4ce32d..eae814a755 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -3,8 +3,11 @@ from pathlib import Path import pytest +import numpy as np + import openmc from openmc.data import decay_photon_energy +from openmc.deplete import Chain import openmc.examples import openmc.model import openmc.stats @@ -712,6 +715,46 @@ def test_avoid_subnormal(run_in_tmpdir): assert mats[0].get_nuclide_atom_densities()['H2'] == 0.0 +def test_material_deplete(): + pristine_material = openmc.Material() + pristine_material.add_nuclide("Ni58", 1.0) + pristine_material.set_density("g/cm3", 7.87) + pristine_material.depletable = True + pristine_material.temperature = 293.6 + pristine_material.volume = 1. + + mg_flux = [0.5e11] * 42 + + chain = Chain.from_xml( + Path(__file__).parents[1] / "chain_ni.xml" + ) + + depleted_material = pristine_material.deplete( + multigroup_flux=mg_flux, + energy_group_structure="VITAMIN-J-42", + timesteps=[10, 70.86], + source_rates=[1e19, 0.0], + timestep_units="d", + chain_file=chain, + ) + + for material in depleted_material: + assert isinstance(material, openmc.Material) + assert len(material.get_nuclides()) > len(pristine_material.get_nuclides()) + + Co58_mat_1_step_0 = depleted_material[0].get_nuclide_atom_densities("Co58")["Co58"] + Co58_mat_1_step_1 = depleted_material[1].get_nuclide_atom_densities("Co58")["Co58"] + Co58_mat_1_step_2 = depleted_material[2].get_nuclide_atom_densities("Co58")["Co58"] + + assert Co58_mat_1_step_0 == 0.0 + + # Check that Co58 is produced in the first step + assert Co58_mat_1_step_1 > 0.0 + + # Check that Co58 is halved in the second step which is one halflife later + assert np.allclose(Co58_mat_1_step_1 * 0.5, Co58_mat_1_step_2) + + def test_mean_free_path(): mat1 = openmc.Material() diff --git a/tests/unit_tests/test_materials.py b/tests/unit_tests/test_materials.py new file mode 100644 index 0000000000..6c8c59fc92 --- /dev/null +++ b/tests/unit_tests/test_materials.py @@ -0,0 +1,63 @@ +from pathlib import Path + +import openmc +from openmc.deplete import Chain + + +def test_materials_deplete(): + pristine_material_1 = openmc.Material() + pristine_material_1.add_nuclide("Ni58", 1.) + pristine_material_1.set_density("g/cm3", 7.87) + pristine_material_1.depletable = True + pristine_material_1.temperature = 293.6 + pristine_material_1.volume = 1. + + pristine_material_2 = openmc.Material() + pristine_material_2.add_nuclide("Ni60", 1.) + pristine_material_2.set_density("g/cm3", 7.87) + pristine_material_2.depletable = True + pristine_material_2.temperature = 293.6 + pristine_material_2.volume = 1. + + pristine_materials = openmc.Materials([pristine_material_1, pristine_material_2]) + + mg_flux = [0.5e11] * 42 + + chain = Chain.from_xml( + Path(__file__).parents[1] / "chain_ni.xml" + ) + + depleted_material = pristine_materials.deplete( + multigroup_fluxes=[mg_flux, mg_flux], + energy_group_structures=["VITAMIN-J-42", "VITAMIN-J-42"], + timesteps=[100, 100], + source_rates=[1e19, 0.0], + timestep_units="d", + chain_file=chain, + ) + + assert list(depleted_material.keys()) == [pristine_material_1.id, pristine_material_2.id] + for mat_id, materials in depleted_material.items(): + for material in materials: + assert isinstance(material, openmc.Material) + assert len(material.get_nuclides()) > 1 + assert mat_id == material.id + + mats = depleted_material[pristine_material_1.id] + Co58_mat_1_step_0 = mats[0].get_nuclide_atom_densities("Co58")["Co58"] + Co58_mat_1_step_1 = mats[1].get_nuclide_atom_densities("Co58")["Co58"] + Co58_mat_1_step_2 = mats[2].get_nuclide_atom_densities("Co58")["Co58"] + + assert Co58_mat_1_step_0 == 0.0 + # Co58 is the main activation product of Ni58 in the first irradiation step. + # It then decays in the second cooling step (flux = 0) + assert Co58_mat_1_step_1 > 0.0 and Co58_mat_1_step_1 > Co58_mat_1_step_2 + + Ni59_mat_1_step_0 = mats[0].get_nuclide_atom_densities("Ni59")["Ni59"] + Ni59_mat_1_step_1 = mats[1].get_nuclide_atom_densities("Ni59")["Ni59"] + Ni59_mat_1_step_2 = mats[2].get_nuclide_atom_densities("Ni59")["Ni59"] + + assert Ni59_mat_1_step_0 == 0.0 + # Ni59 is one of the main activation product of Ni60 in the first irradiation + # step. It then decays in the second cooling step (flux = 0) + assert Ni59_mat_1_step_1 > 0.0 and Ni59_mat_1_step_1 > Ni59_mat_1_step_2