Move most lattice summary.h5 writing to C++

This commit is contained in:
Sterling Harper 2018-05-10 00:12:22 -04:00
parent e8303c863d
commit 42e1148fa6
6 changed files with 93 additions and 94 deletions

View file

@ -198,7 +198,6 @@ module geometry_header
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
end type RectLattice
!===============================================================================
@ -208,7 +207,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
end type HexLattice
!===============================================================================

View file

@ -98,8 +98,16 @@ write_int(hid_t group_id, char const *name, int32_t buffer)
H5Dclose(dataset);
}
template<std::size_t array_len> void
write_int(hid_t group_id, char const *name,
const std::array<int, array_len> &buffer, bool indep)
{
hsize_t dims[1] {array_len};
write_dataset(group_id, 1, dims, name, H5T_NATIVE_INT, buffer.data(), indep);
}
template<std::size_t array_len> inline void
template<std::size_t array_len> void
write_double_1D(hid_t group_id, char const *name,
const std::array<double, array_len> &buffer)
{

View file

@ -1284,15 +1284,6 @@ contains
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 <lower_left> must be the same &
&as the number of entries on <dimension>.")
end if
allocate(lat % lower_left(n))
call get_node_array(node_lat, "lower_left", lat % lower_left)
! 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 &
@ -1365,19 +1356,6 @@ contains
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 <n_axial> must have &
&<center> specified by 3 numbers.")
else if ((.not. lat % is_3d) .and. n /= 2) then
call fatal_error("A hexagonal lattice without <n_axial> must have &
&<center> specified by 2 numbers.")
end if
allocate(lat % center(n))
call get_node_array(node_lat, "center", lat % center)
! Copy number of dimensions
n_rings = lat % n_rings
n_z = lat % n_axial

View file

@ -272,11 +272,57 @@ RectLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const
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_1D(lat_group, "pitch", pitch);
write_double_1D(lat_group, "lower_left", lower_left);
write_int(lat_group, "dimension", n_cells, false);
} else {
std::array<double, 2> out {{pitch[0], pitch[1]}};
write_double_1D(lat_group, "pitch", out);
std::array<double, 2> pitch_short {{pitch[0], pitch[1]}};
write_double_1D(lat_group, "pitch", pitch_short);
std::array<double, 2> ll_short {{lower_left[0], lower_left[1]}};
write_double_1D(lat_group, "lower_left", ll_short);
std::array<int, 2> 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<hsize_t>(n_cells[0])};
hsize_t ny {static_cast<hsize_t>(n_cells[1])};
hsize_t nz {static_cast<hsize_t>(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] = universes_c[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<hsize_t>(n_cells[0])};
hsize_t ny {static_cast<hsize_t>(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] = universes_c[universes[indx1]]->id;
}
}
hsize_t dims[3] {1, ny, nx};
write_int(lat_group, 3, dims, "universes", out, false);
}
}
@ -670,12 +716,45 @@ HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const
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_1D(lat_group, "pitch", pitch);
write_double_1D(lat_group, "center", center);
} else {
std::array<double, 1> out {{pitch[0]}};
write_double_1D(lat_group, "pitch", out);
std::array<double, 1> pitch_short {{pitch[0]}};
write_double_1D(lat_group, "pitch", pitch_short);
std::array<double, 2> center_short {{center[0], center[1]}};
write_double_1D(lat_group, "center", center_short);
}
// Write the universe ids.
hsize_t nx {static_cast<hsize_t>(2*n_rings - 1)};
hsize_t ny {static_cast<hsize_t>(2*n_rings - 1)};
hsize_t nz {static_cast<hsize_t>(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] = universes_c[universes[indx]]->id;
}
}
}
}
hsize_t dims[3] {nz, ny, nx};
write_int(lat_group, 3, dims, "universes", out, false);
}
//==============================================================================

View file

@ -41,7 +41,6 @@ class Lattice
public:
int32_t id; //! Universe ID number
std::string name; //! User-defined name
//std::vector<double> pitch; //! Pitch along each basis
std::vector<int32_t> universes; //! Universes filling each lattice tile
int32_t outer{NO_OUTER_UNIVERSE}; //! Universe tiled outside the lattice
//std::vector<int32_t> offset; //! Distribcell offsets

View file

@ -160,8 +160,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
integer, allocatable :: cell_materials(:)
integer, allocatable :: cell_ids(:)
real(8), allocatable :: cell_temperatures(:)
@ -311,68 +310,6 @@ contains
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%get([j-1, lat%n_cells(2)-k-1, m-1])+1)%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%get([j-1,k-1,m-1])+1)%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