From 5eab6dbed0d2d3d0f564781c51f5c41a6382fc1e Mon Sep 17 00:00:00 2001 From: Adam G Nelson Date: Wed, 5 Dec 2018 21:02:59 -0500 Subject: [PATCH] Got URR working, now just need a few simplifications of the Fortran code --- include/openmc/constants.h | 14 ++-- include/openmc/nuclide.h | 121 ++++++++++++++-------------- include/openmc/urr.h | 7 +- src/nuclide.cpp | 153 +++++++++++++++++++++++++++++++++++- src/nuclide_header.F90 | 156 +++---------------------------------- src/physics.cpp | 2 +- src/urr.cpp | 43 ++++++---- 7 files changed, 262 insertions(+), 234 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 2064252037..dbb5ed3147 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -257,12 +257,12 @@ constexpr int LIBRARY_PHOTON {3}; constexpr int LIBRARY_MULTIGROUP {4}; // Probability table parameters -constexpr int URR_CUM_PROB {1}; -constexpr int URR_TOTAL {2}; -constexpr int URR_ELASTIC {3}; -constexpr int URR_FISSION {4}; -constexpr int URR_N_GAMMA {5}; -constexpr int URR_HEATING {6}; +constexpr int URR_CUM_PROB {0}; +constexpr int URR_TOTAL {1}; +constexpr int URR_ELASTIC {2}; +constexpr int URR_FISSION {3}; +constexpr int URR_N_GAMMA {4}; +constexpr int URR_HEATING {5}; // Maximum number of partial fission reactions constexpr int PARTIAL_FISSION_MAX {4}; @@ -424,7 +424,7 @@ constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE // Interpolation rules enum class Interpolation { - histogram, lin_lin, lin_log, log_lin, log_log + histogram = 1, lin_lin = 2, lin_log = 3, log_lin = 4, log_log = 5 }; // Run modes diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 72e15a02f1..662a6fa7c2 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -24,63 +24,6 @@ namespace openmc { constexpr double CACHE_INVALID {-1.0}; -//=============================================================================== -// Data for a nuclide -//=============================================================================== - -class Nuclide { -public: - // Types, aliases - using EmissionMode = ReactionProduct::EmissionMode; - struct EnergyGrid { - std::vector grid_index; - std::vector energy; - }; - - // Constructors - Nuclide(hid_t group, const double* temperature, int n); - - // Methods - double nu(double E, EmissionMode mode, int group=0) const; - void calculate_elastic_xs(int i_nuclide) const; - - //! Determines the microscopic 0K elastic cross section at a trial relative - //! energy used in resonance scattering - double elastic_xs_0K(double E) const; - - // Data members - std::string name_; //! Name of nuclide, e.g. "U235" - int Z_; //! Atomic number - int A_; //! Mass number - int metastable_; //! Metastable state - double awr_; //! Atomic weight ratio - std::vector kTs_; //! temperatures in eV (k*T) - std::vector grid_; //! Energy grid at each temperature - - bool fissionable_ {false}; //! Whether nuclide is fissionable - bool has_partial_fission_ {false}; //! has partial fission reactions? - std::vector fission_rx_; //! Fission reactions - int n_precursor_ {0}; //! Number of delayed neutron precursors - std::unique_ptr total_nu_; //! Total neutron yield - - // Resonance scattering information - bool resonant_ {false}; - std::vector energy_0K_; - std::vector elastic_0K_; - std::vector xs_cdf_; - - // Unresolved resonance range information - bool urr_present_ {false}; - int urr_inelastic_ {C_NONE}; - std::vector urr_data_; - - std::vector> reactions_; //! Reactions - std::vector index_inelastic_scatter_; - -private: - void create_derived(); -}; - //=============================================================================== //! Cached microscopic cross sections for a particular nuclide at the current //! energy @@ -138,6 +81,68 @@ struct MaterialMacroXS { double pair_production; //!< macroscopic pair production xs }; +//=============================================================================== +// Data for a nuclide +//=============================================================================== + +class Nuclide { +public: + // Types, aliases + using EmissionMode = ReactionProduct::EmissionMode; + struct EnergyGrid { + std::vector grid_index; + std::vector energy; + }; + + // Constructors + Nuclide(hid_t group, const double* temperature, int n, int i_nuclide); + + // Methods + double nu(double E, EmissionMode mode, int group=0) const; + void calculate_elastic_xs() const; + + //! Determines the microscopic 0K elastic cross section at a trial relative + //! energy used in resonance scattering + double elastic_xs_0K(double E) const; + + //! \brief Determines cross sections in the unresolved resonance range + //! from probability tables. + void calculate_urr_xs(const int i_temp, const double E); + + // Data members + std::string name_; //! Name of nuclide, e.g. "U235" + int Z_; //! Atomic number + int A_; //! Mass number + int metastable_; //! Metastable state + double awr_; //! Atomic weight ratio + std::vector kTs_; //! temperatures in eV (k*T) + std::vector grid_; //! Energy grid at each temperature + int i_nuclide_; //! Index in the nuclides array + + bool fissionable_ {false}; //! Whether nuclide is fissionable + bool has_partial_fission_ {false}; //! has partial fission reactions? + std::vector fission_rx_; //! Fission reactions + int n_precursor_ {0}; //! Number of delayed neutron precursors + std::unique_ptr total_nu_; //! Total neutron yield + + // Resonance scattering information + bool resonant_ {false}; + std::vector energy_0K_; + std::vector elastic_0K_; + std::vector xs_cdf_; + + // Unresolved resonance range information + bool urr_present_ {false}; + int urr_inelastic_ {C_NONE}; + std::vector urr_data_; + + std::vector> reactions_; //! Reactions + std::vector index_inelastic_scatter_; + +private: + void create_derived(); +}; + //============================================================================== // Global variables //============================================================================== @@ -170,6 +175,8 @@ extern "C" void set_micro_xs(); extern "C" bool nuclide_wmp_present(int i_nuclide); extern "C" double nuclide_wmp_emin(int i_nuclide); extern "C" double nuclide_wmp_emax(int i_nuclide); +extern "C" void nuclide_calculate_urr_xs(const int i_nuclide, + const int i_temp, const double E); } // namespace openmc diff --git a/include/openmc/urr.h b/include/openmc/urr.h index 224c79a358..acfe4d0590 100644 --- a/include/openmc/urr.h +++ b/include/openmc/urr.h @@ -5,6 +5,7 @@ #include "xtensor/xtensor.hpp" +#include "openmc/constants.h" #include "openmc/hdf5_interface.h" namespace openmc { @@ -15,16 +16,14 @@ namespace openmc { class UrrData{ public: - long unsigned int n_energy_; // # of incident energies - long unsigned int n_prob_; // # of probabilities - int interp_; // interpolation type (2=lin-lin, 5=log-log) + Interpolation interp_; // interpolation type int inelastic_flag_; // inelastic competition flag int absorption_flag_; // other absorption flag bool multiply_smooth_; // multiply by smooth cross section? xt::xtensor energy_; // incident energies xt::xtensor prob_; // Actual probability tables - //! Load the URR data from the provided HDF5 group + //! \brief Load the URR data from the provided HDF5 group void from_hdf5(hid_t group_id); }; diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 2f00e612f3..64c983fd39 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -5,6 +5,7 @@ #include "openmc/error.h" #include "openmc/hdf5_interface.h" #include "openmc/message_passing.h" +#include "openmc/random_lcg.h" #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/string_utils.h" @@ -33,7 +34,7 @@ MaterialMacroXS material_xs; // Nuclide implementation //============================================================================== -Nuclide::Nuclide(hid_t group, const double* temperature, int n) +Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide) { // Get name of nuclide from group, removing leading '/' name_ = object_name(group).substr(1); @@ -42,6 +43,7 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n) read_attribute(group, "A", A_); read_attribute(group, "metastable", metastable_); read_attribute(group, "atomic_weight_ratio", awr_); + i_nuclide_ = i_nuclide; // Determine temperatures available hid_t kT_group = open_group(group, "kTs"); @@ -365,10 +367,10 @@ double Nuclide::nu(double E, EmissionMode mode, int group) const } } -void Nuclide::calculate_elastic_xs(int i_nuclide) const +void Nuclide::calculate_elastic_xs() const { // Get temperature index, grid index, and interpolation factor - auto& micro = simulation::micro_xs[i_nuclide]; + auto& micro = simulation::micro_xs[i_nuclide_]; int i_temp = micro.index_temp - 1; int i_grid = micro.index_grid - 1; double f = micro.interp_factor; @@ -402,6 +404,142 @@ double Nuclide::elastic_xs_0K(double E) const return (1.0 - f)*elastic_0K_[i_grid] + f*elastic_0K_[i_grid + 1]; } +void Nuclide::calculate_urr_xs(const int i_temp, const double E) +{ + auto& micro = simulation::micro_xs[i_nuclide_]; + micro.use_ptable = true; + + // Create a shorthand for the URR data + UrrData* urr = &(urr_data_[i_temp]); + + // Determine the energy table + int i_energy = 0; + while (true) { + if (E < urr->energy_(i_energy + 1)) {break;} + i_energy++; + } + + // Sample the probability table using the cumulative distribution + + // Random nmbers for the xs calculation are sampled from a separate stream. + // This guarantees the rnadomness and, at the same time, makes sure we + // reuse random numbers for the same nuclide at different temperatures, + // therefore preserving correlation of temperature in probability tables. + prn_set_stream(STREAM_URR_PTABLE); + //TODO: to maintain the same random number stream as the Fortran code this + //replaces, the seed is set with i_nuclide_ + 1 instead of i_nuclide_ + double r = future_prn(static_cast(i_nuclide_ + 1)); + prn_set_stream(STREAM_TRACKING); + + int i_low = 0; + while (true) { + if (urr->prob_(i_energy, URR_CUM_PROB, i_low) > r) {break;} + i_low++; + } + int i_up = 0; + while (true) { + if (urr->prob_(i_energy + 1, URR_CUM_PROB, i_up) > r) {break;} + i_up++; + } + + // Determine elastic, fission, and capture cross sections from the + // probability table + double elastic = 0.; + double fission = 0.; + double capture = 0.; + double f; + if (urr->interp_ == Interpolation::lin_lin) { + // Determine the interpolation factor on the table + f = (E - urr->energy_(i_energy)) / + (urr->energy_(i_energy + 1) - urr->energy_(i_energy)); + + elastic = (1. - f) * urr->prob_(i_energy, URR_ELASTIC, i_low) + + f * urr->prob_(i_energy + 1, URR_ELASTIC, i_up); + fission = (1. - f) * urr->prob_(i_energy, URR_FISSION, i_low) + + f * urr->prob_(i_energy + 1, URR_FISSION, i_up); + capture = (1. - f) * urr->prob_(i_energy, URR_N_GAMMA, i_low) + + f * urr->prob_(i_energy + 1, URR_N_GAMMA, i_up); + } else if (urr->interp_ == Interpolation::log_log) { + // Determine interpolation factor on the table + f = std::log(E / urr->energy_(i_energy)) / + std::log(urr->energy_(i_energy + 1) / urr->energy_(i_energy)); + + // Calculate the elastic cross section/factor + if ((urr->prob_(i_energy, URR_ELASTIC, i_low) > 0.) && + (urr->prob_(i_energy + 1, URR_ELASTIC, i_up) > 0.)) { + elastic = + std::exp((1. - f) * + std::log(urr->prob_(i_energy, URR_ELASTIC, i_low)) + + f * std::log(urr->prob_(i_energy + 1, URR_ELASTIC, i_up))); + } else { + elastic = 0.; + } + + // Calculate the fission cross section/factor + if ((urr->prob_(i_energy, URR_FISSION, i_low) > 0.) && + (urr->prob_(i_energy + 1, URR_FISSION, i_up) > 0.)) { + fission = + std::exp((1. - f) * + std::log(urr->prob_(i_energy, URR_FISSION, i_low)) + + f * std::log(urr->prob_(i_energy + 1, URR_FISSION, i_up))); + } else { + fission = 0.; + } + + // Calculate the capture cross section/factor + if ((urr->prob_(i_energy, URR_N_GAMMA, i_low) > 0.) && + (urr->prob_(i_energy + 1, URR_N_GAMMA, i_up) > 0.)) { + capture = + std::exp((1. - f) * + std::log(urr->prob_(i_energy, URR_N_GAMMA, i_low)) + + f * std::log(urr->prob_(i_energy + 1, URR_N_GAMMA, i_up))); + } else { + capture = 0.; + } + } + + // Determine the treatment of inelastic scattering + double inelastic = 0.; + if (urr->inelastic_flag_ != C_NONE) { + // get interpolation factor + f = micro.interp_factor; + + // Determine inelastic scattering cross section + Reaction* rx = reactions_[urr_inelastic_].get(); + int xs_index = micro.index_grid - rx->xs_[i_temp].threshold; + if (xs_index >= 0) { + inelastic = (1. - f) * rx->xs_[i_temp].value[xs_index] + + f * rx->xs_[i_temp].value[xs_index + 1]; + } + } + + // Multiply by smooth cross-section if needed + if (urr->multiply_smooth_) { + calculate_elastic_xs(); + elastic *= micro.elastic; + capture *= (micro.absorption - micro.fission); + fission *= micro.fission; + } + + // Check for negative values + if (elastic < 0.) {elastic = 0.;} + if (fission < 0.) {fission = 0.;} + if (capture < 0.) {capture = 0.;} + + // Set elastic, absorption, fission, and total x/s. Note that the total x/s + // is calculated as a sum of partials instead of the table-provided value + micro.elastic = elastic; + micro.absorption = capture + fission; + micro.fission = fission; + micro.total = elastic + inelastic + capture + fission; + + // Determine nu-fission cross-section + if (fissionable_) { + micro.nu_fission = nu(E, EmissionMode::total) * micro.fission; + } + +} + //============================================================================== // Fortran compatibility functions //============================================================================== @@ -415,7 +553,8 @@ set_particle_energy_bounds(int particle, double E_min, double E_max) extern "C" Nuclide* nuclide_from_hdf5_c(hid_t group, const double* temperature, int n) { - data::nuclides.push_back(std::make_unique(group, temperature, n)); + data::nuclides.push_back(std::make_unique(group, temperature, n, + data::nuclides.size())); return data::nuclides.back().get(); } @@ -436,4 +575,10 @@ void set_micro_xs() } } +extern "C" void +nuclide_calculate_urr_xs(const int i_nuclide, const int i_temp, const double E) +{ + data::nuclides[i_nuclide - 1]->calculate_urr_xs(i_temp - 1, E); +} + } // namespace openmc diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index becb865240..e05a703130 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -198,6 +198,14 @@ module nuclide_header character(kind=C_CHAR), intent(in) :: name(*) type(C_PTR) :: path end function + + subroutine nuclide_calculate_urr_xs_c(i_nuclide, i_temp, E) & + bind(C, name='nuclide_calculate_urr_xs') + import C_INT, C_DOUBLE + integer(C_INT), value, intent(in) :: i_nuclide + integer(C_INT), value, intent(in) :: i_temp + real(C_DOUBLE), value, intent(in) :: E + end subroutine nuclide_calculate_urr_xs_c end interface contains @@ -980,7 +988,7 @@ contains if (urr_ptables_on .and. this % urr_present .and. .not. use_mp) then if (E > this % urr_data(i_temp) % energy(1) .and. E < this % & urr_data(i_temp) % energy(this % urr_data(i_temp) % n_energy)) then - call calculate_urr_xs(this, i_temp, E, micro_xs) + call nuclide_calculate_urr_xs_c(this % i_nuclide, i_temp, E) end if end if @@ -1246,152 +1254,6 @@ contains end if end subroutine multipole_deriv_eval -!=============================================================================== -! CALCULATE_URR_XS determines cross sections in the unresolved resonance range -! from probability tables -!=============================================================================== - - subroutine calculate_urr_xs(this, i_temp, E, micro_xs) - class(Nuclide), intent(in) :: this ! Nuclide object - integer, intent(in) :: i_temp ! temperature index - real(8), intent(in) :: E ! energy - type(NuclideMicroXS), intent(inout) :: micro_xs ! Cross section cache - - integer :: i_energy ! index for energy - integer :: i_low ! band index at lower bounding energy - integer :: i_up ! band index at upper bounding energy - integer :: threshold ! threshold energy index - real(8) :: f ! interpolation factor - real(8) :: r ! pseudo-random number - real(8) :: elastic ! elastic cross section - real(8) :: capture ! (n,gamma) cross section - real(8) :: fission ! fission cross section - real(8) :: inelastic ! inelastic cross section - - micro_xs % use_ptable = .true. - - associate (urr => this % urr_data(i_temp)) - ! determine energy table - i_energy = 1 - do - if (E < urr % energy(i_energy + 1)) exit - i_energy = i_energy + 1 - end do - - ! determine interpolation factor on table - f = (E - urr % energy(i_energy)) / & - (urr % energy(i_energy + 1) - urr % energy(i_energy)) - - ! sample probability table using the cumulative distribution - - ! Random numbers for xs calculation are sampled from a separated stream. - ! This guarantees the randomness and, at the same time, makes sure we reuse - ! random number for the same nuclide at different temperatures, therefore - ! preserving correlation of temperature in probability tables. - call prn_set_stream(STREAM_URR_PTABLE) - r = future_prn(int(this % i_nuclide, 8)) - call prn_set_stream(STREAM_TRACKING) - - i_low = 1 - do - if (urr % prob(i_energy, URR_CUM_PROB, i_low) > r) exit - i_low = i_low + 1 - end do - i_up = 1 - do - if (urr % prob(i_energy + 1, URR_CUM_PROB, i_up) > r) exit - i_up = i_up + 1 - end do - - ! determine elastic, fission, and capture cross sections from probability - ! table - if (urr % interp == LINEAR_LINEAR) then - elastic = (ONE - f) * urr % prob(i_energy, URR_ELASTIC, i_low) + & - f * urr % prob(i_energy + 1, URR_ELASTIC, i_up) - fission = (ONE - f) * urr % prob(i_energy, URR_FISSION, i_low) + & - f * urr % prob(i_energy + 1, URR_FISSION, i_up) - capture = (ONE - f) * urr % prob(i_energy, URR_N_GAMMA, i_low) + & - f * urr % prob(i_energy + 1, URR_N_GAMMA, i_up) - elseif (urr % interp == LOG_LOG) then - ! Get logarithmic interpolation factor - f = log(E / urr % energy(i_energy)) / & - log(urr % energy(i_energy + 1) / urr % energy(i_energy)) - - ! Calculate elastic cross section/factor - elastic = ZERO - if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then - elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, & - i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, & - i_up))) - end if - - ! Calculate fission cross section/factor - fission = ZERO - if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then - fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, & - i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, & - i_up))) - end if - - ! Calculate capture cross section/factor - capture = ZERO - if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then - capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, & - i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, & - i_up))) - end if - end if - - ! Determine treatment of inelastic scattering - inelastic = ZERO - if (urr % inelastic_flag > 0) then - ! Get index on energy grid and interpolation factor - i_energy = micro_xs % index_grid - f = micro_xs % interp_factor - - ! Determine inelastic scattering cross section - associate (rx => this % reactions(this % urr_inelastic)) - threshold = rx % xs_threshold(i_temp) - if (i_energy >= threshold) then - inelastic = (ONE - f) * rx % xs(i_temp, i_energy - threshold + 1) + & - f * rx % xs(i_temp, i_energy - threshold + 2) - end if - end associate - end if - - ! Multiply by smooth cross-section if needed - if (urr % multiply_smooth) then - call this % calculate_elastic_xs(micro_xs) - elastic = elastic * micro_xs % elastic - capture = capture * (micro_xs % absorption - micro_xs % fission) - fission = fission * micro_xs % fission - end if - - ! Check for negative values - if (elastic < ZERO) elastic = ZERO - if (fission < ZERO) fission = ZERO - if (capture < ZERO) capture = ZERO - - ! Set elastic, absorption, fission, and total cross sections. Note that the - ! total cross section is calculated as sum of partials rather than using the - ! table-provided value - micro_xs % elastic = elastic - micro_xs % absorption = capture + fission - micro_xs % fission = fission - micro_xs % total = elastic + inelastic + capture + fission - - ! Determine nu-fission cross section - if (this % fissionable) then - micro_xs % nu_fission = this % nu(E, EMISSION_TOTAL) * & - micro_xs % fission - end if - end associate - - end subroutine calculate_urr_xs - !=============================================================================== ! CHECK_DATA_VERSION checks for the right version of nuclear data within HDF5 ! files diff --git a/src/physics.cpp b/src/physics.cpp index 43fc2871b4..df2f519134 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -627,7 +627,7 @@ void scatter(Particle* p, int i_nuclide, int i_nuc_mat) // Calculate elastic cross section if it wasn't precalculated if (micro.elastic == CACHE_INVALID) { - nuc->calculate_elastic_xs(i_nuclide); + nuc->calculate_elastic_xs(); } double prob = micro.elastic - micro.thermal; diff --git a/src/urr.cpp b/src/urr.cpp index 7a33a20e8f..1ecc8ca8c6 100644 --- a/src/urr.cpp +++ b/src/urr.cpp @@ -1,25 +1,45 @@ #include "openmc/urr.h" +#include + namespace openmc { void UrrData::from_hdf5(hid_t group_id) { // Read interpolation and other flags - read_attribute(group_id, "interpolation", interp_); + int interp_temp; + read_attribute(group_id, "interpolation", interp_temp); + switch (interp_temp) { + case static_cast(Interpolation::histogram): + interp_ = Interpolation::histogram; + break; + case static_cast(Interpolation::lin_lin): + interp_ = Interpolation::lin_lin; + break; + case static_cast(Interpolation::lin_log): + interp_ = Interpolation::lin_log; + break; + case static_cast(Interpolation::log_lin): + interp_ = Interpolation::log_lin; + break; + case static_cast(Interpolation::log_log): + interp_ = Interpolation::log_log; + } + read_attribute(group_id, "inelastic", inelastic_flag_); read_attribute(group_id, "absorption", absorption_flag_); - int i; - read_attribute(group_id, "multiply_smooth", i); - multiply_smooth_ = (i == 1); + int temp_multiply_smooth; + read_attribute(group_id, "multiply_smooth", temp_multiply_smooth); + multiply_smooth_ = (temp_multiply_smooth == 1); // read the enrgies at which tables exist hid_t dset = open_dataset(group_id, "energy"); hsize_t dims[1]; get_shape(dset, dims); close_dataset(dset); - n_energy_ = static_cast(dims[0]); - energy_ = xt::xtensor({n_energy_}, 0.); + // n_energy_ = static_cast(dims[0]); + energy_ = xt::xtensor({dims[0]}, 0.); read_dataset_as_shape(group_id, "energy", energy_); // Read URR tables @@ -27,14 +47,9 @@ UrrData::from_hdf5(hid_t group_id) hsize_t dims3[3]; get_shape(dset, dims3); close_dataset(dset); - n_prob_ = static_cast(dims3[0]); - xt::xarray temp({n_energy_, 6, n_prob_}); - read_dataset(group_id, "table", temp); - - prob_ = xt::xtensor({n_energy_, 6, n_prob_}, 0.); - prob_ = temp; - //TODO: swap 1st and last indices? - + xt::xarray temp_arr({dims3[0], dims3[1], dims3[2]}, 0.); + read_dataset(group_id, "table", temp_arr); + prob_ = temp_arr; } } \ No newline at end of file