mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 21:25:36 -04:00
Add ThermalData::sample routine (based on sab_scatter)
This commit is contained in:
Paul Romano
parent
5f3022989c
commit
aab7b7ca41
8 changed files with 384 additions and 44 deletions
|
|
@ -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<int, 6> 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
|
||||
|
|
|
|||
67
src/nuclide.h
Normal file
67
src/nuclide.h
Normal file
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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<double>(*root, "temperature_range");
|
||||
temperature_range[0] = range[0];
|
||||
temperature_range[1] = range[1];
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
245
src/thermal.cpp
245
src/thermal.cpp
|
|
@ -1,7 +1,7 @@
|
|||
#include "thermal.h"
|
||||
|
||||
#include <algorithm> // for sort, move
|
||||
#include <cmath> // for round
|
||||
#include <algorithm> // for sort, move, min, max
|
||||
#include <cmath> // for round, sqrt, fabs
|
||||
#include <sstream> // for stringstream
|
||||
|
||||
#include "xtensor/xarray.hpp"
|
||||
|
|
@ -37,12 +37,13 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector<double>& 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<double>& 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<double> 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
|
||||
|
|
|
|||
|
|
@ -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<double> kTs_; // temperatures in eV (k*T)
|
||||
std::vector<std::string> nuclides_; // List of valid nuclides
|
||||
int secondary_mode_; // secondary mode (equal/skewed/continuous)
|
||||
|
||||
// cross sections and distributions at each temperature
|
||||
std::vector<ThermalData> data_;
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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 <typename T>
|
||||
std::vector<T> get_node_array(pugi::xml_node node, const char* name)
|
||||
{
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue