subroutine ImageDataVtk_mpi(qunit)
! output for vti format to paraview, non-binary version output
! parallel, uses calc_grid to compute nwauxio variables
! allows renormalizing using convert factors
! allows skipping of w_write selected variables
! implementation such that length of ASCII output is identical when
! run on 1 versus multiple CPUs (however, the order of the vtu pieces can differ)
use mod_forest, only: Morton_start, Morton_stop, tree_node_ptr,&
igrid_to_node, sfc_to_igrid
use mod_global_parameters
use mod_calculate_xw
integer, intent(in) :: qunit
double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,&
ixMlo3-1:ixMhi3,ndim) :: xC_TMP,xC_TMP_recv
double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo3:ixMhi3,&
ndim) :: xCC_TMP,xCC_TMP_recv
double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,&
ixMlo3-1:ixMhi3,ndim) :: xC
double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo3:ixMhi3,&
ndim) :: xCC
double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,&
ixMlo3-1:ixMhi3,nw+nwauxio) :: wC_TMP,wC_TMP_recv
double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo3:ixMhi3,&
nw+nwauxio) :: wCC_TMP,wCC_TMP_recv
double precision, dimension(0:nw+nwauxio) :: normconv
integer:: igrid,iigrid,level,ixCmin1,ixCmin2,ixCmin3,ixCmax1,&
ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3
integer:: NumGridsOnLevel(1:nlevelshi)
integer :: nx1,nx2,nx3
character(len=80):: filename
integer :: filenr
integer, allocatable :: intstatus(:,:)
logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
character(len=name_len) :: wnamei(1:nw+nwauxio),&
xandwnamei(1:ndim+nw+nwauxio)
character(len=1024) :: outfilehead
logical :: fileopen
integer :: itag,ipe,Morton_no,Morton_length
integer :: ixrvCmin1,ixrvCmin2,ixrvCmin3,ixrvCmax1,ixrvCmax2,ixrvCmax3,&
ixrvCCmin1,ixrvCCmin2,ixrvCCmin3,ixrvCCmax1,ixrvCCmax2,ixrvCCmax3,&
siz_ind, ind_send(5*3), ind_recv(5*3)
double precision :: origin(1:3), spacing(1:3)
integer :: wholeExtent(1:6), ig1,ig2,ig3
type(tree_node_ptr) :: tree
if(levmin/=levmax) call mpistop&
('ImageData can only be used when levmin=levmax')
normconv(0) = length_convert_factor
normconv(1:nw) = w_convert_factor
siz_ind=5*3
Morton_length=Morton_stop(npe-1)-Morton_start(0)+1
allocate(Morton_aim(Morton_start(0):Morton_stop(npe-1)))
allocate(Morton_aim_p(Morton_start(0):Morton_stop(npe-1)))
Morton_aim=.false.
Morton_aim_p=.false.
do Morton_no=Morton_start(mype),Morton_stop(mype)
igrid=sfc_to_igrid(Morton_no)
level=node(plevel_,igrid)
! we can clip parts of the grid away, select variables, levels etc.
if(writelevel(level)) then
! only output a grid when fully within clipped region selected
! by writespshift array
if((rnode(rpxmin1_,igrid)>=xprobmin1+&
(xprobmax1-xprobmin1)*writespshift(1,1).and.rnode(rpxmin2_,&
igrid)>=xprobmin2+(xprobmax2-xprobmin2)*writespshift(2,&
1).and.rnode(rpxmin3_,igrid)>=xprobmin3+&
(xprobmax3-xprobmin3)*writespshift(3,1)).and.(rnode(rpxmax1_,&
igrid)<=xprobmax1-(xprobmax1-xprobmin1)*writespshift(1,&
2).and.rnode(rpxmax2_,igrid)<=xprobmax2-(xprobmax2-&
xprobmin2)*writespshift(2,2).and.rnode(rpxmax3_,&
igrid)<=xprobmax3-(xprobmax3-xprobmin3)*writespshift(3,2))) then
Morton_aim_p(Morton_no)=.true.
end if
end if
end do
call MPI_ALLREDUCE(Morton_aim_p,Morton_aim,Morton_length,MPI_LOGICAL,&
MPI_LOR,icomm,ierrmpi)
if(mype /= 0) then
do Morton_no=Morton_start(mype),Morton_stop(mype)
if(.not. Morton_aim(Morton_no)) cycle
igrid=sfc_to_igrid(Morton_no)
call calc_x(igrid,xC,xCC)
call calc_grid(qunit,igrid,xC,xCC,xC_TMP,xCC_TMP,wC_TMP,wCC_TMP,&
normconv,ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,&
ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3,.true.)
tree%node => igrid_to_node(igrid, mype)%node
ig1 = tree%node%ig1; ig2 = tree%node%ig2; ig3 = tree%node%ig3;
itag=Morton_no
ind_send=(/ ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,&
ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3, ig1,ig2,ig3 /)
call mpi_send_wrapper(ind_send,siz_ind,MPI_INTEGER, 0,itag,icomm,ierrmpi)
call mpi_send_wrapper(wC_TMP,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
call mpi_send_wrapper(wCC_TMP,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
end do
else
inquire(qunit,opened=fileopen)
if(.not.fileopen)then
! generate filename
filenr=snapshotini
if (autoconvert) filenr=snapshotnext
write(filename,'(a,i4.4,a)') TRIM(base_filename),filenr,".vti"
! Open the file for the header part
open(qunit,file=filename,status='unknown',form='formatted')
end if
call getheadernames(wnamei,xandwnamei,outfilehead)
! number of cells per grid.
nx1=ixMhi1-ixMlo1+1;nx2=ixMhi2-ixMlo2+1;nx3=ixMhi3-ixMlo3+1;
origin = 0
origin(1) = xprobmin1*normconv(0); origin(2) = xprobmin2*normconv(0)
origin(3) = xprobmin3*normconv(0);
spacing = zero
spacing(1) = dxlevel(1)*normconv(0)
spacing(2) = dxlevel(2)*normconv(0)
spacing(3) = dxlevel(3)*normconv(0);
wholeExtent = 0
! if we use writespshift, the whole extent has to be calculated:
wholeExtent(1*2-1) = nx1 * ceiling(((xprobmax1-xprobmin1)*writespshift(1,&
1)) /(nx1*dxlevel(1)))
wholeExtent(2*2-1) = nx2 * ceiling(((xprobmax2-xprobmin2)*writespshift(2,&
1)) /(nx2*dxlevel(2)))
wholeExtent(3*2-1) = nx3 * ceiling(((xprobmax3-xprobmin3)*writespshift(3,&
1)) /(nx3*dxlevel(3)))
wholeExtent(1*2) = nx1 * floor(((xprobmax1-xprobmin1)*(1.0d0-&
writespshift(1,2))) /(nx1*dxlevel(1)))
wholeExtent(2*2) = nx2 * floor(((xprobmax2-xprobmin2)*(1.0d0-&
writespshift(2,2))) /(nx2*dxlevel(2)))
wholeExtent(3*2) = nx3 * floor(((xprobmax3-xprobmin3)*(1.0d0-&
writespshift(3,2))) /(nx3*dxlevel(3)))
! generate xml header
write(qunit,'(a)')'<?xml version="1.0"?>'
write(qunit,'(a)',advance='no') '<VTKFile type="ImageData"'
write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
write(qunit,'(a,3(1pe14.6),a,6(i10),a,3(1pe14.6),a)')&
' <ImageData Origin="',origin,'" WholeExtent="',wholeExtent,&
'" Spacing="',spacing,'">'
write(qunit,'(a)')'<FieldData>'
write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
'NumberOfTuples="1" format="ascii">'
write(qunit,*) real(global_time*time_convert_factor)
write(qunit,'(a)')'</DataArray>'
write(qunit,'(a)')'</FieldData>'
! write the data from proc 0
do Morton_no=Morton_start(0),Morton_stop(0)
if(.not. Morton_aim(Morton_no)) cycle
igrid=sfc_to_igrid(Morton_no)
tree%node => igrid_to_node(igrid, 0)%node
ig1 = tree%node%ig1; ig2 = tree%node%ig2; ig3 = tree%node%ig3;
call calc_x(igrid,xC,xCC)
call calc_grid(qunit,igrid,xC,xCC,xC_TMP,xCC_TMP,wC_TMP,wCC_TMP,&
normconv,ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,&
ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3,.true.)
call write_vti(qunit,ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,ixCmin1,&
ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,&
ixCCmax1,ixCCmax2,ixCCmax3,ig1,ig2,ig3,nx1,nx2,nx3,normconv,wnamei,&
wC_TMP,wCC_TMP)
end do
if(npe>1)then
allocate(intstatus(MPI_STATUS_SIZE,1))
do ipe=1, npe-1
do Morton_no=Morton_start(ipe),Morton_stop(ipe)
if(.not. Morton_aim(Morton_no)) cycle
itag=Morton_no
call mpi_recv_wrapper(ind_recv,siz_ind, MPI_INTEGER, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
ixrvCmin1=ind_recv(1);ixrvCmin2=ind_recv(2);ixrvCmin3=ind_recv(3)
ixrvCmax1=ind_recv(3+1);ixrvCmax2=ind_recv(3+2)
ixrvCmax3=ind_recv(3+3);
ixrvCCmin1=ind_recv(2*3+1);ixrvCCmin2=ind_recv(2*3+2)
ixrvCCmin3=ind_recv(2*3+3);ixrvCCmax1=ind_recv(3*3+1)
ixrvCCmax2=ind_recv(3*3+2);ixrvCCmax3=ind_recv(3*3+3);
ig1=ind_recv(4*3+1);ig2=ind_recv(4*3+2);ig3=ind_recv(4*3+3);
call mpi_recv_wrapper(wC_TMP,1,type_block_wc_io, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
call mpi_recv_wrapper(wCC_TMP,1,type_block_wcc_io, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
call write_vti(qunit,ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
ixrvCmin1,ixrvCmin2,ixrvCmin3,ixrvCmax1,ixrvCmax2,ixrvCmax3,&
ixrvCCmin1,ixrvCCmin2,ixrvCCmin3,ixrvCCmax1,ixrvCCmax2,&
ixrvCCmax3,ig1,ig2,ig3,nx1,nx2,nx3,normconv,wnamei,wC_TMP,&
wCC_TMP)
end do
end do
end if
write(qunit,'(a)')'</ImageData>'
write(qunit,'(a)')'</VTKFile>'
close(qunit)
if(npe>1) deallocate(intstatus)
end if
deallocate(Morton_aim,Morton_aim_p)
if (npe>1) then
call MPI_BARRIER(icomm,ierrmpi)
endif
end subroutine ImageDataVtk_mpi