diff --git a/CMakeLists.txt b/CMakeLists.txt
index cf5b06953bb2fc2e2841c38edd40f2532cbcbe67..91e19f2c418ce6b46e450c54e9cdbefb984e43fe 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -174,6 +174,9 @@ include_directories(src/output)
 include_directories(src/utils)
 include_directories(src/scheme)
 
+# Pastis generated sources
+include_directories(${PASTIS_BINARY_DIR}/src/utils)
+
 # Pastis tests
 
 set(CATCH_MODULE_PATH "${PASTIS_SOURCE_DIR}/packages/Catch2")
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 276bad33d5624f018b9d512743f4f7695a6213e5..d1d6d5cdd044c62fa661769dccd55bbb972ae621 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -9,6 +9,7 @@
 #include <IConnectivity.hpp>
 
 #include <ItemValue.hpp>
+#include <Messenger.hpp>
 
 class VTKWriter
 {
@@ -18,6 +19,24 @@ class VTKWriter
   double m_last_time;
   const double m_time_period;
 
+  std::string _getFilenamePVTU()
+  {
+    std::ostringstream sout;
+    sout << m_base_filename;
+    sout << '.' << std::setfill('0') << std::setw(4) << m_file_number << ".pvtu";
+    return sout.str();
+  }
+  std::string _getFilenameVTU(const int& rank_number) const
+  {
+    std::ostringstream sout;
+    sout << m_base_filename;
+    if (parallel::size() > 1) {
+      sout << '-' << std::setfill('0') << std::setw(4) << rank_number;
+    }
+    sout << '.' << std::setfill('0') << std::setw(4) << m_file_number << ".vtu";
+    return sout.str();
+  }
+
  public:
   template <typename MeshType>
   void write(const MeshType& mesh,
@@ -30,110 +49,135 @@ class VTKWriter
     } else {
       return;
     }
-    std::ostringstream sout;
-    sout << m_base_filename << '.' << std::setfill('0') << std::setw(4) << m_file_number << ".vtu" << std::ends;
-    std::ofstream fout(sout.str());
-    fout << "<?xml version=\"1.0\"?>\n";
-    fout << "<VTKFile type=\"UnstructuredGrid\">\n";
-    fout << "<UnstructuredGrid>\n";
-    fout << "<Piece NumberOfPoints=\""<< mesh.numberOfNodes()
-         << "\" NumberOfCells=\"" << mesh.numberOfCells() << "\">\n";
-
-    fout << "<Points>\n";
-    fout << "<DataArray Name=\"Positions\" NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">\n";
-    using Rd = TinyVector<MeshType::dimension>;
-    const NodeValue<const Rd>& xr = mesh.xr();
-    if constexpr(MeshType::dimension == 1) {
-      for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
-        for (unsigned short i=0; i<1; ++i) {
-          fout << xr[r][i] << ' ';
-        }
-        fout << "0 0 "; // VTK requires 3 components
-      }
-    } else if constexpr (MeshType::dimension == 2) {
-      for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
-        for (unsigned short i=0; i<2; ++i) {
-          fout << xr[r][i] << ' ';
-        }
-        fout << "0 "; // VTK requires 3 components
-      }
-    } else {
-      for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
-        for (unsigned short i=0; i<3; ++i) {
-          fout << xr[r][i] << ' ';
-        }
-      }
-    }
-    fout << '\n';
-    fout << "</DataArray>\n";
-    fout << "</Points>\n";
 
-    fout << "<Cells>\n";
+    if (parallel::rank() == 0)
+    { // write PVTK file
+      std::ofstream fout(_getFilenamePVTU());
+      fout << "<?xml version=\"1.0\"?>\n";
+      fout << "<VTKFile type=\"PUnstructuredGrid\">\n";
+      fout << "<PUnstructuredGrid GhostLevel=\"0\">\n";
 
-    fout << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n";
+      fout << "<PPoints>\n";
+      fout << "<PDataArray Name=\"Positions\" NumberOfComponents=\"3\" type=\"Float64\"/>\n";
+      fout << "</PPoints>\n";
 
-    const auto& cell_to_node_matrix
-        = mesh.connectivity().cellToNodeMatrix();
+      fout << "<PCells>\n";
+      fout << "<PDataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\"/>\n";
+      fout << "<PDataArray type=\"UInt32\" Name=\"offsets\" NumberOfComponents=\"1\"/>\n";
+      fout << "<PDataArray type=\"Int8\" Name=\"types\" NumberOfComponents=\"1\"/>\n";
+      fout << "</PCells>\n";
 
-    for (CellId j=0; j<mesh.numberOfCells(); ++j) {
-      const auto& cell_nodes = cell_to_node_matrix[j];
-      for (unsigned short r=0; r<cell_nodes.size(); ++r) {
-        fout << cell_nodes[r] << ' ';
+      for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
+        fout << "<Piece Source=\""<< _getFilenameVTU(i_rank) << "\"/>\n";
       }
-    }
-    fout << '\n';
-    fout << "</DataArray>\n";
+      fout << "</PUnstructuredGrid>\n";
+      fout << "</VTKFile>\n";
 
-    fout << "<DataArray type=\"UInt32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">\n";
-    {
-      unsigned int offset=0;
-      for (CellId j=0; j<mesh.numberOfCells(); ++j) {
-        const auto& cell_nodes = cell_to_node_matrix[j];
-        offset += cell_nodes.size();
-        fout << offset << ' ';
-      }
     }
-    fout << '\n';
-    fout << "</DataArray>\n";
-
-    fout << "<DataArray type=\"Int8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">\n";
-    for (CellId j=0; j<mesh.numberOfCells(); ++j) {
-      const auto& cell_nodes = cell_to_node_matrix[j];
-      switch (cell_nodes.size()) {
-        case 2: {
-          fout << "3 ";
-          break;
+
+    { // write VTK files
+      std::ofstream fout(_getFilenameVTU(parallel::rank()));
+      fout << "<?xml version=\"1.0\"?>\n";
+      fout << "<VTKFile type=\"UnstructuredGrid\">\n";
+      fout << "<UnstructuredGrid>\n";
+      fout << "<Piece NumberOfPoints=\""<< mesh.numberOfNodes()
+           << "\" NumberOfCells=\"" << mesh.numberOfCells() << "\">\n";
+
+      fout << "<Points>\n";
+      fout << "<DataArray Name=\"Positions\" NumberOfComponents=\"3\" type=\"Float64\" format=\"ascii\">\n";
+      using Rd = TinyVector<MeshType::dimension>;
+      const NodeValue<const Rd>& xr = mesh.xr();
+      if constexpr(MeshType::dimension == 1) {
+        for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
+          for (unsigned short i=0; i<1; ++i) {
+            fout << xr[r][i] << ' ';
+          }
+          fout << "0 0 "; // VTK requires 3 components
         }
-        case 3: {
-          fout << "5 ";
-          break;
+      } else if constexpr (MeshType::dimension == 2) {
+        for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
+          for (unsigned short i=0; i<2; ++i) {
+            fout << xr[r][i] << ' ';
+          }
+          fout << "0 "; // VTK requires 3 components
         }
-        case 4: {
-          if (mesh.meshDimension() == 3) {
-            fout << "10 ";
-          } else {
-            fout << "9 ";
+      } else {
+        for (NodeId r=0; r<mesh.numberOfNodes(); ++r) {
+          for (unsigned short i=0; i<3; ++i) {
+            fout << xr[r][i] << ' ';
           }
-          break;
         }
-        case 8: {
-          fout << "12 ";
-          break;
+      }
+      fout << '\n';
+      fout << "</DataArray>\n";
+      fout << "</Points>\n";
+
+      fout << "<Cells>\n";
+
+      fout << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n";
+
+      const auto& cell_to_node_matrix
+          = mesh.connectivity().cellToNodeMatrix();
+
+      for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+        for (unsigned short r=0; r<cell_nodes.size(); ++r) {
+          fout << cell_nodes[r] << ' ';
         }
-        default: {
-          fout << "7 ";
-          break;
+      }
+      fout << '\n';
+      fout << "</DataArray>\n";
+
+      fout << "<DataArray type=\"UInt32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">\n";
+      {
+        unsigned int offset=0;
+        for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+          const auto& cell_nodes = cell_to_node_matrix[j];
+          offset += cell_nodes.size();
+          fout << offset << ' ';
         }
       }
-    }
-    fout << '\n';
-    fout << "</DataArray>\n";
+      fout << '\n';
+      fout << "</DataArray>\n";
 
-    fout << "</Cells>\n";
-    fout << "</Piece>\n";
-    fout << "</UnstructuredGrid>\n";
-    fout << "</VTKFile>\n";
+      fout << "<DataArray type=\"Int8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">\n";
+      for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+        switch (cell_nodes.size()) {
+          case 2: {
+            fout << "3 ";
+            break;
+          }
+          case 3: {
+            fout << "5 ";
+            break;
+          }
+          case 4: {
+            if (mesh.meshDimension() == 3) {
+              fout << "10 ";
+            } else {
+              fout << "9 ";
+            }
+            break;
+          }
+          case 8: {
+            fout << "12 ";
+            break;
+          }
+          default: {
+            fout << "7 ";
+            break;
+          }
+        }
+      }
+      fout << '\n';
+      fout << "</DataArray>\n";
 
+      fout << "</Cells>\n";
+      fout << "</Piece>\n";
+      fout << "</UnstructuredGrid>\n";
+      fout << "</VTKFile>\n";
+    }
     m_file_number++;
   }
   VTKWriter(const std::string& base_filename,