From 9da5e36daf5d9a6e489794c625d62283fc27d330 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 21 May 2018 16:32:44 -0400 Subject: [PATCH] Move distribcell offsets to C++ --- src/cell.cpp | 2 + src/cell.h | 2 + src/geometry.F90 | 92 +++--------------------- src/geometry_aux.cpp | 68 +++++++++++------- src/geometry_aux.h | 66 +++++++++++++++++ src/geometry_header.F90 | 49 +++++++++---- src/input_xml.F90 | 44 ++++-------- src/lattice.cpp | 48 +++++++++++++ src/lattice.h | 40 +++++++---- src/summary.F90 | 5 -- src/tallies/tally_filter_distribcell.F90 | 34 ++++----- 11 files changed, 265 insertions(+), 185 deletions(-) create mode 100644 src/geometry_aux.h diff --git a/src/cell.cpp b/src/cell.cpp index 6d7e598eba..1d630b3d4f 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -485,6 +485,8 @@ extern "C" { 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) diff --git a/src/cell.h b/src/cell.h index fff8b749c3..2b03e1cd0a 100644 --- a/src/cell.h +++ b/src/cell.h @@ -71,6 +71,8 @@ public: 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); diff --git a/src/geometry.F90 b/src/geometry.F90 index d053625877..c3b647a822 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -26,7 +26,7 @@ module geometry end function cell_contains_c function count_universe_instances(search_univ, target_univ_id) bind(C) & - result(count) + 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 @@ -174,19 +174,19 @@ contains 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) + 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, & - p % coord(k + 1) % lattice_y, & - p % coord(k + 1) % lattice_z) + offset = offset + lattices(p % coord(k + 1) % lattice) % obj & + % offset(distribcell_index - 1, & + [p % coord(k + 1) % lattice_x - 1, & + p % coord(k + 1) % lattice_y - 1, & + p % coord(k + 1) % lattice_z - 1]) end if end if end do @@ -534,82 +534,6 @@ 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) - - 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 :: i ! index over cells - integer :: j, k, m ! indices in lattice - integer :: offset ! total offset for a given cell - integer :: cell_index ! index in cells array - type(Cell), pointer :: c ! pointer to current cell - class(Lattice), pointer :: lat ! pointer to current lattice - - offset = 0 - - do i = 1, size(univ % cells) - - 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 - c % offset(map) = offset - offset = offset + count_universe_instances(c % fill - 1, univ_id) - - ! ==================================================================== - ! 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) - 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 - offset = offset + count_universe_instances(& - lat % get([j-1, k-1, m-1]), univ_id) - 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 (lat % are_valid_indices([j, k, m])) then - lat % offset(map, j, k, m) = offset - offset = offset + count_universe_instances(& - lat % get([j-1, k-1, m-1]), univ_id) - end if - end do - end do - end do - end select - - end if - end do - - end subroutine calc_offsets - !=============================================================================== ! MAXIMUM_LEVELS determines the maximum number of nested coordinate levels in ! the geometry diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 911e53c281..0685542409 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -1,5 +1,4 @@ -//! \file geometry_aux.cpp -//! Auxilary functions for geometry initialization and general data handling. +#include "geometry_aux.h" #include #include @@ -14,11 +13,9 @@ namespace openmc { -//============================================================================== -//! Replace Universe, Lattice, and Material IDs with indices. //============================================================================== -extern "C" void +void adjust_indices_c() { // Adjust material/fill idices. @@ -65,15 +62,9 @@ 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 +int32_t find_root_universe() { // Find all the universes listed as a cell fill. @@ -115,15 +106,24 @@ find_root_universe() } //============================================================================== -//! 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). + +void +allocate_offset_tables(int n_maps) +{ + for (Cell *c : cells_c) { + if (c->type != FILL_MATERIAL) { + c->offset.resize(n_maps, C_NONE); + } + } + + for (Lattice *lat : lattices_c) { + lat->allocate_offset_table(n_maps); + } +} + //============================================================================== -extern "C" void +void count_cell_instances(int32_t univ_indx) { for (int32_t cell_indx : universes_c[univ_indx]->cells) { @@ -144,14 +144,9 @@ count_cell_instances(int32_t univ_indx) } } -//============================================================================== -//! Recursively search through universes and count the number of instances of -//! the target universe in the geometry tree. -//! @param search_univ The index of the universe to begin searching from. -//! @param target_univ_id The ID of the universe to be counted. //============================================================================== -extern "C" int +int count_universe_instances(int32_t search_univ, int32_t target_univ_id) { // If this is the target, it can't contain itself. @@ -179,4 +174,27 @@ count_universe_instances(int32_t search_univ, int32_t target_univ_id) return count; } +//============================================================================== + +void +fill_offset_tables(int32_t target_univ_id, int map) +{ + for (Universe *univ : universes_c) { + int32_t offset {0}; // TODO: is this a bug? It matches F90 implementation. + for (int32_t cell_indx : univ->cells) { + Cell &c = *cells_c[cell_indx]; + + if (c.type == FILL_UNIVERSE) { + c.offset[map] = offset; + int32_t search_univ = c.fill - 1; // TODO: off-by-one + offset += count_universe_instances(search_univ, target_univ_id); + + } else if (c.type == FILL_LATTICE) { + Lattice &lat = *lattices_c[c.fill - 1]; // TODO: off-by-one + offset = lat.fill_offset_table(offset, target_univ_id, map); + } + } + } +} + } // namespace openmc diff --git a/src/geometry_aux.h b/src/geometry_aux.h new file mode 100644 index 0000000000..0a0500794b --- /dev/null +++ b/src/geometry_aux.h @@ -0,0 +1,66 @@ +//! \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); + +} // namespace openmc +#endif // GEOMETRY_AUX_H diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 012e976afe..fc70a8eeca 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -68,7 +68,7 @@ module geometry_header end subroutine cell_set_universe_c function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') & - result(n_instances) + result(n_instances) import C_PTR, C_INT32_T implicit none type(C_PTR), intent(in), value :: cell_ptr @@ -94,9 +94,16 @@ module geometry_header 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 - implicit none type(C_PTR), intent(in), value :: cell_ptr integer(HID_T), intent(in), value :: group end subroutine cell_to_hdf5_c @@ -104,14 +111,12 @@ module geometry_header function lattice_pointer_c(lat_ind) bind(C, name='lattice_pointer') & result(ptr) import C_PTR, C_INT32_T - implicit none 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 - implicit none type(C_PTR), intent(in), value :: lat_ptr integer(C_INT32_T) :: id end function lattice_id_c @@ -119,7 +124,6 @@ module geometry_header 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 - implicit none type(C_PTR), intent(in), value :: lat_ptr integer(C_INT), intent(in) :: i_xyz(3) logical(C_BOOL) :: is_valid @@ -128,7 +132,6 @@ module geometry_header 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 - implicit none type(C_PTR), intent(in), value :: lat_ptr real(C_DOUBLE), intent(in) :: xyz(3) real(C_DOUBLE), intent(in) :: uvw(3) @@ -140,7 +143,6 @@ module geometry_header subroutine lattice_get_indices_c(lat_ptr, xyz, i_xyz) & bind(C, name='lattice_get_indices') import C_PTR, C_INT, C_DOUBLE - implicit none type(C_PTR), intent(in), value :: lat_ptr real(C_DOUBLE), intent(in) :: xyz(3) integer(C_INT), intent(out) :: i_xyz(3) @@ -149,7 +151,6 @@ module geometry_header 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 - implicit none type(C_PTR), intent(in), value :: lat_ptr real(C_DOUBLE), intent(in) :: global_xyz(3) integer(C_INT), intent(in) :: i_xyz(3) @@ -158,11 +159,19 @@ module geometry_header subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5') import HID_T, C_PTR - implicit none 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 @@ -173,7 +182,6 @@ module geometry_header function lattice_universe_c(lat_ptr, i_xyz) & bind(C, name='lattice_universe') result(univ) import C_PTR, C_INT32_t, C_INT - implicit none type(C_PTR), intent(in), value :: lat_ptr integer(C_INT), intent(in) :: i_xyz(3) integer(C_INT32_T) :: univ @@ -181,7 +189,6 @@ module geometry_header subroutine extend_cells_c(n) bind(C) import C_INT32_t - implicit none integer(C_INT32_T), intent(in), value :: n end subroutine extend_cells_c end interface @@ -207,7 +214,6 @@ module geometry_header type(C_PTR) :: ptr logical :: is_3d ! Lattice has cells on z axis - integer, allocatable :: offset(:,:,:,:) ! Distribcell offsets contains procedure :: id => lattice_id @@ -216,6 +222,7 @@ module geometry_header 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 @@ -256,8 +263,6 @@ module geometry_header 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 :: distribcell_index ! Index corresponding to this cell in @@ -282,6 +287,7 @@ module geometry_header 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 @@ -350,6 +356,14 @@ contains 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 @@ -422,6 +436,13 @@ contains call cell_distance_c(this % ptr, xyz, uvw, on_surface, min_dist, i_surf) end subroutine cell_distance + 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 + subroutine cell_to_hdf5(this, group) class(Cell), intent(in) :: this integer(HID_T), intent(in) :: group diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 44a161f470..c68ed4a8c0 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -11,7 +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, neighbor_lists + use geometry, only: maximum_levels, neighbor_lists use geometry_header use hdf5_interface use list_header, only: ListChar, ListInt, ListReal @@ -50,10 +50,21 @@ module input_xml 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 + end subroutine count_cell_instances subroutine read_surfaces(node_ptr) bind(C) import C_PTR @@ -3884,12 +3895,11 @@ contains ! Allocate offset maps at each level in the geometry call allocate_offsets(univ_list) + call allocate_offset_tables(n_maps) ! 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)) - end do + call fill_offset_tables(univ_list(i), i-1) end do end subroutine prepare_distribcell @@ -3952,30 +3962,6 @@ contains 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 - 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() diff --git a/src/lattice.cpp b/src/lattice.cpp index 024e4cb483..7bb5fb6bb7 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -7,6 +7,7 @@ #include "cell.h" #include "constants.h" #include "error.h" +#include "geometry_aux.h" #include "hdf5_interface.h" #include "xml_interface.h" @@ -98,6 +99,26 @@ Lattice::adjust_indices() //============================================================================== +void +Lattice::allocate_offset_table(int n_maps) +{ + offsets.resize(n_maps * universes.size(), C_NONE); +} + +//============================================================================== + +int32_t +Lattice::fill_offset_table(int32_t offset, int32_t target_univ_id, int map) +{ + for (auto 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 { @@ -316,6 +337,17 @@ RectLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const //============================================================================== +int32_t& +RectLattice::offset(int map, const int i_xyz[3]) +{ + int nx = n_cells[0]; + int ny = n_cells[1]; + int nz = n_cells[2]; + return offsets[nx*ny*nz*map + nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]]; +} + +//============================================================================== + void RectLattice::to_hdf5_inner(hid_t lat_group) const { @@ -763,6 +795,17 @@ HexLattice::is_valid_index(int indx) const //============================================================================== +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]]; +} + +//============================================================================== + void HexLattice::to_hdf5_inner(hid_t lat_group) const { @@ -874,6 +917,11 @@ extern "C" { 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);} diff --git a/src/lattice.h b/src/lattice.h index 938294ae4d..9b55025b58 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -41,11 +41,11 @@ class LatticeIter; 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 offset; //! Distribcell offsets + 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); @@ -59,6 +59,11 @@ public: //! Convert internal universe values from IDs to indices using universe_dict. void adjust_indices(); + //! Allocate offset table for distribcell. + void allocate_offset_table(int n_maps); + + int32_t fill_offset_table(int32_t offset, int32_t target_univ_id, int map); + //! 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 @@ -96,6 +101,13 @@ public: return (indx > 0) && (indx < universes.size()); } + //! 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; + //! Write all information needed to reconstruct the lattice to an HDF5 group. //! @param group_id An HDF5 group id. void to_hdf5(hid_t group_id) const; @@ -166,12 +178,14 @@ public: 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]); + void to_hdf5_inner(hid_t group_id) const; protected: - 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 + 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 }; //============================================================================== @@ -200,13 +214,15 @@ public: bool is_valid_index(int indx) const; + int32_t& offset(int map, const int i_xyz[3]); + void to_hdf5_inner(hid_t group_id) const; protected: - 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 + 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 diff --git a/src/summary.F90 b/src/summary.F90 index 964ec34a35..1d33491230 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -226,11 +226,6 @@ 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 if (allocated(c%translation)) then call write_dataset(cell_group, "translation", c%translation) diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index d804fe0ea0..ce2549f58a 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -59,19 +59,19 @@ contains 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) + 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, & - p % coord(i + 1) % lattice_y, & - p % coord(i + 1) % lattice_z) + offset = offset + lattices(p % coord(i + 1) % lattice) % obj & + % offset(distribcell_index - 1, & + [p % coord(i + 1) % lattice_x - 1, & + p % coord(i + 1) % lattice_y - 1, & + p % coord(i + 1) % lattice_z - 1]) end if end if if (this % cell == p % coord(i) % cell) then @@ -211,12 +211,12 @@ contains ! Two cases, lattice or fill cell if (c % type() == FILL_UNIVERSE) then - temp_offset = c % offset(map) + temp_offset = c % offset(map-1) ! Get the offset of the first lattice location else lat => lattices(c % fill) % obj - temp_offset = lat % offset(map, 1, 1, 1) + temp_offset = lat % offset(map-1, [0, 0, 0]) end if ! If the final offset is in the range of offset - temp_offset+offset @@ -248,7 +248,7 @@ contains if (c % type() == FILL_UNIVERSE) then ! Enter this cell to update the current offset - offset = c % offset(map) + offset + offset = c % offset(map-1) + offset next_univ => universes(c % fill) call find_offset(i_cell, next_univ, target_offset, offset, path) @@ -282,10 +282,11 @@ contains do l = 1, n_y do k = 1, n_x - if (target_offset >= lat % offset(map, k, l, m) + offset) then + if (target_offset >= lat % offset(map-1, [k-1, l-1, m-1]) & + + 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) + lat_offset = lat % offset(map-1, [k-1, l-1, m-1]) offset = offset + lat_offset next_univ => universes(lat % get([k-1, l-1, m-1])+1) if (lat % is_3d) then @@ -306,7 +307,7 @@ contains end if else ! Target is at this lattice position - lat_offset = lat % offset(map, old_k, old_l, old_m) + lat_offset = lat % offset(map-1, [old_k-1, old_l-1, old_m-1]) offset = offset + lat_offset next_univ => universes(lat % get([old_k-1, old_l-1, old_m-1])+1) if (lat % is_3d) then @@ -352,10 +353,11 @@ contains cycle end if - if (target_offset >= lat % offset(map, k, l, m) + offset) then + if (target_offset >= lat % offset(map-1, [k-1, l-1, m-1]) & + + 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) + lat_offset = lat % offset(map-1, [k-1, l-1, m-1]) offset = offset + lat_offset next_univ => universes(lat % get([k-1, l-1, m-1])+1) if (lat % is_3d) then @@ -378,7 +380,7 @@ contains end if else ! Target is at this lattice position - lat_offset = lat % offset(map, old_k, old_l, old_m) + lat_offset = lat % offset(map-1, [old_k-1, old_l-1, old_m-1]) offset = offset + lat_offset next_univ => universes(lat % get([old_k-1, old_l-1, old_m-1])+1) if (lat % is_3d) then