subroutine unstructuredvtk_mpi(qunit)
! output for vtu 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, sfc_to_igrid
use mod_global_parameters
use mod_calculate_xw
integer, intent(in) :: qunit
double precision :: x_VTK(1:3)
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,nxC1,nxC2,nxC3,nodesonlevel,elemsonlevel,nc,np,ix1,&
ix2,ix3
character(len=80):: filename
integer :: filenr
integer, allocatable :: intstatus(:,:)
character(len=name_len) :: wnamei(1:nw+nwauxio),&
xandwnamei(1:ndim+nw+nwauxio)
character(len=1024) :: outfilehead
logical :: fileopen,conv_grid,cond_grid_recv
integer :: itag,ipe,Morton_no,siz_ind
integer :: ind_send(4*3),ind_recv(4*3)
integer :: levmin_recv,levmax_recv,level_recv,igrid_recv,ixrvCmin1,&
ixrvCmin2,ixrvCmin3,ixrvCmax1,ixrvCmax2,ixrvCmax3,ixrvCCmin1,ixrvCCmin2,&
ixrvCCmin3,ixrvCCmax1,ixrvCCmax2,ixrvCCmax3
if(mype==0) then
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,".vtu"
! Open the file for the header part
open(qunit,file=filename,status='unknown',form='formatted')
end if
! generate xml header
write(qunit,'(a)')'<?xml version="1.0"?>'
write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
write(qunit,'(a)')'<UnstructuredGrid>'
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>'
end if
call getheadernames(wnamei,xandwnamei,outfilehead)
! number of cells, number of corner points, per grid.
nx1=ixMhi1-ixMlo1+1;nx2=ixMhi2-ixMlo2+1;nx3=ixMhi3-ixMlo3+1;
nxC1=nx1+1;nxC2=nx2+1;nxC3=nx3+1;
nc=nx1*nx2*nx3
np=nxC1*nxC2*nxC3
! all slave processors send their minmal/maximal levels
if (mype/=0) then
if (Morton_stop(mype)==0) call mpistop("nultag")
itag=1000*Morton_stop(mype)
!print *,'ype,itag for levmin=',mype,itag,levmin
call mpi_send_wrapper(levmin,1,MPI_INTEGER, 0,itag,icomm,ierrmpi)
itag=2000*Morton_stop(mype)
!print *,'mype,itag for levmax=',mype,itag,levmax
call mpi_send_wrapper(levmax,1,MPI_INTEGER, 0,itag,icomm,ierrmpi)
end if
! Note: using the w_write, writelevel, writespshift
! we can clip parts of the grid away, select variables, levels etc.
do level=levmin,levmax
if (.not.writelevel(level)) cycle
do Morton_no=Morton_start(mype),Morton_stop(mype)
igrid=sfc_to_igrid(Morton_no)
if (mype/=0)then
itag=Morton_no
call mpi_send_wrapper(igrid,1,MPI_INTEGER, 0,itag,icomm,ierrmpi)
itag=igrid
call mpi_send_wrapper(node(plevel_,igrid),1,MPI_INTEGER, 0,itag,icomm,&
ierrmpi)
end if
if (node(plevel_,igrid)/=level) cycle
! only output a grid when fully within clipped region selected
! by writespshift array
conv_grid=(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))
if (mype/=0)then
call mpi_send_wrapper(conv_grid,1,MPI_LOGICAL,0,itag,icomm,ierrmpi)
end if
if (.not.conv_grid) cycle
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.)
if(mype/=0) then
itag=Morton_no
ind_send=(/ ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,&
ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3 /)
siz_ind=4*3
call mpi_send_wrapper(ind_send,siz_ind,MPI_INTEGER, 0,itag,icomm,ierrmpi)
call mpi_send_wrapper(normconv,nw+nwauxio+1,MPI_DOUBLE_PRECISION, 0,itag,&
icomm,ierrmpi)
call mpi_send_wrapper(wC_TMP,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
call mpi_send_wrapper(xC_TMP,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
itag=igrid
call mpi_send_wrapper(wCC_TMP,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
call mpi_send_wrapper(xCC_TMP,1,type_block_xcc_io, 0,itag,icomm,ierrmpi)
else
call write_vtk(qunit,ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,&
ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3,igrid,nc,np,nx1,nx2,nx3,nxC1,&
nxC2,nxC3,normconv,wnamei,xC_TMP,xCC_TMP,wC_TMP,wCC_TMP)
end if
end do ! Morton_no loop
end do ! level loop
if(mype==0) then
allocate(intstatus(MPI_STATUS_SIZE,1))
if(npe>1)then
do ipe=1,npe-1
itag=1000*Morton_stop(ipe)
call mpi_recv_wrapper(levmin_recv,1,MPI_INTEGER, ipe,itag,icomm,intstatus(:,&
1),ierrmpi)
!!print *,'mype RECEIVES,itag for levmin=',mype,itag,levmin_recv
itag=2000*Morton_stop(ipe)
call mpi_recv_wrapper(levmax_recv,1,MPI_INTEGER, ipe,itag,icomm,intstatus(:,&
1),ierrmpi)
!!print *,'mype RECEIVES itag for levmax=',mype,itag,levmax_recv
do level=levmin_recv,levmax_recv
if (.not.writelevel(level)) cycle
do Morton_no=Morton_start(ipe),Morton_stop(ipe)
itag=Morton_no
call mpi_recv_wrapper(igrid_recv,1,MPI_INTEGER, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
itag=igrid_recv
call mpi_recv_wrapper(level_recv,1,MPI_INTEGER, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
if (level_recv/=level) cycle
call mpi_recv_wrapper(cond_grid_recv,1,MPI_LOGICAL, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
if(.not.cond_grid_recv)cycle
itag=Morton_no
siz_ind=4*3
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);
call mpi_recv_wrapper(normconv,nw+nwauxio+1, MPI_DOUBLE_PRECISION,ipe,&
itag,icomm,intstatus(:,1),ierrmpi)
call mpi_recv_wrapper(wC_TMP_recv,1,type_block_wc_io, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
call mpi_recv_wrapper(xC_TMP_recv,1,type_block_xc_io, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
itag=igrid_recv
call mpi_recv_wrapper(wCC_TMP_recv,1,type_block_wcc_io, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
call mpi_recv_wrapper(xCC_TMP_recv,1,type_block_xcc_io, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
call write_vtk(qunit,ixGlo1,ixGlo2,ixGlo3,ixGhi1,ixGhi2,ixGhi3,&
ixrvCmin1,ixrvCmin2,ixrvCmin3,ixrvCmax1,ixrvCmax2,ixrvCmax3,&
ixrvCCmin1,ixrvCCmin2,ixrvCCmin3,ixrvCCmax1,ixrvCCmax2,&
ixrvCCmax3,igrid_recv,nc,np,nx1,nx2,nx3,nxC1,nxC2,nxC3,&
normconv,wnamei,xC_TMP_recv,xCC_TMP_recv,wC_TMP_recv,&
wCC_TMP_recv)
end do ! Morton_no loop
end do ! level loop
end do ! processor loop
end if ! multiple processors
write(qunit,'(a)')'</UnstructuredGrid>'
write(qunit,'(a)')'</VTKFile>'
close(qunit)
end if
if (npe>1) then
call MPI_BARRIER(icomm,ierrmpi)
if(mype==0)deallocate(intstatus)
end if
end subroutine unstructuredvtk_mpi