m_data_input.f90

Go to the documentation of this file.

!>
!! @file
!> @brief Contains module m_data_input
 
!> @brief Reads raw simulation grid and conservative-variable data for a given time-step and fills buffer regions
module m_data_input
 
#ifdef MFC_MPI
    use mpi
#endif
 
    use m_derived_types
    use m_global_parameters
    use m_mpi_proxy
    use m_mpi_common
    use m_compile_specific
    use m_boundary_common
    use m_boundary_io
    use m_helper
 
    implicit none
 
    private; public :: s_initialize_data_input_module, s_read_data_files, s_read_serial_data_files, s_read_parallel_data_files, &
        & s_finalize_data_input_module
 
    abstract interface
 
        !> Subroutine for reading data files
        impure subroutine s_read_abstract_data_files(t_step)
 
            implicit none
 
            integer, intent(in) :: t_step
 
        end subroutine s_read_abstract_data_files
    end interface
 
    type(scalar_field), allocatable, dimension(:), public    :: q_cons_vf  !< Conservative variables
    type(scalar_field), allocatable, dimension(:), public    :: q_cons_temp
    type(scalar_field), allocatable, dimension(:), public    :: q_prim_vf  !< Primitive variables
    type(integer_field), allocatable, dimension(:,:), public :: bc_type    !< Boundary condition identifiers
    type(scalar_field), public                               :: q_t_sf     !< Temperature field
    ! type(scalar_field), public :: ib_markers !<
    type(integer_field), public :: ib_markers
 
    procedure(s_read_abstract_data_files), pointer :: s_read_data_files => null()
 
contains
 
    !> Helper subroutine to read grid data files for a given direction
    impure subroutine s_read_grid_data_direction(t_step_dir, direction, cb_array, d_array, cc_array, size_dim)
 
        character(len=*), intent(in)             :: t_step_dir
        character(len=1), intent(in)             :: direction
        real(wp), dimension(-1:), intent(out)    :: cb_array
        real(wp), dimension(0:), intent(out)     :: d_array
        real(wp), dimension(0:), intent(out)     :: cc_array
        integer, intent(in)                      :: size_dim
        character(LEN=len_trim(t_step_dir) + 10) :: file_loc
        logical                                  :: file_check
 
        file_loc = trim(t_step_dir) // '/' // direction // '_cb.dat'
        inquire (file=trim(file_loc), exist=file_check)
 
        if (file_check) then
            open (1, file=trim(file_loc), form='unformatted', status='old', action='read')
            read (1) cb_array(-1:size_dim)
            close (1)
        else
            call s_mpi_abort('File ' // direction // '_cb.dat is missing in ' // trim(t_step_dir) // '. Exiting.')
        end if
 
        d_array(0:size_dim) = cb_array(0:size_dim) - cb_array(-1:size_dim - 1)
        cc_array(0:size_dim) = cb_array(-1:size_dim - 1) + d_array(0:size_dim)/2._wp
 
    end subroutine s_read_grid_data_direction
 
#ifdef MFC_MPI
    !> Helper subroutine to setup MPI data I/O parameters
    impure subroutine s_setup_mpi_io_params(data_size, m_MOK, n_MOK, p_MOK, WP_MOK, MOK, str_MOK, NVARS_MOK)
 
        integer, intent(out)                       :: data_size
        integer(KIND=MPI_OFFSET_KIND), intent(out) :: m_mok, n_mok, p_mok
        integer(KIND=MPI_OFFSET_KIND), intent(out) :: wp_mok, mok, str_mok, nvars_mok
 
        if (ib) then
            call s_initialize_mpi_data(q_cons_vf, ib_markers)
        else
            call s_initialize_mpi_data(q_cons_vf)
        end if
 
        data_size = (m + 1)*(n + 1)*(p + 1)
 
        m_mok = int(m_glb + 1, mpi_offset_kind)
        n_mok = int(n_glb + 1, mpi_offset_kind)
        p_mok = int(p_glb + 1, mpi_offset_kind)
        wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
        mok = int(1._wp, mpi_offset_kind)
        str_mok = int(name_len, mpi_offset_kind)
        nvars_mok = int(sys_size, mpi_offset_kind)
 
    end subroutine s_setup_mpi_io_params
#endif
 
    !> Helper subroutine to read IB data files
    impure subroutine s_read_ib_data_files(file_loc_base, t_step)
 
        character(len=*), intent(in)                :: file_loc_base
        integer, intent(in), optional               :: t_step
        character(LEN=len_trim(file_loc_base) + 20) :: file_loc
        logical                                     :: file_exist
        integer                                     :: ifile, ierr, data_size
 
#ifdef MFC_MPI
        integer, dimension(MPI_STATUS_SIZE) :: status
        integer(KIND=MPI_OFFSET_KIND)       :: disp
        integer(KIND=MPI_OFFSET_KIND)       :: m_mok, n_mok, p_mok, mok, wp_mok, var_mok
        integer                             :: save_index
#endif
 
        if (.not. ib) return
 
        if (parallel_io) then
            write (file_loc, '(A)') trim(file_loc_base) // 'ib.dat'
        else
            write (file_loc, '(A)') trim(file_loc_base) // '/ib_data.dat'
        end if
        inquire (file=trim(file_loc), exist=file_exist)
 
        if (file_exist) then
            if (parallel_io) then
#ifdef MFC_MPI
                call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
 
                m_mok = int(m_glb + 1, mpi_offset_kind)
                n_mok = int(n_glb + 1, mpi_offset_kind)
                p_mok = int(p_glb + 1, mpi_offset_kind)
                mok = int(1._wp, mpi_offset_kind)
                wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
                save_index = t_step/t_step_save  ! get the number of saves done to this point
 
                data_size = (m + 1)*(n + 1)*(p + 1)
                var_mok = int(sys_size + 1, mpi_offset_kind)
                if (t_step == 0) then
                    disp = 0
                else
                    disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1 + int(save_index, mpi_offset_kind))
                end if
 
                call mpi_file_set_view(ifile, disp, mpi_integer, mpi_io_ib_data%view, 'native', mpi_info_int, ierr)
                call mpi_file_read(ifile, mpi_io_ib_data%var%sf, data_size, mpi_integer, status, ierr)
 
                call mpi_file_close(ifile, ierr)
#endif
            else
                open (2, file=trim(file_loc), form='unformatted', action='read', status='old')
                read (2) ib_markers%sf(0:m,0:n,0:p)
                close (2)
            end if
        else
            call s_mpi_abort('File ' // trim(file_loc) // ' is missing. Exiting.')
        end if
 
    end subroutine s_read_ib_data_files
 
    !> Helper subroutine to allocate field arrays for given dimensionality
    impure subroutine s_allocate_field_arrays(local_start_idx, end_x, end_y, end_z)
 
        integer, intent(in) :: local_start_idx, end_x, end_y, end_z
        integer             :: i
 
        do i = 1, sys_size
            allocate (q_cons_vf(i)%sf(local_start_idx:end_x,local_start_idx:end_y,local_start_idx:end_z))
            allocate (q_prim_vf(i)%sf(local_start_idx:end_x,local_start_idx:end_y,local_start_idx:end_z))
        end do
 
        if (ib) then
            allocate (ib_markers%sf(local_start_idx:end_x,local_start_idx:end_y,local_start_idx:end_z))
        end if
 
        if (chemistry) then
            allocate (q_t_sf%sf(local_start_idx:end_x,local_start_idx:end_y,local_start_idx:end_z))
        end if
 
    end subroutine s_allocate_field_arrays
 
    !> Read the raw data files present in the corresponding time-step directory and to populate the associated grid and conservative
    !! variables.
    impure subroutine s_read_serial_data_files(t_step)
 
        integer, intent(in)                                      :: t_step
        character(LEN=len_trim(case_dir) + 2*name_len)           :: t_step_dir
        character(LEN=len_trim(case_dir) + 3*name_len)           :: file_loc
        character(len=int(floor(log10(real(sys_size, wp)))) + 1) :: file_num
        logical                                                  :: dir_check
        logical                                                  :: file_check
        integer                                                  :: i
 
        write (t_step_dir, '(A,I0,A,I0)') '/p_all/p', proc_rank, '/', t_step
        t_step_dir = trim(case_dir) // trim(t_step_dir)
 
        file_loc = trim(t_step_dir) // '/.'
        call my_inquire(file_loc, dir_check)
 
        if (dir_check .neqv. .true.) then
            call s_mpi_abort('Time-step folder ' // trim(t_step_dir) // ' is missing. Exiting.')
        end if
 
        if (bc_io) then
            call s_read_serial_boundary_condition_files(t_step_dir, bc_type)
        else
            call s_assign_default_bc_type(bc_type)
        end if
 
        ! Pass explicit slices so the dummy `dimension(-1:)` / `dimension(0:)` arguments map to the correct interior indices of the
        ! actual arrays. Without slicing, when offset_x%beg or buff_size > 0 (i.e. format=1 parallel 3D ranks), Fortran's
        ! assumed-shape re-mapping shifts the read by that many slots and leaves the last interior cells uninitialized - corrupting
        ! downstream ghost-cell extrapolation.
        call s_read_grid_data_direction(t_step_dir, 'x', x_cb(-1:m), dx(0:m), x_cc(0:m), m)
 
        if (n > 0) then
            call s_read_grid_data_direction(t_step_dir, 'y', y_cb(-1:n), dy(0:n), y_cc(0:n), n)
 
            if (p > 0) then
                call s_read_grid_data_direction(t_step_dir, 'z', z_cb(-1:p), dz(0:p), z_cc(0:p), p)
            end if
        end if
 
        do i = 1, sys_size
            write (file_num, '(I0)') i
            file_loc = trim(t_step_dir) // '/q_cons_vf' // trim(file_num) // '.dat'
            inquire (file=trim(file_loc), exist=file_check)
 
            if (file_check) then
                open (1, file=trim(file_loc), form='unformatted', status='old', action='read')
                read (1) q_cons_vf(i)%sf(0:m,0:n,0:p)
                close (1)
            else
                call s_mpi_abort('File q_cons_vf' // trim(file_num) // '.dat is missing in ' // trim(t_step_dir) // '. Exiting.')
            end if
        end do
 
        call s_read_ib_data_files(t_step_dir)
 
    end subroutine s_read_serial_data_files
 
    !> Parallel-read the raw data files present in the corresponding time-step directory and to populate the associated grid and
    !! conservative variables.
    impure subroutine s_read_parallel_data_files(t_step)
 
        integer, intent(in) :: t_step
 
#ifdef MFC_MPI
        real(wp), allocatable, dimension(:)  :: x_cb_glb, y_cb_glb, z_cb_glb
        integer                              :: ifile, ierr, data_size, filetype, stride
        integer, dimension(MPI_STATUS_SIZE)  :: status
        integer(KIND=MPI_OFFSET_KIND)        :: disp
        integer(KIND=MPI_OFFSET_KIND)        :: m_mok, n_mok, p_mok
        integer(KIND=MPI_OFFSET_KIND)        :: wp_mok, var_mok, str_mok
        integer(KIND=MPI_OFFSET_KIND)        :: nvars_mok
        integer(KIND=MPI_OFFSET_KIND)        :: mok
        integer(kind=MPI_OFFSET_KIND)        :: offset
        character(LEN=path_len + 2*name_len) :: file_loc
        logical                              :: file_exist
        character(len=10)                    :: t_step_string
        integer                              :: i
 
        allocate (x_cb_glb(-1:m_glb))
        allocate (y_cb_glb(-1:n_glb))
        allocate (z_cb_glb(-1:p_glb))
 
        if (down_sample) then
            stride = 3
        else
            stride = 1
        end if
 
        file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // 'x_cb.dat'
        inquire (file=trim(file_loc), exist=file_exist)
 
        if (file_exist) then
            data_size = m_glb + 2
            call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
 
            call mpi_type_vector(data_size, 1, stride, mpi_p, filetype, ierr)
            call mpi_type_commit(filetype, ierr)
 
            offset = 0
            call mpi_file_set_view(ifile, offset, mpi_p, filetype, 'native', mpi_info_int, ierr)
 
            call mpi_file_read(ifile, x_cb_glb, data_size, mpi_p, status, ierr)
            call mpi_file_close(ifile, ierr)
        else
            call s_mpi_abort('File ' // trim(file_loc) // ' is missing. Exiting.')
        end if
 
        ! Bitwise-consistent grid distribution from the global file
        call s_apply_grid_from_global_dim(x_cb_glb, m_glb, m, start_idx(1), bc_x%beg, bc_x%end, offset_x%beg, offset_x%end, &
                                          & buff_size, buff_size, x_cb, x_cc, dx)
 
        if (n > 0) then
            file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // 'y_cb.dat'
            inquire (file=trim(file_loc), exist=file_exist)
 
            if (file_exist) then
                data_size = n_glb + 2
                call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
 
                call mpi_type_vector(data_size, 1, stride, mpi_p, filetype, ierr)
                call mpi_type_commit(filetype, ierr)
 
                offset = 0
                call mpi_file_set_view(ifile, offset, mpi_p, filetype, 'native', mpi_info_int, ierr)
 
                call mpi_file_read(ifile, y_cb_glb, data_size, mpi_p, status, ierr)
                call mpi_file_close(ifile, ierr)
            else
                call s_mpi_abort('File ' // trim(file_loc) // ' is missing. Exiting.')
            end if
 
            call s_apply_grid_from_global_dim(y_cb_glb, n_glb, n, start_idx(2), bc_y%beg, bc_y%end, offset_y%beg, offset_y%end, &
                                              & buff_size, buff_size, y_cb, y_cc, dy)
 
            if (p > 0) then
                file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // 'z_cb.dat'
                inquire (file=trim(file_loc), exist=file_exist)
 
                if (file_exist) then
                    data_size = p_glb + 2
                    call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
 
                    call mpi_type_vector(data_size, 1, stride, mpi_p, filetype, ierr)
                    call mpi_type_commit(filetype, ierr)
 
                    offset = 0
                    call mpi_file_set_view(ifile, offset, mpi_p, filetype, 'native', mpi_info_int, ierr)
 
                    call mpi_file_read(ifile, z_cb_glb, data_size, mpi_p, status, ierr)
                    call mpi_file_close(ifile, ierr)
                else
                    call s_mpi_abort('File ' // trim(file_loc) // ' is missing. Exiting.')
                end if
 
                call s_apply_grid_from_global_dim(z_cb_glb, p_glb, p, start_idx(3), bc_z%beg, bc_z%end, offset_z%beg, &
                                                  & offset_z%end, buff_size, buff_size, z_cb, z_cc, dz)
            end if
        end if
 
        call s_read_parallel_conservative_data(t_step, m_mok, n_mok, p_mok, wp_mok, mok, str_mok, nvars_mok)
 
        deallocate (x_cb_glb, y_cb_glb, z_cb_glb)
 
        if (bc_io) then
            call s_read_parallel_boundary_condition_files(bc_type)
        else
            call s_assign_default_bc_type(bc_type)
        end if
#endif
 
    end subroutine s_read_parallel_data_files
 
#ifdef MFC_MPI
    !> Helper subroutine to read parallel conservative variable data
    impure subroutine s_read_parallel_conservative_data(t_step, m_MOK, n_MOK, p_MOK, WP_MOK, MOK, str_MOK, NVARS_MOK)
 
        integer, intent(in)                          :: t_step
        integer(KIND=MPI_OFFSET_KIND), intent(inout) :: m_mok, n_mok, p_mok
        integer(KIND=MPI_OFFSET_KIND), intent(inout) :: wp_mok, mok, str_mok, nvars_mok
        integer                                      :: ifile, ierr, data_size
        integer, dimension(MPI_STATUS_SIZE)          :: status
        integer(KIND=MPI_OFFSET_KIND)                :: disp, var_mok
        character(LEN=path_len + 2*name_len)         :: file_loc
        logical                                      :: file_exist
        character(len=10)                            :: t_step_string
        integer                                      :: i
 
        if (file_per_process) then
            call s_int_to_str(t_step, t_step_string)
            write (file_loc, '(I0,A1,I7.7,A)') t_step, '_', proc_rank, '.dat'
            file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // trim(mpiiofs) // trim(file_loc)
            inquire (file=trim(file_loc), exist=file_exist)
 
            if (file_exist) then
                call mpi_file_open(mpi_comm_self, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
 
                if (down_sample) then
                    call s_initialize_mpi_data_ds(q_cons_temp)
                else
                    if (ib) then
                        call s_initialize_mpi_data(q_cons_vf, ib_markers)
                    else
                        call s_initialize_mpi_data(q_cons_vf)
                    end if
                end if
 
                if (down_sample) then
                    data_size = (m + 3)*(n + 3)*(p + 3)
                else
                    data_size = (m + 1)*(n + 1)*(p + 1)
                end if
 
                m_mok = int(m_glb + 1, mpi_offset_kind)
                n_mok = int(n_glb + 1, mpi_offset_kind)
                p_mok = int(p_glb + 1, mpi_offset_kind)
                wp_mok = int(storage_size(0._stp)/8, mpi_offset_kind)
                mok = int(1._wp, mpi_offset_kind)
                str_mok = int(name_len, mpi_offset_kind)
                nvars_mok = int(sys_size, mpi_offset_kind)
 
                if (bubbles_euler .or. elasticity .or. mhd) then
                    do i = 1, sys_size
                        var_mok = int(i, mpi_offset_kind)
                        call mpi_file_read_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                    end do
                else
                    do i = 1, sys_size
                        var_mok = int(i, mpi_offset_kind)
                        call mpi_file_read_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                    end do
                end if
 
                call s_mpi_barrier()
                call mpi_file_close(ifile, ierr)
 
                if (down_sample) then
                    do i = 1, sys_size
                        q_cons_vf(i)%sf(0:m,0:n,0:p) = q_cons_temp(i)%sf(0:m,0:n,0:p)
                    end do
                end if
 
                call s_read_ib_data_files(trim(case_dir) // '/restart_data' // trim(mpiiofs), t_step)
            else
                call s_mpi_abort('File ' // trim(file_loc) // ' is missing. Exiting.')
            end if
        else
            write (file_loc, '(I0,A)') t_step, '.dat'
            file_loc = trim(case_dir) // '/restart_data' // trim(mpiiofs) // trim(file_loc)
            inquire (file=trim(file_loc), exist=file_exist)
 
            if (file_exist) then
                call mpi_file_open(mpi_comm_world, file_loc, mpi_mode_rdonly, mpi_info_int, ifile, ierr)
 
                call s_setup_mpi_io_params(data_size, m_mok, n_mok, p_mok, wp_mok, mok, str_mok, nvars_mok)
 
                do i = 1, sys_size
                    var_mok = int(i, mpi_offset_kind)
 
                    disp = m_mok*max(mok, n_mok)*max(mok, p_mok)*wp_mok*(var_mok - 1)
 
                    call mpi_file_set_view(ifile, disp, mpi_p, mpi_io_data%view(i), 'native', mpi_info_int, ierr)
                    call mpi_file_read_all(ifile, mpi_io_data%var(i)%sf, data_size*mpi_io_type, mpi_io_p, status, ierr)
                end do
 
                call s_mpi_barrier()
                call mpi_file_close(ifile, ierr)
 
                call s_read_ib_data_files(trim(case_dir) // '/restart_data' // trim(mpiiofs), t_step)
            else
                call s_mpi_abort('File ' // trim(file_loc) // ' is missing. Exiting.')
            end if
        end if
 
    end subroutine s_read_parallel_conservative_data
#endif
 
    !> Computation of parameters, allocation procedures, and/or any other tasks needed to properly setup the module
    impure subroutine s_initialize_data_input_module
 
        integer :: i
 
        allocate (q_cons_vf(1:sys_size))
        allocate (q_prim_vf(1:sys_size))
        allocate (q_cons_temp(1:sys_size))
 
        if (n > 0) then
            if (p > 0) then
                call s_allocate_field_arrays(-buff_size, m + buff_size, n + buff_size, p + buff_size)
                if (down_sample) then
                    do i = 1, sys_size
                        allocate (q_cons_temp(i)%sf(-1:m + 1,-1:n + 1,-1:p + 1))
                    end do
                end if
            else
                call s_allocate_field_arrays(-buff_size, m + buff_size, n + buff_size, 0)
            end if
        else
            call s_allocate_field_arrays(-buff_size, m + buff_size, 0, 0)
        end if
 
        allocate (bc_type(1:num_dims,1:2))
 
        allocate (bc_type(1, 1)%sf(0:0,0:n,0:p))
        allocate (bc_type(1, 2)%sf(0:0,0:n,0:p))
        if (n > 0) then
            allocate (bc_type(2, 1)%sf(-buff_size:m + buff_size,0:0,0:p))
            allocate (bc_type(2, 2)%sf(-buff_size:m + buff_size,0:0,0:p))
            if (p > 0) then
                allocate (bc_type(3, 1)%sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0))
                allocate (bc_type(3, 2)%sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0))
            end if
        end if
 
        if (parallel_io .neqv. .true.) then
            s_read_data_files => s_read_serial_data_files
        else
            s_read_data_files => s_read_parallel_data_files
        end if
 
    end subroutine s_initialize_data_input_module
 
    !> Deallocation procedures for the module
    impure subroutine s_finalize_data_input_module
 
        integer :: i
 
        do i = 1, sys_size
            deallocate (q_cons_vf(i)%sf)
            deallocate (q_prim_vf(i)%sf)
            if (down_sample) then
                deallocate (q_cons_temp(i)%sf)
            end if
        end do
 
        deallocate (q_cons_vf)
        deallocate (q_prim_vf)
        deallocate (q_cons_temp)
 
        if (ib) then
            deallocate (ib_markers%sf)
        end if
 
        if (chemistry) then
            deallocate (q_t_sf%sf)
        end if
 
        deallocate (bc_type(1, 1)%sf, bc_type(1, 2)%sf)
        if (n > 0) then
            deallocate (bc_type(2, 1)%sf, bc_type(2, 2)%sf)
            if (p > 0) then
                deallocate (bc_type(3, 1)%sf, bc_type(3, 2)%sf)
            end if
        end if
 
        deallocate (bc_type)
 
        s_read_data_files => null()
 
    end subroutine s_finalize_data_input_module
 
end module m_data_input