diff --git a/README.md b/README.md index a1f3e3e..e06a553 100644 --- a/README.md +++ b/README.md @@ -12,23 +12,26 @@ following features are already implemented: * hydrodynamic and magnetohydrodynamic set of equations (HD and MHD), * both classical and special relativity cases for the above equations, -* Cartesian coordinate system, +* Cartesian coordinate system so far, * uniform and adaptive mesh generation and update, -* 2nd to 4th order time integration using Strong Stability Preserving - Runge-Kutta methods, -* 2nd order TVD interpolation with number of limiters and higher order - reconstructions, +* a number of time integration methods, from 2nd to 5th order Runge-Kutta + methods: Strong Stability Preserving and Embedded (with the error control), +* high order reconstructions: from 2nd to 9th order WENO and MP, both explicit + and compact methods, the 2nd order TVD interpolation has a number of limiters + supported, * Riemann solvers of Roe- and HLL-types (HLL, HLLC, and HLLD), * standard boundary conditions: periodic, open, reflective, hydrostatic, etc. * turbulence driving using Alvelius or Ornstein–Uhlenbeck methods, * viscous and resistive source terms, -* support for passive scalars (up to 100), -* data stored in internal XML+binary or HDF5 format, -* support for Zstandard, LZ4, and LZMA compression in XML+binary format, -* Python interface to read snapshots in both formats, +* support for passive scalars, +* data stored in an internal XML+binary or the HDF5 format, +* data integrity of the XML+binary format guaranteed by the XXH64 or XXH3 hashes; +* support for Zstandard, LZ4, and LZMA compressions in the XML+binary format, +* support for Deflate, SZIP, Zstandard, and ZFP compressions in the HDF5 format, +* easy and consistend Python interface to read snapshots in both formats, * MPI parallelization, * completely written in Fortran 2008, -* simple Makefile or CMake for executable building, +* simple Makefile or CMake for building the code executable, * minimum requirements, only Fortran compiler and Python are required to prepare, run, and analyze your simulations. @@ -62,63 +65,59 @@ Requirements compiler version 9.0 or newer. - [NVIDIA HPC](https://developer.nvidia.com/hpc-sdk) compiler version 21.9. Warning: I could not make it run with the included MPI libraries. -* Optional, but recommended, [OpenMPI](https://www.open-mpi.org/) for parallel - runs, tested with version 1.8 or newer. -* Optional support for XML-binary format compression requires: +* Recommended, although optional, [OpenMPI](https://www.open-mpi.org/) for + parallel runs, tested with version 1.8 or newer. +* Optional [CMake](https://cmake.org) version 3.16 or newer, for advanced + compilation option selection. +* Optionally, the XML-binary format compression requires: [LZ4 library](https://lz4.github.io), [Zstandard library](http://facebook.github.io/zstd/), or [LZMA library](https://tukaani.org/xz/) + [XXHASH library](http://www.xxhash.com/). * Optional [HDF5 libraries](https://www.hdfgroup.org/solutions/hdf5/), tested with version 1.10 or newer. The code now uses the new XML-binary snapshot format. However, if you still want to use older HDF5 snapshot format, you will need these libraries. -* Optional [CMake](https://cmake.org) version 3.16 or newer, for managing the - build process. - - -Environment Variables -===================== - -If you need to use the HDF5 libraries and they are not installed in the default -location, i.e. in the system directory **/usr**, make sure that the environment -variable _HDF5DIR_ is set in your **~/.bashrc** (or **~/.cshrc**) and pointing -to the location where the HDF5 libraries have been installed. +* Deflate compression is natively supported in HDF5 libraries, however, + optionally these compression formats are supported through filters: + [SZIP](https://support.hdfgroup.org/doc_resource/SZIP/) + [HDF5Plugin-Zstandard](https://github.com/gkowal/HDF5Plugin-Zstandard), + [H5Z-ZFP](https://github.com/LLNL/H5Z-ZFP). Recommended compilation (using CMake) ===================================== 1. Clone the AMUN source code: - - from Bitbucket: - `git clone https://grzegorz_kowal@bitbucket.org/amunteam/amun-code.git`, - from GitLab: `git clone https://gitlab.com/gkowal/amun-code.git` + - from Bitbucket: + `git clone https://grzegorz_kowal@bitbucket.org/amunteam/amun-code.git`, - or unpack the archive downloaded from page [Downloads](https://bitbucket.org/amunteam/amun-code/downloads/). -2. Create a directory for compilation in any location, - e.g. `mkdir cmake-build && cmake-build`. +2. Create the build directory, e.g. `mkdir amun-build && cd amun-build`. 3. Call `ccmake `, e.g. `ccmake ..`, and press 'c' once. - Configure available options. Press 'c' once again, and 'g' to - generate makefiles. Alternatively, just call `ccmake ` - for default options. + Set available options, if necessary. Press 'c' once again, and 'g' to + generate makefiles. 4. Compile the code using `make`. The executable file **amun.x** should be - created. + available in a few moments. -Alternative compilation (using `make.conf`) +Alternative compilation (using `make`) =========================================== 1. Clone the AMUN source code: - - from Bitbucket: - `git clone https://grzegorz_kowal@bitbucket.org/amunteam/amun-code.git`, - from GitLab: `git clone https://gitlab.com/gkowal/amun-code.git` + - from Bitbucket: + `git clone https://grzegorz_kowal@bitbucket.org/amunteam/amun-code.git`, - or unpack the archive downloaded from page [Downloads](https://bitbucket.org/amunteam/amun-code/downloads/). -2. Go to directory **build/hosts/** and copy file **default** to a new file named - exactly as your host name, i.e. `cp default $HOSTNAME`. +2. Go to directory **build/hosts/** and copy file **default** to a new file + named exactly as your host name, i.e. `cp default $HOSTNAME`. 3. Customize your compiler and compilation options in your new host file. -4. Go up to directory **build/** and copy file **make.default** to **make.config**. +4. Go up to the directory **build/** and copy file **make.default** to + **make.config**. 5. Customize compilation time options in **make.config**. 6. Compile sources by typing `make` in directory **build/**. The executable file **amun.x** should be created there. @@ -147,23 +146,24 @@ where N is the number of processors to use. Reading data ============ -By default, the code uses new XML+binary snapshot data format. It can also be -forced by setting parameter **snapshot_format** to **xml**. +By default, the code uses the new XML+binary snapshot data format. Parameter +**snapshot_format** set to either **xml** or **h5** controls which file format +is used. -In order to read produced data in this format, you will need to install the -provided Python module. Simply change to **python/** directory and run - `python setup.py install --user` +In order to read the data produced in this format, you will need to install the +Python module AmunPy included in subdirectory **python/amunpy**. Simply go to +this directory and run + `python ./setup.py install --user` to install the module in your home directory. Import the module in your python script using `from amunpy import *`, -and then initiate the interface using +and then initiate the interface to the XML+binary snapshots using `snapshot = AmunXML()` -and read desired variable using +or to the HDF5 files using + `snapshot = AmunH5()` +and read desired variables using function `var = snapshot.dataset()`. -The function **dataset()** returns rescaled uniform mesh variable as NumPy -array. - -If you want to read data from HDF5 snapshot, just use - `var = amun_dataset(, )`. +The function **dataset()** returns the requested variable mapped on the uniform +mesh as a NumPy array. diff --git a/python/amunpy/src/amunpy/amunxml.py b/python/amunpy/src/amunpy/amunxml.py index 5bdbfda..78d4902 100644 --- a/python/amunpy/src/amunpy/amunxml.py +++ b/python/amunpy/src/amunpy/amunxml.py @@ -208,14 +208,18 @@ class AmunXML(Amun): return self.__swap__(dset) - def __check_digest(self, filename, digest, data): + def __check_digest(self, filename, hash_type, digest, data): ''' - Verifies if the provided digest matches the XXH64 hash of data - stored in the given filename. + Verifies if the provided digest matches the data. ''' import xxhash - if digest.lower() != xxhash.xxh64(data).hexdigest(): + failed = False + if hash_type == 'xxh64': + failed = digest.lower() != xxhash.xxh64(data).hexdigest() + elif hash_type == 'xxh3': + failed = digest.lower() != xxhash.xxh3_64(data).hexdigest() + if failed: print("File '{}' seems to be corrupted! Proceeding anyway...".format(filename)) @@ -235,7 +239,9 @@ class AmunXML(Amun): if 'compression_format' in self.binaries[dataset]: if 'compressed_digest' in self.binaries[dataset]: - self.__check_digest(fname, self.binaries[dataset]['compressed_digest'], stream) + htype = self.binaries[dataset]['digest_type'] + dhash = self.binaries[dataset]['compressed_digest'] + self.__check_digest(fname, htype, dhash, stream) comp = self.binaries[dataset]['compression_format'] if comp == 'zstd': @@ -248,12 +254,16 @@ class AmunXML(Amun): raise Exception("Binary file '{}' compressed in unsupported format {}!".format(fname, comp)) if 'digest' in self.binaries[dataset]: - self.__check_digest(fname, self.binaries[dataset]['digest'], data) + htype = self.binaries[dataset]['digest_type'] + dhash = self.binaries[dataset]['digest'] + self.__check_digest(fname, htype, dhash, data) return numpy.frombuffer(data, dtype=dtype) else: if 'digest' in self.binaries[dataset]: - self.__check_digest(fname, self.binaries[dataset]['digest'], stream) + htype = self.binaries[dataset]['digest_type'] + dhash = self.binaries[dataset]['digest'] + self.__check_digest(fname, htype, dhash, stream) return numpy.frombuffer(stream, dtype=dtype) else: @@ -276,7 +286,9 @@ class AmunXML(Amun): if 'compression_format' in self.chunks[chunk_number][dataset_name]: if 'compressed_digest' in self.chunks[chunk_number][dataset_name]: - self.__check_digest(fname, self.chunks[chunk_number][dataset_name]['compressed_digest'], stream) + htype = self.chunks[chunk_number][dataset_name]['digest_type'] + dhash = self.chunks[chunk_number][dataset_name]['compressed_digest'] + self.__check_digest(fname, htype, dhash, stream) comp = self.chunks[chunk_number][dataset_name]['compression_format'] if comp == 'zstd': @@ -290,12 +302,16 @@ class AmunXML(Amun): raise Exception("Binary file '{}' compressed in unsupported format {}!".format(fname, comp)) if 'digest' in self.chunks[chunk_number][dataset_name]: - self.__check_digest(fname, self.chunks[chunk_number][dataset_name]['digest'], data) + htype = self.chunks[chunk_number][dataset_name]['digest_type'] + dhash = self.chunks[chunk_number][dataset_name]['digest'] + self.__check_digest(fname, htype, dhash, data) return numpy.frombuffer(data, dtype=dtype) else: if 'digest' in self.chunks[chunk_number][dataset_name]: - self.__check_digest(fname, self.chunks[chunk_number][dataset_name]['digest'], stream) + htype = self.chunks[chunk_number][dataset_name]['digest_type'] + dhash = self.chunks[chunk_number][dataset_name]['digest'] + self.__check_digest(fname, htype, dhash, stream) return numpy.frombuffer(stream, dtype=dtype) else: diff --git a/sources/compression.F90 b/sources/compression.F90 index 081cfe5..967022f 100644 --- a/sources/compression.F90 +++ b/sources/compression.F90 @@ -36,15 +36,15 @@ module compression #ifdef ZSTD interface - integer(c_size_t) function zstd_bound(srcSize) & + integer(c_size_t) function zstd_bound(srcSize) & bind(C, name="ZSTD_compressBound") use iso_c_binding, only: c_size_t implicit none integer(kind=c_size_t), value :: srcSize end function zstd_bound - integer(c_size_t) function zstd_compress(dst, dstCapacity, & - src, srcSize, level) & + integer(c_size_t) function zstd_compress(dst, dstCapacity, & + src, srcSize, level) & bind(C, name="ZSTD_compress") use iso_c_binding, only: c_size_t, c_int, c_ptr implicit none @@ -53,6 +53,12 @@ module compression integer(kind=c_int) , value :: level end function zstd_compress + integer function zstd_iserror(code) bind(C, name="ZSTD_isError") + use iso_c_binding, only: c_size_t + implicit none + integer(kind=c_size_t), value :: code + end function zstd_iserror + end interface #endif /* ZSTD */ #ifdef LZ4 @@ -75,6 +81,12 @@ module compression type(c_ptr) , value :: src, dst, prefsPtr end function lz4_compress + integer function lz4_iserror(code) bind(C, name="LZ4F_isError") + use iso_c_binding, only: c_size_t + implicit none + integer(kind=c_size_t), value :: code + end function lz4_iserror + end interface #endif /* LZ4 */ #ifdef LZMA @@ -94,21 +106,25 @@ module compression end interface #endif /* LZMA */ -! compression parameters -! - integer, save :: compression_format = 0 - integer, save :: compression_level = 0 - ! supported compression formats ! - integer, parameter :: compression_none = 0 - integer, parameter :: compression_zstd = 1 - integer, parameter :: compression_lz4 = 2 - integer, parameter :: compression_lzma = 3 + enum, bind(c) + enumerator compression_none + enumerator compression_zstd + enumerator compression_lz4 + enumerator compression_lzma + end enum + +! compression parameters +! + integer(kind(compression_none)), save :: compression_format = 0 + integer , save :: compression_level = 0 + character(len=4) , save :: compression_suffix = '' private - public :: set_compression, get_compression, compress + public :: set_compression, get_compression, compression_bound, compress + public :: compression_suffix !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -131,17 +147,15 @@ module compression ! ! cformat - the compression format string; ! clevel - the compression level; -! suffix - the compressed file suffix; ! !=============================================================================== ! - subroutine set_compression(cformat, clevel, suffix) + subroutine set_compression(cformat, clevel) implicit none character(len=*) , intent(inout) :: cformat integer , intent(in) :: clevel - character(len=8) , intent(out) :: suffix !------------------------------------------------------------------------------- ! @@ -151,27 +165,27 @@ module compression cformat = "zstd" compression_format = compression_zstd compression_level = max(0, min(19, clevel)) - suffix = ".zst" + compression_suffix = ".zst" #endif /* ZSTD */ #ifdef LZ4 case("lz4", "LZ4") cformat = "lz4" compression_format = compression_lz4 compression_level = max(1, min(12, clevel)) - suffix = ".lz4" + compression_suffix = ".lz4" #endif /* LZ4 */ #ifdef LZMA case("lzma", "LZMA", "xz", "XZ") cformat = "lzma" compression_format = compression_lzma compression_level = max(0, min(9, clevel)) - suffix = ".xz" + compression_suffix = ".xz" #endif /* LZMA */ case default cformat = "none" compression_format = compression_none compression_level = clevel - suffix = "" + compression_suffix = "" end select !------------------------------------------------------------------------------- @@ -203,6 +217,56 @@ module compression ! !=============================================================================== ! +! function COMPRESSION_BOUND: +! -------------------------- +! +! Function returns the minimum buffer size required to perform +! the compression. +! +! Arguments: +! +! ilen - the length of the sequence of bytes to compress; +! +!=============================================================================== +! + integer(kind=8) function compression_bound(ilen) + + use iso_c_binding, only: c_loc + + implicit none + + integer(kind=8), intent(in) :: ilen + +#ifdef LZ4 + integer, dimension(14), target :: prefs = [5, 0, 1, 0, 0, 0, 0, 0, & + 0, 0, 0, 0, 0, 0] +#endif /* LZ4 */ + +!------------------------------------------------------------------------------- +! + select case(compression_format) +#ifdef ZSTD + case(compression_zstd) + compression_bound = zstd_bound(ilen) +#endif /* ZSTD */ +#ifdef LZ4 + case(compression_lz4) + prefs(5:6) = transfer(ilen, [ 0_4 ]) + prefs(9) = compression_level + compression_bound = lz4_bound(ilen, c_loc(prefs)) +#endif /* LZ4 */ + case default + compression_bound = ilen + end select + + return + +!------------------------------------------------------------------------------- +! + end function compression_bound +! +!=============================================================================== +! ! subroutine COMPRESS: ! ------------------- ! @@ -210,80 +274,65 @@ module compression ! ! Arguments: ! -! input - the input sequence of bytes; +! input - the pointer to the input sequence of bytes; +! ilen - the length of input; +! buffer - the buffer where the compression is done; +! clen - the length of buffer in bytes; once the compression was +! successful, it returns the length of the compressed stream; ! !=============================================================================== ! - subroutine compress(input, output, csize) + subroutine compress(input, ilen, buffer, clen) - use iso_c_binding, only: c_int, c_loc + use iso_c_binding, only : c_int, c_loc, c_ptr #ifdef LZ4 - use iso_c_binding, only: c_null_ptr + use iso_c_binding, only : c_null_ptr #endif /* LZ4 */ implicit none - integer(kind=1), dimension(:), target, intent(in) :: input - integer(kind=1), dimension(:), target, intent(out) :: output - integer(kind=8) , target, intent(out) :: csize + type(c_ptr) , intent(in) :: input + integer(kind=8), intent(in) :: ilen + type(c_ptr) , intent(inout) :: buffer + integer(kind=8), intent(inout) :: clen #ifdef LZ4 - integer, dimension(14), target :: prefs = [5, 0, 1, 0, 0, 0, 0, 0, & + integer, dimension(14), target :: prefs = [5, 0, 1, 0, 0, 0, 0, 0, & 0, 0, 0, 0, 0, 0] #endif /* LZ4 */ #ifdef LZMA - integer :: ret + integer(kind=8), target :: lsize + integer :: ret #endif /* LZMA */ - integer(kind=1), dimension(:), allocatable, target :: buffer - !------------------------------------------------------------------------------- ! - csize = min(size(input, kind=8), size(output, kind=8)) select case(compression_format) #ifdef ZSTD case(compression_zstd) - allocate(buffer(zstd_bound(size(input, kind=8)))) - csize = zstd_compress(c_loc(buffer), size(buffer, kind=8), & - c_loc(input), size(input, kind=8), & - compression_level) - if (csize > 0 .and. csize <= size(output, kind=8)) then - output(1:csize) = buffer(1:csize) - else - csize = -1 - end if - deallocate(buffer) + clen = zstd_compress(buffer, clen, input, ilen, compression_level) + if (zstd_iserror(clen) /= 0) clen = 0 #endif /* ZSTD */ #ifdef LZ4 case(compression_lz4) - prefs(5:6) = transfer(size(input, kind=8), [ 0_4 ]) + prefs(5:6) = transfer(ilen, [ 0_4 ]) prefs(9) = compression_level - allocate(buffer(lz4_bound(size(input, kind=8), c_loc(prefs)))) - csize = lz4_compress(c_loc(buffer), size(buffer, kind=8), & - c_loc(input), size(input, kind=8), c_loc(prefs)) - if (csize > 0 .and. csize <= size(output, kind=8)) then - output(1:csize) = buffer(1:csize) - else - csize = -1 - end if - deallocate(buffer) + clen = lz4_compress(buffer, clen, input, ilen, c_loc(prefs)) + if (lz4_iserror(clen) /= 0) clen = 0 #endif /* LZ4 */ #ifdef LZMA case(compression_lzma) - csize = 0 - allocate(buffer(size(input))) - ret = lzma_compress(compression_level, 4, c_null_ptr, & - c_loc(input), size(input, kind=8), & - c_loc(buffer), c_loc(csize), size(buffer, kind=8)) - if (ret == 0 .and. csize <= size(output, kind=8)) then - output(1:csize) = buffer(1:csize) + lsize = 0 + ret = lzma_compress(compression_level, 4, c_null_ptr, & + input, ilen, buffer, c_loc(lsize), clen) + if (ret > 0) then + clen = 0 else - csize = -1 + clen = lsize end if - deallocate(buffer) #endif /* LZMA */ case default - output(1:csize) = input(1:csize) + clen = 0 end select !------------------------------------------------------------------------------- diff --git a/sources/forcing.F90 b/sources/forcing.F90 index 080e001..1978ff8 100644 --- a/sources/forcing.F90 +++ b/sources/forcing.F90 @@ -113,7 +113,7 @@ module forcing ! array for driving mode coefficients ! - complex(kind=8), dimension(:,:), allocatable :: fcoefs + complex(kind=8), dimension(:,:), allocatable, target :: fcoefs ! velocity Fourier coefficients in Alfvelius method ! diff --git a/sources/hash.F90 b/sources/hash.F90 index 385e812..7157341 100644 --- a/sources/hash.F90 +++ b/sources/hash.F90 @@ -24,9 +24,10 @@ !! !! module: HASH !! -!! This module provides 64-bit version of the xxHash64 by Yann Collet. -!! This is a Fortran implementation based on the XXH64 specification -!! published at +!! This module is an interface to the XXH functions by Yann Collet provided +!! by the library libxxhash. If this library is not available, an internal +!! Fortran implementation of the 64-bit version of the xxHash64 is used. +!! The Fortran implementation is based on the XXH64 specification published at !! https://github.com/Cyan4973/xxHash/blob/dev/doc/xxhash_spec.md !! !! For additional info, see @@ -39,58 +40,56 @@ module hash implicit none #ifdef XXHASH -! interfaces to functions XXH64() and XXH3_64bits() provided -! by the systems library libxxhash +! interfaces to functions XXH64() and XXH3_64bits() provided by +! the library libxxhash ! interface - - integer(c_int64_t) function xxh64_system(input, length, seed) & - bind(C, name="XXH64") + integer(c_int64_t) function xxh64_lib(input, length, seed) & + bind(C, name="XXH64") use iso_c_binding, only: c_ptr, c_size_t, c_int64_t implicit none type(c_ptr) , value :: input integer(kind=c_size_t), value :: length integer(c_int64_t) , value :: seed - end function xxh64_system + end function xxh64_lib - integer(c_int64_t) function xxh3_system(input, length) & - bind(C, name="XXH3_64bits") + integer(c_int64_t) function xxh3_lib(input, length) & + bind(C, name="XXH3_64bits") use iso_c_binding, only: c_ptr, c_size_t, c_int64_t implicit none type(c_ptr) , value :: input integer(kind=c_size_t), value :: length - end function xxh3_system - + end function xxh3_lib end interface #else /* XXHASH */ ! hash parameters ! - integer(kind=8), parameter :: prime1 = -7046029288634856825_8, & - prime2 = -4417276706812531889_8, & - prime3 = 1609587929392839161_8, & - prime4 = -8796714831421723037_8, & - prime5 = 2870177450012600261_8, & + integer(kind=8), parameter :: prime1 = -7046029288634856825_8, & + prime2 = -4417276706812531889_8, & + prime3 = 1609587929392839161_8, & + prime4 = -8796714831421723037_8, & + prime5 = 2870177450012600261_8, & prime6 = 6983438078262162902_8 #endif /* XXHASH */ ! supported hash types ! - integer, parameter :: hash_xxh64 = 1 + enum, bind(c) + enumerator hash_none + enumerator hash_xxh64 #ifdef XXHASH - integer, parameter :: hash_xxh3 = 2 + enumerator hash_xxh3 #endif /* XXHASH */ + end enum private - public :: get_hash, check_hash, hash_string, hash_xxh64 -#ifdef XXHASH - public :: hash_xxh3 -#endif /* XXHASH */ + public :: hash_info, hash_name, digest, check_digest + public :: digest_string, digest_integer !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! contains ! -! !=============================================================================== !! !!*** PUBLIC SUBROUTINES ***************************************************** @@ -99,122 +98,242 @@ module hash ! !=============================================================================== ! -! function HASH_STRING: +! subroutine HASH_INFO: ! -------------------- ! -! Function returns the hash type string. +! Subroutine returns the hash ID and the length in bytes (characters) +! by the provided hash name. ! !=============================================================================== ! - character(len=8) function hash_string(hash_type) + subroutine hash_info(hash_name, hash_id, hash_length) + + use helpers, only : print_message + + implicit none + + character(len=*), intent(in) :: hash_name + integer , intent(out) :: hash_id, hash_length + + character(len=*), parameter :: loc = "HASH::hash_info()" + +!------------------------------------------------------------------------------- +! + select case(trim(hash_name)) + case("xxh64", "XXH64") + hash_id = hash_xxh64 + hash_length = 16 +#ifdef XXHASH + case("xxh3", "XXH3") + hash_id = hash_xxh3 + hash_length = 16 +#endif /* XXHASH */ + case("none") + hash_id = hash_none + hash_length = 0 + case default + call print_message(loc, & + "Hash '" // trim(hash_name) // "' is not supported!") + end select + + return + +!------------------------------------------------------------------------------- +! + end subroutine hash_info +! +!=============================================================================== +! +! function HASH_NAME: +! ------------------ +! +! Function returns the hash name by the provided hash ID. +! +!=============================================================================== +! + character(len=8) function hash_name(hash_type) implicit none integer, intent(in) :: hash_type + integer(kind(hash_none)) :: htype + !------------------------------------------------------------------------------- ! + htype = hash_type + select case(htype) + case(hash_xxh64) + hash_name = "xxh64" #ifdef XXHASH - if (hash_type == hash_xxh3) then - hash_string = 'xxh3' - return - end if + case(hash_xxh3) + hash_name = "xxh3" #endif /* XXHASH */ - hash_string = 'xxh64' + case default + hash_name = "none" + end select + return !------------------------------------------------------------------------------- ! - end function hash_string + end function hash_name ! !=============================================================================== ! -! function GET_HASH: -! ----------------- +! function DIGEST: +! --------------- ! -! Function calculates the hash for a given sequence of bytes. +! Function calculates the digest for a given sequence of bytes. ! ! Arguments: ! -! input - the input sequence of bytes; -! hash_type - the number corresponding to the hash type; +! buffer - the buffer pointer; +! length - the buffer length; +! hash_id - the hash ID; ! !=============================================================================== ! - integer(kind=8) function get_hash(input, hash_type) result(hash) + integer(kind=8) function digest(buffer, length, hash_id) result(hash) -#ifdef XXHASH - use iso_c_binding, only: c_loc -#endif /* XXHASH */ + use iso_c_binding, only : c_ptr implicit none - integer(kind=1), dimension(:), target, intent(in) :: input - integer , intent(in) :: hash_type - -#ifdef XXHASH - integer(kind=8) :: length -#endif /* XXHASH */ + type(c_ptr) , intent(in) :: buffer + integer(kind=8), intent(in) :: length + integer , intent(in) :: hash_id !------------------------------------------------------------------------------- ! -#ifdef XXHASH - hash = 0 - length = size(input, kind=8) - - select case(hash_type) - case(hash_xxh3) - hash = xxh3_system(c_loc(input), length) - case default - hash = xxh64_system(c_loc(input), length, hash) - end select + select case(hash_id) + case(hash_xxh64) +#ifndef XXHASH + hash = xxh64(buffer, length) #else /* XXHASH */ - hash = xxh64(input) + hash = xxh64_lib(buffer, length, 0_8) + case(hash_xxh3) + hash = xxh3_lib(buffer, length) #endif /* XXHASH */ + case(hash_none) + hash = 0 + end select return !------------------------------------------------------------------------------- ! - end function get_hash + end function digest ! !=============================================================================== ! -! subroutine CHECK_HASH: -! --------------------- +! subroutine CHECK_DIGEST: +! ----------------------- ! ! Subroutine checks if the provided digest matches the digest of ! the input data. ! ! Arguments: ! -! loc - the location of check; -! fname - the file name; -! input - the input sequence of bytes; -! digest - the data digest to check; -! hash_type - the number corresponding to the hash type; +! loc - the location of check; +! fname - the file name; +! buffer - the buffer pointer; +! length - the buffer length; +! bdigest - the buffer digest to check; +! hash_id - the hash ID; ! !=============================================================================== ! - subroutine check_hash(loc, fname, input, digest, hash_type) + subroutine check_digest(loc, fname, buffer, length, bdigest, hash_id) + + use helpers , only : print_message + use iso_c_binding, only : c_ptr + + implicit none + + character(len=*), intent(in) :: loc, fname + type(c_ptr) , intent(in) :: buffer + integer(kind=8) , intent(in) :: length + integer(kind=8) , intent(in) :: bdigest + integer , intent(in) :: hash_id + +!------------------------------------------------------------------------------- +! + if (hash_id == hash_none) return + + if (bdigest /= digest(buffer, length, hash_id)) & + call print_message(loc, trim(fname) // " seems to be corrupted!") + +!------------------------------------------------------------------------------- +! + end subroutine check_digest +! +!=============================================================================== +! +! subroutine DIGEST_STRING: +! ------------------------ +! +! Subroutine converts the integer digest to string. +! +! Arguments: +! +! idigest - the digest as integer; +! sdigest - the digest as string; +! +!=============================================================================== +! + subroutine digest_string(idigest, sdigest) use helpers, only : print_message implicit none - character(len=*) , intent(in) :: loc, fname - integer(kind=1), dimension(:), intent(in) :: input - integer(kind=8) , intent(in) :: digest - integer , intent(in) :: hash_type + integer(kind=8) , intent(in) :: idigest + character(len=*), intent(inout) :: sdigest + + character(len=*), parameter :: loc = "HASH::digest_string()" !------------------------------------------------------------------------------- ! - if (digest /= get_hash(input, hash_type)) & - call print_message(loc, trim(fname) // " seems to be corrupted!") + if (len(sdigest) >= 16) then + write(sdigest,"(1z16.16)") idigest + else + call print_message(loc, & + "The string is too short to contain the whole digest!") + end if !------------------------------------------------------------------------------- ! - end subroutine check_hash + end subroutine digest_string +! +!=============================================================================== +! +! subroutine DIGEST_INTEGER: +! ------------------------- +! +! Subroutine converts the string digest to its integer representation. +! +! Arguments: +! +! sdigest - the digest as string; +! idigest - the digest as integer; +! +!=============================================================================== +! + subroutine digest_integer(sdigest, idigest) + + implicit none + + character(len=*), intent(in) :: sdigest + integer(kind=8) , intent(out) :: idigest + +!------------------------------------------------------------------------------- +! + read(sdigest, fmt="(1z16)") idigest + +!------------------------------------------------------------------------------- +! + end subroutine digest_integer ! !=============================================================================== !! @@ -232,24 +351,30 @@ module hash ! ! Arguments: ! -! input - the input sequence of bytes; +! buffer - the buffer pointer; +! length - the buffer length; ! !=============================================================================== ! - integer(kind=8) function xxh64(input) result(hash) + integer(kind=8) function xxh64(buffer, length) result(hash) + + use iso_c_binding, only : c_ptr, c_f_pointer implicit none - integer(kind=1), dimension(:), target, intent(in) :: input + type(c_ptr) , intent(in) :: buffer + integer(kind=8), intent(in) :: length - integer(kind=8) :: length integer(kind=8) :: remaining, offset integer(kind=8), dimension(4) :: lane, chunk + integer(kind=1), dimension(:), pointer :: input + !------------------------------------------------------------------------------- ! - length = size(input, kind=8) + call c_f_pointer(buffer, input, [ length ]) + hash = 0_8 offset = 1_8 remaining = length diff --git a/sources/io.F90 b/sources/io.F90 index a9e3ec6..07f2233 100644 --- a/sources/io.F90 +++ b/sources/io.F90 @@ -116,28 +116,33 @@ module io character(len=255), save :: cformat = "none" ! compression format integer , save :: clevel = 3 ! compression level -! the suffix of binary files in the XML+binary format +! the type of digest to use and its length ! - character(len=8) , save :: binary_file_suffix = ".bin" - -! the type of digest to use -! - integer , save :: hash_type = 1 - character(len=8) , save :: hash_str = 'xxh64' + integer, save :: hash_type = 0 + integer, save :: hash_length = 0 #ifdef HDF5 -! compression type +! supported compression types ! - integer , parameter :: H5Z_DEFLATE = 1, H5Z_ZSTANDARD = 32015 + integer, parameter :: H5Z_DEFLATE = 1 + integer, parameter :: H5Z_SZIP = 4 + integer, parameter :: H5Z_ZFP = 32013 + integer, parameter :: H5Z_ZSTANDARD = 32015 -! compression type (0 for no compressions, 1 for deflate, 32015 for zstandard) +! used compression type and level ! - integer , save :: compression = 0, hclevel = 3 + integer, save :: hcformat = 0, hclevel = 20 + +! ZFP compressor parameters +! + character(len=32), save :: zfpmode = "reversible" + integer , save :: zfpprec = 64 + real(kind=8) , save :: zfprate = 6.4d+01 + real(kind=8) , save :: zfpaccu = 0.0d+00 ! HDF5 property object identifier ! - integer(hid_t) , save :: prp_id - + integer(hid_t), save :: prp_id #endif /* HDF5 */ ! array of pointer used during job restart @@ -179,10 +184,7 @@ module io subroutine initialize_io(verbose, status) use compression, only : set_compression, get_compression - use hash , only : hash_xxh64 -#ifdef XXHASH - use hash , only : hash_xxh3 -#endif /* XXHASH */ + use hash , only : hash_info use helpers , only : print_message use mpitools , only : nproc use parameters , only : get_parameter @@ -198,12 +200,10 @@ module io character(len=255) :: precise = "off" character(len=255) :: ghosts = "on" character(len=255) :: xdmf = "off" - character(len=255) :: suffix = "" ! compression file suffix character(len=8) :: dtype = "xxh64" #ifdef HDF5 - logical :: cmpstatus = .false. - integer(hsize_t) :: cd_nelmts = 1 - integer, dimension(1) :: cd_values = 3 + integer(hsize_t) :: cd_nelmts = 6 + integer, dimension(6) :: cd_values = 0 #endif /* HDF5 */ #ifdef HDF5 @@ -279,21 +279,10 @@ module io call get_parameter("compression_format", cformat) call get_parameter("compression_level" , clevel) - call set_compression(cformat, clevel, suffix) - if (get_compression() > 0) & - binary_file_suffix = ".bin" // trim(adjustl(suffix)) + call set_compression(cformat, clevel) call get_parameter("digest_type", dtype) - select case(dtype) -#ifdef XXHASH - case('xxh3', 'XXH3') - hash_type = hash_xxh3 - hash_str = 'xxh3' -#endif /* XXHASH */ - case default - hash_type = hash_xxh64 - hash_str = 'xxh64' - end select + call hash_info(dtype, hash_type, hash_length) if (status == 0) then @@ -323,31 +312,79 @@ module io "Cannot create the compression property for datasets!") else - cmpstatus = .false. - if (.not. cmpstatus) then - call h5zfilter_avail_f(H5Z_ZSTANDARD, cmpstatus, status) - if (cmpstatus) compression = H5Z_ZSTANDARD - end if - if (.not. cmpstatus) then - call h5zfilter_avail_f(H5Z_DEFLATE, cmpstatus, status) - if (cmpstatus) compression = H5Z_DEFLATE - end if + call get_parameter("compression_format", sformat) + call get_parameter("compression_level" , hclevel) + call get_parameter("zfp_mode" , zfpmode) + call get_parameter("zfp_rate" , zfprate) + call get_parameter("zfp_precision" , zfpprec) + call get_parameter("zfp_accuracy" , zfpaccu) - call get_parameter("compression_level", hclevel) + select case(sformat) + case("deflate", "gzip") + call h5zfilter_avail_f(H5Z_DEFLATE, test, status) + if (status == 0) then + if (test) then + hcformat = H5Z_DEFLATE + hclevel = max(1, min(9, hclevel)) + call h5pset_deflate_f(prp_id, hclevel, status) + end if + else + call print_message(loc, & + "Could not check if the filter is available!") + end if + case("szip") + call h5zfilter_avail_f(H5Z_FILTER_SZIP_F, test, status) + if (status == 0) then + if (test) then + hcformat = H5Z_SZIP + call h5pset_szip_f(prp_id, 32, 32, status) + end if + else + call print_message(loc, & + "Could not check if the filter is available!") + end if + case("zstd", "zstandard") + call h5zfilter_avail_f(H5Z_ZSTANDARD, test, status) + if (status == 0) then + if (test) then + hcformat = H5Z_ZSTANDARD + hclevel = max(1, min(20, hclevel)) + cd_values(:) = hclevel + call h5pset_filter_f(prp_id, H5Z_ZSTANDARD, & + H5Z_FLAG_OPTIONAL_F, cd_nelmts, cd_values, status) + end if + else + call print_message(loc, & + "Could not check if the filter is available!") + end if + case("zfp") + call h5zfilter_avail_f(H5Z_ZFP, test, status) + if (status == 0) then + if (test) then + hcformat = H5Z_ZFP + select case(trim(zfpmode)) + case('rate') + cd_values(1) = 1 + cd_values(3:4) = transfer(zfprate, [0_4]) + case('precision') + cd_values(1) = 2 + cd_values(3) = zfpprec + case('accuracy') + cd_values(1) = 3 + cd_values(3:4) = transfer(zfpaccu, [0_4]) + case('reversible') + cd_values(1) = 5 + end select + call h5pset_filter_f(prp_id, H5Z_ZFP, 0, & + cd_nelmts, cd_values, status) + end if + else + call print_message(loc, & + "Could not check if the filter is available!") + end if + case default + end select - if (status == 0) then - select case(compression) - case(H5Z_ZSTANDARD) - hclevel = max(1, min(20, hclevel)) - cd_values(:) = hclevel - call h5pset_filter_f(prp_id, H5Z_ZSTANDARD, H5Z_FLAG_OPTIONAL_F, & - cd_nelmts, cd_values, status) - case(H5Z_DEFLATE) - hclevel = max(1, min(9, hclevel)) - call h5pset_deflate_f(prp_id, hclevel, status) - case default - end select - end if end if end if #endif /* HDF5 */ @@ -417,6 +454,7 @@ module io ! subroutine print_io(verbose) + use hash , only : hash_name use helpers, only : print_section, print_parameter implicit none @@ -445,22 +483,35 @@ module io #endif /* HDF5 */ case default call print_parameter(verbose, "restart snapshot format", "XML+binary") - call print_parameter(verbose, "digest type", hash_str) + call print_parameter(verbose, "digest type", hash_name(hash_type)) call print_parameter(verbose, "compression format", cformat) call print_parameter(verbose, "compression level", clevel) end select call print_parameter(verbose, "precise snapshot intervals", & precise_snapshots, "on") #ifdef HDF5 - select case(compression) - case(H5Z_ZSTANDARD) - call print_parameter(verbose, "HDF5 compression" , "zstd" ) - call print_parameter(verbose, "compression level", hclevel ) + select case(hcformat) case(H5Z_DEFLATE) - call print_parameter(verbose, "HDF5 compression" , "deflate") - call print_parameter(verbose, "compression level", hclevel ) + call print_parameter(verbose, "HDF5 compression", "deflate") + call print_parameter(verbose, "compression level", hclevel) + case(H5Z_SZIP) + call print_parameter(verbose, "HDF5 compression", "szip") + case(H5Z_ZSTANDARD) + call print_parameter(verbose, "HDF5 compression", "zstd") + call print_parameter(verbose, "compression level", hclevel) + case(H5Z_ZFP) + call print_parameter(verbose, "HDF5 compression", "zfp") + call print_parameter(verbose, "ZFP mode", zfpmode) + select case(trim(zfpmode)) + case('rate') + call print_parameter(verbose, "ZFP rate" , zfprate) + case('precision') + call print_parameter(verbose, "ZFP precision", zfpprec) + case('accuracy') + call print_parameter(verbose, "ZFP accuracy" , zfpaccu) + end select case default - call print_parameter(verbose, "HDF5 compression" , "none" ) + call print_parameter(verbose, "HDF5 compression" , "none") end select call print_parameter(verbose, "generate XDMF files", with_xdmf, "on") #endif /* HDF5 */ @@ -996,31 +1047,32 @@ module io ! subroutine read_restart_snapshot_xml(status) - use blocks , only : block_meta, block_data, pointer_meta, list_meta - use blocks , only : ns => nsides, nc => nchildren, nregs - use blocks , only : append_metablock, append_datablock, link_blocks - use blocks , only : get_mblocks - use blocks , only : set_last_id, get_last_id - use blocks , only : metablock_set_id, metablock_set_process - use blocks , only : metablock_set_refinement - use blocks , only : metablock_set_configuration - use blocks , only : metablock_set_level, metablock_set_position - use blocks , only : metablock_set_coordinates, metablock_set_bounds - use blocks , only : metablock_set_leaf - use blocks , only : change_blocks_process - use coordinates, only : nn => bcells, ncells, nghosts - use coordinates, only : xmin, xmax, ymin, ymax, zmin, zmax - use equations , only : cmax, cmax2 - use evolution , only : step, time, dt, dth, dte - use evolution , only : niterations, nrejections, errs - use forcing , only : nmodes, fcoefs, einj - use hash , only : check_hash, hash_xxh64 - use helpers , only : print_message + use blocks , only : block_meta, block_data, pointer_meta, list_meta + use blocks , only : ns => nsides, nc => nchildren, nregs + use blocks , only : append_metablock, append_datablock, link_blocks + use blocks , only : get_mblocks + use blocks , only : set_last_id, get_last_id + use blocks , only : metablock_set_id, metablock_set_process + use blocks , only : metablock_set_refinement + use blocks , only : metablock_set_configuration + use blocks , only : metablock_set_level, metablock_set_position + use blocks , only : metablock_set_coordinates, metablock_set_bounds + use blocks , only : metablock_set_leaf + use blocks , only : change_blocks_process + use coordinates , only : nn => bcells, ncells, nghosts + use coordinates , only : xmin, xmax, ymin, ymax, zmin, zmax + use equations , only : cmax, cmax2 + use evolution , only : step, time, dt, dth, dte + use evolution , only : niterations, nrejections, errs + use forcing , only : nmodes, fcoefs, einj + use hash , only : hash_info, check_digest, digest_integer + use helpers , only : print_message + use iso_c_binding, only : c_loc #ifdef MPI - use mesh , only : redistribute_blocks + use mesh , only : redistribute_blocks #endif /* MPI */ - use mpitools , only : nprocs, nproc - use random , only : gentype, set_seeds + use mpitools , only : nprocs, nproc + use random , only : gentype, set_seeds implicit none @@ -1036,31 +1088,34 @@ module io type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata + integer :: dtype, dlen + integer(kind=4) :: lun = 104 - integer(kind=8) :: digest, bytes, pbytes, ubytes + integer(kind=8) :: bytes, pbytes, ubytes - character(len=16) :: hfield, hchild, hface, hedge, hcorner, hbound - character(len=16) :: hids, harray, hseed, hforce + integer(kind=8) :: hfield, hchild, hface, hedge, hcorner, hbound + integer(kind=8) :: hids, hseed, hforce - integer(kind=8), dimension(:) , allocatable :: hprim, hcons - integer(kind=4), dimension(:) , allocatable :: ids - integer(kind=4), dimension(:,:) , allocatable :: fields - integer(kind=4), dimension(:,:) , allocatable :: children + integer(kind=8), dimension(:) , allocatable :: hprim, hcons + integer(kind=4), dimension(:) , allocatable, target :: ids + integer(kind=4), dimension(:,:) , allocatable, target :: fields + integer(kind=4), dimension(:,:) , allocatable, target :: children #if NDIMS == 2 - integer(kind=4), dimension(:,:,:,:) , allocatable :: edges - integer(kind=4), dimension(:,:,:) , allocatable :: corners + integer(kind=4), dimension(:,:,:,:) , allocatable, target :: edges + integer(kind=4), dimension(:,:,:) , allocatable, target :: corners #endif /* NDIMS == 2 */ #if NDIMS == 3 - integer(kind=4), dimension(:,:,:,:,:), allocatable :: faces - integer(kind=4), dimension(:,:,:,:,:), allocatable :: edges - integer(kind=4), dimension(:,:,:,:) , allocatable :: corners + integer(kind=4), dimension(:,:,:,:,:), allocatable, target :: faces + integer(kind=4), dimension(:,:,:,:,:), allocatable, target :: edges + integer(kind=4), dimension(:,:,:,:) , allocatable, target :: corners #endif /* NDIMS == 3 */ - integer(kind=8), dimension(:,:) , allocatable :: seeds - real(kind=8) , dimension(:,:,:) , allocatable :: bounds - real(kind=8) , dimension(:,:,:,:,:), allocatable :: array + integer(kind=8), dimension(:,:) , allocatable, target :: seeds + real(kind=8) , dimension(:,:,:) , allocatable, target :: bounds + real(kind=8) , dimension(:,:,:,:,:), allocatable, target :: array + complex(kind=8), dimension(:,:) , allocatable, target :: lfcoefs character(len=*), parameter :: loc = 'IO::read_restart_snapshot_xml()' - character(len=*), parameter :: fmt = "(a,a,'_',i6.6,'.',a)" + character(len=*), parameter :: sfmt = "(a,a,'_',i6.6,'.',a)" !------------------------------------------------------------------------------- ! @@ -1069,9 +1124,9 @@ module io write(dname, "(a,'restart-',i5.5)") trim(respath), nrest #ifdef __INTEL_COMPILER - inquire(directory = dname, exist = test) + inquire(directory=dname, exist=test) #else /* __INTEL_COMPILER */ - inquire(file = dname, exist = test) + inquire(file=dname, exist=test) #endif /* __INTEL_COMPILER */ if (.not. test) then call print_message(loc, trim(dname) // " does not exist!") @@ -1081,15 +1136,15 @@ module io dname = trim(dname) // "/" write(fname,"(a,'metadata.xml')") trim(dname) - inquire(file = fname, exist = test) + inquire(file=fname, exist=test) if (.not. test) then call print_message(loc, trim(fname) // " does not exist!") status = 121 return end if - open(newunit = lun, file = fname, status = 'old') -10 read(lun, fmt = "(a)", end = 20) line + open(newunit=lun, file=fname, status='old') +10 read(lun, fmt="(a)", end=20) line if (index(line, ' 0) then il = index(line, 'name="') if (il > 0) then @@ -1097,99 +1152,103 @@ module io iu = index(line(il:), '"') + il - 2 write(sname,*) line(il:iu) il = index(line, '>') + 1 - iu = index(line, '<', back = .true.) - 1 + iu = index(line, '<', back=.true.) - 1 write(svalue,*) line(il:iu) select case(trim(adjustl(sname))) case('ndims') - read(svalue, fmt = *) lndims + read(svalue, fmt=*) lndims case('nprocs') - read(svalue, fmt = *) lnprocs + read(svalue, fmt=*) lnprocs case('nproc') - read(svalue, fmt = *) lnproc + read(svalue, fmt=*) lnproc case('mblocks') - read(svalue, fmt = *) lmblocks + read(svalue, fmt=*) lmblocks case('dblocks') - read(svalue, fmt = *) ldblocks + read(svalue, fmt=*) ldblocks case('nleafs') - read(svalue, fmt = *) lnleafs + read(svalue, fmt=*) lnleafs case('last_id') - read(svalue, fmt = *) llast_id + read(svalue, fmt=*) llast_id case('ncells') - read(svalue, fmt = *) lncells + read(svalue, fmt=*) lncells case('nghosts') - read(svalue, fmt = *) lnghosts + read(svalue, fmt=*) lnghosts case('nseeds') - read(svalue, fmt = *) lnseeds + read(svalue, fmt=*) lnseeds case('step') - read(svalue, fmt = *) step + read(svalue, fmt=*) step case('isnap') - read(svalue, fmt = *) isnap + read(svalue, fmt=*) isnap case('nvars') - read(svalue, fmt = *) nv + read(svalue, fmt=*) nv case('nmodes') - read(svalue, fmt = *) lnmodes + read(svalue, fmt=*) lnmodes case('xmin') - read(svalue, fmt = *) xmin + read(svalue, fmt=*) xmin case('xmax') - read(svalue, fmt = *) xmax + read(svalue, fmt=*) xmax case('ymin') - read(svalue, fmt = *) ymin + read(svalue, fmt=*) ymin case('ymax') - read(svalue, fmt = *) ymax + read(svalue, fmt=*) ymax case('zmin') - read(svalue, fmt = *) zmin + read(svalue, fmt=*) zmin case('zmax') - read(svalue, fmt = *) zmax + read(svalue, fmt=*) zmax case('time') - read(svalue, fmt = *) time + read(svalue, fmt=*) time case('dt') - read(svalue, fmt = *) dt + read(svalue, fmt=*) dt case('dth') - read(svalue, fmt = *) dth + read(svalue, fmt=*) dth case('dte') - read(svalue, fmt = *) dte + read(svalue, fmt=*) dte case('cmax') - read(svalue, fmt = *) cmax + read(svalue, fmt=*) cmax cmax2 = cmax * cmax case('niterations') - read(svalue, fmt = *) niterations + read(svalue, fmt=*) niterations case('nrejections') - read(svalue, fmt = *) nrejections + read(svalue, fmt=*) nrejections case('errs(1)') - read(svalue, fmt = *) errs(1) + read(svalue, fmt=*) errs(1) case('errs(2)') - read(svalue, fmt = *) errs(2) + read(svalue, fmt=*) errs(2) case('errs(3)') - read(svalue, fmt = *) errs(3) + read(svalue, fmt=*) errs(3) case('fields') + il = index(line, 'digest_type="') + 13 + iu = index(line(il:), '"') + il - 2 + call hash_info(line(il:iu), dtype, dlen) + il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hfield + call digest_integer(line(il:iu), hfield) case('children') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hchild + call digest_integer(line(il:iu), hchild) case('faces') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hface + call digest_integer(line(il:iu), hface) case('edges') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hedge + call digest_integer(line(il:iu), hedge) case('corners') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hcorner + call digest_integer(line(il:iu), hcorner) case('bounds') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hbound + call digest_integer(line(il:iu), hbound) case('forcing') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hforce + call digest_integer(line(il:iu), hforce) end select end if end if @@ -1226,66 +1285,48 @@ module io #else /* NDIMS == 3 */ edges(NDIMS,ns,ns,nm), corners(ns,ns,nm), & #endif /* NDIMS == 3 */ - block_array(nx), stat = status) + block_array(nx), stat=status) if (status == 0) then - bytes = size(transfer(fields, [ 0_1 ]), kind=8) write(fname,"(a,'metablock_fields.bin')") trim(dname) - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) fields close(lun) - read(hfield, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(fields, 1_1, bytes), & - digest, hash_xxh64) + bytes = size(fields, kind=8) * kind(fields) + call check_digest(loc, fname, c_loc(fields), bytes, hfield, dtype) write(fname,"(a,'metablock_children.bin')") trim(dname) - bytes = size(transfer(children, [ 0_1 ]), kind=8) - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) children close(lun) - read(hchild, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(children, 1_1, bytes), & - digest, hash_xxh64) + bytes = size(children, kind=8) * kind(children) + call check_digest(loc, fname, c_loc(children), bytes, hchild, dtype) #if NDIMS == 3 - bytes = size(transfer(faces, [ 0_1 ]), kind=8) write(fname,"(a,'metablock_faces.bin')") trim(dname) - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) faces close(lun) - read(hface, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(faces, 1_1, bytes), & - digest, hash_xxh64) + bytes = size(faces, kind=8) * kind(faces) + call check_digest(loc, fname, c_loc(faces), bytes, hface, dtype) #endif /* NDIMS == 3 */ - bytes = size(transfer(edges, [ 0_1 ]), kind=8) write(fname,"(a,'metablock_edges.bin')") trim(dname) - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) edges close(lun) - read(hedge, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(edges, 1_1, bytes), & - digest, hash_xxh64) - bytes = size(transfer(corners, [ 0_1 ]), kind=8) + bytes = size(edges, kind=8) * kind(edges) + call check_digest(loc, fname, c_loc(edges), bytes, hedge, dtype) write(fname,"(a,'metablock_corners.bin')") trim(dname) - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) corners close(lun) - read(hcorner, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(corners, 1_1, bytes), & - digest, hash_xxh64) - bytes = size(transfer(bounds, [ 0_1 ]), kind=8) + bytes = size(corners, kind=8) * kind(corners) + call check_digest(loc, fname, c_loc(corners), bytes, hcorner, dtype) write(fname,"(a,'metablock_bounds.bin')") trim(dname) - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) bounds close(lun) - read(hbound, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(bounds, 1_1, bytes), & - digest, hash_xxh64) + bytes = size(bounds, kind=8) * kind(bounds) + call check_digest(loc, fname, c_loc(bounds), bytes, hbound, dtype) l = 0 pmeta => list_meta @@ -1361,21 +1402,31 @@ module io deallocate(fields, children, bounds, edges, corners, stat=status) #endif /* NDIMS == 3 */ if (status /= 0) & - call print_message(loc, "Could not deallocate space of metablocks!") + call print_message(loc, "Could not release space of metablocks!") + else + call print_message(loc, "Could not allocate space of metablocks!") end if if (lnmodes == nmodes) then if (lnmodes > 0) then - bytes = size(transfer(fcoefs, [ 0_1 ]), kind=8) - write(fname,"(a,'forcing_coefficients.bin')") trim(dname) - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') - read(lun) fcoefs - close(lun) - read(hforce, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(fcoefs, 1_1, bytes), & - digest, hash_xxh64) + allocate(lfcoefs(lnmodes,lndims), stat=status) + if (status == 0) then + write(fname,"(a,'forcing_coefficients.bin')") trim(dname) + open(newunit=lun, file=fname, form='unformatted', access='stream') + read(lun) lfcoefs + close(lun) + bytes = size(lfcoefs, kind=8) * kind(lfcoefs) + call check_digest(loc, fname, c_loc(lfcoefs), bytes, hforce, dtype) + fcoefs = lfcoefs + deallocate(lfcoefs, stat=status) + if (status /= 0) & + call print_message(loc, & + "Could not release space of Fourier coefficients!") + else + call print_message(loc, & + "Could not allocate space of Fourier coefficients!") + end if end if else call print_message(loc, "The number of forcing modes does not match!") @@ -1385,16 +1436,16 @@ module io if (nproc < lnprocs) then - write(fname,fmt) trim(dname), "datablocks", nproc, "xml" - inquire(file = fname, exist = test) + write(fname,sfmt) trim(dname), "datablocks", nproc, "xml" + inquire(file=fname, exist=test) if (.not. test) then write(*,*) trim(fname) // " does not exist!" status = 121 return end if - open(newunit = lun, file = fname, status = 'old') -30 read(lun, fmt = "(a)", end = 40) line + open(newunit=lun, file=fname, status='old') +30 read(lun, fmt="(a)", end=40) line if (index(line, ' 0) then il = index(line, 'name="') if (il > 0) then @@ -1402,43 +1453,42 @@ module io iu = index(line(il:), '"') + il - 2 write(sname,*) line(il:iu) il = index(line, '>') + 1 - iu = index(line, '<', back = .true.) - 1 + iu = index(line, '<', back=.true.) - 1 write(svalue,*) line(il:iu) select case(trim(adjustl(sname))) case('dblocks') - read(svalue, fmt = *) nd + read(svalue, fmt=*) nd if (nd > 0) then - allocate(hprim(nd), hcons(nd), stat = status) + allocate(hprim(nd), hcons(nd), stat=status) + if (status /= 0) & + call print_message(loc, & + "Could not allocate space for hashes!") end if case('nregs') - read(svalue, fmt = *) nr + read(svalue, fmt=*) nr case('einj') - read(svalue, fmt = *) einj + read(svalue, fmt=*) einj case('ids') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hids - case('arrays') - il = index(line, 'digest="') + 8 - iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) harray + call digest_integer(line(il:iu), hids) case('seeds') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hseed + call digest_integer(line(il:iu), hseed) end select if (index(sname, 'prim') > 0) then - read(sname(7:), fmt = *) l + read(sname(7:), fmt=*) l il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = "(1z16)") hprim(l) + call digest_integer(line(il:iu), hprim(l)) end if if (index(sname, 'cons') > 0) then - read(sname(7:), fmt = *) l + read(sname(7:), fmt=*) l il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = "(1z16)") hcons(l) + call digest_integer(line(il:iu), hcons(l)) end if end if end if @@ -1457,38 +1507,34 @@ module io if (nd > 0) then #if NDIMS == 3 - allocate(ids(nd), array(nv,nm,nm,nm,nr), stat = status) + allocate(ids(nd), array(nv,nm,nm,nm,nr), stat=status) #else /* NDIMS == 3 */ - allocate(ids(nd), array(nv,nm,nm, 1,nr), stat = status) + allocate(ids(nd), array(nv,nm,nm, 1,nr), stat=status) #endif /* NDIMS == 3 */ if (status == 0) then - bytes = size(transfer(ids, [ 0_1 ]), kind=8) - write(fname, fmt) trim(dname), "datablock_ids", nproc, "bin" - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + write(fname,sfmt) trim(dname), "datablock_ids", nproc, "bin" + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) ids close(lun) - read(hids, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(ids, 1_1, bytes), & - digest, hash_xxh64) + bytes = size(ids, kind=8) * kind(ids) + call check_digest(loc, fname, c_loc(ids), bytes, hids, dtype) - pbytes = size(transfer(array(:,:,:,:,1), [ 0_1 ]), kind=8) - ubytes = size(transfer(array(:,:,:,:,:), [ 0_1 ]), kind=8) + ubytes = size(array, kind=8) * kind(array) + pbytes = ubytes / nr do l = 1, nd call append_datablock(pdata, status) call link_blocks(block_array(ids(l))%ptr, pdata) write(fname,"(a,'datablock_prim_',i6.6,'_',i6.6,'.bin')") & - trim(dname), nproc, l - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + trim(dname), nproc, l + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) array(:,:,:,:,1) close(lun) - call check_hash(loc, fname, & - transfer(array(:,:,:,:,1), 1_1, pbytes), hprim(l), hash_xxh64) + call check_digest(loc, fname, c_loc(array), & + pbytes, hprim(l), dtype) if (lnghosts >= nghosts) then #if NDIMS == 3 @@ -1505,13 +1551,12 @@ module io end if write(fname,"(a,'datablock_cons_',i6.6,'_',i6.6,'.bin')") & - trim(dname), nproc, l - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + trim(dname), nproc, l + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) array close(lun) - call check_hash(loc, fname, transfer(array, 1_1, ubytes), & - hcons(l), hash_xxh64) + call check_digest(loc, fname, c_loc(array), & + ubytes, hcons(l), dtype) p = min(nregs, nr) if (lnghosts >= nghosts) then @@ -1531,43 +1576,43 @@ module io deallocate(ids, array, hprim, hcons, stat=status) if (status /= 0) & - call print_message(loc, "Could not release memory!") + call print_message(loc, "Could not release space of datablocks!") else - call print_message(loc, "Could not allocate memory!") + call print_message(loc, "Could not allocate space for datablocks!") end if end if - allocate(seeds(4,lnseeds), stat = status) + allocate(seeds(4,lnseeds), stat=status) if (status == 0) then - bytes = size(transfer(seeds, [ 0_1 ]), kind=8) - write(fname, fmt) trim(dname), "random_seeds", nproc, "bin" - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + write(fname,sfmt) trim(dname), "random_seeds", nproc, "bin" + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) seeds close(lun) - read(hseed, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(seeds, 1_1, bytes), & - digest, hash_xxh64) + bytes = size(seeds, kind=8) * kind(seeds) + call check_digest(loc, fname, c_loc(seeds), bytes, hseed, dtype) call set_seeds(lnseeds, seeds(:,:), .false.) - if (allocated(seeds)) deallocate(seeds) - + deallocate(seeds, stat=status) + if (status /= 0) & + call print_message(loc, "Could not release space of seeds!") + else + call print_message(loc, "Could not allocate space for seeds!") end if else ! nproc < lnprocs - write(fname,fmt) trim(dname), "datablocks", 0, "xml" - inquire(file = fname, exist = test) + write(fname,sfmt) trim(dname), "datablocks", 0, "xml" + inquire(file=fname, exist=test) if (.not. test) then write(*,*) trim(fname) // " does not exist!" status = 121 return end if - open(newunit = lun, file = fname, status = 'old') -50 read(lun, fmt = "(a)", end = 60) line + open(newunit=lun, file=fname, status='old') +50 read(lun, fmt="(a)", end=60) line if (index(line, ' 0) then il = index(line, 'name="') if (il > 0) then @@ -1575,14 +1620,14 @@ module io iu = index(line(il:), '"') + il - 2 write(sname,*) line(il:iu) il = index(line, '>') + 1 - iu = index(line, '<', back = .true.) - 1 + iu = index(line, '<', back=.true.) - 1 write(svalue,*) line(il:iu) select case(trim(adjustl(sname))) case('seeds') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hseed + call digest_integer(line(il:iu), hseed) end select end if end if @@ -1593,24 +1638,24 @@ module io ! if (trim(gentype) == "same") then - allocate(seeds(4,lnseeds), stat = status) + allocate(seeds(4,lnseeds), stat=status) if (status == 0) then - bytes = size(transfer(seeds, [ 0_1 ]), kind=8) - write(fname, fmt) trim(dname), "random_seeds", 0, "bin" - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + write(fname,sfmt) trim(dname), "random_seeds", 0, "bin" + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) seeds close(lun) - read(hseed, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(seeds, 1_1, bytes), & - digest, hash_xxh64) + bytes = size(seeds, kind=8) * kind(seeds) + call check_digest(loc, fname, c_loc(seeds), bytes, hseed, dtype) call set_seeds(lnseeds, seeds(:,:), .false.) - if (allocated(seeds)) deallocate(seeds) - - end if ! allocation + deallocate(seeds, stat=status) + if (status /= 0) & + call print_message(loc, "Could not release space of seeds!") + else + call print_message(loc, "Could not allocate space for seeds!") + end if end if ! gentype == "same" @@ -1618,41 +1663,32 @@ module io else ! nprocs < lnprocs -! divide files between processes +! divide files between processes and update the block process accordingly ! nl = 0 - i = mod(lnprocs, nprocs) - j = lnprocs / nprocs - do p = 0, nprocs - k = 0 - do n = 0, p - nl = k - if (n < i) then - nu = k + j - else - nu = k + j - 1 - end if - k = nu + 1 - end do - do n = nl, nu - call change_blocks_process(n, p) - end do - end do - k = 0 - do n = 0, nproc - nl = k - if (n < i) then - nu = k + j + nd = lnprocs / nprocs + nr = mod(lnprocs, nprocs) + do n = 0, nprocs - 1 + if (n < nr) then + il = n * (nd + 1) + iu = il + nd else - nu = k + j - 1 + il = n * nd + nr + iu = il + nd - 1 + end if + do i = il, iu + call change_blocks_process(i, n) + end do + if (n == nproc) then + nl = il + nu = iu end if - k = nu + 1 end do do n = nl, nu - write(fname,fmt) trim(dname), "datablocks", n, "xml" - inquire(file = fname, exist = test) + write(fname,sfmt) trim(dname), "datablocks", n, "xml" + inquire(file=fname, exist=test) if (.not. test) then write(*,*) trim(fname) // " does not exist!" status = 121 @@ -1661,8 +1697,8 @@ module io ! read attributes from the metadata file ! - open(newunit = lun, file = fname, status = 'old') -70 read(lun, fmt = "(a)", end = 80) line + open(newunit=lun, file=fname, status='old') +70 read(lun, fmt="(a)", end=80) line if (index(line, ' 0) then il = index(line, 'name="') if (il > 0) then @@ -1670,29 +1706,43 @@ module io iu = index(line(il:), '"') + il - 2 write(sname,*) line(il:iu) il = index(line, '>') + 1 - iu = index(line, '<', back = .true.) - 1 + iu = index(line, '<', back=.true.) - 1 write(svalue,*) line(il:iu) select case(trim(adjustl(sname))) case('dblocks') - read(svalue, fmt = *) nd + read(svalue, fmt=*) nd + if (nd > 0) then + allocate(hprim(nd), hcons(nd), stat=status) + if (status /= 0) & + call print_message(loc, & + "Could not allocate space for hashes!") + end if case('nregs') - read(svalue, fmt = *) nr + read(svalue, fmt=*) nr case('einj') - read(svalue, fmt = *) deinj + read(svalue, fmt=*) deinj case('ids') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hids - case('arrays') - il = index(line, 'digest="') + 8 - iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) harray + call digest_integer(line(il:iu), hids) case('seeds') il = index(line, 'digest="') + 8 iu = index(line(il:), '"') + il - 2 - read(line(il:iu), fmt = *) hseed + call digest_integer(line(il:iu), hseed) end select + if (index(sname, 'prim') > 0) then + read(sname(7:), fmt=*) l + il = index(line, 'digest="') + 8 + iu = index(line(il:), '"') + il - 2 + call digest_integer(line(il:iu), hprim(l)) + end if + if (index(sname, 'cons') > 0) then + read(sname(7:), fmt=*) l + il = index(line, 'digest="') + 8 + iu = index(line(il:), '"') + il - 2 + call digest_integer(line(il:iu), hcons(l)) + end if end if end if go to 70 @@ -1711,37 +1761,34 @@ module io if (nd > 0) then #if NDIMS == 3 - allocate(ids(nd), array(nv,nm,nm,nm,nr), stat = status) + allocate(ids(nd), array(nv,nm,nm,nm,nr), stat=status) #else /* NDIMS == 3 */ - allocate(ids(nd), array(nv,nm,nm, 1,nr), stat = status) + allocate(ids(nd), array(nv,nm,nm, 1,nr), stat=status) #endif /* NDIMS == 3 */ if (status == 0) then - bytes = size(transfer(ids, [ 0_1 ]), kind=8) - write(fname, fmt) trim(dname), "datablock_ids", nproc, "bin" - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + write(fname,sfmt) trim(dname), "datablock_ids", n, "bin" + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) ids close(lun) - read(hids, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(ids, 1_1, bytes), & - digest, hash_xxh64) + bytes = size(ids, kind=8) * kind(ids) + call check_digest(loc, fname, c_loc(ids), bytes, hids, dtype) - bytes = size(transfer(array(:,:,:,:,1), [ 0_1 ]), kind=8) + ubytes = size(array, kind=8) * kind(array) + pbytes = ubytes / nr do l = 1, nd call append_datablock(pdata, status) call link_blocks(block_array(ids(l))%ptr, pdata) write(fname,"(a,'datablock_prim_',i6.6,'_',i6.6,'.bin')") & - trim(dname), nproc, l - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + trim(dname), n, l + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) array(:,:,:,:,1) close(lun) - call check_hash(loc, fname, & - transfer(array(:,:,:,:,1), 1_1, bytes), hprim(l), hash_xxh64) + call check_digest(loc, fname, c_loc(array), & + pbytes, hprim(l), dtype) if (lnghosts >= nghosts) then #if NDIMS == 3 @@ -1758,13 +1805,12 @@ module io end if write(fname,"(a,'datablock_cons_',i6.6,'_',i6.6,'.bin')") & - trim(dname), nproc, l - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + trim(dname), n, l + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) array close(lun) - call check_hash(loc, fname, transfer(array, 1_1, bytes), & - hcons(l), hash_xxh64) + call check_digest(loc, fname, c_loc(array), & + ubytes, hcons(l), dtype) p = min(nregs, nr) if (lnghosts >= nghosts) then @@ -1784,53 +1830,32 @@ module io deallocate(ids, array, hprim, hcons, stat=status) if (status /= 0) & - call print_message(loc, "Could not release memory!") + call print_message(loc, "Could not release space of datablocks!") else - call print_message(loc, "Could not allocate memory!") + call print_message(loc, "Could not allocate space for datablocks!") end if end if - allocate(seeds(4,lnseeds), stat = status) - - if (status == 0) then - - bytes = size(transfer(seeds, [ 0_1 ]), kind=8) - write(fname, fmt) trim(dname), "random_seeds", nproc, "bin" - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') - read(lun) seeds - close(lun) - read(hseed, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(seeds, 1_1, bytes), & - digest, hash_xxh64) - call set_seeds(lnseeds, seeds(:,:), .false.) - - if (allocated(seeds)) deallocate(seeds) - - end if - end do ! n = nl, nu -! restore seeds -! - allocate(seeds(4,lnseeds), stat = status) + allocate(seeds(4,lnseeds), stat=status) if (status == 0) then - bytes = size(transfer(seeds, [ 0_1 ]), kind=8) - write(fname, fmt) trim(dname), "random_seeds", nproc, "bin" - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream') + write(fname,sfmt) trim(dname), "random_seeds", nproc, "bin" + open(newunit=lun, file=fname, form='unformatted', access='stream') read(lun) seeds close(lun) - read(hseed, fmt = "(1z16)") digest - call check_hash(loc, fname, transfer(seeds, 1_1, bytes), & - digest, hash_xxh64) + bytes = size(seeds, kind=8) * kind(seeds) + call check_digest(loc, fname, c_loc(seeds), bytes, hseed, dtype) call set_seeds(lnseeds, seeds(:,:), .false.) - if (allocated(seeds)) deallocate(seeds) - - end if ! allocation + deallocate(seeds, stat=status) + if (status /= 0) & + call print_message(loc, "Could not release space of seeds!") + else + call print_message(loc, "Could not allocate space for seeds!") + end if end if ! nprocs >= lnprocs @@ -1863,25 +1888,26 @@ module io ! subroutine store_restart_snapshot_xml(problem, nrun, status) - use blocks , only : block_meta, block_data, list_meta, list_data - use blocks , only : get_mblocks, get_dblocks, get_nleafs - use blocks , only : get_last_id - use blocks , only : ns => nsides, nc => nchildren, nr => nregs - use coordinates, only : nn => bcells, ncells, nghosts, minlev, maxlev - use coordinates, only : xmin, xmax, ymin, ymax + use blocks , only : block_meta, block_data, list_meta, list_data + use blocks , only : get_mblocks, get_dblocks, get_nleafs + use blocks , only : get_last_id + use blocks , only : ns => nsides, nc => nchildren, nr => nregs + use coordinates , only : nn => bcells, ncells, nghosts, minlev, maxlev + use coordinates , only : xmin, xmax, ymin, ymax #if NDIMS == 3 - use coordinates, only : zmin, zmax + use coordinates , only : zmin, zmax #endif /* NDIMS == 3 */ - use coordinates, only : bdims => domain_base_dims - use equations , only : eqsys, eos, nv, cmax - use evolution , only : step, time, dt, dth, dte, cfl, glm_alpha, errs - use evolution , only : atol, rtol, mrej, niterations, nrejections - use forcing , only : nmodes, fcoefs, einj - use hash , only : hash_xxh64, hash_string - use helpers , only : print_message - use mpitools , only : nprocs, nproc - use parameters , only : get_parameter_file - use random , only : gentype, nseeds, get_seeds + use coordinates , only : bdims => domain_base_dims + use equations , only : eqsys, eos, nv, cmax + use evolution , only : step, time, dt, dth, dte, cfl, glm_alpha, errs + use evolution , only : atol, rtol, mrej, niterations, nrejections + use forcing , only : nmodes, fcoefs, einj + use hash , only : hash_info + use helpers , only : print_message + use iso_c_binding, only : c_loc + use mpitools , only : nprocs, nproc + use parameters , only : get_parameter_file + use random , only : gentype, nseeds, get_seeds implicit none @@ -1890,9 +1916,10 @@ module io integer , intent(out) :: status logical :: test - character(len=64) :: dname, fname, aname, hname + character(len=64) :: dname, fname, aname integer(kind=8) :: digest, bytes integer(kind=4) :: lun = 103 + integer :: dtype, dlen integer :: nd, nl, nm, nx, i, j, l, n, p #if NDIMS == 3 integer :: k @@ -1901,29 +1928,29 @@ module io type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata - integer(kind=4), dimension(:) , allocatable :: ids - integer(kind=4), dimension(:,:) , allocatable :: fields - integer(kind=4), dimension(:,:) , allocatable :: children + integer(kind=4), dimension(:) , allocatable, target :: ids + integer(kind=4), dimension(:,:) , allocatable, target :: fields + integer(kind=4), dimension(:,:) , allocatable, target :: children #if NDIMS == 2 - integer(kind=4), dimension(:,:,:,:) , allocatable :: edges - integer(kind=4), dimension(:,:,:) , allocatable :: corners + integer(kind=4), dimension(:,:,:,:) , allocatable, target :: edges + integer(kind=4), dimension(:,:,:) , allocatable, target :: corners #endif /* NDIMS == 2 */ #if NDIMS == 3 - integer(kind=4), dimension(:,:,:,:,:) , allocatable :: faces - integer(kind=4), dimension(:,:,:,:,:) , allocatable :: edges - integer(kind=4), dimension(:,:,:,:) , allocatable :: corners + integer(kind=4), dimension(:,:,:,:,:), allocatable, target :: faces + integer(kind=4), dimension(:,:,:,:,:), allocatable, target :: edges + integer(kind=4), dimension(:,:,:,:) , allocatable, target :: corners #endif /* NDIMS == 3 */ - integer(kind=8), dimension(:,:) , allocatable :: seeds - real(kind=8) , dimension(:,:,:) , allocatable :: bounds + integer(kind=8), dimension(:,:) , allocatable, target :: seeds + real(kind=8) , dimension(:,:,:) , allocatable, target :: bounds character(len=*), parameter :: loc = "IO::store_restart_snapshot_xml()" - character(len=*), parameter :: fmt = "(a,a,'_',i6.6,'.',a)" + character(len=*), parameter :: sfmt = "(a,a,'_',i6.6,'.',a)" !------------------------------------------------------------------------------- ! status = 0 - hname = hash_string(hash_xxh64) + call hash_info("xxh64", dtype, dlen) write(dname, "('restart-',i5.5)") nrun @@ -2119,46 +2146,52 @@ module io end do write(fname,"(a,'.bin')") "metablock_fields" - call write_binary_xml(dname, fname, transfer(fields, [ 0_1 ]), & - hash_xxh64, bytes, digest) - call write_attribute_xml(lun, "fields", fname, hname, bytes, digest) + bytes = size(fields, kind=8) * kind(fields) + call write_binary_xml(dname, fname, c_loc(fields), bytes, dtype, digest) + call write_attribute_xml(lun, "fields", fname, 'int32', & + shape(fields), bytes, dtype, digest) write(fname,"(a,'.bin')") "metablock_children" - call write_binary_xml(dname, fname, transfer(children, [ 0_1 ]), & - hash_xxh64, bytes, digest) - call write_attribute_xml(lun, "children", fname, hname, bytes, digest) + bytes = size(children, kind=8) * kind(children) + call write_binary_xml(dname, fname, c_loc(children), & + bytes, dtype, digest) + call write_attribute_xml(lun, "children", fname, 'int32', & + shape(children), bytes, dtype, digest) #if NDIMS == 3 write(fname,"(a,'.bin')") "metablock_faces" - call write_binary_xml(dname, fname, transfer(faces, [ 0_1 ]), & - hash_xxh64, bytes, digest) - call write_attribute_xml(lun, "faces", fname, hname, bytes, digest) + bytes = size(faces, kind=8) * kind(faces) + call write_binary_xml(dname, fname, c_loc(faces), bytes, dtype, digest) + call write_attribute_xml(lun, "faces", fname, 'int32', & + shape(faces), bytes, dtype, digest) #endif /* NDIMS == 3 */ write(fname,"(a,'.bin')") "metablock_edges" - call write_binary_xml(dname, fname, transfer(edges, [ 0_1 ]), & - hash_xxh64, bytes, digest) - call write_attribute_xml(lun, "edges", fname, hname, bytes, digest) + bytes = size(edges, kind=8) * kind(edges) + call write_binary_xml(dname, fname, c_loc(edges), bytes, dtype, digest) + call write_attribute_xml(lun, "edges", fname, 'int32', & + shape(edges), bytes, dtype, digest) write(fname,"(a,'.bin')") "metablock_corners" - call write_binary_xml(dname, fname, transfer(corners, [ 0_1 ]), & - hash_xxh64, bytes, digest) - call write_attribute_xml(lun, "corners", trim(fname), trim(hname), & - bytes, digest) + bytes = size(corners, kind=8) * kind(corners) + call write_binary_xml(dname, fname, c_loc(corners), & + bytes, dtype, digest) + call write_attribute_xml(lun, "corners", fname, 'int32', & + shape(corners), bytes, dtype, digest) write(fname,"(a,'.bin')") "metablock_bounds" - call write_binary_xml(trim(dname), trim(fname), & - transfer(bounds, [ 0_1 ]), hash_xxh64, bytes, digest) - call write_attribute_xml(lun, "bounds", trim(fname), trim(hname), & - bytes, digest) + bytes = size(bounds, kind=8) * kind(bounds) + call write_binary_xml(dname, fname, c_loc(bounds), bytes, dtype, digest) + call write_attribute_xml(lun, "bounds", fname, 'float64', & + shape(bounds), bytes, dtype, digest) if (nmodes > 0) then write(fname,"(a,'.bin')") "forcing_coefficients" - call write_binary_xml(trim(dname), trim(fname), & - transfer(fcoefs, [ 0_1 ]), & - hash_xxh64, bytes, digest) - call write_attribute_xml(lun, "forcing", trim(fname), trim(hname), & - bytes, digest) + bytes = size(fcoefs, kind=8) * kind(fcoefs) + call write_binary_xml(dname, fname, c_loc(fcoefs), & + bytes, dtype, digest) + call write_attribute_xml(lun, "forcing", fname, 'complex64', & + shape(fcoefs), bytes, dtype, digest) end if #if NDIMS == 3 @@ -2181,7 +2214,7 @@ module io end if - write(fname,fmt) trim(dname), "datablocks", nproc, "xml" + write(fname,sfmt) trim(dname), "datablocks", nproc, "xml" open(newunit = lun, file = fname, status = 'replace') write(lun,"(a)") "" write(lun,"(a)") '' @@ -2219,28 +2252,28 @@ module io write(aname,"('_',i6.6)") l write(fname,"('datablock_prim_',i6.6,a,'.bin')") nproc, trim(aname) - call write_binary_xml(trim(dname), trim(fname), & - transfer(pdata%q, [ 0_1 ]), & - hash_xxh64, bytes, digest) - call write_attribute_xml(lun, "prim" // trim(aname), trim(fname), & - trim(hname), bytes, digest) + bytes = size(pdata%q, kind=8) * kind(pdata%q) + call write_binary_xml(dname, fname, c_loc(pdata%q), & + bytes, dtype, digest) + call write_attribute_xml(lun, "prim" // trim(aname), fname, & + 'float64', shape(pdata%q), bytes, dtype, digest) write(fname,"('datablock_cons_',i6.6,a,'.bin')") nproc, trim(aname) - call write_binary_xml(trim(dname), trim(fname), & - transfer(pdata%uu, [ 0_1 ]), & - hash_xxh64, bytes, digest) - call write_attribute_xml(lun, "cons" // trim(aname), trim(fname), & - trim(hname), bytes, digest) + bytes = size(pdata%uu, kind=8) * kind(pdata%uu) + call write_binary_xml(dname, fname, c_loc(pdata%uu), & + bytes, dtype, digest) + call write_attribute_xml(lun, "cons" // trim(aname), fname, & + 'float64', shape(pdata%uu), bytes, dtype, digest) pdata => pdata%next end do write(fname,"(a,'_',a,'_',i6.6,'.bin')") "datablock", "ids", nproc - call write_binary_xml(trim(dname), trim(fname), & - transfer(ids, [ 0_1 ]), hash_xxh64, bytes, digest) - call write_attribute_xml(lun, "ids", trim(fname), trim(hname), & - bytes, digest) + bytes = size(ids, kind=8) * kind(ids) + call write_binary_xml(dname, fname, c_loc(ids), bytes, dtype, digest) + call write_attribute_xml(lun, "ids", fname, 'int32', & + shape(ids), bytes, dtype, digest) if (allocated(ids)) deallocate(ids) @@ -2259,10 +2292,10 @@ module io call get_seeds(seeds(:,:)) write(fname,"(a,'_',a,'_',i6.6,'.bin')") "random", "seeds", nproc - call write_binary_xml(trim(dname), trim(fname), & - transfer(seeds, [ 0_1 ]), hash_xxh64, bytes, digest) - call write_attribute_xml(lun, "seeds", trim(fname), trim(hname), & - bytes, digest) + bytes = size(seeds, kind=8) * kind(seeds) + call write_binary_xml(dname, fname, c_loc(seeds), bytes, dtype, digest) + call write_attribute_xml(lun, "seeds", fname, 'int64', & + shape(seeds), bytes, dtype, digest) if (allocated(seeds)) deallocate(seeds) @@ -2299,20 +2332,21 @@ module io ! subroutine store_snapshot_xml(problem, status) - use blocks , only : block_meta, block_data, list_meta, list_data - use blocks , only : get_dblocks, get_nleafs - use coordinates, only : nn => bcells, ncells, nghosts, minlev, maxlev - use coordinates, only : xmin, xmax, ymin, ymax + use blocks , only : block_meta, block_data, list_meta, list_data + use blocks , only : get_dblocks, get_nleafs + use coordinates , only : nn => bcells, ncells, nghosts, minlev, maxlev + use coordinates , only : xmin, xmax, ymin, ymax #if NDIMS == 3 - use coordinates, only : zmin, zmax + use coordinates , only : zmin, zmax #endif /* NDIMS == 3 */ - use coordinates, only : bdims => domain_base_dims - use equations , only : eqsys, eos, nv, pvars, adiabatic_index, csnd - use evolution , only : step, time, dt, cfl, glm_alpha - use helpers , only : print_message - use mpitools , only : nprocs, nproc - use parameters , only : get_parameter_file - use sources , only : viscosity, resistivity + use coordinates , only : bdims => domain_base_dims + use equations , only : eqsys, eos, nv, pvars, adiabatic_index, csnd + use evolution , only : step, time, dt, cfl, glm_alpha + use helpers , only : print_message + use iso_c_binding, only : c_loc + use mpitools , only : nprocs, nproc + use parameters , only : get_parameter_file + use sources , only : viscosity, resistivity implicit none @@ -2330,16 +2364,14 @@ module io type(block_meta), pointer :: pmeta type(block_data), pointer :: pdata - integer(kind=4), dimension(:) , allocatable :: ids - integer(kind=4), dimension(:,:) , allocatable :: fields - real(kind=8) , dimension(:,:,:) , allocatable :: bounds - real(kind=8) , dimension(:,:,:,:), allocatable :: array + integer(kind=4), dimension(:) , allocatable, target :: ids + integer(kind=4), dimension(:,:) , allocatable, target :: fields + real(kind=8) , dimension(:,:,:) , allocatable, target :: bounds + real(kind=8) , dimension(:,:,:,:), allocatable, target :: array -! local parameters -! character(len=*), parameter :: loc = "IO::store_snapshot_xml()" - character(len=*), parameter :: fmt = "(a,a,'_',i6.6,'.',a)" -! + character(len=*), parameter :: sfmt = "(a,a,'_',i6.6,'.',a)" + !------------------------------------------------------------------------------- ! status = 0 @@ -2467,17 +2499,19 @@ module io pmeta => pmeta%next end do - write(fname,"(a,a)") "metablock_fields", trim(binary_file_suffix) - call write_binary_xml(dname, fname, transfer(fields, [ 0_1 ]), & - hash_type, dbytes, ddigest, cbytes, cdigest) - call write_attribute_xml(lun, "fields", fname, hash_str, & - dbytes, ddigest, cbytes, cdigest) + write(fname,"(a,'.bin')") "metablock_fields" + dbytes = size(fields, kind=8) * kind(fields) + call write_binary_xml(dname, fname, c_loc(fields), & + dbytes, hash_type, ddigest, cbytes, cdigest) + call write_attribute_xml(lun, "fields", fname, 'int32', & + shape(fields), dbytes, hash_type, ddigest, cbytes, cdigest) - write(fname,"(a,a)") "metablock_bounds", trim(binary_file_suffix) - call write_binary_xml(dname, fname, transfer(bounds, [ 0_1 ]), & - hash_type, dbytes, ddigest, cbytes, cdigest) - call write_attribute_xml(lun, "bounds", fname, hash_str, & - dbytes, ddigest, cbytes, cdigest) + write(fname,"(a,'.bin')") "metablock_bounds" + dbytes = size(bounds, kind=8) * kind(bounds) + call write_binary_xml(dname, fname, c_loc(bounds), & + dbytes, hash_type, ddigest, cbytes, cdigest) + call write_attribute_xml(lun, "bounds", fname, 'float64', & + shape(bounds), dbytes, hash_type, ddigest, cbytes, cdigest) if (allocated(fields)) deallocate(fields) if (allocated(bounds)) deallocate(bounds) @@ -2496,7 +2530,7 @@ module io ! prepare and store data block info ! - write(fname,fmt) trim(dname), "datablocks", nproc, "xml" + write(fname,sfmt) trim(dname), "datablocks", nproc, "xml" open(newunit = lun, file = fname, status = 'replace') write(lun,"(a)") "" write(lun,"(a)") '' @@ -2532,12 +2566,14 @@ module io pdata => pdata%next end do ! data blocks - write(fname,"(a,'_',a,'_',i6.6,a)") "datablock", "ids", & - nproc, trim(binary_file_suffix) - call write_binary_xml(dname, fname, transfer(ids, [ 0_1 ]), & - hash_type, dbytes, ddigest, cbytes, cdigest) - call write_attribute_xml(lun, "ids", fname, hash_str, & - dbytes, ddigest, cbytes, cdigest) + write(fname,"(a,'_',a,'_',i6.6,'.bin')") "datablock", "ids", nproc + dbytes = size(ids, kind=8) * kind(ids) + call write_binary_xml(dname, fname, c_loc(ids), & + dbytes, hash_type, ddigest, cbytes, cdigest) + call write_attribute_xml(lun, "ids", fname, 'int32', shape(ids), & + dbytes, hash_type, ddigest, cbytes, cdigest) + + dbytes = size(array, kind=8) * kind(array) do p = 1, nv @@ -2550,13 +2586,12 @@ module io pdata => pdata%next end do - write(fname,"(a,'_',a,'_',i6.6,a)") "datablock", trim(pvars(p)), & - nproc, trim(binary_file_suffix) - call write_binary_xml(dname, fname, & - transfer(array(:,:,:,:), [ 0_1 ]), & - hash_type, dbytes, ddigest, cbytes, cdigest) - call write_attribute_xml(lun, pvars(p), fname, hash_str, & - dbytes, ddigest, cbytes, cdigest) + write(fname,"(a,'_',a,'_',i6.6,'.bin')") "datablock", & + trim(pvars(p)), nproc + call write_binary_xml(dname, fname, c_loc(array), & + dbytes, hash_type, ddigest, cbytes, cdigest) + call write_attribute_xml(lun, pvars(p), fname, 'float64', & + shape(array), dbytes, hash_type, ddigest, cbytes, cdigest) end do if (allocated(ids)) deallocate(ids) @@ -2720,6 +2755,8 @@ module io ! lun - the file handler to write to; ! aname - the file attribute name; ! filename - the file name; +! data_type - the data type of the input data; +! data_shape - the shape of the input data; ! digest_string - the digest type string; ! data_bytes - the file size in bytes; ! data_digest - the digest of the file content; @@ -2728,51 +2765,59 @@ module io ! !=============================================================================== ! - subroutine write_attribute_xml_file(lun, aname, filename, & - hash_string, data_bytes, data_digest, & + subroutine write_attribute_xml_file(lun, aname, filename, & + data_type, data_shape, & + data_bytes, digest_type, data_digest, & compressed_bytes, compressed_digest) + use compression, only : compression_suffix + use hash , only : hash_name, digest_string + implicit none -! input and output arguments -! integer , intent(in) :: lun - character(len=*) , intent(in) :: aname, filename, hash_string + character(len=*) , intent(in) :: aname, filename, data_type + integer , intent(in) :: digest_type + integer, dimension(:) , intent(in) :: data_shape integer(kind=8) , intent(in) :: data_bytes, data_digest - integer(kind=8) , optional, intent(in) :: compressed_bytes, compressed_digest + integer(kind=8) , optional, intent(in) :: compressed_bytes + integer(kind=8) , optional, intent(in) :: compressed_digest -! local variables -! character(len=256) :: fname - character(len=32) :: digest_string, bytes_string character(len=1024) :: string - integer :: l -! + character(len=128) :: str + integer :: n + !------------------------------------------------------------------------------- ! fname = filename string = ' 0) then - write(bytes_string,"(1i32)") compressed_bytes - string = trim(string) // & - ' compression_format="' // trim(adjustl(cformat)) // '"' // & - ' compressed_size="' // trim(adjustl(bytes_string)) // '"' + fname = trim(fname) // trim(compression_suffix) + write(str,"(1i0)") compressed_bytes + string = trim(string) // & + ' compression_format="' // trim(adjustl(cformat)) // & + '" compressed_size="' // trim(adjustl(str)) // '"' if (present(compressed_digest)) then if (compressed_digest /= 0) then - write(digest_string,"(1z0.16)") compressed_digest - string = trim(string) // ' compressed_digest="' // & - trim(adjustl(digest_string)) // '"' + call digest_string(compressed_digest, str) + string = trim(string) // & + ' compressed_digest="' // trim(adjustl(str)) // '"' end if end if - else - l = index(fname, '.bin') + 3 - fname = filename(1:l) end if end if string = trim(string) // '>' // trim(adjustl(fname)) // '' @@ -2793,76 +2838,86 @@ module io ! Arguments: ! ! path, name - the path and name where the array should be written to; -! array - the array of bytes to be written; -! hash_type - the type of digest to hash the data; +! array_ptr - the pointer to the array to store; ! array_bytes - the size of the array in bytes; -! array_digest - the digest of the input array; +! digest_type - the type of digest to hash the data; +! array_digest - the digest of the original array; ! compressed_bytes - the size of the compressed array in bytes; ! compressed_digest - the digest of the compressed array; ! !=============================================================================== ! - subroutine write_binary_xml(path, name, array, hash_type, & - array_bytes, array_digest, & - compressed_bytes, compressed_digest) + subroutine write_binary_xml(path, name, array_ptr, array_bytes, digest_type, & + array_digest, compressed_bytes, compressed_digest) - use compression, only : get_compression, compress - use hash , only : get_hash + use compression , only : get_compression, compression_bound, compress + use compression , only : compression_suffix + use hash , only : digest + use iso_c_binding, only : c_ptr, c_loc, c_f_pointer implicit none -! input and output arguments -! - character(len=*) , intent(in) :: path, name - integer(kind=1), dimension(:), intent(in) :: array - integer , intent(in) :: hash_type - integer(kind=8), optional , intent(out) :: array_bytes, compressed_bytes - integer(kind=8), optional , intent(out) :: array_digest, compressed_digest + character(len=*), intent(in) :: path, name + type(c_ptr) , intent(in) :: array_ptr + integer(kind=8) , intent(in) :: array_bytes + + integer , intent(in) :: digest_type + integer(kind=8), optional , intent(out) :: compressed_bytes + integer(kind=8), optional , intent(out) :: array_digest + integer(kind=8), optional , intent(out) :: compressed_digest -! local variables -! character(len=512) :: fname integer :: lun = 123 logical :: written - integer :: l, status + integer :: status + + integer(kind=1), dimension(:), pointer :: array + + integer(kind=1), dimension(:), allocatable, target :: buffer + type(c_ptr) :: buffer_ptr -! compression buffer -! - integer(kind=1), dimension(:), allocatable :: buffer -! !------------------------------------------------------------------------------- ! status = 0 written = .false. - array_bytes = size(array, kind=8) - if (present(array_digest)) array_digest = get_hash(array, hash_type) - write(fname,"(a,'/',a)") trim(path), trim(name) + if (present(array_digest)) & + array_digest = digest(array_ptr, array_bytes, digest_type) -! try to compress array and write it +! try to compress the array and store it if compression was successful ! if (present(compressed_bytes) .and. get_compression() > 0) then - allocate(buffer(array_bytes), stat = status) + compressed_bytes = compression_bound(array_bytes) + allocate(buffer(compressed_bytes), stat = status) + buffer_ptr = c_loc(buffer) if (status == 0) then - call compress(array, buffer, compressed_bytes) - if (compressed_bytes > 0) then - open(newunit = lun, file = fname, form = 'unformatted', & - access = 'stream', status = 'replace') + call compress(array_ptr, array_bytes, buffer_ptr, compressed_bytes) + if (compressed_bytes > 0 .and. compressed_bytes < array_bytes) then + write(fname,"(a,'/',a,a)") trim(path), trim(name), & + trim(compression_suffix) + open(newunit=lun, file=fname, form='unformatted', & + access='stream', status='replace') write(lun) buffer(1:compressed_bytes) close(lun) written = .true. if (present(compressed_digest)) & - compressed_digest = get_hash(buffer(1:compressed_bytes), hash_type) + compressed_digest = digest(buffer_ptr, & + compressed_bytes, digest_type) + else + compressed_bytes = 0 + compressed_digest = 0 end if deallocate(buffer) end if end if -! compression failed of no compression is used, so writhe the uncompressed array +! compression failed or no compression was used, so write the original array ! if (.not. written) then - l = index(fname, '.bin') + 3 - open(newunit = lun, file = fname(:l), form = 'unformatted', & - access = 'stream', status = 'replace') + call c_f_pointer(array_ptr, array, [ array_bytes ]) + + write(fname,"(a,'/',a)") trim(path), trim(name) + open(newunit=lun, file=fname, form='unformatted', & + access='stream', status='replace') write(lun) array close(lun) end if @@ -3117,7 +3172,7 @@ module io #ifdef MPI use mesh , only : redistribute_blocks #endif /* MPI */ - use mpitools, only : nprocs, npmax, nproc + use mpitools, only : nprocs, nproc implicit none @@ -3125,7 +3180,7 @@ module io character(len=255) :: fname integer(hid_t) :: file_id, grp_id - integer :: nfiles, last_id, n + integer :: nfiles, last_id, n, i, nd, nr, nl, nu, il, iu logical :: flag real(kind=8) :: deinj @@ -3176,7 +3231,7 @@ module io ! larger or equal to the number of files, and when we have less processors than ! files ! - if (nproc < nfiles) then + if (nfiles <= nprocs .and. nproc < nfiles) then write(fname, "(a,'r',i6.6,'_',i5.5,'.h5')") trim(respath), nrest, nproc inquire(file=fname, exist=flag) @@ -3216,70 +3271,85 @@ module io call print_message(loc, "Could not close '" // trim(fname) // "'!") return end if + end if ! nproc < nfiles -#ifdef MPI -! if there are more files than processes, read the remaining files by -! the last process and redistribute blocks after each processed file, -! otherwise only redistribute blocks +! if there are more files than processes, divide the files equally between +! processes ! if (nprocs < nfiles) then - do n = nprocs, nfiles - 1 - call change_blocks_process(n, npmax) + nl = 0 + nd = nfiles / nprocs + nr = mod(nfiles, nprocs) + do n = 0, nprocs - 1 + if (n < nr) then + il = n * (nd + 1) + iu = il + nd + else + il = n * nd + nr + iu = il + nd - 1 + end if + do i = il, iu + call change_blocks_process(i, n) + end do + if (n == nproc) then + nl = il + nu = iu + end if + end do - if (nproc == npmax) then - write(fname, "(a,'r',i6.6,'_',i5.5,'.h5')") trim(respath), nrest, n - inquire(file=fname, exist=flag) - if (.not. flag) then - call print_message(loc, & - "Restart snapshot '" // trim(fname) // "' not found!") - status = 1 - return - end if + do n = nl, nu - call h5fopen_f(fname, H5F_ACC_RDONLY_F, file_id, status) - if (status /= 0) then - call print_message(loc, "Could not open '" // trim(fname) // "'!") - return - end if - - call restore_datablocks_h5(file_id, status) - if (status /= 0) & - call print_message(loc, & - "Could not restore datablocks from '" // trim(fname) // "'!") - - call h5gopen_f(file_id, 'attributes', grp_id, status) - if (status /= 0) then - call print_message(loc, "Could not open group 'attributes'!") - return - end if - call restore_attribute_h5(grp_id, 'einj', & - H5T_NATIVE_DOUBLE, 1, deinj, status) - if (status /= 0) & - call print_message(loc, "Could not get the injected energy!") - einj = einj + deinj - call h5gclose_f(grp_id, status) - if (status /= 0) & - call print_message(loc, "Could not close group 'attributes'!") - - call h5fclose_f(file_id, status) - if (status /= 0) then - call print_message(loc, "Could not close '" // trim(fname) // "'!") - return - end if + write(fname, "(a,'r',i6.6,'_',i5.5,'.h5')") trim(respath), nrest, n + inquire(file=fname, exist=flag) + if (.not. flag) then + call print_message(loc, & + "Restart snapshot '" // trim(fname) // "' not found!") + status = 1 + return end if - call redistribute_blocks(status) + call h5fopen_f(fname, H5F_ACC_RDONLY_F, file_id, status) + if (status /= 0) then + call print_message(loc, "Could not open '" // trim(fname) // "'!") + return + end if + + call restore_datablocks_h5(file_id, status) + if (status /= 0) & + call print_message(loc, & + "Could not restore datablocks from '" // trim(fname) // "'!") + + call h5gopen_f(file_id, 'attributes', grp_id, status) + if (status /= 0) then + call print_message(loc, "Could not open group 'attributes'!") + return + end if + call restore_attribute_h5(grp_id, 'einj', & + H5T_NATIVE_DOUBLE, 1, deinj, status) + if (status /= 0) & + call print_message(loc, "Could not get the injected energy!") + einj = einj + deinj + call h5gclose_f(grp_id, status) + if (status /= 0) & + call print_message(loc, "Could not close group 'attributes'!") + + call h5fclose_f(file_id, status) + if (status /= 0) then + call print_message(loc, "Could not close '" // trim(fname) // "'!") + return + end if end do - else - call redistribute_blocks(status) end if -#endif /* MPI */ if (allocated(block_array)) deallocate(block_array) +#ifdef MPI + call redistribute_blocks(status) +#endif /* MPI */ + !------------------------------------------------------------------------------- ! end subroutine read_restart_snapshot_h5 @@ -3779,13 +3849,13 @@ module io allocate(array(nmodes,NDIMS), stat=status) if (status == 0) then array = real(fcoefs) - call store_dataset_h5(grp_id, 'fcoefs_real', H5T_NATIVE_DOUBLE, & - dims, array, status) + call store_dataset_h5(grp_id, 'fcoefs_real', H5T_NATIVE_DOUBLE, & + dims, array, .false., status) if (status < 0) & call print_message(loc, "Could not store real(fcoefs)!") array = aimag(fcoefs) - call store_dataset_h5(grp_id, 'fcoefs_imag', H5T_NATIVE_DOUBLE, & - dims, array, status) + call store_dataset_h5(grp_id, 'fcoefs_imag', H5T_NATIVE_DOUBLE, & + dims, array, .false., status) if (status < 0) & call print_message(loc, "Could not store imag(fcoefs)!") deallocate(array, stat=status) @@ -3813,7 +3883,7 @@ module io if (status == 0) then call get_seeds(seeds) call store_dataset_h5(grp_id, 'seeds', H5T_STD_I64LE, & - dims, seeds, status) + dims, seeds, .false., status) if (status < 0) & call print_message(loc, "Could not store seeds!") deallocate(seeds, stat=status) @@ -3926,7 +3996,7 @@ module io #if NDIMS == 3 l = rank(faces) dims(1:l) = shape(faces) - call store_dataset_h5(grp_id, 'faces', & + call read_dataset_h5(grp_id, 'faces', & H5T_NATIVE_INTEGER, dims(1:l), faces, status) #endif /* NDIMS == 3 */ l = rank(edges) @@ -4190,30 +4260,30 @@ module io l = rank(fields) dims(1:l) = shape(fields) - call store_dataset_h5(grp_id, 'fields', & - H5T_NATIVE_INTEGER, dims(1:l), fields, status) + call store_dataset_h5(grp_id, 'fields', H5T_NATIVE_INTEGER, & + dims(1:l), fields, .false., status) l = rank(children) dims(1:l) = shape(children) - call store_dataset_h5(grp_id, 'children', & - H5T_NATIVE_INTEGER, dims(1:l), children, status) + call store_dataset_h5(grp_id, 'children', H5T_NATIVE_INTEGER, & + dims(1:l), children, .false., status) #if NDIMS == 3 l = rank(faces) dims(1:l) = shape(faces) - call store_dataset_h5(grp_id, 'faces', & - H5T_NATIVE_INTEGER, dims(1:l), faces, status) + call store_dataset_h5(grp_id, 'faces', H5T_NATIVE_INTEGER, & + dims(1:l), faces, .false., status) #endif /* NDIMS == 3 */ l = rank(edges) dims(1:l) = shape(edges) - call store_dataset_h5(grp_id, 'edges', & - H5T_NATIVE_INTEGER, dims(1:l), edges, status) + call store_dataset_h5(grp_id, 'edges', H5T_NATIVE_INTEGER, & + dims(1:l), edges, .false., status) l = rank(corners) dims(1:l) = shape(corners) - call store_dataset_h5(grp_id, 'corners', & - H5T_NATIVE_INTEGER, dims(1:l), corners, status) + call store_dataset_h5(grp_id, 'corners', H5T_NATIVE_INTEGER, & + dims(1:l), corners, .false., status) l = rank(bounds) dims(1:l) = shape(bounds) - call store_dataset_h5(grp_id, 'bounds', & - H5T_NATIVE_DOUBLE, dims(1:l), bounds, status) + call store_dataset_h5(grp_id, 'bounds', H5T_NATIVE_DOUBLE, & + dims(1:l), bounds, .false., status) #if NDIMS == 3 deallocate(fields, children, bounds, faces, & @@ -4463,13 +4533,13 @@ module io cdims = shape(pdata%uu) call store_dataset_h5(blk_id, 'primitive_variables', & - H5T_NATIVE_DOUBLE, pdims, pdata%q, status) + H5T_NATIVE_DOUBLE, pdims, pdata%q, .false., status) if (status /= 0) & call print_message(loc, & "Could not store the primitive variables in " // & trim(blk_name) // "!") call store_dataset_h5(blk_id, 'conservative_variables', & - H5T_NATIVE_DOUBLE, cdims, pdata%uu, status) + H5T_NATIVE_DOUBLE, cdims, pdata%uu, .false., status) if (status /= 0) & call print_message(loc, & "Could not store the conservative variables in " // & @@ -4555,10 +4625,13 @@ module io end if am(1) = toplev - call store_dataset_h5(grp_id, 'dx', H5T_NATIVE_DOUBLE, am, adx, status) - call store_dataset_h5(grp_id, 'dy', H5T_NATIVE_DOUBLE, am, ady, status) + call store_dataset_h5(grp_id, 'dx', H5T_NATIVE_DOUBLE, & + am, adx, .false., status) + call store_dataset_h5(grp_id, 'dy', H5T_NATIVE_DOUBLE, & + am, ady, .false., status) #if NDIMS == 3 - call store_dataset_h5(grp_id, 'dz', H5T_NATIVE_DOUBLE, am, adz, status) + call store_dataset_h5(grp_id, 'dz', H5T_NATIVE_DOUBLE, & + am, adz, .false., status) #endif /* NDIMS == 3 */ if (get_dblocks() > 0) then @@ -4597,14 +4670,14 @@ module io pdata => pdata%next end do - call store_dataset_h5(grp_id, 'ids', & - H5T_NATIVE_INTEGER, im, ids, status) - call store_dataset_h5(grp_id, 'levels', & - H5T_NATIVE_INTEGER, im, levels, status) - call store_dataset_h5(grp_id, 'coords', & - H5T_NATIVE_INTEGER, cm, coords, status) - call store_dataset_h5(grp_id, 'bounds', & - H5T_NATIVE_DOUBLE, bm, bounds, status) + call store_dataset_h5(grp_id, 'ids', H5T_NATIVE_INTEGER, & + im, ids, .false., status) + call store_dataset_h5(grp_id, 'levels', H5T_NATIVE_INTEGER, & + im, levels, .false., status) + call store_dataset_h5(grp_id, 'coords', H5T_NATIVE_INTEGER, & + cm, coords, .false., status) + call store_dataset_h5(grp_id, 'bounds', H5T_NATIVE_DOUBLE, & + bm, bounds, .false., status) deallocate(ids, levels, coords, bounds, stat=status) if (status > 0) & @@ -4691,7 +4764,7 @@ module io end do call store_dataset_h5(grp_id, trim(pvars(p)), & - H5T_NATIVE_DOUBLE, dims, array, status) + H5T_NATIVE_DOUBLE, dims, array, .true., status) end do deallocate(array, stat=status) @@ -5260,16 +5333,18 @@ module io ! ! Arguments: ! -! loc_id - the location in which the dataset is stored; -! name - the dataset name; -! type_id - the dataset type; -! dims - the dataset dimensions; -! buffer - the dataset buffer to store; -! status - the subroutine call status; +! loc_id - the location in which the dataset is stored; +! name - the dataset name; +! type_id - the dataset type; +! dims - the dataset dimensions; +! buffer - the dataset buffer to store; +! compress - the logical flag inficating is dataset should be compressed; +! status - the subroutine call status; ! !=============================================================================== ! - subroutine store_dataset_h5(loc_id, name, type_id, dims, buffer, status) + subroutine store_dataset_h5(loc_id, name, type_id, dims, & + buffer, compress, status) use helpers , only : print_message use iso_c_binding, only : c_loc @@ -5280,16 +5355,21 @@ module io character(len=*) , intent(in) :: name integer(hsize_t), dimension(:), intent(in) :: dims type(*), target, dimension(..), intent(in) :: buffer + logical , intent(in) :: compress integer , intent(out) :: status integer :: rank integer(hid_t) :: space_id, dset_id + integer(hsize_t), dimension(size(dims)) :: cdims + character(len=*), parameter :: loc = 'IO::store_dataset_h5()' !------------------------------------------------------------------------------- ! - rank = size(dims) + rank = size(dims) + cdims = dims + if (compress .and. hcformat .eq. H5Z_ZFP) cdims(rank) = 1 call h5screate_simple_f(rank, dims, space_id, status) if (status /= 0) then @@ -5297,14 +5377,18 @@ module io "Could not create the dataspace for dataset '" // trim(name) // "'!") return end if - call h5pset_chunk_f(prp_id, rank, dims, status) + call h5pset_chunk_f(prp_id, rank, cdims, status) if (status /= 0) then call print_message(loc, & "Could not set the chunk size for dataset '" // trim(name) // "'!") go to 1000 end if - call h5dcreate_f(loc_id, name, type_id, space_id, dset_id, status, prp_id) + if (compress) then + call h5dcreate_f(loc_id, name, type_id, space_id, dset_id, status, prp_id) + else + call h5dcreate_f(loc_id, name, type_id, space_id, dset_id, status) + end if if (status /= 0) then call print_message(loc, & "Could not create dataset '" // trim(name) // "'!")