From d6f94f7ce8e760033be662d5166148efbbeff10f Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Sun, 2 Sep 2018 18:59:28 -0400 Subject: [PATCH] added specialized fission-getting methods --- include/openmc/mgxs.h | 10 +- include/openmc/xsdata.h | 36 ++ src/input_xml.F90 | 2 +- src/mgxs.cpp | 10 +- src/xsdata.cpp | 611 ++++++++---------- .../regression_tests/mg_benchmark/__init__.py | 0 .../mg_benchmark/inputs_true.dat | 33 + .../mg_benchmark/results_true.dat | 2 + tests/regression_tests/mg_benchmark/test.py | 167 +++++ .../mg_benchmark_delayed/__init__.py | 0 .../mg_benchmark_delayed/test.py | 175 +++++ 11 files changed, 678 insertions(+), 368 deletions(-) create mode 100644 tests/regression_tests/mg_benchmark/__init__.py create mode 100644 tests/regression_tests/mg_benchmark/inputs_true.dat create mode 100644 tests/regression_tests/mg_benchmark/results_true.dat create mode 100644 tests/regression_tests/mg_benchmark/test.py create mode 100644 tests/regression_tests/mg_benchmark_delayed/__init__.py create mode 100644 tests/regression_tests/mg_benchmark_delayed/test.py diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index e80355bee7..cd312ea170 100644 --- a/include/openmc/mgxs.h +++ b/include/openmc/mgxs.h @@ -159,7 +159,7 @@ class Mgxs { //! @param dg delayed group index; use nullptr if irrelevant. //! @return Requested cross section value. double - get_xs(int xstype, int gin, int* gout, double* mu, int* dg); + get_xs(const int xstype, const int gin, int* gout, double* mu, int* dg); //! \brief Samples the fission neutron energy and if prompt or delayed. //! @@ -167,7 +167,7 @@ class Mgxs { //! @param dg Sampled delayed group index. //! @param gout Sampled outgoing energy group. void - sample_fission_energy(int gin, int& dg, int& gout); + sample_fission_energy(const int gin, int& dg, int& gout); //! \brief Samples the outgoing energy and angle from a scatter event. //! @@ -176,7 +176,7 @@ class Mgxs { //! @param mu Sampled cosine of the change-in-angle. //! @param wgt Weight of the particle to be adjusted. void - sample_scatter(int gin, int& gout, double& mu, double& wgt); + sample_scatter(const int gin, int& gout, double& mu, double& wgt); //! \brief Calculates cross section quantities needed for tracking. //! @@ -187,14 +187,14 @@ class Mgxs { //! @param abs_xs Resultant absorption cross section. //! @param nu_fiss_xs Resultant nu-fission cross section. void - calculate_xs(int gin, double sqrtkT, const double uvw[3], + calculate_xs(const int gin, const 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 sqrtkT Temperature of the material. void - set_temperature_index(double sqrtkT); + set_temperature_index(const double sqrtkT); //! \brief Sets the angle index in cache given a direction //! diff --git a/include/openmc/xsdata.h b/include/openmc/xsdata.h index d156de9d95..780855ecfd 100644 --- a/include/openmc/xsdata.h +++ b/include/openmc/xsdata.h @@ -33,6 +33,42 @@ class XsData { fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, size_t delayed_groups, bool is_isotropic); + //! \brief Reads fission data formatted as chi and nu-fission vectors from + // the HDF5 file when beta is provided. + void + fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups); + + //! \brief Reads fission data formatted as chi and nu-fission vectors from + // the HDF5 file when beta is not provided. + void + fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups); + + //! \brief Reads fission data formatted as chi and nu-fission vectors from + // the HDF5 file when no delayed data is provided. + void + fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups); + + //! \brief Reads fission data formatted as a nu-fission matrix from + // the HDF5 file when beta is provided. + void + fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups); + + //! \brief Reads fission data formatted as a nu-fission matrix from + // the HDF5 file when beta is not provided. + void + fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups); + + //! \brief Reads fission data formatted as a nu-fission matrix from + // the HDF5 file when no delayed data is provided. + void + fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups); + public: // The following quantities have the following dimensions: diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 010e09973e..289488c344 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -183,7 +183,7 @@ contains ! Assign temperatures to cells that don't have temperatures already assigned call assign_temperatures() - ! Determine desired txemperatures for each nuclide and S(a,b) table + ! Determine desired temperatures for each nuclide and S(a,b) table call get_temperatures(nuc_temps, sab_temps) ! Check to make sure there are not too many nested coordinate levels in the diff --git a/src/mgxs.cpp b/src/mgxs.cpp index ed68a12e91..93b7d43b77 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -423,7 +423,7 @@ Mgxs::combine(const std::vector& micros, const std::vector& scala //============================================================================== double -Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg) +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 #ifdef _OPENMP @@ -535,7 +535,7 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg) //============================================================================== void -Mgxs::sample_fission_energy(int gin, int& dg, int& gout) +Mgxs::sample_fission_energy(const int gin, int& dg, int& gout) { // This method assumes that the temperature and angle indices are set #ifdef _OPENMP @@ -599,7 +599,7 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout) //============================================================================== void -Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt) +Mgxs::sample_scatter(const int gin, int& gout, double& mu, double& wgt) { // This method assumes that the temperature and angle indices are set // Sample the data @@ -614,7 +614,7 @@ Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt) //============================================================================== void -Mgxs::calculate_xs(int gin, double sqrtkT, const double uvw[3], +Mgxs::calculate_xs(const int gin, const double sqrtkT, const double uvw[3], double& total_xs, double& abs_xs, double& nu_fiss_xs) { // Set our indices @@ -649,7 +649,7 @@ Mgxs::equiv(const Mgxs& that) //============================================================================== void -Mgxs::set_temperature_index(double sqrtkT) +Mgxs::set_temperature_index(const double sqrtkT) { // See if we need to find the new index #ifdef _OPENMP diff --git a/src/xsdata.cpp b/src/xsdata.cpp index d30bc65f83..18e474987c 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -121,6 +121,241 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, //============================================================================== +void +XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups) +{ + // Data is provided as nu-fission and chi with a beta for delayed info + + // Get chi + xt::xtensor temp_chi({n_ang, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "chi", temp_chi, true); + + // Normalize chi by summing over the outgoing groups for each incoming angle + temp_chi = temp_chi / xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis()); + + // Now every incoming group in prompt_chi and delayed_chi is the normalized + // chi we just made + for (size_t a = 0; a < n_ang; a++) { + for (size_t gin = 0; gin < energy_groups; gin++) { + for (size_t gout = 0; gout < energy_groups; gout++) { + chi_prompt(a, gin, gout) = temp_chi(a, gout); + for (size_t dg = 0; dg < delayed_groups; dg++) { + chi_delayed(a, gin, gout, dg) = temp_chi(a, gout); + } + } + } + } + + // Get nu-fission + xt::xtensor temp_nufiss({n_ang, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "nu-fission", temp_nufiss, true); + + // Get beta + xt::xtensor temp_beta({n_ang, delayed_groups}, 0.); + read_nd_vector(xsdata_grp, "beta", temp_beta, true); + + // Set prompt_nu_fission = (1. - beta_total)*nu_fission + prompt_nu_fission = temp_nufiss * (1. - xt::sum(temp_beta, {1})); + + // Set delayed_nu_fission as beta * nu_fission + for (size_t a = 0; a < n_ang; a++) { + for (size_t gin = 0; gin < energy_groups; gin++) { + for (size_t dg = 0; dg < delayed_groups; dg++) { + delayed_nu_fission(a, gin, dg) = temp_beta(a, dg) * temp_nufiss(a, gin); + } + } + } +} + +void +XsData::fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups) +{ + // Data is provided separately as prompt + delayed nu-fission and chi + + // Get chi-prompt + xt::xtensor temp_chi_p({n_ang, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "chi-prompt", temp_chi_p, true); + + // Normalize chi by summing over the outgoing groups for each incoming angle + temp_chi_p = temp_chi_p / xt::view(xt::sum(temp_chi_p, {1}), xt::all(), xt::newaxis()); + + // Get chi-delayed + xt::xtensor temp_chi_d({n_ang, energy_groups, delayed_groups}, 0.); + read_nd_vector(xsdata_grp, "chi-delayed", temp_chi_d, true); + + // Normalize chi by summing over the outgoing groups for each incoming angle + temp_chi_d = temp_chi_d / xt::view(xt::sum(temp_chi_d, {1}), xt::all(), + xt::newaxis(), xt::all()); + + // Now assign the prompt and delayed chis by replicating for each incoming group + for (size_t a = 0; a < n_ang; a++) { + for (size_t gin = 0; gin < energy_groups; gin++) { + for (size_t gout = 0; gout < energy_groups; gout++) { + chi_prompt(a, gin, gout) = temp_chi_p(a, gout); + for (size_t dg = 0; dg < delayed_groups; dg++) { + chi_delayed(a, gin, gout, dg) = temp_chi_d(a, gout, dg); + } + } + } + } + + // Get prompt and delayed nu-fission directly + read_nd_vector(xsdata_grp, "prompt-nu-fission", prompt_nu_fission, true); + read_nd_vector(xsdata_grp, "delayed-nu-fission", delayed_nu_fission, true); +} + +void +XsData::fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups) +{ + // No beta is provided and there is no prompt/delay distinction. + // Therefore, the code only considers the data as prompt. + + // Get chi + xt::xtensor temp_chi({n_ang, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "chi", temp_chi, true); + + // Normalize chi by summing over the outgoing groups for each incoming angle + temp_chi = temp_chi / xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis()); + + // Now every incoming group in self.chi is the normalized chi we just made + for (size_t a = 0; a < n_ang; a++) { + for (size_t gin = 0; gin < energy_groups; gin++) { + for (size_t gout = 0; gout < energy_groups; gout++) { + chi_prompt(a, gin, gout) = temp_chi(a, gout); + } + } + } + + // Get nu-fission directly + if (object_exists(xsdata_grp, "prompt-nu-fission")) { + read_nd_vector(xsdata_grp, "prompt-nu-fission", prompt_nu_fission, true); + } else { + read_nd_vector(xsdata_grp, "nu-fission", prompt_nu_fission, true); + } +} + +//============================================================================== + +void +XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups) +{ + // Data is provided as nu-fission and chi with a beta for delayed info + + // Get nu-fission matrix + xt::xtensor temp_matrix({n_ang, energy_groups, energy_groups}); + read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true); + + // Get beta + xt::xtensor temp_beta({n_ang, energy_groups, delayed_groups}, 0.); + read_nd_vector(xsdata_grp, "beta", temp_beta, true); + + // prompt_nu_fission is the sum of this matrix over outgoing groups and + // multiplied by (1 - beta_tot) + prompt_nu_fission = xt::sum(temp_matrix, {2}) * (1. - xt::sum(temp_beta, {2})); + + // delayed_nu_fission is the sum of this matrix over outgoing groups and + // multiplied by beta + for (size_t a = 0; a < n_ang; a++) { + for (size_t gin = 0; gin < energy_groups; gin++) { + double out_sum = 0.; + for (size_t gout = 0; gout < energy_groups; gout++) { + out_sum += temp_matrix(a, gin, gout); + } + for (size_t dg = 0; dg < delayed_groups; dg++) { + delayed_nu_fission(a, gin, dg) = temp_beta(a, dg) * out_sum; + } + } + } + + // Store chi-prompt + for (size_t a = 0; a < n_ang; a++) { + for (size_t gin = 0; gin < energy_groups; gin++) { + for (size_t gout = 0; gout < energy_groups; gout++) { + double beta = 0.; + for (size_t dg = 0; dg < delayed_groups; dg++) { + beta += temp_beta(a, gin, dg); + } + chi_prompt(a, gin, gout) = (1.0 - beta) * temp_matrix(a, gin, gout); + } + } + } + + // Store chi-delayed + for (size_t a = 0; a < n_ang; a++) { + for (size_t gin = 0; gin < energy_groups; gin++) { + for (size_t gout = 0; gout < energy_groups; gout++) { + for (size_t dg = 0; dg < delayed_groups; dg++) { + chi_delayed(a, gin, gout, dg) = temp_beta(a, gin, dg) * temp_matrix(a, gin, gout); + } + } + } + } + + //Normalize both + chi_prompt = chi_prompt / xt::view(xt::sum(chi_prompt, {2}), xt::all(), xt::newaxis()); + chi_delayed = chi_delayed / xt::view(xt::sum(chi_delayed, {2}), xt::all(), + xt::newaxis(), xt::all()); +} + +void +XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups) +{ + // Data is provided separately as prompt + delayed nu-fission and chi + + // Get the prompt nu-fission matrix + xt::xtensor temp_matrix_p({n_ang, energy_groups, energy_groups}); + read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_matrix_p, true); + + // prompt_nu_fission is the sum over outgoing groups + prompt_nu_fission = xt::sum(temp_matrix_p, {2}); + + // chi_prompt is this matrix but normalized over outgoing groups, which we + // have already stored in prompt_nu_fission + chi_prompt = temp_matrix_p / prompt_nu_fission; + + // Get the delayed nu-fission matrix + xt::xtensor temp_matrix_d({n_ang, energy_groups, energy_groups, + delayed_groups}); + read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_matrix_d, true); + + // delayed_nu_fission is the sum over outgoing groups + delayed_nu_fission = xt::sum(temp_matrix_d, {2}); + + // chi_prompt is this matrix but normalized over outgoing groups, which we + // have already stored in prompt_nu_fission + chi_delayed = temp_matrix_d / delayed_nu_fission; +} + +void +XsData::fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups) +{ + // No beta is provided and there is no prompt/delay distinction. + // Therefore, the code only considers the data as prompt. + + // Get nu-fission matrix + xt::xtensor temp_matrix({n_ang, energy_groups, energy_groups}); + if (object_exists(xsdata_grp, "prompt-nu-fission")) { + read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_matrix, true); + } else { + read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true); + } + + // prompt_nu_fission is the sum over outgoing groups + prompt_nu_fission = xt::sum(temp_matrix, {2}); + + // chi_prompt is this matrix but normalized over outgoing groups, which we + // have already stored in prompt_nu_fission + chi_prompt = temp_matrix / prompt_nu_fission; +} + +//============================================================================== + void XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, size_t delayed_groups, bool is_isotropic) @@ -129,369 +364,31 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, read_nd_vector(xsdata_grp, "fission", fission); read_nd_vector(xsdata_grp, "kappa-fission", kappa_fission); - // Set/get beta - xt::xtensor temp_beta({n_ang, energy_groups, delayed_groups}, 0.); - if (object_exists(xsdata_grp, "beta")) { - hid_t xsdata = open_dataset(xsdata_grp, "beta"); - size_t ndims = dataset_ndims(xsdata); - - // raise ndims to make the isotropic ndims the same as angular - if (is_isotropic) ndims += 2; - - if (ndims == 3) { - // Beta is input as [delayed group] - xt::xtensor temp_arr({n_ang * delayed_groups}, 0.); - read_nd_vector(xsdata_grp, "beta", temp_arr); - - // Broadcast to all incoming groups - size_t temp_idx = 0; - for (size_t a = 0; a < n_ang; a++) { - for (size_t dg = 0; dg < delayed_groups; dg++) { - // Set the first group index and copy the rest - temp_beta(a, 0, dg) = temp_arr[temp_idx++]; - for (size_t gin = 1; gin < energy_groups; gin++) { - temp_beta(a, gin, dg) = temp_beta(a, 0, dg); - } - } - } - } else if (ndims == 4) { - // Beta is input as [in group][delayed group] - read_nd_vector(xsdata_grp, "beta", temp_beta); - } else { - fatal_error("beta must be provided as a 3D or 4D array!"); - } - } - - // If chi is provided, set chi-prompt and chi-delayed + // Get the data; the strategy for doing so depends on if the data is provided + // as a nu-fission matrix or a set of chi and nu-fission vectors if (object_exists(xsdata_grp, "chi")) { - xt::xtensor temp_arr ({n_ang, energy_groups}, 0.); - read_nd_vector(xsdata_grp, "chi", temp_arr); - - for (size_t a = 0; a < n_ang; a++) { - // First set the first group - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_prompt(a, 0, gout) = temp_arr(a, gout); - } - - // Now normalize this data - double chi_sum = 0.; - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_sum += chi_prompt(a, 0, gout); - } - - if (chi_sum <= 0.) { - fatal_error("Encountered chi for a group that is <= 0!"); - } - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_prompt(a, 0, gout) /= chi_sum; - } - - // And extend to the remaining incoming groups - for (size_t gin = 1; gin < energy_groups; gin++) { - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_prompt(a, gin, gout) = chi_prompt(a, 0, gout); - } - } - - // Finally set chi-delayed equal to chi-prompt - // Set chi-delayed to chi-prompt - for(size_t gin = 0; gin < energy_groups; gin++) { - for (size_t gout = 0; gout < energy_groups; gout++) { - for (size_t dg = 0; dg < delayed_groups; dg++) { - chi_delayed(a, gin, gout, dg) = chi_prompt(a, gin, gout); - } - } - } - } - } - - // If nu-fission is provided, set prompt- and delayed-nu-fission; - // if nu-fission is a matrix, set chi-prompt and chi-delayed. - if (object_exists(xsdata_grp, "nu-fission")) { - hid_t xsdata = open_dataset(xsdata_grp, "nu-fission"); - size_t ndims = dataset_ndims(xsdata); - // raise ndims to make the isotropic ndims the same as angular - if (is_isotropic) ndims += 2; - - if (ndims == 3) { - // nu-fission is a 3-d array - read_nd_vector(xsdata_grp, "nu-fission", prompt_nu_fission); - - // set delayed-nu-fission and correct prompt-nu-fission with beta - for (size_t a = 0; a < n_ang; a++) { - for (size_t gin = 0; gin < energy_groups; gin++) { - for (size_t dg = 0; dg < delayed_groups; dg++) { - delayed_nu_fission(a, gin, dg) = - temp_beta(a, gin, dg) * prompt_nu_fission(a, gin); - } - - // Correct the prompt-nu-fission using the delayed neutron fraction - if (delayed_groups > 0) { - double beta_sum = 0.; - for (size_t gin = 0; gin < energy_groups; gin++) { - beta_sum += temp_beta(a, gin); - } - - prompt_nu_fission(a, gin) *= (1. - beta_sum); - } - } - } - - } else if (ndims == 4) { - // nu-fission is a matrix - read_nd_vector(xsdata_grp, "nu-fission", chi_prompt); - - // Normalize the chi info so the CDF is 1. - for (size_t a = 0; a < n_ang; a++) { - for (size_t gin = 0; gin < energy_groups; gin++) { - double chi_sum = 0.; - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_sum += chi_prompt(a, gin, gout); - } - - // Set the vector nu-fission from the matrix nu-fission - prompt_nu_fission(a, gin) = chi_sum; - - if (chi_sum >= 0.) { - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_prompt(a, gin, gout) /= chi_sum; - } - } else { - fatal_error("Encountered chi for a group that is <= 0!"); - } - } - - // set all of chi-delayed to chi-prompt - for (size_t gin = 0; gin < energy_groups; gin++) { - for (size_t gout = 0; gout < energy_groups; gout++) { - for (size_t dg = 0; dg < delayed_groups; dg++) { - chi_delayed(a, gin, gout, dg) = chi_prompt(a, gin, gout); - } - } - } - - // Set the delayed-nu-fission and correct prompt-nu-fission with beta - for (size_t gin = 0; gin < energy_groups; gin++) { - for (size_t dg = 0; dg < delayed_groups; dg++) { - delayed_nu_fission(a, gin, dg) = temp_beta(a, gin, dg) * - prompt_nu_fission(a, gin); - } - - // Correct prompt-nu-fission using the delayed neutron fraction - if (delayed_groups > 0) { - double beta_sum = 0.; - for (size_t dg = 0; dg < delayed_groups; dg++) { - beta_sum += temp_beta(a, gin, dg); - } - prompt_nu_fission(a, gin) *= (1. - beta_sum); - } - } - } + if (delayed_groups == 0) { + fission_vector_no_delayed_from_hdf5(xsdata_grp, n_ang, energy_groups); } else { - fatal_error("nu-fission must be provided as a 3D or 4D array!"); - } - - close_dataset(xsdata); - } - - // If chi-prompt is provided, set chi-prompt - if (object_exists(xsdata_grp, "chi-prompt")) { - xt::xtensor temp_arr({n_ang, energy_groups}, 0.); - read_nd_vector(xsdata_grp, "chi-prompt", temp_arr); - - for (size_t a = 0; a < n_ang; a++) { - for (size_t gin = 0; gin < energy_groups; gin++) { - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_prompt(a, gin, gout) = temp_arr(a, gout); - } - - // Normalize chi so its CDF goes to 1 - double chi_sum = 0.; - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_sum += chi_prompt(a, gin, gout); - } - - if (chi_sum >= 0.) { - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_prompt(a, gin, gout) /= chi_sum; - } - } else { - fatal_error("Encountered chi-prompt for a group that is <= 0.!"); - } + if (object_exists(xsdata_grp, "beta")) { + fission_vector_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, + delayed_groups); + } else { + fission_vector_no_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, + delayed_groups); } } - } - - // If chi-delayed is provided, set chi-delayed - if (object_exists(xsdata_grp, "chi-delayed")) { - hid_t xsdata = open_dataset(xsdata_grp, "chi-delayed"); - size_t ndims = dataset_ndims(xsdata); - // raise ndims to make the isotropic ndims the same as angular - if (is_isotropic) ndims += 2; - close_dataset(xsdata); - - if (ndims == 3) { - // chi-delayed is a [in group] vector - xt::xtensor temp_arr({n_ang, energy_groups}, 0.); - read_nd_vector(xsdata_grp, "chi-delayed", temp_arr); - - for (size_t a = 0; a < n_ang; a++) { - // normalize the chi CDF to 1 - double chi_sum = 0.; - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_sum += temp_arr(a, gout); - } - - if (chi_sum <= 0.) { - fatal_error("Encountered chi-delayed for a group that is <= 0!"); - } - - // set chi-delayed - for (size_t gin = 0; gin < energy_groups; gin++) { - for (size_t gout = 0; gout < energy_groups; gout++) { - for (size_t dg = 0; dg < delayed_groups; dg++) { - chi_delayed(a, gin, gout, dg) = temp_arr(a, gout) / chi_sum; - } - } - } - } - } else if (ndims == 4) { - // chi_delayed is a matrix - read_nd_vector(xsdata_grp, "chi-delayed", chi_delayed); - - // Normalize the chi info so the CDF is 1. - for (size_t a = 0; a < n_ang; a++) { - for (size_t dg = 0; dg < delayed_groups; dg++) { - for (size_t gin = 0; gin < energy_groups; gin++) { - double chi_sum = 0.; - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_sum += chi_delayed(a, gin, gout, dg); - } - - if (chi_sum > 0.) { - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_delayed(a, gin, gout, dg) /= chi_sum; - } - } else { - fatal_error("Encountered chi-delayed for a group that is <= 0!"); - } - } - } - } + } else { + if (delayed_groups == 0) { + fission_matrix_no_delayed_from_hdf5(xsdata_grp, n_ang, energy_groups); } else { - fatal_error("chi-delayed must be provided as a 3D or 4D array!"); - } - } - - // Get prompt-nu-fission, if present - if (object_exists(xsdata_grp, "prompt-nu-fission")) { - hid_t xsdata = open_dataset(xsdata_grp, "prompt-nu-fission"); - size_t ndims = dataset_ndims(xsdata); - // raise ndims to make the isotropic ndims the same as angular - if (is_isotropic) ndims += 2; - close_dataset(xsdata); - - if (ndims == 3) { - // prompt-nu-fission is a [in group] vector - read_nd_vector(xsdata_grp, "prompt-nu-fission", prompt_nu_fission); - } else if (ndims == 4) { - // prompt nu fission is a matrix, - // so set prompt_nu_fiss & chi_prompt - xt::xtensor temp_arr({n_ang, energy_groups, energy_groups}, 0.); - read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_arr); - - // The prompt_nu_fission vector from the matrix form - for (size_t a = 0; a < n_ang; a++) { - for (size_t gin = 0; gin < energy_groups; gin++) { - double prompt_sum = 0.; - for (size_t gout = 0; gout < energy_groups; gout++) { - prompt_sum += temp_arr(a, gin, gout); - } - - prompt_nu_fission(a, gin) = prompt_sum; - } - - // The chi_prompt data is just the normalized fission matrix - for (size_t gin= 0; gin < energy_groups; gin++) { - if (prompt_nu_fission(a, gin) > 0.) { - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_prompt(a, gin, gout) = - temp_arr(a, gin, gout) / prompt_nu_fission(a, gin); - } - } else { - fatal_error("Encountered chi-prompt for a group that is <= 0!"); - } - } + if (object_exists(xsdata_grp, "beta")) { + fission_matrix_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, + delayed_groups); + } else { + fission_matrix_no_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, + delayed_groups); } - - } else { - fatal_error("prompt-nu-fission must be provided as a 3D or 4D array!"); - } - } - - // Get delayed-nu-fission, if present - if (object_exists(xsdata_grp, "delayed-nu-fission")) { - hid_t xsdata = open_dataset(xsdata_grp, "delayed-nu-fission"); - size_t ndims = dataset_ndims(xsdata); - close_dataset(xsdata); - // raise ndims to make the isotropic ndims the same as angular - if (is_isotropic) ndims += 2; - - if (ndims == 3) { - // delayed-nu-fission is an [in group] vector - if (temp_beta(0, 0, 0) == 0.) { - fatal_error("cannot set delayed-nu-fission with a 1D array if " - "beta is not provided"); - } - xt::xtensor temp_arr({n_ang, energy_groups}, 0.); - read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr); - - for (size_t a = 0; a < n_ang; a++) { - for (size_t gin = 0; gin < energy_groups; gin++) { - for (size_t dg = 0; dg < delayed_groups; dg++) { - // Set delayed-nu-fission using beta - delayed_nu_fission(a, gin, dg) = - temp_beta(a, gin, dg) * temp_arr(a, gin); - } - } - } - - } else if (ndims == 4) { - read_nd_vector(xsdata_grp, "delayed-nu-fission", - delayed_nu_fission); - - } else if (ndims == 5) { - // This will contain delayed-nu-fission and chi-delayed data - xt::xtensor temp_arr({n_ang, energy_groups, energy_groups, - delayed_groups}, 0.); - read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr); - - // Set the 3D delayed-nu-fission matrix and 4D chi-delayed matrix - // from the 4D delayed-nu-fission matrix - for (size_t a = 0; a < n_ang; a++) { - for (size_t dg = 0; dg < delayed_groups; dg++) { - for (size_t gin = 0; gin < energy_groups; gin++) { - double gout_sum = 0.; - for (size_t gout = 0; gout < energy_groups; gout++) { - gout_sum += temp_arr(a, gin, gout, dg); - chi_delayed(a, gin, gout, dg) = temp_arr(a, gin, gout, dg); - } - delayed_nu_fission(a, gin, dg) = gout_sum; - // Normalize chi-delayed - if (gout_sum > 0.) { - for (size_t gout = 0; gout < energy_groups; gout++) { - chi_delayed(a, gin, gout, dg) /= gout_sum; - } - } else { - fatal_error("Encountered chi-delayed for a group that is <= 0!"); - } - } - } - } - - } else { - fatal_error("prompt-nu-fission must be provided as a 3D, 4D, or 5D " - "array!"); } } diff --git a/tests/regression_tests/mg_benchmark/__init__.py b/tests/regression_tests/mg_benchmark/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/mg_benchmark/inputs_true.dat b/tests/regression_tests/mg_benchmark/inputs_true.dat new file mode 100644 index 0000000000..8856770ddb --- /dev/null +++ b/tests/regression_tests/mg_benchmark/inputs_true.dat @@ -0,0 +1,33 @@ + + + + + + + + + 2g.h5 + + + + + + + + eigenvalue + 1000 + 10 + 5 + + + -929.45 -1e+50 -1e+50 929.45 1e+50 1e+50 + + + + false + + multi-group + + false + + diff --git a/tests/regression_tests/mg_benchmark/results_true.dat b/tests/regression_tests/mg_benchmark/results_true.dat new file mode 100644 index 0000000000..33e4d68f79 --- /dev/null +++ b/tests/regression_tests/mg_benchmark/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +9.949396E-01 2.047218E-02 diff --git a/tests/regression_tests/mg_benchmark/test.py b/tests/regression_tests/mg_benchmark/test.py new file mode 100644 index 0000000000..cc8f5ba34a --- /dev/null +++ b/tests/regression_tests/mg_benchmark/test.py @@ -0,0 +1,167 @@ +import os + +import numpy as np + +import openmc + +from tests.testing_harness import PyAPITestHarness + + +def create_library(): + # Instantiate the energy group data and file object + groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6]) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + + # Make the base, isotropic data + nu = [2.50, 2.50] + fiss = np.array([0.002817, 0.097]) + capture = [0.008708, 0.02518] + absorption = np.add(capture, fiss) + scatter = np.array( + [[[0.31980, 0.06694], [0.004555, -0.0003972]], + [[0.00000, 0.00000], [0.424100, 0.05439000]]]) + total = [0.33588, 0.54628] + chi = [1., 0.] + + mat_1 = openmc.XSdata('mat_1', groups) + mat_1.order = 1 + mat_1.set_nu_fission(np.multiply(nu, fiss)) + mat_1.set_absorption(absorption) + mat_1.set_scatter_matrix(scatter) + mat_1.set_total(total) + mat_1.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_1) + + # Make a version of mat-1 which has a tabular representation of the + # scattering vice Legendre with 33 points + mat_2 = mat_1.convert_scatter_format('tabular', 33) + mat_2.name = 'mat_2' + mg_cross_sections_file.add_xsdata(mat_2) + + # Make a version of mat-1 which has a histogram representation of the + # scattering vice Legendre with 33 bins + mat_3 = mat_1.convert_scatter_format('histogram', 33) + mat_3.name = 'mat_3' + mg_cross_sections_file.add_xsdata(mat_3) + + # Make a version which uses a fission matrix vice chi & nu-fission + mat_4 = openmc.XSdata('mat_4', groups) + mat_4.order = 1 + mat_4.set_nu_fission(np.outer(np.multiply(nu, fiss), chi)) + mat_4.set_absorption(absorption) + mat_4.set_scatter_matrix(scatter) + mat_4.set_total(total) + mg_cross_sections_file.add_xsdata(mat_4) + + # Make an angle-dependent version of mat_1 with 2 polar and 2 azim. angles + mat_5 = mat_1.convert_representation('angle', 2, 2) + mat_5.name = 'mat_5' + mg_cross_sections_file.add_xsdata(mat_5) + + # Make a copy of mat_1 for testing microscopic cross sections + mat_6 = openmc.XSdata('mat_6', groups) + mat_6.order = 1 + mat_6.set_nu_fission(np.multiply(nu, fiss)) + mat_6.set_absorption(absorption) + mat_6.set_scatter_matrix(scatter) + mat_6.set_total(total) + mat_6.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_6) + + # Write the file + mg_cross_sections_file.export_to_hdf5('2g.h5') + + +def create_model(): + create_library() + + # # Make Materials + materials_file = openmc.Materials() + + mat_names = ['base leg', 'base tab', 'base hist', 'base matrix', 'base ang'] + macros = [] + mats = [] + for i in range(len(mat_names)): + macros.append(openmc.Macroscopic('mat_' + str(i + 1))) + mats.append(openmc.Material(name=mat_names[i])) + mats[-1].set_density('macro', 1.0) + mats[-1].add_macroscopic(macros[-1]) + + # Add in the microscopic data + mats.append(openmc.Material(name='micro')) + mats[-1].set_density("sum") + mats[-1].add_nuclide("mat_1", 0.5) + mats[-1].add_nuclide("mat_6", 0.5) + + materials_file += mats + + materials_file.cross_sections = '2g.h5' + + # # Make Geometry + rad_outer = 929.45 + # Set a cell boundary to exist for every material above (exclude the 0) + rads = np.linspace(0., rad_outer, len(mats) + 1, endpoint=True)[1:] + + # Instantiate Universe + root = openmc.Universe(universe_id=0, name='root universe') + cells = [] + + surfs = [] + surfs.append(openmc.XPlane(x0=0., boundary_type='reflective')) + for r, rad in enumerate(rads): + if r == len(rads) - 1: + surfs.append(openmc.XPlane(x0=rad, boundary_type='vacuum')) + else: + surfs.append(openmc.XPlane(x0=rad)) + + # Instantiate Cells + cells = [] + for c in range(len(surfs) - 1): + cells.append(openmc.Cell()) + cells[-1].region = (+surfs[c] & -surfs[c + 1]) + cells[-1].fill = mats[c] + + # Register Cells with Universe + root.add_cells(cells) + + # Instantiate a Geometry, register the root Universe, and export to XML + geometry_file = openmc.Geometry(root) + + # # Make Settings + # Instantiate a Settings object, set all runtime parameters + settings_file = openmc.Settings() + settings_file.energy_mode = "multi-group" + settings_file.tabular_legendre = {'enable': False} + settings_file.batches = 10 + settings_file.inactive = 5 + settings_file.particles = 1000 + + # Build source distribution + INF = 1000. + bounds = [0., -INF, -INF, rads[0], INF, INF] + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:]) + settings_file.source = openmc.source.Source(space=uniform_dist) + + settings_file.output = {'summary': False} + + model = openmc.model.Model() + model.geometry = geometry_file + model.materials = materials_file + model.settings = settings_file + + return model + + +class MGXSTestHarness(PyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = '2g .h5' + if os.path.exists(f): + os.remove(f) + + +def test_mg_benchmark(): + model = create_model() + harness = PyAPITestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/mg_benchmark_delayed/__init__.py b/tests/regression_tests/mg_benchmark_delayed/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/mg_benchmark_delayed/test.py b/tests/regression_tests/mg_benchmark_delayed/test.py new file mode 100644 index 0000000000..f9f04110ef --- /dev/null +++ b/tests/regression_tests/mg_benchmark_delayed/test.py @@ -0,0 +1,175 @@ +import os + +import numpy as np + +import openmc + +from tests.testing_harness import PyAPITestHarness + + +def create_library(): + # Instantiate the energy group data and file object + groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6]) + n_dg = 2 + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + mg_cross_sections_file.num_delayed_groups = n_dg + + beta = np.array([0.003, 0.003]) + one_m_beta = 1. - np.sum(beta) + nu = [2.50, 2.50] + fiss = np.array([0.002817, 0.097]) + capture = [0.008708, 0.02518] + absorption = np.add(capture, fiss) + scatter = np.array( + [[[0.31980, 0.06694], [0.004555, -0.0003972]], + [[0.00000, 0.00000], [0.424100, 0.05439000]]]) + total = [0.33588, 0.54628] + chi = [1., 0.] + + # Make the base data that uses chi & nu-fission vectors with a beta + mat_1 = openmc.XSdata('mat_1', groups) + mat_1.order = 1 + mat_1.num_delayed_groups = 2 + mat_1.set_beta(beta) + mat_1.set_nu_fission(np.multiply(nu, fiss)) + mat_1.set_absorption(absorption) + mat_1.set_scatter_matrix(scatter) + mat_1.set_total(total) + mat_1.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_1) + + # Make a version that uses prompt and delayed version of nufiss and chi + mat_2 = openmc.XSdata('mat_2', groups) + mat_2.order = 1 + mat_2.num_delayed_groups = 2 + mat_2.set_prompt_nu_fission(one_m_beta * np.multiply(nu, fiss)) + delay_nu_fiss = np.zeros((n_dg, groups.num_groups)) + for dg in range(n_dg): + for g in range(groups.num_groups): + delay_nu_fiss[dg, g] = beta[dg] * nu[g] * fiss[g] + mat_2.set_delayed_nu_fission(delay_nu_fiss) + mat_2.set_absorption(absorption) + mat_2.set_scatter_matrix(scatter) + mat_2.set_total(total) + mat_2.set_chi_prompt(chi) + mat_2.set_chi_delayed(np.stack([chi] * n_dg)) + mg_cross_sections_file.add_xsdata(mat_2) + + # Make a version that uses a nu-fission matrix with a beta + mat_3 = openmc.XSdata('mat_3', groups) + mat_3.order = 1 + mat_3.num_delayed_groups = 2 + mat_3.set_beta(beta) + mat_3.set_nu_fission(np.outer(np.multiply(nu, fiss), chi)) + mat_3.set_absorption(absorption) + mat_3.set_scatter_matrix(scatter) + mat_3.set_total(total) + mg_cross_sections_file.add_xsdata(mat_3) + + # Make a version that uses prompt and delayed version of the nufiss matrix + mat_4 = openmc.XSdata('mat_4', groups) + mat_4.order = 1 + mat_4.num_delayed_groups = 2 + mat_4.set_prompt_nu_fission(one_m_beta * np.outer(np.multiply(nu, fiss), chi)) + delay_nu_fiss = np.zeros((n_dg, groups.num_groups, groups.num_groups)) + for dg in range(n_dg): + for g in range(groups.num_groups): + for go in range(groups.num_groups): + delay_nu_fiss[dg, g, go] = beta[dg] * nu[g] * fiss[g] * chi[go] + mat_4.set_delayed_nu_fission(delay_nu_fiss) + mat_4.set_absorption(absorption) + mat_4.set_scatter_matrix(scatter) + mat_4.set_total(total) + mg_cross_sections_file.add_xsdata(mat_4) + + # Write the file + mg_cross_sections_file.export_to_hdf5('2g.h5') + + +def create_model(): + create_library() + + # # Make Materials + materials_file = openmc.Materials() + + mat_names = ['vec beta', 'vec no beta', 'matrix beta', 'matrix no beta'] + macros = [] + mats = [] + for i in range(len(mat_names)): + macros.append(openmc.Macroscopic('mat_' + str(i + 1))) + mats.append(openmc.Material(name=mat_names[i])) + mats[-1].set_density('macro', 1.0) + mats[-1].add_macroscopic(macros[-1]) + + materials_file += mats + + materials_file.cross_sections = '2g.h5' + + # # Make Geometry + rad_outer = 929.45 + # Set a cell boundary to exist for every material above (exclude the 0) + rads = np.linspace(0., rad_outer, len(mats) + 1, endpoint=True)[1:] + + # Instantiate Universe + root = openmc.Universe(universe_id=0, name='root universe') + cells = [] + + surfs = [] + surfs.append(openmc.XPlane(x0=0., boundary_type='reflective')) + for r, rad in enumerate(rads): + if r == len(rads) - 1: + surfs.append(openmc.XPlane(x0=rad, boundary_type='vacuum')) + else: + surfs.append(openmc.XPlane(x0=rad)) + + # Instantiate Cells + cells = [] + for c in range(len(surfs) - 1): + cells.append(openmc.Cell()) + cells[-1].region = (+surfs[c] & -surfs[c + 1]) + cells[-1].fill = mats[c] + + # Register Cells with Universe + root.add_cells(cells) + + # Instantiate a Geometry, register the root Universe, and export to XML + geometry_file = openmc.Geometry(root) + + # # Make Settings + # Instantiate a Settings object, set all runtime parameters + settings_file = openmc.Settings() + settings_file.energy_mode = "multi-group" + settings_file.tabular_legendre = {'enable': False} + settings_file.batches = 10 + settings_file.inactive = 5 + settings_file.particles = 1000 + + # Build source distribution + INF = 1000. + bounds = [0., -INF, -INF, rads[0], INF, INF] + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:]) + settings_file.source = openmc.source.Source(space=uniform_dist) + + settings_file.output = {'summary': False} + + model = openmc.model.Model() + model.geometry = geometry_file + model.materials = materials_file + model.settings = settings_file + + return model + + +class MGXSTestHarness(PyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = '2g .h5' + if os.path.exists(f): + os.remove(f) + + +def test_mg_benchmark(): + model = create_model() + harness = PyAPITestHarness('statepoint.10.h5', model) + harness.main()