Clean up lattice code

This commit is contained in:
Sterling Harper 2018-04-29 21:21:39 -04:00
parent 98009a268d
commit f217e2fc74
7 changed files with 58 additions and 97 deletions

View file

@ -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

View file

@ -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

View file

@ -75,7 +75,7 @@ write_int(hid_t group_id, char const *name, int32_t buffer)
template<std::size_t array_len> inline void
write_double_1D(hid_t group_id, char const *name,
std::array<double, array_len> &buffer)
const std::array<double, array_len> &buffer)
{
hsize_t dims[1]{array_len};
hid_t dataspace = H5Screate_simple(1, dims, NULL);

View file

@ -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 <pitch> must be the same as &
&the number of entries on <dimension>.")
end if
else if (node_word_count(node_lat, "pitch") /= n) then
call fatal_error("Number of entries on <pitch> must be the same as &
&the number of entries on <dimension>.")
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 <n_axial> must have <pitch> &
&specified by 2 numbers.")
else if ((.not. lat % is_3d) .and. n /= 1) then
call fatal_error("A hexagonal lattice without <n_axial> must have &
&<pitch> 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()

View file

@ -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<double, 2> 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<double, 1> out {{pitch[0]}};
write_double_1D(lat_group, "pitch", out);
}
}
//==============================================================================
// Non-method functions
//==============================================================================

View file

@ -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<double, 3>
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<int, 3> n_cells; //! Number of cells along each axis
std::array<double, 3> lower_left; //! Global lower-left corner of the lattice
@ -142,6 +146,8 @@ public:
std::array<double, 3>
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

View file

@ -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