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 4dafdea5c2c..1e5ce7f2919 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -84,28 +84,15 @@ 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(); -//============================================================================== -//! 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 a740740c1e6..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,109 +606,66 @@ 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(); - } - } - } - - // 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); - } - } - } + // 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; } -} -//============================================================================== + DistribcellBuilder builder(is_target); + builder.build(); -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); + for (auto idx : distribcells) { + Cell& c = *model::cells[idx]; + c.distribcell_index_ = builder.map(c.universe_); } } //============================================================================== -int count_universe_instances(int32_t search_univ, int32_t target_univ_id, - std::unordered_map& univ_count_memo) +void count_universe_instances() { - // 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); + // 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; }); } - - // Remember the number of instances in this universe. - univ_count_memo[search_univ] = count; - - return count; } //============================================================================== @@ -530,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. @@ -570,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}; 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