crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Move neighbor lists to C++

Browse source
This commit is contained in:
Sterling Harper 2018-08-20 23:42:51 -04:00
parent c98142dd23
commit bc1c1e4916
13 changed files with 79 additions and 268 deletions
Download patch file Download diff file Expand all files Collapse all files

18
src/cell.cpp
View file

@ -17,17 +17,6 @@
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
// TODO: Convert to enum
constexpr int32_t OP_LEFT_PAREN {std::numeric_limits<int32_t>::max()};
constexpr int32_t OP_RIGHT_PAREN {std::numeric_limits<int32_t>::max() - 1};
constexpr int32_t OP_COMPLEMENT {std::numeric_limits<int32_t>::max() - 2};
constexpr int32_t OP_INTERSECTION {std::numeric_limits<int32_t>::max() - 3};
constexpr int32_t OP_UNION {std::numeric_limits<int32_t>::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)
{

7
src/cell.h
View file

@ -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<int32_t>::max()};
constexpr int32_t OP_RIGHT_PAREN {std::numeric_limits<int32_t>::max() - 1};
constexpr int32_t OP_COMPLEMENT {std::numeric_limits<int32_t>::max() - 2};
constexpr int32_t OP_INTERSECTION {std::numeric_limits<int32_t>::max() - 3};
constexpr int32_t OP_UNION {std::numeric_limits<int32_t>::max() - 4};
//==============================================================================
// Global variables
//==============================================================================

92
src/geometry.F90
View file

@ -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

27
src/geometry.cpp
View file

@ -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

2
src/geometry.h
View file

@ -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.

28
src/geometry_aux.cpp
View file

@ -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)
{

6
src/geometry_aux.h
View file

@ -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.
//==============================================================================

6
src/geometry_header.F90
View file

@ -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(:,:)

44
src/input_xml.F90
View file

@ -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

92
src/string.F90
View file

@ -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 &
&region 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 &
&region 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

5
src/surface.h
View file

@ -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<int> neighbor_pos; //!< List of cells on positive side
std::vector<int> neighbor_neg; //!< List of cells on negative side
explicit Surface(pugi::xml_node surf_node);
virtual ~Surface() {}

3
src/surface_header.F90
View file

@ -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

17
src/tracking.F90
View file

@ -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