From ce4176deb5b8f64a45c46ad87063c76140ce77dd Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 17 Jun 2024 16:14:12 -0500 Subject: [PATCH] Hexagonal lattice iterators (#2921) Co-authored-by: Paul Romano --- include/openmc/lattice.h | 35 +++--- src/geometry_aux.cpp | 3 +- src/lattice.cpp | 46 +++++-- .../cell_instances/test_hex_multilattice.py | 119 ++++++++++++++++++ .../test_rect_multilattice.py} | 22 ++-- 5 files changed, 192 insertions(+), 33 deletions(-) create mode 100644 tests/unit_tests/cell_instances/test_hex_multilattice.py rename tests/unit_tests/{test_cell_instance.py => cell_instances/test_rect_multilattice.py} (84%) diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index 11dda4a6b8..429d81ddda 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -56,13 +56,14 @@ public: virtual ~Lattice() {} - virtual int32_t const& operator[](array const& i_xyz) = 0; + virtual const int32_t& operator[](const array& i_xyz) = 0; virtual LatticeIter begin(); - LatticeIter end(); + virtual LatticeIter end(); + virtual int32_t& back(); virtual ReverseLatticeIter rbegin(); - ReverseLatticeIter rend(); + virtual ReverseLatticeIter rend(); //! Convert internal universe values from IDs to indices using universe_map. void adjust_indices(); @@ -81,7 +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(array const& i_xyz) const = 0; + virtual bool are_valid_indices(const array& i_xyz) const = 0; //! \brief Find the next lattice surface crossing //! \param r A 3D Cartesian coordinate. @@ -125,7 +126,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, array const& i_xyz) = 0; + virtual int32_t& offset(int map, const array& i_xyz) = 0; //! \brief Get the distribcell offset for a lattice tile. //! \param The map index for the target cell. @@ -167,12 +168,12 @@ public: LatticeIter& operator++() { - while (indx_ < lat_.universes_.size()) { + while (indx_ < lat_.end().indx_) { ++indx_; if (lat_.is_valid_index(indx_)) return *this; } - indx_ = lat_.universes_.size(); + indx_ = lat_.end().indx_; return *this; } @@ -190,7 +191,7 @@ public: ReverseLatticeIter& operator++() { - while (indx_ > -1) { + while (indx_ > lat_.begin().indx_ - 1) { --indx_; if (lat_.is_valid_index(indx_)) return *this; @@ -206,9 +207,9 @@ class RectLattice : public Lattice { public: explicit RectLattice(pugi::xml_node lat_node); - int32_t const& operator[](array const& i_xyz) override; + const int32_t& operator[](const array& i_xyz) override; - bool are_valid_indices(array const& i_xyz) const override; + bool are_valid_indices(const array& i_xyz) const override; std::pair> distance( Position r, Direction u, const array& i_xyz) const override; @@ -221,7 +222,7 @@ public: Position get_local_position( Position r, const array& i_xyz) const override; - int32_t& offset(int map, array const& i_xyz) override; + int32_t& offset(int map, const array& i_xyz) override; int32_t offset(int map, int indx) const override; @@ -241,13 +242,19 @@ class HexLattice : public Lattice { public: explicit HexLattice(pugi::xml_node lat_node); - int32_t const& operator[](array const& i_xyz) override; + const int32_t& operator[](const array& i_xyz) override; LatticeIter begin() override; ReverseLatticeIter rbegin() override; - bool are_valid_indices(array const& i_xyz) const override; + LatticeIter end() override; + + int32_t& back() override; + + ReverseLatticeIter rend() override; + + bool are_valid_indices(const array& i_xyz) const override; std::pair> distance( Position r, Direction u, const array& i_xyz) const override; @@ -262,7 +269,7 @@ public: bool is_valid_index(int indx) const override; - int32_t& offset(int map, array const& i_xyz) override; + int32_t& offset(int map, const array& i_xyz) override; int32_t offset(int map, int indx) const override; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 10b239be83..050d4db968 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -530,7 +530,8 @@ std::string distribcell_path_inner(int32_t target_cell, int32_t map, if (c.type_ != Fill::MATERIAL) { int32_t temp_offset; if (c.type_ == Fill::UNIVERSE) { - temp_offset = offset + c.offset_[map]; + temp_offset = + offset + c.offset_[map]; // TODO: should also apply to lattice fills? } else { Lattice& lat = *model::lattices[c.fill_]; int32_t indx = lat.universes_.size() * map + lat.begin().indx_; diff --git a/src/lattice.cpp b/src/lattice.cpp index fa2e2828ec..73005ad096 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -58,6 +58,11 @@ LatticeIter Lattice::end() return LatticeIter(*this, universes_.size()); } +int32_t& Lattice::back() +{ + return universes_.back(); +} + ReverseLatticeIter Lattice::rbegin() { return ReverseLatticeIter(*this, universes_.size() - 1); @@ -106,9 +111,10 @@ int32_t Lattice::fill_offset_table(int32_t offset, int32_t target_univ_id, // 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()] != C_NONE) { - int last_offset = offsets_[(map + 1) * universes_.size() - 1]; - int last_univ = universes_.back(); + 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); } @@ -117,6 +123,7 @@ int32_t Lattice::fill_offset_table(int32_t offset, int32_t target_univ_id, offsets_[map * universes_.size() + it.indx_] = offset; offset += count_universe_instances(*it, target_univ_id, univ_count_memo); } + return offset; } @@ -225,14 +232,14 @@ RectLattice::RectLattice(pugi::xml_node lat_node) : Lattice {lat_node} //============================================================================== -int32_t const& RectLattice::operator[](array const& i_xyz) +const int32_t& RectLattice::operator[](const array& i_xyz) { return universes_[get_flat_index(i_xyz)]; } //============================================================================== -bool RectLattice::are_valid_indices(array const& i_xyz) const +bool RectLattice::are_valid_indices(const array& i_xyz) const { return ((i_xyz[0] >= 0) && (i_xyz[0] < n_cells_[0]) && (i_xyz[1] >= 0) && (i_xyz[1] < n_cells_[1]) && (i_xyz[2] >= 0) && @@ -354,7 +361,7 @@ Position RectLattice::get_local_position( //============================================================================== -int32_t& RectLattice::offset(int map, array const& i_xyz) +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] + @@ -676,13 +683,19 @@ void HexLattice::fill_lattice_y(const vector& univ_words) //============================================================================== -int32_t const& HexLattice::operator[](array const& i_xyz) +const int32_t& HexLattice::operator[](const array& i_xyz) { return universes_[get_flat_index(i_xyz)]; } //============================================================================== +// The HexLattice iterators need their own versions b/c the universes array is +// "square", meaning that it is allocated with entries that are intentionally +// left empty. As such, the iterator indices need to skip the empty entries to +// get cell instances and geometry paths correct. See the image in the Theory +// and Methodology section on "Hexagonal Lattice Indexing" for a visual of where +// the empty positions are. LatticeIter HexLattice::begin() { return LatticeIter(*this, n_rings_ - 1); @@ -693,9 +706,24 @@ ReverseLatticeIter HexLattice::rbegin() return ReverseLatticeIter(*this, universes_.size() - n_rings_); } +int32_t& HexLattice::back() +{ + return universes_[universes_.size() - n_rings_]; +} + +LatticeIter HexLattice::end() +{ + return LatticeIter(*this, universes_.size() - n_rings_ + 1); +} + +ReverseLatticeIter HexLattice::rend() +{ + return ReverseLatticeIter(*this, n_rings_ - 2); +} + //============================================================================== -bool HexLattice::are_valid_indices(array const& i_xyz) const +bool HexLattice::are_valid_indices(const array& 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) && @@ -992,7 +1020,7 @@ bool HexLattice::is_valid_index(int indx) const //============================================================================== -int32_t& HexLattice::offset(int map, array const& i_xyz) +int32_t& HexLattice::offset(int map, const array& i_xyz) { int nx {2 * n_rings_ - 1}; int ny {2 * n_rings_ - 1}; diff --git a/tests/unit_tests/cell_instances/test_hex_multilattice.py b/tests/unit_tests/cell_instances/test_hex_multilattice.py new file mode 100644 index 0000000000..3f503fe9b8 --- /dev/null +++ b/tests/unit_tests/cell_instances/test_hex_multilattice.py @@ -0,0 +1,119 @@ +from math import sqrt + +import pytest +import numpy as np +import openmc +import openmc.lib + +from tests import cdtemp + + +@pytest.fixture(scope='module', autouse=True) +def double_hex_lattice_model(): + openmc.reset_auto_ids() + radius = 0.9 + pin_lattice_pitch = 2.0 + # make the hex prism a little larger to make sure test + # locations are definitively in the model + hex_prism_edge = 1.2 * pin_lattice_pitch + + model = openmc.Model() + + # materials + nat_u = openmc.Material() + nat_u.set_density('g/cm3', 12.0) + nat_u.add_element('U', 1.0) + + graphite = openmc.Material() + graphite.set_density('g/cm3', 1.1995) + graphite.add_element('C', 1.0) + + # zplanes to define lower and upper region + z_low = openmc.ZPlane(-10, boundary_type='vacuum') + z_mid = openmc.ZPlane(0) + z_high = openmc.ZPlane(10, boundary_type='vacuum') + hex_prism = openmc.model.HexagonalPrism( + edge_length=hex_prism_edge, boundary_type='reflective') + + # geometry + cyl = openmc.ZCylinder(r=radius) + univ = openmc.model.pin([cyl], [nat_u, graphite]) + + # create a hexagonal lattice of compacts + hex_lattice = openmc.HexLattice() + hex_lattice.orientation = 'y' + hex_lattice.pitch = (pin_lattice_pitch,) + hex_lattice.center = (0., 0.) + center = [univ] + ring = [univ, univ, univ, univ, univ, univ] + hex_lattice.universes = [ring, center] + lower_hex_cell = openmc.Cell(fill=hex_lattice, region=-hex_prism & +z_low & -z_mid) + upper_hex_cell = openmc.Cell(fill=hex_lattice, region=-hex_prism & +z_mid & -z_high) + hex_cells = [lower_hex_cell, upper_hex_cell] + model.geometry = openmc.Geometry(hex_cells) + + # moderator + cell = next(iter(univ.get_all_cells().values())) + tally = openmc.Tally(tally_id=1) + filter = openmc.DistribcellFilter(cell) + tally.filters = [filter] + tally.scores = ['flux'] + model.tallies = [tally] + + # settings + # source definition. fission source given bounding box of graphite active region + system_LL = (-pin_lattice_pitch*sqrt(3)/2, -pin_lattice_pitch, -5) + system_UR = (pin_lattice_pitch*sqrt(3)/2, pin_lattice_pitch, 5) + source_dist = openmc.stats.Box(system_LL, system_UR) + model.settings.source = openmc.IndependentSource(space=source_dist) + model.settings.particles = 100 + model.settings.inactive = 2 + model.settings.batches = 10 + + with cdtemp(): + model.export_to_xml() + openmc.lib.init() + yield + openmc.lib.finalize() + + +# Lower cell instances +# 6 +# 5 4 +# 3 +# 2 1 +# 0 +# Upper cell instances +# 13 +# 12 11 +# 10 +# 9 8 +# 7 +hex_expected_results = [ + ((0.0, -2.0, -5.0), 0), + ((1.732, -1.0, -5.0), 1), + ((-1.732, -1.0, -5.0), 2), + ((0.0, 0.0, -0.1), 3), + ((1.732, 1.0, -5.0), 4), + ((-1.732, 1.0, -5.0), 5), + ((0.0, 2.0, -0.1), 6), + ((0.0, -2.0, 5.0), 7), + ((1.732, -1.0, 5.0), 8), + ((-1.732, -1.0, 5.0), 9), + ((0.0, 0.0, 5.0), 10), + ((1.732, 1.0, 5.0), 11), + ((-1.732, 1.0, 5.0), 12), + ((0.0, 2.0, 5.0), 13), +] + + +@pytest.mark.parametrize("r,expected_cell_instance", hex_expected_results, ids=str) +def test_cell_instance_hex_multilattice(r, expected_cell_instance): + _, cell_instance = openmc.lib.find_cell(r) + assert cell_instance == expected_cell_instance + + +def test_cell_instance_multilattice_results(): + openmc.lib.run() + tally_results = openmc.lib.tallies[1].mean + assert (tally_results != 0.0).all() diff --git a/tests/unit_tests/test_cell_instance.py b/tests/unit_tests/cell_instances/test_rect_multilattice.py similarity index 84% rename from tests/unit_tests/test_cell_instance.py rename to tests/unit_tests/cell_instances/test_rect_multilattice.py index 25c20cfef2..aaecb3bdac 100644 --- a/tests/unit_tests/test_cell_instance.py +++ b/tests/unit_tests/cell_instances/test_rect_multilattice.py @@ -1,5 +1,5 @@ -import numpy as np import pytest +import numpy as np import openmc import openmc.lib @@ -8,7 +8,7 @@ from tests import cdtemp @pytest.fixture(scope='module', autouse=True) -def double_lattice_model(): +def double_rect_lattice_model(): openmc.reset_auto_ids() model = openmc.Model() @@ -40,8 +40,9 @@ def double_lattice_model(): cell_with_lattice2.translation = (2., 0., 0.) model.geometry = openmc.Geometry([cell_with_lattice1, cell_with_lattice2]) - tally = openmc.Tally() - tally.filters = [openmc.DistribcellFilter(c)] + tally = openmc.Tally(tally_id=1) + dcell_filter = openmc.DistribcellFilter(c) + tally.filters = [dcell_filter] tally.scores = ['flux'] model.tallies = [tally] @@ -50,7 +51,8 @@ def double_lattice_model(): bbox[0][2] = -0.5 bbox[1][2] = 0.5 space = openmc.stats.Box(*bbox) - model.settings.source = openmc.IndependentSource(space=space) + source = openmc.IndependentSource(space=space) + model.settings.source = source # Add necessary settings and export model.settings.batches = 10 @@ -63,14 +65,13 @@ def double_lattice_model(): yield openmc.lib.finalize() - # This shows the expected cell instance numbers for each lattice position: # ┌─┬─┬─┬─┐ # │2│3│6│7│ # ├─┼─┼─┼─┤ # │0│1│4│5│ # └─┴─┴─┴─┘ -expected_results = [ +rect_expected_results = [ ((0.5, 0.5, 0.0), 0), ((1.5, 0.5, 0.0), 1), ((0.5, 1.5, 0.0), 2), @@ -80,13 +81,16 @@ expected_results = [ ((2.5, 1.5, 0.0), 6), ((3.5, 1.5, 0.0), 7), ] -@pytest.mark.parametrize("r,expected_cell_instance", expected_results) -def test_cell_instance_multilattice(r, expected_cell_instance): + + +@pytest.mark.parametrize("r,expected_cell_instance", rect_expected_results, ids=lambda p : f'{p}') +def test_cell_instance_rect_multilattice(r, expected_cell_instance): _, cell_instance = openmc.lib.find_cell(r) assert cell_instance == expected_cell_instance def test_cell_instance_multilattice_results(): + openmc.run() openmc.lib.run() tally_results = openmc.lib.tallies[1].mean assert (tally_results != 0.0).all()