Replacing PyNE with radioactivedecay for radionuclide decay calculati…

…ons (#1813)

* Added rd decay calculations to decay.py

* minor tweaks after review

* Ran black formatter

* removed pyne references from tardis files

* revert tqdm changes to util base file

* Add util to check if str is valid elem or nuclide

* fixed how inventories are updated, units to grams

* Updated decay test with more accurate decay data

* Fixed convert_abundances_format util test error

* Updated astropy, set radioactivedecay version

* revert astropy fixed version

Co-authored-by: Thomas Lemoine <[email protected]>

tardis/__init__.py

@@ -12,7 +12,6 @@
 # ----------------------------------------------------------------------------
 
 import sys
-import pyne.data
 
 # ----------------------------------------------------------------------------
 

tardis/io/decay.py

@@ -1,5 +1,6 @@
 import pandas as pd
-from pyne import nucname, material
+from radioactivedecay import Nuclide, Inventory
+from radioactivedecay.utils import Z_to_elem
 from astropy import units as u
 
 
@@ -20,56 +21,57 @@ def __init__(self, *args, **kwargs):
     def _constructor(self):
         return IsotopeAbundances
 
-    def _update_material(self):
+    def _update_inventory(self):
         self.comp_dicts = [dict() for i in range(len(self.columns))]
         for (atomic_number, mass_number), abundances in self.iterrows():
-            nuclear_symbol = f"{nucname.name(atomic_number)}{mass_number}"
+            nuclear_symbol = f"{Z_to_elem(atomic_number)}{mass_number}"
             for i in range(len(self.columns)):
                 self.comp_dicts[i][nuclear_symbol] = abundances[i]
 
     @classmethod
-    def from_materials(cls, materials):
+    def from_inventories(cls, inventories):
         multi_index_tuples = set([])
-        for material in materials:
+        for inventory in inventories:
             multi_index_tuples.update(
-                [cls.id_to_tuple(key) for key in material.keys()]
+                [cls.id_to_tuple(key) for key in inventory.contents.keys()]
             )
 
         index = pd.MultiIndex.from_tuples(
             multi_index_tuples, names=["atomic_number", "mass_number"]
         )
 
         abundances = pd.DataFrame(
-            data=0.0, index=index, columns=range(len(materials))
+            data=0.0, index=index, columns=range(len(inventories))
         )
 
-        for i, material in enumerate(materials):
-            for key, value in material.items():
-                abundances.loc[cls.id_to_tuple(key), i] = value
+        for i, inventory in enumerate(inventories):
+            for nuclide, abundance in inventory.masses("g").items():
+                abundances.loc[cls.id_to_tuple(nuclide), i] = abundance
 
         return cls(abundances)
 
     @staticmethod
     def id_to_tuple(atomic_id):
-        return nucname.znum(atomic_id), nucname.anum(atomic_id)
+        nuclide = Nuclide(atomic_id)
+        return nuclide.Z, nuclide.A
 
-    def to_materials(self):
+    def to_inventories(self):
         """
-        Convert DataFrame to a list of materials interpreting the MultiIndex as
+        Convert DataFrame to a list of inventories interpreting the MultiIndex as
         atomic_number and mass_number
 
         Returns
         -------
         list
-            list of pyne Materials
+            list of radioactivedecay Inventories
         """
 
         comp_dicts = [dict() for i in range(len(self.columns))]
         for (atomic_number, mass_number), abundances in self.iterrows():
-            nuclear_symbol = f"{nucname.name(atomic_number):s}{mass_number:d}"
+            nuclear_symbol = f"{Z_to_elem(atomic_number)}{mass_number}"
             for i in range(len(self.columns)):
                 comp_dicts[i][nuclear_symbol] = abundances[i]
-        return [material.Material(comp_dict) for comp_dict in comp_dicts]
+        return [Inventory(comp_dict, "g") for comp_dict in comp_dicts]
 
     def decay(self, t):
         """
@@ -86,14 +88,12 @@ def decay(self, t):
             Decayed abundances
         """
 
-        materials = self.to_materials()
+        inventories = self.to_inventories()
         t_second = (
             u.Quantity(t, u.day).to(u.s).value - self.time_0.to(u.s).value
         )
-        decayed_materials = [item.decay(t_second) for item in materials]
-        for i in range(len(materials)):
-            materials[i].update(decayed_materials[i])
-        df = IsotopeAbundances.from_materials(materials)
+        decayed_inventories = [item.decay(t_second) for item in inventories]
+        df = IsotopeAbundances.from_inventories(decayed_inventories)
         df.sort_index(inplace=True)
         return df
 

tardis/io/logger/logger.py

@@ -1,12 +1,8 @@
 import logging
 import sys
-import warnings
-import pyne.data
 
 from tardis.io.logger.colored_logger import ColoredFormatter, formatter_message
 
-warnings.filterwarnings("ignore", category=pyne.utils.QAWarning)
-
 logging.captureWarnings(True)
 logger = logging.getLogger("tardis")
 

tardis/io/model_reader.py

@@ -5,14 +5,15 @@
 from numpy import recfromtxt, genfromtxt
 import pandas as pd
 from astropy import units as u
-from pyne import nucname
+from radioactivedecay import Nuclide
+from radioactivedecay.utils import Z_DICT, elem_to_Z
 
 import logging
 
 # Adding logging support
 logger = logging.getLogger(__name__)
 
-from tardis.util.base import parse_quantity
+from tardis.util.base import parse_quantity, is_valid_nuclide_or_elem
 
 
 class ConfigurationError(Exception):
@@ -165,14 +166,15 @@ def read_uniform_abundances(abundances_section, no_of_shells):
         if element_symbol_string == "type":
             continue
         try:
-            if element_symbol_string in nucname.name_zz:
-                z = nucname.name_zz[element_symbol_string]
+            if element_symbol_string in Z_DICT.values():
+                z = elem_to_Z(element_symbol_string)
                 abundance.loc[z] = float(
                     abundances_section[element_symbol_string]
                 )
             else:
-                mass_no = nucname.anum(element_symbol_string)
-                z = nucname.znum(element_symbol_string)
+                nuc = Nuclide(element_symbol_string)
+                mass_no = nuc.A
+                z = nuc.Z
                 isotope_abundance.loc[(z, mass_no), :] = float(
                     abundances_section[element_symbol_string]
                 )
@@ -469,12 +471,13 @@ def read_csv_isotope_abundances(
     )
 
     for element_symbol_string in df.index[skip_columns:]:
-        if element_symbol_string in nucname.name_zz:
-            z = nucname.name_zz[element_symbol_string]
+        if element_symbol_string in Z_DICT.values():
+            z = elem_to_Z(element_symbol_string)
             abundance.loc[z, :] = df.loc[element_symbol_string].tolist()
         else:
-            z = nucname.znum(element_symbol_string)
-            mass_no = nucname.anum(element_symbol_string)
+            nuc = Nuclide(element_symbol_string)
+            z = nuc.Z
+            mass_no = nuc.A
             isotope_abundance.loc[(z, mass_no), :] = df.loc[
                 element_symbol_string
             ].tolist()
@@ -498,9 +501,7 @@ def parse_csv_abundances(csvy_data):
     """
 
     abundance_col_names = [
-        name
-        for name in csvy_data.columns
-        if nucname.iselement(name) or nucname.isnuclide(name)
+        name for name in csvy_data.columns if is_valid_nuclide_or_elem(name)
     ]
     df = csvy_data.loc[:, abundance_col_names]
 
@@ -520,12 +521,13 @@ def parse_csv_abundances(csvy_data):
     )
 
     for element_symbol_string in df.index[0:]:
-        if element_symbol_string in nucname.name_zz:
-            z = nucname.name_zz[element_symbol_string]
+        if element_symbol_string in Z_DICT.values():
+            z = elem_to_Z(element_symbol_string)
             abundance.loc[z, :] = df.loc[element_symbol_string].tolist()
         else:
-            z = nucname.znum(element_symbol_string)
-            mass_no = nucname.anum(element_symbol_string)
+            nuc = Nuclide(element_symbol_string)
+            z = nuc.Z
+            mass_no = nuc.A
             isotope_abundance.loc[(z, mass_no), :] = df.loc[
                 element_symbol_string
             ].tolist()

tardis/io/tests/test_decay.py

@@ -15,12 +15,12 @@ def simple_abundance_model():
 
 def test_simple_decay(simple_abundance_model):
     decayed_abundance = simple_abundance_model.decay(100)
-    assert_almost_equal(decayed_abundance.loc[26, 56][0], 0.55752)
-    assert_almost_equal(decayed_abundance.loc[26, 56][1], 0.55752)
-    assert_almost_equal(decayed_abundance.loc[27, 56][0], 0.4423791)
-    assert_almost_equal(decayed_abundance.loc[27, 56][1], 0.4423791)
-    assert_almost_equal(decayed_abundance.loc[28, 56][0], 1.1086e-05)
-    assert_almost_equal(decayed_abundance.loc[28, 56][1], 1.1086e-05)
+    assert_almost_equal(decayed_abundance.loc[26, 56][0], 0.55754786)
+    assert_almost_equal(decayed_abundance.loc[26, 56][1], 0.55754786)
+    assert_almost_equal(decayed_abundance.loc[27, 56][0], 0.44235126)
+    assert_almost_equal(decayed_abundance.loc[27, 56][1], 0.44235126)
+    assert_almost_equal(decayed_abundance.loc[28, 56][0], 1.10859709e-05)
+    assert_almost_equal(decayed_abundance.loc[28, 56][1], 1.10859709e-05)
 
 
 @pytest.fixture

tardis/model/base.py

@@ -5,7 +5,7 @@
 from astropy import units as u
 from tardis import constants
 
-from tardis.util.base import quantity_linspace
+from tardis.util.base import quantity_linspace, is_valid_nuclide_or_elem
 from tardis.io.parsers.csvy import load_csvy
 from tardis.io.model_reader import (
     read_density_file,
@@ -18,7 +18,6 @@
 from tardis.io.util import HDFWriterMixin
 from tardis.io.decay import IsotopeAbundances
 from tardis.model.density import HomologousDensity
-from pyne import nucname
 
 logger = logging.getLogger(__name__)
 
@@ -224,7 +223,7 @@ def __init__(
             self._dilution_factor = 0.5 * (
                 1
                 - np.sqrt(
-                    1 - (self.r_inner[0] ** 2 / self.r_middle ** 2).to(1).value
+                    1 - (self.r_inner[0] ** 2 / self.r_middle**2).to(1).value
                 )
             )
         else:
@@ -389,7 +388,7 @@ def abundance(self):
 
     @property
     def volume(self):
-        return ((4.0 / 3) * np.pi * (self.r_outer ** 3 - self.r_inner ** 3)).cgs
+        return ((4.0 / 3) * np.pi * (self.r_outer**3 - self.r_inner**3)).cgs
 
     @property
     def no_of_shells(self):
@@ -655,7 +654,7 @@ def from_csvy(cls, config):
                 [
                     name
                     for name in csvy_model_data.columns
-                    if nucname.iselement(name) or nucname.isnuclide(name)
+                    if is_valid_nuclide_or_elem(name)
                 ]
             )
             unsupported_columns = (
@@ -675,10 +674,10 @@ def from_csvy(cls, config):
             ), "CSVY field descriptions exist without corresponding csv data"
             if unsupported_columns != set():
                 logger.warning(
-                "The following columns are "
-                "specified in the csvy model file,"
-                 f" but are IGNORED by TARDIS: {str(unsupported_columns)}"
-            )
+                    "The following columns are "
+                    "specified in the csvy model file,"
+                    f" but are IGNORED by TARDIS: {str(unsupported_columns)}"
+                )
 
         time_explosion = config.supernova.time_explosion.cgs
 

tardis/util/base.py

@@ -9,7 +9,8 @@
 import yaml
 from tardis import constants
 from astropy import units as u
-from pyne import nucname
+from radioactivedecay import Nuclide, DEFAULTDATA
+from radioactivedecay.utils import parse_nuclide, Z_DICT
 
 import tardis
 from tardis.io.util import get_internal_data_path
@@ -180,7 +181,7 @@ def calculate_luminosity(
     )
 
     flux_density = np.trapz(flux, wavelength) * (flux_unit * wavelength_unit)
-    luminosity = (flux_density * 4 * np.pi * distance ** 2).to("erg/s")
+    luminosity = (flux_density * 4 * np.pi * distance**2).to("erg/s")
 
     return luminosity.value, wavelength.min(), wavelength.max()
 
@@ -305,7 +306,7 @@ def intensity_black_body(nu, T):
         Returns the intensity of the black body
     """
     beta_rad = 1 / (k_B_cgs * T)
-    coefficient = 2 * h_cgs / c_cgs ** 2
+    coefficient = 2 * h_cgs / c_cgs**2
     intensity = ne.evaluate(
         "coefficient * nu**3 / " "(exp(h_cgs * nu * beta_rad) -1 )"
     )
@@ -491,6 +492,33 @@ def reformat_element_symbol(element_string):
     return element_string[0].upper() + element_string[1:].lower()
 
 
+def is_valid_nuclide_or_elem(input_nuclide):
+    """
+    Parses nuclide string into symbol - mass number format and returns
+    whether the nuclide is either contained in the decay dataset or is a
+    raw element string.
+
+    Parameters
+    ----------
+    input_nuclide : str or int
+        Nuclide name string or element string.
+
+    Returns
+    -------
+    bool
+        Bool indicating if the input nuclide is contained in the decay dataset
+        or is a valid element.
+    """
+
+    try:
+        parse_nuclide(input_nuclide, DEFAULTDATA.nuclides, "ICRP-107")
+        is_nuclide = True
+    except:
+        is_nuclide = True if input_nuclide in Z_DICT.values() else False
+
+    return is_nuclide
+
+
 def quantity_linspace(start, stop, num, **kwargs):
     """
     Essentially the same input parameters as linspace, but
@@ -546,7 +574,7 @@ def convert_abundances_format(fname, delimiter=r"\s+"):
     # Drop shell index column
     df.drop(df.columns[0], axis=1, inplace=True)
     # Assign header row
-    df.columns = [nucname.name(i) for i in range(1, df.shape[1] + 1)]
+    df.columns = [Z_DICT[i] for i in range(1, df.shape[1] + 1)]
     return df
 
 

tardis/visualization/widgets/custom_abundance.py

@@ -6,11 +6,12 @@
 import ipywidgets as ipw
 import plotly.graph_objects as go
 from astropy import units as u
-from pyne import nucname
+from radioactivedecay import Nuclide
+from radioactivedecay.utils import Z_DICT, elem_to_Z
 from pathlib import Path
 
 import tardis
-from tardis.util.base import quantity_linspace
+from tardis.util.base import quantity_linspace, is_valid_nuclide_or_elem
 from tardis.io.config_reader import Configuration
 from tardis.model import Radial1DModel
 from tardis.model.density import (
@@ -999,9 +1000,7 @@ def input_symb_eventhandler(self, obj):
             return
 
         try:
-            if nucname.iselement(element_symbol_string) or nucname.isnuclide(
-                element_symbol_string
-            ):
+            if is_valid_nuclide_or_elem(element_symbol_string):
                 self.symb_warning.layout.visibility = "hidden"
                 self.btn_add_element.disabled = False
                 return
@@ -1024,12 +1023,13 @@ def on_btn_add_element(self, obj):
         """
         element_symbol_string = self.input_symb.value.capitalize()
 
-        if element_symbol_string in nucname.name_zz:
-            z = nucname.name_zz[element_symbol_string]
+        if element_symbol_string in Z_DICT.values():
+            z = elem_to_Z(element_symbol_string)
             self.data.abundance.loc[(z, ""), :] = 0
         else:
-            mass_no = nucname.anum(element_symbol_string)
-            z = nucname.znum(element_symbol_string)
+            nuc = Nuclide(element_symbol_string)
+            mass_no = nuc.A
+            z = nuc.Z
             self.data.abundance.loc[(z, mass_no), :] = 0
 
         self.data.abundance.sort_index(inplace=True)