diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 12dc21b712..0b66f355c1 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -16,11 +16,6 @@ typedef std::vector double_1dvec; typedef std::vector > double_2dvec; typedef std::vector > > double_3dvec; typedef std::vector > > > double_4dvec; -typedef std::vector > > > > double_5dvec; -typedef std::vector > > > > > double_6dvec; -typedef std::vector int_1dvec; -typedef std::vector > int_2dvec; -typedef std::vector > > int_3dvec; // ============================================================================ // VERSIONING NUMBERS diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 48ce41d0e9..74f5ae5810 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -46,37 +46,36 @@ hid_t file_open(const std::string& filename, char mode, bool parallel=false); void write_string(hid_t group_id, const char* name, const std::string& buffer, bool indep); +void +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, + bool must_have = false); + void read_nd_vector(hid_t obj_id, const char* name, std::vector& result, bool must_have = false); void -read_nd_vector(hid_t obj_id, const char* name, - std::vector >& result, +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, bool must_have = false); void -read_nd_vector(hid_t obj_id, const char* name, - std::vector >& result, bool must_have = false); - -void -read_nd_vector(hid_t obj_id, const char* name, - std::vector > >& result, +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, bool must_have = false); void -read_nd_vector(hid_t obj_id, const char* name, - std::vector > >& result, +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, bool must_have = false); void -read_nd_vector(hid_t obj_id, const char* name, - std::vector > > >& result, +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, bool must_have = false); void -read_nd_vector(hid_t obj_id, const char* name, - std::vector > > > >& result, +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, + bool must_have = false); + +void +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, bool must_have = false); std::vector attribute_shape(hid_t obj_id, const char* name); diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index 40746c291d..e80355bee7 100644 --- a/include/openmc/mgxs.h +++ b/include/openmc/mgxs.h @@ -7,6 +7,8 @@ #include #include +#include "xtensor/xtensor.hpp" + #include "openmc/constants.h" #include "openmc/hdf5_interface.h" #include "openmc/xsdata.h" @@ -35,7 +37,7 @@ struct CacheData { class Mgxs { private: - double_1dvec kTs; // temperature in eV (k * T) + xt::xtensor 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 @@ -44,8 +46,8 @@ class Mgxs { bool is_isotropic; // used to skip search for angle indices if isotropic int n_pol; int n_azi; - double_1dvec polar; - double_1dvec azimuthal; + std::vector polar; + std::vector azimuthal; //! \brief Initializes the Mgxs object metadata //! @@ -62,10 +64,10 @@ class Mgxs { //! @param in_polar Polar angle grid. //! @param in_azimuthal Azimuthal angle grid. void - init(const std::string& in_name, double in_awr, const double_1dvec& in_kTs, + init(const std::string& in_name, double in_awr, const std::vector& in_kTs, bool in_fissionable, int in_scatter_format, int in_num_groups, int in_num_delayed_groups, bool in_is_isotropic, - const double_1dvec& in_polar, const double_1dvec& in_azimuthal); + const std::vector& in_polar, const std::vector& in_azimuthal); //! \brief Initializes the Mgxs object metadata from the HDF5 file //! @@ -80,8 +82,8 @@ class Mgxs { //! @param method Method of choosing nearest temperatures. void metadata_from_hdf5(hid_t xs_id, int in_num_groups, - int in_num_delayed_groups, const double_1dvec& temperature, - double tolerance, int_1dvec& temps_to_read, int& order_dim, + int in_num_delayed_groups, const std::vector& temperature, + double tolerance, std::vector& temps_to_read, int& order_dim, int& method); //! \brief Performs the actual act of combining the microscopic data for a @@ -93,8 +95,8 @@ class Mgxs { //! corresponds to the temperature of interest. //! @param this_t The temperature index of the macroscopic object. void - combine(const std::vector& micros, const double_1dvec& scalars, - const int_1dvec& micro_ts, int this_t); + combine(const std::vector& micros, const std::vector& scalars, + const std::vector& micro_ts, int this_t); //! \brief Checks to see if this and that are able to be combined //! @@ -128,7 +130,7 @@ class Mgxs { //! provides the number of points to use in the tabular representation. //! @param method Method of choosing nearest temperatures. Mgxs(hid_t xs_id, int energy_groups, - int delayed_groups, const double_1dvec& temperature, double tolerance, + int delayed_groups, const std::vector& temperature, double tolerance, int max_order, bool legendre_to_tabular, int legendre_to_tabular_points, int& method); @@ -141,8 +143,8 @@ class Mgxs { //! @param atom_densities Atom densities of those microscopic quantities. //! @param tolerance Tolerance of temperature selection method. //! @param method Method of choosing nearest temperatures. - Mgxs(const std::string& in_name, const double_1dvec& mat_kTs, - const std::vector& micros, const double_1dvec& atom_densities, + Mgxs(const std::string& in_name, const std::vector& mat_kTs, + const std::vector& micros, const std::vector& atom_densities, double tolerance, int& method); //! \brief Provides a cross section value given certain parameters diff --git a/include/openmc/scattdata.h b/include/openmc/scattdata.h index 980238aeb5..7e73224b90 100644 --- a/include/openmc/scattdata.h +++ b/include/openmc/scattdata.h @@ -6,6 +6,8 @@ #include +#include "xtensor/xtensor.hpp" + #include "openmc/constants.h" namespace openmc { @@ -25,23 +27,25 @@ class ScattData { protected: //! \brief Initializes the attributes of the base class. void - base_init(int order, const int_1dvec& in_gmin, const int_1dvec& in_gmax, - const double_2dvec& in_energy, const double_2dvec& in_mult); + base_init(int order, const xt::xtensor& in_gmin, + const xt::xtensor& in_gmax, const double_2dvec& in_energy, + const double_2dvec& in_mult); //! \brief Combines microscopic ScattDatas into a macroscopic one. void base_combine(int max_order, const std::vector& those_scatts, - const double_1dvec& scalars, int_1dvec& in_gmin, int_1dvec& in_gmax, - double_2dvec& sparse_mult, double_3dvec& sparse_scatter); + const double_1dvec& scalars, xt::xtensor& in_gmin, + xt::xtensor& in_gmax, double_2dvec& sparse_mult, + double_3dvec& sparse_scatter); public: - double_2dvec energy; // Normalized p0 matrix for sampling Eout - double_2dvec mult; // nu-scatter multiplication (nu-scatt/scatt) - double_3dvec dist; // Angular distribution - int_1dvec gmin; // minimum outgoing group - int_1dvec gmax; // maximum outgoing group - double_1dvec scattxs; // Isotropic Sigma_{s,g_{in}} + double_2dvec energy; // Normalized p0 matrix for sampling Eout + double_2dvec mult; // nu-scatter multiplication (nu-scatt/scatt) + double_3dvec dist; // Angular distribution + xt::xtensor gmin; // minimum outgoing group + xt::xtensor gmax; // maximum outgoing group + xt::xtensor scattxs; // Isotropic Sigma_{s,g_{in}} //! \brief Calculates the value of normalized f(mu). //! @@ -72,7 +76,7 @@ class ScattData { //! @param in_mult Input sparse multiplicity matrix //! @param coeffs Input sparse scattering matrix virtual void - init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, + init(const xt::xtensor& in_gmin, const xt::xtensor& in_gmax, const double_2dvec& in_mult, const double_3dvec& coeffs) = 0; //! \brief Combines the microscopic data. @@ -97,7 +101,7 @@ class ScattData { //! @param max_order If Legendre this is the maximum value of "n" in "Pn" //! requested; ignored otherwise. //! @return The dense scattering matrix. - virtual double_3dvec + virtual xt::xtensor get_matrix(int max_order) = 0; //! \brief Samples the outgoing energy from the ScattData info. @@ -142,7 +146,7 @@ class ScattDataLegendre: public ScattData { public: void - init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, + init(const xt::xtensor& in_gmin, const xt::xtensor& in_gmax, const double_2dvec& in_mult, const double_3dvec& coeffs); void @@ -163,7 +167,7 @@ class ScattDataLegendre: public ScattData { int get_order() {return dist[0][0].size() - 1;}; - double_3dvec + xt::xtensor get_matrix(int max_order); }; @@ -176,14 +180,14 @@ class ScattDataHistogram: public ScattData { protected: - double_1dvec mu; // Angle distribution mu bin boundaries - double dmu; // Quick storage of the spacing between the mu bin points - double_3dvec fmu; // The angular distribution histogram + xt::xtensor mu; // Angle distribution mu bin boundaries + double dmu; // Quick storage of the mu spacing + double_3dvec fmu; // The angular distribution histogram public: void - init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, + init(const xt::xtensor& in_gmin, const xt::xtensor& in_gmax, const double_2dvec& in_mult, const double_3dvec& coeffs); void @@ -199,7 +203,7 @@ class ScattDataHistogram: public ScattData { int get_order() {return dist[0][0].size();}; - double_3dvec + xt::xtensor get_matrix(int max_order); }; @@ -212,9 +216,9 @@ class ScattDataTabular: public ScattData { protected: - double_1dvec mu; // Angle distribution mu grid points - double dmu; // Quick storage of the spacing between the mu points - double_3dvec fmu; // The angular distribution function + xt::xtensor mu; // Angle distribution mu grid points + double dmu; // Quick storage of the mu spacing + double_3dvec fmu; // The angular distribution function // Friend convert_legendre_to_tabular so it has access to protected // parameters @@ -225,7 +229,7 @@ class ScattDataTabular: public ScattData { public: void - init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, + init(const xt::xtensor& in_gmin, const xt::xtensor& in_gmax, const double_2dvec& in_mult, const double_3dvec& coeffs); void @@ -241,7 +245,7 @@ class ScattDataTabular: public ScattData { int get_order() {return dist[0][0].size();}; - double_3dvec get_matrix(int max_order); + xt::xtensor get_matrix(int max_order); }; //============================================================================== diff --git a/include/openmc/xsdata.h b/include/openmc/xsdata.h index 156708c786..9a602c8b73 100644 --- a/include/openmc/xsdata.h +++ b/include/openmc/xsdata.h @@ -7,6 +7,8 @@ #include #include +#include "xtensor/xtensor.hpp" + #include "openmc/hdf5_interface.h" #include "openmc/scattdata.h" @@ -35,26 +37,26 @@ class XsData { // The following quantities have the following dimensions: // [angle][incoming group] - double_2dvec total; - double_2dvec absorption; - double_2dvec nu_fission; - double_2dvec prompt_nu_fission; - double_2dvec kappa_fission; - double_2dvec fission; - double_2dvec inverse_velocity; + xt::xtensor total; + xt::xtensor absorption; + xt::xtensor nu_fission; + xt::xtensor prompt_nu_fission; + xt::xtensor kappa_fission; + xt::xtensor fission; + xt::xtensor inverse_velocity; // decay_rate has the following dimensions: // [angle][delayed group] - double_2dvec decay_rate; + xt::xtensor decay_rate; // delayed_nu_fission has the following dimensions: // [angle][incoming group][delayed group] - double_3dvec delayed_nu_fission; + xt::xtensor delayed_nu_fission; // chi_prompt has the following dimensions: // [angle][incoming group][outgoing group] - double_3dvec chi_prompt; + xt::xtensor chi_prompt; // chi_delayed has the following dimensions: // [angle][incoming group][outgoing group][delayed group] - double_4dvec chi_delayed; + xt::xtensor chi_delayed; // scatter has the following dimensions: [angle] std::vector > scatter; diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index 5056e5f78b..db987cc02e 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -5,6 +5,9 @@ #include #include +#include "xtensor/xtensor.hpp" +#include "xtensor/xarray.hpp" + #include "hdf5.h" #include "hdf5_hl.h" #ifdef OPENMC_MPI @@ -474,6 +477,106 @@ read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id, if (name) H5Dclose(dset); } +// //***************************************************************************** + +// void +// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, +// bool must_have) +// { +// if (object_exists(obj_id, name)) { +// read_double(obj_id, name, result.data(), true); +// } else if (must_have) { +// fatal_error(std::string("Must provide " + std::string(name) + "!")); +// } +// } + + +// void +// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, +// bool must_have) +// { +// if (object_exists(obj_id, name)) { +// xt::xarray temp; +// read_double(obj_id, name, temp.data(), true); + +// result = xt::adapt(temp, result.shape()); +// } else if (must_have) { +// fatal_error(std::string("Must provide " + std::string(name) + "!")); +// } +// } + +// void +// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, +// bool must_have) +// { +// if (object_exists(obj_id, name)) { +// xt::xarray temp; +// read_double(obj_id, name, temp.data(), true); + +// result = xt::adapt(temp, result.shape()); +// } else if (must_have) { +// fatal_error(std::string("Must provide " + std::string(name) + "!")); +// } +// } + +// void +// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, +// bool must_have) +// { +// if (object_exists(obj_id, name)) { +// xt::xarray temp; +// read_double(obj_id, name, temp.data(), true); + +// result = xt::adapt(temp, result.shape()); +// } else if (must_have) { +// fatal_error(std::string("Must provide " + std::string(name) + "!")); +// } +// } + +// void +// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, +// bool must_have) +// { +// if (object_exists(obj_id, name)) { +// xt::xarray temp; +// read_double(obj_id, name, temp.data(), true); + +// result = xt::adapt(temp, result.shape()); +// } else if (must_have) { +// fatal_error(std::string("Must provide " + std::string(name) + "!")); +// } +// } + + +// void +// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, +// bool must_have) +// { +// if (object_exists(obj_id, name)) { +// xt::xarray temp; +// read_int(obj_id, name, temp.data(), true); + +// result = xt::adapt(temp, result.shape()); +// } else if (must_have) { +// fatal_error(std::string("Must provide " + std::string(name) + "!")); +// } +// } + +// void +// read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, +// bool must_have) +// { +// if (object_exists(obj_id, name)) { +// xt::xarray temp; +// read_int(obj_id, name, temp.data(), true); + +// result = xt::adapt(temp, result.shape()); +// } else if (must_have) { +// fatal_error(std::string("Must provide " + std::string(name) + "!")); +// } +// } + +//***************************************************************************** void read_double(hid_t obj_id, const char* name, double* buffer, bool indep) @@ -545,19 +648,38 @@ read_nd_vector(hid_t obj_id, const char* name, std::vector& result, void -read_nd_vector(hid_t obj_id, const char* name, - std::vector >& result, bool must_have) +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, + bool must_have) { if (object_exists(obj_id, name)) { - int dim1 = result.size(); - int dim2 = result[0].size(); + int dim1 = result.shape()[0]; + double temp_arr[dim1]; + read_double(obj_id, name, temp_arr, true); + + int temp_idx = 0; + for (int i = 0; i < dim1; i++) { + result(i) = temp_arr[temp_idx++]; + } + } else if (must_have) { + fatal_error(std::string("Must provide " + std::string(name) + "!")); + } +} + + +void +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, + bool must_have) +{ + if (object_exists(obj_id, name)) { + int dim1 = result.shape()[0]; + int dim2 = result.shape()[1]; double temp_arr[dim1 * dim2]; read_double(obj_id, name, temp_arr, true); int temp_idx = 0; for (int i = 0; i < dim1; i++) { for (int j = 0; j < dim2; j++) { - result[i][j] = temp_arr[temp_idx++]; + result(i, j) = temp_arr[temp_idx++]; } } } else if (must_have) { @@ -567,19 +689,19 @@ read_nd_vector(hid_t obj_id, const char* name, void -read_nd_vector(hid_t obj_id, const char* name, - std::vector >& result, bool must_have) +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, + bool must_have) { if (object_exists(obj_id, name)) { - int dim1 = result.size(); - int dim2 = result[0].size(); + int dim1 = result.shape()[0]; + int dim2 = result.shape()[1]; int temp_arr[dim1 * dim2]; read_int(obj_id, name, temp_arr, true); int temp_idx = 0; for (int i = 0; i < dim1; i++) { for (int j = 0; j < dim2; j++) { - result[i][j] = temp_arr[temp_idx++]; + result(i, j) = temp_arr[temp_idx++]; } } } else if (must_have) { @@ -589,14 +711,13 @@ read_nd_vector(hid_t obj_id, const char* name, void -read_nd_vector(hid_t obj_id, const char* name, - std::vector > >& result, +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, bool must_have) { if (object_exists(obj_id, name)) { - int dim1 = result.size(); - int dim2 = result[0].size(); - int dim3 = result[0][0].size(); + int dim1 = result.shape()[0]; + int dim2 = result.shape()[1]; + int dim3 = result.shape()[2]; double temp_arr[dim1 * dim2 * dim3]; read_double(obj_id, name, temp_arr, true); @@ -604,7 +725,7 @@ read_nd_vector(hid_t obj_id, const char* name, for (int i = 0; i < dim1; i++) { for (int j = 0; j < dim2; j++) { for (int k = 0; k < dim3; k++) { - result[i][j][k] = temp_arr[temp_idx++]; + result(i, j, k) = temp_arr[temp_idx++]; } } } @@ -614,14 +735,13 @@ read_nd_vector(hid_t obj_id, const char* name, } void -read_nd_vector(hid_t obj_id, const char* name, - std::vector > >& result, +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, bool must_have) { if (object_exists(obj_id, name)) { - int dim1 = result.size(); - int dim2 = result[0].size(); - int dim3 = result[0][0].size(); + int dim1 = result.shape()[0]; + int dim2 = result.shape()[1]; + int dim3 = result.shape()[2]; int temp_arr[dim1 * dim2 * dim3]; read_int(obj_id, name, temp_arr, true); @@ -629,7 +749,7 @@ read_nd_vector(hid_t obj_id, const char* name, for (int i = 0; i < dim1; i++) { for (int j = 0; j < dim2; j++) { for (int k = 0; k < dim3; k++) { - result[i][j][k] = temp_arr[temp_idx++]; + result(i, j, k) = temp_arr[temp_idx++]; } } } @@ -639,15 +759,14 @@ read_nd_vector(hid_t obj_id, const char* name, } void -read_nd_vector(hid_t obj_id, const char* name, - std::vector > > >& result, +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, bool must_have) { if (object_exists(obj_id, name)) { - int dim1 = result.size(); - int dim2 = result[0].size(); - int dim3 = result[0][0].size(); - int dim4 = result[0][0][0].size(); + int dim1 = result.shape()[0]; + int dim2 = result.shape()[1]; + int dim3 = result.shape()[2]; + int dim4 = result.shape()[3]; double temp_arr[dim1 * dim2 * dim3 * dim4]; read_double(obj_id, name, temp_arr, true); @@ -656,7 +775,7 @@ read_nd_vector(hid_t obj_id, const char* name, for (int j = 0; j < dim2; j++) { for (int k = 0; k < dim3; k++) { for (int l = 0; l < dim4; l++) { - result[i][j][k][l] = temp_arr[temp_idx++]; + result(i, j, k, l) = temp_arr[temp_idx++]; } } } @@ -667,16 +786,15 @@ read_nd_vector(hid_t obj_id, const char* name, } void -read_nd_vector(hid_t obj_id, const char* name, - std::vector > > > >& result, +read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, bool must_have) { if (object_exists(obj_id, name)) { - int dim1 = result.size(); - int dim2 = result[0].size(); - int dim3 = result[0][0].size(); - int dim4 = result[0][0][0].size(); - int dim5 = result[0][0][0][0].size(); + int dim1 = result.shape()[0]; + int dim2 = result.shape()[1]; + int dim3 = result.shape()[2]; + int dim4 = result.shape()[3]; + int dim5 = result.shape()[4]; double temp_arr[dim1 * dim2 * dim3 * dim4 * dim5]; read_double(obj_id, name, temp_arr, true); @@ -686,7 +804,7 @@ read_nd_vector(hid_t obj_id, const char* name, for (int k = 0; k < dim3; k++) { for (int l = 0; l < dim4; l++) { for (int m = 0; m < dim5; m++) { - result[i][j][k][l][m] = temp_arr[temp_idx++]; + result(i, j, k, l, m) = temp_arr[temp_idx++]; } } } diff --git a/src/mgxs.cpp b/src/mgxs.cpp index eb13a55f61..ed68a12e91 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -3,12 +3,17 @@ #include #include #include -#include +#include #ifdef _OPENMP #include #endif +#include "xtensor/xmath.hpp" +#include "xtensor/xsort.hpp" +#include "xtensor/xadapt.hpp" +#include "xtensor/xview.hpp" + #include "openmc/error.h" #include "openmc/math_functions.h" #include "openmc/random_lcg.h" @@ -28,14 +33,15 @@ std::vector macro_xs; void Mgxs::init(const std::string& in_name, double in_awr, - const double_1dvec& in_kTs, bool in_fissionable, int in_scatter_format, + const std::vector& in_kTs, bool in_fissionable, int in_scatter_format, int in_num_groups, int in_num_delayed_groups, bool in_is_isotropic, - const double_1dvec& in_polar, const double_1dvec& in_azimuthal) + const std::vector& in_polar, const std::vector& in_azimuthal) { // Set the metadata name = in_name; awr = in_awr; - kTs = in_kTs; + //TODO: Remove adapt when in_KTs is an xtensor + kTs = xt::adapt(in_kTs); fissionable = in_fissionable; scatter_format = in_scatter_format; num_groups = in_num_groups; @@ -61,8 +67,8 @@ Mgxs::init(const std::string& in_name, double in_awr, void Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups, - int in_num_delayed_groups, const double_1dvec& temperature, - double tolerance, int_1dvec& temps_to_read, int& order_dim, int& method) + int in_num_delayed_groups, const std::vector& temperature, + double tolerance, std::vector& temps_to_read, int& order_dim, int& method) { // get name char char_name[MAX_WORD_LEN]; @@ -87,7 +93,7 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups, dset_names[i] = new char[151]; } get_datasets(kT_group, dset_names); - double_1dvec available_temps(num_temps); + xt::xarray available_temps(num_temps); for (int i = 0; i < num_temps; i++) { read_double(kT_group, dset_names[i], &available_temps[i], true); @@ -110,24 +116,20 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups, switch(method) { case TEMPERATURE_NEAREST: - // Find the minimum difference - for (int i = 0; i < temperature.size(); i++) { - std::valarray temp_diff(available_temps.data(), - available_temps.size()); - temp_diff = std::abs(temp_diff - temperature[i]); - int i_closest = std::min_element(std::begin(temp_diff), std::end(temp_diff)) - - std::begin(temp_diff); + // Determine actual temperatures to read + for (const auto& T : temperature) { + auto i_closest = xt::argmin(xt::abs(available_temps - T))[0]; double temp_actual = available_temps[i_closest]; - - if (std::abs(temp_actual - temperature[i]) < tolerance) { - if (std::find(temps_to_read.begin(), temps_to_read.end(), - std::round(temp_actual)) == temps_to_read.end()) { + if (std::fabs(temp_actual - T) < tolerance) { + if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temp_actual)) + == temps_to_read.end()) { temps_to_read.push_back(std::round(temp_actual)); - } else { - fatal_error("MGXS Library does not contain cross section for " + - in_name + " at or near " + - std::to_string(std::round(temperature[i])) + " K."); } + } else { + std::stringstream msg; + msg << "MGXS library does not contain cross sections for " + << in_name << " at or near " << std::round(T) << " K."; + fatal_error(msg); } } break; @@ -160,7 +162,7 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups, // Get the library's temperatures int n_temperature = temps_to_read.size(); - double_1dvec in_kTs(n_temperature); + std::vector in_kTs(n_temperature); for (int i = 0; i < n_temperature; i++) { std::string temp_str(std::to_string(temps_to_read[i]) + "K"); @@ -254,12 +256,12 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups, } // Set the angular bins to use equally-spaced bins - double_1dvec in_polar(in_n_pol); + std::vector in_polar(in_n_pol); double dangle = PI / in_n_pol; for (int p = 0; p < in_n_pol; p++) { in_polar[p] = (p + 0.5) * dangle; } - double_1dvec in_azimuthal(in_n_azi); + std::vector in_azimuthal(in_n_azi); dangle = 2. * PI / in_n_azi; for (int a = 0; a < in_n_azi; a++) { in_azimuthal[a] = (a + 0.5) * dangle - PI; @@ -274,12 +276,12 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups, //============================================================================== Mgxs::Mgxs(hid_t xs_id, int energy_groups, int delayed_groups, - const double_1dvec& temperature, double tolerance, int max_order, + const std::vector& temperature, double tolerance, int max_order, bool legendre_to_tabular, int legendre_to_tabular_points, int& method) { // Call generic data gathering routine (will populate the metadata) int order_data; - int_1dvec temps_to_read; + std::vector temps_to_read; metadata_from_hdf5(xs_id, energy_groups, delayed_groups, temperature, tolerance, temps_to_read, order_data, method); @@ -310,8 +312,8 @@ Mgxs::Mgxs(hid_t xs_id, int energy_groups, int delayed_groups, //============================================================================== -Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs, - const std::vector& micros, const double_1dvec& atom_densities, +Mgxs::Mgxs(const std::string& in_name, const std::vector& mat_kTs, + const std::vector& micros, const std::vector& atom_densities, double tolerance, int& method) { // Get the minimum data needed to initialize: @@ -328,8 +330,8 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs, int in_num_groups = micros[0]->num_groups; int in_num_delayed_groups = micros[0]->num_delayed_groups; bool in_is_isotropic = micros[0]->is_isotropic; - double_1dvec in_polar = micros[0]->polar; - double_1dvec in_azimuthal = micros[0]->azimuthal; + std::vector in_polar = micros[0]->polar; + std::vector in_azimuthal = micros[0]->azimuthal; init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format, in_num_groups, in_num_delayed_groups, in_is_isotropic, in_polar, @@ -345,33 +347,27 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs, // Create the list of temperature indices and interpolation factors for // each microscopic data at the material temperature - int_1dvec micro_t(micros.size(), 0); - double_1dvec micro_t_interp(micros.size(), 0.); + std::vector micro_t(micros.size(), 0); + std::vector micro_t_interp(micros.size(), 0.); for (int m = 0; m < micros.size(); m++) { switch(method) { case TEMPERATURE_NEAREST: { - // Find the nearest temperature - std::valarray temp_diff(micros[m]->kTs.data(), - micros[m]->kTs.size()); - temp_diff = std::abs(temp_diff - temp_desired); - micro_t[m] = std::min_element(std::begin(temp_diff), - std::end(temp_diff)) - - std::begin(temp_diff); - double temp_actual = micros[m]->kTs[micro_t[m]]; + micro_t[m] = xt::argmin(xt::abs(micros[m]->kTs - temp_desired))[0]; + auto temp_actual = micros[m]->kTs[micro_t[m]]; if (std::abs(temp_actual - temp_desired) >= K_BOLTZMANN * tolerance) { - fatal_error("MGXS Library does not contain cross section for " + - name + " at or near " + - std::to_string(std::round(temp_desired / K_BOLTZMANN)) - + " K."); + std::stringstream msg; + msg << "MGXS Library does not contain cross section for " << name + << " at or near " << std::round(temp_desired / K_BOLTZMANN) << "K."; + fatal_error(msg); } } break; case TEMPERATURE_INTERPOLATION: // Get a list of bounding temperatures for each actual temperature // present in the model - for (int k = 0; k < micros[m]->kTs.size() - 1; k++) { + for (int k = 0; k < micros[m]->kTs.shape()[0] - 1; k++) { if ((micros[m]->kTs[k] <= temp_desired) && (temp_desired < micros[m]->kTs[k + 1])) { micro_t[m] = k; @@ -394,8 +390,8 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs, int num_interp_points = 2; if (method == TEMPERATURE_NEAREST) num_interp_points = 1; for (int interp_point = 0; interp_point < num_interp_points; interp_point++) { - double_1dvec interp(micros.size()); - double_1dvec temp_indices(micros.size()); + std::vector interp(micros.size()); + std::vector temp_indices(micros.size()); for (int m = 0; m < micros.size(); m++) { interp[m] = (1. - micro_t_interp[m]) * atom_densities[m]; temp_indices[m] = micro_t[m] + interp_point; @@ -409,8 +405,8 @@ Mgxs::Mgxs(const std::string& in_name, const double_1dvec& mat_kTs, //============================================================================== void -Mgxs::combine(const std::vector& micros, const double_1dvec& scalars, - const int_1dvec& micro_ts, int this_t) +Mgxs::combine(const std::vector& micros, const std::vector& scalars, + const std::vector& micro_ts, int this_t) { // Build the vector of pointers to the xs objects within micros std::vector those_xs(micros.size()); @@ -441,19 +437,19 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg) double val; switch(xstype) { case MG_GET_XS_TOTAL: - val = xs_t->total[a][gin]; + val = xs_t->total(a, gin); break; case MG_GET_XS_NU_FISSION: - val = fissionable ? xs_t->nu_fission[a][gin] : 0.; + val = fissionable ? xs_t->nu_fission(a, gin) : 0.; break; case MG_GET_XS_ABSORPTION: - val = xs_t->absorption[a][gin]; + val = xs_t->absorption(a, gin);; break; case MG_GET_XS_FISSION: - val = fissionable ? xs_t->fission[a][gin] : 0.; + val = fissionable ? xs_t->fission(a, gin) : 0.; break; case MG_GET_XS_KAPPA_FISSION: - val = fissionable ? xs_t->kappa_fission[a][gin] : 0.; + val = fissionable ? xs_t->kappa_fission(a, gin) : 0.; break; case MG_GET_XS_SCATTER: case MG_GET_XS_SCATTER_MULT: @@ -462,16 +458,16 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg) val = xs_t->scatter[a]->get_xs(xstype, gin, gout, mu); break; case MG_GET_XS_PROMPT_NU_FISSION: - val = fissionable ? xs_t->prompt_nu_fission[a][gin] : 0.; + val = fissionable ? xs_t->prompt_nu_fission(a, gin) : 0.; break; case MG_GET_XS_DELAYED_NU_FISSION: if (fissionable) { if (dg != nullptr) { - val = xs_t->delayed_nu_fission[a][gin][*dg]; + val = xs_t->delayed_nu_fission(a, gin, *dg); } else { val = 0.; - for (auto& num : xs_t->delayed_nu_fission[a][gin]) { - val += num; + for (int d = 0; d < xs_t->delayed_nu_fission.shape()[2]; d++) { + val += xs_t->delayed_nu_fission(a, gin, d); } } } else { @@ -481,12 +477,12 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg) case MG_GET_XS_CHI_PROMPT: if (fissionable) { if (gout != nullptr) { - val = xs_t->chi_prompt[a][gin][*gout]; + val = xs_t->chi_prompt(a, gin, *gout); } else { // provide an outgoing group-wise sum val = 0.; - for (auto& num : xs_t->chi_prompt[a][gin]) { - val += num; + for (int g = 0; g < xs_t->chi_prompt.shape()[2]; g++) { + val += xs_t->chi_prompt(a, gin, g); } } } else { @@ -497,21 +493,21 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg) if (fissionable) { if (gout != nullptr) { if (dg != nullptr) { - val = xs_t->chi_delayed[a][gin][*gout][*dg]; + val = xs_t->chi_delayed(a, gin, *gout, *dg); } else { - val = xs_t->chi_delayed[a][gin][*gout][0]; + val = xs_t->chi_delayed(a, gin, *gout, 0); } } else { if (dg != nullptr) { val = 0.; - for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) { - val += xs_t->chi_delayed[a][gin][i][*dg]; + for (int g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) { + val += xs_t->delayed_nu_fission(a, gin, g, *dg); } } else { val = 0.; - for (int i = 0; i < xs_t->chi_delayed[a][gin].size(); i++) { - for (auto& num : xs_t->chi_delayed[a][gin][i]) { - val += num; + for (int g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) { + for (int d = 0; d < xs_t->delayed_nu_fission.shape()[3]; d++) { + val += xs_t->delayed_nu_fission(a, gin, g, d); } } } @@ -521,13 +517,13 @@ Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg) } break; case MG_GET_XS_INVERSE_VELOCITY: - val = xs_t->inverse_velocity[a][gin]; + val = xs_t->inverse_velocity(a, gin); break; case MG_GET_XS_DECAY_RATE: if (dg != nullptr) { - val = xs_t->decay_rate[a][*dg + 1]; + val = xs_t->decay_rate(a, *dg + 1); } else { - val = xs_t->decay_rate[a][0]; + val = xs_t->decay_rate(a, 0); } break; default: @@ -548,11 +544,11 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout) int tid = 0; #endif XsData* xs_t = &xs[cache[tid].t]; - double nu_fission = xs_t->nu_fission[cache[tid].a][gin]; + double nu_fission = xs_t->nu_fission(cache[tid].a, gin); // Find the probability of having a prompt neutron double prob_prompt = - xs_t->prompt_nu_fission[cache[tid].a][gin]; + xs_t->prompt_nu_fission(cache[tid].a, gin); // sample random numbers double xi_pd = prn() * nu_fission; @@ -568,10 +564,10 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout) // sample the outgoing energy group gout = 0; double prob_gout = - xs_t->chi_prompt[cache[tid].a][gin][gout]; + xs_t->chi_prompt(cache[tid].a, gin, gout); while (prob_gout < xi_gout) { gout++; - prob_gout += xs_t->chi_prompt[cache[tid].a][gin][gout]; + prob_gout += xs_t->chi_prompt(cache[tid].a, gin, gout); } } else { @@ -582,7 +578,7 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout) while (xi_pd >= prob_prompt) { dg++; prob_prompt += - xs_t->delayed_nu_fission[cache[tid].a][gin][dg]; + xs_t->delayed_nu_fission(cache[tid].a, gin, dg); } // adjust dg in case of round-off error @@ -591,11 +587,11 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout) // sample the outgoing energy group gout = 0; double prob_gout = - xs_t->chi_delayed[cache[tid].a][gin][gout][dg]; + xs_t->chi_delayed(cache[tid].a, gin, gout, dg); while (prob_gout < xi_gout) { gout++; prob_gout += - xs_t->chi_delayed[cache[tid].a][gin][gout][dg]; + xs_t->chi_delayed(cache[tid].a, gin, gout, dg); } } } @@ -630,10 +626,10 @@ Mgxs::calculate_xs(int gin, double sqrtkT, const double uvw[3], set_temperature_index(sqrtkT); set_angle_index(uvw); XsData* xs_t = &xs[cache[tid].t]; - total_xs = xs_t->total[cache[tid].a][gin]; - abs_xs = xs_t->absorption[cache[tid].a][gin]; + total_xs = xs_t->total(cache[tid].a, gin); + abs_xs = xs_t->absorption(cache[tid].a, gin); - nu_fiss_xs = fissionable ? xs_t->nu_fission[cache[tid].a][gin] : 0.; + nu_fiss_xs = fissionable ? xs_t->nu_fission(cache[tid].a, gin) : 0.; } //============================================================================== @@ -662,17 +658,7 @@ Mgxs::set_temperature_index(double sqrtkT) int tid = 0; #endif if (sqrtkT != cache[tid].sqrtkT) { - double kT = sqrtkT * sqrtkT; - - // initialize vector for storage of the differences - std::valarray temp_diff(kTs.data(), kTs.size()); - - // Find the minimum difference of kT and kTs - temp_diff = std::abs(temp_diff - kT); - cache[tid].t = std::min_element(std::begin(temp_diff), std::end(temp_diff)) - - std::begin(temp_diff); - - // store this temperature as the last one used + cache[tid].t = xt::argmin(xt::abs(kTs - sqrtkT * sqrtkT))[0]; cache[tid].sqrtkT = sqrtkT; } } diff --git a/src/scattdata.cpp b/src/scattdata.cpp index 3e18c168f2..69eb20d8e7 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -4,6 +4,8 @@ #include #include +#include "xtensor/xbuilder.hpp" + #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/math_functions.h" @@ -16,8 +18,8 @@ namespace openmc { //============================================================================== void -ScattData::base_init(int order, const int_1dvec& in_gmin, - const int_1dvec& in_gmax, const double_2dvec& in_energy, +ScattData::base_init(int order, const xt::xtensor& in_gmin, + const xt::xtensor& in_gmax, const double_2dvec& in_energy, const double_2dvec& in_mult) { int groups = in_energy.size(); @@ -53,16 +55,15 @@ ScattData::base_init(int order, const int_1dvec& in_gmin, void ScattData::base_combine(int max_order, const std::vector& those_scatts, const double_1dvec& scalars, - int_1dvec& in_gmin, int_1dvec& in_gmax, double_2dvec& sparse_mult, + xt::xtensor& in_gmin, xt::xtensor& in_gmax, double_2dvec& sparse_mult, double_3dvec& sparse_scatter) { int groups = those_scatts[0] -> energy.size(); // Now allocate and zero our storage spaces - double_3dvec this_matrix = double_3dvec(groups, double_2dvec(groups, - double_1dvec(max_order, 0.))); - double_2dvec mult_numer(groups, double_1dvec(groups, 0.)); - double_2dvec mult_denom(groups, double_1dvec(groups, 0.)); + xt::xtensor this_matrix({groups, groups, max_order}, 0.); + xt::xtensor mult_numer({groups, groups}, 0.); + xt::xtensor mult_denom({groups, groups}, 0.); // Build the dense scattering and multiplicity matrices // Get the multiplicity_matrix @@ -80,26 +81,26 @@ ScattData::base_combine(int max_order, ScattData* that = those_scatts[i]; // Build the dense matrix for that object - double_3dvec that_matrix = that->get_matrix(max_order); + xt::xtensor that_matrix = that->get_matrix(max_order); // Now add that to this for the scattering and multiplicity for (int gin = 0; gin < groups; gin++) { // Only spend time adding that's gmin to gmax data since the rest will // be zeros int i_gout = 0; - for (int gout = that->gmin[gin]; gout <= that->gmax[gin]; gout++) { + for (int gout = that->gmin(gin); gout <= that->gmax(gin); gout++) { // Do the scattering matrix for (int l = 0; l < max_order; l++) { - this_matrix[gin][gout][l] += scalars[i] * that_matrix[gin][gout][l]; + this_matrix(gin, gout, l) += scalars[i] * that_matrix(gin, gout, l); } // Incorporate that's contribution to the multiplicity matrix data - double nuscatt = that->scattxs[gin] * that->energy[gin][i_gout]; - mult_numer[gin][gout] += scalars[i] * nuscatt; + double nuscatt = that->scattxs(gin) * that->energy[gin][i_gout]; + mult_numer(gin, gout) += scalars[i] * nuscatt; if (that->mult[gin][i_gout] > 0.) { - mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][i_gout]; + mult_denom(gin, gout) += scalars[i] * nuscatt / that->mult[gin][i_gout]; } else { - mult_denom[gin][gout] += scalars[i]; + mult_denom(gin, gout) += scalars[i]; } i_gout++; } @@ -107,16 +108,14 @@ ScattData::base_combine(int max_order, } // Combine mult_numer and mult_denom into the combined multiplicity matrix - double_2dvec this_mult(groups, double_1dvec(groups, 1.)); + xt::xtensor this_mult = xt::ones({groups, groups}); for (int gin = 0; gin < groups; gin++) { for (int gout = 0; gout < groups; gout++) { - if (mult_denom[gin][gout] > 0.) { - this_mult[gin][gout] = mult_numer[gin][gout] / mult_denom[gin][gout]; + if (mult_denom(gin, gout) > 0.) { + this_mult(gin, gout) = mult_numer(gin, gout) / mult_denom(gin, gout); } } } - mult_numer.clear(); - mult_denom.clear(); // We have the data, now we need to convert to a jagged array and then use // the initialize function to store it on the object. @@ -125,8 +124,8 @@ ScattData::base_combine(int max_order, int gmin_; for (gmin_ = 0; gmin_ < groups; gmin_++) { bool non_zero = false; - for (int l = 0; l < this_matrix[gin][gmin_].size(); l++) { - if (this_matrix[gin][gmin_][l] != 0.) { + for (int l = 0; l < this_matrix.shape()[2]; l++) { + if (this_matrix(gin, gmin_, l) != 0.) { non_zero = true; break; } @@ -136,8 +135,8 @@ ScattData::base_combine(int max_order, int gmax_; for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) { bool non_zero = false; - for (int l = 0; l < this_matrix[gin][gmax_].size(); l++) { - if (this_matrix[gin][gmax_][l] != 0.) { + for (int l = 0; l < this_matrix.shape()[2]; l++) { + if (this_matrix(gin, gmax_, l) != 0.) { non_zero = true; break; } @@ -160,8 +159,11 @@ ScattData::base_combine(int max_order, sparse_mult[gin].resize(gmax_ - gmin_ + 1); int i_gout = 0; for (int gout = gmin_; gout <= gmax_; gout++) { - sparse_scatter[gin][i_gout] = this_matrix[gin][gout]; - sparse_mult[gin][i_gout] = this_mult[gin][gout]; + sparse_scatter[gin][i_gout].resize(this_matrix.shape()[2]); + for (int l = 0; l < this_matrix.shape()[2]; l++) { + sparse_scatter[gin][i_gout][l] = this_matrix(gin, gout, l); + } + sparse_mult[gin][i_gout] = this_mult(gin, gout); i_gout++; } } @@ -241,8 +243,9 @@ ScattData::get_xs(int xstype, int gin, const int* gout, const double* mu) //============================================================================== void -ScattDataLegendre::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, - const double_2dvec& in_mult, const double_3dvec& coeffs) +ScattDataLegendre::init(const xt::xtensor& in_gmin, + const xt::xtensor& in_gmax, const double_2dvec& in_mult, + const double_3dvec& coeffs) { int groups = coeffs.size(); int order = coeffs[0][0].size(); @@ -252,10 +255,9 @@ ScattDataLegendre::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, // Get the scattering cross section value by summing the un-normalized P0 // coefficient in the variable matrix over all outgoing groups. - scattxs.resize(groups); + scattxs = xt::zeros({groups}); for (int gin = 0; gin < groups; gin++) { int num_groups = in_gmax[gin] - in_gmin[gin] + 1; - scattxs[gin] = 0.; for (int i_gout = 0; i_gout < num_groups; i_gout++) { scattxs[gin] += matrix[gin][i_gout][0]; } @@ -401,8 +403,8 @@ ScattDataLegendre::combine(const std::vector& those_scatts, int groups = those_scatts[0] -> energy.size(); - int_1dvec in_gmin(groups); - int_1dvec in_gmax(groups); + xt::xtensor in_gmin({groups}); + xt::xtensor in_gmax({groups}); double_3dvec sparse_scatter(groups); double_2dvec sparse_mult(groups); @@ -418,20 +420,19 @@ ScattDataLegendre::combine(const std::vector& those_scatts, //============================================================================== -double_3dvec +xt::xtensor ScattDataLegendre::get_matrix(int max_order) { // Get the sizes and initialize the data to 0 int groups = energy.size(); int order_dim = max_order + 1; - double_3dvec matrix = double_3dvec(groups, double_2dvec(groups, - double_1dvec(order_dim, 0.))); + xt::xtensor matrix({groups, groups, order_dim}, 0.); for (int gin = 0; gin < groups; gin++) { for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) { int gout = i_gout + gmin[gin]; for (int l = 0; l < order_dim; l++) { - matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] * + matrix(gin, gout, l) = scattxs[gin] * energy[gin][i_gout] * dist[gin][i_gout][l]; } } @@ -444,8 +445,9 @@ ScattDataLegendre::get_matrix(int max_order) //============================================================================== void -ScattDataHistogram::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, - const double_2dvec& in_mult, const double_3dvec& coeffs) +ScattDataHistogram::init(const xt::xtensor& in_gmin, + const xt::xtensor& in_gmax, const double_2dvec& in_mult, + const double_3dvec& coeffs) { int groups = coeffs.size(); int order = coeffs[0][0].size(); @@ -455,9 +457,8 @@ ScattDataHistogram::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, // Get the scattering cross section value by summing the distribution // over all the histogram bins in angle and outgoing energy groups - scattxs.resize(groups); + scattxs = xt::zeros({groups}); for (int gin = 0; gin < groups; gin++) { - scattxs[gin] = 0.; for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) { scattxs[gin] += std::accumulate(matrix[gin][i_gout].begin(), matrix[gin][i_gout].end(), 0.); @@ -484,12 +485,8 @@ ScattDataHistogram::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, ScattData::base_init(order, in_gmin, in_gmax, in_energy, in_mult); // Build the angular distribution mu values - mu = double_1dvec(order); + mu = xt::linspace(-1., 1., order + 1); dmu = 2. / order; - mu[0] = -1.; - for (int imu = 1; imu < order; imu++) { - mu[imu] = -1. + imu * dmu; - } // Calculate f(mu) and integrate it so we can avoid rejection sampling fmu.resize(groups); @@ -581,21 +578,20 @@ ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt) //============================================================================== -double_3dvec +xt::xtensor ScattDataHistogram::get_matrix(int max_order) { // Get the sizes and initialize the data to 0 int groups = energy.size(); // We ignore the requested order for Histogram and Tabular representations int order_dim = get_order(); - double_3dvec matrix = double_3dvec(groups, double_2dvec(groups, - double_1dvec(order_dim, 0.))); + xt::xtensor matrix = xt::zeros({groups, groups, order_dim}); for (int gin = 0; gin < groups; gin++) { for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) { int gout = i_gout + gmin[gin]; for (int l = 0; l < order_dim; l++) { - matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] * + matrix(gin, gout, l) = scattxs[gin] * energy[gin][i_gout] * fmu[gin][i_gout][l]; } } @@ -624,8 +620,8 @@ ScattDataHistogram::combine(const std::vector& those_scatts, int groups = those_scatts[0] -> energy.size(); - int_1dvec in_gmin(groups); - int_1dvec in_gmax(groups); + xt::xtensor in_gmin({groups}); + xt::xtensor in_gmax({groups}); double_3dvec sparse_scatter(groups); double_2dvec sparse_mult(groups); @@ -633,7 +629,7 @@ ScattDataHistogram::combine(const std::vector& those_scatts, // so we use a base class method to sum up xs and create new energy and mult // matrices ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax, - sparse_mult, sparse_scatter); + sparse_mult, sparse_scatter); // Got everything we need, store it. init(in_gmin, in_gmax, sparse_mult, sparse_scatter); @@ -644,8 +640,9 @@ ScattDataHistogram::combine(const std::vector& those_scatts, //============================================================================== void -ScattDataTabular::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, - const double_2dvec& in_mult, const double_3dvec& coeffs) +ScattDataTabular::init(const xt::xtensor& in_gmin, + const xt::xtensor& in_gmax, const double_2dvec& in_mult, + const double_3dvec& coeffs) { int groups = coeffs.size(); int order = coeffs[0][0].size(); @@ -654,19 +651,13 @@ ScattDataTabular::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, double_3dvec matrix = coeffs; // Build the angular distribution mu values - mu = double_1dvec(order); + mu = xt::linspace(-1., 1., order); dmu = 2. / (order - 1); - mu[0] = -1.; - for (int imu = 1; imu < order - 1; imu++) { - mu[imu] = -1. + imu * dmu; - } - mu[order - 1] = 1.; // Get the scattering cross section value by integrating the distribution // over all mu points and then combining over all outgoing groups - scattxs.resize(groups); + scattxs = xt::zeros({groups}); for (int gin = 0; gin < groups; gin++) { - scattxs[gin] = 0.; for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) { for (int imu = 1; imu < order; imu++) { scattxs[gin] += 0.5 * dmu * (matrix[gin][i_gout][imu - 1] + @@ -743,7 +734,7 @@ ScattDataTabular::calc_f(int gin, int gout, double mu) int imu; if (mu == 1.) { // use size -2 to have the index one before the end - imu = this->mu.size() - 2; + imu = this->mu.shape()[0] - 2; } else { imu = std::floor((mu + 1.) / dmu + 1.) - 1; } @@ -764,7 +755,7 @@ ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt) sample_energy(gin, gout, i_gout); // Determine the outgoing cosine bin - int NP = this->mu.size(); + int NP = this->mu.shape()[0]; double xi = prn(); double c_k = dist[gin][i_gout][0]; @@ -804,21 +795,20 @@ ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt) //============================================================================== -double_3dvec +xt::xtensor ScattDataTabular::get_matrix(int max_order) { // Get the sizes and initialize the data to 0 int groups = energy.size(); // We ignore the requested order for Histogram and Tabular representations int order_dim = get_order(); - double_3dvec matrix = double_3dvec(groups, double_2dvec(groups, - double_1dvec(order_dim, 0.))); + xt::xtensor matrix({groups, groups, order_dim}, 0.); for (int gin = 0; gin < groups; gin++) { for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) { int gout = i_gout + gmin[gin]; for (int l = 0; l < order_dim; l++) { - matrix[gin][gout][l] = scattxs[gin] * energy[gin][i_gout] * + matrix(gin, gout, l) = scattxs[gin] * energy[gin][i_gout] * fmu[gin][i_gout][l]; } } @@ -847,8 +837,8 @@ ScattDataTabular::combine(const std::vector& those_scatts, int groups = those_scatts[0] -> energy.size(); - int_1dvec in_gmin(groups); - int_1dvec in_gmax(groups); + xt::xtensor in_gmin({groups}); + xt::xtensor in_gmax({groups}); double_3dvec sparse_scatter(groups); double_2dvec sparse_mult(groups); @@ -856,7 +846,7 @@ ScattDataTabular::combine(const std::vector& those_scatts, // so we use a base class method to sum up xs and create new energy and mult // matrices ScattData::base_combine(max_order, those_scatts, scalars, in_gmin, in_gmax, - sparse_mult, sparse_scatter); + sparse_mult, sparse_scatter); // Got everything we need, store it. init(in_gmin, in_gmax, sparse_mult, sparse_scatter); @@ -885,13 +875,8 @@ convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab, tab.scattxs = leg.scattxs; // Build mu and dmu - tab.mu = double_1dvec(n_mu); + tab.mu = xt::linspace(-1., 1., n_mu); tab.dmu = 2. / (n_mu - 1); - tab.mu[0] = -1.; - for (int imu = 1; imu < n_mu - 1; imu++) { - tab.mu[imu] = -1. + imu * tab.dmu; - } - tab.mu[n_mu - 1] = 1.; // Calculate f(mu) and integrate it so we can avoid rejection sampling int groups = tab.energy.size(); diff --git a/src/xsdata.cpp b/src/xsdata.cpp index d9863694a7..0a48dce2b5 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -5,6 +5,9 @@ #include #include +#include "xtensor/xbuilder.hpp" +#include "xtensor/xview.hpp" + #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/math_functions.h" @@ -20,7 +23,7 @@ namespace openmc { XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable, int scatter_format, int n_pol, int n_azi) { - int n_ang = n_pol * n_azi; + size_t n_ang = n_pol * n_azi; // check to make sure scatter format is OK before we allocate if (scatter_format != ANGLE_HISTOGRAM && scatter_format != ANGLE_TABULAR && @@ -28,31 +31,33 @@ XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable, fatal_error("Invalid scatter_format!"); } // allocate all [temperature][phi][theta][in group] quantities - total = double_2dvec(n_ang, double_1dvec(energy_groups, 0.)); - absorption = double_2dvec(n_ang, double_1dvec(energy_groups, 0.)); - inverse_velocity = double_2dvec(n_ang, double_1dvec(energy_groups, 0.)); + std::vector shape = {n_ang, energy_groups}; + total = xt::zeros(shape); + absorption = xt::zeros(shape); + inverse_velocity = xt::zeros(shape); if (fissionable) { - fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.)); - nu_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.)); - prompt_nu_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.)); - kappa_fission = double_2dvec(n_ang, double_1dvec(energy_groups, 0.)); + fission = xt::zeros(shape); + nu_fission = xt::zeros(shape); + prompt_nu_fission = xt::zeros(shape); + kappa_fission = xt::zeros(shape); } // allocate decay_rate; [temperature][phi][theta][delayed group] - decay_rate = double_2dvec(n_ang, double_1dvec(num_delayed_groups, 0.)); + shape[1] = num_delayed_groups; + decay_rate = xt::zeros(shape); if (fissionable) { + shape = {n_ang, energy_groups, num_delayed_groups}; // allocate delayed_nu_fission; [temperature][phi][theta][in group][delay group] - delayed_nu_fission = double_3dvec(n_ang, double_2dvec(energy_groups, - double_1dvec(num_delayed_groups, 0.))); + delayed_nu_fission = xt::zeros(shape); - // chi_prompt; [temperature][phi][theta][in group][delayed group] - chi_prompt = double_3dvec(n_ang, double_2dvec(energy_groups, - double_1dvec(energy_groups, 0.))); + // chi_prompt; [temperature][phi][theta][in group][out group] + shape = {n_ang, energy_groups, energy_groups}; + chi_prompt = xt::zeros(shape); // chi_delayed; [temperature][phi][theta][in group][out group][delay group] - chi_delayed = double_4dvec(n_ang, double_3dvec(energy_groups, - double_2dvec(energy_groups, double_1dvec(num_delayed_groups, 0.)))); + shape = {n_ang, energy_groups, energy_groups, num_delayed_groups}; + chi_delayed = xt::zeros(shape); } @@ -79,8 +84,8 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, { // Reconstruct the dimension information so it doesn't need to be passed int n_ang = n_pol * n_azi; - int energy_groups = total[0].size(); - int delayed_groups = decay_rate[0].size(); + int energy_groups = total.shape()[1]; + int delayed_groups = decay_rate.shape()[1]; // Set the fissionable-specific data if (fissionable) { @@ -100,7 +105,7 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, // denominator in tally methods for (int a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { - if (absorption[a][gin] == 0.) absorption[a][gin] = 1.e-10; + if (absorption(a, gin) == 0.) absorption(a, gin) = 1.e-10; } } @@ -110,7 +115,7 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, } else { for (int a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { - total[a][gin] = absorption[a][gin] + scatter[a]->scattxs[gin]; + total(a, gin) = absorption(a, gin) + scatter[a]->scattxs[gin]; } } } @@ -118,7 +123,7 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, // Fix if total is 0, since it is in the denominator when tallying for (int a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { - if (total[a][gin] == 0.) total[a][gin] = 1.e-10; + if (total(a, gin) == 0.) total(a, gin) = 1.e-10; } } } @@ -129,14 +134,13 @@ void XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups, int delayed_groups, bool is_isotropic) { - int n_ang = n_pol * n_azi; + size_t n_ang = n_pol * n_azi; // Get the fission and kappa_fission data xs; these are optional read_nd_vector(xsdata_grp, "fission", fission); read_nd_vector(xsdata_grp, "kappa-fission", kappa_fission); // Set/get beta - double_3dvec temp_beta =double_3dvec(n_ang, double_2dvec(energy_groups, - double_1dvec(delayed_groups, 0.))); + 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"); int ndims = dataset_ndims(xsdata); @@ -146,7 +150,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, if (ndims == 3) { // Beta is input as [delayed group] - double_1dvec temp_arr(n_pol * n_azi * delayed_groups); + std::vector temp_arr({n_pol * n_azi * delayed_groups}); read_nd_vector(xsdata_grp, "beta", temp_arr); // Broadcast to all incoming groups @@ -154,9 +158,9 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, for (int a = 0; a < n_ang; a++) { for (int dg = 0; dg < delayed_groups; dg++) { // Set the first group index and copy the rest - temp_beta[a][0][dg] = temp_arr[temp_idx++]; + temp_beta(a, 0, dg) = temp_arr[temp_idx++]; for (int gin = 1; gin < energy_groups; gin++) { - temp_beta[a][gin] = temp_beta[a][0]; + temp_beta(a, gin, dg) = temp_beta(a, 0, dg); } } } @@ -170,29 +174,33 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, // If chi is provided, set chi-prompt and chi-delayed if (object_exists(xsdata_grp, "chi")) { - double_2dvec temp_arr(n_ang, double_1dvec(energy_groups)); + xt::xtensor temp_arr ({n_ang, energy_groups}); read_nd_vector(xsdata_grp, "chi", temp_arr); for (int a = 0; a < n_ang; a++) { // First set the first group for (int gout = 0; gout < energy_groups; gout++) { - chi_prompt[a][0][gout] = temp_arr[a][gout]; + chi_prompt(a, 0, gout) = temp_arr(a, gout); } // Now normalize this data - double chi_sum = std::accumulate(chi_prompt[a][0].begin(), - chi_prompt[a][0].end(), - 0.); + double chi_sum = 0.; + for (int 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 (int gout = 0; gout < energy_groups; gout++) { - chi_prompt[a][0][gout] /= chi_sum; + chi_prompt(a, 0, gout) /= chi_sum; } // And extend to the remaining incoming groups for (int gin = 1; gin < energy_groups; gin++) { - chi_prompt[a][gin] = chi_prompt[a][0]; + for (int gout = 0; gout < energy_groups; gout++) { + chi_prompt(a, gin, gout) = chi_prompt(a, 0, gout); + } } // Finally set chi-delayed equal to chi-prompt @@ -200,8 +208,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, for(int gin = 0; gin < energy_groups; gin++) { for (int gout = 0; gout < energy_groups; gout++) { for (int dg = 0; dg < delayed_groups; dg++) { - chi_delayed[a][gin][gout][dg] = - chi_prompt[a][gin][gout]; + chi_delayed(a, gin, gout, dg) = chi_prompt(a, gin, gout); } } } @@ -224,15 +231,18 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, for (int a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { for (int dg = 0; dg < delayed_groups; dg++) { - delayed_nu_fission[a][gin][dg] = - temp_beta[a][gin][dg] * prompt_nu_fission[a][gin]; + 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 = std::accumulate(temp_beta[a][gin].begin(), - temp_beta[a][gin].end(), 0.); - prompt_nu_fission[a][gin] *= (1. - beta_sum); + double beta_sum = 0.; + for (int gin = 0; gin < energy_groups; gin++) { + beta_sum += temp_beta(a, gin); + } + + prompt_nu_fission(a, gin) *= (1. - beta_sum); } } } @@ -244,14 +254,17 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, // Normalize the chi info so the CDF is 1. for (int a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { - double chi_sum = std::accumulate(chi_prompt[a][gin].begin(), - chi_prompt[a][gin].end(), 0.); + double chi_sum = 0.; + for (int 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; + prompt_nu_fission(a, gin) = chi_sum; if (chi_sum >= 0.) { for (int gout = 0; gout < energy_groups; gout++) { - chi_prompt[a][gin][gout] /= chi_sum; + chi_prompt(a, gin, gout) /= chi_sum; } } else { fatal_error("Encountered chi for a group that is <= 0!"); @@ -262,8 +275,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, for (int gin = 0; gin < energy_groups; gin++) { for (int gout = 0; gout < energy_groups; gout++) { for (int dg = 0; dg < delayed_groups; dg++) { - chi_delayed[a][gin][gout][dg] = - chi_prompt[a][gin][gout]; + chi_delayed(a, gin, gout, dg) = chi_prompt(a, gin, gout); } } } @@ -271,16 +283,17 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, // Set the delayed-nu-fission and correct prompt-nu-fission with beta for (int gin = 0; gin < energy_groups; gin++) { for (int dg = 0; dg < delayed_groups; dg++) { - delayed_nu_fission[a][gin][dg] = - temp_beta[a][gin][dg] * - prompt_nu_fission[a][gin]; + 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 = std::accumulate(temp_beta[a][gin].begin(), - temp_beta[a][gin].end(), 0.); - prompt_nu_fission[a][gin] *= (1. - beta_sum); + double beta_sum = 0.; + for (int dg = 0; dg < delayed_groups; dg++) { + beta_sum += temp_beta(a, gin, dg); + } + prompt_nu_fission(a, gin) *= (1. - beta_sum); } } } @@ -293,21 +306,24 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, // If chi-prompt is provided, set chi-prompt if (object_exists(xsdata_grp, "chi-prompt")) { - double_2dvec temp_arr(n_ang, double_1dvec(energy_groups)); + xt::xtensor temp_arr({n_ang, energy_groups}); read_nd_vector(xsdata_grp, "chi-prompt", temp_arr); for (int a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { for (int gout = 0; gout < energy_groups; gout++) { - chi_prompt[a][gin][gout] = temp_arr[a][gout]; + chi_prompt(a, gin, gout) = temp_arr(a, gout); } // Normalize chi so its CDF goes to 1 - double chi_sum = std::accumulate(chi_prompt[a][gin].begin(), - chi_prompt[a][gin].end(), 0.); + double chi_sum = 0.; + for (int gout = 0; gout < energy_groups; gout++) { + chi_sum += chi_prompt(a, gin, gout); + } + if (chi_sum >= 0.) { for (int gout = 0; gout < energy_groups; gout++) { - chi_prompt[a][gin][gout] /= chi_sum; + chi_prompt(a, gin, gout) /= chi_sum; } } else { fatal_error("Encountered chi-prompt for a group that is <= 0.!"); @@ -326,13 +342,16 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, if (ndims == 3) { // chi-delayed is a [in group] vector - double_2dvec temp_arr(n_ang, double_1dvec(energy_groups)); + xt::xtensor temp_arr({n_ang, energy_groups}); read_nd_vector(xsdata_grp, "chi-delayed", temp_arr); for (int a = 0; a < n_ang; a++) { // normalize the chi CDF to 1 - double chi_sum = std::accumulate(temp_arr[a].begin(), - temp_arr[a].end(), 0.); + double chi_sum = 0.; + for (int 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!"); } @@ -341,7 +360,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, for (int gin = 0; gin < energy_groups; gin++) { for (int gout = 0; gout < energy_groups; gout++) { for (int dg = 0; dg < delayed_groups; dg++) { - chi_delayed[a][gin][gout][dg] = temp_arr[a][gout] / chi_sum; + chi_delayed(a, gin, gout, dg) = temp_arr(a, gout) / chi_sum; } } } @@ -356,12 +375,12 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, for (int gin = 0; gin < energy_groups; gin++) { double chi_sum = 0.; for (int gout = 0; gout < energy_groups; gout++) { - chi_sum += chi_delayed[a][gin][gout][dg]; + chi_sum += chi_delayed(a, gin, gout, dg); } if (chi_sum > 0.) { for (int gout = 0; gout < energy_groups; gout++) { - chi_delayed[a][gin][gout][dg] /= chi_sum; + chi_delayed(a, gin, gout, dg) /= chi_sum; } } else { fatal_error("Encountered chi-delayed for a group that is <= 0!"); @@ -384,29 +403,30 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, if (ndims == 3) { // prompt-nu-fission is a [in group] vector - read_nd_vector(xsdata_grp, "prompt-nu-fission", - prompt_nu_fission); + 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 - double_3dvec temp_arr(n_ang, double_2dvec(energy_groups, - double_1dvec(energy_groups))); + xt::xtensor temp_arr({n_ang, energy_groups, energy_groups}); read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_arr); // The prompt_nu_fission vector from the matrix form for (int a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { - double prompt_sum = std::accumulate(temp_arr[a][gin].begin(), - temp_arr[a][gin].end(), 0.); - prompt_nu_fission[a][gin] = prompt_sum; + double prompt_sum = 0.; + for (int 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 (int gin= 0; gin < energy_groups; gin++) { - if (prompt_nu_fission[a][gin] > 0.) { + if (prompt_nu_fission(a, gin) > 0.) { for (int gout = 0; gout < energy_groups; gout++) { - chi_prompt[a][gin][gout] = - temp_arr[a][gin][gout] / prompt_nu_fission[a][gin]; + 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!"); @@ -429,19 +449,19 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, if (ndims == 3) { // delayed-nu-fission is an [in group] vector - if (temp_beta[0][0][0] == 0.) { + if (temp_beta(0, 0, 0) == 0.) { fatal_error("cannot set delayed-nu-fission with a 1D array if " "beta is not provided"); } - double_2dvec temp_arr(n_ang, double_1dvec(energy_groups)); + xt::xtensor temp_arr({n_ang, energy_groups}); read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr); for (int a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { for (int 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]; + delayed_nu_fission(a, gin, dg) = + temp_beta(a, gin, dg) * temp_arr(a, gin); } } } @@ -451,9 +471,9 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, delayed_nu_fission); } else if (ndims == 5) { - // This will contain delayed-nu-fision and chi-delayed data - double_4dvec temp_arr(n_ang, double_3dvec(energy_groups, - double_2dvec(energy_groups, double_1dvec(delayed_groups)))); + // This will contain delayed-nu-fission and chi-delayed data + xt::xtensor temp_arr({n_ang, energy_groups, energy_groups, + delayed_groups}); read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr); // Set the 3D delayed-nu-fission matrix and 4D chi-delayed matrix @@ -463,14 +483,14 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, for (int gin = 0; gin < energy_groups; gin++) { double gout_sum = 0.; for (int 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]; + 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; + delayed_nu_fission(a, gin, dg) = gout_sum; // Normalize chi-delayed if (gout_sum > 0.) { for (int gout = 0; gout < energy_groups; gout++) { - chi_delayed[a][gin][gout][dg] /= gout_sum; + chi_delayed(a, gin, gout, dg) /= gout_sum; } } else { fatal_error("Encountered chi-delayed for a group that is <= 0!"); @@ -488,10 +508,10 @@ XsData::fission_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 a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { - nu_fission[a][gin] = - std::accumulate(delayed_nu_fission[a][gin].begin(), - delayed_nu_fission[a][gin].end(), - prompt_nu_fission[a][gin]); + nu_fission(a, gin) = prompt_nu_fission(a, gin); + for (int dg = 0; dg < delayed_groups; dg++) { + nu_fission(a, gin) += delayed_nu_fission(a, gin, dg); + } } } } @@ -503,94 +523,92 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups, int scatter_format, int final_scatter_format, int order_data, int max_order, int legendre_to_tabular_points) { - int n_ang = n_pol * n_azi; + size_t n_ang = n_pol * n_azi; if (!object_exists(xsdata_grp, "scatter_data")) { fatal_error("Must provide scatter_data group!"); } hid_t scatt_grp = open_group(xsdata_grp, "scatter_data"); // Get the outgoing group boundary indices - int_2dvec gmin(n_ang, int_1dvec(energy_groups)); + xt::xtensor gmin({n_ang, energy_groups}); read_nd_vector(scatt_grp, "g_min", gmin, true); - int_2dvec gmax(n_ang, int_1dvec(energy_groups)); + xt::xtensor gmax({n_ang, energy_groups}); read_nd_vector(scatt_grp, "g_max", gmax, true); // Make gmin and gmax start from 0 vice 1 as they do in the library - for (int a = 0; a < n_ang; a++) { - for (int gin = 0; gin < energy_groups; gin++) { - gmin[a][gin] -= 1; - gmax[a][gin] -= 1; - } - } + gmin -= 1; + gmax -= 1; // Now use this info to find the length of a vector to hold the flattened // data. int length = 0; for (int a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { - length += order_data * (gmax[a][gin] - gmin[a][gin] + 1); + length += order_data * (gmax(a, gin) - gmin(a, gin) + 1); } } - double_1dvec temp_arr(length); - read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true); - // Compare the number of orders given with the max order of the problem; - // strip off the superfluous orders if needed - int order_dim; - if (scatter_format == ANGLE_LEGENDRE) { - order_dim = std::min(order_data - 1, max_order) + 1; - } else { - order_dim = order_data; - } - - // convert the flattened temp_arr to a jagged array for passing to - // scatt data double_4dvec input_scatt(n_ang, double_3dvec(energy_groups)); + //temp_arr scope + { + std::vector temp_arr(length); + read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true); - int temp_idx = 0; - for (int a = 0; a < n_ang; a++) { - for (int gin = 0; gin < energy_groups; gin++) { - input_scatt[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1); - for (int i_gout = 0; i_gout < input_scatt[a][gin].size(); i_gout++) { - input_scatt[a][gin][i_gout].resize(order_dim); - for (int l = 0; l < order_dim; l++) { - input_scatt[a][gin][i_gout][l] = temp_arr[temp_idx++]; + // Compare the number of orders given with the max order of the problem; + // strip off the superfluous orders if needed + int order_dim; + if (scatter_format == ANGLE_LEGENDRE) { + order_dim = std::min(order_data - 1, max_order) + 1; + } else { + order_dim = order_data; + } + + // convert the flattened temp_arr to a jagged array for passing to + // scatt data + int temp_idx = 0; + for (int a = 0; a < n_ang; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + input_scatt[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1); + for (int i_gout = 0; i_gout < input_scatt[a][gin].size(); i_gout++) { + input_scatt[a][gin][i_gout].resize(order_dim); + for (int l = 0; l < order_dim; l++) { + input_scatt[a][gin][i_gout][l] = temp_arr[temp_idx++]; + } + // Adjust index for the orders we didnt take + temp_idx += (order_data - order_dim); } - // Adjust index for the orders we didnt take - temp_idx += (order_data - order_dim); } } } - temp_arr.clear(); // Get multiplication matrix double_3dvec temp_mult(n_ang, double_2dvec(energy_groups)); if (object_exists(scatt_grp, "multiplicity_matrix")) { - temp_arr.resize(length / order_data); + std::vector temp_arr(length / order_data); read_nd_vector(scatt_grp, "multiplicity_matrix", temp_arr); // convert the flat temp_arr to a jagged array for passing to scatt data int temp_idx = 0; for (int a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { - temp_mult[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1); + temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1); for (int i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) { temp_mult[a][gin][i_gout] = temp_arr[temp_idx++]; } } } + temp_arr.clear(); } else { // Use a default: multiplicities are 1.0. for (int a = 0; a < n_ang; a++) { for (int gin = 0; gin < energy_groups; gin++) { - temp_mult[a][gin].resize(gmax[a][gin] - gmin[a][gin] + 1); + temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1); for (int i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) { temp_mult[a][gin][i_gout] = 1.; } } } } - temp_arr.clear(); close_group(scatt_grp); // Finally, convert the Legendre data to tabular, if needed @@ -598,7 +616,10 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, final_scatter_format == ANGLE_TABULAR) { for (int a = 0; a < n_ang; a++) { ScattDataLegendre legendre_scatt; - legendre_scatt.init(gmin[a], gmax[a], temp_mult[a], input_scatt[a]); + xt::xtensor in_gmin = xt::view(gmin, a, xt::all()); + xt::xtensor in_gmax = xt::view(gmax, a, xt::all()); + legendre_scatt.init(in_gmin, in_gmax, + temp_mult[a], input_scatt[a]); // Now create a tabular version of legendre_scatt convert_legendre_to_tabular(legendre_scatt, @@ -611,7 +632,9 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, // We are sticking with the current representation // Initialize the ScattData object with this data for (int a = 0; a < n_ang; a++) { - scatter[a]->init(gmin[a], gmax[a], temp_mult[a], input_scatt[a]); + scatter[a]->init(xt::view(gmin, a, xt::all()), + xt::view(gmax, a, xt::all()), + temp_mult[a], input_scatt[a]); } } } @@ -627,70 +650,25 @@ XsData::combine(const std::vector& those_xs, XsData* that = those_xs[i]; if (!equiv(*that)) fatal_error("Cannot combine the XsData objects!"); double scalar = scalars[i]; - for (int a = 0; a < total.size(); a++) { - for (int gin = 0; gin < total[a].size(); gin++) { - total[a][gin] += scalar * that->total[a][gin]; - absorption[a][gin] += scalar * that->absorption[a][gin]; - if (i == 0) { - inverse_velocity[a][gin] = that->inverse_velocity[a][gin]; - } - if (that->prompt_nu_fission.size() > 0) { - nu_fission[a][gin] += scalar * that->nu_fission[a][gin]; - prompt_nu_fission[a][gin] += - scalar * that->prompt_nu_fission[a][gin]; - kappa_fission[a][gin] += scalar * that->kappa_fission[a][gin]; - fission[a][gin] += scalar * that->fission[a][gin]; - - for (int dg = 0; dg < delayed_nu_fission[a][gin].size(); dg++) { - delayed_nu_fission[a][gin][dg] += - scalar * that->delayed_nu_fission[a][gin][dg]; - } - - for (int gout = 0; gout < chi_prompt[a][gin].size(); gout++) { - chi_prompt[a][gin][gout] += - scalar * that->chi_prompt[a][gin][gout]; - - for (int dg = 0; dg < chi_delayed[a][gin][gout].size(); dg++) { - chi_delayed[a][gin][gout][dg] += - scalar * that->chi_delayed[a][gin][gout][dg]; - } - } - } - } - - for (int dg = 0; dg < decay_rate[a].size(); dg++) { - decay_rate[a][dg] += scalar * that->decay_rate[a][dg]; - } - - // Normalize chi - if (chi_prompt.size() > 0) { - for (int gin = 0; gin < chi_prompt[a].size(); gin++) { - double norm = std::accumulate(chi_prompt[a][gin].begin(), - chi_prompt[a][gin].end(), 0.); - if (norm > 0.) { - for (int gout = 0; gout < chi_prompt[a][gin].size(); gout++) { - chi_prompt[a][gin][gout] /= norm; - } - } - - for (int dg = 0; dg < chi_delayed[a][gin][0].size(); dg++) { - norm = 0.; - for (int gout = 0; gout < chi_delayed[a][gin].size(); gout++) { - norm += chi_delayed[a][gin][gout][dg]; - } - if (norm > 0.) { - for (int gout = 0; gout < chi_delayed[a][gin].size(); gout++) { - chi_delayed[a][gin][gout][dg] /= norm; - } - } - } - } - } + total += scalar * that->total; + absorption += scalar * that->absorption; + if (i == 0) { + inverse_velocity = that->inverse_velocity; } + if (that->prompt_nu_fission.shape()[0] > 0) { + nu_fission += scalar * that->nu_fission; + prompt_nu_fission += scalar * that->prompt_nu_fission; + kappa_fission += scalar * that->kappa_fission; + fission += scalar * that->fission; + delayed_nu_fission += scalar * that->delayed_nu_fission; + chi_prompt += scalar * that->chi_prompt; + chi_delayed += scalar * that->chi_delayed; + } + decay_rate += scalar * that->decay_rate; } // Allow the ScattData object to combine itself - for (int a = 0; a < total.size(); a++) { + for (int a = 0; a < total.shape()[0]; a++) { // Build vector of the scattering objects to incorporate std::vector those_scatts(those_xs.size()); for (int i = 0; i < those_xs.size(); i++) { @@ -707,8 +685,7 @@ XsData::combine(const std::vector& those_xs, bool XsData::equiv(const XsData& that) { - return ((absorption.size() == that.absorption.size()) && - (absorption[0].size() == that.absorption[0].size())); + return (absorption.shape() == that.absorption.shape()); } } //namespace openmc diff --git a/tests/regression_tests/mgxs_library_ce_to_mg/test.py b/tests/regression_tests/mgxs_library_ce_to_mg/test.py index 72f052e68e..4735961986 100644 --- a/tests/regression_tests/mgxs_library_ce_to_mg/test.py +++ b/tests/regression_tests/mgxs_library_ce_to_mg/test.py @@ -71,10 +71,11 @@ class MGXSTestHarness(PyAPITestHarness): openmc.run(openmc_exec=config['exe']) def _cleanup(self): - super()._cleanup() - f = 'mgxs.h5' - if os.path.exists(f): - os.remove(f) + pass + # super()._cleanup() + # f = 'mgxs.h5' + # if os.path.exists(f): + # os.remove(f) def test_mgxs_library_ce_to_mg():