subroutine convert_dat_generic(nwc, dataset_names, file_suffix,&
convert_vars)
use mod_forest
use mod_global_parameters
use mod_input_output_helper, only: count_ix, create_output_file,&
snapshot_write_header1, block_shape_io
use mod_functions_forest, only: write_forest
integer :: file_handle, igrid, Morton_no, iwrite
integer :: ipe, ix_buffer(2*ndim+1), n_values
integer :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,&
ixImax3, 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(:)
double precision, allocatable:: converted_vars(:,:,:,:)
integer :: idim, itag
integer, allocatable :: block_ig(:, :)
integer, allocatable :: block_lvl(:)
integer(kind=MPI_OFFSET_KIND), allocatable :: block_offset(:)
integer, intent(in) :: nwc
character(len=*), intent(in) :: dataset_names(:)
character(len=*), intent(in) :: file_suffix
interface
function convert_vars(ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,w,x,&
nwc) result(wres)
use mod_global_parameters
integer, intent(in) :: ixImin1,ixImin2,ixImin3,ixImax1,ixImax2,ixImax3,&
ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3, nwc
double precision, intent(in) :: x(ixImin1:ixImax1,ixImin2:ixImax2,&
ixImin3:ixImax3,1:ndim)
double precision, intent(in) :: w(ixImin1:ixImax1,ixImin2:ixImax2,&
ixImin3:ixImax3,1:nw)
double precision :: wres(ixOmin1:ixOmax1,ixOmin2:ixOmax2,&
ixOmin3:ixOmax3,1:nwc)
end function convert_vars
end interface
call MPI_BARRIER(icomm, ierrmpi)
! Allocate send/receive buffer
n_values = count_ix(ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3) * nw
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", file_suffix)
! Don't know offsets yet, we will write header again later
offset_tree_info = -1
offset_block_data = -1
call snapshot_write_header1(file_handle, offset_tree_info,&
offset_block_data, dataset_names, nwc)
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(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
block=>ps(igrid)
! this might be used in convert function,
! it was not used when the output is already computed vars (write_snapshot)
dxlevel(1)=rnode(rpdx1_,igrid);dxlevel(2)=rnode(rpdx2_,igrid)
dxlevel(3)=rnode(rpdx3_,igrid);
! start copied from block_shape_io,
! because nwc is needed as parameter
! TODO check if this will be used elsewehere and put it in separate subroutine
n_ghost(:) = 0
if(save_physical_boundary) then
do idim=1,ndim
! Include ghost cells on lower boundary
if(ps(igrid)%is_physical_boundary(2*idim-1)) &
n_ghost(idim)=nghostcells
! Include ghost cells on upper boundary
if(ps(igrid)%is_physical_boundary(2*idim)) &
n_ghost(ndim+idim)=nghostcells
end do
end if
ixOmin1 = ixMlo1 - n_ghost(1)
ixOmin2 = ixMlo2 - n_ghost(2)
ixOmin3 = ixMlo3 - n_ghost(3)
ixOmax1 = ixMhi1 + n_ghost(ndim+1)
ixOmax2 = ixMhi2 + n_ghost(ndim+2)
ixOmax3 = ixMhi3 + n_ghost(ndim+3)
n_values = count_ix(ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,&
ixOmax3) * nwc
! end copied from block_shape_io
ixImin1 = ixGlo1
ixImin2 = ixGlo2
ixImin3 = ixGlo3
ixImax1 = ixGhi1
ixImax2 = ixGhi2
ixImax3 = ixGhi3
w_buffer(1:n_values) = pack(convert_vars(ixImin1,ixImin2,ixImin3,ixImax1,&
ixImax2,ixImax3, ixOmin1,ixOmin2,ixOmin3,ixOmax1,ixOmax2,ixOmax3,&
ps(igrid)%w(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3, 1:nw),&
ps(igrid)%x(ixImin1:ixImax1,ixImin2:ixImax2,ixImin3:ixImax3, 1:ndim),&
nwc), .true.)
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(file_handle, ix_buffer(2:), 2*ndim, MPI_INTEGER,&
istatus, ierrmpi)
call MPI_FILE_WRITE(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(file_handle, ix_buffer(2:), 2*ndim, MPI_INTEGER,&
istatus, ierrmpi)
call MPI_FILE_WRITE(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(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_header1(file_handle, offset_tree_info,&
offset_block_data, dataset_names, nwc)
call MPI_FILE_CLOSE(file_handle, ierrmpi)
end if
call MPI_BARRIER(icomm, ierrmpi)
end subroutine convert_dat_generic