crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Move thermal scattering reaction data into ThermalData::Reaction type

Browse source
This commit is contained in:
Paul Romano 2019-05-21 13:56:44 -05:00
parent 3e4fe61673
commit fd1f01a911
8 changed files with 69 additions and 354 deletions
Download patch file Download diff file Expand all files Collapse all files

6
include/openmc/secondary_thermal.h
View file

@ -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:

45
include/openmc/thermal.h
View file

@ -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<std::string, int> 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<double, 1> e_out; //!< Outgoing energies
xt::xtensor<double, 1> e_out_pdf; //!< Probability density function
xt::xtensor<double, 1> e_out_cdf; //!< Cumulative distribution function
xt::xtensor<double, 2> mu; //!< Equiprobable angles at each outgoing energy
struct Reaction {
// Default constructor
Reaction() { }
// Data members
std::unique_ptr<Function1D> xs; //!< Cross section
std::unique_ptr<AngleEnergy> 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<double> inelastic_e_in_; //!< incoming E grid for inelastic
std::vector<double> inelastic_sigma_; //!< inelastic scattering cross section
// The following are used only for equal/skewed distributions
xt::xtensor<double, 2> inelastic_e_out_;
xt::xtensor<double, 3> 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<DistEnergySab> 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<double> elastic_e_in_; //!< incoming E grid for elastic
std::vector<double> elastic_P_; //!< elastic scattering cross section
xt::xtensor<double, 2> elastic_mu_; //!< equi-probable angles at each incoming E
Reaction elastic_;
Reaction inelastic_;
// ThermalScattering needs access to private data members
friend class ThermalScattering;

2
openmc/data/thermal.py
View file

@ -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

2
openmc/data/thermal_angle_energy.py
View file

@ -138,4 +138,4 @@ class IncoherentInelasticAEDiscrete(AngleEnergy):
class IncoherentInelasticAEContinuous(CorrelatedAngleEnergy):
_name = 'incoherent_inelastic_continuous'
_name = 'incoherent_inelastic'

4
src/endf.cpp
View file

@ -87,6 +87,10 @@ read_function(hid_t group, const char* name)
func = std::make_unique<Tabulated1D>(dset);
} else if (func_type == "Polynomial") {
func = std::make_unique<Polynomial>(dset);
} else if (func_type == "CoherentElastic") {
func = std::make_unique<CoherentElasticXS>(dset);
} else if (func_type == "IncoherentElastic") {
func = std::make_unique<IncoherentElasticXS>(dset);
} else {
throw std::runtime_error{"Unknown function type " + func_type +
" for dataset " + object_name(dset)};

10
src/reaction_product.cpp
View file

@ -3,6 +3,7 @@
#include <memory> // for unique_ptr
#include <string> // 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<Function1D>{new Tabulated1D{yield}};
} else if (temp == "Polynomial") {
yield_ = std::unique_ptr<Function1D>{new Polynomial{yield}};
}
close_dataset(yield);
yield_ = read_function(group, "yield");
int n;
read_attribute(group, "n_distribution", n);

6
src/secondary_thermal.cpp
View file

@ -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;

348
src/thermal.cpp
View file

@ -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<double>& 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<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);
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<double> 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<CoherentElasticXS*>(elastic_.xs.get());
elastic_.distribution = std::make_unique<CoherentElasticAE>(*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<Tabulated1D*>(elastic_.xs.get());
if (temp == "incoherent_elastic") {
elastic_.distribution = std::make_unique<IncoherentElasticAE>(dgroup);
} else if (temp == "incoherent_elastic_discrete") {
elastic_.distribution = std::make_unique<IncoherentElasticAEDiscrete>(
dgroup, xs->x()
);
}
// Set elastic threshold
threshold_elastic_ = elastic_e_in_.back();
// Read angle distribution
if (elastic_mode_ == SAB_ELASTIC_INCOHERENT) {
xt::xarray<double> 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<double> 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<Tabulated1D*>(inelastic_.xs.get());
threshold_inelastic_ = xs->x().back();
if (secondary_mode != SAB_SECONDARY_CONT) {
// Read energy distribution
xt::xarray<double> 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<double> 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<Tabular*>(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<double>({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<IncoherentInelasticAE>(dgroup);
} else if (temp == "incoherent_inelastic_discrete") {
inelastic_.distribution = std::make_unique<IncoherentInelasticAEDiscrete>(
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()