Using VTK's parallel XML writer for parallel output

Question

I am currently using VTK to write simulation data on the subdomains of a structured grid to .vts files with the vtkXMLStructuredGridWriter in C++. So far, this works great for me. Now I also want to write the parallel file, which contains the information on all the pieces and this part is painfully undocumented.

My target is to collect all consecutive files of a single process in a dedicated folder, named after the rank of the process, e.g '0/mydata_<timestep>.vts'. The parallel file then should be located in the parent folder and be named something like 'mydata_<timestep>.pvts'. In principle, I could write a parallel file with any xml utility to do that job. But I would prefer to use the vtk library for consistency. I found this and this example to get some idea how that is supposed to work, with not a lot of success. The first example is also mentioned on Stackoverflow.

It does write a parallel file with the information which contains:

• number of pieces
• data sets in the pieces
• the location of the pieces in the folders

But it also does:

• create additional '*.vts' files in the same folder as the 'test_<timestep>.pvts' file
• sets the extent of the whole domain and of each piece to the extent of the root subdomain (process with rank 0)

I am also confused by how this parallel writer is supposed to work. If it is only called on the root process as in the examples, how would it know about the extent of the total subdomain? Or if I should call it on each process in parallel, how would it avoid concurrent access to the file?

So finally: Does anyone know, how this parallel XML writer is supposed to be used? Or maybe knows a good example to learn from?

I use following code to write the files:

#include <vtkNew.h>
#include <vtkStructuredGrid.h>
#include <vtkXMLStructuredGridWriter.h>
#include <vtkXMLPStructuredGridWriter.h>
class MyPwriter : public vtkXMLPStructuredGridWriter {
public:
    static MyPwriter* New();

    void WritePPieceAttributes(int index) {
        std::string filename = std::to_string(index);
        filename += "/test.vts";
        this->WriteStringAttribute("Source", filename.c_str());
    }
};
vtkStandardNewMacro(MyPwriter);

int main(int argc, char* argv[]) {
MPI_Init(&argc, &argv);
int myrank, nproc;
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
MPI_Comm_size(MPI_COMM_WORLD, &nproc);
vtkNew<vtkStructuredGrid> grid;
...  //populating the grid with values
vtkNew<vtkXMLStructuredGridWriter> writer;
writer->SetInputData(grid);
std::string fname = std::to_string(myrank);
fname += "/test_1." + writer->GetDefaultFileExtension();   
writer->SetFileName(fname.c_str());
writer->Write();
// Until here, everything works as expected
// Now I want to write the parallel file
if(myrank == 0) {
    vtkNew<MyPwriter> pwriter;
    std::string pfname = "test_1.";
    pfname += pwriter->GetDefaultFileExtension();
    pwriter->SetFileName(pfname.c_str());
    pwriter->SetNumberOfPieces(nproc);
    pwriter->SetStartPiece(0);
    pwriter->SetEndPiece(nproc-1);
    pwriter->SetInputData(grid);
    pwriter->Update();
}
MPI_Finalize();
return 0;
}

Answer 1

I found following post Composing VTK file from multiple MPI outputs concerning using VTK's parallel xml writer for structured grids. However, this post does not deal yet with halo and/or ghost cells.

To make sure that those are excluded properly when loading the parallel file, I took this approach to override the attributes of the displayed pieces. This has the advantage, that the ghost/halo cells of each piece are still there, if required.

An example with two parallel process then looks like this: Example

Note: I also included the option to save the separate pieces in a subfolder

The code:

#include <fstream>
#include <iostream>
#include <string>
#include <mpi.h>


// VTK Library
#include <vtkCellData.h>
#include <vtkFloatArray.h>
#include <vtkMPIController.h>
#include <vtkProgrammableFilter.h>
#include <vtkNew.h>
#include <vtkStructuredGrid.h>
#include <vtkXMLPStructuredGridWriter.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 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);
}

class CustomvtkXMLPStructuredGridWriter : public vtkXMLPStructuredGridWriter {
public:
    static CustomvtkXMLPStructuredGridWriter* New();

    void SetPPieceExtent(const int extent_l[6]) {
        int nranks{0};
        MPI_Comm_size(MPI_COMM_WORLD, &nranks);
        extent_array.resize(nranks * 6);
        MPI_Allgather(extent_l, 6, MPI_INT, extent_array.data(), 6, MPI_INT, MPI_COMM_WORLD);
    }

    void WritePPieceAttributes(int index) {
        int offset = index * 6;
        std::ostringstream os;
        os << std::to_string(extent_array[offset]);
        for (int n{1}; n < 6; n++)
            os << " " << std::to_string(extent_array[offset + n]);

        std::string piecename(this->CreatePieceFileName(index));
        this->WriteStringAttribute("Source", piecename.c_str());
        this->WriteStringAttribute("Extent", os.str().c_str());
    }

private:
    std::vector<int> extent_array;  // list of the extents of each piece
};

vtkStandardNewMacro(CustomvtkXMLPStructuredGridWriter);

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

{
    // Create and Initialize vtkMPIController
    vtkNew<vtkMPIController> contr;
    contr->Initialize(&argc, &argv, 1);
    int nranks = contr->GetNumberOfProcesses();
    int rank = contr->GetLocalProcessId();

    int ghost_levels = 2;
    // local dimensions of the process's grid
    int lx{10}, ly{10}, lz{10};
    int global_extent[6] = {0, nranks * lx + 2* ghost_levels, 0, ly + 2* ghost_levels, 0, lz + 2* ghost_levels};
    int local_extent[6] = { myrank * lx, (myrank + 1) * lx + 2 * ghost_levels,
                            0, ly + 2 * ghost_levels,
                            0, lz + 2 * ghost_levels };
    double offset[3] = {static_cast<double>(rank * lx + ghost_levels), 0., 0.};
    lx += 2 * ghost_levels;
    ly += 2 * ghost_levels;
    lz += 2 * ghost_levels;
    int dims[3] = {lx + 1, ly + 1, lz + 1};

    // Create grid points, allocate memory and Insert them
    vtkNew<vtkPoints>
        points;
    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) {
                int index = i + j * dims[0] + k * dims[0] * dims[1];
                double x = offset[0] + i;
                double y = offset[1] + j;
                double z = offset[2] + k;
                points->InsertPoint(index, x, y, z);
            }

    // 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
    vtkNew<vtkFloatArray> density;
    density->SetNumberOfComponents(1);
    density->SetNumberOfTuples((dims[0] - 1) * (dims[1] - 1) * (dims[2] - 1));
    density->SetName("density");
    int index;
    for (int k = 0; k < lx; ++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, rank);
            }

    // Create a vtkProgrammableFilter
    vtkNew<vtkProgrammableFilter> pf;

    // 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
    vtkNew<vtkStructuredGrid> structuredGrid;
    structuredGrid->SetExtent(global_extent);
    pf->SetInputData(structuredGrid);
    structuredGrid->SetPoints(points);
    structuredGrid->GetCellData()->AddArray(density);

    // remove the halo cell layers from the local extents, so they are not displayed in the root file
    int ppiece_extent[6];
    ppiece_extent[0] = rank == 0 ? 0 : local_extent[0] + ghost_levels;
    ppiece_extent[1] = rank == (nranks - 1) ? local_extent[1] : local_extent[1] - ghost_levels;
    for (int n = 2; n < 6; n++) ppiece_extent[n] = local_extent[n];

    // Create the parallel writer and call some functions
    vtkNew<CustomvtkXMLPStructuredGridWriter> parallel_writer;
    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->SetUseSubdirectory(true);
    parallel_writer->SetPPieceExtent(ppiece_extent);
    parallel_writer->SetGhostLevel(ghost_levels);
    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;
}