From 1b3bbba41e30dc92a58aafefc06e41a47409b5bf Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 15 Aug 2018 11:49:17 -0400 Subject: [PATCH] Start implementing C++ find_cell --- src/geometry.F90 | 62 ++++++++------------------------------- src/geometry.cpp | 59 +++++++++++++++++++++++++++++++++++-- src/geometry.h | 1 + src/input_xml.F90 | 5 ---- src/settings.h | 4 ++- src/simulation_header.F90 | 6 +--- 6 files changed, 74 insertions(+), 63 deletions(-) diff --git a/src/geometry.F90 b/src/geometry.F90 index 5d715e98b..4a033c323 100644 --- a/src/geometry.F90 +++ b/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)) diff --git a/src/geometry.cpp b/src/geometry.cpp index 7c0dd28bc..4c31fe85d 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -3,8 +3,10 @@ #include #include "cell.h" +#include "constants.h" #include "error.h" #include "particle.h" +#include "settings.h" //TODO: remove this include #include @@ -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 diff --git a/src/geometry.h b/src/geometry.h index 7d345b102..e02a17b98 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -8,6 +8,7 @@ namespace openmc { extern "C" int openmc_root_universe; +//TODO: free this memory extern std::vector overlap_check_count; } // namespace openmc diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e81bb92f7..e4302478f 100644 --- a/src/input_xml.F90 +++ b/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) diff --git a/src/settings.h b/src/settings.h index 2bdd3bb5b..ff6b1cf3b 100644 --- a/src/settings.h +++ b/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 diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 60bed6425..340569402 100644 --- a/src/simulation_header.F90 +++ b/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)