Skip to content

Commit

Permalink
Hexagonal lattice iterators (#2921)
Browse files Browse the repository at this point in the history
Co-authored-by: Paul Romano <[email protected]>
  • Loading branch information
pshriwise and paulromano authored Jun 17, 2024
1 parent 89d4daf commit ce4176d
Show file tree
Hide file tree
Showing 5 changed files with 192 additions and 33 deletions.
35 changes: 21 additions & 14 deletions include/openmc/lattice.h
Original file line number Diff line number Diff line change
Expand Up @@ -56,13 +56,14 @@ class Lattice {

virtual ~Lattice() {}

virtual int32_t const& operator[](array<int, 3> const& i_xyz) = 0;
virtual const int32_t& operator[](const array<int, 3>& 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();
Expand All @@ -81,7 +82,7 @@ class Lattice {
//! \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<int, 3> const& i_xyz) const = 0;
virtual bool are_valid_indices(const array<int, 3>& i_xyz) const = 0;

//! \brief Find the next lattice surface crossing
//! \param r A 3D Cartesian coordinate.
Expand Down Expand Up @@ -125,7 +126,7 @@ 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, array<int, 3> const& i_xyz) = 0;
virtual int32_t& offset(int map, const array<int, 3>& i_xyz) = 0;

//! \brief Get the distribcell offset for a lattice tile.
//! \param The map index for the target cell.
Expand Down Expand Up @@ -167,12 +168,12 @@ class LatticeIter {

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;
}

Expand All @@ -190,7 +191,7 @@ class ReverseLatticeIter : public LatticeIter {

ReverseLatticeIter& operator++()
{
while (indx_ > -1) {
while (indx_ > lat_.begin().indx_ - 1) {
--indx_;
if (lat_.is_valid_index(indx_))
return *this;
Expand All @@ -206,9 +207,9 @@ class RectLattice : public Lattice {
public:
explicit RectLattice(pugi::xml_node lat_node);

int32_t const& operator[](array<int, 3> const& i_xyz) override;
const int32_t& operator[](const array<int, 3>& i_xyz) override;

bool are_valid_indices(array<int, 3> const& i_xyz) const override;
bool are_valid_indices(const array<int, 3>& i_xyz) const override;

std::pair<double, array<int, 3>> distance(
Position r, Direction u, const array<int, 3>& i_xyz) const override;
Expand All @@ -221,7 +222,7 @@ class RectLattice : public Lattice {
Position get_local_position(
Position r, const array<int, 3>& i_xyz) const override;

int32_t& offset(int map, array<int, 3> const& i_xyz) override;
int32_t& offset(int map, const array<int, 3>& i_xyz) override;

int32_t offset(int map, int indx) const override;

Expand All @@ -241,13 +242,19 @@ class HexLattice : public Lattice {
public:
explicit HexLattice(pugi::xml_node lat_node);

int32_t const& operator[](array<int, 3> const& i_xyz) override;
const int32_t& operator[](const array<int, 3>& i_xyz) override;

LatticeIter begin() override;

ReverseLatticeIter rbegin() override;

bool are_valid_indices(array<int, 3> const& i_xyz) const override;
LatticeIter end() override;

int32_t& back() override;

ReverseLatticeIter rend() override;

bool are_valid_indices(const array<int, 3>& i_xyz) const override;

std::pair<double, array<int, 3>> distance(
Position r, Direction u, const array<int, 3>& i_xyz) const override;
Expand All @@ -262,7 +269,7 @@ class HexLattice : public Lattice {

bool is_valid_index(int indx) const override;

int32_t& offset(int map, array<int, 3> const& i_xyz) override;
int32_t& offset(int map, const array<int, 3>& i_xyz) override;

int32_t offset(int map, int indx) const override;

Expand Down
3 changes: 2 additions & 1 deletion src/geometry_aux.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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_;
Expand Down
46 changes: 37 additions & 9 deletions src/lattice.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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);
}
Expand All @@ -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;
}

Expand Down Expand Up @@ -225,14 +232,14 @@ RectLattice::RectLattice(pugi::xml_node lat_node) : Lattice {lat_node}

//==============================================================================

int32_t const& RectLattice::operator[](array<int, 3> const& i_xyz)
const int32_t& RectLattice::operator[](const array<int, 3>& i_xyz)
{
return universes_[get_flat_index(i_xyz)];
}

//==============================================================================

bool RectLattice::are_valid_indices(array<int, 3> const& i_xyz) const
bool RectLattice::are_valid_indices(const array<int, 3>& 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) &&
Expand Down Expand Up @@ -354,7 +361,7 @@ Position RectLattice::get_local_position(

//==============================================================================

int32_t& RectLattice::offset(int map, array<int, 3> const& i_xyz)
int32_t& RectLattice::offset(int map, const array<int, 3>& i_xyz)
{
return offsets_[n_cells_[0] * n_cells_[1] * n_cells_[2] * map +
n_cells_[0] * n_cells_[1] * i_xyz[2] +
Expand Down Expand Up @@ -676,13 +683,19 @@ void HexLattice::fill_lattice_y(const vector<std::string>& univ_words)

//==============================================================================

int32_t const& HexLattice::operator[](array<int, 3> const& i_xyz)
const int32_t& HexLattice::operator[](const array<int, 3>& 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);
Expand All @@ -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<int, 3> const& i_xyz) const
bool HexLattice::are_valid_indices(const array<int, 3>& 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) &&
Expand Down Expand Up @@ -992,7 +1020,7 @@ bool HexLattice::is_valid_index(int indx) const

//==============================================================================

int32_t& HexLattice::offset(int map, array<int, 3> const& i_xyz)
int32_t& HexLattice::offset(int map, const array<int, 3>& i_xyz)
{
int nx {2 * n_rings_ - 1};
int ny {2 * n_rings_ - 1};
Expand Down
119 changes: 119 additions & 0 deletions tests/unit_tests/cell_instances/test_hex_multilattice.py
Original file line number Diff line number Diff line change
@@ -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()
Loading

0 comments on commit ce4176d

Please sign in to comment.