From d6a4a3ecb9d14f90e3232c55360d5998beccd343 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 16 Nov 2022 15:10:11 -0600 Subject: [PATCH 1/3] Produce light particles from decay --- openmc/deplete/chain.py | 22 +++++++--- openmc/deplete/independent_operator.py | 6 +-- openmc/deplete/stepresult.py | 26 ++++++----- .../unit_tests/test_deplete_decay_products.py | 43 +++++++++++++++++++ 4 files changed, 79 insertions(+), 18 deletions(-) create mode 100644 tests/unit_tests/test_deplete_decay_products.py diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 33ac8a342..f439cedc3 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -630,15 +630,27 @@ class Chain: # Gain from radioactive decay if nuc.n_decay_modes != 0: - for _, target, branching_ratio in nuc.decay_modes: - # Allow for total annihilation for debug purposes - if target is not None: - branch_val = branching_ratio * decay_constant + for decay_type, target, branching_ratio in nuc.decay_modes: + branch_val = branching_ratio * decay_constant - if branch_val != 0.0: + # Allow for total annihilation for debug purposes + if branch_val != 0.0: + if target is not None: k = self.nuclide_dict[target] matrix[k, i] += branch_val + # Produce alphas and protons from decay + if 'alpha' in decay_type: + k = self.nuclide_dict.get('He4') + if k is not None: + count = decay_type.count('alpha') + matrix[k, i] += count * branch_val + elif 'p' in decay_type: + k = self.nuclide_dict.get('H1') + if k is not None: + count = decay_type.count('p') + matrix[k, i] += count * branch_val + if nuc.name in rates.index_nuc: # Extract all reactions for this nuclide in this cell nuc_ind = rates.index_nuc[nuc.name] diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index f628a7db1..a8249a98e 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -150,7 +150,7 @@ class IndependentOperator(OpenMCOperator): @classmethod def from_nuclides(cls, volume, nuclides, micro_xs, - chain_file, + chain_file=None, nuc_units='atom/b-cm', keff=None, normalization_mode='fission-q', @@ -170,9 +170,9 @@ class IndependentOperator(OpenMCOperator): values. micro_xs : MicroXS One-group microscopic cross sections. - chain_file : str + chain_file : str, optional Path to the depletion chain XML file. - nuc_units : {'atom/cm3', 'atom/b-cm'} + nuc_units : {'atom/cm3', 'atom/b-cm'}, optional Units for nuclide concentration. keff : 2-tuple of float, optional keff eigenvalue and uncertainty from transport calculation. diff --git a/openmc/deplete/stepresult.py b/openmc/deplete/stepresult.py index c8b35a1fa..042fa46a1 100644 --- a/openmc/deplete/stepresult.py +++ b/openmc/deplete/stepresult.py @@ -318,10 +318,11 @@ class StepResult: chunks=(1, 1, n_mats, n_nuc_number), dtype='float64') - handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn), - maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn), - chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn), - dtype='float64') + if n_nuc_rxn > 0 and n_rxn > 0: + handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn), + maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn), + chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn), + dtype='float64') handle.create_dataset("eigenvalues", (1, n_stages, 2), maxshape=(None, n_stages, 2), dtype='float64') @@ -358,7 +359,9 @@ class StepResult: # Grab handles number_dset = handle["/number"] - rxn_dset = handle["/reaction rates"] + has_reactions = ("reaction rates" in handle) + if has_reactions: + rxn_dset = handle["/reaction rates"] eigenvalues_dset = handle["/eigenvalues"] time_dset = handle["/time"] source_rate_dset = handle["/source_rate"] @@ -375,9 +378,10 @@ class StepResult: number_shape[0] = new_shape number_dset.resize(number_shape) - rxn_shape = list(rxn_dset.shape) - rxn_shape[0] = new_shape - rxn_dset.resize(rxn_shape) + if has_reactions: + rxn_shape = list(rxn_dset.shape) + rxn_shape[0] = new_shape + rxn_dset.resize(rxn_shape) eigenvalues_shape = list(eigenvalues_dset.shape) eigenvalues_shape[0] = new_shape @@ -407,7 +411,8 @@ class StepResult: high = max(inds) for i in range(n_stages): number_dset[index, i, low:high+1] = self.data[i] - rxn_dset[index, i, low:high+1] = self.rates[i] + if has_reactions: + rxn_dset[index, i, low:high+1] = self.rates[i] if comm.rank == 0: eigenvalues_dset[index, i] = self.k[i] if comm.rank == 0: @@ -483,7 +488,8 @@ class StepResult: for i in range(results.n_stages): rate = ReactionRates(results.index_mat, rxn_nuc_to_ind, rxn_to_ind, True) - rate[:] = handle["/reaction rates"][step, i, :, :, :] + if "reaction rates" in handle: + rate[:] = handle["/reaction rates"][step, i, :, :, :] results.rates.append(rate) return results diff --git a/tests/unit_tests/test_deplete_decay_products.py b/tests/unit_tests/test_deplete_decay_products.py new file mode 100644 index 000000000..34ca1b067 --- /dev/null +++ b/tests/unit_tests/test_deplete_decay_products.py @@ -0,0 +1,43 @@ +import openmc.deplete +import pytest + + +def test_deplete_decay_products(run_in_tmpdir): + # Create chain file with H1, He4, and Li5 + with open('test_chain.xml', 'w') as chain_file: + chain_file.write(""" + + + + + + + + + """) + + # Create depletion operator with no reactions + micro_xs = openmc.deplete.MicroXS() + op = openmc.deplete.IndependentOperator.from_nuclides( + volume=1.0, + nuclides={'Li5': 1.0}, + micro_xs=micro_xs, + chain_file='test_chain.xml', + normalization_mode='source-rate' + ) + + # Create time-integrator and integrate + integrator = openmc.deplete.PredictorIntegrator( + op, timesteps=[1.0], source_rates=[0.0], timestep_units='d' + ) + integrator.integrate(final_step=False) + + # Get concentration of H1 and He4 + results = openmc.deplete.Results('depletion_results.h5') + _, h1 = results.get_atoms("1", "H1") + _, he4 = results.get_atoms("1", "He4") + + # Since we started with 1e24 atoms of Li5, we should have 1e24 atoms of both + # H1 and He4 + assert h1[1] == pytest.approx(1e24) + assert he4[1] == pytest.approx(1e24) From ea0892f0f3e53169a8ea80e424a5002c3d643935 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 21 Nov 2022 15:20:12 -0600 Subject: [PATCH 2/3] Updated docstring to better explain chain_file and micro_xs --- openmc/deplete/independent_operator.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index a8249a98e..d11730240 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -31,6 +31,9 @@ class IndependentOperator(OpenMCOperator): passed to an integrator class, such as :class:`openmc.deplete.CECMIntegrator`. + Note that passing an empty :class:`~openmc.deplete.MicroXS` instance to the + ``micro_xs`` argument allows a decay-only calculation to be run. + .. versionadded:: 0.13.1 Parameters @@ -38,9 +41,11 @@ class IndependentOperator(OpenMCOperator): materials : openmc.Materials Materials to deplete. micro_xs : MicroXS - One-group microscopic cross sections in [b] . + One-group microscopic cross sections in [b]. If the + :class:`~openmc.deplete.MicroXS` is empty, a decay-only calculation will + be run. chain_file : str - Path to the depletion chain XML file. Defaults to + Path to the depletion chain XML file. Defaults to ``openmc.config['chain_file']``. keff : 2-tuple of float, optional keff eigenvalue and uncertainty from transport calculation. @@ -169,9 +174,12 @@ class IndependentOperator(OpenMCOperator): Dictionary with nuclide names as keys and nuclide concentrations as values. micro_xs : MicroXS - One-group microscopic cross sections. + One-group microscopic cross sections in [b]. If the + :class:`~openmc.deplete.MicroXS` is empty, a decay-only calculation + will be run. chain_file : str, optional - Path to the depletion chain XML file. + Path to the depletion chain XML file. Defaults to + ``openmc.config['chain_file']``. nuc_units : {'atom/cm3', 'atom/b-cm'}, optional Units for nuclide concentration. keff : 2-tuple of float, optional From 102288f94a99f4fd702ab43e9f2999b688c0f803 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 21 Nov 2022 15:36:34 -0600 Subject: [PATCH 3/3] Apply @yardasol suggestions from code review Co-authored-by: Olek <45364492+yardasol@users.noreply.github.com> --- openmc/deplete/independent_operator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index d11730240..0c950b31b 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -42,7 +42,7 @@ class IndependentOperator(OpenMCOperator): Materials to deplete. micro_xs : MicroXS One-group microscopic cross sections in [b]. If the - :class:`~openmc.deplete.MicroXS` is empty, a decay-only calculation will + :class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation will be run. chain_file : str Path to the depletion chain XML file. Defaults to @@ -175,7 +175,7 @@ class IndependentOperator(OpenMCOperator): values. micro_xs : MicroXS One-group microscopic cross sections in [b]. If the - :class:`~openmc.deplete.MicroXS` is empty, a decay-only calculation + :class:`~openmc.deplete.MicroXS` object is empty, a decay-only calculation will be run. chain_file : str, optional Path to the depletion chain XML file. Defaults to