diff --git a/CMakeLists.txt b/CMakeLists.txt index dc5558de18..53aac68c3f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -432,9 +432,12 @@ set(LIBOPENMC_FORTRAN_SRC src/tallies/trigger_header.F90 ) 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/math_functions.cpp src/message_passing.cpp src/plot.cpp diff --git a/docs/source/devguide/styleguide.rst b/docs/source/devguide/styleguide.rst index 89d5673943..ccf6bf5217 100644 --- a/docs/source/devguide/styleguide.rst +++ b/docs/source/devguide/styleguide.rst @@ -219,8 +219,8 @@ Curly braces For a function definition, the opening and closing braces should each be on their own lines. This helps distinguish function code from the argument list. -If the entire function fits on one line, then the braces can be on the same -line. e.g.: +If the entire function fits on one or two lines, then the braces can be on the +same line. e.g.: .. code-block:: C++ @@ -238,6 +238,9 @@ line. e.g.: int return_one() {return 1;} + int return_one() + {return 1;} + For a conditional, the opening brace should be on the same line as the end of the conditional statement. If there is a following ``else if`` or ``else`` statement, the closing brace should be on the same line as that following diff --git a/src/api.F90 b/src/api.F90 index 1740fa2abd..03b19cee3c 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -207,7 +207,7 @@ contains if (found) then if (rtype == 1) then - id = cells(p % coord(p % n_coord) % cell) % id + id = cells(p % coord(p % n_coord) % cell) % id() elseif (rtype == 2) then if (p % material == MATERIAL_VOID) then id = 0 diff --git a/src/cell.cpp b/src/cell.cpp new file mode 100644 index 0000000000..e4388e7012 --- /dev/null +++ b/src/cell.cpp @@ -0,0 +1,523 @@ +#include "cell.h" + +#include +#include +#include +#include + +#include "constants.h" +#include "error.h" +#include "hdf5_interface.h" +#include "lattice.h" +#include "surface.h" +#include "xml_interface.h" + + +namespace openmc { + +//============================================================================== +// Constants +//============================================================================== + +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}; + +extern "C" double FP_PRECISION; + +//============================================================================== +// Global variables +//============================================================================== + +int32_t n_cells {0}; + +std::vector global_cells; +std::unordered_map cell_map; + +std::vector global_universes; +std::unordered_map universe_map; + +//============================================================================== +//! Convert region specification string to integer tokens. +//! +//! The characters (, ), |, and ~ count as separate tokens since they represent +//! operators. +//============================================================================== + +std::vector +tokenize(const std::string region_spec) { + // Check for an empty region_spec first. + std::vector tokens; + if (region_spec.empty()) { + return tokens; + } + + // Parse all halfspaces and operators except for intersection (whitespace). + for (int i = 0; i < region_spec.size(); ) { + if (region_spec[i] == '(') { + tokens.push_back(OP_LEFT_PAREN); + i++; + + } else if (region_spec[i] == ')') { + tokens.push_back(OP_RIGHT_PAREN); + i++; + + } else if (region_spec[i] == '|') { + tokens.push_back(OP_UNION); + i++; + + } else if (region_spec[i] == '~') { + tokens.push_back(OP_COMPLEMENT); + i++; + + } else if (region_spec[i] == '-' || region_spec[i] == '+' + || std::isdigit(region_spec[i])) { + // This is the start of a halfspace specification. Iterate j until we + // find the end, then push-back everything between i and j. + int j = i + 1; + while (j < region_spec.size() && std::isdigit(region_spec[j])) {j++;} + tokens.push_back(std::stoi(region_spec.substr(i, j-i))); + i = j; + + } else if (std::isspace(region_spec[i])) { + i++; + + } else { + std::stringstream err_msg; + err_msg << "Region specification contains invalid character, \"" + << region_spec[i] << "\""; + fatal_error(err_msg); + } + } + + // Add in intersection operators where a missing operator is needed. + int i = 0; + while (i < tokens.size()-1) { + bool left_compat {(tokens[i] < OP_UNION) || (tokens[i] == OP_RIGHT_PAREN)}; + bool right_compat {(tokens[i+1] < OP_UNION) + || (tokens[i+1] == OP_LEFT_PAREN) + || (tokens[i+1] == OP_COMPLEMENT)}; + if (left_compat && right_compat) { + tokens.insert(tokens.begin()+i+1, OP_INTERSECTION); + } + i++; + } + + return tokens; +} + +//============================================================================== +//! Convert infix region specification to Reverse Polish Notation (RPN) +//! +//! This function uses the shunting-yard algorithm. +//============================================================================== + +std::vector +generate_rpn(int32_t cell_id, std::vector infix) +{ + std::vector rpn; + std::vector stack; + + for (int32_t token : infix) { + if (token < OP_UNION) { + // If token is not an operator, add it to output + rpn.push_back(token); + + } else if (token < OP_RIGHT_PAREN) { + // Regular operators union, intersection, complement + while (stack.size() > 0) { + int32_t op = stack.back(); + + if (op < OP_RIGHT_PAREN && + ((token == OP_COMPLEMENT && token < op) || + (token != OP_COMPLEMENT && token <= op))) { + // While there is an operator, op, on top of the stack, if the token + // is left-associative and its precedence is less than or equal to + // that of op or if the token is right-associative and its precedence + // is less than that of op, move op to the output queue and push the + // token on to the stack. Note that only complement is + // right-associative. + rpn.push_back(op); + stack.pop_back(); + } else { + break; + } + } + + stack.push_back(token); + + } else if (token == OP_LEFT_PAREN) { + // If the token is a left parenthesis, push it onto the stack + stack.push_back(token); + + } else { + // If the token is a right parenthesis, move operators from the stack to + // the output queue until reaching the left parenthesis. + for (auto it = stack.rbegin(); *it != OP_LEFT_PAREN; it++) { + // If we run out of operators without finding a left parenthesis, it + // means there are mismatched parentheses. + if (it == stack.rend()) { + std::stringstream err_msg; + err_msg << "Mismatched parentheses in region specification for cell " + << cell_id; + fatal_error(err_msg); + } + + rpn.push_back(stack.back()); + stack.pop_back(); + } + + // Pop the left parenthesis. + stack.pop_back(); + } + } + + while (stack.size() > 0) { + int32_t op = stack.back(); + + // If the operator is a parenthesis it is mismatched. + if (op >= OP_RIGHT_PAREN) { + std::stringstream err_msg; + err_msg << "Mismatched parentheses in region specification for cell " + << cell_id; + fatal_error(err_msg); + } + + rpn.push_back(stack.back()); + stack.pop_back(); + } + + return rpn; +} + +//============================================================================== +// Cell implementation +//============================================================================== + +Cell::Cell(pugi::xml_node cell_node) +{ + if (check_for_node(cell_node, "id")) { + id = stoi(get_node_value(cell_node, "id")); + } else { + fatal_error("Must specify id of cell in geometry XML file."); + } + + //TODO: don't automatically lowercase cell and surface names + if (check_for_node(cell_node, "name")) { + name = get_node_value(cell_node, "name"); + } + + if (check_for_node(cell_node, "universe")) { + universe = stoi(get_node_value(cell_node, "universe")); + } else { + universe = 0; + } + + if (check_for_node(cell_node, "fill")) { + fill = stoi(get_node_value(cell_node, "fill")); + } else { + 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(); + } + + // Make sure that either material or fill was specified. + if ((material[0] == C_NONE) && (fill == C_NONE)) { + std::stringstream err_msg; + err_msg << "Neither material nor fill was specified for cell " << id; + fatal_error(err_msg); + } + + // Make sure that material and fill haven't been specified simultaneously. + if ((material[0] != C_NONE) && (fill != C_NONE)) { + std::stringstream err_msg; + err_msg << "Cell " << id << " has both a material and a fill specified; " + << "only one can be specified per cell"; + fatal_error(err_msg); + } + + // Read the region specification. + std::string region_spec; + if (check_for_node(cell_node, "region")) { + region_spec = get_node_value(cell_node, "region"); + } + + // Get a tokenized representation of the region specification. + region = tokenize(region_spec); + region.shrink_to_fit(); + + // Convert user IDs to surface indices. + for (auto &r : region) { + if (r < OP_UNION) { + r = copysign(surface_map[abs(r)] + 1, r); + } + } + + // Convert the infix region spec to RPN. + rpn = generate_rpn(id, region); + rpn.shrink_to_fit(); + + // Check if this is a simple cell. + simple = true; + for (int32_t token : rpn) { + if ((token == OP_COMPLEMENT) || (token == OP_UNION)) { + simple = false; + break; + } + } +} + +//============================================================================== + +bool +Cell::contains(const double xyz[3], const double uvw[3], + int32_t on_surface) const +{ + if (simple) { + return contains_simple(xyz, uvw, on_surface); + } else { + return contains_complex(xyz, uvw, on_surface); + } +} + +//============================================================================== + +std::pair +Cell::distance(const double xyz[3], const double uvw[3], + int32_t on_surface) const +{ + double min_dist {INFTY}; + int32_t i_surf {std::numeric_limits::max()}; + + for (int32_t token : rpn) { + // Ignore this token if it corresponds to an operator rather than a region. + if (token >= OP_UNION) {continue;} + + // Calculate the distance to this surface. + // Note the off-by-one indexing + bool coincident {token == on_surface}; + double d {surfaces_c[abs(token)-1]->distance(xyz, uvw, coincident)}; + + // Check if this distance is the new minimum. + if (d < min_dist) { + if (std::abs(d - min_dist) / min_dist >= FP_PRECISION) { + min_dist = d; + i_surf = -token; + } + } + } + + return {min_dist, i_surf}; +} + +//============================================================================== + +void +Cell::to_hdf5(hid_t cell_group) const +{ + if (!name.empty()) { + write_string(cell_group, "name", name, false); + } + + //TODO: Fix the off-by-one indexing. + write_int(cell_group, 0, nullptr, "universe", + &global_universes[universe-1]->id, false); + + // Write the region specification. + if (!region.empty()) { + std::stringstream region_spec {}; + for (int32_t token : region) { + if (token == OP_LEFT_PAREN) { + region_spec << " ("; + } else if (token == OP_RIGHT_PAREN) { + region_spec << " )"; + } else if (token == OP_COMPLEMENT) { + region_spec << " ~"; + } else if (token == OP_INTERSECTION) { + } else if (token == OP_UNION) { + region_spec << " |"; + } else { + // Note the off-by-one indexing + region_spec << " " << copysign(surfaces_c[abs(token)-1]->id, token); + } + } + write_string(cell_group, "region", region_spec.str(), false); + } +} + +//============================================================================== + +bool +Cell::contains_simple(const double xyz[3], const double uvw[3], + int32_t on_surface) const +{ + for (int32_t token : rpn) { + if (token < OP_UNION) { + // If the token is not an operator, evaluate the sense of particle with + // respect to the surface and see if the token matches the sense. If the + // particle's surface attribute is set and matches the token, that + // overrides the determination based on sense(). + if (token == on_surface) { + } else if (-token == on_surface) { + return false; + } else { + // Note the off-by-one indexing + bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw); + if (sense != (token > 0)) {return false;} + } + } + } + return true; +} + +//============================================================================== + +bool +Cell::contains_complex(const double xyz[3], const double uvw[3], + int32_t on_surface) const +{ + // Make a stack of booleans. We don't know how big it needs to be, but we do + // know that rpn.size() is an upper-bound. + bool stack[rpn.size()]; + int i_stack = -1; + + for (int32_t token : rpn) { + // If the token is a binary operator (intersection/union), apply it to + // the last two items on the stack. If the token is a unary operator + // (complement), apply it to the last item on the stack. + if (token == OP_UNION) { + stack[i_stack-1] = stack[i_stack-1] || stack[i_stack]; + i_stack --; + } else if (token == OP_INTERSECTION) { + stack[i_stack-1] = stack[i_stack-1] && stack[i_stack]; + i_stack --; + } else if (token == OP_COMPLEMENT) { + stack[i_stack] = !stack[i_stack]; + } else { + // If the token is not an operator, evaluate the sense of particle with + // respect to the surface and see if the token matches the sense. If the + // particle's surface attribute is set and matches the token, that + // overrides the determination based on sense(). + i_stack ++; + if (token == on_surface) { + stack[i_stack] = true; + } else if (-token == on_surface) { + stack[i_stack] = false; + } else { + // Note the off-by-one indexing + bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw);; + stack[i_stack] = (sense == (token > 0)); + } + } + } + + if (i_stack == 0) { + // The one remaining bool on the stack indicates whether the particle is + // in the cell. + return stack[i_stack]; + } else { + // This case occurs if there is no region specification since i_stack will + // still be -1. + return true; + } +} + +//============================================================================== +// Non-method functions +//============================================================================== + +extern "C" void +read_cells(pugi::xml_node *node) +{ + // Count the number of cells. + for (pugi::xml_node cell_node: node->children("cell")) {n_cells++;} + if (n_cells == 0) { + fatal_error("No cells found in geometry.xml!"); + } + + // Allocate the vector of Cells. + global_cells.reserve(n_cells); + + // Loop over XML cell elements and populate the array. + for (pugi::xml_node cell_node: node->children("cell")) { + global_cells.push_back(new Cell(cell_node)); + } + + // Populate the Universe vector and map. + for (int i = 0; i < global_cells.size(); i++) { + int32_t uid = global_cells[i]->universe; + auto it = universe_map.find(uid); + if (it == universe_map.end()) { + global_universes.push_back(new Universe()); + global_universes.back()->id = uid; + global_universes.back()->cells.push_back(i); + universe_map[uid] = global_universes.size() - 1; + } else { + global_universes[it->second]->cells.push_back(i); + } + } +} + +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +extern "C" { + Cell* cell_pointer(int32_t cell_ind) {return global_cells[cell_ind];} + + int32_t cell_id(Cell *c) {return c->id;} + + 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;} + + int32_t cell_fill(Cell *c) {return c->fill;} + + int32_t* cell_fill_ptr(Cell *c) {return &c->fill;} + + 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) + {return c->contains(xyz, uvw, on_surface);} + + void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface, + double *min_dist, int32_t *i_surf) + { + std::pair out = c->distance(xyz, uvw, on_surface); + *min_dist = out.first; + *i_surf = out.second; + } + + int32_t cell_offset(Cell *c, int map) {return c->offset[map];} + + void cell_to_hdf5(Cell *c, hid_t group) {c->to_hdf5(group);} + + void extend_cells_c(int32_t n) + { + global_cells.reserve(global_cells.size() + n); + for (int32_t i = 0; i < n; i++) { + global_cells.push_back(new Cell()); + } + n_cells = global_cells.size(); + } +} + + +} // namespace openmc diff --git a/src/cell.h b/src/cell.h new file mode 100644 index 0000000000..a2f72c8337 --- /dev/null +++ b/src/cell.h @@ -0,0 +1,118 @@ +#ifndef CELL_H +#define CELL_H + +#include +#include +#include +#include + +#include "hdf5.h" +#include "pugixml/pugixml.hpp" + + +namespace openmc { + +//============================================================================== +// Constants +//============================================================================== + +extern "C" int FILL_MATERIAL; +extern "C" int FILL_UNIVERSE; +extern "C" int FILL_LATTICE; + +//============================================================================== +// Global variables +//============================================================================== + +extern "C" int32_t n_cells; + +class Cell; +extern std::vector global_cells; +extern std::unordered_map cell_map; + +class Universe; +extern std::vector global_universes; +extern std::unordered_map universe_map; + +//============================================================================== +//! A geometry primitive that fills all space and contains cells. +//============================================================================== + +class Universe +{ +public: + int32_t id; //!< Unique ID + std::vector cells; //!< Cells within this universe + //double x0, y0, z0; //!< Translation coordinates. +}; + +//============================================================================== +//! A geometry primitive that links surfaces, universes, and materials +//============================================================================== + +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 + int32_t n_instances{0}; //!< Number of instances of this cell + + //! \brief Material(s) 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 + std::vector rpn; + bool simple; //!< Does the region contain only intersections? + + std::vector offset; //!< Distribcell offset table + + Cell() {}; + + explicit Cell(pugi::xml_node cell_node); + + //! \brief Determine if a cell contains the particle at a given location. + //! + //! The bounds of the cell are detemined by a logical expression involving + //! surface half-spaces. At initialization, the expression was converted + //! to RPN notation. + //! + //! The function is split into two cases, one for simple cells (those + //! involving only the intersection of half-spaces) and one for complex cells. + //! Simple cells can be evaluated with short circuit evaluation, i.e., as soon + //! as we know that one half-space is not satisfied, we can exit. This + //! provides a performance benefit for the common case. In + //! contains_complex, we evaluate the RPN expression using a stack, similar to + //! how a RPN calculator would work. + //! @param xyz[3] The 3D Cartesian coordinate to check. + //! @param uvw[3] A direction used to "break ties" the coordinates are very + //! close to a surface. + //! @param on_surface The signed index of a surface that the coordinate is + //! known to be on. This index takes precedence over surface sense + //! calculations. + bool + contains(const double xyz[3], const double uvw[3], int32_t on_surface) const; + + //! Find the oncoming boundary of this cell. + std::pair + distance(const double xyz[3], const double uvw[3], int32_t on_surface) const; + + //! \brief Write cell information to an HDF5 group. + //! @param group_id An HDF5 group id. + void to_hdf5(hid_t group_id) const; + +protected: + bool contains_simple(const double xyz[3], const double uvw[3], + int32_t on_surface) const; + bool contains_complex(const double xyz[3], const double uvw[3], + int32_t on_surface) const; +}; + +} // namespace openmc +#endif // CELL_H diff --git a/src/constants.F90 b/src/constants.F90 index d1893bf9e1..d75ee4e54a 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 @@ -418,6 +421,7 @@ module constants ! indicates that an array index hasn't been set integer, parameter :: NONE = 0 + integer, parameter :: C_NONE = -1 ! Codes for read errors -- better hope these numbers are never used in an ! input file! diff --git a/src/constants.h b/src/constants.h new file mode 100644 index 0000000000..685fbcfd72 --- /dev/null +++ b/src/constants.h @@ -0,0 +1,16 @@ +#ifndef CONSTANTS_H +#define CONSTANTS_H + +#include + + +namespace openmc{ + +extern "C" double FP_COINCIDENT; +extern "C" double FP_PRECISION; +constexpr double INFTY {std::numeric_limits::max()}; +constexpr int C_NONE {-1}; + +} // namespace openmc + +#endif // CONSTANTS_H diff --git a/src/geometry.F90 b/src/geometry.F90 index 651005631e..457fdf3786 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -14,125 +14,36 @@ 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 + integer(C_INT32_T), intent(in), value :: search_univ + integer(C_INT32_T), intent(in), value :: target_univ_id + integer(C_INT) :: count + end function + end interface + contains -!=============================================================================== -! CELL_CONTAINS determines if a cell contains the particle at a given -! location. The bounds of the cell are detemined by a logical expression -! involving surface half-spaces. At initialization, the expression was converted -! to RPN notation. -! -! The function is split into two cases, one for simple cells (those involving -! only the intersection of half-spaces) and one for complex cells. Simple cells -! can be evaluated with short circuit evaluation, i.e., as soon as we know that -! one half-space is not satisfied, we can exit. This provides a performance -! benefit for the common case. In complex_cell_contains, we evaluate the RPN -! expression using a stack, similar to how a RPN calculator would work. -!=============================================================================== - - pure function cell_contains(c, p) result(in_cell) + function cell_contains(c, p) result(in_cell) type(Cell), intent(in) :: c type(Particle), intent(in) :: p logical :: in_cell - - if (c%simple) then - in_cell = simple_cell_contains(c, p) - else - in_cell = complex_cell_contains(c, p) - end if + 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 - pure function simple_cell_contains(c, p) result(in_cell) - type(Cell), intent(in) :: c - type(Particle), intent(in) :: p - logical :: in_cell - - integer :: i - integer :: token - logical :: actual_sense ! sense of particle wrt surface - - in_cell = .true. - do i = 1, size(c%rpn) - token = c%rpn(i) - if (token < OP_UNION) then - ! If the token is not an operator, evaluate the sense of particle with - ! respect to the surface and see if the token matches the sense. If the - ! particle's surface attribute is set and matches the token, that - ! overrides the determination based on sense(). - if (token == p%surface) then - cycle - elseif (-token == p%surface) then - in_cell = .false. - exit - else - actual_sense = surfaces(abs(token)) % sense( & - p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw) - if (actual_sense .neqv. (token > 0)) then - in_cell = .false. - exit - end if - end if - end if - end do - end function simple_cell_contains - - pure function complex_cell_contains(c, p) result(in_cell) - type(Cell), intent(in) :: c - type(Particle), intent(in) :: p - logical :: in_cell - - integer :: i - integer :: token - integer :: i_stack - logical :: actual_sense ! sense of particle wrt surface - logical :: stack(size(c%rpn)) - - i_stack = 0 - do i = 1, size(c%rpn) - token = c%rpn(i) - - ! If the token is a binary operator (intersection/union), apply it to - ! the last two items on the stack. If the token is a unary operator - ! (complement), apply it to the last item on the stack. - select case (token) - case (OP_UNION) - stack(i_stack - 1) = stack(i_stack - 1) .or. stack(i_stack) - i_stack = i_stack - 1 - case (OP_INTERSECTION) - stack(i_stack - 1) = stack(i_stack - 1) .and. stack(i_stack) - i_stack = i_stack - 1 - case (OP_COMPLEMENT) - stack(i_stack) = .not. stack(i_stack) - case default - ! If the token is not an operator, evaluate the sense of particle with - ! respect to the surface and see if the token matches the sense. If the - ! particle's surface attribute is set and matches the token, that - ! overrides the determination based on sense(). - i_stack = i_stack + 1 - if (token == p%surface) then - stack(i_stack) = .true. - elseif (-token == p%surface) then - stack(i_stack) = .false. - else - actual_sense = surfaces(abs(token)) % sense( & - p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw) - stack(i_stack) = (actual_sense .eqv. (token > 0)) - end if - end select - - end do - - if (i_stack == 1) then - ! The one remaining logical on the stack indicates whether the particle is - ! in the cell. - in_cell = stack(i_stack) - else - ! This case occurs if there is no region specification since i_stack will - ! still be zero. - in_cell = .true. - end if - end function complex_cell_contains - !=============================================================================== ! CHECK_CELL_OVERLAP checks for overlapping cells at the current particle's ! position using cell_contains and the LocalCoord's built up by find_cell @@ -166,8 +77,8 @@ contains ! the particle should only be contained in one cell per level if (index_cell /= p % coord(j) % cell) then call fatal_error("Overlapping cells detected: " & - &// trim(to_str(cells(index_cell) % id)) // ", " & - &// trim(to_str(cells(p % coord(j) % cell) % id)) & + &// trim(to_str(cells(index_cell) % id())) // ", " & + &// trim(to_str(cells(p % coord(j) % cell) % id())) & &// " on universe " // trim(to_str(univ % id))) end if @@ -223,7 +134,7 @@ contains 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 + if (cells(i_cell) % universe() /= i_universe) cycle else i_cell = universes(i_universe) % cells(i) end if @@ -236,7 +147,7 @@ contains ! Show cell information on trace if (verbosity >= 10 .or. trace) then call write_message(" Entering cell " // trim(to_str(& - cells(i_cell) % id))) + cells(i_cell) % id()))) end if found = .true. @@ -246,7 +157,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 @@ -262,20 +173,20 @@ contains distribcell_index = c % distribcell_index offset = 0 do k = 1, p % n_coord - 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 + if (cells(p % coord(k) % cell) % type() == FILL_UNIVERSE) then + offset = offset + cells(p % coord(k) % cell) & + % offset(distribcell_index-1) + 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, & p % coord(k + 1) % lattice_y, & p % coord(k + 1) % lattice_z])) then - offset = offset + lattices(p % coord(k + 1) % lattice) % obj % & - offset(distribcell_index, & - p % coord(k + 1) % lattice_x, & + offset = offset + lattices(p % coord(k + 1) % lattice) % obj & + % offset(distribcell_index - 1, & + [p % coord(k + 1) % lattice_x, & p % coord(k + 1) % lattice_y, & - p % coord(k + 1) % lattice_z) + p % coord(k + 1) % lattice_z]) end if end if end do @@ -300,7 +211,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 @@ -311,7 +222,7 @@ contains ! Move particle to next level and set universe j = j + 1 p % n_coord = j - p % coord(j) % universe = c % fill + p % coord(j) % universe = c % fill() + 1 ! Apply translation if (allocated(c % translation)) then @@ -328,11 +239,11 @@ 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 - associate (lat => lattices(c % fill) % obj) + associate (lat => lattices(c % fill() + 1) % obj) ! Determine lattice indices i_xyz = lat % get_indices(p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw) @@ -341,7 +252,7 @@ contains p % coord(j + 1) % uvw = p % coord(j) % uvw ! set particle lattice indices - p % coord(j + 1) % lattice = c % fill + p % coord(j + 1) % lattice = c % fill() + 1 p % coord(j + 1) % lattice_x = i_xyz(1) p % coord(j + 1) % lattice_y = i_xyz(2) p % coord(j + 1) % lattice_z = i_xyz(3) @@ -350,18 +261,18 @@ contains if (lat % are_valid_indices(i_xyz)) then ! Particle is inside the lattice. p % coord(j + 1) % universe = & - lat % universes(i_xyz(1), i_xyz(2), i_xyz(3)) + lat % get([i_xyz(1), i_xyz(2), i_xyz(3)]) + 1 else ! Particle is outside the lattice. - if (lat % outer == NO_OUTER_UNIVERSE) then + if (lat % outer() == NO_OUTER_UNIVERSE) then call warning("Particle " // trim(to_str(p %id)) & - // " is outside lattice " // trim(to_str(lat % id)) & + // " is outside lattice " // trim(to_str(lat % id())) & // " but the lattice has no defined outer universe.") found = .false. return else - p % coord(j + 1) % universe = lat % outer + p % coord(j + 1) % universe = lat % outer() + 1 end if end if end associate @@ -396,7 +307,7 @@ contains lat => lattices(p % coord(j) % lattice) % obj if (verbosity >= 10 .or. trace) then - call write_message(" Crossing lattice " // trim(to_str(lat % id)) & + call write_message(" Crossing lattice " // trim(to_str(lat % id())) & &// ". Current position (" // trim(to_str(p % coord(j) % lattice_x)) & &// "," // trim(to_str(p % coord(j) % lattice_y)) // "," & &// trim(to_str(p % coord(j) % lattice_z)) // ")") @@ -428,7 +339,8 @@ contains else OUTSIDE_LAT ! Find cell in next lattice element - p % coord(j) % universe = lat % universes(i_xyz(1), i_xyz(2), i_xyz(3)) + p % coord(j) % universe = & + lat % get([i_xyz(1), i_xyz(2), i_xyz(3)]) + 1 call find_cell(p, found) if (.not. found) then @@ -465,26 +377,15 @@ contains integer, intent(out) :: lattice_translation(3) integer, intent(out) :: next_level - integer :: i ! index for surface in cell integer :: j - integer :: index_surf ! index in surfaces array (with sign) integer :: i_xyz(3) ! lattice indices integer :: level_surf_cross ! surface crossed on current level integer :: level_lat_trans(3) ! lattice translation on current level - real(8) :: x,y,z ! particle coordinates real(8) :: xyz_t(3) ! local particle coordinates - real(8) :: beta, gama ! skewed particle coordiantes - real(8) :: u,v,w ! particle directions - real(8) :: beta_dir ! skewed particle direction - real(8) :: gama_dir ! skewed particle direction - real(8) :: edge ! distance to oncoming edge - real(8) :: d ! evaluated distance real(8) :: d_lat ! distance to lattice boundary real(8) :: d_surf ! distance to surface - real(8) :: x0,y0,z0 ! coefficients for surface real(8) :: xyz_cross(3) ! coordinates at projected surface crossing real(8) :: surf_uvw(3) ! surface normal direction - logical :: coincident ! is particle on surface? type(Cell), pointer :: c class(Lattice), pointer :: lat @@ -502,35 +403,11 @@ contains ! get pointer to cell on this level c => cells(p % coord(j) % cell) - ! copy directional cosines - u = p % coord(j) % uvw(1) - v = p % coord(j) % uvw(2) - w = p % coord(j) % uvw(3) - ! ======================================================================= ! FIND MINIMUM DISTANCE TO SURFACE IN THIS CELL - SURFACE_LOOP: do i = 1, size(c % region) - index_surf = c % region(i) - coincident = (index_surf == p % surface) - - ! ignore this token if it corresponds to an operator rather than a - ! region. - index_surf = abs(index_surf) - if (index_surf >= OP_UNION) cycle - - ! Calculate distance to surface - d = surfaces(index_surf) % distance(p % coord(j) % xyz, & - p % coord(j) % uvw, logical(coincident, kind=C_BOOL)) - - ! Check if calculated distance is new minimum - if (d < d_surf) then - if (abs(d - d_surf)/d_surf >= FP_PRECISION) then - d_surf = d - level_surf_cross = -c % region(i) - end if - end if - end do SURFACE_LOOP + call c % distance(p % coord(j) % xyz, p % coord(j) % uvw, p % surface, & + d_surf, level_surf_cross) ! ======================================================================= ! FIND MINIMUM DISTANCE TO LATTICE SURFACES @@ -538,185 +415,22 @@ contains LAT_COORD: if (p % coord(j) % lattice /= NONE) then lat => lattices(p % coord(j) % lattice) % obj + i_xyz(1) = p % coord(j) % lattice_x + i_xyz(2) = p % coord(j) % lattice_y + i_xyz(3) = p % coord(j) % lattice_z + LAT_TYPE: select type(lat) type is (RectLattice) - ! copy local coordinates - x = p % coord(j) % xyz(1) - y = p % coord(j) % xyz(2) - z = p % coord(j) % xyz(3) - - ! determine oncoming edge - x0 = sign(lat % pitch(1) * HALF, u) - y0 = sign(lat % pitch(2) * HALF, v) - - ! left and right sides - if (abs(x - x0) < FP_PRECISION) then - d = INFINITY - elseif (u == ZERO) then - d = INFINITY - else - d = (x0 - x)/u - end if - - d_lat = d - if (u > 0) then - level_lat_trans(:) = [1, 0, 0] - else - level_lat_trans(:) = [-1, 0, 0] - end if - - ! front and back sides - if (abs(y - y0) < FP_PRECISION) then - d = INFINITY - elseif (v == ZERO) then - d = INFINITY - else - d = (y0 - y)/v - end if - - if (d < d_lat) then - d_lat = d - if (v > 0) then - level_lat_trans(:) = [0, 1, 0] - else - level_lat_trans(:) = [0, -1, 0] - end if - end if - - if (lat % is_3d) then - z0 = sign(lat % pitch(3) * HALF, w) - - ! top and bottom sides - if (abs(z - z0) < FP_PRECISION) then - d = INFINITY - elseif (w == ZERO) then - d = INFINITY - else - d = (z0 - z)/w - end if - - if (d < d_lat) then - d_lat = d - if (w > 0) then - level_lat_trans(:) = [0, 0, 1] - else - level_lat_trans(:) = [0, 0, -1] - end if - end if - end if + call lat % distance(p % coord(j) % xyz, p % coord(j) % uvw, & + i_xyz, d_lat, level_lat_trans) type is (HexLattice) LAT_TYPE - ! Copy local coordinates. - z = p % coord(j) % xyz(3) - i_xyz(1) = p % coord(j) % lattice_x - i_xyz(2) = p % coord(j) % lattice_y - i_xyz(3) = p % coord(j) % lattice_z - - ! Compute velocities along the hexagonal axes. - beta_dir = u*sqrt(THREE)/TWO + v/TWO - gama_dir = u*sqrt(THREE)/TWO - v/TWO - - ! Note that hexagonal lattice distance calculations are performed - ! using the particle's coordinates relative to the neighbor lattice - ! cells, not relative to the particle's current cell. This is done - ! because there is significant disagreement between neighboring cells - ! on where the lattice boundary is due to the worse finite precision - ! of hex lattices. - - ! Upper right and lower left sides. - edge = -sign(lat % pitch(1)/TWO, beta_dir) ! Oncoming edge - if (beta_dir > ZERO) then - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, 0, 0]) - else - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 0, 0]) - end if - beta = xyz_t(1)*sqrt(THREE)/TWO + xyz_t(2)/TWO - if (abs(beta - edge) < FP_PRECISION) then - d = INFINITY - else if (beta_dir == ZERO) then - d = INFINITY - else - d = (edge - beta)/beta_dir - end if - - d_lat = d - if (beta_dir > 0) then - level_lat_trans(:) = [1, 0, 0] - else - level_lat_trans(:) = [-1, 0, 0] - end if - - ! Lower right and upper left sides. - edge = -sign(lat % pitch(1)/TWO, gama_dir) ! Oncoming edge - if (gama_dir > ZERO) then - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[1, -1, 0]) - else - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[-1, 1, 0]) - end if - gama = xyz_t(1)*sqrt(THREE)/TWO - xyz_t(2)/TWO - if (abs(gama - edge) < FP_PRECISION) then - d = INFINITY - else if (gama_dir == ZERO) then - d = INFINITY - else - d = (edge - gama)/gama_dir - end if - - if (d < d_lat) then - d_lat = d - if (gama_dir > 0) then - level_lat_trans(:) = [1, -1, 0] - else - level_lat_trans(:) = [-1, 1, 0] - end if - end if - - ! Upper and lower sides. - edge = -sign(lat % pitch(1)/TWO, v) ! Oncoming edge - if (v > ZERO) then - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, 1, 0]) - else - xyz_t = lat % get_local_xyz(p % coord(j - 1) % xyz, i_xyz+[0, -1, 0]) - end if - if (abs(xyz_t(2) - edge) < FP_PRECISION) then - d = INFINITY - else if (v == ZERO) then - d = INFINITY - else - d = (edge - xyz_t(2))/v - end if - - if (d < d_lat) then - d_lat = d - if (v > 0) then - level_lat_trans(:) = [0, 1, 0] - else - level_lat_trans(:) = [0, -1, 0] - end if - end if - - ! Top and bottom sides. - if (lat % is_3d) then - z0 = sign(lat % pitch(2) * HALF, w) - - if (abs(z - z0) < FP_PRECISION) then - d = INFINITY - elseif (w == ZERO) then - d = INFINITY - else - d = (z0 - z)/w - end if - - if (d < d_lat) then - d_lat = d - if (w > 0) then - level_lat_trans(:) = [0, 0, 1] - else - level_lat_trans(:) = [0, 0, -1] - end if - end if - end if + xyz_t(1) = p % coord(j-1) % xyz(1) + xyz_t(2) = p % coord(j-1) % xyz(2) + xyz_t(3) = p % coord(j) % xyz(3) + call lat % distance(xyz_t, p % coord(j) % uvw, & + i_xyz, d_lat, level_lat_trans) end select LAT_TYPE if (d_lat < ZERO) then @@ -738,7 +452,7 @@ contains ! positive half-space were given in the region specification. Thus, we ! have to explicitly check which half-space the particle would be ! traveling into if the surface is crossed - if (.not. c % simple) then + if (.not. c % simple()) then xyz_cross(:) = p % coord(j) % xyz + d_surf*p % coord(j) % uvw call surfaces(abs(level_surf_cross)) % normal(xyz_cross, surf_uvw) if (dot_product(p % coord(j) % uvw, surf_uvw) > ZERO) then @@ -820,388 +534,4 @@ contains end subroutine neighbor_lists -!=============================================================================== -! CALC_OFFSETS calculates and stores the offsets in all fill cells. This -! routine is called once upon initialization. -!=============================================================================== - - subroutine calc_offsets(univ_id, map, univ, counts, found) - - integer, intent(in) :: univ_id ! target universe ID - integer, intent(in) :: map ! map index in vector of maps - type(Universe), intent(in) :: univ ! universe searching in - integer, intent(inout) :: counts(:,:) ! target count - logical, intent(inout) :: found(:,:) ! target found - - integer :: i ! index over cells - integer :: j, k, m ! indices in lattice - integer :: n ! number of cells to search - integer :: offset ! total offset for a given cell - integer :: cell_index ! index in cells array - type(Cell), pointer :: c ! pointer to current cell - type(Universe), pointer :: next_univ ! next universe to cycle through - class(Lattice), pointer :: lat ! pointer to current lattice - - n = size(univ % cells) - offset = 0 - - do i = 1, n - - cell_index = univ % cells(i) - - ! get pointer to cell - c => cells(cell_index) - - ! ==================================================================== - ! AT LOWEST UNIVERSE, TERMINATE SEARCH - if (c % type == FILL_MATERIAL) then - - ! ==================================================================== - ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - elseif (c % type == FILL_UNIVERSE) then - ! Set offset for the cell on this level - c % offset(map) = offset - - ! Count contents of this cell - next_univ => universes(c % fill) - offset = offset + count_target(next_univ, counts, found, univ_id, map) - - ! Move into the next universe - next_univ => universes(c % fill) - c => cells(cell_index) - - ! ==================================================================== - ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then - - ! Set current lattice - lat => lattices(c % fill) % obj - - select type (lat) - - type is (RectLattice) - - ! Loop over lattice coordinates - do m = 1, lat % n_cells(3) - do k = 1, lat % n_cells(2) - do j = 1, lat % n_cells(1) - lat % offset(map, j, k, m) = offset - next_univ => universes(lat % universes(j, k, m)) - offset = offset + & - count_target(next_univ, counts, found, univ_id, map) - 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 - ! This array location is never used - if (j + k < lat % n_rings + 1) then - cycle - ! This array location is never used - else if (j + k > 3*lat % n_rings - 1) then - cycle - else - lat % offset(map, j, k, m) = offset - next_univ => universes(lat % universes(j, k, m)) - offset = offset + & - count_target(next_univ, counts, found, univ_id, map) - end if - end do - end do - end do - end select - - end if - end do - - end subroutine calc_offsets - -!=============================================================================== -! COUNT_TARGET recursively totals the numbers of occurances of a given -! universe ID beginning with the universe given. -!=============================================================================== - - recursive function count_target(univ, counts, found, univ_id, map) result(count) - - type(Universe), intent(in) :: univ ! universe to search through - integer, intent(inout) :: counts(:,:) ! target count - logical, intent(inout) :: found(:,:) ! target found - integer, intent(in) :: univ_id ! target universe ID - integer, intent(in) :: map ! current map - - integer :: i ! index over cells - integer :: j, k, m ! indices in lattice - integer :: n ! number of cells to search - integer :: cell_index ! index in cells array - integer :: count ! number of times target located - type(Cell), pointer :: c ! pointer to current cell - type(Universe), pointer :: next_univ ! next univ to loop through - class(Lattice), pointer :: lat ! pointer to current lattice - - ! Don't research places already checked - if (found(universe_dict % get(univ % id), map)) then - count = counts(universe_dict % get(univ % id), map) - return - end if - - ! If this is the target, it can't contain itself. - ! Count = 1, then quit - if (univ % id == univ_id) then - count = 1 - counts(universe_dict % get(univ % id), map) = 1 - found(universe_dict % get(univ % id), map) = .true. - return - end if - - count = 0 - n = size(univ % cells) - - do i = 1, n - - cell_index = univ % cells(i) - - ! get pointer to cell - c => cells(cell_index) - - ! ==================================================================== - ! AT LOWEST UNIVERSE, TERMINATE SEARCH - if (c % type == FILL_MATERIAL) then - - ! ==================================================================== - ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - elseif (c % type == FILL_UNIVERSE) then - - next_univ => universes(c % fill) - - ! Found target - stop since target cannot contain itself - if (next_univ % id == univ_id) then - count = count + 1 - return - end if - - count = count + count_target(next_univ, counts, found, univ_id, map) - c => cells(cell_index) - - ! ==================================================================== - ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then - - ! Set current lattice - lat => lattices(c % fill) % obj - - select type (lat) - - type is (RectLattice) - - ! Loop over lattice coordinates - do m = 1, lat % n_cells(3) - do k = 1, lat % n_cells(2) - do j = 1, lat % n_cells(1) - next_univ => universes(lat % universes(j, k, m)) - - ! Found target - stop since target cannot contain itself - if (next_univ % id == univ_id) then - count = count + 1 - cycle - end if - - count = count + & - count_target(next_univ, counts, found, univ_id, map) - - 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 - ! This array location is never used - if (j + k < lat % n_rings + 1) then - cycle - ! This array location is never used - else if (j + k > 3*lat % n_rings - 1) then - cycle - else - next_univ => universes(lat % universes(j, k, m)) - - ! Found target - stop since target cannot contain itself - if (next_univ % id == univ_id) then - count = count + 1 - cycle - end if - - count = count + & - count_target(next_univ, counts, found, univ_id, map) - end if - end do - end do - end do - - end select - - end if - end do - - counts(universe_dict % get(univ % id), map) = count - found(universe_dict % get(univ % id), map) = .true. - - 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 m = 1, lat % n_cells(3) - do k = 1, lat % n_cells(2) - do j = 1, lat % n_cells(1) - call count_instance(universes(lat % universes(j, k, m))) - 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 - ! This array location is never used - if (j + k < lat % n_rings + 1) then - cycle - ! This array location is never used - else if (j + k > 3*lat % n_rings - 1) then - cycle - else - call count_instance(universes(lat % universes(j, k, m))) - 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 -!=============================================================================== - - recursive function maximum_levels(univ) result(levels) - - type(Universe), intent(in) :: univ ! universe to search through - integer :: levels ! maximum number of levels for this universe - - integer :: i ! index over cells - integer :: j, k, m ! indices in lattice - integer :: levels_below ! max levels below this universe - type(Cell), pointer :: c ! pointer to current cell - type(Universe), pointer :: next_univ ! next universe to loop through - class(Lattice), pointer :: lat ! pointer to current lattice - - levels_below = 0 - do i = 1, size(univ % cells) - c => cells(univ % cells(i)) - - ! ==================================================================== - ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - 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 - - ! Set current lattice - lat => lattices(c % fill) % obj - - select type (lat) - - type is (RectLattice) - - ! Loop over lattice coordinates - do m = 1, lat % n_cells(3) - do k = 1, lat % n_cells(2) - do j = 1, lat % n_cells(1) - next_univ => universes(lat % universes(j, k, m)) - levels_below = max(levels_below, maximum_levels(next_univ)) - 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 - ! This array location is never used - if (j + k < lat % n_rings + 1) then - cycle - ! This array location is never used - else if (j + k > 3*lat % n_rings - 1) then - cycle - else - next_univ => universes(lat % universes(j, k, m)) - levels_below = max(levels_below, maximum_levels(next_univ)) - end if - end do - end do - end do - - end select - - end if - end do - - levels = 1 + levels_below - - end function maximum_levels - end module geometry diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp new file mode 100644 index 0000000000..3f36cf299a --- /dev/null +++ b/src/geometry_aux.cpp @@ -0,0 +1,341 @@ +#include "geometry_aux.h" + +#include // for std::max +#include +#include + +#include "cell.h" +#include "constants.h" +#include "error.h" +#include "lattice.h" + + +namespace openmc { + +//============================================================================== + +void +adjust_indices_c() +{ + // Adjust material/fill idices. + for (Cell *c : global_cells) { + if (c->material[0] == C_NONE) { + int32_t id = c->fill; + auto search_univ = universe_map.find(id); + auto search_lat = lattice_map.find(id); + if (search_univ != universe_map.end()) { + c->type = FILL_UNIVERSE; + c->fill = search_univ->second; + } else if (search_lat != lattice_map.end()) { + c->type = FILL_LATTICE; + c->fill = search_lat->second; + } 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 : global_cells) { + auto search = universe_map.find(c->universe); + if (search != universe_map.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(); + } +} + +//============================================================================== + +int32_t +find_root_universe() +{ + // Find all the universes listed as a cell fill. + std::unordered_set fill_univ_ids; + for (Cell *c : global_cells) { + 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 < global_universes.size(); i++) { + auto search = fill_univ_ids.find(global_universes[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; +} + +//============================================================================== + +void +allocate_offset_tables(int n_maps) +{ + for (Cell *c : global_cells) { + if (c->type != FILL_MATERIAL) { + c->offset.resize(n_maps, C_NONE); + } + } + + for (Lattice *lat : lattices_c) { + lat->allocate_offset_table(n_maps); + } +} + +//============================================================================== + +void +count_cell_instances(int32_t univ_indx) +{ + for (int32_t cell_indx : global_universes[univ_indx]->cells) { + Cell &c = *global_cells[cell_indx]; + ++c.n_instances; + + if (c.type == FILL_UNIVERSE) { + // This cell contains another universe. Recurse into that universe. + count_cell_instances(c.fill); + + } else if (c.type == FILL_LATTICE) { + // This cell contains a lattice. Recurse into the lattice universes. + Lattice &lat = *lattices_c[c.fill]; + for (auto it = lat.begin(); it != lat.end(); ++it) { + count_cell_instances(*it); + } + } + } +} + +//============================================================================== + +int +count_universe_instances(int32_t search_univ, int32_t target_univ_id) +{ + // If this is the target, it can't contain itself. + if (global_universes[search_univ]->id == target_univ_id) { + return 1; + } + + int count {0}; + for (int32_t cell_indx : global_universes[search_univ]->cells) { + Cell &c = *global_cells[cell_indx]; + + if (c.type == FILL_UNIVERSE) { + int32_t next_univ = c.fill; + count += count_universe_instances(next_univ, target_univ_id); + + } else if (c.type == FILL_LATTICE) { + Lattice &lat = *lattices_c[c.fill]; + for (auto it = lat.begin(); it != lat.end(); ++it) { + int32_t next_univ = *it; + count += count_universe_instances(next_univ, target_univ_id); + } + } + } + + return count; +} + +//============================================================================== + +void +fill_offset_tables(int32_t target_univ_id, int map) +{ + for (Universe *univ : global_universes) { + int32_t offset {0}; // TODO: is this a bug? It matches F90 implementation. + for (int32_t cell_indx : univ->cells) { + Cell &c = *global_cells[cell_indx]; + + if (c.type == FILL_UNIVERSE) { + c.offset[map] = offset; + int32_t search_univ = c.fill; + offset += count_universe_instances(search_univ, target_univ_id); + + } else if (c.type == FILL_LATTICE) { + Lattice &lat = *lattices_c[c.fill]; + offset = lat.fill_offset_table(offset, target_univ_id, map); + } + } + } +} + +//============================================================================== + +std::string +distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset, + const Universe &search_univ, int32_t offset) +{ + std::stringstream path; + + path << "u" << search_univ.id << "->"; + + // Check to see if this universe directly contains the target cell. If so, + // write to the path and return. + for (int32_t cell_indx : search_univ.cells) { + if ((cell_indx == target_cell) && (offset == target_offset)) { + Cell &c = *global_cells[cell_indx]; + path << "c" << c.id; + return path.str(); + } + } + + // The target must be further down the geometry tree and contained in a fill + // cell or lattice cell in this universe. Find which cell contains the + // target. + std::vector::const_reverse_iterator cell_it + {search_univ.cells.crbegin()}; + for (; cell_it != search_univ.cells.crend(); ++cell_it) { + Cell &c = *global_cells[*cell_it]; + + // Material cells don't contain other cells so ignore them. + if (c.type != FILL_MATERIAL) { + int32_t temp_offset; + if (c.type == FILL_UNIVERSE) { + temp_offset = offset + c.offset[map]; + } else { + Lattice &lat = *lattices_c[c.fill]; + int32_t indx = lat.universes.size()*map + lat.begin().indx; + temp_offset = offset + lat.offsets[indx]; + } + + // The desired cell is the first cell that gives an offset smaller or + // equal to the target offset. + if (temp_offset <= target_offset) break; + } + } + + // Add the cell to the path string. + Cell &c = *global_cells[*cell_it]; + path << "c" << c.id << "->"; + + if (c.type == FILL_UNIVERSE) { + // Recurse into the fill cell. + offset += c.offset[map]; + path << distribcell_path_inner(target_cell, map, target_offset, + *global_universes[c.fill], offset); + return path.str(); + } else { + // Recurse into the lattice cell. + Lattice &lat = *lattices_c[c.fill]; + path << "l" << lat.id; + for (ReverseLatticeIter it = lat.rbegin(); it != lat.rend(); ++it) { + int32_t indx = lat.universes.size()*map + it.indx; + int32_t temp_offset = offset + lat.offsets[indx]; + if (temp_offset <= target_offset) { + offset = temp_offset; + path << "(" << lat.index_to_string(it.indx) << ")->"; + path << distribcell_path_inner(target_cell, map, target_offset, + *global_universes[*it], offset); + return path.str(); + } + } + } +} + +//============================================================================== + +int +distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset, + int32_t root_univ) +{ + Universe &root = *global_universes[root_univ]; + std::string path_ {distribcell_path_inner(target_cell, map, target_offset, + root, 0)}; + return path_.size() + 1; +} + +//============================================================================== + +void +distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset, + int32_t root_univ, char *path) +{ + Universe &root = *global_universes[root_univ]; + std::string path_ {distribcell_path_inner(target_cell, map, target_offset, + root, 0)}; + path_.copy(path, path_.size()); + path[path_.size()] = '\0'; +} + +//============================================================================== + +int +maximum_levels(int32_t univ) +{ + int levels_below {0}; + + for (int32_t cell_indx : global_universes[univ]->cells) { + Cell &c = *global_cells[cell_indx]; + if (c.type == FILL_UNIVERSE) { + int32_t next_univ = c.fill; + levels_below = std::max(levels_below, maximum_levels(next_univ)); + } else if (c.type == FILL_LATTICE) { + Lattice &lat = *lattices_c[c.fill]; + for (auto it = lat.begin(); it != lat.end(); ++it) { + int32_t next_univ = *it; + levels_below = std::max(levels_below, maximum_levels(next_univ)); + } + } + } + + ++levels_below; + return levels_below; +} + +//============================================================================== + +void +free_memory_geometry_c() +{ + for (Cell *c : global_cells) {delete c;} + global_cells.clear(); + cell_map.clear(); + n_cells = 0; + + for (Universe *u : global_universes) {delete u;} + global_universes.clear(); + universe_map.clear(); + + for (Lattice *lat : lattices_c) {delete lat;} + lattices_c.clear(); + lattice_map.clear(); +} + +} // namespace openmc diff --git a/src/geometry_aux.h b/src/geometry_aux.h new file mode 100644 index 0000000000..b998ad32c9 --- /dev/null +++ b/src/geometry_aux.h @@ -0,0 +1,112 @@ +//! \file geometry_aux.h +//! Auxilary functions for geometry initialization and general data handling. + +#ifndef GEOMETRY_AUX_H +#define GEOMETRY_AUX_H + +#include + + +namespace openmc { + +//============================================================================== +//! Replace Universe, Lattice, and Material IDs with indices. +//============================================================================== + +extern "C" void adjust_indices_c(); + +//============================================================================== +//! Figure out which Universe is the root universe. +//! +//! This function looks for a universe that is not listed in a Cell::fill or in +//! a Lattice. +//! @return The index of the root universe. +//============================================================================== + +extern "C" int32_t find_root_universe(); + +//============================================================================== +//! Allocate storage in Lattice and Cell objects for distribcell offset tables. +//============================================================================== + +extern "C" void allocate_offset_tables(int n_maps); + +//============================================================================== +//! Recursively search through the geometry and count cell instances. +//! +//! This function will update the Cell::n_instances value for each cell in the +//! geometry. +//! @param univ_indx The index of the universe to begin searching from (probably +//! the root universe). +//============================================================================== + +extern "C" void count_cell_instances(int32_t univ_indx); + +//============================================================================== +//! Recursively search through universes and count universe instances. +//! @param search_univ The index of the universe to begin searching from. +//! @param target_univ_id The ID of the universe to be counted. +//! @return The number of instances of target_univ_id in the geometry tree under +//! search_univ. +//============================================================================== + +extern "C" int +count_universe_instances(int32_t search_univ, int32_t target_univ_id); + +//============================================================================== +//! Populate Cell and Lattice distribcell offset tables. +//! @param target_univ_id The ID of the universe to be counted. +//! @param map The index of the distribcell map that defines the offsets for the +//! target universe. +//============================================================================== + +extern "C" void fill_offset_tables(int32_t target_univ_id, int map); + +//============================================================================== +//! Find the length necessary for a string to contain a distribcell path. +//! @param target_cell The index of the Cell in the global Cell array. +//! @param map The index of the distribcell mapping corresponding to the target +//! cell. +//! @param target_offset An instance number for a distributed cell. +//! @param root_univ The index of the root Universe in the global Universe +//! array. +//! @return The size of a character array needed to fit the distribcell path. +//============================================================================== + +extern "C" int +distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset, + int32_t root_univ); + +//============================================================================== +//! Build a character array representing the path to a distribcell instance. +//! @param target_cell The index of the Cell in the global Cell array. +//! @param map The index of the distribcell mapping corresponding to the target +//! cell. +//! @param target_offset An instance number for a distributed cell. +//! @param root_univ The index of the root Universe in the global Universe +//! array. +//! @param[out] path The unique traversal through the geometry tree that leads +//! to the desired instance of the target cell. +//============================================================================== + +extern "C" void +distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset, + int32_t root_univ, char *path); + +//============================================================================== +//! Determine the maximum number of nested coordinate levels in the geometry. +//! @param univ The index of the universe to begin seraching from (probably the +//! root universe). +//! @return The number of coordinate levels. +//============================================================================== + +extern "C" int maximum_levels(int32_t univ); + +//============================================================================== +//! Deallocates global vectors and maps for cells, universes, and lattices. +//============================================================================== + +extern "C" void free_memory_geometry_c(); + +} // namespace openmc +#endif // GEOMETRY_AUX_H diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index f6206adf99..758fc54da7 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -6,6 +6,7 @@ module geometry_header use constants, only: HALF, TWO, THREE, INFINITY, K_BOLTZMANN, & MATERIAL_VOID, NONE use dict_header, only: DictCharInt, DictIntInt + use hdf5_interface, only: HID_T use material_header, only: Material, materials, material_dict, n_materials use nuclide_header use sab_header @@ -14,13 +15,192 @@ module geometry_header implicit none + interface + function cell_pointer_c(cell_ind) bind(C, name='cell_pointer') result(ptr) + import C_PTR, C_INT32_T + integer(C_INT32_T), intent(in), value :: cell_ind + type(C_PTR) :: ptr + end function cell_pointer_c + + function cell_id_c(cell_ptr) bind(C, name='cell_id') result(id) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT32_T) :: id + end function cell_id_c + + subroutine cell_set_id_c(cell_ptr, id) bind(C, name='cell_set_id') + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + 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 + 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 + 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 + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT32_T) :: universe + end function cell_universe_c + + subroutine cell_set_universe_c(cell_ptr, universe) & + bind(C, name='cell_set_universe') + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT32_T), intent(in), value :: universe + end subroutine cell_set_universe_c + + function cell_fill_c(cell_ptr) bind(C, name="cell_fill") result(fill) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT32_T) :: fill + end function cell_fill_c + + function cell_fill_ptr(cell_ptr) bind(C) result(fill_ptr) + import C_PTR + type(C_PTR), intent(in), value :: cell_ptr + type(C_PTR) :: fill_ptr + end function cell_fill_ptr + + function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') & + result(n_instances) + import C_PTR, C_INT32_T + 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 + type(C_PTR), intent(in), value :: cell_ptr + logical(C_BOOL) :: simple + end function cell_simple_c + + subroutine cell_distance_c(cell_ptr, xyz, uvw, on_surface, min_dist, & + i_surf) bind(C, name="cell_distance") + import C_PTR, C_INT32_T, C_DOUBLE + 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 + real(C_DOUBLE), intent(out) :: min_dist + integer(C_INT32_T), intent(out) :: i_surf + end subroutine cell_distance_c + + function cell_offset_c(cell_ptr, map) bind(C, name="cell_offset") & + result(offset) + import C_PTR, C_INT, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT), intent(in), value :: map + integer(C_INT32_T) :: offset + end function cell_offset_c + + subroutine cell_to_hdf5_c(cell_ptr, group) bind(C, name='cell_to_hdf5') + import HID_T, C_PTR + type(C_PTR), intent(in), value :: cell_ptr + integer(HID_T), intent(in), value :: group + end subroutine cell_to_hdf5_c + + function lattice_pointer_c(lat_ind) bind(C, name='lattice_pointer') & + result(ptr) + import C_PTR, C_INT32_T + integer(C_INT32_T), intent(in), value :: lat_ind + type(C_PTR) :: ptr + end function lattice_pointer_c + + function lattice_id_c(lat_ptr) bind(C, name='lattice_id') result(id) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT32_T) :: id + end function lattice_id_c + + function lattice_are_valid_indices_c(lat_ptr, i_xyz) & + bind(C, name='lattice_are_valid_indices') result (is_valid) + import C_PTR, C_INT, C_BOOL + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT), intent(in) :: i_xyz(3) + logical(C_BOOL) :: is_valid + end function lattice_are_valid_indices_c + + subroutine lattice_distance_c(lat_ptr, xyz, uvw, i_xyz, d, lattice_trans) & + bind(C, name='lattice_distance') + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), intent(in), value :: lat_ptr + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + integer(C_INT), intent(in) :: i_xyz(3) + real(C_DOUBLE), intent(out) :: d + integer(C_INT), intent(out) :: lattice_trans(3) + end subroutine lattice_distance_c + + subroutine lattice_get_indices_c(lat_ptr, xyz, i_xyz) & + bind(C, name='lattice_get_indices') + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), intent(in), value :: lat_ptr + real(C_DOUBLE), intent(in) :: xyz(3) + integer(C_INT), intent(out) :: i_xyz(3) + end subroutine lattice_get_indices_c + + subroutine lattice_get_local_xyz_c(lat_ptr, global_xyz, i_xyz, local_xyz) & + bind(C, name='lattice_get_local_xyz') + import C_PTR, C_INT, C_DOUBLE + type(C_PTR), intent(in), value :: lat_ptr + real(C_DOUBLE), intent(in) :: global_xyz(3) + integer(C_INT), intent(in) :: i_xyz(3) + real(C_DOUBLE), intent(out) :: local_xyz(3) + end subroutine lattice_get_local_xyz_c + + subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5') + import HID_T, C_PTR + type(C_PTR), intent(in), value :: lat_ptr + integer(HID_T), intent(in), value :: group + end subroutine lattice_to_hdf5_c + + function lattice_offset_c(lat_ptr, map, i_xyz) & + bind(C, name='lattice_offset') result(offset) + import C_PTR, C_INT, C_INT32_T + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT), intent(in), value :: map + integer(C_INT), intent(in) :: i_xyz(3) + integer(C_INT32_T) :: offset + end function lattice_offset_c + + function lattice_outer_c(lat_ptr) bind(C, name='lattice_outer') & + result(outer) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT32_T) :: outer + end function lattice_outer_c + + function lattice_universe_c(lat_ptr, i_xyz) & + bind(C, name='lattice_universe') result(univ) + import C_PTR, C_INT32_t, C_INT + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT), intent(in) :: i_xyz(3) + integer(C_INT32_T) :: univ + end function lattice_universe_c + + subroutine extend_cells_c(n) bind(C) + import C_INT32_t + integer(C_INT32_T), intent(in), value :: n + end subroutine extend_cells_c + end interface + !=============================================================================== ! UNIVERSE defines a geometry that fills all phase space !=============================================================================== type Universe integer :: id ! Unique ID - integer :: type ! Type integer, allocatable :: cells(:) ! List of cells within real(8) :: x0 ! Translation in x-coordinate real(8) :: y0 ! Translation in y-coordinate @@ -32,68 +212,24 @@ module geometry_header !=============================================================================== type, abstract :: Lattice - integer :: id ! Universe number for lattice - character(len=104) :: name = "" ! User-defined name - real(8), allocatable :: pitch(:) ! Pitch along each axis - integer, allocatable :: universes(:,:,:) ! Specified universes - integer :: outside ! Material to fill area outside - integer :: outer ! universe to tile outside the lat - logical :: is_3d ! Lattice has cells on z axis - integer, allocatable :: offset(:,:,:,:) ! Distribcell offsets + type(C_PTR) :: ptr contains - procedure(lattice_are_valid_indices_), deferred :: are_valid_indices - procedure(lattice_get_indices_), deferred :: get_indices - procedure(lattice_get_local_xyz_), deferred :: get_local_xyz + procedure :: id => lattice_id + procedure :: are_valid_indices => lattice_are_valid_indices + procedure :: distance => lattice_distance + procedure :: get => lattice_get + procedure :: get_indices => lattice_get_indices + procedure :: get_local_xyz => lattice_get_local_xyz + procedure :: offset => lattice_offset + procedure :: outer => lattice_outer + procedure :: to_hdf5 => lattice_to_hdf5 end type Lattice - abstract interface - -!=============================================================================== -! ARE_VALID_INDICES returns .true. if the given lattice indices fit within the -! bounds of the lattice. Returns false otherwise. - - function lattice_are_valid_indices_(this, i_xyz) result(is_valid) - import Lattice - class(Lattice), intent(in) :: this - integer, intent(in) :: i_xyz(3) - logical :: is_valid - end function lattice_are_valid_indices_ - -!=============================================================================== -! GET_INDICES returns the indices in a lattice for the given global xyz. - - function lattice_get_indices_(this, global_xyz) result(i_xyz) - import Lattice - class(Lattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer :: i_xyz(3) - end function lattice_get_indices_ - -!=============================================================================== -! GET_LOCAL_XYZ returns the translated local version of the given global xyz. - - function lattice_get_local_xyz_(this, global_xyz, i_xyz) result(local_xyz) - import Lattice - class(Lattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer, intent(in) :: i_xyz(3) - real(8) :: local_xyz(3) - end function lattice_get_local_xyz_ - end interface - !=============================================================================== ! RECTLATTICE extends LATTICE for rectilinear arrays. !=============================================================================== type, extends(Lattice) :: RectLattice - integer :: n_cells(3) ! Number of cells along each axis - real(8), allocatable :: lower_left(:) ! Global lower-left corner of lat - - contains - - procedure :: are_valid_indices => valid_inds_rect - procedure :: get_indices => get_inds_rect - procedure :: get_local_xyz => get_local_rect end type RectLattice !=============================================================================== @@ -101,15 +237,6 @@ module geometry_header !=============================================================================== type, extends(Lattice) :: HexLattice - integer :: n_rings ! Number of radial ring cell positoins - integer :: n_axial ! Number of axial cell positions - real(8), allocatable :: center(:) ! Global center of lattice - - contains - - procedure :: are_valid_indices => valid_inds_hex - procedure :: get_indices => get_inds_hex - procedure :: get_local_xyz => get_local_hex end type HexLattice !=============================================================================== @@ -125,25 +252,13 @@ module geometry_header !=============================================================================== type Cell - integer :: id ! Unique ID - character(len=104) :: name = "" ! User-defined name - integer :: type ! Type of cell (normal, universe, - ! lattice) - integer :: universe ! universe # this cell is in - integer :: fill ! universe # filling this cell - integer :: instances ! number of instances of this cell in - ! the geom + type(C_PTR) :: ptr + integer, allocatable :: material(:) ! Material within cell. Multiple ! materials for distribcell ! instances. 0 signifies a universe - integer, allocatable :: offset(:) ! Distribcell offset for tally - ! counter integer, allocatable :: region(:) ! Definition of spatial region as ! Boolean expression of half-spaces - integer, allocatable :: rpn(:) ! Reverse Polish notation for region - ! expression - logical :: simple ! Is the region simple (intersections - ! only) integer :: distribcell_index ! Index corresponding to this cell in ! distribcell arrays real(8), allocatable :: sqrtkT(:) ! Square root of k_Boltzmann * @@ -154,6 +269,22 @@ module geometry_header real(8), allocatable :: translation(:) real(8), allocatable :: rotation(:) real(8), allocatable :: rotation_matrix(:,:) + + contains + + 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 :: fill => cell_fill + procedure :: n_instances => cell_n_instances + procedure :: simple => cell_simple + procedure :: distance => cell_distance + procedure :: offset => cell_offset + procedure :: to_hdf5 => cell_to_hdf5 + end type Cell ! array index of the root universe @@ -161,7 +292,6 @@ module geometry_header integer(C_INT32_T), bind(C) :: n_cells ! # of cells integer(C_INT32_T), bind(C) :: n_universes ! # of universes - integer(C_INT32_T), bind(C) :: n_lattices ! # of lattices type(Cell), allocatable, target :: cells(:) type(Universe), allocatable, target :: universes(:) @@ -174,172 +304,150 @@ module geometry_header contains -!=============================================================================== + function lattice_id(this) result(id) + class(Lattice), intent(in) :: this + integer(C_INT32_T) :: id + id = lattice_id_c(this % ptr) + end function lattice_id - function valid_inds_rect(this, i_xyz) result(is_valid) - class(RectLattice), intent(in) :: this - integer, intent(in) :: i_xyz(3) - logical :: is_valid + function lattice_are_valid_indices(this, i_xyz) result (is_valid) + class(Lattice), intent(in) :: this + integer(C_INT), intent(in) :: i_xyz(3) + logical(C_BOOL) :: is_valid + is_valid = lattice_are_valid_indices_c(this % ptr, i_xyz) + end function lattice_are_valid_indices - is_valid = all(i_xyz > 0 .and. i_xyz <= this % n_cells) - end function valid_inds_rect + subroutine lattice_distance(this, xyz, uvw, i_xyz, d, lattice_trans) + class(Lattice), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + integer(C_INT), intent(in) :: i_xyz(3) + real(C_DOUBLE), intent(out) :: d + integer(C_INT), intent(out) :: lattice_trans(3) + call lattice_distance_c(this % ptr, xyz, uvw, i_xyz, d, lattice_trans) + end subroutine lattice_distance + + function lattice_get(this, i_xyz) result(univ) + class(Lattice), intent(in) :: this + integer(C_INT), intent(in) :: i_xyz(3) + integer(C_INT32_T) :: univ + univ = lattice_universe_c(this % ptr, i_xyz) + end function lattice_get + + function lattice_get_indices(this, xyz) result(i_xyz) + class(Lattice), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3) + integer(C_INT) :: i_xyz(3) + call lattice_get_indices_c(this % ptr, xyz, i_xyz) + end function lattice_get_indices + + function lattice_get_local_xyz(this, global_xyz, i_xyz) & + result(local_xyz) + class(Lattice), intent(in) :: this + real(C_DOUBLE), intent(in) :: global_xyz(3) + integer(C_INT), intent(in) :: i_xyz(3) + real(C_DOUBLE) :: local_xyz(3) + call lattice_get_local_xyz_c(this % ptr, global_xyz, i_xyz, local_xyz) + end function lattice_get_local_xyz + + function lattice_offset(this, map, i_xyz) result(offset) + class(Lattice), intent(in) :: this + integer(C_INT), intent(in) :: map + integer(C_INT), intent(in) :: i_xyz(3) + integer(C_INT32_T) :: offset + offset = lattice_offset_c(this % ptr, map, i_xyz) + end function lattice_offset + + function lattice_outer(this) result(outer) + class(Lattice), intent(in) :: this + integer(C_INT32_T) :: outer + outer = lattice_outer_c(this % ptr) + end function lattice_outer + + subroutine lattice_to_hdf5(this, group) + class(Lattice), intent(in) :: this + integer(HID_T), intent(in) :: group + call lattice_to_hdf5_c(this % ptr, group) + end subroutine lattice_to_hdf5 !=============================================================================== - function valid_inds_hex(this, i_xyz) result(is_valid) - class(HexLattice), intent(in) :: this - integer, intent(in) :: i_xyz(3) - logical :: is_valid + function cell_id(this) result(id) + class(Cell), intent(in) :: this + integer(C_INT32_T) :: id + id = cell_id_c(this % ptr) + end function cell_id - is_valid = (all(i_xyz > 0) .and. & - &i_xyz(1) < 2*this % n_rings .and. & - &i_xyz(2) < 2*this % n_rings .and. & - &i_xyz(1) + i_xyz(2) > this % n_rings .and. & - &i_xyz(1) + i_xyz(2) < 3*this % n_rings .and. & - &i_xyz(3) <= this % n_axial) - end function valid_inds_hex + subroutine cell_set_id(this, id) + class(Cell), intent(in) :: this + integer(C_INT32_T), intent(in) :: id + 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 - function get_inds_rect(this, global_xyz) result(i_xyz) - class(RectLattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer :: i_xyz(3) + 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 - real(8) :: xyz(3) ! global_xyz alias + function cell_universe(this) result(universe) + class(Cell), intent(in) :: this + integer(C_INT32_T) :: universe + universe = cell_universe_c(this % ptr) + end function cell_universe - xyz = global_xyz + subroutine cell_set_universe(this, universe) + class(Cell), intent(in) :: this + integer(C_INT32_T), intent(in) :: universe + call cell_set_universe_c(this % ptr, universe) + end subroutine cell_set_universe - i_xyz(1) = ceiling((xyz(1) - this % lower_left(1))/this % pitch(1)) - i_xyz(2) = ceiling((xyz(2) - this % lower_left(2))/this % pitch(2)) - if (this % is_3d) then - i_xyz(3) = ceiling((xyz(3) - this % lower_left(3))/this % pitch(3)) - else - i_xyz(3) = 1 - end if - end function get_inds_rect + function cell_fill(this) result(fill) + class(Cell), intent(in) :: this + integer(C_INT32_T) :: fill + fill = cell_fill_c(this % ptr) + end function cell_fill -!=============================================================================== + 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 get_inds_hex(this, global_xyz) result(i_xyz) - class(HexLattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer :: i_xyz(3) + function cell_simple(this) result(simple) + class(Cell), intent(in) :: this + logical(C_BOOL) :: simple + simple = cell_simple_c(this % ptr) + end function cell_simple - real(8) :: xyz(3) ! global xyz relative to the center - real(8) :: alpha ! Skewed coord axis - real(8) :: xyz_t(3) ! Local xyz - real(8) :: d, d_min ! Squared distance from cell centers - integer :: i, j, k ! Iterators - integer :: k_min ! Minimum distance index + subroutine cell_distance(this, xyz, uvw, on_surface, min_dist, i_surf) + class(Cell), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + integer(C_INT32_T), intent(in) :: on_surface + real(C_DOUBLE), intent(out) :: min_dist + integer(C_INT32_T), intent(out) :: i_surf + call cell_distance_c(this % ptr, xyz, uvw, on_surface, min_dist, i_surf) + end subroutine cell_distance - xyz(1) = global_xyz(1) - this % center(1) - xyz(2) = global_xyz(2) - this % center(2) + function cell_offset(this, map) result(offset) + class(Cell), intent(in) :: this + integer(C_INT), intent(in) :: map + integer(C_INT32_T) :: offset + offset = cell_offset_c(this % ptr, map) + end function cell_offset - ! Index z direction. - if (this % is_3d) then - xyz(3) = global_xyz(3) - this % center(3) - i_xyz(3) = ceiling(xyz(3)/this % pitch(2) + HALF*this % n_axial) - else - xyz(3) = global_xyz(3) - i_xyz(3) = 1 - end if - - ! Convert coordinates into skewed bases. The (x, alpha) basis is used to - ! find the index of the global coordinates to within 4 cells. - alpha = xyz(2) - xyz(1) / sqrt(THREE) - i_xyz(1) = floor(xyz(1) / (sqrt(THREE) / TWO * this % pitch(1))) - i_xyz(2) = floor(alpha / this % pitch(1)) - - ! Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but - ! the array is offset so that the indices never go below 1). - i_xyz(1) = i_xyz(1) + this % n_rings - i_xyz(2) = i_xyz(2) + this % n_rings - - ! Calculate the (squared) distance between the particle and the centers of - ! the four possible cells. Regular hexagonal tiles form a centroidal - ! Voronoi tessellation so the global xyz should be in the hexagonal cell - ! that it is closest to the center of. This method is used over a - ! method that uses the remainders of the floor divisions above because it - ! provides better finite precision performance. Squared distances are - ! used becasue they are more computationally efficient than normal - ! distances. - k = 1 - d_min = INFINITY - do i = 0, 1 - do j = 0, 1 - xyz_t = this % get_local_xyz(global_xyz, i_xyz + [j, i, 0]) - d = xyz_t(1)**2 + xyz_t(2)**2 - if (d < d_min) then - d_min = d - k_min = k - end if - k = k + 1 - end do - end do - - ! Select the minimum squared distance which corresponds to the cell the - ! coordinates are in. - if (k_min == 2) then - i_xyz(1) = i_xyz(1) + 1 - else if (k_min == 3) then - i_xyz(2) = i_xyz(2) + 1 - else if (k_min == 4) then - i_xyz(1) = i_xyz(1) + 1 - i_xyz(2) = i_xyz(2) + 1 - end if - end function get_inds_hex - -!=============================================================================== - - function get_local_rect(this, global_xyz, i_xyz) result(local_xyz) - class(RectLattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer, intent(in) :: i_xyz(3) - real(8) :: local_xyz(3) - - real(8) :: xyz(3) ! global_xyz alias - - xyz = global_xyz - - local_xyz(1) = xyz(1) - (this % lower_left(1) + & - (i_xyz(1) - HALF)*this % pitch(1)) - local_xyz(2) = xyz(2) - (this % lower_left(2) + & - (i_xyz(2) - HALF)*this % pitch(2)) - if (this % is_3d) then - local_xyz(3) = xyz(3) - (this % lower_left(3) + & - (i_xyz(3) - HALF)*this % pitch(3)) - else - local_xyz(3) = xyz(3) - end if - end function get_local_rect - -!=============================================================================== - - function get_local_hex(this, global_xyz, i_xyz) result(local_xyz) - class(HexLattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer, intent(in) :: i_xyz(3) - real(8) :: local_xyz(3) - - real(8) :: xyz(3) ! global_xyz alias - - xyz = global_xyz - - ! x_l = x_g - (center + pitch_x*cos(30)*index_x) - local_xyz(1) = xyz(1) - (this % center(1) + & - sqrt(THREE) / TWO * (i_xyz(1) - this % n_rings) * this % pitch(1)) - ! y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) - local_xyz(2) = xyz(2) - (this % center(2) + & - (i_xyz(2) - this % n_rings) * this % pitch(1) + & - (i_xyz(1) - this % n_rings) * this % pitch(1) / TWO) - if (this % is_3d) then - local_xyz(3) = xyz(3) - this % center(3) & - + (HALF*this % n_axial - i_xyz(3) + HALF) * this % pitch(2) - else - local_xyz(3) = xyz(3) - end if - end function get_local_hex + subroutine cell_to_hdf5(this, group) + class(Cell), intent(in) :: this + integer(HID_T), intent(in) :: group + call cell_to_hdf5_c(this % ptr, group) + end subroutine cell_to_hdf5 !=============================================================================== ! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table @@ -408,10 +516,15 @@ contains !=============================================================================== subroutine free_memory_geometry() + interface + subroutine free_memory_geometry_c() bind(C) + end subroutine free_memory_geometry_c + end interface + + call free_memory_geometry_c() n_cells = 0 n_universes = 0 - n_lattices = 0 if (allocated(cells)) deallocate(cells) if (allocated(universes)) deallocate(universes) @@ -432,6 +545,7 @@ contains integer(C_INT32_T), value, intent(in) :: n integer(C_INT32_T), optional, intent(out) :: index_start integer(C_INT32_T), optional, intent(out) :: index_end + integer(C_INT32_T) :: i integer(C_INT) :: err type(Cell), allocatable :: temp(:) ! temporary cells array @@ -453,7 +567,12 @@ contains ! Return indices in cells array if (present(index_start)) index_start = n_cells + 1 if (present(index_end)) index_end = n_cells + n - n_cells = n_cells + n + + ! Extend the C++ cells array and get pointers to the C++ objects + call extend_cells_c(n) + do i = n_cells - n, n_cells + cells(i) % ptr = cell_pointer_c(i - 1) + end do err = 0 end function openmc_extend_cells @@ -490,14 +609,14 @@ 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) indices = C_LOC(c % material(1)) case (FILL_UNIVERSE, FILL_LATTICE) n = 1 - indices = C_LOC(c % fill) + indices = cell_fill_ptr(c % ptr) end select end associate else @@ -514,7 +633,7 @@ contains integer(C_INT) :: err if (index >= 1 .and. index <= size(cells)) then - id = cells(index) % id + id = cells(index) % id() err = 0 else err = E_OUT_OF_BOUNDS @@ -541,7 +660,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 @@ -553,9 +672,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 @@ -573,7 +692,7 @@ contains integer(C_INT) :: err if (index >= 1 .and. index <= n_cells) then - cells(index) % id = id + call cells(index) % set_id(id) call cell_dict % set(id, index) err = 0 else @@ -609,7 +728,7 @@ contains if (index >= 1 .and. index <= size(cells)) then ! error if the cell is filled with another universe - if (cells(index) % fill /= NONE) then + if (cells(index) % fill() /= C_NONE) then err = E_GEOMETRY call set_errmsg("Cannot set temperature on a cell filled & &with a universe.") diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h index e38a31e997..353d0137af 100644 --- a/src/hdf5_interface.h +++ b/src/hdf5_interface.h @@ -85,12 +85,21 @@ extern "C" void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_ const double* results); template void -write_double_1D(hid_t group_id, char const *name, - std::array &buffer) +write_int(hid_t group_id, char const *name, + const std::array &buffer, bool indep) +{ + hsize_t dims[1] {array_len}; + write_dataset(group_id, 1, dims, name, H5T_NATIVE_INT, buffer.data(), indep); +} + + +template void +write_double(hid_t group_id, char const *name, + const std::array &buffer, bool indep) { hsize_t dims[1] {array_len}; write_dataset(group_id, 1, dims, name, H5T_NATIVE_DOUBLE, - buffer.data(), false); + buffer.data(), indep); } } // namespace openmc diff --git a/src/input_xml.F90 b/src/input_xml.F90 index f8cf50dd3c..e24f3fe253 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: neighbor_lists use geometry_header use hdf5_interface use list_header, only: ListChar, ListInt, ListReal @@ -48,11 +47,50 @@ module input_xml save interface - subroutine read_surfaces(node_ptr) bind(C, name='read_surfaces') - use ISO_C_BINDING - implicit none + subroutine adjust_indices_c() bind(C) + end subroutine adjust_indices_c + + subroutine allocate_offset_tables(n_maps) bind(C) + import C_INT + integer(C_INT), intent(in), value :: n_maps + end subroutine allocate_offset_tables + + subroutine fill_offset_tables(target_univ_id, map) bind(C) + import C_INT32_T, C_INT + integer(C_INT32_T), intent(in), value :: target_univ_id + integer(C_INT), intent(in), value :: map + end subroutine fill_offset_tables + + subroutine count_cell_instances(univ_indx) bind(C) + import C_INT32_T + integer(C_INT32_T), intent(in), value :: univ_indx + end subroutine count_cell_instances + + subroutine read_surfaces(node_ptr) bind(C) + import C_PTR type(C_PTR) :: node_ptr end subroutine read_surfaces + + subroutine read_cells(node_ptr) bind(C) + import C_PTR + type(C_PTR) :: node_ptr + end subroutine read_cells + + subroutine read_lattices(node_ptr) bind(C) + import C_PTR + type(C_PTR) :: node_ptr + end subroutine read_lattices + + function find_root_universe() bind(C) result(root) + import C_INT32_T + integer(C_INT32_T) :: root + end function find_root_universe + + function maximum_levels(univ) bind(C) result(n) + import C_INT32_T, C_INT + integer(C_INT32_T), intent(in), value :: univ + integer(C_INT) :: n + end function maximum_levels end interface contains @@ -122,7 +160,7 @@ contains ! Perform some final operations to set up the geometry call adjust_indices() - call count_instance(universes(root_universe)) + call count_cell_instances(root_universe-1) ! After reading input and basic geometry setup is complete, build lists of ! neighboring cells for efficient tracking @@ -137,7 +175,7 @@ contains ! Check to make sure there are not too many nested coordinate levels in the ! geometry since the coordinate list is statically allocated for performance ! reasons - if (maximum_levels(universes(root_universe)) > MAX_COORD) then + if (maximum_levels(root_universe - 1) > MAX_COORD) then call fatal_error("Too many nested coordinate levels in the geometry. & &Try increasing the maximum number of coordinate levels by & &providing the CMake -Dmaxcoord= option.") @@ -913,12 +951,11 @@ contains subroutine read_geometry_xml() - integer :: i, j, k, m, i_x, i_a, input_index - integer :: n, n_mats, n_x, n_y, n_z, n_rings, n_rlats, n_hlats + integer :: i, j, k + integer :: n, n_mats, n_rlats, n_hlats integer :: id integer :: univ_id integer :: n_cells_in_univ - integer, allocatable :: temp_int_array(:) real(8) :: phi, theta, psi logical :: file_exists logical :: boundary_exists @@ -935,8 +972,6 @@ contains type(XMLNode), allocatable :: node_rlat_list(:) type(XMLNode), allocatable :: node_hlat_list(:) type(VectorInt) :: tokens - type(VectorInt) :: rpn - type(VectorInt) :: fill_univ_ids ! List of fill universe IDs type(VectorInt) :: univ_ids ! List of all universe IDs type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each ! universe contains @@ -988,6 +1023,8 @@ contains ! ========================================================================== ! READ CELLS FROM GEOMETRY.XML + call read_cells(root % ptr) + ! Get pointer to list of XML call get_node_list(root, "cell", node_cell_list) @@ -1011,42 +1048,18 @@ contains do i = 1, n_cells c => cells(i) + 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 node_cell = node_cell_list(i) - ! Copy data into cells - if (check_for_node(node_cell, "id")) then - call get_node_value(node_cell, "id", c % id) - else - call fatal_error("Must specify id of cell in geometry XML file.") - end if - - ! Copy cell name - if (check_for_node(node_cell, "name")) then - call get_node_value(node_cell, "name", c % name) - end if - - if (check_for_node(node_cell, "universe")) then - call get_node_value(node_cell, "universe", c % universe) - else - c % universe = 0 - end if - if (check_for_node(node_cell, "fill")) then - call get_node_value(node_cell, "fill", c % fill) - if (find(fill_univ_ids, c % fill) == -1) & - call fill_univ_ids % push_back(c % fill) - else - c % fill = NONE - end if - ! Check to make sure 'id' hasn't been used - if (cell_dict % has(c % id)) then + if (cell_dict % has(c % id())) then call fatal_error("Two or more cells use the same unique ID: " & - // to_str(c % id)) + // to_str(c % id())) end if ! Read material @@ -1068,7 +1081,7 @@ contains ! Check for error if (c % material(j) == ERROR_INT) then call fatal_error("Invalid material specified on cell " & - // to_str(c % id)) + // to_str(c % id())) end if end select end do @@ -1085,18 +1098,6 @@ contains c % material(1) = NONE end if - ! Check to make sure that either material or fill was specified - if (c % material(1) == NONE .and. c % fill == NONE) then - call fatal_error("Neither material nor fill was specified for cell " & - // trim(to_str(c % id))) - end if - - ! Check to make sure that both material and fill haven't been - ! specified simultaneously - if (c % material(1) /= NONE .and. c % fill /= NONE) then - call fatal_error("Cannot specify material and fill simultaneously") - 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 & @@ -1126,42 +1127,28 @@ contains end if end do - ! Use shunting-yard algorithm to determine RPN for surface algorithm - call generate_rpn(c%id, tokens, rpn) - ! Copy region spec and RPN form to cell arrays allocate(c % region(tokens%size())) - allocate(c % rpn(rpn%size())) c % region(:) = tokens%data(1:tokens%size()) - c % rpn(:) = rpn%data(1:rpn%size()) call tokens%clear() - call rpn%clear() end if if (.not. allocated(c%region)) allocate(c%region(0)) - if (.not. allocated(c%rpn)) allocate(c%rpn(0)) - - ! Check if this is a simple cell - if (any(c%rpn == OP_COMPLEMENT) .or. any(c%rpn == OP_UNION)) then - c%simple = .false. - else - c%simple = .true. - end if ! Rotation matrix if (check_for_node(node_cell, "rotation")) then ! Rotations can only be applied to cells that are being filled with ! another universe - if (c % fill == NONE) then + if (c % fill() == C_NONE) then call fatal_error("Cannot apply a rotation to cell " // trim(to_str(& - &c % id)) // " because it is not filled with another universe") + &c % id())) // " because it is not filled with another universe") end if ! Read number of rotation parameters n = node_word_count(node_cell, "rotation") if (n /= 3) then call fatal_error("Incorrect number of rotation parameters on cell " & - // to_str(c % id)) + // to_str(c % id())) end if ! Copy rotation angles in x,y,z directions @@ -1187,9 +1174,9 @@ contains if (check_for_node(node_cell, "translation")) then ! Translations can only be applied to cells that are being filled with ! another universe - if (c % fill == NONE) then + if (c % fill() == C_NONE) then call fatal_error("Cannot apply a translation to cell " & - // trim(to_str(c % id)) // " because it is not filled with & + // trim(to_str(c % id())) // " because it is not filled with & &another universe") end if @@ -1197,7 +1184,7 @@ contains n = node_word_count(node_cell, "translation") if (n /= 3) then call fatal_error("Incorrect number of translation parameters on & - &cell " // to_str(c % id)) + &cell " // to_str(c % id())) end if ! Copy translation vector @@ -1213,7 +1200,7 @@ contains if (n > 0) then ! Make sure this is a "normal" cell. if (c % material(1) == NONE) call fatal_error("Cell " & - // trim(to_str(c % id)) // " was specified with a temperature & + // trim(to_str(c % id())) // " was specified with a temperature & &but no material. Temperature specification is only valid for & &cells filled with a material.") @@ -1224,7 +1211,7 @@ contains ! Make sure all temperatues are positive do j = 1, size(c % sqrtkT) if (c % sqrtkT(j) < ZERO) call fatal_error("Cell " & - // trim(to_str(c % id)) // " was specified with a negative & + // trim(to_str(c % id())) // " was specified with a negative & &temperature. All cell temperatures must be non-negative.") end do @@ -1240,12 +1227,12 @@ contains end if ! Add cell to dictionary - call cell_dict % set(c % id, i) + call cell_dict % set(c % id(), i) ! For cells, we also need to check if there's a new universe -- ! also for every cell add 1 to the count of cells for the ! specified universe - univ_id = c % universe + univ_id = c % universe() if (.not. cells_in_univ_dict % has(univ_id)) then n_universes = n_universes + 1 n_cells_in_univ = 1 @@ -1261,6 +1248,8 @@ contains ! ========================================================================== ! READ LATTICES FROM GEOMETRY.XML + call read_lattices(root % ptr) + ! Get pointer to list of XML call get_node_list(root, "lattice", node_rlat_list) call get_node_list(root, "hex_lattice", node_hlat_list) @@ -1268,137 +1257,20 @@ contains ! Allocate lattices array n_rlats = size(node_rlat_list) n_hlats = size(node_hlat_list) - n_lattices = n_rlats + n_hlats - allocate(lattices(n_lattices)) + allocate(lattices(n_rlats + n_hlats)) RECT_LATTICES: do i = 1, n_rlats allocate(RectLattice::lattices(i) % obj) lat => lattices(i) % obj + lat % ptr = lattice_pointer_c(i - 1) select type(lat) type is (RectLattice) ! Get pointer to i-th lattice node_lat = node_rlat_list(i) - ! ID of lattice - if (check_for_node(node_lat, "id")) then - call get_node_value(node_lat, "id", lat % id) - else - call fatal_error("Must specify id of lattice in geometry XML file.") - end if - - ! Check to make sure 'id' hasn't been used - if (lattice_dict % has(lat % id)) then - call fatal_error("Two or more lattices use the same unique ID: " & - // to_str(lat % id)) - end if - - ! Copy lattice name - if (check_for_node(node_lat, "name")) then - call get_node_value(node_lat, "name", lat % name) - end if - - ! Read number of lattice cells in each dimension - n = node_word_count(node_lat, "dimension") - if (n == 2) then - call get_node_array(node_lat, "dimension", lat % n_cells(1:2)) - lat % n_cells(3) = 1 - lat % is_3d = .false. - else if (n == 3) then - call get_node_array(node_lat, "dimension", lat % n_cells) - lat % is_3d = .true. - else - call fatal_error("Rectangular lattice must be two or three dimensions.") - end if - - ! Read lattice lower-left location - if (node_word_count(node_lat, "lower_left") /= n) then - call fatal_error("Number of entries on must be the same & - &as the number of entries on .") - end if - - allocate(lat % lower_left(n)) - call get_node_array(node_lat, "lower_left", lat % lower_left) - - ! Read lattice pitches. - ! TODO: Remove this deprecation warning in a future release. - if (check_for_node(node_lat, "width")) then - call warning("The use of 'width' is deprecated and will be disallowed & - &in a future release. Use 'pitch' instead. The utility openmc/& - &src/utils/update_inputs.py can be used to automatically update & - &geometry.xml files.") - if (node_word_count(node_lat, "width") /= n) then - call fatal_error("Number of entries on must be the same as & - &the number of entries on .") - end if - - else if (node_word_count(node_lat, "pitch") /= n) then - call fatal_error("Number of entries on must be the same as & - &the number of entries on .") - end if - - allocate(lat % pitch(n)) - ! TODO: Remove the 'width' code in a future release. - if (check_for_node(node_lat, "width")) then - call get_node_array(node_lat, "width", lat % pitch) - else - call get_node_array(node_lat, "pitch", lat % pitch) - end if - - ! TODO: Remove deprecation warning in a future release. - if (check_for_node(node_lat, "type")) then - call warning("The use of 'type' is no longer needed. The utility & - &openmc/src/utils/update_inputs.py can be used to automatically & - &update geometry.xml files.") - end if - - ! Copy number of dimensions - n_x = lat % n_cells(1) - n_y = lat % n_cells(2) - n_z = lat % n_cells(3) - allocate(lat % universes(n_x, n_y, n_z)) - - ! Check that number of universes matches size - n = node_word_count(node_lat, "universes") - if (n /= n_x*n_y*n_z) then - call fatal_error("Number of universes on does not match & - &size of lattice " // trim(to_str(lat % id)) // ".") - end if - - allocate(temp_int_array(n)) - call get_node_array(node_lat, "universes", temp_int_array) - - ! Read universes - do m = 1, n_z - do k = 0, n_y - 1 - do j = 1, n_x - lat % universes(j, n_y - k, m) = & - temp_int_array(j + n_x*k + n_x*n_y*(m-1)) - if (find(fill_univ_ids, lat % universes(j, n_y - k, m)) == -1) & - call fill_univ_ids % push_back(lat % universes(j, n_y - k, m)) - end do - end do - end do - deallocate(temp_int_array) - - ! Read outer universe for area outside lattice. - lat % outer = NO_OUTER_UNIVERSE - if (check_for_node(node_lat, "outer")) then - call get_node_value(node_lat, "outer", lat % outer) - if (find(fill_univ_ids, lat % outer) == -1) & - call fill_univ_ids % push_back(lat % outer) - end if - - ! Check for 'outside' nodes which are no longer supported. - if (check_for_node(node_lat, "outside")) then - call fatal_error("The use of 'outside' in lattices is no longer & - &supported. Instead, use 'outer' which defines a universe rather & - &than a material. The utility openmc/src/utils/update_inputs.py & - &can be used automatically replace 'outside' with 'outer'.") - end if - ! Add lattice to dictionary - call lattice_dict % set(lat % id, i) + call lattice_dict % set(lat % id(), i) end select end do RECT_LATTICES @@ -1406,186 +1278,15 @@ contains HEX_LATTICES: do i = 1, n_hlats allocate(HexLattice::lattices(n_rlats + i) % obj) lat => lattices(n_rlats + i) % obj + lat % ptr = lattice_pointer_c(n_rlats + i - 1) select type (lat) type is (HexLattice) ! Get pointer to i-th lattice node_lat = node_hlat_list(i) - ! ID of lattice - if (check_for_node(node_lat, "id")) then - call get_node_value(node_lat, "id", lat % id) - else - call fatal_error("Must specify id of lattice in geometry XML file.") - end if - - ! Check to make sure 'id' hasn't been used - if (lattice_dict % has(lat % id)) then - call fatal_error("Two or more lattices use the same unique ID: " & - // to_str(lat % id)) - end if - - ! Copy lattice name - if (check_for_node(node_lat, "name")) then - call get_node_value(node_lat, "name", lat % name) - end if - - ! Read number of lattice cells in each dimension - call get_node_value(node_lat, "n_rings", lat % n_rings) - if (check_for_node(node_lat, "n_axial")) then - call get_node_value(node_lat, "n_axial", lat % n_axial) - lat % is_3d = .true. - else - lat % n_axial = 1 - lat % is_3d = .false. - end if - - ! Read lattice lower-left location - n = node_word_count(node_lat, "center") - if (lat % is_3d .and. n /= 3) then - call fatal_error("A hexagonal lattice with must have & - &
specified by 3 numbers.") - else if ((.not. lat % is_3d) .and. n /= 2) then - call fatal_error("A hexagonal lattice without must have & - &
specified by 2 numbers.") - end if - - allocate(lat % center(n)) - call get_node_array(node_lat, "center", lat % center) - - ! Read lattice pitches - n = node_word_count(node_lat, "pitch") - if (lat % is_3d .and. n /= 2) then - call fatal_error("A hexagonal lattice with must have & - &specified by 2 numbers.") - else if ((.not. lat % is_3d) .and. n /= 1) then - call fatal_error("A hexagonal lattice without must have & - & specified by 1 number.") - end if - - allocate(lat % pitch(n)) - call get_node_array(node_lat, "pitch", lat % pitch) - - ! Copy number of dimensions - n_rings = lat % n_rings - n_z = lat % n_axial - allocate(lat % universes(2*n_rings - 1, 2*n_rings - 1, n_z)) - - ! Check that number of universes matches size - n = node_word_count(node_lat, "universes") - if (n /= (3*n_rings**2 - 3*n_rings + 1)*n_z) then - call fatal_error("Number of universes on does not match & - &size of lattice " // trim(to_str(lat % id)) // ".") - end if - - allocate(temp_int_array(n)) - call get_node_array(node_lat, "universes", temp_int_array) - - ! Read universes - ! Universes in hexagonal lattices are stored in a manner that represents - ! a skewed coordinate system: (x, alpha) rather than (x, y). There is - ! no obvious, direct relationship between the order of universes in the - ! input and the order that they will be stored in the skewed array so - ! the following code walks a set of index values across the skewed array - ! in a manner that matches the input order. Note that i_x = 0, i_a = 0 - ! corresponds to the center of the hexagonal lattice. - - input_index = 1 - do m = 1, n_z - ! Initialize lattice indecies. - i_x = 1 - i_a = n_rings - 1 - - ! Map upper triangular region of hexagonal lattice. - do k = 1, n_rings-1 - ! Walk index to lower-left neighbor of last row start. - i_x = i_x - 1 - do j = 1, k - ! Place universe in array. - lat % universes(i_x + n_rings, i_a + n_rings, m) = & - temp_int_array(input_index) - if (find(fill_univ_ids, temp_int_array(input_index)) == -1) & - call fill_univ_ids % push_back(temp_int_array(input_index)) - ! Walk index to closest non-adjacent right neighbor. - i_x = i_x + 2 - i_a = i_a - 1 - ! Increment XML array index. - input_index = input_index + 1 - end do - ! Return lattice index to start of current row. - i_x = i_x - 2*k - i_a = i_a + k - end do - - ! Map middle square region of hexagonal lattice. - do k = 1, 2*n_rings - 1 - if (mod(k, 2) == 1) then - ! Walk index to lower-left neighbor of last row start. - i_x = i_x - 1 - else - ! Walk index to lower-right neighbor of last row start - i_x = i_x + 1 - i_a = i_a - 1 - end if - do j = 1, n_rings - mod(k-1, 2) - ! Place universe in array. - lat % universes(i_x + n_rings, i_a + n_rings, m) = & - temp_int_array(input_index) - if (find(fill_univ_ids, temp_int_array(input_index)) == -1) & - call fill_univ_ids % push_back(temp_int_array(input_index)) - ! Walk index to closest non-adjacent right neighbor. - i_x = i_x + 2 - i_a = i_a - 1 - ! Increment XML array index. - input_index = input_index + 1 - end do - ! Return lattice index to start of current row. - i_x = i_x - 2*(n_rings - mod(k-1, 2)) - i_a = i_a + n_rings - mod(k-1, 2) - end do - - ! Map lower triangular region of hexagonal lattice. - do k = 1, n_rings-1 - ! Walk index to lower-right neighbor of last row start. - i_x = i_x + 1 - i_a = i_a - 1 - do j = 1, n_rings - k - ! Place universe in array. - lat % universes(i_x + n_rings, i_a + n_rings, m) = & - temp_int_array(input_index) - if (find(fill_univ_ids, temp_int_array(input_index)) == -1) & - call fill_univ_ids % push_back(temp_int_array(input_index)) - ! Walk index to closest non-adjacent right neighbor. - i_x = i_x + 2 - i_a = i_a - 1 - ! Increment XML array index. - input_index = input_index + 1 - end do - ! Return lattice index to start of current row. - i_x = i_x - 2*(n_rings - k) - i_a = i_a + n_rings - k - end do - end do - deallocate(temp_int_array) - - ! Read outer universe for area outside lattice. - lat % outer = NO_OUTER_UNIVERSE - if (check_for_node(node_lat, "outer")) then - call get_node_value(node_lat, "outer", lat % outer) - if (find(fill_univ_ids, lat % outer) == -1) & - call fill_univ_ids % push_back(lat % outer) - end if - - ! Check for 'outside' nodes which are no longer supported. - if (check_for_node(node_lat, "outside")) then - call fatal_error("The use of 'outside' in lattices is no longer & - &supported. Instead, use 'outer' which defines a universe rather & - &than a material. The utility openmc/src/utils/update_inputs.py & - &can be used automatically replace 'outside' with 'outer'.") - end if - ! Add lattice to dictionary - call lattice_dict % set(lat % id, n_rlats + i) + call lattice_dict % set(lat % id(), n_rlats + i) end select end do HEX_LATTICES @@ -1603,23 +1304,13 @@ contains n_cells_in_univ = cells_in_univ_dict % get(u % id) allocate(u % cells(n_cells_in_univ)) u % cells(:) = 0 - - ! Check whether universe is a fill universe - if (find(fill_univ_ids, u % id) == -1) then - if (root_universe > 0) then - call 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_universe = i - end if - end if end associate end do + root_universe = find_root_universe() + 1 do i = 1, n_cells ! Get index in universes array - j = universe_dict % get(cells(i) % universe) + j = universe_dict % get(cells(i) % universe()) ! Set the first zero entry in the universe cells array to the index in the ! global cells array @@ -3743,93 +3434,6 @@ contains end subroutine read_mg_cross_sections_header -!=============================================================================== -! GENERATE_RPN implements the shunting-yard algorithm to generate a Reverse -! Polish notation (RPN) expression for the region specification of a cell given -! the infix notation. -!=============================================================================== - - subroutine generate_rpn(cell_id, tokens, output) - integer, intent(in) :: cell_id - type(VectorInt), intent(in) :: tokens ! infix notation - type(VectorInt), intent(inout) :: output ! RPN notation - - integer :: i - integer :: token - integer :: op - type(VectorInt) :: stack - - do i = 1, tokens%size() - token = tokens%data(i) - - if (token < OP_UNION) then - ! If token is not an operator, add it to output - call output%push_back(token) - - elseif (token < OP_RIGHT_PAREN) then - ! Regular operators union, intersection, complement - do while (stack%size() > 0) - op = stack%data(stack%size()) - - if (op < OP_RIGHT_PAREN .and. & - ((token == OP_COMPLEMENT .and. token < op) .or. & - (token /= OP_COMPLEMENT .and. token <= op))) then - ! While there is an operator, op, on top of the stack, if the token - ! is left-associative and its precedence is less than or equal to - ! that of op or if the token is right-associative and its precedence - ! is less than that of op, move op to the output queue and push the - ! token on to the stack. Note that only complement is - ! right-associative. - call output%push_back(op) - call stack%pop_back() - else - exit - end if - end do - - call stack%push_back(token) - - elseif (token == OP_LEFT_PAREN) then - ! If the token is a left parenthesis, push it onto the stack - call stack%push_back(token) - - else - ! If the token is a right parenthesis, move operators from the stack to - ! the output queue until reaching the left parenthesis. - do - ! If we run out of operators without finding a left parenthesis, it - ! means there are mismatched parentheses. - if (stack%size() == 0) then - call fatal_error('Mimatched parentheses in region specification & - &for cell ' // trim(to_str(cell_id)) // '.') - end if - - op = stack%data(stack%size()) - if (op == OP_LEFT_PAREN) exit - call output%push_back(op) - call stack%pop_back() - end do - - ! Pop the left parenthesis. - call stack%pop_back() - end if - end do - - ! While there are operators on the stack, move them to the output queue - do while (stack%size() > 0) - op = stack%data(stack%size()) - - ! If the operator is a parenthesis, it is mismatched - if (op >= OP_RIGHT_PAREN) then - call fatal_error('Mimatched parentheses in region specification & - &for cell ' // trim(to_str(cell_id)) // '.') - end if - - call output%push_back(op) - call stack%pop_back() - end do - end subroutine generate_rpn - !=============================================================================== ! NORMALIZE_AO Normalize the nuclide atom percents !=============================================================================== @@ -4150,117 +3754,31 @@ contains integer :: i ! index for various purposes integer :: j ! index for various purposes - integer :: k ! loop index for lattices - integer :: m ! loop index for lattices - integer :: lid ! lattice IDs integer :: id ! user-specified id - class(Lattice), pointer :: lat => null() + + call adjust_indices_c() do i = 1, n_cells - ! ======================================================================= - ! ADJUST UNIVERSE INDEX FOR EACH CELL associate (c => cells(i)) - id = c % universe - if (universe_dict % has(id)) then - c % universe = universe_dict % get(id) - else - call fatal_error("Could not find universe " // trim(to_str(id)) & - &// " specified on cell " // trim(to_str(c % id))) - end if - ! ======================================================================= ! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL - 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 + if (c % material(1) /= NONE) then 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)) & - // " specified on cell " // trim(to_str(c % id))) + // " specified on cell " // trim(to_str(c % id()))) end if end do end if end associate end do - ! ========================================================================== - ! ADJUST UNIVERSE INDICES FOR EACH LATTICE - - do i = 1, n_lattices - lat => lattices(i) % obj - select type (lat) - - type is (RectLattice) - do m = 1, lat % n_cells(3) - do k = 1, lat % n_cells(2) - do j = 1, lat % n_cells(1) - id = lat % universes(j,k,m) - if (universe_dict % has(id)) then - lat % universes(j,k,m) = universe_dict % get(id) - else - call fatal_error("Invalid universe number " & - &// trim(to_str(id)) // " specified on lattice " & - &// trim(to_str(lat % id))) - end if - end do - end do - end do - - type is (HexLattice) - do m = 1, lat % n_axial - do k = 1, 2*lat % n_rings - 1 - do j = 1, 2*lat % n_rings - 1 - if (j + k < lat % n_rings + 1) then - cycle - else if (j + k > 3*lat % n_rings - 1) then - cycle - end if - id = lat % universes(j, k, m) - if (universe_dict % has(id)) then - lat % universes(j, k, m) = universe_dict % get(id) - else - call fatal_error("Invalid universe number " & - &// trim(to_str(id)) // " specified on lattice " & - &// trim(to_str(lat % id))) - end if - end do - end do - end do - - end select - - if (lat % outer /= NO_OUTER_UNIVERSE) then - if (universe_dict % has(lat % outer)) then - lat % outer = universe_dict % get(lat % outer) - else - call fatal_error("Invalid universe number " & - &// trim(to_str(lat % outer)) & - &// " specified on lattice " // trim(to_str(lat % id))) - end if - end if - - end do - end subroutine adjust_indices !=============================================================================== @@ -4270,94 +3788,11 @@ contains subroutine prepare_distribcell() - integer :: i, j ! Tally, filter loop counters - logical :: distribcell_active ! Does simulation use distribcell? - integer, allocatable :: univ_list(:) ! Target offsets - integer, allocatable :: counts(:,:) ! Target count - logical, allocatable :: found(:,:) ! Target found + integer :: i, j, k + type(SetInt) :: cell_list ! distribcells to track + type(ListInt) :: univ_list ! universes containing distribcells - ! Assume distribcell is not needed until proven otherwise. - distribcell_active = .false. - - ! We need distribcell if any tallies have distribcell filters. - do i = 1, n_tallies - do j = 1, size(tallies(i) % obj % filter) - select type(filt => filters(tallies(i) % obj % filter(j)) % obj) - type is (DistribcellFilter) - distribcell_active = .true. - end select - end do - end do - - ! We also need distribcell if any distributed materials or distributed - ! temperatues are present. - if (.not. distribcell_active) then - do i = 1, n_cells - if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then - distribcell_active = .true. - exit - end if - end do - end if - - ! If distribcell isn't used in this simulation then no more work left to do. - if (.not. distribcell_active) return - - ! Make sure the number of materials and temperatures matches the number of - ! cell instances. - do i = 1, n_cells - associate (c => cells(i)) - if (size(c % material) > 1) then - if (size(c % material) /= c % 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)) & - // " 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 - 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)) & - // " distributed instances. The number of temperatures must & - &equal one or the number of instances.") - end if - end if - end associate - end do - - ! Allocate offset maps at each level in the geometry - call allocate_offsets(univ_list, counts, found) - - ! Calculate offsets for each target distribcell - do i = 1, n_maps - do j = 1, n_universes - call calc_offsets(univ_list(i), i, universes(j), counts, found) - end do - end do - - end subroutine prepare_distribcell - -!=============================================================================== -! ALLOCATE_OFFSETS determines the number of maps needed and allocates required -! memory for distribcell offset tables -!=============================================================================== - - recursive subroutine allocate_offsets(univ_list, counts, found) - - integer, intent(out), allocatable :: univ_list(:) ! Target offsets - integer, intent(out), allocatable :: counts(:,:) ! Target count - logical, intent(out), allocatable :: found(:,:) ! Target found - - integer :: i, j, k ! Loop counters - type(SetInt) :: cell_list ! distribells to track - - ! Begin gathering list of cells in distribcell tallies - n_maps = 0 - - ! List all cells referenced in distribcell filters. + ! Find all cells listed in a distribcell filter. do i = 1, n_tallies do j = 1, size(tallies(i) % obj % filter) select type(filt => filters(tallies(i) % obj % filter(j)) % obj) @@ -4367,35 +3802,38 @@ contains end do end do - ! List all cells with multiple (distributed) materials or temperatures. + ! Find all cells with multiple (distributed) materials or temperatures. do i = 1, n_cells if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then call cell_list % add(i) end if end do - ! Compute the number of unique universes containing these distribcells - ! to determine the number of offset tables to allocate - do i = 1, n_universes - do j = 1, size(universes(i) % cells) - if (cell_list % contains(universes(i) % cells(j))) then - n_maps = n_maps + 1 + ! Make sure the number of distributed materials and temperatures matches the + ! number of respective cell instances. + do i = 1, n_cells + associate (c => cells(i)) + if (size(c % material) > 1) 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 % n_instances())) & + // " distributed instances. The number of materials must & + &equal one or the number of instances.") + end if end if - end do + if (size(c % sqrtkT) > 1) 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 % n_instances()))& + // " distributed instances. The number of temperatures must & + &equal one or the number of instances.") + end if + end if + end associate end do - ! Allocate the list of offset tables for each unique universe - allocate(univ_list(n_maps)) - - ! Allocate list to accumulate target distribcell counts in each universe - allocate(counts(n_universes, n_maps)) - counts(:,:) = 0 - - ! Allocate list to track if target distribcells are found in each universe - allocate(found(n_universes, n_maps)) - found(:,:) = .false. - - ! Search through universes for distributed cells and assign each one a ! unique distribcell array index. k = 1 @@ -4403,39 +3841,18 @@ contains do j = 1, size(universes(i) % cells) if (cell_list % contains(universes(i) % cells(j))) then cells(universes(i) % cells(j)) % distribcell_index = k - univ_list(k) = universes(i) % id + call univ_list % append(universes(i) % id) k = k + 1 end if end do end do - ! Allocate the offset tables for lattices - do i = 1, n_lattices - associate(lat => lattices(i) % obj) - select type(lat) - - type is (RectLattice) - allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), & - lat % n_cells(3))) - type is (HexLattice) - allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, & - 2 * lat % n_rings - 1, lat % n_axial)) - end select - - lat % offset(:, :, :, :) = 0 - end associate + ! Allocate and fill cell and lattice offset tables. + call allocate_offset_tables(univ_list % size()) + do i = 1, univ_list % size() + call fill_offset_tables(univ_list % get_item(i), i-1) end do - ! Allocate offset table for fill cells - do i = 1, n_cells - if (cells(i) % type /= FILL_MATERIAL) then - allocate(cells(i) % offset(n_maps)) - end if - end do - - ! Free up memory - call cell_list % clear() - - end subroutine allocate_offsets + end subroutine prepare_distribcell end module input_xml diff --git a/src/lattice.cpp b/src/lattice.cpp new file mode 100644 index 0000000000..959ab932e4 --- /dev/null +++ b/src/lattice.cpp @@ -0,0 +1,946 @@ +#include "lattice.h" + +#include +#include +#include + +#include "cell.h" +#include "error.h" +#include "geometry_aux.h" +#include "hdf5_interface.h" +#include "string_utils.h" +#include "xml_interface.h" + + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +std::vector lattices_c; + +std::unordered_map lattice_map; + +//============================================================================== +// Lattice implementation +//============================================================================== + +Lattice::Lattice(pugi::xml_node lat_node) +{ + if (check_for_node(lat_node, "id")) { + id = stoi(get_node_value(lat_node, "id")); + } else { + fatal_error("Must specify id of lattice in geometry XML file."); + } + + if (check_for_node(lat_node, "name")) { + name = get_node_value(lat_node, "name"); + } + + if (check_for_node(lat_node, "outer")) { + outer = stoi(get_node_value(lat_node, "outer")); + } +} + +//============================================================================== + +LatticeIter Lattice::begin() +{return LatticeIter(*this, 0);} + +LatticeIter Lattice::end() +{return LatticeIter(*this, universes.size());} + +ReverseLatticeIter Lattice::rbegin() +{return ReverseLatticeIter(*this, universes.size()-1);} + +ReverseLatticeIter Lattice::rend() +{return ReverseLatticeIter(*this, -1);} + +//============================================================================== + +void +Lattice::adjust_indices() +{ + // Adjust the indices for the universes array. + for (LatticeIter it = begin(); it != end(); ++it) { + int uid = *it; + auto search = universe_map.find(uid); + if (search != universe_map.end()) { + *it = search->second; + } else { + std::stringstream err_msg; + err_msg << "Invalid universe number " << uid << " specified on " + "lattice " << id; + fatal_error(err_msg); + } + } + + // Adjust the index for the outer universe. + if (outer != NO_OUTER_UNIVERSE) { + auto search = universe_map.find(outer); + if (search != universe_map.end()) { + outer = search->second; + } else { + std::stringstream err_msg; + err_msg << "Invalid universe number " << outer << " specified on " + "lattice " << id; + fatal_error(err_msg); + } + } +} + +//============================================================================== + +int32_t +Lattice::fill_offset_table(int32_t offset, int32_t target_univ_id, int map) +{ + for (LatticeIter it = begin(); it != end(); ++it) { + offsets[map * universes.size() + it.indx] = offset; + offset += count_universe_instances(*it, target_univ_id); + } + return offset; +} + +//============================================================================== + +void +Lattice::to_hdf5(hid_t lattices_group) const +{ + // Make a group for the lattice. + std::string group_name {"lattice "}; + group_name += std::to_string(id); + hid_t lat_group = create_group(lattices_group, group_name); + + // Write the name and outer universe. + if (!name.empty()) { + write_string(lat_group, "name", name, false); + } + + if (outer != NO_OUTER_UNIVERSE) { + int32_t outer_id = global_universes[outer]->id; + write_int(lat_group, 0, nullptr, "outer", &outer_id, false); + } else { + write_int(lat_group, 0, nullptr, "outer", &outer, false); + } + + // Call subclass-overriden function to fill in other details. + to_hdf5_inner(lat_group); + + close_group(lat_group); +} + +//============================================================================== +// RectLattice implementation +//============================================================================== + +RectLattice::RectLattice(pugi::xml_node lat_node) + : Lattice {lat_node} +{ + // Read the number of lattice cells in each dimension. + std::string dimension_str {get_node_value(lat_node, "dimension")}; + std::vector dimension_words {split(dimension_str)}; + if (dimension_words.size() == 2) { + n_cells[0] = stoi(dimension_words[0]); + n_cells[1] = stoi(dimension_words[1]); + n_cells[2] = 1; + is_3d = false; + } else if (dimension_words.size() == 3) { + n_cells[0] = stoi(dimension_words[0]); + n_cells[1] = stoi(dimension_words[1]); + n_cells[2] = stoi(dimension_words[2]); + is_3d = true; + } else { + fatal_error("Rectangular lattice must be two or three dimensions."); + } + + // Read the lattice lower-left location. + std::string ll_str {get_node_value(lat_node, "lower_left")}; + std::vector ll_words {split(ll_str)}; + if (ll_words.size() != dimension_words.size()) { + fatal_error("Number of entries on must be the same as the " + "number of entries on ."); + } + lower_left[0] = stod(ll_words[0]); + lower_left[1] = stod(ll_words[1]); + if (is_3d) {lower_left[2] = stod(ll_words[2]);} + + // Read the lattice pitches. + std::string pitch_str {get_node_value(lat_node, "pitch")}; + std::vector pitch_words {split(pitch_str)}; + if (pitch_words.size() != dimension_words.size()) { + fatal_error("Number of entries on must be the same as the " + "number of entries on ."); + } + pitch[0] = stod(pitch_words[0]); + pitch[1] = stod(pitch_words[1]); + if (is_3d) {pitch[2] = stod(pitch_words[2]);} + + // Read the universes and make sure the correct number was specified. + std::string univ_str {get_node_value(lat_node, "universes")}; + std::vector univ_words {split(univ_str)}; + if (univ_words.size() != nx*ny*nz) { + std::stringstream err_msg; + err_msg << "Expected " << nx*ny*nz + << " universes for a rectangular lattice of size " + << nx << "x" << ny << "x" << nz << " but " << univ_words.size() + << " were specified."; + fatal_error(err_msg); + } + + // Parse the universes. + universes.resize(nx*ny*nz, C_NONE); + for (int iz = 0; iz < nz; iz++) { + for (int iy = ny-1; iy > -1; iy--) { + for (int ix = 0; ix < nx; ix++) { + int indx1 = nx*ny*iz + nx*(ny-iy-1) + ix; + int indx2 = nx*ny*iz + nx*iy + ix; + universes[indx1] = stoi(univ_words[indx2]); + } + } + } +} + +//============================================================================== + +int32_t& +RectLattice::operator[](const int i_xyz[3]) +{ + int indx = nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]; + return universes[indx]; +} + +//============================================================================== + +bool +RectLattice::are_valid_indices(const int i_xyz[3]) const +{ + return ( (i_xyz[0] >= 0) && (i_xyz[0] < n_cells[0]) + && (i_xyz[1] >= 0) && (i_xyz[1] < n_cells[1]) + && (i_xyz[2] >= 0) && (i_xyz[2] < n_cells[2])); +} + +//============================================================================== + +std::pair> +RectLattice::distance(const double xyz[3], const double uvw[3], + const int i_xyz[3]) const +{ + // Get short aliases to the coordinates. + double x {xyz[0]}; + double y {xyz[1]}; + double z {xyz[2]}; + double u {uvw[0]}; + double v {uvw[1]}; + + // Determine the oncoming edge. + double x0 {copysign(0.5 * pitch[0], u)}; + double y0 {copysign(0.5 * pitch[1], v)}; + + // Left and right sides + double d {INFTY}; + std::array lattice_trans; + if ((std::abs(x - x0) > FP_PRECISION) && u != 0) { + d = (x0 - x) / u; + if (u > 0) { + lattice_trans = {1, 0, 0}; + } else { + lattice_trans = {-1, 0, 0}; + } + } + + // Front and back sides + if ((std::abs(y - y0) > FP_PRECISION) && v != 0) { + double this_d = (y0 - y) / v; + if (this_d < d) { + d = this_d; + if (v > 0) { + lattice_trans = {0, 1, 0}; + } else { + lattice_trans = {0, -1, 0}; + } + } + } + + // Top and bottom sides + if (is_3d) { + double w {uvw[2]}; + double z0 {copysign(0.5 * pitch[2], w)}; + if ((std::abs(z - z0) > FP_PRECISION) && w != 0) { + double this_d = (z0 - z) / w; + if (this_d < d) { + d = this_d; + if (w > 0) { + lattice_trans = {0, 0, 1}; + } else { + lattice_trans = {0, 0, -1}; + } + } + } + } + + return {d, lattice_trans}; +} + +//============================================================================== + +std::array +RectLattice::get_indices(const double xyz[3]) const +{ + int ix {static_cast(std::ceil((xyz[0] - lower_left[0]) / pitch[0]))-1}; + int iy {static_cast(std::ceil((xyz[1] - lower_left[1]) / pitch[1]))-1}; + int iz; + if (is_3d) { + iz = static_cast(std::ceil((xyz[2] - lower_left[2]) / pitch[2]))-1; + } else { + iz = 0; + } + return {ix, iy, iz}; +} + +//============================================================================== + +std::array +RectLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const +{ + std::array local_xyz; + local_xyz[0] = global_xyz[0] - (lower_left[0] + (i_xyz[0] + 0.5)*pitch[0]); + local_xyz[1] = global_xyz[1] - (lower_left[1] + (i_xyz[1] + 0.5)*pitch[1]); + if (is_3d) { + local_xyz[2] = global_xyz[2] - (lower_left[2] + (i_xyz[2] + 0.5)*pitch[2]); + } else { + local_xyz[2] = global_xyz[2]; + } + return local_xyz; +} + +//============================================================================== + +int32_t& +RectLattice::offset(int map, const int i_xyz[3]) +{ + return offsets[nx*ny*nz*map + nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]]; +} + +//============================================================================== + +std::string +RectLattice::index_to_string(int indx) const +{ + int iz {indx / (nx * ny)}; + int iy {(indx - nx * ny * iz) / nx}; + int ix {indx - nx * ny * iz - nx * iy}; + std::string out {std::to_string(ix)}; + out += ','; + out += std::to_string(iy); + if (is_3d) { + out += ','; + out += std::to_string(iz); + } + return out; +} + +//============================================================================== + +void +RectLattice::to_hdf5_inner(hid_t lat_group) const +{ + // Write basic lattice information. + write_string(lat_group, "type", "rectangular", false); + if (is_3d) { + write_double(lat_group, "pitch", pitch, false); + write_double(lat_group, "lower_left", lower_left, false); + write_int(lat_group, "dimension", n_cells, false); + } else { + std::array pitch_short {{pitch[0], pitch[1]}}; + write_double(lat_group, "pitch", pitch_short, false); + std::array ll_short {{lower_left[0], lower_left[1]}}; + write_double(lat_group, "lower_left", ll_short, false); + std::array nc_short {{n_cells[0], n_cells[1]}}; + write_int(lat_group, "dimension", nc_short, false); + } + + // Write the universe ids. The convention here is to switch the ordering on + // the y-axis to match the way universes are input in a text file. + if (is_3d) { + hsize_t nx {static_cast(n_cells[0])}; + hsize_t ny {static_cast(n_cells[1])}; + hsize_t nz {static_cast(n_cells[2])}; + int out[nx*ny*nz]; + + for (int m = 0; m < nz; m++) { + for (int k = 0; k < ny; k++) { + for (int j = 0; j < nx; j++) { + int indx1 = nx*ny*m + nx*k + j; + int indx2 = nx*ny*m + nx*(ny-k-1) + j; + out[indx2] = global_universes[universes[indx1]]->id; + } + } + } + + hsize_t dims[3] {nz, ny, nx}; + write_int(lat_group, 3, dims, "universes", out, false); + + } else { + hsize_t nx {static_cast(n_cells[0])}; + hsize_t ny {static_cast(n_cells[1])}; + int out[nx*ny]; + + for (int k = 0; k < ny; k++) { + for (int j = 0; j < nx; j++) { + int indx1 = nx*k + j; + int indx2 = nx*(ny-k-1) + j; + out[indx2] = global_universes[universes[indx1]]->id; + } + } + + hsize_t dims[3] {1, ny, nx}; + write_int(lat_group, 3, dims, "universes", out, false); + } +} + +//============================================================================== +// HexLattice implementation +//============================================================================== + +HexLattice::HexLattice(pugi::xml_node lat_node) + : Lattice {lat_node} +{ + // Read the number of lattice cells in each dimension. + n_rings = stoi(get_node_value(lat_node, "n_rings")); + if (check_for_node(lat_node, "n_axial")) { + n_axial = stoi(get_node_value(lat_node, "n_axial")); + is_3d = true; + } else { + n_axial = 1; + is_3d = false; + } + + // Read the lattice center. + std::string center_str {get_node_value(lat_node, "center")}; + std::vector center_words {split(center_str)}; + if (is_3d && (center_words.size() != 3)) { + fatal_error("A hexagonal lattice with must have
" + "specified by 3 numbers."); + } else if (!is_3d && center_words.size() != 2) { + fatal_error("A hexagonal lattice without must have
" + "specified by 2 numbers."); + } + center[0] = stod(center_words[0]); + center[1] = stod(center_words[1]); + if (is_3d) {center[2] = stod(center_words[2]);} + + // Read the lattice pitches. + std::string pitch_str {get_node_value(lat_node, "pitch")}; + std::vector pitch_words {split(pitch_str)}; + if (is_3d && (pitch_words.size() != 2)) { + fatal_error("A hexagonal lattice with must have " + "specified by 2 numbers."); + } else if (!is_3d && (pitch_words.size() != 1)) { + fatal_error("A hexagonal lattice without must have " + "specified by 1 number."); + } + pitch[0] = stod(pitch_words[0]); + if (is_3d) {pitch[1] = stod(pitch_words[1]);} + + // Read the universes and make sure the correct number was specified. + int n_univ = (3*n_rings*n_rings - 3*n_rings + 1) * n_axial; + std::string univ_str {get_node_value(lat_node, "universes")}; + std::vector univ_words {split(univ_str)}; + if (univ_words.size() != n_univ) { + std::stringstream err_msg; + err_msg << "Expected " << n_univ + << " universes for a hexagonal lattice with " << n_rings + << " rings and " << n_axial << " axial levels" << " but " + << univ_words.size() << " were specified."; + fatal_error(err_msg); + } + + // Parse the universes. + // Universes in hexagonal lattices are stored in a manner that represents + // a skewed coordinate system: (x, alpha) rather than (x, y). There is + // no obvious, direct relationship between the order of universes in the + // input and the order that they will be stored in the skewed array so + // the following code walks a set of index values across the skewed array + // in a manner that matches the input order. Note that i_x = 0, i_a = 0 + // corresponds to the center of the hexagonal lattice. + + universes.resize((2*n_rings-1) * (2*n_rings-1) * n_axial, C_NONE); + int input_index = 0; + for (int m = 0; m < n_axial; m++) { + // Initialize lattice indecies. + int i_x = 1; + int i_a = n_rings - 1; + + // Map upper triangular region of hexagonal lattice which is found in the + // first n_rings-1 rows of the input. + for (int k = 0; k < n_rings-1; k++) { + // Walk the index to lower-left neighbor of last row start. + i_x -= 1; + + // Iterate over the input columns. + for (int j = 0; j < k+1; j++) { + int indx = (2*n_rings-1)*(2*n_rings-1) * m + + (2*n_rings-1) * (i_a+n_rings-1) + + (i_x+n_rings-1); + universes[indx] = stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } + + // Return the lattice index to the start of the current row. + i_x -= 2 * (k+1); + i_a += (k+1); + } + + // Map the middle square region of the hexagonal lattice which is found in + // the next 2*n_rings-1 rows of the input. + for (int k = 0; k < 2*n_rings-1; k++) { + if ((k % 2) == 0) { + // Walk the index to the lower-left neighbor of the last row start. + i_x -= 1; + } else { + // Walk the index to the lower-right neighbor of the last row start. + i_x += 1; + i_a -= 1; + } + + // Iterate over the input columns. + for (int j = 0; j < n_rings - (k % 2); j++) { + int indx = (2*n_rings-1)*(2*n_rings-1) * m + + (2*n_rings-1) * (i_a+n_rings-1) + + (i_x+n_rings-1); + universes[indx] = stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } + + // Return the lattice index to the start of the current row. + i_x -= 2*(n_rings - (k % 2)); + i_a += n_rings - (k % 2); + } + + // Map the lower triangular region of the hexagonal lattice. + for (int k = 0; k < n_rings-1; k++) { + // Walk the index to the lower-right neighbor of the last row start. + i_x += 1; + i_a -= 1; + + // Iterate over the input columns. + for (int j = 0; j < n_rings-k-1; j++) { + int indx = (2*n_rings-1)*(2*n_rings-1) * m + + (2*n_rings-1) * (i_a+n_rings-1) + + (i_x+n_rings-1); + universes[indx] = stoi(univ_words[input_index]); + input_index++; + // Walk the index to the right neighbor (which is not adjacent). + i_x += 2; + i_a -= 1; + } + + // Return lattice index to start of current row. + i_x -= 2*(n_rings - k - 1); + i_a += n_rings - k - 1; + } + } +} + +//============================================================================== + +int32_t& +HexLattice::operator[](const int i_xyz[3]) +{ + int indx = (2*n_rings-1)*(2*n_rings-1) * i_xyz[2] + + (2*n_rings-1) * i_xyz[1] + + i_xyz[0]; + return universes[indx]; +} + +//============================================================================== + +LatticeIter HexLattice::begin() +{return LatticeIter(*this, n_rings-1);} + +ReverseLatticeIter HexLattice::rbegin() +{return ReverseLatticeIter(*this, universes.size()-n_rings);} + +//============================================================================== + +bool +HexLattice::are_valid_indices(const int i_xyz[3]) const +{ + return ((i_xyz[0] >= 0) && (i_xyz[1] >= 0) && (i_xyz[2] >= 0) + && (i_xyz[0] < 2*n_rings-1) && (i_xyz[1] < 2*n_rings-1) + && (i_xyz[0] + i_xyz[1] > n_rings-2) + && (i_xyz[0] + i_xyz[1] < 3*n_rings-2) + && (i_xyz[2] < n_axial)); +} + +//============================================================================== + +std::pair> +HexLattice::distance(const double xyz[3], const double uvw[3], + const int i_xyz[3]) const +{ + // Compute the direction on the hexagonal basis. + double beta_dir = uvw[0] * std::sqrt(3.0) / 2.0 + uvw[1] / 2.0; + double gamma_dir = uvw[0] * std::sqrt(3.0) / 2.0 - uvw[1] / 2.0; + + // Note that hexagonal lattice distance calculations are performed + // using the particle's coordinates relative to the neighbor lattice + // cells, not relative to the particle's current cell. This is done + // because there is significant disagreement between neighboring cells + // on where the lattice boundary is due to finite precision issues. + + // Upper-right and lower-left sides. + double d {INFTY}; + std::array lattice_trans; + double edge = -copysign(0.5*pitch[0], beta_dir); // Oncoming edge + std::array xyz_t; + if (beta_dir > 0) { + const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } else { + const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } + double beta = xyz_t[0] * std::sqrt(3.0) / 2.0 + xyz_t[1] / 2.0; + if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { + d = (edge - beta) / beta_dir; + if (beta_dir > 0) { + lattice_trans = {1, 0, 0}; + } else { + lattice_trans = {-1, 0, 0}; + } + } + + // Lower-right and upper-left sides. + edge = -copysign(0.5*pitch[0], gamma_dir); + if (gamma_dir > 0) { + const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } else { + const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } + double gamma = xyz_t[0] * std::sqrt(3.0) / 2.0 - xyz_t[1] / 2.0; + if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { + double this_d = (edge - gamma) / gamma_dir; + if (this_d < d) { + if (gamma_dir > 0) { + lattice_trans = {1, -1, 0}; + } else { + lattice_trans = {-1, 1, 0}; + } + d = this_d; + } + } + + // Upper and lower sides. + edge = -copysign(0.5*pitch[0], uvw[1]); + if (uvw[1] > 0) { + const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } else { + const int i_xyz_t[3] {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; + xyz_t = get_local_xyz(xyz, i_xyz_t); + } + if ((std::abs(xyz_t[1] - edge) > FP_PRECISION) && uvw[1] != 0) { + double this_d = (edge - xyz_t[1]) / uvw[1]; + if (this_d < d) { + if (uvw[1] > 0) { + lattice_trans = {0, 1, 0}; + } else { + lattice_trans = {0, -1, 0}; + } + d = this_d; + } + } + + // Top and bottom sides + if (is_3d) { + double z {xyz[2]}; + double w {uvw[2]}; + double z0 {copysign(0.5 * pitch[1], w)}; + if ((std::abs(z - z0) > FP_PRECISION) && w != 0) { + double this_d = (z0 - z) / w; + if (this_d < d) { + d = this_d; + if (w > 0) { + lattice_trans = {0, 0, 1}; + } else { + lattice_trans = {0, 0, -1}; + } + d = this_d; + } + } + } + + return {d, lattice_trans}; +} + +//============================================================================== + +std::array +HexLattice::get_indices(const double xyz[3]) const +{ + // Offset the xyz by the lattice center. + double xyz_o[3] {xyz[0] - center[0], xyz[1] - center[1], xyz[2]}; + if (is_3d) {xyz_o[2] -= center[2];} + + // Index the z direction. + std::array out; + if (is_3d) { + out[2] = static_cast(std::ceil(xyz_o[2] / pitch[1] + 0.5 * n_axial))-1; + } else { + out[2] = 0; + } + + // Convert coordinates into skewed bases. The (x, alpha) basis is used to + // find the index of the global coordinates to within 4 cells. + double alpha = xyz_o[1] - xyz_o[0] / std::sqrt(3.0); + out[0] = static_cast(std::floor(xyz_o[0] + / (0.5*std::sqrt(3.0) * pitch[0]))); + out[1] = static_cast(std::floor(alpha / pitch[0])); + + // Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but + // the array is offset so that the indices never go below 0). + out[0] += n_rings-1; + out[1] += n_rings-1; + + // Calculate the (squared) distance between the particle and the centers of + // the four possible cells. Regular hexagonal tiles form a Voronoi + // tessellation so the xyz should be in the hexagonal cell that it is closest + // to the center of. This method is used over a method that uses the + // remainders of the floor divisions above because it provides better finite + // precision performance. Squared distances are used becasue they are more + // computationally efficient than normal distances. + int k {1}; + int k_min {1}; + double d_min {INFTY}; + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 2; j++) { + int i_xyz[3] {out[0] + j, out[1] + i, 0}; + std::array xyz_t = get_local_xyz(xyz, i_xyz); + double d = xyz_t[0]*xyz_t[0] + xyz_t[1]*xyz_t[1]; + if (d < d_min) { + d_min = d; + k_min = k; + } + k++; + } + } + + // Select the minimum squared distance which corresponds to the cell the + // coordinates are in. + if (k_min == 2) { + out[0] += 1; + } else if (k_min == 3) { + out[1] += 1; + } else if (k_min == 4) { + out[0] += 1; + out[1] += 1; + } + + return out; +} + +//============================================================================== + +std::array +HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const +{ + std::array local_xyz; + + // x_l = x_g - (center + pitch_x*cos(30)*index_x) + local_xyz[0] = global_xyz[0] - (center[0] + + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]); + // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) + local_xyz[1] = global_xyz[1] - (center[1] + + (i_xyz[1] - n_rings + 1) * pitch[0] + + (i_xyz[0] - n_rings + 1) * pitch[0] / 2.0); + if (is_3d) { + local_xyz[2] = global_xyz[2] - center[2] + + (0.5 * n_axial - i_xyz[2] - 0.5) * pitch[1]; + } else { + local_xyz[2] = global_xyz[2]; + } + + return local_xyz; +} + +//============================================================================== + +bool +HexLattice::is_valid_index(int indx) const +{ + int nx {2*n_rings - 1}; + int ny {2*n_rings - 1}; + int nz {n_axial}; + int iz = indx / (nx * ny); + int iy = (indx - nx*ny*iz) / nx; + int ix = indx - nx*ny*iz - nx*iy; + int i_xyz[3] {ix, iy, iz}; + return are_valid_indices(i_xyz); +} + +//============================================================================== + +int32_t& +HexLattice::offset(int map, const int i_xyz[3]) +{ + int nx {2*n_rings - 1}; + int ny {2*n_rings - 1}; + int nz {n_axial}; + return offsets[nx*ny*nz*map + nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]]; +} + +//============================================================================== + +std::string +HexLattice::index_to_string(int indx) const +{ + int nx {2*n_rings - 1}; + int ny {2*n_rings - 1}; + int iz {indx / (nx * ny)}; + int iy {(indx - nx * ny * iz) / nx}; + int ix {indx - nx * ny * iz - nx * iy}; + std::string out {std::to_string(ix - n_rings + 1)}; + out += ','; + out += std::to_string(iy - n_rings + 1); + if (is_3d) { + out += ','; + out += std::to_string(iz); + } + return out; +} + +//============================================================================== + +void +HexLattice::to_hdf5_inner(hid_t lat_group) const +{ + // Write basic lattice information. + write_string(lat_group, "type", "hexagonal", false); + write_int(lat_group, 0, nullptr, "n_rings", &n_rings, false); + write_int(lat_group, 0, nullptr, "n_axial", &n_axial, false); + if (is_3d) { + write_double(lat_group, "pitch", pitch, false); + write_double(lat_group, "center", center, false); + } else { + std::array pitch_short {{pitch[0]}}; + write_double(lat_group, "pitch", pitch_short, false); + std::array center_short {{center[0], center[1]}}; + write_double(lat_group, "center", center_short, false); + } + + // Write the universe ids. + hsize_t nx {static_cast(2*n_rings - 1)}; + hsize_t ny {static_cast(2*n_rings - 1)}; + hsize_t nz {static_cast(n_axial)}; + int out[nx*ny*nz]; + + for (int m = 0; m < nz; m++) { + for (int k = 0; k < ny; k++) { + for (int j = 0; j < nx; j++) { + int indx = nx*ny*m + nx*k + j; + if (j + k < n_rings - 1) { + // This array position is never used; put a -1 to indicate this. + out[indx] = -1; + } else if (j + k > 3*n_rings - 3) { + // This array position is never used; put a -1 to indicate this. + out[indx] = -1; + } else { + out[indx] = global_universes[universes[indx]]->id; + } + } + } + } + + hsize_t dims[3] {nz, ny, nx}; + write_int(lat_group, 3, dims, "universes", out, false); +} + +//============================================================================== +// Non-method functions +//============================================================================== + +extern "C" void +read_lattices(pugi::xml_node *node) +{ + for (pugi::xml_node lat_node : node->children("lattice")) { + lattices_c.push_back(new RectLattice(lat_node)); + } + for (pugi::xml_node lat_node : node->children("hex_lattice")) { + lattices_c.push_back(new HexLattice(lat_node)); + } + + // Fill the lattice map. + for (int i_lat = 0; i_lat < lattices_c.size(); i_lat++) { + int id = lattices_c[i_lat]->id; + auto in_map = lattice_map.find(id); + if (in_map == lattice_map.end()) { + lattice_map[id] = i_lat; + } else { + std::stringstream err_msg; + err_msg << "Two or more lattices use the same unique ID: " << id; + fatal_error(err_msg); + } + } +} + +//============================================================================== +// Fortran compatibility functions +//============================================================================== + +extern "C" { + Lattice* lattice_pointer(int lat_ind) {return lattices_c[lat_ind];} + + int32_t lattice_id(Lattice *lat) {return lat->id;} + + bool lattice_are_valid_indices(Lattice *lat, const int i_xyz[3]) + {return lat->are_valid_indices(i_xyz);} + + void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3], + const int i_xyz[3], double *d, int lattice_trans[3]) + { + std::pair> ld {lat->distance(xyz, uvw, i_xyz)}; + *d = ld.first; + lattice_trans[0] = ld.second[0]; + lattice_trans[1] = ld.second[1]; + lattice_trans[2] = ld.second[2]; + } + + void lattice_get_indices(Lattice *lat, const double xyz[3], int i_xyz[3]) + { + std::array inds = lat->get_indices(xyz); + i_xyz[0] = inds[0]; + i_xyz[1] = inds[1]; + i_xyz[2] = inds[2]; + } + + void lattice_get_local_xyz(Lattice *lat, const double global_xyz[3], + const int i_xyz[3], double local_xyz[3]) + { + std::array xyz = lat->get_local_xyz(global_xyz, i_xyz); + local_xyz[0] = xyz[0]; + local_xyz[1] = xyz[1]; + local_xyz[2] = xyz[2]; + } + + int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3]) + {return lat->offset(map, i_xyz);} + + int32_t lattice_outer(Lattice *lat) {return lat->outer;} + + void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);} + + int32_t lattice_universe(Lattice *lat, const int i_xyz[3]) + {return (*lat)[i_xyz];} +} + +} // namespace openmc diff --git a/src/lattice.h b/src/lattice.h new file mode 100644 index 0000000000..19749ce4f0 --- /dev/null +++ b/src/lattice.h @@ -0,0 +1,259 @@ +#ifndef LATTICE_H +#define LATTICE_H + +#include +#include +#include +#include +#include + +#include "constants.h" +#include "hdf5.h" +#include "pugixml/pugixml.hpp" + + +namespace openmc { + +//============================================================================== +// Module constants +//============================================================================== + +constexpr int32_t NO_OUTER_UNIVERSE{-1}; + +//============================================================================== +// Global variables +//============================================================================== + +class Lattice; +extern std::vector lattices_c; + +extern std::unordered_map lattice_map; + +//============================================================================== +//! \class Lattice +//! \brief Abstract type for ordered array of universes. +//============================================================================== + +class LatticeIter; +class ReverseLatticeIter; + +class Lattice +{ +public: + int32_t id; //!< Universe ID number + std::string name; //!< User-defined name + std::vector universes; //!< Universes filling each lattice tile + int32_t outer {NO_OUTER_UNIVERSE}; //!< Universe tiled outside the lattice + std::vector offsets; //!< Distribcell offset table + + explicit Lattice(pugi::xml_node lat_node); + + virtual ~Lattice() {} + + virtual int32_t& operator[](const int i_xyz[3]) = 0; + + virtual LatticeIter begin(); + LatticeIter end(); + + virtual ReverseLatticeIter rbegin(); + ReverseLatticeIter rend(); + + //! Convert internal universe values from IDs to indices using universe_map. + void adjust_indices(); + + //! Allocate offset table for distribcell. + void allocate_offset_table(int n_maps) + {offsets.resize(n_maps * universes.size(), C_NONE);} + + //! Populate the distribcell offset tables. + int32_t fill_offset_table(int32_t offset, int32_t target_univ_id, int map); + + //! \brief Check lattice indices. + //! @param i_xyz[3] The indices for a lattice tile. + //! @return true if the given indices fit within the lattice bounds. False + //! otherwise. + virtual bool are_valid_indices(const int i_xyz[3]) const = 0; + + //! \brief Find the next lattice surface crossing + //! @param xyz[3] A 3D Cartesian coordinate. + //! @param uvw[3] A 3D Cartesian direction. + //! @param i_xyz[3] The indices for a lattice tile. + //! @return The distance to the next crossing and an array indicating how the + //! lattice indices would change after crossing that boundary. + virtual std::pair> + distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const + = 0; + + //! \brief Find the lattice tile indices for a given point. + //! @param xyz[3] A 3D Cartesian coordinate. + //! @return An array containing the indices of a lattice tile. + virtual std::array get_indices(const double xyz[3]) const = 0; + + //! \brief Get coordinates local to a lattice tile. + //! @param global_xyz[3] A 3D Cartesian coordinate. + //! @param i_xyz[3] The indices for a lattice tile. + //! @return Local 3D Cartesian coordinates. + virtual std::array + get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const = 0; + + //! \brief Check flattened lattice index. + //! @param indx The index for a lattice tile. + //! @return true if the given index fit within the lattice bounds. False + //! otherwise. + virtual bool is_valid_index(int indx) const + {return (indx >= 0) && (indx < universes.size());} + + //! \brief Get the distribcell offset for a lattice tile. + //! @param The map index for the target cell. + //! @param i_xyz[3] The indices for a lattice tile. + //! @return Distribcell offset i.e. the largest instance number for the target + //! cell found in the geometry tree under this lattice tile. + virtual int32_t& offset(int map, const int i_xyz[3]) = 0; + + //! \brief Convert an array index to a useful human-readable string. + //! @param indx The index for a lattice tile. + //! @return A string representing the lattice tile. + virtual std::string index_to_string(int indx) const = 0; + + //! \brief Write lattice information to an HDF5 group. + //! @param group_id An HDF5 group id. + void to_hdf5(hid_t group_id) const; + +protected: + bool is_3d; //!< Has divisions along the z-axis? + + virtual void to_hdf5_inner(hid_t group_id) const = 0; +}; + +//============================================================================== +//! An iterator over lattice universes. +//============================================================================== + +class LatticeIter +{ +public: + int indx; //!< An index to a Lattice universes or offsets array. + + LatticeIter(Lattice &lat_, int indx_) + : lat(lat_), + indx(indx_) + {} + + bool operator==(const LatticeIter &rhs) {return (indx == rhs.indx);} + + bool operator!=(const LatticeIter &rhs) {return !(*this == rhs);} + + int32_t& operator*() {return lat.universes[indx];} + + LatticeIter& operator++() + { + while (indx < lat.universes.size()) { + ++indx; + if (lat.is_valid_index(indx)) return *this; + } + indx = lat.universes.size(); + return *this; + } + +protected: + Lattice ⪫ +}; + +//============================================================================== +//! A reverse iterator over lattice universes. +//============================================================================== + +class ReverseLatticeIter : public LatticeIter +{ +public: + ReverseLatticeIter(Lattice &lat_, int indx_) + : LatticeIter {lat_, indx_} + {} + + ReverseLatticeIter& operator++() + { + while (indx > -1) { + --indx; + if (lat.is_valid_index(indx)) return *this; + } + indx = -1; + return *this; + } +}; + +//============================================================================== + +class RectLattice : public Lattice +{ +public: + explicit RectLattice(pugi::xml_node lat_node); + + int32_t& operator[](const int i_xyz[3]); + + bool are_valid_indices(const int i_xyz[3]) const; + + std::pair> + distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const; + + std::array get_indices(const double xyz[3]) const; + + std::array + get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const; + + int32_t& offset(int map, const int i_xyz[3]); + + std::string index_to_string(int indx) const; + + void to_hdf5_inner(hid_t group_id) const; + +private: + std::array n_cells; //!< Number of cells along each axis + std::array lower_left; //!< Global lower-left corner of the lattice + std::array pitch; //!< Lattice tile width along each axis + + // Convenience aliases + int &nx {n_cells[0]}; + int &ny {n_cells[1]}; + int &nz {n_cells[2]}; +}; + +//============================================================================== + +class HexLattice : public Lattice +{ +public: + explicit HexLattice(pugi::xml_node lat_node); + + int32_t& operator[](const int i_xyz[3]); + + LatticeIter begin(); + + ReverseLatticeIter rbegin(); + + bool are_valid_indices(const int i_xyz[3]) const; + + std::pair> + distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const; + + std::array get_indices(const double xyz[3]) const; + + std::array + get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const; + + bool is_valid_index(int indx) const; + + int32_t& offset(int map, const int i_xyz[3]); + + std::string index_to_string(int indx) const; + + void to_hdf5_inner(hid_t group_id) const; + +private: + int n_rings; //!< Number of radial tile positions + int n_axial; //!< Number of axial tile positions + std::array center; //!< Global center of lattice + std::array pitch; //!< Lattice tile width and height +}; + +} // namespace openmc +#endif // LATTICE_H diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index ebabe5cdbc..17880f6e1b 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -10,7 +10,6 @@ module nuclide_header use endf_header, only: Function1D, Polynomial, Tabulated1D use error use hdf5_interface - use list_header, only: ListInt use math, only: faddeeva, w_derivative, & broaden_wmp_polynomials use multipole_header, only: FORM_RM, FORM_MLBW, MP_EA, RM_RT, RM_RA, & diff --git a/src/output.F90 b/src/output.F90 index 23d9512eee..b1671df984 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -234,7 +234,7 @@ contains ! Print cell for this level if (p % coord(i) % cell /= NONE) then c => cells(p % coord(i) % cell) - write(ou,*) ' Cell = ' // trim(to_str(c % id)) + write(ou,*) ' Cell = ' // trim(to_str(c % id())) end if ! Print universe for this level @@ -246,7 +246,7 @@ contains ! Print information on lattice if (p % coord(i) % lattice /= NONE) then l => lattices(p % coord(i) % lattice) % obj - write(ou,*) ' Lattice = ' // trim(to_str(l % id)) + write(ou,*) ' Lattice = ' // trim(to_str(l % id())) write(ou,*) ' Lattice position = (' // trim(to_str(& p % coord(i) % lattice_x)) // ',' // trim(to_str(& p % coord(i) % lattice_y)) // ')' @@ -258,7 +258,7 @@ contains end do ! Print surface - if (p % surface /= NONE) then + if (p % surface /= ERROR_INT) then write(ou,*) ' Surface = ' // to_str(sign(surfaces(i)%id(), p % surface)) end if @@ -633,7 +633,7 @@ contains write(ou,100) 'Cell ID','No. Overlap Checks' do i = 1, n_cells - write(ou,101) cells(i) % id, overlap_check_cnt(i) + write(ou,101) cells(i) % id(), overlap_check_cnt(i) if (overlap_check_cnt(i) < 10) num_sparse = num_sparse + 1 end do write(ou,*) @@ -643,7 +643,7 @@ contains do i = 1, n_cells if (overlap_check_cnt(i) < 10) then j = j + 1 - write(ou,'(1X,A8)', advance='no') trim(to_str(cells(i) % id)) + write(ou,'(1X,A8)', advance='no') trim(to_str(cells(i) % id())) if (modulo(j,8) == 0) write(ou,*) end if end do diff --git a/src/plot.F90 b/src/plot.F90 index f17ac33511..789ad0ee9b 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 @@ -101,7 +101,7 @@ contains else if (pl % color_by == PLOT_COLOR_CELLS) then ! Assign color based on cell rgb = pl % colors(p % coord(j) % cell) % rgb - id = cells(p % coord(j) % cell) % id + id = cells(p % coord(j) % cell) % id() else rgb = 0 id = -1 diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 8fe5e2da1f..806bcc65e8 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -72,9 +72,6 @@ module simulation_header logical :: trace - ! Number of distribcell maps - integer :: n_maps - !$omp threadprivate(trace, thread_id, current_work) contains diff --git a/src/string_utils.h b/src/string_utils.h new file mode 100644 index 0000000000..e2cb773ab1 --- /dev/null +++ b/src/string_utils.h @@ -0,0 +1,35 @@ +#ifndef STRING_UTILS_H +#define STRING_UTILS_H + +#include +#include + + +namespace openmc { + +std::vector +split(const std::string &in) +{ + std::vector out; + + for (int i = 0; i < in.size(); ) { + // Increment i until we find a non-whitespace character. + if (std::isspace(in[i])) { + i++; + + } else { + // Find the next whitespace character at j. + int j = i + 1; + while (j < in.size() && std::isspace(in[j]) == 0) {j++;} + + // Push-back everything between i and j. + out.push_back(in.substr(i, j-i)); + i = j + 1; // j is whitespace so leapfrog to j+1 + } + } + + return out; +} + +} // namespace openmc +#endif // STRING_UTILS_H diff --git a/src/summary.F90 b/src/summary.F90 index 45c92062d5..940e3128ab 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -15,7 +15,6 @@ module summary use surface_header use string, only: to_str use tally_header, only: TallyObject - use tally_filter_distribcell, only: find_offset implicit none private @@ -160,8 +159,7 @@ contains subroutine write_geometry(file_id) integer(HID_T), intent(in) :: file_id - integer :: i, j, k, m - integer, allocatable :: lattice_universes(:,:,:) + integer :: i, j, cell_fill integer, allocatable :: cell_materials(:) integer, allocatable :: cell_ids(:) real(8), allocatable :: cell_temperatures(:) @@ -169,8 +167,7 @@ contains integer(HID_T) :: cells_group, cell_group integer(HID_T) :: surfaces_group integer(HID_T) :: universes_group, univ_group - integer(HID_T) :: lattices_group, lattice_group - character(:), allocatable :: region_spec + integer(HID_T) :: lattices_group type(Cell), pointer :: c class(Lattice), pointer :: lat @@ -179,7 +176,7 @@ contains call write_attribute(geom_group, "n_cells", n_cells) call write_attribute(geom_group, "n_surfaces", n_surfaces) call write_attribute(geom_group, "n_universes", n_universes) - call write_attribute(geom_group, "n_lattices", n_lattices) + call write_attribute(geom_group, "n_lattices", size(lattices)) ! ========================================================================== ! WRITE INFORMATION ON CELLS @@ -190,16 +187,12 @@ contains ! Write information on each cell CELL_LOOP: do i = 1, n_cells c => cells(i) - cell_group = create_group(cells_group, "cell " // trim(to_str(c%id))) + cell_group = create_group(cells_group, "cell " // trim(to_str(c%id()))) - ! Write name for this cell - call write_dataset(cell_group, "name", c%name) - - ! Write universe for this cell - call write_dataset(cell_group, "universe", universes(c%universe)%id) + 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") @@ -231,12 +224,7 @@ contains case (FILL_UNIVERSE) call write_dataset(cell_group, "fill_type", "universe") - call write_dataset(cell_group, "fill", universes(c%fill)%id) - if (allocated(c%offset)) then - if (size(c%offset) > 0) then - call write_dataset(cell_group, "offset", c%offset) - end if - end if + call write_dataset(cell_group, "fill", universes(c%fill()+1)%id) if (allocated(c%translation)) then call write_dataset(cell_group, "translation", c%translation) @@ -247,30 +235,12 @@ contains case (FILL_LATTICE) call write_dataset(cell_group, "fill_type", "lattice") - call write_dataset(cell_group, "lattice", lattices(c%fill)%obj%id) + ! Do not access the 'lattices' array with 'c % fill() + 1' directly; it + ! causes a segfault in GCC 7.3.0 + cell_fill = c % fill() + 1 + call write_dataset(cell_group, "lattice", lattices(cell_fill)%obj%id()) end select - ! Write list of bounding surfaces - region_spec = "" - do j = 1, size(c%region) - k = c%region(j) - select case(k) - case (OP_LEFT_PAREN) - region_spec = trim(region_spec) // " (" - case (OP_RIGHT_PAREN) - region_spec = trim(region_spec) // " )" - case (OP_COMPLEMENT) - region_spec = trim(region_spec) // " ~" - case (OP_INTERSECTION) - case (OP_UNION) - region_spec = trim(region_spec) // " |" - case default - region_spec = trim(region_spec) // " " // to_str(& - sign(surfaces(abs(k))%id(), k)) - end select - end do - call write_dataset(cell_group, "region", adjustl(region_spec)) - call close_group(cell_group) end do CELL_LOOP @@ -305,7 +275,7 @@ contains if (size(u % cells) > 0) then allocate(cell_ids(size(u % cells))) do j = 1, size(u % cells) - cell_ids(j) = cells(u % cells(j)) % id + cell_ids(j) = cells(u % cells(j)) % id() end do call write_dataset(univ_group, "cells", cell_ids) deallocate(cell_ids) @@ -320,87 +290,12 @@ contains ! ========================================================================== ! WRITE INFORMATION ON LATTICES - ! Create lattices group (nothing directly written here) then close lattices_group = create_group(geom_group, "lattices") - ! Write information on each lattice - LATTICE_LOOP: do i = 1, n_lattices + do i = 1, size(lattices) lat => lattices(i)%obj - lattice_group = create_group(lattices_group, "lattice " // trim(to_str(lat%id))) - - ! Write name, pitch, and outer universe - call write_dataset(lattice_group, "name", lat%name) - call write_dataset(lattice_group, "pitch", lat%pitch) - if (lat % outer > 0) then - call write_dataset(lattice_group, "outer", universes(lat % outer) % id) - else - call write_dataset(lattice_group, "outer", lat % outer) - end if - - select type (lat) - type is (RectLattice) - ! Write lattice type. - call write_dataset(lattice_group, "type", "rectangular") - - ! Write lattice dimensions, lower left corner, and pitch - if (lat % is_3d) then - call write_dataset(lattice_group, "dimension", lat % n_cells) - call write_dataset(lattice_group, "lower_left", lat % lower_left) - else - call write_dataset(lattice_group, "dimension", lat % n_cells(1:2)) - call write_dataset(lattice_group, "lower_left", lat % lower_left) - end if - - ! Write lattice universes. - allocate(lattice_universes(lat%n_cells(1), lat%n_cells(2), & - &lat%n_cells(3))) - do j = 1, lat%n_cells(1) - do k = 0, lat%n_cells(2) - 1 - do m = 1, lat%n_cells(3) - lattice_universes(j, k+1, m) = & - universes(lat%universes(j, lat%n_cells(2) - k, m))%id - end do - end do - end do - - type is (HexLattice) - ! Write lattice type. - call write_dataset(lattice_group, "type", "hexagonal") - - ! Write number of lattice cells. - call write_dataset(lattice_group, "n_rings", lat%n_rings) - call write_dataset(lattice_group, "n_axial", lat%n_axial) - - ! Write lattice center - call write_dataset(lattice_group, "center", lat%center) - - ! Write lattice universes. - allocate(lattice_universes(2*lat%n_rings - 1, 2*lat%n_rings - 1, & - &lat%n_axial)) - do m = 1, lat%n_axial - do k = 1, 2*lat%n_rings - 1 - do j = 1, 2*lat%n_rings - 1 - if (j + k < lat%n_rings + 1) then - ! This array position is never used; put a -1 to indicate this - lattice_universes(j,k,m) = -1 - cycle - else if (j + k > 3*lat%n_rings - 1) then - ! This array position is never used; put a -1 to indicate this - lattice_universes(j,k,m) = -1 - cycle - end if - lattice_universes(j,k,m) = universes(lat%universes(j,k,m))%id - end do - end do - end do - end select - - ! Write lattice universes - call write_dataset(lattice_group, "universes", lattice_universes) - deallocate(lattice_universes) - - call close_group(lattice_group) - end do LATTICE_LOOP + call lat % to_hdf5(lattices_group) + end do call close_group(lattices_group) call close_group(geom_group) diff --git a/src/surface.cpp b/src/surface.cpp index 2c3279604a..56e57a854f 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -11,12 +11,25 @@ namespace openmc { +//============================================================================== +// Module constant definitions +//============================================================================== + +extern "C" const int BC_TRANSMIT {0}; +extern "C" const int BC_VACUUM {1}; +extern "C" const int BC_REFLECT {2}; +extern "C" const int BC_PERIODIC {3}; + //============================================================================== // Global variables //============================================================================== int32_t n_surfaces; +Surface **surfaces_c; + +std::map surface_map; + //============================================================================== // Helper functions for reading the "coeffs" node of an XML surface element //============================================================================== @@ -305,7 +318,7 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-plane", false); std::array coeffs {{x0}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -370,7 +383,7 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "y-plane", false); std::array coeffs {{y0}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -436,7 +449,7 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "z-plane", false); std::array coeffs {{z0}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -497,7 +510,7 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "plane", false); std::array coeffs {{A, B, C, D}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3], @@ -629,7 +642,7 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-cylinder", false); std::array coeffs {{y0, z0, r}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -663,7 +676,7 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "y-cylinder", false); std::array coeffs {{x0, z0, r}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -697,7 +710,7 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "z-cylinder", false); std::array coeffs {{x0, y0, r}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -768,7 +781,7 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "sphere", false); std::array coeffs {{x0, y0, z0, r}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -885,7 +898,7 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-cone", false); std::array coeffs {{x0, y0, z0, r_sq}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -919,7 +932,7 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "y-cone", false); std::array coeffs {{x0, y0, z0, r_sq}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -953,7 +966,7 @@ void SurfaceZCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "z-cone", false); std::array coeffs {{x0, y0, z0, r_sq}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -1047,7 +1060,7 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "quadric", false); std::array coeffs {{A, B, C, D, E, F, G, H, J, K}}; - write_double_1D(group_id, "coefficients", coeffs); + write_double(group_id, "coefficients", coeffs, false); } //============================================================================== @@ -1116,12 +1129,12 @@ read_surfaces(pugi::xml_node *node) } } - // Fill the surface dictionary. + // Fill the surface map. for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { int id = surfaces_c[i_surf]->id; - auto in_dict = surface_dict.find(id); - if (in_dict == surface_dict.end()) { - surface_dict[id] = i_surf; + auto in_map = surface_map.find(id); + if (in_map == surface_map.end()) { + surface_map[id] = i_surf; } else { std::stringstream err_msg; err_msg << "Two or more surfaces use the same unique ID: " << id; @@ -1211,7 +1224,7 @@ read_surfaces(pugi::xml_node *node) } } else { // Convert the surface id to an index. - surf->i_periodic = surface_dict[surf->i_periodic]; + surf->i_periodic = surface_map[surf->i_periodic]; } } else { // This is a SurfacePlane. We won't try to find it's partner if the @@ -1223,7 +1236,7 @@ read_surfaces(pugi::xml_node *node) fatal_error(err_msg); } else { // Convert the surface id to an index. - surf->i_periodic = surface_dict[surf->i_periodic]; + surf->i_periodic = surface_map[surf->i_periodic]; } } @@ -1242,43 +1255,36 @@ read_surfaces(pugi::xml_node *node) // Fortran compatibility functions //============================================================================== -extern "C" Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];} +extern "C" { + Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];} -extern "C" int surface_id(Surface *surf) {return surf->id;} + int surface_id(Surface *surf) {return surf->id;} -extern "C" int surface_bc(Surface *surf) {return surf->bc;} + int surface_bc(Surface *surf) {return surf->bc;} -extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3]) -{return surf->sense(xyz, uvw);} + void surface_reflect(Surface *surf, double xyz[3], double uvw[3]) + {surf->reflect(xyz, uvw);} -extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3]) -{surf->reflect(xyz, uvw);} + void surface_normal(Surface *surf, double xyz[3], double uvw[3]) + {return surf->normal(xyz, uvw);} -extern "C" double -surface_distance(Surface *surf, double xyz[3], double uvw[3], bool coincident) -{return surf->distance(xyz, uvw, coincident);} + void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);} -extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3]) -{return surf->normal(xyz, uvw);} + int surface_i_periodic(PeriodicSurface *surf) {return surf->i_periodic;} -extern "C" void surface_to_hdf5(Surface *surf, hid_t group) -{surf->to_hdf5(group);} + bool + surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3], + double uvw[3]) + {return surf->periodic_translate(other, xyz, uvw);} -extern "C" int surface_i_periodic(PeriodicSurface *surf) -{return surf->i_periodic;} - -extern "C" bool -surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3], - double uvw[3]) -{return surf->periodic_translate(other, xyz, uvw);} - -extern "C" void free_memory_surfaces_c() -{ - for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];} - delete surfaces_c; - surfaces_c = nullptr; - n_surfaces = 0; - surface_dict.clear(); + void free_memory_surfaces_c() + { + for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];} + delete surfaces_c; + surfaces_c = nullptr; + n_surfaces = 0; + surface_map.clear(); + } } } // namespace openmc diff --git a/src/surface.h b/src/surface.h index 15ac0d5bdd..1d8cc6ac7b 100644 --- a/src/surface.h +++ b/src/surface.h @@ -8,25 +8,19 @@ #include "hdf5.h" #include "pugixml/pugixml.hpp" +#include "constants.h" + namespace openmc { //============================================================================== -// Module constants +// Module constant declarations (defined in .cpp) //============================================================================== -extern "C" const int BC_TRANSMIT {0}; -extern "C" const int BC_VACUUM {1}; -extern "C" const int BC_REFLECT {2}; -extern "C" const int BC_PERIODIC {3}; - -//============================================================================== -// Constants that should eventually be moved out of this file -//============================================================================== - -extern "C" double FP_COINCIDENT; -constexpr double INFTY {std::numeric_limits::max()}; -constexpr int C_NONE {-1}; +extern "C" const int BC_TRANSMIT; +extern "C" const int BC_VACUUM; +extern "C" const int BC_REFLECT; +extern "C" const int BC_PERIODIC; //============================================================================== // Global variables @@ -35,9 +29,9 @@ constexpr int C_NONE {-1}; extern "C" int32_t n_surfaces; class Surface; -Surface **surfaces_c; +extern Surface **surfaces_c; -std::map surface_dict; +extern std::map surface_map; //============================================================================== //! Coordinates for an axis-aligned cube that bounds a geometric object. @@ -64,7 +58,7 @@ public: //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::string name; //!< User-defined name explicit Surface(pugi::xml_node surf_node); @@ -107,6 +101,7 @@ public: //! Write all information needed to reconstruct the surface to an HDF5 group. //! @param group_id An HDF5 group id. + //TODO: this probably needs to include i_periodic for PeriodicSurface void to_hdf5(hid_t group_id) const; protected: @@ -383,19 +378,21 @@ public: // Fortran compatibility functions //============================================================================== -extern "C" Surface* surface_pointer(int surf_ind); -extern "C" int surface_id(Surface *surf); -extern "C" int surface_bc(Surface *surf); -extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3]); -extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3]); -extern "C" double surface_distance(Surface *surf, double xyz[3], double uvw[3], - bool coincident); -extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3]); -extern "C" void surface_to_hdf5(Surface *surf, hid_t group); -extern "C" int surface_i_periodic(PeriodicSurface *surf); -extern "C" bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, - double xyz[3], double uvw[3]); -extern "C" void free_memory_surfaces_c(); +extern "C" { + Surface* surface_pointer(int surf_ind); + int surface_id(Surface *surf); + int surface_bc(Surface *surf); + bool surface_sense(Surface *surf, double xyz[3], double uvw[3]); + void surface_reflect(Surface *surf, double xyz[3], double uvw[3]); + double surface_distance(Surface *surf, double xyz[3], double uvw[3], + bool coincident); + void surface_normal(Surface *surf, double xyz[3], double uvw[3]); + void surface_to_hdf5(Surface *surf, hid_t group); + int surface_i_periodic(PeriodicSurface *surf); + bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, + double xyz[3], double uvw[3]); + void free_memory_surfaces_c(); +} } // namespace openmc #endif // SURFACE_H diff --git a/src/surface_header.F90 b/src/surface_header.F90 index ab6babb5cb..916f9c11c9 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -30,16 +30,6 @@ module surface_header integer(C_INT) :: bc end function surface_bc_c - pure function surface_sense_c(surf_ptr, xyz, uvw) & - bind(C, name='surface_sense') result(sense) - use ISO_C_BINDING - implicit none - type(C_PTR), intent(in), value :: surf_ptr - real(C_DOUBLE), intent(in) :: xyz(3) - real(C_DOUBLE), intent(in) :: uvw(3) - logical(C_BOOL) :: sense - end function surface_sense_c - pure subroutine surface_reflect_c(surf_ptr, xyz, uvw) & bind(C, name='surface_reflect') use ISO_C_BINDING @@ -49,17 +39,6 @@ module surface_header real(C_DOUBLE), intent(inout) :: uvw(3); end subroutine surface_reflect_c - pure function surface_distance_c(surf_ptr, xyz, uvw, coincident) & - bind(C, name='surface_distance') result(d) - use ISO_C_BINDING - implicit none - type(C_PTR), intent(in), value :: surf_ptr - real(C_DOUBLE), intent(in) :: xyz(3); - real(C_DOUBLE), intent(in) :: uvw(3); - logical(C_BOOL), intent(in), value :: coincident; - real(C_DOUBLE) :: d; - end function surface_distance_c - pure subroutine surface_normal_c(surf_ptr, xyz, uvw) & bind(C, name='surface_normal') use ISO_C_BINDING @@ -115,9 +94,7 @@ module surface_header procedure :: id => surface_id procedure :: bc => surface_bc - procedure :: sense => surface_sense procedure :: reflect => surface_reflect - procedure :: distance => surface_distance procedure :: normal => surface_normal procedure :: to_hdf5 => surface_to_hdf5 procedure :: i_periodic => surface_i_periodic @@ -146,14 +123,6 @@ contains bc = surface_bc_c(this % ptr) end function surface_bc - pure function surface_sense(this, xyz, uvw) result(sense) - class(Surface), intent(in) :: this - real(C_DOUBLE), intent(in) :: xyz(3) - real(C_DOUBLE), intent(in) :: uvw(3) - logical(C_BOOL) :: sense - sense = surface_sense_c(this % ptr, xyz, uvw) - end function surface_sense - pure subroutine surface_reflect(this, xyz, uvw) class(Surface), intent(in) :: this real(C_DOUBLE), intent(in) :: xyz(3); @@ -161,15 +130,6 @@ contains call surface_reflect_c(this % ptr, xyz, uvw) end subroutine surface_reflect - pure function surface_distance(this, xyz, uvw, coincident) result(d) - class(Surface), intent(in) :: this - real(C_DOUBLE), intent(in) :: xyz(3); - real(C_DOUBLE), intent(in) :: uvw(3); - logical(C_BOOL), intent(in) :: coincident; - real(C_DOUBLE) :: d; - d = surface_distance_c(this % ptr, xyz, uvw, coincident) - end function surface_distance - pure subroutine surface_normal(this, xyz, uvw) class(Surface), intent(in) :: this real(C_DOUBLE), intent(in) :: xyz(3); diff --git a/src/tallies/tally_filter_cell.F90 b/src/tallies/tally_filter_cell.F90 index 99a9892a77..be37467073 100644 --- a/src/tallies/tally_filter_cell.F90 +++ b/src/tallies/tally_filter_cell.F90 @@ -80,7 +80,7 @@ contains allocate(cell_ids(size(this % cells))) do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id + cell_ids(i) = cells(this % cells(i)) % id() end do call write_dataset(filter_group, "bins", cell_ids) end subroutine to_statepoint_cell @@ -114,7 +114,7 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - label = "Cell " // to_str(cells(this % cells(bin)) % id) + label = "Cell " // to_str(cells(this % cells(bin)) % id()) end function text_label_cell !=============================================================================== diff --git a/src/tallies/tally_filter_cellborn.F90 b/src/tallies/tally_filter_cellborn.F90 index 0d7461ef82..a373c2d4e1 100644 --- a/src/tallies/tally_filter_cellborn.F90 +++ b/src/tallies/tally_filter_cellborn.F90 @@ -74,7 +74,7 @@ contains call write_dataset(filter_group, "n_bins", this % n_bins) allocate(cell_ids(size(this % cells))) do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id + cell_ids(i) = cells(this % cells(i)) % id() end do call write_dataset(filter_group, "bins", cell_ids) end subroutine to_statepoint_cellborn @@ -108,7 +108,7 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - label = "Birth Cell " // to_str(cells(this % cells(bin)) % id) + label = "Birth Cell " // to_str(cells(this % cells(bin)) % id()) end function text_label_cellborn end module tally_filter_cellborn diff --git a/src/tallies/tally_filter_cellfrom.F90 b/src/tallies/tally_filter_cellfrom.F90 index 2c54021226..83a3179d26 100644 --- a/src/tallies/tally_filter_cellfrom.F90 +++ b/src/tallies/tally_filter_cellfrom.F90 @@ -64,7 +64,7 @@ contains allocate(cell_ids(size(this % cells))) do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id + cell_ids(i) = cells(this % cells(i)) % id() end do call write_dataset(filter_group, "bins", cell_ids) end subroutine to_statepoint_cell_from @@ -74,7 +74,7 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - label = "Cell from " // to_str(cells(this % cells(bin)) % id) + label = "Cell from " // to_str(cells(this % cells(bin)) % id()) end function text_label_cell_from end module tally_filter_cellfrom diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index ddf7f4d922..c9564083da 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -14,7 +14,6 @@ module tally_filter_distribcell implicit none private - public :: find_offset !=============================================================================== ! DISTRIBCELLFILTER specifies which distributed geometric cells tally events @@ -58,20 +57,20 @@ 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 - offset = offset + cells(p % coord(i) % cell) % & - offset(distribcell_index) - elseif (cells(p % coord(i) % cell) % type == FILL_LATTICE) then + if (cells(p % coord(i) % cell) % type() == FILL_UNIVERSE) then + offset = offset + cells(p % coord(i) % cell) & + % offset(distribcell_index-1) + 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, & p % coord(i + 1) % lattice_y, & p % coord(i + 1) % lattice_z])) then - offset = offset + lattices(p % coord(i + 1) % lattice) % obj % & - offset(distribcell_index, & - p % coord(i + 1) % lattice_x, & + offset = offset + lattices(p % coord(i + 1) % lattice) % obj & + % offset(distribcell_index - 1, & + [p % coord(i + 1) % lattice_x, & p % coord(i + 1) % lattice_y, & - p % coord(i + 1) % lattice_z) + p % coord(i + 1) % lattice_z]) end if end if if (this % cell == p % coord(i) % cell) then @@ -88,7 +87,7 @@ contains call write_dataset(filter_group, "type", "distribcell") call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", cells(this % cell) % id) + call write_dataset(filter_group, "bins", cells(this % cell) % id()) end subroutine to_statepoint_distribcell subroutine initialize_distribcell(this) @@ -102,7 +101,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.") @@ -114,13 +113,8 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - integer :: offset - type(Universe), pointer :: univ - - univ => universes(root_universe) - offset = 0 label = '' - call find_offset(this % cell, univ, bin-1, offset, label) + call find_offset(this % cell, bin-1, label) label = "Distributed Cell " // label end function text_label_distribcell @@ -130,279 +124,46 @@ contains ! the target cell with the given offset !=============================================================================== - recursive subroutine find_offset(i_cell, univ, target_offset, offset, path) - + subroutine find_offset(i_cell, target_offset, path) integer, intent(in) :: i_cell ! The target cell index - type(Universe), intent(in) :: univ ! Universe to begin search - integer, intent(in) :: target_offset ! Target offset - integer, intent(inout) :: offset ! Current offset - character(*), intent(inout) :: path ! Path to offset + integer, intent(in) :: target_offset ! Target offset + character(*), intent(inout) :: path ! Path to offset integer :: map ! Index in maps vector - integer :: i, j ! Index over cells - integer :: k, l, m ! Indices in lattice - integer :: old_k, old_l, old_m ! Previous indices in lattice - integer :: n_x, n_y, n_z ! Lattice cell array dimensions - integer :: n ! Number of cells to search - integer :: cell_index ! Index in cells array - integer :: lat_offset ! Offset from lattice - integer :: temp_offset ! Looped sum of offsets - integer :: i_univ ! index in universes array - logical :: this_cell = .false. ! Advance in this cell? - logical :: later_cell = .false. ! Fill cells after this one? - type(Cell), pointer :: c ! Pointer to current cell - type(Universe), pointer :: next_univ ! Next universe to loop through - class(Lattice), pointer :: lat ! Pointer to current lattice + integer :: i ! Index over cells + + integer(C_INT) :: path_len + character(kind=C_CHAR), allocatable, target :: path_c(:) + + interface + function distribcell_path_len(target_cell, map, target_offset, root_univ)& + bind(C) result(len) + import C_INT32_T, C_INT + integer(C_INT32_T), intent(in), value :: target_cell, map, & + target_offset, root_univ + integer(C_INT) :: len + end function distribcell_path_len + + subroutine distribcell_path(target_cell, map, target_offset, root_univ, & + path) bind(C) + import C_INT32_T, C_CHAR + integer(C_INT32_T), intent(in), value :: target_cell, map, & + target_offset, root_univ + character(kind=C_CHAR), intent(out) :: path(*) + end subroutine distribcell_path + end interface ! Get the distribcell index for this cell map = cells(i_cell) % distribcell_index - n = size(univ % cells) - - ! Write to the geometry stack - i_univ = universe_dict % get(univ % id) - if (i_univ == root_universe) then - path = trim(path) // "u" // to_str(univ%id) - else - path = trim(path) // "->u" // to_str(univ%id) - end if - - ! Look through all cells in this universe - do i = 1, n - ! If the cell matches the goal and the offset matches final, write to the - ! geometry stack - if (univ % cells(i) == i_cell .and. offset == target_offset) then - c => cells(univ % cells(i)) - path = trim(path) // "->c" // to_str(c % id) - return - end if - end do - - ! Find the fill cell or lattice cell that we need to enter - do i = 1, n - - later_cell = .false. - - c => cells(univ % cells(i)) - - this_cell = .false. - - ! If we got here, we still think the target is in this universe - ! or further down, but it's not this exact cell. - ! Compare offset to next cell to see if we should enter this cell - if (i /= n) then - - do j = i+1, n - - c => cells(univ % cells(j)) - - ! Skip normal cells which do not have offsets - 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 - - ! 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 - temp_offset = c % offset(map) - - ! Get the offset of the first lattice location - else - lat => lattices(c % fill) % obj - temp_offset = lat % offset(map, 1, 1, 1) - end if - - ! If the final offset is in the range of offset - temp_offset+offset - ! then the goal is in this cell - if (target_offset < temp_offset + offset) then - this_cell = .true. - end if - end if - end if - - if (n == 1 .and. c % type /= FILL_MATERIAL) then - this_cell = .true. - end if - - if (.not. later_cell) then - this_cell = .true. - end if - - ! Get pointer to THIS cell because target must be in this cell - if (this_cell) then - - cell_index = univ % cells(i) - c => cells(cell_index) - - path = trim(path) // "->c" // to_str(c%id) - - ! ==================================================================== - ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - if (c % type == FILL_UNIVERSE) then - - ! Enter this cell to update the current offset - offset = c % offset(map) + offset - - next_univ => universes(c % fill) - call find_offset(i_cell, next_univ, target_offset, offset, path) - return - - ! ==================================================================== - ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then - - ! Set current lattice - lat => lattices(c % fill) % obj - - select type (lat) - - ! ================================================================== - ! RECTANGULAR LATTICES - type is (RectLattice) - - ! Write to the geometry stack - path = trim(path) // "->l" // to_str(lat%id) - - n_x = lat % n_cells(1) - n_y = lat % n_cells(2) - n_z = lat % n_cells(3) - old_m = 1 - old_l = 1 - old_k = 1 - - ! Loop over lattice coordinates - do m = 1, n_z - do l = 1, n_y - do k = 1, n_x - - if (target_offset >= lat % offset(map, k, l, m) + offset) then - if (k == n_x .and. l == n_y .and. m == n_z) then - ! This is last lattice cell, so target must be here - lat_offset = lat % offset(map, k, l, m) - offset = offset + lat_offset - next_univ => universes(lat % universes(k, l, m)) - if (lat % is_3d) then - path = trim(path) // "(" // trim(to_str(k-1)) // & - "," // trim(to_str(l-1)) // "," // & - trim(to_str(m-1)) // ")" - else - path = trim(path) // "(" // trim(to_str(k-1)) // & - "," // trim(to_str(l-1)) // ")" - end if - call find_offset(i_cell, next_univ, target_offset, offset, path) - return - else - old_m = m - old_l = l - old_k = k - cycle - end if - else - ! Target is at this lattice position - lat_offset = lat % offset(map, old_k, old_l, old_m) - offset = offset + lat_offset - next_univ => universes(lat % universes(old_k, old_l, old_m)) - if (lat % is_3d) then - path = trim(path) // "(" // trim(to_str(old_k-1)) // & - "," // trim(to_str(old_l-1)) // "," // & - trim(to_str(old_m-1)) // ")" - else - path = trim(path) // "(" // trim(to_str(old_k-1)) // & - "," // trim(to_str(old_l-1)) // ")" - end if - call find_offset(i_cell, next_univ, target_offset, offset, path) - return - end if - - end do - end do - end do - - ! ================================================================== - ! HEXAGONAL LATTICES - type is (HexLattice) - - ! Write to the geometry stack - path = trim(path) // "->l" // to_str(lat%id) - - n_z = lat % n_axial - n_y = 2 * lat % n_rings - 1 - n_x = 2 * lat % n_rings - 1 - old_m = 1 - old_l = 1 - old_k = 1 - - ! Loop over lattice coordinates - do m = 1, n_z - do l = 1, n_y - do k = 1, n_x - - ! This array position is never used - if (k + l < lat % n_rings + 1) then - cycle - ! This array position is never used - else if (k + l > 3*lat % n_rings - 1) then - cycle - end if - - if (target_offset >= lat % offset(map, k, l, m) + offset) then - if (k == lat % n_rings .and. l == n_y .and. m == n_z) then - ! This is last lattice cell, so target must be here - lat_offset = lat % offset(map, k, l, m) - offset = offset + lat_offset - next_univ => universes(lat % universes(k, l, m)) - if (lat % is_3d) then - path = trim(path) // "(" // & - trim(to_str(k - lat % n_rings)) // "," // & - trim(to_str(l - lat % n_rings)) // "," // & - trim(to_str(m - 1)) // ")" - else - path = trim(path) // "(" // & - trim(to_str(k - lat % n_rings)) // "," // & - trim(to_str(l - lat % n_rings)) // ")" - end if - call find_offset(i_cell, next_univ, target_offset, offset, path) - return - else - old_m = m - old_l = l - old_k = k - cycle - end if - else - ! Target is at this lattice position - lat_offset = lat % offset(map, old_k, old_l, old_m) - offset = offset + lat_offset - next_univ => universes(lat % universes(old_k, old_l, old_m)) - if (lat % is_3d) then - path = trim(path) // "(" // & - trim(to_str(old_k - lat % n_rings)) // "," // & - trim(to_str(old_l - lat % n_rings)) // "," // & - trim(to_str(old_m - 1)) // ")" - else - path = trim(path) // "(" // & - trim(to_str(old_k - lat % n_rings)) // "," // & - trim(to_str(old_l - lat % n_rings)) // ")" - end if - call find_offset(i_cell, next_univ, target_offset, offset, path) - return - end if - - end do - end do - end do - - end select - - end if - end if + path_len = distribcell_path_len(i_cell-1, map-1, target_offset, & + root_universe-1) + allocate(path_c(path_len)) + call distribcell_path(i_cell-1, map-1, target_offset, root_universe-1, & + path_c) + do i = 1, min(path_len, MAX_LINE_LEN) + if (path_c(i) == C_NULL_CHAR) exit + path(i:i) = path_c(i) end do end subroutine find_offset diff --git a/src/tracking.F90 b/src/tracking.F90 index 3fe6a065b7..25121d7f44 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -174,7 +174,7 @@ contains p % coord(p % n_coord) % cell = NONE if (any(lattice_translation /= 0)) then ! Particle crosses lattice boundary - p % surface = NONE + p % surface = ERROR_INT call cross_lattice(p, lattice_translation) p % event = EVENT_LATTICE else @@ -205,7 +205,7 @@ contains call score_surface_tally(p, active_meshsurf_tallies) ! Clear surface component - p % surface = NONE + p % surface = ERROR_INT if (run_CE) then call collision(p) @@ -475,7 +475,7 @@ contains ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS ! Remove lower coordinate levels and assignment of surface - p % surface = NONE + p % surface = ERROR_INT p % n_coord = 1 call find_cell(p, found) diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index 8264797a7a..f793d4689e 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -205,7 +205,8 @@ contains elseif (this % domain_type == FILTER_CELL) THEN do level = 1, p % n_coord do i_domain = 1, size(this % domain_id) - if (cells(p % coord(level) % cell) % id == this % domain_id(i_domain)) then + if (cells(p % coord(level) % cell) % id() & + == this % domain_id(i_domain)) then i_material = p % material call check_hit(i_domain, i_material, indices, hits, n_mat) end if diff --git a/src/xml_interface.cpp b/src/xml_interface.cpp index aea1c2c17b..b4fcc524c6 100644 --- a/src/xml_interface.cpp +++ b/src/xml_interface.cpp @@ -2,9 +2,7 @@ #include // for std::transform #include -#include -#include "pugixml/pugixml.hpp" #include "error.h" diff --git a/src/xml_interface.h b/src/xml_interface.h index 52aa6003f4..b6f9e3246d 100644 --- a/src/xml_interface.h +++ b/src/xml_interface.h @@ -2,6 +2,7 @@ #define XML_INTERFACE_H #include +#include #include "pugixml/pugixml.hpp"