diff --git a/include/openmc/constants.h b/include/openmc/constants.h index ab31781e2..9c96258fb 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -11,16 +11,9 @@ namespace openmc { -// TODO: Replace with xtensor/other library? -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; +using double_2dvec = std::vector>; +using double_3dvec = std::vector>>; +using double_4dvec = std::vector>>>; // ============================================================================ // VERSIONING NUMBERS diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index b53703fc9..468c25980 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -14,6 +14,7 @@ #include "xtensor/xarray.hpp" #include "openmc/position.h" +#include "openmc/error.h" namespace openmc { @@ -46,39 +47,6 @@ 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, 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); - -void -read_nd_vector(hid_t obj_id, const char* name, - std::vector > > > >& result, - bool must_have = false); - std::vector attribute_shape(hid_t obj_id, const char* name); std::vector dataset_names(hid_t group_id); void ensure_exists(hid_t group_id, const char* name); @@ -236,7 +204,7 @@ read_attribute(hid_t obj_id, const char* name, std::vector& vec) } //============================================================================== -// Templates/overloads for read_dataset +// Templates/overloads for read_dataset and related methods //============================================================================== template @@ -294,6 +262,40 @@ void read_dataset(hid_t obj_id, const char* name, xt::xarray& arr, bool indep close_dataset(dset); } + +template +void read_dataset_as_shape(hid_t obj_id, const char* name, + xt::xtensor& arr, bool indep=false) +{ + hid_t dset = open_dataset(obj_id, name); + + // Allocate new array to read data into + std::size_t size = 1; + for (const auto x : arr.shape()) + size *= x; + T* buffer = new T[size]; + + // Read data from attribute + read_dataset(dset, nullptr, H5TypeMap::type_id, buffer, indep); + + // Adapt into xarray + arr = xt::adapt(buffer, size, xt::acquire_ownership(), arr.shape()); + + close_dataset(dset); +} + + +template +void read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, + bool must_have=false) +{ + if (object_exists(obj_id, name)) { + read_dataset_as_shape(obj_id, name, result, true); + } else if (must_have) { + fatal_error(std::string("Must provide " + std::string(name) + "!")); + } +} + //============================================================================== // Templates/overloads for write_attribute //============================================================================== diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index 40746c291..e80355bee 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 980238aeb..8ea81f252 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); + base_combine(size_t max_order, const std::vector& those_scatts, + const std::vector& 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. @@ -81,7 +85,7 @@ class ScattData { //! @param scalars Scalars to multiply the microscopic data by. virtual void combine(const std::vector& those_scatts, - const double_1dvec& scalars) = 0; + const std::vector& scalars) = 0; //! \brief Getter for the dimensionality of the scattering order. //! @@ -89,7 +93,7 @@ class ScattData { //! of points, and for Histogram this is the number of bins. //! //! @return The order. - virtual int + virtual size_t get_order() = 0; //! \brief Builds a dense scattering matrix from the constituent parts @@ -97,8 +101,8 @@ 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 - get_matrix(int max_order) = 0; + virtual xt::xtensor + get_matrix(size_t max_order) = 0; //! \brief Samples the outgoing energy from the ScattData info. //! @@ -142,12 +146,12 @@ 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 combine(const std::vector& those_scatts, - const double_1dvec& scalars); + const std::vector& scalars); //! \brief Find the maximal value of the angular distribution to use as a // bounding box with rejection sampling. @@ -160,11 +164,11 @@ class ScattDataLegendre: public ScattData { void sample(int gin, int& gout, double& mu, double& wgt); - int + size_t get_order() {return dist[0][0].size() - 1;}; - double_3dvec - get_matrix(int max_order); + xt::xtensor + get_matrix(size_t max_order); }; //============================================================================== @@ -176,19 +180,19 @@ 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 combine(const std::vector& those_scatts, - const double_1dvec& scalars); + const std::vector& scalars); double calc_f(int gin, int gout, double mu); @@ -196,11 +200,11 @@ class ScattDataHistogram: public ScattData { void sample(int gin, int& gout, double& mu, double& wgt); - int + size_t get_order() {return dist[0][0].size();}; - double_3dvec - get_matrix(int max_order); + xt::xtensor + get_matrix(size_t 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,12 +229,12 @@ 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 combine(const std::vector& those_scatts, - const double_1dvec& scalars); + const std::vector& scalars); double calc_f(int gin, int gout, double mu); @@ -238,10 +242,11 @@ class ScattDataTabular: public ScattData { void sample(int gin, int& gout, double& mu, double& wgt); - int + size_t get_order() {return dist[0][0].size();}; - double_3dvec get_matrix(int max_order); + xt::xtensor + get_matrix(size_t max_order); }; //============================================================================== diff --git a/include/openmc/xsdata.h b/include/openmc/xsdata.h index 156708c78..fa75b10aa 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" @@ -22,41 +24,77 @@ class XsData { private: //! \brief Reads scattering data from the HDF5 file void - scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups, + scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, int scatter_format, int final_scatter_format, int order_data, int max_order, int legendre_to_tabular_points); //! \brief Reads fission data from the HDF5 file void - fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, int energy_groups, - int delayed_groups, bool is_isotropic); + fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, + size_t delayed_groups, bool is_isotropic); + + //! \brief Reads fission data formatted as chi and nu-fission vectors from + // the HDF5 file when beta is provided. + void + fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups, bool is_isotropic); + + //! \brief Reads fission data formatted as chi and nu-fission vectors from + // the HDF5 file when beta is not provided. + void + fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups); + + //! \brief Reads fission data formatted as chi and nu-fission vectors from + // the HDF5 file when no delayed data is provided. + void + fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups); + + //! \brief Reads fission data formatted as a nu-fission matrix from + // the HDF5 file when beta is provided. + void + fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups, bool is_isotropic); + + //! \brief Reads fission data formatted as a nu-fission matrix from + // the HDF5 file when beta is not provided. + void + fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups); + + //! \brief Reads fission data formatted as a nu-fission matrix from + // the HDF5 file when no delayed data is provided. + void + fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups); public: // The following quantities have the following dimensions: // [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; + std::vector> scatter; XsData() = default; @@ -68,7 +106,7 @@ class XsData { //! @param scatter_format The scattering representation of the file. //! @param n_pol Number of polar angles. //! @param n_azi Number of azimuthal angles. - XsData(int num_groups, int num_delayed_groups, bool fissionable, + XsData(size_t num_groups, size_t num_delayed_groups, bool fissionable, int scatter_format, int n_pol, int n_azi); //! \brief Loads the XsData object from the HDF5 file @@ -101,7 +139,7 @@ class XsData { //! @param micros Microscopic objects to combine. //! @param scalars Scalars to multiply the microscopic data by. void - combine(const std::vector& those_xs, const double_1dvec& scalars); + combine(const std::vector& those_xs, const std::vector& scalars); //! \brief Checks to see if this and that are able to be combined //! diff --git a/openmc/examples.py b/openmc/examples.py index d48d26839..a5138377e 100644 --- a/openmc/examples.py +++ b/openmc/examples.py @@ -1,3 +1,5 @@ +from numbers import Integral + import numpy as np import openmc @@ -538,20 +540,20 @@ def pwr_assembly(): return model -def slab_mg(reps=None, as_macro=True): - """Create a one-group, 1D slab model. +def slab_mg(num_regions=1, mat_names=None, mgxslib_name='2g.h5'): + """Create a 1D slab model. Parameters ---------- - reps : list, optional - List of angular representations. Each item corresponds to materials and - dictates the angular representation of the multi-group cross - sections---isotropic ('iso') or angle-dependent ('ang'), and if Legendre - scattering or tabular scattering ('mu') is used. Thus, items can be - 'ang', 'ang_mu', 'iso', or 'iso_mu'. + num_regions : int, optional + Number of regions in the problem, each with a unique MGXS dataset. + Defaults to 1. - as_macro : bool, optional - Whether :class:`openmc.Macroscopic` is used + mat_names : Iterable of str, optional + List of the material names to use; defaults to ['mat_1', 'mat_2',...]. + + mgxslib_name : str, optional + MGXS Library file to use; defaults to '2g.h5'. Returns ------- @@ -559,71 +561,82 @@ def slab_mg(reps=None, as_macro=True): One-group, 1D slab model """ + + openmc.check_type('num_regions', num_regions, Integral) + openmc.check_greater_than('num_regions', num_regions, 0) + if mat_names is not None: + openmc.check_length('mat_names', mat_names, num_regions) + openmc.check_iterable_type('mat_names', mat_names, str) + else: + mat_names = [] + for i in range(num_regions): + mat_names.append('mat_' + str(i + 1)) + + # # Make Materials + materials_file = openmc.Materials() + macros = [] + mats = [] + for i in range(len(mat_names)): + macros.append(openmc.Macroscopic('mat_' + str(i + 1))) + mats.append(openmc.Material(name=mat_names[i])) + mats[-1].set_density('macro', 1.0) + mats[-1].add_macroscopic(macros[-1]) + + materials_file += mats + + materials_file.cross_sections = mgxslib_name + + # # Make Geometry + rad_outer = 929.45 + # Set a cell boundary to exist for every material above (exclude the 0) + rads = np.linspace(0., rad_outer, len(mats) + 1, endpoint=True)[1:] + + # Instantiate Universe + root = openmc.Universe(universe_id=0, name='root universe') + cells = [] + + surfs = [] + surfs.append(openmc.XPlane(x0=0., boundary_type='reflective')) + for r, rad in enumerate(rads): + if r == len(rads) - 1: + surfs.append(openmc.XPlane(x0=rad, boundary_type='vacuum')) + else: + surfs.append(openmc.XPlane(x0=rad)) + + # Instantiate Cells + cells = [] + for c in range(len(surfs) - 1): + cells.append(openmc.Cell()) + cells[-1].region = (+surfs[c] & -surfs[c + 1]) + cells[-1].fill = mats[c] + + # Register Cells with Universe + root.add_cells(cells) + + # Instantiate a Geometry, register the root Universe, and export to XML + geometry_file = openmc.Geometry(root) + + # # Make Settings + # Instantiate a Settings object, set all runtime parameters + settings_file = openmc.Settings() + settings_file.energy_mode = "multi-group" + settings_file.tabular_legendre = {'enable': False} + settings_file.batches = 10 + settings_file.inactive = 5 + settings_file.particles = 1000 + + # Build source distribution + INF = 1000. + bounds = [0., -INF, -INF, rads[0], INF, INF] + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:]) + settings_file.source = openmc.source.Source(space=uniform_dist) + + settings_file.output = {'summary': False} + model = openmc.model.Model() - - # Define materials needed for 1D/1G slab problem - mat_names = ['uo2', 'clad', 'lwtr'] - mgxs_reps = ['ang', 'ang_mu', 'iso', 'iso_mu'] - - if reps is None: - reps = mgxs_reps - - xs = [] - i = 0 - for mat in mat_names: - for rep in reps: - i += 1 - name = mat + '_' + rep - xs.append(name) - if as_macro: - m = openmc.Material(name=str(i)) - m.set_density('macro', 1.) - m.add_macroscopic(name) - else: - m = openmc.Material(name=str(i)) - m.set_density('atom/b-cm', 1.) - m.add_nuclide(name, 1.0, 'ao') - model.materials.append(m) - - # Define the materials file - model.xs_data = xs - model.materials.cross_sections = "../../1d_mgxs.h5" - - # Define surfaces. - # Assembly/Problem Boundary - left = openmc.XPlane(x0=0.0, boundary_type='reflective') - right = openmc.XPlane(x0=10.0, boundary_type='reflective') - bottom = openmc.YPlane(y0=0.0, boundary_type='reflective') - top = openmc.YPlane(y0=10.0, boundary_type='reflective') - - # for each material add a plane - planes = [openmc.ZPlane(z0=0.0, boundary_type='reflective')] - dz = round(5. / float(len(model.materials)), 4) - for i in range(len(model.materials) - 1): - planes.append(openmc.ZPlane(z0=dz * float(i + 1))) - planes.append(openmc.ZPlane(z0=5.0, boundary_type='reflective')) - - # Define cells for each material - model.geometry.root_universe = openmc.Universe(name='root universe') - xy = +left & -right & +bottom & -top - for i, mat in enumerate(model.materials): - c = openmc.Cell(fill=mat, region=xy & +planes[i] & -planes[i + 1]) - model.geometry.root_universe.add_cell(c) - - model.settings.batches = 10 - model.settings.inactive = 5 - model.settings.particles = 100 - model.settings.source = openmc.Source(space=openmc.stats.Box( - [0.0, 0.0, 0.0], [10.0, 10.0, 5.])) - model.settings.energy_mode = "multi-group" - - plot = openmc.Plot() - plot.filename = 'mat' - plot.origin = (5.0, 5.0, 2.5) - plot.width = (2.5, 2.5) - plot.basis = 'xz' - plot.pixels = (3000, 3000) - plot.color_by = 'material' - model.plots.append(plot) + model.geometry = geometry_file + model.materials = materials_file + model.settings = settings_file + model.xs_data = macros return model diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index 5056e5f78..4badc0441 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -5,13 +5,15 @@ #include #include +#include "xtensor/xtensor.hpp" +#include "xtensor/xarray.hpp" + #include "hdf5.h" #include "hdf5_hl.h" #ifdef OPENMC_MPI #include "mpi.h" #include "openmc/message_passing.h" #endif -#include "openmc/error.h" namespace openmc { @@ -532,172 +534,6 @@ read_complex(hid_t obj_id, const char* name, std::complex* buffer, bool } -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.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, - std::vector >& result, bool must_have) -{ - if (object_exists(obj_id, name)) { - int dim1 = result.size(); - int dim2 = result[0].size(); - 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++]; - } - } - } 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(); - 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++]; - } - } - } 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(); - int dim3 = result[0][0].size(); - double temp_arr[dim1 * dim2 * dim3]; - 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++) { - 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)) { - int dim1 = result.size(); - int dim2 = result[0].size(); - int dim3 = result[0][0].size(); - int temp_arr[dim1 * dim2 * dim3]; - 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++) { - 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)) { - int dim1 = result.size(); - int dim2 = result[0].size(); - int dim3 = result[0][0].size(); - int dim4 = result[0][0][0].size(); - double temp_arr[dim1 * dim2 * dim3 * dim4]; - 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++) { - 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)) { - 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(); - double temp_arr[dim1 * dim2 * dim3 * dim4 * dim5]; - 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++) { - 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/input_xml.F90 b/src/input_xml.F90 index e756a6aa5..7361dd85b 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -181,7 +181,7 @@ contains ! Assign temperatures to cells that don't have temperatures already assigned call assign_temperatures() - ! Determine desired txemperatures for each nuclide and S(a,b) table + ! Determine desired temperatures for each nuclide and S(a,b) table call get_temperatures(nuc_temps, sab_temps) ! Check to make sure there are not too many nested coordinate levels in the diff --git a/src/mgxs.cpp b/src/mgxs.cpp index eb13a55f6..12c73fab4 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, *dg, gin); } 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, d, gin); } } } 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, *dg, gin, *gout); } else { - val = xs_t->chi_delayed[a][gin][*gout][0]; + val = xs_t->chi_delayed(a, 0, gin, *gout); } } 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, *dg, gin, g); } } 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, d, gin, g); } } } @@ -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, dg, gin); } // 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, dg, gin, gout); 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, dg, gin, gout); } } } @@ -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/mgxs_interface.cpp b/src/mgxs_interface.cpp index f601f3c71..56f3392ff 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -19,7 +19,7 @@ add_mgxs_c(hid_t file_id, const char* name, int energy_groups, int& method) { // Convert temps to a vector for the from_hdf5 function - double_1dvec temperature(temps, temps + n_temps); + std::vector temperature(temps, temps + n_temps); write_message("Loading " + std::string(name) + " data...", 6); @@ -60,10 +60,10 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[], { if (n_temps > 0) { // // Convert temps to a vector - double_1dvec temperature(temps, temps + n_temps); + std::vector temperature(temps, temps + n_temps); // Convert atom_densities to a vector - double_1dvec atom_densities_vec(atom_densities, + std::vector atom_densities_vec(atom_densities, atom_densities + n_nuclides); // Build array of pointers to nuclides_MG's Mgxs objects needed for this diff --git a/src/scattdata.cpp b/src/scattdata.cpp index 3e18c168f..8d7138624 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,11 +18,11 @@ 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(); + size_t groups = in_energy.size(); gmin = in_gmin; gmax = in_gmax; @@ -51,18 +53,17 @@ 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, +ScattData::base_combine(size_t max_order, + const std::vector& those_scatts, const std::vector& scalars, + xt::xtensor& in_gmin, xt::xtensor& in_gmax, double_2dvec& sparse_mult, double_3dvec& sparse_scatter) { - int groups = those_scatts[0] -> energy.size(); + size_t 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,8 @@ 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.)); - 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(); + xt::xtensor this_mult({groups, groups}, 1.); + this_mult = xt::nan_to_num(mult_numer / mult_denom); // 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 +118,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 +129,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 +153,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,21 +237,21 @@ 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(); + size_t groups = coeffs.size(); + size_t order = coeffs[0][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); + 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]; } @@ -301,7 +297,7 @@ ScattDataLegendre::init(const int_1dvec& in_gmin, const int_1dvec& in_gmax, void ScattDataLegendre::update_max_val() { - int groups = max_val.size(); + size_t groups = max_val.size(); // Step through the polynomial with fixed number of points to identify the // maximal value int Nmu = 1001; @@ -384,25 +380,25 @@ ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt) void ScattDataLegendre::combine(const std::vector& those_scatts, - const double_1dvec& scalars) + const std::vector& scalars) { // Find the max order in the data set and make sure we can combine the sets - int max_order = 0; + size_t max_order = 0; 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 (!that) { fatal_error("Cannot combine the ScattData objects!"); } - int that_order = that->get_order(); + size_t that_order = that->get_order(); if (that_order > max_order) max_order = that_order; } max_order++; // Add one since this is a Legendre - int groups = those_scatts[0] -> energy.size(); + size_t groups = those_scatts[0] -> energy.size(); - int_1dvec in_gmin(groups); - int_1dvec in_gmax(groups); + xt::xtensor in_gmin({groups}, 0); + xt::xtensor in_gmax({groups}, 0); double_3dvec sparse_scatter(groups); double_2dvec sparse_mult(groups); @@ -418,20 +414,19 @@ ScattDataLegendre::combine(const std::vector& those_scatts, //============================================================================== -double_3dvec -ScattDataLegendre::get_matrix(int max_order) +xt::xtensor +ScattDataLegendre::get_matrix(size_t 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.))); + size_t groups = energy.size(); + size_t order_dim = max_order + 1; + 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,20 +439,20 @@ 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(); + size_t groups = coeffs.size(); + size_t order = coeffs[0][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); + 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 +479,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); @@ -534,7 +525,7 @@ ScattDataHistogram::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; } @@ -560,7 +551,6 @@ ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt) if (xi < dist[gin][i_gout][0]) { imu = 0; } 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(); @@ -581,21 +571,20 @@ ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt) //============================================================================== -double_3dvec -ScattDataHistogram::get_matrix(int max_order) +xt::xtensor +ScattDataHistogram::get_matrix(size_t max_order) { // Get the sizes and initialize the data to 0 - int groups = energy.size(); + size_t 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.))); + size_t order_dim = get_order(); + 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]; } } @@ -607,10 +596,10 @@ ScattDataHistogram::get_matrix(int max_order) void ScattDataHistogram::combine(const std::vector& those_scatts, - const double_1dvec& scalars) + const std::vector& scalars) { // Find the max order in the data set and make sure we can combine the sets - int max_order = those_scatts[0]->get_order(); + size_t max_order = those_scatts[0]->get_order(); for (int i = 0; i < those_scatts.size(); i++) { // Lets also make sure these items are combineable ScattDataHistogram* that = dynamic_cast(those_scatts[i]); @@ -622,10 +611,10 @@ ScattDataHistogram::combine(const std::vector& those_scatts, } } - int groups = those_scatts[0] -> energy.size(); + size_t groups = those_scatts[0] -> energy.size(); - int_1dvec in_gmin(groups); - int_1dvec in_gmax(groups); + xt::xtensor in_gmin({groups}, 0); + xt::xtensor in_gmax({groups}, 0); double_3dvec sparse_scatter(groups); double_2dvec sparse_mult(groups); @@ -633,7 +622,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,29 +633,24 @@ 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(); + size_t groups = coeffs.size(); + size_t order = coeffs[0][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); + 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 +727,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 +748,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 +788,20 @@ ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt) //============================================================================== -double_3dvec -ScattDataTabular::get_matrix(int max_order) +xt::xtensor +ScattDataTabular::get_matrix(size_t max_order) { // Get the sizes and initialize the data to 0 - int groups = energy.size(); + size_t 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.))); + size_t order_dim = get_order(); + 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]; } } @@ -830,10 +813,10 @@ ScattDataTabular::get_matrix(int max_order) void ScattDataTabular::combine(const std::vector& those_scatts, - const double_1dvec& scalars) + const std::vector& scalars) { // Find the max order in the data set and make sure we can combine the sets - int max_order = those_scatts[0]->get_order(); + size_t max_order = those_scatts[0]->get_order(); for (int i = 0; i < those_scatts.size(); i++) { // Lets also make sure these items are combineable ScattDataTabular* that = dynamic_cast(those_scatts[i]); @@ -845,10 +828,10 @@ ScattDataTabular::combine(const std::vector& those_scatts, } } - int groups = those_scatts[0] -> energy.size(); + size_t groups = those_scatts[0] -> energy.size(); - int_1dvec in_gmin(groups); - int_1dvec in_gmax(groups); + xt::xtensor in_gmin({groups}, 0); + xt::xtensor in_gmax({groups}, 0); double_3dvec sparse_scatter(groups); double_2dvec sparse_mult(groups); @@ -856,7 +839,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,16 +868,11 @@ 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(); + size_t 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; diff --git a/src/xsdata.cpp b/src/xsdata.cpp index d9863694a..8abc1edfd 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -5,6 +5,11 @@ #include #include +#include "xtensor/xview.hpp" +#include "xtensor/xindex_view.hpp" +#include "xtensor/xmath.hpp" +#include "xtensor/xbuilder.hpp" + #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/math_functions.h" @@ -17,54 +22,53 @@ namespace openmc { // XsData class methods //============================================================================== -XsData::XsData(int energy_groups, int num_delayed_groups, bool fissionable, +XsData::XsData(size_t energy_groups, size_t 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 && scatter_format != ANGLE_LEGENDRE) { 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.)); + // allocate all [temperature][angle][in group] quantities + 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.)); + // allocate decay_rate; [temperature][angle][delayed group] + shape[1] = num_delayed_groups; + decay_rate = xt::zeros(shape); if (fissionable) { - // 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.))); + shape = {n_ang, num_delayed_groups, energy_groups}; + // allocate delayed_nu_fission; [temperature][angle][delay group][in group] + 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][angle][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.)))); + // chi_delayed; [temperature][angle][delay group][in group][out group] + shape = {n_ang, num_delayed_groups, energy_groups, energy_groups}; + chi_delayed = xt::zeros(shape); } for (int a = 0; a < n_ang; a++) { if (scatter_format == ANGLE_HISTOGRAM) { - // scatter[a] = std::make_unique(ScattDataHistogram); scatter.emplace_back(new ScattDataHistogram); } else if (scatter_format == ANGLE_TABULAR) { - // scatter[a] = std::make_unique(ScattDataTabular); scatter.emplace_back(new ScattDataTabular); } else if (scatter_format == ANGLE_LEGENDRE) { - // scatter[a] = std::make_unique(ScattDataLegendre); scatter.emplace_back(new ScattDataLegendre); } } @@ -78,13 +82,13 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, int legendre_to_tabular_points, bool is_isotropic, int n_pol, int n_azi) { // Reconstruct the dimension information so it doesn't need to be passed - int n_ang = n_pol * n_azi; - int energy_groups = total[0].size(); - int delayed_groups = decay_rate[0].size(); + size_t n_ang = n_pol * n_azi; + size_t energy_groups = total.shape()[1]; + size_t delayed_groups = decay_rate.shape()[1]; // Set the fissionable-specific data if (fissionable) { - fission_from_hdf5(xsdata_grp, n_pol, n_azi, energy_groups, delayed_groups, + fission_from_hdf5(xsdata_grp, n_ang, energy_groups, delayed_groups, is_isotropic); } // Get the non-fission-specific data @@ -93,445 +97,343 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, 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, + scatter_from_hdf5(xsdata_grp, n_ang, 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 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; - } - } + xt::filtration(absorption, xt::equal(absorption, 0.)) = 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 a = 0; a < n_ang; a++) { - for (int gin = 0; gin < energy_groups; gin++) { - total[a][gin] = absorption[a][gin] + scatter[a]->scattxs[gin]; + for (size_t a = 0; a < n_ang; a++) { + for (size_t gin = 0; gin < energy_groups; gin++) { + total(a, gin) = absorption(a, gin) + scatter[a]->scattxs[gin]; } } } // 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; - } - } + xt::filtration(total, xt::equal(total, 0.)) = 1.e-10; } //============================================================================== void -XsData::fission_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, - int energy_groups, int delayed_groups, bool is_isotropic) +XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups, bool is_isotropic) +{ + // Data is provided as nu-fission and chi with a beta for delayed info + + // Get chi + xt::xtensor temp_chi({n_ang, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "chi", temp_chi, true); + + // Normalize chi by summing over the outgoing groups for each incoming angle + temp_chi /= xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis()); + + // Now every incoming group in prompt_chi and delayed_chi is the normalized + // chi we just made + chi_prompt = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all()); + chi_delayed = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::newaxis(), + xt::all()); + + // Get nu-fission + xt::xtensor temp_nufiss({n_ang, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "nu-fission", temp_nufiss, true); + + // Get beta (strategy will depend upon the number of dimensions in beta) + hid_t beta_dset = open_dataset(xsdata_grp, "beta"); + int beta_ndims = dataset_ndims(beta_dset); + close_dataset(beta_dset); + int ndim_target = 1; + if (!is_isotropic) ndim_target += 2; + if (beta_ndims == ndim_target) { + xt::xtensor temp_beta({n_ang, delayed_groups}, 0.); + read_nd_vector(xsdata_grp, "beta", temp_beta, true); + + // Set prompt_nu_fission = (1. - beta_total)*nu_fission + prompt_nu_fission = temp_nufiss * (1. - xt::sum(temp_beta, {1})); + + // Set delayed_nu_fission as beta * nu_fission + delayed_nu_fission = + xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis()) * + xt::view(temp_nufiss, xt::all(), xt::newaxis(), xt::all()); + } else if (beta_ndims == ndim_target + 1) { + xt::xtensor temp_beta({n_ang, delayed_groups, energy_groups}, + 0.); + read_nd_vector(xsdata_grp, "beta", temp_beta, true); + + // Set prompt_nu_fission = (1. - beta_total)*nu_fission + prompt_nu_fission = temp_nufiss * (1. - xt::sum(temp_beta, {1})); + + // Set delayed_nu_fission as beta * nu_fission + delayed_nu_fission = temp_beta * + xt::view(temp_nufiss, xt::all(), xt::newaxis(), xt::all()); + } +} + +void +XsData::fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups) +{ + // Data is provided separately as prompt + delayed nu-fission and chi + + // Get chi-prompt + xt::xtensor temp_chi_p({n_ang, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "chi-prompt", temp_chi_p, true); + + // Normalize chi by summing over the outgoing groups for each incoming angle + temp_chi_p /= xt::view(xt::sum(temp_chi_p, {1}), xt::all(), xt::newaxis()); + + // Get chi-delayed + xt::xtensor temp_chi_d({n_ang, delayed_groups, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "chi-delayed", temp_chi_d, true); + + // Normalize chi by summing over the outgoing groups for each incoming angle + temp_chi_d /= xt::view(xt::sum(temp_chi_d, {2}), + xt::all(), xt::all(), xt::newaxis()); + + // Now assign the prompt and delayed chis by replicating for each incoming group + chi_prompt = xt::view(temp_chi_p, xt::all(), xt::newaxis(), xt::all()); + chi_delayed = xt::view(temp_chi_d, xt::all(), xt::all(), xt::newaxis(), + xt::all()); + + // Get prompt and delayed nu-fission directly + read_nd_vector(xsdata_grp, "prompt-nu-fission", prompt_nu_fission, + true); + read_nd_vector(xsdata_grp, "delayed-nu-fission", + delayed_nu_fission, true); +} + +void +XsData::fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups) +{ + // No beta is provided and there is no prompt/delay distinction. + // Therefore, the code only considers the data as prompt. + + // Get chi + xt::xtensor temp_chi({n_ang, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "chi", temp_chi, true); + + // Normalize chi by summing over the outgoing groups for each incoming angle + temp_chi /= xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis()); + + // Now every incoming group in self.chi is the normalized chi we just made + chi_prompt = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all()); + + // Get nu-fission directly + read_nd_vector(xsdata_grp, "nu-fission", prompt_nu_fission, true); +} + +//============================================================================== + +void +XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups, bool is_isotropic) +{ + // Data is provided as nu-fission and chi with a beta for delayed info + + // Get nu-fission matrix + xt::xtensor temp_matrix({n_ang, energy_groups, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true); + + // Get beta (strategy will depend upon the number of dimensions in beta) + hid_t beta_dset = open_dataset(xsdata_grp, "beta"); + int beta_ndims = dataset_ndims(beta_dset); + close_dataset(beta_dset); + int ndim_target = 1; + if (!is_isotropic) ndim_target += 2; + if (beta_ndims == ndim_target) { + xt::xtensor temp_beta({n_ang, delayed_groups}, 0.); + read_nd_vector(xsdata_grp, "beta", temp_beta, true); + + xt::xtensor temp_beta_sum({n_ang}, 0.); + temp_beta_sum = xt::sum(temp_beta, {1}); + + // prompt_nu_fission is the sum of this matrix over outgoing groups and + // multiplied by (1 - beta_sum) + prompt_nu_fission = xt::sum(temp_matrix, {2}) * (1. - temp_beta_sum); + + // Store chi-prompt + chi_prompt = xt::view(1.0 - temp_beta_sum, xt::all(), xt::newaxis(), + xt::newaxis()) * temp_matrix; + + // delayed_nu_fission is the sum of this matrix over outgoing groups and + // multiplied by beta + delayed_nu_fission = + xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis()) * + xt::view(xt::sum(temp_matrix, {2}), xt::all(), xt::newaxis(), xt::all()); + + // Store chi-delayed + chi_delayed = + xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis(), xt::newaxis()) * + xt::view(temp_matrix, xt::all(), xt::newaxis(), xt::all(), xt::all()); + + } else if (beta_ndims == ndim_target + 1) { + xt::xtensor temp_beta({n_ang, delayed_groups, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "beta", temp_beta, true); + + xt::xtensor temp_beta_sum({n_ang, energy_groups}, 0.); + temp_beta_sum = xt::sum(temp_beta, {1}); + + // prompt_nu_fission is the sum of this matrix over outgoing groups and + // multiplied by (1 - beta_sum) + prompt_nu_fission = xt::sum(temp_matrix, {2}) * (1. - temp_beta_sum); + + // Store chi-prompt + chi_prompt = xt::view(1.0 - temp_beta_sum, xt::all(), xt::all(), + xt::newaxis()) * temp_matrix; + + // delayed_nu_fission is the sum of this matrix over outgoing groups and + // multiplied by beta + delayed_nu_fission = temp_beta * + xt::view(xt::sum(temp_matrix, {2}), xt::all(), xt::newaxis(), xt::all()); + + // Store chi-delayed + chi_delayed = + xt::view(temp_beta, xt::all(), xt::all(), xt::all(), xt::newaxis()) * + xt::view(temp_matrix, xt::all(), xt::newaxis(), xt::all(), xt::all()); + } + + //Normalize both chis + chi_prompt /= xt::view(xt::sum(chi_prompt, {2}), + xt::all(), xt::all(), xt::newaxis()); + + chi_delayed /= xt::view(xt::sum(chi_delayed, {3}), + xt::all(), xt::all(), xt::all(), xt::newaxis()); +} + +void +XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups, size_t delayed_groups) +{ + // Data is provided separately as prompt + delayed nu-fission and chi + + // Get the prompt nu-fission matrix + xt::xtensor temp_matrix_p({n_ang, energy_groups, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_matrix_p, true); + + // prompt_nu_fission is the sum over outgoing groups + prompt_nu_fission = xt::sum(temp_matrix_p, {2}); + + // chi_prompt is this matrix but normalized over outgoing groups, which we + // have already stored in prompt_nu_fission + chi_prompt = temp_matrix_p / + xt::view(prompt_nu_fission, xt::all(), xt::all(), xt::newaxis()); + + // Get the delayed nu-fission matrix + xt::xtensor temp_matrix_d({n_ang, delayed_groups, energy_groups, + energy_groups}, 0.); + read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_matrix_d, true); + + // delayed_nu_fission is the sum over outgoing groups + delayed_nu_fission = xt::sum(temp_matrix_d, {3}); + + // chi_prompt is this matrix but normalized over outgoing groups, which we + // have already stored in prompt_nu_fission + chi_delayed = temp_matrix_d / + xt::view(delayed_nu_fission, xt::all(), xt::all(), xt::all(), xt::newaxis()); +} + +void +XsData::fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, + size_t energy_groups) +{ + // No beta is provided and there is no prompt/delay distinction. + // Therefore, the code only considers the data as prompt. + + // Get nu-fission matrix + xt::xtensor temp_matrix({n_ang, energy_groups, energy_groups}, 0.); + read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true); + + // prompt_nu_fission is the sum over outgoing groups + prompt_nu_fission = xt::sum(temp_matrix, {2}); + + // chi_prompt is this matrix but normalized over outgoing groups, which we + // have already stored in prompt_nu_fission + chi_prompt = temp_matrix / xt::view(prompt_nu_fission, xt::all(), xt::all(), + xt::newaxis()); +} + +//============================================================================== + +void +XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, + size_t delayed_groups, bool is_isotropic) { - int 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.))); - if (object_exists(xsdata_grp, "beta")) { - hid_t xsdata = open_dataset(xsdata_grp, "beta"); - int ndims = dataset_ndims(xsdata); - - // raise ndims to make the isotropic ndims the same as angular - if (is_isotropic) ndims += 2; - - if (ndims == 3) { - // Beta is input as [delayed group] - double_1dvec temp_arr(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 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++]; - for (int gin = 1; gin < energy_groups; gin++) { - temp_beta[a][gin] = temp_beta[a][0]; - } - } - } - } else if (ndims == 4) { - // Beta is input as [in group][delayed group] - read_nd_vector(xsdata_grp, "beta", temp_beta); + // Get the data; the strategy for doing so depends on if the data is provided + // as a nu-fission matrix or a set of chi and nu-fission vectors + if (object_exists(xsdata_grp, "chi") || + object_exists(xsdata_grp, "chi-prompt")) { + if (delayed_groups == 0) { + fission_vector_no_delayed_from_hdf5(xsdata_grp, n_ang, energy_groups); } else { - fatal_error("beta must be provided as a 3D or 4D array!"); - } - } - - // 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)); - 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]; - } - - // Now normalize this data - double chi_sum = std::accumulate(chi_prompt[a][0].begin(), - chi_prompt[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[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]; - } - - // 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[a][gin][gout][dg] = - chi_prompt[a][gin][gout]; - } - } + if (object_exists(xsdata_grp, "beta")) { + fission_vector_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, + delayed_groups, is_isotropic); + } else { + fission_vector_no_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, + delayed_groups); } } - } - - // 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); - // raise ndims to make the isotropic ndims the same as angular - if (is_isotropic) ndims += 2; - - if (ndims == 3) { - // nu-fission is a 3-d array - read_nd_vector(xsdata_grp, "nu-fission", prompt_nu_fission); - - // set delayed-nu-fission and correct prompt-nu-fission with beta - for (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]; - } - - // 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); - } - } - } - - } 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 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.); - // Set the vector nu-fission from the matrix nu-fission - 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; - } - } 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[a][gin][gout][dg] = - chi_prompt[a][gin][gout]; - } - } - } - - // 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]; - } - - // 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); - } - } - } + } else { + if (delayed_groups == 0) { + fission_matrix_no_delayed_from_hdf5(xsdata_grp, n_ang, energy_groups); } else { - fatal_error("nu-fission must be provided as a 3D or 4D array!"); - } - - close_dataset(xsdata); - } - - // If chi-prompt is provided, set chi-prompt - if (object_exists(xsdata_grp, "chi-prompt")) { - double_2dvec temp_arr(n_ang, double_1dvec(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]; - } - - // Normalize chi so its CDF goes to 1 - double chi_sum = std::accumulate(chi_prompt[a][gin].begin(), - chi_prompt[a][gin].end(), 0.); - if (chi_sum >= 0.) { - for (int gout = 0; gout < energy_groups; gout++) { - chi_prompt[a][gin][gout] /= chi_sum; - } - } else { - fatal_error("Encountered chi-prompt for a group that is <= 0.!"); - } + if (object_exists(xsdata_grp, "beta")) { + fission_matrix_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, + delayed_groups, is_isotropic); + } else { + fission_matrix_no_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, + delayed_groups); } } } - // If chi-delayed is provided, set chi-delayed - if (object_exists(xsdata_grp, "chi-delayed")) { - hid_t xsdata = open_dataset(xsdata_grp, "chi-delayed"); - int ndims = dataset_ndims(xsdata); - // raise ndims to make the isotropic ndims the same as angular - if (is_isotropic) ndims += 2; - close_dataset(xsdata); - - if (ndims == 3) { - // chi-delayed is a [in group] vector - double_2dvec temp_arr(n_ang, double_1dvec(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.); - 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[a][gin][gout][dg] = temp_arr[a][gout] / chi_sum; - } - } - } - } - } else if (ndims == 4) { - // chi_delayed is a matrix - read_nd_vector(xsdata_grp, "chi-delayed", chi_delayed); - - // Normalize the chi info so the CDF is 1. - for (int a = 0; a < n_ang; a++) { - 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[a][gin][gout][dg]; - } - - if (chi_sum > 0.) { - for (int gout = 0; gout < energy_groups; gout++) { - chi_delayed[a][gin][gout][dg] /= chi_sum; - } - } else { - fatal_error("Encountered chi-delayed for a group that is <= 0!"); - } - } - } - } - } else { - 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); - // raise ndims to make the isotropic ndims the same as angular - if (is_isotropic) ndims += 2; - close_dataset(xsdata); - - if (ndims == 3) { - // prompt-nu-fission is a [in group] vector - read_nd_vector(xsdata_grp, "prompt-nu-fission", - prompt_nu_fission); - } else if (ndims == 4) { - // prompt nu fission is a matrix, - // so set prompt_nu_fiss & chi_prompt - double_3dvec temp_arr(n_ang, 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_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; - } - - // 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.) { - for (int gout = 0; gout < energy_groups; gout++) { - chi_prompt[a][gin][gout] = - temp_arr[a][gin][gout] / prompt_nu_fission[a][gin]; - } - } else { - fatal_error("Encountered chi-prompt for a group that is <= 0!"); - } - } - } - - } 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); - close_dataset(xsdata); - // raise ndims to make the isotropic ndims the same as angular - if (is_isotropic) ndims += 2; - - if (ndims == 3) { - // delayed-nu-fission is an [in group] vector - if (temp_beta[0][0][0] == 0.) { - fatal_error("cannot set delayed-nu-fission with a 1D array if " - "beta is not provided"); - } - double_2dvec temp_arr(n_ang, double_1dvec(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]; - } - } - } - - } else if (ndims == 4) { - 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_4dvec temp_arr(n_ang, double_3dvec(energy_groups, - double_2dvec(energy_groups, double_1dvec(delayed_groups)))); - read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_arr); - - // Set the 3D delayed-nu-fission matrix and 4D chi-delayed matrix - // from the 4D delayed-nu-fission matrix - for (int a = 0; a < n_ang; 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[a][gin][gout][dg]; - chi_delayed[a][gin][gout][dg] = temp_arr[a][gin][gout][dg]; - } - delayed_nu_fission[a][gin][dg] = gout_sum; - // Normalize chi-delayed - if (gout_sum > 0.) { - for (int gout = 0; gout < energy_groups; gout++) { - chi_delayed[a][gin][gout][dg] /= gout_sum; - } - } else { - fatal_error("Encountered chi-delayed for a group that is <= 0!"); - } - } - } - } - - } else { - fatal_error("prompt-nu-fission must be provided as a 3D, 4D, or 5D " - "array!"); - } - } - // 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]); - } + if (delayed_groups == 0) { + nu_fission = prompt_nu_fission; + } else { + nu_fission = prompt_nu_fission + xt::sum(delayed_nu_fission, {1}); } } //============================================================================== 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) +XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, + int scatter_format, int final_scatter_format, int order_data, + int max_order, int legendre_to_tabular_points) { - int n_ang = n_pol * n_azi; if (!object_exists(xsdata_grp, "scatter_data")) { 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}, 0.); 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}, 0.); 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); - } - } - double_1dvec temp_arr(length); + size_t length = order_data * xt::sum(gmax - gmin + 1)(); + + double_4dvec input_scatt(n_ang, double_3dvec(energy_groups)); + xt::xtensor temp_arr({length}, 0.); read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true); // Compare the number of orders given with the max order of the problem; @@ -545,15 +447,13 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, // convert the flattened temp_arr to a jagged array for passing to // scatt data - double_4dvec input_scatt(n_ang, double_3dvec(energy_groups)); - - 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++) { + size_t temp_idx = 0; + for (size_t a = 0; a < n_ang; a++) { + for (size_t gin = 0; gin < energy_groups; gin++) { + input_scatt[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1); + for (size_t 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++) { + for (size_t 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 @@ -561,44 +461,46 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, } } } - 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); + temp_arr.resize({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); - for (int i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) { + size_t temp_idx = 0; + for (size_t a = 0; a < n_ang; a++) { + for (size_t gin = 0; gin < energy_groups; gin++) { + temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1); + for (size_t i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) { temp_mult[a][gin][i_gout] = temp_arr[temp_idx++]; } } } } 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); - for (int i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) { + for (size_t a = 0; a < n_ang; a++) { + for (size_t gin = 0; gin < energy_groups; gin++) { + temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1); + for (size_t 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 if (scatter_format == ANGLE_LEGENDRE && final_scatter_format == ANGLE_TABULAR) { - for (int a = 0; a < n_ang; a++) { + for (size_t 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, @@ -610,8 +512,10 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, } else { // 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]); + for (size_t a = 0; a < n_ang; a++) { + xt::xtensor in_gmin = xt::view(gmin, a, xt::all()); + xt::xtensor in_gmax = xt::view(gmax, a, xt::all()); + scatter[a]->init(in_gmin, in_gmax, temp_mult[a], input_scatt[a]); } } } @@ -620,80 +524,35 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, int n_pol, int n_azi, void XsData::combine(const std::vector& those_xs, - const double_1dvec& scalars) + const std::vector& scalars) { // Combine the non-scattering data - for (int i = 0; i < those_xs.size(); i++) { + for (size_t 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 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 (size_t 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++) { + for (size_t i = 0; i < those_xs.size(); i++) { those_scatts[i] = those_xs[i]->scatter[a].get(); } @@ -707,8 +566,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/1d_mgxs.h5 b/tests/1d_mgxs.h5 deleted file mode 100644 index 0f747345a..000000000 Binary files a/tests/1d_mgxs.h5 and /dev/null differ diff --git a/tests/regression_tests/mg_basic/inputs_true.dat b/tests/regression_tests/mg_basic/inputs_true.dat index a0efdbde0..4f2fd3f0b 100644 --- a/tests/regression_tests/mg_basic/inputs_true.dat +++ b/tests/regression_tests/mg_basic/inputs_true.dat @@ -1,97 +1,64 @@ - - - - - - - - - - - - + + + + + + - - - - - - - - - - - - - - - - + + + + + + - ../../1d_mgxs.h5 - + 2g.h5 + - + - + - + - + - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + eigenvalue - 100 + 1000 10 5 - 0.0 0.0 0.0 10.0 10.0 5.0 + 0.0 -1000.0 -1000.0 154.90833333333333 1000.0 1000.0 + + false + multi-group + + false + diff --git a/tests/regression_tests/mg_basic/results_true.dat b/tests/regression_tests/mg_basic/results_true.dat index ddb57d00b..dede97007 100644 --- a/tests/regression_tests/mg_basic/results_true.dat +++ b/tests/regression_tests/mg_basic/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.073147E+00 1.602384E-02 +1.005345E+00 1.109180E-02 diff --git a/tests/regression_tests/mg_basic/test.py b/tests/regression_tests/mg_basic/test.py index 3c22e6040..7b456bb9d 100644 --- a/tests/regression_tests/mg_basic/test.py +++ b/tests/regression_tests/mg_basic/test.py @@ -1,9 +1,97 @@ +import os + +import numpy as np + +import openmc from openmc.examples import slab_mg from tests.testing_harness import PyAPITestHarness +def create_library(): + # Instantiate the energy group data and file object + groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6]) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + + # Make the base, isotropic data + nu = [2.50, 2.50] + fiss = np.array([0.002817, 0.097]) + capture = [0.008708, 0.02518] + absorption = np.add(capture, fiss) + scatter = np.array( + [[[0.31980, 0.06694], [0.004555, -0.0003972]], + [[0.00000, 0.00000], [0.424100, 0.05439000]]]) + total = [0.33588, 0.54628] + chi = [1., 0.] + + mat_1 = openmc.XSdata('mat_1', groups) + mat_1.order = 1 + mat_1.set_nu_fission(np.multiply(nu, fiss)) + mat_1.set_absorption(absorption) + mat_1.set_scatter_matrix(scatter) + mat_1.set_total(total) + mat_1.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_1) + + # Make a version of mat-1 which has a tabular representation of the + # scattering vice Legendre with 33 points + mat_2 = mat_1.convert_scatter_format('tabular', 33) + mat_2.name = 'mat_2' + mg_cross_sections_file.add_xsdata(mat_2) + + # Make a version of mat-1 which has a histogram representation of the + # scattering vice Legendre with 33 bins + mat_3 = mat_1.convert_scatter_format('histogram', 33) + mat_3.name = 'mat_3' + mg_cross_sections_file.add_xsdata(mat_3) + + # Make a version which uses a fission matrix vice chi & nu-fission + mat_4 = openmc.XSdata('mat_4', groups) + mat_4.order = 1 + mat_4.set_nu_fission(np.outer(np.multiply(nu, fiss), chi)) + mat_4.set_absorption(absorption) + mat_4.set_scatter_matrix(scatter) + mat_4.set_total(total) + mg_cross_sections_file.add_xsdata(mat_4) + + # Make an angle-dependent version of mat_1 with 2 polar and 2 azim. angles + mat_5 = mat_1.convert_representation('angle', 2, 2) + mat_5.name = 'mat_5' + mg_cross_sections_file.add_xsdata(mat_5) + + # Make a copy of mat_1 for testing microscopic cross sections + mat_6 = openmc.XSdata('mat_6', groups) + mat_6.order = 1 + mat_6.set_nu_fission(np.multiply(nu, fiss)) + mat_6.set_absorption(absorption) + mat_6.set_scatter_matrix(scatter) + mat_6.set_total(total) + mat_6.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_6) + + # Write the file + mg_cross_sections_file.export_to_hdf5('2g.h5') + + +class MGXSTestHarness(PyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = '2g.h5' + if os.path.exists(f): + os.remove(f) + + def test_mg_basic(): - model = slab_mg() + create_library() + mat_names = ['base leg', 'base tab', 'base hist', 'base matrix', + 'base ang', 'micro'] + model = slab_mg(num_regions=6, mat_names=mat_names) + # Modify the last material to be a microscopic combination of nuclides + model.materials[-1] = openmc.Material(name='micro', material_id=6) + model.materials[-1].set_density("sum") + model.materials[-1].add_nuclide("mat_1", 0.5) + model.materials[-1].add_nuclide("mat_6", 0.5) + harness = PyAPITestHarness('statepoint.10.h5', model) harness.main() diff --git a/tests/regression_tests/mg_nuclide/__init__.py b/tests/regression_tests/mg_basic_delayed/__init__.py similarity index 100% rename from tests/regression_tests/mg_nuclide/__init__.py rename to tests/regression_tests/mg_basic_delayed/__init__.py diff --git a/tests/regression_tests/mg_basic_delayed/inputs_true.dat b/tests/regression_tests/mg_basic_delayed/inputs_true.dat new file mode 100644 index 000000000..9fb7afe7e --- /dev/null +++ b/tests/regression_tests/mg_basic_delayed/inputs_true.dat @@ -0,0 +1,63 @@ + + + + + + + + + + + + + + + + + + + 2g.h5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 10 + 5 + + + 0.0 -1000.0 -1000.0 154.90833333333333 1000.0 1000.0 + + + + false + + multi-group + + false + + diff --git a/tests/regression_tests/mg_basic_delayed/results_true.dat b/tests/regression_tests/mg_basic_delayed/results_true.dat new file mode 100644 index 000000000..85312b8e1 --- /dev/null +++ b/tests/regression_tests/mg_basic_delayed/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.003463E+00 2.173155E-02 diff --git a/tests/regression_tests/mg_basic_delayed/test.py b/tests/regression_tests/mg_basic_delayed/test.py new file mode 100644 index 000000000..f0474a567 --- /dev/null +++ b/tests/regression_tests/mg_basic_delayed/test.py @@ -0,0 +1,131 @@ +import os + +import numpy as np + +import openmc +from openmc.examples import slab_mg + +from tests.testing_harness import PyAPITestHarness + + +def create_library(): + # Instantiate the energy group data and file object + groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6]) + n_dg = 2 + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + mg_cross_sections_file.num_delayed_groups = n_dg + + beta = np.array([0.003, 0.003]) + one_m_beta = 1. - np.sum(beta) + nu = [2.50, 2.50] + fiss = np.array([0.002817, 0.097]) + capture = [0.008708, 0.02518] + absorption = np.add(capture, fiss) + scatter = np.array( + [[[0.31980, 0.06694], [0.004555, -0.0003972]], + [[0.00000, 0.00000], [0.424100, 0.05439000]]]) + total = [0.33588, 0.54628] + chi = [1., 0.] + + # Make the base data that uses chi & nu-fission vectors with a beta + mat_1 = openmc.XSdata('mat_1', groups) + mat_1.order = 1 + mat_1.num_delayed_groups = 2 + mat_1.set_beta(beta) + mat_1.set_nu_fission(np.multiply(nu, fiss)) + mat_1.set_absorption(absorption) + mat_1.set_scatter_matrix(scatter) + mat_1.set_total(total) + mat_1.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_1) + + # Make a version that uses prompt and delayed version of nufiss and chi + mat_2 = openmc.XSdata('mat_2', groups) + mat_2.order = 1 + mat_2.num_delayed_groups = 2 + mat_2.set_prompt_nu_fission(one_m_beta * np.multiply(nu, fiss)) + delay_nu_fiss = np.zeros((n_dg, groups.num_groups)) + for dg in range(n_dg): + for g in range(groups.num_groups): + delay_nu_fiss[dg, g] = beta[dg] * nu[g] * fiss[g] + mat_2.set_delayed_nu_fission(delay_nu_fiss) + mat_2.set_absorption(absorption) + mat_2.set_scatter_matrix(scatter) + mat_2.set_total(total) + mat_2.set_chi_prompt(chi) + mat_2.set_chi_delayed(np.stack([chi] * n_dg)) + mg_cross_sections_file.add_xsdata(mat_2) + + # Make a version that uses a nu-fission matrix with a beta + mat_3 = openmc.XSdata('mat_3', groups) + mat_3.order = 1 + mat_3.num_delayed_groups = 2 + mat_3.set_beta(beta) + mat_3.set_nu_fission(np.outer(np.multiply(nu, fiss), chi)) + mat_3.set_absorption(absorption) + mat_3.set_scatter_matrix(scatter) + mat_3.set_total(total) + mg_cross_sections_file.add_xsdata(mat_3) + + # Make a version that uses prompt and delayed version of the nufiss matrix + mat_4 = openmc.XSdata('mat_4', groups) + mat_4.order = 1 + mat_4.num_delayed_groups = 2 + mat_4.set_prompt_nu_fission(one_m_beta * + np.outer(np.multiply(nu, fiss), chi)) + delay_nu_fiss = np.zeros((n_dg, groups.num_groups, groups.num_groups)) + for dg in range(n_dg): + for g in range(groups.num_groups): + for go in range(groups.num_groups): + delay_nu_fiss[dg, g, go] = beta[dg] * nu[g] * fiss[g] * chi[go] + mat_4.set_delayed_nu_fission(delay_nu_fiss) + mat_4.set_absorption(absorption) + mat_4.set_scatter_matrix(scatter) + mat_4.set_total(total) + mg_cross_sections_file.add_xsdata(mat_4) + + # Make the base data that uses chi & nu-fiss vectors with a group-wise beta + mat_5 = openmc.XSdata('mat_5', groups) + mat_5.order = 1 + mat_5.num_delayed_groups = 2 + mat_5.set_beta(np.stack([beta] * groups.num_groups)) + mat_5.set_nu_fission(np.multiply(nu, fiss)) + mat_5.set_absorption(absorption) + mat_5.set_scatter_matrix(scatter) + mat_5.set_total(total) + mat_5.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_5) + + # Make a version that uses a nu-fission matrix with a group-wise beta + mat_6 = openmc.XSdata('mat_6', groups) + mat_6.order = 1 + mat_6.num_delayed_groups = 2 + mat_6.set_beta(np.stack([beta] * groups.num_groups)) + mat_6.set_nu_fission(np.outer(np.multiply(nu, fiss), chi)) + mat_6.set_absorption(absorption) + mat_6.set_scatter_matrix(scatter) + mat_6.set_total(total) + mg_cross_sections_file.add_xsdata(mat_6) + + # Write the file + mg_cross_sections_file.export_to_hdf5('2g.h5') + + +class MGXSTestHarness(PyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = '2g.h5' + if os.path.exists(f): + os.remove(f) + + +def test_mg_basic_delayed(): + create_library() + model = slab_mg(num_regions=6, mat_names=['vec beta', 'vec no beta', + 'matrix beta', 'matrix no beta', + 'vec group beta', + 'matrix group beta']) + + harness = PyAPITestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/mg_convert/test.py b/tests/regression_tests/mg_convert/test.py index 1ace10c80..4cc6b656a 100755 --- a/tests/regression_tests/mg_convert/test.py +++ b/tests/regression_tests/mg_convert/test.py @@ -17,7 +17,7 @@ def build_mgxs_library(convert): # Instantiate the energy group data groups = openmc.mgxs.EnergyGroups(group_edges=[1e-5, 0.625, 20.0e6]) - # Instantiate the 7-group (C5G7) cross section data + # Instantiate the 2-group (C5G7) cross section data uo2_xsdata = openmc.XSdata('UO2', groups) uo2_xsdata.order = 2 uo2_xsdata.set_total([2., 2.]) diff --git a/tests/regression_tests/mg_legendre/inputs_true.dat b/tests/regression_tests/mg_legendre/inputs_true.dat index 754808095..ad3b434e6 100644 --- a/tests/regression_tests/mg_legendre/inputs_true.dat +++ b/tests/regression_tests/mg_legendre/inputs_true.dat @@ -1,44 +1,31 @@ - - - + - - - - - - - + - ../../1d_mgxs.h5 - + 2g.h5 + - - - - - - - - - + eigenvalue - 100 + 1000 10 5 - 0.0 0.0 0.0 10.0 10.0 5.0 + 0.0 -1000.0 -1000.0 929.45 1000.0 1000.0 + + false + multi-group false diff --git a/tests/regression_tests/mg_legendre/results_true.dat b/tests/regression_tests/mg_legendre/results_true.dat index d0f8c319e..4a9d98237 100644 --- a/tests/regression_tests/mg_legendre/results_true.dat +++ b/tests/regression_tests/mg_legendre/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.110122E+00 2.549637E-02 +9.934975E-01 2.679669E-02 diff --git a/tests/regression_tests/mg_legendre/test.py b/tests/regression_tests/mg_legendre/test.py index 5a57f758e..b5a05c706 100644 --- a/tests/regression_tests/mg_legendre/test.py +++ b/tests/regression_tests/mg_legendre/test.py @@ -1,10 +1,54 @@ +import os + +import numpy as np + +import openmc from openmc.examples import slab_mg from tests.testing_harness import PyAPITestHarness +def create_library(): + # Instantiate the energy group data and file object + groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6]) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + + # Make the base, isotropic data + nu = [2.50, 2.50] + fiss = np.array([0.002817, 0.097]) + capture = [0.008708, 0.02518] + absorption = np.add(capture, fiss) + scatter = np.array( + [[[0.31980, 0.06694], [0.004555, -0.0003972]], + [[0.00000, 0.00000], [0.424100, 0.05439000]]]) + total = [0.33588, 0.54628] + chi = [1., 0.] + + mat_1 = openmc.XSdata('mat_1', groups) + mat_1.order = 1 + mat_1.set_nu_fission(np.multiply(nu, fiss)) + mat_1.set_absorption(absorption) + mat_1.set_scatter_matrix(scatter) + mat_1.set_total(total) + mat_1.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_1) + + # Write the file + mg_cross_sections_file.export_to_hdf5('2g.h5') + + +class MGXSTestHarness(PyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = '2g.h5' + if os.path.exists(f): + os.remove(f) + + def test_mg_legendre(): - model = slab_mg(reps=['iso']) + create_library() + model = slab_mg() model.settings.tabular_legendre = {'enable': False} harness = PyAPITestHarness('statepoint.10.h5', model) diff --git a/tests/regression_tests/mg_max_order/inputs_true.dat b/tests/regression_tests/mg_max_order/inputs_true.dat index 023d468d4..2ac83852c 100644 --- a/tests/regression_tests/mg_max_order/inputs_true.dat +++ b/tests/regression_tests/mg_max_order/inputs_true.dat @@ -1,44 +1,34 @@ - - - + - - - - - - - + - ../../1d_mgxs.h5 - + 2g.h5 + - - - - - - - - - + eigenvalue - 100 + 1000 10 5 - 0.0 0.0 0.0 10.0 10.0 5.0 + 0.0 -1000.0 -1000.0 929.45 1000.0 1000.0 + + false + multi-group 1 + + false + diff --git a/tests/regression_tests/mg_max_order/results_true.dat b/tests/regression_tests/mg_max_order/results_true.dat index adfcd44a8..4a9d98237 100644 --- a/tests/regression_tests/mg_max_order/results_true.dat +++ b/tests/regression_tests/mg_max_order/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.074551E+00 1.871525E-02 +9.934975E-01 2.679669E-02 diff --git a/tests/regression_tests/mg_max_order/test.py b/tests/regression_tests/mg_max_order/test.py index 20cc4f805..97d3f57d7 100644 --- a/tests/regression_tests/mg_max_order/test.py +++ b/tests/regression_tests/mg_max_order/test.py @@ -1,10 +1,55 @@ +import os + +import numpy as np + +import openmc from openmc.examples import slab_mg from tests.testing_harness import PyAPITestHarness +def create_library(): + # Instantiate the energy group data and file object + groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6]) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + + # Make the base, isotropic data + nu = [2.50, 2.50] + fiss = np.array([0.002817, 0.097]) + capture = [0.008708, 0.02518] + absorption = np.add(capture, fiss) + scatter = np.array( + [[[0.31980, 0.06694, 0.003], [0.004555, -0.0003972, 0.00002]], + [[0.00000, 0.00000, 0.000], [0.424100, 0.05439000, 0.0025]]]) + total = [0.33588, 0.54628] + chi = [1., 0.] + + mat_1 = openmc.XSdata('mat_1', groups) + mat_1.order = 2 + mat_1.set_nu_fission(np.multiply(nu, fiss)) + mat_1.set_absorption(absorption) + mat_1.set_scatter_matrix(scatter) + mat_1.set_total(total) + mat_1.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_1) + + # Write the file + mg_cross_sections_file.export_to_hdf5('2g.h5') + + +class MGXSTestHarness(PyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = '2g.h5' + if os.path.exists(f): + os.remove(f) + + def test_mg_max_order(): - model = slab_mg(reps=['iso']) + create_library() + model = slab_mg() model.settings.max_order = 1 + harness = PyAPITestHarness('statepoint.10.h5', model) harness.main() diff --git a/tests/regression_tests/mg_nuclide/inputs_true.dat b/tests/regression_tests/mg_nuclide/inputs_true.dat deleted file mode 100644 index e11b9e3f0..000000000 --- a/tests/regression_tests/mg_nuclide/inputs_true.dat +++ /dev/null @@ -1,97 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ../../1d_mgxs.h5 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - eigenvalue - 100 - 10 - 5 - - - 0.0 0.0 0.0 10.0 10.0 5.0 - - - multi-group - diff --git a/tests/regression_tests/mg_nuclide/results_true.dat b/tests/regression_tests/mg_nuclide/results_true.dat deleted file mode 100644 index ddb57d00b..000000000 --- a/tests/regression_tests/mg_nuclide/results_true.dat +++ /dev/null @@ -1,2 +0,0 @@ -k-combined: -1.073147E+00 1.602384E-02 diff --git a/tests/regression_tests/mg_nuclide/test.py b/tests/regression_tests/mg_nuclide/test.py deleted file mode 100644 index 44206ef28..000000000 --- a/tests/regression_tests/mg_nuclide/test.py +++ /dev/null @@ -1,9 +0,0 @@ -from openmc.examples import slab_mg - -from tests.testing_harness import PyAPITestHarness - - -def test_mg_nuclide(): - model = slab_mg(as_macro=False) - harness = PyAPITestHarness('statepoint.10.h5', model) - harness.main() diff --git a/tests/regression_tests/mg_survival_biasing/inputs_true.dat b/tests/regression_tests/mg_survival_biasing/inputs_true.dat index 4bc79d48e..5ece3ce9f 100644 --- a/tests/regression_tests/mg_survival_biasing/inputs_true.dat +++ b/tests/regression_tests/mg_survival_biasing/inputs_true.dat @@ -1,98 +1,34 @@ - - - - - - - - - - - - + - - - - - - - - - - - - - - - - + - ../../1d_mgxs.h5 - + 2g.h5 + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + eigenvalue - 100 + 1000 10 5 - 0.0 0.0 0.0 10.0 10.0 5.0 + 0.0 -1000.0 -1000.0 929.45 1000.0 1000.0 + + false + multi-group true + + false + diff --git a/tests/regression_tests/mg_survival_biasing/results_true.dat b/tests/regression_tests/mg_survival_biasing/results_true.dat index b20d63288..ebe98679f 100644 --- a/tests/regression_tests/mg_survival_biasing/results_true.dat +++ b/tests/regression_tests/mg_survival_biasing/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.080832E+00 1.336780E-02 +9.979905E-01 6.207495E-03 diff --git a/tests/regression_tests/mg_survival_biasing/test.py b/tests/regression_tests/mg_survival_biasing/test.py index 3c6c77a37..5d75611a9 100644 --- a/tests/regression_tests/mg_survival_biasing/test.py +++ b/tests/regression_tests/mg_survival_biasing/test.py @@ -1,10 +1,55 @@ +import os + +import numpy as np + +import openmc from openmc.examples import slab_mg from tests.testing_harness import PyAPITestHarness +def create_library(): + # Instantiate the energy group data and file object + groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6]) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + + # Make the base, isotropic data + nu = [2.50, 2.50] + fiss = np.array([0.002817, 0.097]) + capture = [0.008708, 0.02518] + absorption = np.add(capture, fiss) + scatter = np.array( + [[[0.31980, 0.06694], [0.004555, -0.0003972]], + [[0.00000, 0.00000], [0.424100, 0.05439000]]]) + total = [0.33588, 0.54628] + chi = [1., 0.] + + mat_1 = openmc.XSdata('mat_1', groups) + mat_1.order = 1 + mat_1.set_nu_fission(np.multiply(nu, fiss)) + mat_1.set_absorption(absorption) + mat_1.set_scatter_matrix(scatter) + mat_1.set_total(total) + mat_1.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_1) + + # Write the file + mg_cross_sections_file.export_to_hdf5('2g.h5') + + +class MGXSTestHarness(PyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = '2g.h5' + if os.path.exists(f): + os.remove(f) + + def test_mg_survival_biasing(): + create_library() model = slab_mg() model.settings.survival_biasing = True + harness = PyAPITestHarness('statepoint.10.h5', model) harness.main() diff --git a/tests/regression_tests/mg_tallies/inputs_true.dat b/tests/regression_tests/mg_tallies/inputs_true.dat index 7b5067014..a9b821c56 100644 --- a/tests/regression_tests/mg_tallies/inputs_true.dat +++ b/tests/regression_tests/mg_tallies/inputs_true.dat @@ -1,112 +1,48 @@ - - - - - - - - - - - - + - - - - - - - - - - - - - - - - + - ../../1d_mgxs.h5 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + 2g.h5 + + + eigenvalue - 100 + 1000 10 5 - 0.0 0.0 0.0 10.0 10.0 5.0 + 0.0 -1000.0 -1000.0 929.45 1000.0 1000.0 + + false + multi-group + + false + - 1 1 10 + 10 1 1 0.0 0.0 0.0 - 10 10 5 + 929.45 1000 1000 1 - 1 2 3 4 5 6 7 8 9 10 11 12 + 1 0.0 20000000.0 @@ -115,10 +51,10 @@ 0.0 20000000.0 - 1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0 + 0.0 0.625 20000000.0 - 1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0 + 0.0 0.625 20000000.0 5 @@ -170,60 +106,60 @@ 5 - uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu + mat_1 total absorption fission nu-fission analog 5 - uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu + mat_1 total absorption fission nu-fission tracklength 6 1 - uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu + mat_1 total absorption fission nu-fission scatter nu-scatter analog 6 1 - uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu + mat_1 total absorption fission nu-fission collision 6 1 - uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu + mat_1 total absorption fission nu-fission tracklength 6 1 2 - uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu + mat_1 scatter nu-scatter nu-fission 6 3 - uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu + mat_1 total absorption fission nu-fission scatter nu-scatter analog 6 3 - uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu + mat_1 total absorption fission nu-fission collision 6 3 - uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu + mat_1 total absorption fission nu-fission tracklength 6 3 4 - uo2_ang uo2_ang_mu uo2_iso uo2_iso_mu clad_ang clad_ang_mu clad_iso clad_iso_mu lwtr_ang lwtr_ang_mu lwtr_iso lwtr_iso_mu + mat_1 scatter nu-scatter nu-fission diff --git a/tests/regression_tests/mg_tallies/results_true.dat b/tests/regression_tests/mg_tallies/results_true.dat index 87470bdf4..78a5883fe 100644 --- a/tests/regression_tests/mg_tallies/results_true.dat +++ b/tests/regression_tests/mg_tallies/results_true.dat @@ -1 +1 @@ -9183f8b191f2e62334f992acd865d29e3f4e3f871a6df498e280fc4e2d91f2d2d20c732fbd75fa88e2e8c576f86e744f7655af6bb9da66e9b28b1009c8742899 \ No newline at end of file +41ea1f6b17c58a8141921af2f1d044eda93f3a9bca9463ee023af2e9865da613ace90fc8a25b42edde128ed827182ea9df0fe09d9b7887282d0ec092692cf717 \ No newline at end of file diff --git a/tests/regression_tests/mg_tallies/test.py b/tests/regression_tests/mg_tallies/test.py index 8952cc4ad..26d53c230 100644 --- a/tests/regression_tests/mg_tallies/test.py +++ b/tests/regression_tests/mg_tallies/test.py @@ -1,23 +1,66 @@ +import os + +import numpy as np + import openmc from openmc.examples import slab_mg from tests.testing_harness import HashedPyAPITestHarness +def create_library(): + # Instantiate the energy group data and file object + groups = openmc.mgxs.EnergyGroups(group_edges=[0.0, 0.625, 20.0e6]) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + + # Make the base, isotropic data + nu = [2.50, 2.50] + fiss = np.array([0.002817, 0.097]) + capture = [0.008708, 0.02518] + absorption = np.add(capture, fiss) + scatter = np.array( + [[[0.31980, 0.06694], [0.004555, -0.0003972]], + [[0.00000, 0.00000], [0.424100, 0.05439000]]]) + total = [0.33588, 0.54628] + chi = [1., 0.] + + mat_1 = openmc.XSdata('mat_1', groups) + mat_1.order = 1 + mat_1.set_nu_fission(np.multiply(nu, fiss)) + mat_1.set_absorption(absorption) + mat_1.set_scatter_matrix(scatter) + mat_1.set_total(total) + mat_1.set_chi(chi) + mg_cross_sections_file.add_xsdata(mat_1) + + # Write the file + mg_cross_sections_file.export_to_hdf5('2g.h5') + + +class MGXSTestHarness(HashedPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = '2g.h5' + if os.path.exists(f): + os.remove(f) + + def test_mg_tallies(): - model = slab_mg(as_macro=False) + create_library() + model = slab_mg() # Instantiate a tally mesh mesh = openmc.Mesh(mesh_id=1) mesh.type = 'regular' - mesh.dimension = [1, 1, 10] + mesh.dimension = [10, 1, 1] mesh.lower_left = [0.0, 0.0, 0.0] - mesh.upper_right = [10, 10, 5] + mesh.upper_right = [929.45, 1000, 1000] # Instantiate some tally filters energy_filter = openmc.EnergyFilter([0.0, 20.0e6]) energyout_filter = openmc.EnergyoutFilter([0.0, 20.0e6]) - energies = [1e-5, 0.0635, 10.0, 1.0e2, 1.0e3, 0.5e6, 1.0e6, 20.0e6] + energies = [0.0, 0.625, 20.0e6] matching_energy_filter = openmc.EnergyFilter(energies) matching_eout_filter = openmc.EnergyoutFilter(energies) mesh_filter = openmc.MeshFilter(mesh)