From c8bf27af4236636d72ca16e839f139d61c053f89 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 22 May 2019 14:15:26 -0500 Subject: [PATCH] Add more thermal scattering unit tests and clean up some HDF5-related methods --- openmc/data/thermal.py | 104 ++++++++----------- openmc/data/thermal_angle_energy.py | 71 +++++++++++++ tests/unit_tests/test_data_thermal.py | 140 +++++++++++++++++++++++++- 3 files changed, 251 insertions(+), 64 deletions(-) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 5ec9b83659..6f39e478c8 100644 --- a/openmc/data/thermal.py +++ b/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 diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index 3537b91fa4..05393e7bb7 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/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']) diff --git a/tests/unit_tests/test_data_thermal.py b/tests/unit_tests/test_data_thermal.py index d8b21c99c6..63bb98f6a7 100644 --- a/tests/unit_tests/test_data_thermal.py +++ b/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'