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

Merge pull request #1136 from nelsonag/urr_cpp

Browse source 
Convert Unresolved Resonance Range code to C++
This commit is contained in:
Paul Romano 2018-12-11 06:30:46 -06:00 • committed by GitHub
commit ca0bcc922e
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
13 changed files with 396 additions and 375 deletions
Download patch file Download diff file Expand all files Collapse all files

2
CMakeLists.txt
View file

@ -361,7 +361,6 @@ add_library(libopenmc SHARED
src/timer_header.F90
src/tracking.F90
src/track_output.F90
src/urr_header.F90
src/vector_header.F90
src/volume_calc.F90
src/volume_header.F90
@ -417,6 +416,7 @@ add_library(libopenmc SHARED
src/position.cpp
src/progress_bar.cpp
src/pugixml/pugixml_c.cpp
src/urr.cpp
src/random_lcg.cpp
src/reaction.cpp
src/reaction_product.cpp

14
include/openmc/constants.h
View file

@ -257,12 +257,12 @@ constexpr int LIBRARY_PHOTON {3};
constexpr int LIBRARY_MULTIGROUP {4};
// Probability table parameters
constexpr int URR_CUM_PROB {1};
constexpr int URR_TOTAL {2};
constexpr int URR_ELASTIC {3};
constexpr int URR_FISSION {4};
constexpr int URR_N_GAMMA {5};
constexpr int URR_HEATING {6};
constexpr int URR_CUM_PROB {0};
constexpr int URR_TOTAL {1};
constexpr int URR_ELASTIC {2};
constexpr int URR_FISSION {3};
constexpr int URR_N_GAMMA {4};
constexpr int URR_HEATING {5};
// Maximum number of partial fission reactions
constexpr int PARTIAL_FISSION_MAX {4};
@ -424,7 +424,7 @@ constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE
// Interpolation rules
enum class Interpolation {
histogram, lin_lin, lin_log, log_lin, log_log
histogram = 1, lin_lin = 2, lin_log = 3, log_lin = 4, log_log = 5
};
// Run modes

27
include/openmc/hdf5_interface.h
View file

@ -272,6 +272,33 @@ void read_dataset(hid_t obj_id, const char* name, xt::xarray<T>& arr, bool indep
}
template <typename T, std::size_t N>
void read_dataset(hid_t obj_id, const char* name, xt::xtensor<T, N>& arr, bool indep=false)
{
// Open dataset and read array
hid_t dset = open_dataset(obj_id, name);
// Get shape of dataset
std::vector<hsize_t> hsize_t_shape = object_shape(dset);
close_dataset(dset);
// cast from hsize_t to size_t
std::vector<size_t> shape(hsize_t_shape.size());
for (int i = 0; i < shape.size(); i++) {
shape[i] = static_cast<size_t>(hsize_t_shape[i]);
}
// Allocate new xarray to read data into
xt::xarray<T> xarr(shape);
// Read data from the dataset
read_dataset(obj_id, name, xarr);
// Copy into xtensor
arr = xarr;
}
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)

118
include/openmc/nuclide.h
View file

@ -14,6 +14,7 @@
#include "openmc/endf.h"
#include "openmc/reaction.h"
#include "openmc/reaction_product.h"
#include "urr.h"
namespace openmc {
@ -23,58 +24,6 @@ namespace openmc {
constexpr double CACHE_INVALID {-1.0};
//===============================================================================
// Data for a nuclide
//===============================================================================
class Nuclide {
public:
// Types, aliases
using EmissionMode = ReactionProduct::EmissionMode;
struct EnergyGrid {
std::vector<int> grid_index;
std::vector<double> energy;
};
// Constructors
Nuclide(hid_t group, const double* temperature, int n);
// Methods
double nu(double E, EmissionMode mode, int group=0) const;
void calculate_elastic_xs(int i_nuclide) const;
//! Determines the microscopic 0K elastic cross section at a trial relative
//! energy used in resonance scattering
double elastic_xs_0K(double E) const;
// Data members
std::string name_; //! Name of nuclide, e.g. "U235"
int Z_; //! Atomic number
int A_; //! Mass number
int metastable_; //! Metastable state
double awr_; //! Atomic weight ratio
std::vector<double> kTs_; //! temperatures in eV (k*T)
std::vector<EnergyGrid> grid_; //! Energy grid at each temperature
bool fissionable_ {false}; //! Whether nuclide is fissionable
bool has_partial_fission_ {false}; //! has partial fission reactions?
std::vector<Reaction*> fission_rx_; //! Fission reactions
int n_precursor_ {0}; //! Number of delayed neutron precursors
std::unique_ptr<Function1D> total_nu_; //! Total neutron yield
// Resonance scattering information
bool resonant_ {false};
std::vector<double> energy_0K_;
std::vector<double> elastic_0K_;
std::vector<double> xs_cdf_;
std::vector<std::unique_ptr<Reaction>> reactions_; //! Reactions
std::vector<int> index_inelastic_scatter_;
private:
void create_derived();
};
//===============================================================================
//! Cached microscopic cross sections for a particular nuclide at the current
//! energy
@ -132,6 +81,68 @@ struct MaterialMacroXS {
double pair_production; //!< macroscopic pair production xs
};
//===============================================================================
// Data for a nuclide
//===============================================================================
class Nuclide {
public:
// Types, aliases
using EmissionMode = ReactionProduct::EmissionMode;
struct EnergyGrid {
std::vector<int> grid_index;
std::vector<double> energy;
};
// Constructors
Nuclide(hid_t group, const double* temperature, int n, int i_nuclide);
// Methods
double nu(double E, EmissionMode mode, int group=0) const;
void calculate_elastic_xs() const;
//! Determines the microscopic 0K elastic cross section at a trial relative
//! energy used in resonance scattering
double elastic_xs_0K(double E) const;
//! \brief Determines cross sections in the unresolved resonance range
//! from probability tables.
void calculate_urr_xs(int i_temp, double E);
// Data members
std::string name_; //! Name of nuclide, e.g. "U235"
int Z_; //! Atomic number
int A_; //! Mass number
int metastable_; //! Metastable state
double awr_; //! Atomic weight ratio
std::vector<double> kTs_; //! temperatures in eV (k*T)
std::vector<EnergyGrid> grid_; //! Energy grid at each temperature
int i_nuclide_; //! Index in the nuclides array
bool fissionable_ {false}; //! Whether nuclide is fissionable
bool has_partial_fission_ {false}; //! has partial fission reactions?
std::vector<Reaction*> fission_rx_; //! Fission reactions
int n_precursor_ {0}; //! Number of delayed neutron precursors
std::unique_ptr<Function1D> total_nu_; //! Total neutron yield
// Resonance scattering information
bool resonant_ {false};
std::vector<double> energy_0K_;
std::vector<double> elastic_0K_;
std::vector<double> xs_cdf_;
// Unresolved resonance range information
bool urr_present_ {false};
int urr_inelastic_ {C_NONE};
std::vector<UrrData> urr_data_;
std::vector<std::unique_ptr<Reaction>> reactions_; //! Reactions
std::vector<int> index_inelastic_scatter_;
private:
void create_derived();
};
//==============================================================================
// Global variables
//==============================================================================
@ -164,7 +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(bool use_mp, int i_nuclide,
int i_temp, double E);
} // namespace openmc

16
include/openmc/physics_mg.h
View file

@ -10,24 +10,18 @@
namespace openmc {
//TODO: Remove material_xs parameters when they reside on
// the C-side this should happen after materials, physics, input, and tallies
// are brought over
//! \brief samples particle behavior after a collision event.
//! \param p Particle to operate on
//! \param material_xs The cross section cache for the current material
extern "C" void
collision_mg(Particle* p, const MaterialMacroXS* material_xs);
collision_mg(Particle* p);
//! \brief samples a reaction type.
//!
//! Note that there is special logic when suvival biasing is turned on since
//! fission and disappearance are treated implicitly.
//! \param p Particle to operate on
//! \param material_xs The cross section cache for the current material
void
sample_reaction(Particle* p, const MaterialMacroXS* material_xs);
sample_reaction(Particle* p);
//! \brief Samples the scattering event
//! \param p Particle to operate on
@ -40,16 +34,14 @@ scatter(Particle* p);
//! \param bank_array The particle bank to populate
//! \param size_bank Number of particles currently in the bank
//! \param bank_array_size Allocated size of the bank
//! \param material_xs The cross section cache for the current material
void
create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
int64_t bank_array_size, const MaterialMacroXS* material_xs);
int64_t bank_array_size);
//! \brief Handles an absorption event
//! \param p Particle to operate on
//! \param material_xs The cross section cache for the current material
void
absorption(Particle* p, const MaterialMacroXS* material_xs);
absorption(Particle* p);
} // namespace openmc
#endif // OPENMC_PHYSICS_MG_H

33
include/openmc/urr.h Normal file
View file

@ -0,0 +1,33 @@
//! \brief UrrData information for the unresolved resonance treatment
#ifndef OPENMC_URR_H
#define OPENMC_URR_H
#include "xtensor/xtensor.hpp"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
namespace openmc {
//==============================================================================
//! UrrData contains probability tables for the unresolved resonance range.
//==============================================================================
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
//! \brief Load the URR data from the provided HDF5 group
explicit UrrData(hid_t group_id);
};
} // namespace openmc
#endif // OPENMC_URR_H

193
src/nuclide.cpp
View file

@ -5,6 +5,7 @@
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/message_passing.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/string_utils.h"
@ -33,7 +34,7 @@ MaterialMacroXS material_xs;
// Nuclide implementation
//==============================================================================
Nuclide::Nuclide(hid_t group, const double* temperature, int n)
Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide)
{
// Get name of nuclide from group, removing leading '/'
name_ = object_name(group).substr(1);
@ -42,6 +43,7 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n)
read_attribute(group, "A", A_);
read_attribute(group, "metastable", metastable_);
read_attribute(group, "atomic_weight_ratio", awr_);
i_nuclide_ = i_nuclide;
// Determine temperatures available
hid_t kT_group = open_group(group, "kTs");
@ -191,6 +193,47 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n)
}
close_group(rxs_group);
// Read unresolved resonance probability tables if present
if (object_exists(group, "urr")) {
urr_present_ = true;
urr_data_.reserve(temps_to_read.size());
for (int i = 0; i < temps_to_read.size(); i++) {
// Get temperature as a string
std::string temp_str {std::to_string(temps_to_read[i]) + "K"};
// Read probability tables for i-th temperature
hid_t urr_group = open_group(group, ("urr/" + temp_str).c_str());
urr_data_.emplace_back(urr_group);
close_group(urr_group);
// Check for negative values
if (xt::any(urr_data_[i].prob_ < 0.) && mpi::master) {
warning("Negative value(s) found on probability table for nuclide " +
name_ + " at " + temp_str);
}
}
// If the inelastic competition flag indicates that the inelastic cross
// section should be determined from a normal reaction cross section, we
// need to get the index of the reaction.
if (temps_to_read.size() > 0) {
if (urr_data_[0].inelastic_flag_ > 0) {
for (int i = 0; i < reactions_.size(); i++) {
if (reactions_[i]->mt_ == urr_data_[0].inelastic_flag_) {
urr_inelastic_ = i;
}
}
// Abort if no corresponding inelastic reaction was found
if (urr_inelastic_ == C_NONE) {
fatal_error("Could no find inelastic reaction specified on "
"unresolved resonance probability table.");
}
}
}
}
// Check for nu-total
if (object_exists(group, "total_nu")) {
// Read total nu data
@ -324,10 +367,10 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const
}
}
void Nuclide::calculate_elastic_xs(int i_nuclide) const
void Nuclide::calculate_elastic_xs() const
{
// Get temperature index, grid index, and interpolation factor
auto& micro = simulation::micro_xs[i_nuclide];
auto& micro = simulation::micro_xs[i_nuclide_];
int i_temp = micro.index_temp - 1;
int i_grid = micro.index_grid - 1;
double f = micro.interp_factor;
@ -361,6 +404,134 @@ 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(int i_temp, double E)
{
auto& micro = simulation::micro_xs[i_nuclide_];
micro.use_ptable = true;
// Create a shorthand for the URR data
const auto& urr = urr_data_[i_temp];
// Determine the energy table
int i_energy = 0;
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 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);
//TODO: to maintain the same random number stream as the Fortran code this
//replaces, the seed is set with i_nuclide_ + 1 instead of i_nuclide_
double r = future_prn(static_cast<int64_t>(i_nuclide_ + 1));
prn_set_stream(STREAM_TRACKING);
int i_low = 0;
while (urr.prob_(i_energy, URR_CUM_PROB, i_low) <= r) {++i_low;};
int i_up = 0;
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
double elastic = 0.;
double fission = 0.;
double capture = 0.;
double f;
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));
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));
// 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.)) {
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)));
} 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.)) {
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)));
} 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.)) {
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)));
} else {
capture = 0.;
}
}
// Determine the treatment of inelastic scattering
double inelastic = 0.;
if (urr.inelastic_flag_ != C_NONE) {
// get interpolation factor
f = micro.interp_factor;
// Determine inelastic scattering cross section
Reaction* rx = reactions_[urr_inelastic_].get();
int xs_index = micro.index_grid - rx->xs_[i_temp].threshold;
if (xs_index >= 0) {
inelastic = (1. - f) * rx->xs_[i_temp].value[xs_index] +
f * rx->xs_[i_temp].value[xs_index + 1];
}
}
// Multiply by smooth cross-section if needed
if (urr.multiply_smooth_) {
calculate_elastic_xs();
elastic *= micro.elastic;
capture *= (micro.absorption - micro.fission);
fission *= micro.fission;
}
// Check for negative values
if (elastic < 0.) {elastic = 0.;}
if (fission < 0.) {fission = 0.;}
if (capture < 0.) {capture = 0.;}
// Set elastic, absorption, fission, and total x/s. Note that the total x/s
// is calculated as a sum of partials instead of the table-provided value
micro.elastic = elastic;
micro.absorption = capture + fission;
micro.fission = fission;
micro.total = elastic + inelastic + capture + fission;
// Determine nu-fission cross-section
if (fissionable_) {
micro.nu_fission = nu(E, EmissionMode::total) * micro.fission;
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
@ -374,7 +545,8 @@ set_particle_energy_bounds(int particle, double E_min, double E_max)
extern "C" Nuclide* nuclide_from_hdf5_c(hid_t group, const double* temperature, int n)
{
data::nuclides.push_back(std::make_unique<Nuclide>(group, temperature, n));
data::nuclides.push_back(std::make_unique<Nuclide>(group, temperature, n,
data::nuclides.size()));
return data::nuclides.back().get();
}
@ -395,4 +567,17 @@ void set_micro_xs()
}
}
extern "C" void
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)) {
if ((E > nuc->urr_data_[i_temp - 1].energy_(0)) &&
(E < nuc->urr_data_[i_temp - 1].energy_(
nuc->urr_data_[i_temp - 1].n_energy_ - 1))) {
nuc->calculate_urr_xs(i_temp - 1, E);
}
}
}
} // namespace openmc

216
src/nuclide_header.F90
View file

@ -23,7 +23,6 @@ module nuclide_header
use stl_vector, only: VectorInt, VectorReal
use string
use simulation_header, only: need_depletion_rx
use urr_header, only: UrrData
implicit none
@ -76,11 +75,6 @@ module nuclide_header
integer, allocatable :: index_fission(:) ! indices in reactions
class(Function1D), allocatable :: total_nu
! Unresolved resonance data
logical :: urr_present = .false.
integer :: urr_inelastic
type(UrrData), allocatable :: urr_data(:)
! Multipole data
logical :: mp_present = .false.
type(MultipoleArray), pointer :: multipole => null()
@ -198,6 +192,14 @@ module nuclide_header
character(kind=C_CHAR), intent(in) :: name(*)
type(C_PTR) :: path
end function
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
end interface
contains
@ -253,7 +255,7 @@ contains
integer :: i
integer :: i_closest
integer :: n_temperature
integer(HID_T) :: urr_group, nu_group
integer(HID_T) :: nu_group
integer(HID_T) :: energy_group, energy_dset
integer(HID_T) :: kT_group
integer(HID_T) :: rxs_group
@ -438,47 +440,6 @@ contains
end do
call close_group(rxs_group)
! Read unresolved resonance probability tables if present
if (object_exists(group_id, 'urr')) then
this % urr_present = .true.
allocate(this % urr_data(n_temperature))
do i = 1, n_temperature
! Get temperature as a string
temp_str = trim(to_str(temps_to_read % data(i))) // "K"
! Read probability tables for i-th temperature
urr_group = open_group(group_id, 'urr/' // trim(temp_str))
call this % urr_data(i) % from_hdf5(urr_group)
call close_group(urr_group)
! Check for negative values
if (any(this % urr_data(i) % prob < ZERO) .and. master) then
call warning("Negative value(s) found on probability table &
&for nuclide " // this % name // " at " // trim(temp_str))
end if
end do
! if the inelastic competition flag indicates that the inelastic cross
! section should be determined from a normal reaction cross section, we
! need to get the index of the reaction
if (n_temperature > 0) then
if (this % urr_data(1) % inelastic_flag > 0) then
do i = 1, size(this % reactions)
if (this % reactions(i) % MT == this % urr_data(1) % inelastic_flag) then
this % urr_inelastic = i
end if
end do
! Abort if no corresponding inelastic reaction was found
if (this % urr_inelastic == NONE) then
call fatal_error("Could not find inelastic reaction specified on &
&unresolved resonance probability table.")
end if
end if
end if
end if
! Check for nu-total
if (object_exists(group_id, 'total_nu')) then
nu_group = open_group(group_id, 'total_nu')
@ -776,7 +737,7 @@ contains
real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b)
type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache
logical :: use_mp ! true if XS can be calculated with windowed multipole
logical(C_BOOL) :: use_mp ! true if XS can be calculated with windowed multipole
integer :: i_temp ! index for temperature
integer :: i_grid ! index on nuclide energy grid
integer :: i_low ! lower logarithmic mapping index
@ -831,7 +792,7 @@ contains
! 1. physics.F90 - scatter - For inelastic scatter.
! 2. physics.F90 - sample_fission - For partial fissions.
! 3. tally.F90 - score_general - For tallying on MTxxx reactions.
! 4. nuclide_header.F90 - calculate_urr_xs - For unresolved purposes.
! 4. nuclide.h - calculate_urr_xs - For unresolved purposes.
! It is worth noting that none of these occur in the resolved
! resonance range, so the value here does not matter. index_temp is
! set to -1 to force a segfault in case a developer messes up and tries
@ -974,15 +935,10 @@ contains
call calculate_sab_xs(this, i_sab, E, sqrtkT, sab_frac, micro_xs)
end if
! If the particle is in the unresolved resonance range and there are
! probability tables, we need to determine cross sections from the table
if (urr_ptables_on .and. this % urr_present .and. .not. use_mp) then
if (E > this % urr_data(i_temp) % energy(1) .and. E < this % &
urr_data(i_temp) % energy(this % urr_data(i_temp) % n_energy)) then
call calculate_urr_xs(this, i_temp, E, micro_xs)
end if
end if
call nuclide_calculate_urr_xs(use_mp, this % i_nuclide, i_temp, E)
micro_xs % last_E = E
micro_xs % last_sqrtkT = sqrtkT
@ -1246,152 +1202,6 @@ contains
end if
end subroutine multipole_deriv_eval
!===============================================================================
! CALCULATE_URR_XS determines cross sections in the unresolved resonance range
! from probability tables
!===============================================================================
subroutine calculate_urr_xs(this, i_temp, E, micro_xs)
class(Nuclide), intent(in) :: this ! Nuclide object
integer, intent(in) :: i_temp ! temperature index
real(8), intent(in) :: E ! energy
type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache
integer :: i_energy ! index for energy
integer :: i_low ! band index at lower bounding energy
integer :: i_up ! band index at upper bounding energy
integer :: threshold ! threshold energy index
real(8) :: f ! interpolation factor
real(8) :: r ! pseudo-random number
real(8) :: elastic ! elastic cross section
real(8) :: capture ! (n,gamma) cross section
real(8) :: fission ! fission cross section
real(8) :: inelastic ! inelastic cross section
micro_xs % use_ptable = .true.
associate (urr => this % urr_data(i_temp))
! determine energy table
i_energy = 1
do
if (E < urr % energy(i_energy + 1)) exit
i_energy = i_energy + 1
end do
! determine interpolation factor on table
f = (E - urr % energy(i_energy)) / &
(urr % energy(i_energy + 1) - urr % energy(i_energy))
! sample probability table using the cumulative distribution
! Random numbers for xs calculation are sampled from a separated stream.
! This guarantees the randomness and, at the same time, makes sure we reuse
! random number for the same nuclide at different temperatures, therefore
! preserving correlation of temperature in probability tables.
call prn_set_stream(STREAM_URR_PTABLE)
r = future_prn(int(this % i_nuclide, 8))
call prn_set_stream(STREAM_TRACKING)
i_low = 1
do
if (urr % prob(i_energy, URR_CUM_PROB, i_low) > r) exit
i_low = i_low + 1
end do
i_up = 1
do
if (urr % prob(i_energy + 1, URR_CUM_PROB, i_up) > r) exit
i_up = i_up + 1
end do
! determine elastic, fission, and capture cross sections from probability
! table
if (urr % interp == LINEAR_LINEAR) then
elastic = (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_low) + &
f * urr % prob(i_energy + 1, URR_ELASTIC, i_up)
fission = (ONE - f) * urr % prob(i_energy, URR_FISSION, i_low) + &
f * urr % prob(i_energy + 1, URR_FISSION, i_up)
capture = (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_low) + &
f * urr % prob(i_energy + 1, URR_N_GAMMA, i_up)
elseif (urr % interp == LOG_LOG) then
! Get logarithmic interpolation factor
f = log(E / urr % energy(i_energy)) / &
log(urr % energy(i_energy + 1) / urr % energy(i_energy))
! Calculate elastic cross section/factor
elastic = ZERO
if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then
elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, &
i_up)))
end if
! Calculate fission cross section/factor
fission = ZERO
if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then
fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, &
i_up)))
end if
! Calculate capture cross section/factor
capture = ZERO
if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then
capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, &
i_up)))
end if
end if
! Determine treatment of inelastic scattering
inelastic = ZERO
if (urr % inelastic_flag > 0) then
! Get index on energy grid and interpolation factor
i_energy = micro_xs % index_grid
f = micro_xs % interp_factor
! Determine inelastic scattering cross section
associate (rx => this % reactions(this % urr_inelastic))
threshold = rx % xs_threshold(i_temp)
if (i_energy >= threshold) then
inelastic = (ONE - f) * rx % xs(i_temp, i_energy - threshold + 1) + &
f * rx % xs(i_temp, i_energy - threshold + 2)
end if
end associate
end if
! Multiply by smooth cross-section if needed
if (urr % multiply_smooth) then
call this % calculate_elastic_xs(micro_xs)
elastic = elastic * micro_xs % elastic
capture = capture * (micro_xs % absorption - micro_xs % fission)
fission = fission * micro_xs % fission
end if
! Check for negative values
if (elastic < ZERO) elastic = ZERO
if (fission < ZERO) fission = ZERO
if (capture < ZERO) capture = ZERO
! Set elastic, absorption, fission, and total cross sections. Note that the
! total cross section is calculated as sum of partials rather than using the
! table-provided value
micro_xs % elastic = elastic
micro_xs % absorption = capture + fission
micro_xs % fission = fission
micro_xs % total = elastic + inelastic + capture + fission
! Determine nu-fission cross section
if (this % fissionable) then
micro_xs % nu_fission = this % nu(E, EMISSION_TOTAL) * &
micro_xs % fission
end if
end associate
end subroutine calculate_urr_xs
!===============================================================================
! CHECK_DATA_VERSION checks for the right version of nuclear data within HDF5
! files

2
src/physics.cpp
View file

@ -627,7 +627,7 @@ void scatter(Particle* p, int i_nuclide, int i_nuc_mat)
// Calculate elastic cross section if it wasn't precalculated
if (micro.elastic == CACHE_INVALID) {
nuc->calculate_elastic_xs(i_nuclide);
nuc->calculate_elastic_xs();
}
double prob = micro.elastic - micro.thermal;

40
src/physics_mg.cpp
View file

@ -22,13 +22,13 @@
namespace openmc {
void
collision_mg(Particle* p, const MaterialMacroXS* material_xs)
collision_mg(Particle* p)
{
// Add to the collision counter for the particle
p->n_collision++;
// Sample the reaction type
sample_reaction(p, material_xs);
sample_reaction(p);
// Display information about collision
if ((settings::verbosity >= 10) || (simulation::trace)) {
@ -39,7 +39,7 @@ collision_mg(Particle* p, const MaterialMacroXS* material_xs)
}
void
sample_reaction(Particle* p, const MaterialMacroXS* material_xs)
sample_reaction(Particle* p)
{
// Create fission bank sites. Note that while a fission reaction is sampled,
// it never actually "happens", i.e. the weight of the particle does not
@ -48,19 +48,20 @@ sample_reaction(Particle* p, const MaterialMacroXS* material_xs)
if (model::materials[p->material - 1]->fissionable) {
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
create_fission_sites(p, simulation::fission_bank.data(), &simulation::n_bank,
simulation::fission_bank.size(), material_xs);
create_fission_sites(
p, simulation::fission_bank.data(), &simulation::n_bank,
simulation::fission_bank.size());
} else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
(settings::create_fission_neutrons)) {
create_fission_sites(p, p->secondary_bank, &(p->n_secondary),
MAX_SECONDARY, material_xs);
MAX_SECONDARY);
}
}
// If survival biasing is being used, the following subroutine adjusts the
// weight of the particle. Otherwise, it checks to see if absorption occurs.
if (material_xs->absorption > 0.) {
absorption(p, material_xs);
if (simulation::material_xs.absorption > 0.) {
absorption(p);
} else {
p->absorb_wgt = 0.;
}
@ -102,7 +103,7 @@ scatter(Particle* p)
void
create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
int64_t bank_array_size, const MaterialMacroXS* material_xs)
int64_t bank_array_size)
{
// TODO: Heat generation from fission
@ -111,8 +112,8 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
// Determine the expected number of neutrons produced
double nu_t = p->wgt / simulation::keff * weight * material_xs->nu_fission /
material_xs->total;
double nu_t = p->wgt / simulation::keff * weight *
simulation::material_xs.nu_fission / simulation::material_xs.total;
// Sample the number of neutrons produced
int nu = static_cast<int>(nu_t);
@ -200,11 +201,12 @@ create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
}
void
absorption(Particle* p, const MaterialMacroXS* material_xs)
absorption(Particle* p)
{
if (settings::survival_biasing) {
// Determine weight absorbed in survival biasing
p->absorb_wgt = p->wgt * material_xs->absorption / material_xs->total;
p->absorb_wgt = p->wgt *
simulation::material_xs.absorption / simulation::material_xs.total;
// Adjust weight of particle by the probability of absorption
p->wgt -= p->absorb_wgt;
@ -212,13 +214,15 @@ absorption(Particle* p, const MaterialMacroXS* material_xs)
// Score implicit absorpion estimate of keff
#pragma omp atomic
global_tally_absorption += p->absorb_wgt * material_xs->nu_fission /
material_xs->absorption;
global_tally_absorption += p->absorb_wgt *
simulation::material_xs.nu_fission /
simulation::material_xs.absorption;
} else {
if (material_xs->absorption > prn() * material_xs->total) {
if (simulation::material_xs.absorption >
prn() * simulation::material_xs.total) {
#pragma omp atomic
global_tally_absorption += p->wgt * material_xs->nu_fission /
material_xs->absorption;
global_tally_absorption += p->wgt * simulation::material_xs.nu_fission /
simulation::material_xs.absorption;
p->alive = false;
p->event = EVENT_ABSORB;
}

7
src/tracking.F90
View file

@ -33,10 +33,9 @@ module tracking
implicit none
interface
subroutine collision_mg(p, material_xs) bind(C)
import Particle, C_DOUBLE, MaterialMacroXS
subroutine collision_mg(p) bind(C)
import Particle, C_DOUBLE
type(Particle), intent(inout) :: p
type(MaterialMacroXS), intent(in) :: material_xs
end subroutine collision_mg
end interface
@ -234,7 +233,7 @@ contains
if (run_CE) then
call collision(p)
else
call collision_mg(p, material_xs)
call collision_mg(p)
end if
! Score collision estimator tallies -- this is done after a collision

31
src/urr.cpp Normal file
View file

@ -0,0 +1,31 @@
#include "openmc/urr.h"
#include <iostream>
namespace openmc {
UrrData::UrrData(hid_t group_id)
{
// Read interpolation and other flags
int interp_temp;
read_attribute(group_id, "interpolation", interp_temp);
interp_ = static_cast<Interpolation>(interp_temp);
// read the metadata
read_attribute(group_id, "inelastic", inelastic_flag_);
read_attribute(group_id, "absorption", absorption_flag_);
int temp_multiply_smooth;
read_attribute(group_id, "multiply_smooth", temp_multiply_smooth);
multiply_smooth_ = (temp_multiply_smooth == 1);
// read the energies at which tables exist
read_dataset(group_id, "energy", energy_);
// Set n_energy_
n_energy_ = energy_.shape()[0];
// Read URR tables
read_dataset(group_id, "table", prob_);
}
}

72
src/urr_header.F90
View file

@ -1,72 +0,0 @@
module urr_header
use hdf5_interface, only: read_attribute, open_dataset, read_dataset, &
close_dataset, get_shape, HID_T, HSIZE_T
implicit none
!===============================================================================
! URRDATA contains probability tables for the unresolved resonance range.
!===============================================================================
type UrrData
integer :: n_energy ! # of incident neutron energies
integer :: n_prob ! # of probabilities
integer :: interp ! inteprolation (2=lin-lin, 5=log-log)
integer :: inelastic_flag ! inelastic competition flag
integer :: absorption_flag ! other absorption flag
logical :: multiply_smooth ! multiply by smooth cross section?
real(8), allocatable :: energy(:) ! incident energies
real(8), allocatable :: prob(:,:,:) ! actual probabibility tables
contains
procedure :: from_hdf5 => urr_from_hdf5
end type UrrData
contains
subroutine urr_from_hdf5(this, group_id)
class(UrrData), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
integer :: i, j, k
integer(HID_T) :: energy
integer(HID_T) :: table
integer(HSIZE_T) :: dims(1)
integer(HSIZE_T) :: dims3(3)
real(8), allocatable :: temp(:,:,:)
! Read interpolation and other flags
call read_attribute(this % interp, group_id, 'interpolation')
call read_attribute(this % inelastic_flag, group_id, 'inelastic')
call read_attribute(this % absorption_flag, group_id, 'absorption')
call read_attribute(i, group_id, 'multiply_smooth')
this % multiply_smooth = (i == 1)
! Read energies at which tables exist
energy = open_dataset(group_id, 'energy')
call get_shape(energy, dims)
this % n_energy = int(dims(1), 4)
allocate(this % energy(this % n_energy))
call read_dataset(this % energy, energy)
call close_dataset(energy)
! Read URR tables
table = open_dataset(group_id, 'table')
call get_shape(table, dims3)
this % n_prob = int(dims3(1), 4)
allocate(temp(this % n_prob, 6, this % n_energy))
call read_dataset(temp, table)
call close_dataset(table)
! Swap first and last indices
allocate(this % prob(this % n_energy, 6, this % n_prob))
do i = 1, this % n_energy
do j = 1, 6
do k = 1, this % n_prob
this % prob(i, j, k) = temp(k, j, i)
end do
end do
end do
end subroutine urr_from_hdf5
end module urr_header