From fbacfa0e5c83e07dea4998898e4ab7188b81a53a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 18 Apr 2019 08:37:26 -0500 Subject: [PATCH 01/19] Add a from_endf method to ThermalScattering that returns a dict --- openmc/data/thermal.py | 92 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 91 insertions(+), 1 deletion(-) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 27ed7aa291..7eeca685a2 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -1,6 +1,8 @@ from collections.abc import Iterable +from collections import namedtuple from difflib import get_close_matches from numbers import Real +from io import StringIO import itertools import os import re @@ -13,7 +15,7 @@ import h5py import openmc.checkvalue as cv from openmc.mixin import EqualityMixin from openmc.stats import Discrete, Tabular -from . import HDF5_VERSION, HDF5_VERSION_MAJOR +from . import HDF5_VERSION, HDF5_VERSION_MAJOR, endf from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, NATURAL_ABUNDANCE from .ace import Table, get_table, Library from .angle_energy import AngleEnergy @@ -656,3 +658,91 @@ class ThermalScattering(EqualityMixin): data.add_temperature_from_ace(table) return data + + @classmethod + def from_endf(cls, filename): + data = {} + + ev = endf.Evaluation(filename) + if (7, 2) in ev.section: + file_obj = StringIO(ev.section[7, 2]) + lhtr = endf.get_head_record(file_obj)[2] + if lhtr == 1: + # coherent elastic + + # Get structure factor at first temperature + params, S = endf.get_tab1_record(file_obj) + t0 = params[0] + n_temps = params[2] + structure_factor = {t0: S} + + # Get structure factor for subsequent temperatures + for _ in range(n_temps): + params, S = endf.get_list_record(file_obj) + t = params[0] + structure_factor[t] = S + data['structure_factor'] = structure_factor + + elif lhtr == 2: + # incoherent elastic + params, W = endf.get_tab1_record(file_obj) + characteristic_xs = params[0] + + data['debye_waller'] = W + data['characteristic_xs'] = characteristic_xs + + if (7, 4) in ev.section: + file_obj = StringIO(ev.section[7, 4]) + params = endf.get_head_record(file_obj) + data['symmetric'] = (params[4] == 0) + + # Get information about principal atom + params, B = endf.get_list_record(file_obj) + data['log'] = bool(params[2]) + data['free_atom_xs'] = B[0] + data['epsilon'] = B[1] + data['A0'] = B[2] + data['e_max'] = B[3] + data['M0'] = B[5] + + # Get information about non-principal atoms + n_non_principal = params[5] + data['non_principal'] = [] + NonPrincipal = namedtuple('NonPrincipal', ['func', 'xs', 'A', 'M']) + for i in range(1, n_non_principal + 1): + func = {0.0: 'SCT', 1.0: 'free gas', 2.0: 'diffusive'}[B[6*i]] + xs = B[6*i + 1] + A = B[6*i + 2] + M = B[6*i + 5] + data['non_principal'].append(NonPrincipal(func, xs, A, M)) + + # Get S(alpha,beta,T) + if data['free_atom_xs'] > 0.0: + params, _ = endf.get_tab2_record(file_obj) + n_beta = params[5] + sab = {'beta': np.empty(n_beta)} + for i in range(n_beta): + params, S = endf.get_tab1_record(file_obj) + t0, beta, lt = params[:3] + if i == 0: + sab['alpha'] = alpha = S.x + sab[t0] = np.empty((alpha.size, n_beta)) + sab['beta'][i] = beta + sab[t0][:, i] = S.y + for _ in range(lt): + params, S = endf.get_list_record(file_obj) + t = params[0] + if i == 0: + sab[t] = np.empty((alpha.size, n_beta)) + sab[t][:, i] = S + data['sab'] = sab + + # Get effective temperature for each atom + _, Teff = endf.get_tab1_record(file_obj) + data['effective_temperature'] = [Teff] + for atom in data['non_principal']: + if atom.func == 'SCT': + _, Teff = endf.get_tab1_record(file_obj) + data['effective_temperature'].append(Teff) + + return data From 3e5741e1bb688af4a7e56b6389a89ee440b4b749 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 18 Apr 2019 09:45:24 -0500 Subject: [PATCH 02/19] Fix CoherentElastic.__call__ --- openmc/data/thermal.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 7eeca685a2..762de3cd3e 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -145,10 +145,15 @@ class CoherentElastic(EqualityMixin): self.factors = factors def __call__(self, E): + idx = np.searchsorted(self.bragg_edges, E) - 1 if isinstance(E, Iterable): E = np.asarray(E) - idx = np.searchsorted(self.bragg_edges, E) - return self.factors[idx] / E + nonzero = idx >= 0 + xs = np.zeros_like(E) + xs[nonzero] = self.factors[idx[nonzero]] / E[nonzero] + return xs + else: + return self.factors[idx] / E if idx >= 0 else 0.0 def __len__(self): return len(self.bragg_edges) From 4ad1aeac3a35ef89ffadd9f700f6eca21b68eb26 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 19 Apr 2019 14:01:19 -0500 Subject: [PATCH 03/19] Fix potential bug in Decay class --- openmc/data/decay.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openmc/data/decay.py b/openmc/data/decay.py index bbb059f4fa..6a8ad39f09 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -423,7 +423,8 @@ class Decay(EqualityMixin): if spectrum['type'] == 'ec/beta+': di['positron_intensity'] = ufloat(*values[4:6]) elif spectrum['type'] == 'gamma': - di['internal_pair'] = ufloat(*values[4:6]) + if len(values) >= 6: + di['internal_pair'] = ufloat(*values[4:6]) if len(values) >= 8: di['total_internal_conversion'] = ufloat(*values[6:8]) if len(values) == 12: From 51c8c2e0bc3725f906e69ba9f94066f97c525079 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 19 Apr 2019 14:01:33 -0500 Subject: [PATCH 04/19] Refactor classes in openmc.data.ThermalScattering --- docs/source/pythonapi/data.rst | 54 ++- openmc/data/__init__.py | 2 +- openmc/data/angle_energy.py | 8 + openmc/data/function.py | 12 +- openmc/data/thermal.py | 609 ++++++++++++++++++++++++--------- 5 files changed, 511 insertions(+), 174 deletions(-) diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 5cf74bce0a..144d263066 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -13,18 +13,38 @@ Core Classes openmc.data.IncidentNeutron openmc.data.Reaction openmc.data.Product - openmc.data.Tabulated1D - openmc.data.FissionEnergyRelease - openmc.data.ThermalScattering - openmc.data.CoherentElastic openmc.data.FissionEnergyRelease openmc.data.DataLibrary - openmc.data.IncidentPhoton - openmc.data.PhotonReaction - openmc.data.AtomicRelaxation openmc.data.Decay openmc.data.FissionProductYields openmc.data.WindowedMultipole + openmc.data.ProbabilityTables + +The following classes are used for storing atomic data (incident photon cross +sections, atomic relaxation): + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + openmc.data.IncidentPhoton + openmc.data.PhotonReaction + openmc.data.AtomicRelaxation + + +The following classes are used for storing thermal neutron scattering data: + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + openmc.data.ThermalScattering + openmc.data.ThermalScatteringReaction + openmc.data.CoherentElastic + openmc.data.IncoherentElastic + Core Functions -------------- @@ -41,6 +61,22 @@ Core Functions openmc.data.water_density openmc.data.zam +One-dimensional Functions +------------------------- + +.. autosummary:: + :toctree: generated + :nosignatures: + :template: myclass.rst + + openmc.data.Function1D + openmc.data.Tabulated1D + openmc.data.Polynomial + openmc.data.Combination + openmc.data.Sum + openmc.data.Regions1D + openmc.data.ResonancesWithBackground + Angle-Energy Distributions -------------------------- @@ -66,6 +102,10 @@ Angle-Energy Distributions openmc.data.DiscretePhoton openmc.data.LevelInelastic openmc.data.ContinuousTabular + openmc.data.CoherentElasticAE + openmc.data.IncoherentElasticAE + openmc.data.IncoherentElasticAEDiscrete + openmc.data.IncoherentInelasticAEDiscrete Resonance Data -------------- diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index f52c3be755..92acc96354 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -16,7 +16,6 @@ from .decay import * from .reaction import * from . import ace from .angle_distribution import * -from .function import * from . import endf from .energy_distribution import * from .product import * @@ -33,3 +32,4 @@ from .resonance import * from .resonance_covariance import * from .multipole import * from .grid import * +from .function import * diff --git a/openmc/data/angle_energy.py b/openmc/data/angle_energy.py index d67cc6b266..19dddfe19e 100644 --- a/openmc/data/angle_energy.py +++ b/openmc/data/angle_energy.py @@ -35,6 +35,14 @@ class AngleEnergy(EqualityMixin, metaclass=ABCMeta): return openmc.data.KalbachMann.from_hdf5(group) elif dist_type == 'nbody': return openmc.data.NBodyPhaseSpace.from_hdf5(group) + elif dist_type == 'coherent_elastic': + return openmc.data.CoherentElasticAE.from_hdf5(group) + elif dist_type == 'incoherent_elastic': + return openmc.data.IncoherentElasticAE.from_hdf5(group) + elif dist_type == 'incoherent_elastic_discrete': + return openmc.data.IncoherentElasticAEDiscrete.from_hdf5(group) + elif dist_type == 'incoherent_inelastic_discrete': + return openmc.data.IncoherentInelasticAEDiscrete.from_hdf5(group) @staticmethod def from_ace(ace, location_dist, location_start, rx=None): diff --git a/openmc/data/function.py b/openmc/data/function.py index c40e2a699c..99ee20900a 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -440,6 +440,14 @@ class Tabulated1D(Function1D): class Polynomial(np.polynomial.Polynomial, Function1D): + """A power series class. + + Parameters + ---------- + coef : Iterable of float + Polynomial coefficients in order of increasing degree + + """ def to_hdf5(self, group, name='xy'): """Write polynomial function to an HDF5 group @@ -583,8 +591,8 @@ class Regions1D(EqualityMixin): Functions which are to be combined in a piecewise fashion breakpoints : Iterable of float The values of the dependent variable that define the domain of - each function. The *i*th and *(i+1)*th values are the limits of the - domain of the *i*th function. Values must be monotonically increasing. + each function. The `i`th and `(i+1)`th values are the limits of the + domain of the `i`th function. Values must be monotonically increasing. Attributes ---------- diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 762de3cd3e..fc3a967251 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -19,7 +19,7 @@ from . import HDF5_VERSION, HDF5_VERSION_MAJOR, endf from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, NATURAL_ABUNDANCE from .ace import Table, get_table, Library from .angle_energy import AngleEnergy -from .function import Tabulated1D +from .function import Tabulated1D, Function1D from .correlated import CorrelatedAngleEnergy from .njoy import make_ace_thermal @@ -121,22 +121,31 @@ def get_thermal_name(name): return 'c_' + name -class CoherentElastic(EqualityMixin): +class CoherentElastic(Function1D): r"""Coherent elastic scattering data from a crystalline material + The integrated cross section for coherent elastic scattering from a + powdered crystalline material may be represented as: + + .. math:: + \sigma(E,T) = \frac{1}{E} \sum\limits_{i=1}^{E_i < E} s_i(T) + + where :math:`s_i(T)` is proportional the structure factor in [eV-b] at + the moderator temperature :math:`T` in Kelvin. + Parameters ---------- bragg_edges : Iterable of float Bragg edge energies in eV factors : Iterable of float - Partial sum of structure factors, :math:`\sum\limits_{i=1}^{E_i".format(self.name) + return "".format(self.name) else: return "" + @staticmethod + def _temperature_str(T): + return "{}K".format(round(T)) + @property def temperatures(self): - return ["{}K".format(int(round(kT / K_BOLTZMANN))) for kT in self.kTs] + return [self._temperature_str(kT / K_BOLTZMANN) for kT in self.kTs] def export_to_hdf5(self, path, mode='a', libver='earliest'): """Export table to an HDF5 file. @@ -297,7 +561,6 @@ class ThermalScattering(EqualityMixin): g = f.create_group(self.name) g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio g.attrs['nuclides'] = np.array(self.nuclides, dtype='S') - g.attrs['secondary_mode'] = np.string_(self.secondary_mode) ktg = g.create_group('kTs') for i, temperature in enumerate(self.temperatures): ktg.create_dataset(temperature, data=self.kTs[i]) @@ -305,25 +568,18 @@ class ThermalScattering(EqualityMixin): for T in self.temperatures: Tg = g.create_group(T) # Write thermal elastic scattering - if self.elastic_xs: + if self.elastic is not None: elastic_group = Tg.create_group('elastic') - - self.elastic_xs[T].to_hdf5(elastic_group, 'xs') - if self.elastic_mu_out: - elastic_group.create_dataset('mu_out', - data=self.elastic_mu_out[T]) + 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_xs: + if self.inelastic is not None: inelastic_group = Tg.create_group('inelastic') - self.inelastic_xs[T].to_hdf5(inelastic_group, 'xs') - if self.secondary_mode in ('equal', 'skewed'): - inelastic_group.create_dataset('energy_out', - data=self.inelastic_e_out[T]) - inelastic_group.create_dataset('mu_out', - data=self.inelastic_mu_out[T]) - elif self.secondary_mode == 'continuous': - self.inelastic_dist[T].to_hdf5(inelastic_group) + self.inelastic.xs[T].to_hdf5(inelastic_group, 'xs') + dgroup = inelastic_group.create_group('distribution') + self.inelastic.distribution[T].to_hdf5(dgroup) f.close() @@ -363,20 +619,14 @@ class ThermalScattering(EqualityMixin): self.kTs += data.kTs # Add inelastic cross section and distributions - if strT in data.inelastic_xs: - self.inelastic_xs[strT] = data.inelastic_xs[strT] - if strT in data.inelastic_e_out: - self.inelastic_e_out[strT] = data.inelastic_e_out[strT] - if strT in data.inelastic_mu_out: - self.inelastic_mu_out[strT] = data.inelastic_mu_out[strT] - if strT in data.inelastic_dist: - self.inelastic_dist[strT] = data.inelastic_dist[strT] + if data.inelastic is not None: + self.inelastic.xs.update(data.inelastic.xs) + self.inelastic.distribution.update(data.inelastic.distribution) # Add elastic cross sectoin and angular distribution - if strT in data.elastic_xs: - self.elastic_xs[strT] = data.elastic_xs[strT] - if strT in data.elastic_mu_out: - self.elastic_mu_out[strT] = data.elastic_mu_out[strT] + if data.elastic is not None: + self.elastic.xs.update(data.elastic.xs) + self.elastic.distribution.update(data.elastic.distribution) @classmethod def from_hdf5(cls, group_or_filename): @@ -419,45 +669,39 @@ class ThermalScattering(EqualityMixin): name = group.name[1:] atomic_weight_ratio = group.attrs['atomic_weight_ratio'] kTg = group['kTs'] - kTs = [] - for temp in kTg: - kTs.append(kTg[temp][()]) - temperatures = [str(int(round(kT / K_BOLTZMANN))) + "K" for kT in kTs] + kTs = [dataset[()] for dataset in kTg.values()] table = cls(name, atomic_weight_ratio, kTs) table.nuclides = [nuc.decode() for nuc in group.attrs['nuclides']] - table.secondary_mode = group.attrs['secondary_mode'].decode() # Read thermal elastic scattering - for T in temperatures: + elastic_xs = {} + elastic_dist = {} + inelastic_xs = {} + inelastic_dist = {} + for T in table.temperatures: Tgroup = group[T] if 'elastic' in Tgroup: elastic_group = Tgroup['elastic'] # Cross section - elastic_xs_type = elastic_group['xs'].attrs['type'].decode() - if elastic_xs_type == 'Tabulated1D': - table.elastic_xs[T] = Tabulated1D.from_hdf5( - elastic_group['xs']) - elif elastic_xs_type == 'bragg': - table.elastic_xs[T] = CoherentElastic.from_hdf5( - elastic_group['xs']) - - # Angular distribution - if 'mu_out' in elastic_group: - table.elastic_mu_out[T] = elastic_group['mu_out'][()] + 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'] - table.inelastic_xs[T] = Tabulated1D.from_hdf5( - inelastic_group['xs']) - if table.secondary_mode in ('equal', 'skewed'): - table.inelastic_e_out[T] = inelastic_group['energy_out'][()] - table.inelastic_mu_out[T] = inelastic_group['mu_out'][()] - elif table.secondary_mode == 'continuous': - table.inelastic_dist[T] = AngleEnergy.from_hdf5( - inelastic_group) + 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) return table @@ -492,41 +736,32 @@ class ThermalScattering(EqualityMixin): # Assign temperature to the running list kTs = [ace.temperature*EV_PER_MEV] - temperatures = [str(int(round(ace.temperature*EV_PER_MEV - / K_BOLTZMANN))) + "K"] table = cls(name, ace.atomic_weight_ratio, kTs) + T = table.temperatures[0] # Incoherent inelastic scattering cross section idx = ace.jxs[1] n_energy = int(ace.xss[idx]) energy = ace.xss[idx+1 : idx+1+n_energy]*EV_PER_MEV xs = ace.xss[idx+1+n_energy : idx+1+2*n_energy] - table.inelastic_xs[temperatures[0]] = Tabulated1D(energy, xs) - - if ace.nxs[7] == 0: - table.secondary_mode = 'equal' - elif ace.nxs[7] == 1: - table.secondary_mode = 'skewed' - elif ace.nxs[7] == 2: - table.secondary_mode = 'continuous' + inelastic_xs = Tabulated1D(energy, xs) + # Incoherent inelastic angle-energy distribution + continuous = (ace.nxs[7] == 2) n_energy_out = ace.nxs[4] - if table.secondary_mode in ('equal', 'skewed'): + if not continuous: n_mu = ace.nxs[3] idx = ace.jxs[3] - table.inelastic_e_out[temperatures[0]] = \ - ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2): - n_mu + 2]*EV_PER_MEV - table.inelastic_e_out[temperatures[0]].shape = \ - (n_energy, n_energy_out) + energy_out = ace.xss[idx:idx + n_energy * n_energy_out * + (n_mu + 2): n_mu + 2]*EV_PER_MEV + energy_out.shape = (n_energy, n_energy_out) - table.inelastic_mu_out[temperatures[0]] = \ - ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2)] - table.inelastic_mu_out[temperatures[0]].shape = \ - (n_energy, n_energy_out, n_mu+2) - table.inelastic_mu_out[temperatures[0]] = \ - table.inelastic_mu_out[temperatures[0]][:, :, 1:] + mu_out = ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2)] + mu_out.shape = (n_energy, n_energy_out, n_mu+2) + mu_out = mu_out[:, :, 1:] + skewed = (ace.nxs[7] == 1) + distribution = IncoherentInelasticAEDiscrete(energy_out, mu_out, skewed) else: n_mu = ace.nxs[3] - 1 idx = ace.jxs[3] @@ -565,10 +800,14 @@ class ThermalScattering(EqualityMixin): # Create correlated angle-energy distribution breakpoints = [n_energy] interpolation = [2] - energy = table.inelastic_xs[temperatures[0]].x - table.inelastic_dist[temperatures[0]] = CorrelatedAngleEnergy( + energy = inelastic_xs.x + distribution = CorrelatedAngleEnergy( breakpoints, interpolation, energy, energy_out, mu_out) + table.inelastic = ThermalScatteringReaction( + {T: inelastic_xs}, {T: distribution} + ) + # Incoherent/coherent elastic scattering cross section idx = ace.jxs[4] n_mu = ace.nxs[6] + 1 @@ -579,22 +818,23 @@ class ThermalScattering(EqualityMixin): if ace.nxs[5] == 4: # Coherent elastic - table.elastic_xs[temperatures[0]] = CoherentElastic( - energy, P*EV_PER_MEV) + xs = CoherentElastic(energy, P*EV_PER_MEV) + distribution = CoherentElasticAE(xs) # Coherent elastic shouldn't have angular distributions listed assert n_mu == 0 else: # Incoherent elastic - table.elastic_xs[temperatures[0]] = Tabulated1D(energy, P) + xs = Tabulated1D(energy, P) # Angular distribution assert n_mu > 0 idx = ace.jxs[6] - table.elastic_mu_out[temperatures[0]] = \ - ace.xss[idx:idx + n_energy * n_mu] - table.elastic_mu_out[temperatures[0]].shape = \ - (n_energy, n_mu) + mu_out = ace.xss[idx:idx + n_energy * n_mu] + mu_out.shape = (n_energy, n_mu) + distribution = IncoherentElasticAEDiscrete(mu_out) + + table.elastic = ThermalScatteringReaction({T: xs}, {T: distribution}) # Get relevant nuclides -- NJOY only allows one to specify three # nuclides that the S(a,b) table applies to. Thus, for all elements @@ -619,7 +859,7 @@ class ThermalScattering(EqualityMixin): @classmethod def from_njoy(cls, filename, filename_thermal, temperatures=None, evaluation=None, evaluation_thermal=None, **kwargs): - """Generate incident neutron data by running NJOY. + """Generate thermal scattering data by running NJOY. Parameters ---------- @@ -665,11 +905,32 @@ class ThermalScattering(EqualityMixin): return data @classmethod - def from_endf(cls, filename): - data = {} + def from_endf(cls, ev_or_filename): + """Generate thermal scattering data from an ENDF file - ev = endf.Evaluation(filename) + Parameters + ---------- + ev_or_filename : openmc.data.endf.Evaluation or str + ENDF evaluation to read from. If given as a string, it is assumed to + be the filename for the ENDF file. + + Returns + ------- + openmc.data.ThermalScattering + Thermal scattering data + + """ + if isinstance(ev_or_filename, Evaluation): + ev = ev_or_filename + else: + ev = Evaluation(ev_or_filename) + + # Read coherent/incoherent elastic data + elastic = None if (7, 2) in ev.section: + xs = {} + distribution = {} + file_obj = StringIO(ev.section[7, 2]) lhtr = endf.get_head_record(file_obj)[2] if lhtr == 1: @@ -677,77 +938,97 @@ class ThermalScattering(EqualityMixin): # Get structure factor at first temperature params, S = endf.get_tab1_record(file_obj) - t0 = params[0] + strT = self._temperature_str(params[0]) n_temps = params[2] - structure_factor = {t0: S} + bragg_edges = S.x + xs[strT] = CoherentElastic(bragg_edges, S.y) + distribution = {strT: CoherentElasticAE(xs[strT])} # Get structure factor for subsequent temperatures for _ in range(n_temps): params, S = endf.get_list_record(file_obj) - t = params[0] - structure_factor[t] = S - data['structure_factor'] = structure_factor + strT = self._temperature_str(params[0]) + xs[strT] = CoherentElastic(bragg_edges, S) + distribution[strT] = CoherentElasticAE(xs[T]) elif lhtr == 2: # incoherent elastic params, W = endf.get_tab1_record(file_obj) characteristic_xs = params[0] + for T, debye_waller in zip(W.x, W.y): + strT = self._temperature_str(T) + xs[strT] = IncoherentElastic(characteristic_xs, debye_waller) + distribution[strT] = IncoherentElasticAE(debye_waller) - data['debye_waller'] = W - data['characteristic_xs'] = characteristic_xs + elastic = ThermalScatteringReaction(xs, distribution) - if (7, 4) in ev.section: - file_obj = StringIO(ev.section[7, 4]) - params = endf.get_head_record(file_obj) - data['symmetric'] = (params[4] == 0) + # Read incoherent inelastic data + assert (7, 4) in ev.section, 'No MF=7, MT=4 found in thermal scattering' + file_obj = StringIO(ev.section[7, 4]) + params = endf.get_head_record(file_obj) + data = {'symmetric': params[4] == 0} - # Get information about principal atom - params, B = endf.get_list_record(file_obj) - data['log'] = bool(params[2]) - data['free_atom_xs'] = B[0] - data['epsilon'] = B[1] - data['A0'] = B[2] - data['e_max'] = B[3] - data['M0'] = B[5] + # Get information about principal atom + params, B = endf.get_list_record(file_obj) + data['log'] = bool(params[2]) + data['free_atom_xs'] = B[0] + data['epsilon'] = B[1] + data['A0'] = awr = B[2] + data['e_max'] = B[3] + data['M0'] = B[5] - # Get information about non-principal atoms - n_non_principal = params[5] - data['non_principal'] = [] - NonPrincipal = namedtuple('NonPrincipal', ['func', 'xs', 'A', 'M']) - for i in range(1, n_non_principal + 1): - func = {0.0: 'SCT', 1.0: 'free gas', 2.0: 'diffusive'}[B[6*i]] - xs = B[6*i + 1] - A = B[6*i + 2] - M = B[6*i + 5] - data['non_principal'].append(NonPrincipal(func, xs, A, M)) + # Get information about non-principal atoms + n_non_principal = params[5] + data['non_principal'] = [] + NonPrincipal = namedtuple('NonPrincipal', ['func', 'xs', 'A', 'M']) + for i in range(1, n_non_principal + 1): + func = {0.0: 'SCT', 1.0: 'free gas', 2.0: 'diffusive'}[B[6*i]] + xs = B[6*i + 1] + A = B[6*i + 2] + M = B[6*i + 5] + data['non_principal'].append(NonPrincipal(func, xs, A, M)) - # Get S(alpha,beta,T) - if data['free_atom_xs'] > 0.0: - params, _ = endf.get_tab2_record(file_obj) - n_beta = params[5] - sab = {'beta': np.empty(n_beta)} - for i in range(n_beta): - params, S = endf.get_tab1_record(file_obj) - t0, beta, lt = params[:3] + # Get S(alpha,beta,T) + kTs = [] + if data['free_atom_xs'] > 0.0: + params, _ = endf.get_tab2_record(file_obj) + n_beta = params[5] + sab = {'beta': np.empty(n_beta)} + for i in range(n_beta): + params, S = endf.get_tab1_record(file_obj) + T0, beta, lt = params[:3] + if i == 0: + sab['alpha'] = alpha = S.x + sab[T0] = np.empty((alpha.size, n_beta)) + kTs.append(K_BOLTZMANN * T0) + sab['beta'][i] = beta + sab[T0][:, i] = S.y + for _ in range(lt): + params, S = endf.get_list_record(file_obj) + T = params[0] if i == 0: - sab['alpha'] = alpha = S.x - sab[t0] = np.empty((alpha.size, n_beta)) - sab['beta'][i] = beta - sab[t0][:, i] = S.y - for _ in range(lt): - params, S = endf.get_list_record(file_obj) - t = params[0] - if i == 0: - sab[t] = np.empty((alpha.size, n_beta)) - sab[t][:, i] = S - data['sab'] = sab + sab[T] = np.empty((alpha.size, n_beta)) + kTs.append(K_BOLTZMANN * T) + sab[T][:, i] = S + data['sab'] = sab - # Get effective temperature for each atom - _, Teff = endf.get_tab1_record(file_obj) - data['effective_temperature'] = [Teff] - for atom in data['non_principal']: - if atom.func == 'SCT': - _, Teff = endf.get_tab1_record(file_obj) - data['effective_temperature'].append(Teff) + # Get effective temperature for each atom + _, Teff = endf.get_tab1_record(file_obj) + data['effective_temperature'] = [Teff] + for atom in data['non_principal']: + if atom.func == 'SCT': + _, Teff = endf.get_tab1_record(file_obj) + data['effective_temperature'].append(Teff) - return data + name = ev.target['zsymam'].strip() + instance = cls(name, awr, kTs) + if elastic is not None: + instance.elastic = elastic + + # Currently we don't have a proper cross section or distribution for + # incoherent inelastic, so we just create an empty object and attach + # all the data as a dictionary + instance.inelastic = ThermalScatteringReaction(None, None) + instance.inelastic.data = data + + return instance From 8ea820619375b87eb970b32e083de07d10d9c6ae Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 19 Apr 2019 14:33:04 -0500 Subject: [PATCH 05/19] Move thermal angle-energy distributions into separate file --- openmc/data/thermal.py | 170 ++++------------------------ openmc/data/thermal_angle_energy.py | 136 ++++++++++++++++++++++ 2 files changed, 157 insertions(+), 149 deletions(-) create mode 100644 openmc/data/thermal_angle_energy.py diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index fc3a967251..bc2e1d9819 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -19,9 +19,12 @@ from . import HDF5_VERSION, HDF5_VERSION_MAJOR, endf from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, NATURAL_ABUNDANCE from .ace import Table, get_table, Library from .angle_energy import AngleEnergy -from .function import Tabulated1D, Function1D from .correlated import CorrelatedAngleEnergy +from .function import Tabulated1D, Function1D from .njoy import make_ace_thermal +from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE, + IncoherentElasticAEDiscrete, + IncoherentInelasticAEDiscrete) _THERMAL_NAMES = { @@ -71,6 +74,12 @@ _THERMAL_NAMES = { } +def _temperature_str(T): + # round() normally returns an int when called with a single argument, but + # numpy floats overload rounding to return another float + return "{}K".format(int(round(T))) + + def get_thermal_name(name): """Get proper S(a,b) table name, e.g. 'HH2O' -> 'c_H_in_H2O' @@ -112,12 +121,12 @@ def get_thermal_name(name): break warn('Thermal scattering material "{}" is not recognized. ' - 'Assigning a name of {}.'.format(name, match)) + 'Assigning a name of {}.'.format(name, match)) return match else: # OK, we give up. Just use the ACE name. warn('Thermal scattering material "{0}" is not recognized. ' - 'Assigning a name of c_{0}.'.format(name)) + 'Assigning a name of c_{0}.'.format(name)) return 'c_' + name @@ -315,7 +324,7 @@ class ThermalScatteringReaction(EqualityMixin): """ for T, xs in self.xs.items(): - Tgroup = group.create_group(str(round(T)) + "K") + Tgroup = group.create_group(_temperature_str(T)) xs.to_hdf5(Tgroup, 'xs') self.distribution[T].to_hdf5(Tgroup) @@ -339,145 +348,12 @@ class ThermalScatteringReaction(EqualityMixin): xs = {} distribution = {} for T in temperatures: - Tgroup = group[str(round(T)) + "K"] + Tgroup = group[_temperature_str(T)] xs[T] = Function1D.from_hdf5(Tgroup) distribution[T] = AngleEnergy.from_hdf5(Tgroup) return cls(xs, distribution) -class CoherentElasticAE(AngleEnergy): - r"""Differential cross section for coherent elastic scattering - - The differential cross section for coherent elastic scattering from a - powdered crystalline material may be represented as: - - .. math:: - \frac{d^2\sigma}{dE'd\Omega} (E\rightarrow E',\mu,T) = \frac{1}{E} \sum - \limits_{i=1}^{E_i < E} s_i(T) \delta(\mu - \mu_i) \delta (E - E') - /(2\pi) - - where :math:`E_i` are the energies of the Bragg edges in [eV], :math:`s_i(T)` - is the structure factor in [eV-b] at the moderator temperature :math:`T` - in [K], and :math:`\mu_i = 1 - 2E_i/E`. - - Parameters - ---------- - coherent_xs : openmc.data.CoherentElastic - Coherent elastic scattering cross section - - Attributes - ---------- - coherent_xs : openmc.data.CoherentElastic - Coherent elastic scattering cross section - - """ - def __init__(self, coherent_xs): - self.coherent_xs = coherent_xs - - def to_hdf5(self, group): - group.attrs['type'] = np.string_('coherent_elastic') - group['coherent_xs'] = group.parent['xs'] - - -class IncoherentElasticAE(AngleEnergy): - r"""Differential cross section for incoherent elastic scattering - - The differential cross section for incoherent elastic scattering may be - represented as: - - .. math:: - \frac{d^2\sigma}{dE'd\Omega} (E\rightarrow E',\mu,T) = \frac{\sigma_b} - {4\pi} e^{-2EW'(T)(1-\mu)} \delta(E - E') - - where :math:`\sigma_b` is the characteristic cross section in [b] and - :math:`W'(T)` is the Debye-Waller integral divided by the atomic mass in - [eV\ :math:`^{-1}`]. - - Parameters - ---------- - debye_waller : float - Debye-Waller integral in [eV\ :math:`^{-1}`] - - Attributes - ---------- - debye_waller : float - Debye-Waller integral in [eV\ :math:`^{-1}`] - - """ - def __init__(self, debye_waller): - self.debye_waller = debye_waller - - def to_hdf5(self, group): - group.attrs['type'] = np.string_('incoherent_elastic') - group.create_dataset('debye_waller', data=self.debye_waller) - - @classmethod - def from_hdf5(cls, group): - return cls(group['debye_waller']) - - -class IncoherentElasticAEDiscrete(AngleEnergy): - """Discrete angle representation of incoherent elastic scattering - - Parameters - ---------- - mu_out : numpy.ndarray - Equi-probable discrete angles at each incoming energy - - """ - def __init__(self, mu_out): - self.mu_out = mu_out - - def to_hdf5(self, group): - group.attrs['type'] = np.string_('incoherent_elastic_discrete') - group.create_dataset('mu_out', data=self.mu_out) - - @classmethod - def from_hdf5(cls, group): - return cls(group['mu_out'][()]) - - -class IncoherentInelasticAEDiscrete(AngleEnergy): - """Discrete angle representation of incoherent inelastic scattering - - Parameters - ---------- - energy_out : numpy.ndarray - Outgoing energies for each incoming energy - mu_out : numpy.ndarray - Discrete angles for each incoming/outgoing energy - skewed : bool - Whether distance angles are equi-probable or have a skewed distribution - - Attributes - ---------- - energy_out : numpy.ndarray - Outgoing energies for each incoming energy - mu_out : numpy.ndarray - Discrete angles for each incoming/outgoing energy - skewed : bool - Whether distance angles are equi-probable or have a skewed distribution - - """ - def __init__(self, energy_out, mu_out, skewed=False): - self.energy_out = energy_out - self.mu_out = mu_out - self.skewed = skewed - - def to_hdf5(self, group): - 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) - group.create_dataset('skewed', data=self.skewed) - - @classmethod - def from_hdf5(cls, group): - energy_out = group['energy_out'][()] - mu_out = group['mu_out'][()] - skewed = bool(group['skewed']) - return cls(energy_out, mu_out, skewed) - - class ThermalScattering(EqualityMixin): """A ThermalScattering object contains thermal scattering data as represented by an S(alpha, beta) table. @@ -529,13 +405,9 @@ class ThermalScattering(EqualityMixin): else: return "" - @staticmethod - def _temperature_str(T): - return "{}K".format(round(T)) - @property def temperatures(self): - return [self._temperature_str(kT / K_BOLTZMANN) for kT in self.kTs] + return [_temperature_str(kT / K_BOLTZMANN) for kT in self.kTs] def export_to_hdf5(self, path, mode='a', libver='earliest'): """Export table to an HDF5 file. @@ -920,10 +792,10 @@ class ThermalScattering(EqualityMixin): Thermal scattering data """ - if isinstance(ev_or_filename, Evaluation): + if isinstance(ev_or_filename, endf.Evaluation): ev = ev_or_filename else: - ev = Evaluation(ev_or_filename) + ev = endf.Evaluation(ev_or_filename) # Read coherent/incoherent elastic data elastic = None @@ -938,7 +810,7 @@ class ThermalScattering(EqualityMixin): # Get structure factor at first temperature params, S = endf.get_tab1_record(file_obj) - strT = self._temperature_str(params[0]) + strT = _temperature_str(params[0]) n_temps = params[2] bragg_edges = S.x xs[strT] = CoherentElastic(bragg_edges, S.y) @@ -947,16 +819,16 @@ class ThermalScattering(EqualityMixin): # Get structure factor for subsequent temperatures for _ in range(n_temps): params, S = endf.get_list_record(file_obj) - strT = self._temperature_str(params[0]) + strT = _temperature_str(params[0]) xs[strT] = CoherentElastic(bragg_edges, S) - distribution[strT] = CoherentElasticAE(xs[T]) + distribution[strT] = CoherentElasticAE(xs[strT]) elif lhtr == 2: # incoherent elastic params, W = endf.get_tab1_record(file_obj) characteristic_xs = params[0] for T, debye_waller in zip(W.x, W.y): - strT = self._temperature_str(T) + strT = _temperature_str(T) xs[strT] = IncoherentElastic(characteristic_xs, debye_waller) distribution[strT] = IncoherentElasticAE(debye_waller) diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py new file mode 100644 index 0000000000..cf9587cbe8 --- /dev/null +++ b/openmc/data/thermal_angle_energy.py @@ -0,0 +1,136 @@ +import numpy as np + +from .angle_energy import AngleEnergy + + +class CoherentElasticAE(AngleEnergy): + r"""Differential cross section for coherent elastic scattering + + The differential cross section for coherent elastic scattering from a + powdered crystalline material may be represented as: + + .. math:: + \frac{d^2\sigma}{dE'd\Omega} (E\rightarrow E',\mu,T) = \frac{1}{E} \sum + \limits_{i=1}^{E_i < E} s_i(T) \delta(\mu - \mu_i) \delta (E - E') + /(2\pi) + + where :math:`E_i` are the energies of the Bragg edges in [eV], :math:`s_i(T)` + is the structure factor in [eV-b] at the moderator temperature :math:`T` + in [K], and :math:`\mu_i = 1 - 2E_i/E`. + + Parameters + ---------- + coherent_xs : openmc.data.CoherentElastic + Coherent elastic scattering cross section + + Attributes + ---------- + coherent_xs : openmc.data.CoherentElastic + Coherent elastic scattering cross section + + """ + def __init__(self, coherent_xs): + self.coherent_xs = coherent_xs + + def to_hdf5(self, group): + group.attrs['type'] = np.string_('coherent_elastic') + group['coherent_xs'] = group.parent['xs'] + + +class IncoherentElasticAE(AngleEnergy): + r"""Differential cross section for incoherent elastic scattering + + The differential cross section for incoherent elastic scattering may be + represented as: + + .. math:: + \frac{d^2\sigma}{dE'd\Omega} (E\rightarrow E',\mu,T) = \frac{\sigma_b} + {4\pi} e^{-2EW'(T)(1-\mu)} \delta(E - E') + + where :math:`\sigma_b` is the characteristic cross section in [b] and + :math:`W'(T)` is the Debye-Waller integral divided by the atomic mass in + [eV\ :math:`^{-1}`]. + + Parameters + ---------- + debye_waller : float + Debye-Waller integral in [eV\ :math:`^{-1}`] + + Attributes + ---------- + debye_waller : float + Debye-Waller integral in [eV\ :math:`^{-1}`] + + """ + def __init__(self, debye_waller): + self.debye_waller = debye_waller + + def to_hdf5(self, group): + group.attrs['type'] = np.string_('incoherent_elastic') + group.create_dataset('debye_waller', data=self.debye_waller) + + @classmethod + def from_hdf5(cls, group): + return cls(group['debye_waller']) + + +class IncoherentElasticAEDiscrete(AngleEnergy): + """Discrete angle representation of incoherent elastic scattering + + Parameters + ---------- + mu_out : numpy.ndarray + Equi-probable discrete angles at each incoming energy + + """ + def __init__(self, mu_out): + self.mu_out = mu_out + + def to_hdf5(self, group): + group.attrs['type'] = np.string_('incoherent_elastic_discrete') + group.create_dataset('mu_out', data=self.mu_out) + + @classmethod + def from_hdf5(cls, group): + return cls(group['mu_out'][()]) + + +class IncoherentInelasticAEDiscrete(AngleEnergy): + """Discrete angle representation of incoherent inelastic scattering + + Parameters + ---------- + energy_out : numpy.ndarray + Outgoing energies for each incoming energy + mu_out : numpy.ndarray + Discrete angles for each incoming/outgoing energy + skewed : bool + Whether distance angles are equi-probable or have a skewed distribution + + Attributes + ---------- + energy_out : numpy.ndarray + Outgoing energies for each incoming energy + mu_out : numpy.ndarray + Discrete angles for each incoming/outgoing energy + skewed : bool + Whether distance angles are equi-probable or have a skewed distribution + + """ + def __init__(self, energy_out, mu_out, skewed=False): + self.energy_out = energy_out + self.mu_out = mu_out + self.skewed = skewed + + def to_hdf5(self, group): + 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) + group.create_dataset('skewed', data=self.skewed) + + @classmethod + def from_hdf5(cls, group): + energy_out = group['energy_out'][()] + mu_out = group['mu_out'][()] + skewed = bool(group['skewed']) + return cls(energy_out, mu_out, skewed) From c65366d7a0dd83c95b352fd82f5fca1349503785 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 19 Apr 2019 14:35:14 -0500 Subject: [PATCH 06/19] Add IncoherentInelasticAEContinuous class stub --- openmc/data/thermal.py | 6 +++--- openmc/data/thermal_angle_energy.py | 4 ++++ 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index bc2e1d9819..0366af1a66 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -19,12 +19,12 @@ from . import HDF5_VERSION, HDF5_VERSION_MAJOR, endf from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, NATURAL_ABUNDANCE from .ace import Table, get_table, Library from .angle_energy import AngleEnergy -from .correlated import CorrelatedAngleEnergy from .function import Tabulated1D, Function1D from .njoy import make_ace_thermal from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE, IncoherentElasticAEDiscrete, - IncoherentInelasticAEDiscrete) + IncoherentInelasticAEDiscrete, + IncoherentInelasticAEContinuous) _THERMAL_NAMES = { @@ -673,7 +673,7 @@ class ThermalScattering(EqualityMixin): breakpoints = [n_energy] interpolation = [2] energy = inelastic_xs.x - distribution = CorrelatedAngleEnergy( + distribution = IncoherentInelasticAEContinuous( breakpoints, interpolation, energy, energy_out, mu_out) table.inelastic = ThermalScatteringReaction( diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index cf9587cbe8..e5b77b597a 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -1,6 +1,7 @@ import numpy as np from .angle_energy import AngleEnergy +from .correlated import CorrelatedAngleEnergy class CoherentElasticAE(AngleEnergy): @@ -134,3 +135,6 @@ class IncoherentInelasticAEDiscrete(AngleEnergy): mu_out = group['mu_out'][()] skewed = bool(group['skewed']) return cls(energy_out, mu_out, skewed) + +class IncoherentInelasticAEContinuous(CorrelatedAngleEnergy): + pass From 3e4fe616730e62e57e4bda38e5a92f39f1fe9cd8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 17 May 2019 21:13:02 -0500 Subject: [PATCH 07/19] Start adding thermal angle-energy classes on C++ side --- include/openmc/endf.h | 25 +- include/openmc/secondary_correlated.h | 2 +- include/openmc/secondary_kalbach.h | 2 +- include/openmc/secondary_nbody.h | 2 +- include/openmc/secondary_thermal.h | 113 ++++++++ include/openmc/secondary_uncorrelated.h | 2 +- openmc/data/angle_energy.py | 2 + openmc/data/correlated.py | 4 +- openmc/data/thermal.py | 10 +- openmc/data/thermal_angle_energy.py | 3 +- src/endf.cpp | 17 ++ src/secondary_thermal.cpp | 327 ++++++++++++++++++++++++ 12 files changed, 495 insertions(+), 14 deletions(-) create mode 100644 include/openmc/secondary_thermal.h create mode 100644 src/secondary_thermal.cpp diff --git a/include/openmc/endf.h b/include/openmc/endf.h index 7ec078e93f..27958db5c5 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -57,7 +57,7 @@ public: //! Evaluate the polynomials //! \param[in] x independent variable //! \return Polynomial evaluated at x - double operator()(double x) const; + double operator()(double x) const override; private: std::vector coef_; //!< Polynomial coefficients }; @@ -77,7 +77,7 @@ public: //! Evaluate the tabulated function //! \param[in] x independent variable //! \return Function evaluated at x - double operator()(double x) const; + double operator()(double x) const override; // Accessors const std::vector& x() const { return x_; } @@ -96,13 +96,32 @@ private: //============================================================================== class CoherentElasticXS : public Function1D { +public: explicit CoherentElasticXS(hid_t dset); - double operator()(double E) const; + + double operator()(double E) const override; + + const std::vector& bragg_edges() const { return bragg_edges_; } + const std::vector& factors() const { return factors_; } private: std::vector bragg_edges_; //!< Bragg edges in [eV] std::vector factors_; //!< Partial sums of structure factors [eV-b] }; +//============================================================================== +//! Incoherent elastic scattering cross section +//============================================================================== + +class IncoherentElasticXS : public Function1D { +public: + explicit IncoherentElasticXS(hid_t dset); + + double operator()(double E) const override; +private: + double bound_xs_; //!< Characteristic bound xs in [b] + double debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1] +}; + //! Read 1D function from HDF5 dataset //! \param[in] group HDF5 group containing dataset //! \param[in] name Name of dataset diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h index 8ba18d3d15..603d8c414b 100644 --- a/include/openmc/secondary_correlated.h +++ b/include/openmc/secondary_correlated.h @@ -38,7 +38,7 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - void sample(double E_in, double& E_out, double& mu) const; + void sample(double E_in, double& E_out, double& mu) const override; // energy property std::vector& energy() { return energy_; } diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h index a898db67c8..0217c7e776 100644 --- a/include/openmc/secondary_kalbach.h +++ b/include/openmc/secondary_kalbach.h @@ -29,7 +29,7 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - void sample(double E_in, double& E_out, double& mu) const; + void sample(double E_in, double& E_out, double& mu) const override; private: //! Outgoing energy/angle at a single incoming energy struct KMTable { diff --git a/include/openmc/secondary_nbody.h b/include/openmc/secondary_nbody.h index 8509b07143..eb90bd8a04 100644 --- a/include/openmc/secondary_nbody.h +++ b/include/openmc/secondary_nbody.h @@ -24,7 +24,7 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - void sample(double E_in, double& E_out, double& mu) const; + void sample(double E_in, double& E_out, double& mu) const override; private: int n_bodies_; //!< Number of particles distributed double mass_ratio_; //!< Total mass of particles [neutron mass] diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h new file mode 100644 index 0000000000..507176c2e3 --- /dev/null +++ b/include/openmc/secondary_thermal.h @@ -0,0 +1,113 @@ +//! \file secondary_thermal.h +//! Angle-energy distributions for thermal scattering + +#ifndef OPENMC_SECONDARY_THERMAL_H +#define OPENMC_SECONDARY_THERMAL_H + +#include "openmc/angle_energy.h" +#include "openmc/endf.h" +#include "openmc/secondary_correlated.h" + +#include +#include "xtensor/xtensor.hpp" + +#include + +namespace openmc { + +//============================================================================== +//! Coherent elastic scattering angle-energy distribution +//============================================================================== + +class CoherentElasticAE : public AngleEnergy { +public: + //! Construct from a coherent elastic scattering cross section + // + //! \param[in] xs Coherent elastic scattering cross section + explicit CoherentElasticAE(const CoherentElasticXS& xs); + + + //! Sample distribution for an angle and energy + //! \param[in] E_in Incoming energy in [eV] + //! \param[out] E_out Outgoing energy in [eV] + //! \param[out] mu Outgoing cosine with respect to current direction + void sample(double E_in, double& E_out, double& mu) const override; +private: + const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section +}; + +//============================================================================== +//! Incoherent elastic scattering angle-energy distribution +//============================================================================== + +class IncoherentElasticAE : public AngleEnergy { +public: + explicit IncoherentElasticAE(hid_t group); + + //! Sample distribution for an angle and energy + //! \param[in] E_in Incoming energy in [eV] + //! \param[out] E_out Outgoing energy in [eV] + //! \param[out] mu Outgoing cosine with respect to current direction + void sample(double E_in, double& E_out, double& mu) const override; +private: + double debye_waller_; +}; + +//============================================================================== +//! Incoherent elastic scattering angle-energy distribution (discrete) +//============================================================================== + +class IncoherentElasticAEDiscrete : public AngleEnergy { +public: + explicit IncoherentElasticAEDiscrete(hid_t group, const std::vector& energy); + + void sample(double E_in, double& E_out, double& mu) const override; +private: + const std::vector& energy_; //!< Incoherent inelastic scattering cross section + xt::xtensor mu_out_; //!< Cosines for each incident energy +}; + +//============================================================================== +//! Incoherent inelastic scattering angle-energy distribution (discrete) +//============================================================================== + +class IncoherentInelasticAEDiscrete : public AngleEnergy { +public: + explicit IncoherentInelasticAEDiscrete(hid_t group, const std::vector& energy); + + void sample(double E_in, double& E_out, double& mu) const override; +private: + const std::vector& energy_; //!< Incident energies + xt::xtensor energy_out_; //!< Outgoing energies for each incident energy + xt::xtensor mu_out_; //!< Outgoing cosines for each incident/outgoing energy + bool skewed_; //!< Whether outgoing energy distribution is skewed +}; + +//============================================================================== +//! Incoherent inelastic scattering angle-energy distribution (discrete) +//============================================================================== + +class IncoherentInelasticAEContinuous : public AngleEnergy { +public: + explicit IncoherentInelasticAEContinuous(hid_t group); + + void sample(double E_in, double& E_out, double& mu) const override; +private: + //! Secondary energy/angle distribution + struct DistEnergySab { + std::size_t n_e_out; //!< Number of outgoing energies + xt::xtensor e_out; //!< Outgoing energies + xt::xtensor e_out_pdf; //!< Probability density function + xt::xtensor e_out_cdf; //!< Cumulative distribution function + xt::xtensor mu; //!< Equiprobable angles at each outgoing energy + }; + + std::vector energy_; //!< Incident energies + std::vector distribution_; //!< Secondary angle-energy at + //!< each incident energy +}; + + +} // namespace openmc + +#endif // OPENMC_SECONDARY_THERMAL_H diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h index e895a17bdc..e430c75c45 100644 --- a/include/openmc/secondary_uncorrelated.h +++ b/include/openmc/secondary_uncorrelated.h @@ -29,7 +29,7 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - void sample(double E_in, double& E_out, double& mu) const; + void sample(double E_in, double& E_out, double& mu) const override; // Accessors AngleDistribution& angle() { return angle_; } diff --git a/openmc/data/angle_energy.py b/openmc/data/angle_energy.py index 19dddfe19e..5fb176f681 100644 --- a/openmc/data/angle_energy.py +++ b/openmc/data/angle_energy.py @@ -43,6 +43,8 @@ class AngleEnergy(EqualityMixin, metaclass=ABCMeta): return openmc.data.IncoherentElasticAEDiscrete.from_hdf5(group) elif dist_type == 'incoherent_inelastic_discrete': return openmc.data.IncoherentInelasticAEDiscrete.from_hdf5(group) + elif dist_type == 'incoherent_inelastic_continuous': + return openmc.data.IncoherentInelasticAEContinuous.from_hdf5(group) @staticmethod def from_ace(ace, location_dist, location_start, rx=None): diff --git a/openmc/data/correlated.py b/openmc/data/correlated.py index b760f2e2fd..567d23abd4 100644 --- a/openmc/data/correlated.py +++ b/openmc/data/correlated.py @@ -44,6 +44,8 @@ class CorrelatedAngleEnergy(AngleEnergy): """ + _name = 'correlated' + def __init__(self, breakpoints, interpolation, energy, energy_out, mu): super().__init__() self.breakpoints = breakpoints @@ -111,7 +113,7 @@ class CorrelatedAngleEnergy(AngleEnergy): HDF5 group to write to """ - group.attrs['type'] = np.string_('correlated') + group.attrs['type'] = np.string_(self._name) dset = group.create_dataset('energy', data=self.energy) dset.attrs['interpolation'] = np.vstack((self.breakpoints, diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 0366af1a66..eaf8b49bce 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -248,14 +248,14 @@ class IncoherentElastic(Function1D): Parameters ---------- xs : float - Characteristic cross section in [b] + Characteristic bound cross section in [b] debye_waller : float Debye-Waller integral in [eV\ :math:`^{-1}`] Attributes ---------- xs : float - Characteristic cross section in [b] + Characteristic bound cross section in [b] debye_waller : float Debye-Waller integral in [eV\ :math:`^{-1}`] @@ -282,7 +282,7 @@ class IncoherentElastic(Function1D): dataset = group.create_dataset(name) dataset.attrs['type'] = np.string_('incoherent') dataset.attrs['debye_waller'] = self.debye_waller - dataset.attrs['characteristic_xs'] = self.xs + dataset.attrs['bound_xs'] = self.xs @classmethod def from_hdf5(cls, dataset): @@ -826,10 +826,10 @@ class ThermalScattering(EqualityMixin): elif lhtr == 2: # incoherent elastic params, W = endf.get_tab1_record(file_obj) - characteristic_xs = params[0] + bound_xs = params[0] for T, debye_waller in zip(W.x, W.y): strT = _temperature_str(T) - xs[strT] = IncoherentElastic(characteristic_xs, debye_waller) + xs[strT] = IncoherentElastic(bound_xs, debye_waller) distribution[strT] = IncoherentElasticAE(debye_waller) elastic = ThermalScatteringReaction(xs, distribution) diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index e5b77b597a..d810044768 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -136,5 +136,6 @@ class IncoherentInelasticAEDiscrete(AngleEnergy): skewed = bool(group['skewed']) return cls(energy_out, mu_out, skewed) + class IncoherentInelasticAEContinuous(CorrelatedAngleEnergy): - pass + _name = 'incoherent_inelastic_continuous' diff --git a/src/endf.cpp b/src/endf.cpp index 464bc974c2..b10bbd5502 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -234,4 +234,21 @@ double CoherentElasticXS::operator()(double E) const } } +//============================================================================== +// IncoherentElasticXS implementation +//============================================================================== + +IncoherentElasticXS::IncoherentElasticXS(hid_t dset) +{ + read_attribute(dset, "bound_xs", bound_xs_); + read_attribute(dset, "debye_waller", debye_waller_); +} + +double IncoherentElasticXS::operator()(double E) const +{ + // Determine cross section using ENDF-102, Eq. (7.5) + double W = debye_waller_; + return bound_xs_ / 2.0 * ((1 - std::exp(-4.0*E*W))/(2.0*E*W)); +} + } // namespace openmc diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp new file mode 100644 index 0000000000..cc2fce961c --- /dev/null +++ b/src/secondary_thermal.cpp @@ -0,0 +1,327 @@ +#include "openmc/secondary_thermal.h" + +#include "openmc/hdf5_interface.h" +#include "openmc/random_lcg.h" +#include "openmc/search.h" + +#include "xtensor/xview.hpp" + +#include // for log, exp + +namespace openmc { + +// Helper function to get index on incident energy grid +void +get_energy_index(const std::vector& energies, double E, int& i, double& f) +{ + // Get index and interpolation factor for elastic grid + i = 0; + f = 0.0; + if (E >= energies.front()) { + i = lower_bound_index(energies.begin(), energies.end(), E); + f = (E - energies[i]) / (energies[i+1] - energies[i]); + } +} + +//============================================================================== +// CoherentElasticAE implementation +//============================================================================== + +CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs) + : xs_{xs} +{ } + +void +CoherentElasticAE::sample(double E_in, double& E_out, double& mu) const +{ + // Get index and interpolation factor for elastic grid + int i; + double f; + const auto& energies {xs_.bragg_edges()}; + get_energy_index(energies, E_in, i, f); + + // Sample a Bragg edge between 1 and i + const auto& factors = xs_.factors(); + double prob = prn() * factors[i+1]; + int k = 0; + if (prob >= factors.front()) { + k = lower_bound_index(factors.begin(), factors.begin() + (i+1), prob); + } + + // Characteristic scattering cosine for this Bragg edge (ENDF-102, Eq. 7-2) + mu = 1.0 - 2.0*energies[k] / E_in; + + // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1) + E_out = E_in; +} + +//============================================================================== +// IncoherentElasticAE implementation +//============================================================================== + +IncoherentElasticAE::IncoherentElasticAE(hid_t group) +{ + read_attribute(group, "debye_waller", debye_waller_); +} + +void +IncoherentElasticAE::sample(double E_in, double& E_out, double& mu) const +{ + // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4 + double c = 2 * E_in * debye_waller_; + mu = std::log(1.0 + prn()*(std::exp(2.0*c) - 1))/c - 1.0; + + // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4) + E_out = E_in; +} + +//============================================================================== +// IncoherentElasticAEDiscrete implementation +//============================================================================== + +IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(hid_t group, + const std::vector& energy) + : energy_{energy} +{ + read_dataset(group, "mu_out", mu_out_); +} + +void +IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu) const +{ + // Get index and interpolation factor for elastic grid + int i; + double f; + get_energy_index(energy_, E_in, i, f); + + // Interpolate between two discrete cosines corresponding to neighboring + // incoming energies. + + // Sample outgoing cosine bin + int k = prn() * mu_out_.shape()[1]; + + // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1] + double mu_ik = mu_out_(i, k); + double mu_i1k = mu_out_(i+1, k); + + // Cosine of angle between incoming and outgoing neutron + mu = (1 - f)*mu_ik + f*mu_i1k; + + // Energy doesn't change in elastic scattering + E_out = E_in; +} + +//============================================================================== +// IncoherentInelasticAEDiscrete implementation +//============================================================================== + +IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(hid_t group, + const std::vector& energy) + : energy_{energy} +{ + read_dataset(group, "energy_out", energy_out_); + read_dataset(group, "mu_out", mu_out_); + read_dataset(group, "skewed", skewed_); +} + +void +IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu) const +{ + // Get index and interpolation factor for inelastic grid + int i; + double f; + get_energy_index(energy_, E_in, i, f); + + // Now that we have an incoming energy bin, we need to determine the outgoing + // energy bin. This will depend on whether the outgoing energy distribution is + // skewed. If it is skewed, then the first two and last two bins have lower + // probabilities than the other bins (0.1 for the first and last bins and 0.4 + // for the second and second to last bins, relative to a normal bin + // probability of 1). Otherwise, each bin is equally probable. + + int j; + int n = energy_out_.shape()[1]; + if (!skewed_) { + // All bins equally likely + j = prn() * n; + } else { + // Distribution skewed away from edge points + double r = prn() * (n - 3); + if (r > 1.0) { + // equally likely N-4 middle bins + j = r + 1; + } else if (r > 0.6) { + // second to last bin has relative probability of 0.4 + j = n - 2; + } else if (r > 0.5) { + // last bin has relative probability of 0.1 + j = n - 1; + } else if (r > 0.1) { + // second bin has relative probability of 0.4 + j = 1; + } else { + // first bin has relative probability of 0.1 + j = 0; + } + } + + // Determine outgoing energy corresponding to E_in[i] and E_in[i+1] + double E_ij = energy_out_(i, j); + double E_i1j = energy_out_(i+1, j); + + // Outgoing energy + E_out = (1 - f)*E_ij + f*E_i1j; + + // Sample outgoing cosine bin + int m = mu_out_.shape()[2]; + int k = prn() * m; + + // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1] + double mu_ijk = mu_out_(i, j, k); + double mu_i1jk = mu_out_(i+1, j, k); + + // Cosine of angle between incoming and outgoing neutron + mu = (1 - f)*mu_ijk + f*mu_i1jk; +} + +//============================================================================== +// IncoherentInelasticAEContinuous implementation +//============================================================================== + +IncoherentInelasticAEContinuous::IncoherentInelasticAEContinuous(hid_t group) +{ + // Read correlated angle-energy distribution + CorrelatedAngleEnergy dist {group}; + + // Copy incident energies + energy_ = dist.energy(); + + // Convert to S(a,b) native format + for (const auto& edist : dist.distribution()) { + // Create temporary distribution + DistEnergySab d; + + // Copy outgoing energy distribution + d.n_e_out = edist.e_out.size(); + d.e_out = edist.e_out; + d.e_out_pdf = edist.p; + d.e_out_cdf = edist.c; + + for (int j = 0; j < d.n_e_out; ++j) { + auto adist = dynamic_cast(edist.angle[j].get()); + if (adist) { + // On first pass, allocate space for angles + if (j == 0) { + auto n_mu = adist->x().size(); + d.mu = xt::empty({d.n_e_out, n_mu}); + } + + // Copy outgoing angles + auto mu_j = xt::view(d.mu, j); + std::copy(adist->x().begin(), adist->x().end(), mu_j.begin()); + } + } + + distribution_.emplace_back(std::move(d)); + } + +} + +void +IncoherentInelasticAEContinuous::sample(double E_in, double& E_out, double& mu) const +{ + // Get index and interpolation factor for inelastic grid + int i; + double f; + get_energy_index(energy_, E_in, i, f); + + // Sample between ith and [i+1]th bin + int l = f > prn() ? i + 1 : i; + + // Determine endpoints on grid i + auto n = distribution_[i].e_out.size(); + double E_i_1 = distribution_[i].e_out(0); + double E_i_J = distribution_[i].e_out(n - 1); + + // Determine endpoints on grid i + 1 + n = distribution_[i + 1].e_out.size(); + double E_i1_1 = distribution_[i + 1].e_out(0); + double E_i1_J = distribution_[i + 1].e_out(n - 1); + + double E_1 = E_i_1 + f * (E_i1_1 - E_i_1); + double E_J = E_i_J + f * (E_i1_J - E_i_J); + + // Determine outgoing energy bin + // (First reset n_energy_out to the right value) + n = distribution_[l].n_e_out; + double r1 = prn(); + double c_j = distribution_[l].e_out_cdf[0]; + double c_j1; + std::size_t j; + for (j = 0; j < n - 1; ++j) { + c_j1 = distribution_[l].e_out_cdf[j + 1]; + if (r1 < c_j1) break; + c_j = c_j1; + } + + // check to make sure j is <= n_energy_out - 2 + j = std::min(j, n - 2); + + // Get the data to interpolate between + double E_l_j = distribution_[l].e_out[j]; + double p_l_j = distribution_[l].e_out_pdf[j]; + + // Next part assumes linear-linear interpolation in standard + double E_l_j1 = distribution_[l].e_out[j + 1]; + double p_l_j1 = distribution_[l].e_out_pdf[j + 1]; + + // Find secondary energy (variable E) + double frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j); + if (frac == 0.0) { + E_out = E_l_j + (r1 - c_j) / p_l_j; + } else { + E_out = E_l_j + (std::sqrt(std::max(0.0, p_l_j*p_l_j + + 2.0*frac*(r1 - c_j))) - p_l_j) / frac; + } + + // Now interpolate between incident energy bins i and i + 1 + if (l == i) { + E_out = E_1 + (E_out - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1); + } else { + E_out = E_1 + (E_out - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1); + } + + // Sample outgoing cosine bin + int n_mu = distribution_[l].mu.shape()[1]; + std::size_t k = prn() * n_mu; + + // Rather than use the sampled discrete mu directly, it is smeared over + // a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the + // discrete mu value itself. + const auto& mu_l = distribution_[l].mu; + f = (r1 - c_j)/(c_j1 - c_j); + + // Determine (k-1)th mu value + double mu_left; + if (k == 0) { + mu_left = -1.0; + } else { + mu_left = mu_l(j, k-1) + f*(mu_l(j+1, k-1) - mu_l(j, k-1)); + } + + // Determine kth mu value + mu = mu_l(j, k) + f*(mu_l(j+1, k) - mu_l(j, k)); + + // Determine (k+1)th mu value + double mu_right; + if (k == n_mu - 1) { + mu_right = 1.0; + } else { + mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1)); + } + + // Smear angle + mu += std::min(mu - mu_left, mu_right - mu)*(prn() - 0.5); +} + +} // namespace openmc From fd1f01a9110c2a07136d311f3b7924878d3be376 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 21 May 2019 13:56:44 -0500 Subject: [PATCH 08/19] Move thermal scattering reaction data into ThermalData::Reaction type --- include/openmc/secondary_thermal.h | 6 +- include/openmc/thermal.h | 45 +--- openmc/data/thermal.py | 2 +- openmc/data/thermal_angle_energy.py | 2 +- src/endf.cpp | 4 + src/reaction_product.cpp | 10 +- src/secondary_thermal.cpp | 6 +- src/thermal.cpp | 348 ++++------------------------ 8 files changed, 69 insertions(+), 354 deletions(-) diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 507176c2e3..998143d0f5 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -84,12 +84,12 @@ private: }; //============================================================================== -//! Incoherent inelastic scattering angle-energy distribution (discrete) +//! Incoherent inelastic scattering angle-energy distribution //============================================================================== -class IncoherentInelasticAEContinuous : public AngleEnergy { +class IncoherentInelasticAE : public AngleEnergy { public: - explicit IncoherentInelasticAEContinuous(hid_t group); + explicit IncoherentInelasticAE(hid_t group); void sample(double E_in, double& E_out, double& mu) const override; private: diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h index b737f005b3..9f8ebaf733 100644 --- a/include/openmc/thermal.h +++ b/include/openmc/thermal.h @@ -9,6 +9,8 @@ #include "xtensor/xtensor.hpp" +#include "openmc/angle_energy.h" +#include "openmc/endf.h" #include "openmc/hdf5_interface.h" #include "openmc/particle.h" @@ -47,20 +49,19 @@ extern std::unordered_map thermal_scatt_map; class ThermalData { public: - ThermalData(hid_t group, int secondary_mode); + ThermalData(hid_t group); // Sample an outgoing energy and angle void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out, double* mu); private: - //! Secondary energy/angle distributions for inelastic thermal scattering - //! collisions which utilize a continuous secondary energy representation. - struct DistEnergySab { - std::size_t n_e_out; //!< Number of outgoing energies - xt::xtensor e_out; //!< Outgoing energies - xt::xtensor e_out_pdf; //!< Probability density function - xt::xtensor e_out_cdf; //!< Cumulative distribution function - xt::xtensor mu; //!< Equiprobable angles at each outgoing energy + struct Reaction { + // Default constructor + Reaction() { } + + // Data members + std::unique_ptr xs; //!< Cross section + std::unique_ptr distribution; //!< Secondary angle-energy distribution }; //! Upper threshold for incoherent inelastic scattering (usually ~4 eV) @@ -69,30 +70,8 @@ private: double threshold_elastic_ {0.0}; // Inelastic scattering data - int inelastic_mode_; //!< distribution type (equal/skewed/continuous) - std::size_t n_inelastic_e_in_; //!< number of incoming E for inelastic - std::size_t n_inelastic_e_out_; //!< number of outgoing E for inelastic - std::size_t n_inelastic_mu_; //!< number of outgoing angles for inelastic - std::vector inelastic_e_in_; //!< incoming E grid for inelastic - std::vector inelastic_sigma_; //!< inelastic scattering cross section - - // The following are used only for equal/skewed distributions - xt::xtensor inelastic_e_out_; - xt::xtensor inelastic_mu_; - - // The following is used only for continuous S(a,b) distributions. The - // different implementation is necessary because the continuous representation - // has a variable number of outgoing energy points for each incoming energy - std::vector inelastic_data_; //!< Secondary angle-energy at - //!< each incoming energy - - // Elastic scattering data - int elastic_mode_; //!< type of elastic (incoherent/coherent) - std::size_t n_elastic_e_in_; //!< number of incoming E for elastic - std::size_t n_elastic_mu_; //!< number of outgoing angles for elastic - std::vector elastic_e_in_; //!< incoming E grid for elastic - std::vector elastic_P_; //!< elastic scattering cross section - xt::xtensor elastic_mu_; //!< equi-probable angles at each incoming E + Reaction elastic_; + Reaction inelastic_; // ThermalScattering needs access to private data members friend class ThermalScattering; diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index eaf8b49bce..cdb5890930 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -280,7 +280,7 @@ class IncoherentElastic(Function1D): """ dataset = group.create_dataset(name) - dataset.attrs['type'] = np.string_('incoherent') + dataset.attrs['type'] = np.string_(type(self).__name__) dataset.attrs['debye_waller'] = self.debye_waller dataset.attrs['bound_xs'] = self.xs diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index d810044768..d6053ed2a4 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -138,4 +138,4 @@ class IncoherentInelasticAEDiscrete(AngleEnergy): class IncoherentInelasticAEContinuous(CorrelatedAngleEnergy): - _name = 'incoherent_inelastic_continuous' + _name = 'incoherent_inelastic' diff --git a/src/endf.cpp b/src/endf.cpp index b10bbd5502..23512c4dbf 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -87,6 +87,10 @@ read_function(hid_t group, const char* name) func = std::make_unique(dset); } else if (func_type == "Polynomial") { func = std::make_unique(dset); + } else if (func_type == "CoherentElastic") { + func = std::make_unique(dset); + } else if (func_type == "IncoherentElastic") { + func = std::make_unique(dset); } else { throw std::runtime_error{"Unknown function type " + func_type + " for dataset " + object_name(dset)}; diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index 4f3994ff2f..34680a7735 100644 --- a/src/reaction_product.cpp +++ b/src/reaction_product.cpp @@ -3,6 +3,7 @@ #include // for unique_ptr #include // for string +#include "openmc/endf.h" #include "openmc/hdf5_interface.h" #include "openmc/random_lcg.h" #include "openmc/secondary_correlated.h" @@ -42,14 +43,7 @@ ReactionProduct::ReactionProduct(hid_t group) read_attribute(group, "decay_rate", decay_rate_); // Read secondary particle yield - hid_t yield = open_dataset(group, "yield"); - read_attribute(yield, "type", temp); - if (temp == "Tabulated1D") { - yield_ = std::unique_ptr{new Tabulated1D{yield}}; - } else if (temp == "Polynomial") { - yield_ = std::unique_ptr{new Polynomial{yield}}; - } - close_dataset(yield); + yield_ = read_function(group, "yield"); int n; read_attribute(group, "n_distribution", n); diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index cc2fce961c..5fdf6201a0 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -185,10 +185,10 @@ IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu) co } //============================================================================== -// IncoherentInelasticAEContinuous implementation +// IncoherentInelasticAE implementation //============================================================================== -IncoherentInelasticAEContinuous::IncoherentInelasticAEContinuous(hid_t group) +IncoherentInelasticAE::IncoherentInelasticAE(hid_t group) { // Read correlated angle-energy distribution CorrelatedAngleEnergy dist {group}; @@ -228,7 +228,7 @@ IncoherentInelasticAEContinuous::IncoherentInelasticAEContinuous(hid_t group) } void -IncoherentInelasticAEContinuous::sample(double E_in, double& E_out, double& mu) const +IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) const { // Get index and interpolation factor for inelastic grid int i; diff --git a/src/thermal.cpp b/src/thermal.cpp index 40bb19e894..cd7a656a5d 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -12,10 +12,12 @@ #include "xtensor/xview.hpp" #include "openmc/constants.h" +#include "openmc/endf.h" #include "openmc/error.h" #include "openmc/random_lcg.h" #include "openmc/search.h" #include "openmc/secondary_correlated.h" +#include "openmc/secondary_thermal.h" #include "openmc/settings.h" namespace openmc { @@ -43,16 +45,6 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector& tem read_attribute(group, "atomic_weight_ratio", awr_); read_attribute(group, "nuclides", nuclides_); - std::string sec_mode; - read_attribute(group, "secondary_mode", sec_mode); - int secondary_mode; - if (sec_mode == "equal") { - secondary_mode = SAB_SECONDARY_EQUAL; - } else if (sec_mode == "skewed") { - secondary_mode = SAB_SECONDARY_SKEWED; - } else if (sec_mode == "continuous") { - secondary_mode = SAB_SECONDARY_CONT; - } // Read temperatures hid_t kT_group = open_group(group, "kTs"); @@ -148,7 +140,7 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector& tem // Open group for temperature i hid_t T_group = open_group(group, temp_str.data()); - data_.emplace_back(T_group, secondary_mode); + data_.emplace_back(T_group); close_group(T_group); } @@ -250,46 +242,39 @@ ThermalScattering::has_nuclide(const char* name) const // ThermalData implementation //============================================================================== -ThermalData::ThermalData(hid_t group, int secondary_mode) - : inelastic_mode_{secondary_mode} +ThermalData::ThermalData(hid_t group) { - // Coherent elastic data + // Coherent/incoherent elastic data if (object_exists(group, "elastic")) { // Read cross section data hid_t elastic_group = open_group(group, "elastic"); // Read elastic cross section - xt::xarray temp; - hid_t dset = open_dataset(elastic_group, "xs"); - read_dataset(dset, temp); + elastic_.xs = read_function(elastic_group, "xs"); - // Get view on energies and cross section/probability values - auto E_in = xt::view(temp, 0); - auto P = xt::view(temp, 1); + // Read angle-energy distribution + hid_t dgroup = open_group(elastic_group, "distribution"); + std::string temp; + read_attribute(dgroup, "type", temp); + if (temp == "coherent_elastic") { + auto xs = dynamic_cast(elastic_.xs.get()); + elastic_.distribution = std::make_unique(*xs); - // Set cross section data and type - std::copy(E_in.begin(), E_in.end(), std::back_inserter(elastic_e_in_)); - std::copy(P.begin(), P.end(), std::back_inserter(elastic_P_)); - n_elastic_e_in_ = elastic_e_in_.size(); + // Set threshold energy + threshold_elastic_ = xs->bragg_edges().back(); - // Determine whether elastic scattering is incoherent or coherent - std::string type; - read_attribute(dset, "type", type); - if (type == "tab1") { - elastic_mode_ = SAB_ELASTIC_INCOHERENT; - } else if (type == "bragg") { - elastic_mode_ = SAB_ELASTIC_COHERENT; - } - close_dataset(dset); + } else { + auto xs = dynamic_cast(elastic_.xs.get()); + if (temp == "incoherent_elastic") { + elastic_.distribution = std::make_unique(dgroup); + } else if (temp == "incoherent_elastic_discrete") { + elastic_.distribution = std::make_unique( + dgroup, xs->x() + ); + } - // Set elastic threshold - threshold_elastic_ = elastic_e_in_.back(); - - // Read angle distribution - if (elastic_mode_ == SAB_ELASTIC_INCOHERENT) { - xt::xarray mu_out; - read_dataset(elastic_group, "mu_out", mu_out); - elastic_mu_ = mu_out; + // Set threshold energy + threshold_elastic_ = xs->x().back(); } close_group(elastic_group); @@ -300,67 +285,23 @@ ThermalData::ThermalData(hid_t group, int secondary_mode) // Read type of inelastic data hid_t inelastic_group = open_group(group, "inelastic"); - // Read cross section data - xt::xarray temp; - read_dataset(inelastic_group, "xs", temp); - - // Get view of inelastic cross section and energy grid - auto E_in = xt::view(temp, 0); - auto xs = xt::view(temp, 1); - - // Set cross section data - std::copy(E_in.begin(), E_in.end(), std::back_inserter(inelastic_e_in_)); - std::copy(xs.begin(), xs.end(), std::back_inserter(inelastic_sigma_)); - n_inelastic_e_in_ = inelastic_e_in_.size(); + // Read inelastic cross section + inelastic_.xs = read_function(inelastic_group, "xs"); // Set inelastic threshold - threshold_inelastic_ = inelastic_e_in_.back(); + auto xs = dynamic_cast(inelastic_.xs.get()); + threshold_inelastic_ = xs->x().back(); - if (secondary_mode != SAB_SECONDARY_CONT) { - // Read energy distribution - xt::xarray E_out; - read_dataset(inelastic_group, "energy_out", E_out); - inelastic_e_out_ = E_out; - n_inelastic_e_out_ = inelastic_e_out_.shape()[1]; - - // Read angle distribution - xt::xarray mu_out; - read_dataset(inelastic_group, "mu_out", mu_out); - inelastic_mu_ = mu_out; - n_inelastic_mu_ = inelastic_mu_.shape()[2]; - } else { - // Read correlated angle-energy distribution - CorrelatedAngleEnergy dist {inelastic_group}; - - // Convert to S(a,b) native format - for (const auto& edist : dist.distribution()) { - // Create temporary distribution - DistEnergySab d; - - // Copy outgoing energy distribution - d.n_e_out = edist.e_out.size(); - d.e_out = edist.e_out; - d.e_out_pdf = edist.p; - d.e_out_cdf = edist.c; - - for (int j = 0; j < d.n_e_out; ++j) { - auto adist = dynamic_cast(edist.angle[j].get()); - if (adist) { - // On first pass, allocate space for angles - if (j == 0) { - auto n_mu = adist->x().size(); - n_inelastic_mu_ = n_mu; - d.mu = xt::empty({d.n_e_out, n_mu}); - } - - // Copy outgoing angles - auto mu_j = xt::view(d.mu, j); - std::copy(adist->x().begin(), adist->x().end(), mu_j.begin()); - } - } - - inelastic_data_.push_back(std::move(d)); - } + // Read angle-energy distribution + hid_t dgroup = open_group(inelastic_group, "distribution"); + std::string temp; + read_attribute(dgroup, "type", temp); + if (temp == "incoherent_inelastic") { + inelastic_.distribution = std::make_unique(dgroup); + } else if (temp == "incoherent_inelastic_discrete") { + inelastic_.distribution = std::make_unique( + dgroup, xs->x() + ); } close_group(inelastic_group); @@ -373,218 +314,15 @@ ThermalData::sample(const NuclideMicroXS& micro_xs, double E, { // Determine whether inelastic or elastic scattering will occur if (prn() < micro_xs.thermal_elastic / micro_xs.thermal) { - // elastic scattering - - // Get index and interpolation factor for elastic grid - int i = 0; - double f = 0.0; - if (E >= elastic_e_in_.front()) { - auto& E_in = elastic_e_in_; - i = lower_bound_index(E_in.begin(), E_in.end(), E); - f = (E - E_in[i]) / (E_in[i+1] - E_in[i]); - } - - // Select treatment based on elastic mode - if (elastic_mode_ == SAB_ELASTIC_INCOHERENT) { - // With this treatment, we interpolate between two discrete cosines - // corresponding to neighboring incoming energies. This is used for - // data derived in the incoherent approximation - - // Sample outgoing cosine bin - int k = prn() * n_elastic_mu_; - - // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1] - double mu_ijk = elastic_mu_(i, k); - double mu_i1jk = elastic_mu_(i+1, k); - - // Cosine of angle between incoming and outgoing neutron - *mu = (1 - f)*mu_ijk + f*mu_i1jk; - - } else if (elastic_mode_ == SAB_ELASTIC_COHERENT) { - // This treatment is used for data derived in the coherent - // approximation, i.e. for crystalline structures that have Bragg - // edges. - - // Sample a Bragg edge between 1 and i - double prob = prn() * elastic_P_[i+1]; - int k = 0; - if (prob >= elastic_P_.front()) { - k = lower_bound_index(elastic_P_.begin(), elastic_P_.begin() + (i+1), prob); - } - - // Characteristic scattering cosine for this Bragg edge - *mu = 1.0 - 2.0*elastic_e_in_[k] / E; - - } - - // Outgoing energy is same as incoming energy - *E_out = E; - + elastic_.distribution->sample(E, *E_out, *mu); } else { - // Perform inelastic calculations - - // Get index and interpolation factor for inelastic grid - int i = 0; - double f = 0.0; - if (E >= inelastic_e_in_.front()) { - auto& E_in = inelastic_e_in_; - i = lower_bound_index(E_in.begin(), E_in.end(), E); - f = (E - E_in[i]) / (E_in[i+1] - E_in[i]); - } - - // Now that we have an incoming energy bin, we need to determine the - // outgoing energy bin. This will depend on the "secondary energy - // mode". If the mode is 0, then the outgoing energy bin is chosen from a - // set of equally-likely bins. If the mode is 1, then the first - // two and last two bins are skewed to have lower probabilities than the - // other bins (0.1 for the first and last bins and 0.4 for the second and - // second to last bins, relative to a normal bin probability of 1). - // Finally, if the mode is 2, then a continuous distribution (with - // accompanying PDF and CDF is utilized) - - if (inelastic_mode_ == SAB_SECONDARY_EQUAL || - inelastic_mode_ == SAB_SECONDARY_SKEWED) { - int j; - if (inelastic_mode_ == SAB_SECONDARY_EQUAL) { - // All bins equally likely - j = prn() * n_inelastic_e_out_; - } else if (inelastic_mode_ == SAB_SECONDARY_SKEWED) { - // Distribution skewed away from edge points - double r = prn() * (n_inelastic_e_out_ - 3); - if (r > 1.0) { - // equally likely N-4 middle bins - j = r + 1; - } else if (r > 0.6) { - // second to last bin has relative probability of 0.4 - j = n_inelastic_e_out_ - 2; - } else if (r > 0.5) { - // last bin has relative probability of 0.1 - j = n_inelastic_e_out_ - 1; - } else if (r > 0.1) { - // second bin has relative probability of 0.4 - j = 1; - } else { - // first bin has relative probability of 0.1 - j = 0; - } - } - - // Determine outgoing energy corresponding to E_in[i] and E_in[i+1] - double E_ij = inelastic_e_out_(i, j); - double E_i1j = inelastic_e_out_(i+1, j); - - // Outgoing energy - *E_out = (1 - f)*E_ij + f*E_i1j; - - // Sample outgoing cosine bin - int k = prn() * n_inelastic_mu_; - - // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1] - double mu_ijk = inelastic_mu_(i, j, k); - double mu_i1jk = inelastic_mu_(i+1, j, k); - - // Cosine of angle between incoming and outgoing neutron - *mu = (1 - f)*mu_ijk + f*mu_i1jk; - - } else if (inelastic_mode_ == SAB_SECONDARY_CONT) { - // Continuous secondary energy - this is to be similar to - // Law 61 interpolation on outgoing energy - - // Sample between ith and [i+1]th bin - int l = f > prn() ? i + 1 : i; - - // Determine endpoints on grid i - auto n = inelastic_data_[i].e_out.size(); - double E_i_1 = inelastic_data_[i].e_out(0); - double E_i_J = inelastic_data_[i].e_out(n - 1); - - // Determine endpoints on grid i + 1 - n = inelastic_data_[i + 1].e_out.size(); - double E_i1_1 = inelastic_data_[i + 1].e_out(0); - double E_i1_J = inelastic_data_[i + 1].e_out(n - 1); - - double E_1 = E_i_1 + f * (E_i1_1 - E_i_1); - double E_J = E_i_J + f * (E_i1_J - E_i_J); - - // Determine outgoing energy bin - // (First reset n_energy_out to the right value) - n = inelastic_data_[l].n_e_out; - double r1 = prn(); - double c_j = inelastic_data_[l].e_out_cdf[0]; - double c_j1; - std::size_t j; - for (j = 0; j < n - 1; ++j) { - c_j1 = inelastic_data_[l].e_out_cdf[j + 1]; - if (r1 < c_j1) break; - c_j = c_j1; - } - - // check to make sure j is <= n_energy_out - 2 - j = std::min(j, n - 2); - - // Get the data to interpolate between - double E_l_j = inelastic_data_[l].e_out[j]; - double p_l_j = inelastic_data_[l].e_out_pdf[j]; - - // Next part assumes linear-linear interpolation in standard - double E_l_j1 = inelastic_data_[l].e_out[j + 1]; - double p_l_j1 = inelastic_data_[l].e_out_pdf[j + 1]; - - // Find secondary energy (variable E) - double frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j); - if (frac == 0.0) { - *E_out = E_l_j + (r1 - c_j) / p_l_j; - } else { - *E_out = E_l_j + (std::sqrt(std::max(0.0, p_l_j*p_l_j + - 2.0*frac*(r1 - c_j))) - p_l_j) / frac; - } - - // Now interpolate between incident energy bins i and i + 1 - if (l == i) { - *E_out = E_1 + (*E_out - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1); - } else { - *E_out = E_1 + (*E_out - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1); - } - - // Sample outgoing cosine bin - std::size_t k = prn() * n_inelastic_mu_; - - // Rather than use the sampled discrete mu directly, it is smeared over - // a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the - // discrete mu value itself. - const auto& mu_l = inelastic_data_[l].mu; - f = (r1 - c_j)/(c_j1 - c_j); - - // Determine (k-1)th mu value - double mu_left; - if (k == 0) { - mu_left = -1.0; - } else { - mu_left = mu_l(j, k-1) + f*(mu_l(j+1, k-1) - mu_l(j, k-1)); - } - - // Determine kth mu value - *mu = mu_l(j, k) + f*(mu_l(j+1, k) - mu_l(j, k)); - - // Determine (k+1)th mu value - double mu_right; - if (k == n_inelastic_mu_ - 1) { - mu_right = 1.0; - } else { - mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1)); - } - - // Smear angle - *mu += std::min(*mu - mu_left, mu_right - *mu)*(prn() - 0.5); - - } // (inelastic secondary energy treatment) - } // (elastic or inelastic) + inelastic_.distribution->sample(E, *E_out, *mu); + } // Because of floating-point roundoff, it may be possible for mu to be // outside of the range [-1,1). In these cases, we just set mu to exactly // -1 or 1 if (std::abs(*mu) > 1.0) *mu = std::copysign(1.0, *mu); - } void free_memory_thermal() From c851030511a9ad43c5100c3be4bd1761db9de30b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 21 May 2019 14:44:15 -0500 Subject: [PATCH 09/19] Refactor ThermalScattering::calculate_xs and add energy_max to class --- CMakeLists.txt | 1 + include/openmc/thermal.h | 17 ++++++++-- openmc/data/thermal.py | 19 +++++++---- src/material.cpp | 2 +- src/thermal.cpp | 69 ++++++++++------------------------------ 5 files changed, 45 insertions(+), 63 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 9f08124811..8d6e622829 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -219,6 +219,7 @@ list(APPEND libopenmc_SOURCES src/secondary_correlated.cpp src/secondary_kalbach.cpp src/secondary_nbody.cpp + src/secondary_thermal.cpp src/secondary_uncorrelated.cpp src/settings.cpp src/simulation.cpp diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h index 9f8ebaf733..a084c18884 100644 --- a/include/openmc/thermal.h +++ b/include/openmc/thermal.h @@ -51,7 +51,19 @@ class ThermalData { public: ThermalData(hid_t group); - // Sample an outgoing energy and angle + //! Calculate the cross section + // + //! \param[in] E Incident neutron energy in [eV] + //! \param[out] elastic Elastic scattering cross section in [b] + //! \param[out] inelastic Inelastic scattering cross section in [b] + void calculate_xs(double E, double* elastic, double* inelastic) const; + + //! Sample an outgoing energy and angle + // + //! \param[in] micro_xs Microscopic cross sections + //! \param[in] E_in Incident neutron energy in [eV] + //! \param[out] E_out Outgoing neutron energy in [eV] + //! \param[out] mu Outgoing scattering angle cosine void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out, double* mu); private: @@ -106,10 +118,9 @@ public: void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out, double* mu); - double threshold() const { return data_[0].threshold_inelastic_; } - std::string name_; //!< name of table, e.g. "c_H_in_H2O" double awr_; //!< weight of nucleus in neutron masses + double energy_max_; //!< maximum energy for thermal scattering in [eV] std::vector kTs_; //!< temperatures in [eV] (k*T) std::vector nuclides_; //!< Valid nuclides diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index cdb5890930..cacc3478c5 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -372,6 +372,8 @@ class ThermalScattering(EqualityMixin): ---------- atomic_weight_ratio : float Atomic mass ratio of the target nuclide. + energy_max : float + Maximum energy for thermal scattering data in [eV] elastic : openmc.data.ThermalScatteringReaction or None Elastic scattering derived in the coherent or incoherent approximation inelastic : openmc.data.ThermalScatteringReaction @@ -391,9 +393,10 @@ class ThermalScattering(EqualityMixin): """ - def __init__(self, name, atomic_weight_ratio, kTs): + def __init__(self, name, atomic_weight_ratio, energy_max, kTs): self.name = name self.atomic_weight_ratio = atomic_weight_ratio + self.energy_max = energy_max self.kTs = kTs self.elastic = None self.inelastic = None @@ -432,6 +435,7 @@ class ThermalScattering(EqualityMixin): # 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): @@ -540,10 +544,11 @@ class ThermalScattering(EqualityMixin): name = group.name[1:] atomic_weight_ratio = group.attrs['atomic_weight_ratio'] + energy_max = group.attrs['energy_max'] kTg = group['kTs'] kTs = [dataset[()] for dataset in kTg.values()] - table = cls(name, atomic_weight_ratio, kTs) + table = cls(name, atomic_weight_ratio, energy_max, kTs) table.nuclides = [nuc.decode() for nuc in group.attrs['nuclides']] # Read thermal elastic scattering @@ -609,15 +614,13 @@ class ThermalScattering(EqualityMixin): # Assign temperature to the running list kTs = [ace.temperature*EV_PER_MEV] - table = cls(name, ace.atomic_weight_ratio, kTs) - T = table.temperatures[0] - # Incoherent inelastic scattering cross section idx = ace.jxs[1] n_energy = int(ace.xss[idx]) energy = ace.xss[idx+1 : idx+1+n_energy]*EV_PER_MEV xs = ace.xss[idx+1+n_energy : idx+1+2*n_energy] inelastic_xs = Tabulated1D(energy, xs) + energy_max = energy[-1] # Incoherent inelastic angle-energy distribution continuous = (ace.nxs[7] == 2) @@ -676,6 +679,8 @@ class ThermalScattering(EqualityMixin): distribution = IncoherentInelasticAEContinuous( breakpoints, interpolation, energy, energy_out, mu_out) + table = cls(name, ace.atomic_weight_ratio, energy_max, kTs) + T = table.temperatures[0] table.inelastic = ThermalScatteringReaction( {T: inelastic_xs}, {T: distribution} ) @@ -846,7 +851,7 @@ class ThermalScattering(EqualityMixin): data['free_atom_xs'] = B[0] data['epsilon'] = B[1] data['A0'] = awr = B[2] - data['e_max'] = B[3] + data['e_max'] = energy_max = B[3] data['M0'] = B[5] # Get information about non-principal atoms @@ -893,7 +898,7 @@ class ThermalScattering(EqualityMixin): data['effective_temperature'].append(Teff) name = ev.target['zsymam'].strip() - instance = cls(name, awr, kTs) + instance = cls(name, awr, energy_max, kTs) if elastic is not None: instance.elastic = elastic diff --git a/src/material.cpp b/src/material.cpp index 5667821fb7..ccbbe57a8d 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -774,7 +774,7 @@ void Material::calculate_neutron_xs(Particle& p) const // If particle energy is greater than the highest energy for the // S(a,b) table, then don't use the S(a,b) table - if (p.E_ > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE; + if (p.E_ > data::thermal_scatt[i_sab]->energy_max_) i_sab = C_NONE; // Increment position in thermal_tables_ ++j; diff --git a/src/thermal.cpp b/src/thermal.cpp index cd7a656a5d..5cc395e13c 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -44,6 +44,7 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector& tem name_ = name_.substr(1); read_attribute(group, "atomic_weight_ratio", awr_); + read_attribute(group, "energy_max", energy_max_); read_attribute(group, "nuclides", nuclides_); // Read temperatures @@ -177,58 +178,8 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp, // Set temperature index *i_temp = i; - // Get pointer to S(a,b) table - auto& sab = data_[i]; - - // Get index and interpolation factor for inelastic grid - int i_grid = 0; - double f = 0.0; - if (E >= sab.inelastic_e_in_.front()) { - auto& E_in = sab.inelastic_e_in_; - i_grid = lower_bound_index(E_in.begin(), E_in.end(), E); - f = (E - E_in[i_grid]) / (E_in[i_grid+1] - E_in[i_grid]); - } - - // Calculate S(a,b) inelastic scattering cross section - auto& xs = sab.inelastic_sigma_; - *inelastic = xs[i_grid] + f * (xs[i_grid + 1] - xs[i_grid]); - - // Check for elastic data - if (!sab.elastic_e_in_.empty()) { - // Determine whether elastic scattering is given in the coherent or - // incoherent approximation. For coherent, the cross section is - // represented as P/E whereas for incoherent, it is simply P - - auto& E_in = sab.elastic_e_in_; - - if (sab.elastic_mode_ == SAB_ELASTIC_COHERENT) { - if (E < E_in.front()) { - // If energy is below that of the lowest Bragg peak, the elastic - // cross section will be zero - *elastic = 0.0; - } else { - i_grid = lower_bound_index(E_in.begin(), E_in.end(), E); - *elastic = sab.elastic_P_[i_grid] / E; - } - } else { - // Determine index on elastic energy grid - if (E < E_in.front()) { - i_grid = 0; - } else { - i_grid = lower_bound_index(E_in.begin(), E_in.end(), E); - } - - // Get interpolation factor for elastic grid - f = (E - E_in[i_grid]) / (E_in[i_grid+1] - E_in[i_grid]); - - // Calculate S(a,b) elastic scattering cross section - auto& xs = sab.elastic_P_; - *elastic = xs[i_grid] + f*(xs[i_grid + 1] - xs[i_grid]); - } - } else { - // No elastic data - *elastic = 0.0; - } + // Calculate cross sections for ith temperature + data_[i].calculate_xs(E, elastic, inelastic); } bool @@ -308,6 +259,20 @@ ThermalData::ThermalData(hid_t group) } } +void +ThermalData::calculate_xs(double E, double* elastic, double* inelastic) const +{ + // Calculate thermal elastic scattering cross section + if (elastic_.xs) { + *elastic = (*elastic_.xs)(E); + } else { + *elastic = 0.0; + } + + // Calculate thermal inelastic scattering cross section + *inelastic = (*inelastic_.xs)(E); +} + void ThermalData::sample(const NuclideMicroXS& micro_xs, double E, double* E_out, double* mu) From c48d5f07c42c561795ef0e84e4899095e3fd59ec Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 22 May 2019 07:33:42 -0500 Subject: [PATCH 10/19] Change how incoherent elastic data is written to HDF5 --- openmc/data/angle_energy.py | 4 ++-- openmc/data/function.py | 4 ++-- openmc/data/thermal.py | 35 ++++++++++++++++++++--------- openmc/data/thermal_angle_energy.py | 6 ++--- src/endf.cpp | 7 ++++-- 5 files changed, 36 insertions(+), 20 deletions(-) diff --git a/openmc/data/angle_energy.py b/openmc/data/angle_energy.py index 5fb176f681..8cf9bcf5e8 100644 --- a/openmc/data/angle_energy.py +++ b/openmc/data/angle_energy.py @@ -43,8 +43,8 @@ class AngleEnergy(EqualityMixin, metaclass=ABCMeta): return openmc.data.IncoherentElasticAEDiscrete.from_hdf5(group) elif dist_type == 'incoherent_inelastic_discrete': return openmc.data.IncoherentInelasticAEDiscrete.from_hdf5(group) - elif dist_type == 'incoherent_inelastic_continuous': - return openmc.data.IncoherentInelasticAEContinuous.from_hdf5(group) + elif dist_type == 'incoherent_inelastic': + return openmc.data.IncoherentInelasticAE.from_hdf5(group) @staticmethod def from_ace(ace, location_dist, location_start, rx=None): diff --git a/openmc/data/function.py b/openmc/data/function.py index 99ee20900a..9175e5e3bf 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -591,8 +591,8 @@ class Regions1D(EqualityMixin): Functions which are to be combined in a piecewise fashion breakpoints : Iterable of float The values of the dependent variable that define the domain of - each function. The `i`th and `(i+1)`th values are the limits of the - domain of the `i`th function. Values must be monotonically increasing. + each function. The `i`\ th and `(i+1)`\ th values are the limits of the + domain of the `i`\ th function. Values must be monotonically increasing. Attributes ---------- diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index cacc3478c5..5ec9b83659 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -24,7 +24,7 @@ from .njoy import make_ace_thermal from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE, IncoherentElasticAEDiscrete, IncoherentInelasticAEDiscrete, - IncoherentInelasticAEContinuous) + IncoherentInelasticAE) _THERMAL_NAMES = { @@ -247,26 +247,26 @@ class IncoherentElastic(Function1D): Parameters ---------- - xs : float + bound_xs : float Characteristic bound cross section in [b] debye_waller : float Debye-Waller integral in [eV\ :math:`^{-1}`] Attributes ---------- - xs : float + bound_xs : float Characteristic bound cross section in [b] debye_waller : float Debye-Waller integral in [eV\ :math:`^{-1}`] """ - def __init__(self, xs, debye_waller): - self.xs = xs + def __init__(self, bound_xs, debye_waller): + self.bound_xs = bound_xs self.debye_waller = debye_waller def __call__(self, E): W = self.debye_waller - return self.xs / 2.0 * (1 - np.exp(-4*E*W)) / (2*E*W) + return self.bound_xs / 2.0 * (1 - np.exp(-4*E*W)) / (2*E*W) def to_hdf5(self, group, name): """Write incoherent elastic scattering to an HDF5 group @@ -279,14 +279,27 @@ class IncoherentElastic(Function1D): Name of the dataset to create """ - dataset = group.create_dataset(name) + data = np.array([self.bound_xs, self.debye_waller]) + dataset = group.create_dataset(name, data=data) dataset.attrs['type'] = np.string_(type(self).__name__) - dataset.attrs['debye_waller'] = self.debye_waller - dataset.attrs['bound_xs'] = self.xs @classmethod def from_hdf5(cls, dataset): - return cls(dataset.attrs['xs'], dataset.attrs['debye_waller']) + """Read incoherent elastic scattering from an HDF5 dataset + + Parameters + ---------- + dataset : h5py.Dataset + HDF5 dataset to read from + + Returns + ------- + openmc.data.IncoherentElastic + Incoherent elastic scattering cross section + + """ + bound_xs, debye_waller = dataset[()] + return cls(bound_xs, debye_waller) class ThermalScatteringReaction(EqualityMixin): @@ -676,7 +689,7 @@ class ThermalScattering(EqualityMixin): breakpoints = [n_energy] interpolation = [2] energy = inelastic_xs.x - distribution = IncoherentInelasticAEContinuous( + distribution = IncoherentInelasticAE( breakpoints, interpolation, energy, energy_out, mu_out) table = cls(name, ace.atomic_weight_ratio, energy_max, kTs) diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index d6053ed2a4..3537b91fa4 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -106,7 +106,7 @@ class IncoherentInelasticAEDiscrete(AngleEnergy): mu_out : numpy.ndarray Discrete angles for each incoming/outgoing energy skewed : bool - Whether distance angles are equi-probable or have a skewed distribution + Whether discrete angles are equi-probable or have a skewed distribution Attributes ---------- @@ -115,7 +115,7 @@ class IncoherentInelasticAEDiscrete(AngleEnergy): mu_out : numpy.ndarray Discrete angles for each incoming/outgoing energy skewed : bool - Whether distance angles are equi-probable or have a skewed distribution + Whether discrete angles are equi-probable or have a skewed distribution """ def __init__(self, energy_out, mu_out, skewed=False): @@ -137,5 +137,5 @@ class IncoherentInelasticAEDiscrete(AngleEnergy): return cls(energy_out, mu_out, skewed) -class IncoherentInelasticAEContinuous(CorrelatedAngleEnergy): +class IncoherentInelasticAE(CorrelatedAngleEnergy): _name = 'incoherent_inelastic' diff --git a/src/endf.cpp b/src/endf.cpp index 23512c4dbf..f01a643113 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -1,6 +1,7 @@ #include "openmc/endf.h" #include // for copy +#include #include // for log, exp #include // for back_inserter #include // for runtime_error @@ -244,8 +245,10 @@ double CoherentElasticXS::operator()(double E) const IncoherentElasticXS::IncoherentElasticXS(hid_t dset) { - read_attribute(dset, "bound_xs", bound_xs_); - read_attribute(dset, "debye_waller", debye_waller_); + std::array tmp; + read_dataset(dset, nullptr, tmp); + bound_xs_ = tmp[0]; + debye_waller_ = tmp[1]; } double IncoherentElasticXS::operator()(double E) const From 3cf937eaaaef8cb3fc5a921e27ef80ecfafd2a44 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 22 May 2019 07:33:56 -0500 Subject: [PATCH 11/19] Update I/O format for thermal neutron scattering --- docs/source/io_formats/nuclear_data.rst | 107 +++++++++++++++++++----- 1 file changed, 87 insertions(+), 20 deletions(-) diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index 957bfa3a3c..ef61604cb4 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -215,11 +215,9 @@ Thermal Neutron Scattering Data **//** :Attributes: - **atomic_weight_ratio** (*double*) -- Mass in units of neutron masses - - **nuclides** (*char[][]*) -- Names of nuclides for which the thermal - scattering data applies to - - **secondary_mode** (*char[]*) -- Indicates how the inelastic - outgoing angle-energy distributions are represented ('equal', - 'skewed', or 'continuous'). + - **energy_max** (*double*) -- Maximum energy in [eV] + - **nuclides** (*char[][]*) -- Names of nuclides for which the + thermal scattering data applies to **//kTs/** @@ -237,11 +235,13 @@ temperature-dependent data set. For example, the data set corresponding to temperature-dependent data set. For example, the data set corresponding to 300 Kelvin would be located at `300K`. -:Datasets: - **xs** (:ref:`tabulated <1d_tabulated>`) -- Thermal inelastic +:Datasets: + - **xs** (:ref:`function <1d_functions>`) -- Thermal elastic scattering cross section for temperature TTT (in Kelvin) - - **mu_out** (*double[][]*) -- Distribution of outgoing energies - and angles for coherent elastic scattering for temperature TTT - (in Kelvin) + +:Groups: + - **distribution** -- Format for angle-energy distributions are + detailed in :ref:`angle_energy`. **//inelastic/K/** @@ -249,18 +249,13 @@ temperature-dependent data set. For example, the data set corresponding to temperature-dependent data set. For example, the data set corresponding to 300 Kelvin would be located at `300K`. -:Datasets: - **xs** (:ref:`tabulated <1d_tabulated>`) -- Thermal inelastic +:Datasets: + - **xs** (:ref:`function <1d_functions>`) -- Thermal inelastic scattering cross section for temperature TTT (in Kelvin) - - **energy_out** (*double[][]*) -- Distribution of outgoing - energies for each incoming energy for temperature TTT (in Kelvin). - Only present if secondary mode is not continuous. - - **mu_out** (*double[][][]*) -- Distribution of scattering cosines - for each pair of incoming and outgoing energies. for temperature - TTT (in Kelvin). Only present if secondary mode is not continuous. -If the secondary mode is continuous, the outgoing energy-angle distribution is -given as a :ref:`correlated angle-energy distribution -`. +:Groups: + - **distribution** -- Format for angle-energy distributions are + detailed in :ref:`angle_energy`. .. _product: @@ -332,7 +327,17 @@ Coherent elastic scattering :Datatype: *double[2][]* :Description: The first row lists Bragg edges and the second row lists structure factor cumulative sums. -:Attributes: - **type** (*char[]*) -- 'bragg' +:Attributes: - **type** (*char[]*) -- 'CoherentElastic' + +Incoherent elastic scattering +----------------------------- + +:Object type: Dataset +:Datatype: *double[2]* +:Description: The first value is the characteristic bound cross section in [b] + and the second value is the Debye-Waller integral in + [eV\ :math:`^{-1}`]. +:Attributes: - **type** (*char[]*) -- 'IncoherentElastic' .. _angle_energy: @@ -434,6 +439,68 @@ N-Body Phase Space target nuclide in neutron masses - **q_value** (*double*) -- Q value for the reaction in eV +Coherent Elastic +---------------- + +This angle-energy distribution is used specifically for coherent elastic thermal +neutron scattering. + +:Object type: Group +:Attributes: - **type** (*char[]*) -- "coherent_elastic" +:Hard link: - **xs** -- Link to the coherent elastic scattering cross section + +Incoherent Elastic +------------------ + +This angle-energy distribution is used specifically for incoherent elastic +thermal neutron scattering (derived from an ENDF file directly). + +:Object type: Group +:Attributes: - **type** (*char[]*) -- "incoherent_elastic" +:Datasets: + - **debye_waller** (*double*) -- Debye-Waller integral in + [eV\ :math:`^{-1}`] + +Incoherent Elastic (Discrete) +----------------------------- + +This angle-energy distribution is used for discretized incoherent elastic +thermal neutron scattering distributions that are present in ACE files. + +:Object type: Group +:Attributes: - **type** (*char[]*) -- "incoherent_elastic_discrete" +:Datasets: + - **mu_out** (*double[][]*) -- Equiprobable discrete outgoing + angles for each incident neutron energy tabulated + +Incoherent Inelastic +-------------------- + +This angle-energy distribution is used specifically for (continuous) incoherent +inelastic thermal neutron scattering. + +:Object type: Group +:Attributes: - **type** (*char[]*) -- "incoherent_inelastic" +:Datasets: The datasets for this angle-energy distribution are the same as for + :ref:`correlated angle-energy distributions + `. + +Incoherent Inelastic (Discrete) +------------------------------- + +This angle-energy distribution is used specifically for incoherent inelastic +thermal neutron scattering where the distributions have been discretized into +equiprobable bins. + +:Object type: Group +:Attributes: - **type** (*char[]*) -- "incoherent_inelastic_discrete" +:Datasets: - **energy_out** (*double[][]*) -- Distribution of outgoing + energies for each incoming energy. + - **mu_out** (*double[][][]*) -- Distribution of scattering cosines + for each pair of incoming and outgoing energies. + - **skewed** (*int8_t*) -- Whether discrete angles are equi-probable + (0) or have a skewed distribution (1). + .. _energy_distribution: -------------------- From 228392eee592d946745a5f212dc538885111fe75 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 22 May 2019 10:15:13 -0500 Subject: [PATCH 12/19] Fix existing unit tests for thermal scattering --- tests/unit_tests/test_data_thermal.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/tests/unit_tests/test_data_thermal.py b/tests/unit_tests/test_data_thermal.py index 7d7d4e964a..d8b21c99c6 100644 --- a/tests/unit_tests/test_data_thermal.py +++ b/tests/unit_tests/test_data_thermal.py @@ -38,14 +38,14 @@ def h2o_njoy(): def test_h2o_attributes(h2o): assert h2o.name == 'c_H_in_H2O' assert h2o.nuclides == ['H1'] - assert h2o.secondary_mode == 'skewed' assert h2o.temperatures == ['294K'] assert h2o.atomic_weight_ratio == pytest.approx(0.999167) + assert h2o.energy_max == pytest.approx(4.46) def test_h2o_xs(h2o): - assert not h2o.elastic_xs - for temperature, func in h2o.inelastic_xs.items(): + assert not h2o.elastic + for temperature, func in h2o.inelastic.xs.items(): assert temperature.endswith('K') assert isinstance(func, Callable) @@ -53,20 +53,20 @@ def test_h2o_xs(h2o): def test_graphite_attributes(graphite): assert graphite.name == 'c_Graphite' assert graphite.nuclides == ['C0', 'C12', 'C13'] - assert graphite.secondary_mode == 'skewed' assert graphite.temperatures == ['296K'] assert graphite.atomic_weight_ratio == pytest.approx(11.898) + assert graphite.energy_max == pytest.approx(4.46) def test_graphite_xs(graphite): - for temperature, func in graphite.elastic_xs.items(): + for temperature, func in graphite.elastic.xs.items(): assert temperature.endswith('K') assert isinstance(func, openmc.data.CoherentElastic) - for temperature, func in graphite.inelastic_xs.items(): + for temperature, func in graphite.inelastic.xs.items(): assert temperature.endswith('K') assert isinstance(func, Callable) - elastic = graphite.elastic_xs['296K'] - assert elastic([1e-3, 1.0]) == pytest.approx([13.47464936, 0.62590156]) + elastic = graphite.elastic.xs['296K'] + assert elastic([1e-3, 1.0]) == pytest.approx([0.0, 0.62586153]) def test_export_to_hdf5(tmpdir, h2o_njoy, graphite): @@ -80,9 +80,9 @@ def test_export_to_hdf5(tmpdir, h2o_njoy, graphite): def test_continuous_dist(h2o_njoy): - for temperature, dist in h2o_njoy.inelastic_dist.items(): + for temperature, dist in h2o_njoy.inelastic.distribution.items(): assert temperature.endswith('K') - assert isinstance(dist, openmc.data.CorrelatedAngleEnergy) + assert isinstance(dist, openmc.data.IncoherentInelasticAE) def test_get_thermal_name(): From 707eefbf65ed06bcf33569105519d7da77fcd4ec Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 22 May 2019 12:58:31 -0500 Subject: [PATCH 13/19] Allow iwt to be specified on ACER runs for thermal scattering --- openmc/data/njoy.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/openmc/data/njoy.py b/openmc/data/njoy.py index ddc1efb1c4..c8d630851b 100644 --- a/openmc/data/njoy.py +++ b/openmc/data/njoy.py @@ -120,7 +120,7 @@ acer / %%%%%%%%%%%%%%%%%%%%%%%% Write out in ACE format %%%%%%%%%%%%%%%%%%%%%%%% '{library}: {zsymam_thermal} processed by NJOY'/ {mat} {temperature} '{data.name}' / {zaids} / -222 64 {mt_elastic} {elastic_type} {data.nmix} {energy_max} 2/ +222 64 {mt_elastic} {elastic_type} {data.nmix} {energy_max} {iwt}/ """ @@ -355,8 +355,8 @@ def make_ace(filename, temperatures=None, ace='ace', xsdir='xsdir', pendf=None, def make_ace_thermal(filename, filename_thermal, temperatures=None, - ace='ace', xsdir='xsdir', error=0.001, evaluation=None, - evaluation_thermal=None, **kwargs): + ace='ace', xsdir='xsdir', error=0.001, iwt=2, + evaluation=None, evaluation_thermal=None, **kwargs): """Generate thermal scattering ACE file from ENDF files Parameters @@ -374,6 +374,8 @@ def make_ace_thermal(filename, filename_thermal, temperatures=None, Path of xsdir file to write error : float, optional Fractional error tolerance for NJOY processing + iwt : int + `iwt` parameter used in NJOR/ACER card 9 evaluation : openmc.data.endf.Evaluation, optional If the ENDF neutron sublibrary file contains multiple material evaluations, this argument indicates which evaluation to use. From c8bf27af4236636d72ca16e839f139d61c053f89 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 22 May 2019 14:15:26 -0500 Subject: [PATCH 14/19] 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' From adce57093d414e10810640855d8e7d53635862be Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 May 2019 07:20:52 -0500 Subject: [PATCH 15/19] Add ability to use ENDF incoherent elastic data, update unit tests --- openmc/data/thermal.py | 22 +++++++++++- tests/unit_tests/test_data_thermal.py | 51 ++++++++++++++++++--------- 2 files changed, 56 insertions(+), 17 deletions(-) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 6f39e478c8..3a4713097e 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -726,7 +726,8 @@ class ThermalScattering(EqualityMixin): @classmethod def from_njoy(cls, filename, filename_thermal, temperatures=None, - evaluation=None, evaluation_thermal=None, **kwargs): + evaluation=None, evaluation_thermal=None, + use_endf_data=True, **kwargs): """Generate thermal scattering data by running NJOY. Parameters @@ -746,6 +747,9 @@ class ThermalScattering(EqualityMixin): If the ENDF thermal scattering sublibrary file contains multiple material evaluations, this argument indicates which evaluation to use. + use_endf_data : bool + If the material has incoherent elastic scattering, the ENDF data + will be used rather than the ACE data. **kwargs Keyword arguments passed to :func:`openmc.data.njoy.make_ace_thermal` @@ -770,6 +774,22 @@ class ThermalScattering(EqualityMixin): for table in lib.tables[1:]: data.add_temperature_from_ace(table) + # Load ENDF data to replace incoherent elastic + if use_endf_data: + data_endf = cls.from_endf(filename_thermal) + if data_endf.elastic is not None: + # Get appropriate temperatures + if temperatures is None: + temperatures = data_endf.temperatures + else: + temperatures = [_temperature_str(t) for t in temperatures] + + # Replace ACE data with ENDF data + rx, rx_endf = data.elastic, data_endf.elastic + for t in temperatures: + rx.xs[t] = rx_endf.xs[t] + rx.distribution[t] = rx_endf.distribution[t] + return data @classmethod diff --git a/tests/unit_tests/test_data_thermal.py b/tests/unit_tests/test_data_thermal.py index 63bb98f6a7..47cbe27b97 100644 --- a/tests/unit_tests/test_data_thermal.py +++ b/tests/unit_tests/test_data_thermal.py @@ -47,11 +47,16 @@ def hzrh(): @pytest.fixture(scope='module') def hzrh_njoy(): - """H in ZrH genertaed using NJOY.""" + """H in ZrH generated 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) + with_endf_data = openmc.data.ThermalScattering.from_njoy( + path_h1, path_hzrh, temperatures=[296.0], iwt=0 + ) + without_endf_data = openmc.data.ThermalScattering.from_njoy( + path_h1, path_hzrh, temperatures=[296.0], use_endf_data=False, iwt=1 + ) + return with_endf_data, without_endf_data @pytest.fixture(scope='module') @@ -117,10 +122,12 @@ def test_export_to_hdf5(tmpdir, h2o_njoy, hzrh_njoy, graphite): 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 + # H in ZrH covers export of incoherent elastic data, and incoherent + # inelastic angle-energy distributions filename = str(tmpdir.join('hzrh.h5')) - hzrh_njoy.export_to_hdf5(filename) + hzrh_njoy[0].export_to_hdf5(filename) + assert os.path.exists(filename) + hzrh_njoy[1].export_to_hdf5(filename, 'w') assert os.path.exists(filename) @@ -167,20 +174,32 @@ def test_hzrh_elastic(hzrh): 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'] + endf, ace = hzrh_njoy - # Check incoherent elastic distribution - d = hzrh.elastic.distribution['296K'] + # First check version using ENDF incoherent elastic data + assert endf.atomic_weight_ratio == pytest.approx(0.999167) + assert endf.energy_max == pytest.approx(1.855) + assert endf.temperatures == ['296K'] + + # Now check version using ACE incoherent elastic data (discretized) + assert ace.atomic_weight_ratio == endf.atomic_weight_ratio + assert ace.energy_max == endf.energy_max + + # Cross sections should be about the same (within 1%) + E = np.linspace(1e-5, endf.energy_max) + xs1 = endf.elastic.xs['296K'](E) + xs2 = ace.elastic.xs['296K'](E) + assert xs1 == pytest.approx(xs2, rel=0.01) + + # Check discrete incoherent elastic distribution + d = ace.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'] + # Check discrete incoherent inelastic distribution + d = endf.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)) + 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*endf.energy_max)) def test_sio2_attributes(sio2): From 7a1afc5869f94272b77991093e2dd2cdc157d710 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 May 2019 09:26:46 -0500 Subject: [PATCH 16/19] Store 0-based threshold indices in HDF5 files --- openmc/data/reaction.py | 7 ++----- src/reaction.cpp | 2 -- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index c5abc73c34..05b6b8d2d6 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -894,10 +894,7 @@ class Reaction(EqualityMixin): Tgroup = group.create_group(T) if self.xs[T] is not None: dset = Tgroup.create_dataset('xs', data=self.xs[T].y) - if hasattr(self.xs[T], '_threshold_idx'): - threshold_idx = self.xs[T]._threshold_idx + 1 - else: - threshold_idx = 1 + threshold_idx = getattr(self.xs[T], '_threshold_idx', 0) dset.attrs['threshold_idx'] = threshold_idx for i, p in enumerate(self.products): pgroup = group.create_group('product_{}'.format(i)) @@ -939,7 +936,7 @@ class Reaction(EqualityMixin): 'at T={} because no corresponding energy grid ' 'exists.'.format(mt, T)) xs = Tgroup['xs'][()] - threshold_idx = Tgroup['xs'].attrs['threshold_idx'] - 1 + threshold_idx = Tgroup['xs'].attrs['threshold_idx'] tabulated_xs = Tabulated1D(energy[T][threshold_idx:], xs) tabulated_xs._threshold_idx = threshold_idx rx.xs[T] = tabulated_xs diff --git a/src/reaction.cpp b/src/reaction.cpp index 222822230d..b733b20283 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -42,8 +42,6 @@ Reaction::Reaction(hid_t group, const std::vector& temperatures) // Get threshold index TemperatureXS xs; read_attribute(dset, "threshold_idx", xs.threshold); - // TODO: change HDF5 format so that threshold_idx is 0-based - --xs.threshold; // Read cross section values read_dataset(dset, xs.value); From c2e612a5031b35d47b5e46d8536c77b35d5d3b51 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 May 2019 09:22:44 -0500 Subject: [PATCH 17/19] Update HDF5 data version to 3.0 --- docs/source/usersguide/cross_sections.rst | 8 ++++---- include/openmc/constants.h | 3 +-- openmc/data/__init__.py | 2 +- tools/ci/download-xs.sh | 2 +- 4 files changed, 7 insertions(+), 8 deletions(-) diff --git a/docs/source/usersguide/cross_sections.rst b/docs/source/usersguide/cross_sections.rst index e6a79815f9..636b44db46 100644 --- a/docs/source/usersguide/cross_sections.rst +++ b/docs/source/usersguide/cross_sections.rst @@ -11,7 +11,7 @@ or multi-group mode. In continuous-energy mode, OpenMC uses a native `HDF5 `_ format (see :ref:`io_nuclear_data`) to store all nuclear data. Pregenerated HDF5 libraries can be found at -https://openmc.mcs.anl.gov; unless you have specific data needs, it is highly +https://openmc.org; unless you have specific data needs, it is highly recommended to use one of the pregenerated libraries. Alternatively, if you have ACE format data that was produced with NJOY_, such as that distributed with MCNP_ or Serpent_, it can be converted to the HDF5 format using the :ref:`using @@ -61,7 +61,7 @@ Using Pregenerated Libraries ---------------------------- Various evaluated nuclear data libraries have been processed into the HDF5 -format required by OpenMC and can be found at https://openmc.mcs.anl.gov. You +format required by OpenMC and can be found at https://openmc.org. You can find both libraries generated by the OpenMC development team as well as libraries based on ACE files distributed elsewhere. To use these libraries, download the archive file, unpack it, and then set your @@ -198,7 +198,7 @@ file, is distributed with OpenMC. The rest is available from the NNDC_, which provides ENDF data from the photo-atomic and atomic relaxation sublibraries of the ENDF/B-VII.1 library. -Most of the pregenerated HDF5 libraries available at https://openmc.mcs.anl.gov +Most of the pregenerated HDF5 libraries available at https://openmc.org already have photon interaction data included. If you are building a data library yourself, it is possible to use the Python API directly to convert photon interaction data from an ENDF or ACE file to an HDF5 file. The @@ -237,7 +237,7 @@ the :class:`openmc.data.DataLibrary` class to register the .h5 files as described in :ref:`create_xs_library`. The `official ENDF/B-VII.1 HDF5 library -`_ includes the windowed +`_ includes the windowed multipole library, so if you are using this library, the windowed multipole data will already be available to you. diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 6b2e6f63b7..394f90740f 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -26,7 +26,7 @@ constexpr bool VERSION_DEV {true}; constexpr std::array VERSION {VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE}; // HDF5 data format -constexpr int HDF5_VERSION[] {2, 0}; +constexpr int HDF5_VERSION[] {3, 0}; // Version numbers for binary files constexpr std::array VERSION_STATEPOINT {17, 0}; @@ -36,7 +36,6 @@ constexpr std::array VERSION_SUMMARY {6, 0}; constexpr std::array VERSION_VOLUME {1, 0}; constexpr std::array VERSION_VOXEL {2, 0}; constexpr std::array VERSION_MGXS_LIBRARY {1, 0}; -constexpr char VERSION_MULTIPOLE[] {"v0.2"}; // ============================================================================ // ADJUSTABLE PARAMETERS diff --git a/openmc/data/__init__.py b/openmc/data/__init__.py index 92acc96354..277c14ca8d 100644 --- a/openmc/data/__init__.py +++ b/openmc/data/__init__.py @@ -1,5 +1,5 @@ # Version of HDF5 nuclear data format -HDF5_VERSION_MAJOR = 2 +HDF5_VERSION_MAJOR = 3 HDF5_VERSION_MINOR = 0 HDF5_VERSION = (HDF5_VERSION_MAJOR, HDF5_VERSION_MINOR) diff --git a/tools/ci/download-xs.sh b/tools/ci/download-xs.sh index 38a5bc4692..07ade9c70f 100755 --- a/tools/ci/download-xs.sh +++ b/tools/ci/download-xs.sh @@ -3,7 +3,7 @@ set -ex # Download HDF5 data if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then - wget -q -O - https://anl.box.com/shared/static/u1g3n8iai0u1n5f6ev3pg2j3ff941bqa.xz | tar -C $HOME -xJ + wget -q -O - https://anl.box.com/shared/static/teaup95cqv8s9nn56hfn7ku8mmelr95p.xz | tar -C $HOME -xJ fi # Download ENDF/B-VII.1 distribution From 97ea3287c2ff2fddce065717ee23616f826cf037 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 25 Jun 2019 09:47:26 -0500 Subject: [PATCH 18/19] Only replace incoherent elastic in ThermalScattering.from_njoy with ENDF data --- openmc/data/thermal.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 3a4713097e..79cad32823 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -787,8 +787,9 @@ class ThermalScattering(EqualityMixin): # Replace ACE data with ENDF data rx, rx_endf = data.elastic, data_endf.elastic for t in temperatures: - rx.xs[t] = rx_endf.xs[t] - rx.distribution[t] = rx_endf.distribution[t] + if isinstance(rx_endf.xs[t], IncoherentElastic): + rx.xs[t] = rx_endf.xs[t] + rx.distribution[t] = rx_endf.distribution[t] return data From d9b27cca5384bade2308be77d0ad7929127c4d65 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 1 Jul 2019 07:09:10 -0500 Subject: [PATCH 19/19] Respond to @nelsonag comments on #1271 --- docs/source/pythonapi/data.rst | 3 +++ include/openmc/secondary_thermal.h | 28 +++++++++++++++++++++++++++- src/secondary_thermal.cpp | 8 ++++---- 3 files changed, 34 insertions(+), 5 deletions(-) diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 144d263066..92cb785481 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -5,6 +5,9 @@ Core Classes ------------ +The following classes are used for incident neutron data, decay data, fission +and product yields. + .. autosummary:: :toctree: generated :nosignatures: diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 998143d0f5..7b876334f8 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -42,6 +42,9 @@ private: class IncoherentElasticAE : public AngleEnergy { public: + //! Construct from HDF5 file + // + //! \param[in] group HDF5 group explicit IncoherentElasticAE(hid_t group); //! Sample distribution for an angle and energy @@ -59,11 +62,19 @@ private: class IncoherentElasticAEDiscrete : public AngleEnergy { public: + //! Construct from HDF5 file + // + //! \param[in] group HDF5 group + //! \param[in] energy Energies at which cosines are tabulated explicit IncoherentElasticAEDiscrete(hid_t group, const std::vector& energy); + //! Sample distribution for an angle and energy + //! \param[in] E_in Incoming energy in [eV] + //! \param[out] E_out Outgoing energy in [eV] + //! \param[out] mu Outgoing cosine with respect to current direction void sample(double E_in, double& E_out, double& mu) const override; private: - const std::vector& energy_; //!< Incoherent inelastic scattering cross section + const std::vector& energy_; //!< Energies at which cosines are tabulated xt::xtensor mu_out_; //!< Cosines for each incident energy }; @@ -73,8 +84,16 @@ private: class IncoherentInelasticAEDiscrete : public AngleEnergy { public: + //! Construct from HDF5 file + // + //! \param[in] group HDF5 group + //! \param[in] energy Incident energies at which distributions are tabulated explicit IncoherentInelasticAEDiscrete(hid_t group, const std::vector& energy); + //! Sample distribution for an angle and energy + //! \param[in] E_in Incoming energy in [eV] + //! \param[out] E_out Outgoing energy in [eV] + //! \param[out] mu Outgoing cosine with respect to current direction void sample(double E_in, double& E_out, double& mu) const override; private: const std::vector& energy_; //!< Incident energies @@ -89,8 +108,15 @@ private: class IncoherentInelasticAE : public AngleEnergy { public: + //! Construct from HDF5 file + // + //! \param[in] group HDF5 group explicit IncoherentInelasticAE(hid_t group); + //! Sample distribution for an angle and energy + //! \param[in] E_in Incoming energy in [eV] + //! \param[out] E_out Outgoing energy in [eV] + //! \param[out] mu Outgoing cosine with respect to current direction void sample(double E_in, double& E_out, double& mu) const override; private: //! Secondary energy/angle distribution diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index 5fdf6201a0..696907a17e 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -240,13 +240,13 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu) const // Determine endpoints on grid i auto n = distribution_[i].e_out.size(); - double E_i_1 = distribution_[i].e_out(0); - double E_i_J = distribution_[i].e_out(n - 1); + double E_i_1 = distribution_[i].e_out[0]; + double E_i_J = distribution_[i].e_out[n - 1]; // Determine endpoints on grid i + 1 n = distribution_[i + 1].e_out.size(); - double E_i1_1 = distribution_[i + 1].e_out(0); - double E_i1_J = distribution_[i + 1].e_out(n - 1); + double E_i1_1 = distribution_[i + 1].e_out[0]; + double E_i1_J = distribution_[i + 1].e_out[n - 1]; double E_1 = E_i_1 + f * (E_i1_1 - E_i_1); double E_J = E_i_J + f * (E_i1_J - E_i_J);