From b1c73918a8f111f3f3b5e7b2e240b57403c1d2d8 Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Wed, 13 Jun 2018 20:22:51 -0400 Subject: [PATCH] It works! Replaced all mgxs_header functionality with the C++ version and a mgxs_interface module to act as the go-between the C++ and Fortran --- CMakeLists.txt | 5 +- src/api.F90 | 2 - src/cmfd_input.F90 | 4 +- src/constants.F90 | 19 ++ src/constants.h | 17 ++ src/input_xml.F90 | 10 +- src/mgxs.cpp | 265 +++++++----------- src/mgxs.h | 47 +--- src/mgxs_data.F90 | 414 +++++++++++++--------------- src/mgxs_header.F90 | 1 - src/mgxs_interface.F90 | 201 ++++++++++++++ src/mgxs_interface.cpp | 247 +++++++++++++++++ src/mgxs_interface.h | 66 +++++ src/output.F90 | 6 +- src/particle_restart.F90 | 2 +- src/physics_mg.F90 | 33 +-- src/scattdata.cpp | 46 ++-- src/scattdata.h | 2 +- src/simulation.F90 | 2 +- src/source.F90 | 2 +- src/state_point.F90 | 11 +- src/string_functions.cpp | 30 ++ src/string_functions.h | 31 +-- src/summary.F90 | 24 +- src/tallies/tally.F90 | 259 +++++++++-------- src/tallies/tally_filter_energy.F90 | 2 +- src/tracking.F90 | 21 +- src/xsdata.cpp | 15 + src/xsdata.h | 1 + 29 files changed, 1111 insertions(+), 674 deletions(-) create mode 100644 src/mgxs_interface.F90 create mode 100644 src/mgxs_interface.cpp create mode 100644 src/mgxs_interface.h create mode 100644 src/string_functions.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 362e98a830..a13f30b575 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -363,7 +363,7 @@ set(LIBOPENMC_FORTRAN_SRC src/mesh_header.F90 src/message_passing.F90 src/mgxs_data.F90 - src/mgxs_header.F90 + src/mgxs_interface.F90 src/multipole_header.F90 src/nuclide_header.F90 src/output.F90 @@ -380,7 +380,6 @@ set(LIBOPENMC_FORTRAN_SRC src/reaction_header.F90 src/relaxng src/sab_header.F90 - src/scattdata_header.F90 src/secondary_correlated.F90 src/secondary_kalbach.F90 src/secondary_nbody.F90 @@ -439,11 +438,13 @@ set(LIBOPENMC_CXX_SRC src/math_functions.cpp src/message_passing.cpp src/mgxs.cpp + src/mgxs_interface.cpp src/plot.cpp src/random_lcg.cpp src/scattdata.cpp src/simulation.cpp src/state_point.cpp + src/string_functions.cpp src/surface.cpp src/xml_interface.cpp src/xsdata.cpp diff --git a/src/api.F90 b/src/api.F90 index da50007d68..f32de6bebb 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -309,7 +309,6 @@ contains subroutine free_memory() use cmfd_header - use mgxs_header use plot_header use sab_header use settings @@ -327,7 +326,6 @@ contains call free_memory_simulation() call free_memory_nuclide() call free_memory_settings() - call free_memory_mgxs() call free_memory_sab() call free_memory_source() call free_memory_mesh() diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 549cff4191..1e4b8f3742 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -3,8 +3,8 @@ module cmfd_input use, intrinsic :: ISO_C_BINDING use cmfd_header - use mesh_header, only: mesh_dict - use mgxs_header, only: energy_bins + use mesh_header, only: mesh_dict + use mgxs_interface, only: energy_bins, num_energy_groups use tally use tally_header use timer_header diff --git a/src/constants.F90 b/src/constants.F90 index d1893bf9e1..430cc27151 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -402,6 +402,25 @@ module constants DIFF_NUCLIDE_DENSITY = 2, & DIFF_TEMPERATURE = 3 + + ! Mgxs::get_xs enumerated types + integer(C_INT), parameter :: & + MG_GET_XS_TOTAL = 0, & + MG_GET_XS_ABSORPTION = 1, & + MG_GET_XS_INVERSE_VELOCITY = 2, & + MG_GET_XS_DECAY_RATE = 3, & + MG_GET_XS_SCATTER = 4, & + MG_GET_XS_SCATTER_MULT = 5, & + MG_GET_XS_SCATTER_FMU_MULT = 6, & + MG_GET_XS_SCATTER_FMU = 7, & + MG_GET_XS_FISSION = 8, & + MG_GET_XS_KAPPA_FISSION = 9, & + MG_GET_XS_PROMPT_NU_FISSION = 10, & + MG_GET_XS_DELAYED_NU_FISSION = 11, & + MG_GET_XS_NU_FISSION = 12, & + MG_GET_XS_CHI_PROMPT = 13, & + MG_GET_XS_CHI_DELAYED = 14 + ! ============================================================================ ! RANDOM NUMBER STREAM CONSTANTS diff --git a/src/constants.h b/src/constants.h index 8290aa8c9c..3ddadb5a10 100644 --- a/src/constants.h +++ b/src/constants.h @@ -63,6 +63,23 @@ constexpr double PI {3.1415926535898}; const double SQRT_PI {std::sqrt(PI)}; +// Mgxs::get_xs enumerated types +constexpr int MG_GET_XS_TOTAL {0}; +constexpr int MG_GET_XS_ABSORPTION {1}; +constexpr int MG_GET_XS_INVERSE_VELOCITY {2}; +constexpr int MG_GET_XS_DECAY_RATE {3}; +constexpr int MG_GET_XS_SCATTER {4}; +constexpr int MG_GET_XS_SCATTER_MULT {5}; +constexpr int MG_GET_XS_SCATTER_FMU_MULT {6}; +constexpr int MG_GET_XS_SCATTER_FMU {7}; +constexpr int MG_GET_XS_FISSION {8}; +constexpr int MG_GET_XS_KAPPA_FISSION {9}; +constexpr int MG_GET_XS_PROMPT_NU_FISSION {10}; +constexpr int MG_GET_XS_DELAYED_NU_FISSION {11}; +constexpr int MG_GET_XS_NU_FISSION {12}; +constexpr int MG_GET_XS_CHI_PROMPT {13}; +constexpr int MG_GET_XS_CHI_DELAYED {14}; + } // namespace openmc diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 507bdf8a9c..189ca7c212 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -19,8 +19,8 @@ module input_xml use material_header use mesh_header use message_passing - use mgxs_data, only: create_macro_xs, read_mgxs, read_mgxs2, create_macro_xs2 - use mgxs_header + use mgxs_data, only: create_macro_xs, read_mgxs + use mgxs_interface use nuclide_header use output, only: title, header, print_plot use plot_header @@ -84,9 +84,7 @@ contains else ! Create material macroscopic data for MGXS call read_mgxs() - call read_mgxs2() call create_macro_xs() - call create_macro_xs2() end if call time_read_xs % stop() end if @@ -3856,7 +3854,7 @@ contains if (run_CE) then awr = nuclides(mat % nuclide(j)) % awr else - awr = nuclides_MG(mat % nuclide(j)) % obj % awr + awr = get_awr_c(mat % nuclide(j)) end if ! if given weight percent, convert all values so that they are divided @@ -3881,7 +3879,7 @@ contains if (run_CE) then awr = nuclides(mat % nuclide(j)) % awr else - awr = nuclides_MG(mat % nuclide(j)) % obj % awr + awr = get_awr_c(mat % nuclide(j)) end if x = mat % atom_density(j) sum_percent = sum_percent + x*awr diff --git a/src/mgxs.cpp b/src/mgxs.cpp index ec711e1307..4f92a64eee 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -2,6 +2,11 @@ namespace openmc { +// Storage for the MGXS data +std::vector nuclides_MG; +std::vector macro_xs; + + //============================================================================== // Mgxs base-class methods //============================================================================== @@ -390,112 +395,149 @@ void Mgxs::combine(std::vector& micros, double_1dvec& scalars, } -double Mgxs::get_xs(const char* xstype, const int gin, int* gout, double* mu, +double Mgxs::get_xs(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; - if (std::strcmp(xstype, "total")) { + switch(xstype) { + case MG_GET_XS_TOTAL: val = xs[index_temp].total[index_pol][index_azi][gin]; - } else if (std::strcmp(xstype, "absorption")) { + break; + case MG_GET_XS_ABSORPTION: val = xs[index_temp].absorption[index_pol][index_azi][gin]; - } else if (std::strcmp(xstype, "inverse-velocity")) { + break; + case MG_GET_XS_INVERSE_VELOCITY: val = xs[index_temp].inverse_velocity[index_pol][index_azi][gin]; - } else if (std::strcmp(xstype, "decay rate")) { + break; + case MG_GET_XS_DECAY_RATE: if (dg != nullptr) { val = xs[index_temp].decay_rate[index_pol][index_azi][*dg + 1]; } else { val = xs[index_temp].decay_rate[index_pol][index_azi][0]; } - } else if ((std::strcmp(xstype, "scatter")) || - (std::strcmp(xstype, "scatter/mult")) || - (std::strcmp(xstype, "scatter*f_mu/mult")) || - (std::strcmp(xstype, "scatter*f_mu"))) { + 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); - } else if (fissionable && std::strcmp(xstype, "fission")) { - val = xs[index_temp].fission[index_pol][index_azi][gin]; - } else if (fissionable && std::strcmp(xstype, "kappa-fission")) { - val = xs[index_temp].kappa_fission[index_pol][index_azi][gin]; - } else if (fissionable && std::strcmp(xstype, "prompt-nu-fission")) { - val = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin]; - } else if (fissionable && std::strcmp(xstype, "delayed-nu-fission")) { - if (dg != nullptr) { - val = xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][*dg]; + break; + case MG_GET_XS_FISSION: + if (fissionable) { + val = xs[index_temp].fission[index_pol][index_azi][gin]; } else { val = 0.; - for (auto& num : xs[index_temp].delayed_nu_fission[index_pol] - [index_azi][gin]) { - val += num; - } } - } else if (fissionable && std::strcmp(xstype, "nu-fission")) { - val = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin]; - for (auto& num : xs[index_temp].delayed_nu_fission[index_pol] - [index_azi][gin]) { - val += num; - } - } else if (fissionable && std::strcmp(xstype, "chi-prompt")) { - if (gout != nullptr) { - val = xs[index_temp].chi_prompt[index_pol][index_azi][gin][*gout]; + break; + case MG_GET_XS_KAPPA_FISSION: + if (fissionable) { + val = xs[index_temp].kappa_fission[index_pol][index_azi][gin]; } else { - // provide an outgoing group-wise sum val = 0.; - for (auto& num : xs[index_temp].chi_prompt[index_pol][index_azi][gin]) { - val += num; - } } - } else if (fissionable && std::strcmp(xstype, "chi-delayed")) { - if (gout != nullptr) { + break; + case MG_GET_XS_PROMPT_NU_FISSION: + if (fissionable) { + val = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin]; + } else { + val = 0.; + } + break; + case MG_GET_XS_DELAYED_NU_FISSION: + if (fissionable) { if (dg != nullptr) { - val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][*dg]; + val = xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][*dg]; } else { - val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][0]; + val = 0.; + for (auto& num : xs[index_temp].delayed_nu_fission[index_pol] + [index_azi][gin]) { + val += num; + } } } else { - if (dg != nullptr) { + val = 0.; + } + break; + case MG_GET_XS_NU_FISSION: + if (fissionable) { + val = xs[index_temp].nu_fission[index_pol][index_azi][gin]; + } else { + val = 0.; + } + break; + case MG_GET_XS_CHI_PROMPT: + if (fissionable) { + if (gout != nullptr) { + val = xs[index_temp].chi_prompt[index_pol][index_azi][gin][*gout]; + } else { + // provide an outgoing group-wise sum 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 (auto& num : xs[index_temp].chi_prompt[index_pol][index_azi][gin]) { + val += num; + } + } + } else { + val = 0.; + } + break; + case MG_GET_XS_CHI_DELAYED: + if (fissionable) { + if (gout != nullptr) { + if (dg != nullptr) { + val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][*dg]; + } else { + val = xs[index_temp].chi_delayed[index_pol][index_azi][gin][*gout][0]; } } 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]) { - val += num; + 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]; + } + } 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]) { + val += num; + } } } } + } else { + val = 0.; } - } else { + break; + default: val = 0.; } return val; } -void Mgxs::sample_fission_energy(const int gin, const double nu_fission, - int& dg, int& gout) +void Mgxs::sample_fission_energy(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]; + // Find the probability of having a prompt neutron double prob_prompt = - xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin] / - nu_fission; + xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin]; // sample random numbers - double xi_pd = prn(); + double xi_pd = prn() * nu_fission; double xi_gout = prn(); // Select whether the neutron is prompt or delayed if (xi_pd <= prob_prompt) { // the neutron is prompt - // set the delayed group for the particle to be 0, indicating prompt - dg = 0; + // set the delayed group for the particle to be -1, indicating prompt + dg = -1; // sample the outgoing energy group gout = 0; @@ -514,12 +556,11 @@ void Mgxs::sample_fission_energy(const int gin, const double nu_fission, while (xi_pd >= prob_prompt) { dg++; prob_prompt += - xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][dg] / - nu_fission; + xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin][dg]; } // adjust dg in case of round-off error - dg = std::min(dg, num_delayed_groups); + dg = std::min(dg, num_delayed_groups - 1); // sample the outgoing energy group gout = 0; @@ -552,11 +593,7 @@ void Mgxs::calculate_xs(const int gin, const double sqrtkT, const double uvw[3], abs_xs = xs[index_temp].absorption[index_pol][index_azi][gin]; if (fissionable) { - // nu-fission is made up of the prompt and all the delayed nu_fission data - nu_fiss_xs = xs[index_temp].prompt_nu_fission[index_pol][index_azi][gin]; - for (auto& val : xs[index_temp].delayed_nu_fission[index_pol][index_azi][gin]) { - nu_fiss_xs += val; - } + nu_fiss_xs = xs[index_temp].nu_fission[index_pol][index_azi][gin]; } else { nu_fiss_xs = 0.; } @@ -580,7 +617,7 @@ bool Mgxs::equiv(const Mgxs& that) } -inline void Mgxs::set_temperature_index(const double sqrtkT) +void Mgxs::set_temperature_index(const double sqrtkT) { // See if we need to find the new index if (sqrtkT != last_sqrtkT) { @@ -600,7 +637,7 @@ inline void Mgxs::set_temperature_index(const double sqrtkT) } -inline void Mgxs::set_angle_index(const double uvw[3]) +void Mgxs::set_angle_index(const double uvw[3]) { // See if we need to find the new index if ((uvw[0] != last_uvw[0]) || (uvw[1] != last_uvw[1]) || @@ -622,98 +659,4 @@ inline void Mgxs::set_angle_index(const double uvw[3]) } } -//============================================================================== -// Mgxs data loading interface methods -//============================================================================== - -void add_mgxs(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) -{ - //!! mgxs_data.F90 will be modified to just create the list of names - //!! in the order needed - // Convert temps to a vector for the from_hdf5 function - double_1dvec temperature; - temperature.assign(temps, temps + n_temps); - - // TODO: C++ replacement for write_message - // write_message("Loading " + std::string(names[i]) + " data...", 6); - - // Check to make sure cross section set exists in the library - hid_t xs_grp; - if (object_exists(file_id, name)) { - xs_grp = open_group(file_id, name); - } else { - fatal_error("Data for " + std::string(name) + " does not exist in " - + "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); - - nuclides_MG.push_back(mg); -} - - -bool query_fissionable(const int n_nuclides, const int i_nuclides[]) -{ - bool result = false; - for (int n = 0; n < n_nuclides; n++) { - if (nuclides_MG[i_nuclides[n] - 1].fissionable) result = true; - } - return result; -} - - -void create_macro_xs(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) -{ - Mgxs macro; - if (n_temps > 0) { - // // Convert temps to a vector - double_1dvec temperature; - temperature.assign(temps, temps + n_temps); - - // Convert atom_densities to a vector - double_1dvec atom_densities_vec; - atom_densities_vec.assign(atom_densities, atom_densities + n_nuclides); - - // Build array of pointers to nuclides_MG's Mgxs objects needed for this - // material - std::vector mgxs_ptr(n_nuclides); - for (int n = 0; n < n_nuclides; n++) { - mgxs_ptr[n] = &nuclides_MG[i_nuclides[n] - 1]; - } - - macro.build_macro(mat_name, temperature, mgxs_ptr, atom_densities_vec, - method, tolerance); - } - macro_xs.push_back(macro); -} - -//============================================================================== -// Mgxs tracking/transport/tallying interface methods -//============================================================================== - -void calculate_xs(const int i_mat, 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, - nu_fiss_xs); -} - -void scatter(const int i_mat, 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); - gout = gout_c + 1; - - rotate_angle_c(uvw, mu, nullptr); -} - } // namespace openmc \ No newline at end of file diff --git a/src/mgxs.h b/src/mgxs.h index 14c4e4e89c..2fc3e2bec2 100644 --- a/src/mgxs.h +++ b/src/mgxs.h @@ -29,29 +29,31 @@ namespace openmc { class Mgxs { private: - std::string name; // name of dataset, e.g., UO2 - double awr; // atomic weight ratio double_1dvec kTs; // temperature in eV (k * T) 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 - int index_temp; // cache of temperature index double last_sqrtkT; // cache of the temperature corresponding to index_temp std::vector xs; // Cross section data int n_pol; int n_azi; - int index_pol; // cache fof the angle indices - int index_azi; double_1dvec polar; double_1dvec azimuthal; - double last_uvw[3]; 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); + bool equiv(const Mgxs& that); public: + std::string name; // name of dataset, e.g., UO2 + 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 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, @@ -66,40 +68,15 @@ class Mgxs { double_1dvec& temperature, int& method, double tolerance, int max_order, bool legendre_to_tabular, int legendre_to_tabular_points); - double get_xs(const char* xstype, const int gin, int* gout, double* mu, + double get_xs(const int xstype, const int gin, int* gout, double* mu, int* dg); - void sample_fission_energy(const int gin, const double nu_fission, - int& dg, int& gout); + 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); - bool equiv(const Mgxs& that); - inline void set_temperature_index(const double sqrtkT); - inline void set_angle_index(const double uvw[3]); + void set_temperature_index(const double sqrtkT); + void set_angle_index(const double uvw[3]); }; -extern "C" void add_mgxs(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" bool query_fissionable(const int n_nuclides, const int i_nuclides[]); - -extern "C" void create_macro_xs(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); - -extern "C" void calculate_xs(const int i_mat, const int gin, - const double sqrtkT, const double uvw[3], double& total_xs, - double& abs_xs, double& nu_fiss_xs); - -extern "C" void scatter(const int i_mat, const int gin, int& gout, double& mu, - double& wgt, double uvw[3]); - - -// Storage for the MGXS data -std::vector nuclides_MG; -std::vector macro_xs; - } // namespace openmc #endif // MGXS_H \ No newline at end of file diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index 6977378801..0278c96eb6 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -9,7 +9,7 @@ module mgxs_data use geometry_header, only: get_temperatures, cells use hdf5_interface use material_header, only: Material, materials, n_materials - use mgxs_header + use mgxs_interface use nuclide_header, only: n_nuclides use set_header, only: SetChar use settings @@ -17,50 +17,6 @@ module mgxs_data use string, only: to_lower implicit none - 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, name='add_mgxs') - use ISO_C_BINDING - import HID_T - implicit none - integer(HID_T), value, intent(in) :: file_id - character(kind=C_CHAR),intent(in) :: name(*) - integer(C_INT), value, intent(in) :: energy_groups - 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 - end subroutine add_mgxs_c - - function query_fissionable_c(n_nuclides, i_nuclides) result(result) & - bind(C, name='query_fissionable') - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: n_nuclides - integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides) - logical(C_BOOL) :: result - end function query_fissionable_c - - subroutine create_macro_xs_c(name, n_nuclides, i_nuclides, n_temps, temps, & - atom_densities, method, tolerance) bind(C, name='create_macro_xs') - use ISO_C_BINDING - implicit none - character(kind=C_CHAR),intent(in) :: name(*) - integer(C_INT), value, intent(in) :: n_nuclides - integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides) - 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 - end subroutine create_macro_xs_c - end interface - contains !=============================================================================== @@ -68,150 +24,150 @@ contains ! nuclides and sab_tables arrays !=============================================================================== + ! subroutine read_mgxs() + ! integer :: i ! index in materials array + ! integer :: j ! index over nuclides in material + ! integer :: i_nuclide ! index in nuclides array + ! character(20) :: name ! name of library to load + ! integer :: representation ! Data representation + ! character(MAX_LINE_LEN) :: temp_str + ! type(Material), pointer :: mat + ! type(SetChar) :: already_read + ! integer(HID_T) :: file_id + ! integer(HID_T) :: xsdata_group + ! logical :: file_exists + ! type(VectorReal), allocatable :: temps(:) + ! character(MAX_WORD_LEN) :: word + ! integer, allocatable :: array(:) + + ! ! Check if MGXS Library exists + ! inquire(FILE=path_cross_sections, EXIST=file_exists) + ! if (.not. file_exists) then + + ! ! Could not find MGXS Library file + ! call fatal_error("Cross sections HDF5 file '" & + ! // trim(path_cross_sections) // "' does not exist!") + ! end if + + ! call write_message("Loading cross section data...", 5) + + ! ! Get temperatures + ! call get_temperatures(temps) + + ! ! Open file for reading + ! file_id = file_open(path_cross_sections, 'r', parallel=.true.) + + ! ! Read filetype + ! call read_attribute(word, file_id, "filetype") + ! if (word /= 'mgxs') then + ! call fatal_error("Provided MGXS Library is not a MGXS Library file.") + ! end if + + ! ! Read revision number for the MGXS Library file and make sure it matches + ! ! with the current version + ! call read_attribute(array, file_id, "version") + ! if (any(array /= VERSION_MGXS_LIBRARY)) then + ! call fatal_error("MGXS Library file version does not match current & + ! &version supported by OpenMC.") + ! end if + + ! ! allocate arrays for MGXS storage and cross section cache + ! allocate(nuclides_MG(n_nuclides)) + + ! ! ========================================================================== + ! ! READ ALL MGXS CROSS SECTION TABLES + + ! ! Loop over all files + ! MATERIAL_LOOP: do i = 1, n_materials + ! mat => materials(i) + + ! NUCLIDE_LOOP: do j = 1, mat % n_nuclides + ! name = mat % names(j) + + ! if (.not. already_read % contains(name)) then + ! i_nuclide = mat % nuclide(j) + + ! call write_message("Loading " // trim(name) // " data...", 6) + + ! ! Check to make sure cross section set exists in the library + ! if (object_exists(file_id, trim(name))) then + ! xsdata_group = open_group(file_id, trim(name)) + ! else + ! call fatal_error("Data for '" // trim(name) // "' does not exist in "& + ! &// trim(path_cross_sections)) + ! end if + + ! ! First find out the data representation + ! if (attribute_exists(xsdata_group, "representation")) then + + ! call read_attribute(temp_str, xsdata_group, "representation") + + ! if (trim(temp_str) == 'isotropic') then + ! representation = MGXS_ISOTROPIC + ! else if (trim(temp_str) == 'angle') then + ! representation = MGXS_ANGLE + ! else + ! call fatal_error("Invalid Data Representation!") + ! end if + ! else + ! ! Default to isotropic representation + ! representation = MGXS_ISOTROPIC + ! end if + + ! ! Now allocate accordingly + ! select case(representation) + + ! case(MGXS_ISOTROPIC) + ! allocate(MgxsIso :: nuclides_MG(i_nuclide) % obj) + + ! case(MGXS_ANGLE) + ! allocate(MgxsAngle :: nuclides_MG(i_nuclide) % obj) + + ! end select + + ! ! Now read in the data specific to the type we just declared + ! call nuclides_MG(i_nuclide) % obj % from_hdf5(xsdata_group, & + ! num_energy_groups, num_delayed_groups, temps(i_nuclide), & + ! temperature_method, temperature_tolerance, max_order, & + ! legendre_to_tabular, legendre_to_tabular_points) + + ! ! Add name to dictionary + ! call already_read % add(name) + + ! call close_group(xsdata_group) + + ! end if + ! end do NUCLIDE_LOOP + ! end do MATERIAL_LOOP + + ! ! Avoid some valgrind leak errors + ! call already_read % clear() + + ! ! Loop around material + ! MATERIAL_LOOP3: do i = 1, n_materials + + ! ! Get material + ! mat => materials(i) + + ! ! Loop around nuclides in material + ! NUCLIDE_LOOP2: do j = 1, mat % n_nuclides + + ! ! Is this fissionable? + ! if (nuclides_MG(mat % nuclide(j)) % obj % fissionable) then + ! mat % fissionable = .true. + ! end if + ! if (mat % fissionable) then + ! exit NUCLIDE_LOOP2 + ! end if + + ! end do NUCLIDE_LOOP2 + ! end do MATERIAL_LOOP3 + + ! call file_close(file_id) + + ! end subroutine read_mgxs + subroutine read_mgxs() - integer :: i ! index in materials array - integer :: j ! index over nuclides in material - integer :: i_nuclide ! index in nuclides array - character(20) :: name ! name of library to load - integer :: representation ! Data representation - character(MAX_LINE_LEN) :: temp_str - type(Material), pointer :: mat - type(SetChar) :: already_read - integer(HID_T) :: file_id - integer(HID_T) :: xsdata_group - logical :: file_exists - type(VectorReal), allocatable :: temps(:) - character(MAX_WORD_LEN) :: word - integer, allocatable :: array(:) - - ! Check if MGXS Library exists - inquire(FILE=path_cross_sections, EXIST=file_exists) - if (.not. file_exists) then - - ! Could not find MGXS Library file - call fatal_error("Cross sections HDF5 file '" & - // trim(path_cross_sections) // "' does not exist!") - end if - - call write_message("Loading cross section data...", 5) - - ! Get temperatures - call get_temperatures(temps) - - ! Open file for reading - file_id = file_open(path_cross_sections, 'r', parallel=.true.) - - ! Read filetype - call read_attribute(word, file_id, "filetype") - if (word /= 'mgxs') then - call fatal_error("Provided MGXS Library is not a MGXS Library file.") - end if - - ! Read revision number for the MGXS Library file and make sure it matches - ! with the current version - call read_attribute(array, file_id, "version") - if (any(array /= VERSION_MGXS_LIBRARY)) then - call fatal_error("MGXS Library file version does not match current & - &version supported by OpenMC.") - end if - - ! allocate arrays for MGXS storage and cross section cache - allocate(nuclides_MG(n_nuclides)) - - ! ========================================================================== - ! READ ALL MGXS CROSS SECTION TABLES - - ! Loop over all files - MATERIAL_LOOP: do i = 1, n_materials - mat => materials(i) - - NUCLIDE_LOOP: do j = 1, mat % n_nuclides - name = mat % names(j) - - if (.not. already_read % contains(name)) then - i_nuclide = mat % nuclide(j) - - call write_message("Loading " // trim(name) // " data...", 6) - - ! Check to make sure cross section set exists in the library - if (object_exists(file_id, trim(name))) then - xsdata_group = open_group(file_id, trim(name)) - else - call fatal_error("Data for '" // trim(name) // "' does not exist in "& - &// trim(path_cross_sections)) - end if - - ! First find out the data representation - if (attribute_exists(xsdata_group, "representation")) then - - call read_attribute(temp_str, xsdata_group, "representation") - - if (trim(temp_str) == 'isotropic') then - representation = MGXS_ISOTROPIC - else if (trim(temp_str) == 'angle') then - representation = MGXS_ANGLE - else - call fatal_error("Invalid Data Representation!") - end if - else - ! Default to isotropic representation - representation = MGXS_ISOTROPIC - end if - - ! Now allocate accordingly - select case(representation) - - case(MGXS_ISOTROPIC) - allocate(MgxsIso :: nuclides_MG(i_nuclide) % obj) - - case(MGXS_ANGLE) - allocate(MgxsAngle :: nuclides_MG(i_nuclide) % obj) - - end select - - ! Now read in the data specific to the type we just declared - call nuclides_MG(i_nuclide) % obj % from_hdf5(xsdata_group, & - num_energy_groups, num_delayed_groups, temps(i_nuclide), & - temperature_method, temperature_tolerance, max_order, & - legendre_to_tabular, legendre_to_tabular_points) - - ! Add name to dictionary - call already_read % add(name) - - call close_group(xsdata_group) - - end if - end do NUCLIDE_LOOP - end do MATERIAL_LOOP - - ! Avoid some valgrind leak errors - call already_read % clear() - - ! Loop around material - MATERIAL_LOOP3: do i = 1, n_materials - - ! Get material - mat => materials(i) - - ! Loop around nuclides in material - NUCLIDE_LOOP2: do j = 1, mat % n_nuclides - - ! Is this fissionable? - if (nuclides_MG(mat % nuclide(j)) % obj % fissionable) then - mat % fissionable = .true. - end if - if (mat % fissionable) then - exit NUCLIDE_LOOP2 - end if - - end do NUCLIDE_LOOP2 - end do MATERIAL_LOOP3 - - call file_close(file_id) - - end subroutine read_mgxs - - subroutine read_mgxs2() integer :: i ! index in materials array integer :: j ! index over nuclides in material integer :: i_nuclide ! index in nuclides array @@ -286,55 +242,55 @@ contains ! Avoid some valgrind leak errors call already_read % clear() - end subroutine read_mgxs2 + end subroutine read_mgxs !=============================================================================== ! CREATE_MACRO_XS generates the macroscopic xs from the microscopic input data !=============================================================================== + ! subroutine create_macro_xs() + ! integer :: i_mat ! index in materials array + ! type(Material), pointer :: mat ! current material + ! type(VectorReal), allocatable :: kTs(:) + + ! allocate(macro_xs(n_materials)) + + ! ! Get temperatures to read for each material + ! call get_mat_kTs(kTs) + + ! ! Force all nuclides in a material to be the same representation. + ! ! Therefore type(nuclides(mat % nuclide(1)) % obj) dictates type(macroxs). + ! ! At the same time, we will find the scattering type, as that will dictate + ! ! how we allocate the scatter object within macroxs.allocate(macro_xs(n_materials)) + ! do i_mat = 1, n_materials + + ! ! Get the material + ! mat => materials(i_mat) + + ! ! Get the scattering type for the first nuclide + ! select type(nuc => nuclides_MG(mat % nuclide(1)) % obj) + ! type is (MgxsIso) + ! allocate(MgxsIso :: macro_xs(i_mat) % obj) + ! type is (MgxsAngle) + ! allocate(MgxsAngle :: macro_xs(i_mat) % obj) + ! end select + + ! ! Do not read materials which we do not actually use in the problem to + ! ! reduce storage + ! if (allocated(kTs(i_mat) % data)) then + ! call macro_xs(i_mat) % obj % combine(kTs(i_mat), mat, nuclides_MG, & + ! num_energy_groups, num_delayed_groups, max_order, & + ! temperature_tolerance, temperature_method) + ! end if + ! end do + + ! end subroutine create_macro_xs + + subroutine create_macro_xs() integer :: i_mat ! index in materials array type(Material), pointer :: mat ! current material type(VectorReal), allocatable :: kTs(:) - - allocate(macro_xs(n_materials)) - - ! Get temperatures to read for each material - call get_mat_kTs(kTs) - - ! Force all nuclides in a material to be the same representation. - ! Therefore type(nuclides(mat % nuclide(1)) % obj) dictates type(macroxs). - ! At the same time, we will find the scattering type, as that will dictate - ! how we allocate the scatter object within macroxs.allocate(macro_xs(n_materials)) - do i_mat = 1, n_materials - - ! Get the material - mat => materials(i_mat) - - ! Get the scattering type for the first nuclide - select type(nuc => nuclides_MG(mat % nuclide(1)) % obj) - type is (MgxsIso) - allocate(MgxsIso :: macro_xs(i_mat) % obj) - type is (MgxsAngle) - allocate(MgxsAngle :: macro_xs(i_mat) % obj) - end select - - ! Do not read materials which we do not actually use in the problem to - ! reduce storage - if (allocated(kTs(i_mat) % data)) then - call macro_xs(i_mat) % obj % combine(kTs(i_mat), mat, nuclides_MG, & - num_energy_groups, num_delayed_groups, max_order, & - temperature_tolerance, temperature_method) - end if - end do - - end subroutine create_macro_xs - - - subroutine create_macro_xs2() - integer :: i_mat ! index in materials array - type(Material), pointer :: mat ! current material - type(VectorReal), allocatable :: kTs(:) character(MAX_WORD_LEN) :: name ! name of material ! Get temperatures to read for each material @@ -360,7 +316,7 @@ contains end if end do - end subroutine create_macro_xs2 + end subroutine create_macro_xs !=============================================================================== diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index be3e4e9d4d..09db7217ea 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -11,7 +11,6 @@ module mgxs_header use math, only: evaluate_legendre use nuclide_header, only: MaterialMacroXS use random_lcg, only: prn - use scattdata_header use string use stl_vector, only: VectorInt, VectorReal diff --git a/src/mgxs_interface.F90 b/src/mgxs_interface.F90 new file mode 100644 index 0000000000..7404af686d --- /dev/null +++ b/src/mgxs_interface.F90 @@ -0,0 +1,201 @@ +module mgxs_interface + + use, intrinsic :: ISO_C_BINDING + + use hdf5_interface + + implicit none + + 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) + use ISO_C_BINDING + import HID_T + implicit none + integer(HID_T), value, intent(in) :: file_id + character(kind=C_CHAR),intent(in) :: name(*) + integer(C_INT), value, intent(in) :: energy_groups + 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 + end subroutine add_mgxs_c + + function query_fissionable_c(n_nuclides, i_nuclides) result(result) bind(C) + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: n_nuclides + integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides) + logical(C_BOOL) :: result + end function query_fissionable_c + + subroutine create_macro_xs_c(name, n_nuclides, i_nuclides, n_temps, temps, & + atom_densities, method, tolerance) bind(C) + use ISO_C_BINDING + implicit none + character(kind=C_CHAR),intent(in) :: name(*) + integer(C_INT), value, intent(in) :: n_nuclides + integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides) + 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 + end subroutine create_macro_xs_c + + 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) :: gin + real(C_DOUBLE), value, intent(in) :: sqrtkT + real(C_DOUBLE), intent(in) :: uvw(1:3) + real(C_DOUBLE), intent(inout) :: total_xs + real(C_DOUBLE), intent(inout) :: abs_xs + 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) + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: i_mat + integer(C_INT), value, intent(in) :: gin + integer(C_INT), intent(inout) :: gout + real(C_DOUBLE), intent(inout) :: mu + real(C_DOUBLE), intent(inout) :: wgt + 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) + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: i_mat + integer(C_INT), value, intent(in) :: gin + integer(C_INT), intent(inout) :: dg + integer(C_INT), intent(inout) :: gout + end subroutine sample_fission_energy_c + + subroutine get_name_c(index, name_len, name) bind(C) + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: index + integer(C_INT), value, intent(in) :: name_len + character(kind=C_CHAR), intent(inout) :: name(name_len) + end subroutine get_name_c + + function get_awr_c(index) result(awr) bind(C) + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: index + real(C_DOUBLE) :: awr + end function get_awr_c + + 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) :: xstype + integer(C_INT), value, intent(in) :: gin + integer(C_INT), optional, intent(in) :: gout + real(C_DOUBLE), optional, intent(in) :: mu + integer(C_INT), optional, intent(in) :: dg + real(C_DOUBLE) :: val + end function get_nuclide_xs_c + + 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) :: xstype + integer(C_INT), value, intent(in) :: gin + integer(C_INT), optional, intent(in) :: gout + real(C_DOUBLE), optional, intent(in) :: mu + integer(C_INT), optional, intent(in) :: dg + 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) + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: index + 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) + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: index + 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) + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: index + 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 interface + + ! Number of energy groups + integer(C_INT) :: num_energy_groups + + ! Number of delayed groups + integer(C_INT) :: num_delayed_groups + + ! Energy group structure with decreasing energy + real(8), allocatable :: energy_bins(:) + + ! Midpoint of the energy group structure + real(8), allocatable :: energy_bin_avg(:) + + ! Energy group structure with increasing energy + real(8), allocatable :: rev_energy_bins(:) + +end module mgxs_interface \ No newline at end of file diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp new file mode 100644 index 0000000000..60afc28487 --- /dev/null +++ b/src/mgxs_interface.cpp @@ -0,0 +1,247 @@ +#include "mgxs_interface.h" + +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) +{ + //!! mgxs_data.F90 will be modified to just create the list of names + //!! in the order needed + // Convert temps to a vector for the from_hdf5 function + double_1dvec temperature; + temperature.assign(temps, temps + n_temps); + + // TODO: C++ replacement for write_message + // write_message("Loading " + std::string(names[i]) + " data...", 6); + + // Check to make sure cross section set exists in the library + hid_t xs_grp; + if (object_exists(file_id, name)) { + xs_grp = open_group(file_id, name); + } else { + fatal_error("Data for " + std::string(name) + " does not exist in " + + "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); + + nuclides_MG.push_back(mg); +} + + +bool query_fissionable_c(const int n_nuclides, const int i_nuclides[]) +{ + bool result = false; + for (int n = 0; n < n_nuclides; n++) { + if (nuclides_MG[i_nuclides[n] - 1].fissionable) result = true; + } + return result; +} + + +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) +{ + Mgxs macro; + if (n_temps > 0) { + // // Convert temps to a vector + double_1dvec temperature; + temperature.assign(temps, temps + n_temps); + + // Convert atom_densities to a vector + double_1dvec atom_densities_vec; + atom_densities_vec.assign(atom_densities, atom_densities + n_nuclides); + + // Build array of pointers to nuclides_MG's Mgxs objects needed for this + // material + std::vector mgxs_ptr(n_nuclides); + for (int n = 0; n < n_nuclides; n++) { + mgxs_ptr[n] = &nuclides_MG[i_nuclides[n] - 1]; + } + + macro.build_macro(mat_name, temperature, mgxs_ptr, atom_densities_vec, + method, tolerance); + } + macro_xs.push_back(macro); +} + +//============================================================================== +// 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) +{ + 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 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); + + // adjust return value for fortran indexing + gout = gout_c + 1; + + // Rotate the angle + rotate_angle_c(uvw, mu, nullptr); +} + + +void sample_fission_energy_c(const int i_mat, 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); + + // adjust return values for fortran indexing + dg = dg_c + 1; + gout = gout_c + 1; +} + + +void get_name_c(const int index, int name_len, char* name) +{ + // First blank out our input string + std::string str(name_len, ' '); + std::strcpy(name, str.c_str()); + + // Now get the data and copy to the C-string + str = nuclides_MG[index - 1].name; + std::strcpy(name, str.c_str()); + + // Finally, remove the null terminator + name[std::strlen(name)] = ' '; +} + + +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 \ No newline at end of file diff --git a/src/mgxs_interface.h b/src/mgxs_interface.h new file mode 100644 index 0000000000..197ecb0bf9 --- /dev/null +++ b/src/mgxs_interface.h @@ -0,0 +1,66 @@ +//! \file mgxs_interface.h +//! A collection of C interfaces to the C++ Mgxs class + +#ifndef MGXS_INTERFACE_H +#define MGXS_INTERFACE_H + +#include "mgxs.h" + + +namespace openmc { + +extern std::vector nuclides_MG; +extern std::vector 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" 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); + +extern "C" 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); + +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_fission_energy_c(const int i_mat, 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_macro_xs_c(const int index, 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 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 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, + const double sqrtkT); + +extern "C" void reset_nuclide_temperature_index_c(const int index, + const int last_temp); + +} // namespace openmc +#endif // MGXS_INTERFACE_H \ No newline at end of file diff --git a/src/output.F90 b/src/output.F90 index 23d9512eee..9c5d9cb451 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -12,7 +12,7 @@ module output use math, only: t_percentile use mesh_header, only: RegularMesh, meshes use message_passing, only: master, n_procs - use mgxs_header, only: nuclides_MG + use mgxs_interface use nuclide_header use particle_header, only: LocalCoord, Particle use plot_header @@ -680,6 +680,7 @@ contains character(36) :: score_name ! names of scoring function ! to be applied at write-time type(TallyFilterMatch), allocatable :: matches(:) + character(MAX_WORD_LEN) :: temp_name ! Skip if there are no tallies if (n_tallies == 0) return @@ -843,8 +844,9 @@ contains write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), & trim(nuclides(i_nuclide) % name) else + call get_name_c(i_nuclide, len(temp_name), temp_name) write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(nuclides_MG(i_nuclide) % obj % name) + trim(temp_name) end if end if diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index c3d25151b6..200773c256 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -6,7 +6,7 @@ module particle_restart use constants use error, only: write_message use hdf5_interface, only: file_open, file_close, read_dataset, HID_T - use mgxs_header, only: energy_bin_avg + use mgxs_interface, only: energy_bin_avg use nuclide_header, only: micro_xs, n_nuclides use output, only: print_particle use particle_header, only: Particle diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 3030fc099b..797514e25f 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -8,13 +8,12 @@ module physics_mg use material_header, only: Material, materials use math, only: rotate_angle use mesh_header, only: meshes - use mgxs_header + use mgxs_interface use message_passing use nuclide_header, only: material_xs use particle_header, only: Particle use physics_common use random_lcg, only: prn - use scattdata_header use settings use simulation_header use string, only: to_str @@ -22,20 +21,6 @@ module physics_mg implicit none - interface - subroutine scatter_c(i_mat, gin, gout, mu, wgt, uvw) & - bind(C, name='scatter') - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: i_mat - integer(C_INT), value, intent(in) :: gin - integer(C_INT), intent(inout) :: gout - real(C_DOUBLE), intent(inout) :: mu - real(C_DOUBLE), intent(inout) :: wgt - real(C_DOUBLE), intent(inout) :: uvw(1:3) - end subroutine scatter_c - end interface - contains !=============================================================================== @@ -157,14 +142,8 @@ contains type(Particle), intent(inout) :: p - ! call macro_xs(p % material) % obj % sample_scatter(p % coord(1) % uvw, & - ! p % last_g, p % g, p % mu, p % wgt) - - ! Convert change in angle (mu) to new direction - ! p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, p % mu) - - call scatter_c(p % material, p % last_g, p % g, p % mu, p % wgt, & - p % coord(1) % uvw) + 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) p % E = energy_bin_avg(p % g) @@ -194,10 +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 - class(Mgxs), pointer :: xs - - ! Get Pointers - xs => macro_xs(p % material) % obj ! TODO: Heat generation from fission @@ -277,7 +252,7 @@ contains ! Sample secondary energy distribution for fission reaction and set energy ! in fission bank - call xs % sample_fission_energy(p % g, bank_array(i) % uvw, 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 aef620ac5d..b12ec1acb7 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -56,35 +56,38 @@ void ScattData::sample_energy(int gin, int& gout, int& i_gout) } -double ScattData::get_xs(const char* xstype, int gin, int* gout, double* mu) +double ScattData::get_xs(const int xstype, int gin, int* gout, double* mu) { // Set the outgoing group offset index as needed int i_gout = 0; if (gout != nullptr) { // short circuit the function if gout is from a zero portion of the // scattering matrix - if ((*gout < gmin[gin]) || (*gout >= gmax[gin])) { // > gmax? + if ((*gout < gmin[gin]) || (*gout > gmax[gin])) { // > gmax? return 0.; } i_gout = *gout - gmin[gin]; } double val = 0.; - if (std::strcmp(xstype, "scatter")) { + switch(xstype) { + case MG_GET_XS_SCATTER: if (gout != nullptr) { val = scattxs[gin] * energy[gin][i_gout]; } else { val = scattxs[gin]; } - } else if (std::strcmp(xstype, "scatter/mult")) { + break; + case MG_GET_XS_SCATTER_MULT: if (gout != nullptr) { val = scattxs[gin] * energy[gin][i_gout] / mult[gin][i_gout]; } else { - val = scattxs[gin] / std::inner_product(mult[gin].begin(), - mult[gin].end(), - energy[gin].begin(), 0.0); + val = scattxs[gin] / + std::inner_product(mult[gin].begin(), mult[gin].end(), + energy[gin].begin(), 0.0); } - } else if (std::strcmp(xstype, "scatter*f_mu/mult")) { + break; + case MG_GET_XS_SCATTER_FMU_MULT: if ((gout != nullptr) && (mu != nullptr)) { val = scattxs[gin] * energy[gin][i_gout] * calc_f(gin, *gout, *mu); } else { @@ -92,7 +95,8 @@ double ScattData::get_xs(const char* xstype, int gin, int* gout, double* mu) // group or mu is not useful fatal_error("Invalid call to get_xs"); } - } else if (std::strcmp(xstype, "scatter*f_mu")) { + break; + case MG_GET_XS_SCATTER_FMU: if ((gout != nullptr) && (mu != nullptr)) { val = scattxs[gin] * energy[gin][i_gout] * calc_f(gin, *gout, *mu) / mult[gin][i_gout]; @@ -101,6 +105,7 @@ double ScattData::get_xs(const char* xstype, int gin, int* gout, double* mu) // group or mu is not useful fatal_error("Invalid call to get_xs"); } + break; } return val; } @@ -205,13 +210,11 @@ void ScattDataLegendre::update_max_val() double ScattDataLegendre::calc_f(int gin, int gout, double mu) { - // TODO: gout >= or gout >? double f; - if ((gout < gmin[gin]) || (gout >= gmax[gin])) { + if ((gout < gmin[gin]) || (gout > gmax[gin])) { f = 0.; } else { - // TODO: size() -1 or just size? - int i_gout = gout - gmin[gin]; //TODO: + 1? + int i_gout = gout - gmin[gin]; f = evaluate_legendre_c(dist[gin][i_gout].size() - 1, dist[gin][i_gout].data(), mu); } @@ -479,19 +482,18 @@ void ScattDataHistogram::init(int_1dvec& in_gmin, int_1dvec& in_gmax, double ScattDataHistogram::calc_f(int gin, int gout, double mu) { - // TODO: gout >= or gout >? double f; - if ((gout < gmin[gin]) || (gout >= gmax[gin])) { + if ((gout < gmin[gin]) || (gout > gmax[gin])) { f = 0.; } else { // Find mu bin - int i_gout = gout - gmin[gin]; //TODO: + 1? + int i_gout = gout - gmin[gin]; int imu; if (mu == 1.) { // use size -2 to have the index one before the end imu = this->mu.size() - 2; } else { - imu = std::floor((mu + 1.) / dmu + 1.); + imu = std::floor((mu + 1.) / dmu + 1.) - 1; } f = fmu[gin][i_gout][imu]; @@ -776,19 +778,18 @@ void ScattDataTabular::init(int_1dvec& in_gmin, int_1dvec& in_gmax, double ScattDataTabular::calc_f(int gin, int gout, double mu) { - // TODO: gout >= or gout >? double f; - if ((gout < gmin[gin]) || (gout >= gmax[gin])) { + if ((gout < gmin[gin]) || (gout > gmax[gin])) { f = 0.; } else { // Find mu bin - int i_gout = gout - gmin[gin]; //TODO: + 1? + int i_gout = gout - gmin[gin]; int imu; if (mu == 1.) { // use size -2 to have the index one before the end imu = this->mu.size() - 2; } else { - imu = std::floor((mu + 1.) / dmu + 1.); + imu = std::floor((mu + 1.) / dmu + 1.) - 1; } double r = (mu - this->mu[imu]) / (this->mu[imu + 1] - this->mu[imu]); @@ -1006,7 +1007,7 @@ void convert_legendre_to_tabular(ScattDataLegendre& leg, tab.dmu = 2. / (n_mu - 1); tab.mu[0] = -1.; for (int imu = 1; imu < n_mu - 1; imu++) { - tab.mu[imu] = -1. + (imu - 1) * tab.dmu; + tab.mu[imu] = -1. + imu * tab.dmu; } tab.mu[n_mu - 1] = 1.; @@ -1032,6 +1033,7 @@ void convert_legendre_to_tabular(ScattDataLegendre& leg, // Now re-normalize for numerical integration issues and to take care of // the above negative fix-up. Also accrue the CDF double norm = 0.; + tab.dist[gin][i_gout][0] = 0.; for (int imu = 1; imu < n_mu; imu++) { norm += 0.5 * tab.dmu * (tab.fmu[gin][i_gout][imu - 1] + tab.fmu[gin][i_gout][imu]); diff --git a/src/scattdata.h b/src/scattdata.h index 4553156e0d..a9a236b2b7 100644 --- a/src/scattdata.h +++ b/src/scattdata.h @@ -38,7 +38,7 @@ class ScattData { virtual void init(int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& in_mult, double_3dvec& coeffs) = 0; void sample_energy(int gin, int& gout, int& i_gout); - double get_xs(const char* xstype, int gin, int* gout, double* mu); + double get_xs(const int xstype, int gin, int* gout, double* mu); void generic_init(int order, int_1dvec in_gmin, int_1dvec in_gmax, double_2dvec in_energy, double_2dvec in_mult); virtual void combine(std::vector& those_scatts, diff --git a/src/simulation.F90 b/src/simulation.F90 index 4f6163429e..ffcfb8bab8 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -20,7 +20,7 @@ module simulation use geometry_header, only: n_cells use material_header, only: n_materials, materials use message_passing - use mgxs_header, only: energy_bins, energy_bin_avg + use mgxs_interface, only: energy_bins, energy_bin_avg use nuclide_header, only: micro_xs, n_nuclides use output, only: header, print_columns, & print_batch_keff, print_generation, print_runtime, & diff --git a/src/source.F90 b/src/source.F90 index 3683611315..a6b6a92ac2 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -14,7 +14,7 @@ module source use hdf5_interface use math use message_passing, only: rank - use mgxs_header, only: rev_energy_bins, num_energy_groups + use mgxs_interface, only: rev_energy_bins, num_energy_groups use output, only: write_message use particle_header, only: Particle use random_lcg, only: prn, set_particle_seed, prn_set_stream diff --git a/src/state_point.F90 b/src/state_point.F90 index 27fdd969b5..ebb473d496 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -22,7 +22,7 @@ module state_point use hdf5_interface use mesh_header, only: RegularMesh, meshes, n_meshes use message_passing - use mgxs_header, only: nuclides_MG + use mgxs_interface use nuclide_header, only: nuclides use output, only: time_stamp use random_lcg, only: openmc_get_seed, openmc_set_seed @@ -73,6 +73,7 @@ contains character(MAX_WORD_LEN), allocatable :: str_array(:) character(C_CHAR), pointer :: string(:) character(len=:, kind=C_CHAR), allocatable :: filename_ + character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name err = 0 if (present(filename)) then @@ -307,11 +308,13 @@ contains str_array(j) = nuclides(tally % nuclide_bins(j)) % name end if else - i_xs = index(nuclides_MG(tally % nuclide_bins(j)) % obj % name, '.') + call get_name_c(tally % nuclide_bins(j), len(temp_name), & + temp_name) + i_xs = index(temp_name, '.') if (i_xs > 0) then - str_array(j) = nuclides_MG(tally % nuclide_bins(j)) % obj % name(1 : i_xs-1) + str_array(j) = trim(temp_name(1 : i_xs-1)) else - str_array(j) = nuclides_MG(tally % nuclide_bins(j)) % obj % name + str_array(j) = trim(temp_name) end if end if else diff --git a/src/string_functions.cpp b/src/string_functions.cpp new file mode 100644 index 0000000000..f41ec40c1e --- /dev/null +++ b/src/string_functions.cpp @@ -0,0 +1,30 @@ +#include "string_functions.h" + +namespace openmc { + +std::string& strtrim(std::string& s) +{ + const char* t = " \t\n\r\f\v"; + s.erase(s.find_last_not_of(t) + 1); + s.erase(0, s.find_first_not_of(t)); + return s; +} + + +char* strtrim(char* c_str) +{ + std::string std_str; + std_str.assign(c_str); + strtrim(std_str); + int length = std_str.copy(c_str, std_str.size()); + c_str[length] = '\0'; + return c_str; +} + + +void to_lower(std::string& str) +{ + for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]); +} + +} // namespace openmc \ No newline at end of file diff --git a/src/string_functions.h b/src/string_functions.h index d44982bbd3..bf30612fe8 100644 --- a/src/string_functions.h +++ b/src/string_functions.h @@ -4,38 +4,15 @@ #ifndef STRING_FUNCTIONS_H #define STRING_FUNCTIONS_H -// for string functions -#include -#include -#include #include namespace openmc { -std::string& strtrim(std::string& s) -{ - const char* t = " \t\n\r\f\v"; - s.erase(s.find_last_not_of(t) + 1); - s.erase(0, s.find_first_not_of(t)); - return s; -} +std::string& strtrim(std::string& s); +char* strtrim(char* c_str); -char* strtrim(char* c_str) -{ - std::string std_str; - std_str.assign(c_str); - strtrim(std_str); - int length = std_str.copy(c_str, std_str.size()); - c_str[length] = '\0'; - return c_str; -} - - -void to_lower(std::string& str) -{ - for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]); -} +void to_lower(std::string& str); } // namespace openmc -#endif // STRING_FUNCTIONS_H \ No newline at end of file +#endif // STRING_FUNCTIONS_H diff --git a/src/summary.F90 b/src/summary.F90 index bd6ef7158d..cae81a88bc 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -8,7 +8,7 @@ module summary use material_header, only: Material, n_materials use mesh_header, only: RegularMesh use message_passing - use mgxs_header, only: nuclides_MG + use mgxs_interface use nuclide_header use output, only: time_stamp use settings, only: run_CE @@ -94,7 +94,7 @@ contains num_nuclides = 0 num_macros = 0 do i = 1, n_nuclides - if (nuclides_MG(i) % obj % awr /= MACROSCOPIC_AWR) then + if (get_awr_c(i) /= MACROSCOPIC_AWR) then num_nuclides = num_nuclides + 1 else num_macros = num_macros + 1 @@ -119,12 +119,14 @@ contains nuc_names(i) = nuclides(i) % name awrs(i) = nuclides(i) % awr else - if (nuclides_MG(i) % obj % awr /= MACROSCOPIC_AWR) then - nuc_names(j) = nuclides_MG(i) % obj % name - awrs(j) = nuclides_MG(i) % obj % awr + if (get_awr_c(i) /= MACROSCOPIC_AWR) then + call get_name_c(i, len(nuc_names(j)), nuc_names(j)) + nuc_names(j) = trim(nuc_names(j)) + awrs(j) = get_awr_c(i) j = j + 1 else - macro_names(k) = nuclides_MG(i) % obj % name + call get_name_c(i, len(macro_names(k)), macro_names(k)) + macro_names(k) = trim(macro_names(k)) k = k + 1 end if end if @@ -458,7 +460,7 @@ contains num_nuclides = 0 num_macros = 0 do j = 1, m % n_nuclides - if (nuclides_MG(m % nuclide(j)) % obj % awr /= MACROSCOPIC_AWR) then + if (get_awr_c(m % nuclide(j)) /= MACROSCOPIC_AWR) then num_nuclides = num_nuclides + 1 else num_macros = num_macros + 1 @@ -484,12 +486,14 @@ contains k = 1 n = 1 do j = 1, m % n_nuclides - if (nuclides_MG(m % nuclide(j)) % obj % awr /= MACROSCOPIC_AWR) then - nuc_names(k) = nuclides_MG(m % nuclide(j)) % obj % name + if (get_awr_c(m % nuclide(j)) /= MACROSCOPIC_AWR) then + call get_name_c(m % nuclide(j), len(nuc_names(k)), nuc_names(k)) + nuc_names(k) = trim(nuc_names(k)) nuc_densities(k) = m % atom_density(j) k = k + 1 else - macro_names(n) = nuclides_MG(m % nuclide(j)) % obj % name + call get_name_c(m % nuclide(j), len(macro_names(n)), macro_names(n)) + macro_names(n) = trim(macro_names(n)) n = n + 1 end if end do diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index ef8b5915c1..4c24717b43 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -10,7 +10,7 @@ module tally use math, only: t_percentile use mesh_header, only: RegularMesh, meshes use message_passing - use mgxs_header + use mgxs_interface use nuclide_header use output, only: header use particle_header, only: LocalCoord, Particle @@ -1229,8 +1229,14 @@ contains real(8) :: p_uvw(3) ! Particle's current uvw integer :: p_g ! Particle group to use for getting info ! to tally with. - class(Mgxs), pointer :: matxs - class(Mgxs), pointer :: nucxs + ! 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) ! Set the direction and group to use with get_xs if (t % estimator == ESTIMATOR_ANALOG .or. & @@ -1268,13 +1274,15 @@ contains ! To significantly reduce de-referencing, point matxs to the ! macroscopic Mgxs for the material of interest - matxs => macro_xs(p % material) % obj + call set_macro_angle_index_c(p % material, p_uvw, last_mat_pol, & + last_mat_azi, last_mat_uvw) ! Do same for nucxs, point it to the microscopic nuclide data of interest if (i_nuclide > 0) then - nucxs => nuclides_MG(i_nuclide) % obj ! And since we haven't calculated this temperature index yet, do so now - call nucxs % find_temperature(p % sqrtkT) + 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) end if i = 0 @@ -1329,13 +1337,13 @@ contains if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('total', p_g, UVW=p_uvw) / & - matxs % get_xs('total', p_g, UVW=p_uvw) * flux + 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 end if else if (i_nuclide > 0) then - score = nucxs % get_xs('total', p_g, UVW=p_uvw) * & + score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_TOTAL, p_g) * & atom_density * flux else score = material_xs % total * flux @@ -1358,19 +1366,23 @@ contains end if if (i_nuclide > 0) then - score = score * nucxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) & - / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux + score = score * get_nuclide_xs_c(i_nuclide, & + MG_GET_XS_INVERSE_VELOCITY, p_g) / & + get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) * flux else - score = score * matxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) & - / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux + score = score * get_macro_xs_c(p % material, & + MG_GET_XS_INVERSE_VELOCITY, p_g) / & + get_macro_xs_c(p % material, MG_GET_XS_ABSORPTION, p_g) * flux end if else if (i_nuclide > 0) then - score = flux * nucxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) + score = flux * get_nuclide_xs_c(i_nuclide, & + MG_GET_XS_INVERSE_VELOCITY, p_g) else - score = flux * matxs % get_xs('inverse-velocity', p_g, UVW=p_uvw) + score = flux * get_macro_xs_c(p % material, & + MG_GET_XS_INVERSE_VELOCITY, p_g) end if end if @@ -1392,21 +1404,23 @@ contains ! adjust the score by the actual probability for that nuclide. if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, & - UVW=p_uvw, MU=p % mu) / & - matxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, & - UVW=p_uvw, MU=p % mu) + 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, 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 * & - nucxs % get_xs('scatter/mult', p_g, UVW=p_uvw) + 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 * & - matxs % get_xs('scatter/mult', p_g, UVW=p_uvw) + get_macro_xs_c(p % material, MG_GET_XS_SCATTER_MULT, & + p_g, MU=p % mu) end if end if @@ -1428,19 +1442,20 @@ contains ! adjust the score by the actual probability for that nuclide. if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('scatter*f_mu', p % last_g, p % g, & - UVW=p_uvw, MU=p % mu) / & - matxs % get_xs('scatter*f_mu', p % last_g, p % g, & - UVW=p_uvw, MU=p % mu) + 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, MG_GET_XS_SCATTER_FMU, & + p % last_g, p % g, MU=p % mu) end if else if (i_nuclide > 0) then - score = nucxs % get_xs('scatter', p_g, UVW=p_uvw) * & - atom_density * flux + score = atom_density * flux * & + get_nuclide_xs_c(i_nuclide, MG_GET_XS_SCATTER, p_g) else ! Get the scattering x/s, which includes multiplication - score = matxs % get_xs('scatter', p_g, UVW=p_uvw) * flux + score = flux * & + get_macro_xs_c(p % material, MG_GET_XS_SCATTER, p_g) end if end if @@ -1460,13 +1475,13 @@ contains end if if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('absorption', p_g, UVW=p_uvw) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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 = nucxs % get_xs('absorption', p_g, UVW=p_uvw) * & - atom_density * flux + score = atom_density * flux * & + get_nuclide_xs_c(i_nuclide, MG_GET_XS_ABSORPTION, p_g) else score = material_xs % absorption * flux end if @@ -1491,19 +1506,19 @@ contains end if if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('fission', p_g, UVW=p_uvw) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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 * & - matxs % get_xs('fission', p_g, UVW=p_uvw) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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 = nucxs % get_xs('fission', p_g, UVW=p_uvw) * & + score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_FISSION, p_g) * & atom_density * flux else - score = matxs % get_xs('fission', p_g, UVW=p_uvw) * flux + score = get_macro_xs_c(p % material, MG_GET_XS_FISSION, p_g) * flux end if end if @@ -1529,12 +1544,12 @@ contains score = p % absorb_wgt * flux if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('nu-fission', p_g, UVW=p_uvw) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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 * & - matxs % get_xs('nu-fission', p_g, UVW=p_uvw) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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 @@ -1547,17 +1562,17 @@ contains score = keff * p % wgt_bank * flux if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('fission', p_g, UVW=p_uvw) / & - matxs % get_xs('fission', p_g, UVW=p_uvw) + 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 = nucxs % get_xs('nu-fission', p_g, UVW=p_uvw) * & + score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_NU_FISSION, p_g) * & atom_density * flux else - score = matxs % get_xs('nu-fission', p_g, UVW=p_uvw) * flux + score = get_macro_xs_c(p % material, MG_GET_XS_NU_FISSION, p_g) * flux end if end if @@ -1583,12 +1598,12 @@ contains score = p % absorb_wgt * flux if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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 * & - matxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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 @@ -1602,17 +1617,17 @@ contains / real(p % n_bank, 8)) * flux if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('fission', p_g, UVW=p_uvw) / & - matxs % get_xs('fission', p_g, UVW=p_uvw) + 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 = nucxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) * & + score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_PROMPT_NU_FISSION, p_g) * & atom_density * flux else - score = matxs % get_xs('prompt-nu-fission', p_g, UVW=p_uvw) * flux + score = get_macro_xs_c(p % material, MG_GET_XS_PROMPT_NU_FISSION, p_g) * flux end if end if @@ -1638,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 (matxs % get_xs('absorption', p_g, UVW=p_uvw) > 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) @@ -1653,13 +1668,13 @@ contains score = p % absorb_wgt * flux if (i_nuclide > 0) then - score = score * nucxs % get_xs('delayed-nu-fission', & - p_g, UVW=p_uvw, dg=d) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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) else - score = score * matxs % get_xs('delayed-nu-fission', & - p_g, UVW=p_uvw, dg=d) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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) end if call score_fission_delayed_dg(t, d_bin, score, score_index) @@ -1669,11 +1684,13 @@ contains else score = p % absorb_wgt * flux if (i_nuclide > 0) then - score = score * nucxs % get_xs('delayed-nu-fission', p_g, & - UVW=p_uvw) / matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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) else - score = score * matxs % get_xs('delayed-nu-fission', p_g, & - UVW=p_uvw) / matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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) end if end if end if @@ -1703,8 +1720,8 @@ contains if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('fission', p_g, UVW=p_uvw) / & - matxs % get_xs('fission', p_g, UVW=p_uvw) + 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) @@ -1715,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 * & - nucxs % get_xs('fission', p_g, UVW=p_uvw) / & - matxs % get_xs('fission', p_g, UVW=p_uvw) + 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 @@ -1735,11 +1752,11 @@ contains d = filt % groups(d_bin) if (i_nuclide > 0) then - score = nucxs % get_xs('delayed-nu-fission', p_g, & - UVW=p_uvw, dg=d) * atom_density * flux + score = atom_density * flux * & + get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g, DG=d) else - score = matxs % get_xs('delayed-nu-fission', p_g, & - UVW=p_uvw, dg=d) * flux + score = flux * & + 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) @@ -1748,11 +1765,12 @@ contains end select else if (i_nuclide > 0) then - score = nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw) & - * atom_density * flux + score = atom_density * flux * & + get_nuclide_xs_c(i_nuclide, MG_GET_XS_DELAYED_NU_FISSION, p_g) + else - score = matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw) & - * flux + score = flux * & + get_macro_xs_c(p % material, MG_GET_XS_DELAYED_NU_FISSION, p_g) end if end if end if @@ -1767,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 (matxs % get_xs('absorption', p_g, UVW=p_uvw) > 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) @@ -1782,17 +1800,15 @@ contains score = p % absorb_wgt * flux if (i_nuclide > 0) then - score = score * nucxs % get_xs('decay rate', p_g, & - UVW=p_uvw, dg=d) * & - nucxs % get_xs('delayed-nu-fission', p_g, & - UVW=p_uvw, dg=d) / matxs % get_xs('absorption', & - p_g, UVW=p_uvw) + 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) else - score = score * matxs % get_xs('decay rate', p_g, & - UVW=p_uvw, dg=d) * & - matxs % get_xs('delayed-nu-fission', p_g, & - UVW=p_uvw, dg=d) / matxs % get_xs('absorption', & - p_g, UVW=p_uvw) + 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) end if call score_fission_delayed_dg(t, d_bin, score, score_index) @@ -1809,15 +1825,15 @@ contains ! for all delayed groups. do d = 1, num_delayed_groups if (i_nuclide > 0) then - score = score + p % absorb_wgt * & - nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & - nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, & - dg=d) / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux + 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) else - score = score + p % absorb_wgt * & - matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & - matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, & - dg=d) / matxs % get_xs('absorption', p_g, UVW=p_uvw) * flux + 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) end if end do end if @@ -1846,13 +1862,13 @@ contains if (i_nuclide > 0) then score = score + keff * atom_density * & fission_bank(n_bank - p % n_bank + k) % wgt * & - nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=g) * & - nucxs % get_xs('fission', p_g, UVW=p_uvw) / & - matxs % get_xs('fission', p_g, UVW=p_uvw) * 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 * & - matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=g) * 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 @@ -1904,13 +1920,13 @@ contains d = filt % groups(d_bin) if (i_nuclide > 0) then - score = nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & - nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, & - dg=d) * atom_density * flux + 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) else - score = matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & - matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, & - dg=d) * flux + 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) end if call score_fission_delayed_dg(t, d_bin, score, score_index) @@ -1927,12 +1943,12 @@ contains do d = 1, num_delayed_groups if (i_nuclide > 0) then score = score + atom_density * flux * & - nucxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & - nucxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, 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 * & - matxs % get_xs('decay rate', p_g, UVW=p_uvw, dg=d) * & - matxs % get_xs('delayed-nu-fission', p_g, UVW=p_uvw, 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 @@ -1957,19 +1973,20 @@ contains end if if (i_nuclide > 0) then score = score * atom_density * & - nucxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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 * & - matxs % get_xs('kappa-fission', p_g, UVW=p_uvw) / & - matxs % get_xs('absorption', p_g, UVW=p_uvw) + 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 = nucxs % get_xs('kappa-fission', p_g, UVW=p_uvw) * & + score = get_nuclide_xs_c(i_nuclide, MG_GET_XS_KAPPA_FISSION, p_g) * & atom_density * flux else - score = matxs % get_xs('kappa-fission', p_g, UVW=p_uvw) * flux + score = flux * & + get_macro_xs_c(p % material, MG_GET_XS_KAPPA_FISSION, p_g) end if end if @@ -1989,7 +2006,15 @@ contains end do SCORE_LOOP - nullify(matxs, nucxs) + ! 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 !=============================================================================== diff --git a/src/tallies/tally_filter_energy.F90 b/src/tallies/tally_filter_energy.F90 index 93da69edda..f230d438e8 100644 --- a/src/tallies/tally_filter_energy.F90 +++ b/src/tallies/tally_filter_energy.F90 @@ -6,7 +6,7 @@ module tally_filter_energy use constants use error use hdf5_interface - use mgxs_header, only: num_energy_groups, rev_energy_bins + use mgxs_interface, only: num_energy_groups, rev_energy_bins use particle_header, only: Particle use settings, only: run_CE use string, only: to_str diff --git a/src/tracking.F90 b/src/tracking.F90 index 2ba62a8dd8..c832ca9ccc 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -9,7 +9,7 @@ module tracking check_cell_overlap use material_header, only: materials, Material use message_passing - use mgxs_header + use mgxs_interface use nuclide_header use particle_header, only: LocalCoord, Particle use physics, only: collision @@ -29,21 +29,6 @@ module tracking implicit none - interface - subroutine calculate_xs_c(i_mat, gin, sqrtkT, uvw, total_xs, abs_xs, & - nu_fiss_xs) bind(C, name='calculate_xs') - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: i_mat - integer(C_INT), value, intent(in) :: gin - real(C_DOUBLE), value, intent(in) :: sqrtkT - real(C_DOUBLE), intent(in) :: uvw(1:3) - real(C_DOUBLE), intent(inout) :: total_xs - real(C_DOUBLE), intent(inout) :: abs_xs - real(C_DOUBLE), intent(inout) :: nu_fiss_xs - end subroutine calculate_xs_c - end interface - contains !=============================================================================== @@ -129,10 +114,6 @@ contains end if else ! Get the MG data - !!TODO: Remove Fortran call - needed until I'm done replacing Fortran - !!with C++ code because it sets index_temp - call macro_xs(p % material) % obj % calculate_xs(p % g, p % sqrtkT, & - p % coord(p % n_coord) % uvw, material_xs) 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 8afdb7f039..b79d02e200 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -24,6 +24,8 @@ XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable, if (fissionable) { fission = double_3dvec(n_pol, double_2dvec(n_azi, double_1dvec(energy_groups, 0.))); + nu_fission = double_3dvec(n_pol, double_2dvec(n_azi, + double_1dvec(energy_groups, 0.))); prompt_nu_fission = double_3dvec(n_pol, double_2dvec(n_azi, double_1dvec(energy_groups, 0.))); kappa_fission = double_3dvec(n_pol, double_2dvec(n_azi, @@ -514,6 +516,17 @@ void XsData::_fissionable_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, } } + // Combine prompt_nu_fission and delayed_nu_fission into nu_fission + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + nu_fission[p][a][gin] = + std::accumulate(delayed_nu_fission[p][a][gin].begin(), + delayed_nu_fission[p][a][gin].end(), + prompt_nu_fission[p][a][gin]); + } + } + } } @@ -668,6 +681,8 @@ void XsData::combine(std::vector those_xs, double_1dvec& scalars) inverse_velocity[p][a][gin] += scalar * that->inverse_velocity[p][a][gin]; if (that->prompt_nu_fission.size() > 0) { + nu_fission[p][a][gin] += + scalar * that->nu_fission[p][a][gin]; prompt_nu_fission[p][a][gin] += scalar * that->prompt_nu_fission[p][a][gin]; kappa_fission[p][a][gin] += diff --git a/src/xsdata.h b/src/xsdata.h index fec28aece1..478f2e7aa1 100644 --- a/src/xsdata.h +++ b/src/xsdata.h @@ -38,6 +38,7 @@ class XsData { // [phi][theta][incoming group] double_3dvec total; double_3dvec absorption; + double_3dvec nu_fission; double_3dvec prompt_nu_fission; double_3dvec kappa_fission; double_3dvec fission;