MPIIO .pvtr format file

ParaView Support
1

Hi,

I have a problem with my .pvtr output. I do not know the structure following the header. It appears that the header part is correct, but the data cannot be read correctly. The error from Paraview: No available piece provides data for the following extents. A simple version of the code is:

program write_vtr_and_pvtr
  use mpi
  implicit none

  integer :: ierr, rank, size
  integer, parameter :: nx = 100, ny = 100, nz = 1
  integer :: local_nx, local_ny, local_nz
  real(8), allocatable :: rhox(:,:,:)
  character(len=256) :: vtr_filename, pvtr_filename
  integer :: i, j, k, start_x, end_x, offset
  integer(kind=MPI_OFFSET_KIND) :: disp

  call MPI_INIT(ierr)
  call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
  call MPI_COMM_SIZE(MPI_COMM_WORLD, size, ierr)

  ! Determine local grid size
  local_nx = nx / size
  if (rank < mod(nx, size)) local_nx = local_nx + 1
  local_ny = ny
  local_nz = nz

  ! Determine the start and end x-coordinates for each rank
  start_x = rank * (nx / size)
  if (rank < mod(nx, size)) then
    start_x = start_x + rank
  else
    start_x = start_x + mod(nx, size)
  end if
  end_x = start_x + local_nx - 1

  ! Allocate the data array
  allocate(rhox(local_nx, local_ny, local_nz))

  ! Initialize the data array
  rhox = 1.0 * rank

  ! Create a filename for each rank
  write(vtr_filename, '(A, I0, A)') 'output_rank', rank, '.vtr'

  ! Open the file for writing
  open(unit=10, file=vtr_filename, status='unknown')

  ! Write the VTK XML header
  write(10, '(A)') '<?xml version="1.0"?>'
  write(10, '(A)') '<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian">'
  write(10, '(A,I0,A,I0,A,I0,A)') '  <RectilinearGrid WholeExtent="0 ', nx-1, ' 0 ', ny-1, ' 0 ', nz-1, '">'
  write(10, '(A,I0,A,I0,A,I0,A,I0,A)') '    <Piece Extent="', rank*local_nx, ' ', (rank+1)*local_nx-1, ' 0 ', local_ny-1, ' 0 ', local_nz-1, '">'
  write(10, '(A)') '      <PointData Scalars="scalars">'
  write(10, '(A)') '        <DataArray type="Float64" Name="rho" format="ascii">'

  ! Write the data array in ASCII format
  do k = 1, local_nz
    do j = 1, local_ny
      do i = 1, local_nx
        write(10, '(F12.8)', advance='no') rhox(i, j, k)
        if (i < local_nx) write(10, '(A)', advance='no') ' '
      end do
      write(10, '(A)')
    end do
  end do

  ! Close the DataArray and PointData sections
  write(10, '(A)') '        </DataArray>'
  write(10, '(A)') '      </PointData>'

  ! Write the coordinates in ASCII format
  write(10, '(A)') '      <Coordinates>'
  write(10, '(A)') '        <DataArray type="Float64" Name="X_COORDINATES" NumberOfComponents="1" format="ascii">'
  do i = start_x, end_x
    write(10, '(F12.8)', advance='no') real(i, 8)
    if (i < end_x) write(10, '(A)', advance='no') ' '
  end do
  write(10, '(A)') '</DataArray>'
  write(10, '(A)') '        <DataArray type="Float64" Name="Y_COORDINATES" NumberOfComponents="1" format="ascii">'
  do j = 0, local_ny-1
    write(10, '(F12.8)', advance='no') real(j, 8)
    if (j < local_ny-1) write(10, '(A)', advance='no') ' '
  end do
  write(10, '(A)') '</DataArray>'
  write(10, '(A)') '        <DataArray type="Float64" Name="Z_COORDINATES" NumberOfComponents="1" format="ascii">'
  write(10, '(A)') '        0'
  write(10, '(A)') '</DataArray>'
  write(10, '(A)') '      </Coordinates>'
  write(10, '(A)') '    </Piece>'
  write(10, '(A)') '  </RectilinearGrid>'
  write(10, '(A)') '</VTKFile>'

  ! Close the file
  close(10)

  ! Wait for all ranks to finish writing their .vtr files
  call MPI_BARRIER(MPI_COMM_WORLD, ierr)

  ! Create the .pvtr file on rank 0
  if (rank == 0) then
    pvtr_filename = 'output.pvtr'
    open(unit=20, file=pvtr_filename, status='unknown')

    write(20, '(A)') '<?xml version="1.0"?>'
    write(20, '(A)') '<VTKFile type="PRectilinearGrid" version="0.1" byte_order="LittleEndian">'
    write(20, '(A)') '  <PRectilinearGrid WholeExtent="0 99 0 99 0 0" GhostLevel="0">'
    write(20, '(A)') '    <PPointData Scalars="scalars">'
    write(20, '(A)') '      <PDataArray type="Float64" Name="rho" format="ascii"/>'
    write(20, '(A)') '    </PPointData>'
    write(20, '(A)') '    <PCellData/>'
    write(20, '(A)') '    <PCoordinates>'
    write(20, '(A)') '      <PDataArray type="Float64" Name="X_COORDINATES" NumberOfComponents="1"/>'
    write(20, '(A)') '      <PDataArray type="Float64" Name="Y_COORDINATES" NumberOfComponents="1"/>'
    write(20, '(A)') '      <PDataArray type="Float64" Name="Z_COORDINATES" NumberOfComponents="1"/>'
    write(20, '(A)') '    </PCoordinates>'

    do i = 0, size-1
      write(vtr_filename, '(A, I0, A)') 'output_rank', i, '.vtr'
      if (i < size - 1) then
        
       write(20, '(A,I0,A,I0,A,I0,A,A,A)') '    <Piece Extent="', i * local_nx, ' ', (i + 1) * local_nx - 1, ' 0 ', ny - 1, ' 0 0" Source="', trim(vtr_filename), '"/>'
      else
       write(20, '(A,I0,A,I0,A,I0,A,A,A)') '    <Piece Extent="', i * local_nx, ' ', nx - 1, ' 0 ', ny - 1, ' 0 0" Source="', trim(vtr_filename), '"/>'
      end if
    end do

    write(20, '(A)') '  </PRectilinearGrid>'
    write(20, '(A)') '</VTKFile>'

    close(20)
  end if

  call MPI_FINALIZE(ierr)
end program write_vtr_and_pvtr
2

I’d suggest to not write XML manually but rely on VTK instead to write it.

Alternatively, you may be interested by the VTKHDF format (which does not support rectilinear grid yet).