From 5419403e10c6490bd7762ea39e5c4cce5acfa934 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 11 May 2018 18:11:08 -0400 Subject: [PATCH] Move cell instance counting to C++ --- src/cell.cpp | 10 ++- src/cell.h | 3 +- src/geometry.F90 | 66 -------------- src/geometry_aux.cpp | 109 ++++++++++++++--------- src/geometry_header.F90 | 17 +++- src/input_xml.F90 | 19 ++-- src/tallies/tally_filter_distribcell.F90 | 2 +- 7 files changed, 101 insertions(+), 125 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index a3f55d5226..6d7e598eba 100644 --- a/src/cell.cpp +++ b/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) diff --git a/src/cell.h b/src/cell.h index 231bb8ea78..fff8b749c3 100644 --- a/src/cell.h +++ b/src/cell.h @@ -43,7 +43,7 @@ class Universe public: int32_t id; //! Unique ID int32_t type; - std::vector cells; //! Cells within this universe + std::vector 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. diff --git a/src/geometry.F90 b/src/geometry.F90 index 05ad98b86b..be2e2be84d 100644 --- a/src/geometry.F90 +++ b/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 diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 2064bda02d..ec4da6f12a 100644 --- a/src/geometry_aux.cpp +++ b/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 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 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 diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 309076bbe3..012e976afe 100644 --- a/src/geometry_header.F90 +++ b/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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 71408defb0..e43906c9af 100644 --- a/src/input_xml.F90 +++ b/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 diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index 0c369e154a..2b44f3503b 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/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.")