Move cell instance counting to C++

src/cell.cpp

@@ -437,16 +437,16 @@ read_cells(pugi::xml_node *node)
   }
 
   // Populate the Universe vector and dictionary.
-  for (Cell *c : cells_c) {
-    int32_t uid = c->universe;
+  for (int i = 0; i < cells_c.size(); i++) {
+    int32_t uid = cells_c[i]->universe;
     auto it = universe_dict.find(uid);
     if (it == universe_dict.end()) {
       universes_c.push_back(new Universe());
       universes_c.back()->id = uid;
-      universes_c.back()->cells.push_back(c);
+      universes_c.back()->cells.push_back(i);
       universe_dict[uid] = universes_c.size() - 1;
     } else {
-      universes_c[it->second]->cells.push_back(c);
+      universes_c[it->second]->cells.push_back(i);
     }
   }
 }
@@ -470,6 +470,8 @@ extern "C" {
 
   void cell_set_universe(Cell *c, int32_t universe) {c->universe = universe;}
 
+  int32_t cell_n_instances(Cell *c) {return c->n_instances;}
+
   bool cell_simple(Cell *c) {return c->simple;}
 
   bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface)

src/cell.h

@@ -43,7 +43,7 @@ class Universe
 public:
   int32_t id;                  //! Unique ID
   int32_t type;
-  std::vector<Cell*> cells;    //! Cells within this universe
+  std::vector<int32_t> cells;  //! Cells within this universe
   double x0, y0, z0;           //! Translation coordinates.
 };
 
@@ -59,6 +59,7 @@ public:
   int type;                  //!< Material, universe, or lattice
   int32_t universe;          //!< Universe # this cell is in
   int32_t fill;              //!< Universe # filling this cell
+  int32_t n_instances{0};    //!< Number of instances of this cell
 
   //! Material within this cell.  May be multiple materials for distribcell.
   //! C_NONE signifies a universe.

src/geometry.F90

@@ -750,72 +750,6 @@ contains
 
   end function count_target
 
-!===============================================================================
-! COUNT_INSTANCE recursively totals the number of occurrences of all cells
-! beginning with the universe given.
-!===============================================================================
-
-  recursive subroutine count_instance(univ)
-
-    type(Universe), intent(in) :: univ  ! universe to search through
-
-    integer :: i        ! index over cells
-    integer :: j, k, m  ! indices in lattice
-    integer :: n        ! number of cells to search
-
-    n = size(univ % cells)
-
-    do i = 1, n
-      associate (c => cells(univ % cells(i)))
-        c % instances = c % instances + 1
-
-        ! ====================================================================
-        ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
-        if (c % type() == FILL_UNIVERSE) then
-
-          call count_instance(universes(c % fill))
-
-        ! ====================================================================
-        ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
-        elseif (c % type() == FILL_LATTICE) then
-
-          ! Set current lattice
-          associate (lat => lattices(c % fill) % obj)
-
-            select type (lat)
-            type is (RectLattice)
-
-              ! Loop over lattice coordinates
-              do j = 1, lat % n_cells(1)
-                do k = 1, lat % n_cells(2)
-                  do m = 1, lat % n_cells(3)
-                    call count_instance(universes(lat % get([j-1, k-1, m-1])+1))
-                  end do
-                end do
-              end do
-
-            type is (HexLattice)
-
-              ! Loop over lattice coordinates
-              do m = 1, lat % n_axial
-                do k = 1, 2*lat % n_rings - 1
-                  do j = 1, 2*lat % n_rings - 1
-                    if (lat % are_valid_indices([j, k, m])) then
-                      call count_instance(universes(lat % get([j-1, k-1, m-1]) &
-                                                    + 1))
-                    end if
-                  end do
-                end do
-              end do
-
-            end select
-          end associate
-        end if
-      end associate
-    end do
-
-  end subroutine count_instance
-
 !===============================================================================
 ! MAXIMUM_LEVELS determines the maximum number of nested coordinate levels in
 ! the geometry

src/geometry_aux.cpp

@@ -11,49 +11,6 @@
 
 namespace openmc {
 
-extern "C" int32_t
-find_root_universe()
-{
-  // Find all the universes listed as a cell fill.
-  std::unordered_set<int32_t> fill_univ_ids;
-  for (Cell *c : cells_c) {
-    fill_univ_ids.insert(c->fill);
-  }
-
-  // Find all the universes contained in a lattice.
-  for (Lattice *lat : lattices_c) {
-    for (auto it = lat->begin(); it != lat->end(); ++it) {
-      fill_univ_ids.insert(*it);
-    }
-    if (lat->outer != NO_OUTER_UNIVERSE) {
-      fill_univ_ids.insert(lat->outer);
-    }
-  }
-
-  // Figure out which universe is not in the set.  This is the root universe.
-  bool root_found {false};
-  int32_t root_univ;
-  for (int32_t i = 0; i < universes_c.size(); i++) {
-    auto search = fill_univ_ids.find(universes_c[i]->id);
-    if (search == fill_univ_ids.end()) {
-      if (root_found) {
-        fatal_error("Two or more universes are not used as fill universes, so "
-                    "it is not possible to distinguish which one is the root "
-                    "universe.");
-      } else {
-        root_found = true;
-        root_univ = i;
-      }
-    }
-  }
-  if (!root_found) fatal_error("Could not find a root universe.  Make sure "
-       "there are no circular dependencies in the geometry.");
-
-  return root_univ;
-}
-
-//==============================================================================
-
 extern "C" void
 adjust_indices_c()
 {
@@ -101,4 +58,70 @@ adjust_indices_c()
   }
 }
 
+//==============================================================================
+
+extern "C" int32_t
+find_root_universe()
+{
+  // Find all the universes listed as a cell fill.
+  std::unordered_set<int32_t> fill_univ_ids;
+  for (Cell *c : cells_c) {
+    fill_univ_ids.insert(c->fill);
+  }
+
+  // Find all the universes contained in a lattice.
+  for (Lattice *lat : lattices_c) {
+    for (auto it = lat->begin(); it != lat->end(); ++it) {
+      fill_univ_ids.insert(*it);
+    }
+    if (lat->outer != NO_OUTER_UNIVERSE) {
+      fill_univ_ids.insert(lat->outer);
+    }
+  }
+
+  // Figure out which universe is not in the set.  This is the root universe.
+  bool root_found {false};
+  int32_t root_univ;
+  for (int32_t i = 0; i < universes_c.size(); i++) {
+    auto search = fill_univ_ids.find(universes_c[i]->id);
+    if (search == fill_univ_ids.end()) {
+      if (root_found) {
+        fatal_error("Two or more universes are not used as fill universes, so "
+                    "it is not possible to distinguish which one is the root "
+                    "universe.");
+      } else {
+        root_found = true;
+        root_univ = i;
+      }
+    }
+  }
+  if (!root_found) fatal_error("Could not find a root universe.  Make sure "
+       "there are no circular dependencies in the geometry.");
+
+  return root_univ;
+}
+
+//==============================================================================
+
+extern "C" void
+count_instances(int32_t univ_indx)
+{
+  for (int32_t cell_indx : universes_c[univ_indx]->cells) {
+    Cell &c {*cells_c[cell_indx]};
+    ++c.n_instances;
+
+    if (c.type == FILL_UNIVERSE) {
+      // This cell contains another universe.  Recurse into that universe.
+      count_instances(c.fill-1);  // TODO: off-by-one
+
+    } else if (c.type == FILL_LATTICE) {
+      // This cell contains a lattice.  Recurse into the lattice universes.
+      Lattice &lat {*lattices_c[c.fill-1]}; // TODO: off-by-one
+      for (auto it = lat.begin(); it != lat.end(); ++it) {
+        count_instances(*it);
+      }
+    }
+  }
+}
+
 } // namespace openmc

src/geometry_header.F90

@@ -67,6 +67,14 @@ module geometry_header
       integer(C_INT32_T), intent(in), value :: universe
     end subroutine cell_set_universe_c
 
+    function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') &
+        result(n_instances)
+      import C_PTR, C_INT32_T
+      implicit none
+      type(C_PTR), intent(in), value :: cell_ptr
+      integer(C_INT32_T)             :: n_instances
+    end function cell_n_instances_c
+
     function cell_simple_c(cell_ptr) bind(C, name='cell_simple') result(simple)
       import C_PTR, C_BOOL
       implicit none
@@ -245,8 +253,6 @@ module geometry_header
     type(C_PTR) :: ptr
 
     integer :: fill                        ! universe # filling this cell
-    integer :: instances                   ! number of instances of this cell in
-                                           !  the geom
     integer, allocatable :: material(:)    ! Material within cell.  Multiple
                                            !  materials for distribcell
                                            !  instances.  0 signifies a universe
@@ -273,6 +279,7 @@ module geometry_header
     procedure :: set_type => cell_set_type
     procedure :: universe => cell_universe
     procedure :: set_universe => cell_set_universe
+    procedure :: n_instances => cell_n_instances
     procedure :: simple => cell_simple
     procedure :: distance => cell_distance
     procedure :: to_hdf5 => cell_to_hdf5
@@ -393,6 +400,12 @@ contains
     call cell_set_universe_c(this % ptr, universe)
   end subroutine cell_set_universe
 
+  function cell_n_instances(this) result(n_instances)
+    class(Cell), intent(in) :: this
+    integer(C_INT32_T)      :: n_instances
+    n_instances = cell_n_instances_c(this % ptr)
+  end function cell_n_instances
+
   function cell_simple(this) result(simple)
     class(Cell), intent(in) :: this
     logical(C_BOOL)         :: simple

src/input_xml.F90

@@ -11,8 +11,7 @@ module input_xml
   use distribution_univariate
   use endf,             only: reaction_name
   use error,            only: fatal_error, warning, write_message, openmc_err_msg
-  use geometry,         only: calc_offsets, maximum_levels, count_instance, &
-                              neighbor_lists
+  use geometry,         only: calc_offsets, maximum_levels, neighbor_lists
   use geometry_header
   use hdf5_interface
   use list_header,      only: ListChar, ListInt, ListReal
@@ -51,6 +50,11 @@ module input_xml
     subroutine adjust_indices_c() bind(C)
     end subroutine adjust_indices_c
 
+    subroutine count_instances_c(univ_indx) bind(C, name='count_instances')
+      import C_INT32_T
+      integer(C_INT32_T), intent(in), value :: univ_indx
+    end subroutine
+
     subroutine read_surfaces(node_ptr) bind(C)
       import C_PTR
       type(C_PTR) :: node_ptr
@@ -139,7 +143,7 @@ contains
 
     ! Perform some final operations to set up the geometry
     call adjust_indices()
-    call count_instance(universes(root_universe))
+    call count_instances_c(root_universe-1)
 
     ! After reading input and basic geometry setup is complete, build lists of
     ! neighboring cells for efficient tracking
@@ -1030,7 +1034,6 @@ contains
       c % ptr = cell_pointer_c(i - 1)
 
       ! Initialize distribcell instances and distribcell index
-      c % instances = 0
       c % distribcell_index = NONE
 
       ! Get pointer to i-th cell node
@@ -3861,19 +3864,19 @@ contains
     do i = 1, n_cells
       associate (c => cells(i))
         if (size(c % material) > 1) then
-          if (size(c % material) /= c % instances) then
+          if (size(c % material) /= c % n_instances()) then
             call fatal_error("Cell " // trim(to_str(c % id())) // " was &
                  &specified with " // trim(to_str(size(c % material))) &
-                 // " materials but has " // trim(to_str(c % instances)) &
+                 // " materials but has " // trim(to_str(c % n_instances())) &
                  // " distributed instances. The number of materials must &
                  &equal one or the number of instances.")
           end if
         end if
         if (size(c % sqrtkT) > 1) then
-          if (size(c % sqrtkT) /= c % instances) then
+          if (size(c % sqrtkT) /= c % n_instances()) then
             call fatal_error("Cell " // trim(to_str(c % id())) // " was &
                  &specified with " // trim(to_str(size(c % sqrtkT))) &
-                 // " temperatures but has " // trim(to_str(c % instances)) &
+                 // " temperatures but has " // trim(to_str(c % n_instances())) &
                  // " distributed instances. The number of temperatures must &
                  &equal one or the number of instances.")
           end if

src/tallies/tally_filter_distribcell.F90

@@ -102,7 +102,7 @@ contains
     val = cell_dict % get(id)
     if (val /= EMPTY) then
       this % cell = val
-      this % n_bins = cells(this % cell) % instances
+      this % n_bins = cells(this % cell) % n_instances()
     else
       call fatal_error("Could not find cell " // trim(to_str(id)) &
            &// " specified on tally filter.")