Lattice offset table should be allocated differently

[pshriwise]

Bug Description

The Lattice::allocate_offset_table method resizes a vector to the appropriate size for holding N mappings of cell instance offsets based on the Lattice's dimensions when determining distribcell offsets. Nothing goes wrong if the offsets are set once, but if this method is called a second time with a different set of distribcells (e.g. during a mapping/re-mapping phase in Cardinal) the way the offset table is reallocated/resized makes it possible for some previous entries to activate the memoization of offset counts for lattices, but with values that are stale and incorrect. The outcome is the possibility of incorrect cell instance values if distribcell setup occurs more than once in an application.

Additional Detail

As mentioned above, Lattice::allocate_offset_table method resizes a vector of offsets to the appropriate size for holding N mappings of cell instance offsets based on the Lattice's dimensions when determining distribcell offsets in openmc::prepare_distribcell. Any additional entries created in the vector are filled with C_NONE (as descried in the signature for std::vector<T,Allocator>::resize( size_type count, const value_type& value )). None of the existing entries are changed, however.

In a scenario where openmc::prepare_distribcell is called more than once with different distribcell sets, another cell may take up the map value (distribcell index) for the Lattice. During the population of the Lattice::offset_ table, there is a block of code that skips a potentially expensive recursive call into count_universe_instances.

openmc/src/lattice.cpp

Lines 114 to 120 in 339d78c

	if (offsets_[map * universes_.size() + this->begin().indx_] != C_NONE) {
	int last_offset =
	offsets_[(map + 1) * universes_.size() - this->begin().indx_ - 1];
	int last_univ = this->back();
	return last_offset +
	count_universe_instances(last_univ, target_univ_id, univ_count_memo);
	}

In the case that a new call has taken up the map of another cell in a previous prepare_distribcell invocation, this line will skip recalculation of the offset table entries for that map, leaving stale and incorrect values in place. As a result, cells will compute incorrect instance values.

Oddly enough, this wasn't a problem for the HexLattice until #2921. The value being checked in the memoization block above would always discover a C_NONE value by nature of checking one of the intentionally empty entries in the offset table and these values were always recomputed, resulting in the correct cell instance values.

Steps to Reproduce

Here is a relatively simple model file and a simple C++ program that queries the cell and instance at a particular location for Cell 9. The invalid offset table results in an instance value of 3 at this location when it should be 0 as the universe it's part of is only applied once in the lattice.

<?xml version='1.0' encoding='UTF-8'?>
<geometry>
  <cell id="1" material="1" universe="1"/>
  <cell id="2" material="5" region="2" universe="4"/>
  <cell fill="2" id="3" region="-1" universe="3"/>
  <cell id="4" material="5" region="1" universe="3"/>
  <cell id="5" material="7" region="-2" universe="4"/>
  <cell id="6" material="5" region="2" universe="6"/>
  <cell fill="2" id="7" region="-1" universe="5"/>
  <cell id="8" material="5" region="1" universe="5"/>
  <cell id="9" material="7" region="-2" universe="6"/>
  <cell id="10" material="5" region="2" universe="8"/>
  <cell fill="2" id="11" region="-1" universe="7"/>
  <cell id="12" material="5" region="1" universe="7"/>
  <cell id="13" material="7" region="-2" universe="8"/>
  <cell fill="9" id="14" region="-3 4 -5 -8 6 7 9 -12" universe="10"/>
  <hex_lattice id="2" n_rings="2">
    <pitch>1.8800000000000001</pitch>
    <center>0.0 0.0</center>
    <universes>
 1
1 1
 1
1 1
 1</universes>
  </hex_lattice>
  <hex_lattice id="9" n_axial="3" n_rings="2" name="Unit cell lattice" orientation="x">
    <pitch>1.8800000000000001 1.6666666666666667</pitch>
    <center>0.0 0.0 2.5</center>
    <universes>
 3 3
3 4 3
 3 3
 5 5
5 6 5
 5 5
 7 7
7 8 7
 7 7</universes>
  </hex_lattice>
  <surface coeffs="0.0 0.0 0.635" id="1" type="z-cylinder"/>
  <surface coeffs="0.0 0.0 0.8" id="2" type="z-cylinder"/>
  <surface boundary="periodic" coeffs="1.6281277591147447" id="3" periodic_surface_id="4" type="y-plane"/>
  <surface boundary="periodic" coeffs="-1.6281277591147447" id="4" periodic_surface_id="3" type="y-plane"/>
  <surface boundary="periodic" coeffs="1.7320508075688772 1.0 0.0 3.2562555182294894" id="5" periodic_surface_id="7" type="plane"/>
  <surface boundary="periodic" coeffs="-1.7320508075688772 1.0 0.0 -3.2562555182294894" id="6" periodic_surface_id="8" type="plane"/>
  <surface boundary="periodic" coeffs="1.7320508075688772 1.0 0.0 -3.2562555182294894" id="7" periodic_surface_id="5" type="plane"/>
  <surface boundary="periodic" coeffs="-1.7320508075688772 1.0 0.0 3.2562555182294894" id="8" periodic_surface_id="6" type="plane"/>
  <surface boundary="vacuum" coeffs="0.0" id="9" type="z-plane"/>
  <surface boundary="vacuum" coeffs="5.0" id="12" type="z-plane"/>
</geometry>

#include "openmc/capi.h"
#include "openmc/cell.h"
#include "openmc/geometry.h"
#include "openmc/geometry_aux.h"
#include "openmc/particle.h"

int main(int argc, char** argv) {

  openmc_init(argc, argv, nullptr);

  openmc::Particle p;
  p.r() = {-0.11266409, -0.0346915, 2.5};
  p.u() = {0, 0, 1};

  std::vector<int32_t> cells = {9, 12, 11, 10, 5, 8, 7, 6, 1, 4, 3, 2};
  openmc::prepare_distribcell(&cells);

  openmc::exhaustive_find_cell(p);

  int32_t level = 1;

  auto cell_index = p.coord(level).cell;

  auto cell_instance = openmc::cell_instance_at_level(p, level);

  std::cout << "Cell: " << openmc::model::cells[cell_index]->id_ << " (index " << cell_index << "), Instance: " << cell_instance << std::endl;

  openmc_finalize();

  return 0;
}

[pshriwise]

I should note that a simple fix for this is a one-liner to reset the values of the table when it's resized. I'll submit it in a bit.