From 14754b60f82f664a9fcca9aac379f0462e545762 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 9 Jul 2018 14:47:08 -0500 Subject: [PATCH] Convert energy distributions --- CMakeLists.txt | 1 + src/distribution_energy.cpp | 337 ++++++++++++++++++++++++++++++++++++ src/distribution_energy.h | 94 ++++++++++ 3 files changed, 432 insertions(+) create mode 100644 src/distribution_energy.cpp create mode 100644 src/distribution_energy.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 20b9d78d01..5cbc6be538 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -393,6 +393,7 @@ add_library(libopenmc SHARED src/tallies/trigger_header.F90 src/cell.cpp src/distribution.cpp + src/distribution_energy.cpp src/distribution_multi.cpp src/distribution_spatial.cpp src/endf.cpp diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp new file mode 100644 index 0000000000..cadeeb3df1 --- /dev/null +++ b/src/distribution_energy.cpp @@ -0,0 +1,337 @@ +#include "distribution_energy.h" + +#include // for max, min, copy, move +#include // for size_t +#include // for back_inserter + +#include "endf.h" +#include "hdf5_interface.h" +#include "math_functions.h" +#include "random_lcg.h" +#include "search.h" +#include "xtensor/xview.hpp" + +namespace openmc { + +//============================================================================== +// DiscretePhoton implementation +//============================================================================== + +DiscretePhoton::DiscretePhoton(hid_t group) +{ + read_attribute(group, "primary_flag", primary_flag_); + read_attribute(group, "energy", energy_); + read_attribute(group, "atomic_weight_ratio", A_); +} + +double DiscretePhoton::sample(double E) const +{ + if (primary_flag_ == 2) { + return energy_ + A_/(A_+ 1)*E; + } else { + return energy_; + } +} + +//============================================================================== +// LevelInelastic implementation +//============================================================================== + +LevelInelastic::LevelInelastic(hid_t group) +{ + read_attribute(group, "threshold", threshold_); + read_attribute(group, "mass_ratio", mass_ratio_); +} + +double LevelInelastic::sample(double E) const +{ + return mass_ratio_*(E - threshold_); +} + +//============================================================================== +// ContinuousTabular implementation +//============================================================================== + +ContinuousTabular::ContinuousTabular(hid_t group) +{ + // Open incoming energy dataset + hid_t dset = open_dataset(group, "energy"); + + // Get interpolation parameters + xt::xarray temp; + read_attribute(dset, "interpolation", temp); + + auto temp_b = xt::view(temp, 0); // view of breakpoints + auto temp_i = xt::view(temp, 1); // view of interpolation parameters + + std::copy(temp_b.begin(), temp_b.end(), std::back_inserter(breakpoints_)); + for (const auto i : temp_i) + interpolation_.push_back(int2interp(i)); + n_region_ = breakpoints_.size(); + + // Get incoming energies + read_dataset(dset, energy_); + std::size_t n_energy = energy_.size(); + close_dataset(dset); + + // Get outgoing energy distribution data + dset = open_dataset(group, "distribution"); + std::vector offsets; + std::vector interp; + std::vector n_discrete; + read_attribute(dset, "offsets", offsets); + read_attribute(dset, "interpolation", interp); + read_attribute(dset, "n_discrete_lines", n_discrete); + + xt::xarray eout; + read_dataset(dset, eout); + close_dataset(dset); + + for (int i = 0; i < n_energy; ++i) { + // Determine number of outgoing energies + int j = offsets[i]; + int n; + if (i < n_energy - 1) { + n = offsets[i+1] - j; + } else { + n = eout.shape()[1] - j; + } + + // Assign interpolation scheme and number of discrete lines + CTTable d; + d.interpolation = int2interp(interp[i]); + d.n_discrete = n_discrete[i]; + + // Copy data + d.e_out = xt::view(eout, 0, xt::range(j, j+n)); + d.p = xt::view(eout, 1, xt::range(j, j+n)); + + // To get answers that match ACE data, for now we still use the tabulated + // CDF values that were passed through to the HDF5 library. At a later + // time, we can remove the CDF values from the HDF5 library and + // reconstruct them using the PDF + if (true) { + d.c = xt::view(eout, 2, xt::range(j, j+n)); + } else { + // Calculate cumulative distribution function -- discrete portion + for (int k = 0; k < d.n_discrete; ++k) { + if (k == 0) { + d.c[k] = d.p[k]; + } else { + d.c[k] = d.c[k-1] + d.p[k]; + } + } + + // Continuous portion + for (int k = d.n_discrete; k < n; ++k) { + if (k == d.n_discrete) { + d.c[k] = d.c[k-1] + d.p[k]; + } else { + if (d.interpolation == Interpolation::histogram) { + d.c[k] = d.c[k-1] + d.p[k-1]*(d.e_out[k] - d.e_out[k-1]); + } else if (d.interpolation == Interpolation::lin_lin) { + d.c[k] = d.c[k-1] + 0.5*(d.p[k-1] + d.p[k]) * + (d.e_out[k] - d.e_out[k-1]); + } + } + } + + // Normalize density and distribution functions + d.p /= d.c[n - 1]; + d.c /= d.c[n - 1]; + } + + distribution_.push_back(std::move(d)); + } // incoming energies +} + +double ContinuousTabular::sample(double E) const +{ + // Read number of interpolation regions and incoming energies + bool histogram_interp; + if (n_region_ == 1) { + histogram_interp = (interpolation_[0] == Interpolation::histogram); + } else { + histogram_interp = false; + } + + // Find energy bin and calculate interpolation factor -- if the energy is + // outside the range of the tabulated energies, choose the first or last bins + auto n_energy_in = energy_.size(); + int i; + double r; + if (E < energy_[0]) { + i = 0; + r = 0.0; + } else if (E > energy_[n_energy_in - 1]) { + i = n_energy_in - 2; + r = 1.0; + } else { + i = lower_bound_index(energy_.begin(), energy_.end(), E); + r = (E - energy_[i]) / (energy_[i+1] - energy_[i]); + } + + // Sample between the ith and [i+1]th bin + int l; + if (histogram_interp) { + l = i; + } else { + l = r > prn() ? i + 1 : i; + } + + // Interpolation for energy E1 and EK + int n_energy_out = distribution_[i].e_out.size(); + double E_i_1 = distribution_[i].e_out[0]; + double E_i_K = distribution_[i].e_out[n_energy_out - 1]; + + n_energy_out = distribution_[i+1].e_out.size(); + double E_i1_1 = distribution_[i+1].e_out[0]; + double E_i1_K = distribution_[i+1].e_out[n_energy_out - 1]; + + double E_1 = E_i_1 + r*(E_i1_1 - E_i_1); + double E_K = E_i_K + r*(E_i1_K - E_i_K); + + // Determine outgoing energy bin + n_energy_out = distribution_[l].e_out.size(); + double r1 = prn(); + double c_k = distribution_[l].c[0]; + double c_k1; + int k; + for (k = 0; k < n_energy_out - 2; ++k) { + c_k1 = distribution_[l].c[k+1]; + if (r1 < c_k1) break; + c_k = c_k1; + } + + // Check to make sure 1 <= k <= NP - 1 + k = std::max(0, std::min(k, n_energy_out - 2)); + + double E_l_k = distribution_[l].e_out[k]; + double p_l_k = distribution_[l].p[k]; + double E_out; + if (distribution_[l].interpolation == Interpolation::histogram) { + // Histogram interpolation + if (p_l_k > 0.0) { + E_out = E_l_k + (r1 - c_k)/p_l_k; + } else { + E_out = E_l_k; + } + + } else if (distribution_[l].interpolation == Interpolation::lin_lin) { + // Linear-linear interpolation + double E_l_k1 = distribution_[l].e_out[k+1]; + double p_l_k1 = distribution_[l].p[k+1]; + + double frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k); + if (frac == 0.0) { + E_out = E_l_k + (r1 - c_k)/p_l_k; + } else { + E_out = E_l_k + (std::sqrt(std::max(0.0, p_l_k*p_l_k + + 2.0*frac*(r1 - c_k))) - p_l_k)/frac; + } + } + + // Now interpolate between incident energy bins i and i + 1 + if (!histogram_interp && n_energy_out > 1) { + if (l == i) { + return E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1); + } else { + return E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1); + } + } else { + return E_out; + } + +} + +//============================================================================== +// MaxwellEnergy implementation +//============================================================================== + +MaxwellEnergy::MaxwellEnergy(hid_t group) +{ + read_attribute(group, "u", u_); + hid_t dset = open_dataset(group, "theta"); + theta_ = Tabulated1D{dset}; + close_dataset(dset); +} + +double MaxwellEnergy::sample(double E) const +{ + // Get temperature corresponding to incoming energy + double theta = theta_(E); + + while (true) { + // Sample maxwell fission spectrum + double E_out = maxwell_spectrum_c(theta); + + // Accept energy based on restriction energy + if (E_out <= E - u_) return E_out; + } +} + +//============================================================================== +// Evaporation implementation +//============================================================================== + +Evaporation::Evaporation(hid_t group) +{ + read_attribute(group, "u", u_); + hid_t dset = open_dataset(group, "theta"); + theta_ = Tabulated1D{dset}; + close_dataset(dset); +} + +double Evaporation::sample(double E) const +{ + // Get temperature corresponding to incoming energy + double theta = theta_(E); + + double y = (E - u_)/theta; + double v = 1.0 - std::exp(-y); + + // Sample outgoing energy based on evaporation spectrum probability + // density function + double x; + while (true) { + x = -std::log((1.0 - v*prn())*(1.0 - v*prn())); + if (x <= y) break; + } + + return x*theta; +} + +//============================================================================== +// WattEnergy implementation +//============================================================================== + +WattEnergy::WattEnergy(hid_t group) +{ + // Read restriction energy + read_attribute(group, "u", u_); + + // Read tabulated functions + hid_t dset = open_dataset(group, "a"); + a_ = Tabulated1D{dset}; + close_dataset(dset); + dset = open_dataset(group, "b"); + b_ = Tabulated1D{dset}; + close_dataset(dset); +} + +double WattEnergy::sample(double E) const +{ + // Determine Watt parameters at incident energy + double a = a_(E); + double b = b_(E); + + while (true) { + // Sample energy-dependent Watt fission spectrum + double E_out = watt_spectrum_c(a, b); + + // Accept energy based on restriction energy + if (E_out <= E - u_) return E_out; + } +} + +} diff --git a/src/distribution_energy.h b/src/distribution_energy.h new file mode 100644 index 0000000000..288585f378 --- /dev/null +++ b/src/distribution_energy.h @@ -0,0 +1,94 @@ +#ifndef OPENMC_DISTRIBUTION_ENERGY_H +#define OPENMC_DISTRIBUTION_ENERGY_H + +#include + +#include "xtensor/xtensor.hpp" +#include "constants.h" +#include "endf.h" +#include "hdf5.h" + +namespace openmc { + +//=============================================================================== +//! Abstract class defining an energy distribution that is a function of the +//! incident energy of a projectile. Each derived type must implement a sample() +//! function that returns a sampled outgoing energy given an incoming energy +//=============================================================================== + +class EnergyDistribution { +public: + virtual double sample(double E) const = 0; + virtual ~EnergyDistribution() = default; +}; + +class DiscretePhoton : public EnergyDistribution { +public: + explicit DiscretePhoton(hid_t group); + double sample(double E) const; +private: + int primary_flag_; + double energy_; + double A_; +}; + +class LevelInelastic : public EnergyDistribution { +public: + explicit LevelInelastic(hid_t group); + double sample(double E) const; +private: + double threshold_; + double mass_ratio_; +}; + +class ContinuousTabular : public EnergyDistribution { +public: + explicit ContinuousTabular(hid_t group); + double sample(double E) const; +private: + struct CTTable { + Interpolation interpolation; + int n_discrete; + xt::xtensor e_out; + xt::xtensor p; + xt::xtensor c; + }; + + int n_region_; + std::vector breakpoints_; + std::vector interpolation_; + std::vector energy_; + std::vector distribution_; +}; + +class Evaporation : public EnergyDistribution { +public: + explicit Evaporation(hid_t group); + double sample(double E) const; +private: + Tabulated1D theta_; + double u_; +}; + +class MaxwellEnergy : public EnergyDistribution { +public: + explicit MaxwellEnergy(hid_t group); + double sample(double E) const; +private: + Tabulated1D theta_; + double u_; +}; + +class WattEnergy : public EnergyDistribution { +public: + explicit WattEnergy(hid_t group); + double sample(double E) const; +private: + Tabulated1D a_; + Tabulated1D b_; + double u_; +}; + +} + +#endif // OPENMC_DISTRIBUTION_ENERGY_H