From d024a3e00842a6535d8d0711a4131e4bee1ee841 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Oct 2021 07:26:43 -0500 Subject: [PATCH 01/10] Implementation of C++ Sum1D function --- include/openmc/endf.h | 20 ++++++++++++++++++++ src/endf.cpp | 26 ++++++++++++++++++++++++++ 2 files changed, 46 insertions(+) diff --git a/include/openmc/endf.h b/include/openmc/endf.h index 7beb8e452..e580874b4 100644 --- a/include/openmc/endf.h +++ b/include/openmc/endf.h @@ -126,6 +126,26 @@ private: debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1] }; +//============================================================================== +//! Sum of multiple 1D functions +//============================================================================== + +class Sum1D : public Function1D { +public: + // Constructors + explicit Sum1D(hid_t group); + + //! Evaluate each function and sum results + //! \param[in] x independent variable + //! \return Function evaluated at x + double operator()(double E) const override; + + const unique_ptr& functions(int i) const { return functions_[i]; } + +private: + vector> functions_; //!< individual functions +}; + //! Read 1D function from HDF5 dataset //! \param[in] group HDF5 group containing dataset //! \param[in] name Name of dataset diff --git a/src/endf.cpp b/src/endf.cpp index b42c8641d..60db3efa4 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -271,4 +271,30 @@ double IncoherentElasticXS::operator()(double E) const return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W)); } +//============================================================================== +// Sum1D implementation +//============================================================================== + +Sum1D::Sum1D(hid_t group) +{ + // Get number of functions + int n; + read_attribute(group, "n", n); + + // Get each function + for (int i = 0; i < n; ++i) { + auto dset_name = fmt::format("func_{}", i + 1); + functions_.push_back(read_function(group, dset_name.c_str())); + } +} + +double Sum1D::operator()(double x) const +{ + double result = 0.0; + for (auto& func : functions_) { + result += (*func)(x); + } + return result; +} + } // namespace openmc From 8bb2002d8f8bebd58f11d84d88fc88df20552e8d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Oct 2021 07:27:56 -0500 Subject: [PATCH 02/10] Add open_object, close_object functions in HDF5 interface --- include/openmc/hdf5_interface.h | 2 ++ src/endf.cpp | 16 ++++++++-------- src/hdf5_interface.cpp | 12 ++++++++++++ 3 files changed, 22 insertions(+), 8 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index e9d9b4f96..0092c08f8 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -61,6 +61,8 @@ void ensure_exists(hid_t obj_id, const char* name, bool attribute = false); vector group_names(hid_t group_id); vector object_shape(hid_t obj_id); std::string object_name(hid_t obj_id); +hid_t open_object(hid_t group_id, const std::string& name); +void close_object(hid_t obj_id); //============================================================================== // Fortran compatibility functions diff --git a/src/endf.cpp b/src/endf.cpp index 60db3efa4..3ccfdbd21 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -90,23 +90,23 @@ bool is_inelastic_scatter(int mt) unique_ptr read_function(hid_t group, const char* name) { - hid_t dset = open_dataset(group, name); + hid_t obj_id = open_object(group, name); std::string func_type; - read_attribute(dset, "type", func_type); + read_attribute(obj_id, "type", func_type); unique_ptr func; if (func_type == "Tabulated1D") { - func = make_unique(dset); + func = make_unique(obj_id); } else if (func_type == "Polynomial") { - func = make_unique(dset); + func = make_unique(obj_id); } else if (func_type == "CoherentElastic") { - func = make_unique(dset); + func = make_unique(obj_id); } else if (func_type == "IncoherentElastic") { - func = make_unique(dset); + func = make_unique(obj_id); } else { throw std::runtime_error {"Unknown function type " + func_type + - " for dataset " + object_name(dset)}; + " for dataset " + object_name(obj_id)}; } - close_dataset(dset); + close_object(obj_id); return func; } diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index d53ad3879..27b750b93 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -118,6 +118,12 @@ void close_group(hid_t group_id) fatal_error("Failed to close group"); } +void close_object(hid_t obj_id) +{ + if (H5Oclose(obj_id) < 0) + fatal_error("Failed to close object"); +} + int dataset_ndims(hid_t dset) { hid_t dspace = H5Dget_space(dset); @@ -394,6 +400,12 @@ hid_t open_group(hid_t group_id, const char* name) return H5Gopen(group_id, name, H5P_DEFAULT); } +hid_t open_object(hid_t group_id, const std::string& name) +{ + ensure_exists(group_id, name.c_str()); + return H5Oopen(group_id, name.c_str(), H5P_DEFAULT); +} + void read_attr(hid_t obj_id, const char* name, hid_t mem_type_id, void* buffer) { hid_t attr = H5Aopen(obj_id, name, H5P_DEFAULT); From e453901c84e57eca983a52124e7dc903b06dfb00 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Oct 2021 07:28:35 -0500 Subject: [PATCH 03/10] Implement to/from_hdf5 methods for Python Sum class --- openmc/data/function.py | 39 +++++++++++++++++++++++++++++++++++++++ src/endf.cpp | 2 ++ 2 files changed, 41 insertions(+) diff --git a/openmc/data/function.py b/openmc/data/function.py index b5aa2117d..7a73987df 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -578,6 +578,45 @@ class Sum(EqualityMixin): cv.check_type('functions', functions, Iterable, Callable) self._functions = functions + def to_hdf5(self, group, name='xy'): + """Write sum of functions to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + name : str + Name of the dataset to create + + """ + sum_group = group.create_group(name) + sum_group.attrs['type'] = np.string_(type(self).__name__) + sum_group.attrs['n'] = len(self.functions) + for i, f in enumerate(self.functions): + f.to_hdf5(sum_group, f'func_{i+1}') + + @classmethod + def from_hdf5(cls, group): + """Generate sum of functions from an HDF5 group + + Parameters + ---------- + group : h5py.Group + Group to read from + + Returns + ------- + openmc.data.Sum + Functions read from the group + + """ + n = group.attrs['n'] + functions = [ + Function1D.from_hdf5(group[f'func_{i+1}']) + for i in range(n) + ] + return cls(functions) + class Regions1D(EqualityMixin): r"""Piecewise composition of multiple functions. diff --git a/src/endf.cpp b/src/endf.cpp index 3ccfdbd21..c0c1d2e7e 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -102,6 +102,8 @@ unique_ptr read_function(hid_t group, const char* name) func = make_unique(obj_id); } else if (func_type == "IncoherentElastic") { func = make_unique(obj_id); + } else if (func_type == "Sum") { + func = make_unique(obj_id); } else { throw std::runtime_error {"Unknown function type " + func_type + " for dataset " + object_name(obj_id)}; From 70c16d4b04a60606742f6b1c6825bfd04adcad6d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Oct 2021 07:30:56 -0500 Subject: [PATCH 04/10] Implement mixed coherent/incoheret thermal elastic --- include/openmc/secondary_thermal.h | 28 +++++++++ openmc/data/angle_energy.py | 2 + openmc/data/thermal.py | 91 ++++++++++++++++++++--------- openmc/data/thermal_angle_energy.py | 57 ++++++++++++++++++ src/secondary_thermal.cpp | 35 +++++++++++ src/thermal.cpp | 25 +++++--- 6 files changed, 203 insertions(+), 35 deletions(-) diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 84ed68f51..81de5c451 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -150,6 +150,34 @@ private: //!< each incident energy }; +//============================================================================== +//! Mixed coherent/incoherent elastic angle-energy distribution +//============================================================================== + +class MixedElasticAE : public AngleEnergy { +public: + //! Construct from HDF5 file + // + //! \param[in] group HDF5 group + explicit MixedElasticAE( + hid_t group, const CoherentElasticXS& coh_xs, const Tabulated1D& incoh_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 + //! \param[inout] seed Pseudorandom number seed pointer + void sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const override; + +private: + CoherentElasticAE coherent_dist_; //!< Coherent distribution + unique_ptr incoherent_dist_; //!< Incoherent distribution + + const CoherentElasticXS& coherent_xs_; //!< Ref. to coherent XS + const Tabulated1D& incoherent_xs_; //!< Ref. to incoherent XS +}; + } // namespace openmc #endif // OPENMC_SECONDARY_THERMAL_H diff --git a/openmc/data/angle_energy.py b/openmc/data/angle_energy.py index 6009dc748..b8fde5478 100644 --- a/openmc/data/angle_energy.py +++ b/openmc/data/angle_energy.py @@ -44,6 +44,8 @@ class AngleEnergy(EqualityMixin, ABC): return openmc.data.IncoherentInelasticAEDiscrete.from_hdf5(group) elif dist_type == 'incoherent_inelastic': return openmc.data.IncoherentInelasticAE.from_hdf5(group) + elif dist_type == 'mixed_elastic': + return openmc.data.MixedElasticAE.from_hdf5(group) @staticmethod def from_ace(ace, location_dist, location_start, rx=None): diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 4e0a30b64..6cf5a4e9e 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, isotopes from .ace import Table, get_table, Library from .angle_energy import AngleEnergy -from .function import Tabulated1D, Function1D +from .function import Tabulated1D, Function1D, Sum from .njoy import make_ace_thermal from .thermal_angle_energy import (CoherentElasticAE, IncoherentElasticAE, IncoherentElasticAEDiscrete, IncoherentInelasticAEDiscrete, - IncoherentInelasticAE) + IncoherentInelasticAE, MixedElasticAE) _THERMAL_NAMES = { @@ -694,29 +694,53 @@ class ThermalScattering(EqualityMixin): # Incoherent/coherent elastic scattering cross section idx = ace.jxs[4] - n_mu = ace.nxs[6] + 1 if idx != 0: - n_energy = int(ace.xss[idx]) - energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV - P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] - - if ace.nxs[5] == 4: + if ace.nxs[5] in (4, 5): # Coherent elastic - xs = CoherentElastic(energy, P*EV_PER_MEV) - distribution = CoherentElasticAE(xs) + n_energy = int(ace.xss[idx]) + energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV + P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] + coherent_xs = CoherentElastic(energy, P*EV_PER_MEV) + coherent_dist = CoherentElasticAE(xs) # Coherent elastic shouldn't have angular distributions listed + n_mu = ace.nxs[6] + 1 assert n_mu == 0 - else: - # Incoherent elastic - xs = Tabulated1D(energy, P) + + if ace.nxs[5] in (3, 5): + # Incoherent elastic scattering -- first determine if both + # incoherent and coherent are present (mixed) + mixed = (ace.nxs[5] == 5) + + # Get cross section values + idx = ace.jxs[7] if mixed else ace.jxs[4] + n_energy = int(ace.xss[idx]) + energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV + values = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] + + incoherent_xs = Tabulated1D(energy, values) # Angular distribution + n_mu = (ace.nxs[8] if mixed else ace.nxs[6]) + 1 assert n_mu > 0 - idx = ace.jxs[6] + idx = ace.jxs[9] if mixed else ace.jxs[6] mu_out = ace.xss[idx:idx + n_energy * n_mu] mu_out.shape = (n_energy, n_mu) - distribution = IncoherentElasticAEDiscrete(mu_out) + incoherent_dist = IncoherentElasticAEDiscrete(mu_out) + + if ace.nxs[5] == 3: + xs = incoherent_xs + dist = incoherent_dist + elif ace.nxs[5] == 4: + xs = coherent_xs + dist = coherent_dist + else: + # Create mixed cross section -- note that coherent must come + # first due to assumption on C++ side + xs = Sum([coherent_xs, incoherent_xs]) + + # Create mixed distribution + distribution = MixedElasticAE(coherent_dist, incoherent_dist) table.elastic = ThermalScatteringReaction({T: xs}, {T: distribution}) @@ -802,7 +826,7 @@ class ThermalScattering(EqualityMixin): # Replace ACE data with ENDF data rx, rx_endf = data.elastic, data_endf.elastic for t in temperatures: - if isinstance(rx_endf.xs[t], IncoherentElastic): + if isinstance(rx_endf.xs[t], (IncoherentElastic, Sum)): rx.xs[t] = rx_endf.xs[t] rx.distribution[t] = rx_endf.distribution[t] @@ -832,20 +856,14 @@ class ThermalScattering(EqualityMixin): # 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: - # coherent elastic - + # Define helper functions to avoid duplication + def get_coherent_elastic(file_obj): # Get structure factor at first temperature params, S = endf.get_tab1_record(file_obj) strT = _temperature_str(params[0]) n_temps = params[2] bragg_edges = S.x - xs[strT] = CoherentElastic(bragg_edges, S.y) + xs = {strT: CoherentElastic(bragg_edges, S.y)} distribution = {strT: CoherentElasticAE(xs[strT])} # Get structure factor for subsequent temperatures @@ -854,15 +872,34 @@ class ThermalScattering(EqualityMixin): strT = _temperature_str(params[0]) xs[strT] = CoherentElastic(bragg_edges, S) distribution[strT] = CoherentElasticAE(xs[strT]) + return xs, distribution - elif lhtr == 2: - # incoherent elastic + def get_incoherent_elastic(file_obj): params, W = endf.get_tab1_record(file_obj) bound_xs = params[0] + xs = {} + distribution = {} for T, debye_waller in zip(W.x, W.y): strT = _temperature_str(T) xs[strT] = IncoherentElastic(bound_xs, debye_waller) distribution[strT] = IncoherentElasticAE(debye_waller) + return xs, distribution + + file_obj = StringIO(ev.section[7, 2]) + lhtr = endf.get_head_record(file_obj)[2] + if lhtr == 1: + # coherent elastic + xs, distribution = get_coherent_elastic(file_obj) + elif lhtr == 2: + # incoherent elastic + xs, distribution = get_incoherent_elastic(file_obj) + elif lhtr == 3: + # mixed coherent / incoherent elastic + xs_c, dist_c = get_coherent_elastic(file_obj) + xs_i, dist_i = get_incoherent_elastic(file_obj) + assert sorted(xs_c) == sorted(xs_i) + xs = {T: Sum([xs_c[T], xs_i[T]]) for T in xs_c} + distribution = {T: MixedElasticAE(dist_c[T], dist_i[T]) for T in dist_c} elastic = ThermalScatteringReaction(xs, distribution) diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index 05393e7bb..6ef874600 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -210,3 +210,60 @@ class IncoherentInelasticAEDiscrete(AngleEnergy): class IncoherentInelasticAE(CorrelatedAngleEnergy): _name = 'incoherent_inelastic' + + +class MixedElasticAE(AngleEnergy): + """Secondary distribution for mixed coherent/incoherent thermal elastic + + Parameters + ---------- + coherent : AngleEnergy + Secondary distribution for coherent elastic scattering + incoherent : AngleEnergy + Secondary distribution for incoherent elastic scattering + + Attributes + ---------- + coherent : AngleEnergy + Secondary distribution for coherent elastic scattering + incoherent : AngleEnergy + Secondary distribution for incoherent elastic scattering + + """ + def __init__(self, coherent, incoherent): + self.coherent = coherent + self.incoherent = incoherent + + def to_hdf5(self, group): + """Write mixed elastic distribution to an HDF5 group + + Parameters + ---------- + group : h5py.Group + HDF5 group to write to + + """ + group.attrs['type'] = np.string_('mixed_elastic') + coherent_group = group.create_group('coherent') + self.coherent.to_hdf5(coherent_group) + incoherent_group = group.create_group('incoherent') + self.incoherent.to_hdf5(incoherent_group) + + @classmethod + def from_hdf5(cls, group): + """Generate mixed thermal elastic distribution from HDF5 data + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + + Returns + ------- + openmc.data.MixedElasticAE + Mixed thermal elastic distribution + + """ + coherent = AngleEnergy.from_hdf5(group['coherent']) + incoherent = AngleEnergy.from_hdf5(group['incoherent']) + return cls(coherent, incoherent) diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index 3677f6ed3..e5e8274da 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -332,4 +332,39 @@ void IncoherentInelasticAE::sample( mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5); } +//============================================================================== +// MixedElasticAE implementation +//============================================================================== + +MixedElasticAE::MixedElasticAE( + hid_t group, const CoherentElasticXS& coh_xs, const Tabulated1D& incoh_xs) + : coherent_dist_(coh_xs), coherent_xs_(coh_xs), incoherent_xs_(incoh_xs) +{ + // Read incoherent elastic distribution + hid_t incoherent_group = open_group(group, "incoherent"); + std::string temp; + read_attribute(incoherent_group, "type", temp); + if (temp == "incoherent_elastic") { + incoherent_dist_ = make_unique(incoherent_group); + } else if (temp == "incoherent_elastic_discrete") { + incoherent_dist_ = + make_unique(incoherent_group, incoh_xs.x()); + } + close_group(incoherent_group); +} + +void MixedElasticAE::sample( + double E_in, double& E_out, double& mu, uint64_t* seed) const +{ + // Evaluate coherent and incoherent elastic cross sections + double xs_coh = coherent_xs_(E_in); + double xs_incoh = incoherent_xs_(E_in); + + if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) { + coherent_dist_.sample(E_in, E_out, mu, seed); + } else { + incoherent_dist_->sample(E_in, E_out, mu, seed); + } +} + } // namespace openmc diff --git a/src/thermal.cpp b/src/thermal.cpp index 1101d5485..303786915 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -210,14 +210,23 @@ ThermalData::ThermalData(hid_t group) if (temp == "coherent_elastic") { auto xs = dynamic_cast(elastic_.xs.get()); elastic_.distribution = make_unique(*xs); - } else { - if (temp == "incoherent_elastic") { - elastic_.distribution = make_unique(dgroup); - } else if (temp == "incoherent_elastic_discrete") { - auto xs = dynamic_cast(elastic_.xs.get()); - elastic_.distribution = - make_unique(dgroup, xs->x()); - } + } else if (temp == "incoherent_elastic") { + elastic_.distribution = make_unique(dgroup); + } else if (temp == "incoherent_elastic_discrete") { + auto xs = dynamic_cast(elastic_.xs.get()); + elastic_.distribution = + make_unique(dgroup, xs->x()); + } else if (temp == "mixed_elastic") { + // Get coherent/incoherent cross sections + auto mixed_xs = dynamic_cast(elastic_.xs.get()); + const auto& coh_xs = + dynamic_cast(mixed_xs->functions(0).get()); + const auto& incoh_xs = + dynamic_cast(mixed_xs->functions(1).get()); + + // Create mixed elastic distribution + elastic_.distribution = + make_unique(dgroup, *coh_xs, *incoh_xs); } close_group(elastic_group); From edf023adb9d48b5b8078ad3ca4886b0176457902 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 19 Oct 2021 17:05:49 -0500 Subject: [PATCH 05/10] Fix HDF5 export of mixed elastic distribution --- openmc/data/thermal_angle_energy.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index 6ef874600..d7dfe92f9 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -33,17 +33,22 @@ class CoherentElasticAE(AngleEnergy): def __init__(self, coherent_xs): self.coherent_xs = coherent_xs - def to_hdf5(self, group): + def to_hdf5(self, group, xs_group=None): """Write coherent elastic distribution to an HDF5 group Parameters ---------- group : h5py.Group HDF5 group to write to + xs_group : h5py.Group, optional + Group containing 'xs' dataset """ group.attrs['type'] = np.string_('coherent_elastic') - group['coherent_xs'] = group.parent['xs'] + if xs_group is not None: + group['coherent_xs'] = xs_group['xs'] + else: + group['coherent_xs'] = group.parent['xs'] class IncoherentElasticAE(AngleEnergy): @@ -245,7 +250,7 @@ class MixedElasticAE(AngleEnergy): """ group.attrs['type'] = np.string_('mixed_elastic') coherent_group = group.create_group('coherent') - self.coherent.to_hdf5(coherent_group) + self.coherent.to_hdf5(coherent_group, group.parent) incoherent_group = group.create_group('incoherent') self.incoherent.to_hdf5(incoherent_group) From 4e611f59b9ac60839fa0224a91be929572f0dbc3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 20 Oct 2021 06:53:20 -0500 Subject: [PATCH 06/10] Make incoherent XS polymorphic in MixedElasticAE --- include/openmc/secondary_thermal.h | 4 ++-- src/secondary_thermal.cpp | 5 +++-- src/thermal.cpp | 3 +-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 81de5c451..5b18902af 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -160,7 +160,7 @@ public: // //! \param[in] group HDF5 group explicit MixedElasticAE( - hid_t group, const CoherentElasticXS& coh_xs, const Tabulated1D& incoh_xs); + hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs); //! Sample distribution for an angle and energy //! \param[in] E_in Incoming energy in [eV] @@ -175,7 +175,7 @@ private: unique_ptr incoherent_dist_; //!< Incoherent distribution const CoherentElasticXS& coherent_xs_; //!< Ref. to coherent XS - const Tabulated1D& incoherent_xs_; //!< Ref. to incoherent XS + const Function1D& incoherent_xs_; //!< Polymorphic ref. to incoherent XS }; } // namespace openmc diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index e5e8274da..0b8e1ab42 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -337,7 +337,7 @@ void IncoherentInelasticAE::sample( //============================================================================== MixedElasticAE::MixedElasticAE( - hid_t group, const CoherentElasticXS& coh_xs, const Tabulated1D& incoh_xs) + hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs) : coherent_dist_(coh_xs), coherent_xs_(coh_xs), incoherent_xs_(incoh_xs) { // Read incoherent elastic distribution @@ -347,8 +347,9 @@ MixedElasticAE::MixedElasticAE( if (temp == "incoherent_elastic") { incoherent_dist_ = make_unique(incoherent_group); } else if (temp == "incoherent_elastic_discrete") { + auto xs = dynamic_cast(&incoh_xs); incoherent_dist_ = - make_unique(incoherent_group, incoh_xs.x()); + make_unique(incoherent_group, xs->x()); } close_group(incoherent_group); } diff --git a/src/thermal.cpp b/src/thermal.cpp index 303786915..1a9592d0a 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -221,8 +221,7 @@ ThermalData::ThermalData(hid_t group) auto mixed_xs = dynamic_cast(elastic_.xs.get()); const auto& coh_xs = dynamic_cast(mixed_xs->functions(0).get()); - const auto& incoh_xs = - dynamic_cast(mixed_xs->functions(1).get()); + const auto& incoh_xs = mixed_xs->functions(1).get(); // Create mixed elastic distribution elastic_.distribution = From e8294d52a0162b484a0961195dc330ef58b5f051 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 21 Jul 2022 14:37:34 -0500 Subject: [PATCH 07/10] Update documentation for mixed thermal elastic --- docs/source/io_formats/nuclear_data.rst | 23 +++++++++++++++++++++++ docs/source/pythonapi/data.rst | 1 + openmc/data/function.py | 4 ++++ openmc/data/thermal_angle_energy.py | 5 +++++ 4 files changed, 33 insertions(+) diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index ef61604cb..8174108e1 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -339,6 +339,16 @@ Incoherent elastic scattering [eV\ :math:`^{-1}`]. :Attributes: - **type** (*char[]*) -- 'IncoherentElastic' +Sum of functions +---------------- + +:Object type: Group +:Attributes: - **type** (*char[]*) -- "Sum" + - **n** (*int*) -- Number of functions +:Datasets: + - ***func_** (:ref:`function <1d_functions>`) -- Dataset for the + i-th function (indexing starts at 1) + .. _angle_energy: -------------------------- @@ -501,6 +511,19 @@ equiprobable bins. - **skewed** (*int8_t*) -- Whether discrete angles are equi-probable (0) or have a skewed distribution (1). +Mixed Elastic +------------- + +This angle-energy distribution is used when an evaluation specifies both +coherent and incoherent elastic thermal neutron scattering. + +:Object type: Group +:Attributes: - **type** (*char[]*) -- "mixed_elastic" +:Groups: - **coherent** -- Distribution for coherent elastic scattering. The + format is given in :ref:`angle_energy`. + - **incoherent** -- Distribution for incoherent elastic scattering. + The format is given in :ref:`angle_energy`. + .. _energy_distribution: -------------------- diff --git a/docs/source/pythonapi/data.rst b/docs/source/pythonapi/data.rst index 287738774..c74769be6 100644 --- a/docs/source/pythonapi/data.rst +++ b/docs/source/pythonapi/data.rst @@ -116,6 +116,7 @@ Angle-Energy Distributions IncoherentElasticAE IncoherentElasticAEDiscrete IncoherentInelasticAEDiscrete + MixedElasticAE Resonance Data -------------- diff --git a/openmc/data/function.py b/openmc/data/function.py index 7a73987df..751f6cdb5 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -581,6 +581,8 @@ class Sum(EqualityMixin): def to_hdf5(self, group, name='xy'): """Write sum of functions to an HDF5 group + .. versionadded:: 0.13.1 + Parameters ---------- group : h5py.Group @@ -599,6 +601,8 @@ class Sum(EqualityMixin): def from_hdf5(cls, group): """Generate sum of functions from an HDF5 group + .. versionadded:: 0.13.1 + Parameters ---------- group : h5py.Group diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index d7dfe92f9..b7ca349f1 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -36,6 +36,9 @@ class CoherentElasticAE(AngleEnergy): def to_hdf5(self, group, xs_group=None): """Write coherent elastic distribution to an HDF5 group + .. versionchanged:: 0.13.1 + The *xs_group* argument was added. + Parameters ---------- group : h5py.Group @@ -220,6 +223,8 @@ class IncoherentInelasticAE(CorrelatedAngleEnergy): class MixedElasticAE(AngleEnergy): """Secondary distribution for mixed coherent/incoherent thermal elastic + .. versionadded:: 0.13.1 + Parameters ---------- coherent : AngleEnergy From 0f03710e367094db1a822f00db7a1e857ea89dcc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 21 Jul 2022 16:09:49 -0500 Subject: [PATCH 08/10] Fix hdf5 write/read for coherent and mixed elastic --- openmc/data/function.py | 2 +- openmc/data/thermal.py | 2 +- openmc/data/thermal_angle_energy.py | 35 ++++++++++++++++++++--------- 3 files changed, 26 insertions(+), 13 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index 751f6cdb5..b0390d19c 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -544,7 +544,7 @@ class Combination(EqualityMixin): self._operations = operations -class Sum(EqualityMixin): +class Sum(Function1D): """Sum of multiple functions. This class allows you to create a callable object which represents the sum diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 6cf5a4e9e..3d5f2df64 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -701,7 +701,7 @@ class ThermalScattering(EqualityMixin): energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] coherent_xs = CoherentElastic(energy, P*EV_PER_MEV) - coherent_dist = CoherentElasticAE(xs) + coherent_dist = CoherentElasticAE(coherent_xs) # Coherent elastic shouldn't have angular distributions listed n_mu = ace.nxs[6] + 1 diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index b7ca349f1..1f37863ea 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -2,6 +2,7 @@ import numpy as np from .angle_energy import AngleEnergy from .correlated import CorrelatedAngleEnergy +import openmc.data class CoherentElasticAE(AngleEnergy): @@ -33,25 +34,37 @@ class CoherentElasticAE(AngleEnergy): def __init__(self, coherent_xs): self.coherent_xs = coherent_xs - def to_hdf5(self, group, xs_group=None): + def to_hdf5(self, group): """Write coherent elastic distribution to an HDF5 group - .. versionchanged:: 0.13.1 - The *xs_group* argument was added. - Parameters ---------- group : h5py.Group HDF5 group to write to - xs_group : h5py.Group, optional - Group containing 'xs' dataset """ group.attrs['type'] = np.string_('coherent_elastic') - if xs_group is not None: - group['coherent_xs'] = xs_group['xs'] - else: - group['coherent_xs'] = group.parent['xs'] + self.coherent_xs.to_hdf5(group, 'coherent_xs') + + @classmethod + def from_hdf5(cls, group): + """Generate coherent elastic distribution from HDF5 data + + .. versionadded:: 0.13.1 + + Parameters + ---------- + group : h5py.Group + HDF5 group to read from + + Returns + ------- + openmc.data.CoherentElasticAE + Coherent elastic distribution + + """ + coherent_xs = openmc.data.CoherentElastic.from_hdf5(group['coherent_xs']) + return cls(coherent_xs) class IncoherentElasticAE(AngleEnergy): @@ -255,7 +268,7 @@ class MixedElasticAE(AngleEnergy): """ group.attrs['type'] = np.string_('mixed_elastic') coherent_group = group.create_group('coherent') - self.coherent.to_hdf5(coherent_group, group.parent) + self.coherent.to_hdf5(coherent_group) incoherent_group = group.create_group('incoherent') self.incoherent.to_hdf5(incoherent_group) From 425dba5d25f6a3d737838f9cb408dd57bf45c24e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Jul 2022 07:46:34 -0500 Subject: [PATCH 09/10] Add test for mixed elastic thermal scattering --- openmc/mixin.py | 7 +- tests/unit_tests/test_data_thermal.py | 105 ++++++++++++++++++++++++++ 2 files changed, 110 insertions(+), 2 deletions(-) diff --git a/openmc/mixin.py b/openmc/mixin.py index 516162464..31c26ec76 100644 --- a/openmc/mixin.py +++ b/openmc/mixin.py @@ -14,8 +14,11 @@ class EqualityMixin: def __eq__(self, other): if isinstance(other, type(self)): for key, value in self.__dict__.items(): - if not np.array_equal(value, other.__dict__.get(key)): - return False + if isinstance(value, np.ndarray): + if not np.array_equal(value, other.__dict__.get(key)): + return False + else: + return value == other.__dict__.get(key) else: return False diff --git a/tests/unit_tests/test_data_thermal.py b/tests/unit_tests/test_data_thermal.py index dc4d24628..61ae61b24 100644 --- a/tests/unit_tests/test_data_thermal.py +++ b/tests/unit_tests/test_data_thermal.py @@ -262,3 +262,108 @@ def test_get_thermal_name(): # Names that don't remotely match anything assert f('boogie_monster') == 'c_boogie_monster' + + +@pytest.fixture +def fake_mixed_elastic(): + fake_tsl = openmc.data.ThermalScattering("c_D_in_7LiD", 1.9968, 4.9, [0.0253]) + fake_tsl.nuclides = ['H2'] + + # Create elastic reaction + bragg_edges = [0.00370672, 0.00494229, 0.00988458, 0.01359131, 0.01482688, + 0.01976918, 0.02347589, 0.02471147, 0.02965376, 0.03336048, + 0.03953834, 0.04324506, 0.04448063, 0.04942292, 0.05312964, + 0.05436522, 0.05930751, 0.06301423, 0.0642498 , 0.06919209, + 0.07289881, 0.07907667, 0.08278339, 0.08401896, 0.08896126, + 0.09266798, 0.09390355, 0.09884584, 0.1025526 , 0.1037882 , + 0.1087305 , 0.1124372 , 0.1186151 , 0.1223218 , 0.1235574 , + 0.1284997 , 0.1322064 , 0.133442 , 0.142091 , 0.1433266 , + 0.1482688 , 0.1519756 , 0.1581534 , 0.1618601 , 0.1630957 , + 0.168038 , 0.1717447 , 0.1729803 , 0.1779226 , 0.1816293 , + 0.1828649 , 0.1878072 , 0.1915139 , 0.1976918 , 0.2026341 , + 0.2075763 , 0.2125186 , 0.2174609 , 0.2224032 , 0.2273455 , + 0.2421724 , 0.2471147 , 0.252057 , 0.2569993 , 0.2619415 , + 0.2668838 , 0.2767684 , 0.2817107 , 0.2915953 , 0.3064222 , + 0.3261913 , 0.366965] + factors = [0.00375735, 0.01386287, 0.02595574, 0.02992438, 0.03549502, + 0.03855745, 0.04058831, 0.04986305, 0.05703106, 0.05855471, + 0.06078031, 0.06212291, 0.06656602, 0.06930339, 0.0697072 , + 0.07201456, 0.07263853, 0.07313129, 0.07465531, 0.07714482, + 0.07759976, 0.077809 , 0.07790282, 0.07927957, 0.08013058, + 0.08026637, 0.08073475, 0.08112202, 0.08123039, 0.08187171, + 0.08213756, 0.08218236, 0.08236572, 0.08240729, 0.08259795, + 0.08297893, 0.08300455, 0.08314566, 0.08315611, 0.08337715, + 0.08350026, 0.08350663, 0.08352815, 0.08353776, 0.0836098 , + 0.08367017, 0.08367361, 0.0837242 , 0.08375069, 0.08375227, + 0.08377006, 0.08381488, 0.08381644, 0.08382698, 0.08386266, + 0.08387756, 0.08388445, 0.08388974, 0.08390341, 0.08391088, + 0.08391695, 0.08392361, 0.08392684, 0.08392818, 0.08393161, + 0.08393546, 0.08393685, 0.08393801, 0.08393976, 0.08394167, + 0.08394288, 0.08394398] + coherent_xs = openmc.data.CoherentElastic(bragg_edges, factors) + incoherent_xs = openmc.data.Tabulated1D([0.00370672, 0.00370672], [0.00370672, 0.00370672]) + elastic_xs = {'294K': openmc.data.Sum((coherent_xs, incoherent_xs))} + coherent_dist = openmc.data.CoherentElasticAE(coherent_xs) + incoherent_dist = openmc.data.IncoherentElasticAEDiscrete([ + [-0.6, -0.18, 0.18, 0.6], [-0.6, -0.18, 0.18, 0.6] + ]) + elastic_dist = {'294K': openmc.data.MixedElasticAE(coherent_dist, incoherent_dist)} + fake_tsl.elastic = openmc.data.ThermalScatteringReaction(elastic_xs, elastic_dist) + + # Create inelastic reaction + inelastic_xs = {'294K': openmc.data.Tabulated1D([1.0e-5, 4.9], [13.4, 3.35])} + breakpoints = [3] + interpolation = [2] + energy = [1.0e-5, 4.3e-2, 4.9] + energy_out = [ + openmc.data.Tabular([0.0002, 0.067, 0.146, 0.366], [0.25, 0.25, 0.25, 0.25]), + openmc.data.Tabular([0.0001, 0.009, 0.137, 0.277], [0.25, 0.25, 0.25, 0.25]), + openmc.data.Tabular([0.0579, 4.555, 4.803, 4.874], [0.25, 0.25, 0.25, 0.25]), + ] + for eout in energy_out: + eout.normalize() + eout.c = eout.cdf() + discrete = openmc.stats.Discrete([-0.9, -0.6, -0.3, -0.1, 0.1, 0.3, 0.6, 0.9], [1/8]*8) + discrete.c = discrete.cdf()[1:] + mu = [[discrete]*4]*3 + inelastic_dist = {'294K': openmc.data.IncoherentInelasticAE( + breakpoints, interpolation, energy, energy_out, mu)} + inelastic = openmc.data.ThermalScatteringReaction(inelastic_xs, inelastic_dist) + fake_tsl.inelastic = inelastic + + return fake_tsl + + +def test_mixed_elastic(fake_mixed_elastic, run_in_tmpdir): + # Write data to HDF5 and then read back + original = fake_mixed_elastic + original.export_to_hdf5('c_D_in_7LiD.h5') + copy = openmc.data.ThermalScattering.from_hdf5('c_D_in_7LiD.h5') + + # Make sure data did not change as a result of HDF5 writing/reading + assert original == copy + + # Create modified cross_sections.xml file that includes the above data + xs = openmc.data.DataLibrary.from_xml() + xs.register_file('c_D_in_7LiD.h5') + xs.export_to_xml('cross_sections_mixed.xml') + + # Create a minimal model that includes the new data and run it + mat = openmc.Material() + mat.add_nuclide('H2', 1.0) + mat.add_nuclide('Li7', 1.0) + mat.set_density('g/cm3', 1.0) + mat.add_s_alpha_beta('c_D_in_7LiD') + sph = openmc.Sphere(r=10.0, boundary_type="vacuum") + cell = openmc.Cell(fill=mat, region=-sph) + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + model.materials = openmc.Materials([mat]) + model.materials.cross_sections = "cross_sections_mixed.xml" + model.settings.particles = 1000 + model.settings.batches = 10 + model.settings.run_mode = 'fixed source' + model.settings.source = openmc.Source( + energy=openmc.stats.Discrete([3.0], [1.0]) # 3 eV source + ) + model.run() From 0118067e2f7cc1fd37f50591f998b2f97ce897ae Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Jul 2022 16:00:50 -0500 Subject: [PATCH 10/10] Fix for IncoherentElasticAE.from_hdf5 --- openmc/data/thermal_angle_energy.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/data/thermal_angle_energy.py b/openmc/data/thermal_angle_energy.py index 1f37863ea..17a560092 100644 --- a/openmc/data/thermal_angle_energy.py +++ b/openmc/data/thermal_angle_energy.py @@ -122,7 +122,7 @@ class IncoherentElasticAE(AngleEnergy): Incoherent elastic distribution """ - return cls(group['debye_waller']) + return cls(group['debye_waller'][()]) class IncoherentElasticAEDiscrete(AngleEnergy):