diff --git a/src/cell.cpp b/src/cell.cpp index ed16a566a3..56d7c45ca2 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -8,6 +8,7 @@ #include "constants.h" #include "error.h" #include "hdf5_interface.h" +#include "lattice.h" #include "surface.h" #include "xml_interface.h" @@ -29,6 +30,19 @@ constexpr int32_t OP_UNION {std::numeric_limits::max() - 4}; extern "C" double FP_PRECISION; +//============================================================================== +// Global variables +//============================================================================== + +// Braces force n_cells to be defined here, not just declared. +extern "C" {int32_t n_cells {0};} + +std::vector cells_c; +std::map cell_dict; + +std::vector universes_c; +std::map universe_dict; + //============================================================================== //! Convert region specification string to integer tokens. //! @@ -444,6 +458,11 @@ adjust_indices_c() fatal_error(err_msg); } } + + // Change all lattice universe values from IDs to indices. + for (Lattice *l : lattices_c) { + l->adjust_indices(); + } } //============================================================================== diff --git a/src/cell.h b/src/cell.h index 99ca2caf4a..5c06d81200 100644 --- a/src/cell.h +++ b/src/cell.h @@ -16,16 +16,15 @@ namespace openmc { // Global variables //============================================================================== -// Braces force n_cells to be defined here, not just declared. -extern "C" {int32_t n_cells {0};} +extern "C" int32_t n_cells; class Cell; -std::vector cells_c; -std::map cell_dict; +extern std::vector cells_c; +extern std::map cell_dict; class Universe; -std::vector universes_c; -std::map universe_dict; +extern std::vector universes_c; +extern std::map universe_dict; //============================================================================== //! A geometry primitive that fills all space and contains cells. diff --git a/src/constants.h b/src/constants.h index b725809fe0..8cc837a735 100644 --- a/src/constants.h +++ b/src/constants.h @@ -5,7 +5,8 @@ namespace openmc{ extern "C" double FP_COINCIDENT; -constexpr double INFTY{std::numeric_limits::max()}; +extern "C" double FP_PRECISION; +constexpr double INFTY {std::numeric_limits::max()}; constexpr int C_NONE {-1}; } // namespace openmc diff --git a/src/geometry.F90 b/src/geometry.F90 index fff18d6c85..3fd5859f7f 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -254,7 +254,7 @@ 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)-1, i_xyz(2)-1, i_xyz(3)-1]) + 1 else ! Particle is outside the lattice. @@ -331,7 +331,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)-1, i_xyz(2)-1, i_xyz(3)-1]) + 1 call find_cell(p, found) if (.not. found) then @@ -591,7 +592,7 @@ contains do k = 1, lat % n_cells(2) do m = 1, lat % n_cells(3) lat % offset(map, j, k, m) = offset - next_univ => universes(lat % universes(j, k, m)) + next_univ => universes(lat % get([j-1, k-1, m-1]) + 1) offset = offset + & count_target(next_univ, counts, found, univ_id, map) end do @@ -612,7 +613,7 @@ contains cycle else lat % offset(map, j, k, m) = offset - next_univ => universes(lat % universes(j, k, m)) + next_univ => universes(lat % get([j-1, k-1, m-1]) + 1) offset = offset + & count_target(next_univ, counts, found, univ_id, map) end if @@ -707,7 +708,7 @@ contains do j = 1, lat % n_cells(1) do k = 1, lat % n_cells(2) do m = 1, lat % n_cells(3) - next_univ => universes(lat % universes(j, k, m)) + next_univ => universes(lat % get([j-1, k-1, m-1]) + 1) ! Found target - stop since target cannot contain itself if (next_univ % id == univ_id) then @@ -735,7 +736,7 @@ contains else if (j + k > 3*lat % n_rings - 1) then cycle else - next_univ => universes(lat % universes(j, k, m)) + next_univ => universes(lat % get([j-1, k-1, m-1]) + 1) ! Found target - stop since target cannot contain itself if (next_univ % id == univ_id) then @@ -799,7 +800,7 @@ contains do j = 1, lat % n_cells(1) do k = 1, lat % n_cells(2) do m = 1, lat % n_cells(3) - call count_instance(universes(lat % universes(j, k, m))) + call count_instance(universes(lat % get([j-1, k-1, m-1])+1)) end do end do end do @@ -817,7 +818,8 @@ contains else if (j + k > 3*lat % n_rings - 1) then cycle else - call count_instance(universes(lat % universes(j, k, m))) + call count_instance(universes(lat % get([j-1, k-1, m-1]) & + + 1)) end if end do end do @@ -874,7 +876,7 @@ contains do j = 1, lat % n_cells(1) do k = 1, lat % n_cells(2) do m = 1, lat % n_cells(3) - next_univ => universes(lat % universes(j, k, m)) + next_univ => universes(lat % get([j-1, k-1, m-1]) + 1) levels_below = max(levels_below, maximum_levels(next_univ)) end do end do @@ -893,7 +895,7 @@ contains else if (j + k > 3*lat % n_rings - 1) then cycle else - next_univ => universes(lat % universes(j, k, m)) + next_univ => universes(lat % get([j-1, k-1, m-1]) + 1) levels_below = max(levels_below, maximum_levels(next_univ)) end if end do diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index af63cad7f4..ddd9e4c2d2 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -148,6 +148,15 @@ module geometry_header implicit none integer(C_INT32_T), intent(in), value :: n end subroutine extend_cells_c + + function lattice_universe_c(lat_ptr, i_xyz) & + bind(C, name='lattice_universe') result(univ) + use ISO_C_BINDING + implicit none + 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 end interface !=============================================================================== @@ -170,9 +179,7 @@ module geometry_header type, abstract :: Lattice type(C_PTR) :: ptr - !integer :: id ! Universe number for lattice 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 @@ -182,6 +189,7 @@ module geometry_header 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 :: to_hdf5 => lattice_to_hdf5 @@ -297,6 +305,13 @@ contains 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) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 6457f587cc..a3bc6b8260 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1331,7 +1331,6 @@ contains 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") @@ -1340,21 +1339,16 @@ contains &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 m = 0, n_z-1 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)) + do j = 0, n_x - 1 + univ_id = lat % get([j, k, m]) + if (find(fill_univ_ids, univ_id) == -1) & + call fill_univ_ids % push_back(univ_id) end do end do end do - deallocate(temp_int_array) ! Read outer universe for area outside lattice. lat % outer = NO_OUTER_UNIVERSE @@ -1427,7 +1421,6 @@ contains ! 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") @@ -1440,91 +1433,15 @@ contains 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 + do m = 0, n_z-1 + do k = 0, 2*n_rings-2 + do j = 0, 2*n_rings-2 + univ_id = lat % get([j, k, m]) + if (find(fill_univ_ids, univ_id) == -1) & + call fill_univ_ids % push_back(univ_id) 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 @@ -4312,46 +4229,6 @@ contains 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 diff --git a/src/lattice.cpp b/src/lattice.cpp index 2902e3e475..9a716f76b2 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -4,6 +4,8 @@ #include #include +#include "cell.h" +#include "constants.h" #include "error.h" #include "hdf5_interface.h" #include "xml_interface.h" @@ -131,6 +133,36 @@ RectLattice::RectLattice(pugi::xml_node lat_node) //============================================================================== +int32_t& +RectLattice::operator[](const int i_xyz[3]) +{ + int nx = n_cells[0]; + int ny = n_cells[1]; + int nz = n_cells[2]; + int indx = nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]; + return universes[indx]; +} + +//============================================================================== + +void +RectLattice::adjust_indices() +{ + int nx = n_cells[0]; + int ny = n_cells[1]; + int nz = n_cells[2]; + for (int iz = 0; iz < nz; iz++) { + for (int iy = 0; iy < ny; iy++) { + for (int ix = 0; ix < nx; ix++) { + int indx = nx*ny*iz + nx*iy + ix; + universes[indx] = universe_dict[universes[indx]]; + } + } + } +} + +//============================================================================== + bool RectLattice::are_valid_indices(const int i_xyz[3]) const { @@ -386,6 +418,37 @@ HexLattice::HexLattice(pugi::xml_node lat_node) //============================================================================== +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]; +} + +//============================================================================== + +void +HexLattice::adjust_indices() +{ + for (int iz = 0; iz < n_axial; iz++) { + for (int ia = 0; ia < 2*n_rings-1; ia++) { + for (int ix = 0; ix < 2*n_rings-1; ix++) { + int i_xyz[3] {ix+1, ia+1, iz+1}; + if (are_valid_indices(i_xyz)) { + int indx = (2*n_rings-1)*(2*n_rings-1) * iz + + (2*n_rings-1) * ia + + ix; + universes[indx] = universe_dict[universes[indx]]; + } + } + } + } +} + +//============================================================================== + bool HexLattice::are_valid_indices(const int i_xyz[3]) const { @@ -655,6 +718,9 @@ extern "C" { } 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 index 9ade739c02..b21dc56040 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -15,15 +15,7 @@ namespace openmc { //============================================================================== -// Constants that should eventually be moved out of this file -//============================================================================== - -extern "C" double FP_PRECISION; -constexpr double INFTY{std::numeric_limits::max()}; -constexpr int C_NONE {-1}; - -//============================================================================== -// Constants +// Module constants //============================================================================== constexpr int32_t NO_OUTER_UNIVERSE{-1}; @@ -58,6 +50,11 @@ public: virtual ~Lattice() {} + virtual int32_t& operator[](const int i_xyz[3]) = 0; + + //! Convert internal universe values from IDs to indices using universe_dict. + virtual void adjust_indices() = 0; + //! 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 @@ -104,6 +101,10 @@ public: virtual ~RectLattice() {} + int32_t& operator[](const int i_xyz[3]); + + void adjust_indices(); + bool are_valid_indices(const int i_xyz[3]) const; std::pair> @@ -127,6 +128,10 @@ public: virtual ~HexLattice() {} + int32_t& operator[](const int i_xyz[3]); + + void adjust_indices(); + bool are_valid_indices(const int i_xyz[3]) const; std::pair> diff --git a/src/summary.F90 b/src/summary.F90 index e2ef1a2c14..834b2661a9 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -288,7 +288,7 @@ contains 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 + universes(lat%get([j-1, lat%n_cells(2)-k-1, m-1])+1)%id end do end do end do @@ -319,7 +319,7 @@ contains lattice_universes(j,k,m) = -1 cycle end if - lattice_universes(j,k,m) = universes(lat%universes(j,k,m))%id + lattice_universes(j,k,m) = universes(lat%get([j-1,k-1,m-1])+1)%id end do end do end do diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index e7a9645ae5..39f8db5af4 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -289,7 +289,7 @@ contains ! 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)) + next_univ => universes(lat % get([k-1, l-1, m-1])+1) if (lat % is_3d) then path = trim(path) // "(" // trim(to_str(k-1)) // & "," // trim(to_str(l-1)) // "," // & @@ -310,7 +310,7 @@ contains ! 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)) + next_univ => universes(lat % get([old_k-1, old_l-1, old_m-1])+1) if (lat % is_3d) then path = trim(path) // "(" // trim(to_str(old_k-1)) // & "," // trim(to_str(old_l-1)) // "," // & @@ -359,7 +359,7 @@ contains ! 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)) + next_univ => universes(lat % get([k-1, l-1, m-1])+1) if (lat % is_3d) then path = trim(path) // "(" // & trim(to_str(k - lat % n_rings)) // "," // & @@ -382,7 +382,7 @@ contains ! 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)) + next_univ => universes(lat % get([old_k-1, old_l-1, old_m-1])+1) if (lat % is_3d) then path = trim(path) // "(" // & trim(to_str(old_k - lat % n_rings)) // "," // &