crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 21:25:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Fixing @paulromano comments

Browse source
This commit is contained in:
Adam G Nelson 2018-12-07 20:38:24 -05:00
parent 999415be6a
commit af54b60c2e
6 changed files with 52 additions and 80 deletions
Download patch file Download diff file Expand all files Collapse all files

1
include/openmc/hdf5_interface.h
View file

@ -271,7 +271,6 @@ void read_dataset(hid_t obj_id, const char* name, xt::xarray<T>& arr, bool indep
close_dataset(dset);
}
template <typename T, std::size_t N>
void read_dataset_as_shape(hid_t obj_id, const char* name,
xt::xtensor<T, N>& arr, bool indep=false)

6
include/openmc/nuclide.h
View file

@ -107,7 +107,7 @@ public:
//! \brief Determines cross sections in the unresolved resonance range
//! from probability tables.
void calculate_urr_xs(const int i_temp, const double E);
void calculate_urr_xs(int i_temp, double E);
// Data members
std::string name_; //! Name of nuclide, e.g. "U235"
@ -175,8 +175,8 @@ extern "C" void set_micro_xs();
extern "C" bool nuclide_wmp_present(int i_nuclide);
extern "C" double nuclide_wmp_emin(int i_nuclide);
extern "C" double nuclide_wmp_emax(int i_nuclide);
extern "C" void nuclide_calculate_urr_xs(const bool use_mp, const int i_nuclide,
const int i_temp, const double E);
extern "C" void nuclide_calculate_urr_xs(bool use_mp, int i_nuclide,
int i_temp, double E);
} // namespace openmc

17
include/openmc/urr.h
View file

@ -16,17 +16,16 @@ namespace openmc {
class UrrData{
public:
Interpolation interp_; // interpolation type
int inelastic_flag_; // inelastic competition flag
int absorption_flag_; // other absorption flag
bool multiply_smooth_; // multiply by smooth cross section?
int n_energy_; // number of energy points
xt::xtensor<double, 1> energy_; // incident energies
xt::xtensor<double, 3> prob_; // Actual probability tables
Interpolation interp_; //!< interpolation type
int inelastic_flag_; //!< inelastic competition flag
int absorption_flag_; //!< other absorption flag
bool multiply_smooth_; //!< multiply by smooth cross section?
int n_energy_; //!< number of energy points
xt::xtensor<double, 1> energy_; //!< incident energies
xt::xtensor<double, 3> prob_; //!< Actual probability tables
//! \brief Load the URR data from the provided HDF5 group
void
from_hdf5(hid_t group_id);
explicit UrrData(hid_t group_id);
};
} // namespace openmc

81
src/nuclide.cpp
View file

@ -196,7 +196,7 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide)
// Read unresolved resonance probability tables if present
if (object_exists(group, "urr")) {
urr_present_ = true;
urr_data_.resize(temps_to_read.size());
urr_data_.reserve(temps_to_read.size());
for (int i = 0; i < temps_to_read.size(); i++) {
// Get temperature as a string
@ -204,7 +204,7 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide)
// Read probability tables for i-th temperature
hid_t urr_group = open_group(group, ("urr/" + temp_str).c_str());
urr_data_[i].from_hdf5(urr_group);
urr_data_.emplace_back(urr_group);
close_group(urr_group);
// Check for negative values
@ -404,25 +404,22 @@ double Nuclide::elastic_xs_0K(double E) const
return (1.0 - f)*elastic_0K_[i_grid] + f*elastic_0K_[i_grid + 1];
}
void Nuclide::calculate_urr_xs(const int i_temp, const double E)
void Nuclide::calculate_urr_xs(int i_temp, double E)
{
auto& micro = simulation::micro_xs[i_nuclide_];
micro.use_ptable = true;
// Create a shorthand for the URR data
UrrData* urr = &(urr_data_[i_temp]);
const auto& urr = urr_data_[i_temp];
// Determine the energy table
int i_energy = 0;
while (true) {
if (E < urr->energy_(i_energy + 1)) {break;}
i_energy++;
}
while(E >= urr.energy_(i_energy + 1)) {++i_energy;};
// Sample the probability table using the cumulative distribution
// Random nmbers for the xs calculation are sampled from a separate stream.
// This guarantees the rnadomness and, at the same time, makes sure we
// This guarantees the randomness and, at the same time, makes sure we
// reuse random numbers for the same nuclide at different temperatures,
// therefore preserving correlation of temperature in probability tables.
prn_set_stream(STREAM_URR_PTABLE);
@ -432,15 +429,10 @@ void Nuclide::calculate_urr_xs(const int i_temp, const double E)
prn_set_stream(STREAM_TRACKING);
int i_low = 0;
while (true) {
if (urr->prob_(i_energy, URR_CUM_PROB, i_low) > r) {break;}
i_low++;
}
while (urr.prob_(i_energy, URR_CUM_PROB, i_low) <= r) {++i_low;};
int i_up = 0;
while (true) {
if (urr->prob_(i_energy + 1, URR_CUM_PROB, i_up) > r) {break;}
i_up++;
}
while (urr.prob_(i_energy + 1, URR_CUM_PROB, i_up) <= r) {++i_up;};
// Determine elastic, fission, and capture cross sections from the
// probability table
@ -448,51 +440,51 @@ void Nuclide::calculate_urr_xs(const int i_temp, const double E)
double fission = 0.;
double capture = 0.;
double f;
if (urr->interp_ == Interpolation::lin_lin) {
if (urr.interp_ == Interpolation::lin_lin) {
// Determine the interpolation factor on the table
f = (E - urr->energy_(i_energy)) /
(urr->energy_(i_energy + 1) - urr->energy_(i_energy));
f = (E - urr.energy_(i_energy)) /
(urr.energy_(i_energy + 1) - urr.energy_(i_energy));
elastic = (1. - f) * urr->prob_(i_energy, URR_ELASTIC, i_low) +
f * urr->prob_(i_energy + 1, URR_ELASTIC, i_up);
fission = (1. - f) * urr->prob_(i_energy, URR_FISSION, i_low) +
f * urr->prob_(i_energy + 1, URR_FISSION, i_up);
capture = (1. - f) * urr->prob_(i_energy, URR_N_GAMMA, i_low) +
f * urr->prob_(i_energy + 1, URR_N_GAMMA, i_up);
} else if (urr->interp_ == Interpolation::log_log) {
elastic = (1. - f) * urr.prob_(i_energy, URR_ELASTIC, i_low) +
f * urr.prob_(i_energy + 1, URR_ELASTIC, i_up);
fission = (1. - f) * urr.prob_(i_energy, URR_FISSION, i_low) +
f * urr.prob_(i_energy + 1, URR_FISSION, i_up);
capture = (1. - f) * urr.prob_(i_energy, URR_N_GAMMA, i_low) +
f * urr.prob_(i_energy + 1, URR_N_GAMMA, i_up);
} else if (urr.interp_ == Interpolation::log_log) {
// Determine interpolation factor on the table
f = std::log(E / urr->energy_(i_energy)) /
std::log(urr->energy_(i_energy + 1) / urr->energy_(i_energy));
f = std::log(E / urr.energy_(i_energy)) /
std::log(urr.energy_(i_energy + 1) / urr.energy_(i_energy));
// Calculate the elastic cross section/factor
if ((urr->prob_(i_energy, URR_ELASTIC, i_low) > 0.) &&
(urr->prob_(i_energy + 1, URR_ELASTIC, i_up) > 0.)) {
if ((urr.prob_(i_energy, URR_ELASTIC, i_low) > 0.) &&
(urr.prob_(i_energy + 1, URR_ELASTIC, i_up) > 0.)) {
elastic =
std::exp((1. - f) *
std::log(urr->prob_(i_energy, URR_ELASTIC, i_low)) +
f * std::log(urr->prob_(i_energy + 1, URR_ELASTIC, i_up)));
std::log(urr.prob_(i_energy, URR_ELASTIC, i_low)) +
f * std::log(urr.prob_(i_energy + 1, URR_ELASTIC, i_up)));
} else {
elastic = 0.;
}
// Calculate the fission cross section/factor
if ((urr->prob_(i_energy, URR_FISSION, i_low) > 0.) &&
(urr->prob_(i_energy + 1, URR_FISSION, i_up) > 0.)) {
if ((urr.prob_(i_energy, URR_FISSION, i_low) > 0.) &&
(urr.prob_(i_energy + 1, URR_FISSION, i_up) > 0.)) {
fission =
std::exp((1. - f) *
std::log(urr->prob_(i_energy, URR_FISSION, i_low)) +
f * std::log(urr->prob_(i_energy + 1, URR_FISSION, i_up)));
std::log(urr.prob_(i_energy, URR_FISSION, i_low)) +
f * std::log(urr.prob_(i_energy + 1, URR_FISSION, i_up)));
} else {
fission = 0.;
}
// Calculate the capture cross section/factor
if ((urr->prob_(i_energy, URR_N_GAMMA, i_low) > 0.) &&
(urr->prob_(i_energy + 1, URR_N_GAMMA, i_up) > 0.)) {
if ((urr.prob_(i_energy, URR_N_GAMMA, i_low) > 0.) &&
(urr.prob_(i_energy + 1, URR_N_GAMMA, i_up) > 0.)) {
capture =
std::exp((1. - f) *
std::log(urr->prob_(i_energy, URR_N_GAMMA, i_low)) +
f * std::log(urr->prob_(i_energy + 1, URR_N_GAMMA, i_up)));
std::log(urr.prob_(i_energy, URR_N_GAMMA, i_low)) +
f * std::log(urr.prob_(i_energy + 1, URR_N_GAMMA, i_up)));
} else {
capture = 0.;
}
@ -500,7 +492,7 @@ void Nuclide::calculate_urr_xs(const int i_temp, const double E)
// Determine the treatment of inelastic scattering
double inelastic = 0.;
if (urr->inelastic_flag_ != C_NONE) {
if (urr.inelastic_flag_ != C_NONE) {
// get interpolation factor
f = micro.interp_factor;
@ -514,7 +506,7 @@ void Nuclide::calculate_urr_xs(const int i_temp, const double E)
}
// Multiply by smooth cross-section if needed
if (urr->multiply_smooth_) {
if (urr.multiply_smooth_) {
calculate_elastic_xs();
elastic *= micro.elastic;
capture *= (micro.absorption - micro.fission);
@ -576,8 +568,7 @@ void set_micro_xs()
}
extern "C" void
nuclide_calculate_urr_xs(const bool use_mp, const int i_nuclide,
const int i_temp, const double E)
nuclide_calculate_urr_xs(bool use_mp, int i_nuclide, int i_temp, double E)
{
Nuclide* nuc = data::nuclides[i_nuclide - 1].get();
if (settings::urr_ptables_on && (nuc->urr_present_ && !use_mp)) {

7
src/nuclide_header.F90
View file

@ -193,14 +193,13 @@ module nuclide_header
type(C_PTR) :: path
end function
subroutine nuclide_calculate_urr_xs_c(use_mp, i_nuclide, i_temp, E) &
bind(C, name='nuclide_calculate_urr_xs')
subroutine nuclide_calculate_urr_xs(use_mp, i_nuclide, i_temp, E) bind(C)
import C_BOOL, C_INT, C_DOUBLE
logical(C_BOOL), value, intent(in) :: use_mp
integer(C_INT), value, intent(in) :: i_nuclide
integer(C_INT), value, intent(in) :: i_temp
real(C_DOUBLE), value, intent(in) :: E
end subroutine nuclide_calculate_urr_xs_c
end subroutine nuclide_calculate_urr_xs
end interface
contains
@ -939,7 +938,7 @@ contains
! If the particle is in the unresolved resonance range and there are
! probability tables, we need to determine cross sections from the table
call nuclide_calculate_urr_xs_c(use_mp, this % i_nuclide, i_temp, E)
call nuclide_calculate_urr_xs(use_mp, this % i_nuclide, i_temp, E)
micro_xs % last_E = E
micro_xs % last_sqrtkT = sqrtkT

20
src/urr.cpp
View file

@ -4,28 +4,12 @@
namespace openmc {
void
UrrData::from_hdf5(hid_t group_id)
UrrData::UrrData(hid_t group_id)
{
// Read interpolation and other flags
int interp_temp;
read_attribute(group_id, "interpolation", interp_temp);
switch (interp_temp) {
case static_cast<int>(Interpolation::histogram):
interp_ = Interpolation::histogram;
break;
case static_cast<int>(Interpolation::lin_lin):
interp_ = Interpolation::lin_lin;
break;
case static_cast<int>(Interpolation::lin_log):
interp_ = Interpolation::lin_log;
break;
case static_cast<int>(Interpolation::log_lin):
interp_ = Interpolation::log_lin;
break;
case static_cast<int>(Interpolation::log_log):
interp_ = Interpolation::log_log;
}
interp_ = static_cast<Interpolation>(interp_temp);
// read the metadata
read_attribute(group_id, "inelastic", inelastic_flag_);