Avoid writing subnormal nuclide densities to XML (#3144)

openmc/material.py

@@ -5,6 +5,7 @@
 from numbers import Real
 from pathlib import Path
 import re
+import sys
 import warnings
 
 import lxml.etree as ET
@@ -26,6 +27,9 @@
 DENSITY_UNITS = ('g/cm3', 'g/cc', 'kg/m3', 'atom/b-cm', 'atom/cm3', 'sum',
                  'macro')
 
+# Smallest normalized floating point number
+_SMALLEST_NORMAL = sys.float_info.min
+
 
 NuclideTuple = namedtuple('NuclideTuple', ['name', 'percent', 'percent_type'])
 
@@ -1339,10 +1343,16 @@ def _get_nuclide_xml(self, nuclide: NuclideTuple) -> ET.Element:
         xml_element = ET.Element("nuclide")
         xml_element.set("name", nuclide.name)
 
+        # Prevent subnormal numbers from being written to XML, which causes an
+        # exception on the C++ side when calling std::stod
+        val = nuclide.percent
+        if abs(val) < _SMALLEST_NORMAL:
+            val = 0.0
+
         if nuclide.percent_type == 'ao':
-            xml_element.set("ao", str(nuclide.percent))
+            xml_element.set("ao", str(val))
         else:
-            xml_element.set("wo", str(nuclide.percent))
+            xml_element.set("wo", str(val))
 
         return xml_element
 

tests/unit_tests/test_material.py

@@ -683,3 +683,17 @@ def intensity(src):
     stable.add_nuclide('Gd156', 1.0)
     stable.volume = 1.0
     assert stable.get_decay_photon_energy() is None
+
+
+def test_avoid_subnormal(run_in_tmpdir):
+    # Write a materials.xml with a material that has a nuclide density that is
+    # represented as a subnormal floating point value
+    mat = openmc.Material()
+    mat.add_nuclide('H1', 1.0)
+    mat.add_nuclide('H2', 1.0e-315)
+    mats = openmc.Materials([mat])
+    mats.export_to_xml()
+
+    # When read back in, the density should be zero
+    mats = openmc.Materials.from_xml()
+    assert mats[0].get_nuclide_atom_densities()['H2'] == 0.0