diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 6f5c006c24..287738774c 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -61,6 +61,7 @@ Core Functions atomic_mass atomic_weight + decay_constant dose_coefficients gnd_name half_life diff --git a/openmc/cell.py b/openmc/cell.py index 103a34fe63..4b419e1c6e 100644 --- a/openmc/cell.py +++ b/openmc/cell.py @@ -214,8 +214,8 @@ class Cell(IDManagerMixin): self._atoms = self._fill.get_nuclide_atom_densities() # Convert to total number of atoms - for key, nuclide in self._atoms.items(): - atom = nuclide[1] * self._volume * 1.0e+24 + for key, atom_per_bcm in self._atoms.items(): + atom = atom_per_bcm * self._volume * 1.0e+24 self._atoms[key] = atom elif self.fill_type == 'distribmat': @@ -224,12 +224,12 @@ class Cell(IDManagerMixin): partial_volume = self.volume / len(self.fill) self._atoms = OrderedDict() for mat in self.fill: - for key, nuclide in mat.get_nuclide_atom_densities().items(): + for key, atom_per_bcm in mat.get_nuclide_atom_densities().items(): # To account for overlap of nuclides between distribmat # we need to append new atoms to any existing value # hence it is necessary to ask for default. atom = self._atoms.setdefault(key, 0) - atom += nuclide[1] * partial_volume * 1.0e+24 + atom += atom_per_bcm * partial_volume * 1.0e+24 self._atoms[key] = atom else: diff --git a/openmc/data/data.py b/openmc/data/data.py index 3a3cf58268..71b93293e7 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -3,7 +3,7 @@ import json import os import re from pathlib import Path -from math import sqrt +from math import sqrt, log from warnings import warn # Isotopic abundances from Meija J, Coplen T B, et al, "Isotopic compositions @@ -202,7 +202,7 @@ _GND_NAME_RE = re.compile(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)') # Used in half_life function as a cache _HALF_LIFE = {} - +_LOG_TWO = log(2.0) def atomic_mass(isotope): """Return atomic mass of isotope in atomic mass units. @@ -283,6 +283,8 @@ def half_life(isotope): Half-life values are from the `ENDF/B-VIII.0 decay sublibrary `_. + .. versionadded:: 0.13.1 + Parameters ---------- isotope : str @@ -302,6 +304,35 @@ def half_life(isotope): return _HALF_LIFE.get(isotope.lower()) + +def decay_constant(isotope): + """Return decay constant of isotope in [s^-1] + + Decay constants are based on half-life values from the + :func:`~openmc.data.half_life` function. When the isotope is stable, a decay + constant of zero is returned. + + .. versionadded:: 0.13.1 + + Parameters + ---------- + isotope : str + Name of isotope, e.g., 'Pu239' + + Returns + ------- + float + Decay constant of isotope in [s^-1] + + See also + -------- + openmc.data.half_life + + """ + t = half_life(isotope) + return _LOG_TWO / t if t else 0.0 + + def water_density(temperature, pressure=0.1013): """Return the density of liquid water at a given temperature and pressure. diff --git a/openmc/data/decay.py b/openmc/data/decay.py index 5ca2b235a5..51a594bc80 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -465,11 +465,10 @@ class Decay(EqualityMixin): @property def decay_constant(self): - if hasattr(self.half_life, 'n'): - return log(2.)/self.half_life - else: - mu, sigma = self.half_life - return ufloat(log(2.)/mu, log(2.)/mu**2*sigma) + if self.half_life.n == 0.0: + name = self.nuclide['name'] + raise ValueError(f"{name} is listed as unstable but has a zero half-life.") + return log(2.)/self.half_life @property def decay_energy(self): diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 327c762e08..7ef874e87f 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -528,9 +528,9 @@ class Operator(TransportOperator): """ mat_id = str(mat.id) - for nuclide, density in mat.get_nuclide_atom_densities().values(): - number = density * 1.0e24 - self.number.set_atom_density(mat_id, nuclide, number) + for nuclide, atom_per_bcm in mat.get_nuclide_atom_densities().items(): + atom_per_cc = atom_per_bcm * 1.0e24 + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) def _set_number_from_results(self, mat, prev_res): """Extracts material nuclides and number densities. @@ -556,16 +556,15 @@ class Operator(TransportOperator): # Merge lists of nuclides, with the same order for every calculation geom_nuc_densities.update(depl_nuc) - for nuclide in geom_nuc_densities.keys(): + for nuclide, atom_per_bcm in geom_nuc_densities.items(): if nuclide in depl_nuc: concentration = prev_res.get_atoms(mat_id, nuclide)[1][-1] volume = prev_res[-1].volume[mat_id] - number = concentration / volume + atom_per_cc = concentration / volume else: - density = geom_nuc_densities[nuclide][1] - number = density * 1.0e24 + atom_per_cc = atom_per_bcm * 1.0e24 - self.number.set_atom_density(mat_id, nuclide, number) + self.number.set_atom_density(mat_id, nuclide, atom_per_cc) def initial_condition(self): """Performs final setup and returns initial condition. diff --git a/openmc/material.py b/openmc/material.py index 3c2a98ffd3..e8285de7c2 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -3,7 +3,6 @@ from collections.abc import Iterable from copy import deepcopy from numbers import Real from pathlib import Path -import math import re import warnings from xml.etree import ElementTree as ET @@ -86,6 +85,9 @@ class Material(IDManagerMixin): fissionable_mass : float Mass of fissionable nuclides in the material in [g]. Requires that the :attr:`volume` attribute is set. + activity : float + Activity of the material in [Bq]. Requires that the :attr:`volume` + attribute is set. """ @@ -142,20 +144,10 @@ class Material(IDManagerMixin): return string - @property def activity(self): """Returns the total activity of the material in Becquerels.""" - - atoms_per_barn_cm = self.get_nuclide_atom_densities() - total_activity = 0 - for key, value in atoms_per_barn_cm.items(): - half_life = openmc.data.half_life(key) - if half_life: - total_activity += value[1] / half_life - total_activity *= math.log(2) * 1e24 * self.volume - - return total_activity + return sum(self.get_nuclide_activity().values()) @property def name(self): @@ -259,10 +251,10 @@ class Material(IDManagerMixin): if self.volume is None: raise ValueError("Volume must be set in order to determine mass.") density = 0.0 - for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values(): + for nuc, atoms_per_bcm in self.get_nuclide_atom_densities().items(): Z = openmc.data.zam(nuc)[0] if Z >= 90: - density += 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \ + density += 1e24 * atoms_per_bcm * openmc.data.atomic_mass(nuc) \ / openmc.data.AVOGADRO return density*self.volume @@ -786,11 +778,15 @@ class Material(IDManagerMixin): """Returns all nuclides in the material and their atomic densities in units of atom/b-cm + .. versionchanged:: 0.13.1 + The values in the dictionary were changed from a tuple containing + the nuclide name and the density to just the density. + Returns ------- nuclides : dict - Dictionary whose keys are nuclide names and values are tuples of - (nuclide, density in atom/b-cm) + Dictionary whose keys are nuclide names and values are densities in + [atom/b-cm] """ @@ -850,10 +846,46 @@ class Material(IDManagerMixin): nuclides = OrderedDict() for n, nuc in enumerate(nucs): - nuclides[nuc] = (nuc, nuc_densities[n]) + nuclides[nuc] = nuc_densities[n] return nuclides + def get_nuclide_activity(self): + """Return activity in [Bq] for each nuclide in the material + + .. versionadded:: 0.13.1 + + Returns + ------- + dict + Dictionary whose keys are nuclide names and values are activity in + [Bq]. + """ + activity = {} + for nuclide, atoms in self.get_nuclide_atoms().items(): + inv_seconds = openmc.data.decay_constant(nuclide) + activity[nuclide] = inv_seconds * atoms + return activity + + def get_nuclide_atoms(self): + """Return number of atoms of each nuclide in the material + + .. versionadded:: 0.13.1 + + Returns + ------- + dict + Dictionary whose keys are nuclide names and values are number of + atoms present in the material. + + """ + if self.volume is None: + raise ValueError("Volume must be set in order to determine atoms.") + atoms = {} + for nuclide, atom_per_bcm in self.get_nuclide_atom_densities().items(): + atoms[nuclide] = 1.0e24 * atom_per_bcm * self.volume + return atoms + def get_mass_density(self, nuclide=None): """Return mass density of one or all nuclides @@ -870,9 +902,9 @@ class Material(IDManagerMixin): """ mass_density = 0.0 - for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values(): + for nuc, atoms_per_bcm in self.get_nuclide_atom_densities().items(): if nuclide is None or nuclide == nuc: - density_i = 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \ + density_i = 1e24 * atoms_per_bcm * openmc.data.atomic_mass(nuc) \ / openmc.data.AVOGADRO mass_density += density_i return mass_density @@ -1092,7 +1124,7 @@ class Material(IDManagerMixin): nuclides_per_cc = defaultdict(float) mass_per_cc = defaultdict(float) for mat, wgt in zip(materials, wgts): - for nuc, atoms_per_bcm in mat.get_nuclide_atom_densities().values(): + for nuc, atoms_per_bcm in mat.get_nuclide_atom_densities().items(): nuc_per_cc = wgt*1.e24*atoms_per_bcm nuclides_per_cc[nuc] += nuc_per_cc mass_per_cc[nuc] += nuc_per_cc*openmc.data.atomic_mass(nuc) / \ diff --git a/openmc/plotter.py b/openmc/plotter.py index 92acaf5bb5..16760868e3 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -542,15 +542,11 @@ def _calculate_cexs_elem_mat(this, types, temperature=294., if isinstance(this, openmc.Material): # Expand elements in to nuclides with atomic densities - nuclides = this.get_nuclide_atom_densities() - # For ease of processing split out the nuclide and its fraction - nuc_fractions = {nuclide[1][0]: nuclide[1][1] - for nuclide in nuclides.items()} + nuc_fractions = this.get_nuclide_atom_densities() # Create a dict of [nuclide name] = nuclide object to carry forward # with a common nuclides format between openmc.Material and # openmc.Element objects - nuclides = {nuclide[1][0]: nuclide[1][0] - for nuclide in nuclides.items()} + nuclides = {nuclide: nuclide for nuclide in nuc_fractions} else: # Expand elements in to nuclides with atomic densities nuclides = this.expand(1., 'ao', enrichment=enrichment, @@ -885,13 +881,13 @@ def _calculate_mgxs_elem_mat(this, types, library, orders=None, # Check to see if we have nuclides/elements or a macroscopic object if this._macroscopic is not None: # We have macroscopics - nuclides = {this._macroscopic: (this._macroscopic, this.density)} + nuclides = {this._macroscopic: this.density} else: # Expand elements in to nuclides with atomic densities nuclides = this.get_nuclide_atom_densities() # For ease of processing split out nuc and nuc_density - nuc_fraction = [nuclide[1][1] for nuclide in nuclides.items()] + nuc_fraction = list(nuclides.values()) else: T = temperature # Expand elements in to nuclides with atomic densities diff --git a/tests/unit_tests/test_data_misc.py b/tests/unit_tests/test_data_misc.py index 6a81fb1e08..c3f7e3cfff 100644 --- a/tests/unit_tests/test_data_misc.py +++ b/tests/unit_tests/test_data_misc.py @@ -1,6 +1,6 @@ #!/usr/bin/env python -from collections.abc import Mapping +from math import log import os from pathlib import Path @@ -126,3 +126,14 @@ def test_zam(): assert openmc.data.zam('Am242_m10') == (95, 242, 10) with pytest.raises(ValueError): openmc.data.zam('garbage') + + +def test_half_life(): + assert openmc.data.half_life('H2') is None + assert openmc.data.half_life('U235') == pytest.approx(2.22102e16) + assert openmc.data.half_life('Am242') == pytest.approx(57672.0) + assert openmc.data.half_life('Am242_m1') == pytest.approx(4449622000.0) + assert openmc.data.decay_constant('H2') == 0.0 + assert openmc.data.decay_constant('U235') == pytest.approx(log(2.0)/2.22102e16) + assert openmc.data.decay_constant('Am242') == pytest.approx(log(2.0)/57672.0) + assert openmc.data.decay_constant('Am242_m1') == pytest.approx(log(2.0)/4449622000.0) diff --git a/tests/unit_tests/test_deplete_activation.py b/tests/unit_tests/test_deplete_activation.py index 0b82a5fbcf..cb3b86b9c8 100644 --- a/tests/unit_tests/test_deplete_activation.py +++ b/tests/unit_tests/test_deplete_activation.py @@ -92,7 +92,7 @@ def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts w = model.geometry.get_materials_by_name('tungsten')[0] atom_densities = w.get_nuclide_atom_densities() - atom_per_cc = 1e24 * atom_densities['W186'][1] # Density in atom/cm^3 + atom_per_cc = 1e24 * atom_densities['W186'] # Density in atom/cm^3 n0 = atom_per_cc * w.volume # Absolute number of atoms # Pick a random irradiation time and then determine necessary source rate to diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index c15a4b9e55..fd0e4f6247 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -100,7 +100,7 @@ def test_add_elements_by_formula(): m.add_elements_from_formula('Li4SiO4') # checking the ratio of elements is 4:1:4 for Li:Si:O elem = defaultdict(float) - for nuclide, adens in m.get_nuclide_atom_densities().values(): + for nuclide, adens in m.get_nuclide_atom_densities().items(): if nuclide.startswith("Li"): elem["Li"] += adens if nuclide.startswith("Si"): @@ -117,7 +117,7 @@ def test_add_elements_by_formula(): 'O16': 0.443386, 'O17': 0.000168} nuc_dens = m.get_nuclide_atom_densities() for nuclide in ref_dens: - assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2) # testing the correct nuclides are added to the Material when enriched m = openmc.Material() @@ -129,7 +129,7 @@ def test_add_elements_by_formula(): 'O16': 0.443386, 'O17': 0.000168} nuc_dens = m.get_nuclide_atom_densities() for nuclide in ref_dens: - assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2) # testing the use of brackets m = openmc.Material() @@ -137,7 +137,7 @@ def test_add_elements_by_formula(): # checking the ratio of elements is 2:2:6 for Mg:N:O elem = defaultdict(float) - for nuclide, adens in m.get_nuclide_atom_densities().values(): + for nuclide, adens in m.get_nuclide_atom_densities().items(): if nuclide.startswith("Mg"): elem["Mg"] += adens if nuclide.startswith("N"): @@ -155,7 +155,7 @@ def test_add_elements_by_formula(): 'O16': 0.599772, 'O17': 0.000227} nuc_dens = m.get_nuclide_atom_densities() for nuclide in ref_dens: - assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2) # testing non integer multiplier results in a value error m = openmc.Material() @@ -289,12 +289,21 @@ def test_get_nuclide_densities(uo2): def test_get_nuclide_atom_densities(uo2): - nucs = uo2.get_nuclide_atom_densities() - for nuc, density in nucs.values(): + for nuc, density in uo2.get_nuclide_atom_densities().items(): assert nuc in ('U235', 'O16') assert density > 0 +def test_get_nuclide_atoms(): + mat = openmc.Material() + mat.add_nuclide('Li6', 1.0) + mat.set_density('atom/cm3', 3.26e20) + mat.volume = 100.0 + + atoms = mat.get_nuclide_atoms() + assert atoms['Li6'] == pytest.approx(mat.density * mat.volume) + + def test_mass(): m = openmc.Material() m.add_nuclide('Zr90', 1.0, 'wo') @@ -341,7 +350,7 @@ def test_borated_water(): 'O16':2.4672e-02} nuc_dens = m.get_nuclide_atom_densities() for nuclide in ref_dens: - assert nuc_dens[nuclide][1] == pytest.approx(ref_dens[nuclide], 1e-2) + assert nuc_dens[nuclide] == pytest.approx(ref_dens[nuclide], 1e-2) assert m.id == 50 # Test the Celsius conversion. diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py index 10417b20ed..9b05ed2a5c 100644 --- a/tests/unit_tests/test_model.py +++ b/tests/unit_tests/test_model.py @@ -391,16 +391,14 @@ def test_py_lib_attributes(run_in_tmpdir, pin_model_attributes, mpi_intracomm): assert openmc.lib.materials[1].get_density('atom/b-cm') == \ pytest.approx(0.06891296988603757, abs=1e-13) mat_a_dens = np.sum( - [v[1] for v in test_model.materials[0]. - get_nuclide_atom_densities().values()]) + list(test_model.materials[0].get_nuclide_atom_densities().values())) assert mat_a_dens == pytest.approx(0.06891296988603757, abs=1e-8) # Change the density test_model.update_densities(['UO2'], 2.) assert openmc.lib.materials[1].get_density('atom/b-cm') == \ pytest.approx(2., abs=1e-13) mat_a_dens = np.sum( - [v[1] for v in test_model.materials[0]. - get_nuclide_atom_densities().values()]) + list(test_model.materials[0].get_nuclide_atom_densities().values())) assert mat_a_dens == pytest.approx(2., abs=1e-8) # Now lets do the cell temperature updates. @@ -441,7 +439,7 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm): test_model = openmc.Model(geom, mats, settings, tals, plots) initial_mat = mats[0].clone() - initial_u = initial_mat.get_nuclide_atom_densities()['U235'][1] + initial_u = initial_mat.get_nuclide_atom_densities()['U235'] # Note that the chain file includes only U-235 fission to a stable Xe136 w/ # a yield of 100%. Thus all the U235 we lose becomes Xe136 @@ -453,8 +451,8 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm): operator_kwargs=op_kwargs, power=1., output=False) # Get the new Xe136 and U235 atom densities - after_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1] - after_u = mats[0].get_nuclide_atom_densities()['U235'][1] + after_xe = mats[0].get_nuclide_atom_densities()['Xe136'] + after_u = mats[0].get_nuclide_atom_densities()['U235'] assert after_xe + after_u == pytest.approx(initial_u, abs=1e-15) assert test_model.is_initialized is False @@ -462,7 +460,7 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm): mats[0].nuclides.clear() densities = initial_mat.get_nuclide_atom_densities() tot_density = 0. - for nuc, density in densities.values(): + for nuc, density in densities.items(): mats[0].add_nuclide(nuc, density) tot_density += density mats[0].set_density('atom/b-cm', tot_density) @@ -473,8 +471,8 @@ def test_deplete(run_in_tmpdir, pin_model_attributes, mpi_intracomm): operator_kwargs=op_kwargs, power=1., output=False) # Get the new Xe136 and U235 atom densities - after_lib_xe = mats[0].get_nuclide_atom_densities()['Xe136'][1] - after_lib_u = mats[0].get_nuclide_atom_densities()['U235'][1] + after_lib_xe = mats[0].get_nuclide_atom_densities()['Xe136'] + after_lib_u = mats[0].get_nuclide_atom_densities()['U235'] assert after_lib_xe + after_lib_u == pytest.approx(initial_u, abs=1e-15) assert test_model.is_initialized is True