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()