From 267acf627feafa2883545fd2d3433654b44c302b Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Thu, 21 Jun 2018 20:45:10 -0400 Subject: [PATCH] resolving @paulromano comments, including intent(inout) items at the end, consistent usage of const qualifiers --- src/constants.h | 2 +- src/mgxs.cpp | 143 ++++++++++++++++++------------ src/mgxs.h | 119 ++++++++++++------------- src/mgxs_data.F90 | 27 +----- src/mgxs_interface.F90 | 36 +++----- src/mgxs_interface.cpp | 84 ++++++++---------- src/mgxs_interface.h | 49 +++++------ src/physics_mg.F90 | 20 +---- src/scattdata.cpp | 45 +++++----- src/scattdata.h | 129 ++++++++++++++------------- src/tallies/tally.F90 | 196 +++++++++++++++++++---------------------- src/tracking.F90 | 12 +-- src/xsdata.cpp | 28 +++--- src/xsdata.h | 28 +++--- 14 files changed, 421 insertions(+), 497 deletions(-) diff --git a/src/constants.h b/src/constants.h index 3b20e68f6..016a6b348 100644 --- a/src/constants.h +++ b/src/constants.h @@ -21,7 +21,7 @@ typedef std::vector int_1dvec; typedef std::vector > int_2dvec; typedef std::vector > > int_3dvec; -int constexpr MAX_SAMPLE {10000}; +constexpr int MAX_SAMPLE {10000}; constexpr std::array VERSION {0, 10, 0}; constexpr std::array VERSION_PARTICLE_RESTART {2, 0}; diff --git a/src/mgxs.cpp b/src/mgxs.cpp index d30417c28..748e7e949 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -12,13 +12,12 @@ std::vector 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& 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& 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& micros, double_1dvec& scalars, - int_1dvec& micro_ts, int this_t) +Mgxs::combine(const std::vector& 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 those_xs(micros.size()); @@ -415,36 +416,42 @@ Mgxs::combine(std::vector& 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))) { diff --git a/src/mgxs.h b/src/mgxs.h index c0cec5941..8ef174269 100644 --- a/src/mgxs.h +++ b/src/mgxs.h @@ -11,6 +11,10 @@ #include #include + #ifdef _OPENMP + # include + #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& micros, double_1dvec& scalars, - int_1dvec& micro_ts, int this_t); + combine(const std::vector& 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 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& 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& 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 diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index ef416125e..1b3ed3511 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -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 diff --git a/src/mgxs_interface.F90 b/src/mgxs_interface.F90 index 34b0fb611..2a579e930 100644 --- a/src/mgxs_interface.F90 +++ b/src/mgxs_interface.F90 @@ -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 diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 45234f454..a96f53008 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -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; } diff --git a/src/mgxs_interface.h b/src/mgxs_interface.h index 785f72b32..2cdd2889a 100644 --- a/src/mgxs_interface.h +++ b/src/mgxs_interface.h @@ -22,67 +22,58 @@ extern std::vector 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 \ No newline at end of file diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index d04d4cb85..a7a183295 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -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 diff --git a/src/scattdata.cpp b/src/scattdata.cpp index 952c05f64..ad6ae21a6 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -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& 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& 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& 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& 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(); diff --git a/src/scattdata.h b/src/scattdata.h index 69ce9e245..888f96052 100644 --- a/src/scattdata.h +++ b/src/scattdata.h @@ -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& those_scatts, + base_combine(int max_order, const std::vector& 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& 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& 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& 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& 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); }; //============================================================================== diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index f6839081d..fa4a644d1 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -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 diff --git a/src/tracking.F90 b/src/tracking.F90 index c92673984..8fe05e3c2 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -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) diff --git a/src/xsdata.cpp b/src/xsdata.cpp index b447ee0ae..98df2d67d 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -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")) { diff --git a/src/xsdata.h b/src/xsdata.h index adf80cad2..36aacc178 100644 --- a/src/xsdata.h +++ b/src/xsdata.h @@ -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. //!