Composing VTK file from multiple MPI outputs

Question

For a Lattice Boltzmann simulation of a lid-driven cavity (CFD) I'm decomposing my cubic domain into (also cubic) 8 subdomains, which are computed independently by 8 ranks. Each MPI rank is producing a VTK file for each timestep and since I'm using ParaView I want to visualize the whole thing as one cube. To be more specific about what I am trying to achieve:

• I have a cube with length 8 (number of elements for each direction) => 8x8x8 = 512 elements.
• Each dimension is distributed to 2 ranks, i.e. every rank handles 4x4x4 = 64 elements.
• Every rank writes it's result to a file called lbm_out_<rank>.<timestep>.vts in a VTK StructuredGrid format.
• I want to produce a .pvts file that collects the *.vts files and combines the files containing the subdomains to a single file that ParaView can treat as whole domain.

Unfortunately I'm facing many issues with that since I feel ParaView and VTK are extremely poorly documented and the error messages from ParaView are totally useless.

I'm having the following *.pvts file, which includes a ghost layer and:

<?xml version="1.0"?>
<VTKFile type="PStructuredGrid" version="0.1" byte_order="LittleEndian">
    <PStructuredGrid WholeExtent="0 7 0 7 0 7 " GhostLevel="1">
        <PPoints>
            <PDataArray NumberOfComponents="3" type="Float32" />
        </PPoints>
        <Piece Extent="0 4 0 4 0 4" Source="lbm_out_0.0.vts"/>
        <Piece Extent="3 7 0 4 0 4" Source="lbm_out_1.0.vts"/>
        <Piece Extent="0 4 3 7 0 4" Source="lbm_out_2.0.vts"/>
        <Piece Extent="3 7 3 7 0 4" Source="lbm_out_3.0.vts"/>
        <Piece Extent="0 4 0 4 3 7" Source="lbm_out_4.0.vts"/>
        <Piece Extent="3 7 0 4 3 7" Source="lbm_out_5.0.vts"/>
        <Piece Extent="0 4 3 7 3 7" Source="lbm_out_6.0.vts"/>
        <Piece Extent="3 7 3 7 3 7" Source="lbm_out_7.0.vts"/>
    </PStructuredGrid>
</VTKFile>

Having that file, which I feel should work correctly (note that there are not parameters yet, just plain geometry information), my domain ranges are totally messed up, although each *.vts file works fine on its own. I have a screenshot attached to make things more clear:

What may be the problem? Is it possible to use legacy VTK files for this tasks? May I be doing something totally wrong? I really don't know how to accomplish this task and the resources I find in google are very limited. Thank you.

Answer 1

Unfortunately there is no example for vtkXMLPStructuredGridWriter class (VTK Classes not used in the Examples). So I decided to write the simplest code to generate *.vts and .pvts files for a structured grid, very similar to the case you are looking for.

The following code uses MPI and VTK to write parallel structured grid files. In this example, we have two processes which create their own .vts files and the vtkXMLPStructuredGridWriter class writes the .pvts file:

// MPI Library
#include <mpi.h>

//VTK Library
#include <vtkXMLPStructuredGridWriter.h>
#include <vtkStructuredGrid.h>
#include <vtkSmartPointer.h>
#include <vtkFloatArray.h>
#include <vtkCellData.h>
#include <vtkMPIController.h>
#include <vtkProgrammableFilter.h>
#include <vtkInformation.h>

struct Args {
  vtkProgrammableFilter* pf;
  int local_extent[6];
};

// function to operate on the point attribute data
void execute (void* arg) {
  Args* args = reinterpret_cast<Args*>(arg);
  auto info = args->pf->GetOutputInformation(0);
  auto output_tmp = args->pf->GetOutput();
  auto input_tmp  = args->pf->GetInput();
  vtkStructuredGrid* output = dynamic_cast<vtkStructuredGrid*>(output_tmp);
  vtkStructuredGrid* input  = dynamic_cast<vtkStructuredGrid*>(input_tmp);
  output->ShallowCopy(input);
  output->SetExtent(args->local_extent);
}

int main (int argc, char *argv[]) {
  MPI_Init (&argc, &argv);
  int myrank;
  MPI_Comm_rank (MPI_COMM_WORLD, &myrank);

  int lx {10}, ly{10}, lz{10};        //local dimensions of the process's grid
  int dims[3] = {lx+1, ly+1, lz+1};
  int global_extent[6] = {0, 2*lx, 0, ly, 0, lz};
  int local_extent[6] = {myrank*lx, (myrank+1)*lx, 0, ly, 0, lz};

  // Create and Initialize vtkMPIController
  auto contr = vtkSmartPointer<vtkMPIController>::New();
  contr->Initialize(&argc, &argv, 1);
  int nranks = contr->GetNumberOfProcesses();
  int rank   = contr->GetLocalProcessId();

  // Create grid points, allocate memory and Insert them
  auto points = vtkSmartPointer<vtkPoints>::New();
  points->Allocate(dims[0]*dims[1]*dims[2]);
  for (int k=0; k<dims[2]; ++k)
    for (int j=0; j<dims[1]; ++j)
      for (int i=0; i<dims[0]; ++i)
        points->InsertPoint(i + j*dims[0] + k*dims[0]*dims[1],
                            i+rank*(dims[0]-1), j, k);

  // Create a density field. Note that the number of cells is always less than
  // number of grid points by an amount of one so we use dims[i]-1
  auto density = vtkSmartPointer<vtkFloatArray>::New();
  density->SetNumberOfComponents(1);
  density->SetNumberOfTuples((dims[0]-1)*(dims[1]-1)*(dims[2]-1));
  density->SetName ("density");
  int index;
  for (int k=0; k<lz; ++k)
    for (int j=0; j<ly; ++j)
      for (int i=0; i<lx; ++i) {
        index = i + j*(dims[0]-1) + k*(dims[0]-1)*(dims[1]-1);
        density->SetValue(index, i+j+k);
      }

  // Create a vtkProgrammableFilter
  auto pf = vtkSmartPointer<vtkProgrammableFilter>::New();

  // Initialize an instance of Args
  Args args;
  args.pf = pf;
  for(int i=0; i<6; ++i) args.local_extent[i] = local_extent[i];

  pf->SetExecuteMethod(execute, &args);

  // Create a structured grid and assign point data and cell data to it
  auto structuredGrid = vtkSmartPointer<vtkStructuredGrid>::New();
  structuredGrid->SetExtent(global_extent);
  pf->SetInputData(structuredGrid);
  structuredGrid->SetPoints(points);
  structuredGrid->GetCellData()->AddArray(density);

  // Create the parallel writer and call some functions
  auto parallel_writer = vtkSmartPointer<vtkXMLPStructuredGridWriter>::New();
  parallel_writer->SetInputConnection(pf->GetOutputPort());
  parallel_writer->SetController(contr);
  parallel_writer->SetFileName("test.pvts");
  parallel_writer->SetNumberOfPieces(nranks);
  parallel_writer->SetStartPiece(rank);
  parallel_writer->SetEndPiece(rank);
  parallel_writer->SetDataModeToBinary();
  parallel_writer->Update();
  parallel_writer->Write();

  contr->Finalize();

  // WARNING: it seems that MPI_Finalize is not necessary since we are using
  // Finalize() function from vtkMPIController class. Uncomment the following
  // line and see what happens.
//   MPI_Finalize ();
  return 0;
}

This code only writes some data (in this case density which is a scalar) into a file and do not render it. For visualization you will need a software like Paraview. To run this code, you can use this CMake file:

cmake_minimum_required(VERSION 2.8)

PROJECT(PXMLStructuredGridWriter)
add_executable(PXMLStructuredGridWriter parallel_vtk.cpp)

find_package(VTK REQUIRED)
include(${VTK_USE_FILE})
target_link_libraries(PXMLStructuredGridWriter ${VTK_LIBRARIES})

find_package(MPI REQUIRED)
include_directories(${MPI_INCLUDE_PATH})
target_link_libraries(PXMLStructuredGridWriter ${MPI_LIBRARIES})

At the end you will find an xml file like this in the same directory that you have the executable:

<?xml version="1.0"?>
<VTKFile type="PStructuredGrid" version="0.1" byte_order="LittleEndian" header_type="UInt32" compressor="vtkZLibDataCompressor">
  <PStructuredGrid WholeExtent="0 20 0 10 0 10" GhostLevel="2">
    <PCellData>
      <PDataArray type="Float32" Name="density"/>
    </PCellData>
    <PPoints>
      <PDataArray type="Float32" Name="Points" NumberOfComponents="3"/>
    </PPoints>
    <Piece Extent="0 10 0 10 0 10" Source="test_0.vts"/>
    <Piece Extent="10 20 0 10 0 10" Source="test_1.vts"/>
  </PStructuredGrid>
</VTKFile>