diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index a5844c506..e7db9016b 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -1,7 +1,6 @@ #ifndef OPENMC_LATTICE_H #define OPENMC_LATTICE_H -#include #include #include #include @@ -9,6 +8,7 @@ #include "hdf5.h" #include "pugixml.hpp" +#include "openmc/array.h" #include "openmc/constants.h" #include "openmc/memory.h" #include "openmc/position.h" @@ -26,7 +26,6 @@ enum class LatticeType { rect, hex }; - //============================================================================== // Global variables //============================================================================== @@ -60,7 +59,7 @@ public: virtual ~Lattice() {} - virtual int32_t& operator[](std::array i_xyz) = 0; + virtual int32_t const& operator[](array const& i_xyz) = 0; virtual LatticeIter begin(); LatticeIter end(); @@ -83,14 +82,7 @@ public: //! \param i_xyz[3] The indices for a lattice tile. //! \return true if the given indices fit within the lattice bounds. False //! otherwise. - virtual bool are_valid_indices(const int i_xyz[3]) const = 0; - - bool - are_valid_indices(std::array i_xyz) const - { - int i_xyz_[3] {i_xyz[0], i_xyz[1], i_xyz[2]}; - return are_valid_indices(i_xyz_); - } + virtual bool are_valid_indices(array const& i_xyz) const = 0; //! \brief Find the next lattice surface crossing //! \param r A 3D Cartesian coordinate. @@ -103,15 +95,17 @@ public: //! \brief Find the lattice tile indices for a given point. //! \param r A 3D Cartesian coordinate. - //! \return An array containing the indices of a lattice tile. - virtual std::array get_indices(Position r, Direction u) const = 0; + //! \param u Direction of a particle + //! \param result resulting indices to save to + virtual void get_indices( + Position r, Direction u, array& result) const = 0; //! \brief Get coordinates local to a lattice tile. //! \param r A 3D Cartesian coordinate. //! \param i_xyz The indices for a lattice tile. //! \return Local 3D Cartesian coordinates. - virtual Position - get_local_position(Position r, const std::array i_xyz) const = 0; + virtual Position get_local_position( + Position r, const array& i_xyz) const = 0; //! \brief Check flattened lattice index. //! \param indx The index for a lattice tile. @@ -125,7 +119,7 @@ public: //! \param i_xyz[3] The indices for a lattice tile. //! \return Distribcell offset i.e. the largest instance number for the target //! cell found in the geometry tree under this lattice tile. - virtual int32_t& offset(int map, const int i_xyz[3]) = 0; + virtual int32_t& offset(int map, array const& i_xyz) = 0; //! \brief Get the distribcell offset for a lattice tile. //! \param The map index for the target cell. @@ -211,19 +205,18 @@ class RectLattice : public Lattice public: explicit RectLattice(pugi::xml_node lat_node); - int32_t& operator[](std::array i_xyz); + int32_t const& operator[](array const& i_xyz); - bool are_valid_indices(const int i_xyz[3]) const; + bool are_valid_indices(array const& i_xyz) const; std::pair> distance( Position r, Direction u, const array& i_xyz) const; - std::array get_indices(Position r, Direction u) const; + void get_indices(Position r, Direction u, array& result) const; - Position - get_local_position(Position r, const std::array i_xyz) const; + Position get_local_position(Position r, const array& i_xyz) const; - int32_t& offset(int map, const int i_xyz[3]); + int32_t& offset(int map, array const& i_xyz); int32_t offset(int map, int indx) const; @@ -235,11 +228,6 @@ private: array n_cells_; //!< Number of cells along each axis Position lower_left_; //!< Global lower-left corner of the lattice Position pitch_; //!< Lattice tile width along each axis - - // Convenience aliases - int &nx {n_cells_[0]}; - int &ny {n_cells_[1]}; - int &nz {n_cells_[2]}; }; //============================================================================== @@ -249,7 +237,7 @@ class HexLattice : public Lattice public: explicit HexLattice(pugi::xml_node lat_node); - int32_t& operator[](std::array i_xyz); + int32_t const& operator[](array const& i_xyz); LatticeIter begin(); @@ -260,14 +248,13 @@ public: std::pair> distance( Position r, Direction u, const array& i_xyz) const; - std::array get_indices(Position r, Direction u) const; + void get_indices(Position r, Direction u, array& result) const; - Position - get_local_position(Position r, const std::array i_xyz) const; + Position get_local_position(Position r, const array& i_xyz) const; bool is_valid_index(int indx) const; - int32_t& offset(int map, const int i_xyz[3]); + int32_t& offset(int map, array const& i_xyz); int32_t offset(int map, int indx) const; diff --git a/src/geometry.cpp b/src/geometry.cpp index 496c510d3..59607f637 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -159,9 +159,8 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) if (c_i.type_ == Fill::UNIVERSE) { offset += c_i.offset_[c.distribcell_index_]; } else if (c_i.type_ == Fill::LATTICE) { - auto& lat {*model::lattices[p.coord_[i + 1].lattice]}; - int i_xyz[3] {p.coord_[i + 1].lattice_x, p.coord_[i + 1].lattice_y, - p.coord_[i + 1].lattice_z}; + auto& lat {*model::lattices[p.coord(i + 1).lattice]}; + const auto& i_xyz {p.coord(i + 1).lattice_i}; if (lat.are_valid_indices(i_xyz)) { offset += lat.offset(c.distribcell_index_, i_xyz); } @@ -226,16 +225,14 @@ find_cell_inner(Particle& p, const NeighborList* neighbor_list) } // Determine lattice indices. - auto i_xyz = lat.get_indices(coord.r, coord.u); + auto& i_xyz {coord.lattice_i}; + lat.get_indices(coord.r, coord.u, i_xyz); // Get local position in appropriate lattice cell coord.r = lat.get_local_position(coord.r, i_xyz); // Set lattice indices. coord.lattice = c.fill_; - coord.lattice_x = i_xyz[0]; - coord.lattice_y = i_xyz[1]; - coord.lattice_z = i_xyz[2]; // Set the lower coordinate level universe. if (lat.are_valid_indices(i_xyz)) { @@ -318,10 +315,9 @@ cross_lattice(Particle& p, const BoundaryInfo& boundary) } // Set the lattice indices. - coord.lattice_x += boundary.lattice_translation[0]; - coord.lattice_y += boundary.lattice_translation[1]; - coord.lattice_z += boundary.lattice_translation[2]; - std::array i_xyz {coord.lattice_x, coord.lattice_y, coord.lattice_z}; + coord.lattice_i[0] += boundary.lattice_translation[0]; + coord.lattice_i[1] += boundary.lattice_translation[1]; + coord.lattice_i[2] += boundary.lattice_translation[2]; // Set the new coordinate position. const auto& upper_coord {p.coord(p.n_coord() - 2)}; @@ -331,9 +327,9 @@ cross_lattice(Particle& p, const BoundaryInfo& boundary) if (!cell->rotation_.empty()) { r = r.rotate(cell->rotation_); } - p.r_local() = lat.get_local_position(r, i_xyz); + p.r_local() = lat.get_local_position(r, coord.lattice_i); - if (!lat.are_valid_indices(i_xyz)) { + if (!lat.are_valid_indices(coord.lattice_i)) { // The particle is outside the lattice. Search for it from the base coords. p.n_coord() = 1; bool found = exhaustive_find_cell(p); @@ -345,7 +341,7 @@ cross_lattice(Particle& p, const BoundaryInfo& boundary) } else { // Find cell in next lattice element. - p.coord_[p.n_coord_-1].universe = lat[i_xyz]; + p.coord(p.n_coord() - 1).universe = lat[coord.lattice_i]; bool found = exhaustive_find_cell(p); if (!found) { @@ -387,13 +383,12 @@ BoundaryInfo distance_to_boundary(Particle& p) // Find the distance to the next lattice tile crossing. if (coord.lattice != C_NONE) { auto& lat {*model::lattices[coord.lattice]}; - std::array i_xyz {coord.lattice_x, coord.lattice_y, coord.lattice_z}; //TODO: refactor so both lattice use the same position argument (which //also means the lat.type attribute can be removed) std::pair> lattice_distance; switch (lat.type_) { case LatticeType::rect: - lattice_distance = lat.distance(r, u, i_xyz); + lattice_distance = lat.distance(r, u, coord.lattice_i); break; case LatticeType::hex: auto& cell_above {model::cells[p.coord(i - 1).cell]}; @@ -403,7 +398,7 @@ BoundaryInfo distance_to_boundary(Particle& p) r_hex = r_hex.rotate(cell_above->rotation_); } r_hex.z = coord.r.z; - lattice_distance = lat.distance(r_hex, u, i_xyz); + lattice_distance = lat.distance(r_hex, u, coord.lattice_i); break; } d_lat = lattice_distance.first; diff --git a/src/lattice.cpp b/src/lattice.cpp index 37d69b3e9..d7f5a941e 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -184,16 +184,19 @@ RectLattice::RectLattice(pugi::xml_node lat_node) if (univ_words.size() != n_cells_[0] * n_cells_[1] * n_cells_[2]) { fatal_error(fmt::format( "Expected {} universes for a rectangular lattice of size {}x{}x{} but {} " - "were specified.", nx*ny*nz, nx, ny, nz, univ_words.size())); + "were specified.", + n_cells_[0] * n_cells_[1] * n_cells_[2], n_cells_[0], n_cells_[1], + n_cells_[2], univ_words.size())); } // Parse the universes. - universes_.resize(nx*ny*nz, C_NONE); - for (int iz = 0; iz < nz; iz++) { - for (int iy = ny-1; iy > -1; iy--) { - for (int ix = 0; ix < nx; ix++) { - int indx1 = nx*ny*iz + nx*(ny-iy-1) + ix; - int indx2 = nx*ny*iz + nx*iy + ix; + universes_.resize(n_cells_[0] * n_cells_[1] * n_cells_[2], C_NONE); + for (int iz = 0; iz < n_cells_[2]; iz++) { + for (int iy = n_cells_[1] - 1; iy > -1; iy--) { + for (int ix = 0; ix < n_cells_[0]; ix++) { + int indx1 = n_cells_[0] * n_cells_[1] * iz + + n_cells_[0] * (n_cells_[1] - iy - 1) + ix; + int indx2 = n_cells_[0] * n_cells_[1] * iz + n_cells_[0] * iy + ix; universes_[indx1] = std::stoi(univ_words[indx2]); } } @@ -202,17 +205,16 @@ RectLattice::RectLattice(pugi::xml_node lat_node) //============================================================================== -int32_t& -RectLattice::operator[](std::array i_xyz) +int32_t const& RectLattice::operator[](array const& i_xyz) { - int indx = nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]; + int indx = + n_cells_[0] * n_cells_[1] * i_xyz[2] + n_cells_[0] * i_xyz[1] + i_xyz[0]; return universes_[indx]; } //============================================================================== -bool -RectLattice::are_valid_indices(const int i_xyz[3]) const +bool RectLattice::are_valid_indices(array const& i_xyz) const { return ( (i_xyz[0] >= 0) && (i_xyz[0] < n_cells_[0]) && (i_xyz[1] >= 0) && (i_xyz[1] < n_cells_[1]) @@ -279,48 +281,44 @@ std::pair> RectLattice::distance( //============================================================================== -std::array -RectLattice::get_indices(Position r, Direction u) const +void RectLattice::get_indices( + Position r, Direction u, array& result) const { // Determine x index, accounting for coincidence double ix_ {(r.x - lower_left_.x) / pitch_.x}; long ix_close {std::lround(ix_)}; - int ix; if (coincident(ix_, ix_close)) { - ix = (u.x > 0) ? ix_close : ix_close - 1; + result[0] = (u.x > 0) ? ix_close : ix_close - 1; } else { - ix = std::floor(ix_); + result[0] = std::floor(ix_); } // Determine y index, accounting for coincidence double iy_ {(r.y - lower_left_.y) / pitch_.y}; long iy_close {std::lround(iy_)}; - int iy; if (coincident(iy_, iy_close)) { - iy = (u.y > 0) ? iy_close : iy_close - 1; + result[1] = (u.y > 0) ? iy_close : iy_close - 1; } else { - iy = std::floor(iy_); + result[1] = std::floor(iy_); } // Determine z index, accounting for coincidence - int iz = 0; + result[2] = 0; if (is_3d_) { double iz_ {(r.z - lower_left_.z) / pitch_.z}; long iz_close {std::lround(iz_)}; if (coincident(iz_, iz_close)) { - iz = (u.z > 0) ? iz_close : iz_close - 1; + result[2] = (u.z > 0) ? iz_close : iz_close - 1; } else { - iz = std::floor(iz_); + result[2] = std::floor(iz_); } } - return {ix, iy, iz}; } //============================================================================== -Position -RectLattice::get_local_position(Position r, const std::array i_xyz) -const +Position RectLattice::get_local_position( + Position r, const array& i_xyz) const { r.x -= (lower_left_.x + (i_xyz[0] + 0.5)*pitch_.x); r.y -= (lower_left_.y + (i_xyz[1] + 0.5)*pitch_.y); @@ -332,10 +330,11 @@ const //============================================================================== -int32_t& -RectLattice::offset(int map, const int i_xyz[3]) +int32_t& RectLattice::offset(int map, array const& i_xyz) { - return offsets_[nx*ny*nz*map + nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]]; + return offsets_[n_cells_[0] * n_cells_[1] * n_cells_[2] * map + + n_cells_[0] * n_cells_[1] * i_xyz[2] + + n_cells_[0] * i_xyz[1] + i_xyz[0]]; } //============================================================================== @@ -343,7 +342,7 @@ RectLattice::offset(int map, const int i_xyz[3]) int32_t RectLattice::offset(int map, int indx) const { - return offsets_[nx*ny*nz*map + indx]; + return offsets_[n_cells_[0] * n_cells_[1] * n_cells_[2] * map + indx]; } //============================================================================== @@ -351,9 +350,9 @@ RectLattice::offset(int map, int indx) const std::string RectLattice::index_to_string(int indx) const { - int iz {indx / (nx * ny)}; - int iy {(indx - nx * ny * iz) / nx}; - int ix {indx - nx * ny * iz - nx * iy}; + int iz {indx / (n_cells_[0] * n_cells_[1])}; + int iy {(indx - n_cells_[0] * n_cells_[1] * iz) / n_cells_[0]}; + int ix {indx - n_cells_[0] * n_cells_[1] * iz - n_cells_[0] * iy}; std::string out {std::to_string(ix)}; out += ','; out += std::to_string(iy); @@ -653,8 +652,7 @@ void HexLattice::fill_lattice_y(const vector& univ_words) //============================================================================== -int32_t& -HexLattice::operator[](std::array i_xyz) +int32_t const& HexLattice::operator[](array const& i_xyz) { int indx = (2*n_rings_-1)*(2*n_rings_-1) * i_xyz[2] + (2*n_rings_-1) * i_xyz[1] @@ -672,8 +670,7 @@ ReverseLatticeIter HexLattice::rbegin() //============================================================================== -bool -HexLattice::are_valid_indices(const int i_xyz[3]) const +bool HexLattice::are_valid_indices(array const& i_xyz) const { // Check if (x, alpha, z) indices are valid, accounting for number of rings return ((i_xyz[0] >= 0) && (i_xyz[1] >= 0) && (i_xyz[2] >= 0) @@ -826,44 +823,43 @@ std::pair> HexLattice::distance( //============================================================================== -std::array -HexLattice::get_indices(Position r, Direction u) const +void HexLattice::get_indices( + Position r, Direction u, array& result) const { // Offset the xyz by the lattice center. Position r_o {r.x - center_.x, r.y - center_.y, r.z}; if (is_3d_) {r_o.z -= center_.z;} // Index the z direction, accounting for coincidence - int iz = 0; + result[2] = 0; if (is_3d_) { double iz_ {r_o.z / pitch_[1] + 0.5 * n_axial_}; long iz_close {std::lround(iz_)}; if (coincident(iz_, iz_close)) { - iz = (u.z > 0) ? iz_close : iz_close - 1; + result[2] = (u.z > 0) ? iz_close : iz_close - 1; } else { - iz = std::floor(iz_); + result[2] = std::floor(iz_); } } - int i1, i2; if (orientation_ == Orientation::y) { // Convert coordinates into skewed bases. The (x, alpha) basis is used to // find the index of the global coordinates to within 4 cells. double alpha = r_o.y - r_o.x / std::sqrt(3.0); - i1 = std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch_[0])); - i2 = std::floor(alpha / pitch_[0]); + result[0] = std::floor(r_o.x / (0.5 * std::sqrt(3.0) * pitch_[0])); + result[1] = std::floor(alpha / pitch_[0]); } else { // Convert coordinates into skewed bases. The (alpha, y) basis is used to // find the index of the global coordinates to within 4 cells. double alpha = r_o.y - r_o.x * std::sqrt(3.0); - i1 = std::floor(-alpha / (std::sqrt(3.0) * pitch_[0])); - i2 = std::floor(r_o.y / (0.5*std::sqrt(3.0) * pitch_[0])); + result[0] = std::floor(-alpha / (std::sqrt(3.0) * pitch_[0])); + result[1] = std::floor(r_o.y / (0.5 * std::sqrt(3.0) * pitch_[0])); } // Add offset to indices (the center cell is (i1, i2) = (0, 0) but // the array is offset so that the indices never go below 0). - i1 += n_rings_-1; - i2 += n_rings_-1; + result[0] += n_rings_ - 1; + result[1] += n_rings_ - 1; // Calculate the (squared) distance between the particle and the centers of // the four possible cells. Regular hexagonal tiles form a Voronoi @@ -882,14 +878,14 @@ HexLattice::get_indices(Position r, Direction u) const // is kept (i.e. the cell with the lowest dot product as the vectors will be // completely opposed if the particle is moving directly toward the center of // the cell). - int i1_chg {}; - int i2_chg {}; + int i1_chg; + int i2_chg; double d_min {INFTY}; double dp_min {INFTY}; for (int i = 0; i < 2; i++) { for (int j = 0; j < 2; j++) { // get local coordinates - const std::array i_xyz {i1 + j, i2 + i, 0}; + const array i_xyz {result[0] + j, result[1] + i, 0}; Position r_t = get_local_position(r, i_xyz); // calculate distance double d = r_t.x*r_t.x + r_t.y*r_t.y; @@ -916,17 +912,14 @@ HexLattice::get_indices(Position r, Direction u) const } // update outgoing indices - i1 += i1_chg; - i2 += i2_chg; - - return {i1, i2, iz}; + result[0] += i1_chg; + result[1] += i2_chg; } //============================================================================== -Position -HexLattice::get_local_position(Position r, const std::array i_xyz) -const +Position HexLattice::get_local_position( + Position r, const array& i_xyz) const { if (orientation_ == Orientation::y) { // x_l = x_g - (center + pitch_x*cos(30)*index_x) @@ -961,14 +954,13 @@ HexLattice::is_valid_index(int indx) const int iz = indx / (nx * ny); int iy = (indx - nx*ny*iz) / nx; int ix = indx - nx*ny*iz - nx*iy; - int i_xyz[3] {ix, iy, iz}; + array i_xyz {ix, iy, iz}; return are_valid_indices(i_xyz); } //============================================================================== -int32_t& -HexLattice::offset(int map, const int i_xyz[3]) +int32_t& HexLattice::offset(int map, array const& i_xyz) { int nx {2*n_rings_ - 1}; int ny {2*n_rings_ - 1}; diff --git a/src/output.cpp b/src/output.cpp index 8f81d9881..beac21ddd 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -181,8 +181,8 @@ void print_particle(Particle& p) if (p.coord(i).lattice != C_NONE) { const Lattice& lat {*model::lattices[p.coord(i).lattice]}; fmt::print(" Lattice = {}\n", lat.id_); - fmt::print(" Lattice position = ({},{},{})\n", p.coord_[i].lattice_x, - p.coord_[i].lattice_y, p.coord_[i].lattice_z); + fmt::print(" Lattice position = ({},{},{})\n", p.coord(i).lattice_i[0], + p.coord(i).lattice_i[1], p.coord(i).lattice_i[2]); } fmt::print(" r = {}\n", p.coord(i).r); diff --git a/src/tallies/filter_distribcell.cpp b/src/tallies/filter_distribcell.cpp index ea6fef944..c33c88cbc 100644 --- a/src/tallies/filter_distribcell.cpp +++ b/src/tallies/filter_distribcell.cpp @@ -48,10 +48,8 @@ DistribcellFilter::get_all_bins(const Particle& p, TallyEstimator estimator, if (c.type_ == Fill::UNIVERSE) { offset += c.offset_[distribcell_index]; } else if (c.type_ == Fill::LATTICE) { - auto& lat {*model::lattices[p.coord_[i+1].lattice]}; - int i_xyz[3] {p.coord_[i+1].lattice_x, - p.coord_[i+1].lattice_y, - p.coord_[i+1].lattice_z}; + auto& lat {*model::lattices[p.coord(i + 1).lattice]}; + const auto& i_xyz {p.coord(i + 1).lattice_i}; if (lat.are_valid_indices(i_xyz)) { offset += lat.offset(distribcell_index, i_xyz); }