diff --git a/.gitmodules b/.gitmodules index f84d09bb1f9..2bc12389436 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,12 +1,6 @@ [submodule "vendor/pugixml"] path = vendor/pugixml url = https://github.com/zeux/pugixml.git -[submodule "vendor/xtensor"] - path = vendor/xtensor - url = https://github.com/xtensor-stack/xtensor.git -[submodule "vendor/xtl"] - path = vendor/xtl - url = https://github.com/xtensor-stack/xtl.git [submodule "vendor/fmt"] path = vendor/fmt url = https://github.com/fmtlib/fmt.git diff --git a/CMakeLists.txt b/CMakeLists.txt index e29c0d81c68..5d6dd2a0ea4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -266,23 +266,6 @@ else() endif() endif() -#=============================================================================== -# xtensor header-only library -#=============================================================================== - -if(OPENMC_FORCE_VENDORED_LIBS) - add_subdirectory(vendor/xtl) - set(xtl_DIR ${CMAKE_CURRENT_BINARY_DIR}/vendor/xtl) - add_subdirectory(vendor/xtensor) -else() - find_package_write_status(xtensor) - if (NOT xtensor_FOUND) - add_subdirectory(vendor/xtl) - set(xtl_DIR ${CMAKE_CURRENT_BINARY_DIR}/vendor/xtl) - add_subdirectory(vendor/xtensor) - endif() -endif() - #=============================================================================== # Catch2 library #=============================================================================== @@ -498,7 +481,7 @@ endif() # target_link_libraries treats any arguments starting with - but not -l as # linker flags. Thus, we can pass both linker flags and libraries together. target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES} - xtensor fmt::fmt ${CMAKE_DL_LIBS}) + fmt::fmt ${CMAKE_DL_LIBS}) if(TARGET pugixml::pugixml) target_link_libraries(libopenmc pugixml::pugixml) diff --git a/cmake/OpenMCConfig.cmake.in b/cmake/OpenMCConfig.cmake.in index 837a39c7833..cc837bba711 100644 --- a/cmake/OpenMCConfig.cmake.in +++ b/cmake/OpenMCConfig.cmake.in @@ -5,8 +5,6 @@ get_filename_component(_OPENMC_PREFIX "${OpenMC_CMAKE_DIR}/../../.." ABSOLUTE) find_package(fmt CONFIG REQUIRED HINTS ${_OPENMC_PREFIX}) find_package(pugixml CONFIG REQUIRED HINTS ${_OPENMC_PREFIX}) -find_package(xtl CONFIG REQUIRED HINTS ${_OPENMC_PREFIX}) -find_package(xtensor CONFIG REQUIRED HINTS ${_OPENMC_PREFIX}) if(@OPENMC_USE_DAGMC@) find_package(DAGMC REQUIRED HINTS @DAGMC_DIR@) endif() diff --git a/docs/source/quickinstall.rst b/docs/source/quickinstall.rst index 0f887940ea5..6cb774b5b07 100644 --- a/docs/source/quickinstall.rst +++ b/docs/source/quickinstall.rst @@ -119,7 +119,7 @@ packages should be installed, for example in Homebrew via: .. code-block:: sh - brew install llvm cmake xtensor hdf5 python libomp libpng + brew install llvm cmake hdf5 python libomp libpng The compiler provided by the above LLVM package should be used in place of the one provisioned by XCode, which does not support the multithreading library used diff --git a/include/openmc/bremsstrahlung.h b/include/openmc/bremsstrahlung.h index 2f7e41bf087..d3b5868317a 100644 --- a/include/openmc/bremsstrahlung.h +++ b/include/openmc/bremsstrahlung.h @@ -3,7 +3,7 @@ #include "openmc/particle.h" -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" namespace openmc { @@ -14,9 +14,9 @@ namespace openmc { class BremsstrahlungData { public: // Data - xt::xtensor pdf; //!< Bremsstrahlung energy PDF - xt::xtensor cdf; //!< Bremsstrahlung energy CDF - xt::xtensor yield; //!< Photon yield + tensor::Tensor pdf; //!< Bremsstrahlung energy PDF + tensor::Tensor cdf; //!< Bremsstrahlung energy CDF + tensor::Tensor yield; //!< Photon yield }; class Bremsstrahlung { @@ -32,9 +32,9 @@ class Bremsstrahlung { namespace data { -extern xt::xtensor +extern tensor::Tensor ttb_e_grid; //! energy T of incident electron in [eV] -extern xt::xtensor +extern tensor::Tensor ttb_k_grid; //! reduced energy W/T of emitted photon } // namespace data diff --git a/include/openmc/distribution_energy.h b/include/openmc/distribution_energy.h index 9b08ed039dc..efacb913196 100644 --- a/include/openmc/distribution_energy.h +++ b/include/openmc/distribution_energy.h @@ -5,7 +5,7 @@ #define OPENMC_DISTRIBUTION_ENERGY_H #include "hdf5.h" -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include "openmc/constants.h" #include "openmc/endf.h" @@ -86,9 +86,9 @@ class ContinuousTabular : public EnergyDistribution { struct CTTable { Interpolation interpolation; //!< Interpolation law int n_discrete; //!< Number of of discrete energies - xt::xtensor e_out; //!< Outgoing energies in [eV] - xt::xtensor p; //!< Probability density - xt::xtensor c; //!< Cumulative distribution + tensor::Tensor e_out; //!< Outgoing energies in [eV] + tensor::Tensor p; //!< Probability density + tensor::Tensor c; //!< Cumulative distribution }; int n_region_; //!< Number of inteprolation regions diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index b456fee21eb..43874767609 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -6,7 +6,7 @@ #include // for int64_t -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include #include "openmc/array.h" @@ -24,7 +24,7 @@ namespace simulation { extern double keff_generation; //!< Single-generation k on each processor extern array k_sum; //!< Used to reduce sum and sum_sq extern vector entropy; //!< Shannon entropy at each generation -extern xt::xtensor source_frac; //!< Source fraction for UFS +extern tensor::Tensor source_frac; //!< Source fraction for UFS } // namespace simulation diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 28b0d2b113d..93e4bfc820a 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -11,8 +11,7 @@ #include "hdf5.h" #include "hdf5_hl.h" -#include "xtensor/xadapt.hpp" -#include "xtensor/xarray.hpp" +#include "openmc/tensor.h" #include "openmc/array.h" #include "openmc/error.h" @@ -166,24 +165,19 @@ void read_attribute(hid_t obj_id, const char* name, vector& vec) read_attr(obj_id, name, H5TypeMap::type_id, vec.data()); } -// Generic array version +// Tensor version template -void read_attribute(hid_t obj_id, const char* name, xt::xarray& arr) +void read_attribute(hid_t obj_id, const char* name, tensor::Tensor& tensor) { - // Get shape of attribute array + // Get shape of attribute auto shape = attribute_shape(obj_id, name); - // Allocate new array to read data into - std::size_t size = 1; - for (const auto x : shape) - size *= x; - vector buffer(size); + // Resize tensor and read data directly + vector tshape(shape.begin(), shape.end()); + tensor.resize(tshape); // Read data from attribute - read_attr(obj_id, name, H5TypeMap::type_id, buffer.data()); - - // Adapt array into xarray - arr = xt::adapt(buffer, shape); + read_attr(obj_id, name, H5TypeMap::type_id, tensor.data()); } // overload for std::string @@ -290,61 +284,32 @@ void read_dataset( } template -void read_dataset(hid_t dset, xt::xarray& arr, bool indep = false) +void read_dataset(hid_t dset, tensor::Tensor& tensor, bool indep = false) { // Get shape of dataset vector shape = object_shape(dset); - // Allocate space in the array to read data into - std::size_t size = 1; - for (const auto x : shape) - size *= x; - arr.resize(shape); + // Resize tensor and read data directly + vector tshape(shape.begin(), shape.end()); + tensor.resize(tshape); - // Read data from attribute + // Read data from dataset read_dataset_lowlevel( - dset, nullptr, H5TypeMap::type_id, H5S_ALL, indep, arr.data()); + dset, nullptr, H5TypeMap::type_id, H5S_ALL, indep, tensor.data()); } template<> void read_dataset( - hid_t dset, xt::xarray>& arr, bool indep); + hid_t dset, tensor::Tensor>& tensor, bool indep); template void read_dataset( - hid_t obj_id, const char* name, xt::xarray& arr, bool indep = false) -{ - // Open dataset and read array - hid_t dset = open_dataset(obj_id, name); - read_dataset(dset, arr, indep); - close_dataset(dset); -} - -template -void read_dataset( - hid_t obj_id, const char* name, xt::xtensor& arr, bool indep = false) + hid_t obj_id, const char* name, tensor::Tensor& tensor, bool indep = false) { - // Open dataset and read array + // Open dataset and read tensor hid_t dset = open_dataset(obj_id, name); - - // Get shape of dataset - vector hsize_t_shape = object_shape(dset); + read_dataset(dset, tensor, indep); close_dataset(dset); - - // cast from hsize_t to size_t - vector shape(hsize_t_shape.size()); - for (int i = 0; i < shape.size(); i++) { - shape[i] = static_cast(hsize_t_shape[i]); - } - - // Allocate new xarray to read data into - xt::xarray xarr(shape); - - // Read data from the dataset - read_dataset(obj_id, name, xarr); - - // Copy into xtensor - arr = xarr; } // overload for Position @@ -358,31 +323,22 @@ inline void read_dataset( r.z = x[2]; } -template +template inline void read_dataset_as_shape( - hid_t obj_id, const char* name, xt::xtensor& arr, bool indep = false) + hid_t obj_id, const char* name, tensor::Tensor& tensor, 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; - vector buffer(size); - - // Read data from attribute + // Read data directly into pre-shaped tensor read_dataset_lowlevel( - dset, nullptr, H5TypeMap::type_id, H5S_ALL, indep, buffer.data()); - - // Adapt into xarray - arr = xt::adapt(buffer, arr.shape()); + dset, nullptr, H5TypeMap::type_id, H5S_ALL, indep, tensor.data()); close_dataset(dset); } -template -inline void read_nd_vector(hid_t obj_id, const char* name, - xt::xtensor& result, bool must_have = false) +template +inline void read_nd_tensor(hid_t obj_id, const char* name, + tensor::Tensor& result, bool must_have = false) { if (object_exists(obj_id, name)) { read_dataset_as_shape(obj_id, name, result, true); @@ -496,12 +452,16 @@ inline void write_dataset( false, buffer.data()); } -// Template for xarray, xtensor, etc. -template -inline void write_dataset( - hid_t obj_id, const char* name, const xt::xcontainer& arr) +// Template for Tensor and StaticTensor2D. A SFINAE guard is used here to +// prevent this template from matching vector/string types that have their own +// overloads above. A generic Container parameter avoids duplicating the body +// for both Tensor and StaticTensor2D. +template>::value>> +inline void write_dataset(hid_t obj_id, const char* name, const Container& arr) { - using T = typename D::value_type; + using T = typename std::decay_t::value_type; auto s = arr.shape(); vector dims {s.cbegin(), s.cend()}; write_dataset_lowlevel(obj_id, dims.size(), dims.data(), name, diff --git a/include/openmc/material.h b/include/openmc/material.h index c10f25551e8..3967c36878f 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -5,8 +5,8 @@ #include #include "openmc/span.h" +#include "openmc/tensor.h" #include "pugixml.hpp" -#include "xtensor/xtensor.hpp" #include #include "openmc/bremsstrahlung.h" @@ -189,7 +189,7 @@ class Material { vector nuclide_; //!< Indices in nuclides vector vector element_; //!< Indices in elements vector NCrystalMat ncrystal_mat_; //!< NCrystal material object - xt::xtensor atom_density_; //!< Nuclide atom density in [atom/b-cm] + tensor::Tensor atom_density_; //!< Nuclide atom density in [atom/b-cm] double density_; //!< Total atom density in [atom/b-cm] double density_gpcc_; //!< Total atom density in [g/cm^3] double charge_density_; //!< Total charge density in [e/b-cm] diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 5dc327fe359..0d8189caa1d 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -8,8 +8,8 @@ #include #include "hdf5.h" +#include "openmc/tensor.h" #include "pugixml.hpp" -#include "xtensor/xtensor.hpp" #include "openmc/bounding_box.h" #include "openmc/error.h" @@ -284,8 +284,8 @@ class Mesh { virtual Position upper_right() const = 0; // Data members - xt::xtensor lower_left_; //!< Lower-left coordinates of mesh - xt::xtensor upper_right_; //!< Upper-right coordinates of mesh + tensor::Tensor lower_left_; //!< Lower-left coordinates of mesh + tensor::Tensor upper_right_; //!< Upper-right coordinates of mesh int id_ {-1}; //!< Mesh ID std::string name_; //!< User-specified name int n_dimension_ {-1}; //!< Number of dimensions @@ -348,7 +348,7 @@ class StructuredMesh : public Mesh { //! \param[in] Pointer to bank sites //! \param[in] Number of bank sites //! \param[out] Whether any bank sites are outside the mesh - xt::xtensor count_sites( + tensor::Tensor count_sites( const SourceSite* bank, int64_t length, bool* outside) const; //! Get bin given mesh indices @@ -419,8 +419,8 @@ class StructuredMesh : public Mesh { //! Get a label for the mesh bin std::string bin_label(int bin) const override; - //! Get shape as xt::xtensor - xt::xtensor get_x_shape() const; + //! Get mesh dimensions as a tensor + tensor::Tensor get_shape_tensor() const; double volume(int bin) const override { @@ -515,7 +515,7 @@ class RegularMesh : public StructuredMesh { //! \param[in] bank Array of bank sites //! \param[out] Whether any bank sites are outside the mesh //! \return Array indicating number of sites in each mesh/energy bin - xt::xtensor count_sites( + tensor::Tensor count_sites( const SourceSite* bank, int64_t length, bool* outside) const; //! Return the volume for a given mesh index @@ -526,7 +526,7 @@ class RegularMesh : public StructuredMesh { // Data members double volume_frac_; //!< Volume fraction of each mesh element double element_volume_; //!< Volume of each mesh element - xt::xtensor width_; //!< Width of each mesh element + tensor::Tensor width_; //!< Width of each mesh element }; class RectilinearMesh : public StructuredMesh { diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index 43505f432d6..5d8ff58c239 100644 --- a/include/openmc/mgxs.h +++ b/include/openmc/mgxs.h @@ -6,7 +6,7 @@ #include -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include "openmc/constants.h" #include "openmc/hdf5_interface.h" @@ -22,7 +22,7 @@ namespace openmc { class Mgxs { private: - xt::xtensor kTs; // temperature in eV (k * T) + tensor::Tensor kTs; // temperature in eV (k * T) AngleDistributionType scatter_format; // flag for if this is legendre, histogram, or tabular int num_groups; // number of energy groups diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 58d83393963..7a8b2acadd9 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -96,7 +96,7 @@ class Nuclide { // Temperature dependent cross section data vector kTs_; //!< temperatures in eV (k*T) vector grid_; //!< Energy grid at each temperature - vector> xs_; //!< Cross sections at each temperature + vector> xs_; //!< Cross sections at each temperature // Multipole data unique_ptr multipole_; diff --git a/include/openmc/photon.h b/include/openmc/photon.h index f6f28a4df1d..93c6dba53f5 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -6,7 +6,7 @@ #include "openmc/particle.h" #include "openmc/vector.h" -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include #include @@ -62,14 +62,14 @@ class PhotonInteraction { int64_t index_; //!< Index in global elements vector // Microscopic cross sections - xt::xtensor energy_; - xt::xtensor coherent_; - xt::xtensor incoherent_; - xt::xtensor photoelectric_total_; - xt::xtensor pair_production_total_; - xt::xtensor pair_production_electron_; - xt::xtensor pair_production_nuclear_; - xt::xtensor heating_; + tensor::Tensor energy_; + tensor::Tensor coherent_; + tensor::Tensor incoherent_; + tensor::Tensor photoelectric_total_; + tensor::Tensor pair_production_total_; + tensor::Tensor pair_production_electron_; + tensor::Tensor pair_production_nuclear_; + tensor::Tensor heating_; // Form factors Tabulated1D incoherent_form_factor_; @@ -81,27 +81,27 @@ class PhotonInteraction { // stored separately to improve memory access pattern when calculating the // total cross section vector shells_; - xt::xtensor cross_sections_; + tensor::Tensor cross_sections_; // Compton profile data - xt::xtensor profile_pdf_; - xt::xtensor profile_cdf_; - xt::xtensor binding_energy_; - xt::xtensor electron_pdf_; + tensor::Tensor profile_pdf_; + tensor::Tensor profile_cdf_; + tensor::Tensor binding_energy_; + tensor::Tensor electron_pdf_; // Map subshells from Compton profile data obtained from Biggs et al, // "Hartree-Fock Compton profiles for the elements" to ENDF/B atomic // relaxation data - xt::xtensor subshell_map_; + tensor::Tensor subshell_map_; // Stopping power data double I_; // mean excitation energy - xt::xtensor n_electrons_; - xt::xtensor ionization_energy_; - xt::xtensor stopping_power_radiative_; + tensor::Tensor n_electrons_; + tensor::Tensor ionization_energy_; + tensor::Tensor stopping_power_radiative_; // Bremsstrahlung scaled DCS - xt::xtensor dcs_; + tensor::Tensor dcs_; // Whether atomic relaxation data is present bool has_atomic_relaxation_ {false}; @@ -137,7 +137,7 @@ void free_memory_photon(); namespace data { -extern xt::xtensor +extern tensor::Tensor compton_profile_pz; //! Compton profile momentum grid //! Photon interaction data for each element diff --git a/include/openmc/plot.h b/include/openmc/plot.h index 3813101c65c..cd6be43b724 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -6,8 +6,8 @@ #include #include +#include "openmc/tensor.h" #include "pugixml.hpp" -#include "xtensor/xarray.hpp" #include "hdf5.h" #include "openmc/cell.h" @@ -90,7 +90,7 @@ const RGBColor BLACK {0, 0, 0}; * visualized. */ -typedef xt::xtensor ImageData; +typedef tensor::Tensor ImageData; class PlottableInterface { public: PlottableInterface() = default; @@ -154,7 +154,7 @@ struct IdData { void set_overlap(size_t y, size_t x); // Members - xt::xtensor data_; //!< 2D array of cell & material ids + tensor::Tensor data_; //!< 2D array of cell & material ids }; struct PropertyData { @@ -166,7 +166,7 @@ struct PropertyData { void set_overlap(size_t y, size_t x); // Members - xt::xtensor data_; //!< 2D array of temperature & density data + tensor::Tensor data_; //!< 2D array of temperature & density data }; //=============================================================================== diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 693093683a9..c40982712eb 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -178,10 +178,10 @@ class FlatSourceDomain { // Volumes for each tally and bin/score combination. This intermediate data // structure is used when tallying quantities that must be normalized by // volume (i.e., flux). The vector is index by tally index, while the inner 2D - // xtensor is indexed by bin index and score index in a similar manner to the + // tensor is indexed by bin index and score index in a similar manner to the // results tensor in the Tally class, though without the third dimension, as // SUM and SUM_SQ do not need to be tracked. - vector> tally_volumes_; + vector> tally_volumes_; }; // class FlatSourceDomain diff --git a/include/openmc/scattdata.h b/include/openmc/scattdata.h index a75ef09d97b..bea881402ac 100644 --- a/include/openmc/scattdata.h +++ b/include/openmc/scattdata.h @@ -4,7 +4,7 @@ #ifndef OPENMC_SCATTDATA_H #define OPENMC_SCATTDATA_H -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include "openmc/constants.h" #include "openmc/vector.h" @@ -26,23 +26,23 @@ class ScattData { protected: //! \brief Initializes the attributes of the base class. - void base_init(int order, const xt::xtensor& in_gmin, - const xt::xtensor& in_gmax, const double_2dvec& in_energy, + void base_init(int order, const tensor::Tensor& in_gmin, + const tensor::Tensor& in_gmax, const double_2dvec& in_energy, const double_2dvec& in_mult); //! \brief Combines microscopic ScattDatas into a macroscopic one. void base_combine(size_t max_order, size_t order_dim, const vector& those_scatts, const vector& scalars, - xt::xtensor& in_gmin, xt::xtensor& in_gmax, + tensor::Tensor& in_gmin, tensor::Tensor& 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 - xt::xtensor gmin; // minimum outgoing group - xt::xtensor gmax; // maximum outgoing group - xt::xtensor scattxs; // Isotropic Sigma_{s,g_{in}} + tensor::Tensor gmin; // minimum outgoing group + tensor::Tensor gmax; // maximum outgoing group + tensor::Tensor scattxs; // Isotropic Sigma_{s,g_{in}} //! \brief Calculates the value of normalized f(mu). //! @@ -72,8 +72,8 @@ class ScattData { //! @param in_gmax List of maximum outgoing groups for every incoming group //! @param in_mult Input sparse multiplicity matrix //! @param coeffs Input sparse scattering matrix - virtual void init(const xt::xtensor& in_gmin, - const xt::xtensor& in_gmax, const double_2dvec& in_mult, + virtual void init(const tensor::Tensor& in_gmin, + const tensor::Tensor& in_gmax, const double_2dvec& in_mult, const double_3dvec& coeffs) = 0; //! \brief Combines the microscopic data. @@ -96,7 +96,7 @@ class ScattData { //! @param max_order If Legendre this is the maximum value of "n" in "Pn" //! requested; ignored otherwise. //! @return The dense scattering matrix. - virtual xt::xtensor get_matrix(size_t max_order) = 0; + virtual tensor::Tensor get_matrix(size_t max_order) = 0; //! \brief Samples the outgoing energy from the ScattData info. //! @@ -135,8 +135,8 @@ class ScattDataLegendre : public ScattData { ScattDataLegendre& leg, ScattDataTabular& tab); public: - void init(const xt::xtensor& in_gmin, - const xt::xtensor& in_gmax, const double_2dvec& in_mult, + void init(const tensor::Tensor& in_gmin, + const tensor::Tensor& in_gmax, const double_2dvec& in_mult, const double_3dvec& coeffs) override; void combine(const vector& those_scatts, @@ -153,7 +153,7 @@ class ScattDataLegendre : public ScattData { size_t get_order() override { return dist[0][0].size() - 1; }; - xt::xtensor get_matrix(size_t max_order) override; + tensor::Tensor get_matrix(size_t max_order) override; }; //============================================================================== @@ -164,13 +164,13 @@ class ScattDataLegendre : public ScattData { class ScattDataHistogram : public ScattData { protected: - xt::xtensor mu; // Angle distribution mu bin boundaries + tensor::Tensor 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 xt::xtensor& in_gmin, - const xt::xtensor& in_gmax, const double_2dvec& in_mult, + void init(const tensor::Tensor& in_gmin, + const tensor::Tensor& in_gmax, const double_2dvec& in_mult, const double_3dvec& coeffs) override; void combine(const vector& those_scatts, @@ -183,7 +183,7 @@ class ScattDataHistogram : public ScattData { size_t get_order() override { return dist[0][0].size(); }; - xt::xtensor get_matrix(size_t max_order) override; + tensor::Tensor get_matrix(size_t max_order) override; }; //============================================================================== @@ -194,7 +194,7 @@ class ScattDataHistogram : public ScattData { class ScattDataTabular : public ScattData { protected: - xt::xtensor mu; // Angle distribution mu grid points + tensor::Tensor mu; // Angle distribution mu grid points double dmu; // Quick storage of the mu spacing double_3dvec fmu; // The angular distribution function @@ -204,8 +204,8 @@ class ScattDataTabular : public ScattData { ScattDataLegendre& leg, ScattDataTabular& tab); public: - void init(const xt::xtensor& in_gmin, - const xt::xtensor& in_gmax, const double_2dvec& in_mult, + void init(const tensor::Tensor& in_gmin, + const tensor::Tensor& in_gmax, const double_2dvec& in_mult, const double_3dvec& coeffs) override; void combine(const vector& those_scatts, @@ -218,7 +218,7 @@ class ScattDataTabular : public ScattData { size_t get_order() override { return dist[0][0].size(); }; - xt::xtensor get_matrix(size_t max_order) override; + tensor::Tensor get_matrix(size_t max_order) override; }; //============================================================================== diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h index 6905c38e369..b4b7f8480f5 100644 --- a/include/openmc/secondary_correlated.h +++ b/include/openmc/secondary_correlated.h @@ -5,7 +5,7 @@ #define OPENMC_SECONDARY_CORRELATED_H #include "hdf5.h" -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include "openmc/angle_energy.h" #include "openmc/distribution.h" @@ -25,9 +25,9 @@ class CorrelatedAngleEnergy : public AngleEnergy { struct CorrTable { int n_discrete; //!< Number of discrete lines Interpolation interpolation; //!< Interpolation law - xt::xtensor e_out; //!< Outgoing energies [eV] - xt::xtensor p; //!< Probability density - xt::xtensor c; //!< Cumulative distribution + tensor::Tensor e_out; //!< Outgoing energies [eV] + tensor::Tensor p; //!< Probability density + tensor::Tensor c; //!< Cumulative distribution vector> angle; //!< Angle distribution }; diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h index 83806d35248..c9c5849bc77 100644 --- a/include/openmc/secondary_kalbach.h +++ b/include/openmc/secondary_kalbach.h @@ -5,7 +5,7 @@ #define OPENMC_SECONDARY_KALBACH_H #include "hdf5.h" -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include "openmc/angle_energy.h" #include "openmc/constants.h" @@ -37,11 +37,11 @@ class KalbachMann : public AngleEnergy { struct KMTable { int n_discrete; //!< Number of discrete lines Interpolation interpolation; //!< Interpolation law - xt::xtensor e_out; //!< Outgoing energies [eV] - xt::xtensor p; //!< Probability density - xt::xtensor c; //!< Cumulative distribution - xt::xtensor r; //!< Pre-compound fraction - xt::xtensor a; //!< Parameterized function + tensor::Tensor e_out; //!< Outgoing energies [eV] + tensor::Tensor p; //!< Probability density + tensor::Tensor c; //!< Cumulative distribution + tensor::Tensor r; //!< Pre-compound fraction + tensor::Tensor a; //!< Parameterized function }; int n_region_; //!< Number of interpolation regions diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 5b18902afbb..4f33c0e763c 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -9,7 +9,7 @@ #include "openmc/secondary_correlated.h" #include "openmc/vector.h" -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include namespace openmc { @@ -83,7 +83,7 @@ class IncoherentElasticAEDiscrete : public AngleEnergy { private: const vector& energy_; //!< Energies at which cosines are tabulated - xt::xtensor mu_out_; //!< Cosines for each incident energy + tensor::Tensor mu_out_; //!< Cosines for each incident energy }; //============================================================================== @@ -109,9 +109,9 @@ class IncoherentInelasticAEDiscrete : public AngleEnergy { private: const vector& energy_; //!< Incident energies - xt::xtensor + tensor::Tensor energy_out_; //!< Outgoing energies for each incident energy - xt::xtensor + tensor::Tensor mu_out_; //!< Outgoing cosines for each incident/outgoing energy bool skewed_; //!< Whether outgoing energy distribution is skewed }; @@ -139,10 +139,10 @@ class IncoherentInelasticAE : public AngleEnergy { //! Secondary energy/angle distribution struct DistEnergySab { std::size_t n_e_out; //!< Number of outgoing energies - xt::xtensor e_out; //!< Outgoing energies - xt::xtensor e_out_pdf; //!< Probability density function - xt::xtensor e_out_cdf; //!< Cumulative distribution function - xt::xtensor mu; //!< Equiprobable angles at each outgoing energy + tensor::Tensor e_out; //!< Outgoing energies + tensor::Tensor e_out_pdf; //!< Probability density function + tensor::Tensor e_out_cdf; //!< Cumulative distribution function + tensor::Tensor mu; //!< Equiprobable angles at each outgoing energy }; vector energy_; //!< Incident energies diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 374daff92a0..6161f823be1 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -8,9 +8,8 @@ #include "openmc/tallies/trigger.h" #include "openmc/vector.h" +#include "openmc/tensor.h" #include "pugixml.hpp" -#include "xtensor/xfixed.hpp" -#include "xtensor/xtensor.hpp" #include #include @@ -55,7 +54,7 @@ class Tally { void set_nuclides(const vector& nuclides); - const xt::xtensor& results() const { return results_; } + const tensor::Tensor& results() const { return results_; } //! returns vector of indices corresponding to the tally this is called on const vector& filters() const { return filters_; } @@ -125,7 +124,7 @@ class Tally { int score_index(const std::string& score) const; //! Tally results reshaped according to filter sizes - xt::xarray get_reshaped_data() const; + tensor::Tensor get_reshaped_data() const; //! A string representing the i-th score on this tally std::string score_name(int score_idx) const; @@ -160,7 +159,7 @@ class Tally { //! combination of filters (e.g. specific cell, specific energy group, etc.) //! and the second dimension of the array is for scores (e.g. flux, total //! reaction rate, fission reaction rate, etc.) - xt::xtensor results_; + tensor::Tensor results_; //! True if this tally should be written to statepoint files bool writable_ {true}; @@ -220,8 +219,7 @@ extern vector time_grid; namespace simulation { //! Global tallies (such as k-effective estimators) -extern xt::xtensor_fixed> - global_tallies; +extern tensor::StaticTensor2D global_tallies; //! Number of realizations for global tallies extern "C" int32_t n_realizations; @@ -257,12 +255,6 @@ double distance_to_time_boundary(double time, double speed); //! Determine which tallies should be active void setup_active_tallies(); -// Alias for the type returned by xt::adapt(...). N is the dimension of the -// multidimensional array -template -using adaptor_type = - xt::xtensor_adaptor, N>; - #ifdef OPENMC_MPI //! Collect all tally results onto master process void reduce_tally_results(); diff --git a/include/openmc/tensor.h b/include/openmc/tensor.h new file mode 100644 index 00000000000..9d428aaebdd --- /dev/null +++ b/include/openmc/tensor.h @@ -0,0 +1,1185 @@ +//! \file tensor.h +//! \brief Multi-dimensional tensor types for OpenMC. +//! +//! Tensor is the primary type: a dynamic-rank owning container that stores +//! elements contiguously in row-major order. View is a lightweight +//! non-owning reference into a Tensor's storage, returned by the slice() +//! method and flat(). StaticTensor2D is a small stack-allocated 2D +//! array used only for simulation::global_tallies. +//! +//! Slicing follows numpy conventions: each axis takes an index (rank-reducing), +//! All (keep entire axis), or Range (keep sub-range). For example, +//! arr.slice(0, all, range(2, 5)) is equivalent to numpy's arr[0, :, 2:5]. +//! +//! View is declared before Tensor because Tensor's methods return View objects. + +#ifndef OPENMC_TENSOR_H +#define OPENMC_TENSOR_H + +#include "openmc/vector.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace openmc { +namespace tensor { + +//============================================================================== +// Forward declarations +//============================================================================== + +template +class Tensor; + +template +class StaticTensor2D; + +//============================================================================== +// Storage type mapping +// +// std::vector is a bit-packed specialization that returns proxy objects +// instead of real references, which breaks generic code. storage_type_map +// redirects bool to unsigned char so that Tensor stores one byte per +// element with normal reference semantics. +//============================================================================== + +template +struct storage_type_map { + using type = T; +}; +template<> +struct storage_type_map { + using type = unsigned char; +}; +template +using storage_type = typename storage_type_map::type; + +//============================================================================== +// Slice argument types +// +// Used with the variadic slice() method on Tensor, View, and StaticTensor2D. +// Each argument corresponds to one axis: a plain integer fixes that axis at +// a single index (rank-reducing), All keeps the entire axis, and Range keeps +// a sub-range. +//============================================================================== + +//! Keep an entire axis (equivalent to numpy's ':' or xtensor's xt::all()) +struct All {}; +constexpr All all {}; + +//! Sub-range along an axis [start, end) +struct Range { + size_t start; + size_t end; // SIZE_MAX means "to end of axis" +}; + +//! Create a Range [start, end) +inline Range range(size_t start, size_t end) +{ + return {start, end}; +} + +//! Create a Range [0, end) +inline Range range(size_t end) +{ + return {0, end}; +} + +namespace detail { + +//! Internal: normalized representation of a per-axis slice argument +struct SliceArg { + enum Kind { INDEX, ALL, RANGE } kind; + size_t start; + size_t end; +}; + +inline SliceArg to_slice_arg(All) +{ + return {SliceArg::ALL, 0, 0}; +} +inline SliceArg to_slice_arg(Range r) +{ + return {SliceArg::RANGE, r.start, r.end}; +} + +template +inline + typename std::enable_if::value || std::is_enum::value, + SliceArg>::type + to_slice_arg(I i) +{ + return {SliceArg::INDEX, static_cast(i), 0}; +} + +//! Result of a slice computation: pointer offset + new shape/strides +struct SliceResult { + size_t ptr_offset; + vector shape; + vector strides; +}; + +//! Compute the result of applying slice arguments to shape/strides +template +SliceResult compute_slice(const vector& shape, + const vector& strides, First first, Rest... rest) +{ + const size_t n = 1 + sizeof...(Rest); + SliceArg args[1 + sizeof...(Rest)] = { + to_slice_arg(first), to_slice_arg(rest)...}; + + size_t offset = 0; + vector new_shape; + vector new_strides; + + for (size_t a = 0; a < n; ++a) { + switch (args[a].kind) { + case SliceArg::INDEX: + offset += args[a].start * strides[a]; + break; + case SliceArg::ALL: + new_shape.push_back(shape[a]); + new_strides.push_back(strides[a]); + break; + case SliceArg::RANGE: { + offset += args[a].start * strides[a]; + size_t end = (args[a].end == SIZE_MAX) ? shape[a] : args[a].end; + new_shape.push_back(end - args[a].start); + new_strides.push_back(strides[a]); + break; + } + } + } + + // Trailing axes not covered by arguments are implicitly All. + // This matches numpy: a[i] on a 2D array returns a 1D row. + for (size_t a = n; a < shape.size(); ++a) { + new_shape.push_back(shape[a]); + new_strides.push_back(strides[a]); + } + + return {offset, std::move(new_shape), std::move(new_strides)}; +} + +} // namespace detail + +//============================================================================== +// View: a non-owning N-dimensional view into a tensor's storage. +// +// Holds a base pointer, shape, and strides (in elements). Supports arbitrary +// rank and multi-axis slicing via the variadic slice() method. +//============================================================================== + +template +class View { +public: + //-------------------------------------------------------------------------- + // Constructors + + View(T* data, vector shape, vector strides) + : data_(data), shape_(std::move(shape)), strides_(std::move(strides)) + {} + + // Explicitly default copy/move constructors (declaring copy assignment + // below would otherwise suppress the implicit move constructor). + View(const View&) = default; + View(View&&) = default; + + //-------------------------------------------------------------------------- + // Indexing + + //! Multi-index element access (1D, 2D, 3D, ...) + template + T& operator()(Indices... indices) + { + const size_t idx[] = {static_cast(indices)...}; + size_t off = 0; + for (size_t d = 0; d < sizeof...(Indices); ++d) + off += idx[d] * strides_[d]; + return data_[off]; + } + + template + const T& operator()(Indices... indices) const + { + const size_t idx[] = {static_cast(indices)...}; + size_t off = 0; + for (size_t d = 0; d < sizeof...(Indices); ++d) + off += idx[d] * strides_[d]; + return data_[off]; + } + + //! Flat logical index (row-major order) + T& operator[](size_t i) { return data_[flat_to_offset(i)]; } + const T& operator[](size_t i) const { return data_[flat_to_offset(i)]; } + + //-------------------------------------------------------------------------- + // Accessors + + size_t size() const + { + size_t s = 1; + for (auto d : shape_) + s *= d; + return s; + } + size_t ndim() const { return shape_.size(); } + size_t shape(size_t axis) const { return shape_[axis]; } + const vector& shape_vec() const { return shape_; } + T* data() { return data_; } + const T* data() const { return data_; } + + //-------------------------------------------------------------------------- + // View accessors + + //! Multi-axis slice. Each argument corresponds to one axis and is either: + //! - an integer (fixes that axis, rank-reducing) + //! - All (keeps entire axis) + //! - Range (keeps sub-range along that axis) + //! Example: v.slice(0, all, range(2, 5)) == numpy v[0, :, 2:5] + template + View slice(First first, Rest... rest) + { + auto r = detail::compute_slice(shape_, strides_, first, rest...); + return {data_ + r.ptr_offset, std::move(r.shape), std::move(r.strides)}; + } + + template + View slice(First first, Rest... rest) const + { + auto r = detail::compute_slice(shape_, strides_, first, rest...); + return {data_ + r.ptr_offset, std::move(r.shape), std::move(r.strides)}; + } + + //-------------------------------------------------------------------------- + // Assignment operators + + //! Copy assignment: element-wise deep copy (writes through data pointer). + //! Without this, the compiler's implicit copy assignment just copies the + //! View metadata (pointer, shape, strides) instead of the viewed data. + View& operator=(const View& other) + { + size_t n = size(); + for (size_t i = 0; i < n; ++i) + data_[flat_to_offset(i)] = other[i]; + return *this; + } + + //! Fill all elements with a scalar + View& operator=(T val) + { + size_t n = size(); + for (size_t i = 0; i < n; ++i) + data_[flat_to_offset(i)] = val; + return *this; + } + + //! Assignment from initializer_list (for 1D views) + View& operator=(std::initializer_list vals) + { + auto it = vals.begin(); + for (size_t i = 0; i < size() && it != vals.end(); ++i, ++it) + data_[flat_to_offset(i)] = *it; + return *this; + } + + //! Assignment from another View + template + View& operator=(const View& other) + { + size_t n = size(); + for (size_t i = 0; i < n; ++i) + data_[flat_to_offset(i)] = other[i]; + return *this; + } + + //! Assignment from Tensor (forward-declared, defined after Tensor) + template + View& operator=(const Tensor& other); + + //! Compound addition from Tensor (forward-declared, defined after Tensor) + template + View& operator+=(const Tensor& o); + + //! Compound multiply by scalar + View& operator*=(T val) + { + size_t n = size(); + for (size_t i = 0; i < n; ++i) + data_[flat_to_offset(i)] *= val; + return *this; + } + + //! Compound divide by scalar + View& operator/=(T val) + { + size_t n = size(); + for (size_t i = 0; i < n; ++i) + data_[flat_to_offset(i)] /= val; + return *this; + } + + //-------------------------------------------------------------------------- + // Reductions + + //! Sum of all elements + T sum() const + { + // remove_const needed so accumulator is mutable when T is const-qualified + std::remove_const_t s = 0; + size_t n = size(); + for (size_t i = 0; i < n; ++i) + s += data_[flat_to_offset(i)]; + return s; + } + + //-------------------------------------------------------------------------- + // Iterators + // + // Lightweight row-major iterator parameterized on pointer type (Ptr). + // Stores a flat logical position and converts to a physical offset on each + // dereference via divmod over shape/strides. For contiguous 1D views (the + // common case) the divmod chain reduces to a single multiply-by-1, which + // the compiler optimizes away. + // + // view_iterator = mutable iterator (from non-const View) + // view_iterator = read-only iterator (from const View) + + template + class view_iterator { + Ptr base_; + size_t count_; + const size_t* shape_; + const size_t* strides_; + size_t ndim_; + + public: + using iterator_category = std::random_access_iterator_tag; + using value_type = std::remove_const_t; + using difference_type = std::ptrdiff_t; + using pointer = Ptr; + using reference = decltype(*std::declval()); + + view_iterator(Ptr base, size_t count, const View* v) + : base_(base), count_(count), shape_(v->shape_.data()), + strides_(v->strides_.data()), ndim_(v->shape_.size()) + {} + + reference operator*() const { return base_[offset()]; } + reference operator[](difference_type n) const + { + return base_[offset_of(count_ + n)]; + } + view_iterator& operator++() + { + ++count_; + return *this; + } + view_iterator operator++(int) + { + auto tmp = *this; + ++count_; + return tmp; + } + view_iterator& operator--() + { + --count_; + return *this; + } + view_iterator operator+(difference_type n) const + { + auto tmp = *this; + tmp.count_ += n; + return tmp; + } + view_iterator operator-(difference_type n) const + { + auto tmp = *this; + tmp.count_ -= n; + return tmp; + } + difference_type operator-(const view_iterator& o) const + { + return static_cast(count_) - + static_cast(o.count_); + } + view_iterator& operator+=(difference_type n) + { + count_ += n; + return *this; + } + view_iterator& operator-=(difference_type n) + { + count_ -= n; + return *this; + } + bool operator==(const view_iterator& o) const { return count_ == o.count_; } + bool operator!=(const view_iterator& o) const { return count_ != o.count_; } + bool operator<(const view_iterator& o) const { return count_ < o.count_; } + bool operator>(const view_iterator& o) const { return count_ > o.count_; } + bool operator<=(const view_iterator& o) const { return count_ <= o.count_; } + bool operator>=(const view_iterator& o) const { return count_ >= o.count_; } + friend view_iterator operator+(difference_type n, const view_iterator& it) + { + return it + n; + } + + private: + size_t offset() const { return offset_of(count_); } + size_t offset_of(size_t flat) const + { + size_t off = 0; + for (int d = static_cast(ndim_) - 1; d >= 0; --d) { + off += (flat % shape_[d]) * strides_[d]; + flat /= shape_[d]; + } + return off; + } + }; + + using iterator = view_iterator; + using const_iterator = view_iterator; + + iterator begin() { return {data_, 0, this}; } + iterator end() { return {data_, size(), this}; } + const_iterator begin() const { return cbegin(); } + const_iterator end() const { return cend(); } + const_iterator cbegin() const { return {data_, 0, this}; } + const_iterator cend() const { return {data_, size(), this}; } + +private: + //! Convert a logical flat index (row-major) to a physical element offset + size_t flat_to_offset(size_t flat) const + { + size_t off = 0; + for (int d = static_cast(shape_.size()) - 1; d >= 0; --d) { + off += (flat % shape_[d]) * strides_[d]; + flat /= shape_[d]; + } + return off; + } + + T* data_; + vector shape_; + vector strides_; +}; + +//============================================================================== +// Tensor: dynamic-rank N-dimensional tensor. +// +// Stores elements in a contiguous row-major vector> +// with a dynamic shape. +//============================================================================== + +template +class Tensor { +public: + using value_type = T; + using stored_type = storage_type; + using iterator = typename vector::iterator; + using const_iterator = typename vector::const_iterator; + + //-------------------------------------------------------------------------- + // Constructors + + Tensor() = default; + + //! Construct with shape (uninitialized for arithmetic types via vector + //! resize) + explicit Tensor(vector shape) + : shape_(std::move(shape)), data_(compute_size()) + {} + + //! Construct with shape and fill value + Tensor(vector shape, T fill) + : shape_(std::move(shape)), data_(compute_size(), fill) + {} + + //! Construct from initializer_list shape + explicit Tensor(std::initializer_list shape) + : shape_(shape), data_(compute_size()) + {} + + //! Construct from initializer_list shape with fill + Tensor(std::initializer_list shape, T fill) + : shape_(shape), data_(compute_size(), fill) + {} + + //! 1D copy from raw pointer + count + Tensor(const T* ptr, size_t count) : shape_({count}), data_(ptr, ptr + count) + {} + + //! Copy from View (preserves view's shape) + template + explicit Tensor(const View& v) : shape_(v.shape_vec()) + { + size_t n = v.size(); + data_.resize(n); + for (size_t i = 0; i < n; ++i) + data_[i] = v[i]; + } + + //-------------------------------------------------------------------------- + // Assignment + + //! Assignment from View + template + Tensor& operator=(const View& v) + { + shape_ = v.shape_vec(); + size_t n = v.size(); + data_.resize(n); + for (size_t i = 0; i < n; ++i) + data_[i] = v[i]; + return *this; + } + + //! Assignment from initializer_list of values (1D) + Tensor& operator=(std::initializer_list vals) + { + shape_ = {vals.size()}; + data_.assign(vals.begin(), vals.end()); + return *this; + } + + //-------------------------------------------------------------------------- + // Accessors + + stored_type* data() { return data_.data(); } + const stored_type* data() const { return data_.data(); } + size_t size() const { return data_.size(); } + const vector& shape() const { return shape_; } + size_t shape(size_t dim) const + { + return dim < shape_.size() ? shape_[dim] : 0; + } + size_t ndim() const { return shape_.size(); } + bool empty() const { return data_.empty(); } + + //-------------------------------------------------------------------------- + // Indexing (row-major) + + template + stored_type& operator()(Indices... indices) + { + const size_t idx[] = {static_cast(indices)...}; + size_t off = 0; + for (size_t d = 0; d < sizeof...(Indices); ++d) + off = off * shape_[d] + idx[d]; + return data_[off]; + } + + template + const stored_type& operator()(Indices... indices) const + { + const size_t idx[] = {static_cast(indices)...}; + size_t off = 0; + for (size_t d = 0; d < sizeof...(Indices); ++d) + off = off * shape_[d] + idx[d]; + return data_[off]; + } + + stored_type& operator[](size_t i) { return data_[i]; } + const stored_type& operator[](size_t i) const { return data_[i]; } + + //! First and last element + stored_type& front() { return data_.front(); } + const stored_type& front() const { return data_.front(); } + stored_type& back() { return data_.back(); } + const stored_type& back() const { return data_.back(); } + + //-------------------------------------------------------------------------- + // Iterators + + iterator begin() { return data_.begin(); } + iterator end() { return data_.end(); } + const_iterator begin() const { return data_.begin(); } + const_iterator end() const { return data_.end(); } + const_iterator cbegin() const { return data_.cbegin(); } + const_iterator cend() const { return data_.cend(); } + + //-------------------------------------------------------------------------- + // Mutation + + void resize(const vector& shape) + { + shape_ = shape; + data_.resize(compute_size()); + } + + void resize(std::initializer_list shape) + { + shape_.assign(shape.begin(), shape.end()); + data_.resize(compute_size()); + } + + void reshape(const vector& new_shape) { shape_ = new_shape; } + + void fill(T val) { std::fill(data_.begin(), data_.end(), val); } + + //-------------------------------------------------------------------------- + // View accessors + + //! Fix one axis at a given index, returning an (N-1)-dimensional view + //! Multi-axis slice. Each argument corresponds to one axis and is either: + //! - an integer (fixes that axis, rank-reducing) + //! - All (keeps entire axis) + //! - Range (keeps sub-range along that axis) + //! Example: t.slice(0, all, range(2, 5)) == numpy t[0, :, 2:5] + template + View slice(First first, Rest... rest) + { + auto strides = compute_strides(); + auto r = detail::compute_slice(shape_, strides, first, rest...); + return { + data_.data() + r.ptr_offset, std::move(r.shape), std::move(r.strides)}; + } + + template + View slice(First first, Rest... rest) const + { + auto strides = compute_strides(); + auto r = detail::compute_slice(shape_, strides, first, rest...); + return { + data_.data() + r.ptr_offset, std::move(r.shape), std::move(r.strides)}; + } + + //! Flat 1D view of all elements + View flat() + { + return {data_.data(), {data_.size()}, {size_t(1)}}; + } + View flat() const + { + return {data_.data(), {data_.size()}, {size_t(1)}}; + } + + //-------------------------------------------------------------------------- + // Reductions and transforms + + //! Sum of all elements + T sum() const + { + T s = T(0); + for (size_t i = 0; i < data_.size(); ++i) + s += data_[i]; + return s; + } + + //! Sum along an axis, reducing rank by 1 (defined out-of-line below) + Tensor sum(size_t axis) const; + + //! Product of all elements + T prod() const + { + T p = T(1); + for (size_t i = 0; i < data_.size(); ++i) + p *= data_[i]; + return p; + } + + //! True if any element is nonzero + bool any() const + { + for (size_t i = 0; i < data_.size(); ++i) + if (data_[i]) + return true; + return false; + } + + //! True if all elements are nonzero + bool all() const + { + for (size_t i = 0; i < data_.size(); ++i) + if (!data_[i]) + return false; + return true; + } + + //! Flat index of the minimum element + size_t argmin() const + { + return static_cast(std::distance(data_.data(), + std::min_element(data_.data(), data_.data() + data_.size()))); + } + + //! Reverse element order along an axis (e.g. flip(0) reverses rows) + Tensor flip(size_t axis) const + { + size_t outer_size = 1; + for (size_t d = 0; d < axis; ++d) + outer_size *= shape_[d]; + size_t axis_size = shape_[axis]; + size_t inner_size = 1; + for (size_t d = axis + 1; d < shape_.size(); ++d) + inner_size *= shape_[d]; + + Tensor r(shape_); + for (size_t o = 0; o < outer_size; ++o) + for (size_t a = 0; a < axis_size; ++a) + for (size_t i = 0; i < inner_size; ++i) + r.data_[(o * axis_size + (axis_size - 1 - a)) * inner_size + i] = + data_[(o * axis_size + a) * inner_size + i]; + return r; + } + + //-------------------------------------------------------------------------- + // Operators + + Tensor& operator+=(T val) + { + for (auto& x : data_) + x += val; + return *this; + } + Tensor& operator-=(T val) + { + for (auto& x : data_) + x -= val; + return *this; + } + Tensor& operator*=(T val) + { + for (auto& x : data_) + x *= val; + return *this; + } + Tensor& operator/=(T val) + { + for (auto& x : data_) + x /= val; + return *this; + } + Tensor& operator+=(const Tensor& o) + { + for (size_t i = 0; i < data_.size(); ++i) + data_[i] += o.data_[i]; + return *this; + } + + Tensor operator+(const Tensor& o) const + { + Tensor r(shape_); + for (size_t i = 0; i < data_.size(); ++i) + r.data_[i] = data_[i] + o.data_[i]; + return r; + } + Tensor operator-(const Tensor& o) const + { + Tensor r(shape_); + for (size_t i = 0; i < data_.size(); ++i) + r.data_[i] = data_[i] - o.data_[i]; + return r; + } + Tensor operator/(const Tensor& o) const + { + Tensor r(shape_); + for (size_t i = 0; i < data_.size(); ++i) + r.data_[i] = data_[i] / o.data_[i]; + return r; + } + + Tensor operator+(T val) const + { + Tensor r(shape_); + for (size_t i = 0; i < data_.size(); ++i) + r.data_[i] = data_[i] + val; + return r; + } + Tensor operator-(T val) const + { + Tensor r(shape_); + for (size_t i = 0; i < data_.size(); ++i) + r.data_[i] = data_[i] - val; + return r; + } + Tensor operator*(T val) const + { + Tensor r(shape_); + for (size_t i = 0; i < data_.size(); ++i) + r.data_[i] = data_[i] * val; + return r; + } + + Tensor operator<=(T val) const + { + Tensor r(shape_); + for (size_t i = 0; i < data_.size(); ++i) + r.data()[i] = data_[i] <= val; + return r; + } + Tensor operator<(T val) const + { + Tensor r(shape_); + for (size_t i = 0; i < data_.size(); ++i) + r.data()[i] = data_[i] < val; + return r; + } + Tensor operator>=(T val) const + { + Tensor r(shape_); + for (size_t i = 0; i < data_.size(); ++i) + r.data()[i] = data_[i] >= val; + return r; + } + Tensor operator>(T val) const + { + Tensor r(shape_); + for (size_t i = 0; i < data_.size(); ++i) + r.data()[i] = data_[i] > val; + return r; + } + Tensor operator<(const Tensor& o) const + { + Tensor r(shape_); + for (size_t i = 0; i < data_.size(); ++i) + r.data()[i] = data_[i] < o.data_[i]; + return r; + } + +private: + size_t compute_size() const + { + size_t s = 1; + for (auto d : shape_) + s *= d; + return s; + } + + //! Compute row-major strides from shape + vector compute_strides() const + { + vector strides(shape_.size()); + if (!shape_.empty()) { + strides.back() = 1; + for (int d = static_cast(shape_.size()) - 2; d >= 0; --d) + strides[d] = strides[d + 1] * shape_[d + 1]; + } + return strides; + } + + //-------------------------------------------------------------------------- + // Data members + + vector shape_; + vector> data_; +}; + +//============================================================================== +// Non-member operators (scalar op tensor) +//============================================================================== + +template +Tensor operator*(T val, const Tensor& arr) +{ + return arr * val; +} + +template +Tensor operator+(T val, const Tensor& arr) +{ + return arr + val; +} + +// Mixed-type arithmetic: Tensor op Tensor -> Tensor +// A SFINAE guard is used here, as without !is_same Tensor * Tensor +// would be ambiguous between the member operator* and this non-member function. +template::value>> +Tensor operator*(const Tensor& a, const Tensor& b) +{ + Tensor r(a.shape()); + for (size_t i = 0; i < a.size(); ++i) + r.data()[i] = + static_cast(a.data()[i]) * static_cast(b.data()[i]); + return r; +} + +// Same SFINAE guard as operator* above. +template::value>> +Tensor operator/(const Tensor& a, const Tensor& b) +{ + Tensor r(a.shape()); + for (size_t i = 0; i < a.size(); ++i) + r.data()[i] = + static_cast(a.data()[i]) / static_cast(b.data()[i]); + return r; +} + +//============================================================================== +// Out-of-line method definitions (require complete types) +//============================================================================== + +template +template +View& View::operator=(const Tensor& other) +{ + size_t n = size(); + for (size_t i = 0; i < n; ++i) + data_[flat_to_offset(i)] = static_cast(other.data()[i]); + return *this; +} + +template +template +View& View::operator+=(const Tensor& o) +{ + size_t n = size(); + for (size_t i = 0; i < n; ++i) + data_[flat_to_offset(i)] += o.data()[i]; + return *this; +} + +template +Tensor Tensor::sum(size_t axis) const +{ + // Build output shape (all dims except the summed axis) + vector out_shape; + for (size_t d = 0; d < shape_.size(); ++d) + if (d != axis) + out_shape.push_back(shape_[d]); + + // Split dimensions into three zones: outer | axis | inner + size_t outer_size = 1; + for (size_t d = 0; d < axis; ++d) + outer_size *= shape_[d]; + size_t axis_size = shape_[axis]; + size_t inner_size = 1; + for (size_t d = axis + 1; d < shape_.size(); ++d) + inner_size *= shape_[d]; + + Tensor result(out_shape, T(0)); + for (size_t o = 0; o < outer_size; ++o) + for (size_t a = 0; a < axis_size; ++a) + for (size_t i = 0; i < inner_size; ++i) + result.data()[o * inner_size + i] += + data_[(o * axis_size + a) * inner_size + i]; + + return result; +} + +//============================================================================== +// StaticTensor2D: compile-time fixed 2D tensor. +//============================================================================== + +template +class StaticTensor2D { +public: + using value_type = T; + + //-------------------------------------------------------------------------- + // Indexing + + //! Templated to accept enum class indices (e.g. GlobalTally, TallyResult) + //! which don't implicitly convert to integer types. + template + T& operator()(I0 i, I1 j) + { + return data_[static_cast(i) * C + static_cast(j)]; + } + template + const T& operator()(I0 i, I1 j) const + { + return data_[static_cast(i) * C + static_cast(j)]; + } + + //-------------------------------------------------------------------------- + // Accessors + + T* data() { return data_; } + const T* data() const { return data_; } + constexpr size_t size() const { return R * C; } + std::array shape() const { return {R, C}; } + + //-------------------------------------------------------------------------- + // Mutation + + void fill(T val) { std::fill(data_, data_ + R * C, val); } + + //-------------------------------------------------------------------------- + // Iterators + + T* begin() { return data_; } + T* end() { return data_ + R * C; } + const T* begin() const { return data_; } + const T* end() const { return data_ + R * C; } + + //-------------------------------------------------------------------------- + // View accessors + + //! Multi-axis slice (same interface as Tensor/View). + template + View slice(First first, Rest... rest) + { + vector sh = {R, C}; + vector st = {C, 1}; + auto r = detail::compute_slice(sh, st, first, rest...); + return {data_ + r.ptr_offset, std::move(r.shape), std::move(r.strides)}; + } + template + View slice(First first, Rest... rest) const + { + vector sh = {R, C}; + vector st = {C, 1}; + auto r = detail::compute_slice(sh, st, first, rest...); + return {data_ + r.ptr_offset, std::move(r.shape), std::move(r.strides)}; + } + + //! Flat view (1D, contiguous) + View flat() { return {data_, {R * C}, {size_t(1)}}; } + View flat() const { return {data_, {R * C}, {size_t(1)}}; } + +private: + //-------------------------------------------------------------------------- + // Data members + + T data_[R * C] = {}; +}; + +//============================================================================== +// Non-member functions +//============================================================================== + +template +Tensor zeros(std::initializer_list shape) +{ + vector s(shape); + return Tensor(std::move(s), T(0)); +} + +template +Tensor zeros(const vector& shape) +{ + return Tensor(shape, T(0)); +} + +template +Tensor ones(std::initializer_list shape) +{ + vector s(shape); + return Tensor(std::move(s), T(1)); +} + +template +Tensor ones(const vector& shape) +{ + return Tensor(shape, T(1)); +} + +template +Tensor zeros_like(const Tensor& o) +{ + return Tensor(o.shape(), T(0)); +} + +template +Tensor full_like(const Tensor& o, V val) +{ + return Tensor(o.shape(), static_cast(val)); +} + +//! Return a 1D tensor of n evenly spaced values from start to stop (inclusive) +template +Tensor linspace(T start, T stop, size_t n) +{ + Tensor result({n}); + if (n < 2) { + result[0] = start; + return result; + } + for (size_t i = 0; i < n; ++i) { + result[i] = + start + static_cast(i) * (stop - start) / static_cast(n - 1); + } + return result; +} + +//! Concatenate two 1D tensors end-to-end +template +Tensor concatenate(const Tensor& a, const Tensor& b) +{ + size_t total = a.size() + b.size(); + Tensor result({total}); + std::copy(a.data(), a.data() + a.size(), result.data()); + std::copy(b.data(), b.data() + b.size(), result.data() + a.size()); + return result; +} + +//! Element-wise natural logarithm +template +Tensor log(const Tensor& a) +{ + Tensor r(a.shape()); + for (size_t i = 0; i < a.size(); ++i) + r.data()[i] = std::log(a.data()[i]); + return r; +} + +//! Element-wise absolute value +template +Tensor abs(const Tensor& a) +{ + Tensor r(a.shape()); + for (size_t i = 0; i < a.size(); ++i) + r.data()[i] = std::abs(a.data()[i]); + return r; +} + +//! Element-wise conditional: select from true_val where cond is true, +//! otherwise use false_val +template +Tensor where( + const Tensor& cond, const Tensor& true_val, V false_val) +{ + Tensor r(cond.shape()); + for (size_t i = 0; i < cond.size(); ++i) + r.data()[i] = + cond.data()[i] ? true_val.data()[i] : static_cast(false_val); + return r; +} + +//! Replace NaN/Inf values with finite substitutes +template +Tensor nan_to_num(const Tensor& a, T nan_val = T(0), + T posinf_val = std::numeric_limits::max(), + T neginf_val = std::numeric_limits::lowest()) +{ + Tensor r(a.shape()); + for (size_t i = 0; i < a.size(); ++i) { + T val = a.data()[i]; + if (std::isnan(val)) + r.data()[i] = nan_val; + else if (std::isinf(val)) + r.data()[i] = val > 0 ? posinf_val : neginf_val; + else + r.data()[i] = val; + } + return r; +} + +//============================================================================== +// Type traits +//============================================================================== + +//! Type trait that is true for Tensor and StaticTensor2D. +//! Used by hdf5_interface.h to select the correct write_dataset overload. +template +struct is_tensor : std::false_type {}; + +template +struct is_tensor> : std::true_type {}; + +template +struct is_tensor> : std::true_type {}; + +} // namespace tensor +} // namespace openmc + +#endif // OPENMC_TENSOR_H diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h index de0767d0af0..86254b92314 100644 --- a/include/openmc/thermal.h +++ b/include/openmc/thermal.h @@ -5,7 +5,7 @@ #include #include -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include "openmc/angle_energy.h" #include "openmc/endf.h" diff --git a/include/openmc/urr.h b/include/openmc/urr.h index 1e603715847..d40c2aff7e0 100644 --- a/include/openmc/urr.h +++ b/include/openmc/urr.h @@ -3,7 +3,7 @@ #ifndef OPENMC_URR_H #define OPENMC_URR_H -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include "openmc/constants.h" #include "openmc/hdf5_interface.h" @@ -40,11 +40,11 @@ class UrrData { * below, obviously, values of the CDF are stored. For the xs_values * variable, the columns line up with the index of cdf_values. */ - xt::xtensor cdf_values_; // Note: must be row major! - xt::xtensor xs_values_; + tensor::Tensor cdf_values_; // Note: must be row major! + tensor::Tensor xs_values_; // Number of points in the CDF - auto n_cdf() const { return cdf_values_.shape()[1]; } + auto n_cdf() const { return cdf_values_.shape(1); } //! \brief Load the URR data from the provided HDF5 group explicit UrrData(hid_t group_id); diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h index fa8d3d65ece..9d3f1d02615 100644 --- a/include/openmc/volume_calc.h +++ b/include/openmc/volume_calc.h @@ -12,8 +12,8 @@ #include "openmc/tallies/trigger.h" #include "openmc/vector.h" +#include "openmc/tensor.h" #include "pugixml.hpp" -#include "xtensor/xtensor.hpp" #ifdef _OPENMP #include #endif diff --git a/include/openmc/weight_windows.h b/include/openmc/weight_windows.h index 0d7435ca2a5..97abcfe82ef 100644 --- a/include/openmc/weight_windows.h +++ b/include/openmc/weight_windows.h @@ -143,10 +143,10 @@ class WeightWindows { const vector& energy_bounds() const { return energy_bounds_; } - void set_bounds(const xt::xtensor& lower_ww_bounds, - const xt::xtensor& upper_bounds); + void set_bounds(const tensor::Tensor& lower_ww_bounds, + const tensor::Tensor& upper_bounds); - void set_bounds(const xt::xtensor& lower_bounds, double ratio); + void set_bounds(const tensor::Tensor& lower_bounds, double ratio); void set_bounds( span lower_bounds, span upper_bounds); @@ -182,11 +182,11 @@ class WeightWindows { const std::unique_ptr& mesh() const { return model::meshes[mesh_idx_]; } - const xt::xtensor& lower_ww_bounds() const { return lower_ww_; } - xt::xtensor& lower_ww_bounds() { return lower_ww_; } + const tensor::Tensor& lower_ww_bounds() const { return lower_ww_; } + tensor::Tensor& lower_ww_bounds() { return lower_ww_; } - const xt::xtensor& upper_ww_bounds() const { return upper_ww_; } - xt::xtensor& upper_ww_bounds() { return upper_ww_; } + const tensor::Tensor& upper_ww_bounds() const { return upper_ww_; } + tensor::Tensor& upper_ww_bounds() { return upper_ww_; } ParticleType particle_type() const { return particle_type_; } @@ -197,9 +197,9 @@ class WeightWindows { int64_t index_; //!< Index into weight windows vector ParticleType particle_type_; //!< Particle type to apply weight windows to vector energy_bounds_; //!< Energy boundaries [eV] - xt::xtensor lower_ww_; //!< Lower weight window bounds (shape: + tensor::Tensor lower_ww_; //!< Lower weight window bounds (shape: //!< energy_bins, mesh_bins (k, j, i)) - xt::xtensor + tensor::Tensor upper_ww_; //!< Upper weight window bounds (shape: energy_bins, mesh_bins) double survival_ratio_ {3.0}; //!< Survival weight ratio double max_lb_ratio_ {1.0}; //!< Maximum lower bound to particle weight ratio diff --git a/include/openmc/wmp.h b/include/openmc/wmp.h index 6a4abd86714..5cc04e59594 100644 --- a/include/openmc/wmp.h +++ b/include/openmc/wmp.h @@ -2,7 +2,7 @@ #define OPENMC_WMP_H #include "hdf5.h" -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include #include @@ -78,9 +78,9 @@ class WindowedMultipole { int fit_order_; //!< Order of the fit bool fissionable_; //!< Is the nuclide fissionable? vector window_info_; // Information about a window - xt::xtensor + tensor::Tensor curvefit_; // Curve fit coefficients (window, poly order, reaction) - xt::xtensor, 2> data_; //!< Poles and residues + tensor::Tensor> data_; //!< Poles and residues // Constant data static constexpr int MAX_POLY_COEFFICIENTS = diff --git a/include/openmc/xml_interface.h b/include/openmc/xml_interface.h index f49613ecde1..17a34e5c77a 100644 --- a/include/openmc/xml_interface.h +++ b/include/openmc/xml_interface.h @@ -5,9 +5,8 @@ #include // for stringstream #include +#include "openmc/tensor.h" #include "pugixml.hpp" -#include "xtensor/xadapt.hpp" -#include "xtensor/xarray.hpp" #include "openmc/position.h" #include "openmc/vector.h" @@ -42,12 +41,11 @@ vector get_node_array( } template -xt::xarray get_node_xarray( +tensor::Tensor get_node_tensor( pugi::xml_node node, const char* name, bool lowercase = false) { vector v = get_node_array(node, name, lowercase); - vector shape = {v.size()}; - return xt::adapt(v, shape); + return tensor::Tensor(v.data(), v.size()); } std::vector get_node_position_array( diff --git a/include/openmc/xsdata.h b/include/openmc/xsdata.h index feafde68dd3..c9dbde986b2 100644 --- a/include/openmc/xsdata.h +++ b/include/openmc/xsdata.h @@ -4,7 +4,7 @@ #ifndef OPENMC_XSDATA_H #define OPENMC_XSDATA_H -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include "openmc/hdf5_interface.h" #include "openmc/memory.h" @@ -69,26 +69,26 @@ class XsData { public: // The following quantities have the following dimensions: // [angle][incoming group] - 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; + tensor::Tensor total; + tensor::Tensor absorption; + tensor::Tensor nu_fission; + tensor::Tensor prompt_nu_fission; + tensor::Tensor kappa_fission; + tensor::Tensor fission; + tensor::Tensor inverse_velocity; // decay_rate has the following dimensions: // [angle][delayed group] - xt::xtensor decay_rate; + tensor::Tensor decay_rate; // delayed_nu_fission has the following dimensions: // [angle][delayed group][incoming group] - xt::xtensor delayed_nu_fission; + tensor::Tensor delayed_nu_fission; // chi_prompt has the following dimensions: // [angle][incoming group][outgoing group] - xt::xtensor chi_prompt; + tensor::Tensor chi_prompt; // chi_delayed has the following dimensions: // [angle][incoming group][outgoing group][delayed group] - xt::xtensor chi_delayed; + tensor::Tensor chi_delayed; // scatter has the following dimensions: [angle] vector> scatter; diff --git a/src/bank.cpp b/src/bank.cpp index 33790379b85..5afdc2ecab7 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -7,6 +7,7 @@ #include "openmc/vector.h" #include +#include namespace openmc { diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index d77066fb0eb..ec1088101f0 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -6,7 +6,7 @@ #include "openmc/search.h" #include "openmc/settings.h" -#include "xtensor/xmath.hpp" +#include "openmc/tensor.h" namespace openmc { @@ -16,8 +16,8 @@ namespace openmc { namespace data { -xt::xtensor ttb_e_grid; -xt::xtensor ttb_k_grid; +tensor::Tensor ttb_e_grid; +tensor::Tensor ttb_k_grid; vector ttb; } // namespace data diff --git a/src/cmfd_solver.cpp b/src/cmfd_solver.cpp index 943042f67e7..714a5bf3ac5 100644 --- a/src/cmfd_solver.cpp +++ b/src/cmfd_solver.cpp @@ -5,7 +5,7 @@ #ifdef _OPENMP #include #endif -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include "openmc/bank.h" #include "openmc/capi.h" @@ -36,7 +36,7 @@ double spectral; int nx, ny, nz, ng; -xt::xtensor indexmap; +tensor::Tensor indexmap; int use_all_threads; @@ -79,15 +79,14 @@ int get_cmfd_energy_bin(const double E) // COUNT_BANK_SITES bins fission sites according to CMFD mesh and energy //============================================================================== -xt::xtensor count_bank_sites( - xt::xtensor& bins, bool* outside) +tensor::Tensor count_bank_sites( + tensor::Tensor& bins, bool* outside) { // Determine shape of array for counts std::size_t cnt_size = cmfd::nx * cmfd::ny * cmfd::nz * cmfd::ng; - vector cnt_shape = {cnt_size}; // Create array of zeros - xt::xarray cnt {cnt_shape, 0.0}; + tensor::Tensor cnt = tensor::zeros({cnt_size}); bool outside_ = false; auto bank_size = simulation::source_bank.size(); @@ -113,29 +112,22 @@ xt::xtensor count_bank_sites( bins[i] = mesh_bin * cmfd::ng + energy_bin; } - // Create copy of count data. Since ownership will be acquired by xtensor, - // std::allocator must be used to avoid Valgrind mismatched free() / delete - // warnings. int total = cnt.size(); - double* cnt_reduced = std::allocator {}.allocate(total); + tensor::Tensor counts = tensor::zeros({cnt_size}); #ifdef OPENMC_MPI // collect values from all processors MPI_Reduce( - cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); + cnt.data(), counts.data(), total, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); // Check if there were sites outside the mesh for any processor MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm); #else - std::copy(cnt.data(), cnt.data() + total, cnt_reduced); + std::copy(cnt.data(), cnt.data() + total, counts.data()); *outside = outside_; #endif - // Adapt reduced values in array back into an xarray - auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), cnt_shape); - xt::xarray counts = arr; - return counts; } @@ -151,19 +143,19 @@ extern "C" void openmc_cmfd_reweight( std::size_t src_size = cmfd::nx * cmfd::ny * cmfd::nz * cmfd::ng; // count bank sites for CMFD mesh, store bins in bank_bins for reweighting - xt::xtensor bank_bins({bank_size}, 0); + tensor::Tensor bank_bins = tensor::zeros({bank_size}); bool sites_outside; - xt::xtensor sourcecounts = + tensor::Tensor sourcecounts = count_bank_sites(bank_bins, &sites_outside); // Compute CMFD weightfactors - xt::xtensor weightfactors = xt::xtensor({src_size}, 1.); + tensor::Tensor weightfactors = tensor::ones({src_size}); if (mpi::master) { if (sites_outside) { fatal_error("Source sites outside of the CMFD mesh"); } - double norm = xt::sum(sourcecounts)() / cmfd::norm; + double norm = sourcecounts.sum() / cmfd::norm; for (int i = 0; i < src_size; i++) { if (sourcecounts[i] > 0 && cmfd_src[i] > 0) { weightfactors[i] = cmfd_src[i] * norm / sourcecounts[i]; @@ -561,7 +553,7 @@ void free_memory_cmfd() cmfd::indices.clear(); cmfd::egrid.clear(); - // Resize xtensors to be empty + // Resize tensors to be empty cmfd::indexmap.resize({0}); // Set pointers to null diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index b1bfde03d13..ec9da5b8abc 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -254,7 +254,7 @@ void read_ce_cross_sections(const vector>& nuc_temps, if (settings::photon_transport && settings::electron_treatment == ElectronTreatment::TTB) { // Take logarithm of energies since they are log-log interpolated - data::ttb_e_grid = xt::log(data::ttb_e_grid); + data::ttb_e_grid = tensor::log(data::ttb_e_grid); } // Show minimum/maximum temperature diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp index 50f1aca112b..3433f269e35 100644 --- a/src/distribution_angle.cpp +++ b/src/distribution_angle.cpp @@ -2,8 +2,7 @@ #include // for abs, copysign -#include "xtensor/xarray.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include "openmc/endf.h" #include "openmc/hdf5_interface.h" @@ -30,7 +29,7 @@ AngleDistribution::AngleDistribution(hid_t group) hid_t dset = open_dataset(group, "mu"); read_attribute(dset, "offsets", offsets); read_attribute(dset, "interpolation", interp); - xt::xarray temp; + tensor::Tensor temp; read_dataset(dset, temp); close_dataset(dset); @@ -41,13 +40,13 @@ AngleDistribution::AngleDistribution(hid_t group) if (i < n_energy - 1) { n = offsets[i + 1] - j; } else { - n = temp.shape()[1] - j; + n = temp.shape(1) - j; } // Create and initialize tabular distribution - auto xs = xt::view(temp, 0, xt::range(j, j + n)); - auto ps = xt::view(temp, 1, xt::range(j, j + n)); - auto cs = xt::view(temp, 2, xt::range(j, j + n)); + tensor::View xs = temp.slice(0, tensor::range(j, j + n)); + tensor::View ps = temp.slice(1, tensor::range(j, j + n)); + tensor::View cs = temp.slice(2, tensor::range(j, j + n)); vector x {xs.begin(), xs.end()}; vector p {ps.begin(), ps.end()}; vector c {cs.begin(), cs.end()}; diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index a4a5ce9e1b6..2f8e6cf1a99 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -4,7 +4,7 @@ #include // for size_t #include // for back_inserter -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include "openmc/endf.h" #include "openmc/hdf5_interface.h" @@ -60,11 +60,11 @@ ContinuousTabular::ContinuousTabular(hid_t group) hid_t dset = open_dataset(group, "energy"); // Get interpolation parameters - xt::xarray temp; + tensor::Tensor temp; read_attribute(dset, "interpolation", temp); - auto temp_b = xt::view(temp, 0); // view of breakpoints - auto temp_i = xt::view(temp, 1); // view of interpolation parameters + tensor::View temp_b = temp.slice(0); // breakpoints + tensor::View temp_i = temp.slice(1); // interpolation parameters std::copy(temp_b.begin(), temp_b.end(), std::back_inserter(breakpoints_)); for (const auto i : temp_i) @@ -85,7 +85,7 @@ ContinuousTabular::ContinuousTabular(hid_t group) read_attribute(dset, "interpolation", interp); read_attribute(dset, "n_discrete_lines", n_discrete); - xt::xarray eout; + tensor::Tensor eout; read_dataset(dset, eout); close_dataset(dset); @@ -96,7 +96,7 @@ ContinuousTabular::ContinuousTabular(hid_t group) if (i < n_energy - 1) { n = offsets[i + 1] - j; } else { - n = eout.shape()[1] - j; + n = eout.shape(1) - j; } // Assign interpolation scheme and number of discrete lines @@ -105,15 +105,15 @@ ContinuousTabular::ContinuousTabular(hid_t group) d.n_discrete = n_discrete[i]; // Copy data - d.e_out = xt::view(eout, 0, xt::range(j, j + n)); - d.p = xt::view(eout, 1, xt::range(j, j + n)); + d.e_out = eout.slice(0, tensor::range(j, j + n)); + d.p = eout.slice(1, tensor::range(j, j + n)); // To get answers that match ACE data, for now we still use the tabulated // CDF values that were passed through to the HDF5 library. At a later // time, we can remove the CDF values from the HDF5 library and // reconstruct them using the PDF if (true) { - d.c = xt::view(eout, 2, xt::range(j, j + n)); + d.c = eout.slice(2, tensor::range(j, j + n)); } else { // Calculate cumulative distribution function -- discrete portion for (int k = 0; k < d.n_discrete; ++k) { diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index dc4fbb420b8..bc8dcc9ea77 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -1,9 +1,6 @@ #include "openmc/eigenvalue.h" -#include "xtensor/xbuilder.hpp" -#include "xtensor/xmath.hpp" -#include "xtensor/xtensor.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include "openmc/array.h" #include "openmc/bank.h" @@ -39,7 +36,7 @@ namespace simulation { double keff_generation; array k_sum; vector entropy; -xt::xtensor source_frac; +tensor::Tensor source_frac; } // namespace simulation @@ -452,7 +449,7 @@ int openmc_get_keff(double* k_combined) const auto& gt = simulation::global_tallies; array kv {}; - xt::xtensor cov = xt::zeros({3, 3}); + tensor::Tensor cov = tensor::zeros({3, 3}); kv[0] = gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n; kv[1] = gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n; kv[2] = gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n; @@ -591,7 +588,7 @@ void shannon_entropy() { // Get source weight in each mesh bin bool sites_outside; - xt::xtensor p = + tensor::Tensor p = simulation::entropy_mesh->count_sites(simulation::fission_bank.data(), simulation::fission_bank.size(), &sites_outside); @@ -603,7 +600,7 @@ void shannon_entropy() if (mpi::master) { // Normalize to total weight of bank sites - p /= xt::sum(p); + p /= p.sum(); // Sum values to obtain Shannon entropy double H = 0.0; @@ -627,7 +624,7 @@ void ufs_count_sites() std::size_t n = simulation::ufs_mesh->n_bins(); double vol_frac = simulation::ufs_mesh->volume_frac_; - simulation::source_frac = xt::xtensor({n}, vol_frac); + simulation::source_frac = tensor::Tensor({n}, vol_frac); } else { // count number of source sites in each ufs mesh cell @@ -649,7 +646,7 @@ void ufs_count_sites() #endif // Normalize to total weight to get fraction of source in each cell - double total = xt::sum(simulation::source_frac)(); + double total = simulation::source_frac.sum(); simulation::source_frac /= total; // Since the total starting weight is not equal to n_particles, we need to diff --git a/src/endf.cpp b/src/endf.cpp index c0c1d2e7e8b..b298ebe01f9 100644 --- a/src/endf.cpp +++ b/src/endf.cpp @@ -5,8 +5,7 @@ #include // for back_inserter #include // for runtime_error -#include "xtensor/xarray.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include "openmc/array.h" #include "openmc/constants.h" @@ -153,11 +152,11 @@ Tabulated1D::Tabulated1D(hid_t dset) for (const auto i : int_temp) int_.push_back(int2interp(i)); - xt::xarray arr; + tensor::Tensor arr; read_dataset(dset, arr); - auto xs = xt::view(arr, 0); - auto ys = xt::view(arr, 1); + tensor::View xs = arr.slice(0); + tensor::View ys = arr.slice(1); std::copy(xs.begin(), xs.end(), std::back_inserter(x_)); std::copy(ys.begin(), ys.end(), std::back_inserter(y_)); @@ -229,12 +228,12 @@ double Tabulated1D::operator()(double x) const CoherentElasticXS::CoherentElasticXS(hid_t dset) { // Read 2D array from dataset - xt::xarray arr; + tensor::Tensor arr; read_dataset(dset, arr); // Get views for Bragg edges and structure factors - auto E = xt::view(arr, 0); - auto s = xt::view(arr, 1); + tensor::View E = arr.slice(0); + tensor::View s = arr.slice(1); // Copy Bragg edges and partial sums of structure factors std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_)); diff --git a/src/finalize.cpp b/src/finalize.cpp index 5ca9d574627..4ac6d09f321 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -30,7 +30,7 @@ #include "openmc/volume_calc.h" #include "openmc/weight_windows.h" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" namespace openmc { @@ -203,7 +203,7 @@ int openmc_reset() // Reset global tallies simulation::n_realizations = 0; - xt::view(simulation::global_tallies, xt::all()) = 0.0; + simulation::global_tallies.fill(0.0); simulation::k_col_abs = 0.0; simulation::k_col_tra = 0.0; diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index c56d485e281..00c6a4399c0 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -4,8 +4,7 @@ #include #include -#include "xtensor/xarray.hpp" -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include #include "hdf5.h" @@ -466,22 +465,19 @@ void read_dataset_lowlevel(hid_t obj_id, const char* name, hid_t mem_type_id, } template<> -void read_dataset(hid_t dset, xt::xarray>& arr, bool indep) +void read_dataset( + hid_t dset, tensor::Tensor>& tensor, bool indep) { // Get shape of dataset vector shape = object_shape(dset); - // Allocate new array to read data into - std::size_t size = 1; - for (const auto x : shape) - size *= x; - vector> buffer(size); + // Resize tensor and read data directly + vector tshape(shape.begin(), shape.end()); + tensor.resize(tshape); - // Read data from attribute - read_complex(dset, nullptr, buffer.data(), indep); - - // Adapt into xarray - arr = xt::adapt(buffer, shape); + // Read data from dataset + read_complex(dset, nullptr, + reinterpret_cast*>(tensor.data()), indep); } void read_double(hid_t obj_id, const char* name, double* buffer, bool indep) diff --git a/src/material.cpp b/src/material.cpp index 072e6decad1..21b11b8b9ce 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -8,9 +8,7 @@ #include #include -#include "xtensor/xbuilder.hpp" -#include "xtensor/xoperation.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include "openmc/capi.h" #include "openmc/container_util.h" @@ -216,7 +214,7 @@ Material::Material(pugi::xml_node node) // allocate arrays in Material object auto n = names.size(); nuclide_.reserve(n); - atom_density_ = xt::empty({n}); + atom_density_ = tensor::Tensor({n}); if (settings::photon_transport) element_.reserve(n); @@ -290,14 +288,14 @@ Material::Material(pugi::xml_node node) // Check to make sure either all atom percents or all weight percents are // given - if (!(xt::all(atom_density_ >= 0.0) || xt::all(atom_density_ <= 0.0))) { + if (!((atom_density_ >= 0.0).all() || (atom_density_ <= 0.0).all())) { fatal_error( "Cannot mix atom and weight percents in material " + std::to_string(id_)); } // Determine density if it is a sum value if (sum_density) - density_ = xt::sum(atom_density_)(); + density_ = atom_density_.sum(); if (check_for_node(node, "temperature")) { temperature_ = std::stod(get_node_value(node, "temperature")); @@ -435,7 +433,7 @@ void Material::normalize_density() // determine normalized atom percents. if given atom percents, this is // straightforward. if given weight percents, the value is w/awr and is // divided by sum(w/awr) - atom_density_ /= xt::sum(atom_density_)(); + atom_density_ /= atom_density_.sum(); // Change density in g/cm^3 to atom/b-cm. Since all values are now in // atom percent, the sum needs to be re-evaluated as 1/sum(x*awr) @@ -641,14 +639,14 @@ void Material::init_bremsstrahlung() bool positron = (particle == 1); // Allocate arrays for TTB data - ttb->pdf = xt::zeros({n_e, n_e}); - ttb->cdf = xt::zeros({n_e, n_e}); - ttb->yield = xt::zeros({n_e}); + ttb->pdf = tensor::zeros({n_e, n_e}); + ttb->cdf = tensor::zeros({n_e, n_e}); + ttb->yield = tensor::zeros({n_e}); // Allocate temporary arrays - xt::xtensor stopping_power_collision({n_e}, 0.0); - xt::xtensor stopping_power_radiative({n_e}, 0.0); - xt::xtensor dcs({n_e, n_k}, 0.0); + auto stopping_power_collision = tensor::zeros({n_e}); + auto stopping_power_radiative = tensor::zeros({n_e}); + auto dcs = tensor::zeros({n_e, n_k}); double Z_eq_sq = 0.0; double sum_density = 0.0; @@ -698,18 +696,18 @@ void Material::init_bremsstrahlung() 1.0595e-3 * std::pow(t, 5) + 7.0568e-5 * std::pow(t, 6) - 1.808e-6 * std::pow(t, 7)); stopping_power_radiative(i) *= r; - auto dcs_i = xt::view(dcs, i, xt::all()); + tensor::View dcs_i = dcs.slice(i); dcs_i *= r; } } // Total material stopping power - xt::xtensor stopping_power = + tensor::Tensor stopping_power = stopping_power_collision + stopping_power_radiative; // Loop over photon energies - xt::xtensor f({n_e}, 0.0); - xt::xtensor z({n_e}, 0.0); + auto f = tensor::zeros({n_e}); + auto z = tensor::zeros({n_e}); for (int i = 0; i < n_e - 1; ++i) { double w = data::ttb_e_grid(i); @@ -797,7 +795,8 @@ void Material::init_bremsstrahlung() } // Use logarithm of number yield since it is log-log interpolated - ttb->yield = xt::where(ttb->yield > 0.0, xt::log(ttb->yield), -500.0); + ttb->yield = + tensor::where(ttb->yield > 0.0, tensor::log(ttb->yield), -500.0); } } @@ -979,7 +978,7 @@ void Material::set_density(double density, const std::string& units) density_ = density; // Determine normalized atom percents - double sum_percent = xt::sum(atom_density_)(); + double sum_percent = atom_density_.sum(); atom_density_ /= sum_percent; // Recalculate nuclide atom densities based on given density @@ -1020,7 +1019,7 @@ void Material::set_densities( if (n != nuclide_.size()) { nuclide_.resize(n); - atom_density_ = xt::zeros({n}); + atom_density_ = tensor::zeros({n}); if (settings::photon_transport) element_.resize(n); } @@ -1181,8 +1180,8 @@ void Material::add_nuclide(const std::string& name, double density) auto n = nuclide_.size(); // Create copy of atom_density_ array with one extra entry - xt::xtensor atom_density = xt::zeros({n}); - xt::view(atom_density, xt::range(0, n - 1)) = atom_density_; + tensor::Tensor atom_density = tensor::zeros({n}); + atom_density.slice(tensor::range(0, n - 1)) = atom_density_; atom_density(n - 1) = density; atom_density_ = atom_density; diff --git a/src/mesh.cpp b/src/mesh.cpp index 789113ccc31..5ab7ac3988b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -6,6 +6,7 @@ #define _USE_MATH_DEFINES // to make M_PI declared in Intel and MSVC compilers #include // for ceil #include // for size_t +#include // for accumulate #include #ifdef _MSC_VER @@ -16,13 +17,7 @@ #include "mpi.h" #endif -#include "xtensor/xadapt.hpp" -#include "xtensor/xbuilder.hpp" -#include "xtensor/xeval.hpp" -#include "xtensor/xmath.hpp" -#include "xtensor/xsort.hpp" -#include "xtensor/xtensor.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include // for fmt #include "openmc/capi.h" @@ -772,11 +767,9 @@ std::string StructuredMesh::bin_label(int bin) const } } -xt::xtensor StructuredMesh::get_x_shape() const +tensor::Tensor StructuredMesh::get_shape_tensor() const { - // because method is const, shape_ is const as well and can't be adapted - auto tmp_shape = shape_; - return xt::adapt(tmp_shape, {n_dimension_}); + return tensor::Tensor(shape_.data(), static_cast(n_dimension_)); } Position StructuredMesh::sample_element( @@ -961,10 +954,11 @@ void UnstructuredMesh::to_hdf5_inner(hid_t mesh_group) const write_dataset(mesh_group, "length_multiplier", length_multiplier_); // write vertex coordinates - xt::xtensor vertices({static_cast(this->n_vertices()), 3}); + tensor::Tensor vertices( + {static_cast(this->n_vertices()), static_cast(3)}); for (int i = 0; i < this->n_vertices(); i++) { auto v = this->vertex(i); - xt::view(vertices, i, xt::all()) = xt::xarray({v.x, v.y, v.z}); + vertices.slice(i) = {v.x, v.y, v.z}; } write_dataset(mesh_group, "vertices", vertices); @@ -972,8 +966,10 @@ void UnstructuredMesh::to_hdf5_inner(hid_t mesh_group) const // write element types and connectivity vector volumes; - xt::xtensor connectivity({static_cast(this->n_bins()), 8}); - xt::xtensor elem_types({static_cast(this->n_bins()), 1}); + tensor::Tensor connectivity( + {static_cast(this->n_bins()), static_cast(8)}); + tensor::Tensor elem_types( + {static_cast(this->n_bins()), static_cast(1)}); for (int i = 0; i < this->n_bins(); i++) { auto conn = this->connectivity(i); @@ -981,21 +977,18 @@ void UnstructuredMesh::to_hdf5_inner(hid_t mesh_group) const // write linear tet element if (conn.size() == 4) { - xt::view(elem_types, i, xt::all()) = - static_cast(ElementType::LINEAR_TET); - xt::view(connectivity, i, xt::all()) = - xt::xarray({conn[0], conn[1], conn[2], conn[3], -1, -1, -1, -1}); + elem_types.slice(i) = static_cast(ElementType::LINEAR_TET); + connectivity.slice(i) = { + conn[0], conn[1], conn[2], conn[3], -1, -1, -1, -1}; // write linear hex element } else if (conn.size() == 8) { - xt::view(elem_types, i, xt::all()) = - static_cast(ElementType::LINEAR_HEX); - xt::view(connectivity, i, xt::all()) = xt::xarray({conn[0], conn[1], - conn[2], conn[3], conn[4], conn[5], conn[6], conn[7]}); + elem_types.slice(i) = static_cast(ElementType::LINEAR_HEX); + connectivity.slice(i) = { + conn[0], conn[1], conn[2], conn[3], conn[4], conn[5], conn[6], conn[7]}; } else { num_elem_skipped++; - xt::view(elem_types, i, xt::all()) = - static_cast(ElementType::UNSUPPORTED); - xt::view(connectivity, i, xt::all()) = -1; + elem_types.slice(i) = static_cast(ElementType::UNSUPPORTED); + connectivity.slice(i) = -1; } } @@ -1096,7 +1089,7 @@ int StructuredMesh::n_surface_bins() const return 4 * n_dimension_ * n_bins(); } -xt::xtensor StructuredMesh::count_sites( +tensor::Tensor StructuredMesh::count_sites( const SourceSite* bank, int64_t length, bool* outside) const { // Determine shape of array for counts @@ -1104,7 +1097,7 @@ xt::xtensor StructuredMesh::count_sites( vector shape = {m}; // Create array of zeros - xt::xarray cnt {shape, 0.0}; + auto cnt = tensor::zeros(shape); bool outside_ = false; for (int64_t i = 0; i < length; i++) { @@ -1123,31 +1116,25 @@ xt::xtensor StructuredMesh::count_sites( cnt(mesh_bin) += site.wgt; } - // Create copy of count data. Since ownership will be acquired by xtensor, - // std::allocator must be used to avoid Valgrind mismatched free() / delete - // warnings. + // Create reduced count data + auto counts = tensor::zeros(shape); int total = cnt.size(); - double* cnt_reduced = std::allocator {}.allocate(total); #ifdef OPENMC_MPI // collect values from all processors MPI_Reduce( - cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); + cnt.data(), counts.data(), total, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); // Check if there were sites outside the mesh for any processor if (outside) { MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm); } #else - std::copy(cnt.data(), cnt.data() + total, cnt_reduced); + std::copy(cnt.data(), cnt.data() + total, counts.data()); if (outside) *outside = outside_; #endif - // Adapt reduced values in array back into an xarray - auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), shape); - xt::xarray counts = arr; - return counts; } @@ -1340,10 +1327,10 @@ void StructuredMesh::surface_bins_crossed( int RegularMesh::set_grid() { - auto shape = xt::adapt(shape_, {n_dimension_}); + tensor::Tensor shape(shape_.data(), static_cast(n_dimension_)); // Check that dimensions are all greater than zero - if (xt::any(shape <= 0)) { + if ((shape <= 0).any()) { set_errmsg("All entries for a regular mesh dimensions " "must be positive."); return OPENMC_E_INVALID_ARGUMENT; @@ -1365,13 +1352,13 @@ int RegularMesh::set_grid() } // Check for negative widths - if (xt::any(width_ < 0.0)) { + if ((width_ < 0.0).any()) { set_errmsg("Cannot have a negative width on a regular mesh."); return OPENMC_E_INVALID_ARGUMENT; } // Set width and upper right coordinate - upper_right_ = xt::eval(lower_left_ + shape * width_); + upper_right_ = lower_left_ + shape * width_; } else if (upper_right_.size() > 0) { @@ -1383,7 +1370,7 @@ int RegularMesh::set_grid() } // Check that upper-right is above lower-left - if (xt::any(upper_right_ < lower_left_)) { + if ((upper_right_ < lower_left_).any()) { set_errmsg( "The upper_right coordinates of a regular mesh must be greater than " "the lower_left coordinates."); @@ -1391,11 +1378,11 @@ int RegularMesh::set_grid() } // Set width - width_ = xt::eval((upper_right_ - lower_left_) / shape); + width_ = (upper_right_ - lower_left_) / shape; } // Set material volumes - volume_frac_ = 1.0 / xt::prod(shape)(); + volume_frac_ = 1.0 / shape.prod(); element_volume_ = 1.0; for (int i = 0; i < n_dimension_; i++) { @@ -1411,7 +1398,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node} fatal_error("Must specify on a regular mesh."); } - xt::xtensor shape = get_node_xarray(node, "dimension"); + tensor::Tensor shape = get_node_tensor(node, "dimension"); int n = n_dimension_ = shape.size(); if (n != 1 && n != 2 && n != 3) { fatal_error("Mesh must be one, two, or three dimensions."); @@ -1421,7 +1408,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node} // Check for lower-left coordinates if (check_for_node(node, "lower_left")) { // Read mesh lower-left corner location - lower_left_ = get_node_xarray(node, "lower_left"); + lower_left_ = get_node_tensor(node, "lower_left"); } else { fatal_error("Must specify on a mesh."); } @@ -1432,11 +1419,11 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node} fatal_error("Cannot specify both and on a mesh."); } - width_ = get_node_xarray(node, "width"); + width_ = get_node_tensor(node, "width"); } else if (check_for_node(node, "upper_right")) { - upper_right_ = get_node_xarray(node, "upper_right"); + upper_right_ = get_node_tensor(node, "upper_right"); } else { fatal_error("Must specify either or on a mesh."); @@ -1454,7 +1441,7 @@ RegularMesh::RegularMesh(hid_t group) : StructuredMesh {group} fatal_error("Must specify on a regular mesh."); } - xt::xtensor shape; + tensor::Tensor shape; read_dataset(group, "dimension", shape); int n = n_dimension_ = shape.size(); if (n != 1 && n != 2 && n != 3) { @@ -1569,13 +1556,13 @@ std::pair, vector> RegularMesh::plot( void RegularMesh::to_hdf5_inner(hid_t mesh_group) const { - write_dataset(mesh_group, "dimension", get_x_shape()); + write_dataset(mesh_group, "dimension", get_shape_tensor()); write_dataset(mesh_group, "lower_left", lower_left_); write_dataset(mesh_group, "upper_right", upper_right_); write_dataset(mesh_group, "width", width_); } -xt::xtensor RegularMesh::count_sites( +tensor::Tensor RegularMesh::count_sites( const SourceSite* bank, int64_t length, bool* outside) const { // Determine shape of array for counts @@ -1583,7 +1570,7 @@ xt::xtensor RegularMesh::count_sites( vector shape = {m}; // Create array of zeros - xt::xarray cnt {shape, 0.0}; + auto cnt = tensor::zeros(shape); bool outside_ = false; for (int64_t i = 0; i < length; i++) { @@ -1602,31 +1589,25 @@ xt::xtensor RegularMesh::count_sites( cnt(mesh_bin) += site.wgt; } - // Create copy of count data. Since ownership will be acquired by xtensor, - // std::allocator must be used to avoid Valgrind mismatched free() / delete - // warnings. + // Create reduced count data + auto counts = tensor::zeros(shape); int total = cnt.size(); - double* cnt_reduced = std::allocator {}.allocate(total); #ifdef OPENMC_MPI // collect values from all processors MPI_Reduce( - cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); + cnt.data(), counts.data(), total, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); // Check if there were sites outside the mesh for any processor if (outside) { MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm); } #else - std::copy(cnt.data(), cnt.data() + total, cnt_reduced); + std::copy(cnt.data(), cnt.data() + total, counts.data()); if (outside) *outside = outside_; #endif - // Adapt reduced values in array back into an xarray - auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), shape); - xt::xarray counts = arr; - return counts; } @@ -2698,7 +2679,7 @@ extern "C" int openmc_regular_mesh_get_params( return err; RegularMesh* m = dynamic_cast(model::meshes[index].get()); - if (m->lower_left_.dimension() == 0) { + if (m->lower_left_.empty()) { set_errmsg("Mesh parameters have not been set."); return OPENMC_E_ALLOCATE; } @@ -2725,17 +2706,17 @@ extern "C" int openmc_regular_mesh_set_params( vector shape = {static_cast(n)}; if (ll && ur) { - m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape); - m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape); - m->width_ = (m->upper_right_ - m->lower_left_) / m->get_x_shape(); + m->lower_left_ = tensor::Tensor(ll, n); + m->upper_right_ = tensor::Tensor(ur, n); + m->width_ = (m->upper_right_ - m->lower_left_) / m->get_shape_tensor(); } else if (ll && width) { - m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape); - m->width_ = xt::adapt(width, n, xt::no_ownership(), shape); - m->upper_right_ = m->lower_left_ + m->get_x_shape() * m->width_; + m->lower_left_ = tensor::Tensor(ll, n); + m->width_ = tensor::Tensor(width, n); + m->upper_right_ = m->lower_left_ + m->get_shape_tensor() * m->width_; } else if (ur && width) { - m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape); - m->width_ = xt::adapt(width, n, xt::no_ownership(), shape); - m->lower_left_ = m->upper_right_ - m->get_x_shape() * m->width_; + m->upper_right_ = tensor::Tensor(ur, n); + m->width_ = tensor::Tensor(width, n); + m->lower_left_ = m->upper_right_ - m->get_shape_tensor() * m->width_; } else { set_errmsg("At least two parameters must be specified."); return OPENMC_E_INVALID_ARGUMENT; @@ -2745,7 +2726,7 @@ extern "C" int openmc_regular_mesh_set_params( // TODO: incorporate this into method in RegularMesh that can be called from // here and from constructor - m->volume_frac_ = 1.0 / xt::prod(m->get_x_shape())(); + m->volume_frac_ = 1.0 / m->get_shape_tensor().prod(); m->element_volume_ = 1.0; for (int i = 0; i < m->n_dimension_; i++) { m->element_volume_ *= m->width_[i]; @@ -2794,7 +2775,7 @@ int openmc_structured_mesh_get_grid_impl(int32_t index, double** grid_x, return err; C* m = dynamic_cast(model::meshes[index].get()); - if (m->lower_left_.dimension() == 0) { + if (m->lower_left_.empty()) { set_errmsg("Mesh parameters have not been set."); return OPENMC_E_ALLOCATE; } diff --git a/src/mgxs.cpp b/src/mgxs.cpp index ea78238b809..65a7a0c55ca 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -5,10 +5,7 @@ #include #include -#include "xtensor/xadapt.hpp" -#include "xtensor/xmath.hpp" -#include "xtensor/xsort.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include #include "openmc/error.h" @@ -33,8 +30,7 @@ void Mgxs::init(const std::string& in_name, double in_awr, // Set the metadata name = in_name; awr = in_awr; - // TODO: Remove adapt when in_KTs is an xtensor - kTs = xt::adapt(in_kTs); + kTs = tensor::Tensor(in_kTs.data(), in_kTs.size()); fissionable = in_fissionable; scatter_format = in_scatter_format; xs.resize(in_kTs.size()); @@ -73,7 +69,7 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector& temperature, } get_datasets(kT_group, dset_names); vector shape = {num_temps}; - xt::xarray temps_available(shape); + tensor::Tensor temps_available(shape); for (int i = 0; i < num_temps; i++) { read_double(kT_group, dset_names[i], &temps_available[i], true); @@ -108,19 +104,7 @@ void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector& temperature, // Determine actual temperatures to read for (const auto& T : temperature) { // Determine the closest temperature value - // NOTE: the below block could be replaced with the following line, - // though this gives a runtime error if using LLVM 20 or newer, - // likely due to a bug in xtensor. - // auto i_closest = xt::argmin(xt::abs(temps_available - T))[0]; - double closest = std::numeric_limits::max(); - int i_closest = 0; - for (int i = 0; i < temps_available.size(); i++) { - double diff = std::abs(temps_available[i] - T); - if (diff < closest) { - closest = diff; - i_closest = i; - } - } + auto i_closest = tensor::abs(temps_available - T).argmin(); double temp_actual = temps_available[i_closest]; if (std::fabs(temp_actual - T) < settings::temperature_tolerance) { @@ -355,7 +339,7 @@ Mgxs::Mgxs(const std::string& in_name, const vector& mat_kTs, for (int m = 0; m < micros.size(); m++) { switch (settings::temperature_method) { case TemperatureMethod::NEAREST: { - micro_t[m] = xt::argmin(xt::abs(micros[m]->kTs - temp_desired))[0]; + micro_t[m] = tensor::abs(micros[m]->kTs - temp_desired).argmin(); auto temp_actual = micros[m]->kTs[micro_t[m]]; if (std::abs(temp_actual - temp_desired) >= @@ -368,7 +352,7 @@ Mgxs::Mgxs(const std::string& in_name, const vector& mat_kTs, case TemperatureMethod::INTERPOLATION: // Get a list of bounding temperatures for each actual temperature // present in the model - for (int k = 0; k < micros[m]->kTs.shape()[0] - 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; @@ -474,7 +458,7 @@ double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu, val = xs_t->delayed_nu_fission(a, *dg, gin); } else { val = 0.; - for (int d = 0; d < xs_t->delayed_nu_fission.shape()[1]; d++) { + for (int d = 0; d < xs_t->delayed_nu_fission.shape(1); d++) { val += xs_t->delayed_nu_fission(a, d, gin); } } @@ -489,7 +473,7 @@ double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu, } else { // provide an outgoing group-wise sum val = 0.; - for (int g = 0; g < xs_t->chi_prompt.shape()[2]; g++) { + for (int g = 0; g < xs_t->chi_prompt.shape(2); g++) { val += xs_t->chi_prompt(a, gin, g); } } @@ -508,13 +492,13 @@ double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu, } else { if (dg != nullptr) { val = 0.; - for (int g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) { + 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 g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) { - for (int d = 0; d < xs_t->delayed_nu_fission.shape()[3]; d++) { + 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); } } @@ -650,7 +634,7 @@ bool Mgxs::equiv(const Mgxs& that) int Mgxs::get_temperature_index(double sqrtkT) const { - return xt::argmin(xt::abs(kTs - sqrtkT * sqrtkT))[0]; + return tensor::abs(kTs - sqrtkT * sqrtkT).argmin(); } //============================================================================== diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 69e603a7c66..17d6e952c39 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -17,8 +17,7 @@ #include -#include "xtensor/xbuilder.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include // for sort, min_element #include @@ -361,8 +360,7 @@ void Nuclide::create_derived( { for (const auto& grid : grid_) { // Allocate and initialize cross section - array shape {grid.energy.size(), 5}; - xs_.emplace_back(shape, 0.0); + xs_.push_back(tensor::zeros({grid.energy.size(), 5})); } reaction_index_.fill(C_NONE); @@ -375,9 +373,8 @@ void Nuclide::create_derived( for (int t = 0; t < kTs_.size(); ++t) { int j = rx->xs_[t].threshold; int n = rx->xs_[t].value.size(); - auto xs = xt::adapt(rx->xs_[t].value); - auto pprod = xt::view(xs_[t], xt::range(j, j + n), XS_PHOTON_PROD); - + auto xs = tensor::Tensor( + rx->xs_[t].value.data(), rx->xs_[t].value.size()); for (const auto& p : rx->products_) { if (p.particle_.is_photon()) { for (int k = 0; k < n; ++k) { @@ -396,7 +393,7 @@ void Nuclide::create_derived( } } - pprod[k] += f * xs[k] * (*p.yield_)(E); + xs_[t](j + k, XS_PHOTON_PROD) += f * xs[k] * (*p.yield_)(E); } } } @@ -406,20 +403,17 @@ void Nuclide::create_derived( continue; // Add contribution to total cross section - auto total = xt::view(xs_[t], xt::range(j, j + n), XS_TOTAL); - total += xs; + xs_[t].slice(tensor::range(j, j + n), XS_TOTAL) += xs; // Add contribution to absorption cross section - auto absorption = xt::view(xs_[t], xt::range(j, j + n), XS_ABSORPTION); if (is_disappearance(rx->mt_)) { - absorption += xs; + xs_[t].slice(tensor::range(j, j + n), XS_ABSORPTION) += xs; } if (is_fission(rx->mt_)) { fissionable_ = true; - auto fission = xt::view(xs_[t], xt::range(j, j + n), XS_FISSION); - fission += xs; - absorption += xs; + xs_[t].slice(tensor::range(j, j + n), XS_FISSION) += xs; + xs_[t].slice(tensor::range(j, j + n), XS_ABSORPTION) += xs; // Keep track of fission reactions if (t == 0) { @@ -510,7 +504,7 @@ void Nuclide::init_grid() double spacing = std::log(E_max / E_min) / M; // Create equally log-spaced energy grid - auto umesh = xt::linspace(0.0, M * spacing, M + 1); + auto umesh = tensor::linspace(0.0, M * spacing, M + 1); for (auto& grid : grid_) { // Resize array for storing grid indices diff --git a/src/output.cpp b/src/output.cpp index ae2daaffc14..8f9802bb9bc 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -17,7 +17,7 @@ #ifdef _OPENMP #include #endif -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include "openmc/capi.h" #include "openmc/cell.h" diff --git a/src/photon.cpp b/src/photon.cpp index 951acb9fbd8..5e2ba688427 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -13,11 +13,7 @@ #include "openmc/search.h" #include "openmc/settings.h" -#include "xtensor/xbuilder.hpp" -#include "xtensor/xmath.hpp" -#include "xtensor/xoperation.hpp" -#include "xtensor/xslice.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include #include @@ -33,7 +29,7 @@ constexpr int PhotonInteraction::MAX_STACK_SIZE; namespace data { -xt::xtensor compton_profile_pz; +tensor::Tensor compton_profile_pz; std::unordered_map element_map; vector> elements; @@ -46,8 +42,6 @@ vector> elements; PhotonInteraction::PhotonInteraction(hid_t group) { - using namespace xt::placeholders; - // Set index of element in global vector index_ = data::elements.size(); @@ -96,7 +90,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) read_dataset(rgroup, "xs", pair_production_electron_); close_group(rgroup); } else { - pair_production_electron_ = xt::zeros_like(energy_); + pair_production_electron_ = tensor::zeros_like(energy_); } // Read pair production @@ -105,7 +99,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) read_dataset(rgroup, "xs", pair_production_nuclear_); close_group(rgroup); } else { - pair_production_nuclear_ = xt::zeros_like(energy_); + pair_production_nuclear_ = tensor::zeros_like(energy_); } // Read photoelectric @@ -119,7 +113,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) read_dataset(rgroup, "xs", heating_); close_group(rgroup); } else { - heating_ = xt::zeros_like(energy_); + heating_ = tensor::zeros_like(energy_); } // Read subshell photoionization cross section and atomic relaxation data @@ -133,7 +127,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) } shells_.resize(n_shell); - cross_sections_ = xt::zeros({energy_.size(), n_shell}); + cross_sections_ = tensor::zeros({energy_.size(), n_shell}); // Create mapping from designator to index std::unordered_map shell_map; @@ -168,15 +162,17 @@ PhotonInteraction::PhotonInteraction(hid_t group) } // Read subshell cross section - xt::xtensor xs; + tensor::Tensor xs; dset = open_dataset(tgroup, "xs"); read_attribute(dset, "threshold_idx", shell.threshold); close_dataset(dset); read_dataset(tgroup, "xs", xs); auto cross_section = - xt::view(cross_sections_, xt::range(shell.threshold, _), i); - cross_section = xt::where(xs > 0, xt::log(xs), 0); + cross_sections_.slice(tensor::range(static_cast(shell.threshold), + cross_sections_.shape(0)), + i); + cross_section = tensor::where(xs > 0, tensor::log(xs), 0); if (object_exists(tgroup, "transitions")) { // Determine dimensions of transitions @@ -186,11 +182,12 @@ PhotonInteraction::PhotonInteraction(hid_t group) int n_transition = dims[0]; if (n_transition > 0) { - xt::xtensor matrix; + tensor::Tensor matrix; read_dataset(tgroup, "transitions", matrix); // Transition probability normalization - double norm = xt::sum(xt::col(matrix, 3))(); + double norm = + tensor::Tensor(matrix.slice(tensor::all, 3)).sum(); shell.transitions.resize(n_transition); for (int j = 0; j < n_transition; ++j) { @@ -220,7 +217,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) // Read electron shell PDF and binding energies read_dataset(rgroup, "num_electrons", electron_pdf_); - electron_pdf_ /= xt::sum(electron_pdf_); + electron_pdf_ /= electron_pdf_.sum(); read_dataset(rgroup, "binding_energy", binding_energy_); // Read Compton profiles @@ -238,7 +235,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) auto is_close = [](double a, double b) { return std::abs(a - b) / a < FP_REL_PRECISION; }; - subshell_map_ = xt::full_like(binding_energy_, -1); + subshell_map_ = tensor::Tensor(binding_energy_.shape(), -1); for (int i = 0; i < binding_energy_.size(); ++i) { double E_b = binding_energy_[i]; if (i < n_shell && is_close(E_b, shells_[i].binding_energy)) { @@ -257,7 +254,7 @@ PhotonInteraction::PhotonInteraction(hid_t group) // Create Compton profile CDF auto n_profile = data::compton_profile_pz.size(); auto n_shell_compton = profile_pdf_.shape(0); - profile_cdf_ = xt::empty({n_shell_compton, n_profile}); + profile_cdf_ = tensor::Tensor({n_shell_compton, n_profile}); for (int i = 0; i < n_shell_compton; ++i) { double c = 0.0; profile_cdf_(i, 0) = 0.0; @@ -276,11 +273,11 @@ PhotonInteraction::PhotonInteraction(hid_t group) // Read bremsstrahlung scaled DCS rgroup = open_group(group, "bremsstrahlung"); read_dataset(rgroup, "dcs", dcs_); - auto n_e = dcs_.shape()[0]; - auto n_k = dcs_.shape()[1]; + auto n_e = dcs_.shape(0); + auto n_k = dcs_.shape(1); // Get energy grids used for bremsstrahlung DCS and for stopping powers - xt::xtensor electron_energy; + tensor::Tensor electron_energy; read_dataset(rgroup, "electron_energy", electron_energy); if (data::ttb_k_grid.size() == 0) { read_dataset(rgroup, "photon_energy", data::ttb_k_grid); @@ -305,12 +302,12 @@ PhotonInteraction::PhotonInteraction(hid_t group) (std::log(E(i_grid + 1)) - std::log(E(i_grid))); // Interpolate bremsstrahlung DCS at the cutoff energy and truncate - xt::xtensor dcs({n_e - i_grid, n_k}); + tensor::Tensor dcs({n_e - i_grid, n_k}); for (int i = 0; i < n_k; ++i) { double y = std::exp( std::log(dcs_(i_grid, i)) + f * (std::log(dcs_(i_grid + 1, i)) - std::log(dcs_(i_grid, i)))); - auto col_i = xt::view(dcs, xt::all(), i); + tensor::View col_i = dcs.slice(tensor::all, i); col_i(0) = y; for (int j = i_grid + 1; j < n_e; ++j) { col_i(j - i_grid) = dcs_(j, i); @@ -318,9 +315,11 @@ PhotonInteraction::PhotonInteraction(hid_t group) } dcs_ = dcs; - xt::xtensor frst {cutoff}; - electron_energy = xt::concatenate(xt::xtuple( - frst, xt::view(electron_energy, xt::range(i_grid + 1, n_e)))); + tensor::Tensor frst({static_cast(1)}); + frst(0) = cutoff; + tensor::Tensor rest(electron_energy.slice( + tensor::range(i_grid + 1, electron_energy.size()))); + electron_energy = tensor::concatenate(frst, rest); } // Set incident particle energy grid @@ -329,7 +328,8 @@ PhotonInteraction::PhotonInteraction(hid_t group) } // Calculate the radiative stopping power - stopping_power_radiative_ = xt::empty({data::ttb_e_grid.size()}); + stopping_power_radiative_ = + tensor::Tensor({data::ttb_e_grid.size()}); for (int i = 0; i < data::ttb_e_grid.size(); ++i) { // Integrate over reduced photon energy double c = 0.0; @@ -354,14 +354,15 @@ PhotonInteraction::PhotonInteraction(hid_t group) // values below exp(-499) we store the log as -900, for which exp(-900) // evaluates to zero. double limit = std::exp(-499.0); - energy_ = xt::log(energy_); - coherent_ = xt::where(coherent_ > limit, xt::log(coherent_), -900.0); - incoherent_ = xt::where(incoherent_ > limit, xt::log(incoherent_), -900.0); - photoelectric_total_ = xt::where( - photoelectric_total_ > limit, xt::log(photoelectric_total_), -900.0); - pair_production_total_ = xt::where( - pair_production_total_ > limit, xt::log(pair_production_total_), -900.0); - heating_ = xt::where(heating_ > limit, xt::log(heating_), -900.0); + energy_ = tensor::log(energy_); + coherent_ = tensor::where(coherent_ > limit, tensor::log(coherent_), -900.0); + incoherent_ = + tensor::where(incoherent_ > limit, tensor::log(incoherent_), -900.0); + photoelectric_total_ = tensor::where( + photoelectric_total_ > limit, tensor::log(photoelectric_total_), -900.0); + pair_production_total_ = tensor::where(pair_production_total_ > limit, + tensor::log(pair_production_total_), -900.0); + heating_ = tensor::where(heating_ > limit, tensor::log(heating_), -900.0); } PhotonInteraction::~PhotonInteraction() @@ -512,7 +513,7 @@ void PhotonInteraction::compton_doppler( c = prn(seed) * c_max; // Determine pz corresponding to sampled cdf value - auto cdf_shell = xt::view(profile_cdf_, shell, xt::all()); + tensor::View cdf_shell = profile_cdf_.slice(shell); int i = lower_bound_index(cdf_shell.cbegin(), cdf_shell.cend(), c); double pz_l = data::compton_profile_pz(i); double pz_r = data::compton_profile_pz(i + 1); @@ -608,8 +609,8 @@ void PhotonInteraction::calculate_xs(Particle& p) const // Calculate microscopic photoelectric cross section xs.photoelectric = 0.0; - const auto& xs_lower = xt::row(cross_sections_, i_grid); - const auto& xs_upper = xt::row(cross_sections_, i_grid + 1); + tensor::View xs_lower = cross_sections_.slice(i_grid); + tensor::View xs_upper = cross_sections_.slice(i_grid + 1); for (int i = 0; i < xs_upper.size(); ++i) if (xs_lower(i) != 0) diff --git a/src/physics.cpp b/src/physics.cpp index 1f72e4fac07..c9737205ced 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -30,9 +30,9 @@ #include +#include "openmc/tensor.h" #include // for max, min, max_element #include // for sqrt, exp, log, abs, copysign -#include namespace openmc { @@ -375,8 +375,9 @@ void sample_photon_reaction(Particle& p) // cross sections int i_grid = micro.index_grid; double f = micro.interp_factor; - const auto& xs_lower = xt::row(element.cross_sections_, i_grid); - const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1); + tensor::View xs_lower = element.cross_sections_.slice(i_grid); + tensor::View xs_upper = + element.cross_sections_.slice(i_grid + 1); for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) { const auto& shell {element.shells_[i_shell]}; diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 58a21a67895..32ae10a6092 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -2,7 +2,7 @@ #include -#include "xtensor/xarray.hpp" +#include "openmc/tensor.h" #include #include "openmc/bank.h" diff --git a/src/plot.cpp b/src/plot.cpp index e3b1e84c80e..ea0555202d5 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -7,8 +7,7 @@ #include #include -#include "xtensor/xmanipulation.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include #include #ifdef USE_LIBPNG @@ -74,7 +73,8 @@ void IdData::set_value(size_t y, size_t x, const GeometryState& p, int level) void IdData::set_overlap(size_t y, size_t x) { - xt::view(data_, y, x, xt::all()) = OVERLAP; + for (size_t k = 0; k < data_.shape(2); ++k) + data_(y, x, k) = OVERLAP; } PropertyData::PropertyData(size_t h_res, size_t v_res) @@ -783,14 +783,14 @@ void output_ppm(const std::string& filename, const ImageData& data) // Write header of << "P6\n"; - of << data.shape()[0] << " " << data.shape()[1] << "\n"; + of << data.shape(0) << " " << data.shape(1) << "\n"; of << "255\n"; of.close(); of.open(fname, std::ios::binary | std::ios::app); // Write color for each pixel - for (int y = 0; y < data.shape()[1]; y++) { - for (int x = 0; x < data.shape()[0]; x++) { + for (int y = 0; y < data.shape(1); y++) { + for (int x = 0; x < data.shape(0); x++) { RGBColor rgb = data(x, y); of << rgb.red << rgb.green << rgb.blue; } @@ -822,8 +822,8 @@ void output_png(const std::string& filename, const ImageData& data) png_init_io(png_ptr, fp); // Write header (8 bit colour depth) - int width = data.shape()[0]; - int height = data.shape()[1]; + int width = data.shape(0); + int height = data.shape(1); png_set_IHDR(png_ptr, info_ptr, width, height, 8, PNG_COLOR_TYPE_RGB, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_BASE); png_write_info(png_ptr, info_ptr); @@ -1024,9 +1024,14 @@ void Plot::create_voxel() const // select only cell/material ID data and flip the y-axis int idx = color_by_ == PlotColorBy::cells ? 0 : 2; - xt::xtensor data_slice = - xt::view(ids.data_, xt::all(), xt::all(), idx); - xt::xtensor data_flipped = xt::flip(data_slice, 0); + // Extract 2D slice at index idx from 3D data + size_t rows = ids.data_.shape(0); + size_t cols = ids.data_.shape(1); + tensor::Tensor data_slice({rows, cols}); + for (size_t r = 0; r < rows; ++r) + for (size_t c = 0; c < cols; ++c) + data_slice(r, c) = ids.data_(r, c, idx); + tensor::Tensor data_flipped = data_slice.flip(0); // Write to HDF5 dataset voxel_write_slice(z, dspace, dset, memspace, data_flipped.data()); @@ -1272,7 +1277,8 @@ ImageData WireframeRayTracePlot::create_image() const // This array marks where the initial wireframe was drawn. We convolve it with // a filter that gets adjusted with the wireframe thickness in order to // thicken the lines. - xt::xtensor wireframe_initial({width, height}, 0); + tensor::Tensor wireframe_initial( + {static_cast(width), static_cast(height)}, 0); /* Holds all of the track segments for the current rendered line of pixels. * old_segments holds a copy of this_line_segments from the previous line. diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 13a1b0163b6..b881ced716b 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -18,6 +18,7 @@ #include "openmc/weight_windows.h" #include +#include namespace openmc { @@ -63,8 +64,7 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) // Create a new 2D tensor with the same size as the first // two dimensions of the 3D tensor - tally_volumes_[i] = - xt::xtensor::from_shape({shape[0], shape[1]}); + tally_volumes_[i] = tensor::Tensor({shape[0], shape[1]}); } } diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 47d901d63a0..44e161a8e63 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -10,6 +10,8 @@ #include "openmc/settings.h" #include "openmc/simulation.h" +#include + #include "openmc/distribution_spatial.h" #include "openmc/random_dist.h" #include "openmc/source.h" diff --git a/src/scattdata.cpp b/src/scattdata.cpp index 21b18cbd923..9aa09956d50 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -4,8 +4,7 @@ #include #include -#include "xtensor/xbuilder.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include "openmc/constants.h" #include "openmc/error.h" @@ -19,8 +18,8 @@ namespace openmc { // ScattData base-class methods //============================================================================== -void ScattData::base_init(int order, const xt::xtensor& in_gmin, - const xt::xtensor& in_gmax, const double_2dvec& in_energy, +void ScattData::base_init(int order, const tensor::Tensor& in_gmin, + const tensor::Tensor& in_gmax, const double_2dvec& in_energy, const double_2dvec& in_mult) { size_t groups = in_energy.size(); @@ -63,23 +62,26 @@ void ScattData::base_init(int order, const xt::xtensor& in_gmin, void ScattData::base_combine(size_t max_order, size_t order_dim, const vector& those_scatts, const vector& scalars, - xt::xtensor& in_gmin, xt::xtensor& in_gmax, + tensor::Tensor& in_gmin, tensor::Tensor& in_gmax, double_2dvec& sparse_mult, double_3dvec& sparse_scatter) { size_t groups = those_scatts[0]->energy.size(); // Now allocate and zero our storage spaces - xt::xtensor this_nuscatt_matrix({groups, groups, order_dim}, 0.); - xt::xtensor this_nuscatt_P0({groups, groups}, 0.); - xt::xtensor this_scatt_P0({groups, groups}, 0.); - xt::xtensor this_mult({groups, groups}, 1.); + tensor::Tensor this_nuscatt_matrix = + tensor::zeros({groups, groups, order_dim}); + tensor::Tensor this_nuscatt_P0 = + tensor::zeros({groups, groups}); + tensor::Tensor this_scatt_P0 = + tensor::zeros({groups, groups}); + tensor::Tensor this_mult = tensor::ones({groups, groups}); // Build the dense scattering and multiplicity matrices for (int i = 0; i < those_scatts.size(); i++) { ScattData* that = those_scatts[i]; // Build the dense matrix for that object - xt::xtensor that_matrix = that->get_matrix(max_order); + tensor::Tensor that_matrix = that->get_matrix(max_order); // Now add that to this for the nu-scatter matrix this_nuscatt_matrix += scalars[i] * that_matrix; @@ -97,7 +99,7 @@ void ScattData::base_combine(size_t max_order, size_t order_dim, // Now we have the dense nuscatt and scatt, we can easily compute the // multiplicity matrix by dividing the two and fixing any nans - this_mult = xt::nan_to_num(this_nuscatt_P0 / this_scatt_P0); + this_mult = tensor::nan_to_num(this_nuscatt_P0 / this_scatt_P0); // 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. @@ -106,7 +108,7 @@ void ScattData::base_combine(size_t max_order, size_t order_dim, int gmin_; for (gmin_ = 0; gmin_ < groups; gmin_++) { bool non_zero = false; - for (int l = 0; l < this_nuscatt_matrix.shape()[2]; l++) { + for (int l = 0; l < this_nuscatt_matrix.shape(2); l++) { if (this_nuscatt_matrix(gin, gmin_, l) != 0.) { non_zero = true; break; @@ -118,7 +120,7 @@ void ScattData::base_combine(size_t max_order, size_t order_dim, int gmax_; for (gmax_ = groups - 1; gmax_ >= 0; gmax_--) { bool non_zero = false; - for (int l = 0; l < this_nuscatt_matrix.shape()[2]; l++) { + for (int l = 0; l < this_nuscatt_matrix.shape(2); l++) { if (this_nuscatt_matrix(gin, gmax_, l) != 0.) { non_zero = true; break; @@ -143,8 +145,8 @@ void ScattData::base_combine(size_t max_order, size_t order_dim, sparse_mult[gin].resize(gmax_ - gmin_ + 1); int i_gout = 0; for (int gout = gmin_; gout <= gmax_; gout++) { - sparse_scatter[gin][i_gout].resize(this_nuscatt_matrix.shape()[2]); - for (int l = 0; l < this_nuscatt_matrix.shape()[2]; l++) { + sparse_scatter[gin][i_gout].resize(this_nuscatt_matrix.shape(2)); + for (int l = 0; l < this_nuscatt_matrix.shape(2); l++) { sparse_scatter[gin][i_gout][l] = this_nuscatt_matrix(gin, gout, l); } sparse_mult[gin][i_gout] = this_mult(gin, gout); @@ -227,8 +229,8 @@ double ScattData::get_xs( // ScattDataLegendre methods //============================================================================== -void ScattDataLegendre::init(const xt::xtensor& in_gmin, - const xt::xtensor& in_gmax, const double_2dvec& in_mult, +void ScattDataLegendre::init(const tensor::Tensor& in_gmin, + const tensor::Tensor& in_gmax, const double_2dvec& in_mult, const double_3dvec& coeffs) { size_t groups = coeffs.size(); @@ -239,7 +241,7 @@ void ScattDataLegendre::init(const xt::xtensor& in_gmin, // Get the scattering cross section value by summing the un-normalized P0 // coefficient in the variable matrix over all outgoing groups. - scattxs = xt::zeros({groups}); + scattxs = tensor::zeros({groups}); for (int gin = 0; gin < groups; gin++) { int num_groups = in_gmax[gin] - in_gmin[gin] + 1; for (int i_gout = 0; i_gout < num_groups; i_gout++) { @@ -386,8 +388,8 @@ void ScattDataLegendre::combine( size_t groups = those_scatts[0]->energy.size(); - xt::xtensor in_gmin({groups}, 0); - xt::xtensor in_gmax({groups}, 0); + tensor::Tensor in_gmin({groups}, 0); + tensor::Tensor in_gmax({groups}, 0); double_3dvec sparse_scatter(groups); double_2dvec sparse_mult(groups); @@ -404,12 +406,13 @@ void ScattDataLegendre::combine( //============================================================================== -xt::xtensor ScattDataLegendre::get_matrix(size_t max_order) +tensor::Tensor ScattDataLegendre::get_matrix(size_t max_order) { // Get the sizes and initialize the data to 0 size_t groups = energy.size(); size_t order_dim = max_order + 1; - xt::xtensor matrix({groups, groups, order_dim}, 0.); + tensor::Tensor matrix = + tensor::zeros({groups, groups, order_dim}); for (int gin = 0; gin < groups; gin++) { for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) { @@ -427,8 +430,8 @@ xt::xtensor ScattDataLegendre::get_matrix(size_t max_order) // ScattDataHistogram methods //============================================================================== -void ScattDataHistogram::init(const xt::xtensor& in_gmin, - const xt::xtensor& in_gmax, const double_2dvec& in_mult, +void ScattDataHistogram::init(const tensor::Tensor& in_gmin, + const tensor::Tensor& in_gmax, const double_2dvec& in_mult, const double_3dvec& coeffs) { size_t groups = coeffs.size(); @@ -439,7 +442,7 @@ void ScattDataHistogram::init(const xt::xtensor& in_gmin, // Get the scattering cross section value by summing the distribution // over all the histogram bins in angle and outgoing energy groups - scattxs = xt::zeros({groups}); + scattxs = tensor::zeros({groups}); for (int gin = 0; gin < groups; gin++) { for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) { scattxs[gin] += std::accumulate( @@ -468,7 +471,7 @@ void ScattDataHistogram::init(const xt::xtensor& in_gmin, ScattData::base_init(order, in_gmin, in_gmax, in_energy, in_mult); // Build the angular distribution mu values - mu = xt::linspace(-1., 1., order + 1); + mu = tensor::linspace(-1., 1., order + 1); dmu = 2. / order; // Calculate f(mu) and integrate it so we can avoid rejection sampling @@ -513,7 +516,7 @@ double 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.shape()[0] - 2; + imu = this->mu.shape(0) - 2; } else { imu = std::floor((mu + 1.) / dmu + 1.) - 1; } @@ -559,13 +562,13 @@ void ScattDataHistogram::sample( //============================================================================== -xt::xtensor ScattDataHistogram::get_matrix(size_t max_order) +tensor::Tensor ScattDataHistogram::get_matrix(size_t max_order) { // Get the sizes and initialize the data to 0 size_t groups = energy.size(); // We ignore the requested order for Histogram and Tabular representations size_t order_dim = get_order(); - xt::xtensor matrix({groups, groups, order_dim}, 0); + tensor::Tensor 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++) { @@ -600,8 +603,8 @@ void ScattDataHistogram::combine( size_t groups = those_scatts[0]->energy.size(); - xt::xtensor in_gmin({groups}, 0); - xt::xtensor in_gmax({groups}, 0); + tensor::Tensor in_gmin({groups}, 0); + tensor::Tensor in_gmax({groups}, 0); double_3dvec sparse_scatter(groups); double_2dvec sparse_mult(groups); @@ -620,8 +623,8 @@ void ScattDataHistogram::combine( // ScattDataTabular methods //============================================================================== -void ScattDataTabular::init(const xt::xtensor& in_gmin, - const xt::xtensor& in_gmax, const double_2dvec& in_mult, +void ScattDataTabular::init(const tensor::Tensor& in_gmin, + const tensor::Tensor& in_gmax, const double_2dvec& in_mult, const double_3dvec& coeffs) { size_t groups = coeffs.size(); @@ -631,12 +634,12 @@ void ScattDataTabular::init(const xt::xtensor& in_gmin, double_3dvec matrix = coeffs; // Build the angular distribution mu values - mu = xt::linspace(-1., 1., order); + mu = tensor::linspace(-1., 1., order); dmu = 2. / (order - 1); // Get the scattering cross section value by integrating the distribution // over all mu points and then combining over all outgoing groups - scattxs = xt::zeros({groups}); + scattxs = tensor::zeros({groups}); for (int gin = 0; gin < groups; gin++) { for (int i_gout = 0; i_gout < matrix[gin].size(); i_gout++) { for (int imu = 1; imu < order; imu++) { @@ -713,7 +716,7 @@ double 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.shape()[0] - 2; + imu = this->mu.shape(0) - 2; } else { imu = std::floor((mu + 1.) / dmu + 1.) - 1; } @@ -734,7 +737,7 @@ void ScattDataTabular::sample( sample_energy(gin, gout, i_gout, seed); // Determine the outgoing cosine bin - int NP = this->mu.shape()[0]; + int NP = this->mu.shape(0); double xi = prn(seed); double c_k = dist[gin][i_gout][0]; @@ -776,13 +779,14 @@ void ScattDataTabular::sample( //============================================================================== -xt::xtensor ScattDataTabular::get_matrix(size_t max_order) +tensor::Tensor ScattDataTabular::get_matrix(size_t max_order) { // Get the sizes and initialize the data to 0 size_t groups = energy.size(); // We ignore the requested order for Histogram and Tabular representations size_t order_dim = get_order(); - xt::xtensor matrix({groups, groups, order_dim}, 0.); + tensor::Tensor matrix = + tensor::zeros({groups, groups, order_dim}); for (int gin = 0; gin < groups; gin++) { for (int i_gout = 0; i_gout < energy[gin].size(); i_gout++) { @@ -816,8 +820,8 @@ void ScattDataTabular::combine( size_t groups = those_scatts[0]->energy.size(); - xt::xtensor in_gmin({groups}, 0); - xt::xtensor in_gmax({groups}, 0); + tensor::Tensor in_gmin({groups}, 0); + tensor::Tensor in_gmax({groups}, 0); double_3dvec sparse_scatter(groups); double_2dvec sparse_mult(groups); @@ -854,7 +858,7 @@ void convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab) tab.scattxs = leg.scattxs; // Build mu and dmu - tab.mu = xt::linspace(-1., 1., n_mu); + tab.mu = tensor::linspace(-1., 1., n_mu); tab.dmu = 2. / (n_mu - 1); // Calculate f(mu) and integrate it so we can avoid rejection sampling diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index e820419b484..cc0ab8af19f 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -5,8 +5,7 @@ #include // for size_t #include // for back_inserter -#include "xtensor/xarray.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include "openmc/endf.h" #include "openmc/hdf5_interface.h" @@ -26,11 +25,11 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) hid_t dset = open_dataset(group, "energy"); // Get interpolation parameters - xt::xarray temp; + tensor::Tensor temp; read_attribute(dset, "interpolation", temp); - auto temp_b = xt::view(temp, 0); // view of breakpoints - auto temp_i = xt::view(temp, 1); // view of interpolation parameters + tensor::View temp_b = temp.slice(0); // breakpoints + tensor::View temp_i = temp.slice(1); // interpolation parameters std::copy(temp_b.begin(), temp_b.end(), std::back_inserter(breakpoints_)); for (const auto i : temp_i) @@ -51,12 +50,12 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) read_attribute(dset, "interpolation", interp); read_attribute(dset, "n_discrete_lines", n_discrete); - xt::xarray eout; + tensor::Tensor eout; read_dataset(dset, eout); close_dataset(dset); // Read angle distributions - xt::xarray mu; + tensor::Tensor mu; read_dataset(group, "mu", mu); for (int i = 0; i < n_energy; ++i) { @@ -66,7 +65,7 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) if (i < n_energy - 1) { n = offsets[i + 1] - j; } else { - n = eout.shape()[1] - j; + n = eout.shape(1) - j; } // Assign interpolation scheme and number of discrete lines @@ -75,9 +74,9 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) d.n_discrete = n_discrete[i]; // Copy data - d.e_out = xt::view(eout, 0, xt::range(j, j + n)); - d.p = xt::view(eout, 1, xt::range(j, j + n)); - d.c = xt::view(eout, 2, xt::range(j, j + n)); + d.e_out = eout.slice(0, tensor::range(j, j + n)); + d.p = eout.slice(1, tensor::range(j, j + n)); + d.c = eout.slice(2, tensor::range(j, j + n)); // To get answers that match ACE data, for now we still use the tabulated // CDF values that were passed through to the HDF5 library. At a later @@ -119,10 +118,10 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) // Determine offset and size of distribution int offset_mu = std::lround(eout(4, offsets[i] + j)); int m; - if (offsets[i] + j + 1 < eout.shape()[1]) { + if (offsets[i] + j + 1 < eout.shape(1)) { m = std::lround(eout(4, offsets[i] + j + 1)) - offset_mu; } else { - m = mu.shape()[1] - offset_mu; + m = mu.shape(1) - offset_mu; } // For incoherent inelastic thermal scattering, the angle distributions @@ -133,9 +132,12 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) interp_mu = 1; auto interp = int2interp(interp_mu); - auto xs = xt::view(mu, 0, xt::range(offset_mu, offset_mu + m)); - auto ps = xt::view(mu, 1, xt::range(offset_mu, offset_mu + m)); - auto cs = xt::view(mu, 2, xt::range(offset_mu, offset_mu + m)); + tensor::View xs = + mu.slice(0, tensor::range(offset_mu, offset_mu + m)); + tensor::View ps = + mu.slice(1, tensor::range(offset_mu, offset_mu + m)); + tensor::View cs = + mu.slice(2, tensor::range(offset_mu, offset_mu + m)); vector x {xs.begin(), xs.end()}; vector p {ps.begin(), ps.end()}; diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index 6ac91e665b9..8470c8c18e6 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -5,8 +5,7 @@ #include // for size_t #include // for back_inserter -#include "xtensor/xarray.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include "openmc/hdf5_interface.h" #include "openmc/math_functions.h" @@ -27,11 +26,11 @@ KalbachMann::KalbachMann(hid_t group) hid_t dset = open_dataset(group, "energy"); // Get interpolation parameters - xt::xarray temp; + tensor::Tensor temp; read_attribute(dset, "interpolation", temp); - auto temp_b = xt::view(temp, 0); // view of breakpoints - auto temp_i = xt::view(temp, 1); // view of interpolation parameters + tensor::View temp_b = temp.slice(0); // breakpoints + tensor::View temp_i = temp.slice(1); // interpolation parameters std::copy(temp_b.begin(), temp_b.end(), std::back_inserter(breakpoints_)); for (const auto i : temp_i) @@ -52,7 +51,7 @@ KalbachMann::KalbachMann(hid_t group) read_attribute(dset, "interpolation", interp); read_attribute(dset, "n_discrete_lines", n_discrete); - xt::xarray eout; + tensor::Tensor eout; read_dataset(dset, eout); close_dataset(dset); @@ -63,7 +62,7 @@ KalbachMann::KalbachMann(hid_t group) if (i < n_energy - 1) { n = offsets[i + 1] - j; } else { - n = eout.shape()[1] - j; + n = eout.shape(1) - j; } // Assign interpolation scheme and number of discrete lines @@ -72,11 +71,11 @@ KalbachMann::KalbachMann(hid_t group) d.n_discrete = n_discrete[i]; // Copy data - d.e_out = xt::view(eout, 0, xt::range(j, j + n)); - d.p = xt::view(eout, 1, xt::range(j, j + n)); - d.c = xt::view(eout, 2, xt::range(j, j + n)); - d.r = xt::view(eout, 3, xt::range(j, j + n)); - d.a = xt::view(eout, 4, xt::range(j, j + n)); + d.e_out = eout.slice(0, tensor::range(j, j + n)); + d.p = eout.slice(1, tensor::range(j, j + n)); + d.c = eout.slice(2, tensor::range(j, j + n)); + d.r = eout.slice(3, tensor::range(j, j + n)); + d.a = eout.slice(4, tensor::range(j, j + n)); // To get answers that match ACE data, for now we still use the tabulated // CDF values that were passed through to the HDF5 library. At a later diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index 030d398aabe..2ab7d8a63e1 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -5,7 +5,7 @@ #include "openmc/random_lcg.h" #include "openmc/search.h" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include #include // for log, exp @@ -85,7 +85,7 @@ void IncoherentElasticAEDiscrete::sample( // incoming energies. // Sample outgoing cosine bin - int n_mu = mu_out_.shape()[1]; + int n_mu = mu_out_.shape(1); int k = prn(seed) * n_mu; // Rather than use the sampled discrete mu directly, it is smeared over @@ -145,7 +145,7 @@ void IncoherentInelasticAEDiscrete::sample( // probability of 1). Otherwise, each bin is equally probable. int j; - int n = energy_out_.shape()[1]; + int n = energy_out_.shape(1); if (!skewed_) { // All bins equally likely j = prn(seed) * n; @@ -178,7 +178,7 @@ void IncoherentInelasticAEDiscrete::sample( E_out = (1 - f) * E_ij + f * E_i1j; // Sample outgoing cosine bin - int m = mu_out_.shape()[2]; + int m = mu_out_.shape(2); int k = prn(seed) * m; // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1] @@ -218,11 +218,11 @@ IncoherentInelasticAE::IncoherentInelasticAE(hid_t group) // On first pass, allocate space for angles if (j == 0) { auto n_mu = adist->x().size(); - d.mu = xt::empty({d.n_e_out, n_mu}); + d.mu = tensor::Tensor({d.n_e_out, n_mu}); } // Copy outgoing angles - auto mu_j = xt::view(d.mu, j); + tensor::View mu_j = d.mu.slice(j); std::copy(adist->x().begin(), adist->x().end(), mu_j.begin()); } } @@ -287,7 +287,7 @@ void IncoherentInelasticAE::sample( } // Sample outgoing cosine bin - int n_mu = distribution_[l].mu.shape()[1]; + int n_mu = distribution_[l].mu.shape(1); std::size_t k = prn(seed) * n_mu; // Rather than use the sampled discrete mu directly, it is smeared over diff --git a/src/simulation.cpp b/src/simulation.cpp index 18e40a8bc73..d0d6f037f9e 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -30,7 +30,7 @@ #ifdef _OPENMP #include #endif -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #ifdef OPENMC_MPI #include @@ -413,7 +413,7 @@ void finalize_batch() // Reset global tally results if (simulation::current_batch <= settings::n_inactive) { - xt::view(simulation::global_tallies, xt::all()) = 0.0; + simulation::global_tallies.fill(0.0); simulation::n_realizations = 0; } diff --git a/src/source.cpp b/src/source.cpp index 8bd9c789352..5300a685f78 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -10,7 +10,7 @@ #include // for dlopen, dlsym, dlclose, dlerror #endif -#include "xtensor/xadapt.hpp" +#include "openmc/tensor.h" #include #include "openmc/bank.h" @@ -400,8 +400,9 @@ SourceSite IndependentSource::sample(uint64_t* seed) const auto p = particle_.transport_index(); auto energy_ptr = dynamic_cast(energy_.get()); if (energy_ptr) { - auto energies = xt::adapt(energy_ptr->x()); - if (xt::any(energies > data::energy_max[p])) { + auto energies = + tensor::Tensor(energy_ptr->x().data(), energy_ptr->x().size()); + if ((energies > data::energy_max[p]).any()) { fatal_error("Source energy above range of energies of at least " "one cross section table"); } diff --git a/src/state_point.cpp b/src/state_point.cpp index 455299529b3..c0d8ab5b2fa 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -4,8 +4,7 @@ #include // for int64_t #include -#include "xtensor/xbuilder.hpp" // for empty_like -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include #include "openmc/bank.h" @@ -277,8 +276,8 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) std::string name = "tally " + std::to_string(tally->id_); hid_t tally_group = open_group(tallies_group, name.c_str()); auto& results = tally->results_; - write_tally_results(tally_group, results.shape()[0], - results.shape()[1], results.shape()[2], results.data()); + write_tally_results(tally_group, results.shape(0), results.shape(1), + results.shape(2), results.data()); close_group(tally_group); } } else { @@ -517,8 +516,8 @@ extern "C" int openmc_statepoint_load(const char* filename) tally->writable_ = false; } else { auto& results = tally->results_; - read_tally_results(tally_group, results.shape()[0], - results.shape()[1], results.shape()[2], results.data()); + read_tally_results(tally_group, results.shape(0), results.shape(1), + results.shape(2), results.data()); read_dataset(tally_group, "n_realizations", tally->n_realizations_); close_group(tally_group); @@ -827,7 +826,7 @@ void write_unstructured_mesh_results() // construct result vectors vector mean_vec(umesh->n_bins()), std_dev_vec(umesh->n_bins()); - for (int j = 0; j < tally->results_.shape()[0]; j++) { + for (int j = 0; j < tally->results_.shape(0); j++) { // get the volume for this bin double volume = umesh->volume(j); // compute the mean @@ -889,7 +888,7 @@ void write_tally_results_nr(hid_t file_id) #ifdef OPENMC_MPI // Reduce global tallies - xt::xtensor gt_reduced = xt::empty_like(gt); + tensor::Tensor gt_reduced({N_GLOBAL_TALLIES, 3}); MPI_Reduce(gt.data(), gt_reduced.data(), gt.size(), MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); @@ -918,13 +917,18 @@ void write_tally_results_nr(hid_t file_id) write_attribute(file_id, "tallies_present", 1); } - // Get view of accumulated tally values - auto values_view = xt::view(t->results_, xt::all(), xt::all(), - xt::range(static_cast(TallyResult::SUM), - static_cast(TallyResult::SUM_SQ) + 1)); - - // Make copy of tally values in contiguous array - xt::xtensor values = values_view; + // Copy the SUM and SUM_SQ columns from the tally results into a + // contiguous array for MPI reduction + const int r_start = static_cast(TallyResult::SUM); + const int r_end = static_cast(TallyResult::SUM_SQ) + 1; + const size_t r_count = r_end - r_start; + const size_t ni = t->results_.shape(0); + const size_t nj = t->results_.shape(1); + tensor::Tensor values({ni, nj, r_count}); + for (size_t i = 0; i < ni; i++) + for (size_t j = 0; j < nj; j++) + for (size_t r = 0; r < r_count; r++) + values(i, j, r) = t->results_(i, j, r_start + r); if (mpi::master) { // Open group for tally @@ -938,19 +942,22 @@ void write_tally_results_nr(hid_t file_id) MPI_SUM, 0, mpi::intracomm); #endif - // At the end of the simulation, store the results back in the - // regular TallyResults array + // At the end of the simulation, store the reduced results back + // into the tally results array if (simulation::current_batch == settings::n_max_batches || simulation::satisfy_triggers) { - values_view = values; + for (size_t i = 0; i < ni; i++) + for (size_t j = 0; j < nj; j++) + for (size_t r = 0; r < r_count; r++) + t->results_(i, j, r_start + r) = values(i, j, r); } - // Put in temporary tally result - xt::xtensor results_copy = xt::zeros_like(t->results_); - auto copy_view = xt::view(results_copy, xt::all(), xt::all(), - xt::range(static_cast(TallyResult::SUM), - static_cast(TallyResult::SUM_SQ) + 1)); - copy_view = values; + // Put reduced values into a full-sized copy for writing to HDF5 + tensor::Tensor results_copy = tensor::zeros_like(t->results_); + for (size_t i = 0; i < ni; i++) + for (size_t j = 0; j < nj; j++) + for (size_t r = 0; r < r_count; r++) + results_copy(i, j, r_start + r) = values(i, j, r); // Write reduced tally results to file auto shape = results_copy.shape(); diff --git a/src/tallies/filter_cell_instance.cpp b/src/tallies/filter_cell_instance.cpp index 316a758d11a..928bfb6c5e7 100644 --- a/src/tallies/filter_cell_instance.cpp +++ b/src/tallies/filter_cell_instance.cpp @@ -9,6 +9,7 @@ #include "openmc/cell.h" #include "openmc/error.h" #include "openmc/geometry.h" +#include "openmc/tensor.h" #include "openmc/xml_interface.h" namespace openmc { @@ -108,7 +109,7 @@ void CellInstanceFilter::to_statepoint(hid_t filter_group) const { Filter::to_statepoint(filter_group); size_t n = cell_instances_.size(); - xt::xtensor data({n, 2}); + tensor::Tensor data({n, 2}); for (int64_t i = 0; i < n; ++i) { const auto& x = cell_instances_[i]; data(i, 0) = model::cells[x.index_cell]->id_; diff --git a/src/tallies/filter_meshmaterial.cpp b/src/tallies/filter_meshmaterial.cpp index 6e1f30380f6..b45bb4164eb 100644 --- a/src/tallies/filter_meshmaterial.cpp +++ b/src/tallies/filter_meshmaterial.cpp @@ -1,5 +1,6 @@ #include "openmc/tallies/filter_meshmaterial.h" +#include #include // for move #include @@ -10,6 +11,7 @@ #include "openmc/error.h" #include "openmc/material.h" #include "openmc/mesh.h" +#include "openmc/tensor.h" #include "openmc/xml_interface.h" namespace openmc { @@ -161,7 +163,7 @@ void MeshMaterialFilter::to_statepoint(hid_t filter_group) const write_dataset(filter_group, "mesh", model::meshes[mesh_]->id_); size_t n = bins_.size(); - xt::xtensor data({n, 2}); + tensor::Tensor data({n, 2}); for (int64_t i = 0; i < n; ++i) { const auto& x = bins_[i]; data(i, 0) = x.index_element; diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 8cdf3a9050c..296c3fb5051 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -35,9 +35,7 @@ #include "openmc/tallies/filter_time.h" #include "openmc/xml_interface.h" -#include "xtensor/xadapt.hpp" -#include "xtensor/xbuilder.hpp" // for empty_like -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include #include // for max, set_union @@ -69,7 +67,7 @@ vector time_grid; } // namespace model namespace simulation { -xt::xtensor_fixed> global_tallies; +tensor::StaticTensor2D global_tallies; int32_t n_realizations {0}; } // namespace simulation @@ -806,9 +804,11 @@ void Tally::init_results() { int n_scores = scores_.size() * nuclides_.size(); if (higher_moments_) { - results_ = xt::empty({n_filter_bins_, n_scores, 5}); + results_ = tensor::Tensor({static_cast(n_filter_bins_), + static_cast(n_scores), size_t {5}}); } else { - results_ = xt::empty({n_filter_bins_, n_scores, 3}); + results_ = tensor::Tensor({static_cast(n_filter_bins_), + static_cast(n_scores), size_t {3}}); } } @@ -816,7 +816,7 @@ void Tally::reset() { n_realizations_ = 0; if (results_.size() != 0) { - xt::view(results_, xt::all()) = 0.0; + results_.fill(0.0); } } @@ -851,9 +851,9 @@ void Tally::accumulate() if (higher_moments_) { #pragma omp parallel for // filter bins (specific cell, energy bins) - for (int i = 0; i < results_.shape()[0]; ++i) { + for (int i = 0; i < results_.shape(0); ++i) { // score bins (flux, total reaction rate, fission reaction rate, etc.) - for (int j = 0; j < results_.shape()[1]; ++j) { + for (int j = 0; j < results_.shape(1); ++j) { double val = results_(i, j, TallyResult::VALUE) * norm; double val2 = val * val; results_(i, j, TallyResult::VALUE) = 0.0; @@ -866,9 +866,9 @@ void Tally::accumulate() } else { #pragma omp parallel for // filter bins (specific cell, energy bins) - for (int i = 0; i < results_.shape()[0]; ++i) { + for (int i = 0; i < results_.shape(0); ++i) { // score bins (flux, total reaction rate, fission reaction rate, etc.) - for (int j = 0; j < results_.shape()[1]; ++j) { + for (int j = 0; j < results_.shape(1); ++j) { double val = results_(i, j, TallyResult::VALUE) * norm; results_(i, j, TallyResult::VALUE) = 0.0; results_(i, j, TallyResult::SUM) += val; @@ -888,18 +888,18 @@ int Tally::score_index(const std::string& score) const return -1; } -xt::xarray Tally::get_reshaped_data() const +tensor::Tensor Tally::get_reshaped_data() const { - std::vector shape; + vector shape; for (auto f : filters()) { shape.push_back(model::tally_filters[f]->n_bins()); } // add number of scores and nuclides to tally - shape.push_back(results_.shape()[1]); - shape.push_back(results_.shape()[2]); + shape.push_back(results_.shape(1)); + shape.push_back(results_.shape(2)); - xt::xarray reshaped_results = results_; + tensor::Tensor reshaped_results = results_; reshaped_results.reshape(shape); return reshaped_results; } @@ -1004,13 +1004,14 @@ void reduce_tally_results() // Skip any tallies that are not active auto& tally {model::tallies[i_tally]}; - // Get view of accumulated tally values - auto values_view = xt::view(tally->results_, xt::all(), xt::all(), - static_cast(TallyResult::VALUE)); + // Extract 2D view of the VALUE column from the 3D results tensor, + // then copy into a contiguous array for MPI reduction + const int val_idx = static_cast(TallyResult::VALUE); + tensor::View val_view = + tally->results_.slice(tensor::all, tensor::all, val_idx); + tensor::Tensor values(val_view); - // Make copy of tally values in contiguous array - xt::xtensor values = values_view; - xt::xtensor values_reduced = xt::empty_like(values); + tensor::Tensor values_reduced(values.shape()); // Reduce contiguous set of tally results MPI_Reduce(values.data(), values_reduced.data(), values.size(), @@ -1018,9 +1019,9 @@ void reduce_tally_results() // Transfer values on master and reset on other ranks if (mpi::master) { - values_view = values_reduced; + val_view = values_reduced; } else { - values_view = 0.0; + val_view = 0.0; } } } @@ -1028,14 +1029,13 @@ void reduce_tally_results() // Note that global tallies are *always* reduced even when no_reduce option // is on. - // Get view of global tally values + // Get reference to global tallies auto& gt = simulation::global_tallies; - auto gt_values_view = - xt::view(gt, xt::all(), static_cast(TallyResult::VALUE)); + const int val_col = static_cast(TallyResult::VALUE); - // Make copy of values in contiguous array - xt::xtensor gt_values = gt_values_view; - xt::xtensor gt_values_reduced = xt::empty_like(gt_values); + // Copy VALUE column into contiguous array for MPI reduction + tensor::Tensor gt_values(gt.slice(tensor::all, val_col)); + tensor::Tensor gt_values_reduced({size_t {N_GLOBAL_TALLIES}}); // Reduce contiguous data MPI_Reduce(gt_values.data(), gt_values_reduced.data(), N_GLOBAL_TALLIES, @@ -1043,9 +1043,9 @@ void reduce_tally_results() // Transfer values on master and reset on other ranks if (mpi::master) { - gt_values_view = gt_values_reduced; + gt.slice(tensor::all, val_col) = gt_values_reduced; } else { - gt_values_view = 0.0; + gt.slice(tensor::all, val_col) = 0.0; } // We also need to determine the total starting weight of particles from the diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index ccc70fc8f14..fb851f78388 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -20,6 +20,7 @@ #include "openmc/tallies/filter_delayedgroup.h" #include "openmc/tallies/filter_energy.h" +#include #include namespace openmc { diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp index f1f83e2982c..6c54edd5e1d 100644 --- a/src/tallies/trigger.cpp +++ b/src/tallies/trigger.cpp @@ -71,7 +71,7 @@ void check_tally_triggers(double& ratio, int& tally_id, int& score) continue; const auto& results = t.results_; - for (auto filter_index = 0; filter_index < results.shape()[0]; + for (auto filter_index = 0; filter_index < results.shape(0); ++filter_index) { // Compute the tally uncertainty metrics. auto uncert_pair = diff --git a/src/thermal.cpp b/src/thermal.cpp index cbe0983ed65..6ed59f68693 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -3,12 +3,7 @@ #include // for sort, move, min, max, find #include // for round, sqrt, abs -#include "xtensor/xarray.hpp" -#include "xtensor/xbuilder.hpp" -#include "xtensor/xmath.hpp" -#include "xtensor/xsort.hpp" -#include "xtensor/xtensor.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include #include "openmc/constants.h" @@ -55,7 +50,7 @@ ThermalScattering::ThermalScattering( // Determine temperatures available auto dset_names = dataset_names(kT_group); auto n = dset_names.size(); - auto temps_available = xt::empty({n}); + auto temps_available = tensor::Tensor({n}); for (int i = 0; i < dset_names.size(); ++i) { // Read temperature value double T; @@ -82,7 +77,7 @@ ThermalScattering::ThermalScattering( // Determine actual temperatures to read for (const auto& T : temperature) { - auto i_closest = xt::argmin(xt::abs(temps_available - T))[0]; + auto i_closest = tensor::abs(temps_available - T).argmin(); auto temp_actual = temps_available[i_closest]; if (std::abs(temp_actual - T) < settings::temperature_tolerance) { if (std::find(temps_to_read.begin(), temps_to_read.end(), diff --git a/src/track_output.cpp b/src/track_output.cpp index e86f774f046..02e09223645 100644 --- a/src/track_output.cpp +++ b/src/track_output.cpp @@ -8,7 +8,7 @@ #include "openmc/simulation.h" #include "openmc/vector.h" -#include "xtensor/xtensor.hpp" +#include "openmc/tensor.h" #include #include diff --git a/src/urr.cpp b/src/urr.cpp index 02cef228099..b791eecb0da 100644 --- a/src/urr.cpp +++ b/src/urr.cpp @@ -25,7 +25,7 @@ UrrData::UrrData(hid_t group_id) // Read URR tables. The HDF5 format is a little // different from how we want it laid out in memory. // This array used to be called "prob_". - xt::xtensor tmp_prob; + tensor::Tensor tmp_prob; read_dataset(group_id, "table", tmp_prob); auto shape = tmp_prob.shape(); @@ -38,7 +38,7 @@ UrrData::UrrData(hid_t group_id) xs_values_.resize({n_energy, n_cdf_values}); // Now fill in the values. Using manual loops here since we might - // not have fancy xtensor slicing code written for GPU tensors. + // not have fancy tensor slicing code written for GPU tensors. // The below enum gives how URR tables are laid out in our HDF5 tables. enum class URRTableParam { CUM_PROB, diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp index 1deffb80488..f675ae78a2e 100644 --- a/src/volume_calc.cpp +++ b/src/volume_calc.cpp @@ -17,8 +17,7 @@ #include "openmc/timer.h" #include "openmc/xml_interface.h" -#include "xtensor/xadapt.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include #include // for copy @@ -242,7 +241,8 @@ vector VolumeCalculation::execute() const // non-zero auto n_nuc = settings::run_CE ? data::nuclides.size() : data::mg.nuclides_.size(); - xt::xtensor atoms({n_nuc, 2}, 0.0); + auto atoms = + tensor::zeros({static_cast(n_nuc), size_t {2}}); #ifdef OPENMC_MPI if (mpi::master) { @@ -452,9 +452,11 @@ void VolumeCalculation::to_hdf5( } // Create array of total # of atoms with uncertainty for each nuclide - xt::xtensor atom_data({n_nuc, 2}); - xt::view(atom_data, xt::all(), 0) = xt::adapt(result.atoms); - xt::view(atom_data, xt::all(), 1) = xt::adapt(result.uncertainty); + tensor::Tensor atom_data({static_cast(n_nuc), size_t {2}}); + for (size_t k = 0; k < static_cast(n_nuc); ++k) { + atom_data(k, 0) = result.atoms[k]; + atom_data(k, 1) = result.uncertainty[k]; + } // Write results write_dataset(group_id, "nuclides", nucnames); diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 5ca4addbbf7..c0516913204 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -6,12 +6,7 @@ #include #include -#include "xtensor/xdynamic_view.hpp" -#include "xtensor/xindex_view.hpp" -#include "xtensor/xio.hpp" -#include "xtensor/xmasked_view.hpp" -#include "xtensor/xnoalias.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include "openmc/error.h" #include "openmc/file_utils.h" @@ -265,8 +260,12 @@ WeightWindows* WeightWindows::from_hdf5( } wws->set_mesh(model::mesh_map[mesh_id]); - wws->lower_ww_ = xt::empty(wws->bounds_size()); - wws->upper_ww_ = xt::empty(wws->bounds_size()); + wws->lower_ww_ = + tensor::Tensor({static_cast(wws->bounds_size()[0]), + static_cast(wws->bounds_size()[1])}); + wws->upper_ww_ = + tensor::Tensor({static_cast(wws->bounds_size()[0]), + static_cast(wws->bounds_size()[1])}); read_dataset(ww_group, "lower_ww_bounds", wws->lower_ww_); read_dataset(ww_group, "upper_ww_bounds", wws->upper_ww_); @@ -301,9 +300,11 @@ void WeightWindows::allocate_ww_bounds() "Size of weight window bounds is zero for WeightWindows {}", id()); warning(msg); } - lower_ww_ = xt::empty(shape); + lower_ww_ = tensor::Tensor( + {static_cast(shape[0]), static_cast(shape[1])}); lower_ww_.fill(-1); - upper_ww_ = xt::empty(shape); + upper_ww_ = tensor::Tensor( + {static_cast(shape[0]), static_cast(shape[1])}); upper_ww_.fill(-1); } @@ -448,8 +449,8 @@ void WeightWindows::check_bounds(const T& bounds) const } } -void WeightWindows::set_bounds(const xt::xtensor& lower_bounds, - const xt::xtensor& upper_bounds) +void WeightWindows::set_bounds(const tensor::Tensor& lower_bounds, + const tensor::Tensor& upper_bounds) { this->check_bounds(lower_bounds, upper_bounds); @@ -460,7 +461,7 @@ void WeightWindows::set_bounds(const xt::xtensor& lower_bounds, } void WeightWindows::set_bounds( - const xt::xtensor& lower_bounds, double ratio) + const tensor::Tensor& lower_bounds, double ratio) { this->check_bounds(lower_bounds); @@ -475,14 +476,16 @@ void WeightWindows::set_bounds( { check_bounds(lower_bounds, upper_bounds); auto shape = this->bounds_size(); - lower_ww_ = xt::empty(shape); - upper_ww_ = xt::empty(shape); - - // set new weight window values - xt::view(lower_ww_, xt::all()) = - xt::adapt(lower_bounds.data(), lower_ww_.shape()); - xt::view(upper_ww_, xt::all()) = - xt::adapt(upper_bounds.data(), upper_ww_.shape()); + lower_ww_ = tensor::Tensor( + {static_cast(shape[0]), static_cast(shape[1])}); + upper_ww_ = tensor::Tensor( + {static_cast(shape[0]), static_cast(shape[1])}); + + // Copy weight window values from input spans into the tensors + std::copy(lower_bounds.data(), lower_bounds.data() + lower_ww_.size(), + lower_ww_.data()); + std::copy(upper_bounds.data(), upper_bounds.data() + upper_ww_.size(), + upper_ww_.data()); } void WeightWindows::set_bounds(span lower_bounds, double ratio) @@ -490,14 +493,16 @@ void WeightWindows::set_bounds(span lower_bounds, double ratio) this->check_bounds(lower_bounds); auto shape = this->bounds_size(); - lower_ww_ = xt::empty(shape); - upper_ww_ = xt::empty(shape); - - // set new weight window values - xt::view(lower_ww_, xt::all()) = - xt::adapt(lower_bounds.data(), lower_ww_.shape()); - xt::view(upper_ww_, xt::all()) = - xt::adapt(lower_bounds.data(), upper_ww_.shape()); + lower_ww_ = tensor::Tensor( + {static_cast(shape[0]), static_cast(shape[1])}); + upper_ww_ = tensor::Tensor( + {static_cast(shape[0]), static_cast(shape[1])}); + + // Copy lower bounds into both arrays, then scale upper by ratio + std::copy(lower_bounds.data(), lower_bounds.data() + lower_ww_.size(), + lower_ww_.data()); + std::copy(lower_bounds.data(), lower_bounds.data() + upper_ww_.size(), + upper_ww_.data()); upper_ww_ *= ratio; } @@ -510,8 +515,8 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value, this->check_tally_update_compatibility(tally); // Dimensions of weight window arrays - int e_bins = lower_ww_.shape()[0]; - int64_t mesh_bins = lower_ww_.shape()[1]; + int e_bins = lower_ww_.shape(0); + int64_t mesh_bins = lower_ww_.shape(1); // Initialize weight window arrays to -1.0 by default #pragma omp parallel for collapse(2) schedule(static) @@ -542,16 +547,16 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value, /////////////////////////// // Extract tally data // - // At the end of this section, the mean and rel_err array - // is a 2D view of tally data (n_e_groups, n_mesh_bins) + // At the end of this section, mean and rel_err are + // 2D tensors of tally data (n_e_groups, n_mesh_bins) // /////////////////////////// - // build a shape for a view of the tally results, this will always be + // build a shape for the tally results, this will always be // dimension 5 (3 filter dimensions, 1 score dimension, 1 results dimension) - // Look for the size of the last dimension of the results array - const auto& results_arr = tally->results(); - const int results_dim = static_cast(results_arr.shape()[2]); + // Look for the size of the last dimension of the results tensor + const auto& results = tally->results(); + const int results_dim = static_cast(results.shape(2)); std::array shape = {1, 1, 1, tally->n_scores(), results_dim}; // set the shape for the filters applied on the tally @@ -588,25 +593,14 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value, std::find(filter_types.begin(), filter_types.end(), FilterType::MESH) - filter_types.begin(); - // get a fully reshaped view of the tally according to tally ordering of - // filters - auto tally_values = xt::reshape_view(results_arr, shape); - - // get a that is (particle, energy, mesh, scores, values) - auto transposed_view = xt::transpose(tally_values, transpose); - - // determine the dimension and index of the particle data + // determine the index of the particle within its filter int particle_idx = 0; if (tally->has_filter(FilterType::PARTICLE)) { - // get the particle filter auto pf = tally->get_filter(); const auto& particles = pf->particles(); - // find the index of the particle that matches these weight windows auto p_it = std::find(particles.begin(), particles.end(), this->particle_type_); - // if the particle filter doesn't have particle data for the particle - // used on this weight windows instance, report an error if (p_it == particles.end()) { auto msg = fmt::format("Particle type '{}' not present on Filter {} for " "Tally {} used to update WeightWindows {}", @@ -614,17 +608,46 @@ void WeightWindows::update_weights(const Tally* tally, const std::string& value, fatal_error(msg); } - // use the index of the particle in the filter to down-select data later particle_idx = p_it - particles.begin(); } - // down-select data based on particle and score - auto sum = xt::dynamic_view( - transposed_view, {particle_idx, xt::all(), xt::all(), score_index, - static_cast(TallyResult::SUM)}); - auto sum_sq = xt::dynamic_view( - transposed_view, {particle_idx, xt::all(), xt::all(), score_index, - static_cast(TallyResult::SUM_SQ)}); + // The tally results array is 3D: (n_filter_combos, n_scores, n_result_types). + // The first dimension is a row-major flattening of up to 3 filter dimensions + // (particle, energy, mesh) whose storage order depends on which filters the + // tally has. We need to map our desired indices (particle, energy, mesh) + // into the correct flat filter combination index. + // + // transpose[i] tells us which storage position holds dimension i: + // i=0 -> particle, i=1 -> energy, i=2 -> mesh + // shape[j] gives the number of bins for filter storage position j. + + // Row-major strides for the 3 filter dimensions + const int stride0 = shape[1] * shape[2]; + const int stride1 = shape[2]; + + tensor::Tensor sum( + {static_cast(e_bins), static_cast(mesh_bins)}); + tensor::Tensor sum_sq( + {static_cast(e_bins), static_cast(mesh_bins)}); + + const int i_sum = static_cast(TallyResult::SUM); + const int i_sum_sq = static_cast(TallyResult::SUM_SQ); + + for (int e = 0; e < e_bins; e++) { + for (int64_t m = 0; m < mesh_bins; m++) { + // Place particle, energy, and mesh indices into their storage positions + std::array idx = {0, 0, 0}; + idx[transpose[0]] = particle_idx; + idx[transpose[1]] = e; + idx[transpose[2]] = static_cast(m); + + // Compute flat filter combination index (row-major over filter dims) + int flat = idx[0] * stride0 + idx[1] * stride1 + idx[2]; + + sum(e, m) = results(flat, score_index, i_sum); + sum_sq(e, m) = results(flat, score_index, i_sum_sq); + } + } int n = tally->n_realizations_; ////////////////////////////////////////////// @@ -1155,7 +1178,8 @@ extern "C" int openmc_weight_windows_set_bounds(int32_t index, return err; const auto& wws = variance_reduction::weight_windows[index]; - wws->set_bounds({lower_bounds, size}, {upper_bounds, size}); + wws->set_bounds(span(lower_bounds, size), + span(upper_bounds, size)); return 0; } diff --git a/src/wmp.cpp b/src/wmp.cpp index e9c4fb5e3e6..6f72e1ba4cc 100644 --- a/src/wmp.cpp +++ b/src/wmp.cpp @@ -33,22 +33,22 @@ WindowedMultipole::WindowedMultipole(hid_t group) // Read the "data" array. Use its shape to figure out the number of poles // and residue types in this data. read_dataset(group, "data", data_); - int n_residues = data_.shape()[1] - 1; + int n_residues = data_.shape(1) - 1; // Check to see if this data includes fission residues. fissionable_ = (n_residues == 3); // Read the "windows" array and use its shape to figure out the number of // windows. - xt::xtensor windows; + tensor::Tensor windows; read_dataset(group, "windows", windows); - int n_windows = windows.shape()[0]; + int n_windows = windows.shape(0); windows -= 1; // Adjust to 0-based indices // Read the "broaden_poly" arrays. - xt::xtensor broaden_poly; + tensor::Tensor broaden_poly; read_dataset(group, "broaden_poly", broaden_poly); - if (n_windows != broaden_poly.shape()[0]) { + if (n_windows != broaden_poly.shape(0)) { fatal_error("broaden_poly array shape is not consistent with the windows " "array shape in WMP library for " + name_ + "."); @@ -56,12 +56,12 @@ WindowedMultipole::WindowedMultipole(hid_t group) // Read the "curvefit" array. read_dataset(group, "curvefit", curvefit_); - if (n_windows != curvefit_.shape()[0]) { + if (n_windows != curvefit_.shape(0)) { fatal_error("curvefit array shape is not consistent with the windows " "array shape in WMP library for " + name_ + "."); } - fit_order_ = curvefit_.shape()[1] - 1; + fit_order_ = curvefit_.shape(1) - 1; // Check the code is compiling to work with sufficiently high fit order if (fit_order_ + 1 > MAX_POLY_COEFFICIENTS) { diff --git a/src/xsdata.cpp b/src/xsdata.cpp index 1929f51e6fa..33d063a7b9a 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -5,10 +5,7 @@ #include #include -#include "xtensor/xbuilder.hpp" -#include "xtensor/xindex_view.hpp" -#include "xtensor/xmath.hpp" -#include "xtensor/xview.hpp" +#include "openmc/tensor.h" #include "openmc/constants.h" #include "openmc/error.h" @@ -37,32 +34,32 @@ XsData::XsData(bool fissionable, AngleDistributionType scatter_format, } // allocate all [temperature][angle][in group] quantities vector shape {n_ang, n_g_}; - total = xt::zeros(shape); - absorption = xt::zeros(shape); - inverse_velocity = xt::zeros(shape); + total = tensor::zeros(shape); + absorption = tensor::zeros(shape); + inverse_velocity = tensor::zeros(shape); if (fissionable) { - fission = xt::zeros(shape); - nu_fission = xt::zeros(shape); - prompt_nu_fission = xt::zeros(shape); - kappa_fission = xt::zeros(shape); + fission = tensor::zeros(shape); + nu_fission = tensor::zeros(shape); + prompt_nu_fission = tensor::zeros(shape); + kappa_fission = tensor::zeros(shape); } // allocate decay_rate; [temperature][angle][delayed group] shape[1] = n_dg_; - decay_rate = xt::zeros(shape); + decay_rate = tensor::zeros(shape); if (fissionable) { shape = {n_ang, n_dg_, n_g_}; // allocate delayed_nu_fission; [temperature][angle][delay group][in group] - delayed_nu_fission = xt::zeros(shape); + delayed_nu_fission = tensor::zeros(shape); // chi_prompt; [temperature][angle][in group][out group] shape = {n_ang, n_g_, n_g_}; - chi_prompt = xt::zeros(shape); + chi_prompt = tensor::zeros(shape); // chi_delayed; [temperature][angle][delay group][in group][out group] shape = {n_ang, n_dg_, n_g_, n_g_}; - chi_delayed = xt::zeros(shape); + chi_delayed = tensor::zeros(shape); } for (int a = 0; a < n_ang; a++) { @@ -85,28 +82,30 @@ void XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, { // Reconstruct the dimension information so it doesn't need to be passed size_t n_ang = n_pol * n_azi; - size_t energy_groups = total.shape()[1]; + size_t energy_groups = total.shape(1); // Set the fissionable-specific data if (fissionable) { fission_from_hdf5(xsdata_grp, n_ang, is_isotropic); } // Get the non-fission-specific data - read_nd_vector(xsdata_grp, "decay-rate", decay_rate); - read_nd_vector(xsdata_grp, "absorption", absorption, true); - read_nd_vector(xsdata_grp, "inverse-velocity", inverse_velocity); + read_nd_tensor(xsdata_grp, "decay-rate", decay_rate); + read_nd_tensor(xsdata_grp, "absorption", absorption, true); + read_nd_tensor(xsdata_grp, "inverse-velocity", inverse_velocity); // Get scattering data scatter_from_hdf5( xsdata_grp, n_ang, scatter_format, final_scatter_format, order_data); - // Check absorption to ensure it is not 0 since it is often the - // denominator in tally methods - xt::filtration(absorption, xt::equal(absorption, 0.)) = 1.e-10; + // Replace zero absorption values with a small number to avoid + // division by zero in tally methods + for (size_t i = 0; i < absorption.size(); i++) + if (absorption.data()[i] == 0.0) + absorption.data()[i] = 1.e-10; // Get or calculate the total x/s if (object_exists(xsdata_grp, "total")) { - read_nd_vector(xsdata_grp, "total", total); + read_nd_tensor(xsdata_grp, "total", total); } else { for (size_t a = 0; a < n_ang; a++) { for (size_t gin = 0; gin < energy_groups; gin++) { @@ -115,8 +114,11 @@ void XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, } } - // Fix if total is 0, since it is in the denominator when tallying - xt::filtration(total, xt::equal(total, 0.)) = 1.e-10; + // Replace zero total cross sections with a small number to avoid + // division by zero in tally methods + for (size_t i = 0; i < total.size(); i++) + if (total.data()[i] == 0.0) + total.data()[i] = 1.e-10; } //============================================================================== @@ -127,21 +129,30 @@ void XsData::fission_vector_beta_from_hdf5( // Data is provided as nu-fission and chi with a beta for delayed info // Get chi - xt::xtensor temp_chi({n_ang, n_g_}, 0.); - read_nd_vector(xsdata_grp, "chi", temp_chi, true); + tensor::Tensor temp_chi = tensor::zeros({n_ang, n_g_}); + read_nd_tensor(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()); + // Normalize chi so it sums to 1 over outgoing groups for each angle + for (size_t a = 0; a < n_ang; a++) { + tensor::View row = temp_chi.slice(a); + row /= row.sum(); + } - // 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()); + // Replicate the energy spectrum across all incoming groups — the + // spectrum is independent of the incoming neutron energy + for (size_t a = 0; a < n_ang; a++) + for (size_t gin = 0; gin < n_g_; gin++) + chi_prompt.slice(a, gin) = temp_chi.slice(a); + + // Same spectrum for delayed neutrons, replicated across delayed groups + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg_; d++) + for (size_t gin = 0; gin < n_g_; gin++) + chi_delayed.slice(a, d, gin) = temp_chi.slice(a); // Get nu-fission - xt::xtensor temp_nufiss({n_ang, n_g_}, 0.); - read_nd_vector(xsdata_grp, "nu-fission", temp_nufiss, true); + tensor::Tensor temp_nufiss = tensor::zeros({n_ang, n_g_}); + read_nd_tensor(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"); @@ -151,26 +162,39 @@ void XsData::fission_vector_beta_from_hdf5( if (!is_isotropic) ndim_target += 2; if (beta_ndims == ndim_target) { - xt::xtensor temp_beta({n_ang, n_dg_}, 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()); + tensor::Tensor temp_beta = tensor::zeros({n_ang, n_dg_}); + read_nd_tensor(xsdata_grp, "beta", temp_beta, true); + + // prompt_nu_fission = (1 - sum_of_beta) * nu_fission + auto beta_sum = temp_beta.sum(1); + for (size_t a = 0; a < n_ang; a++) + for (size_t g = 0; g < n_g_; g++) + prompt_nu_fission(a, g) = temp_nufiss(a, g) * (1.0 - beta_sum(a)); + + // Delayed nu-fission is the outer product of the delayed neutron + // fraction (beta) and the fission production rate (nu-fission) + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg_; d++) + for (size_t g = 0; g < n_g_; g++) + delayed_nu_fission(a, d, g) = temp_beta(a, d) * temp_nufiss(a, g); } else if (beta_ndims == ndim_target + 1) { - xt::xtensor temp_beta({n_ang, n_dg_, n_g_}, 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()); + tensor::Tensor temp_beta = + tensor::zeros({n_ang, n_dg_, n_g_}); + read_nd_tensor(xsdata_grp, "beta", temp_beta, true); + + // prompt_nu_fission = (1 - sum_of_beta) * nu_fission + // Here beta is energy-dependent, so sum over delayed groups (axis 1) + auto beta_sum = temp_beta.sum(1); + for (size_t a = 0; a < n_ang; a++) + for (size_t g = 0; g < n_g_; g++) + prompt_nu_fission(a, g) = temp_nufiss(a, g) * (1.0 - beta_sum(a, g)); + + // Delayed nu-fission: beta is already energy-dependent [n_ang, n_dg, n_g], + // so scale each delayed group's beta by the total nu-fission for that group + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg_; d++) + for (size_t g = 0; g < n_g_; g++) + delayed_nu_fission(a, d, g) = temp_beta(a, d, g) * temp_nufiss(a, g); } } @@ -179,29 +203,42 @@ void XsData::fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang) // Data is provided separately as prompt + delayed nu-fission and chi // Get chi-prompt - xt::xtensor temp_chi_p({n_ang, n_g_}, 0.); - read_nd_vector(xsdata_grp, "chi-prompt", temp_chi_p, true); + tensor::Tensor temp_chi_p = tensor::zeros({n_ang, n_g_}); + read_nd_tensor(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()); + // Normalize prompt chi so it sums to 1 over outgoing groups for each angle + for (size_t a = 0; a < n_ang; a++) { + tensor::View row = temp_chi_p.slice(a); + row /= row.sum(); + } // Get chi-delayed - xt::xtensor temp_chi_d({n_ang, n_dg_, n_g_}, 0.); - read_nd_vector(xsdata_grp, "chi-delayed", temp_chi_d, true); + tensor::Tensor temp_chi_d = + tensor::zeros({n_ang, n_dg_, n_g_}); + read_nd_tensor(xsdata_grp, "chi-delayed", temp_chi_d, true); + + // Normalize delayed chi so it sums to 1 over outgoing groups for each + // angle and delayed group + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg_; d++) { + tensor::View row = temp_chi_d.slice(a, d); + row /= row.sum(); + } - // 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()); + // Replicate the prompt spectrum across all incoming groups + for (size_t a = 0; a < n_ang; a++) + for (size_t gin = 0; gin < n_g_; gin++) + chi_prompt.slice(a, gin) = temp_chi_p.slice(a); - // 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()); + // Replicate the delayed spectrum across all incoming groups + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg_; d++) + for (size_t gin = 0; gin < n_g_; gin++) + chi_delayed.slice(a, d, gin) = temp_chi_d.slice(a, d); // 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); + read_nd_tensor(xsdata_grp, "prompt-nu-fission", prompt_nu_fission, true); + read_nd_tensor(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) @@ -210,17 +247,22 @@ void XsData::fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang) // Therefore, the code only considers the data as prompt. // Get chi - xt::xtensor temp_chi({n_ang, n_g_}, 0.); - read_nd_vector(xsdata_grp, "chi", temp_chi, true); + tensor::Tensor temp_chi = tensor::zeros({n_ang, n_g_}); + read_nd_tensor(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()); + // Normalize chi so it sums to 1 over outgoing groups for each angle + for (size_t a = 0; a < n_ang; a++) { + tensor::View row = temp_chi.slice(a); + row /= row.sum(); + } - // 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()); + // Replicate the energy spectrum across all incoming groups + for (size_t a = 0; a < n_ang; a++) + for (size_t gin = 0; gin < n_g_; gin++) + chi_prompt.slice(a, gin) = temp_chi.slice(a); // Get nu-fission directly - read_nd_vector(xsdata_grp, "nu-fission", prompt_nu_fission, true); + read_nd_tensor(xsdata_grp, "nu-fission", prompt_nu_fission, true); } //============================================================================== @@ -231,8 +273,9 @@ void XsData::fission_matrix_beta_from_hdf5( // Data is provided as nu-fission and chi with a beta for delayed info // Get nu-fission matrix - xt::xtensor temp_matrix({n_ang, n_g_, n_g_}, 0.); - read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true); + tensor::Tensor temp_matrix = + tensor::zeros({n_ang, n_g_, n_g_}); + read_nd_tensor(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"); @@ -242,65 +285,92 @@ void XsData::fission_matrix_beta_from_hdf5( if (!is_isotropic) ndim_target += 2; if (beta_ndims == ndim_target) { - xt::xtensor temp_beta({n_ang, n_dg_}, 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()); + tensor::Tensor temp_beta = tensor::zeros({n_ang, n_dg_}); + read_nd_tensor(xsdata_grp, "beta", temp_beta, true); + + auto beta_sum = temp_beta.sum(1); + auto matrix_gout_sum = temp_matrix.sum(2); + + // prompt_nu_fission = sum_gout(matrix) * (1 - beta_total) + for (size_t a = 0; a < n_ang; a++) + for (size_t g = 0; g < n_g_; g++) + prompt_nu_fission(a, g) = matrix_gout_sum(a, g) * (1.0 - beta_sum(a)); + + // chi_prompt = (1 - beta_total) * nu-fission matrix (unnormalized) + for (size_t a = 0; a < n_ang; a++) + for (size_t gin = 0; gin < n_g_; gin++) + for (size_t gout = 0; gout < n_g_; gout++) + chi_prompt(a, gin, gout) = + (1.0 - beta_sum(a)) * temp_matrix(a, gin, gout); + + // Delayed nu-fission is the outer product of the delayed neutron + // fraction (beta) and the total fission rate summed over outgoing groups + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg_; d++) + for (size_t g = 0; g < n_g_; g++) + delayed_nu_fission(a, d, g) = temp_beta(a, d) * matrix_gout_sum(a, g); + + // chi_delayed = beta * nu-fission matrix, expanded across delayed groups + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg_; d++) + for (size_t gin = 0; gin < n_g_; gin++) + for (size_t gout = 0; gout < n_g_; gout++) + chi_delayed(a, d, gin, gout) = + temp_beta(a, d) * temp_matrix(a, gin, gout); } else if (beta_ndims == ndim_target + 1) { - xt::xtensor temp_beta({n_ang, n_dg_, n_g_}, 0.); - read_nd_vector(xsdata_grp, "beta", temp_beta, true); - - xt::xtensor temp_beta_sum({n_ang, n_g_}, 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()); + tensor::Tensor temp_beta = + tensor::zeros({n_ang, n_dg_, n_g_}); + read_nd_tensor(xsdata_grp, "beta", temp_beta, true); + + auto beta_sum = temp_beta.sum(1); + auto matrix_gout_sum = temp_matrix.sum(2); + + // prompt_nu_fission = sum_gout(matrix) * (1 - beta_total) + // Here beta is energy-dependent, so beta_sum is 2D [n_ang, n_g] + for (size_t a = 0; a < n_ang; a++) + for (size_t g = 0; g < n_g_; g++) + prompt_nu_fission(a, g) = + matrix_gout_sum(a, g) * (1.0 - beta_sum(a, g)); + + // chi_prompt = (1 - beta_sum) * nu-fission matrix (unnormalized) + for (size_t a = 0; a < n_ang; a++) + for (size_t gin = 0; gin < n_g_; gin++) + for (size_t gout = 0; gout < n_g_; gout++) + chi_prompt(a, gin, gout) = + (1.0 - beta_sum(a, gin)) * temp_matrix(a, gin, gout); + + // Delayed nu-fission: beta is energy-dependent [n_ang, n_dg, n_g], + // scale by total fission rate summed over outgoing groups + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg_; d++) + for (size_t g = 0; g < n_g_; g++) + delayed_nu_fission(a, d, g) = + temp_beta(a, d, g) * matrix_gout_sum(a, g); + + // chi_delayed = beta * nu-fission matrix, expanded across delayed groups + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg_; d++) + for (size_t gin = 0; gin < n_g_; gin++) + for (size_t gout = 0; gout < n_g_; gout++) + chi_delayed(a, d, gin, gout) = + temp_beta(a, d, gin) * temp_matrix(a, gin, gout); } - // Normalize both chis - chi_prompt /= - xt::view(xt::sum(chi_prompt, {2}), xt::all(), xt::all(), xt::newaxis()); + // Normalize chi_prompt so it sums to 1 over outgoing groups + for (size_t a = 0; a < n_ang; a++) + for (size_t gin = 0; gin < n_g_; gin++) { + tensor::View row = chi_prompt.slice(a, gin); + row /= row.sum(); + } - chi_delayed /= xt::view( - xt::sum(chi_delayed, {3}), xt::all(), xt::all(), xt::all(), xt::newaxis()); + // Normalize chi_delayed so it sums to 1 over outgoing groups + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg_; d++) + for (size_t gin = 0; gin < n_g_; gin++) { + tensor::View row = chi_delayed.slice(a, d, gin); + row /= row.sum(); + } } void XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang) @@ -308,28 +378,36 @@ void XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang) // Data is provided separately as prompt + delayed nu-fission and chi // Get the prompt nu-fission matrix - xt::xtensor temp_matrix_p({n_ang, n_g_, n_g_}, 0.); - read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_matrix_p, true); + tensor::Tensor temp_matrix_p = + tensor::zeros({n_ang, n_g_, n_g_}); + read_nd_tensor(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}); + prompt_nu_fission = temp_matrix_p.sum(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()); + // chi_prompt is the nu-fission matrix normalized over outgoing groups + for (size_t a = 0; a < n_ang; a++) + for (size_t gin = 0; gin < n_g_; gin++) + for (size_t gout = 0; gout < n_g_; gout++) + chi_prompt(a, gin, gout) = + temp_matrix_p(a, gin, gout) / prompt_nu_fission(a, gin); // Get the delayed nu-fission matrix - xt::xtensor temp_matrix_d({n_ang, n_dg_, n_g_, n_g_}, 0.); - read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_matrix_d, true); + tensor::Tensor temp_matrix_d = + tensor::zeros({n_ang, n_dg_, n_g_, n_g_}); + read_nd_tensor(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()); + delayed_nu_fission = temp_matrix_d.sum(3); + + // chi_delayed is the delayed nu-fission matrix normalized over outgoing + // groups + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg_; d++) + for (size_t gin = 0; gin < n_g_; gin++) + for (size_t gout = 0; gout < n_g_; gout++) + chi_delayed(a, d, gin, gout) = + temp_matrix_d(a, d, gin, gout) / delayed_nu_fission(a, d, gin); } void XsData::fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang) @@ -338,16 +416,19 @@ void XsData::fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang) // Therefore, the code only considers the data as prompt. // Get nu-fission matrix - xt::xtensor temp_matrix({n_ang, n_g_, n_g_}, 0.); - read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true); + tensor::Tensor temp_matrix = + tensor::zeros({n_ang, n_g_, n_g_}); + read_nd_tensor(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()); + prompt_nu_fission = temp_matrix.sum(2); + + // chi_prompt is the nu-fission matrix normalized over outgoing groups + for (size_t a = 0; a < n_ang; a++) + for (size_t gin = 0; gin < n_g_; gin++) + for (size_t gout = 0; gout < n_g_; gout++) + chi_prompt(a, gin, gout) = + temp_matrix(a, gin, gout) / prompt_nu_fission(a, gin); } //============================================================================== @@ -356,8 +437,8 @@ void XsData::fission_from_hdf5( hid_t xsdata_grp, size_t n_ang, bool is_isotropic) { // 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); + read_nd_tensor(xsdata_grp, "fission", fission); + read_nd_tensor(xsdata_grp, "kappa-fission", kappa_fission); // 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 @@ -388,7 +469,7 @@ void XsData::fission_from_hdf5( if (n_dg_ == 0) { nu_fission = prompt_nu_fission; } else { - nu_fission = prompt_nu_fission + xt::sum(delayed_nu_fission, {1}); + nu_fission = prompt_nu_fission + delayed_nu_fission.sum(1); } } @@ -404,10 +485,10 @@ void XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, hid_t scatt_grp = open_group(xsdata_grp, "scatter_data"); // Get the outgoing group boundary indices - xt::xtensor gmin({n_ang, n_g_}, 0.); - read_nd_vector(scatt_grp, "g_min", gmin, true); - xt::xtensor gmax({n_ang, n_g_}, 0.); - read_nd_vector(scatt_grp, "g_max", gmax, true); + tensor::Tensor gmin = tensor::zeros({n_ang, n_g_}); + read_nd_tensor(scatt_grp, "g_min", gmin, true); + tensor::Tensor gmax = tensor::zeros({n_ang, n_g_}); + read_nd_tensor(scatt_grp, "g_max", gmax, true); // Make gmin and gmax start from 0 vice 1 as they do in the library gmin -= 1; @@ -415,11 +496,11 @@ void XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, // Now use this info to find the length of a vector to hold the flattened // data. - size_t length = order_data * xt::sum(gmax - gmin + 1)(); + size_t length = order_data * (gmax - gmin + 1).sum(); double_4dvec input_scatt(n_ang, double_3dvec(n_g_)); - xt::xtensor temp_arr({length}, 0.); - read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true); + tensor::Tensor temp_arr = tensor::zeros({length}); + read_nd_tensor(scatt_grp, "scatter_matrix", temp_arr, true); // Compare the number of orders given with the max order of the problem; // strip off the superfluous orders if needed @@ -451,7 +532,7 @@ void XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, double_3dvec temp_mult(n_ang, double_2dvec(n_g_)); if (object_exists(scatt_grp, "multiplicity_matrix")) { temp_arr.resize({length / order_data}); - read_nd_vector(scatt_grp, "multiplicity_matrix", temp_arr); + read_nd_tensor(scatt_grp, "multiplicity_matrix", temp_arr); // convert the flat temp_arr to a jagged array for passing to scatt data size_t temp_idx = 0; @@ -481,8 +562,8 @@ void XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, final_scatter_format == AngleDistributionType::TABULAR) { for (size_t a = 0; a < n_ang; a++) { ScattDataLegendre legendre_scatt; - xt::xtensor in_gmin = xt::view(gmin, a, xt::all()); - xt::xtensor in_gmax = xt::view(gmax, a, xt::all()); + tensor::Tensor in_gmin(gmin.slice(a)); + tensor::Tensor in_gmax(gmax.slice(a)); legendre_scatt.init(in_gmin, in_gmax, temp_mult[a], input_scatt[a]); @@ -496,8 +577,8 @@ void XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, // We are sticking with the current representation // Initialize the ScattData object with this data 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()); + tensor::Tensor in_gmin(gmin.slice(a)); + tensor::Tensor in_gmax(gmax.slice(a)); scatter[a]->init(in_gmin, in_gmax, temp_mult[a], input_scatt[a]); } } @@ -519,33 +600,67 @@ void XsData::combine( if (i == 0) { inverse_velocity = that->inverse_velocity; } - if (that->prompt_nu_fission.shape()[0] > 0) { + if (!that->prompt_nu_fission.empty()) { 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 * - xt::view(xt::sum(that->prompt_nu_fission, {1}), xt::all(), - xt::newaxis(), xt::newaxis()) * - that->chi_prompt; - chi_delayed += scalar * - xt::view(xt::sum(that->delayed_nu_fission, {2}), xt::all(), - xt::all(), xt::newaxis(), xt::newaxis()) * - that->chi_delayed; + // Accumulate chi_prompt weighted by total prompt nu-fission + // (summed over energy groups) for this constituent + { + auto pnf_sum = that->prompt_nu_fission.sum(1); + size_t n_ang = chi_prompt.shape(0); + size_t n_g = chi_prompt.shape(1); + for (size_t a = 0; a < n_ang; a++) + for (size_t gin = 0; gin < n_g; gin++) + for (size_t gout = 0; gout < n_g; gout++) + chi_prompt(a, gin, gout) += + scalar * pnf_sum(a) * that->chi_prompt(a, gin, gout); + } + // Accumulate chi_delayed weighted by total delayed nu-fission + // (summed over energy groups) for this constituent + { + auto dnf_sum = that->delayed_nu_fission.sum(2); + size_t n_ang = chi_delayed.shape(0); + size_t n_dg = chi_delayed.shape(1); + size_t n_g = chi_delayed.shape(2); + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg; d++) + for (size_t gin = 0; gin < n_g; gin++) + for (size_t gout = 0; gout < n_g; gout++) + chi_delayed(a, d, gin, gout) += + scalar * dnf_sum(a, d) * that->chi_delayed(a, d, gin, gout); + } } decay_rate += scalar * that->decay_rate; } - // Ensure the chi_prompt and chi_delayed are normalized to 1 for each - // azimuthal angle and delayed group (for chi_delayed) - 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()); + // Normalize chi_prompt so it sums to 1 over outgoing groups + { + size_t n_ang = chi_prompt.shape(0); + size_t n_g = chi_prompt.shape(1); + for (size_t a = 0; a < n_ang; a++) + for (size_t gin = 0; gin < n_g; gin++) { + tensor::View row = chi_prompt.slice(a, gin); + row /= row.sum(); + } + } + // Normalize chi_delayed so it sums to 1 over outgoing groups + { + size_t n_ang = chi_delayed.shape(0); + size_t n_dg = chi_delayed.shape(1); + size_t n_g = chi_delayed.shape(2); + for (size_t a = 0; a < n_ang; a++) + for (size_t d = 0; d < n_dg; d++) + for (size_t gin = 0; gin < n_g; gin++) { + tensor::View row = chi_delayed.slice(a, d, gin); + row /= row.sum(); + } + } // Allow the ScattData object to combine itself - for (size_t a = 0; a < total.shape()[0]; a++) { + for (size_t a = 0; a < total.shape(0); a++) { // Build vector of the scattering objects to incorporate vector those_scatts(those_xs.size()); for (size_t i = 0; i < those_xs.size(); i++) { diff --git a/tests/cpp_unit_tests/CMakeLists.txt b/tests/cpp_unit_tests/CMakeLists.txt index 20341b420f8..ce7e539ea57 100644 --- a/tests/cpp_unit_tests/CMakeLists.txt +++ b/tests/cpp_unit_tests/CMakeLists.txt @@ -7,6 +7,7 @@ set(TEST_NAMES test_mcpl_stat_sum test_mesh test_region + test_tensor # Add additional unit test files here ) diff --git a/tests/cpp_unit_tests/test_tensor.cpp b/tests/cpp_unit_tests/test_tensor.cpp new file mode 100644 index 00000000000..6eae1cea6a4 --- /dev/null +++ b/tests/cpp_unit_tests/test_tensor.cpp @@ -0,0 +1,987 @@ +#include +#include + +#include +#include + +#include "openmc/tensor.h" + +using namespace openmc; +using namespace openmc::tensor; + +// ============================================================================ +// Tensor constructors +// ============================================================================ + +TEST_CASE("Tensor default constructor") +{ + Tensor t; + REQUIRE(t.size() == 0); + REQUIRE(t.empty()); + REQUIRE(t.shape().empty()); +} + +TEST_CASE("Tensor shape constructor") +{ + Tensor t1({5}); + REQUIRE(t1.size() == 5); + REQUIRE(t1.shape().size() == 1); + REQUIRE(t1.shape(0) == 5); + + Tensor t2({3, 4}); + REQUIRE(t2.size() == 12); + REQUIRE(t2.shape().size() == 2); + REQUIRE(t2.shape(0) == 3); + REQUIRE(t2.shape(1) == 4); + + Tensor t3({2, 3, 4}); + REQUIRE(t3.size() == 24); + REQUIRE(t3.shape().size() == 3); +} + +TEST_CASE("Tensor shape + fill constructor") +{ + Tensor t({2, 3}, 7.0); + REQUIRE(t.size() == 6); + for (size_t i = 0; i < t.size(); ++i) + REQUIRE(t[i] == 7.0); +} + +TEST_CASE("Tensor pointer constructor") +{ + double vals[] = {1.0, 2.0, 3.0, 4.0}; + Tensor t(vals, 4); + REQUIRE(t.size() == 4); + REQUIRE(t.shape(0) == 4); + REQUIRE(t[0] == 1.0); + REQUIRE(t[1] == 2.0); + REQUIRE(t[2] == 3.0); + REQUIRE(t[3] == 4.0); +} + +TEST_CASE("Tensor copy and move") +{ + Tensor a({2, 3}, 5.0); + Tensor b(a); + REQUIRE(b.size() == 6); + REQUIRE(b(0, 0) == 5.0); + // Modifying copy doesn't affect original + b(0, 0) = 99.0; + REQUIRE(a(0, 0) == 5.0); + + Tensor c(std::move(b)); + REQUIRE(c(0, 0) == 99.0); + REQUIRE(c.size() == 6); +} + +// ============================================================================ +// Tensor indexing +// ============================================================================ + +TEST_CASE("Tensor 1D indexing") +{ + Tensor t({4}, 0); + t[0] = 10; + t[1] = 20; + t[2] = 30; + t[3] = 40; + REQUIRE(t(0) == 10); + REQUIRE(t(1) == 20); + REQUIRE(t(2) == 30); + REQUIRE(t(3) == 40); +} + +TEST_CASE("Tensor 2D indexing (row-major)") +{ + // Layout: [[1, 2, 3], [4, 5, 6]] + Tensor t({2, 3}, 0); + int val = 1; + for (size_t i = 0; i < 2; ++i) + for (size_t j = 0; j < 3; ++j) + t(i, j) = val++; + + REQUIRE(t(0, 0) == 1); + REQUIRE(t(0, 2) == 3); + REQUIRE(t(1, 0) == 4); + REQUIRE(t(1, 2) == 6); + // Flat index should match row-major order + REQUIRE(t[0] == 1); + REQUIRE(t[3] == 4); + REQUIRE(t[5] == 6); +} + +TEST_CASE("Tensor 3D indexing") +{ + // 2x3x4 tensor + Tensor t({2, 3, 4}, 0); + t(1, 2, 3) = 42; + // Flat index: 1*12 + 2*4 + 3 = 23 + REQUIRE(t[23] == 42); + REQUIRE(t(1, 2, 3) == 42); +} + +// ============================================================================ +// Tensor assignment +// ============================================================================ + +TEST_CASE("Tensor initializer_list assignment") +{ + Tensor t; + t = {1.0, 2.0, 3.0}; + REQUIRE(t.size() == 3); + REQUIRE(t.shape(0) == 3); + REQUIRE(t[0] == 1.0); + REQUIRE(t[2] == 3.0); +} + +// ============================================================================ +// Tensor mutation +// ============================================================================ + +TEST_CASE("Tensor resize") +{ + Tensor t({2, 3}, 1.0); + REQUIRE(t.size() == 6); + t.resize({4, 5}); + REQUIRE(t.size() == 20); + REQUIRE(t.shape(0) == 4); + REQUIRE(t.shape(1) == 5); +} + +TEST_CASE("Tensor reshape") +{ + Tensor t({12}, 0); + for (size_t i = 0; i < 12; ++i) + t[i] = static_cast(i); + + t.reshape({3, 4}); + REQUIRE(t.shape(0) == 3); + REQUIRE(t.shape(1) == 4); + REQUIRE(t.size() == 12); + // Data unchanged, just reinterpreted + REQUIRE(t(0, 0) == 0); + REQUIRE(t(1, 0) == 4); // row 1, col 0 = flat index 4 + REQUIRE(t(2, 3) == 11); // row 2, col 3 = flat index 11 +} + +TEST_CASE("Tensor fill") +{ + Tensor t({3, 3}, 0.0); + t.fill(42.0); + for (size_t i = 0; i < t.size(); ++i) + REQUIRE(t[i] == 42.0); +} + +// ============================================================================ +// Tensor iterators +// ============================================================================ + +TEST_CASE("Tensor iterators") +{ + Tensor t({4}, 0); + t = {10, 20, 30, 40}; + int sum = 0; + for (auto val : t) + sum += val; + REQUIRE(sum == 100); +} + +// ============================================================================ +// Tensor reductions +// ============================================================================ + +TEST_CASE("Tensor sum (full)") +{ + Tensor t({3}, 0.0); + t = {1.0, 2.0, 3.0}; + REQUIRE(t.sum() == 6.0); +} + +TEST_CASE("Tensor sum (axis) on 2D") +{ + // [[1, 2, 3], + // [4, 5, 6]] + Tensor t({2, 3}, 0); + int v = 1; + for (size_t i = 0; i < 2; ++i) + for (size_t j = 0; j < 3; ++j) + t(i, j) = v++; + + // Sum along axis 0 -> [5, 7, 9] + Tensor s0 = t.sum(0); + REQUIRE(s0.size() == 3); + REQUIRE(s0[0] == 5); + REQUIRE(s0[1] == 7); + REQUIRE(s0[2] == 9); + + // Sum along axis 1 -> [6, 15] + Tensor s1 = t.sum(1); + REQUIRE(s1.size() == 2); + REQUIRE(s1[0] == 6); + REQUIRE(s1[1] == 15); +} + +TEST_CASE("Tensor sum (axis) on 3D") +{ + // 2x3x2 tensor filled with sequential values 1..12 + Tensor t({2, 3, 2}, 0); + int v = 1; + for (size_t i = 0; i < 2; ++i) + for (size_t j = 0; j < 3; ++j) + for (size_t k = 0; k < 2; ++k) + t(i, j, k) = v++; + + // Sum along axis 1 (middle) -> 2x2, each sums 3 values + // [0,0]: t(0,0,0)+t(0,1,0)+t(0,2,0) = 1+3+5 = 9 + // [0,1]: t(0,0,1)+t(0,1,1)+t(0,2,1) = 2+4+6 = 12 + // [1,0]: t(1,0,0)+t(1,1,0)+t(1,2,0) = 7+9+11 = 27 + // [1,1]: t(1,0,1)+t(1,1,1)+t(1,2,1) = 8+10+12 = 30 + Tensor s = t.sum(1); + REQUIRE(s.shape(0) == 2); + REQUIRE(s.shape(1) == 2); + REQUIRE(s(0, 0) == 9); + REQUIRE(s(0, 1) == 12); + REQUIRE(s(1, 0) == 27); + REQUIRE(s(1, 1) == 30); +} + +TEST_CASE("Tensor prod") +{ + Tensor t({4}, 0); + t = {1, 2, 3, 4}; + REQUIRE(t.prod() == 24); +} + +TEST_CASE("Tensor any and all") +{ + Tensor t({4}, false); + REQUIRE(!t.any()); + REQUIRE(!t.all()); + + // Set one element true + t.data()[0] = true; + REQUIRE(t.any()); + REQUIRE(!t.all()); + + // Set all true + for (size_t i = 0; i < t.size(); ++i) + t.data()[i] = true; + REQUIRE(t.any()); + REQUIRE(t.all()); +} + +TEST_CASE("Tensor argmin") +{ + Tensor t({5}, 0.0); + t = {3.0, 1.0, 4.0, 0.5, 2.0}; + REQUIRE(t.argmin() == 3); +} + +TEST_CASE("Tensor flip") +{ + Tensor t({5}, 0); + t = {1, 2, 3, 4, 5}; + Tensor f = t.flip(0); + REQUIRE(f[0] == 5); + REQUIRE(f[1] == 4); + REQUIRE(f[2] == 3); + REQUIRE(f[3] == 2); + REQUIRE(f[4] == 1); +} + +TEST_CASE("Tensor flip 2D") +{ + // [[1, 2], [3, 4], [5, 6]] + Tensor t({3, 2}, 0); + t(0, 0) = 1; + t(0, 1) = 2; + t(1, 0) = 3; + t(1, 1) = 4; + t(2, 0) = 5; + t(2, 1) = 6; + + // Flip axis 0 reverses rows -> [[5,6],[3,4],[1,2]] + Tensor f = t.flip(0); + REQUIRE(f(0, 0) == 5); + REQUIRE(f(0, 1) == 6); + REQUIRE(f(1, 0) == 3); + REQUIRE(f(2, 0) == 1); +} + +// ============================================================================ +// Tensor operators +// ============================================================================ + +TEST_CASE("Tensor scalar compound assignment") +{ + Tensor t({3}, 0.0); + t = {2.0, 4.0, 6.0}; + + t += 1.0; + REQUIRE(t[0] == 3.0); + REQUIRE(t[1] == 5.0); + + t -= 1.0; + REQUIRE(t[0] == 2.0); + + t *= 3.0; + REQUIRE(t[0] == 6.0); + REQUIRE(t[1] == 12.0); + + t /= 2.0; + REQUIRE(t[0] == 3.0); + REQUIRE(t[1] == 6.0); +} + +TEST_CASE("Tensor element-wise arithmetic") +{ + Tensor a({3}, 0.0); + Tensor b({3}, 0.0); + a = {1.0, 2.0, 3.0}; + b = {4.0, 5.0, 6.0}; + + Tensor c = a + b; + REQUIRE(c[0] == 5.0); + REQUIRE(c[1] == 7.0); + REQUIRE(c[2] == 9.0); + + c = a - b; + REQUIRE(c[0] == -3.0); + + c = a / b; + REQUIRE(c[0] == 0.25); +} + +TEST_CASE("Tensor scalar arithmetic") +{ + Tensor a({3}, 0.0); + a = {1.0, 2.0, 3.0}; + + Tensor b = a + 10.0; + REQUIRE(b[0] == 11.0); + REQUIRE(b[2] == 13.0); + + b = a - 1.0; + REQUIRE(b[0] == 0.0); + + b = a * 2.0; + REQUIRE(b[0] == 2.0); + REQUIRE(b[2] == 6.0); + + // Non-member scalar * tensor (commutativity) + b = 2.0 * a; + REQUIRE(b[0] == 2.0); + REQUIRE(b[2] == 6.0); + + // Non-member scalar + tensor + b = 10.0 + a; + REQUIRE(b[0] == 11.0); +} + +TEST_CASE("Tensor compound addition with tensor") +{ + Tensor a({3}, 0.0); + Tensor b({3}, 0.0); + a = {1.0, 2.0, 3.0}; + b = {10.0, 20.0, 30.0}; + a += b; + REQUIRE(a[0] == 11.0); + REQUIRE(a[1] == 22.0); + REQUIRE(a[2] == 33.0); +} + +TEST_CASE("Tensor comparison operators") +{ + Tensor t({4}, 0.0); + t = {1.0, 2.0, 3.0, 4.0}; + + Tensor r = t < 3.0; + REQUIRE(r.data()[0] == true); + REQUIRE(r.data()[1] == true); + REQUIRE(r.data()[2] == false); + REQUIRE(r.data()[3] == false); + + r = t >= 3.0; + REQUIRE(r.data()[0] == false); + REQUIRE(r.data()[2] == true); + REQUIRE(r.data()[3] == true); + + r = t <= 2.0; + REQUIRE(r.data()[0] == true); + REQUIRE(r.data()[1] == true); + REQUIRE(r.data()[2] == false); + + r = t > 3.0; + REQUIRE(r.data()[0] == false); + REQUIRE(r.data()[3] == true); +} + +TEST_CASE("Tensor element-wise comparison") +{ + Tensor a({3}, 0.0); + Tensor b({3}, 0.0); + a = {1.0, 5.0, 3.0}; + b = {2.0, 4.0, 3.0}; + + Tensor r = a < b; + REQUIRE(r.data()[0] == true); + REQUIRE(r.data()[1] == false); + REQUIRE(r.data()[2] == false); +} + +TEST_CASE("Tensor mixed-type multiply") +{ + Tensor a({3}, 0); + Tensor b({3}, 0.0); + a = {2, 3, 4}; + b = {1.5, 2.5, 3.5}; + + Tensor c = a * b; + REQUIRE(c[0] == 3.0); + REQUIRE(c[1] == 7.5); + REQUIRE(c[2] == 14.0); +} + +TEST_CASE("Tensor mixed-type divide") +{ + Tensor a({3}, 0.0); + Tensor b({3}, 0); + a = {10.0, 20.0, 30.0}; + b = {2, 4, 5}; + + Tensor c = a / b; + REQUIRE(c[0] == 5.0); + REQUIRE(c[1] == 5.0); + REQUIRE(c[2] == 6.0); +} + +// ============================================================================ +// Tensor bool specialization +// ============================================================================ + +TEST_CASE("Tensor storage") +{ + // Tensor uses unsigned char internally to avoid std::vector proxy + Tensor t({4}, false); + t.data()[0] = true; + t.data()[2] = true; + REQUIRE(t.any()); + REQUIRE(!t.all()); + REQUIRE(t.data()[0] == true); + REQUIRE(t.data()[1] == false); +} + +// ============================================================================ +// View (via Tensor accessors) +// ============================================================================ + +TEST_CASE("Tensor slice axis 0 (2D)") +{ + // [[1, 2, 3], [4, 5, 6]] + Tensor t({2, 3}, 0); + int v = 1; + for (size_t i = 0; i < 2; ++i) + for (size_t j = 0; j < 3; ++j) + t(i, j) = v++; + + auto r0 = t.slice(0); + REQUIRE(r0.size() == 3); + REQUIRE(r0[0] == 1); + REQUIRE(r0[1] == 2); + REQUIRE(r0[2] == 3); + + auto r1 = t.slice(1); + REQUIRE(r1[0] == 4); + REQUIRE(r1[1] == 5); + REQUIRE(r1[2] == 6); + + // Writing through view modifies the tensor + r0[1] = 99; + REQUIRE(t(0, 1) == 99); +} + +TEST_CASE("Tensor slice axis 1 (2D)") +{ + // [[1, 2], [3, 4], [5, 6]] + Tensor t({3, 2}, 0); + t(0, 0) = 1; + t(0, 1) = 2; + t(1, 0) = 3; + t(1, 1) = 4; + t(2, 0) = 5; + t(2, 1) = 6; + + auto c0 = t.slice(all, 0); + REQUIRE(c0.size() == 3); + REQUIRE(c0[0] == 1); + REQUIRE(c0[1] == 3); + REQUIRE(c0[2] == 5); + + auto c1 = t.slice(all, 1); + REQUIRE(c1[0] == 2); + REQUIRE(c1[1] == 4); + REQUIRE(c1[2] == 6); + + // Write through column view + c1[0] = 77; + REQUIRE(t(0, 1) == 77); +} + +TEST_CASE("Tensor slice with range") +{ + Tensor t({6}, 0); + t = {10, 20, 30, 40, 50, 60}; + + // range(start, end) + auto s = t.slice(range(1, 4)); + REQUIRE(s.size() == 3); + REQUIRE(s[0] == 20); + REQUIRE(s[1] == 30); + REQUIRE(s[2] == 40); + + // range(end) from start — range(3) means [0, 3) + auto s2 = t.slice(range(3)); + REQUIRE(s2.size() == 3); + REQUIRE(s2[0] == 10); + REQUIRE(s2[2] == 30); + + // range(start, SIZE_MAX) to end + auto s3 = t.slice(range(3, 6)); + REQUIRE(s3.size() == 3); + REQUIRE(s3[0] == 40); + REQUIRE(s3[2] == 60); + + // Write through slice + s[0] = 99; + REQUIRE(t[1] == 99); +} + +TEST_CASE("Tensor flat view") +{ + Tensor t({2, 3}, 0); + int v = 1; + for (size_t i = 0; i < 2; ++i) + for (size_t j = 0; j < 3; ++j) + t(i, j) = v++; + + auto f = t.flat(); + REQUIRE(f.size() == 6); + REQUIRE(f[0] == 1); + REQUIRE(f[5] == 6); +} + +TEST_CASE("Tensor slice on 3D") +{ + // 2x3x4 tensor + Tensor t({2, 3, 4}, 0); + int v = 0; + for (size_t i = 0; i < 2; ++i) + for (size_t j = 0; j < 3; ++j) + for (size_t k = 0; k < 4; ++k) + t(i, j, k) = v++; + + // slice(1) -> fix axis 0 at 1 -> 3x4 view + auto s = t.slice(1); + REQUIRE(s.size() == 12); + // t(1,0,0) = 12, t(1,0,1) = 13, ... + REQUIRE(s(0, 0) == 12); + REQUIRE(s(0, 1) == 13); + REQUIRE(s(2, 3) == 23); + + // slice(all, 2) -> fix axis 1 at 2 -> 2x4 view + auto s2 = t.slice(all, 2); + REQUIRE(s2.size() == 8); + // t(0,2,0)=8, t(0,2,1)=9, t(1,2,0)=20 + REQUIRE(s2(0, 0) == 8); + REQUIRE(s2(0, 1) == 9); + REQUIRE(s2(1, 0) == 20); +} + +TEST_CASE("Tensor multi-axis slice") +{ + // 2x3x4 tensor with sequential values + Tensor t({2, 3, 4}, 0); + int v = 0; + for (size_t i = 0; i < 2; ++i) + for (size_t j = 0; j < 3; ++j) + for (size_t k = 0; k < 4; ++k) + t(i, j, k) = v++; + + // slice(1, 2) -> fix axes 0 and 1 -> 1D view of 4 elements + // Equivalent to numpy t[1, 2, :] -> t(1,2,0..3) = [20, 21, 22, 23] + auto s = t.slice(1, 2); + REQUIRE(s.size() == 4); + REQUIRE(s[0] == 20); + REQUIRE(s[1] == 21); + REQUIRE(s[3] == 23); + + // slice(all, 1, range(1, 3)) -> keep axis 0, fix axis 1, range on axis 2 + // Equivalent to numpy t[:, 1, 1:3] + // t(0,1,1)=5, t(0,1,2)=6, t(1,1,1)=17, t(1,1,2)=18 + auto s2 = t.slice(all, 1, range(1, 3)); + REQUIRE(s2.size() == 4); + REQUIRE(s2.ndim() == 2); + REQUIRE(s2(0, 0) == 5); + REQUIRE(s2(0, 1) == 6); + REQUIRE(s2(1, 0) == 17); + REQUIRE(s2(1, 1) == 18); + + // slice(0, range(0, 2)) -> fix axis 0 at 0, range on axis 1 + // Equivalent to numpy t[0, 0:2, :] -> shape (2, 4) + auto s3 = t.slice(0, range(0, 2)); + REQUIRE(s3.ndim() == 2); + REQUIRE(s3.shape(0) == 2); + REQUIRE(s3.shape(1) == 4); + REQUIRE(s3(0, 0) == 0); // t(0,0,0) + REQUIRE(s3(1, 3) == 7); // t(0,1,3) +} + +// ============================================================================ +// View assignment and arithmetic +// ============================================================================ + +TEST_CASE("View scalar assignment (fill)") +{ + Tensor t({2, 3}, 0.0); + auto r = t.slice(0); + r = 7.0; + REQUIRE(t(0, 0) == 7.0); + REQUIRE(t(0, 1) == 7.0); + REQUIRE(t(0, 2) == 7.0); + REQUIRE(t(1, 0) == 0.0); // Other row unchanged +} + +TEST_CASE("View initializer_list assignment") +{ + Tensor t({2, 3}, 0.0); + auto r = t.slice(1); + r = {10.0, 20.0, 30.0}; + REQUIRE(t(1, 0) == 10.0); + REQUIRE(t(1, 1) == 20.0); + REQUIRE(t(1, 2) == 30.0); +} + +TEST_CASE("View copy assignment (deep copy)") +{ + Tensor t({2, 3}, 0.0); + t.slice(0) = {1.0, 2.0, 3.0}; + t.slice(1) = {4.0, 5.0, 6.0}; + + // Copy row 0 into row 1 + t.slice(1) = t.slice(0); + REQUIRE(t(1, 0) == 1.0); + REQUIRE(t(1, 1) == 2.0); + REQUIRE(t(1, 2) == 3.0); +} + +TEST_CASE("View compound operators") +{ + Tensor t({2, 3}, 0.0); + t.slice(0) = {1.0, 2.0, 3.0}; + + t.slice(0) *= 2.0; + REQUIRE(t(0, 0) == 2.0); + REQUIRE(t(0, 1) == 4.0); + + t.slice(0) /= 2.0; + REQUIRE(t(0, 0) == 1.0); + REQUIRE(t(0, 1) == 2.0); +} + +TEST_CASE("View assignment from tensor") +{ + Tensor t({2, 3}, 0.0); + Tensor vals({3}, 0.0); + vals = {7.0, 8.0, 9.0}; + + t.slice(1) = vals; + REQUIRE(t(1, 0) == 7.0); + REQUIRE(t(1, 1) == 8.0); + REQUIRE(t(1, 2) == 9.0); +} + +TEST_CASE("View compound addition from tensor") +{ + Tensor t({2, 3}, 0.0); + t.slice(0) = {1.0, 2.0, 3.0}; + Tensor vals({3}, 0.0); + vals = {10.0, 20.0, 30.0}; + + t.slice(0) += vals; + REQUIRE(t(0, 0) == 11.0); + REQUIRE(t(0, 1) == 22.0); + REQUIRE(t(0, 2) == 33.0); +} + +TEST_CASE("View sum") +{ + Tensor t({2, 3}, 0.0); + t.slice(0) = {1.0, 2.0, 3.0}; + t.slice(1) = {4.0, 5.0, 6.0}; + + REQUIRE(t.slice(0).sum() == 6.0); + REQUIRE(t.slice(1).sum() == 15.0); +} + +TEST_CASE("View iteration") +{ + Tensor t({2, 3}, 0); + t.slice(0) = {1, 2, 3}; + + int sum = 0; + for (auto val : t.slice(0)) + sum += val; + REQUIRE(sum == 6); +} + +TEST_CASE("View sub-slice") +{ + Tensor t({6}, 0); + t = {10, 20, 30, 40, 50, 60}; + + auto s = t.slice(range(1, 5)); // [20, 30, 40, 50] + auto ss = s.slice(range(1, 3)); // [30, 40] + REQUIRE(ss.size() == 2); + REQUIRE(ss[0] == 30); + REQUIRE(ss[1] == 40); +} + +TEST_CASE("Tensor from View") +{ + Tensor t({2, 3}, 0.0); + t.slice(0) = {1.0, 2.0, 3.0}; + + // Construct a new tensor from a view (copies data) + Tensor t2(t.slice(0)); + REQUIRE(t2.size() == 3); + REQUIRE(t2[0] == 1.0); + REQUIRE(t2[2] == 3.0); + + // Modifying the new tensor doesn't affect the original + t2[0] = 99.0; + REQUIRE(t(0, 0) == 1.0); +} + +// ============================================================================ +// Const View +// ============================================================================ + +TEST_CASE("Const tensor produces const views") +{ + Tensor t({2, 3}, 0.0); + int v = 1; + for (size_t i = 0; i < 2; ++i) + for (size_t j = 0; j < 3; ++j) + t(i, j) = v++; + + const Tensor& ct = t; + auto r = ct.slice(0); // View + REQUIRE(r[0] == 1.0); + REQUIRE(r[2] == 3.0); + + auto c = ct.slice(all, 1); + REQUIRE(c[0] == 2.0); + REQUIRE(c[1] == 5.0); +} + +// ============================================================================ +// StaticTensor2D +// ============================================================================ + +TEST_CASE("StaticTensor2D basics") +{ + StaticTensor2D t; + REQUIRE(t.size() == 12); + REQUIRE(t.shape()[0] == 3); + REQUIRE(t.shape()[1] == 4); + + // Default-initialized to zero + REQUIRE(t(0, 0) == 0.0); + + t(1, 2) = 42.0; + REQUIRE(t(1, 2) == 42.0); + // Flat data: row 1, col 2 = index 1*4 + 2 = 6 + REQUIRE(t.data()[6] == 42.0); +} + +TEST_CASE("StaticTensor2D fill") +{ + StaticTensor2D t; + t.fill(5); + for (size_t i = 0; i < t.size(); ++i) + REQUIRE(t.data()[i] == 5); +} + +TEST_CASE("StaticTensor2D iteration") +{ + StaticTensor2D t; + t.fill(1); + int sum = 0; + for (auto val : t) + sum += val; + REQUIRE(sum == 6); +} + +TEST_CASE("StaticTensor2D slice") +{ + StaticTensor2D t; + t(0, 0) = 1; + t(0, 1) = 2; + t(1, 0) = 3; + t(1, 1) = 4; + t(2, 0) = 5; + t(2, 1) = 6; + + // slice(1) = row 1 (fix axis 0 at 1) + auto r1 = t.slice(1); + REQUIRE(r1.size() == 2); + REQUIRE(r1[0] == 3); + REQUIRE(r1[1] == 4); + + // slice(all, 0) = column 0 (fix axis 1 at 0) + auto c0 = t.slice(all, 0); + REQUIRE(c0.size() == 3); + REQUIRE(c0[0] == 1); + REQUIRE(c0[1] == 3); + REQUIRE(c0[2] == 5); +} + +TEST_CASE("StaticTensor2D flat view") +{ + StaticTensor2D t; + t(0, 0) = 1.0; + t(0, 1) = 2.0; + t(1, 0) = 3.0; + t(1, 1) = 4.0; + + auto f = t.flat(); + REQUIRE(f.size() == 4); + f = 0.0; + REQUIRE(t(0, 0) == 0.0); + REQUIRE(t(1, 1) == 0.0); +} + +// ============================================================================ +// Non-member functions +// ============================================================================ + +TEST_CASE("zeros") +{ + auto t = zeros({3, 4}); + REQUIRE(t.size() == 12); + for (size_t i = 0; i < t.size(); ++i) + REQUIRE(t[i] == 0.0); +} + +TEST_CASE("zeros_like") +{ + Tensor a({2, 5}, 7.0); + auto b = zeros_like(a); + REQUIRE(b.size() == 10); + REQUIRE(b.shape(0) == 2); + REQUIRE(b.shape(1) == 5); + for (size_t i = 0; i < b.size(); ++i) + REQUIRE(b[i] == 0.0); +} + +TEST_CASE("full_like") +{ + Tensor a({4}, 0); + auto b = full_like(a, 42); + REQUIRE(b.size() == 4); + for (size_t i = 0; i < b.size(); ++i) + REQUIRE(b[i] == 42); +} + +TEST_CASE("linspace") +{ + auto t = linspace(0.0, 1.0, 5); + REQUIRE(t.size() == 5); + REQUIRE(t[0] == 0.0); + REQUIRE(t[4] == 1.0); + REQUIRE_THAT(t[1], Catch::Matchers::WithinRel(0.25, 1e-12)); + REQUIRE_THAT(t[2], Catch::Matchers::WithinRel(0.5, 1e-12)); + REQUIRE_THAT(t[3], Catch::Matchers::WithinRel(0.75, 1e-12)); +} + +TEST_CASE("concatenate") +{ + Tensor a({3}, 0); + Tensor b({2}, 0); + a = {1, 2, 3}; + b = {4, 5}; + + auto c = concatenate(a, b); + REQUIRE(c.size() == 5); + REQUIRE(c[0] == 1); + REQUIRE(c[2] == 3); + REQUIRE(c[3] == 4); + REQUIRE(c[4] == 5); +} + +TEST_CASE("log") +{ + Tensor t({3}, 0.0); + t = {1.0, std::exp(1.0), std::exp(2.0)}; + + auto r = log(t); + REQUIRE_THAT(r[0], Catch::Matchers::WithinAbs(0.0, 1e-12)); + REQUIRE_THAT(r[1], Catch::Matchers::WithinAbs(1.0, 1e-12)); + REQUIRE_THAT(r[2], Catch::Matchers::WithinAbs(2.0, 1e-12)); +} + +TEST_CASE("abs") +{ + Tensor t({4}, 0.0); + t = {-3.0, -1.0, 0.0, 2.0}; + + auto r = abs(t); + REQUIRE(r[0] == 3.0); + REQUIRE(r[1] == 1.0); + REQUIRE(r[2] == 0.0); + REQUIRE(r[3] == 2.0); +} + +TEST_CASE("where") +{ + Tensor cond({4}, false); + cond.data()[0] = true; + cond.data()[2] = true; + + Tensor vals({4}, 0.0); + vals = {10.0, 20.0, 30.0, 40.0}; + + auto r = where(cond, vals, -1.0); + REQUIRE(r[0] == 10.0); + REQUIRE(r[1] == -1.0); + REQUIRE(r[2] == 30.0); + REQUIRE(r[3] == -1.0); +} + +TEST_CASE("nan_to_num") +{ + Tensor t({4}, 0.0); + t[0] = 1.0; + t[1] = std::nan(""); + t[2] = std::numeric_limits::infinity(); + t[3] = -std::numeric_limits::infinity(); + + auto r = nan_to_num(t); + REQUIRE(r[0] == 1.0); + REQUIRE(r[1] == 0.0); // NaN -> 0 + REQUIRE(r[2] == std::numeric_limits::max()); // +inf -> max + REQUIRE(r[3] == std::numeric_limits::lowest()); // -inf -> lowest +} + +// ============================================================================ +// is_tensor trait +// ============================================================================ + +TEST_CASE("is_tensor trait") +{ + REQUIRE(is_tensor>::value); + REQUIRE(is_tensor>::value); + REQUIRE(is_tensor>::value); + REQUIRE(!is_tensor::value); + REQUIRE(!is_tensor>::value); +} diff --git a/tests/regression_tests/cpp_driver/driver.cpp b/tests/regression_tests/cpp_driver/driver.cpp index a99c97b64e4..a6c3e651037 100644 --- a/tests/regression_tests/cpp_driver/driver.cpp +++ b/tests/regression_tests/cpp_driver/driver.cpp @@ -2,6 +2,8 @@ #include #endif +#include + #include "openmc/capi.h" #include "openmc/cell.h" #include "openmc/error.h" diff --git a/vendor/xtensor b/vendor/xtensor deleted file mode 160000 index 3634f2ded19..00000000000 --- a/vendor/xtensor +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 3634f2ded19e0cf38208c8b86cea9e1d7c8e397d diff --git a/vendor/xtl b/vendor/xtl deleted file mode 160000 index a7c1c5444df..00000000000 --- a/vendor/xtl +++ /dev/null @@ -1 +0,0 @@ -Subproject commit a7c1c5444dfc57f76620391af4c94785ff82c8d6