`StackLattice` implementation

[yardasol]

This PR adds a new type of Lattice, StackLattice, meant to simplify the creation of universes repeated in one dimenson.

Specifically this PR:

• In lattice.py:
  • Adds the StackLattice class
  • Adds StackLattice as one of the possible classes in Lattice.from_hdf5
• In geometry.py
  • Adds flow control for StackLattice in geometry.py::from_xml
• In include/openmc/lattice.h
  • Add class diagram for StackLattice
• In src/lattice.cpp
  • Implement StackLattice (WIP)
• In src/geometry.cpp
  • Added StackLattice case to the switch statement in distance_to_boundary
• Adds functions for writing booleans to HDF5 files:
  • write_bool in hdf5_interface.cpp and hdf5_interface.h
• In tests/
  • Added slat_u and slat_nu objects, associated lines in test functions to unit_tests/test_lattice.py
  • Added stack_lat regression test to regression_tests. The tests cover both uniform and non-uniform stack lattices

This PR also adds a section on Stack Lattices to the User's Guide, a section on StackLattice tiling to the Theory Manual, as well as a reference to StackLattice in the Python API docs.

I tried my best on the C++ side but I'm much less familiar with that area of OpenMC (and rustier on my C++), so I'm sure there are some things I did inefficiently.

[yardasol]

This PR isn't quite ready yet, but it's getting there. I still need to do a bunch of testing though. Unfortunately, I'm running into a compiler error locally:

[  0%] Building CXX object vendor/fmt/CMakeFiles/fmt.dir/src/format.cc.o
.
.
.
[100%] Linking CXX executable bin/openmc
/usr/bin/ld: lib/libopenmc.so: undefined reference to `openmc::H5TypeMap<float>::type_id'
collect2: error: ld returned 1 exit status
make[2]: *** [CMakeFiles/openmc.dir/build.make:113: bin/openmc] Error 1
make[1]: *** [CMakeFiles/Makefile2:235: CMakeFiles/openmc.dir/all] Error 2
make: *** [Makefile:149: all] Error 2

Looks like the CI is also running into this. I double checked and haven't modified any variable of that type or name, so I'm not sure how to interpret this. Maybe I need to define the StackLattice C++ class somewhere else?

My C++ experience -- especially with things like shared libraries -- is limited in comparison to Python. Any ideas? @pshriwise @paulromano

[yardasol]

Took another look at this today and realized that float is indeed not a defined type. There was one variable I added in lattice.h that used this type.

[yardasol]

I have started to write a regression test. Right now the test is failing due to hitting the maximum number of lost particles due to many of the them having a negative distance to the lattice boundary. I looked at my implementation of distance function but couldn't see anything that stood out. I'll take another look tomorrow when I'm fresh, but if anyone wants to check out just that function, for anything obvious that I'm missing I'd appreciate it.

UPDATE: This isn't an obvious bug, but I think what must be happening is that the lattice index never updates even when a particle crosses a lattice boundary. I'll see if I can figure out where this is happening.

[yardasol]

I've done some pretty in depth debugging and everything seems to be working fine until we actually cross a lattice boundary for the first time. Something is going wrong in the call to Geometry::exhaustive_find_cell in Geometry::cross_lattice.

On the first call to exhaustive_find_cell in cross_lattice, we should be able to find the cell in the next lattice element without stepping into the if-statement at line 350 in geometry.cpp, as a StackLattice does not have any lattice tile corners. However, we do indeed step into the if statement.

Digging deeper, it looks like we fail to find the cell in the next lattice element due to the repeated universe not containing that cell. This is because the particle's local coordinate in the z- direction is -0.75, which is indeed outside any cell in the lattice. For some reason the lattice crossing occurs at z=0.75, but it should occur at z=1.5 (the absolute position of the lattice boundary). So it's entirely possible there is an issue with the model I'm using for this regression test, but I'm not ruling out an implementation error somewhere.

[yardasol]

I'm working on a regression test for nonuniform stack lattices. It's more or less a copy of the one for uniform stack lattices, but with the nonuniform bits added in. Nothing else about the model has changed, but now I get a bunch of warning messages like the following:

 WARNING: Particle 1 is outside lattice 3 but the lattice has no defined outer
          universe.
 WARNING: Particle 1 is outside lattice 3 but the lattice has no defined outer
          universe.
 WARNING: Couldn't find particle after reflecting from surface 7.

I tried to plot the geometry to see if my model is bugged but I get a similar error. Adding in an outer universe, running openmc --plot works, but the plot it produces is a bit wonky:

I suspect there may be something wrong with how the lattice indices are being found. Need to investigate further.

UPDATE: I fixed the "particle outside lattice" error by adding an outer universe, but I'm still getting the "Couldn't find particle after reflecting" error. I'm also getting errors like

 WARNING: After particle 2 crossed surface 5 it could not be located in any cell
          and it did not leak.

I'll have to take another look at my model, but maybe something is happening internally? Any ideas @paulromano?

[yardasol]

layer_boundaries_ for uniform stack lattice with 2 elements of pitch = 1.5:

(gdb) p layer_boundaries_
$1 = std::vector of length 3, capacity 4 = {0, 1.5, 3}

layer_boundaries_ for nonunifrom stack lattice with 2 elements of 'pitch' [1.0, 2.0]:

(gdb) p layer_boundaries_
$1 = std::vector of length 3, capacity 4 = {0, 1, 2}

The uniform case is correct, but the nonuniform case should be {0, 1, 3}. The same happens in the Python API so I'll need to adjust the implementation for creating the _layer_boundaries array there as well.

Note that this is because my current implementation does not catch the case where a non-iterable pitch is set for a non-uniform stack lattice. I'm working on a solution to this.

UPDATE: Commit 2bff4cc fixes the implementation and now catches errors when setting pitch for a StackLattice object.

[yardasol]

I did some more poking around today and realized that the get_local_position function doesn't handle cases where the lattice index is invalid, which I speculate is a component of telling OpenMC that the particle is no longer in the lattice. Regardless, I added control flow to handle these cases in 6789365 and the particle not found error has disappeared.

[yardasol]

This is a pretty big PR so I think 2 or 3 reviewers would be ideal @paulromano @pshriwise @anyone else who wants to review

[paulromano]

Thanks @yardasol for this PR! I'm definitely planning on reviewing this and would appreciate getting @pshriwise's input too. I unfortunately have to plead with you to be patient. We're in the midst of preparing for our 4-day OpenMC course next week, hence the pileup of unreviewed PRs. Once things are all set for the course, we should be able to start handling a lot of the PRs, including this one.

[yardasol]

Hi @paulromano, I apologize for the frequent pinging. I don't expect anyone to review PRs immediately when I tag them, so I'm a bit embarrassed that I came off in that way! I'll try to be more suave about that. I pinged you and Patrick because this PR has been sitting in draft mode for weeks now and I just wanted to make sure you knew it existed :). Good luck with the course next week!

[paulromano]

@yardasol I apologize that I haven't gotten you a full review on this yet. This is quite an impressive PR! But, it's also a big one that has some design implications so I need a bit of time to work through it and also try out the feature as well so I can better understand it. I just wanted you to know I haven't forgotten it :smile: If you happen to be free tomorrow at 9am CT and are up for it, it would be great if you could join our OpenMC monthly meet-up and walk through an example with the stack lattice feature.

[yardasol]

Hi @paulromano. Yes this PR is quite large. I appreciate you taking the time to go through it. What fortuitous timing for the monthly meet-up; I'd be happy to swing by and go through a notebook example with everyone.

Update: So, I decided to go with bit more sophisticated of a model for the example I want to show (a toy BWR fuel rod). Turns out that this contribution isn't ready at all in the way I've implemented it. The toy BWR fuel rod model is giving me bunch of particle not found errors again. Here's my hypothesis: StackLattice tiles are not explicitly bounded in two out of three dimensions, so if a model uses StackLattice but doesn't assign a cell to every possible point in each tile, there will be errors. The reason my integration tests work is because I happened to stumble upon this exact corner case where the current implementation works.

I'm thinking about how to allow users to specify a bound on their lattice in the unbounded dimensions. Perhaps an attribute tile_boundary that contains a Surface object (in the nonuniform case it would be an array of Surface objects) defining the boundary of each tile in the uniform case? This way we could easily compute the distance to the surface (on both the C++ backend and the Python frontend) and then if that surface is crossed, we are no longer in the lattice.

Even though I realize this is now not working as intended, I think I'll still show up to the meeting tomorrow anyways so I can get ideas and feedback from folks.

[yardasol]

Thanks for proposing this feature @yardasol! I've had a chance to look into the implementation a little bit more and do feel that the right path forward is to try to extend our current RectLattice class to support a one-dimensional lattice (as we discussed on the last monthly meeting call). The only case that won't cover is the non-uniform case, but I think that is much less common and raises some thorny design issues. Let me know what you think and if there's any support I can provide in moving toward the 1D RectLattice capability.

Hi @paulromano, I generally agree with your conclusion that extending RectLattice to have true 1D lattices would be the best path forwards.

We'll still need to figure out how to implement the 1D tiling structure without needing to define every region in space off to infinity within the pitch-defined bounds of the tile. I'll try to think about some good ways to do this.