From 110ff69cac8d75181f0f2bef8a861f9f4f868c9f Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 11 May 2018 17:20:25 -0400 Subject: [PATCH] Move cell%type to C++ --- CMakeLists.txt | 1 + src/cell.cpp | 82 +++--------------- src/cell.h | 13 +++ src/constants.F90 | 3 + src/constants.h | 2 + src/geometry.F90 | 30 +++---- src/geometry_aux.cpp | 104 +++++++++++++++++++++++ src/geometry_header.F90 | 38 +++++++-- src/input_xml.F90 | 10 +-- src/plot.F90 | 2 +- src/summary.F90 | 2 +- src/tallies/tally_filter_distribcell.F90 | 16 ++-- 12 files changed, 194 insertions(+), 109 deletions(-) create mode 100644 src/geometry_aux.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index d2a4348fd..dff755999 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -435,6 +435,7 @@ set(LIBOPENMC_CXX_SRC src/cell.cpp src/initialize.cpp src/finalize.cpp + src/geometry_aux.cpp src/hdf5_interface.cpp src/lattice.cpp src/message_passing.cpp diff --git a/src/cell.cpp b/src/cell.cpp index 5bfd9104d..a3f55d522 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -4,7 +4,6 @@ #include #include #include -#include #include "constants.h" #include "error.h" @@ -226,6 +225,15 @@ Cell::Cell(pugi::xml_node cell_node) fill = C_NONE; } + if (check_for_node(cell_node, "material")) { + //TODO: read material ids. + material.push_back(C_NONE+1); + material.shrink_to_fit(); + } else { + material.push_back(C_NONE); + material.shrink_to_fit(); + } + std::string region_spec {""}; if (check_for_node(cell_node, "region")) { region_spec = get_node_value(cell_node, "region"); @@ -443,74 +451,6 @@ read_cells(pugi::xml_node *node) } } -//============================================================================== - -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() -{ - // Change cell.universe values from IDs to indices. - for (Cell *c : cells_c) { - auto it = universe_dict.find(c->universe); - if (it != universe_dict.end()) { - //TODO: Remove this off-by-one indexing. - c->universe = it->second + 1; - } else { - std::stringstream err_msg; - err_msg << "Could not find universe " << c->universe - << " specified on cell " << c->id; - fatal_error(err_msg); - } - } - - // Change all lattice universe values from IDs to indices. - for (Lattice *l : lattices_c) { - l->adjust_indices(); - } -} - //============================================================================== // Fortran compatibility functions //============================================================================== @@ -522,6 +462,10 @@ extern "C" { void cell_set_id(Cell *c, int32_t id) {c->id = id;} + int cell_type(Cell *c) {return c->type;} + + void cell_set_type(Cell *c, int type) {c->type = type;} + int32_t cell_universe(Cell *c) {return c->universe;} void cell_set_universe(Cell *c, int32_t universe) {c->universe = universe;} diff --git a/src/cell.h b/src/cell.h index 5c06d8120..231bb8ea7 100644 --- a/src/cell.h +++ b/src/cell.h @@ -12,6 +12,14 @@ namespace openmc { +//============================================================================== +// Constants +//============================================================================== + +extern "C" int FILL_MATERIAL; +extern "C" int FILL_UNIVERSE; +extern "C" int FILL_LATTICE; + //============================================================================== // Global variables //============================================================================== @@ -48,9 +56,14 @@ class Cell public: int32_t id; //!< Unique ID std::string name{""}; //!< User-defined name + int type; //!< Material, universe, or lattice int32_t universe; //!< Universe # this cell is in int32_t fill; //!< Universe # filling this cell + //! Material within this cell. May be multiple materials for distribcell. + //! C_NONE signifies a universe. + std::vector material; + //! Definition of spatial region as Boolean expression of half-spaces std::vector region; //! Reverse Polish notation for region expression diff --git a/src/constants.F90 b/src/constants.F90 index d1893bf9e..dbd4181da 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -113,6 +113,9 @@ module constants FILL_MATERIAL = 1, & ! Cell with a specified material FILL_UNIVERSE = 2, & ! Cell filled by a separate universe FILL_LATTICE = 3 ! Cell filled with a lattice + integer(C_INT), bind(C, name='FILL_MATERIAL') :: FILL_MATERIAL_C = FILL_MATERIAL + integer(C_INT), bind(C, name='FILL_UNIVERSE') :: FILL_UNIVERSE_C = FILL_UNIVERSE + integer(C_INT), bind(C, name='FILL_LATTICE') :: FILL_LATTICE_C = FILL_LATTICE ! Void material integer, parameter :: MATERIAL_VOID = -1 diff --git a/src/constants.h b/src/constants.h index 8cc837a73..685fbcfd7 100644 --- a/src/constants.h +++ b/src/constants.h @@ -1,6 +1,8 @@ #ifndef CONSTANTS_H #define CONSTANTS_H +#include + namespace openmc{ diff --git a/src/geometry.F90 b/src/geometry.F90 index 35e4d3ca5..05ad98b86 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -150,7 +150,7 @@ contains if (found) then associate(c => cells(i_cell)) - CELL_TYPE: if (c % type == FILL_MATERIAL) then + CELL_TYPE: if (c % type() == FILL_MATERIAL) then ! ====================================================================== ! AT LOWEST UNIVERSE, TERMINATE SEARCH @@ -166,10 +166,10 @@ contains distribcell_index = c % distribcell_index offset = 0 do k = 1, p % n_coord - if (cells(p % coord(k) % cell) % type == FILL_UNIVERSE) then + if (cells(p % coord(k) % cell) % type() == FILL_UNIVERSE) then offset = offset + cells(p % coord(k) % cell) % & offset(distribcell_index) - elseif (cells(p % coord(k) % cell) % type == FILL_LATTICE) then + elseif (cells(p % coord(k) % cell) % type() == FILL_LATTICE) then if (lattices(p % coord(k + 1) % lattice) % obj & % are_valid_indices([& p % coord(k + 1) % lattice_x, & @@ -204,7 +204,7 @@ contains p % sqrtkT = c % sqrtkT(1) end if - elseif (c % type == FILL_UNIVERSE) then CELL_TYPE + elseif (c % type() == FILL_UNIVERSE) then CELL_TYPE ! ====================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL @@ -232,7 +232,7 @@ contains call find_cell(p, found) j = p % n_coord - elseif (c % type == FILL_LATTICE) then CELL_TYPE + elseif (c % type() == FILL_LATTICE) then CELL_TYPE ! ====================================================================== ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL @@ -561,11 +561,11 @@ contains ! ==================================================================== ! AT LOWEST UNIVERSE, TERMINATE SEARCH - if (c % type == FILL_MATERIAL) then + if (c % type() == FILL_MATERIAL) then ! ==================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - elseif (c % type == FILL_UNIVERSE) then + elseif (c % type() == FILL_UNIVERSE) then ! Set offset for the cell on this level c % offset(map) = offset @@ -579,7 +579,7 @@ contains ! ==================================================================== ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then + elseif (c % type() == FILL_LATTICE) then ! Set current lattice lat => lattices(c % fill) % obj @@ -671,11 +671,11 @@ contains ! ==================================================================== ! AT LOWEST UNIVERSE, TERMINATE SEARCH - if (c % type == FILL_MATERIAL) then + if (c % type() == FILL_MATERIAL) then ! ==================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - elseif (c % type == FILL_UNIVERSE) then + elseif (c % type() == FILL_UNIVERSE) then next_univ => universes(c % fill) @@ -690,7 +690,7 @@ contains ! ==================================================================== ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then + elseif (c % type() == FILL_LATTICE) then ! Set current lattice lat => lattices(c % fill) % obj @@ -771,13 +771,13 @@ contains ! ==================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - if (c % type == FILL_UNIVERSE) then + if (c % type() == FILL_UNIVERSE) then call count_instance(universes(c % fill)) ! ==================================================================== ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then + elseif (c % type() == FILL_LATTICE) then ! Set current lattice associate (lat => lattices(c % fill) % obj) @@ -840,14 +840,14 @@ contains ! ==================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - if (c % type == FILL_UNIVERSE) then + if (c % type() == FILL_UNIVERSE) then next_univ => universes(c % fill) levels_below = max(levels_below, maximum_levels(next_univ)) ! ==================================================================== ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then + elseif (c % type() == FILL_LATTICE) then ! Set current lattice lat => lattices(c % fill) % obj diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp new file mode 100644 index 000000000..2064bda02 --- /dev/null +++ b/src/geometry_aux.cpp @@ -0,0 +1,104 @@ +#include +#include + +#include "cell.h" +#include "constants.h" +#include "error.h" +#include "lattice.h" + +#include //TODO: remove this + + +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() +{ + // Adjust material/fill idices. + for (Cell *c : cells_c) { + if (c->material[0] == C_NONE) { + int32_t id = c->fill; + auto search_univ = universe_dict.find(id); + auto search_lat = lattice_dict.find(id); + if (search_univ != universe_dict.end()) { + c->type = FILL_UNIVERSE; + c->fill = search_univ->second + 1; //TODO: off-by-one + } else if (search_lat != lattice_dict.end()) { + c->type = FILL_LATTICE; + c->fill = search_lat->second + 1; //TODO: off-by-one + } else { + std::stringstream err_msg; + err_msg << "Specified fill " << id << " on cell " << c->id + << " is neither a universe nor a lattice."; + fatal_error(err_msg); + } + } else { + //TODO: materials + c->type = FILL_MATERIAL; + } + } + + // Change cell.universe values from IDs to indices. + for (Cell *c : cells_c) { + auto search = universe_dict.find(c->universe); + if (search != universe_dict.end()) { + //TODO: Remove this off-by-one indexing. + c->universe = search->second + 1; + } else { + std::stringstream err_msg; + err_msg << "Could not find universe " << c->universe + << " specified on cell " << c->id; + fatal_error(err_msg); + } + } + + // Change all lattice universe values from IDs to indices. + for (Lattice *l : lattices_c) { + l->adjust_indices(); + } +} + +} // namespace openmc diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index ebce31c45..309076bbe 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -37,6 +37,20 @@ module geometry_header integer(C_INT32_T), intent(in), value :: id end subroutine cell_set_id_c + function cell_type_c(cell_ptr) bind(C, name='cell_type') result(type) + import C_PTR, C_INT + implicit none + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT) :: type + end function cell_type_c + + subroutine cell_set_type_c(cell_ptr, type) bind(C, name='cell_set_type') + import C_PTR, C_INT + implicit none + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT), intent(in), value :: type + end subroutine cell_set_type_c + function cell_universe_c(cell_ptr) bind(C, name='cell_universe') & result(universe) import C_PTR, C_INT32_T @@ -230,8 +244,6 @@ module geometry_header type Cell type(C_PTR) :: ptr - integer :: type ! Type of cell (normal, universe, - ! lattice) integer :: fill ! universe # filling this cell integer :: instances ! number of instances of this cell in ! the geom @@ -257,6 +269,8 @@ module geometry_header procedure :: id => cell_id procedure :: set_id => cell_set_id + procedure :: type => cell_type + procedure :: set_type => cell_set_type procedure :: universe => cell_universe procedure :: set_universe => cell_set_universe procedure :: simple => cell_simple @@ -355,6 +369,18 @@ contains call cell_set_id_c(this % ptr, id) end subroutine cell_set_id + function cell_type(this) result(type) + class(Cell), intent(in) :: this + integer(C_INT) :: type + type = cell_type_c(this % ptr) + end function cell_type + + subroutine cell_set_type(this, type) + class(Cell), intent(in) :: this + integer(C_INT), intent(in) :: type + call cell_set_type_c(this % ptr, type) + end subroutine cell_set_type + function cell_universe(this) result(universe) class(Cell), intent(in) :: this integer(C_INT32_T) :: universe @@ -544,7 +570,7 @@ contains err = 0 if (index >= 1 .and. index <= size(cells)) then associate (c => cells(index)) - type = c % type + type = c % type() select case (type) case (FILL_MATERIAL) n = size(c % material) @@ -595,7 +621,7 @@ contains if (allocated(c % material)) deallocate(c % material) allocate(c % material(n)) - c % type = FILL_MATERIAL + call c % set_type(FILL_MATERIAL) do i = 1, n j = indices(i) if ((j >= 1 .and. j <= n_materials) .or. j == MATERIAL_VOID) then @@ -607,9 +633,9 @@ contains end if end do case (FILL_UNIVERSE) - c % type = FILL_UNIVERSE + call c % set_type(FILL_UNIVERSE) case (FILL_LATTICE) - c % type = FILL_LATTICE + call c % set_type(FILL_LATTICE) end select end associate else diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 09da2d1f4..71408defb 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3794,24 +3794,16 @@ contains if (c % material(1) == NONE) then id = c % fill if (universe_dict % has(id)) then - c % type = FILL_UNIVERSE c % fill = universe_dict % get(id) elseif (lattice_dict % has(id)) then lid = lattice_dict % get(id) - c % type = FILL_LATTICE c % fill = lid - else - call fatal_error("Specified fill " // trim(to_str(id)) // " on cell "& - // trim(to_str(c % id())) // " is neither a universe nor a & - &lattice.") end if else do j = 1, size(c % material) id = c % material(j) if (id == MATERIAL_VOID) then - c % type = FILL_MATERIAL else if (material_dict % has(id)) then - c % type = FILL_MATERIAL c % material(j) = material_dict % get(id) else call fatal_error("Could not find material " // trim(to_str(id)) & @@ -3989,7 +3981,7 @@ contains ! Allocate offset table for fill cells do i = 1, n_cells - if (cells(i) % type /= FILL_MATERIAL) then + if (cells(i) % type() /= FILL_MATERIAL) then allocate(cells(i) % offset(n_maps)) end if end do diff --git a/src/plot.F90 b/src/plot.F90 index 2613acd9d..789ad0ee9 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -85,7 +85,7 @@ contains if (pl % color_by == PLOT_COLOR_MATS) then ! Assign color based on material associate (c => cells(p % coord(j) % cell)) - if (c % type == FILL_UNIVERSE) then + if (c % type() == FILL_UNIVERSE) then ! If we stopped on a middle universe level, treat as if not found rgb = pl % not_found % rgb id = -1 diff --git a/src/summary.F90 b/src/summary.F90 index fb934726e..964ec34a3 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -193,7 +193,7 @@ contains call c % to_hdf5(cell_group) ! Write information on what fills this cell - select case (c%type) + select case (c%type()) case (FILL_MATERIAL) call write_dataset(cell_group, "fill_type", "material") diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index dd12e808c..0c369e154 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -58,10 +58,10 @@ contains distribcell_index = cells(this % cell) % distribcell_index offset = 0 do i = 1, p % n_coord - if (cells(p % coord(i) % cell) % type == FILL_UNIVERSE) then + if (cells(p % coord(i) % cell) % type() == FILL_UNIVERSE) then offset = offset + cells(p % coord(i) % cell) % & offset(distribcell_index) - elseif (cells(p % coord(i) % cell) % type == FILL_LATTICE) then + elseif (cells(p % coord(i) % cell) % type() == FILL_LATTICE) then if (lattices(p % coord(i + 1) % lattice) % obj & % are_valid_indices([& p % coord(i + 1) % lattice_x, & @@ -197,20 +197,20 @@ contains c => cells(univ % cells(j)) ! Skip normal cells which do not have offsets - if (c % type == FILL_MATERIAL) cycle + if (c % type() == FILL_MATERIAL) cycle ! Break loop once we've found the next cell with an offset exit end do ! Ensure we didn't just end the loop by iteration - if (c % type /= FILL_MATERIAL) then + if (c % type() /= FILL_MATERIAL) then ! There are more cells in this universe that it could be in later_cell = .true. ! Two cases, lattice or fill cell - if (c % type == FILL_UNIVERSE) then + if (c % type() == FILL_UNIVERSE) then temp_offset = c % offset(map) ! Get the offset of the first lattice location @@ -227,7 +227,7 @@ contains end if end if - if (n == 1 .and. c % type /= FILL_MATERIAL) then + if (n == 1 .and. c % type() /= FILL_MATERIAL) then this_cell = .true. end if @@ -245,7 +245,7 @@ contains ! ==================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - if (c % type == FILL_UNIVERSE) then + if (c % type() == FILL_UNIVERSE) then ! Enter this cell to update the current offset offset = c % offset(map) + offset @@ -256,7 +256,7 @@ contains ! ==================================================================== ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then + elseif (c % type() == FILL_LATTICE) then ! Set current lattice lat => lattices(c % fill) % obj