write_snapshot Subroutine

public subroutine write_snapshot()

Uses

  • proc~~write_snapshot~~UsesGraph proc~write_snapshot write_snapshot mod_physics mod_physics proc~write_snapshot->mod_physics module~mod_forest mod_forest proc~write_snapshot->module~mod_forest module~mod_functions_forest mod_functions_forest proc~write_snapshot->module~mod_functions_forest module~mod_global_parameters mod_global_parameters proc~write_snapshot->module~mod_global_parameters module~mod_input_output_helper mod_input_output_helper proc~write_snapshot->module~mod_input_output_helper 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

Arguments

None

Calls

proc~~write_snapshot~~CallsGraph proc~write_snapshot write_snapshot mpi_barrier mpi_barrier proc~write_snapshot->mpi_barrier mpi_file_close mpi_file_close proc~write_snapshot->mpi_file_close mpi_file_get_position mpi_file_get_position proc~write_snapshot->mpi_file_get_position mpi_file_seek mpi_file_seek proc~write_snapshot->mpi_file_seek mpi_file_write mpi_file_write proc~write_snapshot->mpi_file_write proc~block_shape_io block_shape_io proc~write_snapshot->proc~block_shape_io proc~count_ix count_ix proc~write_snapshot->proc~count_ix proc~create_output_file create_output_file proc~write_snapshot->proc~create_output_file proc~mpi_file_write_wrapper mpi_file_write_wrapper proc~write_snapshot->proc~mpi_file_write_wrapper proc~mpi_recv_wrapper mpi_recv_wrapper proc~write_snapshot->proc~mpi_recv_wrapper proc~mpi_send_wrapper mpi_send_wrapper proc~write_snapshot->proc~mpi_send_wrapper proc~snapshot_write_header snapshot_write_header proc~write_snapshot->proc~snapshot_write_header proc~write_forest write_forest proc~write_snapshot->proc~write_forest proc~block_shape_io->proc~count_ix mpi_file_open mpi_file_open proc~create_output_file->mpi_file_open proc~mpistop mpistop proc~create_output_file->proc~mpistop proc~mpi_file_write_wrapper->mpi_file_write mpi_recv mpi_recv proc~mpi_recv_wrapper->mpi_recv mpi_send mpi_send proc~mpi_send_wrapper->mpi_send proc~snapshot_write_header1 snapshot_write_header1 proc~snapshot_write_header->proc~snapshot_write_header1 proc~write_forest->mpi_file_write mpi_abort mpi_abort proc~mpistop->mpi_abort proc~snapshot_write_header1->mpi_file_write phys_write_info phys_write_info proc~snapshot_write_header1->phys_write_info

Called by

proc~~write_snapshot~~CalledByGraph proc~write_snapshot write_snapshot proc~saveamrfile saveamrfile proc~saveamrfile->proc~write_snapshot proc~timeintegration timeintegration proc~timeintegration->proc~saveamrfile proc~main main proc~main->proc~timeintegration program~agile agile program~agile->proc~main

Source Code

  subroutine write_snapshot
    use mod_forest
    use mod_global_parameters
    use mod_physics
    use mod_input_output_helper, only: count_ix,block_shape_io,&
       create_output_file
    use mod_functions_forest, only: write_forest

    integer                       :: file_handle, igrid, Morton_no, iwrite
    integer                       :: ipe, ix_buffer(2*ndim+1), n_values
    integer                       :: ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,&
       ixOmax3, n_ghost(2*ndim)
    integer                       :: ixOsmin1,ixOsmin2,ixOsmin3,ixOsmax1,&
       ixOsmax2,ixOsmax3,n_values_stagger
    integer                       :: iorecvstatus(MPI_STATUS_SIZE)
    integer                       :: ioastatus(MPI_STATUS_SIZE)
    integer                       :: igrecvstatus(MPI_STATUS_SIZE)
    integer                       :: istatus(MPI_STATUS_SIZE)
    type(tree_node), pointer      :: pnode
    integer(kind=MPI_OFFSET_KIND) :: offset_tree_info
    integer(kind=MPI_OFFSET_KIND) :: offset_block_data
    integer(kind=MPI_OFFSET_KIND) :: offset_offsets
    double precision, allocatable :: w_buffer(:)

    integer, allocatable                       :: block_ig(:, :)
    integer, allocatable                       :: block_lvl(:)
    integer(kind=MPI_OFFSET_KIND), allocatable :: block_offset(:)

    call MPI_BARRIER(icomm, ierrmpi)

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

    ! Allocate arrays with information about grid blocks
    allocate(block_ig(ndim, nleafs))
    allocate(block_lvl(nleafs))
    allocate(block_offset(nleafs+1))

    ! master processor
    if (mype==0) then
      call create_output_file(file_handle, snapshotnext, ".dat")

      ! Don't know offsets yet, we will write header again later
      offset_tree_info = -1
      offset_block_data = -1
      call snapshot_write_header(file_handle, offset_tree_info,&
          offset_block_data)

      call MPI_File_get_position(file_handle, offset_tree_info, ierrmpi)

      call write_forest(file_handle)

      ! Collect information about the spatial index (ig^D) and refinement level
      ! of leaves
      do Morton_no = Morton_start(0), Morton_stop(npe-1)
         igrid = sfc(1, Morton_no)
         ipe = sfc(2, Morton_no)
         pnode => igrid_to_node(igrid, ipe)%node

         block_ig(:, Morton_no) = [ pnode%ig1,pnode%ig2,pnode%ig3 ]
         block_lvl(Morton_no) = pnode%level
         block_offset(Morton_no) = 0 ! Will be determined later
      end do

      call MPI_FILE_WRITE(file_handle, block_lvl, size(block_lvl), MPI_INTEGER,&
          istatus, ierrmpi)

      call MPI_FILE_WRITE(file_handle, block_ig, size(block_ig), MPI_INTEGER,&
          istatus, ierrmpi)

      ! Block offsets are currently unknown, but will be overwritten later
      call MPI_File_get_position(file_handle, offset_offsets, ierrmpi)
      call mpi_file_write_wrapper(file_handle, block_offset(1:nleafs), nleafs,&
          MPI_OFFSET, istatus, ierrmpi)

      call MPI_File_get_position(file_handle, offset_block_data, ierrmpi)

      ! Check whether data was written as expected
      if (offset_block_data - offset_tree_info /= (nleafs + nparents) * &
         size_logical + nleafs * ((1+ndim) * size_int + 2 * size_int)) then
        if (mype == 0) then
          print *, "Warning: MPI_OFFSET type /= 8 bytes"
          print *, "This *could* cause problems when reading .dat files"
        end if
      end if

      block_offset(1) = offset_block_data
      iwrite = 0
    end if

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

      call block_shape_io(igrid, n_ghost, ixOmin1,ixOmin2,ixOmin3,ixOmax1,&
         ixOmax2,ixOmax3, n_values)
      if(stagger_grid) then
        w_buffer(1:n_values) = pack(ps(igrid)%w(ixOmin1:ixOmax1,&
           ixOmin2:ixOmax2,ixOmin3:ixOmax3, 1:nw), .true.)
        ixOsmin1 = ixOmin1 -1
        ixOsmin2 = ixOmin2 -1
        ixOsmin3 = ixOmin3 -1
        ixOsmax1 = ixOmax1
        ixOsmax2 = ixOmax2
        ixOsmax3 = ixOmax3
        n_values_stagger= count_ix(ixOsmin1,ixOsmin2,ixOsmin3,ixOsmax1,&
           ixOsmax2,ixOsmax3)*nws
        w_buffer(n_values+1:n_values+n_values_stagger) = &
           pack(ps(igrid)%ws(ixOsmin1:ixOsmax1,ixOsmin2:ixOsmax2,&
           ixOsmin3:ixOsmax3, 1:nws), .true.)
        n_values=n_values+n_values_stagger
      else
        w_buffer(1:n_values) = pack(ps(igrid)%w(ixOmin1:ixOmax1,&
           ixOmin2:ixOmax2,ixOmin3:ixOmax3, 1:nw), .true.)
      end if
      ix_buffer(1) = n_values
      ix_buffer(2:) = n_ghost

      if (mype /= 0) then
        call mpi_send_wrapper(ix_buffer, 2*ndim+1, MPI_INTEGER, 0, itag, icomm,&
            ierrmpi)
        call mpi_send_wrapper(w_buffer, n_values, MPI_DOUBLE_PRECISION, 0, itag, icomm,&
            ierrmpi)
      else
        iwrite = iwrite+1
        call mpi_file_write_wrapper(file_handle, ix_buffer(2:), 2*ndim, MPI_INTEGER,&
            istatus, ierrmpi)
        call mpi_file_write_wrapper(file_handle, w_buffer, n_values,&
            MPI_DOUBLE_PRECISION, istatus, ierrmpi)

        ! Set offset of next block
        block_offset(iwrite+1) = block_offset(iwrite) + int(n_values,&
            MPI_OFFSET_KIND) * size_double + 2 * ndim * size_int
      end if
    end do

    ! Write data communicated from other processors
    if (mype == 0) then
      do ipe = 1, npe-1
        do Morton_no=Morton_start(ipe), Morton_stop(ipe)
          iwrite=iwrite+1
          itag=Morton_no

          call mpi_recv_wrapper(ix_buffer, 2*ndim+1, MPI_INTEGER, ipe, itag, icomm,&
             igrecvstatus, ierrmpi)
          n_values = ix_buffer(1)

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

          call mpi_file_write_wrapper(file_handle, ix_buffer(2:), 2*ndim, MPI_INTEGER,&
              istatus, ierrmpi)
          call mpi_file_write_wrapper(file_handle, w_buffer, n_values,&
              MPI_DOUBLE_PRECISION, istatus, ierrmpi)

          ! Set offset of next block
          block_offset(iwrite+1) = block_offset(iwrite) + int(n_values,&
              MPI_OFFSET_KIND) * size_double + 2 * ndim * size_int
        end do
      end do

      ! Write block offsets (now we know them)
      call MPI_FILE_SEEK(file_handle, offset_offsets, MPI_SEEK_SET, ierrmpi)
      call mpi_file_write_wrapper(file_handle, block_offset(1:nleafs), nleafs,&
          MPI_OFFSET, istatus, ierrmpi)

      ! Write header again, now with correct offsets
      call MPI_FILE_SEEK(file_handle, 0_MPI_OFFSET_KIND, MPI_SEEK_SET,&
          ierrmpi)
      call snapshot_write_header(file_handle, offset_tree_info,&
          offset_block_data)

      call MPI_FILE_CLOSE(file_handle, ierrmpi)
    end if

    call MPI_BARRIER(icomm, ierrmpi)

  end subroutine write_snapshot