diff --git a/src/constants.h b/src/constants.h new file mode 100644 index 000000000..b14984a14 --- /dev/null +++ b/src/constants.h @@ -0,0 +1,60 @@ + +#ifndef CONSTANTS_H +#define CONSTANTS_H + +#include +#include + +namespace openmc { + +typedef std::array dir_arr; +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; + +int constexpr MAX_SAMPLE {10000}; + +constexpr std::array VERSION {0, 10, 0}; +constexpr std::array VERSION_PARTICLE_RESTART {2, 0}; + +// Maximum number of words in a single line, length of line, and length of +// single word +constexpr int MAX_WORDS {500}; +constexpr int MAX_LINE_LEN {250}; +constexpr int MAX_WORD_LEN {150}; +constexpr int MAX_FILE_LEN {255}; + +// Physical Constants +constexpr double K_BOLTZMANN {8.6173303e-5}; // Boltzmann constant in eV/K + +// Angular distribution type +constexpr int ANGLE_ISOTROPIC {1}; +constexpr int ANGLE_32_EQUI {2}; +constexpr int ANGLE_TABULAR {3}; +constexpr int ANGLE_LEGENDRE {4}; +constexpr int ANGLE_HISTOGRAM {5}; + +// MGXS Table Types +constexpr int MGXS_ISOTROPIC {1}; // Isotroically weighted data +constexpr int MGXS_ANGLE {2}; // Data by angular bins + +// Flag to denote this was a macroscopic data object +constexpr double MACROSCOPIC_AWR {-2.}; + +// Number of mu bins to use when converting Legendres to tabular type +constexpr int DEFAULT_NMU {33}; + +// Temperature treatment method +constexpr int TEMPERATURE_NEAREST {1}; +constexpr int TEMPERATURE_INTERPOLATION {2}; + + +} // namespace openmc + +#endif // CONSTANTS_H \ No newline at end of file diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index 8a8391bb5..3d2f3ee7b 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -447,6 +447,152 @@ read_complex(hid_t obj_id, const char* name, double _Complex* buffer, bool indep } +void +read_nd_vector(hid_t obj_id, const char* name, std::vector& result, + bool must_have) +{ + if (object_exists(obj_id, name)) { + read_double(obj_id, name, &result[0], 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, + std::vector >& result, bool must_have) +{ + if (object_exists(obj_id, name)) { + int dim1 = result.size(); + int dim2 = result[0].size(); + std::vector temp_arr = std::vector(dim1 * dim2); + read_double(obj_id, name, &temp_arr[0], 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++]; + } + } + } else if (must_have) { + fatal_error(std::string("Must provide " + std::string(name) + "!")); + } +} + + +void +read_nd_vector(hid_t obj_id, const char* name, + std::vector > >& result, + bool must_have) +{ + if (object_exists(obj_id, name)) { + dim1 = result.size(); + dim2 = result[0].size(); + dim3 = result[0][0].size(); + std::vector temp_arr = std::vector(dim1 * dim2 * dim3); + read_double(obj_id, name, &temp_arr[0], true); + + int temp_idx = 0; + 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++]; + } + } + } + } else if (must_have) { + fatal_error(std::string("Must provide " + std::string(name) + "!")); + } +} + +void +read_nd_vector(hid_t obj_id, const char* name, + std::vector > >& result, + bool must_have) +{ + if (object_exists(obj_id, name)) { + dim1 = result.size(); + dim2 = result[0].size(); + dim3 = result[0][0].size(); + std::vector temp_arr = std::vector(dim1 * dim2 * dim3); + read_int(obj_id, name, &temp_arr[0], true); + + int temp_idx = 0; + 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++]; + } + } + } + } else if (must_have) { + fatal_error(std::string("Must provide " + std::string(name) + "!")); + } +} + +void +read_nd_vector(hid_t obj_id, const char* name, + std::vector > > >& result, + bool must_have) +{ + if (object_exists(obj_id, name)) { + dim1 = result.size(); + dim2 = result[0].size(); + dim3 = result[0][0].size(); + dim4 = result[0][0][0].size(); + std::vector temp_arr = std::vector( + dim1 * dim2 * dim3 * dim4); + read_double(obj_id, name, &temp_arr[0], true); + + int temp_idx = 0; + for (int i = 0; i < dim1; i++) { + 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++]; + } + } + } + } + } else if (must_have) { + fatal_error(std::string("Must provide " + std::string(name) + "!")); + } +} + +void +read_nd_vector(hid_t obj_id, const char* name, + std::vector > > > >& result, + bool must_have) +{ + if (object_exists(obj_id, name)) { + dim1 = result.size(); + dim2 = result[0].size(); + dim3 = result[0][0].size(); + dim4 = result[0][0][0].size(); + dim5 = result[0][0][0][0].size(); + std::vector temp_arr = std::vector( + dim1 * dim2 * dim3 * dim4 * dim5); + read_double(obj_id, name, &temp_arr[0], true); + + int temp_idx = 0; + for (int i = 0; i < dim1; i++) { + for (int j = 0; j < dim2; j++) { + 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++]; + } + } + } + } + } + } else if (must_have) { + fatal_error(std::string("Must provide " + std::string(name) + "!")); + } +} + + void read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results) { diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h index e38a31e99..734cf0f91 100644 --- a/src/hdf5_interface.h +++ b/src/hdf5_interface.h @@ -55,6 +55,35 @@ extern "C" void read_string(hid_t obj_id, const char* name, size_t slen, extern "C" void read_complex(hid_t obj_id, const char* name, double _Complex* buffer, bool indep); +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, + 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, + 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, + bool must_have = false); + extern "C" void read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, double* results); diff --git a/src/mgxs.h b/src/mgxs.h new file mode 100644 index 000000000..8ed4ba99e --- /dev/null +++ b/src/mgxs.h @@ -0,0 +1,95 @@ +//! \file mgxs.h +//! A collection of classes for Multi-Group Cross Section data + +#ifndef MGXS_H +#define MGXS_H + +#include +#include +#include +#include +#include +#include + +#include "constants.h" +#include "hdf5_interface.h" +#include "math_functions.h" +#include "random_lcg.h" +#include "scattdata.h" +#include "string_functions.h" +#include "xsdata.h" + + +namespace openmc { + + +//============================================================================== +// MGXS contains the mgxs data for a nuclide/material +//============================================================================== + +class Mgxs { + private: + std::string name; // name of dataset, e.g., UO2 + double awr; // atomic weight ratio + double_1dvec kTs; // temperature in eV (k * T) + int scatter_format; // flag for if this is legendre, histogram, or tabular + int num_delayed_groups; // number of delayed neutron groups + int num_groups; // number of energy groups + int index_temp; // cache of temperature index + double last_sqrtkT; // cache of the temperature corresponding to index_temp + std::vector xs; // Cross section data + int n_pol; + int n_azi; + int index_pol; // cache fof the angle indices + int index_azi; + double_1dvec polar; + double_1dvec azimuthal; + dir_arr last_uvw; + void _metadata_from_hdf5(hid_t xs_id, int in_num_groups, + int in_num_delayed_groups, double_1dvec temperature, int& method, + double tolerance, double_1dvec& temps_to_read, int& order_dim); + + public: + bool fissionable; // Is this fissionable + void init(const std::string& in_name, double in_awr, double_1dvec& in_kTs, + bool in_fissionable, int in_scatter_format, int in_num_groups, + int in_num_delayed_groups, double_1dvec& in_polar, + double_1dvec& in_azimuthal); + void build_macro(const std::string& in_name, double_1dvec& mat_kTs, + std::vector& micros, double_1dvec& atom_densities, + int& method, double tolerance); + void combine(std::vector& micros, double_1dvec& scalars, + int_1dvec& micro_ts, int this_t); + void from_hdf5(hid_t xs_id, int energy_groups, int delayed_groups, + double_1dvec temperature, int& method, + double tolerance, int max_order, + bool legendre_to_tabular, + int legendre_to_tabular_points); + double get_xs(const char* xstype, int gin, int* gout, double* mu, + int* dg); + void sample_fission_energy(int gin, double nu_fission, int& dg, int& gout); + void sample_scatter(dir_arr& uvw, int gin, int& gout, double& mu, + double& wgt); + void calculate_xs(int gin, double sqrtkT, dir_arr& uvw, double& total_xs, + double& abs_xs, double& nu_fiss_xs); + bool equiv(const Mgxs& that); + inline void set_temperature_index(double sqrtkT); + inline void set_angle_index(dir_arr& uvw); +}; + +extern "C" void read_mgxs_library(hid_t file_id, int n_nuclides, char** names, + int energy_groups, int delayed_groups, int n_temps, double temps[], + int& method, double tolerance, int max_order, bool legendre_to_tabular, + int legendre_to_tabular_points); +extern "C" bool query_fissionable(const int i_nuclides[], const int n_nuclides); +void create_macro_xs(int n_materials, double_2dvec& mat_kTs, + std::vector& mat_names, double_1dvec& atom_densities, + int& method, double tolerance); + + +// Storage for the MGXS data +std::vector nuclides_MG; +std::vector macro_xs; + +} // namespace openmc +#endif // MGXS_H \ No newline at end of file diff --git a/src/scattdata.cpp b/src/scattdata.cpp new file mode 100644 index 000000000..8283febc8 --- /dev/null +++ b/src/scattdata.cpp @@ -0,0 +1,831 @@ +#include "scattdata.h" + +namespace openmc { + +//============================================================================== +// Methods for use by all extended types +//============================================================================== + + + +//============================================================================== +// ScattData base-class methods +//============================================================================== + +void ScattData::generic_init(int order, int_1dvec in_gmin, + int_1dvec in_gmax, double_2dvec in_energy, double_2dvec in_mult) +{ + int groups = in_energy.size(); + + gmin = in_gmin; + gmax = in_gmax; + energy.resize(groups); + mult.resize(groups); + dist.resize(groups); + + for (int gin = 0; gin < groups; gin++) { + // Make sure the energy is normalized + double norm = std::accumulate(in_energy[gin].begin(), + in_energy[gin].end(), 0.); + + if (norm != 0.) { + for (auto& n : in_energy[gin]) n /= norm; + } + + // Store the inputted data + energy[gin] = in_energy[gin]; + mult[gin] = in_mult[gin]; + + // Initialize the distribution data + dist[gin].resize(in_gmax[gin] - in_gmin[gin] + 1); + for (auto& v : dist[gin]) { + v.resize(order); + for (auto& n : v) n = 0.; + } + } +} + + +void ScattData::sample_energy(int gin, int& gout, int& i_gout) +{ + // Sample the outgoing group + double xi = prn(); + i_gout = 0; //TODO: + 1? + gout = gmin[gin]; + double prob = energy[gin][i_gout]; + while((prob < xi) && (gout < gmax[gin])) { + gout++; + i_gout++; + prob += energy[gin][i_gout]; + } +} + + +double ScattData::get_xs(const char* xstype, int gin, int* gout, double* mu) +{ + // Set the outgoing group offset index as needed + int i_gout = 0; + if (gout != nullptr) { + // short circuit the function if gout is from a zero portion of the + // scattering matrix + if ((*gout < gmin[gin]) || (*gout >= gmax[gin])) { // > gmax? + return 0.; + } + i_gout = *gout - gmin[gin]; + } + + double val = 0.; + if (std::strcmp(xstype, "scatter")) { + if (gout != nullptr) { + val = scattxs[gin] * energy[gin][i_gout]; + } else { + val = scattxs[gin]; + } + } else if (std::strcmp(xstype, "scatter/mult")) { + if (gout != nullptr) { + val = scattxs[gin] * energy[gin][i_gout] / mult[gin][i_gout]; + } else { + val = scattxs[gin] / std::inner_product(mult[gin].begin(), + mult[gin].end(), + energy[gin].begin(), 0.0); + } + } else if (std::strcmp(xstype, "scatter*f_mu/mult")) { + if ((gout != nullptr) && (mu != nullptr)) { + val = scattxs[gin] * energy[gin][i_gout] * calc_f(gin, *gout, *mu); + } else { + // This is not an expected path (asking for f_mu without asking for a + // group or mu is not useful + fatal_error("Invalid call to get_xs"); + } + } else if (std::strcmp(xstype, "scatter*f_mu")) { + if ((gout != nullptr) && (mu != nullptr)) { + val = scattxs[gin] * energy[gin][i_gout] * calc_f(gin, *gout, *mu) / + mult[gin][i_gout]; + } else { + // This is not an expected path (asking for f_mu without asking for a + // group or mu is not useful + fatal_error("Invalid call to get_xs"); + } + } + return val; +} + + +//============================================================================== +// ScattDataLegendre methods +//============================================================================== + +void ScattDataLegendre::init(int_1dvec in_gmin, int_1dvec in_gmax, + double_2dvec in_mult, double_3dvec coeffs) +{ + int groups = coeffs.size(); + int order = coeffs[0].size(); + + // make a copy of coeffs that we can use to both extract data and normalize + double_3dvec matrix = coeffs; + + // Get the scattering cross section value by summing the un-normalized P0 + // coefficient in the variable matrix over all outgoing groups. + scattxs.resize(groups); + for (int gin = 0; gin < groups; gin++) { + int num_groups = gmax[gin] - gmin[gin] + 1; + scattxs[gin] = 0.; + for (int i_gout = 0; i_gout < num_groups; i_gout++) { + scattxs[gin] = std::accumulate(matrix[gin][i_gout].begin(), + matrix[gin][i_gout].end(), + scattxs[gin]); + } + } + + // Build the energy transfer matrix from data in the variable matrix while + // also normalizing the variable matrix itself + // (forcing the CDF of f(mu=1) == 1) + double_2dvec in_energy; + in_energy.resize(groups); + for (int gin = 0; gin < groups; gin++) { + int num_groups = gmax[gin] - gmin[gin] + 1; + in_energy[gin].resize(num_groups); + for (int i_gout = 0; i_gout < num_groups; i_gout++) { + double norm = matrix[gin][i_gout][0]; + in_energy[gin][i_gout] = norm; + if (norm != 0.) { + for (auto& n : matrix[gin][i_gout]) n /= norm; + } + } + } + + // Initialize the base class attributes + ScattData::generic_init(order, in_gmin, in_gmax, in_energy, in_mult); + + // Set the distribution (sdata.dist) values and initialize max_val + max_val.resize(groups); + for (int gin = 0; gin < groups; gin++) { + int num_groups = gmax[gin] - gmin[gin] + 1; + for (int i_gout = 0; i_gout < num_groups; i_gout++) { + dist[gin][i_gout] = matrix[gin][i_gout]; + } + max_val[gin].resize(num_groups); + for (auto& n : max_val[gin]) n = 0.; + } + + // Now update the maximum value + update_max_val(); +} + + +void ScattDataLegendre::update_max_val() +{ + int groups = max_val.size(); + // Step through the polynomial with fixed number of points to identify the + // maximal value + int Nmu = 1001; + double dmu = 2. / (Nmu - 1); + for (int gin = 0; gin < groups; gin++) { + int num_groups = gmax[gin] - gmin[gin] + 1; + for (int i_gout = 0; i_gout < num_groups; i_gout++) { + for (int imu = 0; imu < Nmu; imu++) { + double mu; + if (imu == 0) { + mu = -1.; + } else if (imu == (Nmu - 1)) { + mu = 1.; + } else { + mu = -1. + (imu - 1) * dmu; + } + + // Calculate probability + double f = evaluate_legendre_c(dist[gin][i_gout].size() - 1, + dist[gin][i_gout].data(), mu); + + // if this is a new maximum, store it + if (f > max_val[gin][i_gout]) max_val[gin][i_gout] = f; + } // end imu loop + + // Since we may not have caught the true max, add 10% margin + max_val[gin][i_gout] *= 1.1; + } + } +} + + +double ScattDataLegendre::calc_f(int gin, int gout, double mu) +{ + // TODO: gout >= or gout >? + double f; + if ((gout < gmin[gin]) || (gout >= gmax[gin])) { + f = 0.; + } else { + // TODO: size() -1 or just size? + int i_gout = gout - gmin[gin]; //TODO: + 1? + f = evaluate_legendre_c(dist[gin][i_gout].size() - 1, + dist[gin][i_gout].data(), mu); + } + return f; +} + + +void ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt) +{ + // Sample the outgoing energy using the base-class method + int i_gout; + sample_energy(gin, gout, i_gout); + + // Now we can sample mu using the scattering kernel using rejection + // sampling from a rectangular bounding box + double M = max_val[gin][i_gout]; + int samples = 0; + + while(true) { + double mu = 2. * prn() - 1.; + double f = calc_f(gin, gout, mu); + if (f > 0.) { + double u = prn() * M; + if (u <= f) break; + } + samples++; + if (samples > MAX_SAMPLE) { + fatal_error("Maximum number of Legendre expansion samples reached"); + } + }; + + // Update the weight to reflect neutron multiplicity + wgt *= mult[gin][i_gout]; +} + + +void ScattDataLegendre::combine(std::vector those_scatts, + double_1dvec& scalars) +{ + int groups = energy.size(); + // Find the maximum order in the data set + int max_order = get_order(); + for (int i = 0; i < those_scatts.size(); i++) { + // Lets also make sure these items are combineable + ScattDataLegendre* that = dynamic_cast(those_scatts[i]); + if (!equiv(*that)) { + fatal_error("Cannot combine the ScattData objects!"); + } + int that_order = that->get_order(); + if (that_order > max_order) max_order = that_order; + } + max_order++; // Add one since this is a Legendre + + // Now allocate and zero our storage spaces + double_3dvec this_matrix = get_matrix(max_order); + double_2dvec mult_numer(groups, double_1dvec(groups, 0.)); + double_2dvec mult_denom(groups, double_1dvec(groups, 0.)); + + // Build the dense scattering and multiplicity matrices + // Get the multiplicity_matrix + // To combine from nuclidic data we need to use the final relationship + // mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / + // sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'})) + // Developed as follows: + // mult_{gg'} = nuScatt{g,g'} / Scatt{g,g'} + // mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / sum(N_i*scatt_{i,g,g'}) + // mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) / + // sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'})) + // nuscatt_{i,g,g'} can be reconstructed from the energy and scattxs member + // variables + for (int i = 0; i < those_scatts.size(); i++) { + ScattDataLegendre* that = dynamic_cast(those_scatts[i]); + + // Build the dense matrix for that object + double_3dvec 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 + 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]; + } + + // Incorporate that's contribution to the multiplicity matrix data + double nuscatt = that->scattxs[gin] * that->energy[gin][gout]; + mult_numer[gin][gout] += scalars[i] * nuscatt; + if (that->mult[gin][gout] > 0.) { + mult_denom[gin][gout] += scalars[i] * nuscatt / that->mult[gin][gout]; + } else { + mult_denom[gin][gout] += scalars[i]; + } + } + } + } + + // Combine mult_numer and mult_denom into the combined multiplicity matrix + double_2dvec this_mult(groups, double_1dvec(groups, 1.)); + 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]; + } + } + } + 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. + int_1dvec in_gmin(groups); + int_1dvec in_gmax(groups); + double_3dvec sparse_scatter(groups); + double_2dvec sparse_mult(groups); + for (int gin = 0; gin < groups; gin++) { + // Find the minimum and maximum group boundaries + int gmin_; + for (gmin_ = 0; gmin_ < groups; gmin_++) { + bool non_zero = std::all_of(this_matrix[gin][gmin_].begin(), + this_matrix[gin][gmin_].end(), + [](double val){return val > 0.;}); + if (non_zero) break; + } + int gmax_; + for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) { + bool non_zero = std::all_of(this_matrix[gin][gmax_].begin(), + this_matrix[gin][gmax_].end(), + [](double val){return val > 0.;}); + if (non_zero) break; + } + + // treat the case of all values being 0 + if (gmin_ > gmax_) { + gmin_ = gin; + gmax_ = gin; + } + + // Store the group bounds + in_gmin[gin] = gmin_; + in_gmax[gin] = gmax_; + + // Store the data in the compressed format + sparse_scatter[gin].resize(gmax_ - gmin_ + 1); + 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]; + } + } + + // Got everything we need, store it. + init(in_gmin, in_gmax, sparse_mult, sparse_scatter); +} + + +bool ScattDataLegendre::equiv(const ScattDataLegendre& that) +{ + // ensure that the number of groups match + return (this->energy.size() == that.energy.size()); +} + + +double_3dvec 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(order_dim, + double_1dvec(order_dim, 0.))); + + for (int gin = 0; gin < groups; gin++) { + for (int i_gout = 0; i_gout < energy[0].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] * + dist[gin][i_gout][l]; + } + } + } + return matrix; +} + +//============================================================================== +// ScattDataHistogram methods +//============================================================================== + +void ScattDataHistogram::init(int_1dvec in_gmin, int_1dvec in_gmax, + double_2dvec in_mult, double_3dvec coeffs) +{ + int groups = coeffs.size(); + int order = coeffs[0].size(); + + // make a copy of coeffs that we can use to both extract data and normalize + double_3dvec matrix = coeffs; + + // Get the scattering cross section value by summing the distribution + // over all the histogram bins in angle and outgoing energy groups + scattxs.resize(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.); + } + } + + // Build the energy transfer matrix from data in the variable matrix + double_2dvec in_energy; + in_energy.resize(groups); + for (int gin = 0; gin < groups; gin++) { + int num_groups = gmax[gin] - gmin[gin] + 1; + in_energy[gin].resize(num_groups); + for (int i_gout = 0; i_gout < num_groups; i_gout++) { + double norm = std::accumulate(matrix[gin][i_gout].begin(), + matrix[gin][i_gout].end(), 0.); + if (norm != 0.) { + for (auto& n : matrix[gin][i_gout]) n /= norm; + } + } + } + + // Initialize the base class attributes + ScattData::generic_init(order, in_gmin, in_gmax, in_energy, + in_mult); + + // Build the angular distributio mu values + mu = double_1dvec(order); + dmu = 2. / order; + mu[0] = -1.; + for (int imu = 1; imu < order; imu++) { + mu[imu] = -1. + (imu - 1) * dmu; + } + + // Calculate f(mu) and integrate it so we can avoid rejection sampling + fmu.resize(groups); + for (int gin = 0; gin < groups; gin++) { + int num_groups = gmax[gin] - gmin[gin] + 1; + fmu[gin].resize(num_groups); + for (int i_gout = 0; i_gout < num_groups; i_gout) { + fmu[gin][i_gout].resize(order); + // The variable matrix contains f(mu); so directly assign it + fmu[gin][i_gout] = matrix[gin][i_gout]; + + // Integrate the histogram + dist[gin][i_gout][0] = dmu * matrix[gin][i_gout][0]; + for (int imu = 1; imu < order; imu++) { + dist[gin][i_gout][imu] = dmu * matrix[gin][i_gout][imu] + + dist[gin][i_gout][imu - 1]; + } + + // Now re-normalize for integral to unity + double norm = dist[gin][i_gout][order - 1]; + if (norm > 0.) { + for (int imu = 0; imu < order; imu++) { + fmu[gin][i_gout][imu] /= norm; + dist[gin][i_gout][imu] /= norm; + } + } + } + } +} + + +double ScattDataHistogram::calc_f(int gin, int gout, double mu) +{ + // TODO: gout >= or gout >? + double f; + if ((gout < gmin[gin]) || (gout >= gmax[gin])) { + f = 0.; + } else { + // Find mu bin + int i_gout = gout - gmin[gin]; //TODO: + 1? + int imu; + if (mu == 1.) { + // use size -2 to have the index one before the end + imu = this->mu.size() - 2; + } else { + imu = std::floor((mu + 1.) / dmu + 1.); + } + + f = fmu[gin][i_gout][imu]; + } + return f; +} + + +void ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt) +{ + // Sample the outgoing energy using the base-class method + int i_gout; + sample_energy(gin, gout, i_gout); + + // Determine the outgoing cosine bin + double xi = prn(); + + int imu; + if (xi < dist[gin][i_gout][0]) { + imu = 1; + } else { + // TODO lower_bound? + 1? + imu = std::upper_bound(dist[gin][i_gout].begin(), + dist[gin][i_gout].end(), xi) - + dist[gin][i_gout].begin(); + } + + // Randomly select mu within the imu bin + mu = prn() * dmu + this->mu[imu]; + + if (mu < -1.) { + mu = -1.; + } else if (mu > 1.) { + mu = 1.; + } + + // Update the weight to reflect neutron multiplicity + wgt *= mult[gin][i_gout]; +} + + +bool ScattDataHistogram::equiv(const ScattDataHistogram& that) +{ + bool match = false; + if (this->energy.size() == that.energy.size() && + this->dmu == that.dmu && + std::equal(this->mu.begin(), this->mu.end(), that.mu.begin())) { + match = true; + } + return match; +} + + +double_3dvec 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(order_dim, + double_1dvec(order_dim, 0.))); + + for (int gin = 0; gin < groups; gin++) { + for (int i_gout = 0; i_gout < energy[0].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] * + fmu[gin][i_gout][l]; + } + } + } + return matrix; +} + +//============================================================================== +// ScattDataTabular methods +//============================================================================== + +void ScattDataTabular::init(int_1dvec in_gmin, int_1dvec in_gmax, + double_2dvec in_mult, double_3dvec coeffs) +{ + int groups = coeffs.size(); + int order = coeffs[0].size(); + + // make a copy of coeffs that we can use to both extract data and normalize + double_3dvec matrix = coeffs; + + // Build the angular distribution mu values + mu = double_1dvec(order); + dmu = 2. / (order - 1); + mu[0] = -1.; + for (int imu = 1; imu < order - 1; imu++) { + mu[imu] = -1. + (imu - 1) * 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); + 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] + + matrix[gin][i_gout][imu]); + } + } + } + + // Build the energy transfer matrix from data in the variable matrix + double_2dvec in_energy; + in_energy.resize(groups); + for (int gin = 0; gin < groups; gin++) { + int num_groups = gmax[gin] - gmin[gin] + 1; + in_energy[gin].resize(num_groups); + for (int i_gout = 0; i_gout < num_groups; i_gout++) { + double norm = 0.; + for (int imu = 1; imu < order; imu++) { + norm += 0.5 * dmu * (matrix[gin][i_gout][imu - 1] + + matrix[gin][i_gout][imu]); + } + if (norm != 0.) { + for (auto& n : matrix[gin][i_gout]) n /= norm; + } + } + } + + // Initialize the base class attributes + ScattData::generic_init(order, in_gmin, in_gmax, in_energy, + in_mult); + + // Calculate f(mu) and integrate it so we can avoid rejection sampling + fmu.resize(groups); + for (int gin = 0; gin < groups; gin++) { + int num_groups = gmax[gin] - gmin[gin] + 1; + fmu[gin].resize(num_groups); + for (int i_gout = 0; i_gout < num_groups; i_gout) { + fmu[gin][i_gout].resize(order); + // The variable matrix contains f(mu); so directly assign it + fmu[gin][i_gout] = matrix[gin][i_gout]; + + // Ensure positivity + for (auto& val : fmu[gin][i_gout]) { + if (val < 0.) val = 0.; + } + + // Now re-normalize for numerical integration issues and to take care of + // the above negative fix-up. Also accrue the CDF + double norm = 0.; + for (int imu = 1; imu < order; imu++) { + norm += 0.5 * dmu * (fmu[gin][i_gout][imu - 1] + + fmu[gin][i_gout][imu]); + // incorporate to the CDF + dist[gin][i_gout][imu] = norm; + } + + // now do the normalization + if (norm > 0.) { + for (int imu = 0; imu < order; imu++) { + fmu[gin][i_gout][imu] /= norm; + dist[gin][i_gout][imu] /= norm; + } + } + } + } +} + + +double ScattDataTabular::calc_f(int gin, int gout, double mu) +{ + // TODO: gout >= or gout >? + double f; + if ((gout < gmin[gin]) || (gout >= gmax[gin])) { + f = 0.; + } else { + // Find mu bin + int i_gout = gout - gmin[gin]; //TODO: + 1? + int imu; + if (mu == 1.) { + // use size -2 to have the index one before the end + imu = this->mu.size() - 2; + } else { + imu = std::floor((mu + 1.) / dmu + 1.); + } + + double r = (mu - this->mu[imu]) / (this->mu[imu + 1] - this->mu[imu]); + f = (1. - r) * fmu[gin][i_gout][imu] + r * fmu[gin][i_gout][imu + 1]; + } + return f; +} + + +void ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt) +{ + // Sample the outgoing energy using the base-class method + int i_gout; + sample_energy(gin, gout, i_gout); + + // Determine the outgoing cosine bin + int NP = this->mu.size(); + double xi = prn(); + + double c_k = dist[gin][i_gout][0]; + int k; + for (k = 0; k < NP - 2; k++) { + double c_k1 = dist[gin][i_gout][k + 1]; + if (xi < c_k1) break; + c_k = c_k1; + } + + // Check to make sure k is <= NP - 1 + k = std::min(k, NP - 1); + + // Find the pdf values we want + double p0 = fmu[gin][i_gout][k]; + double mu0 = this -> mu[k]; + double p1 = fmu[gin][i_gout][k + 1]; + double mu1 = this -> mu[k + 1]; + + if (p0 == p1) { + mu = mu0 + (xi - c_k) / p0; + } else { + double frac = (p1 - p0) / (mu1 - mu0); + mu = mu0 + (std::sqrt(std::max(0., p0 * p0 + 2. * frac * (xi - c_k))) + - p0) / frac; + } + + if (mu < -1.) { + mu = -1.; + } else if (mu > 1.) { + mu = 1.; + } + + // Update the weight to reflect neutron multiplicity + wgt *= mult[gin][i_gout]; +} + + +bool ScattDataTabular::equiv(const ScattDataTabular& that) +{ + bool match = false; + if (this->energy.size() == that.energy.size() && + this->dmu == that.dmu && + std::equal(this->mu.begin(), this->mu.end(), that.mu.begin())) { + match = true; + } + return match; +} + + +double_3dvec 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(order_dim, + double_1dvec(order_dim, 0.))); + + for (int gin = 0; gin < groups; gin++) { + for (int i_gout = 0; i_gout < energy[0].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] * + fmu[gin][i_gout][l]; + } + } + } + return matrix; +} + + +void convert_legendre_to_tabular(ScattDataLegendre& leg, + ScattDataTabular& tab, int n_mu) +{ + // Copy the obvious data + tab.energy = leg.energy; + tab.mult = leg.mult; + tab.gmin = leg.gmin; + tab.gmax = leg.gmax; + + // Build mu and dmu + tab.mu = double_1dvec(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 - 1) * 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(); + tab.fmu.resize(groups); + for (int gin = 0; gin < groups; gin++) { + int num_groups = tab.gmax[gin] - tab.gmin[gin] + 1; + tab.fmu[gin].resize(num_groups); + for (int i_gout = 0; i_gout < num_groups; i_gout) { + tab.fmu[gin][i_gout].resize(n_mu); + for (int imu = 0; imu < n_mu; imu++) { + tab.fmu[gin][i_gout][imu] = + evaluate_legendre_c(leg.dist[gin][i_gout].size() - 1, + leg.dist[gin][i_gout].data(), tab.mu[imu]); + } + + // Ensure positivity + for (auto& val : tab.fmu[gin][i_gout]) { + if (val < 0.) val = 0.; + } + + // Now re-normalize for numerical integration issues and to take care of + // the above negative fix-up. Also accrue the CDF + double norm = 0.; + for (int imu = 1; imu < n_mu; imu++) { + norm += 0.5 * tab.dmu * (tab.fmu[gin][i_gout][imu - 1] + + tab.fmu[gin][i_gout][imu]); + // incorporate to the CDF + tab.dist[gin][i_gout][imu] = norm; + } + + // now do the normalization + if (norm > 0.) { + for (int imu = 0; imu < n_mu; imu++) { + tab.fmu[gin][i_gout][imu] /= norm; + tab.dist[gin][i_gout][imu] /= norm; + } + } + } + } +} + +} // namespace openmc diff --git a/src/scattdata.h b/src/scattdata.h new file mode 100644 index 000000000..787a8f005 --- /dev/null +++ b/src/scattdata.h @@ -0,0 +1,106 @@ +//! \file scattdata.h +//! A collection of multi-group scattering data classes + +#ifndef SCATTDATA_H +#define SCATTDATA_H + +#include +#include +#include +#include + +#include "constants.h" +#include "math_functions.h" +#include "random_lcg.h" +#include "error.h" + +namespace openmc { + +//============================================================================== +// SCATTDATA contains all the data needed to describe the scattering energy and +// angular distribution data +//============================================================================== +// temporary declaations so we can name our friend functions +class ScattDataLegendre; +class ScattDataTabular; + +class ScattData { + protected: + 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 + public: + double_1dvec scattxs; // Isotropic Sigma_{s,g_{in}} + virtual double calc_f(int gin, int gout, double mu) = 0; + virtual void sample(int gin, int& gout, double& mu, double& wgt) = 0; + virtual void init(int_1dvec in_gmin, int_1dvec in_gmax, + double_2dvec in_mult, double_3dvec coeffs) = 0; + void sample_energy(int gin, int& gout, int& i_gout); + double get_xs(const char* xstype, int gin, int* gout, double* mu); + void generic_init(int order, int_1dvec in_gmin, int_1dvec in_gmax, + double_2dvec in_energy, double_2dvec in_mult); + virtual void combine(std::vector those_scatts, + double_1dvec& scalars) = 0; + virtual int get_order() = 0; + virtual double_3dvec get_matrix(int max_order) = 0; +}; + +class ScattDataLegendre: public ScattData { + protected: + // Maximal value for rejection sampling from a rectangle + double_2dvec max_val; + friend void convert_legendre_to_tabular(ScattDataLegendre& leg, + ScattDataTabular& tab, int n_mu); + public: + void init(int_1dvec in_gmin, int_1dvec in_gmax, double_2dvec in_mult, + double_3dvec coeffs); + void update_max_val(); + double calc_f(int gin, int gout, double mu); + void sample(int gin, int& gout, double& mu, double& wgt); + bool equiv(const ScattDataLegendre& that); + void combine(std::vector those_scatts, double_1dvec& scalars); + int get_order() {return dist[0][0].size() - 1;}; + double_3dvec get_matrix(int max_order); +}; + +class ScattDataHistogram: public ScattData { + protected: + double_1dvec mu; + double dmu; + double_3dvec fmu; + public: + void init(int_1dvec in_gmin, int_1dvec in_gmax, double_2dvec in_mult, + double_3dvec coeffs); + double calc_f(int gin, int gout, double mu); + void sample(int gin, int& gout, double& mu, double& wgt); + void combine(std::vector those_scatts, double_1dvec& scalars); + bool equiv(const ScattDataHistogram& that); + int get_order() {return dist[0][0].size();}; + double_3dvec get_matrix(int max_order); +}; + +class ScattDataTabular: public ScattData { + protected: + double_1dvec mu; + double dmu; + double_3dvec fmu; + friend void convert_legendre_to_tabular(ScattDataLegendre& leg, + ScattDataTabular& tab, int n_mu); + public: + void init(int_1dvec in_gmin, int_1dvec in_gmax, double_2dvec in_mult, + double_3dvec coeffs); + double calc_f(int gin, int gout, double mu); + void sample(int gin, int& gout, double& mu, double& wgt); + void combine(std::vector those_scatts, double_1dvec& scalars); + bool equiv(const ScattDataTabular& that); + int get_order() {return dist[0][0].size();}; + double_3dvec get_matrix(int max_order); +}; + +void convert_legendre_to_tabular(ScattDataLegendre& leg, + ScattDataTabular& tab, int n_mu); + +} // namespace openmc +#endif // SCATTDATA_H \ No newline at end of file diff --git a/src/string_functions.h b/src/string_functions.h new file mode 100644 index 000000000..94d4ce627 --- /dev/null +++ b/src/string_functions.h @@ -0,0 +1,54 @@ +//! \file string_functions.h +//! A collection of helper routines for C-strings and STL strings + +#ifndef STRING_FUNCTIONS_H +#define STRING_FUNCTIONS_H + +// for string functions +#include +#include +#include +#include + +namespace openmc { + +void strtrim(char* str) +{ + int start = 0; // number of leading spaces + char* buffer = str; + + while (*str && *str++ == ' ') ++start; + + while (*str++); // move to end of string + + // backup over trailing spaces + int end = str - buffer - 1; + while (end > 0 && buffer[end - 1] == ' ') --end; + buffer[end] = 0; // remove trailing spaces + + // exit if no leading spaces or string is now empty + if (end <= start || start == 0) return; + str = buffer + start; + + while ((*buffer++ = *str++)); // remove leading spaces: K&R +} + +std::string strtrim(std::string in_str) +{ + std::string str = in_str; + // perform the left trim + str.erase(str.begin(), std::find_if(str.begin(), str.end(), + std::not1(std::ptr_fun(std::isspace)))); + // perform the right trim + str.erase(std::find_if(str.rbegin(), str.rend(), + std::not1(std::ptr_fun(std::isspace))).base(), + str.end()); +} + +void to_lower(std::string& str) +{ + for (int i = 0; i < str.size(); i++) str[i] = std::tolower(str[i]); +} + +} // namespace openmc +#endif // STRING_FUNCTIONS_H \ No newline at end of file diff --git a/src/xsdata.cpp b/src/xsdata.cpp new file mode 100644 index 000000000..a027b18dd --- /dev/null +++ b/src/xsdata.cpp @@ -0,0 +1,782 @@ +#include "xsdata.h" + +namespace openmc { + +//============================================================================== +// XsData class methods +//============================================================================== + +XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable, + int scatter_format, int n_pol, int n_azi) +{ + // check to make sure scatter format is OK before we allocate + if (scatter_format != ANGLE_HISTOGRAM && scatter_format != ANGLE_TABULAR && + scatter_format != ANGLE_LEGENDRE) { + fatal_error("Invalid scatter_format!"); + } + // allocate all [temperature][phi][theta][in group] quantities + total = double_3dvec(n_pol, double_2dvec(n_azi, + double_1dvec(energy_groups, 0.))); + absorption = double_3dvec(n_pol, double_2dvec(n_azi, + double_1dvec(energy_groups, 0.))); + inverse_velocity = double_3dvec(n_pol, + double_2dvec(n_azi, double_1dvec(energy_groups, 0.))); + if (fissionable) { + fission = double_3dvec(n_pol, double_2dvec(n_azi, + double_1dvec(energy_groups, 0.))); + prompt_nu_fission = double_3dvec(n_pol, double_2dvec(n_azi, + double_1dvec(energy_groups, 0.))); + kappa_fission = double_3dvec(n_pol, double_2dvec(n_azi, + double_1dvec(energy_groups, 0.))); + } + + // allocate decay_rate; [temperature][phi][theta][delayed group] + decay_rate = double_3dvec(n_pol, double_2dvec(n_azi, + double_1dvec(num_delayed_groups, 0.))); + + if (fissionable) { + // allocate delayed_nu_fission; [temperature][phi][theta][in group][out group] + delayed_nu_fission = double_4dvec(n_pol, double_3dvec(n_azi, + double_2dvec(energy_groups, double_1dvec(energy_groups, 0.)))); + + // chi_prompt; [temperature][phi][theta][in group][delayed group] + chi_prompt = double_4dvec(n_pol, double_3dvec(n_azi, + double_2dvec(energy_groups, double_1dvec(num_delayed_groups, 0.)))); + + // chi_delayed; [temperature][phi][theta][in group][out group][delay group] + chi_delayed = double_5dvec(n_pol, double_4dvec(n_azi, + double_3dvec(energy_groups, double_2dvec(energy_groups, + double_1dvec(num_delayed_groups, 0.))))); + } + + scatter.resize(n_pol); + for (int p = 0; p < n_pol; p++) { + scatter[p].resize(n_azi); + for (int a = 0; a < n_azi; a++) { + if (scatter_format == ANGLE_HISTOGRAM) { + scatter[p][a] = new ScattDataHistogram; + } else if (scatter_format == ANGLE_TABULAR) { + scatter[p][a] = new ScattDataTabular; + } else if (scatter_format == ANGLE_LEGENDRE) { + scatter[p][a] = new ScattDataLegendre; + } + } + } +} + +XsData::~XsData() +{ + for (int p = 0; p < scatter.size(); p++) { + for (int a = 0; a < scatter[p].size(); a++) delete scatter[p][a]; + scatter[p].clear(); + } + scatter.clear(); +} + + +void XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, + int final_scatter_format, int order_data, int max_order, + int legendre_to_tabular_points) +{ + // Reconstruct the dimension information so it doesn't need to be passed + int n_pol = total.size(); + int n_azi = total[0].size(); + int energy_groups = total[0][0].size(); + int delayed_groups = decay_rate[0][0].size(); + + // Set the fissionable-specific data + if (fissionable) { + _fissionable_from_hdf5(xsdata_grp, n_pol, n_azi, energy_groups, + delayed_groups); + } + // Get the non-fission-specific data + read_nd_vector(xsdata_grp, "decay_rate", decay_rate); + read_nd_vector(xsdata_grp, "absorption", absorption, true); + read_nd_vector(xsdata_grp, "inverse-velocity", inverse_velocity); + + // Get scattering data + _scatter_from_hdf5(xsdata_grp, n_pol, n_azi, energy_groups, scatter_format, + final_scatter_format, order_data, max_order, + legendre_to_tabular_points); + + // Check absorption to ensure it is not 0 since it is often the + // denominator in tally methods + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + if (absorption[gin][p][a] == 0.) { + absorption[p][a][gin] = 1.e-10; + } + } + } + } + + // Get or calculate the total x/s + if (object_exists(xsdata_grp, "total")) { + read_nd_vector(xsdata_grp, "total", total); + } else { + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + total[p][a][gin] = absorption[p][a][gin] + + scatter[p][a]->scattxs[gin]; + } + } + } + } + + // Check total to ensure it is not 0 since it is often the denominator in + // tally methods + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + if (total[p][a][gin] == 0.) total[p][a][gin] = 1.e-10; + } + } + } +} + + +void XsData::_fissionable_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, + int energy_groups, int delayed_groups) +{ + double_4dvec temp_beta = + double_4dvec(n_pol, double_3dvec(n_azi, + double_2dvec(energy_groups, double_1dvec(delayed_groups, 0.)))); + + // Set/get beta + if (object_exists(xsdata_grp, "beta")) { + hid_t xsdata = open_dataset(xsdata_grp, "beta"); + int ndims = dataset_ndims(xsdata); + + if (ndims == 3) { + // Beta is input as [delayed group] + double_1dvec temp_arr = double_1dvec(n_pol * n_azi * delayed_groups); + read_nd_vector(xsdata_grp, "beta", temp_arr); + + // Broadcast to all incoming groups + int temp_idx = 0; + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int dg = 0; dg < delayed_groups; dg++) { + // Set the first group index and copy the rest + temp_beta[p][a][0][dg] = temp_arr[temp_idx++]; + for (int gin = 1; gin < energy_groups; gin++) { + temp_beta[p][a][gin] = temp_beta[p][a][0]; + } + } + } + } + } 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 1D or 2D array!"); + } + } + + // If chi is provided, set chi-prompt and chi-delayed + if (object_exists(xsdata_grp, "chi")) { + double_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi, + double_1dvec(energy_groups))); + read_nd_vector(xsdata_grp, "chi", temp_arr); + + int temp_idx = 0; + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + // First set the first group + for (int gout = 0; gout < energy_groups; gout++) { + chi_prompt[p][a][0][gout] = temp_arr[p][a][temp_idx++]; + } + + // Now normalize this data + double chi_sum = std::accumulate(chi_prompt[p][a][0].begin(), + chi_prompt[p][a][0].end(), + 0.); + if (chi_sum <= 0.) { + fatal_error("Encountered chi for a group that is <= 0!"); + } + for (int gout = 0; gout < energy_groups; gout++) { + chi_prompt[p][a][0][gout] /= chi_sum; + } + + // And extend to the remaining incoming groups + for (int gin = 1; gin < energy_groups; gin++) { + chi_prompt[p][a][gin] = chi_prompt[p][a][0]; + } + + // Finally set chi-delayed equal to chi-prompt + // Set chi-delayed to chi-prompt + 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[p][a][gin][gout][dg] = + chi_prompt[p][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"); + int ndims = dataset_ndims(xsdata); + + 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 (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + for (int dg = 0; dg < delayed_groups; dg++) { + delayed_nu_fission[p][a][gin][dg] = + temp_beta[p][a][gin][dg] * + prompt_nu_fission[p][a][gin]; + } + + // Correct the prompt-nu-fission using the delayed neutron fraction + if (delayed_groups > 0) { + double beta_sum = std::accumulate(temp_beta[p][a][gin].begin(), + temp_beta[p][a][gin].end(), 0.); + prompt_nu_fission[p][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 (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + double chi_sum = std::accumulate(chi_prompt[p][a][gin].begin(), + chi_prompt[p][a][gin].end(), 0.); + if (chi_sum >= 0.) { + for (int gout = 0; gout < energy_groups; gout++) { + chi_prompt[p][a][gin][gout] /= chi_sum; + } + } else { + fatal_error("Encountered chi for a group that is <= 0!"); + } + } + + // set chi-delayed to chi-prompt + 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[p][a][gin][gout][dg] = + chi_prompt[p][a][gin][gout]; + } + } + } + + // Set the vector nu-fission from the matrix nu-fission + for (int gin = 0; gin < energy_groups; gin++) { + double sum = std::accumulate(chi_prompt[p][a][gin].begin(), + chi_prompt[p][a][gin].end(), 0.); + prompt_nu_fission[p][a][gin] = sum; + } + + // 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[p][a][gin][dg] = + temp_beta[p][a][gin][dg] * + prompt_nu_fission[p][a][gin]; + } + + // Correct prompt-nu-fission using the delayed neutron fraction + if (delayed_groups > 0) { + double beta_sum = std::accumulate(temp_beta[p][a][gin].begin(), + temp_beta[p][a][gin].end(), 0.); + prompt_nu_fission[p][a][gin] *= (1. - beta_sum); + } + } + } + } + } else { + fatal_error("beta 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")) { + double_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi, + double_1dvec(energy_groups))); + read_nd_vector(xsdata_grp, "chi-prompt", temp_arr); + + for (int a = 0; a < n_azi; a++) { + for (int p = 0; p < n_pol; p++) { + for (int gin = 0; gin < energy_groups; gin++) { + for (int gout = 0; gout < energy_groups; gout++) { + chi_prompt[p][a][gin][gout] = temp_arr[p][a][gout]; + } + + // Normalize chi so its CDF goes to 1 + double chi_sum = std::accumulate(chi_prompt[p][a][gin].begin(), + chi_prompt[p][a][gin].end(), 0.); + if (chi_sum >= 0.) { + for (int gout = 0; gout < energy_groups; gout++) { + chi_prompt[p][a][gin][gout] /= chi_sum; + } + } else { + fatal_error("Encountered chi-prompt for a group that is <= 0.!"); + } + } + } + } + } + + // If chi-delayed is provided, set chi-delayed + if (object_exists(xsdata_grp, "chi-delayed")) { + hid_t xsdata = open_dataset(xsdata_grp, "chi-delayed"); + int ndims = dataset_ndims(xsdata); + close_dataset(xsdata); + + if (ndims == 3) { + // chi-delayed is a [in group] vector + double_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi, + double_1dvec(energy_groups))); + read_nd_vector(xsdata_grp, "chi-delayed", temp_arr); + + for (int a = 0; a < n_azi; a++) { + for (int p = 0; p < n_pol; p++) { + // normalize the chi CDF to 1 + double chi_sum = std::accumulate(temp_arr[p][a].begin(), + temp_arr[p][a].end(), 0.); + if (chi_sum <= 0.) { + fatal_error("Encountered chi-delayed for a group that is <= 0!"); + } + + // set chi-delayed + 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[p][a][gin][gout][dg] = + temp_arr[p][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 (int a = 0; a < n_azi; a++) { + for (int p = 0; p < n_pol; p++) { + for (int dg = 0; dg < delayed_groups; dg++) { + for (int gin = 0; gin < energy_groups; gin++) { + double chi_sum = 0.; + for (int gout = 0; gout < energy_groups; gout++) { + chi_sum += chi_delayed[p][a][gin][gout][dg]; + } + + if (chi_sum > 0.) { + for (int gout = 0; gout < energy_groups; gout++) { + chi_delayed[p][a][gin][gout][dg] /= chi_sum; + } + } else { + fatal_error("Encountered chi-delayed for a group that is <= 0!"); + } + } + } + } + } + } 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"); + int ndims = dataset_ndims(xsdata); + 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 + double_4dvec temp_arr = double_4dvec(n_pol, double_3dvec(n_azi, + double_2dvec(energy_groups, double_1dvec(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_azi; a++) { + for (int p = 0; p < n_pol; p++) { + for (int gin = 0; gin < energy_groups; gin++) { + double prompt_sum = std::accumulate(temp_arr[p][a][gin].begin(), + temp_arr[p][a][gin].end(), 0.); + prompt_nu_fission[p][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[p][a][gin] > 0.) { + for (int gout = 0; gout < energy_groups; gout++) { + chi_prompt[p][a][gin][gout] = + temp_arr[p][a][gin][gout] / + prompt_nu_fission[p][a][gin]; + } + } else { + fatal_error("Encountered chi-prompt for a group that is <= 0!"); + } + } + } + } + + } 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"); + int ndims = dataset_ndims(xsdata); + + if (ndims == 3) { + // delayed-nu-fission is a [in group] vector + if (temp_beta[0][0][0][0] == 0.) { + fatal_error("cannot set delayed-nu-fission with a 1D array if " + "beta is not provided"); + } + double_3dvec temp_arr = double_3dvec(n_pol, double_2dvec(n_azi, + double_1dvec(energy_groups))); + read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr); + + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; 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[p][a][gin][dg] = + temp_beta[p][a][gin][dg] * temp_arr[p][a][gin]; + } + } + } + } + + } else if (ndims == 4) { + // delayed nu fission is a [pol][azi][energy_group][delayed_group] matrix; + // matrix use this to set delayed-nu-fission separately for each + // delayed group + std::vector dims(ndims); + get_shape(xsdata, &dims[0]); + + if (dims[2] != delayed_groups) { + fatal_error("The delayed-nu-fission matrix was input with a 1st " + "dimension not equal to the number of delayed groups"); + } + if (dims[3] != energy_groups) { + fatal_error("The delayed-nu-fission matrix was input with a 2nd " + "dimension not equal to the number of energy groups"); + } + if (delayed_groups == energy_groups) { + warning("delayed-nu-fission was input as a dimension-4 matrix " + "with the same number of delayed groups and energy " + "groups. OpenMC assumes the dimensions in the matrix " + "are [delayed_groups][energy_groups]. Currently, " + "delayed-nu-fission cannot be set as a group-by-group " + "matrix"); + } + + read_nd_vector(xsdata_grp, "delayed-nu-fission", + delayed_nu_fission); + + } else if (ndims == 5) { + // This will contain delayed-nu-fision and chi-delayed data + double_5dvec temp_arr = double_5dvec(n_pol, double_4dvec(n_azi, + double_3dvec(energy_groups, double_2dvec(energy_groups, + double_1dvec(delayed_groups))))); + read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr); + + // Set the 4D delayed-nu-fission matrix and 5D chi-delayed matrix + // from the 5D delayed-nu-fission matrix + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int dg = 0; dg < delayed_groups; dg++) { + for (int gin = 0; gin < energy_groups; gin++) { + double gout_sum = 0.; + for (int gout = 0; gout < energy_groups; gout++) { + gout_sum += temp_arr[p][a][gin][gout][dg]; + chi_delayed[p][a][gin][gout][dg] = + temp_arr[p][a][gin][gout][dg]; + } + delayed_nu_fission[p][a][gin][dg] = gout_sum; + // Normalize chi-delayed + if (gout_sum > 0.) { + for (int gout = 0; gout < energy_groups; gout++) { + chi_delayed[p][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!"); + } + close_dataset(xsdata); + } + + // Get the fission and kappa_fission data xs + read_nd_vector(xsdata_grp, "fission", fission); + read_nd_vector(xsdata_grp, "kappa-fission", kappa_fission); +} + + +void 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) +{ + 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_3dvec gmin = int_3dvec(n_pol, int_2dvec(n_azi, + int_1dvec(energy_groups))); + read_nd_vector(scatt_grp, "g_min", gmin, true); + int_3dvec gmax = int_3dvec(n_pol, int_2dvec(n_azi, + int_1dvec(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 p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + gmin[p][a][gin] -= 1; + gmax[p][a][gin] -= 1; + } + } + } + + // Now use this info to find the length of a vector to hold the flattened + // data. + int length = 0; + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + length += order_data * (gmax[p][a][gin] - gmin[p][a][gin] + 1); + } + } + } + double_1dvec temp_arr = double_1dvec(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_5dvec input_scatt = + double_5dvec(n_pol, double_4dvec(n_azi, double_3dvec(energy_groups))); + + int temp_idx = 0; + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + input_scatt[p][a][gin].resize(gmax[p][a][gin] - gmin[p][a][gin] + 1); + for (int i_gout = 0; i_gout < input_scatt[p][a][gin].size(); i_gout++) { + input_scatt[p][a][gin][i_gout].resize(order_dim); + for (int l = 0; l < order_dim; l++) { + input_scatt[p][a][gin][i_gout][l] = temp_arr[temp_idx++]; + } + // Adjust index for the orders we didnt take + temp_idx += (order_data - order_dim); + } + } + } + } + temp_arr.clear(); + + // Get multiplication matrix + double_4dvec temp_mult = double_4dvec(n_pol, double_3dvec(n_azi, + double_2dvec(energy_groups))); + if (object_exists(scatt_grp, "multiplicity_matrix")) { + temp_arr.resize(length); + 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 p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + temp_mult[p][a][gin].resize(gmax[p][a][gin] - gmin[p][a][gin] + 1); + for (int i_gout = 0; i_gout < temp_mult[p][a][gin].size(); i_gout++) { + temp_mult[p][a][gin][i_gout] = temp_arr[temp_idx++]; + } + } + } + } + } else { + // Use a default: multiplicities are 1.0. + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + for (int gin = 0; gin < energy_groups; gin++) { + temp_mult[p][a][gin].resize(gmax[p][a][gin] - gmin[p][a][gin] + 1); + for (int i_gout = 0; i_gout < temp_mult[p][a][gin].size(); i_gout++) { + temp_mult[p][a][gin][i_gout] = 1.; + } + } + } + } + } + close_group(scatt_grp); + + // Finally, convert the Legendre data to tabular, if needed + if (scatter_format == ANGLE_LEGENDRE && + final_scatter_format == ANGLE_TABULAR) { + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + ScattDataLegendre legendre_scatt; + legendre_scatt.init(gmin[p][a], gmax[p][a], temp_mult[p][a], + input_scatt[p][a]); + + // Now create a tabular version of legendre_scatt + convert_legendre_to_tabular(legendre_scatt, + *static_cast(scatter[p][a]), + legendre_to_tabular_points); + + scatter_format = final_scatter_format; + } + } + } else { + // We are sticking with the current representation + // Initialize the ScattData object with this data + for (int p = 0; p < n_pol; p++) { + for (int a = 0; a < n_azi; a++) { + scatter[p][a]->init(gmin[p][a], gmax[p][a], temp_mult[p][a], + input_scatt[p][a]); + } + } + } +} + + +void XsData::combine(std::vector those_xs, double_1dvec& scalars) +{ + // Combine the non-scattering data + for (int i = 0; i < those_xs.size(); i++) { + XsData* that = those_xs[i]; + if (!equiv(*that)) fatal_error("Cannot combine the XsData objects!"); + double scalar = scalars[i]; + for (int p = 0; p < total.size(); p++) { + for (int a = 0; a < total[p].size(); a++) { + for (int gin = 0; gin < total[p][a].size(); gin++) { + total[p][a][gin] += scalar * that->total[p][a][gin]; + absorption[p][a][gin] += scalar * that->absorption[p][a][gin]; + inverse_velocity[p][a][gin] += + scalar * that->inverse_velocity[p][a][gin]; + + prompt_nu_fission[p][a][gin] += + scalar * that->prompt_nu_fission[p][a][gin]; + kappa_fission[p][a][gin] += + scalar * that->kappa_fission[p][a][gin]; + fission[p][a][gin] += + scalar * that->fission[p][a][gin]; + + for (int dg = 0; dg < delayed_nu_fission[p][a][gin].size(); dg++) { + delayed_nu_fission[p][a][gin][dg] += + scalar * that->delayed_nu_fission[p][a][gin][dg]; + } + + for (int gout = 0; gout < chi_prompt[p][a][gin].size(); gout++) { + chi_prompt[p][a][gin][gout] += + scalar * that->chi_prompt[p][a][gin][gout]; + + for (int dg = 0; dg < chi_delayed[p][a][gin][gout].size(); dg++) { + chi_delayed[p][a][gin][gout][dg] += + scalar * that->chi_delayed[p][a][gin][gout][dg]; + } + } + } + + for (int dg = 0; dg < decay_rate[p][a].size(); dg++) { + decay_rate[p][a][dg] += scalar * that->decay_rate[p][a][dg]; + } + + // Normalize chi + for (int gin = 0; gin < chi_prompt[p][a].size(); gin++) { + double norm = std::accumulate(chi_prompt[p][a][gin].begin(), + chi_prompt[p][a][gin].end(), 0.); + if (norm > 0.) { + for (int gout = 0; gout < chi_prompt[p][a][gin].size(); gout++) { + chi_prompt[p][a][gin][gout] /= norm; + } + } + + for (int dg = 0; dg < chi_delayed[p][a][gin][0].size(); dg++) { + norm = 0.; + for (int gout = 0; gout < chi_delayed[p][a][gin].size(); gout++) { + norm += chi_delayed[p][a][gin][gout][dg]; + } + if (norm > 0.) { + for (int gout = 0; gout < chi_delayed[p][a][gin].size(); gout++) { + chi_delayed[p][a][gin][gout][dg] /= norm; + } + } + } + } + } + } + } + + // Allow the ScattData object to combine itself + for (int p = 0; p < total.size(); p++) { + for (int a = 0; a < total[p].size(); 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++) { + those_scatts[i] = those_xs[i]->scatter[p][a]; + } + + // Now combine these guys + scatter[p][a]->combine(those_scatts, scalars); + } + } +} + + +bool XsData::equiv(const XsData& that) +{ + bool match = false; + // check n_pol (total.size()), n_azi (total[0].size()), and + // groups (total[0][0].size()) + // This assumes correct initializatino of the remaining cross sections + if ((total.size() == that.total.size()) && + (total[0].size() == that.total[0].size()) && + (total[0][0].size() == that.total[0][0].size())) { + match = true; + } + return match; +} + +} //namespace openmc diff --git a/src/xsdata.h b/src/xsdata.h new file mode 100644 index 000000000..1c1f8241b --- /dev/null +++ b/src/xsdata.h @@ -0,0 +1,72 @@ +//! \file xsdata.h +//! A collection of classes for containing the Multi-Group Cross Section data + +#ifndef XSDATA_H +#define XSDATA_H + +#include +#include +#include +#include +#include +#include +#include + +#include "constants.h" +#include "hdf5_interface.h" +#include "math_functions.h" +#include "random_lcg.h" +#include "scattdata.h" + + +namespace openmc { + +//============================================================================== +// XSDATA contains the temperature-independent cross section data for an MGXS +//============================================================================== + +class XsData { + private: + void _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); + void _fissionable_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, + int energy_groups, int delayed_groups); + public: + // The following quantities have the following dimensions: + // [phi][theta][incoming group] + double_3dvec total; + double_3dvec absorption; + double_3dvec prompt_nu_fission; + double_3dvec kappa_fission; + double_3dvec fission; + double_3dvec inverse_velocity; + // decay_rate has the following dimensions: + // [phi][theta][delayed group] + double_3dvec decay_rate; + // delayed_nu_fission has the following dimensions: + // [phi][theta][incoming group][delayed group] + double_4dvec delayed_nu_fission; + // chi_prompt has the following dimensions: + // [phi][theta][incoming group][outgoing group] + double_4dvec chi_prompt; + // chi_delayed has the following dimensions: + // [phi][theta][incoming group][outgoing group][delayed group] + double_5dvec chi_delayed; + // scatter has the following dimensions: [phi][theta] + std::vector > scatter; + + XsData() = default; + XsData(int num_groups, int num_delayed_groups, bool fissionable, + int scatter_format, int n_pol, int n_azi); + ~XsData(); + void from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, + int final_scatter_format, int order_data, int max_order, + int legendre_to_tabular_points); + void combine(std::vector those_xs, double_1dvec& scalars); + bool equiv(const XsData& that); +}; + + +} //namespace openmc +#endif // XSDATA_H \ No newline at end of file