Start implementing C++ find_cell

src/geometry.F90

@@ -37,6 +37,15 @@ module geometry
       import Particle
       type(Particle), intent(in) :: p
     end subroutine check_cell_overlap
+
+    function find_cell_c(p, n_search_cells, search_cells) &
+         bind(C, name="find_cell") result(found)
+      import Particle, C_INT, C_BOOL
+      type(Particle),  intent(in)           :: p
+      integer(C_INT),  intent(in), value    :: n_search_cells
+      integer(C_INT),  intent(in), optional :: search_cells(n_search_cells)
+      logical(C_BOOL)                       :: found
+    end function find_cell_c
   end interface
 
 contains
@@ -60,64 +69,19 @@ contains
     type(Particle), intent(inout) :: p
     logical,        intent(inout) :: found
     integer,        optional      :: search_cells(:)
-    integer :: i                    ! index over cells
     integer :: j, k                 ! coordinate level index
     integer :: offset               ! instance # of a distributed cell
     integer :: distribcell_index
     integer :: i_xyz(3)             ! indices in lattice
-    integer :: n                    ! number of cells to search
     integer :: i_cell               ! index in cells array
-    integer :: i_universe           ! index in universes array
-    logical :: use_search_cells     ! use cells provided as argument
 
-    do j = p % n_coord + 1, MAX_COORD
-      call reset_coord(p % coord(j))
-    end do
-    j = p % n_coord
-
-    ! Determine universe (if not yet set, use root universe)
-    i_universe = p % coord(j) % universe
-    if (i_universe == C_NONE) then
-      p % coord(j) % universe = root_universe
-      i_universe = root_universe
-    end if
-
-    ! set size of list to search
     if (present(search_cells)) then
-      use_search_cells = .true.
-      n = size(search_cells)
+      found = find_cell_c(p, size(search_cells), search_cells-1)
     else
-      use_search_cells = .false.
-      n = size(universes(i_universe) % cells)
+      found = find_cell_c(p, 0)
     end if
-
-    found = .false.
-    CELL_LOOP: do i = 1, n
-      ! select cells based on whether we are searching a universe or a provided
-      ! list of cells (this would be for lists of neighbor cells)
-      if (use_search_cells) then
-        i_cell = search_cells(i)
-        ! check to make sure search cell is in same universe
-        if (cells(i_cell) % universe() /= i_universe) cycle
-      else
-        i_cell = universes(i_universe) % cells(i)
-      end if
-
-      ! Move on to the next cell if the particle is not inside this cell
-      if (cell_contains(cells(i_cell), p)) then
-        ! Set cell on this level
-        p % coord(j) % cell = i_cell
-
-        ! Show cell information on trace
-        if (verbosity >= 10 .or. trace) then
-          call write_message("    Entering cell " // trim(to_str(&
-               cells(i_cell) % id())))
-        end if
-
-        found = .true.
-        exit
-      end if
-    end do CELL_LOOP
+    j = p % n_coord
+    i_cell = p % coord(j) % cell
 
     if (found) then
       associate(c => cells(i_cell))

src/geometry.cpp

@@ -3,8 +3,10 @@
 #include <sstream>
 
 #include "cell.h"
+#include "constants.h"
 #include "error.h"
 #include "particle.h"
+#include "settings.h"
 
 //TODO: remove this include
 #include <iostream>
@@ -20,14 +22,13 @@ check_cell_overlap(Particle* p) {
 
   // loop through each coordinate level
   for (int j = 0; j < n_coord; j++) {
-    //p->n_coord = j + 1;
-    Universe& univ {*global_universes[p->coord[j].universe - 1]};
+    Universe& univ = *global_universes[p->coord[j].universe - 1];
     int n = univ.cells.size();
 
     // loop through each cell on this level
     for (int i = 0; i < n; i++) {
       int index_cell = univ.cells[i];
-      Cell& c {*global_cells[index_cell]};
+      Cell& c = *global_cells[index_cell];
 
       if (c.contains(p->coord[j].xyz, p->coord[j].uvw, p->surface)) {
         //TODO: off-by-one indexing
@@ -46,4 +47,56 @@ check_cell_overlap(Particle* p) {
   return false;
 }
 
+//==============================================================================
+
+extern "C" bool
+find_cell(Particle* p, int n_search_cells, int* search_cells) {
+  for (int i = p->n_coord; i < MAX_COORD; i++) {
+    p->coord[i].reset();
+  }
+
+  // Determine universe (if not yet set, use root universe)
+  int i_universe = p->coord[p->n_coord-1].universe;
+  if (i_universe == C_NONE) {
+    p->coord[p->n_coord-1].universe = openmc_root_universe;
+    i_universe = openmc_root_universe;
+  }
+  //TODO: off-by-one indexing
+  --i_universe;
+
+  // If not given a set of search cells, search all cells in the uninverse.
+  if (n_search_cells == 0) {
+    search_cells = global_universes[i_universe]->cells.data();
+    n_search_cells = global_universes[i_universe]->cells.size();
+  }
+
+  // Find which cell of this universe the particle is in.
+  bool found = false;
+  for (int i = 0; i < n_search_cells; i++) {
+    int32_t i_cell = search_cells[i];
+
+    // Make sure the search cell is in the same universe
+    //TODO: off-by-one indexing
+    if (global_cells[i_cell]->universe - 1 != i_universe) continue;
+
+    Position r {p->coord[p->n_coord-1].xyz};
+    Direction u {p->coord[p->n_coord-1].uvw};
+    int32_t surf = p->surface;
+    if (global_cells[i_cell]->contains(r, u, surf)) {
+      //TODO: off-by-one indexing
+      p->coord[p->n_coord-1].cell = i_cell + 1;
+
+      if (openmc_verbosity >= 10 || openmc_trace) {
+        std::stringstream msg;
+        msg << "    Entering cell " << global_cells[i_cell]->id;
+        write_message(msg, 1);
+      }
+      found = true;
+      break;
+    }
+  }
+
+  return found;
+}
+
 } // namespace openmc

src/geometry.h

@@ -8,6 +8,7 @@ namespace openmc {
 
 extern "C" int openmc_root_universe;
 
+//TODO: free this memory
 extern std::vector<int64_t> overlap_check_count;
 
 } // namespace openmc

src/input_xml.F90

@@ -1086,11 +1086,6 @@ contains
     ! Allocate cells array
     allocate(cells(n_cells))
 
-    if (check_overlaps) then
-      allocate(overlap_check_cnt(n_cells))
-      overlap_check_cnt = 0
-    end if
-
     n_universes = 0
     do i = 1, n_cells
       c => cells(i)

src/settings.h

@@ -18,8 +18,10 @@ namespace openmc {
 extern "C" bool openmc_check_overlaps;
 extern "C" bool openmc_particle_restart_run;
 extern "C" bool openmc_restart_run;
-extern "C" bool openmc_write_all_tracks;
+extern "C" bool openmc_trace;
 extern "C" int openmc_verbosity;
+extern "C" bool openmc_write_all_tracks;
+#pragma omp threadprivate(openmc_trace)
 
 // Defined in .cpp
 // TODO: Make strings instead of char* once Fortran is gone

src/simulation_header.F90

@@ -67,10 +67,7 @@ module simulation_header
 
   integer :: restart_batch
 
-  ! Flag for enabling cell overlap checking during transport
-  integer(8), allocatable  :: overlap_check_cnt(:)
-
-  logical :: trace
+  logical(C_BOOL), bind(C, name='openmc_trace') :: trace
 
 !$omp threadprivate(trace, thread_id, current_work)
 
@@ -90,7 +87,6 @@ contains
 !===============================================================================
 
   subroutine free_memory_simulation()
-    if (allocated(overlap_check_cnt)) deallocate(overlap_check_cnt)
     if (allocated(entropy_p)) deallocate(entropy_p)
     if (allocated(source_frac)) deallocate(source_frac)
     if (allocated(work_index)) deallocate(work_index)