Adding material depletion function (#3420)

Co-authored-by: Jon Shimwell <jon@proximafusion.com>
Co-authored-by: Micah Gale <mgale@fastmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
Jonathan Shimwell 2025-07-18 16:37:57 +02:00 committed by GitHub
parent 659e43af7d
commit 8be65513b7
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
4 changed files with 281 additions and 1 deletions

View file

@ -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)

View file

@ -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

View file

@ -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()

View file

@ -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