Merge pull request #2079 from shimwell/adding_activity_to_material

Adding activity to material

docs/source/pythonapi/data.rst

@@ -63,6 +63,7 @@ Core Functions
     atomic_weight
     dose_coefficients
     gnd_name
+    half_life
     isotopes
     linearize
     thin

openmc/data/data.py

@@ -1,10 +1,11 @@
 import itertools
-from math import sqrt
+import json
 import os
 import re
+from pathlib import Path
+from math import sqrt
 from warnings import warn
 
-
 # Isotopic abundances from Meija J, Coplen T B, et al, "Isotopic compositions
 # of the elements 2013 (IUPAC Technical Report)", Pure. Appl. Chem. 88 (3),
 # pp. 293-306 (2013). The "representative isotopic abundance" values from
@@ -199,6 +200,9 @@ _ATOMIC_MASS = {}
 # Regex for GND nuclide names (used in zam function)
 _GND_NAME_RE = re.compile(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)')
 
+# Used in half_life function as a cache
+_HALF_LIFE = {}
+
 
 def atomic_mass(isotope):
     """Return atomic mass of isotope in atomic mass units.
@@ -273,6 +277,31 @@ def atomic_weight(element):
                          .format(element))
 
 
+def half_life(isotope):
+    """Return half-life of isotope in seconds or None if isotope is stable
+
+    Half-life values are from the `ENDF/B-VIII.0 decay sublibrary
+    <https://www.nndc.bnl.gov/endf-b8.0/download.html>`_.
+
+    Parameters
+    ----------
+    isotope : str
+        Name of isotope, e.g., 'Pu239'
+
+    Returns
+    -------
+    float
+        Half-life of isotope in [s]
+
+    """
+    global _HALF_LIFE
+    if not _HALF_LIFE:
+        # Load ENDF/B-VIII.0 data from JSON file
+        half_life_path = Path(__file__).with_name('half_life.json')
+        _HALF_LIFE = json.loads(half_life_path.read_text())
+
+    return _HALF_LIFE.get(isotope.lower())
+
 def water_density(temperature, pressure=0.1013):
     """Return the density of liquid water at a given temperature and pressure.
 

openmc/material.py

@@ -3,6 +3,7 @@ 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
@@ -141,6 +142,21 @@ 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
+
     @property
     def name(self):
         return self._name

setup.py

@@ -32,7 +32,7 @@ kwargs = {
     # Data files and libraries
     'package_data': {
         'openmc.lib': ['libopenmc.{}'.format(suffix)],
-        'openmc.data': ['mass16.txt', 'BREMX.DAT', '*.h5'],
+        'openmc.data': ['mass16.txt', 'BREMX.DAT', 'half_life.json', '*.h5'],
         'openmc.data.effective_dose': ['*.txt']
     },
 

tests/unit_tests/test_material.py

@@ -404,3 +404,30 @@ def test_mix_materials():
     assert m3.density == pytest.approx(dens3)
     assert m4.density == pytest.approx(dens4)
     assert m5.density == pytest.approx(dens5)
+
+
+def test_activity_of_stable():
+    """Creates a material with stable isotopes to checks the activity is 0"""
+    m1 = openmc.Material()
+    m1.add_element("Fe", 1)
+    m1.set_density('g/cm3', 1)
+    m1.volume = 1
+    assert m1.activity == 0
+
+
+def test_activity_of_tritium():
+    """Checks that 1g of tritium has the correct activity"""
+    m1 = openmc.Material()
+    m1.add_nuclide("H3", 1)
+    m1.set_density('g/cm3', 1)
+    m1.volume = 1
+    assert pytest.approx(m1.activity) == 3.559778e14 
+
+
+def test_activity_of_metastable():
+    """Checks that 1 mol of a Tc99_m1 nuclides has the correct activity"""
+    m1 = openmc.Material()
+    m1.add_nuclide("Tc99_m1", 1)
+    m1.set_density('g/cm3', 1)
+    m1.volume = 98.9
+    assert pytest.approx(m1.activity, rel=0.001) == 1.93e19