mit-crpg/openmc — pull request #313 by smharper (2014-09-06 22:05:39+00:00)
Title: Hex lattice

This pull request gives OpenMC the ability to simulate geometries with hexagonal lattices. In addition to adding the new hexagonal lattice features, this PR makes significant changes to the existing code. Here are the top-level changes:
* Users can specify a new element, hex_lattice, in the geometry.xml file
* The documentation has been updated, and two new tests have been added.
* Lattices in the geometry.xml file no longer need the `type` element/attribute.
* `pitch` is now used rather than `width` when specifying lattices in the geometry.xml file.
* In order to handle the changes with `type` and `width`, I added deprecation warnings to input.F90. I know that we don't technically need them since we are still in 0.x.y releases, but I think it's more polite to our users. I also added a Python utility to automatically update geometry.xml files.

Low-level changes:
* The lattice code in geometry.F90 is now a little more modular. Various lattice calculations were moved into the functions `get_lat_trans`, `get_lat_indices`, and `is_valid_lat_index`.
* Lattices use a polymorphic data structure---a cool feature of the Fortran 2003 standard. There is an abstract `Lattice` class. `Lattice` objects can be passed around by most functions without worrying about whether they are actually instances of the `RectLattice` or `HexLattice` types. Unfortunately, it also means that the global `lattices` array needs to be accessed like `lattices(i) % obj` rather than just `lattices(i)`. Sometime in the future I would like to convert the lattice primitives like `get_lat_trans` into type-bound procedures.
* Removed `lattice_edge` calculations from `find_cell`. This had no effect on the test suite, and I have not seen any damage caused by it's removal in my larger test calculations.
* `distance_to_boundary` stores lattice distances and surface distances in different variables in order to address issue #258

Design decisions that may affect performance or readability:
* Hexagonal lattice cell indexing has an algorithm that uses the particle's distance from various cell centers. It could probably be replaced with an algorithm that uses division remainders, but all of my attempted remainder algorithms suffered from finite precision issues that caused lost particles.
* The hexagonal lattice portion of `distance_to_boundary` uses a lot of expensive calculations to measure distances relative to the neighbor cells rather than the current cell. This algorithm is also used due to it's superior finite precision performance.
* There is a check in `distance_to_boundary` for negative lattice distances. I needed this for debugging, but I am considering leaving it in so that no future developers accidentally make changes that allow negative distances.
* `cross_lattice` now uses `get_lat_trans` which is a little bit more expensive than the older method of subtracting cell widths from the particle position, but it's agnostic to the lattice type and provides better finite precision performance (I had a lot of finite precision issues, can you tell?).

Verification and validation:
* I tested a large 2d hex lattice model against a model of manually placed cylinders. I got a disagreement of 3 +/- 12 pcm.
* I tested a 3d hex lattice model against a model using a 2d hex lattice placed in a 3d rectangular lattice. I got a disagreement of 17 +/- 22 pcm.
* I used an earlier version of hex_lattice to model various cores in the LR-0 reactor. The calculations differed from MCNP5 calculations by roughly 30 +/ 10 pcm and agreed with the experiment within uncertainties.

I can provide my V&V models to anyone who wants them. I'll also try to rerun them after we've gotten through a lot of the review process to make sure no errors have sneaked in. Also, I should mention that none of the regression test results are different (despite the changes to the handling of rectangular lattices) with the exception of score_flux_yn. I haven't updated that test yet because I am investigating it to see if I can find the reason for the change.