diff --git a/src/geometry.F90 b/src/geometry.F90 index 3fd5859f7f..a98ae8b95d 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -605,13 +605,7 @@ contains 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 + if (lat % are_valid_indices([j, k, m])) then lat % offset(map, j, k, m) = offset next_univ => universes(lat % get([j-1, k-1, m-1]) + 1) offset = offset + & @@ -729,13 +723,7 @@ contains 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 + if (lat % are_valid_indices([j, k, m])) then next_univ => universes(lat % get([j-1, k-1, m-1]) + 1) ! Found target - stop since target cannot contain itself @@ -811,13 +799,7 @@ contains 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 + if (lat % are_valid_indices([j, k, m])) then call count_instance(universes(lat % get([j-1, k-1, m-1]) & + 1)) end if @@ -844,6 +826,7 @@ contains integer :: levels ! maximum number of levels for this universe integer :: i ! index over cells + integer :: nx, ny, nz ! lattice shape integer :: j, k, m ! indices in lattice integer :: levels_below ! max levels below this universe type(Cell), pointer :: c ! pointer to current cell @@ -871,39 +854,29 @@ contains select type (lat) type is (RectLattice) - - ! Loop over lattice coordinates - do j = 1, lat % n_cells(1) - do k = 1, lat % n_cells(2) - do m = 1, lat % n_cells(3) - 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 - end do + nx = lat % n_cells(1) + ny = lat % n_cells(2) + nz = lat % n_cells(3) 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 % get([j-1, k-1, m-1]) + 1) - levels_below = max(levels_below, maximum_levels(next_univ)) - end if - end do - end do - end do + nx = 2 * lat % n_rings - 1 + ny = 2 * lat % n_rings - 1 + nz = lat % n_axial end select + ! Loop over lattice coordinates + do j = 1, nx + do k = 1, ny + do m = 1, nz + if (lat % are_valid_indices([j, k, m])) then + 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 + end do + end do + end if end do diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index ddd9e4c2d2..117d3feaf4 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -179,7 +179,6 @@ module geometry_header type, abstract :: Lattice type(C_PTR) :: ptr - real(8), allocatable :: pitch(:) ! Pitch along each axis integer :: outside ! Material to fill area outside integer :: outer ! universe to tile outside the lat logical :: is_3d ! Lattice has cells on z axis diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h index f5a3028551..0b2e030c48 100644 --- a/src/hdf5_interface.h +++ b/src/hdf5_interface.h @@ -75,7 +75,7 @@ write_int(hid_t group_id, char const *name, int32_t buffer) template inline void write_double_1D(hid_t group_id, char const *name, - std::array &buffer) + const std::array &buffer) { hsize_t dims[1]{array_len}; hid_t dataspace = H5Screate_simple(1, dims, NULL); diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a3bc6b8260..26e2ff0a3b 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -931,12 +931,11 @@ contains subroutine read_geometry_xml() - integer :: i, j, k, m, i_x, i_a, input_index + integer :: i, j, k, m integer :: n, n_mats, n_x, n_y, n_z, n_rings, 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 @@ -1295,31 +1294,6 @@ contains 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 & @@ -1405,19 +1379,6 @@ contains 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 @@ -1429,9 +1390,6 @@ 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 = 0, n_z-1 do k = 0, 2*n_rings-2 @@ -4179,8 +4137,6 @@ 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() diff --git a/src/lattice.cpp b/src/lattice.cpp index 9a716f76b2..cda5dff372 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -55,6 +55,8 @@ Lattice::to_hdf5(hid_t lat_group) const if (!name.empty()) { write_string(lat_group, "name", name); } + + to_hdf5_inner(lat_group); } //============================================================================== @@ -265,6 +267,19 @@ RectLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const return local_xyz; } +//============================================================================== + +void +RectLattice::to_hdf5_inner(hid_t lat_group) const +{ + if (is_3d) { + write_double_1D(lat_group, "pitch", pitch); + } else { + std::array out {{pitch[0], pitch[1]}}; + write_double_1D(lat_group, "pitch", out); + } +} + //============================================================================== // HexLattice implementation //============================================================================== @@ -650,6 +665,19 @@ HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const return local_xyz; } +//============================================================================== + +void +HexLattice::to_hdf5_inner(hid_t lat_group) const +{ + if (is_3d) { + write_double_1D(lat_group, "pitch", pitch); + } else { + std::array out {{pitch[0]}}; + write_double_1D(lat_group, "pitch", out); + } +} + //============================================================================== // Non-method functions //============================================================================== diff --git a/src/lattice.h b/src/lattice.h index b21dc56040..78d83a2cbd 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -89,6 +89,8 @@ public: protected: bool is_3d; //! Has divisions along the z-axis + + virtual void to_hdf5_inner(hid_t group_id) const = 0; }; //============================================================================== @@ -115,6 +117,8 @@ public: std::array get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const; + 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 @@ -142,6 +146,8 @@ public: std::array get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const; + 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 diff --git a/src/summary.F90 b/src/summary.F90 index 834b2661a9..160f3bfbda 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -260,7 +260,6 @@ contains call lat % to_hdf5(lattice_group) ! Write name, pitch, and outer universe - call write_dataset(lattice_group, "pitch", lat%pitch) if (lat % outer > 0) then call write_dataset(lattice_group, "outer", universes(lat % outer) % id) else