From bc1c1e491628741fb8457d5a1ecc69d9c259b541 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 20 Aug 2018 23:42:51 -0400 Subject: [PATCH] Move neighbor lists to C++ --- src/cell.cpp | 18 -------- src/cell.h | 7 ++++ src/geometry.F90 | 92 +++++------------------------------------ src/geometry.cpp | 27 ++++++++---- src/geometry.h | 2 +- src/geometry_aux.cpp | 28 +++++++++++++ src/geometry_aux.h | 6 +++ src/geometry_header.F90 | 6 --- src/input_xml.F90 | 44 +------------------- src/string.F90 | 92 ----------------------------------------- src/surface.h | 5 ++- src/surface_header.F90 | 3 -- src/tracking.F90 | 17 +------- 13 files changed, 79 insertions(+), 268 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 9786a89fc6..c22813ea71 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -17,17 +17,6 @@ namespace openmc { -//============================================================================== -// Constants -//============================================================================== - -// TODO: Convert to enum -constexpr int32_t OP_LEFT_PAREN {std::numeric_limits::max()}; -constexpr int32_t OP_RIGHT_PAREN {std::numeric_limits::max() - 1}; -constexpr int32_t OP_COMPLEMENT {std::numeric_limits::max() - 2}; -constexpr int32_t OP_INTERSECTION {std::numeric_limits::max() - 3}; -constexpr int32_t OP_UNION {std::numeric_limits::max() - 4}; - //============================================================================== // Global variables //============================================================================== @@ -705,13 +694,6 @@ extern "C" { bool cell_simple(Cell* c) {return c->simple;} - bool cell_contains(Cell* c, double xyz[3], double uvw[3], int32_t on_surface) - { - Position r {xyz}; - Direction u {uvw}; - return c->contains(r, u, on_surface); - } - void cell_distance(Cell* c, double xyz[3], double uvw[3], int32_t on_surface, double* min_dist, int32_t* i_surf) { diff --git a/src/cell.h b/src/cell.h index 4f3a8978f9..a10c407f0b 100644 --- a/src/cell.h +++ b/src/cell.h @@ -25,6 +25,13 @@ extern "C" int FILL_MATERIAL; extern "C" int FILL_UNIVERSE; extern "C" int FILL_LATTICE; +// TODO: Convert to enum +constexpr int32_t OP_LEFT_PAREN {std::numeric_limits::max()}; +constexpr int32_t OP_RIGHT_PAREN {std::numeric_limits::max() - 1}; +constexpr int32_t OP_COMPLEMENT {std::numeric_limits::max() - 2}; +constexpr int32_t OP_INTERSECTION {std::numeric_limits::max() - 3}; +constexpr int32_t OP_UNION {std::numeric_limits::max() - 4}; + //============================================================================== // Global variables //============================================================================== diff --git a/src/geometry.F90 b/src/geometry.F90 index 09f7066bef..dcd44360be 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -15,16 +15,6 @@ module geometry implicit none interface - function cell_contains_c(cell_ptr, xyz, uvw, on_surface) & - bind(C, name="cell_contains") result(in_cell) - import C_PTR, C_DOUBLE, C_INT32_T, C_BOOL - type(C_PTR), intent(in), value :: cell_ptr - real(C_DOUBLE), intent(in) :: xyz(3) - real(C_DOUBLE), intent(in) :: uvw(3) - integer(C_INT32_T), intent(in), value :: on_surface - logical(C_BOOL) :: in_cell - end function cell_contains_c - function count_universe_instances(search_univ, target_univ_id) bind(C) & result(count) import C_INT32_T, C_INT @@ -38,12 +28,11 @@ module geometry type(Particle), intent(in) :: p end subroutine check_cell_overlap - function find_cell_c(p, n_search_cells, search_cells) & + function find_cell_c(p, search_surf) & bind(C, name="find_cell") result(found) import Particle, C_INT, C_BOOL type(Particle), intent(inout) :: p - integer(C_INT), intent(in), value :: n_search_cells - integer(C_INT), intent(in), optional :: search_cells(n_search_cells) + integer(C_INT), intent(in), value :: search_surf logical(C_BOOL) :: found end function find_cell_c @@ -53,31 +42,26 @@ module geometry type(Particle), intent(inout) :: p integer(C_INT), intent(in) :: lattice_translation(3) end subroutine cross_lattice + + subroutine neighbor_lists() bind(C) + end subroutine neighbor_lists end interface contains - function cell_contains(c, p) result(in_cell) - type(Cell), intent(in) :: c - type(Particle), intent(in) :: p - logical :: in_cell - in_cell = cell_contains_c(c%ptr, p%coord(p%n_coord)%xyz, & - p%coord(p%n_coord)%uvw, p%surface) - end function cell_contains - !=============================================================================== ! FIND_CELL determines what cell a source particle is in within a particular ! universe. If the base universe is passed, the particle should be found as long ! as it's within the geometry !=============================================================================== - recursive subroutine find_cell(p, found, search_cells) - type(Particle), intent(inout) :: p - logical, intent(inout) :: found - integer, optional :: search_cells(:) + subroutine find_cell(p, found, search_surf) + type(Particle), intent(inout) :: p + logical, intent(inout) :: found + integer, optional, intent(in) :: search_surf - if (present(search_cells)) then - found = find_cell_c(p, size(search_cells), search_cells-1) + if (present(search_surf)) then + found = find_cell_c(p, search_surf) else found = find_cell_c(p, 0) end if @@ -201,58 +185,4 @@ contains end subroutine distance_to_boundary -!=============================================================================== -! NEIGHBOR_LISTS builds a list of neighboring cells to each surface to speed up -! searches when a cell boundary is crossed. -!=============================================================================== - - subroutine neighbor_lists() - - integer :: i ! index in cells/surfaces array - integer :: j ! index in region specification - integer :: k ! surface half-space spec - integer :: n ! size of vector - type(VectorInt), allocatable :: neighbor_pos(:) - type(VectorInt), allocatable :: neighbor_neg(:) - - call write_message("Building neighboring cells lists for each surface...", & - 6) - - allocate(neighbor_pos(n_surfaces)) - allocate(neighbor_neg(n_surfaces)) - - do i = 1, n_cells - do j = 1, size(cells(i)%region) - ! Get token from region specification and skip any tokens that - ! correspond to operators rather than regions - k = cells(i)%region(j) - if (abs(k) >= OP_UNION) cycle - - ! Add this cell ID to neighbor list for k-th surface - if (k > 0) then - call neighbor_pos(abs(k))%push_back(i) - else - call neighbor_neg(abs(k))%push_back(i) - end if - end do - end do - - do i = 1, n_surfaces - ! Copy positive neighbors to Surface instance - n = neighbor_pos(i)%size() - if (n > 0) then - allocate(surfaces(i)%neighbor_pos(n)) - surfaces(i)%neighbor_pos(:) = neighbor_pos(i)%data(1:n) - end if - - ! Copy negative neighbors to Surface instance - n = neighbor_neg(i)%size() - if (n > 0) then - allocate(surfaces(i)%neighbor_neg(n)) - surfaces(i)%neighbor_neg(:) = neighbor_neg(i)%data(1:n) - end if - end do - - end subroutine neighbor_lists - end module geometry diff --git a/src/geometry.cpp b/src/geometry.cpp index 85cefb7e2d..87bf0abfda 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -8,6 +8,7 @@ #include "error.h" #include "lattice.h" #include "settings.h" +#include "surface.h" namespace openmc { @@ -50,7 +51,7 @@ check_cell_overlap(Particle* p) { //============================================================================== extern "C" bool -find_cell(Particle* p, int n_search_cells, int* search_cells) { +find_cell(Particle* p, int search_surf) { for (int i = p->n_coord; i < MAX_COORD; i++) { p->coord[i].reset(); } @@ -64,8 +65,18 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) { //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) { + // If a surface was indicated, only search cells from the neighbor list of + // that surface. + int* search_cells; + int n_search_cells; + if (search_surf > 0) { + search_cells = global_surfaces[search_surf-1]->neighbor_pos.data(); + n_search_cells = global_surfaces[search_surf-1]->neighbor_pos.size(); + } else if (search_surf < 0) { + search_cells = global_surfaces[-search_surf-1]->neighbor_neg.data(); + n_search_cells = global_surfaces[-search_surf-1]->neighbor_neg.size(); + } else { + // No surface was indicated, search all cells in the universe. search_cells = global_universes[i_universe]->cells.data(); n_search_cells = global_universes[i_universe]->cells.size(); } @@ -190,7 +201,7 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) { // Update the coordinate level and recurse. ++p->n_coord; - find_cell(p, 0, nullptr); + find_cell(p, 0); } else if (c.type == FILL_LATTICE) { //======================================================================== @@ -237,7 +248,7 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) { // Update the coordinate level and recurse. ++p->n_coord; - find_cell(p, 0, nullptr); + find_cell(p, 0); } } @@ -277,7 +288,7 @@ cross_lattice(Particle* p, int lattice_translation[3]) if (!lat.are_valid_indices(i_xyz)) { // The particle is outside the lattice. Search for it from the base coords. p->n_coord = 1; - bool found = find_cell(p, 0, nullptr); + bool found = find_cell(p, 0); if (!found && p->alive) { std::stringstream err_msg; err_msg << "Could not locate particle " << p->id @@ -288,13 +299,13 @@ cross_lattice(Particle* p, int lattice_translation[3]) } else { // Find cell in next lattice element. p->coord[p->n_coord-1].universe = lat[i_xyz] + 1; - bool found = find_cell(p, 0, nullptr); + bool found = find_cell(p, 0); if (!found) { // A particle crossing the corner of a lattice tile may not be found. In // this case, search for it from the base coords. p->n_coord = 1; - bool found = find_cell(p, 0, nullptr); + bool found = find_cell(p, 0); if (!found && p->alive) { std::stringstream err_msg; err_msg << "Could not locate particle " << p->id diff --git a/src/geometry.h b/src/geometry.h index 47162c73e1..9a13ae2209 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -24,7 +24,7 @@ check_cell_overlap(Particle* p); //============================================================================== extern "C" bool -find_cell(Particle* p, int n_search_cells, int* search_cells); +find_cell(Particle* p, int search_surf); //============================================================================== //! Move a particle into a new lattice tile. diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 6d7d5badb5..6ca3b6d3b3 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -11,6 +11,7 @@ #include "lattice.h" #include "material.h" #include "settings.h" +#include "surface.h" namespace openmc { @@ -152,6 +153,33 @@ find_root_universe() //============================================================================== +void +neighbor_lists() +{ + write_message("Building neighboring cells lists for each surface...", 6); + + for (int i = 0; i < global_cells.size(); i++) { + for (auto token : global_cells[i]->region) { + // Skip operator tokens. + if (std::abs(token) >= OP_UNION) continue; + + // This token is a surface index. Add the cell to the surface's list. + if (token > 0) { + global_surfaces[std::abs(token)-1]->neighbor_pos.push_back(i); + } else { + global_surfaces[std::abs(token)-1]->neighbor_neg.push_back(i); + } + } + } + + for (Surface* surf : global_surfaces) { + surf->neighbor_pos.shrink_to_fit(); + surf->neighbor_neg.shrink_to_fit(); + } +} + +//============================================================================== + void prepare_distribcell(int32_t* filter_cell_list, int n) { diff --git a/src/geometry_aux.h b/src/geometry_aux.h index e3fd034c69..3511a84bee 100644 --- a/src/geometry_aux.h +++ b/src/geometry_aux.h @@ -31,6 +31,12 @@ extern "C" void assign_temperatures(); extern "C" int32_t find_root_universe(); +//!============================================================================= +//! Build a list of neighboring cells to each surface to speed up tracking. +//!============================================================================= + +extern "C" void neighbor_lists(); + //============================================================================== //! Populate all data structures needed for distribcells. //============================================================================== diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 23700d8291..0223a63ff8 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -212,9 +212,6 @@ module geometry_header type Universe integer :: id ! Unique ID integer, allocatable :: cells(:) ! List of cells within - real(8) :: x0 ! Translation in x-coordinate - real(8) :: y0 ! Translation in y-coordinate - real(8) :: z0 ! Translation in z-coordinate end type Universe !=============================================================================== @@ -263,9 +260,6 @@ module geometry_header type Cell type(C_PTR) :: ptr - integer, allocatable :: region(:) ! Definition of spatial region as - ! Boolean expression of half-spaces - ! Rotation matrix and translation vector real(8), allocatable :: rotation(:) real(8), allocatable :: rotation_matrix(:,:) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 27b8be9467..bea0ed6379 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -31,7 +31,7 @@ module input_xml use source_header use stl_vector, only: VectorInt, VectorReal, VectorChar use string, only: to_lower, to_str, str_to_int, str_to_real, & - starts_with, ends_with, tokenize, split_string, & + starts_with, ends_with, split_string, & zero_padded, to_c_string use summary, only: write_summary use tally @@ -1002,16 +1002,14 @@ contains subroutine read_geometry_xml() - integer :: i, j, k + integer :: i, j integer :: n, n_rlats, n_hlats - integer :: id integer :: univ_id integer :: n_cells_in_univ real(8) :: phi, theta, psi logical :: file_exists logical :: boundary_exists character(MAX_LINE_LEN) :: filename - character(:), allocatable :: region_spec type(Cell), pointer :: c class(Lattice), pointer :: lat type(XMLDocument) :: doc @@ -1021,7 +1019,6 @@ contains type(XMLNode), allocatable :: node_cell_list(:) type(XMLNode), allocatable :: node_rlat_list(:) type(XMLNode), allocatable :: node_hlat_list(:) - type(VectorInt) :: tokens type(VectorInt) :: univ_ids ! List of all universe IDs type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each ! universe contains @@ -1104,43 +1101,6 @@ contains // to_str(c % id())) end if - ! Check for region specification (also under deprecated name surfaces) - if (check_for_node(node_cell, "surfaces")) then - call warning("The use of 'surfaces' is deprecated and will be & - &disallowed in a future release. Use 'region' instead. The & - &openmc-update-inputs utility can be used to automatically & - &update geometry.xml files.") - region_spec = node_value_string(node_cell, "surfaces") - call get_node_value(node_cell, "surfaces", region_spec) - elseif (check_for_node(node_cell, "region")) then - region_spec = node_value_string(node_cell, "region") - else - region_spec = '' - end if - - if (len_trim(region_spec) > 0) then - ! Create surfaces array from string - call tokenize(region_spec, tokens) - - ! Convert user IDs to surface indices - do j = 1, tokens % size() - id = tokens % data(j) - if (id < OP_UNION) then - if (surface_dict % has(abs(id))) then - k = surface_dict % get(abs(id)) - tokens % data(j) = sign(k, id) - end if - end if - end do - - ! Copy region spec and RPN form to cell arrays - allocate(c % region(tokens%size())) - c % region(:) = tokens%data(1:tokens%size()) - - call tokens%clear() - end if - if (.not. allocated(c%region)) allocate(c%region(0)) - ! Rotation matrix if (check_for_node(node_cell, "rotation")) then ! Rotations can only be applied to cells that are being filled with diff --git a/src/string.F90 b/src/string.F90 index 71d9d6e966..90b51d7aa0 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -64,98 +64,6 @@ contains end subroutine split_string -!=============================================================================== -! TOKENIZE takes a string that includes logical expressions for a list of -! bounding surfaces in a cell and splits it into separate tokens. The characters -! (, ), |, and ~ count as separate tokens since they represent operators. -!=============================================================================== - - subroutine tokenize(string, tokens) - character(*), intent(in) :: string - type(VectorInt), intent(inout) :: tokens - - integer :: i ! current index - integer :: i_start ! starting index of word - integer :: token - character(len=len_trim(string)) :: string_ - - ! Remove leading blanks - string_ = adjustl(string) - - i_start = 0 - i = 1 - do while (i <= len_trim(string_)) - ! Check for special characters - if (index('()|~ ', string_(i:i)) > 0) then - ! If the special character appears immediately after a non-operator, - ! create a token with the surface half-space - if (i_start > 0) then - call tokens%push_back(int(str_to_int(& - string_(i_start:i - 1)), 4)) - end if - - select case (string_(i:i)) - case ('(') - call tokens%push_back(OP_LEFT_PAREN) - case (')') - if (tokens%size() > 0) then - token = tokens%data(tokens%size()) - if (token >= OP_UNION .and. token < OP_RIGHT_PAREN) then - call fatal_error("Right parentheses cannot follow an operator in & - ®ion specification: " // trim(string)) - end if - end if - call tokens%push_back(OP_RIGHT_PAREN) - case ('|') - if (tokens%size() > 0) then - token = tokens%data(tokens%size()) - if (.not. (token < OP_UNION .or. token == OP_RIGHT_PAREN)) then - call fatal_error("Union cannot follow an operator in region & - &specification: " // trim(string)) - end if - end if - call tokens%push_back(OP_UNION) - case ('~') - call tokens%push_back(OP_COMPLEMENT) - case (' ') - ! Find next non-space character - do while (string_(i+1:i+1) == ' ') - i = i + 1 - end do - - ! If previous token is a halfspace or right parenthesis and next token - ! is not a left parenthese or union operator, that implies that the - ! whitespace is to be interpreted as an intersection operator - if (i_start > 0 .or. tokens%data(tokens%size()) == OP_RIGHT_PAREN) then - if (index(')|', string_(i+1:i+1)) == 0) then - call tokens%push_back(OP_INTERSECTION) - end if - end if - end select - - i_start = 0 - else - ! Check for invalid characters - if (index('-+0123456789', string_(i:i)) == 0) then - call fatal_error("Invalid character '" // string_(i:i) // "' in & - ®ion specification.") - end if - - ! If we haven't yet reached the start of a word, start a new word - if (i_start == 0) i_start = i - end if - - i = i + 1 - end do - - ! If we've reached the end and we're still in a word, create a token from it - ! and add it to the list - if (i_start > 0) then - call tokens%push_back(int(str_to_int(& - string_(i_start:len_trim(string_))), 4)) - end if - end subroutine tokenize - !=============================================================================== ! CONCATENATE takes an array of words and concatenates them together in one ! string with a single space between words diff --git a/src/surface.h b/src/surface.h index 746db5c3b5..c9d8648d8d 100644 --- a/src/surface.h +++ b/src/surface.h @@ -58,11 +58,12 @@ class Surface { public: int id; //!< Unique ID - //int neighbor_pos[], //!< List of cells on positive side - // neighbor_neg[]; //!< List of cells on negative side int bc; //!< Boundary condition std::string name; //!< User-defined name + std::vector neighbor_pos; //!< List of cells on positive side + std::vector neighbor_neg; //!< List of cells on negative side + explicit Surface(pugi::xml_node surf_node); virtual ~Surface() {} diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 4a5834ec2c..401ac90aae 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -84,9 +84,6 @@ module surface_header !=============================================================================== type :: Surface - integer, allocatable :: & - neighbor_pos(:), & ! List of cells on positive side - neighbor_neg(:) ! List of cells on negative side type(C_PTR) :: ptr contains diff --git a/src/tracking.F90 b/src/tracking.F90 index dd0c26b433..e9568fe54a 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -467,21 +467,8 @@ contains ! ========================================================================== ! SEARCH NEIGHBOR LISTS FOR NEXT CELL - if (p % surface > 0 .and. allocated(surf%neighbor_pos)) then - ! If coming from negative side of surface, search all the neighboring - ! cells on the positive side - - call find_cell(p, found, surf%neighbor_pos) - if (found) return - - elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then - ! If coming from positive side of surface, search all the neighboring - ! cells on the negative side - - call find_cell(p, found, surf%neighbor_neg) - if (found) return - - end if + call find_cell(p, found, p % surface) + if (found) return ! ========================================================================== ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS