From 35511121075047d73c7cf588531b03e2f7f6ffa5 Mon Sep 17 00:00:00 2001 From: Will Boyd Date: Tue, 25 Nov 2014 14:45:51 -0500 Subject: [PATCH] Began refactorization of summary/statepoint to conform to API --- src/ace.F90 | 1 - src/geometry.F90 | 6 +- src/geometry_header.F90 | 3 +- src/hdf5_interface.F90 | 74 +++++++++++++++++++++++++ src/hdf5_summary.F90 | 65 +++++++++++++++++++++- src/input_xml.F90 | 5 +- src/mpiio_interface.F90 | 2 + src/output_interface.F90 | 85 +++++++++++++++++++++++++++++ src/tracking.F90 | 2 +- tests/test_infinite_cell/results.py | 7 +-- 10 files changed, 235 insertions(+), 15 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index 584254149d..4d13bfc1df 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -190,7 +190,6 @@ contains ! Deallocate temporary arrays for names of nuclides and S(a,b) tables if (allocated(mat % names)) deallocate(mat % names) - if (allocated(mat % sab_names)) deallocate(mat % sab_names) end do MATERIAL_LOOP2 diff --git a/src/geometry.F90 b/src/geometry.F90 index 2d9f3355a8..3868412277 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -209,9 +209,9 @@ contains end if ! Apply rotation - if (allocated(c % rotation)) then - p % coord % xyz = matmul(c % rotation, p % coord % xyz) - p % coord % uvw = matmul(c % rotation, p % coord % uvw) + if (allocated(c % rotation_matrix)) then + p % coord % xyz = matmul(c % rotation_matrix, p % coord % xyz) + p % coord % uvw = matmul(c % rotation_matrix, p % coord % uvw) p % coord % rotated = .true. end if diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 93863865c4..0cfb357b81 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -65,8 +65,9 @@ module geometry_header ! here too ! Rotation matrix and translation vector - real(8), allocatable :: rotation(:,:) real(8), allocatable :: translation(:) + real(8), allocatable :: rotation(:) + real(8), allocatable :: rotation_matrix(:,:) end type Cell ! array index of universe 0 diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index c61880780a..125fdaa0f0 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -38,6 +38,7 @@ module hdf5_interface module procedure hdf5_write_integer_4Darray module procedure hdf5_write_long module procedure hdf5_write_string + module procedure hdf5_write_string_1Darray #ifdef MPI module procedure hdf5_write_double_parallel module procedure hdf5_write_double_1Darray_parallel @@ -51,6 +52,7 @@ module hdf5_interface module procedure hdf5_write_integer_4Darray_parallel module procedure hdf5_write_long_parallel module procedure hdf5_write_string_parallel +! module procedure hdf5_write_string_1Darray_parallel #endif end interface hdf5_write_data @@ -68,6 +70,7 @@ module hdf5_interface module procedure hdf5_read_integer_4Darray module procedure hdf5_read_long module procedure hdf5_read_string + module procedure hdf5_read_string_1Darray #ifdef MPI module procedure hdf5_read_double_parallel module procedure hdf5_read_double_1Darray_parallel @@ -81,6 +84,7 @@ module hdf5_interface module procedure hdf5_read_integer_4Darray_parallel module procedure hdf5_read_long_parallel module procedure hdf5_read_string_parallel +! module procedure hdf5_read_string_1Darray_parallel #endif end interface hdf5_read_data @@ -786,6 +790,76 @@ contains end subroutine hdf5_read_string +!=============================================================================== +! HDF5_WRITE_STRING_1DARRAY writes string data +!=============================================================================== + + subroutine hdf5_write_string_1Darray(group, name, buffer, length, rank) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer, intent(in) :: length ! length of strings + integer, intent(in) :: rank ! number of strings + character(*), intent(in) :: buffer(length,rank) ! data to write + + integer :: dims ! number of strings + + ! Insert null character at end of string when writing + call h5tset_strpad_f(H5T_STRING, H5T_STR_NULLPAD_F, hdf5_err) + + ! Create the dataspace and dataset + call h5screate_simple_f(1, dims1, dspace, hdf5_err) + call h5dcreate_f(group, name, H5T_STRING, dspace, dset, hdf5_err) + + ! Set up dimesnions of string to write + dims2 = (/length, rank/) ! full array of strings to write + dims = 1 ! length of string + + ! Write the variable dataset + !call h5ltmake_dataset_string_f(dset, H5T_STRING, buffer, dims2, dims1, hdf5_err, & + ! mem_space_id=dspace) + + !call h5ltmake_dataset_f(group, name, dims, dims2, H5T_STRING, buffer, hdf5_err) + + ! Close all + call h5dclose_f(dset, hdf5_err) + call h5sclose_f(dspace, hdf5_err) + + end subroutine hdf5_write_string_1Darray + +!=============================================================================== +! HDF5_READ_STRING_1DARRAY reads string data +!=============================================================================== + + subroutine hdf5_read_string_1Darray(group, name, buffer, length, rank) + + integer(HID_T), intent(in) :: group ! name of group + character(*), intent(in) :: name ! name of data + integer, intent(in) :: length ! length of strings + integer, intent(in) :: rank ! number of strings + character(*), intent(inout) :: buffer(length,rank) ! read data to here + + ! Open dataset + call h5dopen_f(group, name, dset, hdf5_err) + + ! Get dataspace to read + call h5dget_space_f(dset, dspace, hdf5_err) + + ! Set dimensions + dims2 = (/length, rank/) + dims1(1) = length + + ! Read in the data + !call h5ltread_dataset_string_f(dset, H5T_STRING, buffer, dims2, dims1, hdf5_err, & + ! mem_space_id=dspace, xfer_prp = plist) + + !call h5ltread_dataset_f(group, name, dims, dims2, H5T_STRING, buffer, hdf5_err) + + ! Close dataset + call h5dclose_f(dset, hdf5_err) + + end subroutine hdf5_read_string_1Darray + !=============================================================================== ! HDF5_WRITE_ATTRIBUTE_STRING writes a string attribute to a variables !=============================================================================== diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90 index c9e367200c..09761c87a1 100644 --- a/src/hdf5_summary.F90 +++ b/src/hdf5_summary.F90 @@ -124,6 +124,10 @@ contains CELL_LOOP: do i = 1, n_cells c => cells(i) + ! Write internal OpenMC index for this cell + call su % write_data(i, "index", & + group="geometry/cells/cell " // trim(to_str(c % id))) + ! Write universe for this cell call su % write_data(universes(c % universe) % id, "universe", & group="geometry/cells/cell " // trim(to_str(c % id))) @@ -140,11 +144,33 @@ contains call su % write_data(materials(c % material) % id, "material", & group="geometry/cells/cell " // trim(to_str(c % id))) end if + case (CELL_FILL) call su % write_data("universe", "fill_type", & group="geometry/cells/cell " // trim(to_str(c % id))) - call su % write_data(universes(c % fill) % id, "material", & + call su % write_data(universes(c % fill) % id, "universe", & group="geometry/cells/cell " // trim(to_str(c % id))) + + if (allocated(c % translation)) then + call su % write_data(1, "translated", & + group="geometry/cells/cell " // trim(to_str(c % id))) + call su % write_data(c % translation, "translation", length=3, & + group="geometry/cells/cell " // trim(to_str(c % id))) + else + call su % write_data(0, "translated", & + group="geometry/cells/cell " // trim(to_str(c % id))) + end if + + if (allocated(c % rotation_matrix)) then + call su % write_data(1, "rotated", & + group="geometry/cells/cell " // trim(to_str(c % id))) + call su % write_data(c % rotation, "rotation", length=3, & + group="geometry/cells/cell " // trim(to_str(c % id))) + else + call su % write_data(0, "rotated", & + group="geometry/cells/cell " // trim(to_str(c % id))) + end if + case (CELL_LATTICE) call su % write_data("lattice", "fill_type", & group="geometry/cells/cell " // trim(to_str(c % id))) @@ -171,6 +197,10 @@ contains SURFACE_LOOP: do i = 1, n_surfaces s => surfaces(i) + ! Write internal OpenMC index for this surface + call su % write_data(i, "index", & + group="geometry/surfaces/surface " // trim(to_str(s % id))) + ! Write surface type select case (s % type) case (SURF_PX) @@ -255,6 +285,10 @@ contains UNIVERSE_LOOP: do i = 1, n_universes u => universes(i) + ! Write internal OpenMC index for this universe + call su % write_data(i, "index", & + group="geometry/universes/universe " // trim(to_str(u % id))) + ! Write list of cells in this universe if (u % n_cells > 0) then call su % write_data(u % cells, "cells", length=u % n_cells, & @@ -274,6 +308,10 @@ contains LATTICE_LOOP: do i = 1, n_lattices lat => lattices(i) + ! Write internal OpenMC index for this lattice + call su % write_data(i, "index", & + group="geometry/lattices/lattice " // trim(to_str(lat % id))) + ! Write lattice type select case(lat % type) case (LATTICE_RECT) @@ -295,6 +333,9 @@ contains length=lat % n_dimension, & group="geometry/lattices/lattice " // trim(to_str(lat % id))) + call su % write_data(lat % outside, "outside", & + group="geometry/lattices/lattice " // trim(to_str(lat % id))) + ! Determine dimensions of lattice n_x = lat % dimension(1) n_y = lat % dimension(2) @@ -303,7 +344,7 @@ contains else n_z = 1 end if - + ! Write lattice universes allocate(lattice_universes(n_x, n_y, n_z)) do j = 1, n_x @@ -340,6 +381,10 @@ contains do i = 1, n_materials m => materials(i) + ! Write internal OpenMC index for this material + call su % write_data(i, "index", & + group="materials/material " // trim(to_str(m % id))) + ! Write atom density with units call su % write_data(m % density, "atom_density", & group="materials/material " // trim(to_str(m % id))) @@ -365,13 +410,21 @@ contains group="materials/material " // trim(to_str(m % id))) ! Write S(a,b) information if present + call su % write_data(m % n_sab, "n_sab", & + group="materials/material " // trim(to_str(m % id))) + if (m % n_sab > 0) then call su % write_data(m % i_sab_nuclides, "i_sab_nuclides", & length=m % n_sab, & group="materials/material " // trim(to_str(m % id))) call su % write_data(m % i_sab_tables, "i_sab_tables", & length=m % n_sab, & - group="materials/material " // trim(to_str(m % id))) + group="materials/material " // trim(to_str(m % id))) + do j = 1, m % n_sab + call su % write_data(m % sab_names(j), to_str(j), & + group="materials/material " // & + trim(to_str(m % id)) // "/sab_tables") + end do end if end do @@ -550,6 +603,10 @@ contains NUCLIDE_LOOP: do i = 1, n_nuclides_total nuc => nuclides(i) + ! Write internal OpenMC index for this nuclide + call su % write_data(i, "index", & + group="nuclides/" // trim(nuc % name)) + ! Determine size of cross-sections size_xs = (5 + nuc % n_reaction) * nuc % n_grid * 8 size_total = size_xs @@ -557,6 +614,8 @@ contains ! Write some basic attributes call su % write_data(nuc % zaid, "zaid", & group="nuclides/" // trim(nuc % name)) + call su % write_data(xs_listings(nuc % listing) % alias, "alias", & + group="nuclides/" // trim(nuc % name)) call su % write_data(nuc % awr, "awr", & group="nuclides/" // trim(nuc % name)) call su % write_data(nuc % kT, "kT", & diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 69677d1757..b76ec39277 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1042,13 +1042,14 @@ contains ! Copy rotation angles in x,y,z directions call get_node_array(node_cell, "rotation", temp_int_array3) + c % rotation = temp_int_array3 phi = -temp_int_array3(1) * PI/180.0_8 theta = -temp_int_array3(2) * PI/180.0_8 psi = -temp_int_array3(3) * PI/180.0_8 ! Calculate rotation matrix based on angles given - allocate(c % rotation(3,3)) - c % rotation = reshape((/ & + allocate(c % rotation_matrix(3,3)) + c % rotation_matrix = reshape((/ & cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), & -cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), & cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), & diff --git a/src/mpiio_interface.F90 b/src/mpiio_interface.F90 index 30093adec9..97ca73fb3c 100644 --- a/src/mpiio_interface.F90 +++ b/src/mpiio_interface.F90 @@ -21,6 +21,7 @@ module mpiio_interface module procedure mpi_write_integer_4Darray module procedure mpi_write_long module procedure mpi_write_string + !module procedure mpi_write_string_1Darray end interface mpi_write_data ! Generic HDF5 read procedure interface @@ -37,6 +38,7 @@ module mpiio_interface module procedure mpi_read_integer_4Darray module procedure mpi_read_long module procedure mpi_read_string + !module procedure mpi_read_string_1Darray end interface mpi_read_data contains diff --git a/src/output_interface.F90 b/src/output_interface.F90 index 799fc4d21a..70f46a8f25 100644 --- a/src/output_interface.F90 +++ b/src/output_interface.F90 @@ -1761,6 +1761,91 @@ contains #endif end subroutine read_string +!=============================================================================== +! WRITE_STRING_1DARRAY writes 1-D string data +!=============================================================================== + + subroutine write_string_1Darray(self, buffer, name, group, length, collect) + + integer, intent(in) :: length ! length of array to write + character(*), intent(in) :: buffer(:) ! data to write + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + logical, intent(in), optional :: collect ! collective I/O + class(BinaryOutput) :: self + + character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name + character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name + integer :: n + logical :: collect_ + + ! Get string length + n = len_trim(buffer(1)) + + ! Set name + name_ = trim(name) + + ! Set group + if (present(group)) then + group_ = trim(group) + end if + + ! Set up collective vs. independent I/O + if (present(collect)) then + collect_ = collect + else + collect_ = .true. + end if + +#ifdef HDF5 + ! Check if HDF5 group should be created/opened + if (present(group)) then + call hdf5_open_group(self % hdf5_fh, group_, self % hdf5_grp) + else + self % hdf5_grp = self % hdf5_fh + endif + if (present(group)) call hdf5_close_group(self % hdf5_grp) +#endif + + end subroutine write_string_1Darray + +!=============================================================================== +! READ_STRING_1DARRAY reads 1-D string data +!=============================================================================== + + subroutine read_string_1Darray(self, buffer, name, group, length, collect) + + integer, intent(in) :: length ! length of array to write + character(*), intent(inout) :: buffer(:) ! data to write + character(*), intent(in) :: name ! name of data + character(*), intent(in), optional :: group ! HDF5 group name + logical, intent(in), optional :: collect ! collective I/O + class(BinaryOutput) :: self + + character(len=MAX_WORD_LEN) :: name_ ! HDF5 dataset name + character(len=MAX_WORD_LEN) :: group_ ! HDF5 group name + integer :: n + logical :: collect_ + + ! Get string length + n = len(buffer(1)) + + ! Set name + name_ = trim(name) + + ! Set group + if (present(group)) then + group_ = trim(group) + end if + + ! Set up collective vs. independent I/O + if (present(collect)) then + collect_ = collect + else + collect_ = .true. + end if + + end subroutine read_string_1Darray !=============================================================================== ! WRITE_ATTRIBUTE_STRING diff --git a/src/tracking.F90 b/src/tracking.F90 index d108325fa0..d96269fc1d 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -184,7 +184,7 @@ contains if (coord % next % rotated) then ! If next level is rotated, apply rotation matrix coord % next % uvw = matmul(cells(coord % cell) % & - rotation, coord % uvw) + rotation_matrix, coord % uvw) else ! Otherwise, copy this level's direction coord % next % uvw = coord % uvw diff --git a/tests/test_infinite_cell/results.py b/tests/test_infinite_cell/results.py index be13ee66f1..e1c7836111 100644 --- a/tests/test_infinite_cell/results.py +++ b/tests/test_infinite_cell/results.py @@ -3,14 +3,13 @@ import sys # import statepoint -sys.path.insert(0, '../../src/utils') -import statepoint +from openmc.statepoint import StatePoint # read in statepoint file if len(sys.argv) > 1: - sp = statepoint.StatePoint(sys.argv[1]) + sp = StatePoint(sys.argv[1]) else: - sp = statepoint.StatePoint('statepoint.10.binary') + sp = StatePoint('statepoint.10.binary') sp.read_results() # set up output string