diff --git a/src/constants.h b/src/constants.h index b6bb82d37..12dc21b71 100644 --- a/src/constants.h +++ b/src/constants.h @@ -130,17 +130,6 @@ constexpr int ANGLE_HISTOGRAM {5}; constexpr int TEMPERATURE_NEAREST {1}; constexpr int TEMPERATURE_INTERPOLATION {2}; -// Secondary energy mode for S(a,b) inelastic scattering -// TODO: Convert to enum -constexpr int SAB_SECONDARY_EQUAL {0}; // Equally-likely outgoing energy bins -constexpr int SAB_SECONDARY_SKEWED {1}; // Skewed outgoing energy bins -constexpr int SAB_SECONDARY_CONT {2}; // Continuous, linear-linear interpolation - -// Elastic mode for S(a,b) elastic scattering -// TODO: Convert to enum -constexpr int SAB_ELASTIC_DISCRETE {3}; // Sample from discrete cosines -constexpr int SAB_ELASTIC_EXACT {4}; // Exact treatment for coherent elastic - // Reaction types // TODO: Convert to enum constexpr int TOTAL_XS {1}; @@ -260,6 +249,8 @@ constexpr int N_AC {849}; constexpr int N_2N0 {875}; constexpr int N_2NC {891}; +constexpr std::array DEPLETION_RX {N_GAMMA, N_P, N_A, N_2N, N_3N, N_4N}; + // Fission neutron emission (nu) type constexpr int NU_NONE {0}; // No nu values (non-fissionable) constexpr int NU_POLYNOMIAL {1}; // Nu values given by polynomial diff --git a/src/nuclide.h b/src/nuclide.h new file mode 100644 index 000000000..8c43d55a0 --- /dev/null +++ b/src/nuclide.h @@ -0,0 +1,67 @@ +#ifndef OPENMC_NUCLIDE_H +#define OPENMC_NUCLIDE_H + +#include "constants.h" + +namespace openmc { + +//=============================================================================== +//! Cached microscopic cross sections for a particular nuclide at the current +//! energy +//=============================================================================== + +struct NuclideMicroXS { + // Microscopic cross sections in barns + double total; //!< total cross section + double absorption; //!< absorption (disappearance) + double fission; //!< fission + double nu_fission; //!< neutron production from fission + + double elastic; //!< If sab_frac is not 1 or 0, then this value is + //!< averaged over bound and non-bound nuclei + double thermal; //!< Bound thermal elastic & inelastic scattering + double thermal_elastic; //!< Bound thermal elastic scattering + double photon_prod; //!< microscopic photon production xs + + // Cross sections for depletion reactions (note that these are not stored in + // macroscopic cache) + double reaction[DEPLETION_RX.size()]; + + // Indicies and factors needed to compute cross sections from the data tables + int index_grid; //!< Index on nuclide energy grid + int index_temp; //!< Temperature index for nuclide + double interp_factor; //!< Interpolation factor on nuc. energy grid + int index_sab {-1}; //!< Index in sab_tables + int index_temp_sab; //!< Temperature index for sab_tables + double sab_frac; //!< Fraction of atoms affected by S(a,b) + bool use_ptable; //!< In URR range with probability tables? + + // Energy and temperature last used to evaluate these cross sections. If + // these values have changed, then the cross sections must be re-evaluated. + double last_E {0.0}; //!< Last evaluated energy + double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant + //!< * temperature (eV)) +}; + +//=============================================================================== +// MATERIALMACROXS contains cached macroscopic cross sections for the material a +// particle is traveling through +//=============================================================================== + + struct MaterialMacroXS { + double total; //!< macroscopic total xs + double absorption; //!< macroscopic absorption xs + double fission; //!< macroscopic fission xs + double nu_fission; //!< macroscopic production xs + double photon_prod; //!< macroscopic photon production xs + + // Photon cross sections + double coherent; //!< macroscopic coherent xs + double incoherent; //!< macroscopic incoherent xs + double photoelectric; //!< macroscopic photoelectric xs + double pair_production; //!< macroscopic pair production xs + }; + +} // namespace openmc + +#endif // OPENMC_NUCLIDE_H diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 8cd225bc9..6faaaec7f 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -126,36 +126,36 @@ module nuclide_header ! (NuclideMicroXS % elastic) real(8), parameter :: CACHE_INVALID = dble(Z"FFE0000000000000") - type NuclideMicroXS + type, bind(C) :: NuclideMicroXS ! Microscopic cross sections in barns - real(8) :: total - real(8) :: absorption ! absorption (disappearance) - real(8) :: fission ! fission - real(8) :: nu_fission ! neutron production from fission + real(C_DOUBLE) :: total + real(C_DOUBLE) :: absorption ! absorption (disappearance) + real(C_DOUBLE) :: fission ! fission + real(C_DOUBLE) :: nu_fission ! neutron production from fission - real(8) :: elastic ! If sab_frac is not 1 or 0, then this value is + real(C_DOUBLE) :: elastic ! If sab_frac is not 1 or 0, then this value is ! averaged over bound and non-bound nuclei - real(8) :: thermal ! Bound thermal elastic & inelastic scattering - real(8) :: thermal_elastic ! Bound thermal elastic scattering - real(8) :: photon_prod ! microscopic photon production xs + real(C_DOUBLE) :: thermal ! Bound thermal elastic & inelastic scattering + real(C_DOUBLE) :: thermal_elastic ! Bound thermal elastic scattering + real(C_DOUBLE) :: photon_prod ! microscopic photon production xs ! Cross sections for depletion reactions (note that these are not stored in ! macroscopic cache) - real(8) :: reaction(size(DEPLETION_RX)) + real(C_DOUBLE) :: reaction(size(DEPLETION_RX)) ! Indicies and factors needed to compute cross sections from the data tables - integer :: index_grid ! Index on nuclide energy grid - integer :: index_temp ! Temperature index for nuclide - real(8) :: interp_factor ! Interpolation factor on nuc. energy grid - integer :: index_sab = NONE ! Index in sab_tables - integer :: index_temp_sab ! Temperature index for sab_tables - real(8) :: sab_frac ! Fraction of atoms affected by S(a,b) - logical :: use_ptable ! In URR range with probability tables? + integer(C_INT) :: index_grid ! Index on nuclide energy grid + integer(C_INT) :: index_temp ! Temperature index for nuclide + real(C_DOUBLE) :: interp_factor ! Interpolation factor on nuc. energy grid + integer(C_INT) :: index_sab = NONE ! Index in sab_tables + integer(C_INT) :: index_temp_sab ! Temperature index for sab_tables + real(C_DOUBLE) :: sab_frac ! Fraction of atoms affected by S(a,b) + logical(C_BOOL) :: use_ptable ! In URR range with probability tables? ! Energy and temperature last used to evaluate these cross sections. If ! these values have changed, then the cross sections must be re-evaluated. - real(8) :: last_E = ZERO ! Last evaluated energy - real(8) :: last_sqrtkT = ZERO ! Last temperature in sqrt(Boltzmann + real(C_DOUBLE) :: last_E = ZERO ! Last evaluated energy + real(C_DOUBLE) :: last_sqrtkT = ZERO ! Last temperature in sqrt(Boltzmann ! constant * temperature (eV)) end type NuclideMicroXS @@ -164,7 +164,7 @@ module nuclide_header ! particle is traveling through !=============================================================================== - type MaterialMacroXS + type, bind(C) :: MaterialMacroXS real(C_DOUBLE) :: total ! macroscopic total xs real(C_DOUBLE) :: absorption ! macroscopic absorption xs real(C_DOUBLE) :: fission ! macroscopic fission xs diff --git a/src/settings.cpp b/src/settings.cpp index 8ca07a527..c705ef8a7 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -72,6 +72,32 @@ void read_settings(pugi::xml_node* root) } } } + + // Get temperature settings + if (check_for_node(*root, "temperature_default")) { + temperature_default = std::stod(get_node_value(*root, "temperature_default")); + } + if (check_for_node(*root, "temperature_method")) { + auto temp_str = get_node_value(*root, "temperature_method", true, true); + if (temp_str == "nearest") { + temperature_method = TEMPERATURE_NEAREST; + } else if (temp_str == "interpolation") { + temperature_method = TEMPERATURE_INTERPOLATION; + } else { + fatal_error("Unknown temperature method: " + temp_str); + } + } + if (check_for_node(*root, "temperature_tolerance")) { + temperature_tolerance = std::stod(get_node_value(*root, "temperature_tolerance")); + } + if (check_for_node(*root, "temperature_multipole")) { + temperature_multipole = get_node_value_bool(*root, "temperature_multipole"); + } + if (check_for_node(*root, "temperature_range")) { + auto range = get_node_array(*root, "temperature_range"); + temperature_range[0] = range[0]; + temperature_range[1] = range[1]; + } } } // namespace openmc diff --git a/src/thermal.cpp b/src/thermal.cpp index b4818d080..29f8c83a4 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -1,7 +1,7 @@ #include "thermal.h" -#include // for sort, move -#include // for round +#include // for sort, move, min, max +#include // for round, sqrt, fabs #include // for stringstream #include "xtensor/xarray.hpp" @@ -37,12 +37,13 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector& tem 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; + secondary_mode = SAB_SECONDARY_EQUAL; } else if (sec_mode == "skewed") { - secondary_mode_ = SAB_SECONDARY_SKEWED; + secondary_mode = SAB_SECONDARY_SKEWED; } else if (sec_mode == "continuous") { - secondary_mode_ = SAB_SECONDARY_CONT; + secondary_mode = SAB_SECONDARY_CONT; } // Read temperatures @@ -138,7 +139,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, secondary_mode); close_group(group); } @@ -202,7 +203,7 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp, auto& E_in = sab.elastic_e_in_; - if (sab.elastic_mode_ == SAB_ELASTIC_EXACT) { + 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 @@ -237,6 +238,7 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp, //============================================================================== ThermalData::ThermalData(hid_t group, int secondary_mode) + : inelastic_mode_{secondary_mode} { // Coherent elastic data if (object_exists(group, "elastic")) { @@ -257,13 +259,13 @@ ThermalData::ThermalData(hid_t group, int secondary_mode) std::copy(P.begin(), P.end(), std::back_inserter(elastic_P_)); n_elastic_e_in_ = elastic_e_in_.size(); - // Determine elastic type + // Determine whether elastic scattering is incoherent or coherent std::string type; read_attribute(dset, "type", type); if (type == "tab1") { - elastic_mode_ = SAB_ELASTIC_DISCRETE; + elastic_mode_ = SAB_ELASTIC_INCOHERENT; } else if (type == "bragg") { - elastic_mode_ = SAB_ELASTIC_EXACT; + elastic_mode_ = SAB_ELASTIC_COHERENT; } close_dataset(dset); @@ -271,7 +273,7 @@ ThermalData::ThermalData(hid_t group, int secondary_mode) threshold_elastic_ = elastic_e_in_.back(); // Read angle distribution - if (elastic_mode_ != SAB_ELASTIC_EXACT) { + if (elastic_mode_ == SAB_ELASTIC_INCOHERENT) { xt::xarray mu_out; read_dataset(elastic_group, "mu_out", mu_out); elastic_mu_ = mu_out; @@ -350,4 +352,225 @@ ThermalData::ThermalData(hid_t group, int secondary_mode) } } +void +ThermalData::sample(const NuclideMicroXS* micro_xs, double E, + double* E_out, double* mu) +{ + // 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; + if (prob < elastic_P_.front()) { + k = 1; + } else { + 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 -- no need to do anything + + } 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); + + // Determine endpoints on grid i + 1 + n = inelastic_data_[i].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); + + 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 - 2; ++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 - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1); + } else { + *E_out = E_1 + (E - 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 - *mu; + } else { + mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1)) - *mu; + } + + // Smear angle + *mu += std::min(*mu - mu_left, mu_right - *mu)*(prn() - 0.5); + + } // (inelastic secondary energy treatment) + } // (elastic or inelastic) + + // 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::fabs(*mu) > 1.0) *mu = std::copysign(1.0, *mu); + +} + } // namespace openmc diff --git a/src/thermal.h b/src/thermal.h index 59f5c97ee..922a897b5 100644 --- a/src/thermal.h +++ b/src/thermal.h @@ -8,12 +8,28 @@ #include "xtensor/xtensor.hpp" #include "hdf5_interface.h" +#include "nuclide.h" namespace openmc { +// Secondary energy mode for S(a,b) inelastic scattering +// TODO: Convert to enum +constexpr int SAB_SECONDARY_EQUAL {0}; // Equally-likely outgoing energy bins +constexpr int SAB_SECONDARY_SKEWED {1}; // Skewed outgoing energy bins +constexpr int SAB_SECONDARY_CONT {2}; // Continuous, linear-linear interpolation + +// Elastic mode for S(a,b) elastic scattering +// TODO: Convert to enum +constexpr int SAB_ELASTIC_INCOHERENT {3}; // Incoherent elastic scattering +constexpr int SAB_ELASTIC_COHERENT {4}; // Coherent elastic scattering (Bragg edges) + class ThermalData { public: ThermalData(hid_t group, int secondary_mode); + + // Sample an outgoing energy and angle + void sample(const NuclideMicroXS* micro_xs, double E_in, + double* E_out, double* mu); private: struct DistEnergySab { std::size_t n_e_out; @@ -28,6 +44,7 @@ private: double threshold_elastic_ {0.0}; // Inelastic scattering data + int inelastic_mode_; // secondary mode (equal/skewed/continuous) std::size_t n_inelastic_e_in_; // # of incoming E for inelastic std::size_t n_inelastic_e_out_; // # of outgoing E for inelastic std::size_t n_inelastic_mu_; // # of outgoing angles for inelastic @@ -67,7 +84,6 @@ public: double awr_; // weight of nucleus in neutron masses std::vector kTs_; // temperatures in eV (k*T) std::vector nuclides_; // List of valid nuclides - int secondary_mode_; // secondary mode (equal/skewed/continuous) // cross sections and distributions at each temperature std::vector data_; diff --git a/src/xml_interface.cpp b/src/xml_interface.cpp index c89e19ef1..5835c8600 100644 --- a/src/xml_interface.cpp +++ b/src/xml_interface.cpp @@ -40,4 +40,19 @@ get_node_value(pugi::xml_node node, const char* name, bool lowercase, return value; } +bool +get_node_value_bool(pugi::xml_node node, const char* name) +{ + if (node.attribute(name)) { + return node.attribute(name).as_bool(); + } else if (node.child(name)) { + return node.child(name).text().as_bool(); + } else { + std::stringstream err_msg; + err_msg << "Node \"" << name << "\" is not a member of the \"" + << node.name() << "\" XML node"; + fatal_error(err_msg); + } +} + } // namespace openmc diff --git a/src/xml_interface.h b/src/xml_interface.h index f10143697..1530333c4 100644 --- a/src/xml_interface.h +++ b/src/xml_interface.h @@ -16,9 +16,11 @@ check_for_node(pugi::xml_node node, const char *name) return node.attribute(name) || node.child(name); } -std::string get_node_value(pugi::xml_node node, const char *name, +std::string get_node_value(pugi::xml_node node, const char* name, bool lowercase=false, bool strip=false); +bool get_node_value_bool(pugi::xml_node node, const char* name); + template std::vector get_node_array(pugi::xml_node node, const char* name) {