mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Got URR working, now just need a few simplifications of the Fortran code
This commit is contained in:
Adam G Nelson
parent
12540cbc5b
commit
5eab6dbed0
7 changed files with 262 additions and 234 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -24,63 +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_;
|
||||
|
||||
// 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();
|
||||
};
|
||||
|
||||
//===============================================================================
|
||||
//! Cached microscopic cross sections for a particular nuclide at the current
|
||||
//! energy
|
||||
|
|
@ -138,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(const int i_temp, const 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
|
||||
//==============================================================================
|
||||
|
|
@ -170,6 +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 int i_nuclide,
|
||||
const int i_temp, const double E);
|
||||
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -5,6 +5,7 @@
|
|||
|
||||
#include "xtensor/xtensor.hpp"
|
||||
|
||||
#include "openmc/constants.h"
|
||||
#include "openmc/hdf5_interface.h"
|
||||
|
||||
namespace openmc {
|
||||
|
|
@ -15,16 +16,14 @@ namespace openmc {
|
|||
|
||||
class UrrData{
|
||||
public:
|
||||
long unsigned int n_energy_; // # of incident energies
|
||||
long unsigned int n_prob_; // # of probabilities
|
||||
int interp_; // interpolation type (2=lin-lin, 5=log-log)
|
||||
Interpolation interp_; // interpolation type
|
||||
int inelastic_flag_; // inelastic competition flag
|
||||
int absorption_flag_; // other absorption flag
|
||||
bool multiply_smooth_; // multiply by smooth cross section?
|
||||
xt::xtensor<double, 1> energy_; // incident energies
|
||||
xt::xtensor<double, 3> prob_; // Actual probability tables
|
||||
|
||||
//! Load the URR data from the provided HDF5 group
|
||||
//! \brief Load the URR data from the provided HDF5 group
|
||||
void
|
||||
from_hdf5(hid_t group_id);
|
||||
};
|
||||
|
|
|
|||
153
src/nuclide.cpp
153
src/nuclide.cpp
|
|
@ -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");
|
||||
|
|
@ -365,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;
|
||||
|
|
@ -402,6 +404,142 @@ 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)
|
||||
{
|
||||
auto& micro = simulation::micro_xs[i_nuclide_];
|
||||
micro.use_ptable = true;
|
||||
|
||||
// Create a shorthand for the URR data
|
||||
UrrData* 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++;
|
||||
}
|
||||
|
||||
// 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
|
||||
// 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 (true) {
|
||||
if (urr->prob_(i_energy, URR_CUM_PROB, i_low) > r) {break;}
|
||||
i_low++;
|
||||
}
|
||||
int i_up = 0;
|
||||
while (true) {
|
||||
if (urr->prob_(i_energy + 1, URR_CUM_PROB, i_up) > r) {break;}
|
||||
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
|
||||
//==============================================================================
|
||||
|
|
@ -415,7 +553,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();
|
||||
}
|
||||
|
||||
|
|
@ -436,4 +575,10 @@ void set_micro_xs()
|
|||
}
|
||||
}
|
||||
|
||||
extern "C" void
|
||||
nuclide_calculate_urr_xs(const int i_nuclide, const int i_temp, const double E)
|
||||
{
|
||||
data::nuclides[i_nuclide - 1]->calculate_urr_xs(i_temp - 1, E);
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -198,6 +198,14 @@ module nuclide_header
|
|||
character(kind=C_CHAR), intent(in) :: name(*)
|
||||
type(C_PTR) :: path
|
||||
end function
|
||||
|
||||
subroutine nuclide_calculate_urr_xs_c(i_nuclide, i_temp, E) &
|
||||
bind(C, name='nuclide_calculate_urr_xs')
|
||||
import C_INT, C_DOUBLE
|
||||
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 interface
|
||||
|
||||
contains
|
||||
|
|
@ -980,7 +988,7 @@ contains
|
|||
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)
|
||||
call nuclide_calculate_urr_xs_c(this % i_nuclide, i_temp, E)
|
||||
end if
|
||||
end if
|
||||
|
||||
|
|
@ -1246,152 +1254,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
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
43
src/urr.cpp
43
src/urr.cpp
|
|
@ -1,25 +1,45 @@
|
|||
#include "openmc/urr.h"
|
||||
|
||||
#include <iostream>
|
||||
|
||||
namespace openmc {
|
||||
|
||||
void
|
||||
UrrData::from_hdf5(hid_t group_id)
|
||||
{
|
||||
// Read interpolation and other flags
|
||||
read_attribute(group_id, "interpolation", interp_);
|
||||
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;
|
||||
}
|
||||
|
||||
read_attribute(group_id, "inelastic", inelastic_flag_);
|
||||
read_attribute(group_id, "absorption", absorption_flag_);
|
||||
int i;
|
||||
read_attribute(group_id, "multiply_smooth", i);
|
||||
multiply_smooth_ = (i == 1);
|
||||
int temp_multiply_smooth;
|
||||
read_attribute(group_id, "multiply_smooth", temp_multiply_smooth);
|
||||
multiply_smooth_ = (temp_multiply_smooth == 1);
|
||||
|
||||
// read the enrgies at which tables exist
|
||||
hid_t dset = open_dataset(group_id, "energy");
|
||||
hsize_t dims[1];
|
||||
get_shape(dset, dims);
|
||||
close_dataset(dset);
|
||||
n_energy_ = static_cast<int>(dims[0]);
|
||||
energy_ = xt::xtensor<double, 1>({n_energy_}, 0.);
|
||||
// n_energy_ = static_cast<int>(dims[0]);
|
||||
energy_ = xt::xtensor<double, 1>({dims[0]}, 0.);
|
||||
read_dataset_as_shape(group_id, "energy", energy_);
|
||||
|
||||
// Read URR tables
|
||||
|
|
@ -27,14 +47,9 @@ UrrData::from_hdf5(hid_t group_id)
|
|||
hsize_t dims3[3];
|
||||
get_shape(dset, dims3);
|
||||
close_dataset(dset);
|
||||
n_prob_ = static_cast<int>(dims3[0]);
|
||||
xt::xarray<double> temp({n_energy_, 6, n_prob_});
|
||||
read_dataset(group_id, "table", temp);
|
||||
|
||||
prob_ = xt::xtensor<double, 3>({n_energy_, 6, n_prob_}, 0.);
|
||||
prob_ = temp;
|
||||
//TODO: swap 1st and last indices?
|
||||
|
||||
xt::xarray<double> temp_arr({dims3[0], dims3[1], dims3[2]}, 0.);
|
||||
read_dataset(group_id, "table", temp_arr);
|
||||
prob_ = temp_arr;
|
||||
}
|
||||
|
||||
}
|
||||
Loading…
Add table
Add a link
Reference in a new issue