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

got it working with OpenMP

Browse source
This commit is contained in:
Adam G Nelson 2018-06-14 20:22:31 -04:00
parent 7fd4cbb140
commit 7c819f02bc
9 changed files with 404 additions and 416 deletions
Download patch file Download diff file Expand all files Collapse all files

156
src/mgxs.cpp
View file

@ -15,7 +15,7 @@ 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 double_1dvec& in_polar,
const double_1dvec& in_azimuthal)
const double_1dvec& in_azimuthal, const int n_threads)
{
name = in_name;
awr = in_awr;
@ -29,20 +29,23 @@ void Mgxs::init(const std::string& in_name, const double in_awr,
azimuthal = in_azimuthal;
n_pol = polar.size();
n_azi = azimuthal.size();
index_temp = 0;
last_sqrtkT = 0.;
index_pol = 0;
index_azi = 0;
last_uvw[0] = 1.;
last_uvw[1] = 0.;
last_uvw[2] = 0.;
cache.resize(n_threads);
for (int thread = 0; thread < n_threads; thread++) {
cache[thread].sqrtkT = 0.;
cache[thread].t = 0;
cache[thread].p = 0;
cache[thread].a = 0;
cache[thread].uvw[0] = 1.;
cache[thread].uvw[1] = 0.;
cache[thread].uvw[2] = 0.;
}
}
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,
bool& is_isotropic)
bool& is_isotropic, const int n_threads)
{
// get name
char char_name[MAX_WORD_LEN];
@ -242,21 +245,23 @@ void 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_polar, in_azimuthal);
in_num_groups, in_num_delayed_groups, in_polar, in_azimuthal,
n_threads);
}
void Mgxs::from_hdf5(hid_t xs_id, int energy_groups, int delayed_groups,
double_1dvec& temperature, int& method, double tolerance,
int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points)
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)
{
// Call generic data gathering routine (will populate the metadata)
int order_data;
int_1dvec temps_to_read;
bool is_isotropic;
_metadata_from_hdf5(xs_id, energy_groups, delayed_groups, temperature,
method, tolerance, temps_to_read, order_data, is_isotropic);
method, tolerance, temps_to_read, order_data, is_isotropic, n_threads);
// Set number of energy and delayed groups
int final_scatter_format = scatter_format;
@ -285,8 +290,8 @@ void Mgxs::from_hdf5(hid_t xs_id, int energy_groups, int delayed_groups,
void Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs,
std::vector<Mgxs*>& micros, double_1dvec& atom_densities,
int& method, double tolerance)
std::vector<Mgxs*>& micros, double_1dvec& atom_densities, int& method,
const double tolerance, const int n_threads)
{
// Get the minimum data needed to initialize:
// Dont need awr, but lets just initialize it anyways
@ -305,7 +310,8 @@ void Mgxs::build_macro(const std::string& in_name, double_1dvec& mat_kTs,
double_1dvec in_azimuthal = micros[0]->azimuthal;
init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format,
in_num_groups, in_num_delayed_groups, in_polar, in_azimuthal);
in_num_groups, in_num_delayed_groups, in_polar, in_azimuthal,
n_threads);
// Create the xs data for each temperature
for (int t = 0; t < mat_kTs.size(); t++) {
@ -395,52 +401,52 @@ void Mgxs::combine(std::vector<Mgxs*>& micros, double_1dvec& scalars,
}
double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
int* dg)
double Mgxs::get_xs(const int tid, const int xstype, const int gin,
int* gout, double* mu, int* dg)
{
// This method assumes that the temperature and angle indices are set
double val;
switch(xstype) {
case MG_GET_XS_TOTAL:
val = xs[index_temp].total[index_pol][index_azi][gin];
val = xs[cache[tid].t].total[cache[tid].p][cache[tid].a][gin];
break;
case MG_GET_XS_ABSORPTION:
val = xs[index_temp].absorption[index_pol][index_azi][gin];
val = xs[cache[tid].t].absorption[cache[tid].p][cache[tid].a][gin];
break;
case MG_GET_XS_INVERSE_VELOCITY:
val = xs[index_temp].inverse_velocity[index_pol][index_azi][gin];
val = xs[cache[tid].t].inverse_velocity[cache[tid].p][cache[tid].a][gin];
break;
case MG_GET_XS_DECAY_RATE:
if (dg != nullptr) {
val = xs[index_temp].decay_rate[index_pol][index_azi][*dg + 1];
val = xs[cache[tid].t].decay_rate[cache[tid].p][cache[tid].a][*dg + 1];
} else {
val = xs[index_temp].decay_rate[index_pol][index_azi][0];
val = xs[cache[tid].t].decay_rate[cache[tid].p][cache[tid].a][0];
}
break;
case MG_GET_XS_SCATTER:
case MG_GET_XS_SCATTER_MULT:
case MG_GET_XS_SCATTER_FMU_MULT:
case MG_GET_XS_SCATTER_FMU:
val = xs[index_temp].scatter[index_pol]
[index_azi]->get_xs(xstype, gin, gout, mu);
val = xs[cache[tid].t].scatter[cache[tid].p]
[cache[tid].a]->get_xs(xstype, gin, gout, mu);
break;
case MG_GET_XS_FISSION:
if (fissionable) {
val = xs[index_temp].fission[index_pol][index_azi][gin];
val = xs[cache[tid].t].fission[cache[tid].p][cache[tid].a][gin];
} else {
val = 0.;
}
break;
case MG_GET_XS_KAPPA_FISSION:
if (fissionable) {
val = xs[index_temp].kappa_fission[index_pol][index_azi][gin];
val = xs[cache[tid].t].kappa_fission[cache[tid].p][cache[tid].a][gin];
} else {
val = 0.;
}
break;
case MG_GET_XS_PROMPT_NU_FISSION:
if (fissionable) {
val = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
val = xs[cache[tid].t].prompt_nu_fission[cache[tid].p][cache[tid].a][gin];
} else {
val = 0.;
}
@ -448,11 +454,11 @@ double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
case MG_GET_XS_DELAYED_NU_FISSION:
if (fissionable) {
if (dg != nullptr) {
val = xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][*dg];
val = xs[cache[tid].t].delayed_nu_fission[cache[tid].p][cache[tid].a][gin][*dg];
} else {
val = 0.;
for (auto& num : xs[index_temp].delayed_nu_fission[index_pol]
[index_azi][gin]) {
for (auto& num : xs[cache[tid].t].delayed_nu_fission[cache[tid].p]
[cache[tid].a][gin]) {
val += num;
}
}
@ -462,7 +468,7 @@ double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
break;
case MG_GET_XS_NU_FISSION:
if (fissionable) {
val = xs[index_temp].nu_fission[index_pol][index_azi][gin];
val = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
} else {
val = 0.;
}
@ -470,11 +476,11 @@ double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
case MG_GET_XS_CHI_PROMPT:
if (fissionable) {
if (gout != nullptr) {
val = xs[index_temp].chi_prompt[index_pol][index_azi][gin][*gout];
val = xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][*gout];
} else {
// provide an outgoing group-wise sum
val = 0.;
for (auto& num : xs[index_temp].chi_prompt[index_pol][index_azi][gin]) {
for (auto& num : xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin]) {
val += num;
}
}
@ -486,23 +492,23 @@ double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
if (fissionable) {
if (gout != nullptr) {
if (dg != nullptr) {
val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][*dg];
val = xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][*gout][*dg];
} else {
val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][0];
val = xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][*gout][0];
}
} else {
if (dg != nullptr) {
val = 0.;
for (int i = 0; i < xs[index_temp].chi_delayed[index_pol]
[index_azi][gin].size(); i++) {
val += xs[index_temp].chi_delayed[index_pol][index_azi][gin][i][*dg];
for (int i = 0; i < xs[cache[tid].t].chi_delayed[cache[tid].p]
[cache[tid].a][gin].size(); i++) {
val += xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][i][*dg];
}
} else {
val = 0.;
for (int i = 0; i < xs[index_temp].chi_delayed[index_pol]
[index_azi][gin].size(); i++) {
for (auto& num : xs[index_temp].chi_delayed[index_pol]
[index_azi][gin][i]) {
for (int i = 0; i < xs[cache[tid].t].chi_delayed[cache[tid].p]
[cache[tid].a][gin].size(); i++) {
for (auto& num : xs[cache[tid].t].chi_delayed[cache[tid].p]
[cache[tid].a][gin][i]) {
val += num;
}
}
@ -519,14 +525,15 @@ double Mgxs::get_xs(const int xstype, const int gin, int* gout, double* mu,
}
void Mgxs::sample_fission_energy(const int gin, int& dg, int& gout)
void Mgxs::sample_fission_energy(const int tid, const int gin, int& dg,
int& gout)
{
// This method assumes that the temperature and angle indices are set
double nu_fission = xs[index_temp].nu_fission[index_pol][index_azi][gin];
double nu_fission = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
// Find the probability of having a prompt neutron
double prob_prompt =
xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin];
xs[cache[tid].t].prompt_nu_fission[cache[tid].p][cache[tid].a][gin];
// sample random numbers
double xi_pd = prn() * nu_fission;
@ -542,10 +549,10 @@ void Mgxs::sample_fission_energy(const int gin, int& dg, int& gout)
// sample the outgoing energy group
gout = 0;
double prob_gout =
xs[index_temp].chi_prompt[index_pol][index_azi][gin][gout];
xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][gout];
while (prob_gout < xi_gout) {
gout++;
prob_gout += xs[index_temp].chi_prompt[index_pol][index_azi][gin][gout];
prob_gout += xs[cache[tid].t].chi_prompt[cache[tid].p][cache[tid].a][gin][gout];
}
} else {
@ -556,7 +563,7 @@ void Mgxs::sample_fission_energy(const int gin, int& dg, int& gout)
while (xi_pd >= prob_prompt) {
dg++;
prob_prompt +=
xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][dg];
xs[cache[tid].t].delayed_nu_fission[cache[tid].p][cache[tid].a][gin][dg];
}
// adjust dg in case of round-off error
@ -565,35 +572,36 @@ void Mgxs::sample_fission_energy(const int gin, int& dg, int& gout)
// sample the outgoing energy group
gout = 0;
double prob_gout =
xs[index_temp].chi_delayed[index_pol][index_azi][gin][gout][dg];
xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][gout][dg];
while (prob_gout < xi_gout) {
gout++;
prob_gout +=
xs[index_temp].chi_delayed[index_pol][index_azi][gin][gout][dg];
xs[cache[tid].t].chi_delayed[cache[tid].p][cache[tid].a][gin][gout][dg];
}
}
}
void Mgxs::sample_scatter(const int gin, int& gout, double& mu, double& wgt)
void Mgxs::sample_scatter(const int tid, const int gin, int& gout, double& mu,
double& wgt)
{
// This method assumes that the temperature and angle indices are set
// Sample the data
xs[index_temp].scatter[index_pol][index_azi]->sample(gin, gout, mu, wgt);
xs[cache[tid].t].scatter[cache[tid].p][cache[tid].a]->sample(gin, gout, mu, wgt);
}
void Mgxs::calculate_xs(const int gin, const double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs)
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)
{
// Set our indices
set_temperature_index(sqrtkT);
set_angle_index(uvw);
total_xs = xs[index_temp].total[index_pol][index_azi][gin];
abs_xs = xs[index_temp].absorption[index_pol][index_azi][gin];
set_temperature_index(tid, sqrtkT);
set_angle_index(tid, uvw);
total_xs = xs[cache[tid].t].total[cache[tid].p][cache[tid].a][gin];
abs_xs = xs[cache[tid].t].absorption[cache[tid].p][cache[tid].a][gin];
if (fissionable) {
nu_fiss_xs = xs[index_temp].nu_fission[index_pol][index_azi][gin];
nu_fiss_xs = xs[cache[tid].t].nu_fission[cache[tid].p][cache[tid].a][gin];
} else {
nu_fiss_xs = 0.;
}
@ -617,10 +625,10 @@ bool Mgxs::equiv(const Mgxs& that)
}
void Mgxs::set_temperature_index(const double sqrtkT)
void Mgxs::set_temperature_index(const int tid, const double sqrtkT)
{
// See if we need to find the new index
if (sqrtkT != last_sqrtkT) {
if (sqrtkT != cache[tid].sqrtkT) {
double kT = sqrtkT * sqrtkT;
// initialize vector for storage of the differences
@ -628,34 +636,34 @@ void Mgxs::set_temperature_index(const double sqrtkT)
// Find the minimum difference of kT and kTs
temp_diff = std::abs(temp_diff - kT);
index_temp = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
cache[tid].t = std::min_element(std::begin(temp_diff), std::end(temp_diff)) -
std::begin(temp_diff);
// store this temperature as the last one used
last_sqrtkT = sqrtkT;
cache[tid].sqrtkT = sqrtkT;
}
}
void Mgxs::set_angle_index(const double uvw[3])
void Mgxs::set_angle_index(const int tid, const double uvw[3])
{
// See if we need to find the new index
if ((uvw[0] != last_uvw[0]) || (uvw[1] != last_uvw[1]) ||
(uvw[2] != last_uvw[2])) {
if ((uvw[0] != cache[tid].uvw[0]) || (uvw[1] != cache[tid].uvw[1]) ||
(uvw[2] != cache[tid].uvw[2])) {
// convert uvw to polar and azimuthal angles
double my_pol = std::acos(uvw[2]);
double my_azi = std::atan2(uvw[1], uvw[0]);
// Find the location, assuming equal-bin angles
double delta_angle = PI / n_pol;
index_pol = std::floor(my_pol / delta_angle);
cache[tid].p = std::floor(my_pol / delta_angle);
delta_angle = 2. * PI / n_azi;
index_azi = std::floor((my_azi + PI) / delta_angle);
cache[tid].a = std::floor((my_azi + PI) / delta_angle);
// store this direction as the last one used
last_uvw[0] = uvw[0];
last_uvw[1] = uvw[1];
last_uvw[2] = uvw[2];
cache[tid].uvw[0] = uvw[0];
cache[tid].uvw[1] = uvw[1];
cache[tid].uvw[2] = uvw[2];
}
}

50
src/mgxs.h
View file

@ -23,6 +23,18 @@
namespace openmc {
//==============================================================================
// Cache contains the cached data for an MGXS object
//==============================================================================
struct CacheData {
double sqrtkT; // last temperature corresponding to t
int t; // temperature index
int p; // polar angle index
int a; // azimuthal angle index
double uvw[3]; // last angle that corresponds to p and a
};
//==============================================================================
// MGXS contains the mgxs data for a nuclide/material
//==============================================================================
@ -33,7 +45,6 @@ class Mgxs {
int scatter_format; // flag for if this is legendre, histogram, or tabular
int num_delayed_groups; // number of delayed neutron groups
int num_groups; // number of energy groups
double last_sqrtkT; // cache of the temperature corresponding to index_temp
std::vector<XsData> xs; // Cross section data
int n_pol;
int n_azi;
@ -42,7 +53,7 @@ class Mgxs {
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, bool& is_isotropic);
int& order_dim, bool& is_isotropic, const int n_threads);
bool equiv(const Mgxs& that);
public:
@ -50,32 +61,31 @@ class Mgxs {
double awr; // atomic weight ratio
bool fissionable; // Is this fissionable
// TODO: The following attributes be private when Fortran is fully replaced
int index_pol; // cache for the angle indices
int index_azi;
double last_uvw[3]; // cache of the angle corresponding to the above indices
int index_temp; // cache of temperature index
std::vector<CacheData> cache; // index and data cache
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 double_1dvec& in_polar,
const double_1dvec& in_azimuthal);
const double_1dvec& in_azimuthal, const int n_threads);
void build_macro(const std::string& in_name, double_1dvec& mat_kTs,
std::vector<Mgxs*>& micros, double_1dvec& atom_densities,
int& method, double tolerance);
int& method, const double tolerance, const int n_threads);
void combine(std::vector<Mgxs*>& micros, double_1dvec& scalars,
int_1dvec& micro_ts, int this_t);
void from_hdf5(hid_t xs_id, int energy_groups, int delayed_groups,
double_1dvec& temperature, int& method, double tolerance,
int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points);
double get_xs(const int xstype, const int gin, int* gout, double* mu,
int* dg);
void sample_fission_energy(const int gin, int& dg, int& gout);
void sample_scatter(const int gin, int& gout, double& mu, double& wgt);
void calculate_xs(const int gin, const double sqrtkT, const double uvw[3],
double& total_xs, double& abs_xs, double& nu_fiss_xs);
void set_temperature_index(const double sqrtkT);
void set_angle_index(const double uvw[3]);
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);
double get_xs(const int tid, const int xstype, const int gin, int* gout,
double* mu, int* dg);
void sample_fission_energy(const int tid, const int gin, int& dg, int& gout);
void sample_scatter(const int tid, const int gin, int& gout, double& mu,
double& wgt);
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);
void set_temperature_index(const int tid, const double sqrtkT);
void set_angle_index(const int tid, const double uvw[3]);
};
} // namespace openmc

23
src/mgxs_data.F90
View file

@ -2,6 +2,10 @@ 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
@ -36,6 +40,13 @@ 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)
@ -84,7 +95,8 @@ contains
num_energy_groups, num_delayed_groups, &
temps(i_nuclide) % size(), temps(i_nuclide) % data, &
temperature_method, temperature_tolerance, max_order, &
logical(legendre_to_tabular, C_BOOL), legendre_to_tabular_points)
logical(legendre_to_tabular, C_BOOL), &
legendre_to_tabular_points, n_threads)
call already_read % add(name)
end if
@ -110,6 +122,13 @@ 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)
@ -130,7 +149,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)
temperature_method, temperature_tolerance, n_threads)
end if
end do

69
src/mgxs_interface.F90
View file

@ -10,7 +10,7 @@ module mgxs_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) bind(C)
legendre_to_tabular_points, n_threads) bind(C)
use ISO_C_BINDING
import HID_T
implicit none
@ -25,6 +25,7 @@ module mgxs_interface
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
end subroutine add_mgxs_c
function query_fissionable_c(n_nuclides, i_nuclides) result(result) bind(C)
@ -36,7 +37,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) bind(C)
atom_densities, method, tolerance, n_threads) bind(C)
use ISO_C_BINDING
implicit none
character(kind=C_CHAR),intent(in) :: name(*)
@ -47,13 +48,15 @@ module mgxs_interface
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
end subroutine create_macro_xs_c
subroutine calculate_xs_c(i_mat, gin, sqrtkT, uvw, total_xs, abs_xs, &
subroutine calculate_xs_c(i_mat, tid, 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)
@ -62,10 +65,11 @@ module mgxs_interface
real(C_DOUBLE), intent(inout) :: nu_fiss_xs
end subroutine calculate_xs_c
subroutine sample_scatter_c(i_mat, gin, gout, mu, wgt, uvw) bind(C)
subroutine sample_scatter_c(i_mat, tid, 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
@ -73,10 +77,11 @@ module mgxs_interface
real(C_DOUBLE), intent(inout) :: uvw(1:3)
end subroutine sample_scatter_c
subroutine sample_fission_energy_c(i_mat, gin, dg, gout) bind(C)
subroutine sample_fission_energy_c(i_mat, tid, 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
@ -97,11 +102,12 @@ module mgxs_interface
real(C_DOUBLE) :: awr
end function get_awr_c
function get_nuclide_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
function get_nuclide_xs_c(index, tid, 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
@ -110,11 +116,12 @@ module mgxs_interface
real(C_DOUBLE) :: val
end function get_nuclide_xs_c
function get_macro_xs_c(index, xstype, gin, gout, mu, dg) result(val) &
function get_macro_xs_c(index, tid, 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
@ -123,63 +130,29 @@ module mgxs_interface
real(C_DOUBLE) :: val
end function get_macro_xs_c
subroutine set_nuclide_angle_index_c(index, uvw, last_pol, last_azi, &
last_uvw) bind(C)
subroutine set_nuclide_angle_index_c(index, tid, 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)
integer(C_INT), intent(inout) :: last_pol
integer(C_INT), intent(inout) :: last_azi
real(C_DOUBLE), intent(inout) :: last_uvw(1:3)
end subroutine set_nuclide_angle_index_c
subroutine reset_nuclide_angle_index_c(index, last_pol, last_azi, &
last_uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: last_pol
integer(C_INT), value, intent(in) :: last_azi
real(C_DOUBLE), intent(in) :: last_uvw(1:3)
end subroutine reset_nuclide_angle_index_c
subroutine set_macro_angle_index_c(index, uvw, last_pol, last_azi, &
last_uvw) bind(C)
subroutine set_macro_angle_index_c(index, tid, 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)
integer(C_INT), intent(inout) :: last_pol
integer(C_INT), intent(inout) :: last_azi
real(C_DOUBLE), intent(inout) :: last_uvw(1:3)
end subroutine set_macro_angle_index_c
subroutine reset_macro_angle_index_c(index, last_pol, last_azi, &
last_uvw) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: last_pol
integer(C_INT), value, intent(in) :: last_azi
real(C_DOUBLE), intent(in) :: last_uvw(1:3)
end subroutine reset_macro_angle_index_c
function set_nuclide_temperature_index_c(index, sqrtkT) result(last_temp) &
bind(C)
subroutine set_nuclide_temperature_index_c(index, tid, 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
integer(C_INT) :: last_temp
end function set_nuclide_temperature_index_c
subroutine reset_nuclide_temperature_index_c(index, last_temp) bind(C)
use ISO_C_BINDING
implicit none
integer(C_INT), value, intent(in) :: index
integer(C_INT), value, intent(in) :: last_temp
end subroutine reset_nuclide_temperature_index_c
end subroutine set_nuclide_temperature_index_c
end interface

221
src/mgxs_interface.cpp
View file

@ -6,10 +6,11 @@ namespace openmc {
// Mgxs data loading interface methods
//==============================================================================
void add_mgxs_c(hid_t file_id, char* name, int energy_groups,
int delayed_groups, int n_temps, double temps[], int& method,
double tolerance, int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points)
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)
{
//!! mgxs_data.F90 will be modified to just create the list of names
//!! in the order needed
@ -32,7 +33,7 @@ void add_mgxs_c(hid_t file_id, char* name, int energy_groups,
Mgxs mg;
mg.from_hdf5(xs_grp, energy_groups, delayed_groups,
temperature, method, tolerance, max_order, legendre_to_tabular,
legendre_to_tabular_points);
legendre_to_tabular_points, n_threads);
nuclides_MG.push_back(mg);
}
@ -50,7 +51,8 @@ bool 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 double atom_densities[], int& method, const double tolerance,
const int n_threads)
{
Mgxs macro;
if (n_temps > 0) {
@ -70,7 +72,7 @@ void create_macro_xs_c(char* mat_name, const int n_nuclides,
}
macro.build_macro(mat_name, temperature, mgxs_ptr, atom_densities_vec,
method, tolerance);
method, tolerance, n_threads);
}
macro_xs.push_back(macro);
}
@ -79,19 +81,20 @@ void create_macro_xs_c(char* mat_name, const int n_nuclides,
// Mgxs tracking/transport/tallying interface methods
//==============================================================================
void calculate_xs_c(const int i_mat, const int gin, const double sqrtkT,
const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs)
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)
{
macro_xs[i_mat - 1].calculate_xs(gin - 1, sqrtkT, uvw, total_xs, abs_xs,
macro_xs[i_mat - 1].calculate_xs(tid, gin - 1, sqrtkT, uvw, total_xs, abs_xs,
nu_fiss_xs);
}
void sample_scatter_c(const int i_mat, const int gin, int& gout, double& mu,
double& wgt, double uvw[3])
void sample_scatter_c(const int i_mat, const int tid, const int gin, int& gout,
double& mu, double& wgt, double uvw[3])
{
int gout_c = gout - 1;
macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt);
macro_xs[i_mat - 1].sample_scatter(tid, gin - 1, gout_c, mu, wgt);
// adjust return value for fortran indexing
gout = gout_c + 1;
@ -101,11 +104,12 @@ void sample_scatter_c(const int i_mat, const int gin, int& gout, double& mu,
}
void sample_fission_energy_c(const int i_mat, const int gin, int& dg, int& gout)
void sample_fission_energy_c(const int i_mat, const int tid, const int gin,
int& dg, int& gout)
{
int dg_c = 0;
int gout_c = 0;
macro_xs[i_mat - 1].sample_fission_energy(gin - 1, dg_c, gout_c);
macro_xs[i_mat - 1].sample_fission_energy(tid, gin - 1, dg_c, gout_c);
// adjust return values for fortran indexing
dg = dg_c + 1;
@ -113,6 +117,82 @@ void sample_fission_energy_c(const int i_mat, const int gin, int& dg, int& gout)
}
double get_nuclide_xs_c(const int index, const int tid, const int xstype,
const int gin, int* gout, double* mu, int* dg)
{
int gout_c;
int* gout_c_p;
int dg_c;
int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
} else {
gout_c_p = gout;
}
if (dg != nullptr) {
dg_c = *dg - 1;
dg_c_p = &dg_c;
} else {
dg_c_p = dg;
}
return nuclides_MG[index - 1].get_xs(tid, 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)
{
int gout_c;
int* gout_c_p;
int dg_c;
int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
} else {
gout_c_p = gout;
}
if (dg != nullptr) {
dg_c = *dg - 1;
dg_c_p = &dg_c;
} else {
dg_c_p = dg;
}
return macro_xs[index - 1].get_xs(tid, 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])
{
// Update the values
nuclides_MG[index - 1].set_angle_index(tid, uvw);
}
void set_macro_angle_index_c(const int index, const int tid,
const double uvw[3])
{
// Update the values
macro_xs[index - 1].set_angle_index(tid, uvw);
}
void set_nuclide_temperature_index_c(const int index, const int tid,
const double sqrtkT)
{
// Update the values
nuclides_MG[index - 1].set_temperature_index(tid, sqrtkT);
}
//==============================================================================
// Mgxs general methods
//==============================================================================
void get_name_c(const int index, int name_len, char* name)
{
// First blank out our input string
@ -133,115 +213,4 @@ double get_awr_c(const int index)
return nuclides_MG[index - 1].awr;
}
double get_nuclide_xs_c(const int index, const int xstype, const int gin,
int* gout, double* mu, int* dg)
{
int gout_c;
int* gout_c_p;
int dg_c;
int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
} else {
gout_c_p = gout;
}
if (dg != nullptr) {
dg_c = *dg - 1;
dg_c_p = &dg_c;
} else {
dg_c_p = dg;
}
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 xstype, const int gin,
int* gout, double* mu, int* dg)
{
int gout_c;
int* gout_c_p;
int dg_c;
int* dg_c_p;
if (gout != nullptr) {
gout_c = *gout - 1;
gout_c_p = &gout_c;
} else {
gout_c_p = gout;
}
if (dg != nullptr) {
dg_c = *dg - 1;
dg_c_p = &dg_c;
} else {
dg_c_p = dg;
}
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 double uvw[3],
int& last_pol, int& last_azi, double last_uvw[3])
{
// Store the old
last_pol = nuclides_MG[index - 1].index_pol;
last_azi = nuclides_MG[index - 1].index_azi;
last_uvw[0] = nuclides_MG[index - 1].last_uvw[0];
last_uvw[1] = nuclides_MG[index - 1].last_uvw[1];
last_uvw[2] = nuclides_MG[index - 1].last_uvw[2];
// Update the values
nuclides_MG[index - 1].set_angle_index(uvw);
}
void reset_nuclide_angle_index_c(const int index, const int last_pol,
const int last_azi, const double last_uvw[3])
{
nuclides_MG[index - 1].index_pol = last_pol;
nuclides_MG[index - 1].index_azi = last_azi;
nuclides_MG[index - 1].last_uvw[0] = last_uvw[0];
nuclides_MG[index - 1].last_uvw[1] = last_uvw[1];
nuclides_MG[index - 1].last_uvw[2] = last_uvw[2];
}
void set_macro_angle_index_c(const int index, const double uvw[3],
int& last_pol, int& last_azi, double last_uvw[3])
{
// Store the old
last_pol = macro_xs[index - 1].index_pol;
last_azi = macro_xs[index - 1].index_azi;
last_uvw[0] = macro_xs[index - 1].last_uvw[0];
last_uvw[1] = macro_xs[index - 1].last_uvw[1];
last_uvw[2] = macro_xs[index - 1].last_uvw[2];
// Update the values
macro_xs[index - 1].set_angle_index(uvw);
}
void reset_macro_angle_index_c(const int index, const int last_pol,
const int last_azi, const double last_uvw[3])
{
macro_xs[index - 1].index_pol = last_pol;
macro_xs[index - 1].index_azi = last_azi;
macro_xs[index - 1].last_uvw[0] = last_uvw[0];
macro_xs[index - 1].last_uvw[1] = last_uvw[1];
macro_xs[index - 1].last_uvw[2] = last_uvw[2];
}
int set_nuclide_temperature_index_c(const int index, const double sqrtkT)
{
int old = nuclides_MG[index - 1].index_temp;
nuclides_MG[index - 1].set_temperature_index(sqrtkT);
return old;
}
void reset_nuclide_temperature_index_c(const int index, const int last_temp)
{
nuclides_MG[index - 1].index_temp = last_temp;
}
} // namespace openmc

49
src/mgxs_interface.h
View file

@ -13,54 +13,47 @@ extern std::vector<Mgxs> nuclides_MG;
extern std::vector<Mgxs> macro_xs;
extern "C" void add_mgxs_c(hid_t file_id, char* name, int energy_groups,
int delayed_groups, int n_temps, double temps[], int& method,
double tolerance, int max_order, bool legendre_to_tabular,
int legendre_to_tabular_points);
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);
extern "C" bool query_fissionable_c(const 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 double atom_densities[], int& method, const double tolerance,
const int n_threads);
extern "C" void calculate_xs_c(const int i_mat, const int gin,
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);
extern "C" void sample_scatter_c(const int i_mat, const int gin, int& gout,
double& mu, double& wgt, double uvw[3]);
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]);
extern "C" void sample_fission_energy_c(const int i_mat, const int gin,
int& dg, int& gout);
extern "C" void sample_fission_energy_c(const int i_mat, const int tid,
const int gin, int& dg, int& gout);
extern "C" void get_name_c(const int index, int name_len, char* name);
extern "C" double get_awr_c(const int index);
extern "C" double get_nuclide_xs_c(const int index, const int xstype,
const int gin, int* gout, double* mu, int* dg);
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);
extern "C" double get_macro_xs_c(const int index, const int xstype,
const 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);
extern "C" void set_nuclide_angle_index_c(const int index, const double uvw[3],
int& last_pol, int& last_azi, double last_uvw[3]);
extern "C" void set_nuclide_angle_index_c(const int index, const int tid,
const double uvw[3]);
extern "C" void reset_nuclide_angle_index_c(const int index, const int last_pol,
const int last_azi, const double last_uvw[3]);
extern "C" void set_macro_angle_index_c(const int index, const int tid,
const double uvw[3]);
extern "C" void set_macro_angle_index_c(const int index, const double uvw[3],
int& last_pol, int& last_azi, double last_uvw[3]);
extern "C" void reset_macro_angle_index_c(const int index, const int last_pol,
const int last_azi, const double last_uvw[3]);
extern "C" int set_nuclide_temperature_index_c(const int index,
extern "C" void set_nuclide_temperature_index_c(const int index, const int tid,
const double sqrtkT);
extern "C" void reset_nuclide_temperature_index_c(const int index,
const int last_temp);
} // namespace openmc
#endif // MGXS_INTERFACE_H

22
src/physics_mg.F90
View file

@ -2,6 +2,10 @@ 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
@ -141,9 +145,15 @@ 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, p % last_g, p % g, p % mu, p % wgt, &
p % coord(1) % uvw)
call sample_scatter_c(p % material, tid, p % last_g, p % g, p % mu, &
p % wgt, p % coord(1) % uvw)
! Update energy value for downstream compatability (in tallying)
p % E = energy_bin_avg(p % g)
@ -173,6 +183,12 @@ 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
@ -252,7 +268,7 @@ contains
! Sample secondary energy distribution for fission reaction and set energy
! in fission bank
call sample_fission_energy_c(p % material, p % g, dg, gout)
call sample_fission_energy_c(p % material, tid, p % g, dg, gout)
bank_array(i) % E = real(gout, 8)
bank_array(i) % delayed_group = dg

218
src/tallies/tally.F90
View file

@ -2,6 +2,10 @@ 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
@ -1229,14 +1233,12 @@ contains
real(8) :: p_uvw(3) ! Particle's current uvw
integer :: p_g ! Particle group to use for getting info
! to tally with.
! Storage of the indices the Mgxs object arrived with for resetting later
integer(C_INT) :: last_nuc_azi
integer(C_INT) :: last_nuc_pol
integer(C_INT) :: last_mat_azi
integer(C_INT) :: last_mat_pol
integer(C_INT) :: last_nuc_temp
real(C_DOUBLE) :: last_mat_uvw(3)
real(C_DOUBLE) :: last_nuc_uvw(3)
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. &
@ -1274,15 +1276,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, p_uvw, last_mat_pol, &
last_mat_azi, last_mat_uvw)
call set_macro_angle_index_c(p % material, tid, 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
last_nuc_temp = set_nuclide_temperature_index_c(i_nuclide, p % sqrtkT)
call set_nuclide_angle_index_c(i_nuclide, p_uvw, last_nuc_pol, &
last_nuc_azi, last_nuc_uvw)
call set_nuclide_temperature_index_c(i_nuclide, tid, p % sqrtkT)
call set_nuclide_angle_index_c(i_nuclide, tid, p_uvw)
end if
i = 0
@ -1336,14 +1336,14 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_TOTAL, p_g) * flux
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)
end if
else
if (i_nuclide > 0) then
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) * &
score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_TOTAL, p_g) * &
atom_density * flux
else
score = material_xs % total * flux
@ -1366,22 +1366,22 @@ contains
end if
if (i_nuclide > 0) then
score = score * get_nuclide_xs_c(i_nuclide, &
score = score * flux * get_nuclide_xs_c(i_nuclide, tid, &
MG_GET_XS_INVERSE_VELOCITY, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) * flux
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
else
score = score * get_macro_xs_c(p % material, &
score = score * flux * get_macro_xs_c(p % material, tid, &
MG_GET_XS_INVERSE_VELOCITY, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) * flux
get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g)
end if
else
if (i_nuclide > 0) then
score = flux * get_nuclide_xs_c(i_nuclide, &
score = flux * get_nuclide_xs_c(i_nuclide, tid, &
MG_GET_XS_INVERSE_VELOCITY, p_g)
else
score = flux * get_macro_xs_c(p % material, &
score = flux * get_macro_xs_c(p % material, tid, &
MG_GET_XS_INVERSE_VELOCITY, p_g)
end if
end if
@ -1404,22 +1404,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, MG_GET_XS_SCATTER_FMU_MULT, &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_SCATTER_FMU_MULT, &
p % last_g, p % g, MU=p % mu) / &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU_MULT, &
get_macro_xs_c(p % material, tid, 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, MG_GET_XS_SCATTER_MULT, &
get_nuclide_xs_c(i_nuclide, tid, 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, MG_GET_XS_SCATTER_MULT, &
get_macro_xs_c(p % material, tid, MG_GET_XS_SCATTER_MULT, &
p_g, MU=p % mu)
end if
end if
@ -1442,20 +1442,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, MG_GET_XS_SCATTER_FMU, &
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_SCATTER_FMU, &
p % last_g, p % g, MU=p % mu) / &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER_FMU, &
get_macro_xs_c(p % material, tid, 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, MG_GET_XS_SCATTER, p_g)
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_SCATTER, p_g)
else
! Get the scattering x/s, which includes multiplication
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_SCATTER, p_g)
get_macro_xs_c(p % material, tid, MG_GET_XS_SCATTER, p_g)
end if
end if
@ -1475,13 +1475,13 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
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)
end if
else
if (i_nuclide > 0) then
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g)
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_ABSORPTION, p_g)
else
score = material_xs % absorption * flux
end if
@ -1506,19 +1506,19 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
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)
else
score = score * &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g)
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)
end if
else
if (i_nuclide > 0) then
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) * &
score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_FISSION, p_g) * &
atom_density * flux
else
score = get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux
score = get_macro_xs_c(p % material, tid, MG_GET_XS_FISSION, p_g) * flux
end if
end if
@ -1544,12 +1544,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * atom_density * &
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)
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)
else
score = score * &
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)
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)
end if
else
! Skip any non-fission events
@ -1562,17 +1562,17 @@ contains
score = keff * p % wgt_bank * flux
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
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)
end if
end if
else
if (i_nuclide > 0) then
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) * &
score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_NU_FISSION, p_g) * &
atom_density * flux
else
score = get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) * flux
score = get_macro_xs_c(p % material, tid, MG_GET_XS_NU_FISSION, p_g) * flux
end if
end if
@ -1598,12 +1598,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * atom_density * &
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)
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)
else
score = score * &
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)
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)
end if
else
! Skip any non-fission events
@ -1617,17 +1617,17 @@ contains
/ real(p % n_bank, 8)) * flux
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
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)
end if
end if
else
if (i_nuclide > 0) then
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) * &
score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_PROMPT_NU_FISSION, p_g) * &
atom_density * flux
else
score = get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux
score = get_macro_xs_c(p % material, tid, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux
end if
end if
@ -1653,7 +1653,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, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (dg_filter > 0) then
select type(filt => filters(t % filter(dg_filter)) % obj)
@ -1669,12 +1669,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * &
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)
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)
else
score = score * &
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)
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)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1685,12 +1685,12 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * &
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)
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)
else
score = score * &
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)
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)
end if
end if
end if
@ -1720,8 +1720,8 @@ contains
if (i_nuclide > 0) then
score = score * atom_density * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
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)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1732,8 +1732,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, MG_GET_XS_FISSION, p_g) / &
get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g)
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)
end if
end if
end if
@ -1753,10 +1753,10 @@ contains
if (i_nuclide > 0) then
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
else
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d)
get_macro_xs_c(p % material, tid, 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 +1766,11 @@ contains
else
if (i_nuclide > 0) then
score = atom_density * flux * &
get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g)
get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g)
else
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g)
get_macro_xs_c(p % material, tid, MG_GET_XS_DELAYED_NU_FISSION, p_g)
end if
end if
end if
@ -1785,7 +1785,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, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (get_macro_xs_c(p % material, tid, MG_GET_XS_ABSORPTION, p_g) > ZERO) then
if (dg_filter > 0) then
select type(filt => filters(t % filter(dg_filter)) % obj)
@ -1801,14 +1801,14 @@ contains
score = p % absorb_wgt * flux
if (i_nuclide > 0) then
score = score * &
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)
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)
else
score = score * &
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)
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)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1826,14 +1826,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, 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)
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)
else
score = score + p % absorb_wgt * flux * &
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)
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)
end if
end do
end if
@ -1862,13 +1862,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, 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
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
else
score = score + keff * &
fission_bank(n_bank - p % n_bank + k) % wgt * &
get_macro_xs_c(p % material, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
get_macro_xs_c(p % material, tid, MG_GET_XS_DECAY_RATE, p_g, DG=d) * flux
end if
! if the delayed group filter is present, tally to corresponding
@ -1921,12 +1921,12 @@ contains
if (i_nuclide > 0) then
score = atom_density * 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_DELAYED_NU_FISSION, p_g, DG=d)
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)
else
score = flux * &
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, 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)
end if
call score_fission_delayed_dg(t, d_bin, score, score_index)
@ -1943,12 +1943,12 @@ contains
do d = 1, num_delayed_groups
if (i_nuclide > 0) then
score = score + atom_density * 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_DELAYED_NU_FISSION, p_g, DG=d)
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)
else
score = score + flux * &
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, 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)
end if
end do
end if
@ -1973,20 +1973,20 @@ contains
end if
if (i_nuclide > 0) then
score = score * atom_density * &
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)
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)
else
score = score * &
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)
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)
end if
else
if (i_nuclide > 0) then
score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) * &
score = get_nuclide_xs_c(i_nuclide, tid, MG_GET_XS_KAPPA_FISSION, p_g) * &
atom_density * flux
else
score = flux * &
get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g)
get_macro_xs_c(p % material, tid, MG_GET_XS_KAPPA_FISSION, p_g)
end if
end if
@ -2005,16 +2005,6 @@ contains
t % results(RESULT_VALUE, score_index, filter_index) + score
end do SCORE_LOOP
! Reset temporary Mgxs indices
call reset_macro_angle_index_c(p % material, last_mat_pol, last_mat_azi, &
last_mat_uvw);
if (i_nuclide > 0) then
call reset_nuclide_temperature_index_c(i_nuclide, last_nuc_temp)
call reset_nuclide_angle_index_c(i_nuclide, last_nuc_pol, last_nuc_azi, &
last_nuc_uvw)
end if
end subroutine score_general_mg
!===============================================================================

12
src/tracking.F90
View file

@ -2,6 +2,10 @@ 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
@ -48,6 +52,12 @@ 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
@ -114,7 +124,7 @@ contains
end if
else
! Get the MG data
call calculate_xs_c(p % material, p % g, p % sqrtkT, &
call calculate_xs_c(p % material, tid, p % g, p % sqrtkT, &
p % coord(p % n_coord) % uvw, material_xs % total, &
material_xs % absorption, material_xs % nu_fission)