terrapower / armi

An open-source nuclear reactor analysis automation framework that helps design teams increase efficiency and quality
https://terrapower.github.io/armi/
Apache License 2.0
231 stars 89 forks source link

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

Closed drewj-usnctech closed 7 months ago

drewj-usnctech commented 1 year ago

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 ```yaml 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 ```python 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 https://github.com/terrapower/armi/blob/328458cd62c1d4c6459ebcee1f2f6af10b6fb859/armi/utils/tests/test_asciimaps.py#L112-L131

https://github.com/terrapower/armi/blob/328458cd62c1d4c6459ebcee1f2f6af10b6fb859/armi/utils/tests/test_asciimaps.py#L312-L318

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 commented 1 year ago

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 commented 1 year ago

@ntouran - please help! Grids are magical

drewj-usnctech commented 1 year ago

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

drewj-usnctech commented 1 year ago

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

drewj-usnctech commented 1 year ago

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 https://github.com/terrapower/armi/blob/56762e20541a38640146370b93fba4bd4307e1d8/armi/reactor/grids.py#L1647-L1652

HexGrid.getSymmetricEquivalents

https://github.com/terrapower/armi/blob/56762e20541a38640146370b93fba4bd4307e1d8/armi/reactor/grids.py#L1628-L1635

HexGrid.overlapsWhichSymmetryLine

https://github.com/terrapower/armi/blob/56762e20541a38640146370b93fba4bd4307e1d8/armi/reactor/grids.py#L1584-L1611

HexGrid.getIndicesFromRingAndPos

https://github.com/terrapower/armi/blob/56762e20541a38640146370b93fba4bd4307e1d8/armi/reactor/grids.py#L1529-L1571

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 commented 1 year ago

I think this is a pretty appealing idea.

drewj-usnctech commented 1 year ago

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 commented 11 months ago

@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 commented 11 months ago

@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 commented 11 months ago

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 commented 8 months ago

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:

https://github.com/terrapower/armi/blob/571e4f58f364229b754f11b56d10ef2c153b92c9/armi/utils/asciimaps.py#L611-L615

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 commented 8 months ago

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

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

PXL_20240214_183657398~2

@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 commented 8 months ago

@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:

https://github.com/terrapower/armi/blob/571e4f58f364229b754f11b56d10ef2c153b92c9/armi/reactor/grids/hexagonal.py#L379

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

hex_grid_flats_up_third_equivalents

And then this "corners up" diagram:

hex_grid_corners_up_third_equivalents

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 commented 8 months ago

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

hex_grids_flats_v_corners

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

john-science commented 8 months ago

@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:

hex_grids_flats_v_corners_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:

https://github.com/terrapower/armi/blob/571e4f58f364229b754f11b56d10ef2c153b92c9/armi/reactor/grids/hexagonal.py#L330

https://github.com/terrapower/armi/blob/571e4f58f364229b754f11b56d10ef2c153b92c9/armi/reactor/grids/hexagonal.py#L333

https://github.com/terrapower/armi/blob/571e4f58f364229b754f11b56d10ef2c153b92c9/armi/reactor/grids/hexagonal.py#L336

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

keckler commented 8 months ago

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 commented 8 months ago

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 commented 8 months ago

@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):

hex_grid_top_left_corner

(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:

https://github.com/terrapower/armi/blob/6163e942abe873129533b84a6f2efa632f260b47/armi/reactor/grids/locations.py#L289

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.