mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-29 06:35:48 -04:00
Merge pull request #1047 from paulromano/cpp-thermal
Move thermal scattering data to C++
This commit is contained in:
Paul Romano
GitHub
commit
cf299aac73
24 changed files with 1153 additions and 1066 deletions
|
|
@ -288,7 +288,6 @@ set_target_properties(faddeeva PROPERTIES
|
|||
|
||||
add_library(libopenmc SHARED
|
||||
src/algorithm.F90
|
||||
src/angleenergy_header.F90
|
||||
src/bank_header.F90
|
||||
src/api.F90
|
||||
src/cmfd_data.F90
|
||||
|
|
@ -338,7 +337,6 @@ add_library(libopenmc SHARED
|
|||
src/reaction_header.F90
|
||||
src/relaxng
|
||||
src/sab_header.F90
|
||||
src/secondary_correlated.F90
|
||||
src/set_header.F90
|
||||
src/settings.F90
|
||||
src/simulation_header.F90
|
||||
|
|
@ -419,6 +417,7 @@ add_library(libopenmc SHARED
|
|||
src/state_point.cpp
|
||||
src/string_functions.cpp
|
||||
src/surface.cpp
|
||||
src/thermal.cpp
|
||||
src/xml_interface.cpp
|
||||
src/xsdata.cpp)
|
||||
set_target_properties(libopenmc PROPERTIES
|
||||
|
|
|
|||
|
|
@ -1,6 +1,7 @@
|
|||
include CMakeLists.txt
|
||||
include LICENSE
|
||||
include schemas.xml
|
||||
include pyproject.toml
|
||||
include openmc/data/reconstruct.pyx
|
||||
include docs/source/_templates/layout.html
|
||||
include docs/sphinxext/LICENSE
|
||||
|
|
|
|||
2
pyproject.toml
Normal file
2
pyproject.toml
Normal file
|
|
@ -0,0 +1,2 @@
|
|||
[build-system]
|
||||
requires = ["setuptools", "wheel", "numpy", "cython"]
|
||||
|
|
@ -1,38 +0,0 @@
|
|||
module angleenergy_header
|
||||
|
||||
use hdf5_interface, only: HID_T
|
||||
|
||||
!===============================================================================
|
||||
! ANGLEENERGY (abstract) defines a correlated or uncorrelated angle-energy
|
||||
! distribution that is a function of incoming energy. Each derived type must
|
||||
! implement a sample() subroutine that returns an outgoing energy and scattering
|
||||
! cosine given an incoming energy.
|
||||
!===============================================================================
|
||||
|
||||
type, abstract :: AngleEnergy
|
||||
contains
|
||||
procedure(angleenergy_sample_), deferred :: sample
|
||||
procedure(angleenergy_from_hdf5_), deferred :: from_hdf5
|
||||
end type AngleEnergy
|
||||
|
||||
abstract interface
|
||||
subroutine angleenergy_sample_(this, E_in, E_out, mu)
|
||||
import AngleEnergy
|
||||
class(AngleEnergy), intent(in) :: this
|
||||
real(8), intent(in) :: E_in
|
||||
real(8), intent(out) :: E_out
|
||||
real(8), intent(out) :: mu
|
||||
end subroutine angleenergy_sample_
|
||||
|
||||
subroutine angleenergy_from_hdf5_(this, group_id)
|
||||
import AngleEnergy, HID_T
|
||||
class(AngleEnergy), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
end subroutine angleenergy_from_hdf5_
|
||||
end interface
|
||||
|
||||
type :: AngleEnergyContainer
|
||||
class(AngleEnergy), allocatable :: obj
|
||||
end type AngleEnergyContainer
|
||||
|
||||
end module angleenergy_header
|
||||
|
|
@ -130,17 +130,6 @@ constexpr int ANGLE_HISTOGRAM {5};
|
|||
constexpr int TEMPERATURE_NEAREST {1};
|
||||
constexpr int TEMPERATURE_INTERPOLATION {2};
|
||||
|
||||
// Secondary energy mode for S(a,b) inelastic scattering
|
||||
// TODO: Convert to enum
|
||||
constexpr int SAB_SECONDARY_EQUAL {0}; // Equally-likely outgoing energy bins
|
||||
constexpr int SAB_SECONDARY_SKEWED {1}; // Skewed outgoing energy bins
|
||||
constexpr int SAB_SECONDARY_CONT {2}; // Continuous, linear-linear interpolation
|
||||
|
||||
// Elastic mode for S(a,b) elastic scattering
|
||||
// TODO: Convert to enum
|
||||
constexpr int SAB_ELASTIC_DISCRETE {3}; // Sample from discrete cosines
|
||||
constexpr int SAB_ELASTIC_EXACT {4}; // Exact treatment for coherent elastic
|
||||
|
||||
// Reaction types
|
||||
// TODO: Convert to enum
|
||||
constexpr int TOTAL_XS {1};
|
||||
|
|
@ -260,6 +249,8 @@ constexpr int N_AC {849};
|
|||
constexpr int N_2N0 {875};
|
||||
constexpr int N_2NC {891};
|
||||
|
||||
constexpr std::array<int, 6> DEPLETION_RX {N_GAMMA, N_P, N_A, N_2N, N_3N, N_4N};
|
||||
|
||||
// Fission neutron emission (nu) type
|
||||
constexpr int NU_NONE {0}; // No nu values (non-fissionable)
|
||||
constexpr int NU_POLYNOMIAL {1}; // Nu values given by polynomial
|
||||
|
|
|
|||
|
|
@ -107,6 +107,10 @@ public:
|
|||
//! Sample a value from the distribution
|
||||
//! \return Sampled value
|
||||
double sample() const;
|
||||
|
||||
// x property
|
||||
std::vector<double>& x() { return x_; }
|
||||
const std::vector<double>& x() const { return x_; }
|
||||
private:
|
||||
std::vector<double> x_; //!< tabulated independent variable
|
||||
std::vector<double> p_; //!< tabulated probability density
|
||||
|
|
|
|||
31
src/endf.cpp
31
src/endf.cpp
|
|
@ -144,4 +144,35 @@ double Tabulated1D::operator()(double x) const
|
|||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// CoherentElasticXS implementation
|
||||
//==============================================================================
|
||||
|
||||
CoherentElasticXS::CoherentElasticXS(hid_t dset)
|
||||
{
|
||||
// Read 2D array from dataset
|
||||
xt::xarray<double> arr;
|
||||
read_dataset(dset, arr);
|
||||
|
||||
// Get views for Bragg edges and structure factors
|
||||
auto E = xt::view(arr, 0);
|
||||
auto s = xt::view(arr, 1);
|
||||
|
||||
// Copy Bragg edges and partial sums of structure factors
|
||||
std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_));
|
||||
std::copy(s.begin(), s.end(), std::back_inserter(factors_));
|
||||
}
|
||||
|
||||
double CoherentElasticXS::operator()(double E) const
|
||||
{
|
||||
if (E < bragg_edges_[0]) {
|
||||
// If energy is below that of the lowest Bragg peak, the elastic cross
|
||||
// section will be zero
|
||||
return 0.0;
|
||||
} else {
|
||||
auto i_grid = lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E);
|
||||
return factors_[i_grid] / E;
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
13
src/endf.h
13
src/endf.h
|
|
@ -73,6 +73,19 @@ private:
|
|||
std::vector<double> y_; //!< values of ordinate
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
//! Coherent elastic scattering data from a crystalline material
|
||||
//==============================================================================
|
||||
|
||||
class CoherentElasticXS : public Function1D {
|
||||
explicit CoherentElasticXS(hid_t dset);
|
||||
double operator()(double E) const;
|
||||
private:
|
||||
std::vector<double> bragg_edges_; //!< Bragg edges in [eV]
|
||||
std::vector<double> factors_; //!< Partial sums of structure factors [eV-b]
|
||||
};
|
||||
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_ENDF_H
|
||||
|
|
|
|||
|
|
@ -327,7 +327,7 @@ get_groups(hid_t group_id, char* name[])
|
|||
}
|
||||
|
||||
std::vector<std::string>
|
||||
group_names(hid_t group_id)
|
||||
member_names(hid_t group_id, H5O_type_t type)
|
||||
{
|
||||
// Determine number of links in the group
|
||||
H5G_info_t info;
|
||||
|
|
@ -341,7 +341,7 @@ group_names(hid_t group_id)
|
|||
// Determine type of object (and skip non-group)
|
||||
H5Oget_info_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC, i, &oinfo,
|
||||
H5P_DEFAULT);
|
||||
if (oinfo.type != H5O_TYPE_GROUP) continue;
|
||||
if (oinfo.type != type) continue;
|
||||
|
||||
// Get size of name
|
||||
size = 1 + H5Lget_name_by_idx(group_id, ".", H5_INDEX_NAME, H5_ITER_INC,
|
||||
|
|
@ -356,6 +356,17 @@ group_names(hid_t group_id)
|
|||
return names;
|
||||
}
|
||||
|
||||
std::vector<std::string>
|
||||
group_names(hid_t group_id)
|
||||
{
|
||||
return member_names(group_id, H5O_TYPE_GROUP);
|
||||
}
|
||||
|
||||
std::vector<std::string>
|
||||
dataset_names(hid_t group_id)
|
||||
{
|
||||
return member_names(group_id, H5O_TYPE_DATASET);
|
||||
}
|
||||
|
||||
bool
|
||||
object_exists(hid_t object_id, const char* name)
|
||||
|
|
|
|||
|
|
@ -74,6 +74,7 @@ read_nd_vector(hid_t obj_id, const char* name,
|
|||
bool must_have = false);
|
||||
|
||||
std::vector<hsize_t> attribute_shape(hid_t obj_id, const char* name);
|
||||
std::vector<std::string> dataset_names(hid_t group_id);
|
||||
void ensure_exists(hid_t group_id, const char* name);
|
||||
std::vector<std::string> group_names(hid_t group_id);
|
||||
std::vector<hsize_t> object_shape(hid_t obj_id);
|
||||
|
|
@ -209,12 +210,31 @@ read_attribute(hid_t obj_id, const char* name, std::string& str)
|
|||
str = std::string{buffer, n};
|
||||
}
|
||||
|
||||
// overload for std::vector<std::string>
|
||||
inline void
|
||||
read_attribute(hid_t obj_id, const char* name, std::vector<std::string>& vec)
|
||||
{
|
||||
auto dims = attribute_shape(obj_id, name);
|
||||
auto m = dims[0];
|
||||
|
||||
// Allocate a C char array to get strings
|
||||
auto n = attribute_typesize(obj_id, name);
|
||||
char buffer[m][n+1];
|
||||
|
||||
// Read char data in attribute
|
||||
read_attr_string(obj_id, name, n, buffer[0]);
|
||||
|
||||
for (int i = 0; i < m; ++i) {
|
||||
vec.emplace_back(&buffer[i][0]);
|
||||
}
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Templates/overloads for read_dataset
|
||||
//==============================================================================
|
||||
|
||||
template<typename T>
|
||||
void read_dataset(hid_t obj_id, const char* name, T buffer, bool indep=false)
|
||||
void read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false)
|
||||
{
|
||||
read_dataset(obj_id, name, H5TypeMap<T>::type_id, &buffer, indep);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -223,7 +223,7 @@ contains
|
|||
found = .false.
|
||||
associate (sab => sab_tables(this % i_sab_tables(k)))
|
||||
FIND_NUCLIDE: do j = 1, size(this % nuclide)
|
||||
if (any(sab % nuclides == nuclides(this % nuclide(j)) % name)) then
|
||||
if (sab % has_nuclide(nuclides(this % nuclide(j)) % name)) then
|
||||
call i_sab_tables % push_back(this % i_sab_tables(k))
|
||||
call i_sab_nuclides % push_back(j)
|
||||
call sab_fracs % push_back(this % sab_fracs(k))
|
||||
|
|
@ -377,7 +377,7 @@ contains
|
|||
|
||||
! If particle energy is greater than the highest energy for the
|
||||
! S(a,b) table, then don't use the S(a,b) table
|
||||
if (p % E > sab_tables(i_sab) % data(1) % threshold_inelastic) then
|
||||
if (p % E > sab_tables(i_sab) % threshold()) then
|
||||
i_sab = 0
|
||||
end if
|
||||
|
||||
|
|
|
|||
67
src/nuclide.h
Normal file
67
src/nuclide.h
Normal file
|
|
@ -0,0 +1,67 @@
|
|||
#ifndef OPENMC_NUCLIDE_H
|
||||
#define OPENMC_NUCLIDE_H
|
||||
|
||||
#include "constants.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//===============================================================================
|
||||
//! Cached microscopic cross sections for a particular nuclide at the current
|
||||
//! energy
|
||||
//===============================================================================
|
||||
|
||||
struct NuclideMicroXS {
|
||||
// Microscopic cross sections in barns
|
||||
double total; //!< total cross section
|
||||
double absorption; //!< absorption (disappearance)
|
||||
double fission; //!< fission
|
||||
double nu_fission; //!< neutron production from fission
|
||||
|
||||
double elastic; //!< If sab_frac is not 1 or 0, then this value is
|
||||
//!< averaged over bound and non-bound nuclei
|
||||
double thermal; //!< Bound thermal elastic & inelastic scattering
|
||||
double thermal_elastic; //!< Bound thermal elastic scattering
|
||||
double photon_prod; //!< microscopic photon production xs
|
||||
|
||||
// Cross sections for depletion reactions (note that these are not stored in
|
||||
// macroscopic cache)
|
||||
double reaction[DEPLETION_RX.size()];
|
||||
|
||||
// Indicies and factors needed to compute cross sections from the data tables
|
||||
int index_grid; //!< Index on nuclide energy grid
|
||||
int index_temp; //!< Temperature index for nuclide
|
||||
double interp_factor; //!< Interpolation factor on nuc. energy grid
|
||||
int index_sab {-1}; //!< Index in sab_tables
|
||||
int index_temp_sab; //!< Temperature index for sab_tables
|
||||
double sab_frac; //!< Fraction of atoms affected by S(a,b)
|
||||
bool use_ptable; //!< In URR range with probability tables?
|
||||
|
||||
// Energy and temperature last used to evaluate these cross sections. If
|
||||
// these values have changed, then the cross sections must be re-evaluated.
|
||||
double last_E {0.0}; //!< Last evaluated energy
|
||||
double last_sqrtkT {0.0}; //!< Last temperature in sqrt(Boltzmann constant
|
||||
//!< * temperature (eV))
|
||||
};
|
||||
|
||||
//===============================================================================
|
||||
// MATERIALMACROXS contains cached macroscopic cross sections for the material a
|
||||
// particle is traveling through
|
||||
//===============================================================================
|
||||
|
||||
struct MaterialMacroXS {
|
||||
double total; //!< macroscopic total xs
|
||||
double absorption; //!< macroscopic absorption xs
|
||||
double fission; //!< macroscopic fission xs
|
||||
double nu_fission; //!< macroscopic production xs
|
||||
double photon_prod; //!< macroscopic photon production xs
|
||||
|
||||
// Photon cross sections
|
||||
double coherent; //!< macroscopic coherent xs
|
||||
double incoherent; //!< macroscopic incoherent xs
|
||||
double photoelectric; //!< macroscopic photoelectric xs
|
||||
double pair_production; //!< macroscopic pair production xs
|
||||
};
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_NUCLIDE_H
|
||||
|
|
@ -126,36 +126,36 @@ module nuclide_header
|
|||
! (NuclideMicroXS % elastic)
|
||||
real(8), parameter :: CACHE_INVALID = dble(Z"FFE0000000000000")
|
||||
|
||||
type NuclideMicroXS
|
||||
type, bind(C) :: NuclideMicroXS
|
||||
! Microscopic cross sections in barns
|
||||
real(8) :: total
|
||||
real(8) :: absorption ! absorption (disappearance)
|
||||
real(8) :: fission ! fission
|
||||
real(8) :: nu_fission ! neutron production from fission
|
||||
real(C_DOUBLE) :: total
|
||||
real(C_DOUBLE) :: absorption ! absorption (disappearance)
|
||||
real(C_DOUBLE) :: fission ! fission
|
||||
real(C_DOUBLE) :: nu_fission ! neutron production from fission
|
||||
|
||||
real(8) :: elastic ! If sab_frac is not 1 or 0, then this value is
|
||||
real(C_DOUBLE) :: elastic ! If sab_frac is not 1 or 0, then this value is
|
||||
! averaged over bound and non-bound nuclei
|
||||
real(8) :: thermal ! Bound thermal elastic & inelastic scattering
|
||||
real(8) :: thermal_elastic ! Bound thermal elastic scattering
|
||||
real(8) :: photon_prod ! microscopic photon production xs
|
||||
real(C_DOUBLE) :: thermal ! Bound thermal elastic & inelastic scattering
|
||||
real(C_DOUBLE) :: thermal_elastic ! Bound thermal elastic scattering
|
||||
real(C_DOUBLE) :: photon_prod ! microscopic photon production xs
|
||||
|
||||
! Cross sections for depletion reactions (note that these are not stored in
|
||||
! macroscopic cache)
|
||||
real(8) :: reaction(size(DEPLETION_RX))
|
||||
real(C_DOUBLE) :: reaction(size(DEPLETION_RX))
|
||||
|
||||
! Indicies and factors needed to compute cross sections from the data tables
|
||||
integer :: index_grid ! Index on nuclide energy grid
|
||||
integer :: index_temp ! Temperature index for nuclide
|
||||
real(8) :: interp_factor ! Interpolation factor on nuc. energy grid
|
||||
integer :: index_sab = NONE ! Index in sab_tables
|
||||
integer :: index_temp_sab ! Temperature index for sab_tables
|
||||
real(8) :: sab_frac ! Fraction of atoms affected by S(a,b)
|
||||
logical :: use_ptable ! In URR range with probability tables?
|
||||
integer(C_INT) :: index_grid ! Index on nuclide energy grid
|
||||
integer(C_INT) :: index_temp ! Temperature index for nuclide
|
||||
real(C_DOUBLE) :: interp_factor ! Interpolation factor on nuc. energy grid
|
||||
integer(C_INT) :: index_sab = NONE ! Index in sab_tables
|
||||
integer(C_INT) :: index_temp_sab ! Temperature index for sab_tables
|
||||
real(C_DOUBLE) :: sab_frac ! Fraction of atoms affected by S(a,b)
|
||||
logical(C_BOOL) :: use_ptable ! In URR range with probability tables?
|
||||
|
||||
! Energy and temperature last used to evaluate these cross sections. If
|
||||
! these values have changed, then the cross sections must be re-evaluated.
|
||||
real(8) :: last_E = ZERO ! Last evaluated energy
|
||||
real(8) :: last_sqrtkT = ZERO ! Last temperature in sqrt(Boltzmann
|
||||
real(C_DOUBLE) :: last_E = ZERO ! Last evaluated energy
|
||||
real(C_DOUBLE) :: last_sqrtkT = ZERO ! Last temperature in sqrt(Boltzmann
|
||||
! constant * temperature (eV))
|
||||
end type NuclideMicroXS
|
||||
|
||||
|
|
@ -164,7 +164,7 @@ module nuclide_header
|
|||
! particle is traveling through
|
||||
!===============================================================================
|
||||
|
||||
type MaterialMacroXS
|
||||
type, bind(C) :: MaterialMacroXS
|
||||
real(C_DOUBLE) :: total ! macroscopic total xs
|
||||
real(C_DOUBLE) :: absorption ! macroscopic absorption xs
|
||||
real(C_DOUBLE) :: fission ! macroscopic fission xs
|
||||
|
|
@ -198,7 +198,7 @@ module nuclide_header
|
|||
type(DictCharInt) :: nuclide_dict
|
||||
|
||||
! Cross section caches
|
||||
type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide
|
||||
type(NuclideMicroXS), allocatable, target :: micro_xs(:) ! Cache for each nuclide
|
||||
type(MaterialMacroXS) :: material_xs ! Cache for current material
|
||||
!$omp threadprivate(micro_xs, material_xs)
|
||||
|
||||
|
|
@ -1096,9 +1096,9 @@ contains
|
|||
real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b)
|
||||
type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache
|
||||
|
||||
integer :: i_temp ! temperature index
|
||||
real(8) :: inelastic ! S(a,b) inelastic cross section
|
||||
real(8) :: elastic ! S(a,b) elastic cross section
|
||||
integer(C_INT) :: i_temp ! temperature index
|
||||
real(C_DOUBLE) :: inelastic ! S(a,b) inelastic cross section
|
||||
real(C_DOUBLE) :: elastic ! S(a,b) elastic cross section
|
||||
|
||||
! Set flag that S(a,b) treatment should be used for scattering
|
||||
micro_xs % index_sab = i_sab
|
||||
|
|
|
|||
261
src/physics.F90
261
src/physics.F90
|
|
@ -884,263 +884,16 @@ contains
|
|||
real(8), intent(inout) :: uvw(3) ! directional cosines
|
||||
real(8), intent(out) :: mu ! scattering cosine
|
||||
|
||||
integer :: i ! incoming energy bin
|
||||
integer :: j ! outgoing energy bin
|
||||
integer :: k ! outgoing cosine bin
|
||||
integer :: i_temp ! temperature index
|
||||
integer :: n_energy_out ! number of outgoing energy bins
|
||||
real(8) :: f ! interpolation factor
|
||||
real(8) :: r ! used for skewed sampling & continuous
|
||||
real(8) :: E_ij ! outgoing energy j for E_in(i)
|
||||
real(8) :: E_i1j ! outgoing energy j for E_in(i+1)
|
||||
real(8) :: mu_ijk ! outgoing cosine k for E_in(i) and E_out(j)
|
||||
real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j)
|
||||
real(8) :: prob ! probability for sampling Bragg edge
|
||||
! Following are needed only for SAB_SECONDARY_CONT scattering
|
||||
integer :: l ! sampled incoming E bin (is i or i + 1)
|
||||
real(8) :: E_i_1, E_i_J ! endpoints on outgoing grid i
|
||||
real(8) :: E_i1_1, E_i1_J ! endpoints on outgoing grid i+1
|
||||
real(8) :: E_1, E_J ! endpoints interpolated between i and i+1
|
||||
real(8) :: E_l_j, E_l_j1 ! adjacent E on outgoing grid l
|
||||
real(8) :: p_l_j, p_l_j1 ! adjacent p on outgoing grid l
|
||||
real(8) :: c_j, c_j1 ! cumulative probability
|
||||
real(8) :: frac ! interpolation factor on outgoing energy
|
||||
real(8) :: r1 ! RNG for outgoing energy
|
||||
real(8) :: mu_left, mu_right ! adjacent mu values
|
||||
real(C_DOUBLE) :: E_out
|
||||
type(C_PTR) :: ptr
|
||||
|
||||
i_temp = micro_xs(i_nuclide) % index_temp_sab
|
||||
! Sample from C++ side
|
||||
ptr = C_LOC(micro_xs(i_nuclide))
|
||||
call sab_tables(i_sab) % sample(ptr, E, E_out, mu)
|
||||
|
||||
! Get pointer to S(a,b) table
|
||||
associate (sab => sab_tables(i_sab) % data(i_temp))
|
||||
|
||||
! Determine whether inelastic or elastic scattering will occur
|
||||
if (prn() < micro_xs(i_nuclide) % thermal_elastic / &
|
||||
micro_xs(i_nuclide) % thermal) then
|
||||
! elastic scattering
|
||||
|
||||
! Get index and interpolation factor for elastic grid
|
||||
if (E < sab % elastic_e_in(1)) then
|
||||
i = 1
|
||||
f = ZERO
|
||||
else
|
||||
i = binary_search(sab % elastic_e_in, sab % n_elastic_e_in, E)
|
||||
f = (E - sab%elastic_e_in(i)) / &
|
||||
(sab%elastic_e_in(i+1) - sab%elastic_e_in(i))
|
||||
end if
|
||||
|
||||
! Select treatment based on elastic mode
|
||||
if (sab % elastic_mode == SAB_ELASTIC_DISCRETE) then
|
||||
! 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
|
||||
k = 1 + int(prn() * sab % n_elastic_mu)
|
||||
|
||||
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
|
||||
mu_ijk = sab % elastic_mu(k,i)
|
||||
mu_i1jk = sab % elastic_mu(k,i+1)
|
||||
|
||||
! Cosine of angle between incoming and outgoing neutron
|
||||
mu = (1 - f)*mu_ijk + f*mu_i1jk
|
||||
|
||||
elseif (sab % elastic_mode == SAB_ELASTIC_EXACT) then
|
||||
! 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
|
||||
prob = prn() * sab % elastic_P(i+1)
|
||||
if (prob < sab % elastic_P(1)) then
|
||||
k = 1
|
||||
else
|
||||
k = binary_search(sab % elastic_P(1:i+1), i+1, prob)
|
||||
end if
|
||||
|
||||
! Characteristic scattering cosine for this Bragg edge
|
||||
mu = ONE - TWO*sab % elastic_e_in(k) / E
|
||||
|
||||
end if
|
||||
|
||||
! Outgoing energy is same as incoming energy -- no need to do anything
|
||||
|
||||
else
|
||||
! Perform inelastic calculations
|
||||
|
||||
! Get index and interpolation factor for inelastic grid
|
||||
if (E < sab % inelastic_e_in(1)) then
|
||||
i = 1
|
||||
f = ZERO
|
||||
else
|
||||
i = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E)
|
||||
f = (E - sab%inelastic_e_in(i)) / &
|
||||
(sab%inelastic_e_in(i+1) - sab%inelastic_e_in(i))
|
||||
end if
|
||||
|
||||
! 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 ((sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_EQUAL) .or. &
|
||||
(sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_SKEWED)) then
|
||||
if (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_EQUAL) then
|
||||
! All bins equally likely
|
||||
|
||||
j = 1 + int(prn() * sab % n_inelastic_e_out)
|
||||
elseif (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_SKEWED) then
|
||||
! Distribution skewed away from edge points
|
||||
|
||||
! Determine number of outgoing energy and angle bins
|
||||
n_energy_out = sab % n_inelastic_e_out
|
||||
|
||||
r = prn() * (n_energy_out - 3)
|
||||
if (r > ONE) then
|
||||
! equally likely N-4 middle bins
|
||||
j = int(r) + 2
|
||||
elseif (r > 0.6_8) then
|
||||
! second to last bin has relative probability of 0.4
|
||||
j = n_energy_out - 1
|
||||
elseif (r > HALF) then
|
||||
! last bin has relative probability of 0.1
|
||||
j = n_energy_out
|
||||
elseif (r > 0.1_8) then
|
||||
! second bin has relative probability of 0.4
|
||||
j = 2
|
||||
else
|
||||
! first bin has relative probability of 0.1
|
||||
j = 1
|
||||
end if
|
||||
end if
|
||||
|
||||
! Determine outgoing energy corresponding to E_in(i) and E_in(i+1)
|
||||
E_ij = sab % inelastic_e_out(j,i)
|
||||
E_i1j = sab % inelastic_e_out(j,i+1)
|
||||
|
||||
! Outgoing energy
|
||||
E = (1 - f)*E_ij + f*E_i1j
|
||||
|
||||
! Sample outgoing cosine bin
|
||||
k = 1 + int(prn() * sab % n_inelastic_mu)
|
||||
|
||||
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
|
||||
mu_ijk = sab % inelastic_mu(k,j,i)
|
||||
mu_i1jk = sab % inelastic_mu(k,j,i+1)
|
||||
|
||||
! Cosine of angle between incoming and outgoing neutron
|
||||
mu = (1 - f)*mu_ijk + f*mu_i1jk
|
||||
|
||||
else if (sab_tables(i_sab) % secondary_mode == SAB_SECONDARY_CONT) then
|
||||
! Continuous secondary energy - this is to be similar to
|
||||
! Law 61 interpolation on outgoing energy
|
||||
|
||||
! Sample between ith and (i+1)th bin
|
||||
r = prn()
|
||||
if (f > r) then
|
||||
l = i + 1
|
||||
else
|
||||
l = i
|
||||
end if
|
||||
|
||||
! Determine endpoints on grid i
|
||||
n_energy_out = sab % inelastic_data(i) % n_e_out
|
||||
E_i_1 = sab % inelastic_data(i) % e_out(1)
|
||||
E_i_J = sab % inelastic_data(i) % e_out(n_energy_out)
|
||||
|
||||
! Determine endpoints on grid i + 1
|
||||
n_energy_out = sab % inelastic_data(i + 1) % n_e_out
|
||||
E_i1_1 = sab % inelastic_data(i + 1) % e_out(1)
|
||||
E_i1_J = sab % inelastic_data(i + 1) % e_out(n_energy_out)
|
||||
|
||||
E_1 = E_i_1 + f * (E_i1_1 - E_i_1)
|
||||
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_energy_out = sab % inelastic_data(l) % n_e_out
|
||||
r1 = prn()
|
||||
c_j = sab % inelastic_data(l) % e_out_cdf(1)
|
||||
do j = 1, n_energy_out - 1
|
||||
c_j1 = sab % inelastic_data(l) % e_out_cdf(j + 1)
|
||||
if (r1 < c_j1) exit
|
||||
c_j = c_j1
|
||||
end do
|
||||
|
||||
! check to make sure k is <= n_energy_out - 1
|
||||
j = min(j, n_energy_out - 1)
|
||||
|
||||
! Get the data to interpolate between
|
||||
E_l_j = sab % inelastic_data(l) % e_out(j)
|
||||
p_l_j = sab % inelastic_data(l) % e_out_pdf(j)
|
||||
|
||||
! Next part assumes linear-linear interpolation in standard
|
||||
E_l_j1 = sab % inelastic_data(l) % e_out(j + 1)
|
||||
p_l_j1 = sab % inelastic_data(l) % e_out_pdf(j + 1)
|
||||
|
||||
! Find secondary energy (variable E)
|
||||
frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j)
|
||||
if (frac == ZERO) then
|
||||
E = E_l_j + (r1 - c_j) / p_l_j
|
||||
else
|
||||
E = E_l_j + (sqrt(max(ZERO, p_l_j * p_l_j + &
|
||||
TWO * frac * (r1 - c_j))) - p_l_j) / frac
|
||||
end if
|
||||
|
||||
! Now interpolate between incident energy bins i and i + 1
|
||||
if (l == i) then
|
||||
E = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1)
|
||||
else
|
||||
E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1)
|
||||
end if
|
||||
|
||||
! Sample outgoing cosine bin
|
||||
k = 1 + int(prn() * sab % 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.
|
||||
associate (mu_l => sab % inelastic_data(l) % mu)
|
||||
f = (r1 - c_j)/(c_j1 - c_j)
|
||||
|
||||
! Determine (k-1)th mu value
|
||||
if (k == 1) then
|
||||
mu_left = -ONE
|
||||
else
|
||||
mu_left = mu_l(k-1, j) + f*(mu_l(k-1, j+1) - mu_l(k-1,j))
|
||||
end if
|
||||
|
||||
! Determine kth mu value
|
||||
mu = mu_l(k, j) + f*(mu_l(k, j+1) - mu_l(k, j))
|
||||
|
||||
! Determine (k+1)th mu value
|
||||
if (k == sab % n_inelastic_mu) then
|
||||
mu_right = ONE - mu
|
||||
else
|
||||
mu_right = mu_l(k+1, j) + f*(mu_l(k+1, j+1) - mu_l(k+1,j)) - mu
|
||||
end if
|
||||
end associate
|
||||
|
||||
! Smear angle
|
||||
mu = mu + min(mu - mu_left, mu_right - mu)*(prn() - HALF)
|
||||
|
||||
end if ! (inelastic secondary energy treatment)
|
||||
end if ! (elastic or inelastic)
|
||||
end associate
|
||||
|
||||
! 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 (abs(mu) > ONE) mu = sign(ONE,mu)
|
||||
|
||||
! change direction of particle
|
||||
! Set energy to outgoing, change direction of particle
|
||||
E = E_out
|
||||
uvw = rotate_angle(uvw, mu)
|
||||
|
||||
end subroutine sab_scatter
|
||||
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -1,468 +1,172 @@
|
|||
module sab_header
|
||||
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
use, intrinsic :: ISO_FORTRAN_ENV
|
||||
|
||||
use algorithm, only: find, sort, binary_search
|
||||
use constants
|
||||
use dict_header, only: DictIntInt, DictCharInt
|
||||
use distribution_univariate, only: Tabular
|
||||
use error, only: warning, fatal_error
|
||||
use dict_header, only: DictCharInt
|
||||
use hdf5_interface
|
||||
use random_lcg, only: prn
|
||||
use secondary_correlated, only: CorrelatedAngleEnergy
|
||||
use settings
|
||||
use stl_vector, only: VectorInt, VectorReal
|
||||
use string, only: to_str, str_to_int
|
||||
use stl_vector, only: VectorReal
|
||||
use string, only: to_c_string
|
||||
|
||||
implicit none
|
||||
private
|
||||
|
||||
!===============================================================================
|
||||
! DISTENERGYSAB contains the secondary energy/angle distributions for inelastic
|
||||
! thermal scattering collisions which utilize a continuous secondary energy
|
||||
! representation.
|
||||
!===============================================================================
|
||||
|
||||
type DistEnergySab
|
||||
integer :: n_e_out
|
||||
real(8), allocatable :: e_out(:)
|
||||
real(8), allocatable :: e_out_pdf(:)
|
||||
real(8), allocatable :: e_out_cdf(:)
|
||||
real(8), allocatable :: mu(:,:)
|
||||
end type DistEnergySab
|
||||
public :: free_memory_sab
|
||||
|
||||
!===============================================================================
|
||||
! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off
|
||||
! of light isotopes such as water, graphite, Be, etc
|
||||
!===============================================================================
|
||||
|
||||
type SabData
|
||||
! threshold for S(a,b) treatment (usually ~4 eV)
|
||||
real(8) :: threshold_inelastic
|
||||
real(8) :: threshold_elastic = ZERO
|
||||
|
||||
! Inelastic scattering data
|
||||
integer :: n_inelastic_e_in ! # of incoming E for inelastic
|
||||
integer :: n_inelastic_e_out ! # of outgoing E for inelastic
|
||||
integer :: n_inelastic_mu ! # of outgoing angles for inelastic
|
||||
real(8), allocatable :: inelastic_e_in(:)
|
||||
real(8), allocatable :: inelastic_sigma(:)
|
||||
! The following are used only if secondary_mode is 0 or 1
|
||||
real(8), allocatable :: inelastic_e_out(:,:)
|
||||
real(8), allocatable :: inelastic_mu(:,:,:)
|
||||
! The following is used only if secondary_mode is 3
|
||||
! The different implementation is necessary because the continuous
|
||||
! representation has a variable number of outgoing energy points for each
|
||||
! incoming energy
|
||||
type(DistEnergySab), allocatable :: inelastic_data(:) ! One for each Ein
|
||||
|
||||
! Elastic scattering data
|
||||
integer :: elastic_mode ! elastic mode (discrete/exact)
|
||||
integer :: n_elastic_e_in ! # of incoming E for elastic
|
||||
integer :: n_elastic_mu ! # of outgoing angles for elastic
|
||||
real(8), allocatable :: elastic_e_in(:)
|
||||
real(8), allocatable :: elastic_P(:)
|
||||
real(8), allocatable :: elastic_mu(:,:)
|
||||
end type SabData
|
||||
|
||||
type SAlphaBeta
|
||||
character(150) :: name ! name of table, e.g. lwtr.10t
|
||||
real(8) :: awr ! weight of nucleus in neutron masses
|
||||
real(8), allocatable :: kTs(:) ! temperatures in eV (k*T)
|
||||
character(10), allocatable :: nuclides(:) ! List of valid nuclides
|
||||
integer :: secondary_mode ! secondary mode (equal/skewed/continuous)
|
||||
|
||||
! cross sections and distributions at each temperature
|
||||
type(SabData), allocatable :: data(:)
|
||||
type, public :: SAlphaBeta
|
||||
type(C_PTR) :: ptr
|
||||
contains
|
||||
procedure :: from_hdf5 => salphabeta_from_hdf5
|
||||
procedure :: calculate_xs => sab_calculate_xs
|
||||
procedure :: from_hdf5
|
||||
procedure :: calculate_xs
|
||||
procedure :: free
|
||||
procedure :: has_nuclide
|
||||
procedure :: sample
|
||||
procedure :: threshold
|
||||
end type SAlphaBeta
|
||||
|
||||
! S(a,b) tables
|
||||
type(SAlphaBeta), allocatable, target :: sab_tables(:)
|
||||
integer(C_INT), bind(C) :: n_sab_tables
|
||||
type(DictCharInt) :: sab_dict
|
||||
type(SAlphaBeta), public, allocatable, target :: sab_tables(:)
|
||||
integer(C_INT), public, bind(C) :: n_sab_tables
|
||||
type(DictCharInt), public :: sab_dict
|
||||
|
||||
interface
|
||||
function sab_from_hdf5(group_id, temperature, n, method, &
|
||||
tolerance, minmax) result(ptr) bind(C)
|
||||
import HID_T, C_DOUBLE, C_INT, C_PTR
|
||||
integer(HID_T), value :: group_id
|
||||
real(C_DOUBLE), intent(in) :: temperature
|
||||
integer(C_INT), value :: n
|
||||
integer(C_INT), value :: method
|
||||
real(C_DOUBLE), value :: tolerance
|
||||
real(C_DOUBLE), intent(in) :: minmax(2)
|
||||
type(C_PTR) :: ptr
|
||||
end function
|
||||
|
||||
subroutine sab_calculate_xs(ptr, E, sqrtkT, i_temp, elastic, &
|
||||
inelastic) bind(C)
|
||||
import C_PTR, C_DOUBLE, C_INT
|
||||
type(C_PTR), value :: ptr
|
||||
real(C_DOUBLE), value :: E
|
||||
real(C_DOUBLE), value :: sqrtkT
|
||||
integer(C_INT), intent(out) :: i_temp
|
||||
real(C_DOUBLE), intent(out) :: elastic
|
||||
real(C_DOUBLE), intent(out) :: inelastic
|
||||
end subroutine
|
||||
|
||||
subroutine sab_free(ptr) bind(C)
|
||||
import C_PTR
|
||||
type(C_PTR), value :: ptr
|
||||
end subroutine
|
||||
|
||||
function sab_has_nuclide(ptr, name) result(val) bind(C)
|
||||
import C_PTR, C_CHAR, C_BOOL
|
||||
type(C_PTR), value :: ptr
|
||||
character(kind=C_CHAR), intent(in) :: name(*)
|
||||
logical(C_BOOL) :: val
|
||||
end function
|
||||
|
||||
subroutine sab_sample(ptr, micro_xs, E_in, E_out, mu) bind(C)
|
||||
import C_PTR, C_INT, C_DOUBLE
|
||||
type(C_PTR), value :: ptr
|
||||
type(C_PTR), value :: micro_xs
|
||||
real(C_DOUBLE), value :: E_in
|
||||
real(C_DOUBLE), intent(out) :: E_out
|
||||
real(C_DOUBLE), intent(out) :: mu
|
||||
end subroutine
|
||||
|
||||
function sab_threshold(ptr) result(threshold) bind(C)
|
||||
import C_PTR, C_double
|
||||
type(C_PTR), value :: ptr
|
||||
real(C_DOUBLE) :: threshold
|
||||
end function
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
subroutine salphabeta_from_hdf5(this, group_id, temperature, method, &
|
||||
subroutine from_hdf5(this, group_id, temperature, method, &
|
||||
tolerance, minmax)
|
||||
class(SAlphaBeta), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
type(VectorReal), intent(in) :: temperature ! list of temperatures
|
||||
integer, intent(in) :: method
|
||||
real(8), intent(in) :: tolerance
|
||||
real(8), intent(in) :: minmax(2)
|
||||
integer(C_INT), intent(in) :: method
|
||||
real(C_DOUBLE), intent(in) :: tolerance
|
||||
real(C_DOUBLE), intent(in) :: minmax(2)
|
||||
|
||||
integer :: i, j
|
||||
integer :: t
|
||||
integer :: n_energy, n_energy_out, n_mu
|
||||
integer :: i_closest
|
||||
integer :: n_temperature
|
||||
integer(HID_T) :: T_group
|
||||
integer(HID_T) :: elastic_group
|
||||
integer(HID_T) :: inelastic_group
|
||||
integer(HID_T) :: dset_id
|
||||
integer(HID_T) :: kT_group
|
||||
integer(HSIZE_T) :: dims2(2)
|
||||
integer(HSIZE_T) :: dims3(3)
|
||||
real(8), allocatable :: temp(:,:)
|
||||
character(20) :: type
|
||||
type(CorrelatedAngleEnergy) :: correlated_dist
|
||||
real(C_DOUBLE) :: dummy
|
||||
integer(C_INT) :: n
|
||||
|
||||
character(MAX_WORD_LEN) :: temp_str
|
||||
character(MAX_WORD_LEN), allocatable :: dset_names(:)
|
||||
real(8), allocatable :: temps_available(:) ! temperatures available
|
||||
real(8) :: temp_desired
|
||||
real(8) :: temp_actual
|
||||
type(VectorInt) :: temps_to_read
|
||||
|
||||
! Get name of table from group
|
||||
this % name = get_name(group_id)
|
||||
|
||||
! Get rid of leading '/'
|
||||
this % name = trim(this % name(2:))
|
||||
|
||||
call read_attribute(this % awr, group_id, 'atomic_weight_ratio')
|
||||
call read_attribute(this % nuclides, group_id, 'nuclides')
|
||||
call read_attribute(type, group_id, 'secondary_mode')
|
||||
select case (type)
|
||||
case ('equal')
|
||||
this % secondary_mode = SAB_SECONDARY_EQUAL
|
||||
case ('skewed')
|
||||
this % secondary_mode = SAB_SECONDARY_SKEWED
|
||||
case ('continuous')
|
||||
this % secondary_mode = SAB_SECONDARY_CONT
|
||||
end select
|
||||
|
||||
! Read temperatures
|
||||
kT_group = open_group(group_id, 'kTs')
|
||||
|
||||
! Determine temperatures available
|
||||
call get_datasets(kT_group, dset_names)
|
||||
allocate(temps_available(size(dset_names)))
|
||||
do i = 1, size(dset_names)
|
||||
! Read temperature value
|
||||
call read_dataset(temps_available(i), kT_group, trim(dset_names(i)))
|
||||
temps_available(i) = temps_available(i) / K_BOLTZMANN
|
||||
end do
|
||||
call sort(temps_available)
|
||||
|
||||
! Determine actual temperatures to read -- start by checking whether a
|
||||
! temperature range was given, in which case all temperatures in the range
|
||||
! are loaded irrespective of what temperatures actually appear in the model
|
||||
if (minmax(2) > ZERO) then
|
||||
do i = 1, size(temps_available)
|
||||
temp_actual = temps_available(i)
|
||||
if (minmax(1) <= temp_actual .and. temp_actual <= minmax(2)) then
|
||||
call temps_to_read % push_back(nint(temp_actual))
|
||||
end if
|
||||
end do
|
||||
n = temperature % size()
|
||||
if (n > 0) then
|
||||
this % ptr = sab_from_hdf5(group_id, temperature % data(1), n, method, tolerance, minmax)
|
||||
else
|
||||
! In this case, temperatures % data(1) doesn't exist, so we just pass a
|
||||
! dummy value
|
||||
this % ptr = sab_from_hdf5(group_id, dummy, n, method, tolerance, minmax)
|
||||
end if
|
||||
|
||||
select case (method)
|
||||
case (TEMPERATURE_NEAREST)
|
||||
! Determine actual temperatures to read
|
||||
do i = 1, temperature % size()
|
||||
temp_desired = temperature % data(i)
|
||||
i_closest = minloc(abs(temps_available - temp_desired), dim=1)
|
||||
temp_actual = temps_available(i_closest)
|
||||
if (abs(temp_actual - temp_desired) < tolerance) then
|
||||
if (find(temps_to_read, nint(temp_actual)) == -1) then
|
||||
call temps_to_read % push_back(nint(temp_actual))
|
||||
end if
|
||||
else
|
||||
call fatal_error("Nuclear data library does not contain cross sections &
|
||||
&for " // trim(this % name) // " at or near " // &
|
||||
trim(to_str(nint(temp_desired))) // " K.")
|
||||
end if
|
||||
end do
|
||||
|
||||
case (TEMPERATURE_INTERPOLATION)
|
||||
! If temperature interpolation or multipole is selected, get a list of
|
||||
! bounding temperatures for each actual temperature present in the model
|
||||
TEMP_LOOP: do i = 1, temperature % size()
|
||||
temp_desired = temperature % data(i)
|
||||
|
||||
do j = 1, size(temps_available) - 1
|
||||
if (temps_available(j) <= temp_desired .and. &
|
||||
temp_desired < temps_available(j + 1)) then
|
||||
if (find(temps_to_read, nint(temps_available(j))) == -1) then
|
||||
call temps_to_read % push_back(nint(temps_available(j)))
|
||||
end if
|
||||
if (find(temps_to_read, nint(temps_available(j + 1))) == -1) then
|
||||
call temps_to_read % push_back(nint(temps_available(j + 1)))
|
||||
end if
|
||||
cycle TEMP_LOOP
|
||||
end if
|
||||
end do
|
||||
|
||||
call fatal_error("Nuclear data library does not contain cross sections &
|
||||
&for " // trim(this % name) // " at temperatures that bound " // &
|
||||
trim(to_str(nint(temp_desired))) // " K.")
|
||||
end do TEMP_LOOP
|
||||
|
||||
end select
|
||||
|
||||
! Sort temperatures to read
|
||||
call sort(temps_to_read)
|
||||
|
||||
n_temperature = temps_to_read % size()
|
||||
allocate(this % kTs(n_temperature))
|
||||
allocate(this % data(n_temperature))
|
||||
|
||||
do t = 1, n_temperature
|
||||
! Get temperature as a string
|
||||
temp_str = trim(to_str(temps_to_read % data(t))) // "K"
|
||||
|
||||
! Read exact temperature value
|
||||
call read_dataset(this % kTs(t), kT_group, temp_str)
|
||||
|
||||
! Open group for temperature i
|
||||
T_group = open_group(group_id, temp_str)
|
||||
|
||||
! Coherent elastic data
|
||||
if (object_exists(T_group, 'elastic')) then
|
||||
! Read cross section data
|
||||
elastic_group = open_group(T_group, 'elastic')
|
||||
dset_id = open_dataset(elastic_group, 'xs')
|
||||
call read_attribute(type, dset_id, 'type')
|
||||
call get_shape(dset_id, dims2)
|
||||
allocate(temp(dims2(1), dims2(2)))
|
||||
call read_dataset(temp, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
! Set cross section data and type
|
||||
this % data(t) % n_elastic_e_in = int(dims2(1), 4)
|
||||
allocate(this % data(t) % elastic_e_in(this % data(t) % n_elastic_e_in))
|
||||
allocate(this % data(t) % elastic_P(this % data(t) % n_elastic_e_in))
|
||||
this % data(t) % elastic_e_in(:) = temp(:, 1)
|
||||
this % data(t) % elastic_P(:) = temp(:, 2)
|
||||
select case (type)
|
||||
case ('tab1')
|
||||
this % data(t) % elastic_mode = SAB_ELASTIC_DISCRETE
|
||||
case ('bragg')
|
||||
this % data(t) % elastic_mode = SAB_ELASTIC_EXACT
|
||||
end select
|
||||
deallocate(temp)
|
||||
|
||||
! Set elastic threshold
|
||||
this % data(t) % threshold_elastic = this % data(t) % elastic_e_in(&
|
||||
this % data(t) % n_elastic_e_in)
|
||||
|
||||
! Read angle distribution
|
||||
if (this % data(t) % elastic_mode /= SAB_ELASTIC_EXACT) then
|
||||
dset_id = open_dataset(elastic_group, 'mu_out')
|
||||
call get_shape(dset_id, dims2)
|
||||
this % data(t) % n_elastic_mu = int(dims2(1), 4)
|
||||
allocate(this % data(t) % elastic_mu(dims2(1), dims2(2)))
|
||||
call read_dataset(this % data(t) % elastic_mu, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
end if
|
||||
|
||||
call close_group(elastic_group)
|
||||
end if
|
||||
|
||||
! Inelastic data
|
||||
if (object_exists(T_group, 'inelastic')) then
|
||||
! Read type of inelastic data
|
||||
inelastic_group = open_group(T_group, 'inelastic')
|
||||
|
||||
! Read cross section data
|
||||
dset_id = open_dataset(inelastic_group, 'xs')
|
||||
call get_shape(dset_id, dims2)
|
||||
allocate(temp(dims2(1), dims2(2)))
|
||||
call read_dataset(temp, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
! Set cross section data
|
||||
this % data(t) % n_inelastic_e_in = int(dims2(1), 4)
|
||||
allocate(this % data(t) % inelastic_e_in(this % data(t) % n_inelastic_e_in))
|
||||
allocate(this % data(t) % inelastic_sigma(this % data(t) % n_inelastic_e_in))
|
||||
this % data(t) % inelastic_e_in(:) = temp(:, 1)
|
||||
this % data(t) % inelastic_sigma(:) = temp(:, 2)
|
||||
deallocate(temp)
|
||||
|
||||
! Set inelastic threshold
|
||||
this % data(t) % threshold_inelastic = this % data(t) % inelastic_e_in(&
|
||||
this % data(t) % n_inelastic_e_in)
|
||||
|
||||
if (this % secondary_mode /= SAB_SECONDARY_CONT) then
|
||||
! Read energy distribution
|
||||
dset_id = open_dataset(inelastic_group, 'energy_out')
|
||||
call get_shape(dset_id, dims2)
|
||||
this % data(t) % n_inelastic_e_out = int(dims2(1), 4)
|
||||
allocate(this % data(t) % inelastic_e_out(dims2(1), dims2(2)))
|
||||
call read_dataset(this % data(t) % inelastic_e_out, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
! Read angle distribution
|
||||
dset_id = open_dataset(inelastic_group, 'mu_out')
|
||||
call get_shape(dset_id, dims3)
|
||||
this % data(t) % n_inelastic_mu = int(dims3(1), 4)
|
||||
allocate(this % data(t) % inelastic_mu(dims3(1), dims3(2), dims3(3)))
|
||||
call read_dataset(this % data(t) % inelastic_mu, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
else
|
||||
! Read correlated angle-energy distribution
|
||||
call correlated_dist % from_hdf5(inelastic_group)
|
||||
|
||||
! Convert to S(a,b) native format
|
||||
n_energy = size(correlated_dist % energy)
|
||||
allocate(this % data(t) % inelastic_data(n_energy))
|
||||
do i = 1, n_energy
|
||||
associate (edist => correlated_dist % distribution(i))
|
||||
! Get number of outgoing energies for incoming energy i
|
||||
n_energy_out = size(edist % e_out)
|
||||
this % data(t) % inelastic_data(i) % n_e_out = n_energy_out
|
||||
allocate(this % data(t) % inelastic_data(i) % e_out(n_energy_out))
|
||||
allocate(this % data(t) % inelastic_data(i) % e_out_pdf(n_energy_out))
|
||||
allocate(this % data(t) % inelastic_data(i) % e_out_cdf(n_energy_out))
|
||||
|
||||
! Copy outgoing energy distribution
|
||||
this % data(t) % inelastic_data(i) % e_out(:) = edist % e_out
|
||||
this % data(t) % inelastic_data(i) % e_out_pdf(:) = edist % p
|
||||
this % data(t) % inelastic_data(i) % e_out_cdf(:) = edist % c
|
||||
|
||||
do j = 1, n_energy_out
|
||||
select type (adist => edist % angle(j) % obj)
|
||||
type is (Tabular)
|
||||
! On first pass, allocate space for angles
|
||||
if (j == 1) then
|
||||
n_mu = size(adist % x)
|
||||
this % data(t) % n_inelastic_mu = n_mu
|
||||
allocate(this % data(t) % inelastic_data(i) % mu(&
|
||||
n_mu, n_energy_out))
|
||||
end if
|
||||
|
||||
! Copy outgoing angles
|
||||
this % data(t) % inelastic_data(i) % mu(:, j) = adist % x
|
||||
end select
|
||||
end do
|
||||
end associate
|
||||
end do
|
||||
|
||||
! Clear data on correlated angle-energy object
|
||||
deallocate(correlated_dist % breakpoints)
|
||||
deallocate(correlated_dist % interpolation)
|
||||
deallocate(correlated_dist % energy)
|
||||
deallocate(correlated_dist % distribution)
|
||||
end if
|
||||
|
||||
call close_group(inelastic_group)
|
||||
end if
|
||||
call close_group(T_group)
|
||||
end do
|
||||
|
||||
call close_group(kT_group)
|
||||
end subroutine salphabeta_from_hdf5
|
||||
end subroutine from_hdf5
|
||||
|
||||
!===============================================================================
|
||||
! SAB_CALCULATE_XS determines the elastic and inelastic scattering
|
||||
! cross-sections in the thermal energy range.
|
||||
!===============================================================================
|
||||
|
||||
subroutine sab_calculate_xs(this, E, sqrtkT, i_temp, elastic, inelastic)
|
||||
subroutine calculate_xs(this, E, sqrtkT, i_temp, elastic, inelastic)
|
||||
class(SAlphaBeta), intent(in) :: this ! S(a,b) object
|
||||
real(8), intent(in) :: E ! energy
|
||||
real(8), intent(in) :: sqrtkT ! temperature
|
||||
real(C_DOUBLE), intent(in) :: E ! energy
|
||||
real(C_DOUBLE), intent(in) :: sqrtkT ! temperature
|
||||
integer, intent(out) :: i_temp ! index in the S(a,b)'s temperature
|
||||
real(8), intent(out) :: elastic ! thermal elastic cross section
|
||||
real(8), intent(out) :: inelastic ! thermal inelastic cross section
|
||||
real(C_DOUBLE), intent(out) :: elastic ! thermal elastic cross section
|
||||
real(C_DOUBLE), intent(out) :: inelastic ! thermal inelastic cross section
|
||||
|
||||
integer :: i_grid ! index on S(a,b) energy grid
|
||||
real(8) :: f ! interp factor on S(a,b) energy grid
|
||||
real(8) :: kT
|
||||
call sab_calculate_xs(this % ptr, E, sqrtkT, i_temp, elastic, inelastic)
|
||||
end subroutine
|
||||
|
||||
! Determine temperature for S(a,b) table
|
||||
kT = sqrtkT**2
|
||||
if (temperature_method == TEMPERATURE_NEAREST) then
|
||||
! If using nearest temperature, do linear search on temperature
|
||||
do i_temp = 1, size(this % kTs)
|
||||
if (abs(this % kTs(i_temp) - kT) < &
|
||||
K_BOLTZMANN*temperature_tolerance) exit
|
||||
end do
|
||||
else
|
||||
! Find temperatures that bound the actual temperature
|
||||
do i_temp = 1, size(this % kTs) - 1
|
||||
if (this % kTs(i_temp) <= kT .and. &
|
||||
kT < this % kTs(i_temp + 1)) exit
|
||||
end do
|
||||
subroutine free(this)
|
||||
class(SAlphaBeta), intent(inout) :: this
|
||||
call sab_free(this % ptr)
|
||||
end subroutine
|
||||
|
||||
! Randomly sample between temperature i and i+1
|
||||
f = (kT - this % kTs(i_temp)) / &
|
||||
(this % kTs(i_temp + 1) - this % kTs(i_temp))
|
||||
if (f > prn()) i_temp = i_temp + 1
|
||||
end if
|
||||
function has_nuclide(this, name) result(val)
|
||||
class(SAlphaBeta), intent(in) :: this
|
||||
character(len=*), intent(in) :: name
|
||||
logical(C_BOOL) :: val
|
||||
|
||||
val = sab_has_nuclide(this % ptr, to_c_string(name))
|
||||
end function
|
||||
|
||||
! Get pointer to S(a,b) table
|
||||
associate (sab => this % data(i_temp))
|
||||
subroutine sample(this, micro_xs, E_in, E_out, mu)
|
||||
class(SAlphaBeta), intent(in) :: this
|
||||
type(C_PTR), value :: micro_xs
|
||||
real(C_DOUBLE), value :: E_in
|
||||
real(C_DOUBLE), intent(out) :: E_out
|
||||
real(C_DOUBLE), intent(out) :: mu
|
||||
|
||||
! Get index and interpolation factor for inelastic grid
|
||||
if (E < sab % inelastic_e_in(1)) then
|
||||
i_grid = 1
|
||||
f = ZERO
|
||||
else
|
||||
i_grid = binary_search(sab % inelastic_e_in, sab % n_inelastic_e_in, E)
|
||||
f = (E - sab%inelastic_e_in(i_grid)) / &
|
||||
(sab%inelastic_e_in(i_grid+1) - sab%inelastic_e_in(i_grid))
|
||||
end if
|
||||
|
||||
! Calculate S(a,b) inelastic scattering cross section
|
||||
inelastic = (ONE - f) * sab % inelastic_sigma(i_grid) + &
|
||||
f * sab % inelastic_sigma(i_grid + 1)
|
||||
|
||||
! Check for elastic data
|
||||
if (E < sab % threshold_elastic) then
|
||||
! Determine whether elastic scattering is given in the coherent or
|
||||
! incoherent approximation. For coherent, the cross section is
|
||||
! represented as P/E whereas for incoherent, it is simply P
|
||||
|
||||
if (sab % elastic_mode == SAB_ELASTIC_EXACT) then
|
||||
if (E < sab % elastic_e_in(1)) then
|
||||
! If energy is below that of the lowest Bragg peak, the elastic
|
||||
! cross section will be zero
|
||||
elastic = ZERO
|
||||
else
|
||||
i_grid = binary_search(sab % elastic_e_in, &
|
||||
sab % n_elastic_e_in, E)
|
||||
elastic = sab % elastic_P(i_grid) / E
|
||||
end if
|
||||
else
|
||||
! Determine index on elastic energy grid
|
||||
if (E < sab % elastic_e_in(1)) then
|
||||
i_grid = 1
|
||||
else
|
||||
i_grid = binary_search(sab % elastic_e_in, &
|
||||
sab % n_elastic_e_in, E)
|
||||
end if
|
||||
|
||||
! Get interpolation factor for elastic grid
|
||||
f = (E - sab%elastic_e_in(i_grid))/(sab%elastic_e_in(i_grid+1) - &
|
||||
sab%elastic_e_in(i_grid))
|
||||
|
||||
! Calculate S(a,b) elastic scattering cross section
|
||||
elastic = (ONE - f) * sab % elastic_P(i_grid) + &
|
||||
f * sab % elastic_P(i_grid + 1)
|
||||
end if
|
||||
else
|
||||
! No elastic data
|
||||
elastic = ZERO
|
||||
end if
|
||||
end associate
|
||||
|
||||
end subroutine sab_calculate_xs
|
||||
call sab_sample(this % ptr, micro_xs, E_in, E_out, mu)
|
||||
end subroutine
|
||||
|
||||
function threshold(this)
|
||||
class(SAlphaBeta), intent(in) :: this
|
||||
real(C_DOUBLE) :: threshold
|
||||
threshold = sab_threshold(this % ptr)
|
||||
end function
|
||||
|
||||
!===============================================================================
|
||||
! FREE_MEMORY_SAB deallocates global arrays defined in this module
|
||||
!===============================================================================
|
||||
|
||||
subroutine free_memory_sab()
|
||||
integer :: i
|
||||
n_sab_tables = 0
|
||||
if (allocated(sab_tables)) deallocate(sab_tables)
|
||||
if (allocated(sab_tables)) then
|
||||
do i = 1, size(sab_tables)
|
||||
call sab_tables(i) % free()
|
||||
end do
|
||||
deallocate(sab_tables)
|
||||
end if
|
||||
call sab_dict % clear()
|
||||
end subroutine free_memory_sab
|
||||
|
||||
|
|
|
|||
|
|
@ -1,302 +0,0 @@
|
|||
module secondary_correlated
|
||||
|
||||
use algorithm, only: binary_search
|
||||
use angleenergy_header, only: AngleEnergy
|
||||
use constants, only: ZERO, ONE, HALF, TWO, HISTOGRAM, LINEAR_LINEAR
|
||||
use distribution_univariate, only: DistributionContainer, Tabular
|
||||
use hdf5_interface
|
||||
use random_lcg, only: prn
|
||||
|
||||
!===============================================================================
|
||||
! CORRELATEDANGLEENERGY represents a correlated angle-energy distribution. This
|
||||
! corresponds to ACE law 61 and ENDF File 6, LAW=1, LANG/=2.
|
||||
!===============================================================================
|
||||
|
||||
type AngleEnergyTable
|
||||
integer :: interpolation
|
||||
integer :: n_discrete
|
||||
real(8), allocatable :: e_out(:)
|
||||
real(8), allocatable :: p(:)
|
||||
real(8), allocatable :: c(:)
|
||||
type(DistributionContainer), allocatable :: angle(:)
|
||||
end type AngleEnergyTable
|
||||
|
||||
type, extends(AngleEnergy) :: CorrelatedAngleEnergy
|
||||
integer :: n_region ! number of interpolation regions
|
||||
integer, allocatable :: breakpoints(:) ! breakpoints of interpolation regions
|
||||
integer, allocatable :: interpolation(:) ! interpolation region codes
|
||||
real(8), allocatable :: energy(:) ! incoming energies
|
||||
type(AngleEnergyTable), allocatable :: distribution(:) ! outgoing E/mu distributions
|
||||
contains
|
||||
procedure :: sample => correlated_sample
|
||||
procedure :: from_hdf5 => correlated_from_hdf5
|
||||
end type CorrelatedAngleEnergy
|
||||
|
||||
contains
|
||||
|
||||
subroutine correlated_sample(this, E_in, E_out, mu)
|
||||
class(CorrelatedAngleEnergy), intent(in) :: this
|
||||
real(8), intent(in) :: E_in ! incoming energy
|
||||
real(8), intent(out) :: E_out ! sampled outgoing energy
|
||||
real(8), intent(out) :: mu ! sapmled scattering cosine
|
||||
|
||||
integer :: i, k, l ! indices
|
||||
integer :: n_energy_in ! number of incoming energies
|
||||
integer :: n_energy_out ! number of outgoing energies
|
||||
real(8) :: r ! interpolation factor on incoming energy
|
||||
real(8) :: r1 ! random number on [0,1)
|
||||
real(8) :: frac ! interpolation factor on outgoing energy
|
||||
real(8) :: E_i_1, E_i_K ! endpoints on outgoing grid i
|
||||
real(8) :: E_i1_1, E_i1_K ! endpoints on outgoing grid i+1
|
||||
real(8) :: E_1, E_K ! endpoints interpolated between i and i+1
|
||||
real(8) :: E_l_k, E_l_k1 ! adjacent E on outgoing grid l
|
||||
real(8) :: p_l_k, p_l_k1 ! adjacent p on outgoing grid l
|
||||
real(8) :: c_k, c_k1 ! cumulative probability
|
||||
|
||||
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
|
||||
! Before the secondary distribution refactor, an isotropic polar cosine was
|
||||
! always sampled but then overwritten with the polar cosine sampled from the
|
||||
! correlated distribution. To preserve the random number stream, we keep
|
||||
! this dummy sampling here but can remove it later (will change answers)
|
||||
mu = TWO*prn() - ONE
|
||||
! <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
|
||||
|
||||
! find energy bin and calculate interpolation factor -- if the energy is
|
||||
! outside the range of the tabulated energies, choose the first or last bins
|
||||
n_energy_in = size(this%energy)
|
||||
if (E_in < this%energy(1)) then
|
||||
i = 1
|
||||
r = ZERO
|
||||
elseif (E_in > this%energy(n_energy_in)) then
|
||||
i = n_energy_in - 1
|
||||
r = ONE
|
||||
else
|
||||
i = binary_search(this%energy, n_energy_in, E_in)
|
||||
r = (E_in - this%energy(i)) / &
|
||||
(this%energy(i+1) - this%energy(i))
|
||||
end if
|
||||
|
||||
! Sample between the ith and (i+1)th bin
|
||||
if (r > prn()) then
|
||||
l = i + 1
|
||||
else
|
||||
l = i
|
||||
end if
|
||||
|
||||
! interpolation for energy E1 and EK
|
||||
n_energy_out = size(this%distribution(i)%e_out)
|
||||
E_i_1 = this%distribution(i)%e_out(1)
|
||||
E_i_K = this%distribution(i)%e_out(n_energy_out)
|
||||
|
||||
n_energy_out = size(this%distribution(i+1)%e_out)
|
||||
E_i1_1 = this%distribution(i+1)%e_out(1)
|
||||
E_i1_K = this%distribution(i+1)%e_out(n_energy_out)
|
||||
|
||||
E_1 = E_i_1 + r*(E_i1_1 - E_i_1)
|
||||
E_K = E_i_K + r*(E_i1_K - E_i_K)
|
||||
|
||||
! determine outgoing energy bin
|
||||
n_energy_out = size(this%distribution(l)%e_out)
|
||||
r1 = prn()
|
||||
c_k = this%distribution(l)%c(1)
|
||||
do k = 1, n_energy_out - 1
|
||||
c_k1 = this%distribution(l)%c(k+1)
|
||||
if (r1 < c_k1) exit
|
||||
c_k = c_k1
|
||||
end do
|
||||
|
||||
! check to make sure k is <= NP - 1
|
||||
k = min(k, n_energy_out - 1)
|
||||
|
||||
E_l_k = this%distribution(l)%e_out(k)
|
||||
p_l_k = this%distribution(l)%p(k)
|
||||
if (this%distribution(l)%interpolation == HISTOGRAM) then
|
||||
! Histogram interpolation
|
||||
if (p_l_k > ZERO) then
|
||||
E_out = E_l_k + (r1 - c_k)/p_l_k
|
||||
else
|
||||
E_out = E_l_k
|
||||
end if
|
||||
|
||||
elseif (this%distribution(l)%interpolation == LINEAR_LINEAR) then
|
||||
! Linear-linear interpolation
|
||||
E_l_k1 = this%distribution(l)%e_out(k+1)
|
||||
p_l_k1 = this%distribution(l)%p(k+1)
|
||||
|
||||
frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k)
|
||||
if (frac == ZERO) then
|
||||
E_out = E_l_k + (r1 - c_k)/p_l_k
|
||||
else
|
||||
E_out = E_l_k + (sqrt(max(ZERO, p_l_k*p_l_k + &
|
||||
TWO*frac*(r1 - c_k))) - p_l_k)/frac
|
||||
end if
|
||||
end if
|
||||
|
||||
! Now interpolate between incident energy bins i and i + 1
|
||||
if (l == i) then
|
||||
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
|
||||
else
|
||||
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
|
||||
end if
|
||||
|
||||
! Find correlated angular distribution for closest outgoing energy bin
|
||||
if (r1 - c_k < c_k1 - r1) then
|
||||
mu = this%distribution(l)%angle(k)%obj%sample()
|
||||
else
|
||||
mu = this%distribution(l)%angle(k + 1)%obj%sample()
|
||||
end if
|
||||
end subroutine correlated_sample
|
||||
|
||||
subroutine correlated_from_hdf5(this, group_id)
|
||||
class(CorrelatedAngleEnergy), intent(inout) :: this
|
||||
integer(HID_T), intent(in) :: group_id
|
||||
|
||||
integer :: i, j, k
|
||||
integer :: n_energy
|
||||
integer :: m, n
|
||||
integer :: offset_mu
|
||||
integer :: interp_mu
|
||||
integer(HID_T) :: dset_id
|
||||
integer(HSIZE_T) :: dims(1), dims2(2)
|
||||
integer, allocatable :: temp(:,:)
|
||||
integer, allocatable :: offsets(:)
|
||||
integer, allocatable :: interp(:)
|
||||
integer, allocatable :: n_discrete(:)
|
||||
real(8), allocatable :: eout(:,:)
|
||||
real(8), allocatable :: mu(:,:)
|
||||
|
||||
! Open incoming energy dataset
|
||||
dset_id = open_dataset(group_id, 'energy')
|
||||
|
||||
! Get interpolation parameters
|
||||
call read_attribute(temp, dset_id, 'interpolation')
|
||||
allocate(this%breakpoints(size(temp, 1)))
|
||||
allocate(this%interpolation(size(temp, 1)))
|
||||
this%breakpoints(:) = temp(:, 1)
|
||||
this%interpolation(:) = temp(:, 2)
|
||||
this%n_region = size(this%breakpoints)
|
||||
deallocate(temp)
|
||||
|
||||
! Get incoming energies
|
||||
call get_shape(dset_id, dims)
|
||||
n_energy = int(dims(1), 4)
|
||||
allocate(this%energy(n_energy))
|
||||
allocate(this%distribution(n_energy))
|
||||
call read_dataset(this%energy, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
! Get outgoing energy distribution data
|
||||
dset_id = open_dataset(group_id, 'energy_out')
|
||||
call read_attribute(offsets, dset_id, 'offsets')
|
||||
call read_attribute(interp, dset_id, 'interpolation')
|
||||
call read_attribute(n_discrete, dset_id, 'n_discrete_lines')
|
||||
call get_shape(dset_id, dims2)
|
||||
allocate(eout(dims2(1), dims2(2)))
|
||||
call read_dataset(eout, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
! Get outgoing angle data
|
||||
dset_id = open_dataset(group_id, 'mu')
|
||||
call get_shape(dset_id, dims2)
|
||||
allocate(mu(dims2(1), dims2(2)))
|
||||
call read_dataset(mu, dset_id)
|
||||
call close_dataset(dset_id)
|
||||
|
||||
do i = 1, n_energy
|
||||
! Determine number of outgoing energies
|
||||
j = offsets(i)
|
||||
if (i < n_energy) then
|
||||
n = offsets(i+1) - j
|
||||
else
|
||||
n = size(eout, 1) - j
|
||||
end if
|
||||
|
||||
associate (d => this % distribution(i))
|
||||
! Assign interpolation scheme and number of discrete lines
|
||||
d % interpolation = interp(i)
|
||||
d % n_discrete = n_discrete(i)
|
||||
|
||||
! Allocate arrays for energies and PDF/CDF
|
||||
allocate(d % e_out(n))
|
||||
allocate(d % p(n))
|
||||
allocate(d % c(n))
|
||||
allocate(d % angle(n))
|
||||
|
||||
! Copy data
|
||||
d % e_out(:) = eout(j+1:j+n, 1)
|
||||
d % p(:) = eout(j+1:j+n, 2)
|
||||
d % c(:) = eout(j+1:j+n, 3)
|
||||
|
||||
! To get answers that match ACE data, for now we still use the tabulated
|
||||
! CDF values that were passed through to the HDF5 library. At a later
|
||||
! time, we can remove the CDF values from the HDF5 library and
|
||||
! reconstruct them using the PDF
|
||||
if (.false.) then
|
||||
! Calculate cumulative distribution function -- discrete portion
|
||||
do k = 1, d % n_discrete
|
||||
if (k == 1) then
|
||||
d % c(k) = d % p(k)
|
||||
else
|
||||
d % c(k) = d % c(k-1) + d % p(k)
|
||||
end if
|
||||
end do
|
||||
|
||||
! Continuous portion
|
||||
do k = d % n_discrete + 1, n
|
||||
if (k == d % n_discrete + 1) then
|
||||
d % c(k) = sum(d % p(1:d % n_discrete))
|
||||
else
|
||||
if (d % interpolation == HISTOGRAM) then
|
||||
d % c(k) = d % c(k-1) + d % p(k-1) * &
|
||||
(d % e_out(k) - d % e_out(k-1))
|
||||
elseif (d % interpolation == LINEAR_LINEAR) then
|
||||
d % c(k) = d % c(k-1) + HALF*(d % p(k-1) + d % p(k)) * &
|
||||
(d % e_out(k) - d % e_out(k-1))
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
|
||||
! Normalize density and distribution functions
|
||||
d % p(:) = d % p(:)/d % c(n)
|
||||
d % c(:) = d % c(:)/d % c(n)
|
||||
end if
|
||||
end associate
|
||||
|
||||
do j = 1, n
|
||||
allocate(Tabular :: this%distribution(i)%angle(j)%obj)
|
||||
select type(mudist => this%distribution(i)%angle(j)%obj)
|
||||
type is (Tabular)
|
||||
! Get interpolation scheme
|
||||
interp_mu = nint(eout(offsets(i)+j, 4))
|
||||
|
||||
! Determine offset and size of distribution
|
||||
offset_mu = nint(eout(offsets(i)+j, 5))
|
||||
if (offsets(i) + j < size(eout, 1)) then
|
||||
m = nint(eout(offsets(i)+j+1, 5)) - offset_mu
|
||||
else
|
||||
m = size(mu, 1) - offset_mu
|
||||
end if
|
||||
|
||||
! To get answers that match ACE data, for now we still use the tabulated
|
||||
! CDF values that were passed through to the HDF5 library. At a later
|
||||
! time, we can remove the CDF values from the HDF5 library and
|
||||
! reconstruct them using the PDF
|
||||
if (.true.) then
|
||||
mudist % interpolation = interp_mu
|
||||
allocate(mudist % x(m))
|
||||
allocate(mudist % p(m))
|
||||
allocate(mudist % c(m))
|
||||
mudist % x(:) = mu(offset_mu+1:offset_mu+m, 1)
|
||||
mudist % p(:) = mu(offset_mu+1:offset_mu+m, 2)
|
||||
mudist % c(:) = mu(offset_mu+1:offset_mu+m, 3)
|
||||
else
|
||||
! Initialize tabular distribution
|
||||
call mudist % initialize(mu(offset_mu+1:offset_mu+m, 1), &
|
||||
mu(offset_mu+1:offset_mu+m, 2), interp_mu)
|
||||
end if
|
||||
end select
|
||||
end do
|
||||
end do
|
||||
end subroutine correlated_from_hdf5
|
||||
|
||||
end module secondary_correlated
|
||||
|
|
@ -21,14 +21,6 @@ namespace openmc {
|
|||
|
||||
class CorrelatedAngleEnergy : public AngleEnergy {
|
||||
public:
|
||||
explicit CorrelatedAngleEnergy(hid_t group);
|
||||
|
||||
//! Sample distribution for an angle and energy
|
||||
//! \param[in] E_in Incoming energy in [eV]
|
||||
//! \param[out] E_out Outgoing energy in [eV]
|
||||
//! \param[out] mu Outgoing cosine with respect to current direction
|
||||
void sample(double E_in, double& E_out, double& mu) const;
|
||||
private:
|
||||
//! Outgoing energy/angle at a single incoming energy
|
||||
struct CorrTable {
|
||||
int n_discrete; //!< Number of discrete lines
|
||||
|
|
@ -39,6 +31,22 @@ private:
|
|||
std::vector<UPtrDist> angle; //!< Angle distribution
|
||||
};
|
||||
|
||||
explicit CorrelatedAngleEnergy(hid_t group);
|
||||
|
||||
//! Sample distribution for an angle and energy
|
||||
//! \param[in] E_in Incoming energy in [eV]
|
||||
//! \param[out] E_out Outgoing energy in [eV]
|
||||
//! \param[out] mu Outgoing cosine with respect to current direction
|
||||
void sample(double E_in, double& E_out, double& mu) const;
|
||||
|
||||
// energy property
|
||||
std::vector<double>& energy() { return energy_; }
|
||||
const std::vector<double>& energy() const { return energy_; }
|
||||
|
||||
// distribution property
|
||||
std::vector<CorrTable>& distribution() { return distribution_; }
|
||||
const std::vector<CorrTable>& distribution() const { return distribution_; }
|
||||
private:
|
||||
int n_region_; //!< Number of interpolation regions
|
||||
std::vector<int> breakpoints_; //!< Breakpoints between regions
|
||||
std::vector<Interpolation> interpolation_; //!< Interpolation laws
|
||||
|
|
|
|||
|
|
@ -1,5 +1,6 @@
|
|||
#include "settings.h"
|
||||
|
||||
#include "constants.h"
|
||||
#include "error.h"
|
||||
#include "openmc.h"
|
||||
#include "string_utils.h"
|
||||
|
|
@ -20,6 +21,12 @@ std::string path_multipole;
|
|||
std::string path_output;
|
||||
std::string path_source;
|
||||
|
||||
int temperature_method {TEMPERATURE_NEAREST};
|
||||
bool temperature_multipole {false};
|
||||
double temperature_tolerance {10.0};
|
||||
double temperature_default {293.6};
|
||||
std::array<double, 2> temperature_range {0.0, 0.0};
|
||||
|
||||
//==============================================================================
|
||||
// Functions
|
||||
//==============================================================================
|
||||
|
|
@ -65,6 +72,32 @@ void read_settings(pugi::xml_node* root)
|
|||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Get temperature settings
|
||||
if (check_for_node(*root, "temperature_default")) {
|
||||
temperature_default = std::stod(get_node_value(*root, "temperature_default"));
|
||||
}
|
||||
if (check_for_node(*root, "temperature_method")) {
|
||||
auto temp_str = get_node_value(*root, "temperature_method", true, true);
|
||||
if (temp_str == "nearest") {
|
||||
temperature_method = TEMPERATURE_NEAREST;
|
||||
} else if (temp_str == "interpolation") {
|
||||
temperature_method = TEMPERATURE_INTERPOLATION;
|
||||
} else {
|
||||
fatal_error("Unknown temperature method: " + temp_str);
|
||||
}
|
||||
}
|
||||
if (check_for_node(*root, "temperature_tolerance")) {
|
||||
temperature_tolerance = std::stod(get_node_value(*root, "temperature_tolerance"));
|
||||
}
|
||||
if (check_for_node(*root, "temperature_multipole")) {
|
||||
temperature_multipole = get_node_value_bool(*root, "temperature_multipole");
|
||||
}
|
||||
if (check_for_node(*root, "temperature_range")) {
|
||||
auto range = get_node_array<double>(*root, "temperature_range");
|
||||
temperature_range[0] = range[0];
|
||||
temperature_range[1] = range[1];
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -4,6 +4,7 @@
|
|||
//! \file settings.h
|
||||
//! \brief Settings for OpenMC
|
||||
|
||||
#include <array>
|
||||
#include <string>
|
||||
|
||||
#include "pugixml.hpp"
|
||||
|
|
@ -31,6 +32,12 @@ extern std::string path_multipole;
|
|||
extern std::string path_output;
|
||||
extern std::string path_source;
|
||||
|
||||
extern int temperature_method;
|
||||
extern bool temperature_multipole;
|
||||
extern double temperature_tolerance;
|
||||
extern double temperature_default;
|
||||
extern std::array<double, 2> temperature_range;
|
||||
|
||||
//==============================================================================
|
||||
//! Read settings from XML file
|
||||
//! \param[in] root XML node for <settings>
|
||||
|
|
@ -40,4 +47,4 @@ extern "C" void read_settings(pugi::xml_node* root);
|
|||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_SETTINGS_H
|
||||
#endif // OPENMC_SETTINGS_H
|
||||
|
|
|
|||
623
src/thermal.cpp
Normal file
623
src/thermal.cpp
Normal file
|
|
@ -0,0 +1,623 @@
|
|||
#include "thermal.h"
|
||||
|
||||
#include <algorithm> // for sort, move, min, max, find
|
||||
#include <cmath> // for round, sqrt, fabs
|
||||
#include <sstream> // for stringstream
|
||||
|
||||
#include "xtensor/xarray.hpp"
|
||||
#include "xtensor/xbuilder.hpp"
|
||||
#include "xtensor/xmath.hpp"
|
||||
#include "xtensor/xsort.hpp"
|
||||
#include "xtensor/xtensor.hpp"
|
||||
#include "xtensor/xview.hpp"
|
||||
|
||||
#include "constants.h"
|
||||
#include "error.h"
|
||||
#include "random_lcg.h"
|
||||
#include "search.h"
|
||||
#include "secondary_correlated.h"
|
||||
#include "settings.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
// ThermalScattering implementation
|
||||
//==============================================================================
|
||||
|
||||
ThermalScattering::ThermalScattering(hid_t group, const std::vector<double>& temperature,
|
||||
int method, double tolerance, const double* minmax)
|
||||
{
|
||||
// Get name of table from group
|
||||
name_ = object_name(group);
|
||||
|
||||
// Get rid of leading '/'
|
||||
name_ = name_.substr(1);
|
||||
|
||||
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");
|
||||
|
||||
// Determine temperatures available
|
||||
auto dset_names = dataset_names(kT_group);
|
||||
auto n = dset_names.size();
|
||||
auto temps_available = xt::empty<double>({n});
|
||||
for (int i = 0; i < dset_names.size(); ++i) {
|
||||
// Read temperature value
|
||||
double T;
|
||||
read_dataset(kT_group, dset_names[i].data(), T);
|
||||
temps_available[i] = T / K_BOLTZMANN;
|
||||
}
|
||||
std::sort(temps_available.begin(), temps_available.end());
|
||||
|
||||
// Determine actual temperatures to read -- start by checking whether a
|
||||
// temperature range was given, in which case all temperatures in the range
|
||||
// are loaded irrespective of what temperatures actually appear in the model
|
||||
std::vector<int> temps_to_read;
|
||||
if (minmax[1] > 0.0) {
|
||||
for (const auto& T : temps_available) {
|
||||
if (minmax[0] <= T && T <= minmax[1]) {
|
||||
temps_to_read.push_back(std::round(T));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
switch (method) {
|
||||
case TEMPERATURE_NEAREST:
|
||||
// Determine actual temperatures to read
|
||||
for (const auto& T : temperature) {
|
||||
|
||||
auto i_closest = xt::argmin(xt::abs(temps_available - T))[0];
|
||||
auto temp_actual = temps_available[i_closest];
|
||||
if (std::fabs(temp_actual - T) < tolerance) {
|
||||
if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temp_actual))
|
||||
== temps_to_read.end()) {
|
||||
temps_to_read.push_back(std::round(temp_actual));
|
||||
}
|
||||
} else {
|
||||
std::stringstream msg;
|
||||
msg << "Nuclear data library does not contain cross sections for "
|
||||
<< name_ << " at or near " << std::round(T) << " K.";
|
||||
fatal_error(msg);
|
||||
}
|
||||
}
|
||||
break;
|
||||
|
||||
case TEMPERATURE_INTERPOLATION:
|
||||
// If temperature interpolation or multipole is selected, get a list of
|
||||
// bounding temperatures for each actual temperature present in the model
|
||||
for (const auto& T : temperature) {
|
||||
bool found = false;
|
||||
for (int j = 0; j < temps_available.size() - 1; ++j) {
|
||||
if (temps_available[j] <= T && T < temps_available[j + 1]) {
|
||||
int T_j = std::round(temps_available[j]);
|
||||
int T_j1 = std::round(temps_available[j + 1]);
|
||||
if (std::find(temps_to_read.begin(), temps_to_read.end(), T_j) == temps_to_read.end()) {
|
||||
temps_to_read.push_back(T_j);
|
||||
}
|
||||
if (std::find(temps_to_read.begin(), temps_to_read.end(), T_j1) == temps_to_read.end()) {
|
||||
temps_to_read.push_back(T_j1);
|
||||
}
|
||||
found = true;
|
||||
}
|
||||
}
|
||||
if (!found) {
|
||||
std::stringstream msg;
|
||||
msg << "Nuclear data library does not contain cross sections for "
|
||||
<< name_ << " at temperatures that bound " << std::round(T) << " K.";
|
||||
fatal_error(msg);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Sort temperatures to read
|
||||
std::sort(temps_to_read.begin(), temps_to_read.end());
|
||||
|
||||
auto n_temperature = temps_to_read.size();
|
||||
kTs_.reserve(n_temperature);
|
||||
data_.reserve(n_temperature);
|
||||
|
||||
for (auto T : temps_to_read) {
|
||||
// Get temperature as a string
|
||||
std::string temp_str = std::to_string(T) + "K";
|
||||
|
||||
// Read exact temperature value
|
||||
double kT;
|
||||
read_dataset(kT_group, temp_str.data(), kT);
|
||||
kTs_.push_back(kT);
|
||||
|
||||
// Open group for temperature i
|
||||
hid_t T_group = open_group(group, temp_str.data());
|
||||
data_.emplace_back(T_group, secondary_mode);
|
||||
close_group(T_group);
|
||||
}
|
||||
|
||||
close_group(kT_group);
|
||||
}
|
||||
|
||||
void
|
||||
ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
|
||||
double* elastic, double* inelastic) const
|
||||
{
|
||||
// Determine temperature for S(a,b) table
|
||||
double kT = sqrtkT*sqrtkT;
|
||||
int i;
|
||||
if (temperature_method == TEMPERATURE_NEAREST) {
|
||||
// If using nearest temperature, do linear search on temperature
|
||||
for (i = 0; i < kTs_.size(); ++i) {
|
||||
if (abs(kTs_[i] - kT) < K_BOLTZMANN*temperature_tolerance) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// Find temperatures that bound the actual temperature
|
||||
for (i = 0; i < kTs_.size() - 1; ++i) {
|
||||
if (kTs_[i] <= kT && kT < kTs_[i+1]) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
// Randomly sample between temperature i and i+1
|
||||
double f = (kT - kTs_[i]) / (kTs_[i+1] - kTs_[i]);
|
||||
if (f > prn()) ++i;
|
||||
}
|
||||
|
||||
// Set temperature index
|
||||
*i_temp = i;
|
||||
|
||||
// Get pointer to S(a,b) table
|
||||
auto& sab = data_[i];
|
||||
|
||||
// Get index and interpolation factor for inelastic grid
|
||||
int i_grid = 0;
|
||||
double f = 0.0;
|
||||
if (E >= sab.inelastic_e_in_.front()) {
|
||||
auto& E_in = sab.inelastic_e_in_;
|
||||
i_grid = lower_bound_index(E_in.begin(), E_in.end(), E);
|
||||
f = (E - E_in[i_grid]) / (E_in[i_grid+1] - E_in[i_grid]);
|
||||
}
|
||||
|
||||
// Calculate S(a,b) inelastic scattering cross section
|
||||
auto& xs = sab.inelastic_sigma_;
|
||||
*inelastic = (1.0 - f) * xs[i_grid] + f * xs[i_grid + 1];
|
||||
|
||||
// Check for elastic data
|
||||
if (E < sab.threshold_elastic_) {
|
||||
// Determine whether elastic scattering is given in the coherent or
|
||||
// incoherent approximation. For coherent, the cross section is
|
||||
// represented as P/E whereas for incoherent, it is simply P
|
||||
|
||||
auto& E_in = sab.elastic_e_in_;
|
||||
|
||||
if (sab.elastic_mode_ == SAB_ELASTIC_COHERENT) {
|
||||
if (E < E_in.front()) {
|
||||
// If energy is below that of the lowest Bragg peak, the elastic
|
||||
// cross section will be zero
|
||||
*elastic = 0.0;
|
||||
} else {
|
||||
i_grid = lower_bound_index(E_in.begin(), E_in.end(), E);
|
||||
*elastic = sab.elastic_P_[i_grid] / E;
|
||||
}
|
||||
} else {
|
||||
// Determine index on elastic energy grid
|
||||
if (E < E_in.front()) {
|
||||
i_grid = 0;
|
||||
} else {
|
||||
i_grid = lower_bound_index(E_in.begin(), E_in.end(), E);
|
||||
}
|
||||
|
||||
// Get interpolation factor for elastic grid
|
||||
f = (E - E_in[i_grid]) / (E_in[i_grid+1] - E_in[i_grid]);
|
||||
|
||||
// Calculate S(a,b) elastic scattering cross section
|
||||
auto& xs = sab.elastic_P_;
|
||||
*elastic = (1.0 - f) * xs[i_grid] + f * xs[i_grid + 1];
|
||||
}
|
||||
} else {
|
||||
// No elastic data
|
||||
*elastic = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
bool
|
||||
ThermalScattering::has_nuclide(const char* name) const
|
||||
{
|
||||
std::string nuc {name};
|
||||
return std::find(nuclides_.begin(), nuclides_.end(), nuc) != nuclides_.end();
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// ThermalData implementation
|
||||
//==============================================================================
|
||||
|
||||
ThermalData::ThermalData(hid_t group, int secondary_mode)
|
||||
: inelastic_mode_{secondary_mode}
|
||||
{
|
||||
// Coherent 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);
|
||||
|
||||
// Get view on energies and cross section/probability values
|
||||
auto E_in = xt::view(temp, 0);
|
||||
auto P = xt::view(temp, 1);
|
||||
|
||||
// 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();
|
||||
|
||||
// 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);
|
||||
|
||||
// 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;
|
||||
}
|
||||
|
||||
close_group(elastic_group);
|
||||
}
|
||||
|
||||
// Inelastic data
|
||||
if (object_exists(group, "inelastic")) {
|
||||
// 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();
|
||||
|
||||
// Set inelastic threshold
|
||||
threshold_inelastic_ = inelastic_e_in_.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));
|
||||
}
|
||||
}
|
||||
|
||||
close_group(inelastic_group);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
ThermalData::sample(const NuclideMicroXS* micro_xs, double E,
|
||||
double* E_out, double* mu)
|
||||
{
|
||||
// Determine whether inelastic or elastic scattering will occur
|
||||
if (prn() < micro_xs->thermal_elastic / micro_xs->thermal) {
|
||||
// elastic scattering
|
||||
|
||||
// Get index and interpolation factor for elastic grid
|
||||
int i = 0;
|
||||
double f = 0.0;
|
||||
if (E >= elastic_e_in_.front()) {
|
||||
auto& E_in = elastic_e_in_;
|
||||
i = lower_bound_index(E_in.begin(), E_in.end(), E);
|
||||
f = (E - E_in[i]) / (E_in[i+1] - E_in[i]);
|
||||
}
|
||||
|
||||
// Select treatment based on elastic mode
|
||||
if (elastic_mode_ == SAB_ELASTIC_INCOHERENT) {
|
||||
// With this treatment, we interpolate between two discrete cosines
|
||||
// corresponding to neighboring incoming energies. This is used for
|
||||
// data derived in the incoherent approximation
|
||||
|
||||
// Sample outgoing cosine bin
|
||||
int k = prn() * n_elastic_mu_;
|
||||
|
||||
// Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
|
||||
double mu_ijk = elastic_mu_(i, k);
|
||||
double mu_i1jk = elastic_mu_(i+1, k);
|
||||
|
||||
// Cosine of angle between incoming and outgoing neutron
|
||||
*mu = (1 - f)*mu_ijk + f*mu_i1jk;
|
||||
|
||||
} else if (elastic_mode_ == SAB_ELASTIC_COHERENT) {
|
||||
// This treatment is used for data derived in the coherent
|
||||
// approximation, i.e. for crystalline structures that have Bragg
|
||||
// edges.
|
||||
|
||||
// Sample a Bragg edge between 1 and i
|
||||
double prob = prn() * elastic_P_[i+1];
|
||||
int k = 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;
|
||||
|
||||
} else {
|
||||
// Perform inelastic calculations
|
||||
|
||||
// Get index and interpolation factor for inelastic grid
|
||||
int i = 0;
|
||||
double f = 0.0;
|
||||
if (E >= inelastic_e_in_.front()) {
|
||||
auto& E_in = inelastic_e_in_;
|
||||
i = lower_bound_index(E_in.begin(), E_in.end(), E);
|
||||
f = (E - E_in[i]) / (E_in[i+1] - E_in[i]);
|
||||
}
|
||||
|
||||
// Now that we have an incoming energy bin, we need to determine the
|
||||
// outgoing energy bin. This will depend on the "secondary energy
|
||||
// mode". If the mode is 0, then the outgoing energy bin is chosen from a
|
||||
// set of equally-likely bins. If the mode is 1, then the first
|
||||
// two and last two bins are skewed to have lower probabilities than the
|
||||
// other bins (0.1 for the first and last bins and 0.4 for the second and
|
||||
// second to last bins, relative to a normal bin probability of 1).
|
||||
// Finally, if the mode is 2, then a continuous distribution (with
|
||||
// accompanying PDF and CDF is utilized)
|
||||
|
||||
if (inelastic_mode_ == SAB_SECONDARY_EQUAL ||
|
||||
inelastic_mode_ == SAB_SECONDARY_SKEWED) {
|
||||
int j;
|
||||
if (inelastic_mode_ == SAB_SECONDARY_EQUAL) {
|
||||
// All bins equally likely
|
||||
j = prn() * n_inelastic_e_out_;
|
||||
} else if (inelastic_mode_ == SAB_SECONDARY_SKEWED) {
|
||||
// Distribution skewed away from edge points
|
||||
double r = prn() * (n_inelastic_e_out_ - 3);
|
||||
if (r > 1.0) {
|
||||
// equally likely N-4 middle bins
|
||||
j = r + 1;
|
||||
} else if (r > 0.6) {
|
||||
// second to last bin has relative probability of 0.4
|
||||
j = n_inelastic_e_out_ - 2;
|
||||
} else if (r > 0.5) {
|
||||
// last bin has relative probability of 0.1
|
||||
j = n_inelastic_e_out_ - 1;
|
||||
} else if (r > 0.1) {
|
||||
// second bin has relative probability of 0.4
|
||||
j = 1;
|
||||
} else {
|
||||
// first bin has relative probability of 0.1
|
||||
j = 0;
|
||||
}
|
||||
}
|
||||
|
||||
// Determine outgoing energy corresponding to E_in[i] and E_in[i+1]
|
||||
double E_ij = inelastic_e_out_(i, j);
|
||||
double E_i1j = inelastic_e_out_(i+1, j);
|
||||
|
||||
// Outgoing energy
|
||||
*E_out = (1 - f)*E_ij + f*E_i1j;
|
||||
|
||||
// Sample outgoing cosine bin
|
||||
int k = prn() * n_inelastic_mu_;
|
||||
|
||||
// Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
|
||||
double mu_ijk = inelastic_mu_(i, j, k);
|
||||
double mu_i1jk = inelastic_mu_(i+1, j, k);
|
||||
|
||||
// Cosine of angle between incoming and outgoing neutron
|
||||
*mu = (1 - f)*mu_ijk + f*mu_i1jk;
|
||||
|
||||
} else if (inelastic_mode_ == SAB_SECONDARY_CONT) {
|
||||
// Continuous secondary energy - this is to be similar to
|
||||
// Law 61 interpolation on outgoing energy
|
||||
|
||||
// Sample between ith and [i+1]th bin
|
||||
int l = f > prn() ? i + 1 : i;
|
||||
|
||||
// Determine endpoints on grid i
|
||||
auto n = inelastic_data_[i].e_out.size();
|
||||
double E_i_1 = inelastic_data_[i].e_out(0);
|
||||
double E_i_J = inelastic_data_[i].e_out(n);
|
||||
|
||||
// Determine endpoints on grid i + 1
|
||||
n = inelastic_data_[i].e_out.size();
|
||||
double E_i1_1 = inelastic_data_[i + 1].e_out(0);
|
||||
double E_i1_J = inelastic_data_[i + 1].e_out(n);
|
||||
|
||||
double E_1 = E_i_1 + f * (E_i1_1 - E_i_1);
|
||||
double E_J = E_i_J + f * (E_i1_J - E_i_J);
|
||||
|
||||
// Determine outgoing energy bin
|
||||
// (First reset n_energy_out to the right value)
|
||||
n = inelastic_data_[l].n_e_out;
|
||||
double r1 = prn();
|
||||
double c_j = inelastic_data_[l].e_out_cdf[0];
|
||||
double c_j1;
|
||||
std::size_t j;
|
||||
for (j = 0; j < n - 2; ++j) {
|
||||
c_j1 = inelastic_data_[l].e_out_cdf[j + 1];
|
||||
if (r1 < c_j1) break;
|
||||
c_j = c_j1;
|
||||
}
|
||||
|
||||
// check to make sure j is <= n_energy_out - 2
|
||||
j = std::min(j, n - 2);
|
||||
|
||||
// Get the data to interpolate between
|
||||
double E_l_j = inelastic_data_[l].e_out[j];
|
||||
double p_l_j = inelastic_data_[l].e_out_pdf[j];
|
||||
|
||||
// Next part assumes linear-linear interpolation in standard
|
||||
double E_l_j1 = inelastic_data_[l].e_out[j + 1];
|
||||
double p_l_j1 = inelastic_data_[l].e_out_pdf[j + 1];
|
||||
|
||||
// Find secondary energy (variable E)
|
||||
double frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j);
|
||||
if (frac == 0.0) {
|
||||
*E_out = E_l_j + (r1 - c_j) / p_l_j;
|
||||
} else {
|
||||
*E_out = E_l_j + (std::sqrt(std::max(0.0, p_l_j*p_l_j +
|
||||
2.0*frac*(r1 - c_j))) - p_l_j) / frac;
|
||||
}
|
||||
|
||||
// Now interpolate between incident energy bins i and i + 1
|
||||
if (l == i) {
|
||||
*E_out = E_1 + (E - E_i_1) * (E_J - E_1) / (E_i_J - E_i_1);
|
||||
} else {
|
||||
*E_out = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1);
|
||||
}
|
||||
|
||||
// Sample outgoing cosine bin
|
||||
std::size_t k = prn() * n_inelastic_mu_;
|
||||
|
||||
// Rather than use the sampled discrete mu directly, it is smeared over
|
||||
// a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
|
||||
// discrete mu value itself.
|
||||
const auto& mu_l = inelastic_data_[l].mu;
|
||||
f = (r1 - c_j)/(c_j1 - c_j);
|
||||
|
||||
// Determine (k-1)th mu value
|
||||
double mu_left;
|
||||
if (k == 0) {
|
||||
mu_left = -1.0;
|
||||
} else {
|
||||
mu_left = mu_l(j, k-1) + f*(mu_l(j+1, k-1) - mu_l(j, k-1));
|
||||
}
|
||||
|
||||
// Determine kth mu value
|
||||
*mu = mu_l(j, k) + f*(mu_l(j+1, k) - mu_l(j, k));
|
||||
|
||||
// Determine (k+1)th mu value
|
||||
double mu_right;
|
||||
if (k == n_inelastic_mu_ - 1) {
|
||||
mu_right = 1.0 - *mu;
|
||||
} else {
|
||||
mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1)) - *mu;
|
||||
}
|
||||
|
||||
// Smear angle
|
||||
*mu += std::min(*mu - mu_left, mu_right - *mu)*(prn() - 0.5);
|
||||
|
||||
} // (inelastic secondary energy treatment)
|
||||
} // (elastic or inelastic)
|
||||
|
||||
// Because of floating-point roundoff, it may be possible for mu to be
|
||||
// outside of the range [-1,1). In these cases, we just set mu to exactly
|
||||
// -1 or 1
|
||||
if (std::fabs(*mu) > 1.0) *mu = std::copysign(1.0, *mu);
|
||||
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Fortran compatibility functions
|
||||
//==============================================================================
|
||||
|
||||
ThermalScattering*
|
||||
sab_from_hdf5(hid_t group, const double* temperature, int n,
|
||||
int method, double tolerance, const double* minmax)
|
||||
{
|
||||
// Convert temperatures to a vector
|
||||
std::vector<double> T {temperature, temperature + n};
|
||||
|
||||
// Create new object and return it
|
||||
return new ThermalScattering{group, T, method, tolerance, minmax};
|
||||
}
|
||||
|
||||
void sab_calculate_xs(ThermalScattering* data, double E, double sqrtkT,
|
||||
int* i_temp, double* elastic, double* inelastic)
|
||||
{
|
||||
// Calculate cross section
|
||||
int t;
|
||||
data->calculate_xs(E, sqrtkT, &t, elastic, inelastic);
|
||||
|
||||
// Fortran needs index plus one
|
||||
*i_temp = t + 1;
|
||||
}
|
||||
|
||||
void sab_free(ThermalScattering* data) { delete data; }
|
||||
|
||||
bool sab_has_nuclide(ThermalScattering* data, const char* name)
|
||||
{
|
||||
return data->has_nuclide(name);
|
||||
}
|
||||
|
||||
void sab_sample(ThermalScattering* data, const NuclideMicroXS* micro_xs,
|
||||
double E_in, double* E_out, double* mu)
|
||||
{
|
||||
int i_temp = micro_xs->index_temp_sab;
|
||||
data->data_[i_temp - 1].sample(micro_xs, E_in, E_out, mu);
|
||||
}
|
||||
|
||||
double sab_threshold(ThermalScattering* data) { return data->threshold(); }
|
||||
|
||||
} // namespace openmc
|
||||
147
src/thermal.h
Normal file
147
src/thermal.h
Normal file
|
|
@ -0,0 +1,147 @@
|
|||
#ifndef OPENMC_THERMAL_SCATTERING_H
|
||||
#define OPENMC_THERMAL_SCATTERING_H
|
||||
|
||||
#include <cstddef>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
#include "xtensor/xtensor.hpp"
|
||||
|
||||
#include "hdf5_interface.h"
|
||||
#include "nuclide.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
// Constants
|
||||
//==============================================================================
|
||||
|
||||
// Secondary energy mode for S(a,b) inelastic scattering
|
||||
// TODO: Convert to enum
|
||||
constexpr int SAB_SECONDARY_EQUAL {0}; // Equally-likely outgoing energy bins
|
||||
constexpr int SAB_SECONDARY_SKEWED {1}; // Skewed outgoing energy bins
|
||||
constexpr int SAB_SECONDARY_CONT {2}; // Continuous, linear-linear interpolation
|
||||
|
||||
// Elastic mode for S(a,b) elastic scattering
|
||||
// TODO: Convert to enum
|
||||
constexpr int SAB_ELASTIC_INCOHERENT {3}; // Incoherent elastic scattering
|
||||
constexpr int SAB_ELASTIC_COHERENT {4}; // Coherent elastic scattering (Bragg edges)
|
||||
|
||||
//==============================================================================
|
||||
//! Secondary angle-energy data for thermal neutron scattering at a single
|
||||
//! temperature
|
||||
//==============================================================================
|
||||
|
||||
class ThermalData {
|
||||
public:
|
||||
ThermalData(hid_t group, int secondary_mode);
|
||||
|
||||
// Sample an outgoing energy and angle
|
||||
void sample(const NuclideMicroXS* micro_xs, double E_in,
|
||||
double* E_out, double* mu);
|
||||
private:
|
||||
//! 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
|
||||
};
|
||||
|
||||
//! Upper threshold for incoherent inelastic scattering (usually ~4 eV)
|
||||
double threshold_inelastic_;
|
||||
//! Upper threshold for coherent/incoherent elastic scattering
|
||||
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
|
||||
|
||||
// ThermalScattering needs access to private data members
|
||||
friend class ThermalScattering;
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
//! Data for thermal neutron scattering, typically off light isotopes in
|
||||
//! moderating materials such as water, graphite, BeO, etc.
|
||||
//==============================================================================
|
||||
|
||||
class ThermalScattering {
|
||||
public:
|
||||
ThermalScattering(hid_t group, const std::vector<double>& temperature, int method,
|
||||
double tolerance, const double* minmax);
|
||||
|
||||
//! Determine inelastic/elastic cross section at given energy
|
||||
//!
|
||||
//! \param[in] E incoming energy in [eV]
|
||||
//! \param[in] sqrtkT square-root of temperature multipled by Boltzmann's constant
|
||||
//! \param[out] i_temp corresponding temperature index
|
||||
//! \param[out] elastic Thermal elastic scattering cross section
|
||||
//! \param[out] inelastic Thermal inelastic scattering cross section
|
||||
void calculate_xs(double E, double sqrtkT, int* i_temp, double* elastic,
|
||||
double* inelastic) const;
|
||||
|
||||
//! Determine whether table applies to a particular nuclide
|
||||
//!
|
||||
//! \param[in] name Name of the nuclide, e.g., "H1"
|
||||
//! \return Whether table applies to the nuclide
|
||||
bool has_nuclide(const char* name) const;
|
||||
|
||||
// Sample an outgoing energy and angle
|
||||
void sample(const NuclideMicroXS* micro_xs, double E_in,
|
||||
double* E_out, double* mu);
|
||||
|
||||
double threshold() const { return data_[0].threshold_inelastic_; }
|
||||
|
||||
std::string name_; //!< name of table, e.g. "c_H_in_H2O"
|
||||
double awr_; //!< weight of nucleus in neutron masses
|
||||
std::vector<double> kTs_; //!< temperatures in [eV] (k*T)
|
||||
std::vector<std::string> nuclides_; //!< Valid nuclides
|
||||
|
||||
//! cross sections and distributions at each temperature
|
||||
std::vector<ThermalData> data_;
|
||||
};
|
||||
|
||||
//==============================================================================
|
||||
// Fortran compatibility functions
|
||||
//==============================================================================
|
||||
|
||||
extern "C" {
|
||||
ThermalScattering* sab_from_hdf5(hid_t group, const double* temperature,
|
||||
int n, int method, double tolerance, const double* minmax);
|
||||
void sab_calculate_xs(ThermalScattering* data, double E, double sqrtkT,
|
||||
int* i_temp, double* elastic, double* inelastic);
|
||||
void sab_free(ThermalScattering* data);
|
||||
bool sab_has_nuclide(ThermalScattering* data, const char* name);
|
||||
void sab_sample(ThermalScattering* data, const NuclideMicroXS* micro_xs,
|
||||
double E_in, double* E_out, double* mu);
|
||||
double sab_threshold(ThermalScattering* data);
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
#endif // OPENMC_THERMAL_SCATTERING_H
|
||||
|
|
@ -40,4 +40,19 @@ get_node_value(pugi::xml_node node, const char* name, bool lowercase,
|
|||
return value;
|
||||
}
|
||||
|
||||
bool
|
||||
get_node_value_bool(pugi::xml_node node, const char* name)
|
||||
{
|
||||
if (node.attribute(name)) {
|
||||
return node.attribute(name).as_bool();
|
||||
} else if (node.child(name)) {
|
||||
return node.child(name).text().as_bool();
|
||||
} else {
|
||||
std::stringstream err_msg;
|
||||
err_msg << "Node \"" << name << "\" is not a member of the \""
|
||||
<< node.name() << "\" XML node";
|
||||
fatal_error(err_msg);
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -16,9 +16,11 @@ check_for_node(pugi::xml_node node, const char *name)
|
|||
return node.attribute(name) || node.child(name);
|
||||
}
|
||||
|
||||
std::string get_node_value(pugi::xml_node node, const char *name,
|
||||
std::string get_node_value(pugi::xml_node node, const char* name,
|
||||
bool lowercase=false, bool strip=false);
|
||||
|
||||
bool get_node_value_bool(pugi::xml_node node, const char* name);
|
||||
|
||||
template <typename T>
|
||||
std::vector<T> get_node_array(pugi::xml_node node, const char* name,
|
||||
bool lowercase=false)
|
||||
|
|
|
|||
|
|
@ -7,10 +7,6 @@ set -ex
|
|||
# Upgrade pip before doing anything else
|
||||
pip install --upgrade pip
|
||||
|
||||
# Running OpenMC's setup.py requires numpy/cython already
|
||||
pip install numpy
|
||||
pip install cython
|
||||
|
||||
# pytest installed by default -- make sure we get latest
|
||||
pip install --upgrade pytest
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue