From a4fb7b15aa0c4eb3db318661f1561e77c78024e6 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 3 Oct 2026 15:07:54 +0000 Subject: [PATCH 1/2] Count universe instances in one pass over the geometry The number of instances of each universe was found by searching the geometry from the root once per universe, so initialization took time proportional to the number of universes times the size of the geometry. Models with many universes, such as TRISO particles in a fine lattice, spent many minutes in this step. The universes are now ordered so that each comes after the universes containing it, and the counts are propagated from the root in that order, visiting each universe once. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014RN7JroEkrbkLKVHbMDap9 --- include/openmc/geometry_aux.h | 7 ++- src/geometry_aux.cpp | 43 ++++++++++++-- tests/unit_tests/test_universe_instances.py | 63 +++++++++++++++++++++ 3 files changed, 106 insertions(+), 7 deletions(-) create mode 100644 tests/unit_tests/test_universe_instances.py diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index 4dafdea5c2c..cd59627bbb1 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -84,10 +84,11 @@ void prepare_distribcell( const std::vector* user_distribcells = nullptr); //============================================================================== -//! Recursively search through the geometry and count universe instances. +//! Count the number of instances of every universe in the geometry. //! -//! This function will update Universe.n_instances_ for each -//! universe in the geometry. +//! This function will update Universe.n_instances_ for each universe in the +//! geometry. The universes are visited once each, in an order where every +//! universe comes after the universes containing it. //============================================================================== void count_universe_instances(); diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index a740740c1e6..71cf79a8bae 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -455,10 +455,45 @@ void prepare_distribcell(const std::vector* user_distribcells) void count_universe_instances() { - for (auto& univ : model::universes) { - std::unordered_map univ_count_memo; - univ->n_instances_ = count_universe_instances( - model::root_universe, univ->id_, univ_count_memo); + // Call a function with each universe filling a cell or lattice element of a + // universe, once per cell or lattice element + auto for_each_fill = [](int32_t i_univ, auto&& f) { + for (int32_t i_cell : model::universes[i_univ]->cells_) { + Cell& c = *model::cells[i_cell]; + if (c.type_ == Fill::UNIVERSE) { + f(c.fill_); + } else if (c.type_ == Fill::LATTICE) { + Lattice& lat = *model::lattices[c.fill_]; + for (auto it = lat.begin(); it != lat.end(); ++it) + f(*it); + } + } + }; + + // Order the universes reachable from the root so that every universe comes + // after all universes that contain it + vector order; + vector visited(model::universes.size(), false); + auto visit = [&](int32_t i_univ, auto& self) -> void { + visited[i_univ] = true; + for_each_fill(i_univ, [&](int32_t next) { + if (!visited[next]) + self(next, self); + }); + order.push_back(i_univ); + }; + visit(model::root_universe, visit); + + // The number of instances of a universe is the sum over the cells and + // lattice elements it fills of the number of instances of the universe + // containing them. Universes not reachable from the root have none. + for (auto& univ : model::universes) + univ->n_instances_ = 0; + model::universes[model::root_universe]->n_instances_ = 1; + for (auto it = order.rbegin(); it != order.rend(); ++it) { + int n = model::universes[*it]->n_instances_; + for_each_fill( + *it, [&](int32_t next) { model::universes[next]->n_instances_ += n; }); } } diff --git a/tests/unit_tests/test_universe_instances.py b/tests/unit_tests/test_universe_instances.py new file mode 100644 index 00000000000..bba7878eca1 --- /dev/null +++ b/tests/unit_tests/test_universe_instances.py @@ -0,0 +1,63 @@ +import openmc +import openmc.lib + + +def test_universe_instances(run_in_tmpdir): + """Number of instances of cells in universes that are used through + nested lattices and directly, compared with the Python API.""" + water = openmc.Material() + water.add_nuclide('H1', 2.0) + water.add_nuclide('O16', 1.0) + water.set_density('g/cm3', 1.0) + + # Pin universes + cyl = openmc.ZCylinder(r=0.4) + fuel_pin = openmc.Universe(cells=[ + openmc.Cell(fill=water, region=-cyl), + openmc.Cell(fill=water, region=+cyl)]) + guide_pin = openmc.Universe(cells=[openmc.Cell(fill=water)]) + + # Two assembly types, used several times in a core lattice + def assembly(pins): + lat = openmc.RectLattice() + lat.lower_left = (-1.5, -1.5) + lat.pitch = (1.0, 1.0) + lat.universes = pins + return openmc.Universe(cells=[openmc.Cell(fill=lat)]) + + f, g = fuel_pin, guide_pin + assembly_a = assembly([[f, f, f], [f, g, f], [f, f, f]]) + assembly_b = assembly([[f, g, f], [g, g, g], [f, g, f]]) + core_lat = openmc.RectLattice() + core_lat.lower_left = (-6.0, -3.0) + core_lat.pitch = (3.0, 3.0) + core_lat.universes = [[assembly_a, assembly_b, assembly_a, assembly_a], + [assembly_b, assembly_a, assembly_b, assembly_a]] + + # The fuel pin is also used directly next to the core + box = openmc.model.RectangularParallelepiped( + -6.0, 7.0, -3.0, 3.0, -1.0, 1.0, boundary_type='vacuum') + x_core = openmc.XPlane(6.0) + root = openmc.Universe(cells=[ + openmc.Cell(fill=core_lat, region=-box & -x_core), + openmc.Cell(fill=fuel_pin, region=-box & +x_core)]) + + model = openmc.Model() + model.geometry = openmc.Geometry(root) + model.materials = openmc.Materials([water]) + model.settings.particles = 10 + model.settings.batches = 1 + model.settings.run_mode = 'fixed source' + model.export_to_model_xml() + + model.geometry.determine_paths() + with openmc.lib.run_in_memory(): + for cell in model.geometry.get_all_cells().values(): + assert openmc.lib.cells[cell.id].num_instances == \ + cell.num_instances + + # Five A and three B assemblies, plus the pin used directly + for cell in fuel_pin.cells.values(): + assert openmc.lib.cells[cell.id].num_instances == 5*8 + 3*4 + 1 + for cell in guide_pin.cells.values(): + assert openmc.lib.cells[cell.id].num_instances == 5*1 + 3*5 From ac50777ce83ee4d9a91caa9a45b8ae264ac7e7a3 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 3 Oct 2026 18:22:07 +0000 Subject: [PATCH 2/2] Store distribcell offsets only where they can be read Every fill cell and lattice tile stored one offset per distribcell map, where a map is a universe containing distributed cells. By default every universe containing a material cell is a map, so the tables grew with the number of such universes times the size of the geometry. In models with many universes, such as TRISO particles in a fine lattice where each tile has its own universe, they could take tens of gigabytes. An offset is only read when the map's universe is in the cell's fill or the tile's universe. Each universe and lattice now keeps the sorted list of maps it contains, and each fill cell and tile stores offsets for those maps only. Maps are numbered in post-order, so the maps in a universe are mostly consecutive and a lookup is usually a single range check. The offsets are found in one visit of each universe and lattice from the root, instead of one pass over the geometry per map. Instance numbers are unchanged. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014RN7JroEkrbkLKVHbMDap9 --- include/openmc/cell.h | 3 +- include/openmc/distribcell_offsets.h | 54 +++++ include/openmc/geometry_aux.h | 14 -- include/openmc/lattice.h | 35 ++- include/openmc/universe.h | 7 + src/geometry_aux.cpp | 323 +++++++++++++++++++-------- src/lattice.cpp | 62 ----- 7 files changed, 301 insertions(+), 197 deletions(-) create mode 100644 include/openmc/distribcell_offsets.h diff --git a/include/openmc/cell.h b/include/openmc/cell.h index c2a09d774b5..4f90853a655 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -13,6 +13,7 @@ #include "openmc/bounding_box.h" #include "openmc/constants.h" +#include "openmc/distribcell_offsets.h" #include "openmc/memory.h" // for unique_ptr #include "openmc/neighbor_list.h" #include "openmc/position.h" @@ -410,7 +411,7 @@ class Cell { //! also present at the end of the vector, making it of length 12. vector rotation_; - vector offset_; //!< Distribcell offset table + DistribcellOffsets offset_; //!< Distribcell offsets // Right now, either CSG or DAGMC cells are used. virtual GeometryType geom_type() const = 0; diff --git a/include/openmc/distribcell_offsets.h b/include/openmc/distribcell_offsets.h new file mode 100644 index 00000000000..d95dcbc90bf --- /dev/null +++ b/include/openmc/distribcell_offsets.h @@ -0,0 +1,54 @@ +#ifndef OPENMC_DISTRIBCELL_OFFSETS_H +#define OPENMC_DISTRIBCELL_OFFSETS_H + +#include + +#include "openmc/vector.h" + +namespace openmc { + +//============================================================================== +//! Distributed cell offsets of a fill cell or lattice tile. +// +//! A map numbers the instances of a universe containing distributed cells. The +//! offset for a map is the number of instances of the map's universe before +//! the fill cell in its universe, or before the tile in its lattice. Offsets +//! are only read for maps whose universe is in the cell's fill or the tile, so +//! only those are stored, in the order of the sorted maps of the fill or tile. +//============================================================================== + +class DistribcellOffsets { +public: + DistribcellOffsets() = default; + + //! \param maps Sorted maps in the fill or tile + //! \param values Offset for each map + //! Both must outlive this object. + DistribcellOffsets(const vector& maps, const int32_t* values); + + //! \return Offset for a map, or 0 if the map's universe isn't contained + int32_t operator[](int32_t map) const + { + // Maps are numbered so that the maps in a universe are mostly consecutive, + // so first look in the consecutive maps at the start of the list + auto d = static_cast(map - first_map_); + if (d < static_cast(n_first_)) + return values_[d]; + return find(map); + } + + //! \return Whether the map's universe is contained + bool contains(int32_t map) const; + +private: + int32_t find(int32_t map) const; + + const int32_t* values_ {nullptr}; + const vector* maps_ {nullptr}; + int32_t first_map_ {0}; //!< First map + int32_t n_first_ {0}; //!< Number of consecutive maps from the first +}; + +} // namespace openmc + +#endif // OPENMC_DISTRIBCELL_OFFSETS_H diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index cd59627bbb1..1e5ce7f2919 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -93,20 +93,6 @@ void prepare_distribcell( void count_universe_instances(); -//============================================================================== -//! Recursively search through universes and count universe instances. -//! \param search_univ The index of the universe to begin searching from. -//! \param target_univ_id The ID of the universe to be counted. -//! \param univ_count_memo Memoized counts that make this function faster for -//! large systems. The first call to this function for each target_univ_id -//! should start with an empty memo. -//! \return The number of instances of target_univ_id in the geometry tree under -//! search_univ. -//============================================================================== - -int count_universe_instances(int32_t search_univ, int32_t target_univ_id, - std::unordered_map& univ_count_memo); - //============================================================================== //! Build a character array representing the path to a distribcell instance. //! \param target_cell The index of the Cell in the global Cell array. diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index ca40bbc2a38..9ddc36095ac 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -10,6 +10,7 @@ #include "openmc/array.h" #include "openmc/constants.h" +#include "openmc/distribcell_offsets.h" #include "openmc/memory.h" #include "openmc/position.h" #include "openmc/vector.h" @@ -50,7 +51,13 @@ class Lattice { LatticeType type_; vector universes_; //!< Universes filling each lattice tile int32_t outer_ {NO_OUTER_UNIVERSE}; //!< Universe tiled outside the lattice - vector offsets_; //!< Distribcell offset table + + //! Sorted distributed cell maps of the universes in the tiles and the outer + //! universe + vector distribcell_maps_; + + vector offset_values_; //!< Distribcell offsets of all tiles + vector offsets_; //!< Distribcell offsets of each tile explicit Lattice(pugi::xml_node lat_node); @@ -68,17 +75,6 @@ class Lattice { //! Convert internal universe values from IDs to indices using universe_map. void adjust_indices(); - //! Allocate offset table for distribcell. - void allocate_offset_table(int n_maps) - { - offsets_.resize(n_maps * universes_.size()); - std::fill(offsets_.begin(), offsets_.end(), C_NONE); - } - - //! Populate the distribcell offset tables. - int32_t fill_offset_table(int32_t target_univ_id, int map, - std::unordered_map& univ_count_memo); - //! \brief Check lattice indices. //! \param i_xyz[3] The indices for a lattice tile. //! \return true if the given indices fit within the lattice bounds. False @@ -135,14 +131,17 @@ class Lattice { //! \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 array& i_xyz) = 0; + int32_t offset(int map, const array& i_xyz) const + { + return offset(map, get_flat_index(i_xyz)); + } //! \brief Get the distribcell offset for a lattice tile. //! \param The map index for the target cell. //! \param indx The index for a lattice tile. //! \return Distribcell offset i.e. the largest instance number for the target //! cell found in the geometry tree for this lattice index. - virtual int32_t offset(int map, int indx) const = 0; + int32_t offset(int map, int indx) const { return offsets_[indx][map]; } //! \brief Convert an array index to a useful human-readable string. //! \param indx The index for a lattice tile. @@ -234,10 +233,6 @@ class RectLattice : public Lattice { Direction get_normal( const array& i_xyz, bool& is_valid) const override; - int32_t& offset(int map, const array& i_xyz) override; - - int32_t offset(int map, int indx) const override; - std::string index_to_string(int indx) const override; void to_hdf5_inner(hid_t group_id) const override; @@ -284,10 +279,6 @@ class HexLattice : public Lattice { bool is_valid_index(int indx) const override; - int32_t& offset(int map, const array& i_xyz) override; - - int32_t offset(int map, int indx) const override; - std::string index_to_string(int indx) const override; void to_hdf5_inner(hid_t group_id) const override; diff --git a/include/openmc/universe.h b/include/openmc/universe.h index b7450224f4f..f073fcf9fe9 100644 --- a/include/openmc/universe.h +++ b/include/openmc/universe.h @@ -31,6 +31,13 @@ class Universe { vector cells_; //!< Cells within this universe int32_t n_instances_; //!< Number of instances of this universe + //! Sorted distributed cell maps of the universes in this universe, + //! including itself + vector distribcell_maps_; + + //! Distributed cell offsets of the fill cells in this universe + vector offset_values_; + //! \brief Write universe information to an HDF5 group. //! \param group_id An HDF5 group id. virtual void to_hdf5(hid_t group_id) const; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 71cf79a8bae..609d3ed95cd 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -329,6 +329,214 @@ int32_t find_root_universe() //============================================================================== +DistribcellOffsets::DistribcellOffsets( + const vector& maps, const int32_t* values) + : values_(values), maps_(&maps) +{ + if (!maps.empty()) { + first_map_ = maps[0]; + n_first_ = 1; + while (n_first_ < maps.size() && maps[n_first_] == first_map_ + n_first_) { + ++n_first_; + } + } +} + +int32_t DistribcellOffsets::find(int32_t map) const +{ + if (!maps_) + return 0; + auto it = std::lower_bound(maps_->begin(), maps_->end(), map); + if (it == maps_->end() || *it != map) + return 0; + return values_[it - maps_->begin()]; +} + +bool DistribcellOffsets::contains(int32_t map) const +{ + return maps_ && std::binary_search(maps_->begin(), maps_->end(), map); +} + +//============================================================================== + +namespace { + +//! Builds the distributed cell offsets, visiting each universe and lattice +//! once, after the universes they contain. +class DistribcellBuilder { +public: + explicit DistribcellBuilder(const vector& is_target) + : is_target_(is_target), univ_map_(model::universes.size(), C_NONE), + univ_visited_(model::universes.size(), false), + lat_visited_(model::lattices.size(), false), + univ_counts_(model::universes.size()), lat_counts_(model::lattices.size()) + {} + + //! Visit the geometry from the root universe. Maps are numbered in + //! post-order, so that the maps in a universe are mostly consecutive. + void build() { visit_universe(model::root_universe); } + + int32_t map(int32_t univ) const { return univ_map_[univ]; } + +private: + //! Find the maps in a universe, store the offsets of its fill cells, and + //! count the instances of each map. + void visit_universe(int32_t univ_indx) + { + if (univ_visited_[univ_indx]) + return; + univ_visited_[univ_indx] = true; + Universe& univ = *model::universes[univ_indx]; + + auto& maps = univ.distribcell_maps_; + maps.clear(); + for (int32_t cell_indx : univ.cells_) { + const Cell& c = *model::cells[cell_indx]; + if (c.type_ == Fill::UNIVERSE) { + visit_universe(c.fill_); + } else if (c.type_ == Fill::LATTICE) { + visit_lattice(c.fill_); + } else { + continue; + } + const auto& fill_maps = this->fill_maps(c); + maps.insert(maps.end(), fill_maps.begin(), fill_maps.end()); + } + sort_unique(maps); + + // The universe's own map comes after the maps it contains + if (is_target_[univ_indx]) { + univ_map_[univ_indx] = n_maps_++; + maps.push_back(univ_map_[univ_indx]); + } + + // The instances counted before a cell are its offsets. The offsets of all + // cells are stored together, so the storage is allocated once. + size_t n_values = 0; + for (int32_t cell_indx : univ.cells_) { + const Cell& c = *model::cells[cell_indx]; + if (c.type_ != Fill::MATERIAL) + n_values += fill_maps(c).size(); + } + univ.offset_values_.clear(); + univ.offset_values_.reserve(n_values); + auto& counts = univ_counts_[univ_indx]; + counts.assign(maps.size(), 0); + for (int32_t cell_indx : univ.cells_) { + Cell& c = *model::cells[cell_indx]; + if (c.type_ == Fill::MATERIAL) + continue; + const auto& fill_maps = this->fill_maps(c); + c.offset_ = store(univ.offset_values_, maps, counts, fill_maps); + add_counts(maps, counts, fill_maps, + c.type_ == Fill::UNIVERSE ? univ_counts_[c.fill_] + : lat_counts_[c.fill_]); + } + if (is_target_[univ_indx]) + counts.back() = 1; + } + + //! Find the maps in a lattice, store the offsets of its tiles, and count the + //! instances of each map. The outer universe is not counted, but its maps are + //! included since a lattice cell's offsets are read for particles in the + //! outer universe. + void visit_lattice(int32_t lat_indx) + { + if (lat_visited_[lat_indx]) + return; + lat_visited_[lat_indx] = true; + Lattice& lat = *model::lattices[lat_indx]; + + auto& maps = lat.distribcell_maps_; + maps.clear(); + std::unordered_set univs; + for (LatticeIter it = lat.begin(); it != lat.end(); ++it) { + if (univs.insert(*it).second) { + visit_universe(*it); + const auto& univ_maps = model::universes[*it]->distribcell_maps_; + maps.insert(maps.end(), univ_maps.begin(), univ_maps.end()); + } + } + if (lat.outer_ != NO_OUTER_UNIVERSE) { + visit_universe(lat.outer_); + const auto& outer_maps = model::universes[lat.outer_]->distribcell_maps_; + maps.insert(maps.end(), outer_maps.begin(), outer_maps.end()); + } + sort_unique(maps); + + // The instances counted before a tile are its offsets. The offsets of all + // tiles are stored together, so the storage is allocated once. + size_t n_values = 0; + for (LatticeIter it = lat.begin(); it != lat.end(); ++it) { + n_values += model::universes[*it]->distribcell_maps_.size(); + } + lat.offset_values_.clear(); + lat.offset_values_.reserve(n_values); + auto& counts = lat_counts_[lat_indx]; + counts.assign(maps.size(), 0); + lat.offsets_.assign(lat.universes_.size(), {}); + for (LatticeIter it = lat.begin(); it != lat.end(); ++it) { + const auto& univ_maps = model::universes[*it]->distribcell_maps_; + lat.offsets_[it.indx_] = + store(lat.offset_values_, maps, counts, univ_maps); + add_counts(maps, counts, univ_maps, univ_counts_[*it]); + } + } + + //! Sorted maps in the fill of a fill cell + static const vector& fill_maps(const Cell& c) + { + return c.type_ == Fill::UNIVERSE + ? model::universes[c.fill_]->distribcell_maps_ + : model::lattices[c.fill_]->distribcell_maps_; + } + + //! Append the counts of a subset of maps to storage with enough capacity + //! \return Offsets of the subset + static DistribcellOffsets store(vector& storage, + const vector& maps, const vector& counts, + const vector& subset) + { + const int32_t* values = storage.data() + storage.size(); + for (int32_t m : subset) { + storage.push_back(counts[position(maps, m)]); + } + return DistribcellOffsets(subset, values); + } + + //! Add the counts of a subset of maps to counts + static void add_counts(const vector& maps, vector& counts, + const vector& subset, const vector& subset_counts) + { + for (int32_t i = 0; i < subset.size(); ++i) { + counts[position(maps, subset[i])] += subset_counts[i]; + } + } + + //! Position of a map in a sorted list of maps that contains it + static int32_t position(const vector& maps, int32_t map) + { + return std::lower_bound(maps.begin(), maps.end(), map) - maps.begin(); + } + + static void sort_unique(vector& v) + { + std::sort(v.begin(), v.end()); + v.erase(std::unique(v.begin(), v.end()), v.end()); + } + + const vector& is_target_; + int32_t n_maps_ {0}; + vector univ_map_; + vector univ_visited_; + vector lat_visited_; + // Number of instances of each map in a universe or lattice + vector> univ_counts_; + vector> lat_counts_; +}; + +} // namespace + void prepare_distribcell(const std::vector* user_distribcells) { write_message("Preparing distributed cell instances...", 5); @@ -398,56 +606,19 @@ void prepare_distribcell(const std::vector* user_distribcells) } } - // Search through universes for material cells and assign each one a - // distribcell array index according to the containing universe. - vector target_univ_ids; - for (const auto& u : model::universes) { - for (auto idx : u->cells_) { - if (distribcells.find(idx) != distribcells.end()) { - if (!contains(target_univ_ids, u->id_)) { - target_univ_ids.push_back(u->id_); - } - model::cells[idx]->distribcell_index_ = - std::find(target_univ_ids.begin(), target_univ_ids.end(), u->id_) - - target_univ_ids.begin(); - } - } + // Each universe containing distributed cells gets a map, which numbers the + // instances of that universe. + vector is_target(model::universes.size(), false); + for (auto idx : distribcells) { + is_target[model::cells[idx]->universe_] = true; } - // Allocate the cell and lattice offset tables. - int n_maps = target_univ_ids.size(); - for (auto& c : model::cells) { - if (c->type_ != Fill::MATERIAL) { - c->offset_.resize(n_maps, C_NONE); - } - } - for (auto& lat : model::lattices) { - lat->allocate_offset_table(n_maps); - } - -// Fill the cell and lattice offset tables. -#pragma omp parallel for - for (int map = 0; map < target_univ_ids.size(); map++) { - auto target_univ_id = target_univ_ids[map]; - std::unordered_map univ_count_memo; - for (const auto& univ : model::universes) { - int32_t offset = 0; - for (int32_t cell_indx : univ->cells_) { - Cell& c = *model::cells[cell_indx]; - - if (c.type_ == Fill::UNIVERSE) { - c.offset_[map] = offset; - int32_t search_univ = c.fill_; - offset += count_universe_instances( - search_univ, target_univ_id, univ_count_memo); - - } else if (c.type_ == Fill::LATTICE) { - c.offset_[map] = offset; - Lattice& lat = *model::lattices[c.fill_]; - offset += lat.fill_offset_table(target_univ_id, map, univ_count_memo); - } - } - } + DistribcellBuilder builder(is_target); + builder.build(); + + for (auto idx : distribcells) { + Cell& c = *model::cells[idx]; + c.distribcell_index_ = builder.map(c.universe_); } } @@ -499,47 +670,6 @@ void count_universe_instances() //============================================================================== -int count_universe_instances(int32_t search_univ, int32_t target_univ_id, - std::unordered_map& univ_count_memo) -{ - // If this is the target, it can't contain itself. - if (model::universes[search_univ]->id_ == target_univ_id) { - return 1; - } - - // If we have already counted the number of instances, reuse that value. - auto search = univ_count_memo.find(search_univ); - if (search != univ_count_memo.end()) { - return search->second; - } - - int count {0}; - for (int32_t cell_indx : model::universes[search_univ]->cells_) { - Cell& c = *model::cells[cell_indx]; - - if (c.type_ == Fill::UNIVERSE) { - int32_t next_univ = c.fill_; - count += - count_universe_instances(next_univ, target_univ_id, univ_count_memo); - - } else if (c.type_ == Fill::LATTICE) { - Lattice& lat = *model::lattices[c.fill_]; - for (auto it = lat.begin(); it != lat.end(); ++it) { - int32_t next_univ = *it; - count += - count_universe_instances(next_univ, target_univ_id, univ_count_memo); - } - } - } - - // Remember the number of instances in this universe. - univ_count_memo[search_univ] = count; - - return count; -} - -//============================================================================== - std::string distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, const Universe& search_univ, int32_t offset) { @@ -565,14 +695,10 @@ std::string distribcell_path_inner(int32_t target_cell, int32_t map, for (; cell_it != search_univ.cells_.crend(); ++cell_it) { Cell& c = *model::cells[*cell_it]; - // Material cells don't contain other cells so ignore them. - if (c.type_ != Fill::MATERIAL) { + // Material cells don't contain other cells and other cells may not contain + // the target's universe, so ignore them. + if (c.type_ != Fill::MATERIAL && c.offset_.contains(map)) { int32_t temp_offset = offset + c.offset_[map]; - if (c.type_ == Fill::LATTICE) { - Lattice& lat = *model::lattices[c.fill_]; - int32_t indx = lat.universes_.size() * map + lat.begin().indx_; - temp_offset += lat.offsets_[indx]; - } // The desired cell is the first cell that gives an offset smaller or // equal to the target offset. @@ -605,8 +731,9 @@ std::string distribcell_path_inner(int32_t target_cell, int32_t map, Lattice& lat = *model::lattices[c.fill_]; path << "l" << lat.id_; for (ReverseLatticeIter it = lat.rbegin(); it != lat.rend(); ++it) { - int32_t indx = lat.universes_.size() * map + it.indx_; - int32_t temp_offset = offset + lat.offsets_[indx] + c.offset_[map]; + if (!lat.offsets_[it.indx_].contains(map)) + continue; + int32_t temp_offset = offset + lat.offset(map, it.indx_) + c.offset_[map]; if (temp_offset <= target_offset) { offset = temp_offset; path << "(" << lat.index_to_string(it.indx_) << ")->"; diff --git a/src/lattice.cpp b/src/lattice.cpp index 9e56c58e8ec..eebce4d7c2c 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -104,33 +104,6 @@ void Lattice::adjust_indices() //============================================================================== -int32_t Lattice::fill_offset_table(int32_t target_univ_id, int map, - std::unordered_map& univ_count_memo) -{ - // If the offsets have already been determined for this "map", don't bother - // recalculating all of them and just return the total offset. Note that the - // offsets_ array doesn't actually include the offset accounting for the last - // universe, so we get the before-last offset for the given map and then - // explicitly add the count for the last universe. - if (offsets_[map * universes_.size() + this->begin().indx_] != C_NONE) { - int last_offset = - offsets_[(map + 1) * universes_.size() - this->begin().indx_ - 1]; - int last_univ = this->back(); - return last_offset + - count_universe_instances(last_univ, target_univ_id, univ_count_memo); - } - - int32_t offset = 0; - for (LatticeIter it = begin(); it != end(); ++it) { - offsets_[map * universes_.size() + it.indx_] = offset; - offset += count_universe_instances(*it, target_univ_id, univ_count_memo); - } - - return offset; -} - -//============================================================================== - void Lattice::to_hdf5(hid_t lattices_group) const { // Make a group for the lattice. @@ -369,22 +342,6 @@ Direction RectLattice::get_normal( //============================================================================== -int32_t& RectLattice::offset(int map, const array& i_xyz) -{ - 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]]; -} - -//============================================================================== - -int32_t RectLattice::offset(int map, int indx) const -{ - return offsets_[n_cells_[0] * n_cells_[1] * n_cells_[2] * map + indx]; -} - -//============================================================================== - std::string RectLattice::index_to_string(int indx) const { int iz {indx / (n_cells_[0] * n_cells_[1])}; @@ -1113,25 +1070,6 @@ bool HexLattice::is_valid_index(int indx) const //============================================================================== -int32_t& HexLattice::offset(int map, const array& i_xyz) -{ - int nx {2 * n_rings_ - 1}; - int ny {2 * n_rings_ - 1}; - int nz {n_axial_}; - return offsets_[nx * ny * nz * map + nx * ny * i_xyz[2] + nx * i_xyz[1] + - i_xyz[0]]; -} - -int32_t HexLattice::offset(int map, int indx) const -{ - int nx {2 * n_rings_ - 1}; - int ny {2 * n_rings_ - 1}; - int nz {n_axial_}; - return offsets_[nx * ny * nz * map + indx]; -} - -//============================================================================== - std::string HexLattice::index_to_string(int indx) const { int nx {2 * n_rings_ - 1};