From 4c919e5e3e97ae9a7ce6afaece411921f0def551 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 12:47:56 -0700 Subject: [PATCH 01/16] Move grids.py to own module armi.reactor.grids Goal is to introduce separate files for each grid type, but still make those importable from armi.reactor.grids --- armi/reactor/{grids.py => grids/__init__.py} | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename armi/reactor/{grids.py => grids/__init__.py} (99%) diff --git a/armi/reactor/grids.py b/armi/reactor/grids/__init__.py similarity index 99% rename from armi/reactor/grids.py rename to armi/reactor/grids/__init__.py index 9ce755450..2ec5b4a24 100644 --- a/armi/reactor/grids.py +++ b/armi/reactor/grids/__init__.py @@ -1,4 +1,4 @@ -# Copyright 2019 TerraPower, LLC +# Copyright 2023 TerraPower, LLC # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. From ddf372bf43edf9f3c49983052f1fc77aeac1a128 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 12:53:04 -0700 Subject: [PATCH 02/16] Add armi.reactor.grids.constants Intended to share some constants and type hints across the `armi.reactor.grids` submodule. Intentionally sparse. Items will be imported into `armi.reactor.grids` to help with back compatibility. --- armi/reactor/grids/__init__.py | 12 ++++++++---- armi/reactor/grids/constants.py | 19 +++++++++++++++++++ 2 files changed, 27 insertions(+), 4 deletions(-) create mode 100644 armi/reactor/grids/constants.py diff --git a/armi/reactor/grids/__init__.py b/armi/reactor/grids/__init__.py index 2ec5b4a24..f6c0dc0a6 100644 --- a/armi/reactor/grids/__init__.py +++ b/armi/reactor/grids/__init__.py @@ -77,10 +77,7 @@ ("unitSteps", "bounds", "unitStepLimits", "offset", "geomType", "symmetry"), ) TAU = math.pi * 2.0 -BOUNDARY_0_DEGREES = 1 -BOUNDARY_60_DEGREES = 2 -BOUNDARY_120_DEGREES = 3 -BOUNDARY_CENTER = 4 + COS30 = math.sqrt(3) / 2.0 SIN30 = 1.0 / 2.0 @@ -96,6 +93,13 @@ ] ) +from .constants import ( + BOUNDARY_CENTER, + BOUNDARY_0_DEGREES, + BOUNDARY_120_DEGREES, + BOUNDARY_60_DEGREES, +) + class LocationBase: """ diff --git a/armi/reactor/grids/constants.py b/armi/reactor/grids/constants.py new file mode 100644 index 000000000..f2fa78bd8 --- /dev/null +++ b/armi/reactor/grids/constants.py @@ -0,0 +1,19 @@ +# Copyright 2023 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +"""Some constants often used in grid manipulation""" + +BOUNDARY_0_DEGREES = 1 +BOUNDARY_60_DEGREES = 2 +BOUNDARY_120_DEGREES = 3 +BOUNDARY_CENTER = 4 From fbca33dc8fb379e8243f6c20c7faaadab4e38bc1 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 12:56:21 -0700 Subject: [PATCH 03/16] Add armi.reactor.grids.locations Provides abstract `LocationBase` as the interface defining base class. Concrete subclasses `IndexLocation`, `MultiIndexLocation`, and `CoordinateLocation` are also defined. There is little to no implementation differences in this change. Only that `LocationBase` is a proper abstract base class, so subclasses that don't override the necessary abstract methods will need to change. Some other key changes - `LocationBase.__eq__` returns `NotImplemented` if the RHS argument is not a tuple nor a `LocationBase` The four main classes defined are imported into `armi.reactor.grids` to maintain back compatibility. --- armi/reactor/grids/__init__.py | 417 +-------------------------- armi/reactor/grids/locations.py | 480 ++++++++++++++++++++++++++++++++ 2 files changed, 486 insertions(+), 411 deletions(-) create mode 100644 armi/reactor/grids/locations.py diff --git a/armi/reactor/grids/__init__.py b/armi/reactor/grids/__init__.py index f6c0dc0a6..54270d5ab 100644 --- a/armi/reactor/grids/__init__.py +++ b/armi/reactor/grids/__init__.py @@ -100,417 +100,12 @@ BOUNDARY_60_DEGREES, ) - -class LocationBase: - """ - A namedtuple-like object for storing location information. - - It's immutable (you can't set things after construction) and has names. - - Notes - ----- - This was originally a namedtuple but there was a somewhat unbelievable - bug in Python 2.7.8 where unpickling a reactor full of these ended up - inexplicably replacing one of them with an AssemblyParameterCollection. - The bug did not show up in Python 3. - - Converting to this class solved that problem. - """ - - __slots__ = ("_i", "_j", "_k", "_grid") - - def __init__(self, i, j, k, grid): - self._i = i - self._j = j - self._k = k - self._grid = grid - - def __repr__(self): - return "<{} @ ({},{:},{})>".format( - self.__class__.__name__, self.i, self.j, self.k - ) - - def __getstate__(self): - """Used in pickling and deepcopy, this detaches the grid.""" - return (self._i, self._j, self._k, None) - - def __setstate__(self, state): - """ - Unpickle a locator, the grid will attach itself if it was also pickled, otherwise this will - be detached. - """ - self.__init__(*state) - - @property - def i(self): - return self._i - - @property - def j(self): - return self._j - - @property - def k(self): - return self._k - - @property - def grid(self): - return self._grid - - def __getitem__(self, index): - return (self.i, self.j, self.k, self.grid)[index] - - def __hash__(self): - """ - Define a hash so we can use these as dict keys w/o having exact object. - - Notes - ----- - Including the ``grid`` attribute may be more robust; however, using only (i, j, k) allows - dictionaries to use IndexLocations and (i,j,k) tuples interchangeably. - """ - return hash((self.i, self.j, self.k)) - - def __eq__(self, other): - if isinstance(other, tuple): - return (self.i, self.j, self.k) == other - return ( - self.i == other.i - and self.j == other.j - and self.k == other.k - and self.grid is other.grid - ) - - def __lt__(self, that): - """ - A Locationbase is less than another if the pseudo-radius is less, or if equal, in order - any index is less. - - Examples - -------- - >>> grid = grids.HexGrid.fromPitch(1.0) - >>> grid[0, 0, 0] < grid[2, 3, 4] # the "radius" is less - True - >>> grid[2, 3, 4] < grid[2, 3, 4] # they are equal - False - >>> grid[2, 3, 4] < grid[-2, 3, 4] # 2 is greater than -2 - False - >>> grid[-2, 3, 4] < grid[2, 3, 4] # -2 is less than 2 - True - >>> grid[1, 3, 4] < grid[-2, 3, 4] # the "radius" is less - True - """ - selfIndices = self.indices - thatIndices = that.indices - # this is not really r, but it is fast and consistent - selfR = abs(selfIndices).sum() - thatR = abs(thatIndices).sum() - - # this cannot be reduced to - # return selfR < thatR or (selfIndices < thatIndices).any() - # because the comparison is not symmetric. - if selfR < thatR: - return True - else: - for lt, eq in zip(selfIndices < thatIndices, selfIndices == thatIndices): - if eq: - continue - - return lt - - return False - - def __len__(self): - """Returns 3, the number of directions.""" - return 3 - - def associate(self, grid): - """Re-assign locator to another Grid.""" - self._grid = grid - - -class IndexLocation(LocationBase): - """ - An immutable location representing one cell in a grid. - - The locator is intimately tied to a grid and together, they represent - a grid cell somewhere in the coordinate system of the grid. - - ``grid`` is not in the constructor (must be added after construction ) because - the extra argument (grid) gives an inconsistency between __init__ and __new__. - Unfortunately this decision makes whipping up IndexLocations on the fly awkward. - But perhaps that's ok because they should only be created by their grids. - """ - - __slots__ = () - - def __add__(self, that): - """ - Enable adding with other objects like this and/or 3-tuples. - - Tuples are needed so we can terminate the recursive additions with a (0,0,0) basis. - """ - # New location is not associated with any particular grid. - return self.__class__( - self[0] + that[0], self[1] + that[1], self[2] + that[2], None - ) - - def __sub__(self, that): - return self.__class__( - self[0] - that[0], self[1] - that[1], self[2] - that[2], None - ) - - def detachedCopy(self): - """ - Make a copy of this locator that is not associated with a grid. - - See Also - -------- - armi.reactor.reactors.detach : uses this - """ - return self.__class__(self.i, self.j, self.k, None) - - @property - def parentLocation(self): - """ - Get the spatialLocator of the ArmiObject that this locator's grid is anchored to. - - For example, if this is one of many spatialLocators in a 2-D grid representing - a reactor, then the ``parentLocation`` is the spatialLocator of the reactor, which - will often be a ``CoordinateLocation``. - """ - grid = self.grid # performance matters a lot here so we remove a dot - # check for None rather than __nonzero__ for speed (otherwise it checks the length) - if ( - grid is not None - and grid.armiObject is not None - and grid.armiObject.parent is not None - ): - return grid.armiObject.spatialLocator - return None - - @property - def indices(self): - """ - Get the non-grid indices (i,j,k) of this locator. - - This strips off the annoying ``grid`` tagalong which is there to ensure proper - equality (i.e. (0,0,0) in a storage rack is not equal to (0,0,0) in a core). - - It is a numpy array for two reasons: - - 1. It can be added and subtracted for the recursive computations - through different coordinate systems - 2. It can be written/read from the database. - - """ - return numpy.array(self[:3]) - - def getCompleteIndices(self) -> Tuple[int, int, int]: - """ - Transform the indices of this object up to the top mesh. - - The top mesh is either the one where there's no more parent (true top) - or when an axis gets added twice. Unlike with coordinates, - you can only add each index axis one time. Thus a *complete* - set of indices is one where an index for each axis has been defined - by a set of 1, 2, or 3 nested grids. - - This is useful for getting the reactor-level (i,j,k) indices of an object - in a multi-layered 2-D(assemblies in core)/1-D(blocks in assembly) mesh - like the one mapping blocks up to reactor in Hex reactors. - - The benefit of that particular mesh over a 3-D one is that different - assemblies can have different axial meshes, a common situation. - - It will just return local indices for pin-meshes inside of blocks. - - A tuple is returned so that it is easy to compare pairs of indices. - """ - parentLocation = self.parentLocation # to avoid evaluating property if's twice - indices = self.indices - if parentLocation is not None: - if parentLocation.grid is not None and addingIsValid( - self.grid, parentLocation.grid - ): - indices += parentLocation.indices - return tuple(indices) - - def getLocalCoordinates(self, nativeCoords=False): - """Return the coordinates of the center of the mesh cell here in cm.""" - if self.grid is None: - raise ValueError( - "Cannot get local coordinates of {} because grid is None.".format(self) - ) - return self.grid.getCoordinates(self.indices, nativeCoords=nativeCoords) - - def getGlobalCoordinates(self, nativeCoords=False): - """Get coordinates in global 3D space of the centroid of this object.""" - parentLocation = self.parentLocation # to avoid evaluating property if's twice - if parentLocation: - return self.getLocalCoordinates( - nativeCoords=nativeCoords - ) + parentLocation.getGlobalCoordinates(nativeCoords=nativeCoords) - return self.getLocalCoordinates(nativeCoords=nativeCoords) - - def getGlobalCellBase(self): - """Return the cell base (i.e. "bottom left"), in global coordinate system.""" - parentLocation = self.parentLocation # to avoid evaluating property if's twice - if parentLocation: - return parentLocation.getGlobalCellBase() + self.grid.getCellBase( - self.indices - ) - return self.grid.getCellBase(self.indices) - - def getGlobalCellTop(self): - """Return the cell top (i.e. "top right"), in global coordinate system.""" - parentLocation = self.parentLocation # to avoid evaluating property if's twice - if parentLocation: - return parentLocation.getGlobalCellTop() + self.grid.getCellTop( - self.indices - ) - return self.grid.getCellTop(self.indices) - - def getRingPos(self): - """Return ring and position of this locator.""" - return self.grid.getRingPos(self.getCompleteIndices()) - - def getSymmetricEquivalents(self): - """ - Get symmetrically-equivalent locations, based on Grid symmetry. - - See Also - -------- - Grid.getSymmetricEquivalents - """ - return self.grid.getSymmetricEquivalents(self.indices) - - def distanceTo(self, other) -> float: - """Return the distance from this locator to another.""" - return math.sqrt( - ( - ( - numpy.array(self.getGlobalCoordinates()) - - numpy.array(other.getGlobalCoordinates()) - ) - ** 2 - ).sum() - ) - - -class MultiIndexLocation(IndexLocation): - """ - A collection of index locations that can be used as a spatialLocator. - - This allows components with multiplicity>1 to have location information - within a parent grid. The implication is that there are multiple - discrete components, each one residing in one of the actual locators - underlying this collection. - - This class contains an implementation that allows a multi-index - location to be used in the ARMI data model similar to a - individual IndexLocation. - """ - - # MIL's cannot be hashed, so we need to scrape off the implementation from - # LocationBase. This raises some interesting questions of substitutability of the - # various Location classes, which should be addressed. - __hash__ = None - - def __init__(self, grid): - IndexLocation.__init__(self, 0, 0, 0, grid) - self._locations = [] - - def __getstate__(self): - """Used in pickling and deepcopy, this detaches the grid.""" - return self._locations - - def __setstate__(self, state): - """ - Unpickle a locator, the grid will attach itself if it was also pickled, otherwise this will - be detached. - """ - self.__init__(None) - self._locations = state - - def __repr__(self): - return "<{} with {} locations>".format( - self.__class__.__name__, len(self._locations) - ) - - def __getitem__(self, index): - return self._locations[index] - - def __setitem__(self, index, obj): - self._locations[index] = obj - - def __iter__(self): - return iter(self._locations) - - def __len__(self): - return len(self._locations) - - def detachedCopy(self): - loc = MultiIndexLocation(None) - loc.extend(self._locations) - return loc - - def associate(self, grid): - self._grid = grid - for loc in self._locations: - loc.associate(grid) - - def getCompleteIndices(self) -> Tuple[int, int, int]: - raise NotImplementedError("Multi locations cannot do this yet.") - - def append(self, location: IndexLocation): - self._locations.append(location) - - def extend(self, locations: List[IndexLocation]): - self._locations.extend(locations) - - def pop(self, location: IndexLocation): - self._locations.pop(location) - - @property - def indices(self): - """ - Return indices for all locations. - - Notes - ----- - Notice that this returns a list of all of the indices, unlike the ``indices()`` - implementation for :py:class:`IndexLocation`. This is intended to make the - behavior of getting the indices from the Locator symmetric with passing a list - of indices to the Grid's ``__getitem__()`` function, which constructs and - returns a ``MultiIndexLocation`` containing those indices. - """ - return [loc.indices for loc in self._locations] - - -class CoordinateLocation(IndexLocation): - """ - A triple representing a point in space. - - This is still associated with a grid. The grid defines the continuous coordinate - space and axes that the location is within. This also links to the composite tree. - """ - - __slots__ = () - - def getLocalCoordinates(self, nativeCoords=False): - """Return x,y,z coordinates in cm within the grid's coordinate system.""" - return self.indices - - def getCompleteIndices(self): - """Top of chain. Stop recursion and return basis.""" - return 0, 0, 0 - - def getGlobalCellBase(self): - return self.indices - - def getGlobalCellTop(self): - return self.indices +from .locations import ( + LocationBase, + IndexLocation, + MultiIndexLocation, + CoordinateLocation, +) class Grid: diff --git a/armi/reactor/grids/locations.py b/armi/reactor/grids/locations.py new file mode 100644 index 000000000..854b1a1ac --- /dev/null +++ b/armi/reactor/grids/locations.py @@ -0,0 +1,480 @@ +# Copyright 2023 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from typing import Optional, TYPE_CHECKING, Union, Hashable, Tuple, List, Iterator +from abc import ABC, abstractmethod +import math + +import numpy + +if TYPE_CHECKING: + # Avoid some circular imports + from . import Grid + + +IJType = Tuple[int, int] +IJKType = Tuple[int, int, int] + + +class LocationBase(ABC): + """ + A namedtuple-like object for storing location information. + + It's immutable (you can't set things after construction) and has names. + + Notes + ----- + This was originally a namedtuple but there was a somewhat unbelievable + bug in Python 2.7.8 where unpickling a reactor full of these ended up + inexplicably replacing one of them with an AssemblyParameterCollection. + The bug did not show up in Python 3. + + Converting to this class solved that problem. + """ + + __slots__ = ("_i", "_j", "_k", "_grid") + + def __init__(self, i: int, j: int, k: int, grid: Optional["Grid"]): + self._i = i + self._j = j + self._k = k + self._grid = grid + + def __repr__(self) -> str: + return "<{} @ ({},{:},{})>".format( + self.__class__.__name__, self.i, self.j, self.k + ) + + def __getstate__(self) -> Hashable: + """ + Used in pickling and deepcopy, this detaches the grid. + """ + return (self._i, self._j, self._k, None) + + def __setstate__(self, state: Hashable): + """ + Unpickle a locator, the grid will attach itself if it was also pickled, otherwise this will + be detached. + """ + self.__init__(*state) + + @property + def i(self) -> int: + return self._i + + @property + def j(self) -> int: + return self._j + + @property + def k(self) -> int: + return self._k + + @property + def grid(self) -> Optional["Grid"]: + return self._grid + + def __getitem__(self, index: int) -> Union[int, "Grid"]: + return (self.i, self.j, self.k, self.grid)[index] + + def __hash__(self) -> Hashable: + """ + Define a hash so we can use these as dict keys w/o having exact object. + + Notes + ----- + Including the ``grid`` attribute may be more robust; however, using only (i, j, k) allows + dictionaries to use IndexLocations and (i,j,k) tuples interchangeably. + """ + return hash((self.i, self.j, self.k)) + + def __eq__(self, other: Union[Tuple[int, int, int], "LocationBase"]) -> bool: + if isinstance(other, tuple): + return (self.i, self.j, self.k) == other + if isinstance(other, LocationBase): + return ( + self.i == other.i + and self.j == other.j + and self.k == other.k + and self.grid is other.grid + ) + return NotImplemented + + def __lt__(self, that: "LocationBase") -> bool: + """ + A Locationbase is less than another if the pseudo-radius is less, or if equal, in order + any index is less. + + Examples + -------- + >>> grid = grids.HexGrid.fromPitch(1.0) + >>> grid[0, 0, 0] < grid[2, 3, 4] # the "radius" is less + True + >>> grid[2, 3, 4] < grid[2, 3, 4] # they are equal + False + >>> grid[2, 3, 4] < grid[-2, 3, 4] # 2 is greater than -2 + False + >>> grid[-2, 3, 4] < grid[2, 3, 4] # -2 is less than 2 + True + >>> grid[1, 3, 4] < grid[-2, 3, 4] # the "radius" is less + True + """ + selfIndices = self.indices + thatIndices = that.indices + # this is not really r, but it is fast and consistent + selfR = abs(selfIndices).sum() + thatR = abs(thatIndices).sum() + + # this cannot be reduced to + # return selfR < thatR or (selfIndices < thatIndices).any() + # because the comparison is not symmetric. + if selfR < thatR: + return True + else: + for lt, eq in zip(selfIndices < thatIndices, selfIndices == thatIndices): + if eq: + continue + + return lt + + return False + + def __len__(self) -> int: + """Returns 3, the number of directions.""" + return 3 + + def associate(self, grid: "Grid"): + """Re-assign locator to another Grid.""" + self._grid = grid + + @property + @abstractmethod + def indices(self) -> numpy.ndarray: + """Get the non-grid indices (i,j,k) of this locator. + + This strips off the annoying ``grid`` tagalong which is there to ensure proper + equality (i.e. (0,0,0) in a storage rack is not equal to (0,0,0) in a core). + + It is a numpy array for two reasons: + + 1. It can be added and subtracted for the recursive computations + through different coordinate systems + 2. It can be written/read from the database. + """ + + +class IndexLocation(LocationBase): + """ + An immutable location representing one cell in a grid. + + The locator is intimately tied to a grid and together, they represent + a grid cell somewhere in the coordinate system of the grid. + + ``grid`` is not in the constructor (must be added after construction ) because + the extra argument (grid) gives an inconsistency between __init__ and __new__. + Unfortunately this decision makes whipping up IndexLocations on the fly awkward. + But perhaps that's ok because they should only be created by their grids. + + TODO Is the above correct still? The constructor has an optional ``Grid`` + + """ + + # TODO Maybe __slots__ = LocationBase.__slots__ + ("parentLocation", ) + # But parentLocation is a property... + # Maybe other parts of ARMI set attributes? + __slots__ = () + + def __add__(self, that: Union[IJKType, "IndexLocation"]) -> "IndexLocation": + """ + Enable adding with other objects like this and/or 3-tuples. + + Tuples are needed so we can terminate the recursive additions with a (0,0,0) basis. + """ + # New location is not associated with any particular grid. + return self.__class__( + self[0] + that[0], self[1] + that[1], self[2] + that[2], None + ) + + def __sub__(self, that: Union[IJKType, "IndexLocation"]) -> "IndexLocation": + return self.__class__( + self[0] - that[0], self[1] - that[1], self[2] - that[2], None + ) + + def detachedCopy(self) -> "IndexLocation": + """ + Make a copy of this locator that is not associated with a grid. + + See Also + -------- + armi.reactor.reactors.detach : uses this + """ + return self.__class__(self.i, self.j, self.k, None) + + @property + def parentLocation(self): + """ + Get the spatialLocator of the ArmiObject that this locator's grid is anchored to. + + For example, if this is one of many spatialLocators in a 2-D grid representing + a reactor, then the ``parentLocation`` is the spatialLocator of the reactor, which + will often be a ``CoordinateLocation``. + """ + grid = self.grid # performance matters a lot here so we remove a dot + # check for None rather than __nonzero__ for speed (otherwise it checks the length) + if ( + grid is not None + and grid.armiObject is not None + and grid.armiObject.parent is not None + ): + return grid.armiObject.spatialLocator + return None + + @property + def indices(self) -> numpy.ndarray: + """ + Get the non-grid indices (i,j,k) of this locator. + + This strips off the annoying ``grid`` tagalong which is there to ensure proper + equality (i.e. (0,0,0) in a storage rack is not equal to (0,0,0) in a core). + + It is a numpy array for two reasons: + + 1. It can be added and subtracted for the recursive computations + through different coordinate systems + 2. It can be written/read from the database. + + """ + return numpy.array(self[:3]) + + def getCompleteIndices(self) -> IJKType: + """ + Transform the indices of this object up to the top mesh. + + The top mesh is either the one where there's no more parent (true top) + or when an axis gets added twice. Unlike with coordinates, + you can only add each index axis one time. Thus a *complete* + set of indices is one where an index for each axis has been defined + by a set of 1, 2, or 3 nested grids. + + This is useful for getting the reactor-level (i,j,k) indices of an object + in a multi-layered 2-D(assemblies in core)/1-D(blocks in assembly) mesh + like the one mapping blocks up to reactor in Hex reactors. + + The benefit of that particular mesh over a 3-D one is that different + assemblies can have different axial meshes, a common situation. + + It will just return local indices for pin-meshes inside of blocks. + + A tuple is returned so that it is easy to compare pairs of indices. + """ + parentLocation = self.parentLocation # to avoid evaluating property if's twice + indices = self.indices + if parentLocation is not None: + if parentLocation.grid is not None and addingIsValid( + self.grid, parentLocation.grid + ): + indices += parentLocation.indices + return tuple(indices) + + def getLocalCoordinates(self, nativeCoords=False): + """Return the coordinates of the center of the mesh cell here in cm.""" + if self.grid is None: + raise ValueError( + "Cannot get local coordinates of {} because grid is None.".format(self) + ) + return self.grid.getCoordinates(self.indices, nativeCoords=nativeCoords) + + def getGlobalCoordinates(self, nativeCoords=False): + """Get coordinates in global 3D space of the centroid of this object.""" + parentLocation = self.parentLocation # to avoid evaluating property if's twice + if parentLocation: + return self.getLocalCoordinates( + nativeCoords=nativeCoords + ) + parentLocation.getGlobalCoordinates(nativeCoords=nativeCoords) + return self.getLocalCoordinates(nativeCoords=nativeCoords) + + def getGlobalCellBase(self): + """Return the cell base (i.e. "bottom left"), in global coordinate system.""" + parentLocation = self.parentLocation # to avoid evaluating property if's twice + if parentLocation: + return parentLocation.getGlobalCellBase() + self.grid.getCellBase( + self.indices + ) + return self.grid.getCellBase(self.indices) + + def getGlobalCellTop(self): + """Return the cell top (i.e. "top right"), in global coordinate system.""" + parentLocation = self.parentLocation # to avoid evaluating property if's twice + if parentLocation: + return parentLocation.getGlobalCellTop() + self.grid.getCellTop( + self.indices + ) + return self.grid.getCellTop(self.indices) + + def getRingPos(self): + """Return ring and position of this locator.""" + return self.grid.getRingPos(self.getCompleteIndices()) + + def getSymmetricEquivalents(self): + """ + Get symmetrically-equivalent locations, based on Grid symmetry. + + See Also + -------- + Grid.getSymmetricEquivalents + """ + return self.grid.getSymmetricEquivalents(self.indices) + + def distanceTo(self, other: "IndexLocation") -> float: + """Return the distance from this locator to another.""" + return math.sqrt( + ( + ( + numpy.array(self.getGlobalCoordinates()) + - numpy.array(other.getGlobalCoordinates()) + ) + ** 2 + ).sum() + ) + + +class MultiIndexLocation(IndexLocation): + """ + A collection of index locations that can be used as a spatialLocator. + + This allows components with multiplicity>1 to have location information + within a parent grid. The implication is that there are multiple + discrete components, each one residing in one of the actual locators + underlying this collection. + + This class contains an implementation that allows a multi-index + location to be used in the ARMI data model similar to a + individual IndexLocation. + """ + + # MIL's cannot be hashed, so we need to scrape off the implementation from + # LocationBase. This raises some interesting questions of substitutability of the + # various Location classes, which should be addressed. + __hash__ = None + + _locations: List[IndexLocation] + + def __init__(self, grid: "Grid"): + IndexLocation.__init__(self, 0, 0, 0, grid) + self._locations = [] + + def __getstate__(self) -> List[IndexLocation]: + """ + Used in pickling and deepcopy, this detaches the grid. + """ + return self._locations + + def __setstate__(self, state: List[IndexLocation]): + """ + Unpickle a locator, the grid will attach itself if it was also pickled, otherwise this will + be detached. + """ + self.__init__(None) + self._locations = state + + def __repr__(self) -> str: + return "<{} with {} locations>".format( + self.__class__.__name__, len(self._locations) + ) + + def __getitem__(self, index: int) -> IndexLocation: + return self._locations[index] + + def __setitem__(self, index: int, obj: IndexLocation): + self._locations[index] = obj + + def __iter__(self) -> Iterator[IndexLocation]: + return iter(self._locations) + + def __len__(self) -> int: + return len(self._locations) + + def detachedCopy(self) -> "MultiIndexLocation": + loc = MultiIndexLocation(None) + loc.extend(self._locations) + return loc + + def associate(self, grid: "Grid"): + self._grid = grid + for loc in self._locations: + loc.associate(grid) + + def getCompleteIndices(self) -> IJKType: + raise NotImplementedError("Multi locations cannot do this yet.") + + def append(self, location: IndexLocation): + self._locations.append(location) + + def extend(self, locations: List[IndexLocation]): + self._locations.extend(locations) + + def pop(self, location: IndexLocation): + self._locations.pop(location) + + @property + def indices(self) -> List[numpy.ndarray]: + """ + Return indices for all locations. + + Notes + ----- + Notice that this returns a list of all of the indices, unlike the ``indices()`` + implementation for :py:class:`IndexLocation`. This is intended to make the + behavior of getting the indices from the Locator symmetric with passing a list + of indices to the Grid's ``__getitem__()`` function, which constructs and + returns a ``MultiIndexLocation`` containing those indices. + """ + return [loc.indices for loc in self._locations] + + +class CoordinateLocation(IndexLocation): + """ + A triple representing a point in space. + + This is still associated with a grid. The grid defines the continuous coordinate + space and axes that the location is within. This also links to the composite tree. + """ + + __slots__ = () + + def getLocalCoordinates(self, nativeCoords=False): + """Return x,y,z coordinates in cm within the grid's coordinate system.""" + return self.indices + + def getCompleteIndices(self) -> IJKType: + """Top of chain. Stop recursion and return basis.""" + return 0, 0, 0 + + def getGlobalCellBase(self): + return self.indices + + def getGlobalCellTop(self): + return self.indices + + +def addingIsValid(myGrid: "Grid", parentGrid: "Grid"): + """ + True if adding a indices from one grid to another is considered valid. + + In ARMI we allow the addition of a 1-D axial grid with a 2-D grid. + We do not allow any other kind of adding. This enables the 2D/1D + grid layout in Assemblies/Blocks but does not allow 2D indexing + in pins to become inconsistent. + """ + return myGrid.isAxialOnly and not parentGrid.isAxialOnly From 8b765f287f0090a1d9e717d768bd4e171fcec044 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 13:00:33 -0700 Subject: [PATCH 04/16] Move reactor/tests/test_grids to grids/tests --- armi/reactor/{ => grids}/tests/test_grids.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename armi/reactor/{ => grids}/tests/test_grids.py (100%) diff --git a/armi/reactor/tests/test_grids.py b/armi/reactor/grids/tests/test_grids.py similarity index 100% rename from armi/reactor/tests/test_grids.py rename to armi/reactor/grids/tests/test_grids.py From e0bb4f0254edb3c63e13a03ab73e4a0943fb84de Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 13:08:45 -0700 Subject: [PATCH 05/16] Add armi.reactor.grids.grid for Grid base class Imported into `armi.reactor.grids` so still back compatible --- armi/reactor/grids/__init__.py | 606 +------------------------------- armi/reactor/grids/grid.py | 624 +++++++++++++++++++++++++++++++++ 2 files changed, 627 insertions(+), 603 deletions(-) create mode 100644 armi/reactor/grids/grid.py diff --git a/armi/reactor/grids/__init__.py b/armi/reactor/grids/__init__.py index 54270d5ab..1766ad1b7 100644 --- a/armi/reactor/grids/__init__.py +++ b/armi/reactor/grids/__init__.py @@ -60,22 +60,16 @@ while the word **local** refers to within the current coordinate system defined by the current grid. """ -import collections import itertools import math -from typing import Tuple, List, Optional, Sequence, Union +from typing import Tuple, Optional import numpy.linalg - from armi.utils import hexagon from armi.reactor import geometry -# data structure for database-serialization of grids -GridParameters = collections.namedtuple( - "GridParameters", - ("unitSteps", "bounds", "unitStepLimits", "offset", "geomType", "symmetry"), -) + TAU = math.pi * 2.0 @@ -107,591 +101,7 @@ CoordinateLocation, ) - -class Grid: - """ - A connected set of cells characterized by indices mapping to space and vice versa. - - The cells may be characterized by any mixture of regular repeating steps and - user-defined steps in any direction. - - For example, a 2-D hex lattice has constant, regular steps whereas a 3-D hex mesh - may have user-defined axial meshes. Similar for Cartesian, RZT, etc. - - Parameters - ---------- - unitSteps : tuple of tuples, optional - Describes the grid spatially as a function on indices. - Each tuple describes how each ``(x,y,or z)`` dimension is influenced by - ``(i,j,k)``. In other words, it is:: - - (dxi, dxj, jxk), (dyi, dyj, dyk), (dzi, dzj, dzk) - - where ``dmn`` is the distance (in cm) that dimension ``m`` will change as a - function of index ``n``. - - Unit steps are used as a generic method for defining repetitive grids in a - variety of geometries, including hexagonal and Cartesian. The tuples are not - vectors in the direction of the translation, but rather grouped by direction. If - the bounds argument is described for a direction, the bounds will be used rather - than the unit step information. The default of (0, 0, 0) makes all dimensions - insensitive to indices since the coordinates are calculated by the dot product - of this and the indices. With this default, any dimension that is desired to - change with indices should be defined with bounds. RZtheta grids are created - exclusively with bounds. - - bounds : 3-tuple - Absolute increasing bounds in cm including endpoints of a non-uniform grid. - Each item represents the boundaries in the associated direction. Use Nones when - unitSteps should be applied instead. Most useful for thetaRZ grids or other - non-uniform grids. - - unitStepLimits : 3-tuple - The limit of the steps in all three directions. This constrains step-defined - grids to be finite so we can populate them with SpatialLocator objects. - - offset : 3-tuple, optional - Offset in cm for each axis. By default the center of the (0,0,0)-th object is in - the center of the grid. Offsets can move it so that the (0,0,0)-th object can - be fully within a quadrant (i.e. in a Cartesian grid). - - armiObject : ArmiObject, optional - The ArmiObject that this grid describes. For example if it's a 1-D assembly - grid, the armiObject is the assembly. Note that ``self.armiObject.spatialGrid`` - is ``self``. - - Examples - -------- - A 2D a rectangular grid with width (x) 2 and height (y) 3 would be:: - - >>> grid = Grid(unitSteps=((2, 0, 0), (0, 3, 0),(0, 0, 0))) - - A regular hex grid with pitch 1 is:: - - >>> grid = Grid(unitSteps= ((sqrt(3)/2, 0.0, 0.0), (0.5, 1.0, 0.0), (0, 0, 0)) - - .. note:: For this unit hex the magnitude of the vector constructed using the - 0th index of each tuple is 1.0. - - Notes - ----- - Each dimension must either be defined through unitSteps or bounds. - The combination of unitSteps with bounds was settled upon after some struggle to - have one unified definition of a grid (i.e. just bounds). A hexagonal grid is - somewhat challenging to represent with bounds because the axes are not orthogonal, - so a unit-direction vector plus bounds would be required. And then the bounds would - be wasted space because they can be derived simply by unit steps. Memory efficiency - is important in this object so the compact representation of - unitSteps-when-possible, bounds-otherwise was settled upon. - - Design considerations include: - - * unitSteps are more intuitive as operations starting from the center of a cell, - particularly with hexagons and rectangles. Otherwise the 0,0 position of a hexagon - in the center of 1/3-symmetric hexagon is at the phantom bottom left of the - hexagon. - - * Users generally prefer to input mesh bounds rather than centers (e.g. starting at - 0.5 instead of 0.0 in a unit mesh is weird). - - * If we store bounds, computing bounds is simple and computing centers takes ~2x the - effort. If we store centers, it's the opposite. - - * Regardless of how we store things, we'll need a Grid that has the lower-left - assembly fully inside the problem (i.e. for full core Cartesian) as well as - another one that has the lower-left assembly half-way or quarter-way sliced off - (for 1/2, 1/4, and 1/8 symmetries). The ``offset`` parameter handles this. - - * Looking up mesh boundaries (to define a mesh in another code) is generally more - common than looking up centers (for plotting or measuring distance). - - * A grid can be anchored to the object that it is in with a backreference. This - gives it the ability to traverse the composite tree and map local to global - locations without having to duplicate the composite pattern on grids. This remains - optional so grids can be used for non-reactor-package reasons. It may seem - slightly cleaner to set the armiObject to the parent's spatialLocator itself - but the major disadvantage of this is that when an object moves, the armiObject - would have to be updated. By anchoring directly to Composite objects, the parent - is always up to date no matter where or how things get moved. - - * Unit step calculations use dot products and must not be polluted by the bound - indices. Thus we reduce the size of the unitSteps tuple accordingly. - - .. impl:: ARMI supports a number of structured mesh options. - :id: IMPL_REACTOR_MESH_0 - :links: REQ_REACTOR_MESH - """ - - def __init__( - self, - unitSteps=(0, 0, 0), - bounds=(None, None, None), - unitStepLimits=((0, 1), (0, 1), (0, 1)), - offset=None, - geomType="", - symmetry="", - armiObject=None, - ): - # these lists contain the indices representing which dimensions for which steps - # are used, or for which bounds are used. index 0 is x direction, etc. - self._boundDims = [] - self._stepDims = [] - for dimensionIndex, bound in enumerate(bounds): - if bound is None: - self._stepDims.append(dimensionIndex) - else: - self._boundDims.append(dimensionIndex) - - # numpy prefers tuples like this to do slicing on arrays - self._boundDims = (tuple(self._boundDims),) - self._stepDims = (tuple(self._stepDims),) - - unitSteps = _tuplify(unitSteps) - - self._bounds = bounds - self._unitStepLimits = _tuplify(unitStepLimits) - - # only represent unit steps in dimensions they're being used so as to not - # pollute the dot product. This may reduce the length of this from 3 to 2 or 1 - self._unitSteps = numpy.array(unitSteps)[self._stepDims] - self._offset = numpy.zeros(3) if offset is None else numpy.array(offset) - - self._locations = {} # indices -> SpatialLocator map - self.armiObject = armiObject - self.buildLocations() # locations are owned by a grid, so the grid builds them. - self._backup = None # for retainState - - (_ii, iLen), (_ji, jLen), (_ki, kLen) = self.getIndexBounds() - # True if only contains k-cells. - self.isAxialOnly = iLen == jLen == 1 and kLen > 1 - - # geometric metadata encapsulated here because it's related to the grid. - # They do not impact the grid object itself. - # Notice that these are stored using their string representations, rather than - # the GridType enum. This avoids the danger of deserializing an enum value from - # an old version of the code that may have had different numeric values. - self._geomType: str = str(geomType) - self._symmetry: str = str(symmetry) - - def reduce(self): - """ - Return the set of arguments used to create this Grid. - - This is very much like the argument tuple from ``__reduce__``, but we do not - implement ``__reduce__`` for real, because we are generally happy with - ``__getstate__`` and ``__setstate__`` for pickling purposes. However, getting - these arguments to ``__init__`` is useful for storing Grids to the database, as - they are more stable (less likely to change) than the actual internal state of - the objects. - - The return value should be hashable, such that a set of these can be created. - """ - offset = None if not self._offset.any() else tuple(self._offset) - - bounds = _tuplify(self._bounds) - - # recreate a constructor-friendly version of `_unitSteps` from live data (may have been reduced from - # length 3 to length 2 or 1 based on mixing the step-based definition and the bounds-based definition - # described in Design Considerations above.) - # We don't just save the original tuple passed in because that may miss transformations that - # occur between instantiation and reduction. - unitSteps = [] - compressedSteps = list(self._unitSteps[:]) - for i in range(3): - # Recall _stepDims are stored as a single-value tuple (for numpy indexing) - # So this just is grabbing the actual data. - if i in self._stepDims[0]: - unitSteps.append(compressedSteps.pop(0)) - else: - # Add dummy value which will never get used (it gets reduced away) - unitSteps.append(0) - unitSteps = _tuplify(unitSteps) - - return GridParameters( - unitSteps, - bounds, - self._unitStepLimits, - offset, - self._geomType, - self._symmetry, - ) - - @property - def geomType(self) -> geometry.GeomType: - return geometry.GeomType.fromStr(self._geomType) - - @geomType.setter - def geomType(self, geomType: Union[str, geometry.GeomType]): - self._geomType = str(geometry.GeomType.fromAny(geomType)) - - @property - def symmetry(self) -> geometry.SymmetryType: - return geometry.SymmetryType.fromStr(self._symmetry) - - @symmetry.setter - def symmetry(self, symmetry: Union[str, geometry.SymmetryType]): - self._symmetry = str(geometry.SymmetryType.fromAny(symmetry)) - - @property - def offset(self): - return self._offset - - @offset.setter - def offset(self, offset): - self._offset = offset - - def __repr__(self): - msg = ( - ["<{} -- {}\nBounds:\n".format(self.__class__.__name__, id(self))] - + [" {}\n".format(b) for b in self._bounds] - + ["Steps:\n"] - + [" {}\n".format(b) for b in self._unitSteps] - + [ - "Anchor: {}\n".format(self.armiObject), - "Offset: {}\n".format(self._offset), - "Num Locations: {}>".format(len(self)), - ] - ) - return "".join(msg) - - def __getstate__(self): - """ - Pickling removes reference to ``armiObject``. - - Removing the ``armiObject`` allows us to pickle an assembly without pickling the entire - reactor. An ``Assembly.spatialLocator.grid.armiObject`` is the reactor, by removing the link - here, we still have spatial orientation, but are not required to pickle the entire reactor - to pickle an assembly. - - This relies on the ``armiObject.__setstate__`` to assign itself. - """ - state = self.__dict__.copy() - state["armiObject"] = None - - return state - - def __setstate__(self, state): - """ - Pickling removes reference to ``armiObject``. - - This relies on the ``ArmiObject.__setstate__`` to assign itself. - """ - self.__dict__.update(state) - - for _indices, locator in self.items(): - locator._grid = self - - def __getitem__(self, ijk: Union[Tuple[int, int, int], List[Tuple[int, int, int]]]): - """ - Get a location by (i, j, k) indices. If it does not exist, create a new one and return it. - - Parameters - ---------- - ijk : tuple of indices or list of the same - If provided a tuple, an IndexLocation will be created (if necessary) and - returned. If provided a list, each element will create a new IndexLocation - (if necessary), and a MultiIndexLocation containing all of the passed - indices will be returned. - - - Notes - ----- - The method is defaultdict-like, in that it will create a new location on the fly. However, - the class itself is not really a dictionary, it is just index-able. For example, there is no - desire to have a ``__setitem__`` method, because the only way to create a location is by - retrieving it or through ``buildLocations``. - """ - try: - return self._locations[ijk] - except (KeyError, TypeError): - pass - - if isinstance(ijk, tuple): - i, j, k = ijk - val = IndexLocation(i, j, k, self) - self._locations[ijk] = val - elif isinstance(ijk, list): - val = MultiIndexLocation(self) - locators = [self[idx] for idx in ijk] - val.extend(locators) - else: - raise TypeError( - "Unsupported index type `{}` for `{}`".format(type(ijk), ijk) - ) - return val - - def __len__(self): - return len(self._locations) - - def items(self): - """Return list of ((i, j, k), IndexLocation) tuples.""" - return self._locations.items() - - def backUp(self): - """Gather internal info that should be restored within a retainState.""" - self._backup = self._unitSteps, self._bounds, self._offset - - def restoreBackup(self): - self._unitSteps, self._bounds, self._offset = self._backup - - def getCoordinates(self, indices, nativeCoords=False) -> numpy.ndarray: - """Return the coordinates of the center of the mesh cell at the given given indices in cm.""" - indices = numpy.array(indices) - return self._evaluateMesh( - indices, self._centroidBySteps, self._centroidByBounds - ) - - def getCellBase(self, indices) -> numpy.ndarray: - """Get the mesh base (lower left) of this mesh cell in cm.""" - indices = numpy.array(indices) - return self._evaluateMesh( - indices, self._meshBaseBySteps, self._meshBaseByBounds - ) - - def getCellTop(self, indices) -> numpy.ndarray: - """Get the mesh top (upper right) of this mesh cell in cm.""" - indices = numpy.array(indices) + 1 - return self._evaluateMesh( - indices, self._meshBaseBySteps, self._meshBaseByBounds - ) - - def locatorInDomain( - self, locator: LocationBase, symmetryOverlap: Optional[bool] = False - ) -> bool: - """ - Return whether the passed locator is in the domain represented by the Grid. - - For instance, if we have a 1/3rd core hex grid, this would return False for - locators that are outside of the first third of the grid. - - Parameters - ---------- - locator : LocationBase - The location to test - symmetryOverlap : bool, optional - Whether grid locations along the symmetry line should be considered "in the - represented domain". This can be useful when assemblies are split along the - domain boundary, with fractions of the assembly on either side. - """ - raise NotImplementedError("Not implemented on base Grid type.") - - def _evaluateMesh(self, indices, stepOperator, boundsOperator) -> numpy.ndarray: - """ - Evaluate some function of indices on this grid. - - Recall from above that steps are mesh centered and bounds are mesh edged. - - Notes - ----- - This method may be able to be simplified. Complications from arbitrary - mixtures of bounds-based and step-based meshing caused it to get bad. - These were separate subclasses first, but in practice almost all cases have some mix - of step-based (hexagons, squares), and bounds based (radial, zeta). - - Improvements welcome! - """ - boundCoords = [] - for ii, bounds in enumerate(self._bounds): - if bounds is not None: - boundCoords.append(boundsOperator(indices[ii], bounds)) - - # limit step operator to the step dimensions - stepCoords = stepOperator(numpy.array(indices)[self._stepDims]) - - # now mix/match bounds coords with step coords appropriately. - result = numpy.zeros(len(indices)) - result[self._stepDims] = stepCoords - result[self._boundDims] = boundCoords - return result + self._offset - - def _centroidBySteps(self, indices): - return numpy.dot(self._unitSteps, indices) - - def _meshBaseBySteps(self, indices): - return ( - self._centroidBySteps(indices - 1) + self._centroidBySteps(indices) - ) / 2.0 - - @staticmethod - def _centroidByBounds(index, bounds): - if index < 0: - # avoid wrap-around - raise IndexError("Bounds-defined indices may not be negative.") - return (bounds[index + 1] + bounds[index]) / 2.0 - - @staticmethod - def _meshBaseByBounds(index, bounds): - if index < 0: - raise IndexError("Bounds-defined indices may not be negative.") - return bounds[index] - - @staticmethod - def getNeighboringCellIndices(i, j=0, k=0): - """Return the indices of the immediate neighbors of a mesh point in the plane.""" - return ((i + 1, j, k), (1, j + 1, k), (i - 1, j, k), (i, j - 1, k)) - - def getSymmetricEquivalents(self, indices): - """ - Return a list of grid indices that contain matching contents based on symmetry. - - The length of the list will depend on the type of symmetry being used, and - potentially the location of the requested indices. E.g., - third-core will return the two sets of indices at the matching location in the - other two thirds of the grid, unless it is the central location, in which case - no indices will be returned. - """ - raise NotImplementedError - - @staticmethod - def getAboveAndBelowCellIndices(indices): - i, j, k = indices - return ((i, j, k + 1), (i, j, k - 1)) - - def cellIndicesContainingPoint(self, x, y=0.0, z=0.0): - """Return the indices of a mesh cell that contains a point.""" - raise NotImplementedError - - def overlapsWhichSymmetryLine(self, indices): - return None - - def getLabel(self, indices): - """ - Get a string label from a 0-based spatial locator. - - Returns a string representing i, j, and k indices of the locator - """ - i, j = indices[:2] - label = f"{i:03d}-{j:03d}" - if len(indices) == 3: - label += f"-{indices[2]:03d}" - return label - - def getIndexBounds(self): - """ - Get min index and number of indices in this grid. - - Step-defined grids would be infinite but for the step limits defined in the constructor. - - Notes - ----- - This produces output that is intended to be passed to a ``range`` statement. - """ - indexBounds = [] - for minMax, bounds in zip(self._unitStepLimits, self._bounds): - if bounds is None: - indexBounds.append(minMax) - else: - indexBounds.append((0, len(bounds))) - return tuple(indexBounds) - - def getBounds( - self, - ) -> Tuple[ - Optional[Sequence[float]], Optional[Sequence[float]], Optional[Sequence[float]] - ]: - """Return the grid bounds for each dimension, if present.""" - return self._bounds - - def getLocatorFromRingAndPos(self, ring, pos, k=0): - """ - Return the location based on ring and position. - - Parameters - ---------- - ring : int - Ring number (1-based indexing) - pos : int - Position number (1-based indexing) - k : int, optional - Axial index (0-based indexing) - - See Also - -------- - getIndicesFromRingAndPos - This implements the transform into i, j indices based on ring and position. - """ - i, j = self.getIndicesFromRingAndPos(ring, pos) - return self[i, j, k] - - @staticmethod - def getIndicesFromRingAndPos(ring, pos): - """ - Return i, j indices given ring and position. - - Note - ---- - This should be implemented as a staticmethod, since no Grids currently in - exsistence actually need any instance data to perform this task, and - staticmethods provide the convenience of calling the method without an instance - of the class in the first place. - """ - raise NotImplementedError("Base Grid does not know about ring/pos") - - def getMinimumRings(self, n: int) -> int: - """ - Return the minimum number of rings needed to fit ``n`` objects. - - Warning - ------- - While this is useful and safe for answering the question of "how many rings do I - need to hold N things?", is generally not safe to use it to answer "I have N - things; within how many rings are they distributed?". This function provides a - lower bound, assuming that objects are densely-packed. If they are not actually - densely packed, this may be unphysical. - """ - raise NotImplementedError("Base grid does not know about rings") - - def getPositionsInRing(self, ring: int) -> int: - """Return the number of positions within a ring.""" - raise NotImplementedError("Base grid does not know about rings") - - def getRingPos(self, indices) -> Tuple[int, int]: - """ - Get ring and position number in this grid. - - For non-hex grids this is just i and j. - - A tuple is returned so that it is easy to compare pairs of indices. - """ - # Regular grids dont really know about ring and position. We can try to see if - # their parent does! - if ( - self.armiObject is not None - and self.armiObject.parent is not None - and self.armiObject.parent.spatialGrid is not None - ): - return self.armiObject.parent.spatialGrid.getRingPos(indices) - - # For compatibility's sake, return __something__. TODO: We may want to just - # throw here, to be honest. - return tuple(i + 1 for i in indices[:2]) - - def getAllIndices(self): - """Get all possible indices in this grid.""" - iBounds, jBounds, kBounds = self.getIndexBounds() - allIndices = tuple( - itertools.product(range(*iBounds), range(*jBounds), range(*kBounds)) - ) - return allIndices - - def buildLocations(self): - """Populate all grid cells with a characteristic SpatialLocator.""" - for i, j, k in self.getAllIndices(): - loc = IndexLocation(i, j, k, self) - self._locations[(i, j, k)] = loc - - @property - def pitch(self): - """ - The pitch of the grid. - - Assumes 2-d unit-step defined (works for cartesian) - """ - # TODO move this to the CartesianGrid - pitch = (self._unitSteps[0][0], self._unitSteps[1][1]) - if pitch[0] == 0: - raise ValueError(f"Grid {self} does not have a defined pitch.") - return pitch +from .grid import Grid, GridParameters, _tuplify class CartesianGrid(Grid): @@ -1477,13 +887,3 @@ def addingIsValid(myGrid, parentGrid): in pins to become inconsistent. """ return myGrid.isAxialOnly and not parentGrid.isAxialOnly - - -def _tuplify(maybeArray): - if isinstance(maybeArray, (numpy.ndarray, list, tuple)): - maybeArray = tuple( - tuple(row) if isinstance(row, (numpy.ndarray, list)) else row - for row in maybeArray - ) - - return maybeArray diff --git a/armi/reactor/grids/grid.py b/armi/reactor/grids/grid.py new file mode 100644 index 000000000..368ae644b --- /dev/null +++ b/armi/reactor/grids/grid.py @@ -0,0 +1,624 @@ +# Copyright 2023 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import itertools +import collections +from typing import Optional, List, Sequence, Union, Tuple + +import numpy + +from armi.reactor import geometry + +from .locations import IndexLocation, LocationBase, MultiIndexLocation + +# data structure for database-serialization of grids +GridParameters = collections.namedtuple( + "GridParameters", + ("unitSteps", "bounds", "unitStepLimits", "offset", "geomType", "symmetry"), +) + + +class Grid: + """ + A connected set of cells characterized by indices mapping to space and vice versa. + + The cells may be characterized by any mixture of regular repeating steps and + user-defined steps in any direction. + + For example, a 2-D hex lattice has constant, regular steps whereas a 3-D hex mesh + may have user-defined axial meshes. Similar for Cartesian, RZT, etc. + + Parameters + ---------- + unitSteps : tuple of tuples, optional + Describes the grid spatially as a function on indices. + Each tuple describes how each ``(x,y,or z)`` dimension is influenced by + ``(i,j,k)``. In other words, it is:: + + (dxi, dxj, jxk), (dyi, dyj, dyk), (dzi, dzj, dzk) + + where ``dmn`` is the distance (in cm) that dimension ``m`` will change as a + function of index ``n``. + + Unit steps are used as a generic method for defining repetitive grids in a + variety of geometries, including hexagonal and Cartesian. The tuples are not + vectors in the direction of the translation, but rather grouped by direction. If + the bounds argument is described for a direction, the bounds will be used rather + than the unit step information. The default of (0, 0, 0) makes all dimensions + insensitive to indices since the coordinates are calculated by the dot product + of this and the indices. With this default, any dimension that is desired to + change with indices should be defined with bounds. RZtheta grids are created + exclusively with bounds. + + bounds : 3-tuple + Absolute increasing bounds in cm including endpoints of a non-uniform grid. + Each item represents the boundaries in the associated direction. Use Nones when + unitSteps should be applied instead. Most useful for thetaRZ grids or other + non-uniform grids. + + unitStepLimits : 3-tuple + The limit of the steps in all three directions. This constrains step-defined + grids to be finite so we can populate them with SpatialLocator objects. + + offset : 3-tuple, optional + Offset in cm for each axis. By default the center of the (0,0,0)-th object is in + the center of the grid. Offsets can move it so that the (0,0,0)-th object can + be fully within a quadrant (i.e. in a Cartesian grid). + + armiObject : ArmiObject, optional + The ArmiObject that this grid describes. For example if it's a 1-D assembly + grid, the armiObject is the assembly. Note that ``self.armiObject.spatialGrid`` + is ``self``. + + Examples + -------- + A 2D a rectangular grid with width (x) 2 and height (y) 3 would be:: + + >>> grid = Grid(unitSteps=((2, 0, 0), (0, 3, 0),(0, 0, 0))) + + A regular hex grid with pitch 1 is:: + + >>> grid = Grid(unitSteps= ((sqrt(3)/2, 0.0, 0.0), (0.5, 1.0, 0.0), (0, 0, 0)) + + .. note:: For this unit hex the magnitude of the vector constructed using the + 0th index of each tuple is 1.0. + + Notes + ----- + Each dimension must either be defined through unitSteps or bounds. + The combination of unitSteps with bounds was settled upon after some struggle to + have one unified definition of a grid (i.e. just bounds). A hexagonal grid is + somewhat challenging to represent with bounds because the axes are not orthogonal, + so a unit-direction vector plus bounds would be required. And then the bounds would + be wasted space because they can be derived simply by unit steps. Memory efficiency + is important in this object so the compact representation of + unitSteps-when-possible, bounds-otherwise was settled upon. + + Design considerations include: + + * unitSteps are more intuitive as operations starting from the center of a cell, + particularly with hexagons and rectangles. Otherwise the 0,0 position of a hexagon + in the center of 1/3-symmetric hexagon is at the phantom bottom left of the + hexagon. + + * Users generally prefer to input mesh bounds rather than centers (e.g. starting at + 0.5 instead of 0.0 in a unit mesh is weird). + + * If we store bounds, computing bounds is simple and computing centers takes ~2x the + effort. If we store centers, it's the opposite. + + * Regardless of how we store things, we'll need a Grid that has the lower-left + assembly fully inside the problem (i.e. for full core Cartesian) as well as + another one that has the lower-left assembly half-way or quarter-way sliced off + (for 1/2, 1/4, and 1/8 symmetries). The ``offset`` parameter handles this. + + * Looking up mesh boundaries (to define a mesh in another code) is generally more + common than looking up centers (for plotting or measuring distance). + + * A grid can be anchored to the object that it is in with a backreference. This + gives it the ability to traverse the composite tree and map local to global + locations without having to duplicate the composite pattern on grids. This remains + optional so grids can be used for non-reactor-package reasons. It may seem + slightly cleaner to set the armiObject to the parent's spatialLocator itself + but the major disadvantage of this is that when an object moves, the armiObject + would have to be updated. By anchoring directly to Composite objects, the parent + is always up to date no matter where or how things get moved. + + * Unit step calculations use dot products and must not be polluted by the bound + indices. Thus we reduce the size of the unitSteps tuple accordingly. + + .. impl:: ARMI supports a number of structured mesh options. + :id: IMPL_REACTOR_MESH_0 + :links: REQ_REACTOR_MESH + """ + + def __init__( + self, + unitSteps=(0, 0, 0), + bounds=(None, None, None), + unitStepLimits=((0, 1), (0, 1), (0, 1)), + offset=None, + geomType="", + symmetry="", + armiObject=None, + ): + # these lists contain the indices representing which dimensions for which steps + # are used, or for which bounds are used. index 0 is x direction, etc. + self._boundDims = [] + self._stepDims = [] + for dimensionIndex, bound in enumerate(bounds): + if bound is None: + self._stepDims.append(dimensionIndex) + else: + self._boundDims.append(dimensionIndex) + + # numpy prefers tuples like this to do slicing on arrays + self._boundDims = (tuple(self._boundDims),) + self._stepDims = (tuple(self._stepDims),) + + unitSteps = _tuplify(unitSteps) + + self._bounds = bounds + self._unitStepLimits = _tuplify(unitStepLimits) + + # only represent unit steps in dimensions they're being used so as to not + # pollute the dot product. This may reduce the length of this from 3 to 2 or 1 + self._unitSteps = numpy.array(unitSteps)[self._stepDims] + self._offset = numpy.zeros(3) if offset is None else numpy.array(offset) + + self._locations = {} # indices -> SpatialLocator map + self.armiObject = armiObject + self.buildLocations() # locations are owned by a grid, so the grid builds them. + self._backup = None # for retainState + + (_ii, iLen), (_ji, jLen), (_ki, kLen) = self.getIndexBounds() + # True if only contains k-cells. + self.isAxialOnly = iLen == jLen == 1 and kLen > 1 + + # geometric metadata encapsulated here because it's related to the grid. + # They do not impact the grid object itself. + # Notice that these are stored using their string representations, rather than + # the GridType enum. This avoids the danger of deserializing an enum value from + # an old version of the code that may have had different numeric values. + self._geomType: str = str(geomType) + self._symmetry: str = str(symmetry) + + def reduce(self): + """ + Return the set of arguments used to create this Grid. + + This is very much like the argument tuple from ``__reduce__``, but we do not + implement ``__reduce__`` for real, because we are generally happy with + ``__getstate__`` and ``__setstate__`` for pickling purposes. However, getting + these arguments to ``__init__`` is useful for storing Grids to the database, as + they are more stable (less likely to change) than the actual internal state of + the objects. + + The return value should be hashable, such that a set of these can be created. + """ + offset = None if not self._offset.any() else tuple(self._offset) + + bounds = _tuplify(self._bounds) + + # recreate a constructor-friendly version of `_unitSteps` from live data (may have been reduced from + # length 3 to length 2 or 1 based on mixing the step-based definition and the bounds-based definition + # described in Design Considerations above.) + # We don't just save the original tuple passed in because that may miss transformations that + # occur between instantiation and reduction. + unitSteps = [] + compressedSteps = list(self._unitSteps[:]) + for i in range(3): + # Recall _stepDims are stored as a single-value tuple (for numpy indexing) + # So this just is grabbing the actual data. + if i in self._stepDims[0]: + unitSteps.append(compressedSteps.pop(0)) + else: + # Add dummy value which will never get used (it gets reduced away) + unitSteps.append(0) + unitSteps = _tuplify(unitSteps) + + return GridParameters( + unitSteps, + bounds, + self._unitStepLimits, + offset, + self._geomType, + self._symmetry, + ) + + @property + def geomType(self) -> geometry.GeomType: + return geometry.GeomType.fromStr(self._geomType) + + @geomType.setter + def geomType(self, geomType: Union[str, geometry.GeomType]): + self._geomType = str(geometry.GeomType.fromAny(geomType)) + + @property + def symmetry(self) -> geometry.SymmetryType: + return geometry.SymmetryType.fromStr(self._symmetry) + + @symmetry.setter + def symmetry(self, symmetry: Union[str, geometry.SymmetryType]): + self._symmetry = str(geometry.SymmetryType.fromAny(symmetry)) + + @property + def offset(self): + return self._offset + + @offset.setter + def offset(self, offset): + self._offset = offset + + def __repr__(self): + msg = ( + ["<{} -- {}\nBounds:\n".format(self.__class__.__name__, id(self))] + + [" {}\n".format(b) for b in self._bounds] + + ["Steps:\n"] + + [" {}\n".format(b) for b in self._unitSteps] + + [ + "Anchor: {}\n".format(self.armiObject), + "Offset: {}\n".format(self._offset), + "Num Locations: {}>".format(len(self)), + ] + ) + return "".join(msg) + + def __getstate__(self): + """ + Pickling removes reference to ``armiObject``. + + Removing the ``armiObject`` allows us to pickle an assembly without pickling the entire + reactor. An ``Assembly.spatialLocator.grid.armiObject`` is the reactor, by removing the link + here, we still have spatial orientation, but are not required to pickle the entire reactor + to pickle an assembly. + + This relies on the ``armiObject.__setstate__`` to assign itself. + """ + state = self.__dict__.copy() + state["armiObject"] = None + + return state + + def __setstate__(self, state): + """ + Pickling removes reference to ``armiObject``. + + This relies on the ``ArmiObject.__setstate__`` to assign itself. + """ + self.__dict__.update(state) + + for _indices, locator in self.items(): + locator._grid = self + + def __getitem__(self, ijk: Union[Tuple[int, int, int], List[Tuple[int, int, int]]]): + """ + Get a location by (i, j, k) indices. If it does not exist, create a new one and return it. + + Parameters + ---------- + ijk : tuple of indices or list of the same + If provided a tuple, an IndexLocation will be created (if necessary) and + returned. If provided a list, each element will create a new IndexLocation + (if necessary), and a MultiIndexLocation containing all of the passed + indices will be returned. + + + Notes + ----- + The method is defaultdict-like, in that it will create a new location on the fly. However, + the class itself is not really a dictionary, it is just index-able. For example, there is no + desire to have a ``__setitem__`` method, because the only way to create a location is by + retrieving it or through ``buildLocations``. + """ + try: + return self._locations[ijk] + except (KeyError, TypeError): + pass + + if isinstance(ijk, tuple): + i, j, k = ijk + val = IndexLocation(i, j, k, self) + self._locations[ijk] = val + elif isinstance(ijk, list): + val = MultiIndexLocation(self) + locators = [self[idx] for idx in ijk] + val.extend(locators) + else: + raise TypeError( + "Unsupported index type `{}` for `{}`".format(type(ijk), ijk) + ) + return val + + def __len__(self): + return len(self._locations) + + def items(self): + """Return list of ((i, j, k), IndexLocation) tuples.""" + return self._locations.items() + + def backUp(self): + """Gather internal info that should be restored within a retainState.""" + self._backup = self._unitSteps, self._bounds, self._offset + + def restoreBackup(self): + self._unitSteps, self._bounds, self._offset = self._backup + + def getCoordinates(self, indices, nativeCoords=False) -> numpy.ndarray: + """Return the coordinates of the center of the mesh cell at the given given indices in cm.""" + indices = numpy.array(indices) + return self._evaluateMesh( + indices, self._centroidBySteps, self._centroidByBounds + ) + + def getCellBase(self, indices) -> numpy.ndarray: + """Get the mesh base (lower left) of this mesh cell in cm.""" + indices = numpy.array(indices) + return self._evaluateMesh( + indices, self._meshBaseBySteps, self._meshBaseByBounds + ) + + def getCellTop(self, indices) -> numpy.ndarray: + """Get the mesh top (upper right) of this mesh cell in cm.""" + indices = numpy.array(indices) + 1 + return self._evaluateMesh( + indices, self._meshBaseBySteps, self._meshBaseByBounds + ) + + def locatorInDomain( + self, locator: LocationBase, symmetryOverlap: Optional[bool] = False + ) -> bool: + """ + Return whether the passed locator is in the domain represented by the Grid. + + For instance, if we have a 1/3rd core hex grid, this would return False for + locators that are outside of the first third of the grid. + + Parameters + ---------- + locator : LocationBase + The location to test + symmetryOverlap : bool, optional + Whether grid locations along the symmetry line should be considered "in the + represented domain". This can be useful when assemblies are split along the + domain boundary, with fractions of the assembly on either side. + """ + raise NotImplementedError("Not implemented on base Grid type.") + + def _evaluateMesh(self, indices, stepOperator, boundsOperator) -> numpy.ndarray: + """ + Evaluate some function of indices on this grid. + + Recall from above that steps are mesh centered and bounds are mesh edged. + + Notes + ----- + This method may be able to be simplified. Complications from arbitrary + mixtures of bounds-based and step-based meshing caused it to get bad. + These were separate subclasses first, but in practice almost all cases have some mix + of step-based (hexagons, squares), and bounds based (radial, zeta). + + Improvements welcome! + """ + boundCoords = [] + for ii, bounds in enumerate(self._bounds): + if bounds is not None: + boundCoords.append(boundsOperator(indices[ii], bounds)) + + # limit step operator to the step dimensions + stepCoords = stepOperator(numpy.array(indices)[self._stepDims]) + + # now mix/match bounds coords with step coords appropriately. + result = numpy.zeros(len(indices)) + result[self._stepDims] = stepCoords + result[self._boundDims] = boundCoords + return result + self._offset + + def _centroidBySteps(self, indices): + return numpy.dot(self._unitSteps, indices) + + def _meshBaseBySteps(self, indices): + return ( + self._centroidBySteps(indices - 1) + self._centroidBySteps(indices) + ) / 2.0 + + @staticmethod + def _centroidByBounds(index, bounds): + if index < 0: + # avoid wrap-around + raise IndexError("Bounds-defined indices may not be negative.") + return (bounds[index + 1] + bounds[index]) / 2.0 + + @staticmethod + def _meshBaseByBounds(index, bounds): + if index < 0: + raise IndexError("Bounds-defined indices may not be negative.") + return bounds[index] + + @staticmethod + def getNeighboringCellIndices(i, j=0, k=0): + """Return the indices of the immediate neighbors of a mesh point in the plane.""" + return ((i + 1, j, k), (1, j + 1, k), (i - 1, j, k), (i, j - 1, k)) + + def getSymmetricEquivalents(self, indices): + """ + Return a list of grid indices that contain matching contents based on symmetry. + + The length of the list will depend on the type of symmetry being used, and + potentially the location of the requested indices. E.g., + third-core will return the two sets of indices at the matching location in the + other two thirds of the grid, unless it is the central location, in which case + no indices will be returned. + """ + raise NotImplementedError + + @staticmethod + def getAboveAndBelowCellIndices(indices): + i, j, k = indices + return ((i, j, k + 1), (i, j, k - 1)) + + def cellIndicesContainingPoint(self, x, y=0.0, z=0.0): + """Return the indices of a mesh cell that contains a point.""" + raise NotImplementedError + + def overlapsWhichSymmetryLine(self, indices): + return None + + def getLabel(self, indices): + """ + Get a string label from a 0-based spatial locator. + + Returns a string representing i, j, and k indices of the locator + """ + i, j = indices[:2] + label = f"{i:03d}-{j:03d}" + if len(indices) == 3: + label += f"-{indices[2]:03d}" + return label + + def getIndexBounds(self): + """ + Get min index and number of indices in this grid. + + Step-defined grids would be infinite but for the step limits defined in the constructor. + + Notes + ----- + This produces output that is intended to be passed to a ``range`` statement. + """ + indexBounds = [] + for minMax, bounds in zip(self._unitStepLimits, self._bounds): + if bounds is None: + indexBounds.append(minMax) + else: + indexBounds.append((0, len(bounds))) + return tuple(indexBounds) + + def getBounds( + self, + ) -> Tuple[ + Optional[Sequence[float]], Optional[Sequence[float]], Optional[Sequence[float]] + ]: + """Return the grid bounds for each dimension, if present.""" + return self._bounds + + def getLocatorFromRingAndPos(self, ring, pos, k=0): + """ + Return the location based on ring and position. + + Parameters + ---------- + ring : int + Ring number (1-based indexing) + pos : int + Position number (1-based indexing) + k : int, optional + Axial index (0-based indexing) + + See Also + -------- + getIndicesFromRingAndPos + This implements the transform into i, j indices based on ring and position. + """ + i, j = self.getIndicesFromRingAndPos(ring, pos) + return self[i, j, k] + + @staticmethod + def getIndicesFromRingAndPos(ring, pos): + """ + Return i, j indices given ring and position. + + Note + ---- + This should be implemented as a staticmethod, since no Grids currently in + exsistence actually need any instance data to perform this task, and + staticmethods provide the convenience of calling the method without an instance + of the class in the first place. + """ + raise NotImplementedError("Base Grid does not know about ring/pos") + + def getMinimumRings(self, n: int) -> int: + """ + Return the minimum number of rings needed to fit ``n`` objects. + + Warning + ------- + While this is useful and safe for answering the question of "how many rings do I + need to hold N things?", is generally not safe to use it to answer "I have N + things; within how many rings are they distributed?". This function provides a + lower bound, assuming that objects are densely-packed. If they are not actually + densely packed, this may be unphysical. + """ + raise NotImplementedError("Base grid does not know about rings") + + def getPositionsInRing(self, ring: int) -> int: + """Return the number of positions within a ring.""" + raise NotImplementedError("Base grid does not know about rings") + + def getRingPos(self, indices) -> Tuple[int, int]: + """ + Get ring and position number in this grid. + + For non-hex grids this is just i and j. + + A tuple is returned so that it is easy to compare pairs of indices. + """ + # Regular grids dont really know about ring and position. We can try to see if + # their parent does! + if ( + self.armiObject is not None + and self.armiObject.parent is not None + and self.armiObject.parent.spatialGrid is not None + ): + return self.armiObject.parent.spatialGrid.getRingPos(indices) + + # For compatibility's sake, return __something__. TODO: We may want to just + # throw here, to be honest. + return tuple(i + 1 for i in indices[:2]) + + def getAllIndices(self): + """Get all possible indices in this grid.""" + iBounds, jBounds, kBounds = self.getIndexBounds() + allIndices = tuple( + itertools.product(range(*iBounds), range(*jBounds), range(*kBounds)) + ) + return allIndices + + def buildLocations(self): + """Populate all grid cells with a characteristic SpatialLocator.""" + for i, j, k in self.getAllIndices(): + loc = IndexLocation(i, j, k, self) + self._locations[(i, j, k)] = loc + + @property + def pitch(self): + """ + The pitch of the grid. + + Assumes 2-d unit-step defined (works for cartesian) + """ + # TODO move this to the CartesianGrid + pitch = (self._unitSteps[0][0], self._unitSteps[1][1]) + if pitch[0] == 0: + raise ValueError(f"Grid {self} does not have a defined pitch.") + return pitch + + +def _tuplify(maybeArray): + if isinstance(maybeArray, (numpy.ndarray, list, tuple)): + maybeArray = tuple( + tuple(row) if isinstance(row, (numpy.ndarray, list)) else row + for row in maybeArray + ) + + return maybeArray From 5e440646482d1195dcdf4689e54fb0aea8a5da65 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 13:11:33 -0700 Subject: [PATCH 06/16] Add armi.reactor.grids.cartesian for CartesianGrid --- armi/reactor/grids/__init__.py | 261 +---------------------------- armi/reactor/grids/cartesian.py | 280 ++++++++++++++++++++++++++++++++ 2 files changed, 281 insertions(+), 260 deletions(-) create mode 100644 armi/reactor/grids/cartesian.py diff --git a/armi/reactor/grids/__init__.py b/armi/reactor/grids/__init__.py index 1766ad1b7..358bb9363 100644 --- a/armi/reactor/grids/__init__.py +++ b/armi/reactor/grids/__init__.py @@ -60,7 +60,6 @@ while the word **local** refers to within the current coordinate system defined by the current grid. """ -import itertools import math from typing import Tuple, Optional @@ -102,265 +101,7 @@ ) from .grid import Grid, GridParameters, _tuplify - - -class CartesianGrid(Grid): - """ - Grid class representing a conformal Cartesian mesh. - - Notes - ----- - In Cartesian, (i, j, k) indices map to (x, y, z) coordinates. - In an axial plane (i, j) are as follows:: - - (-1, 1) ( 0, 1) ( 1, 1) - (-1, 0) ( 0, 0) ( 1, 0) - (-1,-1) ( 0,-1) ( 1,-1) - - The concepts of ring and position are a bit tricker in Cartesian grids than in Hex, - because unlike in the Hex case, there is no guaranteed center location. For example, - when using a CartesianGrid to lay out assemblies in a core, there is only a single - central location if the number of assemblies in the core is odd-by-odd; in an - even-by-even case, there are four center-most assemblies. Therefore, the number of - locations per ring will vary depending on the "through center" nature of - ``symmetry``. - - Furthermore, notice that in the "through center" (odd-by-odd) case, the central - index location, (0,0) is typically centered at the origin (0.0, 0.0), whereas with - the "not through center" (even-by-even) case, the (0,0) index location is offset, - away from the origin. - - These concepts are illustrated in the example drawings below. - - .. figure:: ../.static/through-center.png - :width: 400px - :align: center - - Grid example where the axes pass through the "center assembly" (odd-by-odd). - Note that ring 1 only has one location in it. - - .. figure:: ../.static/not-through-center.png - :width: 400px - :align: center - - Grid example where the axes lie between the "center assemblies" (even-by-even). - Note that ring 1 has four locations, and that the center of the (0, 0)-index - location is offset from the origin. - - .. impl:: ARMI supports a Cartesian mesh. - :id: IMPL_REACTOR_MESH_1 - :links: REQ_REACTOR_MESH - """ - - @classmethod - def fromRectangle( - cls, width, height, numRings=5, symmetry="", isOffset=False, armiObject=None - ): - """ - Build a finite step-based 2-D Cartesian grid based on a width and height in cm. - - Parameters - ---------- - width : float - Width of the unit rectangle - height : float - Height of the unit rectangle - numRings : int - Number of rings that the grid should span - symmetry : str - The symmetry condition (see :py:mod:`armi.reactor.geometry`) - isOffset : bool - If True, the origin of the Grid's coordinate system will be placed at the - bottom-left corner of the center-most cell. Otherwise, the origin will be - placed at the center of the center-most cell. - armiObject : ArmiObject - An object in a Composite model that the Grid should be bound to. - """ - unitSteps = ((width, 0.0, 0.0), (0.0, height, 0.0), (0, 0, 0)) - offset = numpy.array((width / 2.0, height / 2.0, 0.0)) if isOffset else None - return cls( - unitSteps=unitSteps, - unitStepLimits=((-numRings, numRings), (-numRings, numRings), (0, 1)), - offset=offset, - armiObject=armiObject, - symmetry=symmetry, - ) - - def getRingPos(self, indices): - """ - Return ring and position from indices. - - Ring is the Manhattan distance from (0, 0) to the passed indices. Position - counts up around the ring counter-clockwise from the quadrant 1 diagonal, like - this:: - - 7 6 5 4 3 2 1 - 8 | 24 - 9 | 23 - 10 -------|------ 22 - 11 | 21 - 12 | 20 - 13 14 15 16 17 18 19 - - Grids that split the central locations have 1 location in in inner-most ring, - whereas grids without split central locations will have 4. - - Notes - ----- - This is needed to support GUI, but should not often be used. - i, j (0-based) indices are much more useful. For example: - - >>> locator = core.spatialGrid[i, j, 0] # 3rd index is 0 for assembly - >>> a = core.childrenByLocator[locator] - - >>> a = core.childrenByLocator[core.spatialGrid[i, j, 0]] # one liner - """ - i, j = indices[0:2] - split = self._isThroughCenter() - - if not split: - i += 0.5 - j += 0.5 - - ring = max(abs(int(i)), abs(int(j))) - - if not split: - ring += 0.5 - - if j == ring: - # region 1 - pos = -i + ring - elif i == -ring: - # region 2 - pos = 3 * ring - j - elif j == -ring: - # region 3 - pos = 5 * ring + i - else: - # region 4 - pos = 7 * ring + j - return (int(ring) + 1, int(pos) + 1) - - @staticmethod - def getIndicesFromRingAndPos(ring, pos): - """Not implemented for Cartesian-see getRingPos notes.""" - raise NotImplementedError( - "Cartesian should not need need ring/pos, use i, j indices." - "See getRingPos doc string notes for more information/example." - ) - - def getMinimumRings(self, n): - """Return the minimum number of rings needed to fit ``n`` objects.""" - numPositions = 0 - ring = 0 - for ring in itertools.count(1): - ringPositions = self.getPositionsInRing(ring) - numPositions += ringPositions - if numPositions >= n: - break - - return ring - - def getPositionsInRing(self, ring): - """ - Return the number of positions within a ring. - - Notes - ----- - The number of positions within a ring will change - depending on whether the central position in the - grid is at origin, or if origin is the point - where 4 positions meet (i.e., the ``_isThroughCenter`` - method returns True). - """ - if ring == 1: - ringPositions = 1 if self._isThroughCenter() else 4 - else: - ringPositions = (ring - 1) * 8 - if not self._isThroughCenter(): - ringPositions += 4 - return ringPositions - - def locatorInDomain(self, locator, symmetryOverlap: Optional[bool] = False): - if self.symmetry.domain == geometry.DomainType.QUARTER_CORE: - return locator.i >= 0 and locator.j >= 0 - else: - return True - - def changePitch(self, xw, yw): - """ - Change the pitch of a Cartesian grid. - - This also scales the offset. - """ - xwOld = self._unitSteps[0][0] - ywOld = self._unitSteps[1][1] - self._unitSteps = numpy.array(((xw, 0.0, 0.0), (0.0, yw, 0.0), (0, 0, 0)))[ - self._stepDims - ] - newOffsetX = self._offset[0] * xw / xwOld - newOffsetY = self._offset[1] * yw / ywOld - self._offset = numpy.array((newOffsetX, newOffsetY, 0.0)) - - def getSymmetricEquivalents(self, indices): - symmetry = self.symmetry # construct the symmetry object once up top - isRotational = symmetry.boundary == geometry.BoundaryType.PERIODIC - - i, j = indices[0:2] - if symmetry.domain == geometry.DomainType.FULL_CORE: - return [] - elif symmetry.domain == geometry.DomainType.QUARTER_CORE: - if symmetry.isThroughCenterAssembly: - # some locations lie on the symmetric boundary - if i == 0 and j == 0: - # on the split corner, so the location is its own symmetric - # equivalent - return [] - elif i == 0: - if isRotational: - return [(j, i), (i, -j), (-j, i)] - else: - return [(i, -j)] - elif j == 0: - if isRotational: - return [(j, i), (-i, j), (j, -i)] - else: - return [(-i, j)] - else: - # Math is a bit easier for the split case, since there is an actual - # center location for (0, 0) - if isRotational: - return [(-j, i), (-i, -j), (j, -i)] - else: - return [(-i, j), (-i, -j), (i, -j)] - else: - # most objects have 3 equivalents. the bottom-left corner of Quadrant I - # is (0, 0), so to reflect, add one and negate each index in - # combination. To rotate, first flip the indices for the Quadrant II and - # Quadrant IV - if isRotational: - # rotational - # QII QIII QIV - return [(-j - 1, i), (-i - 1, -j - 1), (j, -i - 1)] - else: - # reflective - # QII QIII QIV - return [(-i - 1, j), (-i - 1, -j - 1), (i, -j - 1)] - - elif symmetry.domain == geometry.DomainType.EIGHTH_CORE: - raise NotImplementedError( - "Eighth-core symmetry isn't fully implemented for grids yet!" - ) - else: - raise NotImplementedError( - "Unhandled symmetry condition for {}: {}".format( - type(self).__name__, symmetry.domain - ) - ) - - def _isThroughCenter(self): - """Return whether the central cells are split through the middle for symmetry.""" - return all(self._offset == [0, 0, 0]) +from .cartesian import CartesianGrid class HexGrid(Grid): diff --git a/armi/reactor/grids/cartesian.py b/armi/reactor/grids/cartesian.py new file mode 100644 index 000000000..99a6ded09 --- /dev/null +++ b/armi/reactor/grids/cartesian.py @@ -0,0 +1,280 @@ +# Copyright 2023 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import itertools +from typing import Optional + +import numpy + +from armi.reactor import geometry + +from .grid import Grid + + +class CartesianGrid(Grid): + """ + Grid class representing a conformal Cartesian mesh. + + Notes + ----- + In Cartesian, (i, j, k) indices map to (x, y, z) coordinates. + In an axial plane (i, j) are as follows:: + + (-1, 1) ( 0, 1) ( 1, 1) + (-1, 0) ( 0, 0) ( 1, 0) + (-1,-1) ( 0,-1) ( 1,-1) + + The concepts of ring and position are a bit tricker in Cartesian grids than in Hex, + because unlike in the Hex case, there is no guaranteed center location. For example, + when using a CartesianGrid to lay out assemblies in a core, there is only a single + central location if the number of assemblies in the core is odd-by-odd; in an + even-by-even case, there are four center-most assemblies. Therefore, the number of + locations per ring will vary depending on the "through center" nature of + ``symmetry``. + + Furthermore, notice that in the "through center" (odd-by-odd) case, the central + index location, (0,0) is typically centered at the origin (0.0, 0.0), whereas with + the "not through center" (even-by-even) case, the (0,0) index location is offset, + away from the origin. + + These concepts are illustrated in the example drawings below. + + .. figure:: ../.static/through-center.png + :width: 400px + :align: center + + Grid example where the axes pass through the "center assembly" (odd-by-odd). + Note that ring 1 only has one location in it. + + .. figure:: ../.static/not-through-center.png + :width: 400px + :align: center + + Grid example where the axes lie between the "center assemblies" (even-by-even). + Note that ring 1 has four locations, and that the center of the (0, 0)-index + location is offset from the origin. + + .. impl:: ARMI supports a Cartesian mesh. + :id: IMPL_REACTOR_MESH_1 + :links: REQ_REACTOR_MESH + """ + + @classmethod + def fromRectangle( + cls, width, height, numRings=5, symmetry="", isOffset=False, armiObject=None + ): + """ + Build a finite step-based 2-D Cartesian grid based on a width and height in cm. + + Parameters + ---------- + width : float + Width of the unit rectangle + height : float + Height of the unit rectangle + numRings : int + Number of rings that the grid should span + symmetry : str + The symmetry condition (see :py:mod:`armi.reactor.geometry`) + isOffset : bool + If True, the origin of the Grid's coordinate system will be placed at the + bottom-left corner of the center-most cell. Otherwise, the origin will be + placed at the center of the center-most cell. + armiObject : ArmiObject + An object in a Composite model that the Grid should be bound to. + """ + unitSteps = ((width, 0.0, 0.0), (0.0, height, 0.0), (0, 0, 0)) + offset = numpy.array((width / 2.0, height / 2.0, 0.0)) if isOffset else None + return cls( + unitSteps=unitSteps, + unitStepLimits=((-numRings, numRings), (-numRings, numRings), (0, 1)), + offset=offset, + armiObject=armiObject, + symmetry=symmetry, + ) + + def getRingPos(self, indices): + """ + Return ring and position from indices. + + Ring is the Manhattan distance from (0, 0) to the passed indices. Position + counts up around the ring counter-clockwise from the quadrant 1 diagonal, like + this:: + + 7 6 5 4 3 2 1 + 8 | 24 + 9 | 23 + 10 -------|------ 22 + 11 | 21 + 12 | 20 + 13 14 15 16 17 18 19 + + Grids that split the central locations have 1 location in in inner-most ring, + whereas grids without split central locations will have 4. + + Notes + ----- + This is needed to support GUI, but should not often be used. + i, j (0-based) indices are much more useful. For example: + + >>> locator = core.spatialGrid[i, j, 0] # 3rd index is 0 for assembly + >>> a = core.childrenByLocator[locator] + + >>> a = core.childrenByLocator[core.spatialGrid[i, j, 0]] # one liner + """ + i, j = indices[0:2] + split = self._isThroughCenter() + + if not split: + i += 0.5 + j += 0.5 + + ring = max(abs(int(i)), abs(int(j))) + + if not split: + ring += 0.5 + + if j == ring: + # region 1 + pos = -i + ring + elif i == -ring: + # region 2 + pos = 3 * ring - j + elif j == -ring: + # region 3 + pos = 5 * ring + i + else: + # region 4 + pos = 7 * ring + j + return (int(ring) + 1, int(pos) + 1) + + @staticmethod + def getIndicesFromRingAndPos(ring, pos): + """Not implemented for Cartesian-see getRingPos notes.""" + raise NotImplementedError( + "Cartesian should not need need ring/pos, use i, j indices." + "See getRingPos doc string notes for more information/example." + ) + + def getMinimumRings(self, n): + """Return the minimum number of rings needed to fit ``n`` objects.""" + numPositions = 0 + ring = 0 + for ring in itertools.count(1): + ringPositions = self.getPositionsInRing(ring) + numPositions += ringPositions + if numPositions >= n: + break + + return ring + + def getPositionsInRing(self, ring): + """ + Return the number of positions within a ring. + + Notes + ----- + The number of positions within a ring will change + depending on whether the central position in the + grid is at origin, or if origin is the point + where 4 positions meet (i.e., the ``_isThroughCenter`` + method returns True). + """ + if ring == 1: + ringPositions = 1 if self._isThroughCenter() else 4 + else: + ringPositions = (ring - 1) * 8 + if not self._isThroughCenter(): + ringPositions += 4 + return ringPositions + + def locatorInDomain(self, locator, symmetryOverlap: Optional[bool] = False): + if self.symmetry.domain == geometry.DomainType.QUARTER_CORE: + return locator.i >= 0 and locator.j >= 0 + else: + return True + + def changePitch(self, xw, yw): + """ + Change the pitch of a Cartesian grid. + + This also scales the offset. + """ + xwOld = self._unitSteps[0][0] + ywOld = self._unitSteps[1][1] + self._unitSteps = numpy.array(((xw, 0.0, 0.0), (0.0, yw, 0.0), (0, 0, 0)))[ + self._stepDims + ] + newOffsetX = self._offset[0] * xw / xwOld + newOffsetY = self._offset[1] * yw / ywOld + self._offset = numpy.array((newOffsetX, newOffsetY, 0.0)) + + def getSymmetricEquivalents(self, indices): + symmetry = self.symmetry # construct the symmetry object once up top + isRotational = symmetry.boundary == geometry.BoundaryType.PERIODIC + + i, j = indices[0:2] + if symmetry.domain == geometry.DomainType.FULL_CORE: + return [] + elif symmetry.domain == geometry.DomainType.QUARTER_CORE: + if symmetry.isThroughCenterAssembly: + # some locations lie on the symmetric boundary + if i == 0 and j == 0: + # on the split corner, so the location is its own symmetric + # equivalent + return [] + elif i == 0: + if isRotational: + return [(j, i), (i, -j), (-j, i)] + else: + return [(i, -j)] + elif j == 0: + if isRotational: + return [(j, i), (-i, j), (j, -i)] + else: + return [(-i, j)] + else: + # Math is a bit easier for the split case, since there is an actual + # center location for (0, 0) + if isRotational: + return [(-j, i), (-i, -j), (j, -i)] + else: + return [(-i, j), (-i, -j), (i, -j)] + else: + # most objects have 3 equivalents. the bottom-left corner of Quadrant I + # is (0, 0), so to reflect, add one and negate each index in + # combination. To rotate, first flip the indices for the Quadrant II and + # Quadrant IV + if isRotational: + # rotational + # QII QIII QIV + return [(-j - 1, i), (-i - 1, -j - 1), (j, -i - 1)] + else: + # reflective + # QII QIII QIV + return [(-i - 1, j), (-i - 1, -j - 1), (i, -j - 1)] + + elif symmetry.domain == geometry.DomainType.EIGHTH_CORE: + raise NotImplementedError( + "Eighth-core symmetry isn't fully implemented for grids yet!" + ) + else: + raise NotImplementedError( + "Unhandled symmetry condition for {}: {}".format( + type(self).__name__, symmetry.domain + ) + ) + + def _isThroughCenter(self): + """Return whether the central cells are split through the middle for symmetry.""" + return all(self._offset == [0, 0, 0]) From b06a78e9988feaab10043fdce46ae6002ced1214 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 13:17:47 -0700 Subject: [PATCH 07/16] Add armi.reactor.grids.hexagonal --- armi/reactor/grids/__init__.py | 404 +----------------------------- armi/reactor/grids/hexagonal.py | 426 ++++++++++++++++++++++++++++++++ 2 files changed, 427 insertions(+), 403 deletions(-) create mode 100644 armi/reactor/grids/hexagonal.py diff --git a/armi/reactor/grids/__init__.py b/armi/reactor/grids/__init__.py index 358bb9363..eff3bd595 100644 --- a/armi/reactor/grids/__init__.py +++ b/armi/reactor/grids/__init__.py @@ -63,29 +63,9 @@ import math from typing import Tuple, Optional -import numpy.linalg - -from armi.utils import hexagon -from armi.reactor import geometry - - TAU = math.pi * 2.0 -COS30 = math.sqrt(3) / 2.0 -SIN30 = 1.0 / 2.0 -# going CCW from "position 1" (top right) -TRIANGLES_IN_HEXAGON = numpy.array( - [ - (+COS30, SIN30), - (+0, 1.0), - (-COS30, SIN30), - (-COS30, -SIN30), - (+0, -1.0), - (+COS30, -SIN30), - ] -) - from .constants import ( BOUNDARY_CENTER, BOUNDARY_0_DEGREES, @@ -102,389 +82,7 @@ from .grid import Grid, GridParameters, _tuplify from .cartesian import CartesianGrid - - -class HexGrid(Grid): - """ - Has 6 neighbors in plane. - - Notes - ----- - In an axial plane (i, j) are as follows (second one is pointedEndUp):: - - - ( 0, 1) - (-1, 1) ( 1, 0) - ( 0, 0) - (-1, 0) ( 1,-1) - ( 0,-1) - - - ( 0, 1) ( 1, 0) - - (-1, 1) ( 0, 0) ( 1,-1) - - (-1, 0) ( 0,-1) - - .. impl:: ARMI supports a Hexagonal mesh. - :id: IMPL_REACTOR_MESH_2 - :links: REQ_REACTOR_MESH - """ - - @staticmethod - def fromPitch(pitch, numRings=25, armiObject=None, pointedEndUp=False, symmetry=""): - """ - Build a finite step-based 2-D hex grid from a hex pitch in cm. - - Parameters - ---------- - pitch : float - Hex pitch (flat-to-flat) in cm - numRings : int, optional - The number of rings in the grid to pre-populate with locatator objects. - Even if positions are not pre-populated, locators will be generated - there on the fly. - armiObject : ArmiObject, optional - The object that this grid is anchored to (i.e. the reactor for a grid of - assemblies) - pointedEndUp : bool, optional - Rotate the hexagons 30 degrees so that the pointed end faces up instead of - the flat. - symmetry : string, optional - A string representation of the symmetry options for the grid. - - Returns - ------- - HexGrid - A functional hexagonal grid object. - """ - side = pitch / math.sqrt(3.0) - if pointedEndUp: - # rotated 30 degrees CCW from normal - # increases in i move you in x and y - # increases in j also move you in x and y - unitSteps = ( - (pitch / 2.0, -pitch / 2.0, 0), - (1.5 * side, 1.5 * side, 0), - (0, 0, 0), - ) - else: - # x direction is only a function of i because j-axis is vertical. - # y direction is a function of both. - unitSteps = ((1.5 * side, 0.0, 0.0), (pitch / 2.0, pitch, 0.0), (0, 0, 0)) - - return HexGrid( - unitSteps=unitSteps, - unitStepLimits=((-numRings, numRings), (-numRings, numRings), (0, 1)), - armiObject=armiObject, - symmetry=symmetry, - ) - - @property - def pitch(self): - """ - Get the hex-pitch of a regular hexagonal array. - - See Also - -------- - armi.reactor.grids.HexGrid.fromPitch - """ - return self._unitSteps[1][1] - - @staticmethod - def indicesToRingPos(i, j): - """ - Convert spatialLocator indices to ring/position. - - One benefit it has is that it never has negative numbers. - - Notes - ----- - Ring, pos index system goes in counterclockwise hex rings. - """ - if i > 0 and j >= 0: - edge = 0 - ring = i + j + 1 - offset = j - elif i <= 0 and j > -i: - edge = 1 - ring = j + 1 - offset = -i - elif i < 0 and j > 0: - edge = 2 - ring = -i + 1 - offset = -j - i - elif i < 0: - edge = 3 - ring = -i - j + 1 - offset = -j - elif i >= 0 and j < -i: - edge = 4 - ring = -j + 1 - offset = i - else: - edge = 5 - ring = i + 1 - offset = i + j - - positionBase = 1 + edge * (ring - 1) - return ring, positionBase + offset - - def getMinimumRings(self, n): - """ - Return the minimum number of rings needed to fit ``n`` objects. - - Notes - ----- - ``self`` is not used because hex grids always behave the same w.r.t. - rings/positions. - """ - return hexagon.numRingsToHoldNumCells(n) - - def getPositionsInRing(self, ring): - """Return the number of positions within a ring.""" - return hexagon.numPositionsInRing(ring) - - def getNeighboringCellIndices(self, i, j=0, k=0): - """ - Return the indices of the immediate neighbors of a mesh point in the plane. - - Note that these neighbors are ordered counter-clockwise beginning from the - 30 or 60 degree direction. Exact direction is dependent on pointedEndUp arg. - """ - return [ - (i + 1, j, k), - (i, j + 1, k), - (i - 1, j + 1, k), - (i - 1, j, k), - (i, j - 1, k), - (i + 1, j - 1, k), - ] - - def getLabel(self, indices): - """ - Hex labels start at 1, and are ring/position based rather than i,j. - - This difference is partially because ring/pos is easier to understand in hex - geometry, and partially because it is used in some codes ARMI originally was focused - on. - """ - ring, pos = self.getRingPos(indices) - if len(indices) == 2: - return Grid.getLabel(self, (ring, pos)) - else: - return Grid.getLabel(self, (ring, pos, indices[2])) - - @staticmethod - def _indicesAndEdgeFromRingAndPos(ring, position): - ring = ring - 1 - pos = position - 1 - - if ring == 0: - if pos != 0: - raise ValueError(f"Position in center ring must be 1, not {position}") - return 0, 0, 0 - - # Edge indicates which edge of the ring in which the hexagon resides. - # Edge 0 is the NE edge, edge 1 is the N edge, etc. - # Offset is (0-based) index of the hexagon in that edge. For instance, - # ring 3, pos 12 resides in edge 5 at index 1; it is the second hexagon - # in ring 3, edge 5. - edge, offset = divmod(pos, ring) # = pos//ring, pos%ring - if edge == 0: - i = ring - offset - j = offset - elif edge == 1: - i = -offset - j = ring - elif edge == 2: - i = -ring - j = -offset + ring - elif edge == 3: - i = -ring + offset - j = -offset - elif edge == 4: - i = offset - j = -ring - elif edge == 5: - i = ring - j = offset - ring - else: - raise ValueError( - "Edge {} is invalid. From ring {}, pos {}".format(edge, ring, pos) - ) - return i, j, edge - - @staticmethod - def getIndicesFromRingAndPos(ring, pos): - i, j, _edge = HexGrid._indicesAndEdgeFromRingAndPos(ring, pos) - return i, j - - def getRingPos(self, indices): - """ - Get 1-based ring and position from normal indices. - - See Also - -------- - getIndicesFromRingAndPos : does the reverse - """ - i, j = indices[:2] - return HexGrid.indicesToRingPos(i, j) - - def overlapsWhichSymmetryLine(self, indices): - """Return a list of which lines of symmetry this is on. - - If none, returns [] - If on a line of symmetry in 1/6 geometry, returns a list containing a 6. - If on a line of symmetry in 1/3 geometry, returns a list containing a 3. - Only the 1/3 core view geometry is actually coded in here right now. - - Being "on" a symmetry line means the line goes through the middle of you. - - """ - i, j = indices[:2] - - if i == 0 and j == 0: - symmetryLine = BOUNDARY_CENTER - elif i > 0 and i == -2 * j: - # edge 1: 1/3 symmetry line (bottom horizontal side in 1/3 core view, theta = 0) - symmetryLine = BOUNDARY_0_DEGREES - elif i == j and i > 0 and j > 0: - # edge 2: 1/6 symmetry line (bisects 1/3 core view, theta = pi/3) - symmetryLine = BOUNDARY_60_DEGREES - elif j == -2 * i and j > 0: - # edge 3: 1/3 symmetry line (left oblique side in 1/3 core view, theta = 2*pi/3) - symmetryLine = BOUNDARY_120_DEGREES - else: - symmetryLine = None - - return symmetryLine - - def getSymmetricEquivalents(self, indices): - if ( - self.symmetry.domain == geometry.DomainType.THIRD_CORE - and self.symmetry.boundary == geometry.BoundaryType.PERIODIC - ): - return HexGrid._getSymmetricIdenticalsThird(indices) - elif self.symmetry.domain == geometry.DomainType.FULL_CORE: - return [] - else: - raise NotImplementedError( - "Unhandled symmetry condition for HexGrid: {}".format( - str(self.symmetry) - ) - ) - - @staticmethod - def _getSymmetricIdenticalsThird(indices): - """This works by rotating the indices by 120 degrees twice, counterclockwise.""" - i, j = indices[:2] - if i == 0 and j == 0: - return [] - identicals = [(-i - j, i), (j, -i - j)] - return identicals - - def triangleCoords(self, indices): - """ - Return 6 coordinate pairs representing the centers of the 6 triangles in a hexagon centered here. - - Ignores z-coordinate and only operates in 2D for now. - """ - xy = self.getCoordinates(indices)[:2] - scale = self.pitch / 3.0 - return xy + scale * TRIANGLES_IN_HEXAGON - - def changePitch(self, newPitchCm): - """Change the hex pitch.""" - side = newPitchCm / math.sqrt(3.0) - self._unitSteps = numpy.array( - ((1.5 * side, 0.0, 0.0), (newPitchCm / 2.0, newPitchCm, 0.0), (0, 0, 0)) - )[self._stepDims] - - def locatorInDomain(self, locator, symmetryOverlap: Optional[bool] = False): - # This will include the "top" 120-degree symmetry lines. This is to support - # adding of edge assemblies. - if self.symmetry.domain == geometry.DomainType.THIRD_CORE: - return self.isInFirstThird(locator, includeTopEdge=symmetryOverlap) - else: - return True - - def isInFirstThird(self, locator, includeTopEdge=False): - """True if locator is in first third of hex grid.""" - ring, pos = self.getRingPos(locator.indices) - if ring == 1: - return True - maxPosTotal = self.getPositionsInRing(ring) - - maxPos1 = ring + ring // 2 - 1 - maxPos2 = maxPosTotal - ring // 2 + 1 - if ring % 2: - # odd ring. Upper edge assem typically not included. - if includeTopEdge: - maxPos1 += 1 - else: - maxPos2 += 1 # make a table to understand this. - if pos <= maxPos1 or pos >= maxPos2: - return True - return False - - def generateSortedHexLocationList(self, nLocs): - """ - Generate a list IndexLocations, sorted based on their distance from the center. - - IndexLocation are taken from a full core. - - Ties between locations with the same distance (e.g. A3001 and A3003) are broken - by ring number then position number. - """ - # first, roughly calculate how many rings need to be created to cover nLocs worth of assemblies - nLocs = int(nLocs) # need to make this an integer - - # next, generate a list of locations and corresponding distances - locs = [] - for ring in range(1, hexagon.numRingsToHoldNumCells(nLocs) + 1): - positions = self.getPositionsInRing(ring) - for position in range(1, positions + 1): - i, j = self.getIndicesFromRingAndPos(ring, position) - locs.append(self[(i, j, 0)]) - # round to avoid differences due to floating point math - locs.sort( - key=lambda loc: ( - round(numpy.linalg.norm(loc.getGlobalCoordinates()), 6), - loc.i, # loc.i=ring - loc.j, - ) - ) # loc.j= pos - return locs[:nLocs] - - # TODO: this is only used by testing and another method that just needs the count of assemblies - # in a ring, not the actual positions - def allPositionsInThird(self, ring, includeEdgeAssems=False): - """ - Returns a list of all the positions in a ring (in the first third). - - Parameters - ---------- - ring : int - The ring to check - includeEdgeAssems : bool, optional - If True, include repeated positions in odd ring numbers. Default: False - - Notes - ----- - Rings start at 1, positions start at 1 - - Returns - ------- - positions : int - """ - positions = [] - for pos in range(1, self.getPositionsInRing(ring) + 1): - i, j = self.getIndicesFromRingAndPos(ring, pos) - loc = IndexLocation(i, j, 0, None) - if self.isInFirstThird(loc, includeEdgeAssems): - positions.append(pos) - return positions +from .hexagonal import HexGrid, COS30, SIN30, TRIANGLES_IN_HEXAGON class ThetaRZGrid(Grid): diff --git a/armi/reactor/grids/hexagonal.py b/armi/reactor/grids/hexagonal.py new file mode 100644 index 000000000..ab7cfd631 --- /dev/null +++ b/armi/reactor/grids/hexagonal.py @@ -0,0 +1,426 @@ +# Copyright 2023 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import math +from typing import Optional + +import numpy + +from armi.utils import hexagon +from armi.reactor import geometry + +from .constants import ( + BOUNDARY_0_DEGREES, + BOUNDARY_120_DEGREES, + BOUNDARY_60_DEGREES, + BOUNDARY_CENTER, +) +from .locations import IndexLocation +from .grid import Grid + +COS30 = math.sqrt(3) / 2.0 +SIN30 = 1.0 / 2.0 +# going CCW from "position 1" (top right) +TRIANGLES_IN_HEXAGON = numpy.array( + [ + (+COS30, SIN30), + (+0, 1.0), + (-COS30, SIN30), + (-COS30, -SIN30), + (+0, -1.0), + (+COS30, -SIN30), + ] +) + + +class HexGrid(Grid): + """ + Has 6 neighbors in plane. + + Notes + ----- + In an axial plane (i, j) are as follows (second one is pointedEndUp):: + + + ( 0, 1) + (-1, 1) ( 1, 0) + ( 0, 0) + (-1, 0) ( 1,-1) + ( 0,-1) + + + ( 0, 1) ( 1, 0) + + (-1, 1) ( 0, 0) ( 1,-1) + + (-1, 0) ( 0,-1) + + .. impl:: ARMI supports a Hexagonal mesh. + :id: IMPL_REACTOR_MESH_2 + :links: REQ_REACTOR_MESH + """ + + @staticmethod + def fromPitch(pitch, numRings=25, armiObject=None, pointedEndUp=False, symmetry=""): + """ + Build a finite step-based 2-D hex grid from a hex pitch in cm. + + Parameters + ---------- + pitch : float + Hex pitch (flat-to-flat) in cm + numRings : int, optional + The number of rings in the grid to pre-populate with locatator objects. + Even if positions are not pre-populated, locators will be generated + there on the fly. + armiObject : ArmiObject, optional + The object that this grid is anchored to (i.e. the reactor for a grid of + assemblies) + pointedEndUp : bool, optional + Rotate the hexagons 30 degrees so that the pointed end faces up instead of + the flat. + symmetry : string, optional + A string representation of the symmetry options for the grid. + + Returns + ------- + HexGrid + A functional hexagonal grid object. + """ + side = pitch / math.sqrt(3.0) + if pointedEndUp: + # rotated 30 degrees CCW from normal + # increases in i move you in x and y + # increases in j also move you in x and y + unitSteps = ( + (pitch / 2.0, -pitch / 2.0, 0), + (1.5 * side, 1.5 * side, 0), + (0, 0, 0), + ) + else: + # x direction is only a function of i because j-axis is vertical. + # y direction is a function of both. + unitSteps = ((1.5 * side, 0.0, 0.0), (pitch / 2.0, pitch, 0.0), (0, 0, 0)) + + return HexGrid( + unitSteps=unitSteps, + unitStepLimits=((-numRings, numRings), (-numRings, numRings), (0, 1)), + armiObject=armiObject, + symmetry=symmetry, + ) + + @property + def pitch(self): + """ + Get the hex-pitch of a regular hexagonal array. + + See Also + -------- + armi.reactor.grids.HexGrid.fromPitch + """ + return self._unitSteps[1][1] + + @staticmethod + def indicesToRingPos(i, j): + """ + Convert spatialLocator indices to ring/position. + + One benefit it has is that it never has negative numbers. + + Notes + ----- + Ring, pos index system goes in counterclockwise hex rings. + """ + if i > 0 and j >= 0: + edge = 0 + ring = i + j + 1 + offset = j + elif i <= 0 and j > -i: + edge = 1 + ring = j + 1 + offset = -i + elif i < 0 and j > 0: + edge = 2 + ring = -i + 1 + offset = -j - i + elif i < 0: + edge = 3 + ring = -i - j + 1 + offset = -j + elif i >= 0 and j < -i: + edge = 4 + ring = -j + 1 + offset = i + else: + edge = 5 + ring = i + 1 + offset = i + j + + positionBase = 1 + edge * (ring - 1) + return ring, positionBase + offset + + def getMinimumRings(self, n): + """ + Return the minimum number of rings needed to fit ``n`` objects. + + Notes + ----- + ``self`` is not used because hex grids always behave the same w.r.t. + rings/positions. + """ + return hexagon.numRingsToHoldNumCells(n) + + def getPositionsInRing(self, ring): + """Return the number of positions within a ring.""" + return hexagon.numPositionsInRing(ring) + + def getNeighboringCellIndices(self, i, j=0, k=0): + """ + Return the indices of the immediate neighbors of a mesh point in the plane. + + Note that these neighbors are ordered counter-clockwise beginning from the + 30 or 60 degree direction. Exact direction is dependent on pointedEndUp arg. + """ + return [ + (i + 1, j, k), + (i, j + 1, k), + (i - 1, j + 1, k), + (i - 1, j, k), + (i, j - 1, k), + (i + 1, j - 1, k), + ] + + def getLabel(self, indices): + """ + Hex labels start at 1, and are ring/position based rather than i,j. + + This difference is partially because ring/pos is easier to understand in hex + geometry, and partially because it is used in some codes ARMI originally was focused + on. + """ + ring, pos = self.getRingPos(indices) + if len(indices) == 2: + return Grid.getLabel(self, (ring, pos)) + else: + return Grid.getLabel(self, (ring, pos, indices[2])) + + @staticmethod + def _indicesAndEdgeFromRingAndPos(ring, position): + ring = ring - 1 + pos = position - 1 + + if ring == 0: + if pos != 0: + raise ValueError(f"Position in center ring must be 1, not {position}") + return 0, 0, 0 + + # Edge indicates which edge of the ring in which the hexagon resides. + # Edge 0 is the NE edge, edge 1 is the N edge, etc. + # Offset is (0-based) index of the hexagon in that edge. For instance, + # ring 3, pos 12 resides in edge 5 at index 1; it is the second hexagon + # in ring 3, edge 5. + edge, offset = divmod(pos, ring) # = pos//ring, pos%ring + if edge == 0: + i = ring - offset + j = offset + elif edge == 1: + i = -offset + j = ring + elif edge == 2: + i = -ring + j = -offset + ring + elif edge == 3: + i = -ring + offset + j = -offset + elif edge == 4: + i = offset + j = -ring + elif edge == 5: + i = ring + j = offset - ring + else: + raise ValueError( + "Edge {} is invalid. From ring {}, pos {}".format(edge, ring, pos) + ) + return i, j, edge + + @staticmethod + def getIndicesFromRingAndPos(ring, pos): + i, j, _edge = HexGrid._indicesAndEdgeFromRingAndPos(ring, pos) + return i, j + + def getRingPos(self, indices): + """ + Get 1-based ring and position from normal indices. + + See Also + -------- + getIndicesFromRingAndPos : does the reverse + """ + i, j = indices[:2] + return HexGrid.indicesToRingPos(i, j) + + def overlapsWhichSymmetryLine(self, indices): + """Return a list of which lines of symmetry this is on. + + If none, returns [] + If on a line of symmetry in 1/6 geometry, returns a list containing a 6. + If on a line of symmetry in 1/3 geometry, returns a list containing a 3. + Only the 1/3 core view geometry is actually coded in here right now. + + Being "on" a symmetry line means the line goes through the middle of you. + + """ + i, j = indices[:2] + + if i == 0 and j == 0: + symmetryLine = BOUNDARY_CENTER + elif i > 0 and i == -2 * j: + # edge 1: 1/3 symmetry line (bottom horizontal side in 1/3 core view, theta = 0) + symmetryLine = BOUNDARY_0_DEGREES + elif i == j and i > 0 and j > 0: + # edge 2: 1/6 symmetry line (bisects 1/3 core view, theta = pi/3) + symmetryLine = BOUNDARY_60_DEGREES + elif j == -2 * i and j > 0: + # edge 3: 1/3 symmetry line (left oblique side in 1/3 core view, theta = 2*pi/3) + symmetryLine = BOUNDARY_120_DEGREES + else: + symmetryLine = None + + return symmetryLine + + def getSymmetricEquivalents(self, indices): + if ( + self.symmetry.domain == geometry.DomainType.THIRD_CORE + and self.symmetry.boundary == geometry.BoundaryType.PERIODIC + ): + return HexGrid._getSymmetricIdenticalsThird(indices) + elif self.symmetry.domain == geometry.DomainType.FULL_CORE: + return [] + else: + raise NotImplementedError( + "Unhandled symmetry condition for HexGrid: {}".format( + str(self.symmetry) + ) + ) + + @staticmethod + def _getSymmetricIdenticalsThird(indices): + """This works by rotating the indices by 120 degrees twice, counterclockwise.""" + i, j = indices[:2] + if i == 0 and j == 0: + return [] + identicals = [(-i - j, i), (j, -i - j)] + return identicals + + def triangleCoords(self, indices): + """ + Return 6 coordinate pairs representing the centers of the 6 triangles in a hexagon centered here. + + Ignores z-coordinate and only operates in 2D for now. + """ + xy = self.getCoordinates(indices)[:2] + scale = self.pitch / 3.0 + return xy + scale * TRIANGLES_IN_HEXAGON + + def changePitch(self, newPitchCm): + """Change the hex pitch.""" + side = newPitchCm / math.sqrt(3.0) + self._unitSteps = numpy.array( + ((1.5 * side, 0.0, 0.0), (newPitchCm / 2.0, newPitchCm, 0.0), (0, 0, 0)) + )[self._stepDims] + + def locatorInDomain(self, locator, symmetryOverlap: Optional[bool] = False): + # This will include the "top" 120-degree symmetry lines. This is to support + # adding of edge assemblies. + if self.symmetry.domain == geometry.DomainType.THIRD_CORE: + return self.isInFirstThird(locator, includeTopEdge=symmetryOverlap) + else: + return True + + def isInFirstThird(self, locator, includeTopEdge=False): + """True if locator is in first third of hex grid.""" + ring, pos = self.getRingPos(locator.indices) + if ring == 1: + return True + maxPosTotal = self.getPositionsInRing(ring) + + maxPos1 = ring + ring // 2 - 1 + maxPos2 = maxPosTotal - ring // 2 + 1 + if ring % 2: + # odd ring. Upper edge assem typically not included. + if includeTopEdge: + maxPos1 += 1 + else: + maxPos2 += 1 # make a table to understand this. + if pos <= maxPos1 or pos >= maxPos2: + return True + return False + + def generateSortedHexLocationList(self, nLocs): + """ + Generate a list IndexLocations, sorted based on their distance from the center. + + IndexLocation are taken from a full core. + + Ties between locations with the same distance (e.g. A3001 and A3003) are broken + by ring number then position number. + """ + # first, roughly calculate how many rings need to be created to cover nLocs worth of assemblies + nLocs = int(nLocs) # need to make this an integer + + # next, generate a list of locations and corresponding distances + locs = [] + for ring in range(1, hexagon.numRingsToHoldNumCells(nLocs) + 1): + positions = self.getPositionsInRing(ring) + for position in range(1, positions + 1): + i, j = self.getIndicesFromRingAndPos(ring, position) + locs.append(self[(i, j, 0)]) + # round to avoid differences due to floating point math + locs.sort( + key=lambda loc: ( + round(numpy.linalg.norm(loc.getGlobalCoordinates()), 6), + loc.i, # loc.i=ring + loc.j, + ) + ) # loc.j= pos + return locs[:nLocs] + + # TODO: this is only used by testing and another method that just needs the count of assemblies + # in a ring, not the actual positions + def allPositionsInThird(self, ring, includeEdgeAssems=False): + """ + Returns a list of all the positions in a ring (in the first third). + + Parameters + ---------- + ring : int + The ring to check + includeEdgeAssems : bool, optional + If True, include repeated positions in odd ring numbers. Default: False + + Notes + ----- + Rings start at 1, positions start at 1 + + Returns + ------- + positions : int + """ + positions = [] + for pos in range(1, self.getPositionsInRing(ring) + 1): + i, j = self.getIndicesFromRingAndPos(ring, pos) + loc = IndexLocation(i, j, 0, None) + if self.isInFirstThird(loc, includeEdgeAssems): + positions.append(pos) + return positions From 38e9180f0d8972e5809996c37b370009c8213276 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 13:21:56 -0700 Subject: [PATCH 08/16] Add armi.reactor.grids.thetarz --- armi/reactor/grids/__init__.py | 110 +--------------------------- armi/reactor/grids/thetarz.py | 126 +++++++++++++++++++++++++++++++++ 2 files changed, 128 insertions(+), 108 deletions(-) create mode 100644 armi/reactor/grids/thetarz.py diff --git a/armi/reactor/grids/__init__.py b/armi/reactor/grids/__init__.py index eff3bd595..7e4c1e648 100644 --- a/armi/reactor/grids/__init__.py +++ b/armi/reactor/grids/__init__.py @@ -60,11 +60,9 @@ while the word **local** refers to within the current coordinate system defined by the current grid. """ -import math from typing import Tuple, Optional -TAU = math.pi * 2.0 - +import numpy from .constants import ( BOUNDARY_CENTER, @@ -83,111 +81,7 @@ from .grid import Grid, GridParameters, _tuplify from .cartesian import CartesianGrid from .hexagonal import HexGrid, COS30, SIN30, TRIANGLES_IN_HEXAGON - - -class ThetaRZGrid(Grid): - """ - A grid characterized by azimuthal, radial, and zeta indices. - - The angular meshes are limited to 0 to 2pi radians. R and Zeta are as in other - meshes. - - See Figure 2.2 in Derstine 1984, ANL. [DIF3D]_. - - .. impl:: ARMI supports an RZTheta mesh. - :id: IMPL_REACTOR_MESH_3 - :links: REQ_REACTOR_MESH - """ - - @staticmethod - def fromGeom(geom, armiObject=None): - """ - Build 2-D R-theta grid based on a Geometry object. - - Parameters - ---------- - geomInfo : list - list of ((indices), assemName) tuples for all positions in core with input - in radians - - See Also - -------- - armi.reactor.systemLayoutInput.SystemLayoutInput.readGeomXML : produces the geomInfo - structure - - Examples - -------- - >>> grid = grids.ThetaRZGrid.fromGeom(geomInfo) - """ - allIndices = [ - indices for indices, _assemName in geom.assemTypeByIndices.items() - ] - - # create ordered lists of all unique theta and R points - thetas, radii = set(), set() - for rad1, rad2, theta1, theta2, _numAzi, _numRadial in allIndices: - radii.add(rad1) - radii.add(rad2) - thetas.add(theta1) - thetas.add(theta2) - radii = numpy.array(sorted(radii), dtype=numpy.float64) - thetaRadians = numpy.array(sorted(thetas), dtype=numpy.float64) - - return ThetaRZGrid( - bounds=(thetaRadians, radii, (0.0, 0.0)), armiObject=armiObject - ) - - def getRingPos(self, indices): - return (indices[1] + 1, indices[0] + 1) - - @staticmethod - def getIndicesFromRingAndPos(ring, pos): - return (pos - 1, ring - 1) - - def getCoordinates(self, indices, nativeCoords=False): - meshCoords = theta, r, z = Grid.getCoordinates(self, indices) - if not 0 <= theta <= TAU: - raise ValueError("Invalid theta value: {}. Check mesh.".format(theta)) - if nativeCoords: - # return Theta, R, Z values directly. - return meshCoords - else: - # return x, y ,z - return numpy.array((r * math.cos(theta), r * math.sin(theta), z)) - - def indicesOfBounds(self, rad0, rad1, theta0, theta1, sigma=1e-4): - """ - Return indices corresponding to upper and lower radial and theta bounds. - - Parameters - ---------- - rad0 : float - inner radius of control volume - rad1 : float - outer radius of control volume - theta0 : float - inner azimuthal location of control volume in radians - theta1 : float - inner azimuthal of control volume in radians - sigma: float - acceptable relative error (i.e. if one of the positions in the mesh are within - this error it'll act the same if it matches a position in the mesh) - - Returns - ------- - tuple : i, j, k of given bounds - """ - i = int(numpy.abs(self._bounds[0] - theta0).argmin()) - j = int(numpy.abs(self._bounds[1] - rad0).argmin()) - - return (i, j, 0) - - def locatorInDomain(self, locator, symmetryOverlap: Optional[bool] = False): - """ - ThetaRZGrids do not check for bounds, though they could if that becomes a - problem. - """ - return True +from .thetarz import ThetaRZGrid, TAU def axialUnitGrid(numCells, armiObject=None): diff --git a/armi/reactor/grids/thetarz.py b/armi/reactor/grids/thetarz.py new file mode 100644 index 000000000..9551f064b --- /dev/null +++ b/armi/reactor/grids/thetarz.py @@ -0,0 +1,126 @@ +# Copyright 2023 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +import math +from typing import Optional + +import numpy + +from .grid import Grid + +TAU = math.pi * 2.0 + + +class ThetaRZGrid(Grid): + """ + A grid characterized by azimuthal, radial, and zeta indices. + + The angular meshes are limited to 0 to 2pi radians. R and Zeta are as in other + meshes. + + See Figure 2.2 in Derstine 1984, ANL. [DIF3D]_. + + .. impl:: ARMI supports an RZTheta mesh. + :id: IMPL_REACTOR_MESH_3 + :links: REQ_REACTOR_MESH + """ + + @staticmethod + def fromGeom(geom, armiObject=None): + """ + Build 2-D R-theta grid based on a Geometry object. + + Parameters + ---------- + geomInfo : list + list of ((indices), assemName) tuples for all positions in core with input + in radians + + See Also + -------- + armi.reactor.systemLayoutInput.SystemLayoutInput.readGeomXML : produces the geomInfo + structure + + Examples + -------- + >>> grid = grids.ThetaRZGrid.fromGeom(geomInfo) + """ + allIndices = [ + indices for indices, _assemName in geom.assemTypeByIndices.items() + ] + + # create ordered lists of all unique theta and R points + thetas, radii = set(), set() + for rad1, rad2, theta1, theta2, _numAzi, _numRadial in allIndices: + radii.add(rad1) + radii.add(rad2) + thetas.add(theta1) + thetas.add(theta2) + radii = numpy.array(sorted(radii), dtype=numpy.float64) + thetaRadians = numpy.array(sorted(thetas), dtype=numpy.float64) + + return ThetaRZGrid( + bounds=(thetaRadians, radii, (0.0, 0.0)), armiObject=armiObject + ) + + def getRingPos(self, indices): + return (indices[1] + 1, indices[0] + 1) + + @staticmethod + def getIndicesFromRingAndPos(ring, pos): + return (pos - 1, ring - 1) + + def getCoordinates(self, indices, nativeCoords=False): + meshCoords = theta, r, z = Grid.getCoordinates(self, indices) + if not 0 <= theta <= TAU: + raise ValueError("Invalid theta value: {}. Check mesh.".format(theta)) + if nativeCoords: + # return Theta, R, Z values directly. + return meshCoords + else: + # return x, y ,z + return numpy.array((r * math.cos(theta), r * math.sin(theta), z)) + + def indicesOfBounds(self, rad0, rad1, theta0, theta1, sigma=1e-4): + """ + Return indices corresponding to upper and lower radial and theta bounds. + + Parameters + ---------- + rad0 : float + inner radius of control volume + rad1 : float + outer radius of control volume + theta0 : float + inner azimuthal location of control volume in radians + theta1 : float + inner azimuthal of control volume in radians + sigma: float + acceptable relative error (i.e. if one of the positions in the mesh are within + this error it'll act the same if it matches a position in the mesh) + + Returns + ------- + tuple : i, j, k of given bounds + """ + i = int(numpy.abs(self._bounds[0] - theta0).argmin()) + j = int(numpy.abs(self._bounds[1] - rad0).argmin()) + + return (i, j, 0) + + def locatorInDomain(self, locator, symmetryOverlap: Optional[bool] = False): + """ + ThetaRZGrids do not check for bounds, though they could if that becomes a + problem. + """ + return True From 62bac65f346bb93d8e2b3d3e578914fec841b70a Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 16:05:02 -0700 Subject: [PATCH 09/16] Remove six.cPickle -> pickle in test_grids.py --- armi/reactor/grids/tests/test_grids.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/armi/reactor/grids/tests/test_grids.py b/armi/reactor/grids/tests/test_grids.py index 79f0543e8..67013eaa8 100644 --- a/armi/reactor/grids/tests/test_grids.py +++ b/armi/reactor/grids/tests/test_grids.py @@ -16,9 +16,9 @@ from io import BytesIO import math import unittest +import pickle from numpy.testing import assert_allclose -from six.moves import cPickle import numpy from armi.reactor import geometry @@ -348,11 +348,11 @@ def test_buildLocations(self): def test_is_pickleable(self): grid = grids.HexGrid.fromPitch(1.0, numRings=3) loc = grid[1, 1, 0] - for protocol in range(cPickle.HIGHEST_PROTOCOL + 1): + for protocol in range(pickle.HIGHEST_PROTOCOL + 1): buf = BytesIO() - cPickle.dump(loc, buf, protocol=protocol) + pickle.dump(loc, buf, protocol=protocol) buf.seek(0) - newLoc = cPickle.load(buf) + newLoc = pickle.load(buf) assert_allclose(loc.indices, newLoc.indices) def test_adjustPitch(self): From 33f6babf7aa392c17093d19d4e3a2ea847673c61 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 16:19:27 -0700 Subject: [PATCH 10/16] Add armi.reactor.grids.axial with AxialGrid class Replaces grids.axialUnitGrid --- armi/reactor/grids/__init__.py | 17 +-------- armi/reactor/grids/axial.py | 49 ++++++++++++++++++++++++++ armi/reactor/grids/tests/test_grids.py | 26 ++++++++++++-- 3 files changed, 73 insertions(+), 19 deletions(-) create mode 100644 armi/reactor/grids/axial.py diff --git a/armi/reactor/grids/__init__.py b/armi/reactor/grids/__init__.py index 7e4c1e648..3e531eb14 100644 --- a/armi/reactor/grids/__init__.py +++ b/armi/reactor/grids/__init__.py @@ -62,8 +62,6 @@ """ from typing import Tuple, Optional -import numpy - from .constants import ( BOUNDARY_CENTER, BOUNDARY_0_DEGREES, @@ -79,25 +77,12 @@ ) from .grid import Grid, GridParameters, _tuplify +from .axial import AxialGrid, axialUnitGrid from .cartesian import CartesianGrid from .hexagonal import HexGrid, COS30, SIN30, TRIANGLES_IN_HEXAGON from .thetarz import ThetaRZGrid, TAU -def axialUnitGrid(numCells, armiObject=None): - """ - Build a 1-D unit grid in the k-direction based on a number of times. Each mesh is 1cm wide. - - numCells + 1 mesh boundaries are added, since one block would require a bottom and a - top. - """ - # need float bounds or else we truncate integers - return Grid( - bounds=(None, None, numpy.arange(numCells + 1, dtype=numpy.float64)), - armiObject=armiObject, - ) - - def locatorLabelToIndices(label: str) -> Tuple[int, int, Optional[int]]: """ Convert a locator label to numerical i,j,k indices. diff --git a/armi/reactor/grids/axial.py b/armi/reactor/grids/axial.py new file mode 100644 index 000000000..f719205d4 --- /dev/null +++ b/armi/reactor/grids/axial.py @@ -0,0 +1,49 @@ +from typing import Optional, TYPE_CHECKING + +import numpy + +from .grid import Grid + +if TYPE_CHECKING: + from armi.reactor.composites import ArmiObject + + +class AxialGrid(Grid): + """1-D grid in the k-direction (z) + + .. note::: + + It is recommended to use :meth:`fromNCells` rather than calling + the ``__init_`` constructor directly + + """ + + @classmethod + def fromNCells( + cls, numCells: int, armiObject: Optional["ArmiObject"] = None + ) -> "AxialGrid": + """Produces an unit grid where each bin is 1-cm tall + + ``numCells + 1`` mesh boundaries are added, since one block would + require a bottom and a top. + + """ + # Need float bounds or else we truncate integers + return cls( + bounds=(None, None, numpy.arange(numCells + 1, dtype=numpy.float64)), + armiObject=armiObject, + ) + + +def axialUnitGrid( + numCells: int, armiObject: Optional["ArmiObject"] = None +) -> AxialGrid: + """ + Build a 1-D unit grid in the k-direction based on a number of times. Each mesh is 1cm wide. + + .. deprecated:: + + Use :class:`AxialUnitGrid` class instead + + """ + return AxialGrid.fromNCells(numCells, armiObject) diff --git a/armi/reactor/grids/tests/test_grids.py b/armi/reactor/grids/tests/test_grids.py index 67013eaa8..45a42eb6f 100644 --- a/armi/reactor/grids/tests/test_grids.py +++ b/armi/reactor/grids/tests/test_grids.py @@ -18,7 +18,7 @@ import unittest import pickle -from numpy.testing import assert_allclose +from numpy.testing import assert_allclose, assert_array_equal import numpy from armi.reactor import geometry @@ -161,7 +161,7 @@ def test_isAxialOnly(self): grid = grids.HexGrid.fromPitch(1.0, numRings=3) self.assertEqual(grid.isAxialOnly, False) - grid2 = grids.axialUnitGrid(10) + grid2 = grids.AxialGrid.fromNCells(10) self.assertEqual(grid2.isAxialOnly, True) def test_lookupFactory(self): @@ -368,7 +368,7 @@ def test_badIndices(self): # this is actually ok because step-defined grids are infinite self.assertEqual(grid.getCoordinates((-100, 2000, 5))[2], 0.0) - grid = grids.axialUnitGrid(10) + grid = grids.AxialGrid.fromNCells(10) with self.assertRaises(IndexError): grid.getCoordinates((0, 5, -1)) @@ -657,3 +657,23 @@ def test_symmetry(self): ) with self.assertRaises(NotImplementedError): grid.getSymmetricEquivalents((5, 6)) + + +class TestAxialGrid(unittest.TestCase): + def test_simpleBounds(self): + N_CELLS = 5 + g = grids.AxialGrid.fromNCells(N_CELLS) + _x, _y, z = g.getBounds() + self.assertEqual(len(z), N_CELLS + 1) + assert_array_equal(z, [0, 1, 2, 3, 4, 5]) + self.assertTrue(g.isAxialOnly) + + def test_getLocations(self): + N_CELLS = 10 + g = grids.AxialGrid.fromNCells(N_CELLS) + for count in range(N_CELLS): + index = g[(0, 0, count)] + x, y, z = index.getLocalCoordinates() + self.assertEqual(x, 0.0) + self.assertEqual(y, 0.0) + self.assertEqual(z, count + 0.5) From 1df8cf1eceb9a919ff4b73b7105e8456ee477ca6 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 16:25:38 -0700 Subject: [PATCH 11/16] Move duplicated grids.addingIsValid to grids.locations --- armi/reactor/grids/__init__.py | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/armi/reactor/grids/__init__.py b/armi/reactor/grids/__init__.py index 3e531eb14..137a3b633 100644 --- a/armi/reactor/grids/__init__.py +++ b/armi/reactor/grids/__init__.py @@ -74,6 +74,7 @@ IndexLocation, MultiIndexLocation, CoordinateLocation, + addingIsValid, ) from .grid import Grid, GridParameters, _tuplify @@ -93,15 +94,3 @@ def locatorLabelToIndices(label: str) -> Tuple[int, int, Optional[int]]: if len(intVals) == 2: intVals = (intVals[0], intVals[1], None) return intVals - - -def addingIsValid(myGrid, parentGrid): - """ - True if adding a indices from one grid to another is considered valid. - - In ARMI we allow the addition of a 1-D axial grid with a 2-D grid. - We do not allow any other kind of adding. This enables the 2D/1D - grid layout in Assemblies/Blocks but does not allow 2D indexing - in pins to become inconsistent. - """ - return myGrid.isAxialOnly and not parentGrid.isAxialOnly From a312849fc95d2f539819d636e56357670c6d916d Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Tue, 1 Aug 2023 16:31:34 -0700 Subject: [PATCH 12/16] Refactor Grid to be ABC, intermediate StructuredGrid Heavier lift, but the classes defined in `armi.reactor.grids` follow this inheritance pattern - `Grid(abc.ABC)` :: abstract interface for _a_ grid, or things that reside in space - `StructuredGrid(Grid)` :: still abstract, but for things that have some structure, known to ARMI as unit steps and bounds - `HexGrid(StructuredGrid)` - `CartesianGrid(StructuredGrid) - `ThetaRZGrid(StructuredGrid)` - though how "structured" is this really? - `AxialGrid(StructuredGrid)` In some cases, the shared interface between structured grids falls apart. For example, `HexGrid.pitch` returns a single float, while `CartesianGrid.pitch` is a 2-ple of float, one for x and one for y. And the concept of ring and position appears to be stressed since the inception of cartesian grids. But parts of ARMI require an interface and so an attempt was made to define that interface. The decision was made to not have `StructuredGrid` the base class, and to introduce `Grid` as the base intentionally. There are cases for some multiphysics codes where a grid may still be "structured" but not with the level of structure the previous `Grid` implementations expect. One example would be an irregular cartesian grid, where the unit cell spacing in one or more dimensions may change. So the cell spacing is a function of the cell position in that dimension. --- armi/reactor/grids/__init__.py | 3 +- armi/reactor/grids/axial.py | 50 +- armi/reactor/grids/cartesian.py | 68 ++- armi/reactor/grids/grid.py | 637 +++++-------------------- armi/reactor/grids/hexagonal.py | 110 +++-- armi/reactor/grids/structuredgrid.py | 514 ++++++++++++++++++++ armi/reactor/grids/tests/test_grids.py | 78 ++- armi/reactor/grids/thetarz.py | 76 ++- 8 files changed, 900 insertions(+), 636 deletions(-) create mode 100644 armi/reactor/grids/structuredgrid.py diff --git a/armi/reactor/grids/__init__.py b/armi/reactor/grids/__init__.py index 137a3b633..598295a65 100644 --- a/armi/reactor/grids/__init__.py +++ b/armi/reactor/grids/__init__.py @@ -77,7 +77,8 @@ addingIsValid, ) -from .grid import Grid, GridParameters, _tuplify +from .grid import Grid +from .structuredgrid import StructuredGrid, GridParameters, _tuplify from .axial import AxialGrid, axialUnitGrid from .cartesian import CartesianGrid from .hexagonal import HexGrid, COS30, SIN30, TRIANGLES_IN_HEXAGON diff --git a/armi/reactor/grids/axial.py b/armi/reactor/grids/axial.py index f719205d4..69b8d35de 100644 --- a/armi/reactor/grids/axial.py +++ b/armi/reactor/grids/axial.py @@ -1,14 +1,16 @@ -from typing import Optional, TYPE_CHECKING +from typing import List, Optional, TYPE_CHECKING, NoReturn +import warnings import numpy -from .grid import Grid +from .locations import IJType, LocationBase +from .structuredgrid import StructuredGrid if TYPE_CHECKING: from armi.reactor.composites import ArmiObject -class AxialGrid(Grid): +class AxialGrid(StructuredGrid): """1-D grid in the k-direction (z) .. note::: @@ -34,6 +36,43 @@ def fromNCells( armiObject=armiObject, ) + @staticmethod + def getSymmetricEquivalents(indices: IJType) -> List[IJType]: + return [] + + @staticmethod + def locatorInDomain( + locator: LocationBase, symmetryOverlap: Optional[bool] = False + ) -> NoReturn: + raise NotImplementedError + + @staticmethod + def getIndicesFromRingAndPos(ring: int, pos: int) -> NoReturn: + raise NotImplementedError + + @staticmethod + def getMinimumRings(n: int) -> NoReturn: + raise NotImplementedError + + @staticmethod + def getPositionsInRing(ring: int) -> NoReturn: + raise NotImplementedError + + @staticmethod + def overlapsWhichSymmetryLine(indices: IJType) -> None: + return None + + @property + def pitch(self) -> float: + """Grid spacing in the z-direction + + Returns + ------- + float + Pitch in cm + + """ + def axialUnitGrid( numCells: int, armiObject: Optional["ArmiObject"] = None @@ -46,4 +85,9 @@ def axialUnitGrid( Use :class:`AxialUnitGrid` class instead """ + warnings.warn( + "Use grids.AxialGrid class rather than function", + PendingDeprecationWarning, + stacklevel=2, + ) return AxialGrid.fromNCells(numCells, armiObject) diff --git a/armi/reactor/grids/cartesian.py b/armi/reactor/grids/cartesian.py index 99a6ded09..5ce795bac 100644 --- a/armi/reactor/grids/cartesian.py +++ b/armi/reactor/grids/cartesian.py @@ -1,29 +1,22 @@ -# Copyright 2023 TerraPower, LLC -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. import itertools -from typing import Optional +from typing import Optional, NoReturn, Tuple import numpy from armi.reactor import geometry -from .grid import Grid +from .locations import IJType +from .structuredgrid import StructuredGrid -class CartesianGrid(Grid): +class CartesianGrid(StructuredGrid): """ - Grid class representing a conformal Cartesian mesh. + Grid class representing a conformal Cartesian mesh + + .. note:: + + It is recommended to call :meth:`fromRectangle` to construct, + rather than directly constructing with ``__init__`` Notes ----- @@ -103,6 +96,17 @@ def fromRectangle( symmetry=symmetry, ) + def overlapsWhichSymmetryLine(self, indices: IJType) -> None: + """Return lines of symmetry position at a given index can be found + + .. warning:: + + This is not really implemented, but parts of ARMI need it to + not fail, so it always returns None. + + """ + return None + def getRingPos(self, indices): """ Return ring and position from indices. @@ -159,14 +163,14 @@ def getRingPos(self, indices): return (int(ring) + 1, int(pos) + 1) @staticmethod - def getIndicesFromRingAndPos(ring, pos): + def getIndicesFromRingAndPos(ring: int, pos: int) -> NoReturn: """Not implemented for Cartesian-see getRingPos notes.""" raise NotImplementedError( "Cartesian should not need need ring/pos, use i, j indices." "See getRingPos doc string notes for more information/example." ) - def getMinimumRings(self, n): + def getMinimumRings(self, n: int) -> int: """Return the minimum number of rings needed to fit ``n`` objects.""" numPositions = 0 ring = 0 @@ -178,10 +182,15 @@ def getMinimumRings(self, n): return ring - def getPositionsInRing(self, ring): + def getPositionsInRing(self, ring: int) -> int: """ Return the number of positions within a ring. + Parameters + ---------- + ring : int + Ring in question + Notes ----- The number of positions within a ring will change @@ -204,7 +213,7 @@ def locatorInDomain(self, locator, symmetryOverlap: Optional[bool] = False): else: return True - def changePitch(self, xw, yw): + def changePitch(self, xw: float, yw: float): """ Change the pitch of a Cartesian grid. @@ -278,3 +287,20 @@ def getSymmetricEquivalents(self, indices): def _isThroughCenter(self): """Return whether the central cells are split through the middle for symmetry.""" return all(self._offset == [0, 0, 0]) + + @property + def pitch(self) -> Tuple[float, float]: + """Grid pitch in the x and y dimension + + Returns + ------- + float + x-pitch (cm) + float + y-pitch (cm) + + """ + pitch = (self._unitSteps[0][0], self._unitSteps[1][1]) + if pitch[0] == 0: + raise ValueError(f"Grid {self} does not have a defined pitch.") + return pitch diff --git a/armi/reactor/grids/grid.py b/armi/reactor/grids/grid.py index 368ae644b..0e32ea52c 100644 --- a/armi/reactor/grids/grid.py +++ b/armi/reactor/grids/grid.py @@ -1,287 +1,91 @@ -# Copyright 2023 TerraPower, LLC -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -import itertools -import collections -from typing import Optional, List, Sequence, Union, Tuple +from abc import ABC, abstractmethod +from typing import Union, Optional, Hashable, TYPE_CHECKING, Dict, Iterable, Tuple, List import numpy from armi.reactor import geometry -from .locations import IndexLocation, LocationBase, MultiIndexLocation +from .locations import LocationBase, IndexLocation, IJType, IJKType -# data structure for database-serialization of grids -GridParameters = collections.namedtuple( - "GridParameters", - ("unitSteps", "bounds", "unitStepLimits", "offset", "geomType", "symmetry"), -) +if TYPE_CHECKING: + from armi.reactor.composites import ArmiObject -class Grid: - """ - A connected set of cells characterized by indices mapping to space and vice versa. +class Grid(ABC): + """Base class that defines the interface for grids - The cells may be characterized by any mixture of regular repeating steps and - user-defined steps in any direction. + Most work will be done with structured grids, e.g., hexagonal grid, Cartesian grids, + but some physics codes accept irregular or unstructured grids. Consider + a Cartesian grid but with variable stepping between cells, where ``dx`` may not be + constant. - For example, a 2-D hex lattice has constant, regular steps whereas a 3-D hex mesh - may have user-defined axial meshes. Similar for Cartesian, RZT, etc. + So here, we define an interface so things that rely on grids can worry less + about how the location data are stored. Parameters ---------- - unitSteps : tuple of tuples, optional - Describes the grid spatially as a function on indices. - Each tuple describes how each ``(x,y,or z)`` dimension is influenced by - ``(i,j,k)``. In other words, it is:: - - (dxi, dxj, jxk), (dyi, dyj, dyk), (dzi, dzj, dzk) - - where ``dmn`` is the distance (in cm) that dimension ``m`` will change as a - function of index ``n``. - - Unit steps are used as a generic method for defining repetitive grids in a - variety of geometries, including hexagonal and Cartesian. The tuples are not - vectors in the direction of the translation, but rather grouped by direction. If - the bounds argument is described for a direction, the bounds will be used rather - than the unit step information. The default of (0, 0, 0) makes all dimensions - insensitive to indices since the coordinates are calculated by the dot product - of this and the indices. With this default, any dimension that is desired to - change with indices should be defined with bounds. RZtheta grids are created - exclusively with bounds. - - bounds : 3-tuple - Absolute increasing bounds in cm including endpoints of a non-uniform grid. - Each item represents the boundaries in the associated direction. Use Nones when - unitSteps should be applied instead. Most useful for thetaRZ grids or other - non-uniform grids. - - unitStepLimits : 3-tuple - The limit of the steps in all three directions. This constrains step-defined - grids to be finite so we can populate them with SpatialLocator objects. - - offset : 3-tuple, optional - Offset in cm for each axis. By default the center of the (0,0,0)-th object is in - the center of the grid. Offsets can move it so that the (0,0,0)-th object can - be fully within a quadrant (i.e. in a Cartesian grid). - - armiObject : ArmiObject, optional - The ArmiObject that this grid describes. For example if it's a 1-D assembly - grid, the armiObject is the assembly. Note that ``self.armiObject.spatialGrid`` - is ``self``. - - Examples - -------- - A 2D a rectangular grid with width (x) 2 and height (y) 3 would be:: - - >>> grid = Grid(unitSteps=((2, 0, 0), (0, 3, 0),(0, 0, 0))) - - A regular hex grid with pitch 1 is:: - - >>> grid = Grid(unitSteps= ((sqrt(3)/2, 0.0, 0.0), (0.5, 1.0, 0.0), (0, 0, 0)) - - .. note:: For this unit hex the magnitude of the vector constructed using the - 0th index of each tuple is 1.0. - - Notes - ----- - Each dimension must either be defined through unitSteps or bounds. - The combination of unitSteps with bounds was settled upon after some struggle to - have one unified definition of a grid (i.e. just bounds). A hexagonal grid is - somewhat challenging to represent with bounds because the axes are not orthogonal, - so a unit-direction vector plus bounds would be required. And then the bounds would - be wasted space because they can be derived simply by unit steps. Memory efficiency - is important in this object so the compact representation of - unitSteps-when-possible, bounds-otherwise was settled upon. - - Design considerations include: - - * unitSteps are more intuitive as operations starting from the center of a cell, - particularly with hexagons and rectangles. Otherwise the 0,0 position of a hexagon - in the center of 1/3-symmetric hexagon is at the phantom bottom left of the - hexagon. - - * Users generally prefer to input mesh bounds rather than centers (e.g. starting at - 0.5 instead of 0.0 in a unit mesh is weird). - - * If we store bounds, computing bounds is simple and computing centers takes ~2x the - effort. If we store centers, it's the opposite. - - * Regardless of how we store things, we'll need a Grid that has the lower-left - assembly fully inside the problem (i.e. for full core Cartesian) as well as - another one that has the lower-left assembly half-way or quarter-way sliced off - (for 1/2, 1/4, and 1/8 symmetries). The ``offset`` parameter handles this. - - * Looking up mesh boundaries (to define a mesh in another code) is generally more - common than looking up centers (for plotting or measuring distance). - - * A grid can be anchored to the object that it is in with a backreference. This - gives it the ability to traverse the composite tree and map local to global - locations without having to duplicate the composite pattern on grids. This remains - optional so grids can be used for non-reactor-package reasons. It may seem - slightly cleaner to set the armiObject to the parent's spatialLocator itself - but the major disadvantage of this is that when an object moves, the armiObject - would have to be updated. By anchoring directly to Composite objects, the parent - is always up to date no matter where or how things get moved. - - * Unit step calculations use dot products and must not be polluted by the bound - indices. Thus we reduce the size of the unitSteps tuple accordingly. - - .. impl:: ARMI supports a number of structured mesh options. - :id: IMPL_REACTOR_MESH_0 - :links: REQ_REACTOR_MESH + geomType : str or armi.reactor.geometry.GeomType + Underlying geometric representation + symmetry : str or armi.reactor.geometry.SymmetryType + Symmetry conditions + armiObject : optional, armi.reactor.composites.ArmiObject + If given, what is this grid attached to or what does it describe? + Something like a :class:`armi.reactor.Core` + """ + _geomType: str + _symmetry: str + armiObject: Optional["ArmiObject"] + def __init__( self, - unitSteps=(0, 0, 0), - bounds=(None, None, None), - unitStepLimits=((0, 1), (0, 1), (0, 1)), - offset=None, - geomType="", - symmetry="", - armiObject=None, + geomType: Union[str, geometry.GeomType] = "", + symmetry: Union[str, geometry.SymmetryType] = "", + armiObject: Optional["ArmiObject"] = None, ): - # these lists contain the indices representing which dimensions for which steps - # are used, or for which bounds are used. index 0 is x direction, etc. - self._boundDims = [] - self._stepDims = [] - for dimensionIndex, bound in enumerate(bounds): - if bound is None: - self._stepDims.append(dimensionIndex) - else: - self._boundDims.append(dimensionIndex) - - # numpy prefers tuples like this to do slicing on arrays - self._boundDims = (tuple(self._boundDims),) - self._stepDims = (tuple(self._stepDims),) - - unitSteps = _tuplify(unitSteps) - - self._bounds = bounds - self._unitStepLimits = _tuplify(unitStepLimits) - - # only represent unit steps in dimensions they're being used so as to not - # pollute the dot product. This may reduce the length of this from 3 to 2 or 1 - self._unitSteps = numpy.array(unitSteps)[self._stepDims] - self._offset = numpy.zeros(3) if offset is None else numpy.array(offset) - - self._locations = {} # indices -> SpatialLocator map - self.armiObject = armiObject - self.buildLocations() # locations are owned by a grid, so the grid builds them. - self._backup = None # for retainState - - (_ii, iLen), (_ji, jLen), (_ki, kLen) = self.getIndexBounds() - # True if only contains k-cells. - self.isAxialOnly = iLen == jLen == 1 and kLen > 1 - # geometric metadata encapsulated here because it's related to the grid. # They do not impact the grid object itself. # Notice that these are stored using their string representations, rather than # the GridType enum. This avoids the danger of deserializing an enum value from # an old version of the code that may have had different numeric values. - self._geomType: str = str(geomType) - self._symmetry: str = str(symmetry) - - def reduce(self): - """ - Return the set of arguments used to create this Grid. - - This is very much like the argument tuple from ``__reduce__``, but we do not - implement ``__reduce__`` for real, because we are generally happy with - ``__getstate__`` and ``__setstate__`` for pickling purposes. However, getting - these arguments to ``__init__`` is useful for storing Grids to the database, as - they are more stable (less likely to change) than the actual internal state of - the objects. - - The return value should be hashable, such that a set of these can be created. - """ - offset = None if not self._offset.any() else tuple(self._offset) - - bounds = _tuplify(self._bounds) - - # recreate a constructor-friendly version of `_unitSteps` from live data (may have been reduced from - # length 3 to length 2 or 1 based on mixing the step-based definition and the bounds-based definition - # described in Design Considerations above.) - # We don't just save the original tuple passed in because that may miss transformations that - # occur between instantiation and reduction. - unitSteps = [] - compressedSteps = list(self._unitSteps[:]) - for i in range(3): - # Recall _stepDims are stored as a single-value tuple (for numpy indexing) - # So this just is grabbing the actual data. - if i in self._stepDims[0]: - unitSteps.append(compressedSteps.pop(0)) - else: - # Add dummy value which will never get used (it gets reduced away) - unitSteps.append(0) - unitSteps = _tuplify(unitSteps) - - return GridParameters( - unitSteps, - bounds, - self._unitStepLimits, - offset, - self._geomType, - self._symmetry, - ) + self.geomType = geomType + self.symmetry = symmetry + self.armiObject = armiObject + self._backup = None @property def geomType(self) -> geometry.GeomType: + """Geometric representation""" return geometry.GeomType.fromStr(self._geomType) @geomType.setter def geomType(self, geomType: Union[str, geometry.GeomType]): - self._geomType = str(geometry.GeomType.fromAny(geomType)) + if geomType: + self._geomType = str(geometry.GeomType.fromAny(geomType)) + else: + self._geomType = "" @property - def symmetry(self) -> geometry.SymmetryType: + def symmetry(self) -> str: + """Symmetry applied to the grid""" return geometry.SymmetryType.fromStr(self._symmetry) @symmetry.setter def symmetry(self, symmetry: Union[str, geometry.SymmetryType]): - self._symmetry = str(geometry.SymmetryType.fromAny(symmetry)) + if symmetry: + self._symmetry = str(geometry.SymmetryType.fromAny(symmetry)) + else: + self._symmetry = "" - @property - def offset(self): - return self._offset - - @offset.setter - def offset(self, offset): - self._offset = offset - - def __repr__(self): - msg = ( - ["<{} -- {}\nBounds:\n".format(self.__class__.__name__, id(self))] - + [" {}\n".format(b) for b in self._bounds] - + ["Steps:\n"] - + [" {}\n".format(b) for b in self._unitSteps] - + [ - "Anchor: {}\n".format(self.armiObject), - "Offset: {}\n".format(self._offset), - "Num Locations: {}>".format(len(self)), - ] - ) - return "".join(msg) - - def __getstate__(self): + def __getstate__(self) -> Dict: """ - Pickling removes reference to ``armiObject``. + Pickling removes reference to ``armiObject`` - Removing the ``armiObject`` allows us to pickle an assembly without pickling the entire - reactor. An ``Assembly.spatialLocator.grid.armiObject`` is the reactor, by removing the link - here, we still have spatial orientation, but are not required to pickle the entire reactor - to pickle an assembly. + Removing the ``armiObject`` allows us to pickle an assembly without pickling + the entire reactor. An ``Assembly.spatialLocator.grid.armiObject`` is the + reactor, by removing the link here, we still have spatial orientation, but are + not required to pickle the entire reactor to pickle an assembly. This relies on the ``armiObject.__setstate__`` to assign itself. """ @@ -290,91 +94,31 @@ def __getstate__(self): return state - def __setstate__(self, state): + def __setstate__(self, state: Dict): """ - Pickling removes reference to ``armiObject``. + Pickling removes reference to ``armiObject`` This relies on the ``ArmiObject.__setstate__`` to assign itself. """ self.__dict__.update(state) - for _indices, locator in self.items(): + for _index, locator in self.items(): locator._grid = self - def __getitem__(self, ijk: Union[Tuple[int, int, int], List[Tuple[int, int, int]]]): - """ - Get a location by (i, j, k) indices. If it does not exist, create a new one and return it. - - Parameters - ---------- - ijk : tuple of indices or list of the same - If provided a tuple, an IndexLocation will be created (if necessary) and - returned. If provided a list, each element will create a new IndexLocation - (if necessary), and a MultiIndexLocation containing all of the passed - indices will be returned. - - - Notes - ----- - The method is defaultdict-like, in that it will create a new location on the fly. However, - the class itself is not really a dictionary, it is just index-able. For example, there is no - desire to have a ``__setitem__`` method, because the only way to create a location is by - retrieving it or through ``buildLocations``. - """ - try: - return self._locations[ijk] - except (KeyError, TypeError): - pass - - if isinstance(ijk, tuple): - i, j, k = ijk - val = IndexLocation(i, j, k, self) - self._locations[ijk] = val - elif isinstance(ijk, list): - val = MultiIndexLocation(self) - locators = [self[idx] for idx in ijk] - val.extend(locators) - else: - raise TypeError( - "Unsupported index type `{}` for `{}`".format(type(ijk), ijk) - ) - return val + @property + @abstractmethod + def isAxialOnly(self) -> bool: + """Indicate to parts of ARMI if this Grid handles only axial cells""" - def __len__(self): - return len(self._locations) + @abstractmethod + def __len__(self) -> int: + """Number of items in the grid""" - def items(self): + @abstractmethod + def items(self) -> Iterable[Tuple[IJKType, IndexLocation]]: """Return list of ((i, j, k), IndexLocation) tuples.""" - return self._locations.items() - - def backUp(self): - """Gather internal info that should be restored within a retainState.""" - self._backup = self._unitSteps, self._bounds, self._offset - - def restoreBackup(self): - self._unitSteps, self._bounds, self._offset = self._backup - - def getCoordinates(self, indices, nativeCoords=False) -> numpy.ndarray: - """Return the coordinates of the center of the mesh cell at the given given indices in cm.""" - indices = numpy.array(indices) - return self._evaluateMesh( - indices, self._centroidBySteps, self._centroidByBounds - ) - - def getCellBase(self, indices) -> numpy.ndarray: - """Get the mesh base (lower left) of this mesh cell in cm.""" - indices = numpy.array(indices) - return self._evaluateMesh( - indices, self._meshBaseBySteps, self._meshBaseByBounds - ) - - def getCellTop(self, indices) -> numpy.ndarray: - """Get the mesh top (upper right) of this mesh cell in cm.""" - indices = numpy.array(indices) + 1 - return self._evaluateMesh( - indices, self._meshBaseBySteps, self._meshBaseByBounds - ) + @abstractmethod def locatorInDomain( self, locator: LocationBase, symmetryOverlap: Optional[bool] = False ) -> bool: @@ -392,65 +136,15 @@ def locatorInDomain( Whether grid locations along the symmetry line should be considered "in the represented domain". This can be useful when assemblies are split along the domain boundary, with fractions of the assembly on either side. - """ - raise NotImplementedError("Not implemented on base Grid type.") - - def _evaluateMesh(self, indices, stepOperator, boundsOperator) -> numpy.ndarray: - """ - Evaluate some function of indices on this grid. - Recall from above that steps are mesh centered and bounds are mesh edged. - - Notes - ----- - This method may be able to be simplified. Complications from arbitrary - mixtures of bounds-based and step-based meshing caused it to get bad. - These were separate subclasses first, but in practice almost all cases have some mix - of step-based (hexagons, squares), and bounds based (radial, zeta). - - Improvements welcome! + Returns + ------- + bool + If the given locator is within the given grid """ - boundCoords = [] - for ii, bounds in enumerate(self._bounds): - if bounds is not None: - boundCoords.append(boundsOperator(indices[ii], bounds)) - - # limit step operator to the step dimensions - stepCoords = stepOperator(numpy.array(indices)[self._stepDims]) - - # now mix/match bounds coords with step coords appropriately. - result = numpy.zeros(len(indices)) - result[self._stepDims] = stepCoords - result[self._boundDims] = boundCoords - return result + self._offset - - def _centroidBySteps(self, indices): - return numpy.dot(self._unitSteps, indices) - - def _meshBaseBySteps(self, indices): - return ( - self._centroidBySteps(indices - 1) + self._centroidBySteps(indices) - ) / 2.0 - - @staticmethod - def _centroidByBounds(index, bounds): - if index < 0: - # avoid wrap-around - raise IndexError("Bounds-defined indices may not be negative.") - return (bounds[index + 1] + bounds[index]) / 2.0 - - @staticmethod - def _meshBaseByBounds(index, bounds): - if index < 0: - raise IndexError("Bounds-defined indices may not be negative.") - return bounds[index] - @staticmethod - def getNeighboringCellIndices(i, j=0, k=0): - """Return the indices of the immediate neighbors of a mesh point in the plane.""" - return ((i + 1, j, k), (1, j + 1, k), (i - 1, j, k), (i, j - 1, k)) - - def getSymmetricEquivalents(self, indices): + @abstractmethod + def getSymmetricEquivalents(self, indices: IJType) -> List[IJType]: """ Return a list of grid indices that contain matching contents based on symmetry. @@ -460,165 +154,82 @@ def getSymmetricEquivalents(self, indices): other two thirds of the grid, unless it is the central location, in which case no indices will be returned. """ - raise NotImplementedError - @staticmethod - def getAboveAndBelowCellIndices(indices): - i, j, k = indices - return ((i, j, k + 1), (i, j, k - 1)) + @abstractmethod + def overlapsWhichSymmetryLine(self, indices: IJType) -> Optional[int]: + """Return lines of symmetry position at a given index can be found - def cellIndicesContainingPoint(self, x, y=0.0, z=0.0): - """Return the indices of a mesh cell that contains a point.""" - raise NotImplementedError + Parameters + ---------- + indices : tuple of [int, int] + Indices for the requested object - def overlapsWhichSymmetryLine(self, indices): - return None + Returns + ------- + None or int + None if not line of symmetry goes through the object at the + requested index. Otherwise, some grid constants like ``BOUNDARY_CENTER`` + will be returned. - def getLabel(self, indices): """ - Get a string label from a 0-based spatial locator. - Returns a string representing i, j, and k indices of the locator - """ - i, j = indices[:2] - label = f"{i:03d}-{j:03d}" - if len(indices) == 3: - label += f"-{indices[2]:03d}" - return label - - def getIndexBounds(self): - """ - Get min index and number of indices in this grid. + @abstractmethod + def getCoordinates( + self, + indices: Union[IJKType, List[IJKType]], + nativeCoords: bool = False, + ) -> numpy.ndarray: + pass - Step-defined grids would be infinite but for the step limits defined in the constructor. + @abstractmethod + def backUp(self): + """Subclasses should modify the internal backup variable""" - Notes - ----- - This produces output that is intended to be passed to a ``range`` statement. - """ - indexBounds = [] - for minMax, bounds in zip(self._unitStepLimits, self._bounds): - if bounds is None: - indexBounds.append(minMax) - else: - indexBounds.append((0, len(bounds))) - return tuple(indexBounds) - - def getBounds( - self, - ) -> Tuple[ - Optional[Sequence[float]], Optional[Sequence[float]], Optional[Sequence[float]] - ]: - """Return the grid bounds for each dimension, if present.""" - return self._bounds + @abstractmethod + def restoreBackup(self): + """Restore state from backup""" - def getLocatorFromRingAndPos(self, ring, pos, k=0): - """ - Return the location based on ring and position. + @abstractmethod + def getCellBase(self, indices: IJKType) -> numpy.ndarray: + """Return the lower left case of this cell in cm""" - Parameters - ---------- - ring : int - Ring number (1-based indexing) - pos : int - Position number (1-based indexing) - k : int, optional - Axial index (0-based indexing) - - See Also - -------- - getIndicesFromRingAndPos - This implements the transform into i, j indices based on ring and position. - """ - i, j = self.getIndicesFromRingAndPos(ring, pos) - return self[i, j, k] + @abstractmethod + def getCellTop(self, indices: IJKType) -> numpy.ndarray: + """Get the upper right of this cell in cm""" @staticmethod - def getIndicesFromRingAndPos(ring, pos): + def getLabel(indices): """ - Return i, j indices given ring and position. - - Note - ---- - This should be implemented as a staticmethod, since no Grids currently in - exsistence actually need any instance data to perform this task, and - staticmethods provide the convenience of calling the method without an instance - of the class in the first place. - """ - raise NotImplementedError("Base Grid does not know about ring/pos") + Get a string label from a 0-based spatial locator. - def getMinimumRings(self, n: int) -> int: + Returns a string representing i, j, and k indices of the locator """ - Return the minimum number of rings needed to fit ``n`` objects. + i, j = indices[:2] + label = f"{i:03d}-{j:03d}" + if len(indices) == 3: + label += f"-{indices[2]:03d}" + return label - Warning - ------- - While this is useful and safe for answering the question of "how many rings do I - need to hold N things?", is generally not safe to use it to answer "I have N - things; within how many rings are they distributed?". This function provides a - lower bound, assuming that objects are densely-packed. If they are not actually - densely packed, this may be unphysical. + @abstractmethod + def reduce(self) -> Tuple[Hashable, ...]: """ - raise NotImplementedError("Base grid does not know about rings") - - def getPositionsInRing(self, ring: int) -> int: - """Return the number of positions within a ring.""" - raise NotImplementedError("Base grid does not know about rings") + Return the set of arguments used to create this Grid. - def getRingPos(self, indices) -> Tuple[int, int]: - """ - Get ring and position number in this grid. + This is very much like the argument tuple from ``__reduce__``, but we do not + implement ``__reduce__`` for real, because we are generally happy with + ``__getstate__`` and ``__setstate__`` for pickling purposes. However, getting + these arguments to ``__init__`` is useful for storing Grids to the database, as + they are more stable (less likely to change) than the actual internal state of + the objects. - For non-hex grids this is just i and j. + The return value should be hashable, such that a set of these can be created. - A tuple is returned so that it is easy to compare pairs of indices. - """ - # Regular grids dont really know about ring and position. We can try to see if - # their parent does! - if ( - self.armiObject is not None - and self.armiObject.parent is not None - and self.armiObject.parent.spatialGrid is not None - ): - return self.armiObject.parent.spatialGrid.getRingPos(indices) - - # For compatibility's sake, return __something__. TODO: We may want to just - # throw here, to be honest. - return tuple(i + 1 for i in indices[:2]) - - def getAllIndices(self): - """Get all possible indices in this grid.""" - iBounds, jBounds, kBounds = self.getIndexBounds() - allIndices = tuple( - itertools.product(range(*iBounds), range(*jBounds), range(*kBounds)) - ) - return allIndices - - def buildLocations(self): - """Populate all grid cells with a characteristic SpatialLocator.""" - for i, j, k in self.getAllIndices(): - loc = IndexLocation(i, j, k, self) - self._locations[(i, j, k)] = loc + The return type should be symmetric such that a similar grid can be + created just with the outputs of ``Grid.reduce``, e.g., + ``type(grid)(*grid.reduce())`` - @property - def pitch(self): - """ - The pitch of the grid. + Notes + ----- + For consistency, the second to last argument **must** be the geomType - Assumes 2-d unit-step defined (works for cartesian) """ - # TODO move this to the CartesianGrid - pitch = (self._unitSteps[0][0], self._unitSteps[1][1]) - if pitch[0] == 0: - raise ValueError(f"Grid {self} does not have a defined pitch.") - return pitch - - -def _tuplify(maybeArray): - if isinstance(maybeArray, (numpy.ndarray, list, tuple)): - maybeArray = tuple( - tuple(row) if isinstance(row, (numpy.ndarray, list)) else row - for row in maybeArray - ) - - return maybeArray diff --git a/armi/reactor/grids/hexagonal.py b/armi/reactor/grids/hexagonal.py index ab7cfd631..8a06d5d7c 100644 --- a/armi/reactor/grids/hexagonal.py +++ b/armi/reactor/grids/hexagonal.py @@ -1,23 +1,10 @@ -# Copyright 2023 TerraPower, LLC -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -import math -from typing import Optional +from math import sqrt +from typing import Tuple, List, Optional import numpy -from armi.utils import hexagon from armi.reactor import geometry +from armi.utils import hexagon from .constants import ( BOUNDARY_0_DEGREES, @@ -25,10 +12,10 @@ BOUNDARY_60_DEGREES, BOUNDARY_CENTER, ) -from .locations import IndexLocation -from .grid import Grid +from .locations import IndexLocation, IJKType, IJType +from .structuredgrid import StructuredGrid -COS30 = math.sqrt(3) / 2.0 +COS30 = sqrt(3) / 2.0 SIN30 = 1.0 / 2.0 # going CCW from "position 1" (top right) TRIANGLES_IN_HEXAGON = numpy.array( @@ -43,10 +30,15 @@ ) -class HexGrid(Grid): +class HexGrid(StructuredGrid): """ Has 6 neighbors in plane. + .. note:: + + It is recommended to use :meth:`fromPitch` rather than + calling the ``__init__`` constructor directly. + Notes ----- In an axial plane (i, j) are as follows (second one is pointedEndUp):: @@ -97,7 +89,7 @@ def fromPitch(pitch, numRings=25, armiObject=None, pointedEndUp=False, symmetry= HexGrid A functional hexagonal grid object. """ - side = pitch / math.sqrt(3.0) + side = hexagon.side(pitch) if pointedEndUp: # rotated 30 degrees CCW from normal # increases in i move you in x and y @@ -120,7 +112,7 @@ def fromPitch(pitch, numRings=25, armiObject=None, pointedEndUp=False, symmetry= ) @property - def pitch(self): + def pitch(self) -> float: """ Get the hex-pitch of a regular hexagonal array. @@ -131,7 +123,7 @@ def pitch(self): return self._unitSteps[1][1] @staticmethod - def indicesToRingPos(i, j): + def indicesToRingPos(i: int, j: int) -> Tuple[int, int]: """ Convert spatialLocator indices to ring/position. @@ -169,7 +161,7 @@ def indicesToRingPos(i, j): positionBase = 1 + edge * (ring - 1) return ring, positionBase + offset - def getMinimumRings(self, n): + def getMinimumRings(self, n: int) -> int: """ Return the minimum number of rings needed to fit ``n`` objects. @@ -180,11 +172,15 @@ def getMinimumRings(self, n): """ return hexagon.numRingsToHoldNumCells(n) - def getPositionsInRing(self, ring): - """Return the number of positions within a ring.""" + def getPositionsInRing(self, ring: int) -> int: + """ + Return the number of positions within a ring. + """ return hexagon.numPositionsInRing(ring) - def getNeighboringCellIndices(self, i, j=0, k=0): + def getNeighboringCellIndices( + self, i: int, j: int = 0, k: int = 0 + ) -> List[IJKType]: """ Return the indices of the immediate neighbors of a mesh point in the plane. @@ -210,9 +206,9 @@ def getLabel(self, indices): """ ring, pos = self.getRingPos(indices) if len(indices) == 2: - return Grid.getLabel(self, (ring, pos)) + return super().getLabel((ring, pos)) else: - return Grid.getLabel(self, (ring, pos, indices[2])) + return super().getLabel((ring, pos, indices[2])) @staticmethod def _indicesAndEdgeFromRingAndPos(ring, position): @@ -255,11 +251,11 @@ def _indicesAndEdgeFromRingAndPos(ring, position): return i, j, edge @staticmethod - def getIndicesFromRingAndPos(ring, pos): + def getIndicesFromRingAndPos(ring: int, pos: int) -> IJType: i, j, _edge = HexGrid._indicesAndEdgeFromRingAndPos(ring, pos) return i, j - def getRingPos(self, indices): + def getRingPos(self, indices: IJKType) -> Tuple[int, int]: """ Get 1-based ring and position from normal indices. @@ -268,17 +264,27 @@ def getRingPos(self, indices): getIndicesFromRingAndPos : does the reverse """ i, j = indices[:2] - return HexGrid.indicesToRingPos(i, j) + return self.indicesToRingPos(i, j) - def overlapsWhichSymmetryLine(self, indices): + def overlapsWhichSymmetryLine(self, indices: IJType) -> Optional[int]: """Return a list of which lines of symmetry this is on. - If none, returns [] - If on a line of symmetry in 1/6 geometry, returns a list containing a 6. - If on a line of symmetry in 1/3 geometry, returns a list containing a 3. - Only the 1/3 core view geometry is actually coded in here right now. + Parameters + ---------- + indices : tuple of [int, int] + Indices for the requested object + + Returns + ------- + None or int + None if not line of symmetry goes through the object at the + requested index. Otherwise, some grid constants like ``BOUNDARY_CENTER`` + will be returned. - Being "on" a symmetry line means the line goes through the middle of you. + Notes + ----- + - Only the 1/3 core view geometry is actually coded in here right now. + - Being "on" a symmetry line means the line goes through the middle of you. """ i, j = indices[:2] @@ -299,12 +305,12 @@ def overlapsWhichSymmetryLine(self, indices): return symmetryLine - def getSymmetricEquivalents(self, indices): + def getSymmetricEquivalents(self, indices: IJKType) -> List[IJType]: if ( self.symmetry.domain == geometry.DomainType.THIRD_CORE and self.symmetry.boundary == geometry.BoundaryType.PERIODIC ): - return HexGrid._getSymmetricIdenticalsThird(indices) + return self._getSymmetricIdenticalsThird(indices) elif self.symmetry.domain == geometry.DomainType.FULL_CORE: return [] else: @@ -315,7 +321,7 @@ def getSymmetricEquivalents(self, indices): ) @staticmethod - def _getSymmetricIdenticalsThird(indices): + def _getSymmetricIdenticalsThird(indices) -> List[IJType]: """This works by rotating the indices by 120 degrees twice, counterclockwise.""" i, j = indices[:2] if i == 0 and j == 0: @@ -323,7 +329,7 @@ def _getSymmetricIdenticalsThird(indices): identicals = [(-i - j, i), (j, -i - j)] return identicals - def triangleCoords(self, indices): + def triangleCoords(self, indices: IJKType) -> numpy.ndarray: """ Return 6 coordinate pairs representing the centers of the 6 triangles in a hexagon centered here. @@ -333,14 +339,14 @@ def triangleCoords(self, indices): scale = self.pitch / 3.0 return xy + scale * TRIANGLES_IN_HEXAGON - def changePitch(self, newPitchCm): + def changePitch(self, newPitchCm: float): """Change the hex pitch.""" - side = newPitchCm / math.sqrt(3.0) + side = hexagon.side(newPitchCm) self._unitSteps = numpy.array( ((1.5 * side, 0.0, 0.0), (newPitchCm / 2.0, newPitchCm, 0.0), (0, 0, 0)) )[self._stepDims] - def locatorInDomain(self, locator, symmetryOverlap: Optional[bool] = False): + def locatorInDomain(self, locator, symmetryOverlap: Optional[bool] = False) -> bool: # This will include the "top" 120-degree symmetry lines. This is to support # adding of edge assemblies. if self.symmetry.domain == geometry.DomainType.THIRD_CORE: @@ -348,7 +354,7 @@ def locatorInDomain(self, locator, symmetryOverlap: Optional[bool] = False): else: return True - def isInFirstThird(self, locator, includeTopEdge=False): + def isInFirstThird(self, locator, includeTopEdge=False) -> bool: """True if locator is in first third of hex grid.""" ring, pos = self.getRingPos(locator.indices) if ring == 1: @@ -367,7 +373,7 @@ def isInFirstThird(self, locator, includeTopEdge=False): return True return False - def generateSortedHexLocationList(self, nLocs): + def generateSortedHexLocationList(self, nLocs: int): """ Generate a list IndexLocations, sorted based on their distance from the center. @@ -380,27 +386,27 @@ def generateSortedHexLocationList(self, nLocs): nLocs = int(nLocs) # need to make this an integer # next, generate a list of locations and corresponding distances - locs = [] + locList = [] for ring in range(1, hexagon.numRingsToHoldNumCells(nLocs) + 1): positions = self.getPositionsInRing(ring) for position in range(1, positions + 1): i, j = self.getIndicesFromRingAndPos(ring, position) - locs.append(self[(i, j, 0)]) + locList.append(self[(i, j, 0)]) # round to avoid differences due to floating point math - locs.sort( + locList.sort( key=lambda loc: ( round(numpy.linalg.norm(loc.getGlobalCoordinates()), 6), loc.i, # loc.i=ring loc.j, ) ) # loc.j= pos - return locs[:nLocs] + return locList[:nLocs] # TODO: this is only used by testing and another method that just needs the count of assemblies # in a ring, not the actual positions def allPositionsInThird(self, ring, includeEdgeAssems=False): """ - Returns a list of all the positions in a ring (in the first third). + Returns a list of all the positions in a ring (in the first third) Parameters ---------- diff --git a/armi/reactor/grids/structuredgrid.py b/armi/reactor/grids/structuredgrid.py new file mode 100644 index 000000000..fa1f32fc6 --- /dev/null +++ b/armi/reactor/grids/structuredgrid.py @@ -0,0 +1,514 @@ +# Copyright 2023 Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from abc import abstractmethod +import collections +import itertools +from typing import Tuple, Union, List, Iterable, Optional, Sequence + +import numpy + +from .locations import IJKType, LocationBase, IndexLocation, MultiIndexLocation +from .grid import Grid + +# data structure for database-serialization of grids +GridParameters = collections.namedtuple( + "GridParameters", + ("unitSteps", "bounds", "unitStepLimits", "offset", "geomType", "symmetry"), +) + + +class StructuredGrid(Grid): + """ + A connected set of cells characterized by indices mapping to space and vice versa. + + The cells may be characterized by any mixture of regular repeating steps and + user-defined steps in any direction. + + For example, a 2-D hex lattice has constant, regular steps whereas a 3-D hex mesh + may have user-defined axial meshes. Similar for Cartesian, RZT, etc. + + Parameters + ---------- + unitSteps : tuple of tuples, optional + Describes the grid spatially as a function on indices. + Each tuple describes how each ``(x,y,or z)`` dimension is influenced by + ``(i,j,k)``. In other words, it is:: + + (dxi, dxj, jxk), (dyi, dyj, dyk), (dzi, dzj, dzk) + + where ``dmn`` is the distance (in cm) that dimension ``m`` will change as a + function of index ``n``. + + Unit steps are used as a generic method for defining repetitive grids in a + variety of geometries, including hexagonal and Cartesian. The tuples are not + vectors in the direction of the translation, but rather grouped by direction. If + the bounds argument is described for a direction, the bounds will be used rather + than the unit step information. The default of (0, 0, 0) makes all dimensions + insensitive to indices since the coordinates are calculated by the dot product + of this and the indices. With this default, any dimension that is desired to + change with indices should be defined with bounds. RZtheta grids are created + exclusively with bounds. + bounds : 3-tuple + Absolute increasing bounds in cm including endpoints of a non-uniform grid. + Each item represents the boundaries in the associated direction. Use Nones when + unitSteps should be applied instead. Most useful for thetaRZ grids or other + non-uniform grids. + unitStepLimits : 3-tuple + The limit of the steps in all three directions. This constrains step-defined + grids to be finite so we can populate them with SpatialLocator objects. + offset : 3-tuple, optional + Offset in cm for each axis. By default the center of the (0,0,0)-th object is in + the center of the grid. Offsets can move it so that the (0,0,0)-th object can + be fully within a quadrant (i.e. in a Cartesian grid). + armiObject : ArmiObject, optional + The ArmiObject that this grid describes. For example if it's a 1-D assembly + grid, the armiObject is the assembly. Note that ``self.armiObject.spatialGrid`` + is ``self``. + + Examples + -------- + A 2D a rectangular grid with width (x) 2 and height (y) 3 would be:: + + >>> grid = Grid(unitSteps=((2, 0, 0), (0, 3, 0),(0, 0, 0))) + + A regular hex grid with pitch 1 is:: + + >>> grid = Grid(unitSteps= ((sqrt(3)/2, 0.0, 0.0), (0.5, 1.0, 0.0), (0, 0, 0)) + + .. note:: For this unit hex the magnitude of the vector constructed using the + 0th index of each tuple is 1.0. + + Notes + ----- + Each dimension must either be defined through unitSteps or bounds. + The combination of unitSteps with bounds was settled upon after some struggle to + have one unified definition of a grid (i.e. just bounds). A hexagonal grid is + somewhat challenging to represent with bounds because the axes are not orthogonal, + so a unit-direction vector plus bounds would be required. And then the bounds would + be wasted space because they can be derived simply by unit steps. Memory efficiency + is important in this object so the compact representation of + unitSteps-when-possible, bounds-otherwise was settled upon. + + Design considerations include: + + * unitSteps are more intuitive as operations starting from the center of a cell, + particularly with hexagons and rectangles. Otherwise the 0,0 position of a hexagon + in the center of 1/3-symmetric hexagon is at the phantom bottom left of the + hexagon. + + * Users generally prefer to input mesh bounds rather than centers (e.g. starting at + 0.5 instead of 0.0 in a unit mesh is weird). + + * If we store bounds, computing bounds is simple and computing centers takes ~2x the + effort. If we store centers, it's the opposite. + + * Regardless of how we store things, we'll need a Grid that has the lower-left + assembly fully inside the problem (i.e. for full core Cartesian) as well as + another one that has the lower-left assembly half-way or quarter-way sliced off + (for 1/2, 1/4, and 1/8 symmetries). The ``offset`` parameter handles this. + + * Looking up mesh boundaries (to define a mesh in another code) is generally more + common than looking up centers (for plotting or measuring distance). + + * A grid can be anchored to the object that it is in with a backreference. This + gives it the ability to traverse the composite tree and map local to global + locations without having to duplicate the composite pattern on grids. This remains + optional so grids can be used for non-reactor-package reasons. It may seem + slightly cleaner to set the armiObject to the parent's spatialLocator itself + but the major disadvantage of this is that when an object moves, the armiObject + would have to be updated. By anchoring directly to Composite objects, the parent + is always up to date no matter where or how things get moved. + + * Unit step calculations use dot products and must not be polluted by the bound + indices. Thus we reduce the size of the unitSteps tuple accordingly. + + .. impl:: ARMI supports a number of structured mesh options. + :id: IMPL_REACTOR_MESH_0 + :links: REQ_REACTOR_MESH + """ + + def __init__( + self, + unitSteps=(0, 0, 0), + bounds=(None, None, None), + unitStepLimits=((0, 1), (0, 1), (0, 1)), + offset=None, + geomType="", + symmetry="", + armiObject=None, + ): + super().__init__(geomType, symmetry, armiObject) + # these lists contain the indices representing which dimensions for which steps + # are used, or for which bounds are used. index 0 is x direction, etc. + self._boundDims = [] + self._stepDims = [] + for dimensionIndex, bound in enumerate(bounds): + if bound is None: + self._stepDims.append(dimensionIndex) + else: + self._boundDims.append(dimensionIndex) + + # numpy prefers tuples like this to do slicing on arrays + self._boundDims = (tuple(self._boundDims),) + self._stepDims = (tuple(self._stepDims),) + + unitSteps = _tuplify(unitSteps) + + self._bounds = bounds + self._unitStepLimits = _tuplify(unitStepLimits) + + # only represent unit steps in dimensions they're being used so as to not + # pollute the dot product. This may reduce the length of this from 3 to 2 or 1 + self._unitSteps = numpy.array(unitSteps)[self._stepDims] + self._offset = numpy.zeros(3) if offset is None else numpy.array(offset) + self._locations = {} + self._buildLocations() # locations are owned by a grid, so the grid builds them. + + (_ii, iLen), (_ji, jLen), (_ki, kLen) = self.getIndexBounds() + # True if only contains k-cells. + self._isAxialOnly = iLen == jLen == 1 and kLen > 1 + + def __len__(self) -> int: + return len(self._locations) + + @property + def isAxialOnly(self) -> bool: + return self._isAxialOnly + + def reduce(self) -> GridParameters: + """Recreate the parameter necessary to create this grid""" + offset = None if not self._offset.any() else tuple(self._offset) + + bounds = _tuplify(self._bounds) + + # recreate a constructor-friendly version of `_unitSteps` from live data (may have been reduced from + # length 3 to length 2 or 1 based on mixing the step-based definition and the bounds-based definition + # described in Design Considerations above.) + # We don't just save the original tuple passed in because that may miss transformations that + # occur between instantiation and reduction. + unitSteps = [] + compressedSteps = list(self._unitSteps[:]) + for i in range(3): + # Recall _stepDims are stored as a single-value tuple (for numpy indexing) + # So this just is grabbing the actual data. + if i in self._stepDims[0]: + unitSteps.append(compressedSteps.pop(0)) + else: + # Add dummy value which will never get used (it gets reduced away) + unitSteps.append(0) + unitSteps = _tuplify(unitSteps) + + return GridParameters( + unitSteps, + bounds, + self._unitStepLimits, + offset, + self._geomType, + self._symmetry, + ) + + @property + def offset(self) -> numpy.ndarray: + """Offset in cm for each axis""" + return self._offset + + @offset.setter + def offset(self, offset: numpy.ndarray): + self._offset = offset + + def __repr__(self) -> str: + msg = ( + ["<{} -- {}\nBounds:\n".format(self.__class__.__name__, id(self))] + + [" {}\n".format(b) for b in self._bounds] + + ["Steps:\n"] + + [" {}\n".format(b) for b in self._unitSteps] + + [ + "Anchor: {}\n".format(self.armiObject), + "Offset: {}\n".format(self._offset), + "Num Locations: {}>".format(len(self)), + ] + ) + return "".join(msg) + + def __getitem__(self, ijk: Union[IJKType, List[IJKType]]) -> LocationBase: + """ + Get a location by (i, j, k) indices. If it does not exist, create a new one and return it. + + Parameters + ---------- + ijk : tuple of indices or list of the same + If provided a tuple, an IndexLocation will be created (if necessary) and + returned. If provided a list, each element will create a new IndexLocation + (if necessary), and a MultiIndexLocation containing all of the passed + indices will be returned. + + Notes + ----- + The method is defaultdict-like, in that it will create a new location on the fly. However, + the class itself is not really a dictionary, it is just index-able. For example, there is no + desire to have a ``__setitem__`` method, because the only way to create a location is by + retrieving it or through ``buildLocations``. + """ + try: + return self._locations[ijk] + except (KeyError, TypeError): + pass + + if isinstance(ijk, tuple): + i, j, k = ijk + val = IndexLocation(i, j, k, self) + self._locations[ijk] = val + elif isinstance(ijk, list): + val = MultiIndexLocation(self) + locators = [self[idx] for idx in ijk] + val.extend(locators) + else: + raise TypeError( + "Unsupported index type `{}` for `{}`".format(type(ijk), ijk) + ) + return val + + def items(self) -> Iterable[Tuple[IJKType, IndexLocation]]: + return self._locations.items() + + def backUp(self): + """Gather internal info that should be restored within a retainState.""" + self._backup = self._unitSteps, self._bounds, self._offset + + def restoreBackup(self): + self._unitSteps, self._bounds, self._offset = self._backup + + def getCoordinates(self, indices, nativeCoords=False) -> numpy.ndarray: + """Return the coordinates of the center of the mesh cell at the given given indices in cm.""" + indices = numpy.array(indices) + return self._evaluateMesh( + indices, self._centroidBySteps, self._centroidByBounds + ) + + def getCellBase(self, indices) -> numpy.ndarray: + """Get the mesh base (lower left) of this mesh cell in cm""" + indices = numpy.array(indices) + return self._evaluateMesh( + indices, self._meshBaseBySteps, self._meshBaseByBounds + ) + + def getCellTop(self, indices) -> numpy.ndarray: + """Get the mesh top (upper right) of this mesh cell in cm""" + indices = numpy.array(indices) + 1 + return self._evaluateMesh( + indices, self._meshBaseBySteps, self._meshBaseByBounds + ) + + def _evaluateMesh(self, indices, stepOperator, boundsOperator) -> numpy.ndarray: + """ + Evaluate some function of indices on this grid. + + Recall from above that steps are mesh centered and bounds are mesh edged. + + Notes + ----- + This method may be able to be simplified. Complications from arbitrary + mixtures of bounds-based and step-based meshing caused it to get bad. + These were separate subclasses first, but in practice almost all cases have some mix + of step-based (hexagons, squares), and bounds based (radial, zeta). + + Improvements welcome! + """ + boundCoords = [] + for ii, bounds in enumerate(self._bounds): + if bounds is not None: + boundCoords.append(boundsOperator(indices[ii], bounds)) + + # limit step operator to the step dimensions + stepCoords = stepOperator(numpy.array(indices)[self._stepDims]) + + # now mix/match bounds coords with step coords appropriately. + result = numpy.zeros(len(indices)) + result[self._stepDims] = stepCoords + result[self._boundDims] = boundCoords + return result + self._offset + + def _centroidBySteps(self, indices): + return numpy.dot(self._unitSteps, indices) + + def _meshBaseBySteps(self, indices): + return ( + self._centroidBySteps(indices - 1) + self._centroidBySteps(indices) + ) / 2.0 + + @staticmethod + def _centroidByBounds(index, bounds): + if index < 0: + # avoid wrap-around + raise IndexError("Bounds-defined indices may not be negative.") + return (bounds[index + 1] + bounds[index]) / 2.0 + + @staticmethod + def _meshBaseByBounds(index, bounds): + if index < 0: + raise IndexError("Bounds-defined indices may not be negative.") + return bounds[index] + + @staticmethod + def getNeighboringCellIndices(i, j=0, k=0): + """Return the indices of the immediate neighbors of a mesh point in the plane.""" + return ((i + 1, j, k), (1, j + 1, k), (i - 1, j, k), (i, j - 1, k)) + + @staticmethod + def getAboveAndBelowCellIndices(indices): + i, j, k = indices + return ((i, j, k + 1), (i, j, k - 1)) + + def getIndexBounds(self): + """ + Get min index and number of indices in this grid. + + Step-defined grids would be infinite but for the step limits defined in the constructor. + + Notes + ----- + This produces output that is intended to be passed to a ``range`` statement. + """ + indexBounds = [] + for minMax, bounds in zip(self._unitStepLimits, self._bounds): + if bounds is None: + indexBounds.append(minMax) + else: + indexBounds.append((0, len(bounds))) + return tuple(indexBounds) + + def getBounds( + self, + ) -> Tuple[ + Optional[Sequence[float]], Optional[Sequence[float]], Optional[Sequence[float]] + ]: + """ + Return the grid bounds for each dimension, if present. + """ + return self._bounds + + def getLocatorFromRingAndPos(self, ring, pos, k=0): + """ + Return the location based on ring and position. + + Parameters + ---------- + ring : int + Ring number (1-based indexing) + pos : int + Position number (1-based indexing) + k : int, optional + Axial index (0-based indexing) + + See Also + -------- + getIndicesFromRingAndPos + This implements the transform into i, j indices based on ring and position. + """ + i, j = self.getIndicesFromRingAndPos(ring, pos) + return self[i, j, k] + + @staticmethod + @abstractmethod + def getIndicesFromRingAndPos(ring: int, pos: int): + """ + Return i, j indices given ring and position. + + Note + ---- + This should be implemented as a staticmethod, since no Grids currently in + exsistence actually need any instance data to perform this task, and + staticmethods provide the convenience of calling the method without an instance + of the class in the first place. + """ + + @abstractmethod + def getMinimumRings(self, n: int) -> int: + """ + Return the minimum number of rings needed to fit ``n`` objects. + + Warning + ------- + While this is useful and safe for answering the question of "how many rings do I + need to hold N things?", is generally not safe to use it to answer "I have N + things; within how many rings are they distributed?". This function provides a + lower bound, assuming that objects are densely-packed. If they are not actually + densely packed, this may be unphysical. + """ + + @abstractmethod + def getPositionsInRing(self, ring: int) -> int: + """Return the number of positions within a ring.""" + + def getRingPos(self, indices) -> Tuple[int, int]: + """ + Get ring and position number in this grid. + + For non-hex grids this is just i and j. + + A tuple is returned so that it is easy to compare pairs of indices. + """ + + # Regular grids dont really know about ring and position. We can try to see if + # their parent does! + if ( + self.armiObject is not None + and self.armiObject.parent is not None + and self.armiObject.parent.spatialGrid is not None + ): + return self.armiObject.parent.spatialGrid.getRingPos(indices) + + # For compatibility's sake, return __something__. + # TODO: We may want to just throw here, to be honest. + return tuple(i + 1 for i in indices[:2]) + + def getAllIndices(self): + """Get all possible indices in this grid.""" + iBounds, jBounds, kBounds = self.getIndexBounds() + allIndices = tuple( + itertools.product(range(*iBounds), range(*jBounds), range(*kBounds)) + ) + return allIndices + + def _buildLocations(self): + """Populate all grid cells with a characteristic SpatialLocator.""" + for i, j, k in self.getAllIndices(): + loc = IndexLocation(i, j, k, self) + self._locations[(i, j, k)] = loc + + @property + @abstractmethod + def pitch(self) -> Union[float, Tuple[float, float]]: + """Grid pitch + + Some implementations may rely on a single pitch, such + as axial or hexagonal grids. Cartesian grids may use + a single pitch between elements or separate pitches + for the x and y dimensions. + + Returns + ------- + float or tuple of (float, float) + Grid spacing in cm + + """ + + +def _tuplify(maybeArray) -> tuple: + if isinstance(maybeArray, (numpy.ndarray, list, tuple)): + maybeArray = tuple( + tuple(row) if isinstance(row, (numpy.ndarray, list)) else row + for row in maybeArray + ) + + return maybeArray diff --git a/armi/reactor/grids/tests/test_grids.py b/armi/reactor/grids/tests/test_grids.py index 45a42eb6f..095ecb604 100644 --- a/armi/reactor/grids/tests/test_grids.py +++ b/armi/reactor/grids/tests/test_grids.py @@ -1,4 +1,4 @@ -# Copyright 2019 TerraPower, LLC +# Copyright 2023 TerraPower, LLC # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -12,14 +12,15 @@ # See the License for the specific language governing permissions and # limitations under the License. -"""Tests for grids.""" +"""Tests for grids""" +# pylint: disable=missing-function-docstring,missing-class-docstring,abstract-method,protected-access,no-self-use,attribute-defined-outside-init from io import BytesIO import math import unittest import pickle -from numpy.testing import assert_allclose, assert_array_equal import numpy +from numpy.testing import assert_allclose, assert_array_equal from armi.reactor import geometry from armi.reactor import grids @@ -49,6 +50,25 @@ def __init__(self, parent=None): self.parent = parent +class MockStructuredGrid(grids.StructuredGrid): + """Need a concrete class to test a lot of inherited methods + + Abstract methods from the parent now raise ``NotImplementedError`` + """ + + +# De-abstract the mock structured grid to test some basic +# properties, but let the abstract methods error +def _throwsNotImplemented(*args, **kwargs): + raise NotImplementedError + + +for f in MockStructuredGrid.__abstractmethods__: + setattr(MockStructuredGrid, f, _throwsNotImplemented) + +MockStructuredGrid.__abstractmethods__ = () + + class TestSpatialLocator(unittest.TestCase): def test_add(self): loc1 = grids.IndexLocation(1, 2, 0, None) @@ -68,13 +88,12 @@ def test_recursion(self): block = MockArmiObject(assem) # build meshes just like how they're used on a regular system. - coreGrid = grids.CartesianGrid.fromRectangle( - 1.0, 1.0, armiObject=core - ) # 2-D grid + # 2-D grid + coreGrid = grids.CartesianGrid.fromRectangle(1.0, 1.0, armiObject=core) + # 1-D z-mesh - assemblyGrid = grids.Grid( - bounds=(None, None, numpy.arange(5)), armiObject=assem - ) + assemblyGrid = grids.AxialGrid.fromNCells(5, armiObject=assem) + # pins sit in this 2-D grid. blockGrid = grids.CartesianGrid.fromRectangle(0.1, 0.1, armiObject=block) @@ -115,9 +134,7 @@ def test_recursionPin(self): # 2-D grid coreGrid = grids.CartesianGrid.fromRectangle(1.0, 1.0, armiObject=core) # 1-D z-mesh - assemblyGrid = grids.Grid( - bounds=(None, None, numpy.arange(5)), armiObject=assem - ) + assemblyGrid = grids.AxialGrid.fromNCells(5, armiObject=assem) # pins sit in this 2-D grid. blockGrid = grids.CartesianGrid.fromRectangle(0.1, 0.1, armiObject=block) @@ -142,19 +159,25 @@ def test_basicPosition(self): Full core Cartesian meshes will want to be shifted to bottom left of 0th cell. """ - grid = grids.Grid(unitSteps=((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))) + grid = MockStructuredGrid( + unitSteps=((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0)) + ) assert_allclose(grid.getCoordinates((1, 1, 1)), (1, 1, 1)) assert_allclose(grid.getCoordinates((0, 0, 0)), (0.0, 0.0, 0.0)) assert_allclose(grid.getCoordinates((0, 0, -1)), (0, 0, -1)) assert_allclose(grid.getCoordinates((1, 0, 0)), (1, 0, 0)) def test_neighbors(self): - grid = grids.Grid(unitSteps=((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))) + grid = MockStructuredGrid( + unitSteps=((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0)) + ) neighbs = grid.getNeighboringCellIndices(0, 0, 0) self.assertEqual(len(neighbs), 4) def test_label(self): - grid = grids.Grid(unitSteps=((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))) + grid = MockStructuredGrid( + unitSteps=((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0)) + ) self.assertEqual(grid.getLabel((1, 1, 2)), "001-001-002") def test_isAxialOnly(self): @@ -190,7 +213,7 @@ def test_getitem(self): class TestHexGrid(unittest.TestCase): - """A set of tests for the Hexagonal Grid. + """A set of tests for the Hexagonal Grid .. test: Tests of the Hexagonal grid. :id: TEST_REACTOR_MESH_0 @@ -431,27 +454,33 @@ def test_indicesAndEdgeFromRingAndPos(self): class TestBoundsDefinedGrid(unittest.TestCase): def test_positions(self): - grid = grids.Grid(bounds=([0, 1, 2, 3, 4], [0, 10, 20, 50], [0, 20, 60, 90])) + grid = MockStructuredGrid( + bounds=([0, 1, 2, 3, 4], [0, 10, 20, 50], [0, 20, 60, 90]) + ) assert_allclose(grid.getCoordinates((1, 1, 1)), (1.5, 15.0, 40.0)) def test_base(self): - grid = grids.Grid(bounds=([0, 1, 2, 3, 4], [0, 10, 20, 50], [0, 20, 60, 90])) + grid = MockStructuredGrid( + bounds=([0, 1, 2, 3, 4], [0, 10, 20, 50], [0, 20, 60, 90]) + ) assert_allclose(grid.getCellBase((1, 1, 1)), (1.0, 10.0, 20.0)) def test_positionsMixedDefinition(self): - grid = grids.Grid( + grid = MockStructuredGrid( unitSteps=((1.0, 0.0), (0.0, 1.0)), bounds=(None, None, [0, 20, 60, 90]) ) assert_allclose(grid.getCoordinates((1, 1, 1)), (1, 1, 40.0)) def test_getIndexBounds(self): - grid = grids.Grid(bounds=([0, 1, 2, 3, 4], [0, 10, 20, 50], [0, 20, 60, 90])) + grid = MockStructuredGrid( + bounds=([0, 1, 2, 3, 4], [0, 10, 20, 50], [0, 20, 60, 90]) + ) boundsIJK = grid.getIndexBounds() self.assertEqual(boundsIJK, ((0, 5), (0, 4), (0, 4))) class TestThetaRZGrid(unittest.TestCase): - """A set of tests for the RZTheta Grid. + """A set of tests for the RZTheta Grid .. test: Tests of the RZTheta grid. :id: TEST_REACTOR_MESH_1 @@ -474,7 +503,7 @@ def test_positions(self): class TestCartesianGrid(unittest.TestCase): - """A set of tests for the Cartesian Grid. + """A set of tests for the Cartesian Grid .. test: Tests of the Cartesian grid. :id: TEST_REACTOR_MESH_2 @@ -677,3 +706,8 @@ def test_getLocations(self): self.assertEqual(x, 0.0) self.assertEqual(y, 0.0) self.assertEqual(z, count + 0.5) + + +if __name__ == "__main__": + # import sys;sys.argv = ["", "TestHexGrid.testPositions"] + unittest.main() diff --git a/armi/reactor/grids/thetarz.py b/armi/reactor/grids/thetarz.py index 9551f064b..e5f228dd1 100644 --- a/armi/reactor/grids/thetarz.py +++ b/armi/reactor/grids/thetarz.py @@ -1,33 +1,30 @@ -# Copyright 2023 TerraPower, LLC -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. import math -from typing import Optional +from typing import TYPE_CHECKING, Optional, NoReturn import numpy -from .grid import Grid +from .locations import IJType, IJKType +from .structuredgrid import StructuredGrid -TAU = math.pi * 2.0 +if TYPE_CHECKING: + # Avoid circular imports + from armi.reactor.composites import ArmiObject +TAU = math.tau -class ThetaRZGrid(Grid): + +class ThetaRZGrid(StructuredGrid): """ A grid characterized by azimuthal, radial, and zeta indices. The angular meshes are limited to 0 to 2pi radians. R and Zeta are as in other meshes. + .. note:: + + It is recommended to call :meth:`fromGeom` to construct, + rather than directly constructing with ``__init__`` + See Figure 2.2 in Derstine 1984, ANL. [DIF3D]_. .. impl:: ARMI supports an RZTheta mesh. @@ -35,8 +32,13 @@ class ThetaRZGrid(Grid): :links: REQ_REACTOR_MESH """ - @staticmethod - def fromGeom(geom, armiObject=None): + def getSymmetricEquivalents(self, indices: IJType) -> NoReturn: + raise NotImplementedError( + f"{self.__class__.__name__} does not support symmetric equivalents" + ) + + @classmethod + def fromGeom(cls, geom, armiObject: Optional["ArmiObject"] = None) -> "ThetaRZGrid": """ Build 2-D R-theta grid based on a Geometry object. @@ -77,11 +79,13 @@ def getRingPos(self, indices): return (indices[1] + 1, indices[0] + 1) @staticmethod - def getIndicesFromRingAndPos(ring, pos): + def getIndicesFromRingAndPos(ring: int, pos: int) -> IJType: return (pos - 1, ring - 1) - def getCoordinates(self, indices, nativeCoords=False): - meshCoords = theta, r, z = Grid.getCoordinates(self, indices) + def getCoordinates(self, indices, nativeCoords=False) -> numpy.ndarray: + meshCoords = theta, r, z = super().getCoordinates( + indices, nativeCoords=nativeCoords + ) if not 0 <= theta <= TAU: raise ValueError("Invalid theta value: {}. Check mesh.".format(theta)) if nativeCoords: @@ -91,7 +95,14 @@ def getCoordinates(self, indices, nativeCoords=False): # return x, y ,z return numpy.array((r * math.cos(theta), r * math.sin(theta), z)) - def indicesOfBounds(self, rad0, rad1, theta0, theta1, sigma=1e-4): + def indicesOfBounds( + self, + rad0: float, + rad1: float, + theta0: float, + theta1: float, + sigma: float = 1e-4, + ) -> IJKType: """ Return indices corresponding to upper and lower radial and theta bounds. @@ -118,9 +129,26 @@ def indicesOfBounds(self, rad0, rad1, theta0, theta1, sigma=1e-4): return (i, j, 0) - def locatorInDomain(self, locator, symmetryOverlap: Optional[bool] = False): + @staticmethod + def locatorInDomain(*args, **kwargs) -> bool: """ ThetaRZGrids do not check for bounds, though they could if that becomes a problem. """ return True + + @staticmethod + def getMinimumRings(n: int) -> NoReturn: + raise NotImplementedError + + @staticmethod + def getPositionsInRing(ring: int) -> NoReturn: + raise NotImplementedError + + @staticmethod + def overlapsWhichSymmetryLine(indices: IJType) -> None: + return None + + @staticmethod + def pitch() -> NoReturn: + raise NotImplementedError() From ffe446a5fa6b67649c9cb38fe1040071aba9e4f5 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Thu, 24 Aug 2023 08:01:50 -0700 Subject: [PATCH 13/16] Add license headers to some grid files --- armi/reactor/grids/grid.py | 13 +++++++++++++ armi/reactor/grids/hexagonal.py | 13 +++++++++++++ armi/reactor/grids/thetarz.py | 13 +++++++++++++ 3 files changed, 39 insertions(+) diff --git a/armi/reactor/grids/grid.py b/armi/reactor/grids/grid.py index 0e32ea52c..3146e122f 100644 --- a/armi/reactor/grids/grid.py +++ b/armi/reactor/grids/grid.py @@ -1,3 +1,16 @@ +# Copyright 2023 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. from abc import ABC, abstractmethod from typing import Union, Optional, Hashable, TYPE_CHECKING, Dict, Iterable, Tuple, List diff --git a/armi/reactor/grids/hexagonal.py b/armi/reactor/grids/hexagonal.py index 8a06d5d7c..1517b39dd 100644 --- a/armi/reactor/grids/hexagonal.py +++ b/armi/reactor/grids/hexagonal.py @@ -1,3 +1,16 @@ +# Copyright 2023 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. from math import sqrt from typing import Tuple, List, Optional diff --git a/armi/reactor/grids/thetarz.py b/armi/reactor/grids/thetarz.py index e5f228dd1..13e05ef39 100644 --- a/armi/reactor/grids/thetarz.py +++ b/armi/reactor/grids/thetarz.py @@ -1,3 +1,16 @@ +# Copyright 2023 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. import math from typing import TYPE_CHECKING, Optional, NoReturn From a01054e21385cabf62e56bb8cffd5582345c8373 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Thu, 24 Aug 2023 08:03:00 -0700 Subject: [PATCH 14/16] Remove note:: callout in grid docstrings --- armi/reactor/grids/cartesian.py | 6 ++---- armi/reactor/grids/hexagonal.py | 6 ++---- armi/reactor/grids/thetarz.py | 6 ++---- 3 files changed, 6 insertions(+), 12 deletions(-) diff --git a/armi/reactor/grids/cartesian.py b/armi/reactor/grids/cartesian.py index 5ce795bac..0ec47cfb3 100644 --- a/armi/reactor/grids/cartesian.py +++ b/armi/reactor/grids/cartesian.py @@ -13,10 +13,8 @@ class CartesianGrid(StructuredGrid): """ Grid class representing a conformal Cartesian mesh - .. note:: - - It is recommended to call :meth:`fromRectangle` to construct, - rather than directly constructing with ``__init__`` + It is recommended to call :meth:`fromRectangle` to construct, + rather than directly constructing with ``__init__`` Notes ----- diff --git a/armi/reactor/grids/hexagonal.py b/armi/reactor/grids/hexagonal.py index 1517b39dd..2d0391de5 100644 --- a/armi/reactor/grids/hexagonal.py +++ b/armi/reactor/grids/hexagonal.py @@ -47,10 +47,8 @@ class HexGrid(StructuredGrid): """ Has 6 neighbors in plane. - .. note:: - - It is recommended to use :meth:`fromPitch` rather than - calling the ``__init__`` constructor directly. + It is recommended to use :meth:`fromPitch` rather than + calling the ``__init__`` constructor directly. Notes ----- diff --git a/armi/reactor/grids/thetarz.py b/armi/reactor/grids/thetarz.py index 13e05ef39..53116e5db 100644 --- a/armi/reactor/grids/thetarz.py +++ b/armi/reactor/grids/thetarz.py @@ -33,10 +33,8 @@ class ThetaRZGrid(StructuredGrid): The angular meshes are limited to 0 to 2pi radians. R and Zeta are as in other meshes. - .. note:: - - It is recommended to call :meth:`fromGeom` to construct, - rather than directly constructing with ``__init__`` + It is recommended to call :meth:`fromGeom` to construct, + rather than directly constructing with ``__init__`` See Figure 2.2 in Derstine 1984, ANL. [DIF3D]_. From 728f4912ba3b1bcb481f15b2be59dbabcf43146a Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Thu, 24 Aug 2023 08:12:06 -0700 Subject: [PATCH 15/16] Add two more missing license files --- armi/reactor/grids/axial.py | 13 +++++++++++++ armi/reactor/grids/cartesian.py | 13 +++++++++++++ 2 files changed, 26 insertions(+) diff --git a/armi/reactor/grids/axial.py b/armi/reactor/grids/axial.py index 69b8d35de..1bfaa997a 100644 --- a/armi/reactor/grids/axial.py +++ b/armi/reactor/grids/axial.py @@ -1,3 +1,16 @@ +# Copyright 2023 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. from typing import List, Optional, TYPE_CHECKING, NoReturn import warnings diff --git a/armi/reactor/grids/cartesian.py b/armi/reactor/grids/cartesian.py index 0ec47cfb3..c36f8923f 100644 --- a/armi/reactor/grids/cartesian.py +++ b/armi/reactor/grids/cartesian.py @@ -1,3 +1,16 @@ +# Copyright 2023 TerraPower, LLC +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. import itertools from typing import Optional, NoReturn, Tuple From 80696833dbc3fc3038031472294c08205f8ba838 Mon Sep 17 00:00:00 2001 From: Drew Johnson Date: Thu, 24 Aug 2023 13:41:08 -0700 Subject: [PATCH 16/16] Remove an abstract Grid from test_grids.py --- armi/reactor/grids/tests/test_grids.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/armi/reactor/grids/tests/test_grids.py b/armi/reactor/grids/tests/test_grids.py index 0daee57ab..924deee3f 100644 --- a/armi/reactor/grids/tests/test_grids.py +++ b/armi/reactor/grids/tests/test_grids.py @@ -213,7 +213,9 @@ def test_getitem(self): def test_ringPosFromIndicesIncorrect(self): """Test the getRingPos fails if there is no armiObect or parent.""" - grid = grids.Grid(unitSteps=((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0))) + grid = MockStructuredGrid( + unitSteps=((1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0)) + ) grid.armiObject = None with self.assertRaises(ValueError):