read_snapshot Subroutine

public subroutine read_snapshot()

Uses

  • proc~~read_snapshot~~UsesGraph proc~read_snapshot read_snapshot mod_amr_solution_node mod_amr_solution_node proc~read_snapshot->mod_amr_solution_node module~mod_forest mod_forest proc~read_snapshot->module~mod_forest module~mod_functions_forest mod_functions_forest proc~read_snapshot->module~mod_functions_forest module~mod_global_parameters mod_global_parameters proc~read_snapshot->module~mod_global_parameters module~mod_input_output_helper mod_input_output_helper proc~read_snapshot->module~mod_input_output_helper module~mod_usr_methods mod_usr_methods proc~read_snapshot->module~mod_usr_methods module~mod_mpi_wrapper mod_mpi_wrapper module~mod_functions_forest->module~mod_mpi_wrapper module~mod_basic_types mod_basic_types module~mod_global_parameters->module~mod_basic_types module~mod_connectivity mod_connectivity module~mod_global_parameters->module~mod_connectivity module~mod_constants mod_constants module~mod_global_parameters->module~mod_constants module~mod_physicaldata mod_physicaldata module~mod_global_parameters->module~mod_physicaldata module~mod_variables mod_variables module~mod_global_parameters->module~mod_variables mpi mpi module~mod_global_parameters->mpi module~mod_mpi_wrapper->mpi module~mod_variables->module~mod_basic_types

Routine to read in snapshots (.dat files). When it cannot recognize the file version, it will automatically try the 'old' reader.

Arguments

None

Calls

proc~~read_snapshot~~CallsGraph proc~read_snapshot read_snapshot alloc_node alloc_node proc~read_snapshot->alloc_node mpi_barrier mpi_barrier proc~read_snapshot->mpi_barrier mpi_bcast mpi_bcast proc~read_snapshot->mpi_bcast mpi_file_close mpi_file_close proc~read_snapshot->mpi_file_close mpi_file_open mpi_file_open proc~read_snapshot->mpi_file_open mpi_file_seek mpi_file_seek proc~read_snapshot->mpi_file_seek proc~count_ix count_ix proc~read_snapshot->proc~count_ix proc~get_snapshot_index get_snapshot_index proc~read_snapshot->proc~get_snapshot_index proc~mpi_file_read_wrapper mpi_file_read_wrapper proc~read_snapshot->proc~mpi_file_read_wrapper proc~mpi_recv_wrapper mpi_recv_wrapper proc~read_snapshot->proc~mpi_recv_wrapper proc~mpi_send_wrapper mpi_send_wrapper proc~read_snapshot->proc~mpi_send_wrapper proc~mpistop mpistop proc~read_snapshot->proc~mpistop proc~read_forest read_forest proc~read_snapshot->proc~read_forest proc~read_snapshot_old read_snapshot_old proc~read_snapshot->proc~read_snapshot_old proc~snapshot_read_header snapshot_read_header proc~read_snapshot->proc~snapshot_read_header mpi_file_read mpi_file_read proc~mpi_file_read_wrapper->mpi_file_read mpi_recv mpi_recv proc~mpi_recv_wrapper->mpi_recv mpi_send mpi_send proc~mpi_send_wrapper->mpi_send mpi_abort mpi_abort proc~mpistop->mpi_abort proc~read_forest->mpi_bcast proc~read_forest->mpi_file_seek proc~read_forest->proc~mpi_file_read_wrapper getnode getnode proc~read_forest->getnode mpi_get_count mpi_get_count proc~read_forest->mpi_get_count proc~add_to_linked_list add_to_linked_list proc~read_forest->proc~add_to_linked_list proc~amr_morton_order amr_Morton_order proc~read_forest->proc~amr_morton_order proc~asign_tree_neighbor asign_tree_neighbor proc~read_forest->proc~asign_tree_neighbor proc~build_connectivity build_connectivity proc~read_forest->proc~build_connectivity proc~get_level_range get_level_range proc~read_forest->proc~get_level_range proc~get_morton_range get_Morton_range proc~read_forest->proc~get_morton_range proc~getigrids getigrids proc~read_forest->proc~getigrids proc~level1_morton_order level1_Morton_order proc~read_forest->proc~level1_morton_order proc~read_snapshot_old->alloc_node proc~read_snapshot_old->mpi_barrier proc~read_snapshot_old->mpi_bcast proc~read_snapshot_old->mpi_file_close proc~read_snapshot_old->mpi_file_open proc~read_snapshot_old->mpi_file_seek proc~read_snapshot_old->proc~mpi_file_read_wrapper proc~read_snapshot_old->proc~mpi_recv_wrapper proc~read_snapshot_old->proc~mpi_send_wrapper proc~read_snapshot_old->proc~mpistop proc~read_snapshot_old->proc~read_forest proc~mpi_file_read_at_wrapper mpi_file_read_at_wrapper proc~read_snapshot_old->proc~mpi_file_read_at_wrapper proc~snapshot_read_header->proc~get_snapshot_index proc~snapshot_read_header->proc~mpi_file_read_wrapper proc~snapshot_read_header->proc~mpistop proc~amr_morton_order->proc~mpistop proc~find_neighbor find_neighbor proc~asign_tree_neighbor->proc~find_neighbor proc~add_to_c nbprocs_info_t%add_to_c proc~build_connectivity->proc~add_to_c proc~add_to_f nbprocs_info_t%add_to_f proc~build_connectivity->proc~add_to_f proc~add_to_srl nbprocs_info_t%add_to_srl proc~build_connectivity->proc~add_to_srl proc~alloc_buffers_c nbprocs_info_t%alloc_buffers_c proc~build_connectivity->proc~alloc_buffers_c proc~alloc_buffers_f nbprocs_info_t%alloc_buffers_f proc~build_connectivity->proc~alloc_buffers_f proc~alloc_buffers_srl nbprocs_info_t%alloc_buffers_srl proc~build_connectivity->proc~alloc_buffers_srl proc~build_connectivity->proc~find_neighbor proc~identifyphysbound_connectivity identifyphysbound_connectivity proc~build_connectivity->proc~identifyphysbound_connectivity proc~reset nbprocs_info_t%reset proc~build_connectivity->proc~reset proc~mortonencode mortonEncode proc~level1_morton_order->proc~mortonencode proc~mpi_file_read_at_wrapper->mpi_file_read proc~add_igrid_to_c nbprocs_info_t%add_igrid_to_c proc~add_to_c->proc~add_igrid_to_c proc~add_igrid_to_f nbprocs_info_t%add_igrid_to_f proc~add_to_f->proc~add_igrid_to_f proc~add_igrid_to_srl nbprocs_info_t%add_igrid_to_srl proc~add_to_srl->proc~add_igrid_to_srl proc~idecode nbprocs_info_t%idecode proc~alloc_buffers_srl->proc~idecode proc~find_root_neighbor find_root_neighbor proc~find_neighbor->proc~find_root_neighbor proc~init_cf_info nbinfo_cf_t%init_cf_info proc~add_igrid_to_c->proc~init_cf_info proc~add_igrid_to_f->proc~init_cf_info proc~iencode nbprocs_info_t%iencode proc~add_igrid_to_srl->proc~iencode proc~init_srl_info init_srl_info proc~add_igrid_to_srl->proc~init_srl_info

Called by

proc~~read_snapshot~~CalledByGraph proc~read_snapshot read_snapshot proc~main main proc~main->proc~read_snapshot program~agile agile program~agile->proc~main

Source Code

  subroutine read_snapshot
    use mod_usr_methods, only: usr_transform_w
    use mod_input_output_helper, only: count_ix
    use mod_forest
    use mod_global_parameters
    use mod_amr_solution_node, only: alloc_node
    use mod_functions_forest, only: read_forest

    integer                       :: ix_buffer(2*ndim+1), n_values,&
        n_values_stagger
    integer                       :: ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,&
       ixOmax3, ixOsmin1,ixOsmin2,ixOsmin3,ixOsmax1,ixOsmax2,ixOsmax3
    integer                       :: file_handle, amode, igrid, Morton_no,&
        iread
    integer                       :: istatus(MPI_STATUS_SIZE)
    integer                       :: iorecvstatus(MPI_STATUS_SIZE)
    integer                       :: ipe,inrecv,nrecv, file_version
    logical                       :: fexist
    integer(MPI_OFFSET_KIND)      :: offset_tree_info
    integer(MPI_OFFSET_KIND)      :: offset_block_data
    double precision, allocatable :: w_buffer(:)
    double precision, dimension(:,:,:,:), allocatable :: w
    double precision :: ws(ixGslo1:ixGshi1,ixGslo2:ixGshi2,ixGslo3:ixGshi3,&
       1:ndim)

    if (mype==0) then
      inquire(file=trim(restart_from_file), exist=fexist)
      if (.not.fexist) call mpistop(trim(restart_from_file)//" not found!")

      call MPI_FILE_OPEN(MPI_COMM_SELF,restart_from_file,MPI_MODE_RDONLY,&
          MPI_INFO_NULL,file_handle,ierrmpi)
      call mpi_file_read_wrapper(file_handle, file_version, 1, MPI_INTEGER, istatus,&
          ierrmpi)
    end if

    call MPI_BCAST(file_version,1,MPI_INTEGER,0,icomm,ierrmpi)

    if (all(compatible_versions /= file_version)) then
      if (mype == 0) print *, "Unknown version, trying old snapshot reader..."
      call MPI_FILE_CLOSE(file_handle,ierrmpi)
      call read_snapshot_old()

      ! Guess snapshotnext from file name if not set
      if (snapshotnext == -1) snapshotnext = &
         get_snapshot_index(trim(restart_from_file)) + 1
      ! Set slicenext and collapsenext if not set
      if (slicenext == -1) slicenext = 0
      if (collapsenext == -1) collapsenext = 0

      ! Still used in convert
      snapshotini = snapshotnext-1

      return ! Leave this routine
    else if (mype == 0) then
      call MPI_FILE_SEEK(file_handle, 0_MPI_OFFSET_KIND, MPI_SEEK_SET,&
          ierrmpi)
      call snapshot_read_header(file_handle, offset_tree_info,&
          offset_block_data)
    end if

    ! Share information about restart file
    call MPI_BCAST(nw_found,1,MPI_INTEGER,0,icomm,ierrmpi)
    call MPI_BCAST(nleafs,1,MPI_INTEGER,0,icomm,ierrmpi)
    call MPI_BCAST(nparents,1,MPI_INTEGER,0,icomm,ierrmpi)
    call MPI_BCAST(it,1,MPI_INTEGER,0,icomm,ierrmpi)
    call MPI_BCAST(global_time,1,MPI_DOUBLE_PRECISION,0,icomm,ierrmpi)

    call MPI_BCAST(snapshotnext,1,MPI_INTEGER,0,icomm,ierrmpi)
    call MPI_BCAST(slicenext,1,MPI_INTEGER,0,icomm,ierrmpi)
    call MPI_BCAST(collapsenext,1,MPI_INTEGER,0,icomm,ierrmpi)
    call MPI_BCAST(stagger_mark_dat,1,MPI_LOGICAL,0,icomm,ierrmpi)

    ! Allocate send/receive buffer
    n_values = count_ix(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3) * nw_found
    if(stagger_mark_dat) then
      n_values = n_values + count_ix(ixGslo1,ixGslo2,ixGslo3,ixGshi1,ixGshi2,&
         ixGshi3) * nws
    end if
    allocate(w_buffer(n_values))
    allocate(w(ixGlo1:ixGhi1,ixGlo2:ixGhi2,ixGlo3:ixGhi3,1:nw_found))

    nleafs_active = nleafs

    if (mype == 0) then
      call MPI_FILE_SEEK(file_handle, offset_tree_info, MPI_SEEK_SET, ierrmpi)
    end if

    call read_forest(file_handle)

    do Morton_no=Morton_start(mype),Morton_stop(mype)
      igrid=sfc_to_igrid(Morton_no)
      call alloc_node(igrid)
    end do

    if (mype==0) then
      call MPI_FILE_SEEK(file_handle, offset_block_data, MPI_SEEK_SET,&
          ierrmpi)

      iread = 0
      do ipe = 0, npe-1
        do Morton_no=Morton_start(ipe),Morton_stop(ipe)
          iread=iread+1
          itag=Morton_no

          call mpi_file_read_wrapper(file_handle,ix_buffer(1:2*ndim), 2*ndim,&
              MPI_INTEGER, istatus,ierrmpi)

          ! Construct ixO^L array from number of ghost cells
          ixOmin1 = ixMlo1 - ix_buffer(1)
          ixOmin2 = ixMlo2 - ix_buffer(2)
          ixOmin3 = ixMlo3 - ix_buffer(3)
          ixOmax1 = ixMhi1 + ix_buffer(ndim+1)
          ixOmax2 = ixMhi2 + ix_buffer(ndim+2)
          ixOmax3 = ixMhi3 + ix_buffer(ndim+3)
          n_values = count_ix(ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,&
             ixOmax3) * nw_found
          if(stagger_mark_dat) then
            ixOsmin1 = ixOmin1 - 1
            ixOsmin2 = ixOmin2 - 1
            ixOsmin3 = ixOmin3 - 1
            ixOsmax1 = ixOmax1
            ixOsmax2 = ixOmax2
            ixOsmax3 = ixOmax3
            n_values_stagger = n_values
            n_values = n_values + count_ix(ixOsmin1,ixOsmin2,ixOsmin3,ixOsmax1,&
               ixOsmax2,ixOsmax3) * nws
          end if

          call mpi_file_read_wrapper(file_handle, w_buffer, n_values,&
              MPI_DOUBLE_PRECISION, istatus, ierrmpi)

          if (mype == ipe) then ! Root task
            igrid=sfc_to_igrid(Morton_no)
            block=>ps(igrid)
            if(stagger_mark_dat) then
              w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                  1:nw_found) = reshape(w_buffer(1:n_values_stagger),&
                  shape(w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                  1:nw_found)))
              if(stagger_grid) ps(igrid)%ws(ixOsmin1:ixOsmax1,&
                 ixOsmin2:ixOsmax2,ixOsmin3:ixOsmax3,&
                 1:nws)=reshape(w_buffer(n_values_stagger+1:n_values),&
                  shape(ws(ixOsmin1:ixOsmax1,ixOsmin2:ixOsmax2,&
                 ixOsmin3:ixOsmax3, 1:nws)))
            else
              w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                  1:nw_found) = reshape(w_buffer(1:n_values),&
                  shape(w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                  1:nw_found)))
            end if
            if (nw_found<nw) then
              if (associated(usr_transform_w)) then
                call usr_transform_w(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
                   ixMlo1,ixMlo2,ixMlo3,ixMhi1,ixMhi2,ixMhi3,nw_found,w,&
                   ps(igrid)%x,ps(igrid)%w)
              else
                ps(igrid)%w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                   1:nw_found)=w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
                   ixOmin3:ixOmax3,1:nw_found)
              end if
            else if (nw_found>nw) then
              if (associated(usr_transform_w)) then
                call usr_transform_w(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
                   ixMlo1,ixMlo2,ixMlo3,ixMhi1,ixMhi2,ixMhi3,nw_found,w,&
                   ps(igrid)%x,ps(igrid)%w)
              else
                ps(igrid)%w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                   1:nw)=w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                   1:nw)
              end if
            else
              ps(igrid)%w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
                 1:nw)=w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,1:nw)
            end if
            
          else
            call mpi_send_wrapper([ ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,&
                n_values ], 2*ndim+1, MPI_INTEGER, ipe, itag, icomm, ierrmpi)
            call mpi_send_wrapper(w_buffer, n_values, MPI_DOUBLE_PRECISION, ipe, itag,&
                icomm, ierrmpi)
          end if
        end do
      end do

      call MPI_FILE_CLOSE(file_handle,ierrmpi)

    else ! mype > 0

      do Morton_no=Morton_start(mype),Morton_stop(mype)
        igrid=sfc_to_igrid(Morton_no)
        block=>ps(igrid)
        itag=Morton_no

        call mpi_recv_wrapper(ix_buffer, 2*ndim+1, MPI_INTEGER, 0, itag, icomm,&
           iorecvstatus, ierrmpi)
        ixOmin1 = ix_buffer(1)
        ixOmin2 = ix_buffer(2)
        ixOmin3 = ix_buffer(3)
        ixOmax1 = ix_buffer(ndim+1)
        ixOmax2 = ix_buffer(ndim+2)
        ixOmax3 = ix_buffer(ndim+3)
        n_values = ix_buffer(2*ndim+1)

        call mpi_recv_wrapper(w_buffer, n_values, MPI_DOUBLE_PRECISION,0, itag, icomm,&
            iorecvstatus, ierrmpi)

        if(stagger_mark_dat) then
          n_values_stagger = count_ix(ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,&
             ixOmax3) * nw_found
          ixOsmin1 = ixOmin1 - 1
          ixOsmin2 = ixOmin2 - 1
          ixOsmin3 = ixOmin3 - 1
          ixOsmax1 = ixOmax1
          ixOsmax2 = ixOmax2
          ixOsmax3 = ixOmax3
          w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
              1:nw_found) = reshape(w_buffer(1:n_values_stagger),&
              shape(w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
              1:nw_found)))
          if(stagger_grid) ps(igrid)%ws(ixOsmin1:ixOsmax1,ixOsmin2:ixOsmax2,&
             ixOsmin3:ixOsmax3,1:nws)=reshape(w_buffer(n_values_stagger+&
             1:n_values), shape(ws(ixOsmin1:ixOsmax1,ixOsmin2:ixOsmax2,&
             ixOsmin3:ixOsmax3, 1:nws)))
        else
          w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
              1:nw_found) = reshape(w_buffer(1:n_values),&
              shape(w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
              1:nw_found)))
        end if
        if (nw_found<nw) then
          if (associated(usr_transform_w)) then
            call usr_transform_w(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
               ixMlo1,ixMlo2,ixMlo3,ixMhi1,ixMhi2,ixMhi3,nw_found,w,&
               ps(igrid)%x,ps(igrid)%w)
          else
            ps(igrid)%w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
               1:nw_found)=w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
               1:nw_found)
          end if
        else if (nw_found>nw) then
          if (associated(usr_transform_w)) then
            call usr_transform_w(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
               ixMlo1,ixMlo2,ixMlo3,ixMhi1,ixMhi2,ixMhi3,nw_found,w,&
               ps(igrid)%x,ps(igrid)%w)
          else
            ps(igrid)%w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
               1:nw)=w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,1:nw)
          end if
        else
          ps(igrid)%w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,&
             1:nw)=w(ixOmin1:ixOmax1,ixOmin2:ixOmax2,ixOmin3:ixOmax3,1:nw)
        end if
        
      end do
    end if

    call MPI_BARRIER(icomm,ierrmpi)

  end subroutine read_snapshot