From a585a7d63dadb7a40eeba077064a4eb5a3184e09 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 26 Sep 2024 21:13:34 -0500
Subject: [PATCH] Avoid writing subnormal nuclide densities to XML

---
 openmc/material.py                | 14 ++++++++++++--
 tests/unit_tests/test_material.py | 14 ++++++++++++++
 2 files changed, 26 insertions(+), 2 deletions(-)

diff --git a/openmc/material.py b/openmc/material.py
index 5b958c75a68..5dff15ccb9a 100644
--- a/openmc/material.py
+++ b/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
@@ -24,6 +25,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'])
 
@@ -1295,10 +1299,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
 
diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py
index 44c0730fbd2..9e82fa4bd51 100644
--- a/tests/unit_tests/test_material.py
+++ b/tests/unit_tests/test_material.py
@@ -661,3 +661,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