Add more thermal scattering unit tests and clean up some HDF5-related methods

openmc/data/thermal.py

@@ -327,30 +327,36 @@ class ThermalScatteringReaction(EqualityMixin):
         self.xs = xs
         self.distribution = distribution
 
-    def to_hdf5(self, group):
+    def to_hdf5(self, group, name):
         """Write thermal scattering reaction to HDF5
 
         Parameters
         ----------
         group : h5py.Group
             HDF5 group to write to
+        name : {'elastic', 'inelastic'}
+            Name of reaction to write
 
         """
         for T, xs in self.xs.items():
-            Tgroup = group.create_group(_temperature_str(T))
-            xs.to_hdf5(Tgroup, 'xs')
-            self.distribution[T].to_hdf5(Tgroup)
+            Tgroup = group.require_group(T)
+            rx_group = Tgroup.create_group(name)
+            xs.to_hdf5(rx_group, 'xs')
+            dgroup = rx_group.create_group('distribution')
+            self.distribution[T].to_hdf5(dgroup)
 
     @classmethod
-    def from_hdf5(cls, group, temperatures):
+    def from_hdf5(cls, group, name, temperatures):
         """Generate thermal scattering reaction data from HDF5
 
         Parameters
         ----------
         group : h5py.Group
             HDF5 group to read from
-        temperatures : Iterable of float
-            Temperatures in [K] to read
+        name : {'elastic', 'inelastic'}
+            Name of the reaction to read
+        temperatures : Iterable of str
+            Temperatures to read
 
         Returns
         -------
@@ -361,9 +367,12 @@ class ThermalScatteringReaction(EqualityMixin):
         xs = {}
         distribution = {}
         for T in temperatures:
-            Tgroup = group[_temperature_str(T)]
-            xs[T] = Function1D.from_hdf5(Tgroup)
-            distribution[T] = AngleEnergy.from_hdf5(Tgroup)
+            rx_group = group[T][name]
+            xs[T] = Function1D.from_hdf5(rx_group['xs'])
+            if isinstance(xs[T], CoherentElastic):
+                distribution[T] = CoherentElasticAE(xs[T])
+            else:
+                distribution[T] = AngleEnergy.from_hdf5(rx_group['distribution'])
         return cls(xs, distribution)
 
 
@@ -441,36 +450,23 @@ class ThermalScattering(EqualityMixin):
 
         """
         # Open file and write version
-        f = h5py.File(str(path), mode, libver=libver)
-        f.attrs['filetype'] = np.string_('data_thermal')
-        f.attrs['version'] = np.array(HDF5_VERSION)
+        with h5py.File(str(path), mode, libver=libver) as f:
+            f.attrs['filetype'] = np.string_('data_thermal')
+            f.attrs['version'] = np.array(HDF5_VERSION)
 
-        # Write basic data
-        g = f.create_group(self.name)
-        g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
-        g.attrs['energy_max'] = self.energy_max
-        g.attrs['nuclides'] = np.array(self.nuclides, dtype='S')
-        ktg = g.create_group('kTs')
-        for i, temperature in enumerate(self.temperatures):
-            ktg.create_dataset(temperature, data=self.kTs[i])
+            # Write basic data
+            g = f.create_group(self.name)
+            g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
+            g.attrs['energy_max'] = self.energy_max
+            g.attrs['nuclides'] = np.array(self.nuclides, dtype='S')
+            ktg = g.create_group('kTs')
+            for i, temperature in enumerate(self.temperatures):
+                ktg.create_dataset(temperature, data=self.kTs[i])
 
-        for T in self.temperatures:
-            Tg = g.create_group(T)
-            # Write thermal elastic scattering
+            # Write elastic/inelastic reaction data
             if self.elastic is not None:
-                elastic_group = Tg.create_group('elastic')
-                self.elastic.xs[T].to_hdf5(elastic_group, 'xs')
-                dgroup = elastic_group.create_group('distribution')
-                self.elastic.distribution[T].to_hdf5(dgroup)
-
-            # Write thermal inelastic scattering
-            if self.inelastic is not None:
-                inelastic_group = Tg.create_group('inelastic')
-                self.inelastic.xs[T].to_hdf5(inelastic_group, 'xs')
-                dgroup = inelastic_group.create_group('distribution')
-                self.inelastic.distribution[T].to_hdf5(dgroup)
-
-        f.close()
+                self.elastic.to_hdf5(g, 'elastic')
+            self.inelastic.to_hdf5(g, 'inelastic')
 
     def add_temperature_from_ace(self, ace_or_filename, name=None):
         """Add data to the ThermalScattering object from an ACE file at a
@@ -565,33 +561,15 @@ class ThermalScattering(EqualityMixin):
         table.nuclides = [nuc.decode() for nuc in group.attrs['nuclides']]
 
         # Read thermal elastic scattering
-        elastic_xs = {}
-        elastic_dist = {}
-        inelastic_xs = {}
-        inelastic_dist = {}
-        for T in table.temperatures:
-            Tgroup = group[T]
-            if 'elastic' in Tgroup:
-                elastic_group = Tgroup['elastic']
+        if 'elastic' in group[table.temperatures[0]]:
+            table.elastic = ThermalScatteringReaction.from_hdf5(
+                group, 'elastic', table.temperatures
+            )
 
-                # Cross section
-                elastic_xs[T] = Function1D.from_hdf5(elastic_group['xs'])
-                if isinstance(elastic_xs[T], CoherentElastic):
-                    elastic_dist[T] = CoherentElasticAE(elastic_xs[T])
-                else:
-                    dgroup = elastic_group['distribution']
-                    elastic_dist[T] = AngleEnergy.from_hdf5(dgroup)
-
-            # Read thermal inelastic scattering
-            if 'inelastic' in Tgroup:
-                inelastic_group = Tgroup['inelastic']
-                inelastic_xs[T] = Function1D.from_hdf5(inelastic_group['xs'])
-                inelastic_dist[T] = AngleEnergy.from_hdf5(
-                    inelastic_group['distribution'])
-
-        if elastic_xs:
-            table.elastic = ThermalScatteringReaction(elastic_xs, elastic_dist)
-        table.inelastic = ThermalScatteringReaction(inelastic_xs, inelastic_dist)
+        # Read thermal inelastic scattering
+        table.inelastic = ThermalScatteringReaction.from_hdf5(
+            group, 'inelastic', table.temperatures
+        )
 
         return table
 

openmc/data/thermal_angle_energy.py

@@ -34,6 +34,14 @@ class CoherentElasticAE(AngleEnergy):
         self.coherent_xs = coherent_xs
 
     def to_hdf5(self, group):
+        """Write coherent elastic distribution to an HDF5 group
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to write to
+
+        """
         group.attrs['type'] = np.string_('coherent_elastic')
         group['coherent_xs'] = group.parent['xs']
 
@@ -67,11 +75,32 @@ class IncoherentElasticAE(AngleEnergy):
         self.debye_waller = debye_waller
 
     def to_hdf5(self, group):
+        """Write incoherent elastic distribution to an HDF5 group
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to write to
+
+        """
         group.attrs['type'] = np.string_('incoherent_elastic')
         group.create_dataset('debye_waller', data=self.debye_waller)
 
     @classmethod
     def from_hdf5(cls, group):
+        """Generate incoherent elastic distribution from HDF5 data
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to read from
+
+        Returns
+        -------
+        openmc.data.IncoherentElasticAE
+            Incoherent elastic distribution
+
+        """
         return cls(group['debye_waller'])
 
 
@@ -88,11 +117,32 @@ class IncoherentElasticAEDiscrete(AngleEnergy):
         self.mu_out = mu_out
 
     def to_hdf5(self, group):
+        """Write discrete incoherent elastic distribution to an HDF5 group
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to write to
+
+        """
         group.attrs['type'] = np.string_('incoherent_elastic_discrete')
         group.create_dataset('mu_out', data=self.mu_out)
 
     @classmethod
     def from_hdf5(cls, group):
+        """Generate discrete incoherent elastic distribution from HDF5 data
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to read from
+
+        Returns
+        -------
+        openmc.data.IncoherentElasticAEDiscrete
+            Discrete incoherent elastic distribution
+
+        """
         return cls(group['mu_out'][()])
 
 
@@ -124,6 +174,14 @@ class IncoherentInelasticAEDiscrete(AngleEnergy):
         self.skewed = skewed
 
     def to_hdf5(self, group):
+        """Write discrete incoherent inelastic distribution to an HDF5 group
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to write to
+
+        """
         group.attrs['type'] = np.string_('incoherent_inelastic_discrete')
         group.create_dataset('energy_out', data=self.energy_out)
         group.create_dataset('mu_out', data=self.mu_out)
@@ -131,6 +189,19 @@ class IncoherentInelasticAEDiscrete(AngleEnergy):
 
     @classmethod
     def from_hdf5(cls, group):
+        """Generate discrete incoherent inelastic distribution from HDF5 data
+
+        Parameters
+        ----------
+        group : h5py.Group
+            HDF5 group to read from
+
+        Returns
+        -------
+        openmc.data.IncoherentInelasticAEDiscrete
+            Discrete incoherent inelastic distribution
+
+        """
         energy_out = group['energy_out'][()]
         mu_out = group['mu_out'][()]
         skewed = bool(group['skewed'])

tests/unit_tests/test_data_thermal.py

@@ -1,7 +1,9 @@
 #!/usr/bin/env python
 
 from collections.abc import Callable
+from math import exp
 import os
+import random
 
 import numpy as np
 import pytest
@@ -29,18 +31,43 @@ def graphite():
 
 @pytest.fixture(scope='module')
 def h2o_njoy():
+    """H in H2O generated using NJOY."""
     path_h1 = os.path.join(_ENDF_DATA, 'neutrons', 'n-001_H_001.endf')
     path_h2o = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-HinH2O.endf')
     return openmc.data.ThermalScattering.from_njoy(
         path_h1, path_h2o, temperatures=[293.6, 500.0])
 
 
+@pytest.fixture(scope='module')
+def hzrh():
+    """H in ZrH thermal scattering data."""
+    filename = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-HinZrH.endf')
+    return openmc.data.ThermalScattering.from_endf(filename)
+
+
+@pytest.fixture(scope='module')
+def hzrh_njoy():
+    """H in ZrH genertaed using NJOY."""
+    path_h1 = os.path.join(_ENDF_DATA, 'neutrons', 'n-001_H_001.endf')
+    path_hzrh = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-HinZrH.endf')
+    return openmc.data.ThermalScattering.from_njoy(
+        path_h1, path_hzrh, temperatures=[296.0], iwt=0)
+
+
+@pytest.fixture(scope='module')
+def sio2():
+    """SiO2 thermal scattering data."""
+    filename = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-SiO2.endf')
+    return openmc.data.ThermalScattering.from_endf(filename)
+
+
 def test_h2o_attributes(h2o):
     assert h2o.name == 'c_H_in_H2O'
     assert h2o.nuclides == ['H1']
     assert h2o.temperatures == ['294K']
     assert h2o.atomic_weight_ratio == pytest.approx(0.999167)
     assert h2o.energy_max == pytest.approx(4.46)
+    assert isinstance(repr(h2o), str)
 
 
 def test_h2o_xs(h2o):
@@ -69,15 +96,33 @@ def test_graphite_xs(graphite):
     assert elastic([1e-3, 1.0]) == pytest.approx([0.0, 0.62586153])
 
 
-def test_export_to_hdf5(tmpdir, h2o_njoy, graphite):
+def test_graphite_njoy():
+    path_c0 = os.path.join(_ENDF_DATA, 'neutrons', 'n-006_C_000.endf')
+    path_gr = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-graphite.endf')
+    graphite = openmc.data.ThermalScattering.from_njoy(
+        path_c0, path_gr, temperatures=[296.0])
+    assert graphite.nuclides == ['C0', 'C12', 'C13']
+    assert graphite.atomic_weight_ratio == pytest.approx(11.898)
+    assert graphite.energy_max == pytest.approx(2.02)
+    assert graphite.temperatures == ['296K']
+
+
+def test_export_to_hdf5(tmpdir, h2o_njoy, hzrh_njoy, graphite):
     filename = str(tmpdir.join('water.h5'))
     h2o_njoy.export_to_hdf5(filename)
     assert os.path.exists(filename)
 
+    # Graphite covers export of coherent elastic data
     filename = str(tmpdir.join('graphite.h5'))
     graphite.export_to_hdf5(filename)
     assert os.path.exists(filename)
 
+    # H in ZrH covers export of incoherent elastic data, and discrete incoherent
+    # inelastic angle-energy distribution
+    filename = str(tmpdir.join('hzrh.h5'))
+    hzrh_njoy.export_to_hdf5(filename)
+    assert os.path.exists(filename)
+
 
 def test_continuous_dist(h2o_njoy):
     for temperature, dist in h2o_njoy.inelastic.distribution.items():
@@ -85,6 +130,95 @@ def test_continuous_dist(h2o_njoy):
         assert isinstance(dist, openmc.data.IncoherentInelasticAE)
 
 
+def test_h2o_endf():
+    filename = os.path.join(_ENDF_DATA, 'thermal_scatt', 'tsl-HinH2O.endf')
+    h2o = openmc.data.ThermalScattering.from_endf(filename)
+    assert not h2o.elastic
+    assert h2o.atomic_weight_ratio == pytest.approx(0.99917)
+    assert h2o.energy_max == pytest.approx(3.99993)
+    assert h2o.temperatures == ['294K', '350K', '400K', '450K', '500K', '550K',
+                                '600K', '650K', '800K']
+
+
+def test_hzrh_attributes(hzrh):
+    assert hzrh.atomic_weight_ratio == pytest.approx(0.99917)
+    assert hzrh.energy_max == pytest.approx(1.9734)
+    assert hzrh.temperatures == ['296K', '400K', '500K', '600K', '700K', '800K',
+                                 '1000K', '1200K']
+
+
+def test_hzrh_elastic(hzrh):
+    rx = hzrh.elastic
+    for temperature, func in rx.xs.items():
+        assert temperature.endswith('K')
+        assert isinstance(func, openmc.data.IncoherentElastic)
+
+    xs = rx.xs['296K']
+    sig_b, W = xs.bound_xs, xs.debye_waller
+    assert sig_b == pytest.approx(81.98006)
+    assert W == pytest.approx(8.486993)
+    for i in range(10):
+        E = random.uniform(0.0, hzrh.energy_max)
+        assert xs(E) == pytest.approx(sig_b/2 * ((1 - exp(-4*E*W))/(2*E*W)))
+
+    for temperature, dist in rx.distribution.items():
+        assert temperature.endswith('K')
+        assert dist.debye_waller > 0.0
+
+
+def test_hzrh_njoy(hzrh_njoy):
+    hzrh = hzrh_njoy
+    assert hzrh.atomic_weight_ratio == pytest.approx(0.999167)
+    assert hzrh.energy_max == pytest.approx(1.855)
+    assert hzrh.temperatures == ['296K']
+
+    # Check incoherent elastic distribution
+    d = hzrh.elastic.distribution['296K']
+    assert np.all((-1.0 <= d.mu_out) & (d.mu_out <= 1.0))
+
+    # Check incoherent inelastic distribution
+    d = hzrh.inelastic.distribution['296K']
+    assert d.skewed
+    assert np.all((-1.0 < d.mu_out) & (d.mu_out < 1.0))
+    assert np.all((0.0 <= d.energy_out) & (d.energy_out < 3*hzrh.energy_max))
+
+
+def test_sio2_attributes(sio2):
+    assert sio2.atomic_weight_ratio == pytest.approx(27.84423)
+    assert sio2.energy_max == pytest.approx(2.46675)
+    assert sio2.temperatures == ['294K', '350K', '400K', '500K', '800K',
+                                 '1000K', '1200K']
+
+
+def test_sio2_elastic(sio2):
+    rx = sio2.elastic
+    for temperature, func in rx.xs.items():
+        assert temperature.endswith('K')
+        assert isinstance(func, openmc.data.CoherentElastic)
+    xs = rx.xs['294K']
+    assert len(xs) == 317
+    assert xs.bragg_edges[0] == pytest.approx(0.000711634)
+    assert xs.factors[0] == pytest.approx(2.6958e-14)
+
+    # Below first bragg edge, cross section should be zero
+    E = xs.bragg_edges[0] / 2.0
+    assert xs(E) == 0.0
+
+    # Between bragg edges, cross section is P/E where P is the factor
+    E = (xs.bragg_edges[0] + xs.bragg_edges[1]) / 2.0
+    P = xs.factors[0]
+    assert xs(E) == pytest.approx(P / E)
+
+    # Check the last Bragg edge
+    E = 1.1 * xs.bragg_edges[-1]
+    P = xs.factors[-1]
+    assert xs(E) == pytest.approx(P / E)
+
+    for temperature, dist in rx.distribution.items():
+        assert temperature.endswith('K')
+        assert dist.coherent_xs is rx.xs[temperature]
+
+
 def test_get_thermal_name():
     f = openmc.data.get_thermal_name
     # Names which are recognized
@@ -97,5 +231,9 @@ def test_get_thermal_name():
         assert f('graphite') == 'c_Graphite'
         assert f('D_in_D2O') == 'c_D_in_D2O'
 
+        # Not in values, but very close
+        assert f('hluci') == 'c_H_in_C5O2H8'
+        assert f('ortho_d') == 'c_ortho_D'
+
         # Names that don't remotely match anything
         assert f('boogie_monster') == 'c_boogie_monster'