diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 092c9431a..a57387f8a 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -10,6 +10,7 @@ stages: script: - echo "Setup environment" - apt-get update + - apt-get install -y git - apt-get install -y cmake - apt-get install -y gcc - apt-get install -y infiniband-diags ibverbs-utils diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 000000000..e162648e3 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "lib/ADIOS2"] + path = lib/ADIOS2 + url = git@github.com:ornladios/ADIOS2.git diff --git a/CMakeLists.txt b/CMakeLists.txt index 8ee17c0e6..f631c44b3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,17 @@ cmake_minimum_required(VERSION 3.10) -project(x3d2 LANGUAGES Fortran) +project(x3d2 LANGUAGES Fortran CXX C) enable_testing() +# TODO switch adios between system install and this installed version + +include(FetchContent) + +FetchContent_Declare( + Adios2 + GIT_REPOSITORY https://github.com/ornladios/ADIOS2.git + GIT_TAG v2.10.0-rc1 +) +FetchContent_MakeAvailable(Adios2) + add_subdirectory(src) add_subdirectory(tests) diff --git a/lib/ADIOS2 b/lib/ADIOS2 new file mode 160000 index 000000000..aab3bf23e --- /dev/null +++ b/lib/ADIOS2 @@ -0,0 +1 @@ +Subproject commit aab3bf23ea19b09b17e6171e11ba0caf67ddef30 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8f3527576..7e2c8ffcb 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,5 @@ set(SRC + adios_io.f90 allocator.f90 backend.f90 common.f90 @@ -67,3 +68,5 @@ find_package(MPI REQUIRED) target_link_libraries(x3d2 PRIVATE MPI::MPI_Fortran) target_link_libraries(xcompact PRIVATE MPI::MPI_Fortran) +target_link_libraries(x3d2 PRIVATE adios2::fortran_mpi) +target_link_libraries(xcompact PRIVATE adios2::fortran_mpi) diff --git a/src/adios_io.f90 b/src/adios_io.f90 new file mode 100644 index 000000000..4a8e8724f --- /dev/null +++ b/src/adios_io.f90 @@ -0,0 +1,330 @@ +module m_adios_io + use adios2 + use mpi + use m_allocator, only: field_t + + use, intrinsic :: iso_fortran_env, only : stdin=>input_unit, & + stdout=>output_unit, & + stderr=>error_unit + + implicit none + + type :: adios_io_t + type(adios2_adios) :: adios_ctx + type(adios2_io) :: aio + integer :: irank, nproc + integer :: io_mode_write = adios2_mode_write + integer :: io_mode_read = adios2_mode_read + contains + procedure :: init + procedure :: deinit + procedure :: handle_fatal_error + procedure :: read_field + procedure :: write_field + procedure :: write_real8 + procedure :: write_integer8 + procedure :: open_file + procedure :: close_file + end type + + type :: adios_file_t + type(adios2_engine) :: engine + end type + +contains + subroutine init(self, comm_in) + class(adios_io_t), intent(inout) :: self + integer, intent(in), optional :: comm_in + + integer :: comm + integer :: ierr + logical :: is_mpi_initialised + + if(present(comm_in)) then + comm = comm_in + else + comm = MPI_COMM_WORLD + endif + + call mpi_initialized(is_mpi_initialised, ierr) + + if(.not. is_mpi_initialised) then + call self%handle_fatal_error("IO: ADIOS must be initialised after MPI!", 0) + endif + + call MPI_Comm_rank(comm, self%irank, ierr) + + call adios2_init(self%adios_ctx, comm, ierr) + if (.not.self%adios_ctx%valid) then + call self%handle_fatal_error("IO: Cannot initialise ADIOS context.", ierr) + endif + + call adios2_declare_io (self%aio, self%adios_ctx, "main_io", ierr) + if (.not.self%aio%valid) then + call self%handle_fatal_error("IO: Cannot create ADIOS2 IO", ierr) + endif + + + end subroutine + + subroutine deinit(self) + class(adios_io_t), intent(inout) :: self + integer :: ierr + + if (self%adios_ctx%valid) then + call adios2_finalize(self%adios_ctx, ierr) + endif + end subroutine + + pure function calc_local_strided_offset(offset, stride) result(local_offset) + ! Returns the first index in the local + ! subregion of a strided, distributed array. + integer(8), intent(in) :: offset ! offset of subregion owned by this rank + integer(8), intent(in) :: stride ! number of grid points between adjacent output indices + integer(8) :: local_offset ! return value: offset into local region respecting global stride + + local_offset = mod(stride - mod(offset, stride), stride) + 1 + end function + + function open_file(self, fpath, mode) result(file) + class(adios_io_t), intent(inout) :: self + character(*), intent(in) :: fpath !! Path to ouptut file + integer :: mode !! Read/write mode. Should be adios2_mode_write or adios2_mode_read + type(adios_file_t) :: file + + integer :: ierr + + call adios2_open(file%engine, self%aio, fpath, mode, ierr) + if (.not.file%engine%valid) then + call self%handle_fatal_error("IO: Cannot create ADIOS2 engine", ierr) + endif + + if(mode == adios2_mode_write) then + call adios2_begin_step(file%engine, adios2_step_mode_append, ierr) + else if(mode == adios2_mode_read) then + call adios2_begin_step(file%engine, adios2_step_mode_read, ierr) + else + error stop "Unsupported mode" + endif + + if (ierr /= 0) then + call self%handle_fatal_error("IO: Cannot begin step", ierr) + endif + + end function + + subroutine close_file(self, file) + class(adios_io_t), intent(inout) :: self + class(adios_file_t), intent(inout) :: file + + integer :: ierr + ! Might be abusing ADIOS here... We tend to redefined variables every time we open a file + call adios2_remove_all_variables(self%aio, ierr) + call adios2_end_step(file%engine, ierr) + + if (file%engine%valid) then + call adios2_close(file%engine, ierr) + endif + end subroutine + + subroutine write_real8(self, var, file, varname) + class(adios_io_t), intent(inout) :: self + real(8), intent(in) :: var + type(adios_file_t), intent(inout) :: file + character(*), intent(in) :: varname !! Name of variable in output file + + type(adios2_variable) :: adios_var + integer :: vartype + + integer :: ierr + + vartype = adios2_type_real8 + + call adios2_define_variable(adios_var, self%aio, varname, vartype, ierr) + call adios2_put(file%engine, adios_var, var, ierr) + end subroutine + + subroutine write_integer8(self, var, file, varname) + class(adios_io_t), intent(inout) :: self + integer(8), intent(in) :: var + type(adios_file_t), intent(inout) :: file + character(*), intent(in) :: varname !! Name of variable in output file + + type(adios2_variable) :: adios_var + integer :: vartype + + integer :: ierr + + vartype = adios2_type_integer8 + + call adios2_define_variable(adios_var, self%aio, varname, vartype, ierr) + call adios2_put(file%engine, adios_var, var, ierr) + end subroutine + + subroutine write_field(self, in_arr, file, varname, icount, ishape, istart, convert_to_sp_in, istride_in) + class(adios_io_t), intent(inout) :: self + class(field_t), pointer, intent(in) :: in_arr !! Field to be outputted + ! TODO should this be a field or a fortran array? + type(adios_file_t), intent(inout) :: file + character(*), intent(in) :: varname !! Name of variable in output file + integer(8), dimension(3), intent(in) :: icount !! Local size of in_arr + integer(8), dimension(3), intent(in) :: ishape !! Global size of in_arr + integer(8), dimension(3), intent(in) :: istart !! Local offset of in_arr + + !> if .true. input array will be converted to single precision before being + !> outputted. + !> Defaults to .false. + logical, intent(in), optional :: convert_to_sp_in + logical :: convert_to_sp = .false. + logical :: should_stride = .false. + !> If set, will coarsen output in each direction by only dumping every + !> `istride` gridpoints. + !> Defaults to (1,1,1). + integer(8), dimension(3), intent(in), optional :: istride_in + integer(8), dimension(3) :: istride = [1,1,1] + + real(4), allocatable :: data_sp(:,:,:) + real(4), allocatable :: data_strided_sp(:,:,:) + real(8), allocatable :: data_strided_dp(:,:,:) + integer(8), dimension(3) :: istart_strided + integer(8), dimension(3) :: icount_strided + + type(adios2_variable) :: adios_var + integer :: vartype + integer :: ierr + + integer :: i + + ! Set our optional inputs + if(present(convert_to_sp_in)) then + convert_to_sp = convert_to_sp_in + else + convert_to_sp = .false. + endif + if(present(istride_in)) istride = istride_in + + if(istride(1) < 1 .or. istride(2) < 1 .or. istride(3) < 1) then + call self%handle_fatal_error("Output stride < 1. Cannot continue.", 0) + endif + + if(convert_to_sp) then + print*, "converting to single precision!" + allocate(data_sp(icount(1), icount(2), icount(3))) + data_sp(:,:,:) = real(in_arr%data(:,:,:)) + vartype = adios2_type_real + else + vartype = adios2_type_dp + endif + + should_stride = any(istride /= [1,1,1]) + + if(should_stride) then + print*, "striding!" + do i = 1, 3 + istart_strided(i) = calc_local_strided_offset(istart(i), istride(i)) + end do + icount_strided = icount/istride + if(convert_to_sp) then + allocate(data_strided_sp(& + icount_strided(1),& + icount_strided(2),& + icount_strided(3)& + )) + data_strided_sp(:,:,:) = data_sp(& + istart_strided(1)::istride(1),& + istart_strided(2)::istride(2),& + istart_strided(3)::istride(3)& + ) + else + allocate(data_strided_dp(& + icount_strided(1),& + icount_strided(2),& + icount_strided(3)& + )) + data_strided_dp(:,:,:) = in_arr%data(& + istart_strided(1)::istride(1),& + istart_strided(2)::istride(2),& + istart_strided(3)::istride(3)& + ) + endif + endif + + if(should_stride) then + call adios2_define_variable(adios_var, self%aio, varname, vartype, & + 3, ishape/istride, istart/istride, icount_strided, adios2_constant_dims, ierr) + else + call adios2_define_variable(adios_var, self%aio, varname, vartype, & + 3, ishape, istart, icount, adios2_constant_dims, ierr) + endif + + if(should_stride) then + if(convert_to_sp) then + call adios2_put(file%engine, adios_var, data_strided_sp, ierr) + else + call adios2_put(file%engine, adios_var, data_strided_dp, ierr) + endif + else + if(convert_to_sp) then + call adios2_put(file%engine, adios_var, data_sp, ierr) + else + call adios2_put(file%engine, adios_var, in_arr%data, ierr) + endif + endif + + if(allocated(data_sp)) then + deallocate(data_sp) + endif + + end subroutine + + subroutine read_field(self, out_arr, file, varname, icount, ishape, istart, idump_in) + class(adios_io_t), intent(inout) :: self + class(field_t), pointer, intent(in) :: out_arr !! Field to be read from file + type(adios_file_t), intent(inout) :: file !! File already opened for reading + character(*), intent(in) :: varname !! Name of variable in input file + integer(kind=8), dimension(3), intent(in) :: icount !! Local size of in_arr + integer(kind=8), dimension(3), intent(in) :: ishape !! Global size of in_arr + integer(kind=8), dimension(3), intent(in) :: istart !! Local offset of in_arr + !> When multiple timesteps have been dumped into one variable, chose which + !> dump to read. + integer(kind=8), intent(in), optional :: idump_in + integer(kind=8) :: idump = 0 + + integer(8), parameter :: n_steps = 1 !! This version reads one dump at a time + + type(adios2_variable) :: adios_var + integer :: ierr + + if (present(idump_in)) idump = idump_in + + call adios2_inquire_variable(adios_var, self%aio, varname, ierr) + if (.not.adios_var%valid) then + call self%handle_fatal_error("Cannot fetch ADIOS2 variable", ierr) + endif + + call adios2_set_step_selection(adios_var, idump, n_steps, ierr) + if (ierr /= 0) then + call self%handle_fatal_error("IO: Cannot set step selection", ierr) + endif + + call adios2_set_selection(adios_var, 3, istart, icount, ierr) + if (ierr /= 0) then + call self%handle_fatal_error("IO: Cannot set selection", ierr) + endif + + call adios2_get(file%engine, adios_var, out_arr%data, ierr) + if (ierr /= 0) then + call self%handle_fatal_error("IO: Cannot fetch data from file", ierr) + endif + end subroutine + + subroutine handle_fatal_error(self, msg, ierr) + class(adios_io_t), intent(in) :: self + integer, intent(in) :: ierr + character(*), intent(in) :: msg + + write(stderr, *) "ADIOS2 Error code:", ierr + write(stderr, *) msg + error stop -1 + end subroutine +end module diff --git a/src/solver.f90 b/src/solver.f90 index 4f9547b3f..5bea1c1fb 100644 --- a/src/solver.f90 +++ b/src/solver.f90 @@ -12,6 +12,7 @@ module m_solver use m_time_integrator, only: time_intg_t use m_vector_calculus, only: vector_calculus_t use m_mesh, only: mesh_t + use m_adios_io, only: adios_io_t, adios_file_t implicit none @@ -46,11 +47,14 @@ module m_solver real(dp) :: dt, nu integer :: n_iters, n_output integer :: ngrid + integer :: steps_between_checkpoints ! checkpoints can be used as restarts + integer :: steps_between_snapshots - class(field_t), pointer :: u, v, w + class(field_t), pointer :: u, v, w, pressure class(base_backend_t), pointer :: backend class(mesh_t), pointer :: mesh + type(adios_io_t) :: io type(time_intg_t) :: time_integrator type(allocator_t), pointer :: host_allocator type(dirps_t), pointer :: xdirps, ydirps, zdirps @@ -63,6 +67,7 @@ module m_solver procedure :: curl procedure :: output procedure :: run + procedure :: write_variable end type solver_t abstract interface @@ -122,6 +127,7 @@ function init(backend, mesh, host_allocator) result(solver) solver%u => solver%backend%allocator%get_block(DIR_X, VERT) solver%v => solver%backend%allocator%get_block(DIR_X, VERT) solver%w => solver%backend%allocator%get_block(DIR_X, VERT) + solver%pressure => solver%backend%allocator%get_block(DIR_Z, VERT) ! set defaults poisson_solver_type = 'FFT' @@ -145,6 +151,8 @@ function init(backend, mesh, host_allocator) result(solver) print *, time_intg//' time integrator instantiated' end if + solver%steps_between_checkpoints = 1 + solver%steps_between_snapshots = 1 solver%dt = dt solver%backend%nu = 1._dp/Re solver%n_iters = n_iters @@ -204,6 +212,8 @@ function init(backend, mesh, host_allocator) result(solver) error stop 'poisson_solver_type is not valid. Use "FFT" or "CG".' end select + call solver%io%init(MPI_COMM_WORLD) + end function init subroutine allocate_tdsops(dirps, backend, der1st_scheme, der2nd_scheme, & @@ -397,6 +407,24 @@ subroutine poisson_fft(self, pressure, div_u) end subroutine poisson_fft + subroutine write_variable(self, file, varname, varfield, icount, ishape, istart) + class(solver_t), intent(in) :: self + type(adios_file_t) :: file + character(*), intent(in) :: varname !! Name of variable in output file + class(field_t), intent(in) :: varfield !! Field to be written out + integer(8), dimension(3), intent(in) :: icount !! Local size of in_arr + integer(8), dimension(3), intent(in) :: ishape !! Global size of in_arr + integer(8), dimension(3), intent(in) :: istart !! Local offset of in_arr + + class(field_t), pointer :: out + + out => self%host_allocator%get_block(DIR_C) + ! Backend should handle moving from device if necessary + call self%backend%get_field_data(out%data, varfield) + call self%io%write_field(out, file, varname, icount, ishape, istart) + call self%host_allocator%release_block(out) + end subroutine write_variable + subroutine poisson_cg(self, pressure, div_u) implicit none @@ -461,17 +489,24 @@ subroutine run(self) class(solver_t), intent(inout) :: self - class(field_t), pointer :: du, dv, dw, div_u, pressure, dpdx, dpdy, dpdz + class(field_t), pointer :: du, dv, dw, div_u, dpdx, dpdy, dpdz class(field_t), pointer :: u_out, v_out, w_out + integer(8), dimension(3) :: icount !! Local size of output array + integer(8), dimension(3) :: ishape !! Global size of output array + integer(8), dimension(3) :: istart !! Local offset of output array + + type(adios_file_t) :: checkpoint_file real(dp) :: t integer :: i, j + integer :: checkpoint_at_i = 1 if (self%mesh%par%is_root()) print *, 'initial conditions' t = 0._dp call self%output(t) if (self%mesh%par%is_root()) print *, 'start run' + checkpoint_at_i = checkpoint_at_i + self%steps_between_checkpoints do i = 1, self%n_iters do j = 1, self%time_integrator%nstage @@ -494,9 +529,7 @@ subroutine run(self) call self%divergence_v2p(div_u, self%u, self%v, self%w) - pressure => self%backend%allocator%get_block(DIR_Z, CELL) - - call self%poisson(pressure, div_u) + call self%poisson(self%pressure, div_u) call self%backend%allocator%release_block(div_u) @@ -504,9 +537,7 @@ subroutine run(self) dpdy => self%backend%allocator%get_block(DIR_X) dpdz => self%backend%allocator%get_block(DIR_X) - call self%gradient_p2v(dpdx, dpdy, dpdz, pressure) - - call self%backend%allocator%release_block(pressure) + call self%gradient_p2v(dpdx, dpdy, dpdz, self%pressure) ! velocity correction call self%backend%vecadd(-1._dp, dpdx, 1._dp, self%u) @@ -518,6 +549,34 @@ subroutine run(self) call self%backend%allocator%release_block(dpdz) end do + if(i == checkpoint_at_i) then + checkpoint_file = self%io%open_file("test_output.bp", self%io%io_mode_write) + + icount = self%mesh%get_dims(self%pressure%data_loc) + ishape = self%mesh%get_global_dims(self%pressure%data_loc) + istart = self%mesh%par%n_offset ! Is this offset the same for pressure and velocity? + + call self%time_integrator%write_checkpoint(checkpoint_file, self%io) + + call self%write_variable(checkpoint_file, "pressure", self%pressure, & + icount, ishape, istart) + + icount = self%mesh%get_dims(self%u%data_loc) + ishape = self%mesh%get_global_dims(self%u%data_loc) + istart = self%mesh%par%n_offset + + call self%write_variable(checkpoint_file, "u", self%u, icount, ishape, istart) + call self%write_variable(checkpoint_file, "v", self%v, icount, ishape, istart) + call self%write_variable(checkpoint_file, "w", self%w, icount, ishape, istart) + + call self%io%write_real8(t, checkpoint_file, "time") + call self%io%write_integer8(1_8, checkpoint_file, "is_valid") + + call self%io%close_file(checkpoint_file) + + checkpoint_at_i = checkpoint_at_i + self%steps_between_checkpoints + end if + if (mod(i, self%n_output) == 0) then t = i*self%dt call self%output(t) diff --git a/src/time_integrator.f90 b/src/time_integrator.f90 index 4dd326af6..ddb85c7f8 100644 --- a/src/time_integrator.f90 +++ b/src/time_integrator.f90 @@ -2,6 +2,7 @@ module m_time_integrator use m_allocator, only: allocator_t, field_t, flist_t use m_base_backend, only: base_backend_t use m_common, only: dp, DIR_X + use m_adios_io, only: adios_io_t, adios_file_t implicit none @@ -24,6 +25,7 @@ module m_time_integrator procedure :: step procedure :: runge_kutta procedure :: adams_bashforth + procedure :: write_checkpoint end type time_intg_t interface time_intg_t @@ -166,6 +168,23 @@ function init(backend, allocator, method, nvars) end function init + subroutine write_checkpoint(self, file, io) + class(time_intg_t), intent(inout) :: self + type(adios_file_t), intent(inout) :: file !! File already opened for writing + class(adios_io_t), intent(inout) :: io + + call io%write_real8(3.14_8, file, "pi") ! TODO save method/order data + + ! if adams_bashforth: + ! io%write(..., self%istep) ! required to ensure startup doesn't happen again + ! for var in olds: + ! for step in olds[var]: ! index 1 is the most recent + ! data = olds[var][step] + ! call io%write(filename, "old_{var}_{step}", data) + ! else if runge_kutta: + ! do nothing (TODO) + end subroutine + subroutine step(self, u, v, w, du, dv, dw, dt) implicit none diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5795f0076..064f290d4 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -2,6 +2,8 @@ set(CMAKE_CTEST_ARGUMENTS "--output-on-failure" CACHE STRING "comma separated list of arguments to pass to ctest") set(CMAKE_CTEST_NPROCS "4" CACHE STRING "number of process used in tests") +find_package(MPI REQUIRED) + function(define_test testfile np backend) get_filename_component(test_name ${testfile} NAME_WE) string(CONCAT test_name ${test_name} _ ${backend} _ ${np}) @@ -20,6 +22,8 @@ function(define_test testfile np backend) target_link_libraries(${test_name} PRIVATE OpenMP::OpenMP_Fortran) endif() target_link_libraries(${test_name} PRIVATE x3d2) + target_link_libraries(${test_name} PRIVATE MPI::MPI_Fortran) + target_link_libraries(${test_name} PRIVATE adios2::fortran_mpi) add_test(NAME ${test_name} COMMAND sh -c "mpirun --oversubscribe -np ${np} ${test_name}") @@ -32,6 +36,7 @@ define_test(test_reordering.f90 2 omp) define_test(test_reordering.f90 4 omp) define_test(test_reorder_map.f90 1 omp) define_test(test_time_integrator.f90 1 omp) +define_test(test_hostside_io.f90 ${CMAKE_CTEST_NPROCS} omp) define_test(omp/test_omp_tridiag.f90 ${CMAKE_CTEST_NPROCS} omp) define_test(omp/test_omp_transeq.f90 ${CMAKE_CTEST_NPROCS} omp) define_test(omp/test_omp_dist_transeq.f90 ${CMAKE_CTEST_NPROCS} omp) diff --git a/tests/test_hostside_io.f90 b/tests/test_hostside_io.f90 new file mode 100644 index 000000000..0a5f91233 --- /dev/null +++ b/tests/test_hostside_io.f90 @@ -0,0 +1,557 @@ +program test_hostside_io + use iso_fortran_env, only: stderr => error_unit + use adios2 + use mpi + + use m_adios_io + use m_omp_backend, only: omp_backend_t + use m_mesh, only: mesh_t + + implicit none + + !> MPI vars + integer :: irank, nproc + !> Error code, used for MPI and ADIOS + integer :: ierr + + call run_tests() + +contains + + subroutine run_tests() + logical :: allpass = .true. !! flag indicating all tests pass + + call MPI_Init(ierr) + call MPI_Comm_rank(MPI_COMM_WORLD, irank, ierr) + call MPI_Comm_size(MPI_COMM_WORLD, nproc, ierr) + + if (irank == 0) print*, 'Run with', nproc, 'ranks' + + if(.not. test_basic_io()) then + write(stderr, '(a)') 'Basic IO test failed.' + allpass = .false. + endif + + if(.not. test_single_precision_io()) then + write(stderr, '(a)') 'Single precision IO test failed.' + allpass = .false. + endif + + if(.not. test_strided_io()) then + write(stderr, '(a)') 'Strided IO test failed.' + allpass = .false. + endif + + if(.not. test_strided_single_precision_io()) then + write(stderr, '(a)') 'Strided + SP IO test failed.' + allpass = .false. + endif + + if (allpass) then + if (irank == 0) write(stderr, '(a)') 'ALL TESTS PASSED SUCCESSFULLY.' + end if + + call MPI_Finalize(ierr) + + if(.not. allpass) then + write(stderr, '(a)') 'SOME TESTS FAILED.' + stop + endif + end subroutine + + function test_basic_io() result(test_passed) + use m_allocator, only: allocator_t, field_t + use m_common, only: dp, DIR_C + use m_omp_common, only: SZ + + class(field_t), pointer :: arr_to_write + class(field_t), pointer :: arr_to_read + + !> Did the test pass? + logical :: test_passed + + type(allocator_t), target :: omp_allocator + class(mesh_t), allocatable :: mesh + + integer, dimension(3) :: dims_global, nproc_dir + real(dp), dimension(3) :: L_global + + type(adios_io_t) :: io + type(adios_file_t) :: write_file + type(adios_file_t) :: read_file + + integer(kind=8), dimension(3) :: ishape, istart, icount + integer, parameter :: nx = 64, ny = 32, nz = 16 + + integer :: i, j, k + + test_passed = .true. + + call io%init(MPI_COMM_WORLD) + + if(irank == 0) write(stderr, '(a)') 'Starting basic IO test.' + + !================ + ! SETUP TEST DATA + !================ + + icount = (/ nx, ny, nz/) ! local size + ishape = (/ nproc*nx, ny, nz/) ! global size + istart = (/ irank*nx, 0, 0/) ! local offset + + dims_global = (/ nx, ny, nz/) + ! Global domain dimensions (unused) + L_global = [1.0, 1.0, 1.0] + nproc_dir = [nproc, 1, 1] + mesh = mesh_t(dims_global, nproc_dir, L_global) + + omp_allocator = allocator_t(mesh, SZ) + print *, 'OpenMP allocator instantiated' + + arr_to_write => omp_allocator%get_block(DIR_C) + + ! Initialise with a simple index to verify read later + do k = 1, nz + do j = 1, ny + do i = 1, nx + arr_to_write%data(i, j, k) = (istart(1) + i-1) + (j-1)*nx*nproc + (k-1)*nx*nproc*ny + end do + end do + end do + + !================ + ! WRITE TEST DATA + !================ + + write_file = io%open_file("test_basic_io.bp", io%io_mode_write) + call io%write_field(arr_to_write, write_file, "TestArr", icount, ishape, istart) + call io%close_file(write_file) + call omp_allocator%release_block(arr_to_write) + + !================ + ! READ AND VERIFY TEST DATA + !================ + + arr_to_read => omp_allocator%get_block(DIR_C) + read_file = io%open_file("test_basic_io.bp", io%io_mode_read) + call io%read_field(arr_to_read, read_file, "TestArr", icount, ishape, istart) + call io%close_file(read_file) + + do k = 1, nz + do j = 1, ny + do i = 1, nx + if(arr_to_read%data(i, j, k) /= (istart(1) + i-1) + (j-1)*nx*nproc + (k-1)*nx*nproc*ny) then + write(stderr, '(a, f8.4, a, i8, a, i5, i5, i5)') & + 'Mismatch between read array(', arr_to_read%data(i, j, k), & + ") and expected index (", (istart(1) + i-1) + (j-1)*nx*nproc + (k-1)*nx*nproc*ny, ") at (i,j,k) = ", i, j, k + test_passed = .false. + end if + end do + end do + end do + + call omp_allocator%release_block(arr_to_read) + + call io%deinit() + + if(irank == 0) write(stderr, '(a)') 'Finishing basic IO test.' + end function + + function test_single_precision_io() result(test_passed) + use m_allocator, only: allocator_t, field_t + use m_common, only: dp, DIR_C + use m_omp_common, only: SZ + + class(field_t), pointer :: arr_to_write + + !> Did the test pass? + logical :: test_passed + + type(allocator_t), target :: omp_allocator + class(mesh_t), allocatable :: mesh + + integer, dimension(3) :: dims_global, nproc_dir + real(dp), dimension(3) :: L_global + type(adios_io_t) :: io + type(adios_file_t) :: write_file + + type(adios2_io) :: aio + type(adios2_engine) :: reader + type(adios2_variable) :: adios_var + integer(kind=8) :: idump = 0 + + integer(kind=8), dimension(3) :: ishape, istart, icount + integer, parameter :: nx = 64, ny = 32, nz = 16 + + real(kind=4) :: arr_to_read(nx, ny, nz) + + integer :: i, j, k + + test_passed = .true. + + call io%init(MPI_COMM_WORLD) + + if(irank == 0) write(stderr, '(a)') 'Starting single precision IO test.' + + !================ + ! SETUP TEST DATA + !================ + + icount = (/ nx, ny, nz/) ! local size + ishape = (/ nproc*nx, ny, nz/) ! global size + istart = (/ irank*nx, 0, 0/) ! local offset + + dims_global = (/ nx, ny, nz/) + ! Global domain dimensions (unused) + L_global = [1.0, 1.0, 1.0] + nproc_dir = [nproc, 1, 1] + mesh = mesh_t(dims_global, nproc_dir, L_global) + + omp_allocator = allocator_t(mesh, SZ) + print *, 'OpenMP allocator instantiated' + + arr_to_write => omp_allocator%get_block(DIR_C) + + ! Initialise with a simple index to verify read later + do k = 1, nz + do j = 1, ny + do i = 1, nx + arr_to_write%data(i, j, k) = (istart(1) + i-1) + (j-1)*nx*nproc + (k-1)*nx*nproc*ny + end do + end do + end do + + !================ + ! WRITE TEST DATA + !================ + + write_file = io%open_file("test_single_precision_io_data.bp", io%io_mode_write) + call io%write_field(arr_to_write, write_file, "TestArr", icount, ishape, istart, .true.) + call io%close_file(write_file) + call omp_allocator%release_block(arr_to_write) + + !================ + ! READ AND VERIFY TEST DATA + !================ + + ! Need to read manually because we haven't implemented a single-precision reader + + call adios2_declare_io (aio, io%adios_ctx, 'read', ierr) + if (.not. aio%valid) then + write(stderr, *) "Cannot create ADIOS2 IO" + test_passed = .false. + return + endif + call adios2_open(reader, aio, "test_single_precision_io_data.bp", adios2_mode_read, MPI_COMM_SELF, ierr) + if (.not. reader%valid) then + write(stderr, *) "Cannot create ADIOS2 reader" + test_passed = .false. + return + endif + call adios2_begin_step(reader, adios2_step_mode_read, ierr) + call adios2_inquire_variable(adios_var, aio, "TestArr", ierr) + if (.not. adios_var%valid) then + write(stderr, *) "Cannot fetch ADIOS2 variable" + test_passed = .false. + return + endif + + call adios2_set_step_selection(adios_var, 0_8, 1_8, ierr) + call adios2_set_selection(adios_var, 3, istart, icount, ierr) + call adios2_get(reader, adios_var, arr_to_read, ierr) + call adios2_end_step(reader, ierr) + + if (reader%valid) then + call adios2_close(reader, ierr) + endif + + do k = 1, nz + do j = 1, ny + do i = 1, nx + if(arr_to_read(i, j, k) /= (istart(1) + i-1) + (j-1)*nx*nproc + (k-1)*nx*nproc*ny) then + write(stderr, '(a, f8.4, a, i8, a, i5, i5, i5)') & + 'Mismatch between read array(', arr_to_read(i, j, k), & + ") and expected index (", (istart(1) + i-1) + (j-1)*nx*nproc + (k-1)*nx*nproc*ny, ") at (i,j,k) = ", i, j, k + test_passed = .false. + stop + end if + end do + end do + end do + + call io%deinit() + + if(irank == 0) write(stderr, '(a)') 'Finishing single precision IO test.' + end function + + function test_strided_io() result(test_passed) + use m_allocator, only: allocator_t, field_t + use m_common, only: dp, DIR_C + use m_omp_common, only: SZ + + class(field_t), pointer :: arr_to_write + real(kind=8), dimension(:,:,:), allocatable :: arr_to_read + + !> Did the test pass? + logical :: test_passed + + type(allocator_t), target :: omp_allocator + class(mesh_t), allocatable :: mesh + + integer, dimension(3) :: dims_global, nproc_dir + real(dp), dimension(3) :: L_global + type(adios_io_t) :: io + type(adios_file_t) :: write_file + + type(adios2_io) :: aio + type(adios2_engine) :: reader + type(adios2_variable) :: adios_var + integer(kind=8) :: idump = 0 + + integer(kind=8), dimension(3) :: ishape, istart, icount, istride, icount_strided, ishape_strided + integer, parameter :: nx = 32, ny = 16, nz = 64 + + integer :: i, j, k, idx + + test_passed = .true. + + call io%init(MPI_COMM_WORLD) + + if(irank == 0) write(stderr, '(a)') 'Starting strided IO test.' + + !================ + ! SETUP TEST DATA + !================ + + icount = (/ nx, ny, nz/) ! local size + ishape = (/ nproc*nx, ny, nz/) ! global size + istart = (/ irank*nx, 0, 0/) ! local offset + + dims_global = (/ nx, ny, nz/) + ! Global domain dimensions (unused) + L_global = [1.0, 1.0, 1.0] + nproc_dir = [nproc, 1, 1] + mesh = mesh_t(dims_global, nproc_dir, L_global) + + omp_allocator = allocator_t(mesh, SZ) + print *, 'OpenMP allocator instantiated' + + arr_to_write => omp_allocator%get_block(DIR_C) + + ! Initialise with a simple index to verify read later + do k = 1, nz + do j = 1, ny + do i = 1, nx + arr_to_write%data(i, j, k) = (istart(1) + i-1) + (j-1)*nx*nproc + (k-1)*nx*nproc*ny + end do + end do + end do + + !================ + ! WRITE TEST DATA + !================ + + istride = [2_8, 3_8, 4_8] + icount_strided = icount / istride + + write_file = io%open_file("test_strided_io_data.bp", io%io_mode_write) + call io%write_field(arr_to_write, write_file, "TestArr", icount, ishape, istart, istride_in=istride) + call io%close_file(write_file) + call omp_allocator%release_block(arr_to_write) + + !================ + ! READ AND VERIFY TEST DATA + !================ + + ! Read only from one rank; let's not overcomplicate this + if(irank == 0) then + ishape_strided = ishape/istride + allocate(arr_to_read(ishape_strided(1), ishape_strided(2), ishape_strided(3))) + + call adios2_declare_io (aio, io%adios_ctx, 'read', ierr) + if (.not. aio%valid) then + write(stderr, *) "Cannot create ADIOS2 IO" + test_passed = .false. + return + endif + call adios2_open(reader, aio, "test_strided_io_data.bp", adios2_mode_read, MPI_COMM_SELF, ierr) + if (.not. reader%valid) then + write(stderr, *) "Cannot create ADIOS2 reader" + test_passed = .false. + return + endif + call adios2_begin_step(reader, adios2_step_mode_read, ierr) + call adios2_inquire_variable(adios_var, aio, "TestArr", ierr) + if (.not. adios_var%valid) then + write(stderr, *) "Cannot fetch ADIOS2 variable" + test_passed = .false. + return + endif + + call adios2_set_step_selection(adios_var, 0_8, 1_8, ierr) + call adios2_set_selection(adios_var, 3, [0_8,0_8,0_8], ishape_strided, ierr) + call adios2_get(reader, adios_var, arr_to_read, ierr) + call adios2_end_step(reader, ierr) + + if (reader%valid) then + call adios2_close(reader, ierr) + endif + + do k = 1, ishape_strided(3) + do j = 1, ishape_strided(2) + do i = 1, ishape_strided(1) + idx = (istart(1) + i-1)*istride(1) + (j-1)*istride(2)*nx*nproc + (k-1)*istride(3)*nx*nproc*ny + if(arr_to_read(i, j, k) /= idx) then + write(stderr, '(a, f8.4, a, i8, a, i5, i5, i5)') & + 'Mismatch between read array(', arr_to_read(i, j, k), & + ") and expected index (", idx, ") at (i,j,k) = ", i, j, k + test_passed = .false. + end if + end do + end do + end do + + deallocate(arr_to_read) + endif + + call io%deinit() + + if(irank == 0) write(stderr, '(a)') 'Finishing strided IO test.' + end function + + function test_strided_single_precision_io() result(test_passed) + use m_allocator, only: allocator_t, field_t + use m_common, only: dp, DIR_C + use m_omp_common, only: SZ + + class(field_t), pointer :: arr_to_write + real(kind=4), dimension(:,:,:), allocatable :: arr_to_read + + !> Did the test pass? + logical :: test_passed + + type(allocator_t), target :: omp_allocator + class(mesh_t), allocatable :: mesh + + integer, dimension(3) :: dims_global, nproc_dir + real(dp), dimension(3) :: L_global + type(adios_io_t) :: io + type(adios_file_t) :: write_file + + type(adios2_io) :: aio + type(adios2_engine) :: reader + type(adios2_variable) :: adios_var + integer(kind=8) :: idump = 0 + + integer(kind=8), dimension(3) :: ishape, istart, icount, istride, icount_strided, ishape_strided + integer, parameter :: nx = 32, ny = 16, nz = 64 + + integer :: i, j, k, idx + + test_passed = .true. + + call io%init(MPI_COMM_WORLD) + + if(irank == 0) write(stderr, '(a)') 'Starting strided IO test.' + + !================ + ! SETUP TEST DATA + !================ + + icount = (/ nx, ny, nz/) ! local size + ishape = (/ nproc*nx, ny, nz/) ! global size + istart = (/ irank*nx, 0, 0/) ! local offset + + dims_global = (/ nx, ny, nz/) + ! Global domain dimensions (unused) + L_global = [1.0, 1.0, 1.0] + nproc_dir = [nproc, 1, 1] + mesh = mesh_t(dims_global, nproc_dir, L_global) + + omp_allocator = allocator_t(mesh, SZ) + print *, 'OpenMP allocator instantiated' + + arr_to_write => omp_allocator%get_block(DIR_C) + + ! Initialise with a simple index to verify read later + do k = 1, nz + do j = 1, ny + do i = 1, nx + arr_to_write%data(i, j, k) = (istart(1) + i-1) + (j-1)*nx*nproc + (k-1)*nx*nproc*ny + end do + end do + end do + + !================ + ! WRITE TEST DATA + !================ + + istride = [2_8, 3_8, 4_8] + icount_strided = icount / istride + + write_file = io%open_file("test_strided_io_sp_data.bp", io%io_mode_write) + call io%write_field(arr_to_write, write_file, "TestArr", icount, ishape, istart, convert_to_sp_in=.true., istride_in=istride) + call io%close_file(write_file) + call omp_allocator%release_block(arr_to_write) + + !================ + ! READ AND VERIFY TEST DATA + !================ + + ! Read only from one rank; let's not overcomplicate this + if(irank == 0) then + ishape_strided = ishape/istride + allocate(arr_to_read(ishape_strided(1), ishape_strided(2), ishape_strided(3))) + + call adios2_declare_io (aio, io%adios_ctx, 'read', ierr) + if (.not. aio%valid) then + write(stderr, *) "Cannot create ADIOS2 IO" + test_passed = .false. + return + endif + call adios2_open(reader, aio, "test_strided_io_sp_data.bp", adios2_mode_read, MPI_COMM_SELF, ierr) + if (.not. reader%valid) then + write(stderr, *) "Cannot create ADIOS2 reader" + test_passed = .false. + return + endif + call adios2_begin_step(reader, adios2_step_mode_read, ierr) + call adios2_inquire_variable(adios_var, aio, "TestArr", ierr) + if (.not. adios_var%valid) then + write(stderr, *) "Cannot fetch ADIOS2 variable" + test_passed = .false. + return + endif + + call adios2_set_step_selection(adios_var, 0_8, 1_8, ierr) + call adios2_set_selection(adios_var, 3, [0_8,0_8,0_8], ishape_strided, ierr) + call adios2_get(reader, adios_var, arr_to_read, ierr) + call adios2_end_step(reader, ierr) + + if (reader%valid) then + call adios2_close(reader, ierr) + endif + + do k = 1, ishape_strided(3) + do j = 1, ishape_strided(2) + do i = 1, ishape_strided(1) + idx = (istart(1) + i-1)*istride(1) + (j-1)*istride(2)*nx*nproc + (k-1)*istride(3)*nx*nproc*ny + if(arr_to_read(i, j, k) /= idx) then + write(stderr, '(a, f8.4, a, i8, a, i5, i5, i5)') & + 'Mismatch between read array(', arr_to_read(i, j, k), & + ") and expected index (", idx, ") at (i,j,k) = ", i, j, k + test_passed = .false. + end if + end do + end do + end do + + deallocate(arr_to_read) + endif + + call io%deinit() + + if(irank == 0) write(stderr, '(a)') 'Finishing strided IO test.' + end function + +end program