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

resolving @paulromano comments, including intent(inout) items at the end, consistent usage of const qualifiers

Browse source
This commit is contained in:
Adam G Nelson 2018-06-21 20:45:10 -04:00
parent b704480f9f
commit 267acf627f
14 changed files with 421 additions and 497 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/constants.h
View file

@ -21,7 +21,7 @@ typedef std::vector<int> int_1dvec;
typedef std::vector<std::vector<int> > int_2dvec;
typedef std::vector<std::vector<std::vector<int> > > int_3dvec;
int constexpr MAX_SAMPLE {10000};
constexpr int MAX_SAMPLE {10000};
constexpr std::array<int, 3> VERSION {0, 10, 0};
constexpr std::array<int, 2> VERSION_PARTICLE_RESTART {2, 0};

143
src/mgxs.cpp
View file

@ -12,13 +12,12 @@ std::vector<Mgxs> macro_xs;
//==============================================================================
void
Mgxs::init(const std::string& in_name, const double in_awr,
const double_1dvec& in_kTs, const bool in_fissionable,
const int in_scatter_format, const int in_num_groups,
const int in_num_delayed_groups, const bool in_is_isotropic,
const double_1dvec& in_polar, const double_1dvec& in_azimuthal,
const int n_threads)
Mgxs::init(const std::string& in_name, double in_awr,
const double_1dvec& in_kTs, bool in_fissionable, int in_scatter_format,
int in_num_groups, int in_num_delayed_groups, bool in_is_isotropic,
const double_1dvec& in_polar, const double_1dvec& in_azimuthal)
{
// Set the metadata
name = in_name;
awr = in_awr;
kTs = in_kTs;
@ -32,6 +31,13 @@ Mgxs::init(const std::string& in_name, const double in_awr,
n_azi = in_azimuthal.size();
polar = in_polar;
azimuthal = in_azimuthal;
// Set the cross section index cache
#ifdef _OPENMP
int n_threads = omp_get_max_threads();
#else
int n_threads = 1;
#endif
cache.resize(n_threads);
for (int thread = 0; thread < n_threads; thread++) {
cache[thread].sqrtkT = 0.;
@ -46,10 +52,9 @@ Mgxs::init(const std::string& in_name, const double in_awr,
//==============================================================================
void
Mgxs::metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
const int in_num_delayed_groups, double_1dvec& temperature, int& method,
const double tolerance, int_1dvec& temps_to_read, int& order_dim,
const int n_threads)
Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups,
int in_num_delayed_groups, const double_1dvec& temperature,
double tolerance, int_1dvec& temps_to_read, int& order_dim, int& method)
{
// get name
char char_name[MAX_WORD_LEN];
@ -252,23 +257,20 @@ Mgxs::metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
// Finally use this data to initialize the MGXS Object
init(in_name, in_awr, in_kTs, in_fissionable, in_scatter_format,
in_num_groups, in_num_delayed_groups, in_is_isotropic, in_polar,
in_azimuthal, n_threads);
in_azimuthal);
}
//==============================================================================
void
Mgxs::from_hdf5(hid_t xs_id, const int energy_groups,
const int delayed_groups, double_1dvec& temperature, int& method,
const double tolerance, const int max_order,
const bool legendre_to_tabular, const int legendre_to_tabular_points,
const int n_threads)
Mgxs::Mgxs(hid_t xs_id, int energy_groups, int delayed_groups,
const double_1dvec& temperature, double tolerance, int max_order,
bool legendre_to_tabular, int legendre_to_tabular_points, int& method)
{
// Call generic data gathering routine (will populate the metadata)
int order_data;
int_1dvec temps_to_read;
metadata_from_hdf5(xs_id, energy_groups, delayed_groups, temperature,
method, tolerance, temps_to_read, order_data, n_threads);
tolerance, temps_to_read, order_data, method);
// Set number of energy and delayed groups
int final_scatter_format = scatter_format;
@ -297,10 +299,9 @@ Mgxs::from_hdf5(hid_t xs_id, const int energy_groups,
//==============================================================================
void
Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs,
std::vector<Mgxs*>& micros, double_1dvec& atom_densities, int& method,
const double tolerance, const int n_threads)
Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
const std::vector<Mgxs*>& micros, const double_1dvec& atom_densities,
double tolerance, int& method)
{
// Get the minimum data needed to initialize:
// Dont need awr, but lets just initialize it anyways
@ -321,7 +322,7 @@ Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs,
init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format,
in_num_groups, in_num_delayed_groups, in_is_isotropic, in_polar,
in_azimuthal, n_threads);
in_azimuthal);
// Create the xs data for each temperature
for (int t = 0; t < mat_kTs.size(); t++) {
@ -397,8 +398,8 @@ Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs,
//==============================================================================
void
Mgxs::combine(std::vector<Mgxs*>& micros, double_1dvec& scalars,
int_1dvec& micro_ts, int this_t)
Mgxs::combine(const std::vector<Mgxs*>& micros, const double_1dvec& scalars,
const int_1dvec& micro_ts, int this_t)
{
// Build the vector of pointers to the xs objects within micros
std::vector<XsData*> those_xs(micros.size());
@ -415,36 +416,42 @@ Mgxs::combine(std::vector<Mgxs*>& micros, double_1dvec& scalars,
//==============================================================================
double
Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
double* mu, int* dg)
Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
{
// This method assumes that the temperature and angle indices are set
#ifdef _OPENMP
int tid = omp_get_thread_num();
XsData* xs_t = &xs[cache[tid].t];
int a = cache[tid].a;
#else
XsData* xs_t = &xs[cache[0].t];
int a = cache[0].a;
#endif
double val;
switch(xstype) {
case MG_GET_XS_TOTAL:
val = xs_t->total[cache[tid].a][gin];
val = xs_t->total[a][gin];
break;
case MG_GET_XS_NU_FISSION:
if (fissionable) {
val = xs_t->nu_fission[cache[tid].a][gin];
val = xs_t->nu_fission[a][gin];
} else {
val = 0.;
}
break;
case MG_GET_XS_ABSORPTION:
val = xs_t->absorption[cache[tid].a][gin];
val = xs_t->absorption[a][gin];
break;
case MG_GET_XS_FISSION:
if (fissionable) {
val = xs_t->fission[cache[tid].a][gin];
val = xs_t->fission[a][gin];
} else {
val = 0.;
}
break;
case MG_GET_XS_KAPPA_FISSION:
if (fissionable) {
val = xs_t->kappa_fission[cache[tid].a][gin];
val = xs_t->kappa_fission[a][gin];
} else {
val = 0.;
}
@ -453,11 +460,11 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
case MG_GET_XS_SCATTER_MULT:
case MG_GET_XS_SCATTER_FMU_MULT:
case MG_GET_XS_SCATTER_FMU:
val = xs_t->scatter[cache[tid].a]->get_xs(xstype, gin, gout, mu);
val = xs_t->scatter[a]->get_xs(xstype, gin, gout, mu);
break;
case MG_GET_XS_PROMPT_NU_FISSION:
if (fissionable) {
val = xs_t->prompt_nu_fission[cache[tid].a][gin];
val = xs_t->prompt_nu_fission[a][gin];
} else {
val = 0.;
}
@ -465,10 +472,10 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
case MG_GET_XS_DELAYED_NU_FISSION:
if (fissionable) {
if (dg != nullptr) {
val = xs_t->delayed_nu_fission[cache[tid].a][gin][*dg];
val = xs_t->delayed_nu_fission[a][gin][*dg];
} else {
val = 0.;
for (auto& num : xs_t->delayed_nu_fission[cache[tid].a][gin]) {
for (auto& num : xs_t->delayed_nu_fission[a][gin]) {
val += num;
}
}
@ -479,11 +486,11 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
case MG_GET_XS_CHI_PROMPT:
if (fissionable) {
if (gout != nullptr) {
val = xs_t->chi_prompt[cache[tid].a][gin][*gout];
val = xs_t->chi_prompt[a][gin][*gout];
} else {
// provide an outgoing group-wise sum
val = 0.;
for (auto& num : xs_t->chi_prompt[cache[tid].a][gin]) {
for (auto& num : xs_t->chi_prompt[a][gin]) {
val += num;
}
}
@ -495,20 +502,20 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
if (fissionable) {
if (gout != nullptr) {
if (dg != nullptr) {
val = xs_t->chi_delayed[cache[tid].a][gin][*gout][*dg];
val = xs_t->chi_delayed[a][gin][*gout][*dg];
} else {
val = xs_t->chi_delayed[cache[tid].a][gin][*gout][0];
val = xs_t->chi_delayed[a][gin][*gout][0];
}
} else {
if (dg != nullptr) {
val = 0.;
for (int i = 0; i < xs_t->chi_delayed[cache[tid].a][gin].size(); i++) {
val += xs_t->chi_delayed[cache[tid].a][gin][i][*dg];
for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) {
val += xs_t->chi_delayed[a][gin][i][*dg];
}
} else {
val = 0.;
for (int i = 0; i < xs_t->chi_delayed[cache[tid].a][gin].size(); i++) {
for (auto& num : xs_t->chi_delayed[cache[tid].a][gin][i]) {
for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) {
for (auto& num : xs_t->chi_delayed[a][gin][i]) {
val += num;
}
}
@ -519,13 +526,13 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
}
break;
case MG_GET_XS_INVERSE_VELOCITY:
val = xs_t->inverse_velocity[cache[tid].a][gin];
val = xs_t->inverse_velocity[a][gin];
break;
case MG_GET_XS_DECAY_RATE:
if (dg != nullptr) {
val = xs_t->decay_rate[cache[tid].a][*dg + 1];
val = xs_t->decay_rate[a][*dg + 1];
} else {
val = xs_t->decay_rate[cache[tid].a][0];
val = xs_t->decay_rate[a][0];
}
break;
default:
@ -537,9 +544,14 @@ Mgxs::get_xs(const int tid, const int xstype, const int gin, int* gout,
//==============================================================================
void
Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, int& gout)
Mgxs::sample_fission_energy(int gin, int& dg, int& gout)
{
// This method assumes that the temperature and angle indices are set
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
XsData* xs_t = &xs[cache[tid].t];
double nu_fission = xs_t->nu_fission[cache[tid].a][gin];
@ -596,23 +608,32 @@ Mgxs::sample_fission_energy(const int tid, const int gin, int& dg, int& gout)
//==============================================================================
void
Mgxs::sample_scatter(const int tid, const int gin, int& gout, double& mu,
double& wgt)
Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt)
{
// This method assumes that the temperature and angle indices are set
// Sample the data
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt);
}
//==============================================================================
void
Mgxs::calculate_xs(const int tid, const int gin, const double sqrtkT,
const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs)
Mgxs::calculate_xs(int gin, double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs)
{
// Set our indices
set_temperature_index(tid, sqrtkT);
set_angle_index(tid, uvw);
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
set_temperature_index(sqrtkT);
set_angle_index(uvw);
XsData* xs_t = &xs[cache[tid].t];
total_xs = xs_t->total[cache[tid].a][gin];
abs_xs = xs_t->absorption[cache[tid].a][gin];
@ -646,9 +667,14 @@ Mgxs::equiv(const Mgxs& that)
//==============================================================================
void
Mgxs::set_temperature_index(const int tid, const double sqrtkT)
Mgxs::set_temperature_index(double sqrtkT)
{
// See if we need to find the new index
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
if (sqrtkT != cache[tid].sqrtkT) {
double kT = sqrtkT * sqrtkT;
@ -668,9 +694,14 @@ Mgxs::set_temperature_index(const int tid, const double sqrtkT)
//==============================================================================
void
Mgxs::set_angle_index(const int tid, const double uvw[3])
Mgxs::set_angle_index(const double uvw[3])
{
// See if we need to find the new index
#ifdef _OPENMP
int tid = omp_get_thread_num();
#else
int tid = 0;
#endif
if (!is_isotropic &&
((uvw[0] != cache[tid].u) || (uvw[1] != cache[tid].v) ||
(uvw[2] != cache[tid].w))) {

119
src/mgxs.h
View file

@ -11,6 +11,10 @@
#include <valarray>
#include <vector>
#ifdef _OPENMP
# include <omp.h>
#endif
#include "constants.h"
#include "hdf5_interface.h"
#include "math_functions.h"
@ -55,24 +59,6 @@ class Mgxs {
double_1dvec polar;
double_1dvec azimuthal;
//! \brief Initializes the Mgxs object metadata from the HDF5 file
//!
//! @param xs_id HDF5 group id for the cross section data.
//! @param in_num_groups Number of energy groups.
//! @param in_num_delayed_groups Number of delayed groups.
//! @param temperature Temperatures to read.
//! @param method Method of choosing nearest temperatures.
//! @param tolerance Tolerance of temperature selection method.
//! @param temps_to_read Resultant list of temperatures in the library
//! to read which correspond to the requested temperatures.
//! @param order_dim Resultant dimensionality of the scattering order.
//! @param n_threads Number of threads at runtime.
void
metadata_from_hdf5(const hid_t xs_id, const int in_num_groups,
const int in_num_delayed_groups, double_1dvec& temperature,
int& method, const double tolerance, int_1dvec& temps_to_read,
int& order_dim, const int n_threads);
//! \brief Initializes the Mgxs object metadata
//!
//! @param in_name Name of the object.
@ -87,14 +73,28 @@ class Mgxs {
//! the incoming particle.
//! @param in_polar Polar angle grid.
//! @param in_azimuthal Azimuthal angle grid.
//! @param n_threads Number of threads at runtime.
void
init(const std::string& in_name, const double in_awr,
const double_1dvec& in_kTs, const bool in_fissionable,
const int in_scatter_format, const int in_num_groups,
const int in_num_delayed_groups, const bool in_is_isotropic,
const double_1dvec& in_polar, const double_1dvec& in_azimuthal,
const int n_threads);
init(const std::string& in_name, double in_awr, const double_1dvec& in_kTs,
bool in_fissionable, int in_scatter_format, int in_num_groups,
int in_num_delayed_groups, bool in_is_isotropic,
const double_1dvec& in_polar, const double_1dvec& in_azimuthal);
//! \brief Initializes the Mgxs object metadata from the HDF5 file
//!
//! @param xs_id HDF5 group id for the cross section data.
//! @param in_num_groups Number of energy groups.
//! @param in_num_delayed_groups Number of delayed groups.
//! @param temperature Temperatures to read.
//! @param tolerance Tolerance of temperature selection method.
//! @param temps_to_read Resultant list of temperatures in the library
//! to read which correspond to the requested temperatures.
//! @param order_dim Resultant dimensionality of the scattering order.
//! @param method Method of choosing nearest temperatures.
void
metadata_from_hdf5(hid_t xs_id, int in_num_groups,
int in_num_delayed_groups, const double_1dvec& temperature,
double tolerance, int_1dvec& temps_to_read, int& order_dim,
int& method);
//! \brief Performs the actual act of combining the microscopic data for a
//! single temperature.
@ -105,8 +105,8 @@ class Mgxs {
//! corresponds to the temperature of interest.
//! @param this_t The temperature index of the macroscopic object.
void
combine(std::vector<Mgxs*>& micros, double_1dvec& scalars,
int_1dvec& micro_ts, int this_t);
combine(const std::vector<Mgxs*>& micros, const double_1dvec& scalars,
const int_1dvec& micro_ts, int this_t);
//! \brief Checks to see if this and that are able to be combined
//!
@ -123,28 +123,14 @@ class Mgxs {
bool fissionable; // Is this fissionable
std::vector<CacheData> cache; // index and data cache
//! \brief Initializes and populates all data to build a macroscopic
//! cross section from microscopic cross section.
//!
//! @param in_name Name of the object.
//! @param mat_kTs temperatures (in units of eV) that data is needed.
//! @param micros Microscopic objects to combine.
//! @param atom_densities Atom densities of those microscopic quantities.
//! @param method Method of choosing nearest temperatures.
//! @param tolerance Tolerance of temperature selection method.
//! @param n_threads Number of threads at runtime.
void
build_macro(const std::string& in_name, double_1dvec& mat_kTs,
std::vector<Mgxs*>& micros, double_1dvec& atom_densities,
int& method, const double tolerance, const int n_threads);
Mgxs() = default;
//! \brief Loads the Mgxs object from the HDF5 file
//! \brief Constructor that loads the Mgxs object from the HDF5 file
//!
//! @param xs_id HDF5 group id for the cross section data.
//! @param energy_groups Number of energy groups.
//! @param delayed_groups Number of delayed groups.
//! @param temperature Temperatures to read.
//! @param method Method of choosing nearest temperatures.
//! @param tolerance Tolerance of temperature selection method.
//! @param max_order Maximum order requested by the user;
//! this is only used for Legendre scattering.
@ -152,13 +138,24 @@ class Mgxs {
//! should be converted to a Tabular representation.
//! @param legendre_to_tabular_points If a conversion is requested, this
//! provides the number of points to use in the tabular representation.
//! @param n_threads Number of threads at runtime.
void
from_hdf5(hid_t xs_id, const int energy_groups,
const int delayed_groups, double_1dvec& temperature, int& method,
const double tolerance, const int max_order,
const bool legendre_to_tabular, const int legendre_to_tabular_points,
const int n_threads);
//! @param method Method of choosing nearest temperatures.
Mgxs(hid_t xs_id, int energy_groups,
int delayed_groups, const double_1dvec& temperature, double tolerance,
int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points, int& method);
//! \brief Constructor that initializes and populates all data to build a
//! macroscopic cross section from microscopic cross section.
//!
//! @param in_name Name of the object.
//! @param mat_kTs temperatures (in units of eV) that data is needed.
//! @param micros Microscopic objects to combine.
//! @param atom_densities Atom densities of those microscopic quantities.
//! @param tolerance Tolerance of temperature selection method.
//! @param method Method of choosing nearest temperatures.
Mgxs(const std::string& in_name, const double_1dvec& mat_kTs,
const std::vector<Mgxs*>& micros, const double_1dvec& atom_densities,
double tolerance, int& method);
//! \brief Provides a cross section value given certain parameters
//!
@ -172,32 +169,27 @@ class Mgxs {
//! @param dg delayed group index; use nullptr if irrelevant.
//! @return Requested cross section value.
double
get_xs(const int tid, const int xstype, const int gin, int* gout,
double* mu, int* dg);
get_xs(int xstype, int gin, int* gout, double* mu, int* dg);
//! \brief Samples the fission neutron energy and if prompt or delayed.
//!
//! @param tid Thread id to use when using the index cache.
//! @param gin Incoming energy group.
//! @param dg Sampled delayed group index.
//! @param gout Sampled outgoing energy group.
void
sample_fission_energy(const int tid, const int gin, int& dg, int& gout);
sample_fission_energy(int gin, int& dg, int& gout);
//! \brief Samples the outgoing energy and angle from a scatter event.
//!
//! @param tid Thread id to use when using the index cache.
//! @param gin Incoming energy group.
//! @param gout Sampled outgoing energy group.
//! @param mu Sampled cosine of the change-in-angle.
//! @param wgt Weight of the particle to be adjusted.
void
sample_scatter(const int tid, const int gin, int& gout, double& mu,
double& wgt);
sample_scatter(int gin, int& gout, double& mu, double& wgt);
//! \brief Calculates cross section quantities needed for tracking.
//!
//! @param tid Thread id to use when using the index cache.
//! @param gin Incoming energy group.
//! @param sqrtkT Temperature of the material.
//! @param uvw Incoming particle direction.
@ -205,23 +197,20 @@ class Mgxs {
//! @param abs_xs Resultant absorption cross section.
//! @param nu_fiss_xs Resultant nu-fission cross section.
void
calculate_xs(const int tid, const int gin, const double sqrtkT,
const double uvw[3], double& total_xs, double& abs_xs,
double& nu_fiss_xs);
calculate_xs(int gin, double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs);
//! \brief Sets the temperature index in cache given a temperature
//!
//! @param tid Thread id to use when setting the index cache.
//! @param sqrtkT Temperature of the material.
void
set_temperature_index(const int tid, const double sqrtkT);
set_temperature_index(double sqrtkT);
//! \brief Sets the angle index in cache given a direction
//!
//! @param tid Thread id to use when setting the index cache.
//! @param uvw Incoming particle direction.
void
set_angle_index(const int tid, const double uvw[3]);
set_angle_index(const double uvw[3]);
};
} // namespace openmc

27
src/mgxs_data.F90
View file

@ -2,10 +2,6 @@ module mgxs_data
use, intrinsic :: ISO_C_BINDING
#ifdef _OPENMP
use omp_lib
#endif
use constants
use algorithm, only: find
use dict_header, only: DictCharInt
@ -40,13 +36,6 @@ contains
type(VectorReal), allocatable, target :: temps(:)
character(MAX_WORD_LEN) :: word
integer, allocatable :: array(:)
integer(C_INT) :: n_threads
#ifdef _OPENMP
n_threads = OMP_GET_MAX_THREADS()
#else
n_threads = 1
#endif
! Check if MGXS Library exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
@ -91,12 +80,11 @@ contains
i_nuclide = mat % nuclide(j)
if (.not. already_read % contains(name)) then
call add_mgxs_c(file_id, name, &
num_energy_groups, num_delayed_groups, &
call add_mgxs_c(file_id, name, num_energy_groups, num_delayed_groups, &
temps(i_nuclide) % size(), temps(i_nuclide) % data, &
temperature_method, temperature_tolerance, max_order, &
temperature_tolerance, max_order, &
logical(legendre_to_tabular, C_BOOL), &
legendre_to_tabular_points, n_threads)
legendre_to_tabular_points, temperature_method)
call already_read % add(name)
end if
@ -122,13 +110,6 @@ contains
type(Material), pointer :: mat ! current material
type(VectorReal), allocatable :: kTs(:)
character(MAX_WORD_LEN) :: name ! name of material
integer(C_INT) :: n_threads
#ifdef _OPENMP
n_threads = OMP_GET_MAX_THREADS()
#else
n_threads = 1
#endif
! Get temperatures to read for each material
call get_mat_kTs(kTs)
@ -149,7 +130,7 @@ contains
if (allocated(kTs(i_mat) % data)) then
call create_macro_xs_c(name, mat % n_nuclides, mat % nuclide, &
kTs(i_mat) % size(), kTs(i_mat) % data, mat % atom_density, &
temperature_method, temperature_tolerance, n_threads)
temperature_tolerance, temperature_method)
end if
end do

36
src/mgxs_interface.F90
View file

@ -9,8 +9,8 @@ module mgxs_interface
interface
subroutine add_mgxs_c(file_id, name, energy_groups, delayed_groups, &
n_temps, temps, method, tolerance, max_order, legendre_to_tabular, &
legendre_to_tabular_points, n_threads) bind(C)
n_temps, temps, tolerance, max_order, legendre_to_tabular, &
legendre_to_tabular_points, method) bind(C)
use ISO_C_BINDING
import HID_T
implicit none
@ -20,12 +20,11 @@ module mgxs_interface
integer(C_INT), value, intent(in) :: delayed_groups
integer(C_INT), value, intent(in) :: n_temps
real(C_DOUBLE), intent(in) :: temps(1:n_temps)
integer(C_INT), intent(inout) :: method
real(C_DOUBLE), value, intent(in) :: tolerance
integer(C_INT), value, intent(in) :: max_order
logical(C_BOOL),value, intent(in) :: legendre_to_tabular
integer(C_INT), value, intent(in) :: legendre_to_tabular_points
integer(C_INT), value, intent(in) :: n_threads
integer(C_INT), intent(inout) :: method
end subroutine add_mgxs_c
function query_fissionable_c(n_nuclides, i_nuclides) result(result) bind(C)
@ -37,7 +36,7 @@ module mgxs_interface
end function query_fissionable_c
subroutine create_macro_xs_c(name, n_nuclides, i_nuclides, n_temps, temps, &
atom_densities, method, tolerance, n_threads) bind(C)
atom_densities, tolerance, method) bind(C)
use ISO_C_BINDING
implicit none
character(kind=C_CHAR),intent(in) :: name(*)
@ -46,17 +45,15 @@ module mgxs_interface
integer(C_INT), value, intent(in) :: n_temps
real(C_DOUBLE), intent(in) :: temps(1:n_temps)
real(C_DOUBLE), intent(in) :: atom_densities(1:n_nuclides)
integer(C_INT), intent(inout) :: method
real(C_DOUBLE), value, intent(in) :: tolerance
integer(C_INT), value, intent(in) :: n_threads
integer(C_INT), intent(inout) :: method
end subroutine create_macro_xs_c
subroutine calculate_xs_c(i_mat, tid, gin, sqrtkT, uvw, total_xs, abs_xs, &
subroutine calculate_xs_c(i_mat, gin, sqrtkT, uvw, total_xs, abs_xs, &
nu_fiss_xs) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: i_mat
integer(C_INT), value, intent(in) :: tid
integer(C_INT), value, intent(in) :: gin
real(C_DOUBLE), value, intent(in) :: sqrtkT
real(C_DOUBLE), intent(in) :: uvw(1:3)
@ -65,11 +62,10 @@ module mgxs_interface
real(C_DOUBLE), intent(inout) :: nu_fiss_xs
end subroutine calculate_xs_c
subroutine sample_scatter_c(i_mat, tid, gin, gout, mu, wgt, uvw) bind(C)
subroutine sample_scatter_c(i_mat, gin, gout, mu, wgt, uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: i_mat
integer(C_INT), value, intent(in) :: tid
integer(C_INT), value, intent(in) :: gin
integer(C_INT), intent(inout) :: gout
real(C_DOUBLE), intent(inout) :: mu
@ -77,11 +73,10 @@ module mgxs_interface
real(C_DOUBLE), intent(inout) :: uvw(1:3)
end subroutine sample_scatter_c
subroutine sample_fission_energy_c(i_mat, tid, gin, dg, gout) bind(C)
subroutine sample_fission_energy_c(i_mat, gin, dg, gout) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: i_mat
integer(C_INT), value, intent(in) :: tid
integer(C_INT), value, intent(in) :: gin
integer(C_INT), intent(inout) :: dg
integer(C_INT), intent(inout) :: gout
@ -102,12 +97,11 @@ module mgxs_interface
real(C_DOUBLE) :: awr
end function get_awr_c
function get_nuclide_xs_c(index, tid, xstype, gin, gout, mu, dg) result(val) &
function get_nuclide_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: tid
integer(C_INT), value, intent(in) :: xstype
integer(C_INT), value, intent(in) :: gin
integer(C_INT), optional, intent(in) :: gout
@ -116,12 +110,11 @@ module mgxs_interface
real(C_DOUBLE) :: val
end function get_nuclide_xs_c
function get_macro_xs_c(index, tid, xstype, gin, gout, mu, dg) result(val) &
function get_macro_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: tid
integer(C_INT), value, intent(in) :: xstype
integer(C_INT), value, intent(in) :: gin
integer(C_INT), optional, intent(in) :: gout
@ -130,27 +123,24 @@ module mgxs_interface
real(C_DOUBLE) :: val
end function get_macro_xs_c
subroutine set_nuclide_angle_index_c(index, tid, uvw) bind(C)
subroutine set_nuclide_angle_index_c(index, uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: tid
real(C_DOUBLE), intent(in) :: uvw(1:3)
end subroutine set_nuclide_angle_index_c
subroutine set_macro_angle_index_c(index, tid, uvw) bind(C)
subroutine set_macro_angle_index_c(index, uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: tid
real(C_DOUBLE), intent(in) :: uvw(1:3)
end subroutine set_macro_angle_index_c
subroutine set_nuclide_temperature_index_c(index, tid, sqrtkT) bind(C)
subroutine set_nuclide_temperature_index_c(index, sqrtkT) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: tid
real(C_DOUBLE), value, intent(in) :: sqrtkT
end subroutine set_nuclide_temperature_index_c

84
src/mgxs_interface.cpp
View file

@ -7,11 +7,10 @@ namespace openmc {
//==============================================================================
void
add_mgxs_c(hid_t file_id, char* name, const int energy_groups,
const int delayed_groups, const int n_temps, double temps[], int& method,
const double tolerance, const int max_order,
const bool legendre_to_tabular, const int legendre_to_tabular_points,
const int n_threads)
add_mgxs_c(hid_t file_id, const char* name, int energy_groups,
int delayed_groups, int n_temps, const double temps[], double tolerance,
int max_order, bool legendre_to_tabular, int legendre_to_tabular_points,
int& method)
{
//!! mgxs_data.F90 will be modified to just create the list of names
//!! in the order needed
@ -30,10 +29,8 @@ add_mgxs_c(hid_t file_id, char* name, const int energy_groups,
+ "provided MGXS Library");
}
Mgxs mg;
mg.from_hdf5(xs_grp, energy_groups, delayed_groups,
temperature, method, tolerance, max_order, legendre_to_tabular,
legendre_to_tabular_points, n_threads);
Mgxs mg(xs_grp, energy_groups, delayed_groups, temperature, tolerance,
max_order, legendre_to_tabular, legendre_to_tabular_points, method);
nuclides_MG.push_back(mg);
}
@ -41,7 +38,7 @@ add_mgxs_c(hid_t file_id, char* name, const int energy_groups,
//==============================================================================
bool
query_fissionable_c(const int n_nuclides, const int i_nuclides[])
query_fissionable_c(int n_nuclides, const int i_nuclides[])
{
bool result = false;
for (int n = 0; n < n_nuclides; n++) {
@ -53,12 +50,10 @@ query_fissionable_c(const int n_nuclides, const int i_nuclides[])
//==============================================================================
void
create_macro_xs_c(char* mat_name, const int n_nuclides,
const int i_nuclides[], const int n_temps, const double temps[],
const double atom_densities[], int& method, const double tolerance,
const int n_threads)
create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
int n_temps, const double temps[], const double atom_densities[],
double tolerance, int& method)
{
Mgxs macro;
if (n_temps > 0) {
// // Convert temps to a vector
double_1dvec temperature;
@ -75,10 +70,14 @@ create_macro_xs_c(char* mat_name, const int n_nuclides,
mgxs_ptr[n] = &nuclides_MG[i_nuclides[n] - 1];
}
macro.build_macro(mat_name, temperature, mgxs_ptr, atom_densities_vec,
method, tolerance, n_threads);
Mgxs macro(mat_name, temperature, mgxs_ptr, atom_densities_vec,
tolerance, method);
macro_xs.push_back(macro);
} else {
// Preserve the ordering of materials by including a blank entry
Mgxs macro;
macro_xs.push_back(macro);
}
macro_xs.push_back(macro);
}
//==============================================================================
@ -86,22 +85,21 @@ create_macro_xs_c(char* mat_name, const int n_nuclides,
//==============================================================================
void
calculate_xs_c(const int i_mat, const int tid, const int gin,
const double sqrtkT, const double uvw[3], double& total_xs, double& abs_xs,
double& nu_fiss_xs)
calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs)
{
macro_xs[i_mat - 1].calculate_xs(tid, gin - 1, sqrtkT, uvw, total_xs, abs_xs,
macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
nu_fiss_xs);
}
//==============================================================================
void
sample_scatter_c(const int i_mat, const int tid, const int gin, int& gout,
double& mu, double& wgt, double uvw[3])
sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
double uvw[3])
{
int gout_c = gout - 1;
macro_xs[i_mat - 1].sample_scatter(tid, gin - 1, gout_c, mu, wgt);
macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt);
// adjust return value for fortran indexing
gout = gout_c + 1;
@ -113,12 +111,11 @@ sample_scatter_c(const int i_mat, const int tid, const int gin, int& gout,
//==============================================================================
void
sample_fission_energy_c(const int i_mat, const int tid, const int gin,
int& dg, int& gout)
sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout)
{
int dg_c = 0;
int gout_c = 0;
macro_xs[i_mat - 1].sample_fission_energy(tid, gin - 1, dg_c, gout_c);
macro_xs[i_mat - 1].sample_fission_energy(gin - 1, dg_c, gout_c);
// adjust return values for fortran indexing
dg = dg_c + 1;
@ -128,8 +125,7 @@ sample_fission_energy_c(const int i_mat, const int tid, const int gin,
//==============================================================================
double
get_nuclide_xs_c(const int index, const int tid, const int xstype,
const int gin, int* gout, double* mu, int* dg)
get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
{
int gout_c;
int* gout_c_p;
@ -147,15 +143,13 @@ get_nuclide_xs_c(const int index, const int tid, const int xstype,
} else {
dg_c_p = dg;
}
return nuclides_MG[index - 1].get_xs(tid, xstype, gin - 1, gout_c_p, mu,
dg_c_p);
return nuclides_MG[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
}
//==============================================================================
double
get_macro_xs_c(const int index, const int tid, const int xstype,
const int gin, int* gout, double* mu, int* dg)
get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
{
int gout_c;
int* gout_c_p;
@ -173,38 +167,34 @@ get_macro_xs_c(const int index, const int tid, const int xstype,
} else {
dg_c_p = dg;
}
return macro_xs[index - 1].get_xs(tid, xstype, gin - 1, gout_c_p, mu,
dg_c_p);
return macro_xs[index - 1].get_xs(xstype, gin - 1, gout_c_p, mu, dg_c_p);
}
//==============================================================================
void
set_nuclide_angle_index_c(const int index, const int tid,
const double uvw[3])
set_nuclide_angle_index_c(int index, const double uvw[3])
{
// Update the values
nuclides_MG[index - 1].set_angle_index(tid, uvw);
nuclides_MG[index - 1].set_angle_index(uvw);
}
//==============================================================================
void
set_macro_angle_index_c(const int index, const int tid,
const double uvw[3])
set_macro_angle_index_c(int index, const double uvw[3])
{
// Update the values
macro_xs[index - 1].set_angle_index(tid, uvw);
macro_xs[index - 1].set_angle_index(uvw);
}
//==============================================================================
void
set_nuclide_temperature_index_c(const int index, const int tid,
const double sqrtkT)
set_nuclide_temperature_index_c(int index, double sqrtkT)
{
// Update the values
nuclides_MG[index - 1].set_temperature_index(tid, sqrtkT);
nuclides_MG[index - 1].set_temperature_index(sqrtkT);
}
//==============================================================================
@ -212,7 +202,7 @@ set_nuclide_temperature_index_c(const int index, const int tid,
//==============================================================================
void
get_name_c(const int index, int name_len, char* name)
get_name_c(int index, int name_len, char* name)
{
// First blank out our input string
std::string str(name_len, ' ');
@ -229,7 +219,7 @@ get_name_c(const int index, int name_len, char* name)
//==============================================================================
double
get_awr_c(const int index)
get_awr_c(int index)
{
return nuclides_MG[index - 1].awr;
}

49
src/mgxs_interface.h
View file

@ -22,67 +22,58 @@ extern std::vector<Mgxs> macro_xs;
//==============================================================================
extern "C" void
add_mgxs_c(hid_t file_id, char* name, const int energy_groups,
const int delayed_groups, const int n_temps, double temps[], int& method,
const double tolerance, const int max_order,
const bool legendre_to_tabular, const int legendre_to_tabular_points,
const int n_threads);
add_mgxs_c(hid_t file_id, const char* name, int energy_groups,
int delayed_groups, int n_temps, const double temps[], double tolerance,
int max_order, bool legendre_to_tabular, int legendre_to_tabular_points,
int& method);
extern "C" bool
query_fissionable_c(const int n_nuclides, const int i_nuclides[]);
query_fissionable_c(int n_nuclides, const int i_nuclides[]);
extern "C" void
create_macro_xs_c(char* mat_name, const int n_nuclides,
const int i_nuclides[], const int n_temps, const double temps[],
const double atom_densities[], int& method, const double tolerance,
const int n_threads);
create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[],
int n_temps, const double temps[], const double atom_densities[],
double tolerance, int& method);
//==============================================================================
// Mgxs tracking/transport/tallying interface methods
//==============================================================================
extern "C" void
calculate_xs_c(const int i_mat, const int tid, const int gin,
const double sqrtkT, const double uvw[3], double& total_xs,
double& abs_xs, double& nu_fiss_xs);
calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs);
extern "C" void
sample_scatter_c(const int i_mat, const int tid, const int gin,
int& gout, double& mu, double& wgt, double uvw[3]);
sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
double uvw[3]);
extern "C" void
sample_fission_energy_c(const int i_mat, const int tid,
const int gin, int& dg, int& gout);
sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout);
extern "C" double
get_nuclide_xs_c(const int index, const int tid,
const int xstype, const int gin, int* gout, double* mu, int* dg);
get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
extern "C" double
get_macro_xs_c(const int index, const int tid,
const int xstype, const int gin, int* gout, double* mu, int* dg);
get_macro_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
extern "C" void
set_nuclide_angle_index_c(const int index, const int tid,
const double uvw[3]);
set_nuclide_angle_index_c(int index, const double uvw[3]);
extern "C" void
set_macro_angle_index_c(const int index, const int tid,
const double uvw[3]);
set_macro_angle_index_c(int index, const double uvw[3]);
extern "C" void
set_nuclide_temperature_index_c(const int index, const int tid,
const double sqrtkT);
set_nuclide_temperature_index_c(int index, double sqrtkT);
//==============================================================================
// General Mgxs methods
//==============================================================================
extern "C" void
get_name_c(const int index, int name_len, char* name);
get_name_c(int index, int name_len, char* name);
extern "C" double
get_awr_c(const int index);
get_awr_c(int index);
} // namespace openmc
#endif // MGXS_INTERFACE_H

20
src/physics_mg.F90
View file

@ -2,10 +2,6 @@ module physics_mg
! This module contains the multi-group specific physics routines so as to not
! hinder performance of the CE versions with multiple if-thens.
#ifdef _OPENMP
use omp_lib
#endif
use bank_header
use constants
use error, only: fatal_error, warning, write_message
@ -145,14 +141,8 @@ contains
subroutine scatter(p)
type(Particle), intent(inout) :: p
integer(C_INT) :: tid
#ifdef _OPENMP
tid = OMP_GET_THREAD_NUM()
#else
tid = 0
#endif
call sample_scatter_c(p % material, tid, p % last_g, p % g, p % mu, &
call sample_scatter_c(p % material, p % last_g, p % g, p % mu, &
p % wgt, p % coord(1) % uvw)
! Update energy value for downstream compatability (in tallying)
@ -183,12 +173,6 @@ contains
real(8) :: mu ! fission neutron angular cosine
real(8) :: phi ! fission neutron azimuthal angle
real(8) :: weight ! weight adjustment for ufs method
integer(C_INT) :: tid
#ifdef _OPENMP
tid = OMP_GET_THREAD_NUM()
#else
tid = 0
#endif
! TODO: Heat generation from fission
@ -268,7 +252,7 @@ contains
! Sample secondary energy distribution for fission reaction and set energy
! in fission bank
call sample_fission_energy_c(p % material, tid, p % g, dg, gout)
call sample_fission_energy_c(p % material, p % g, dg, gout)
bank_array(i) % E = real(gout, 8)
bank_array(i) % delayed_group = dg

45
src/scattdata.cpp
View file

@ -7,8 +7,9 @@ namespace openmc {
//==============================================================================
void
ScattData::base_init(int order, int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& in_energy, double_2dvec& in_mult)
ScattData::base_init(int order, const int_1dvec& in_gmin,
const int_1dvec& in_gmax, const double_2dvec& in_energy,
const double_2dvec& in_mult)
{
int groups = in_energy.size();
@ -42,7 +43,7 @@ ScattData::base_init(int order, int_1dvec& in_gmin, int_1dvec& in_gmax,
//==============================================================================
void
ScattData::base_combine(const int max_order,
ScattData::base_combine(int max_order,
const std::vector<ScattData*>& those_scatts, const double_1dvec& scalars,
int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& sparse_mult,
double_3dvec& sparse_scatter)
@ -179,8 +180,7 @@ ScattData::sample_energy(int gin, int& gout, int& i_gout)
//==============================================================================
double
ScattData::get_xs(const int xstype, const int gin, const int* gout,
const double* mu)
ScattData::get_xs(int xstype, int gin, const int* gout, const double* mu)
{
// Set the outgoing group offset index as needed
int i_gout = 0;
@ -239,8 +239,8 @@ ScattData::get_xs(const int xstype, const int gin, const int* gout,
//==============================================================================
void
ScattDataLegendre::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& in_mult, double_3dvec& coeffs)
ScattDataLegendre::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
@ -334,7 +334,7 @@ ScattDataLegendre::update_max_val()
//==============================================================================
double
ScattDataLegendre::calc_f(const int gin, const int gout, const double mu)
ScattDataLegendre::calc_f(int gin, int gout, double mu)
{
double f;
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
@ -350,7 +350,7 @@ ScattDataLegendre::calc_f(const int gin, const int gout, const double mu)
//==============================================================================
void
ScattDataLegendre::sample(const int gin, int& gout, double& mu, double& wgt)
ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt)
{
// Sample the outgoing energy using the base-class method
int i_gout;
@ -408,7 +408,7 @@ ScattDataLegendre::combine(const std::vector<ScattData*>& those_scatts,
// so we use a base class method to sum up xs and create new energy and mult
// matrices
ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax,
sparse_mult, sparse_scatter);
sparse_mult, sparse_scatter);
// Got everything we need, store it.
init(in_gmin, in_gmax, sparse_mult, sparse_scatter);
@ -417,7 +417,7 @@ ScattDataLegendre::combine(const std::vector<ScattData*>& those_scatts,
//==============================================================================
double_3dvec
ScattDataLegendre::get_matrix(const int max_order)
ScattDataLegendre::get_matrix(int max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();
@ -442,8 +442,8 @@ ScattDataLegendre::get_matrix(const int max_order)
//==============================================================================
void
ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& in_mult, double_3dvec& coeffs)
ScattDataHistogram::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
@ -479,8 +479,7 @@ ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
}
// Initialize the base class attributes
ScattData::base_init(order, in_gmin, in_gmax, in_energy,
in_mult);
ScattData::base_init(order, in_gmin, in_gmax, in_energy, in_mult);
// Build the angular distribution mu values
mu = double_1dvec(order);
@ -522,7 +521,7 @@ ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
//==============================================================================
double
ScattDataHistogram::calc_f(const int gin, const int gout, const double mu)
ScattDataHistogram::calc_f(int gin, int gout, double mu)
{
double f;
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
@ -546,7 +545,7 @@ ScattDataHistogram::calc_f(const int gin, const int gout, const double mu)
//==============================================================================
void
ScattDataHistogram::sample(const int gin, int& gout, double& mu, double& wgt)
ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt)
{
// Sample the outgoing energy using the base-class method
int i_gout;
@ -581,7 +580,7 @@ ScattDataHistogram::sample(const int gin, int& gout, double& mu, double& wgt)
//==============================================================================
double_3dvec
ScattDataHistogram::get_matrix(const int max_order)
ScattDataHistogram::get_matrix(int max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();
@ -643,8 +642,8 @@ ScattDataHistogram::combine(const std::vector<ScattData*>& those_scatts,
//==============================================================================
void
ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& in_mult, double_3dvec& coeffs)
ScattDataTabular::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs)
{
int groups = coeffs.size();
int order = coeffs[0][0].size();
@ -732,7 +731,7 @@ ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax,
//==============================================================================
double
ScattDataTabular::calc_f(const int gin, const int gout, const double mu)
ScattDataTabular::calc_f(int gin, int gout, double mu)
{
double f;
if ((gout < gmin[gin]) || (gout > gmax[gin])) {
@ -757,7 +756,7 @@ ScattDataTabular::calc_f(const int gin, const int gout, const double mu)
//==============================================================================
void
ScattDataTabular::sample(const int gin, int& gout, double& mu, double& wgt)
ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt)
{
// Sample the outgoing energy using the base-class method
int i_gout;
@ -805,7 +804,7 @@ ScattDataTabular::sample(const int gin, int& gout, double& mu, double& wgt)
//==============================================================================
double_3dvec
ScattDataTabular::get_matrix(const int max_order)
ScattDataTabular::get_matrix(int max_order)
{
// Get the sizes and initialize the data to 0
int groups = energy.size();

129
src/scattdata.h
View file

@ -29,18 +29,17 @@ class ScattData {
protected:
//! \brief Initializes the attributes of the base class.
void
base_init(int order, int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& in_energy, double_2dvec& in_mult);
base_init(int order, const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_energy, const double_2dvec& in_mult);
//! \brief Combines microscopic ScattDatas into a macroscopic one.
void
base_combine(const int max_order,
const std::vector<ScattData*>& those_scatts,
base_combine(int max_order, const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars, int_1dvec& in_gmin, int_1dvec& in_gmax,
double_2dvec& sparse_mult, double_3dvec& sparse_scatter);
public:
double_2dvec energy; // Normalized p0 matrix for sampling Eout
double_2dvec mult; // nu-scatter multiplication (nu-scatt/scatt)
double_3dvec dist; // Angular distribution
@ -51,15 +50,15 @@ class ScattData {
//! \brief Calculates the value of normalized f(mu).
//!
//! The value of f(mu) is normalized as in the integral of f(mu)dmu across
//! [-1,1] is 1.
//! [-1,1] is 1.
//!
//! @param gin Incoming energy group of interest.
//! @param gout Outgoing energy group of interest.
//! @param mu Cosine of the change-in-angle of interest.
//! @return The value of f(mu).
virtual double
calc_f(const int gin, const int gout, const double mu) = 0;
calc_f(int gin, int gout, double mu) = 0;
//! \brief Samples the outgoing energy and angle from the ScattData info.
//!
//! @param gin Incoming energy group.
@ -67,8 +66,8 @@ class ScattData {
//! @param mu Sampled cosine of the change-in-angle.
//! @param wgt Weight of the particle to be adjusted.
virtual void
sample(const int gin, int& gout, double& mu, double& wgt) = 0;
sample(int gin, int& gout, double& mu, double& wgt) = 0;
//! \brief Initializes the ScattData object from a given scatter and
//! multiplicity matrix.
//!
@ -77,9 +76,9 @@ class ScattData {
//! @param in_mult Input sparse multiplicity matrix
//! @param coeffs Input sparse scattering matrix
virtual void
init(int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& in_mult,
double_3dvec& coeffs) = 0;
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs) = 0;
//! \brief Combines the microscopic data.
//!
//! @param those_scatts Microscopic objects to combine.
@ -87,7 +86,7 @@ class ScattData {
virtual void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars) = 0;
//! \brief Getter for the dimensionality of the scattering order.
//!
//! If Legendre this is the "n" in "Pn"; for Tabular, this is the number
@ -96,14 +95,14 @@ class ScattData {
//! @return The order.
virtual int
get_order() = 0;
//! \brief Builds a dense scattering matrix from the constituent parts
//!
//! @param max_order If Legendre this is the maximum value of "n" in "Pn"
//! requested; ignored otherwise.
//! @return The dense scattering matrix.
virtual double_3dvec
get_matrix(const int max_order) = 0;
get_matrix(int max_order) = 0;
//! \brief Samples the outgoing energy from the ScattData info.
//!
@ -124,7 +123,7 @@ class ScattData {
//! use nullptr if irrelevant.
//! @return Requested cross section value.
double
get_xs(const int xstype, const int gin, const int* gout, const double* mu);
get_xs(int xstype, int gin, const int* gout, const double* mu);
};
//==============================================================================
@ -132,44 +131,44 @@ class ScattData {
//==============================================================================
class ScattDataLegendre: public ScattData {
protected:
// Maximal value for rejection sampling from a rectangle
double_2dvec max_val;
// Friend convert_legendre_to_tabular so it has access to protected
// parameters
friend void
convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
int n_mu);
convert_legendre_to_tabular(ScattDataLegendre& leg,
ScattDataTabular& tab, int n_mu);
public:
void
init(int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& in_mult,
double_3dvec& coeffs);
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
const double_1dvec& scalars);
//! \brief Find the maximal value of the angular distribution to use as a
// bounding box with rejection sampling.
void
update_max_val();
double
calc_f(const int gin, const int gout, const double mu);
calc_f(int gin, int gout, double mu);
void
sample(const int gin, int& gout, double& mu, double& wgt);
sample(int gin, int& gout, double& mu, double& wgt);
int
get_order() {return dist[0][0].size() - 1;};
double_3dvec
get_matrix(const int max_order);
get_matrix(int max_order);
};
//==============================================================================
@ -180,32 +179,32 @@ class ScattDataLegendre: public ScattData {
class ScattDataHistogram: public ScattData {
protected:
double_1dvec mu; // Angle distribution mu bin boundaries
double dmu; // Quick storage of the spacing between the mu bin points
double_3dvec fmu; // The angular distribution histogram
public:
void
init(int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& in_mult,
double_3dvec& coeffs);
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
const double_1dvec& scalars);
double
calc_f(const int gin, const int gout, const double mu);
calc_f(int gin, int gout, double mu);
void
sample(const int gin, int& gout, double& mu, double& wgt);
sample(int gin, int& gout, double& mu, double& wgt);
int
get_order() {return dist[0][0].size();};
double_3dvec
get_matrix(const int max_order);
get_matrix(int max_order);
};
//==============================================================================
@ -214,9 +213,9 @@ class ScattDataHistogram: public ScattData {
//==============================================================================
class ScattDataTabular: public ScattData {
protected:
double_1dvec mu; // Angle distribution mu grid points
double dmu; // Quick storage of the spacing between the mu points
double_3dvec fmu; // The angular distribution function
@ -224,29 +223,29 @@ class ScattDataTabular: public ScattData {
// Friend convert_legendre_to_tabular so it has access to protected
// parameters
friend void
convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
int n_mu);
convert_legendre_to_tabular(ScattDataLegendre& leg,
ScattDataTabular& tab, int n_mu);
public:
void
init(int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& in_mult,
double_3dvec& coeffs);
init(const int_1dvec& in_gmin, const int_1dvec& in_gmax,
const double_2dvec& in_mult, const double_3dvec& coeffs);
void
combine(const std::vector<ScattData*>& those_scatts,
const double_1dvec& scalars);
const double_1dvec& scalars);
double
calc_f(const int gin, const int gout, const double mu);
calc_f(int gin, int gout, double mu);
void
sample(const int gin, int& gout, double& mu, double& wgt);
sample(int gin, int& gout, double& mu, double& wgt);
int
get_order() {return dist[0][0].size();};
double_3dvec get_matrix(const int max_order);
double_3dvec get_matrix(int max_order);
};
//==============================================================================

196
src/tallies/tally.F90
View file

@ -2,10 +2,6 @@ module tally
use, intrinsic :: ISO_C_BINDING
#ifdef _OPENMP
use omp_lib
#endif
use algorithm, only: binary_search
use constants
use dict_header, only: EMPTY
@ -1233,12 +1229,6 @@ contains
real(8) :: p_uvw(3) ! Particle's current uvw
integer :: p_g ! Particle group to use for getting info
! to tally with.
integer(C_INT) :: tid
#ifdef _OPENMP
tid = OMP_GET_THREAD_NUM()
#else
tid = 0
#endif
! Set the direction and group to use with get_xs
if (t % estimator == ESTIMATOR_ANALOG .or. &
@ -1276,13 +1266,13 @@ contains
! To significantly reduce de-referencing, point matxs to the
! macroscopic Mgxs for the material of interest
call set_macro_angle_index_c(p % material, tid, p_uvw)
call set_macro_angle_index_c(p % material, p_uvw)
! Do same for nucxs, point it to the microscopic nuclide data of interest
if (i_nuclide > 0) then
! And since we haven't calculated this temperature index yet, do so now
call set_nuclide_temperature_index_c(i_nuclide, tid, p % sqrtkT)
call set_nuclide_angle_index_c(i_nuclide, tid, p_uvw)
call set_nuclide_temperature_index_c(i_nuclide, p % sqrtkT)
call set_nuclide_angle_index_c(i_nuclide, p_uvw)
end if
i = 0
@ -1337,13 +1327,13 @@ contains
if (i_nuclide > 0) then
score = score * flux * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_TOTAL, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_TOTAL, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_TOTAL, p_g)
end if
else
if (i_nuclide > 0) then
score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_TOTAL, p_g) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) * &
atom_density * flux
else
score = material_xs % total * flux
@ -1366,22 +1356,22 @@ contains
end if
if (i_nuclide > 0) then
score = score * flux * get_nuclide_xs_c(i_nuclide, tid, &
score = score * flux * get_nuclide_xs_c(i_nuclide, &
MG_GET_XS_INVERSE_VELOCITY, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * flux * get_macro_xs_c(p % material, tid, &
score = score * flux * get_macro_xs_c(p % material, &
MG_GET_XS_INVERSE_VELOCITY, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = flux * get_nuclide_xs_c(i_nuclide, tid, &
score = flux * get_nuclide_xs_c(i_nuclide, &
MG_GET_XS_INVERSE_VELOCITY, p_g)
else
score = flux * get_macro_xs_c(p % material, tid, &
score = flux * get_macro_xs_c(p % material, &
MG_GET_XS_INVERSE_VELOCITY, p_g)
end if
end if
@ -1404,22 +1394,22 @@ contains
! adjust the score by the actual probability for that nuclide.
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_SCATTER_FMU_MULT, &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_FMU_MULT, &
p % last_g, p % g, MU=p % mu) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_SCATTER_FMU_MULT, &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU_MULT, &
p % last_g, p % g, MU=p % mu)
end if
else
if (i_nuclide > 0) then
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_SCATTER_MULT, &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_MULT, &
p_g, MU=p % mu)
else
! Get the scattering x/s and take away
! the multiplication baked in to sigS
score = flux * &
get_macro_xs_c(p % material, tid, MG_GET_XS_SCATTER_MULT, &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER_MULT, &
p_g, MU=p % mu)
end if
end if
@ -1442,20 +1432,20 @@ contains
! adjust the score by the actual probability for that nuclide.
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_SCATTER_FMU, &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER_FMU, &
p % last_g, p % g, MU=p % mu) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_SCATTER_FMU, &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU, &
p % last_g, p % g, MU=p % mu)
end if
else
if (i_nuclide > 0) then
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_SCATTER, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER, p_g)
else
! Get the scattering x/s, which includes multiplication
score = flux * &
get_macro_xs_c(p % material, tid, MG_GET_XS_SCATTER, p_g)
get_macro_xs_c(p % material, MG_GET_XS_SCATTER, p_g)
end if
end if
@ -1475,13 +1465,13 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_ABSORPTION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_ABSORPTION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g)
else
score = material_xs % absorption * flux
end if
@ -1506,19 +1496,19 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) * &
atom_density * flux
else
score = get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g) * flux
score = get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
end if
end if
@ -1544,12 +1534,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
get_macro_xs_c(p % material, tid, MG_GET_XS_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
! Skip any non-fission events
@ -1562,17 +1552,17 @@ contains
score = keff * p % wgt_bank * flux
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
end if
else
if (i_nuclide > 0) then
score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_NU_FISSION, p_g) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) * &
atom_density * flux
else
score = get_macro_xs_c(p % material, tid, MG_GET_XS_NU_FISSION, p_g) * flux
score = get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) * flux
end if
end if
@ -1598,12 +1588,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
get_macro_xs_c(p % material, tid, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
! Skip any non-fission events
@ -1617,17 +1607,17 @@ contains
/ real(p % n_bank, 8)) * flux
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
end if
else
if (i_nuclide > 0) then
score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_PROMPT_NU_FISSION, p_g) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) * &
atom_density * flux
else
score = get_macro_xs_c(p % material, tid, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux
score = get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux
end if
end if
@ -1653,7 +1643,7 @@ contains
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! nu-fission
if (get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (dg_filter > 0) then
select type(filt => filters(t % filter(dg_filter)) % obj)
@ -1669,12 +1659,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1685,12 +1675,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
end if
end if
@ -1720,8 +1710,8 @@ contains
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1732,8 +1722,8 @@ contains
score = keff * p % wgt_bank / p % n_bank * sum(p % n_delayed_bank) * flux
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
end if
end if
end if
@ -1753,10 +1743,10 @@ contains
if (i_nuclide > 0) then
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
else
score = flux * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1766,11 +1756,11 @@ contains
else
if (i_nuclide > 0) then
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g)
else
score = flux * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g)
end if
end if
end if
@ -1785,7 +1775,7 @@ contains
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! nu-fission
if (get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (dg_filter > 0) then
select type(filt => filters(t % filter(dg_filter)) % obj)
@ -1801,14 +1791,14 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1826,14 +1816,14 @@ contains
do d = 1, num_delayed_groups
if (i_nuclide > 0) then
score = score + p % absorb_wgt * flux * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score + p % absorb_wgt * flux * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
end do
end if
@ -1862,13 +1852,13 @@ contains
if (i_nuclide > 0) then
score = score + keff * atom_density * &
fission_bank(n_bank - p % n_bank + k) % wgt * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g) * flux
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
else
score = score + keff * &
fission_bank(n_bank - p % n_bank + k) % wgt * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
end if
! if the delayed group filter is present, tally to corresponding
@ -1921,12 +1911,12 @@ contains
if (i_nuclide > 0) then
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
else
score = flux * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1943,12 +1933,12 @@ contains
do d = 1, num_delayed_groups
if (i_nuclide > 0) then
score = score + atom_density * flux * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
else
score = score + flux * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
end if
end do
end if
@ -1973,20 +1963,20 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_KAPPA_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
else
score = score * &
get_macro_xs_c(p % material, tid, MG_GET_XS_KAPPA_FISSION, p_g) / &
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_KAPPA_FISSION, p_g) * &
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) * &
atom_density * flux
else
score = flux * &
get_macro_xs_c(p % material, tid, MG_GET_XS_KAPPA_FISSION, p_g)
get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g)
end if
end if

12
src/tracking.F90
View file

@ -2,10 +2,6 @@ module tracking
use, intrinsic :: ISO_C_BINDING
#ifdef _OPENMP
use omp_lib
#endif
use constants
use error, only: warning, write_message
use geometry_header, only: cells
@ -52,12 +48,6 @@ contains
real(8) :: d_collision ! sampled distance to collision
real(8) :: distance ! distance particle travels
logical :: found_cell ! found cell which particle is in?
integer(C_INT) :: tid
#ifdef _OPENMP
tid = OMP_GET_THREAD_NUM()
#else
tid = 0
#endif
! Display message if high verbosity or trace is on
if (verbosity >= 9 .or. trace) then
@ -124,7 +114,7 @@ contains
end if
else
! Get the MG data
call calculate_xs_c(p % material, tid, p % g, p % sqrtkT, &
call calculate_xs_c(p % material, p % g, p % sqrtkT, &
p % coord(p % n_coord) % uvw, material_xs % total, &
material_xs % absorption, material_xs % nu_fission)

28
src/xsdata.cpp
View file

@ -6,9 +6,8 @@ namespace openmc {
// XsData class methods
//==============================================================================
XsData::XsData(const int energy_groups, const int num_delayed_groups,
const bool fissionable, const int scatter_format,
const int n_pol, const int n_azi)
XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable,
int scatter_format, int n_pol, int n_azi)
{
int n_ang = n_pol * n_azi;
@ -60,11 +59,9 @@ XsData::XsData(const int energy_groups, const int num_delayed_groups,
//==============================================================================
void
XsData::from_hdf5(const hid_t xsdata_grp, const bool fissionable,
const int scatter_format, const int final_scatter_format,
const int order_data, const int max_order,
const int legendre_to_tabular_points, const bool is_isotropic,
const int n_pol, const int n_azi)
XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
int final_scatter_format, int order_data, int max_order,
int legendre_to_tabular_points, bool is_isotropic, int n_pol, int n_azi)
{
// Reconstruct the dimension information so it doesn't need to be passed
int n_ang = n_pol * n_azi;
@ -83,8 +80,7 @@ XsData::from_hdf5(const hid_t xsdata_grp, const bool fissionable,
// Get scattering data
scatter_from_hdf5(xsdata_grp, n_pol, n_azi, energy_groups, scatter_format,
final_scatter_format, order_data, max_order,
legendre_to_tabular_points);
final_scatter_format, order_data, max_order, legendre_to_tabular_points);
// Check absorption to ensure it is not 0 since it is often the
// denominator in tally methods
@ -116,9 +112,8 @@ XsData::from_hdf5(const hid_t xsdata_grp, const bool fissionable,
//==============================================================================
void
XsData::fission_from_hdf5(const hid_t xsdata_grp, const int n_pol,
const int n_azi, const int energy_groups, const int delayed_groups,
const bool is_isotropic)
XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
int energy_groups, int delayed_groups, bool is_isotropic)
{
int n_ang = n_pol * n_azi;
// Get the fission and kappa_fission data xs; these are optional
@ -485,10 +480,9 @@ XsData::fission_from_hdf5(const hid_t xsdata_grp, const int n_pol,
//==============================================================================
void
XsData::scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol,
const int n_azi, const int energy_groups, int scatter_format,
const int final_scatter_format, const int order_data, const int max_order,
const int legendre_to_tabular_points)
XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi,
int energy_groups, int scatter_format, int final_scatter_format,
int order_data, int max_order, int legendre_to_tabular_points)
{
int n_ang = n_pol * n_azi;
if (!object_exists(xsdata_grp, "scatter_data")) {

28
src/xsdata.h
View file

@ -27,19 +27,17 @@ class XsData {
private:
//! \brief Reads scattering data from the HDF5 file
void
scatter_from_hdf5(const hid_t xsdata_grp, const int n_pol, const int n_azi,
const int energy_groups, int scatter_format,
const int final_scatter_format, const int order_data,
const int max_order, const int legendre_to_tabular_points);
scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups,
int scatter_format, int final_scatter_format, int order_data,
int max_order, int legendre_to_tabular_points);
//! \brief Reads fission data from the HDF5 file
void
fission_from_hdf5(const hid_t xsdata_grp, const int n_pol, const int n_azi,
const int energy_groups, const int delayed_groups,
const bool is_isotropic);
fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups,
int delayed_groups, bool is_isotropic);
public:
// The following quantities have the following dimensions:
// [angle][incoming group]
double_2dvec total;
@ -75,9 +73,8 @@ class XsData {
//! @param scatter_format The scattering representation of the file.
//! @param n_pol Number of polar angles.
//! @param n_azi Number of azimuthal angles.
XsData(const int num_groups, const int num_delayed_groups,
const bool fissionable, const int scatter_format, const int n_pol,
const int n_azi);
XsData(int num_groups, int num_delayed_groups, bool fissionable,
int scatter_format, int n_pol, int n_azi);
//! \brief Loads the XsData object from the HDF5 file
//!
@ -99,11 +96,10 @@ class XsData {
//! @param n_pol Number of polar angles.
//! @param n_azi Number of azimuthal angles.
void
from_hdf5(const hid_t xsdata_grp, const bool fissionable,
const int scatter_format, const int final_scatter_format,
const int order_data, const int max_order,
const int legendre_to_tabular_points, const bool is_isotropic,
const int n_pol, const int n_azi);
from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
int final_scatter_format, int order_data, int max_order,
int legendre_to_tabular_points, bool is_isotropic, int n_pol,
int n_azi);
//! \brief Combines the microscopic data to a macroscopic object.
//!