Added parallel timeseries writer

dune/vtk/utility/uid.hh

+#pragma once
+
+#include <cstdlib>
+#include <cstring>
+#include <ctime>
+#include <string>
+
+namespace Dune
+{
+  inline std::string uid (std::size_t len = 8)
+  {
+    static const auto digits = "0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ";
+    static const int N = std::strlen(digits);
+
+    std::string id(len,' ');
+    for (std::size_t i = 0; i < len; ++i)
+      id[i] = digits[std::rand()%N];
+
+    return id;
+  }
+
+} // end namespace Dune

dune/vtk/vtktimeserieswriter.hh

@@ -6,10 +6,17 @@
 
 #include <dune/vtk/filewriter.hh>
 #include <dune/vtk/vtktypes.hh>
+#include <dune/vtk/utility/filesystem.hh>
+#include <dune/vtk/utility/uid.hh>
 
 namespace Dune
 {
-  /// File-Writer for Vtk .vtu files
+  /// File-Writer for Vtk timeseries .vtu files
+  /**
+   * \tparam VtkWriter  Type of a FileWriter derived from \ref VtkWriterInterface that
+   *                    additionally supports writeTimeseriesSerialFile() and writeTimeseriesParallelFile(),
+   *                    e.g. \ref VtkUnstructuredGridWriter.
+   **/
   template <class VtkWriter>
   class VtkTimeseriesWriter
       : public FileWriter
@@ -18,6 +25,14 @@ namespace Dune
     using Self = VtkTimeseriesWriter;
     using pos_type = typename std::ostream::pos_type;
 
+    template <class W>
+    using HasWriteTimeseriesSerialFile = decltype(&W::writeTimeseriesSerialFile);
+    static_assert(Std::is_detected<HasWriteTimeseriesSerialFile, VtkWriter>::value, "");
+
+    template <class W>
+    using HasWriteTimeseriesParallelFile = decltype(&W::writeTimeseriesParallelFile);
+    static_assert(Std::is_detected<HasWriteTimeseriesParallelFile, VtkWriter>::value, "");
+
   public:
     /// Constructor, stores the gridView
     template <class... Args, disableCopyMove<Self, Args...> = 0>
@@ -25,6 +40,16 @@ namespace Dune
       : vtkWriter_{std::forward<Args>(args)...}
     {
       assert(vtkWriter_.format_ != Vtk::ASCII && "Timeseries writer requires APPENDED mode");
+      std::srand(std::time(nullptr));
+      // put temporary file to /tmp directory
+      tmpDir_ = filesystem::path("/tmp/vtktimeserieswriter_" + uid(10) + "/");
+      assert( filesystem::exists("/tmp") );
+      filesystem::create_directories(tmpDir_);
+    }
+
+    ~VtkTimeseriesWriter ()
+    {
+      std::remove(tmpDir_.string().c_str());
     }
 
     /// Write the attached data to the file with \ref Vtk::FormatTypes and \ref Vtk::DataTypes
@@ -52,6 +77,7 @@ namespace Dune
 
   protected:
     VtkWriter vtkWriter_;
+    filesystem::path tmpDir_;
 
     bool initialized_ = false;
 

dune/vtk/vtktimeserieswriter.impl.hh

@@ -26,22 +26,32 @@ template <class W>
 void VtkTimeseriesWriter<W>
   ::writeTimestep (double time, std::string const& fn)
 {
-  auto p = filesystem::path(fn);
-  auto name = p.stem();
-  p.remove_filename();
-  p /= name.string();
+  auto name = filesystem::path(fn).stem();
+  auto tmp = tmpDir_;
+  tmp /= name.string();
 
   vtkWriter_.dataCollector_.update();
 
+  std::string filenameBase = tmp.string();
+
+  int rank = 0;
+  int num_ranks = 1;
+  #ifdef HAVE_MPI
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    MPI_Comm_size(MPI_COMM_WORLD, &num_ranks);
+    if (num_ranks > 1)
+      filenameBase = tmp.string() + "_p" + std::to_string(rank);
+  #endif
+
   if (!initialized_) {
     // write points and cells only once
-    filenameMesh_ = p.string() + ".mesh.vtkdata";
+    filenameMesh_ = filenameBase + ".mesh.vtkdata";
     std::ofstream out(filenameMesh_, std::ios_base::ate | std::ios::binary);
     vtkWriter_.writeGridAppended(out, blocks_);
     initialized_ = true;
   }
 
-  std::string filenameData = p.string() + "_t" + std::to_string(timesteps_.size()) + ".vtkdata";
+  std::string filenameData = filenameBase + "_t" + std::to_string(timesteps_.size()) + ".vtkdata";
   std::ofstream out(filenameData, std::ios_base::ate | std::ios::binary);
   vtkWriter_.writeDataAppended(out, blocks_);
   timesteps_.emplace_back(time, filenameData);
@@ -59,8 +69,25 @@ void VtkTimeseriesWriter<W>
   auto name = p.stem();
   p.remove_filename();
   p /= name.string();
-  std::string filename = p.string() + ".ts." + vtkWriter_.getFileExtension();
-  vtkWriter_.writeTimeseriesFile(filename, filenameMesh_, timesteps_, blocks_);
+
+  std::string filenameBase = p.string() + "_ts";
+
+  int rank = 0;
+  int num_ranks = 1;
+#ifdef HAVE_MPI
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    MPI_Comm_size(MPI_COMM_WORLD, &num_ranks);
+    if (num_ranks > 1)
+      filenameBase = p.string() + "_ts_p" + std::to_string(rank);
+#endif
+
+  std::string filename = filenameBase + "." + vtkWriter_.getFileExtension();
+  vtkWriter_.writeTimeseriesSerialFile(filename, filenameMesh_, timesteps_, blocks_);
+
+#ifdef HAVE_MPI
+  if (num_ranks > 1 && rank == 0)
+    vtkWriter_.writeTimeseriesParallelFile(p.string() + "_ts", num_ranks, timesteps_);
+#endif
 
   // remove all temporary data files
   int ec = std::remove(filenameMesh_.c_str());

dune/vtk/vtkwriterinterface.impl.hh

@@ -30,21 +30,21 @@ void VtkWriterInterface<GV,DC>
 
   std::string filename = p.string() + "." + fileExtension();
 
+  int rank = 0;
+  int num_ranks = 1;
+
 #ifdef HAVE_MPI
-  int rank = -1;
-  int num_ranks = -1;
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
   MPI_Comm_size(MPI_COMM_WORLD, &num_ranks);
-  if (num_ranks > 1) {
+  if (num_ranks > 1)
     filename = p.string() + "_p" + std::to_string(rank) + "." + fileExtension();
+#endif
 
-    writeSerialFile(filename);
-    if (rank == 0) {
-      writeParallelFile(p.string(), num_ranks);
-    }
-  } else {
-    writeSerialFile(filename);
-  }
+  writeSerialFile(filename);
+
+#ifdef HAVE_MPI
+  if (num_ranks > 1 && rank == 0)
+    writeParallelFile(p.string(), num_ranks);
 #endif
 }
 

dune/vtk/writers/vtkunstructuredgridwriter.hh

@@ -47,10 +47,22 @@ namespace Dune
     virtual void writeParallelFile (std::string const& pfilename, int size) const override;
 
     /// Write a series of timesteps in one file
-    void writeTimeseriesFile (std::string const& filename,
-                              std::string const& filenameMesh,
-                              std::vector<std::pair<double, std::string>> const& timesteps,
-                              std::vector<std::uint64_t> const& blocksize) const;
+    /**
+     * \param filename      The name of the output file
+     * \param filenameMesh  The name of a file where the mesh is stored. Must exist.
+     * \param timesteps     A vector of pairs (timestep, filename) where the filename indicates
+     *                      a file where the data of the timestep is stored.
+     * \param blocks        A list of block sizes of the binary data stored in the files.
+     *                      Order: (points, cells, pointdata[0], celldata[0], pointdata[1], celldata[1],...)
+     **/
+    void writeTimeseriesSerialFile (std::string const& filename,
+                                    std::string const& filenameMesh,
+                                    std::vector<std::pair<double, std::string>> const& timesteps,
+                                    std::vector<std::uint64_t> const& blocks) const;
+
+    /// Write parallel VTK file for series of timesteps
+    void writeTimeseriesParallelFile (std::string const& pfilename, int size,
+                                      std::vector<std::pair<double, std::string>> const& timesteps) const;
 
     virtual std::string fileExtension () const override
     {

dune/vtk/writers/vtkunstructuredgridwriter.impl.hh

@@ -136,10 +136,10 @@ void VtkUnstructuredGridWriter<GV,DC>
 
 template <class GV, class DC>
 void VtkUnstructuredGridWriter<GV,DC>
-  ::writeTimeseriesFile (std::string const& filename,
-                         std::string const& filenameMesh,
-                         std::vector<std::pair<double, std::string>> const& timesteps,
-                         std::vector<std::uint64_t> const& blocks) const
+  ::writeTimeseriesSerialFile (std::string const& filename,
+                               std::string const& filenameMesh,
+                               std::vector<std::pair<double, std::string>> const& timesteps,
+                               std::vector<std::uint64_t> const& blocks) const
 {
   std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
   assert(out.is_open());
@@ -206,15 +206,16 @@ void VtkUnstructuredGridWriter<GV,DC>
   out << "</Piece>\n";
   out << "</UnstructuredGrid>\n";
 
-  pos_type appended_pos = 0;
   out << "<AppendedData encoding=\"raw\">\n_";
-  appended_pos = out.tellp();
+  pos_type appended_pos = out.tellp();
 
-  std::ifstream file_mesh(filenameMesh, std::ios_base::in | std::ios_base::binary);
-  out << file_mesh.rdbuf();
-  file_mesh.close();
-  assert( std::uint64_t(out.tellp()) == std::accumulate(blocks.begin(), std::next(blocks.begin(),shift), std::uint64_t(appended_pos)) );
+  { // write grid (points, cells)
+    std::ifstream file_mesh(filenameMesh, std::ios_base::in | std::ios_base::binary);
+    out << file_mesh.rdbuf();
+    assert( std::uint64_t(out.tellp()) == std::accumulate(blocks.begin(), std::next(blocks.begin(),shift), std::uint64_t(appended_pos)) );
+  }
 
+  // write point-data and cell-data
   for (auto const& timestep : timesteps) {
     std::ifstream file(timestep.second, std::ios_base::in | std::ios_base::binary);
     out << file.rdbuf();
@@ -250,6 +251,78 @@ void VtkUnstructuredGridWriter<GV,DC>
 }
 
 
+template <class GV, class DC>
+void VtkUnstructuredGridWriter<GV,DC>
+  ::writeTimeseriesParallelFile (std::string const& pfilename, int size, std::vector<std::pair<double, std::string>> const& timesteps) const
+{
+  std::string filename = pfilename + ".pvtu";
+  std::ofstream out(filename, std::ios_base::ate | std::ios::binary);
+  assert(out.is_open());
+
+  out << "<VTKFile"
+      << " type=\"PUnstructuredGrid\""
+      << " version=\"1.0\""
+      << " byte_order=\"" << this->getEndian() << "\""
+      << " header_type=\"UInt64\""
+      << (format_ == Vtk::COMPRESSED ? " compressor=\"vtkZLibDataCompressor\"" : "")
+      << ">\n";
+
+  out << "<PUnstructuredGrid GhostLevel=\"0\""
+      << " TimeValues=\"";
+  {
+    std::size_t i = 0;
+    for (auto const& timestep : timesteps)
+      out << timestep.first << (++i % 6 != 0 ? ' ' : '\n');
+  }
+  out << "\">\n";
+
+  // Write points
+  out << "<PPoints>\n";
+  out << "<PDataArray"
+      << " type=\"" << to_string(datatype_) << "\""
+      << " NumberOfComponents=\"3\""
+      << " />\n";
+  out << "</PPoints>\n";
+
+  // Write data associated with grid points
+  out << "<PPointData" << this->getNames(pointData_) << ">\n";
+  for (std::size_t i = 0; i < timesteps.size(); ++i) {
+    for (auto const& v : pointData_) {
+      out << "<PDataArray"
+          << " Name=\"" << v.name() << "\""
+          << " type=\"" << to_string(v.type()) << "\""
+          << " NumberOfComponents=\"" << v.ncomps() << "\""
+          << " TimeStep=\"" << i << "\""
+          << " />\n";
+    }
+  }
+  out << "</PPointData>\n";
+
+  // Write data associated with grid cells
+  out << "<PCellData" << this->getNames(cellData_) << ">\n";
+  for (std::size_t i = 0; i < timesteps.size(); ++i) {
+    for (auto const& v : cellData_) {
+      out << "<PDataArray"
+          << " Name=\"" << v.name() << "\""
+          << " type=\"" << to_string(v.type()) << "\""
+          << " NumberOfComponents=\"" << v.ncomps() << "\""
+          << " TimeStep=\"" << i << "\""
+          << " />\n";
+    }
+  }
+  out << "</PCellData>\n";
+
+  // Write piece file references
+  for (int p = 0; p < size; ++p) {
+    std::string piece_source = pfilename + "_p" + std::to_string(p) + "." + this->fileExtension();
+    out << "<Piece Source=\"" << piece_source << "\" />\n";
+  }
+
+  out << "</PUnstructuredGrid>\n";
+  out << "</VTKFile>";
+}
+
+
 template <class GV, class DC>
 void VtkUnstructuredGridWriter<GV,DC>
   ::writeCells (std::ofstream& out, std::vector<pos_type>& offsets,

src/timeserieswriter.cc

@@ -49,7 +49,7 @@ int main (int argc, char** argv)
 
   using GridType = YaspGrid<3>;
   FieldVector<double,3> upperRight; upperRight = 1.0;
-  auto numElements = filledArray<3,int>(4);
+  auto numElements = filledArray<3,int>(8);
   GridType grid(upperRight, numElements);
   write("yasp", grid.leafGridView());
 }
\ No newline at end of file