Location of pins in tips-up hexagonal grid is off by 60 degrees

[drewj-usnctech]

I have a tips up grid with fuel pins and one non-fuel pin on a tips-up hexagonal grid. The lattice map is

  pins:
    geom: hex_corners_up
    symmetry: full
    lattice map: |
      - - O F F
       - F F F F
        F F F F F
         F F F F
          F F F

where O is the non-fuel and F is the fuel pin.

As I read it, the non-fuel pin should be in the top left of the hexagonal grid. But, the locator attached to that component has a getLocalCoordinates of [-2.0, 0, 0], which implies it is located at due-left, 60 degrees off of where it should be.

Blueprint

blocks:
  fuel: &block_fuel
    grid name: pins
    duct:
      shape: Hexagon
      material: HT9
      Tinput: 450.0
      Thot: 450.0
      ip: 16
      op: 16.5
    coolant:
      shape: DerivedShape
      material: Sodium
      Tinput: 600.0
      Thot: 600.0
    fuel:
      shape: Circle
      material: UZr
      Tinput: 900
      Thot: 900
      od: 1
      latticeIDs: [F]
    other:
      shape: Circle
      material: HT9
      od: 1
      Tinput: 900
      Thot: 900
      latticeIDs: [O]
assemblies:
  fuel:
    specifier: F
    height: [10]
    axial mesh points: [1]
    blocks: [*block_fuel]
    xs types: [A]
systems:
  core:
    grid name: core
    origin:
      x: 0.0
      y: 0.0
      z: 0.0
grids:
  core:
    geom: hex
    symmetry: full
    lattice map: |
      F
  pins:
    geom: hex_corners_up
    symmetry: full
    lattice map: |
      - - O F F
       - F F F F
        F F F F F
         F F F F
          F F F
nuclide flags:
  NA: &nuc_flags {burn: false, xs: true}
  U: *nuc_flags
  ZR: *nuc_flags
  MO: *nuc_flags
  MN: *nuc_flags
  SI: *nuc_flags
  CR: *nuc_flags
  V: *nuc_flags
  C: *nuc_flags
  FE: *nuc_flags
  W: *nuc_flags
  NI: *nuc_flags

Processing script

import armi

armi.configure()

o = armi.init(fName="sample.yaml")

b = o.r.core[0][0]

for child in b:
    locator = child.spatialLocator
    if locator is None:
        continue
    if len(locator) == 1:
        print(child, locator[0].getLocalCoordinates())
        break

The asciimap testing shows it's being read in properly

armi/armi/utils/tests/test_asciimaps.py

Lines 112 to 131 in 328458c

	HEX_FULL_MAP = """- - - - - - - - - 1 1 1 1 1 1 1 1 1 4
	- - - - - - - - 1 1 1 1 1 1 1 1 1 1 1
	- - - - - - - 1 8 1 1 1 1 1 1 1 1 1 1
	- - - - - - 1 1 1 1 1 1 1 1 1 1 1 1 1
	- - - - - 1 1 1 1 1 1 1 1 1 1 1 1 1 1
	- - - - 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
	- - - 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
	- - 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
	- 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
	7 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1
	1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1
	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
	1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
	1 1 1 1 1 1 1 1 1 1 1 1 1 1
	1 1 1 1 1 1 1 1 1 3 1 1 1
	1 1 1 1 1 1 1 1 1 1 1 1
	1 6 1 1 1 1 1 1 1 1 1
	1 1 1 1 1 1 1 1 1 1
	"""

armi/armi/utils/tests/test_asciimaps.py

Lines 312 to 318 in 328458c

	self.assertEqual(asciimap[0, 0], "0")
	self.assertEqual(asciimap[0, -1], "2")
	self.assertEqual(asciimap[0, -8], "6")
	self.assertEqual(asciimap[0, 9], "4")
	self.assertEqual(asciimap[-9, 0], "7")
	self.assertEqual(asciimap[-8, 7], "8")
	self.assertEqual(asciimap[6, -6], "3")

with asciimap[-9, 0] == "7" being the location due-left (or negative x axis). But this doesn't seem to be sent to the grid properly

[keckler]

I am not very deep on this part of ARMI, so I can't speculate much right now about what's going on here. But I just wanted to chime in to add that this is a feature that TP cares about, so if it is broken, it will need to be fixed.

[jakehader]

@ntouran - please help! Grids are magical

[drewj-usnctech]

Just wanted to bring this back up. We have some work that is negatively impacted by this

[drewj-usnctech]

I'm working on a fix. No ETA on when it'll be done and shippable

[drewj-usnctech]

Looking through the code, it seems the spatial locators attached to these components have correct (i, j) indices. But the problem falls through HexGrid.getRingPos as that seems (after minor investigation) to assume a flats-up orientation.

Looking through more of the HexGrid implementation, making HexGrid.getRingPos more aware of the orientation seems to be the tip of the iceberg. It seems like multiple methods make this assumption

HexGrid.changePitch

Converts back to a flats up orientation

armi/armi/reactor/grids.py

Lines 1647 to 1652 in 56762e2

	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]

HexGrid.getSymmetricEquivalents

armi/armi/reactor/grids.py

Lines 1628 to 1635 in 56762e2

	@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

HexGrid.overlapsWhichSymmetryLine

armi/armi/reactor/grids.py

Lines 1584 to 1611 in 56762e2

	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

HexGrid.getIndicesFromRingAndPos

armi/armi/reactor/grids.py

Lines 1529 to 1571 in 56762e2

	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

Proposal

Would it be acceptable to introduce subclasses for these grid orientations because their implementations are so different? I think there's only so much that looking at unit steps can get you. Something like

class HexGrid(...):
    @classmethod
    def fromPitch(..., pointedEndsUp: bool=False) -> HexGrid:
        if pointedEndsUp:
            return TipsUpHexGrid.fromPitch(...)
       return FlatsUpHexGrid.fromPitch(...)

this has a nice side benefit of allowing grid orientation to be known simply while also not breaking calls to isinstance(grid, HexGrid)

cc @keckler @john-science @jakehader @ntouran

[keckler]

I think this is a pretty appealing idea.

[drewj-usnctech]

Maybe there's a confusion in the hex grid docs that is making this harder to understand / resolve

hex tips up ascii map has

     I axis is pure horizontal here 
     J axis is 60 degrees up. (upper right corner) 

But the HexGrid doc string implies the i axis is 60 degrees counterclockwise from x axis (upper right corner) and j axis is 120 degrees counter clockwise from x-axis (upper left corner)

                ( 0, 1) ( 1, 0)

            (-1, 1) ( 0, 0) ( 1,-1)

                (-1, 0) ( 0,-1)

[drewj-usnctech]

@john-science @keckler do you have any insight into which of the two tips up basis vectors are correct or preferred? I am still committed to getting this resolved, however slowly I can manage. But I don't want to take a gamble and miss and need to re-do the effort.

My gut is to keep the HexGrid orientation consistent, where j axis points to the upper left, because the code base seems to use a tips up hex grid only for auto-creating pin lattice grids. So keeping TipsUpHexGrid (ring, pos) <-> (i, j) consistent might not be as big of a pain.

edit after some thought: The whole point of this ticket though is that hex grid (ring, pos) <-> (i, j) is broken for tips up grids though. So things are likely going to break (towards a more correct thing!) regardless 🙃

[keckler]

@drewj-usnctech Sorry for the slow reply.

Unfortunately I really don't have much insight on which way is preferred to resolve this. As you noted, it is actively broken, so I think the most important part is that it gets fixed. For one reason or another, nobody has come across this error yet, so I wouldn't expect your fix to cause any huge issues for our current workflows.

My expectation is that whatever basis direction you choose to fix this will be fine. As long as it is consistent throughout the framework and fixes the broken functionality, we will be happy.

[drewj-usnctech]

No worries @keckler from the guy who's also being slow to respond. but thanks for responding nonetheless!

I've decided to use a 120 degree angle between the ij basis vectors, with i being in the positive x axis, and j being 120 degrees rotated counter clockwise. This is what is alluded to in the grids package.

I'm going to try to align the ring,pos orientation as much as possible with the flats up hex grid. I think the flats up orientation (position 1 at 60 degrees counter clockwise from x-axis, increasing counter clockwise) is linked to DIF3D. I can't find anything to verify that orientation in what little DIF3D docs I can find on the web. Not sure if DIF3D would even be able to support a tips up hexagonal core so this might be less of an issue.

Anything you or @ntouran may be able to provide on the (ring, pos) layout for a tips up hex grid?

[john-science]

I am working through this Issue. It will seem random and unrelated, but it's not; I just found a bug in how hex corners-up grids calculate their grid offsets:

armi/armi/utils/asciimaps.py

Lines 611 to 615 in 571e4f5

	def _makeOffsets(self):
	"""Full hex tips-up grids have linearly incrementing offset."""
	AsciiMap._makeOffsets(self)
	for li, _line in enumerate(self.asciiLines):
	self.asciiOffsets.append(li)

Line 613 there should read:

     self.asciiOffsets = []

I will continue working through all the other issues and questions in the ticket. Just an update; we are back to full power again.

[john-science]

Progress Update: I fixed changePitch(). It was indeed broken.

Here is my little drawing that our original unitSteps were correct:

@drewj-usnctech I know you suggested using new subclasses of HexGrid. But that turns out to be problematic for people downstream who use blueprint-reading code to directly generate grids. It would just mean a lot of downstream changes. So, right now, I am just added as self.cornersUp attribute to HexGrid, that way the change can be 100% ARMI, and not mean changes downstream. (TBD, but that's how it looks right now.)

[john-science]

@drewj-usnctech I believe the HexGrid._getSymmetricIdenticalsThird() method is already correct, and does not need any changes. It currently provides the same solution for "flats up" and "corners up" grids:

armi/armi/reactor/grids/hexagonal.py

Line 379 in 571e4f5

identicals = [(-i - j, i), (j, -i - j)]

But I believe that is correct. See this "flats up" diagram:

And then this "corners up" diagram:

Just do the math (ignore the center hex): the Dark Green rotate to the Red, and the Light Green rotate to each other.

[john-science]

By the way, if anyone wants a good plot of the "flats up" vs "corners up" hex grids, here they are (thanks @ntouran ):

I may try to highlight this more in the docs, as part of this PR.

[john-science]

@ntouran @keckler I was looking at the HexGrid.overlapsWhichSymmetryLine() method, it does the same thing for "corners up" and "flats up" hex grids. Here is what it identifies as the symmetry lines:

What I can't tell is... is this right or wrong? I can't tell. The only reference I have is these comments in the code:

armi/armi/reactor/grids/hexagonal.py

Line 330 in 571e4f5

# edge 1: 1/3 symmetry line (bottom horizontal side in 1/3 core view, theta = 0)

armi/armi/reactor/grids/hexagonal.py

Line 333 in 571e4f5

# edge 2: 1/6 symmetry line (bisects 1/3 core view, theta = pi/3)

armi/armi/reactor/grids/hexagonal.py

Line 336 in 571e4f5

# edge 3: 1/3 symmetry line (left oblique side in 1/3 core view, theta = 2*pi/3)

On balance, I think this might be correct as it is now.

[keckler]

I don't immediately see anything wrong with it.

But it seems like it must only make sense with respect to some "definition" of how to define 1/3 of a hex grid. I could also imagine that someone would draw the 1/3 boundaries extending from the center to the corners of the grid.

I'm not the best person to answer this question. @mgjarrett might have more insight.

[mgjarrett]

I think there is nothing wrong with the current definition of the symmetry lines. As Chris said above, you could just define the hex grid with the symmetry lines rotated by 30 degrees, but this current configuration makes the most sense for our applications.

I don't think either way is inherently more "correct" than the other; the only thing that matters is that it is treated consistently throughout ARMI.

[john-science]

@drewj-usnctech In response to your original concern:

As I read it, the non-fuel pin should be in the top left of the hexagonal grid. But, the locator attached to that component has a getLocalCoordinates of [-2.0, 0, 0], which implies it is located at due-left, 60 degrees off of where it should be.

So, if you I made the following little change to the example you show (notice the 99 on the top-left):

HEX_FULL_MAP = """- - - - - - - - - 99 1 1 1 1 1 1 1 1 4 
  - - - - - - - - 1 1 1 1 1 1 1 1 1 1 1 
   - - - - - - - 1 8 1 1 1 1 1 1 1 1 1 1 
    - - - - - - 1 1 1 1 1 1 1 1 1 1 1 1 1 
     - - - - - 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
      - - - - 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
       - - - 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
        - - 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
         - 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
          7 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 
           1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 
            1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
             1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
              1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 
               1 1 1 1 1 1 1 1 1 1 1 1 1 1 
                1 1 1 1 1 1 1 1 1 3 1 1 1 
                 1 1 1 1 1 1 1 1 1 1 1 1 
                  1 6 1 1 1 1 1 1 1 1 1 
                   1 1 1 1 1 1 1 1 1 1 
 """ 

Then the unit test line you link would show:

self.assertEqual(asciimap[-9, 9], "99") 

And, so, that's at indices for (-9,9). But above I gave a plot of our (I, J) indexes):

(Okay, that plot only goes out to "6" and not "9", but the text was too small to read with a grid that went to "9". The idea is the same.)

So, perhaps the (I,J) mapping is 60 degrees off from where you thought it should be. This must be your primary complaint?

So, if I alter your script to have this pin grid:

  pins:
    geom: hex_corners_up
    symmetry: full
    lattice map: |
      - - - - - - - - - O F F F F F F F F F
       - - - - - - - - F F F F F F F F F F F
        - - - - - - - F F F F F F F F F F F F
         - - - - - - F F F F F F F F F F F F F
          - - - - - F F F F F F F F F F F F F F
           - - - - F F F F F F F F F F F F F F F
            - - - F F F F F F F F F F F F F F F F
             - - F F F F F F F F F F F F F F F F F
              - F F F F F F F F F F F F F F F F F F
               F F F F F F F F F F F F F F F F F F F
                F F F F F F F F F F F F F F F F F F
                 F F F F F F F F F F F F F F F F F
                  F F F F F F F F F F F F F F F F
                   F F F F F F F F F F F F F F F
                    F F F F F F F F F F F F F F
                     F F F F F F F F F F F F F
                      F F F F F F F F F F F F
                       F F F F F F F F F F F
                        F F F F F F F F F F

And I run your script, with an extra print, we get this:

print(child, locator[0].indices, locator[0].getLocalCoordinates())

# <Circle: other>    [-9  9  0]    [-9 0 0]

So, the locator[0].indices matches our (-9, 9) . The getLocalCoordinates() are mean to be:

armi/armi/reactor/grids/locations.py

Line 289 in 6163e94

"""Return the coordinates of the center of the mesh cell here in cm."""

So, for the moment, I'm going to leave the "position in centimeters" of the center of a hex grid cell alone. And we can focus on this 60-degree rotation of the hex grid you don't like. The most important thing is all these grids are self-consistent, and that appears to be true. But I guess your question is:

Can we rotate all the hex grids back to where they would match the graphics in the blueprint?

Or something sort of like that.