From bbd77569309e032b1577067d41e2c38c21b69e6a Mon Sep 17 00:00:00 2001 From: Olek <45364492+yardasol@users.noreply.github.com> Date: Wed, 13 Sep 2023 20:00:52 -0500 Subject: [PATCH] Allow decay only nuclides to be specified in `Material` for a depletion simulation (#2616) Co-authored-by: Jonathan Shimwell Co-authored-by: Paul Romano --- openmc/deplete/coupled_operator.py | 2 +- openmc/deplete/openmc_operator.py | 30 +++-- openmc/material.py | 40 +++++-- tests/chain_simple_decay.xml | 66 +++++++++++ .../deplete_decay_only/__init__.py | 0 .../deplete_decay_only/test.py | 103 ++++++++++++++++++ tests/unit_tests/test_material.py | 12 ++ 7 files changed, 233 insertions(+), 20 deletions(-) create mode 100644 tests/chain_simple_decay.xml create mode 100644 tests/regression_tests/deplete_decay_only/__init__.py create mode 100644 tests/regression_tests/deplete_decay_only/test.py diff --git a/openmc/deplete/coupled_operator.py b/openmc/deplete/coupled_operator.py index ab7314ac786..54f8f9462bf 100644 --- a/openmc/deplete/coupled_operator.py +++ b/openmc/deplete/coupled_operator.py @@ -416,7 +416,7 @@ def _generate_materials_xml(self): for mat in self.materials: mat._nuclides.sort(key=lambda x: nuclides.index(x[0])) - self.materials.export_to_xml() + self.materials.export_to_xml(nuclides_to_ignore=self._decay_nucs) def __call__(self, vec, source_rate): """Runs a simulation. diff --git a/openmc/deplete/openmc_operator.py b/openmc/deplete/openmc_operator.py index 319d9870122..30203d23fb5 100644 --- a/openmc/deplete/openmc_operator.py +++ b/openmc/deplete/openmc_operator.py @@ -6,6 +6,7 @@ """ from abc import abstractmethod +from warnings import warn from typing import List, Tuple, Dict import numpy as np @@ -133,6 +134,18 @@ def __init__( # This nuclides variables contains every nuclides # for which there is an entry in the micro_xs parameter openmc.reset_auto_ids() + + self.nuclides_with_data = self._get_nuclides_with_data( + self.cross_sections) + + # Select nuclides with data that are also in the chain + self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides + if nuc.name in self.nuclides_with_data] + + # Select nuclides without data that are also in the chain + self._decay_nucs = [nuc.name for nuc in self.chain.nuclides + if nuc.name not in self.nuclides_with_data] + self.burnable_mats, volumes, all_nuclides = self._get_burnable_mats() self.local_mats = _distribute(self.burnable_mats) @@ -142,13 +155,6 @@ def __init__( if self.prev_res is not None: self._load_previous_results() - self.nuclides_with_data = self._get_nuclides_with_data( - self.cross_sections) - - # Select nuclides with data that are also in the chain - self._burnable_nucs = [nuc.name for nuc in self.chain.nuclides - if nuc.name in self.nuclides_with_data] - # Extract number densities from the geometry / previous depletion run self._extract_number(self.local_mats, volumes, @@ -171,7 +177,7 @@ def _get_burnable_mats(self) -> Tuple[List[str], Dict[str, float], List[str]]: Returns ------- burnable_mats : list of str - List of burnable material IDs + list of burnable material IDs volume : dict of str to float Volume of each material in [cm^3] nuclides : list of str @@ -188,7 +194,13 @@ def _get_burnable_mats(self) -> Tuple[List[str], Dict[str, float], List[str]]: # Iterate once through the geometry to get dictionaries for mat in self.materials: for nuclide in mat.get_nuclides(): - model_nuclides.add(nuclide) + if nuclide in self.nuclides_with_data or self._decay_nucs: + model_nuclides.add(nuclide) + else: + msg = (f"Nuclilde {nuclide} in material {mat.id} is not " + "present in the depletion chain and has no cross " + "section data.") + raise warn(msg) if mat.depletable: burnable_mats.add(str(mat.id)) if mat.volume is None: diff --git a/openmc/material.py b/openmc/material.py index f0e9117e63d..4aa5c9f0ffd 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -1250,7 +1250,7 @@ def clone(self, memo: Optional[dict] = None) -> Material: return memo[self] - def _get_nuclide_xml(self, nuclide: str) -> ET.Element: + def _get_nuclide_xml(self, nuclide: NuclideTuple) -> ET.Element: xml_element = ET.Element("nuclide") xml_element.set("name", nuclide.name) @@ -1267,15 +1267,28 @@ def _get_macroscopic_xml(self, macroscopic: str) -> ET.Element: return xml_element - def _get_nuclides_xml(self, nuclides: typing.Iterable[str]) -> List[ET.Element]: + def _get_nuclides_xml( + self, nuclides: typing.Iterable[NuclideTuple], + nuclides_to_ignore: Optional[typing.Iterable[str]] = None)-> List[ET.Element]: xml_elements = [] - for nuclide in nuclides: - xml_elements.append(self._get_nuclide_xml(nuclide)) + + # Remove any nuclides to ignore from the XML export + if nuclides_to_ignore: + nuclides = [nuclide for nuclide in nuclides if nuclide.name not in nuclides_to_ignore] + + xml_elements = [self._get_nuclide_xml(nuclide) for nuclide in nuclides] + return xml_elements - def to_xml_element(self) -> ET.Element: + def to_xml_element( + self, nuclides_to_ignore: Optional[typing.Iterable[str]] = None) -> ET.Element: """Return XML representation of the material + Parameters + ---------- + nuclides_to_ignore : list of str + Nuclides to ignore when exporting to XML. + Returns ------- element : lxml.etree._Element @@ -1320,7 +1333,8 @@ def to_xml_element(self) -> ET.Element: if self._macroscopic is None: # Create nuclide XML subelements - subelements = self._get_nuclides_xml(self._nuclides) + subelements = self._get_nuclides_xml(self._nuclides, + nuclides_to_ignore=nuclides_to_ignore) for subelement in subelements: element.append(subelement) else: @@ -1576,7 +1590,8 @@ def make_isotropic_in_lab(self): for material in self: material.make_isotropic_in_lab() - def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_indent=True): + def _write_xml(self, file, header=True, level=0, spaces_per_level=2, + trailing_indent=True, nuclides_to_ignore=None): """Writes XML content of the materials to an open file handle. Parameters @@ -1591,6 +1606,8 @@ def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_in Number of spaces per indentation trailing_indentation : bool Whether or not to write a trailing indentation for the materials element + nuclides_to_ignore : list of str + Nuclides to ignore when exporting to XML. """ indentation = level*spaces_per_level*' ' @@ -1611,7 +1628,7 @@ def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_in # Write the elements. for material in sorted(self, key=lambda x: x.id): - element = material.to_xml_element() + element = material.to_xml_element(nuclides_to_ignore=nuclides_to_ignore) clean_indentation(element, level=level+1) element.tail = element.tail.strip(' ') file.write((level+1)*spaces_per_level*' ') @@ -1626,13 +1643,16 @@ def _write_xml(self, file, header=True, level=0, spaces_per_level=2, trailing_in if trailing_indent: file.write(indentation) - def export_to_xml(self, path: PathLike = 'materials.xml'): + def export_to_xml(self, path: PathLike = 'materials.xml', + nuclides_to_ignore: Optional[typing.Iterable[str]] = None): """Export material collection to an XML file. Parameters ---------- path : str Path to file to write. Defaults to 'materials.xml'. + nuclides_to_ignore : list of str + Nuclides to ignore when exporting to XML. """ # Check if path is a directory @@ -1645,7 +1665,7 @@ def export_to_xml(self, path: PathLike = 'materials.xml'): # one go. with open(str(p), 'w', encoding='utf-8', errors='xmlcharrefreplace') as fh: - self._write_xml(fh) + self._write_xml(fh, nuclides_to_ignore=nuclides_to_ignore) @classmethod def from_xml_element(cls, elem) -> Material: diff --git a/tests/chain_simple_decay.xml b/tests/chain_simple_decay.xml new file mode 100644 index 00000000000..1c8c5e4eecc --- /dev/null +++ b/tests/chain_simple_decay.xml @@ -0,0 +1,66 @@ + + + + + + + 3696.125 4095.822 4477.27 5097.122 29452.1 29781.3 33566.5 33629.4 33865.1 33878.5 34395.3 34408.0 34486.2 34488.2 112780.0 113150.0 162650.0 165740.0 184490.0 197190.0 220502.0 229720.0 247500.0 254740.0 264260.0 288451.0 290270.0 304910.0 305830.0 326000.0 333600.0 342520.0 361850.0 403030.0 414830.0 417633.0 429930.0 433741.0 451630.0 530800.0 546557.0 575970.0 588280.0 616900.0 649850.0 656090.0 679220.0 684600.0 690130.0 707920.0 785480.0 795500.0 797710.0 807200.0 836804.0 960290.0 961430.0 971960.0 972620.0 995090.0 1038760.0 1096860.0 1101580.0 1124000.0 1131511.0 1151510.0 1159900.0 1169040.0 1225600.0 1240470.0 1254800.0 1260409.0 1315770.0 1334800.0 1343660.0 1367890.0 1441800.0 1448350.0 1457560.0 1502790.0 1521990.0 1543700.0 1566410.0 1678027.0 1706459.0 1791196.0 1830690.0 1845300.0 1927300.0 1948490.0 2045880.0 2112400.0 2151500.0 2189400.0 2255457.0 2408650.0 2466070.0 2477100.0 9.714352819815078e-10 7.460941651551526e-09 6.047882056201745e-09 8.510389107205747e-10 3.7979729633684727e-08 7.033747061198817e-08 6.602815946851458e-09 1.2800909198245071e-08 6.513060244589454e-11 8.894667887850525e-11 1.3843783302275766e-09 2.700005498135763e-09 1.2141115256573815e-11 1.654893875618862e-11 3.7007705885806645e-09 2.0186021392258173e-09 2.859686363903241e-09 9.16781804898392e-09 6.896890642354875e-09 9.588360161322631e-09 5.130613770532286e-07 7.06510748729036e-08 8.410842246774237e-09 6.7286737974193905e-09 5.3829390379355124e-08 9.083709626516178e-07 8.915492781580693e-08 9.251926471451662e-09 2.783988783682273e-08 6.728673797419391e-10 1.093409492080651e-08 2.5232526740322717e-10 5.467047460403255e-08 6.812782219887132e-08 8.83138435911295e-08 1.0345335963532313e-06 8.915492781580693e-08 1.623292553627428e-07 9.251926471451663e-08 9.251926471451662e-09 2.0942997194467853e-06 3.784879011048407e-08 1.513951604419363e-08 1.093409492080651e-08 1.337323917237104e-07 2.186818984161302e-08 1.5980600268871052e-08 6.7286737974193905e-09 3.784879011048407e-08 1.9344937167580748e-07 4.457746390790346e-08 6.7286737974193905e-09 5.0465053480645434e-08 1.3457347594838781e-08 1.9597262434983976e-06 1.0093010696129087e-08 4.289529545854862e-08 2.607361096500014e-07 3.53255374364518e-07 4.541854813258089e-08 2.329803302356464e-06 2.607361096500014e-08 4.7100716581935733e-07 1.059766123093554e-06 6.6193328482113254e-06 4.205421123387119e-10 3.027903208838726e-08 2.565306885266143e-07 1.2616263370161358e-08 2.649415307733885e-07 3.3643368987096953e-09 8.410842246774237e-06 1.934493716758075e-08 9.251926471451662e-09 2.2709274066290446e-08 1.7830985563161385e-07 5.046505348064544e-09 9.251926471451663e-08 2.5400743585258203e-06 3.154065842540339e-07 1.093409492080651e-08 7.569758022096815e-09 3.784879011048407e-07 2.8008104681758214e-06 1.2027504412887161e-06 2.26251656438227e-06 1.6989901338483962e-07 1.6821684493548476e-09 8.663167514177466e-08 1.8503852942903323e-08 2.5568960430193683e-07 2.0186021392258175e-08 6.560456952483906e-09 3.784879011048407e-09 1.7999202408096872e-07 2.800810468175822e-07 2.1027105616935597e-08 4.205421123387119e-10 + + + + + + 0.0 10000.0 20000.0 50000.0 100000.0 200000.0 300000.0 400000.0 600000.0 800000.0 1000000.0 1220000.0 1.1612176249914943e-11 0.0 3.1976524142990123e-11 0.0 6.418466676129901e-13 1.894551923065874e-10 5.099949011192198e-13 3.170644340540729e-13 2.756798470867507e-12 6.67627508515641e-14 3.711746246444946e-15 0.0 + + + + + + + + + + + + + + + + + + + + 2.53000e-02 + + Gd157 Gd156 I135 Xe135 Xe136 Cs135 + 1.093250e-04 2.087260e-04 2.780820e-02 6.759540e-03 2.392300e-02 4.356330e-05 + + + + + + + 3065.349 12960.11 13197.49 16125.43 19185.05 19590.0 31600.0 34700.0 41400.0 41960.0 51220.0 54100.0 54250.0 64350.0 72700.0 75020.0 76198.0 90330.0 93795.0 96090.0 105278.0 106074.0 106608.0 106771.0 108948.0 109154.0 109160.0 109395.0 109433.0 115450.0 120350.0 136550.0 140760.0 142400.0 143760.0 150930.0 163330.0 173300.0 182610.0 185715.0 194940.0 198900.0 202110.0 205311.0 215280.0 221380.0 228780.0 233500.0 240870.0 246840.0 266450.0 275129.0 275430.0 281420.0 282920.0 289560.0 291650.0 301700.0 317100.0 343500.0 345900.0 356030.0 387820.0 410290.0 428710.0 448400.0 1.4211820389290126e-18 3.695705148646686e-18 4.752472005970374e-19 4.0825252294602915e-18 8.998127221991115e-19 1.9035285506819052e-21 5.305446177665699e-21 1.1547147563154756e-20 9.362552078233586e-21 1.86835843588509e-20 1.0610892355331398e-20 2.736290732002608e-22 4.776811520523089e-21 3.977552295559136e-21 3.432935762018982e-20 1.872510415646717e-20 2.4966805541956234e-21 1.0265454754703584e-18 1.7575878976520235e-18 2.839974130397521e-20 2.1057468800839102e-19 4.1342252420362975e-19 7.098565272539688e-21 8.20050332279016e-21 5.276915360632628e-20 1.0602918581811435e-19 4.806110066826575e-19 2.0521431485853304e-21 2.375669880669523e-21 9.362552078233586e-21 8.267152210184412e-21 3.745020831293435e-21 6.865871524037964e-20 1.5604253463722645e-21 3.420452359248004e-18 2.4966805541956232e-20 1.5853921519142206e-18 1.8725104156467174e-21 1.0610892355331397e-19 1.7851265962498703e-17 1.966135936429053e-19 1.3107572909527022e-20 3.370518748164091e-19 1.5635461970650089e-18 9.050467008959134e-21 3.745020831293434e-20 2.18459548492117e-21 9.050467008959134e-21 2.3406380195583967e-20 1.6540508671546e-20 1.8725104156467174e-21 1.622842360227155e-20 2.18459548492117e-21 1.8725104156467174e-21 1.8725104156467174e-21 2.18459548492117e-21 1.2483402770978116e-20 1.5604253463722645e-21 3.120850692744529e-22 9.362552078233587e-22 1.2483402770978116e-20 1.5604253463722645e-21 1.2483402770978116e-20 9.362552078233587e-22 3.120850692744529e-22 3.120850692744529e-22 + + + 2.53000e-02 + + Gd157 Gd156 I135 Xe135 Xe136 Cs135 + 6.142710e-5 1.483250e-04 0.0292737 0.002566345 0.0219242 4.9097e-6 + + + + + + + 3061.32 12959.8 13440.07 16150.05 19148.83 49550.0 90330.0 93795.0 105278.0 106074.0 106608.0 106771.0 108948.0 109154.0 109395.0 109433.0 113500.0 5.3130501476983077e-20 1.4936931167066328e-19 1.943509007980312e-20 1.8224649880365157e-19 4.0452400597373603e-20 3.146222282132249e-21 3.300026341588747e-23 5.444173877278485e-23 6.3486292118621375e-24 1.2400176384733938e-23 2.124537722121886e-25 2.4542156071495764e-25 1.5792541400719878e-24 3.1731991718126067e-24 6.141568457919309e-26 7.109808533300877e-26 5.014291762148272e-22 + + + 2.53000e-02 + + Gd157 Gd156 I135 Xe135 Xe136 Cs135 + 4.141120e-04 7.605360e-04 0.0135457 0.00026864 0.0024432 3.7100E-07 + + + + diff --git a/tests/regression_tests/deplete_decay_only/__init__.py b/tests/regression_tests/deplete_decay_only/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/deplete_decay_only/test.py b/tests/regression_tests/deplete_decay_only/test.py new file mode 100644 index 00000000000..cb716a5ad2d --- /dev/null +++ b/tests/regression_tests/deplete_decay_only/test.py @@ -0,0 +1,103 @@ +""" Transport-free depletion test suite """ + +from pathlib import Path +import shutil + +import numpy as np +import pytest +import openmc +import openmc.deplete +from openmc.deplete import CoupledOperator, IndependentOperator, MicroXS + + +@pytest.fixture(scope="module") +def model(): + fuel = openmc.Material(name="uo2") + fuel.add_element("U", 1, percent_type="ao", enrichment=4.25) + fuel.add_element("O", 2) + fuel.add_nuclide("Xe135_m1", 1) + fuel.add_nuclide("Cs135_m1", 1) + fuel.set_density("g/cc", 10.4) + + clad = openmc.Material(name="clad") + clad.add_element("Zr", 1) + clad.set_density("g/cc", 6) + + water = openmc.Material(name="water") + water.add_element("O", 1) + water.add_element("H", 2) + water.set_density("g/cc", 1.0) + water.add_s_alpha_beta("c_H_in_H2O") + + radii = [0.42, 0.45] + fuel.volume = np.pi * radii[0] ** 2 + + materials = openmc.Materials([fuel, clad, water]) + + pin_surfaces = [openmc.ZCylinder(r=r) for r in radii] + pin_univ = openmc.model.pin(pin_surfaces, materials) + bound_box = openmc.rectangular_prism(1.24, 1.24, boundary_type="reflective") + root_cell = openmc.Cell(fill=pin_univ, region=bound_box) + geometry = openmc.Geometry([root_cell]) + + settings = openmc.Settings() + settings.particles = 1000 + settings.inactive = 5 + settings.batches = 10 + + return openmc.Model(geometry, materials, settings) + +@pytest.fixture(scope="module") +def micro_xs(): + micro_xs_file = Path(__file__).parents[2] / 'micro_xs_simple.csv' + return MicroXS.from_csv(micro_xs_file) + + +@pytest.fixture(scope="module") +def chain_file(): + return Path(__file__).parents[2] / 'chain_simple_decay.xml' + + +@pytest.mark.parametrize("operator_type", ["coupled", "independent"]) +def test_decay_only(run_in_tmpdir, operator_type, model, micro_xs, chain_file): + """Transport free system test suite. + + """ + # Create operator + if operator_type == "coupled": + op = CoupledOperator(model, chain_file=chain_file) + else: + op = IndependentOperator(openmc.Materials([model.materials[0]]), + [1e15], + [micro_xs], + chain_file) + + # Power and timesteps + dt = [917.4, 2262.6] # one Xe135_m1 half life and one Cs135_m1 half life + + # Perform simulation using the predictor algorithm + openmc.deplete.PredictorIntegrator(op, + dt, + power=0.0, + timestep_units='s').integrate() + + # Get path to test and reference results + path_test = op.output_dir / 'depletion_results.h5' + + # Load the reference/test results + res_test = openmc.deplete.Results(path_test) + + _, xe135m1_atoms = res_test.get_atoms('1', 'Xe135_m1') + _, xe135_atoms = res_test.get_atoms('1', 'Xe135') + _, cs135m1_atoms = res_test.get_atoms('1', 'Cs135_m1') + _, cs135_atoms = res_test.get_atoms('1', 'Cs135') + + tol = 1.0e-14 + assert xe135m1_atoms[0] == pytest.approx(xe135m1_atoms[1] * 2, rel=tol) + + # WARNING: this is generally not true as Xe135_m1 has two + # decay modes, and Xe135 will also decay, but we've modified the depletion chain so + # that Xe135_m1 only decays to Xe135, and that Xe135 has has no decay modes + assert xe135_atoms[1] == pytest.approx(xe135m1_atoms[1], rel=tol) + assert cs135m1_atoms[0] == pytest.approx(cs135m1_atoms[2] * 2, rel=tol) + assert cs135_atoms[2] == pytest.approx(cs135m1_atoms[2], rel=tol) diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index 3b93b6879a0..bf6b19ee5e7 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -74,6 +74,18 @@ def test_add_components(): with pytest.raises(ValueError): m.add_components({'H1': 1.0}, percent_type = 'oa') +def test_nuclides_to_ignore(run_in_tmpdir): + """Test nuclides_to_ignore when exporting a material to XML""" + m = openmc.Material() + m.add_nuclide('U235', 1.0) + m.add_nuclide('H1', 1.0) + m.add_nuclide('O16', 1.0) + + mats = openmc.Materials([m]) + mats.export_to_xml(nuclides_to_ignore=['H1']) + + test_mats = openmc.Materials.from_xml() + assert 'H1' not in test_mats[0].get_nuclides() def test_remove_nuclide(): """Test removing nuclides."""