Merge branch 'mpiio-gf-dev' of github.com:mfem/mfem into mpiio-gf-dev

examples/ex1p.cpp

@@ -70,6 +70,7 @@ int main(int argc, char *argv[])
    bool pa = false;
    const char *device_config = "cpu";
    bool visualization = true;
+   int nfiles = 1;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -86,6 +87,7 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&nfiles, "-nf", "--num-files", "Number of files to write.");
    args.Parse();
    if (!args.Good())
    {
@@ -158,7 +160,7 @@ int main(int argc, char *argv[])
    {
       fec = new H1_FECollection(order = 1, dim);
    }
-   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec, 1, 0);
    HYPRE_Int size = fespace->GlobalTrueVSize();
    if (myid == 0)
    {
@@ -237,20 +239,52 @@ int main(int argc, char *argv[])
    //     local finite element solution on each processor.
    a->RecoverFEMSolution(X, *b, x);
 
+   std::string filename("nranks_");
+   filename += to_string(num_procs);
+   filename += ".gf";
+   {
+      double t1;
+      t1 = MPI_Wtime();
+      x.Save(filename.c_str(), nfiles);
+      double t2 = MPI_Wtime();
+
+      if (myid == 0)
+      {
+         err << "elapsed write time: " << t2 - t1 << endl;
+      }
+   }
+   {
+      double t1;
+      t1 = MPI_Wtime();
+      ParGridFunction new_x(fespace, filename.c_str());
+      double t2 = MPI_Wtime();
+      if (myid == 0)
+      {
+         err << "elapsed read time: " << t2 - t1 << endl;
+      }
+      // new_x -= x;
+      // out << "GF difference: " << new_x.Norml1() << endl;
+   }
+
    // 15. Save the refined mesh and the solution in parallel. This output can
    //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
    {
       ostringstream mesh_name, sol_name;
-      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
-      sol_name << "sol." << setfill('0') << setw(6) << myid;
-
-      ofstream mesh_ofs(mesh_name.str().c_str());
-      mesh_ofs.precision(8);
-      pmesh->Print(mesh_ofs);
+      //mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << num_procs << setfill('0') << setw(6) << myid;
 
+      //ofstream mesh_ofs(mesh_name.str().c_str());
+      //mesh_ofs.precision(8);
+      //pmesh->Print(mesh_ofs);
+      double t1 = MPI_Wtime();
       ofstream sol_ofs(sol_name.str().c_str());
       sol_ofs.precision(8);
       x.Save(sol_ofs);
+      double t2 = MPI_Wtime();
+      if (myid == 0)
+      {
+         err << t2 - t1 << endl;
+      }
    }
 
    // 16. Send the solution by socket to a GLVis server.

examples/ex2p.cpp

@@ -274,6 +274,13 @@ int main(int argc, char *argv[])
       pmesh->SetNodalFESpace(fespace);
    }
 
+   {
+      x.Save("ex2p.gf", 1);
+      ParGridFunction new_x(fespace, "ex2p.gf");
+      new_x -= x;
+      out << "GF difference: " << new_x.Norml1() << endl;
+   }
+
    // 16. Save in parallel the displaced mesh and the inverted solution (which
    //     gives the backward displacements to the original grid). This output
    //     can be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".

fem/pgridfunc.cpp

@@ -16,6 +16,7 @@
 #include "fem.hpp"
 #include <iostream>
 #include <limits>
+#include <string>
 #include "../general/forall.hpp"
 using namespace std;
 
@@ -78,6 +79,229 @@ ParGridFunction::ParGridFunction(ParMesh *pmesh, std::istream &input)
    fes = pfes;
 }
 
+ParGridFunction::ParGridFunction(ParFiniteElementSpace *pf,
+                                 const char *_filename)
+   : GridFunction(pf), pfes(pf)
+{
+   MPI_Comm fes_comm;
+   int fes_rank, n_fes_ranks;
+   fes_comm = pfes->GetComm();
+   MPI_Comm_size(fes_comm, &n_fes_ranks);
+   MPI_Comm_rank(fes_comm, &fes_rank);
+
+   std::string filename(_filename);
+   std::string file_prefix;
+   std::string file_ext;
+   {
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos)
+      {
+         file_prefix = (filename.substr(0, i));
+         file_ext = (filename.substr(i, filename.length() - i));
+      }
+   }
+
+   int nfiles = 1;
+   if (fes_rank == 0)
+   {
+      int n_rfes_ranks;
+      int tmp[2];
+      std::string mpi_filename;
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos)
+      {
+         mpi_filename = file_prefix + to_string(0) + file_ext;
+      }
+      else
+      {
+         mpi_filename = filename + to_string(0);
+      }
+      
+      MPI_File fh;
+      MPI_File_open(MPI_COMM_SELF, mpi_filename.c_str(), MPI_MODE_RDONLY,
+                                                      MPI_INFO_NULL, &fh);
+      MPI_File_read_at(fh, 0, tmp, 2, MPI_INT, MPI_STATUS_IGNORE);
+      MPI_File_close(&fh);
+      
+      n_rfes_ranks = tmp[0];
+      nfiles = tmp[1];
+
+      MFEM_ASSERT(n_fes_ranks == n_rfes_ranks,
+         "ParGridFunction::ParGridFunction(ParFiniteElementSpace *pf,"
+               " const char *_filename):\n"
+               "\tThe number of MPI ranks used to save the GridFunction is\n"
+               "\tnot the same as the number used to load it!");
+   }
+   MPI_Bcast(&nfiles, 1, MPI_INT, 0, fes_comm);
+
+   int color = fes_rank * nfiles / n_fes_ranks;
+   
+   MPI_Comm file_comm;
+   MPI_Comm_split(fes_comm, color, fes_rank, &file_comm);
+
+   int file_rank, n_file_ranks;
+   MPI_Comm_size(file_comm, &n_file_ranks);
+   MPI_Comm_rank(file_comm, &file_rank);
+
+   std::string mpi_filename;
+   {
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos) {
+         mpi_filename = file_prefix + std::to_string(color) + file_ext;
+      }
+      else
+      {
+         mpi_filename = filename + std::to_string(color);
+      }
+   }
+
+   MPI_File fh;
+   MPI_File_open(file_comm, mpi_filename.c_str(), MPI_MODE_RDONLY,
+                                                  MPI_INFO_NULL, &fh);
+
+   int *dof_counts = new int[5*n_file_ranks];
+   int **nv = new int*[n_file_ranks];
+   int **nvdofs = new int*[n_file_ranks];
+   int **nedofs = new int*[n_file_ranks];
+   int **nfdofs = new int*[n_file_ranks];
+   int **nrdofs = new int*[n_file_ranks];
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      nv[i]     = &dof_counts[i*5+0];
+      nvdofs[i] = &dof_counts[i*5+1];
+      nedofs[i] = &dof_counts[i*5+2];
+      nfdofs[i] = &dof_counts[i*5+3];
+      nrdofs[i] = &dof_counts[i*5+4];
+   }
+
+   *nv[file_rank] = pfes->GetVSize();
+   *nvdofs[file_rank] = pfes->GetNVDofs();
+   *nedofs[file_rank] = pfes->GetNEDofs();
+   *nfdofs[file_rank] = pfes->GetNFDofs();
+
+   int vdim = pfes->GetVDim();
+   *nrdofs[file_rank] = *nv[file_rank] / vdim - *nvdofs[file_rank] -
+                                       *nedofs[file_rank] - *nfdofs[file_rank];
+
+   MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &dof_counts[0], 5,
+                 MPI_INT, file_comm);
+
+   double *data_ = HostWrite();
+
+   MPI_Offset header_offset = 0;
+   header_offset += 2 * sizeof(int);
+   MPI_Offset v_offset, e_offset, f_offset, r_offset;
+
+   int total_vdofs = 0, total_edofs = 0, total_fdofs = 0, total_rdofs = 0;
+   int total_scalar_dofs = 0;
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      total_vdofs += *nvdofs[i];
+      total_edofs += *nedofs[i];
+      total_fdofs += *nfdofs[i];
+      total_rdofs += *nrdofs[i];
+      total_scalar_dofs += *nv[i];
+   }
+
+   total_scalar_dofs /= vdim;
+
+   if (pfes->GetOrdering() == Ordering::byNODES)
+   {
+      for (int d = 0; d < vdim; ++d)
+      {
+         int v_data_offset = 0 + *nv[file_rank] * d / vdim ;
+         int e_data_offset = v_data_offset + *nvdofs[file_rank];
+         int f_data_offset = e_data_offset + *nedofs[file_rank];
+         int r_data_offset = f_data_offset + *nfdofs[file_rank];
+
+         v_offset = header_offset;
+         e_offset = header_offset;
+         f_offset = header_offset;
+         r_offset = header_offset;
+         
+         v_offset += total_scalar_dofs * d * sizeof(double);
+         e_offset += (total_vdofs + total_scalar_dofs * d) * sizeof(double);
+         f_offset += (total_vdofs + total_edofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+         r_offset += (total_vdofs + total_edofs + total_fdofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+
+
+         for (int i = 0; i < file_rank; ++i)
+         {
+            v_offset += *nvdofs[i] * sizeof(double);
+            e_offset += *nedofs[i] * sizeof(double);
+            f_offset += *nfdofs[i] * sizeof(double);
+            r_offset += *nrdofs[i] * sizeof(double);
+         }
+
+         MPI_File_read_at_all(fh, v_offset, &data_[v_data_offset],
+                              *nvdofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+         MPI_File_read_at_all(fh, e_offset, &data_[e_data_offset],
+                              *nedofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+         MPI_File_read_at_all(fh, f_offset, &data_[f_data_offset],
+                              *nfdofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+         MPI_File_read_at_all(fh, r_offset, &data_[r_data_offset],
+                              *nrdofs[file_rank], MPI_DOUBLE,
+                              MPI_STATUS_IGNORE);
+      }
+   }
+   else
+   {
+      v_offset = header_offset;
+      e_offset = v_offset + total_vdofs * vdim * sizeof(double);
+      f_offset = e_offset + total_edofs * vdim * sizeof(double);
+      r_offset = f_offset + total_fdofs * vdim * sizeof(double);
+
+      for (int i = 0; i < file_rank; ++i)
+      {
+         v_offset += *nvdofs[i] * sizeof(double) * vdim;
+         e_offset += *nedofs[i] * sizeof(double) * vdim;
+         f_offset += *nfdofs[i] * sizeof(double) * vdim;
+         r_offset += *nrdofs[i] * sizeof(double) * vdim;
+      }
+
+      int v_data_offset = 0;
+      int e_data_offset = v_data_offset + *nvdofs[file_rank] * vdim;
+      int f_data_offset = e_data_offset + *nedofs[file_rank] * vdim;
+      int r_data_offset = f_data_offset + *nfdofs[file_rank] * vdim;
+      
+      MPI_File_read_at_all(fh, v_offset, &data_[v_data_offset],
+                           *nvdofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+      MPI_File_read_at_all(fh, e_offset, &data_[e_data_offset],
+                           *nedofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+      MPI_File_read_at_all(fh, f_offset, &data_[f_data_offset],
+                           *nfdofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+      MPI_File_read_at_all(fh, r_offset, &data_[r_data_offset],
+                           *nrdofs[file_rank] * vdim, MPI_DOUBLE,
+                           MPI_STATUS_IGNORE);
+   }
+
+   MPI_File_close(&fh);
+   MPI_Comm_free(&file_comm);
+
+   for (int i = 0; i < size; i++)
+   {
+      if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
+   }
+
+   delete[] dof_counts;
+   delete[] nv;
+   delete[] nvdofs;
+   delete[] nedofs;
+   delete[] nfdofs;
+   delete[] nrdofs;
+}
+
+
 void ParGridFunction::Update()
 {
    face_nbr_data.Destroy();
@@ -518,6 +742,201 @@ void ParGridFunction::Save(adios2stream &out,
 }
 #endif
 
+void ParGridFunction::Save(const char *_filename, const int nfiles)
+{
+   MPI_Comm fes_comm;
+   int fes_rank, n_fes_ranks;
+   fes_comm = pfes->GetComm();
+
+   MPI_Comm_size(fes_comm, &n_fes_ranks);
+   MPI_Comm_rank(fes_comm, &fes_rank);
+
+   int color = fes_rank * nfiles / n_fes_ranks;
+   
+   MPI_Comm file_comm;
+   MPI_Comm_split(fes_comm, color, fes_rank, &file_comm);
+
+   int file_rank, n_file_ranks;
+   MPI_Comm_size(file_comm, &n_file_ranks);
+   MPI_Comm_rank(file_comm, &file_rank);
+
+   std::string filename(_filename);
+   std::string file_prefix;
+   std::string file_ext;
+   std::string mpi_filename;
+   {
+      size_t i = filename.rfind('.', filename.length());
+      if (i != string::npos)
+      {
+         file_prefix = (filename.substr(0, i));
+         file_ext = (filename.substr(i, filename.length() - i));
+         mpi_filename = file_prefix + std::to_string(color) + file_ext;
+      }
+      else
+      {
+         mpi_filename = filename + std::to_string(color);
+      }
+   }
+
+   MPI_File fh;
+   MPI_File_open(file_comm, mpi_filename.c_str(), MPI_MODE_CREATE |
+                                                  MPI_MODE_WRONLY,
+                                                  MPI_INFO_NULL, &fh);
+
+   int *dof_counts = new int[5*n_file_ranks];
+   int **nv = new int*[n_file_ranks];
+   int **nvdofs = new int*[n_file_ranks];
+   int **nedofs = new int*[n_file_ranks];
+   int **nfdofs = new int*[n_file_ranks];
+   int **nrdofs = new int*[n_file_ranks];
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      nv[i]     = &dof_counts[i*5+0];
+      nvdofs[i] = &dof_counts[i*5+1];
+      nedofs[i] = &dof_counts[i*5+2];
+      nfdofs[i] = &dof_counts[i*5+3];
+      nrdofs[i] = &dof_counts[i*5+4];
+   }
+
+   *nv[file_rank] = pfes->GetVSize();
+   *nvdofs[file_rank] = pfes->GetNVDofs();
+   *nedofs[file_rank] = pfes->GetNEDofs();
+   *nfdofs[file_rank] = pfes->GetNFDofs();
+
+   int vdim = pfes->GetVDim();
+   *nrdofs[file_rank] = *nv[file_rank] / vdim - *nvdofs[file_rank] -
+                                       *nedofs[file_rank] - *nfdofs[file_rank];
+
+   MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &dof_counts[0], 5,
+                 MPI_INT, file_comm);
+
+   double *data_  = const_cast<double*>(HostRead());
+   for (int i = 0; i < size; i++)
+   {
+      if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
+   }
+
+   MPI_Offset header_offset = 0;
+
+   if (file_rank == 0)
+   {
+      int tmp[] = {n_fes_ranks, nfiles};
+      MPI_File_write_at(fh, header_offset, &tmp, 2, MPI_INT,
+                        MPI_STATUS_IGNORE);
+   }
+
+   header_offset += 2 * sizeof(int);
+
+   MPI_Offset v_offset, e_offset, f_offset, r_offset;
+
+   int total_vdofs = 0, total_edofs = 0, total_fdofs = 0, total_rdofs = 0;
+   int total_scalar_dofs = 0;
+
+   for (int i = 0; i < n_file_ranks; ++i)
+   {
+      total_vdofs += *nvdofs[i];
+      total_edofs += *nedofs[i];
+      total_fdofs += *nfdofs[i];
+      total_rdofs += *nrdofs[i];
+      total_scalar_dofs += *nv[i];
+   }
+
+   total_scalar_dofs /= vdim;
+
+   if (pfes->GetOrdering() == Ordering::byNODES)
+   {
+      for (int d = 0; d < vdim; ++d)
+      {
+         int v_data_offset = 0 + *nv[file_rank] * d / vdim ;
+         int e_data_offset = v_data_offset + *nvdofs[file_rank];
+         int f_data_offset = e_data_offset + *nedofs[file_rank];
+         int r_data_offset = f_data_offset + *nfdofs[file_rank];
+
+         v_offset = header_offset;
+         e_offset = header_offset;
+         f_offset = header_offset;
+         r_offset = header_offset;
+         
+         v_offset += total_scalar_dofs * d * sizeof(double);
+         e_offset += (total_vdofs + total_scalar_dofs * d) * sizeof(double);
+         f_offset += (total_vdofs + total_edofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+         r_offset += (total_vdofs + total_edofs + total_fdofs +
+                                    total_scalar_dofs * d) * sizeof(double);
+
+         for (int i = 0; i < file_rank; ++i)
+         {
+            v_offset += *nvdofs[i] * sizeof(double);
+            e_offset += *nedofs[i] * sizeof(double);
+            f_offset += *nfdofs[i] * sizeof(double);
+            r_offset += *nrdofs[i] * sizeof(double);
+         }
+
+         MPI_File_write_at_all(fh, v_offset, &data_[v_data_offset],
+                               *nvdofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+         MPI_File_write_at_all(fh, e_offset, &data_[e_data_offset],
+                               *nedofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+         MPI_File_write_at_all(fh, f_offset, &data_[f_data_offset],
+                               *nfdofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+         MPI_File_write_at_all(fh, r_offset, &data_[r_data_offset],
+                               *nrdofs[file_rank], MPI_DOUBLE,
+                               MPI_STATUS_IGNORE);
+      }
+   }
+   else
+   {
+      v_offset = header_offset;
+      e_offset = v_offset + total_vdofs * vdim * sizeof(double);
+      f_offset = e_offset + total_edofs * vdim * sizeof(double);
+      r_offset = f_offset + total_fdofs * vdim * sizeof(double);
+
+      for (int i = 0; i < file_rank; ++i)
+      {
+         v_offset += *nvdofs[i] * sizeof(double) * vdim;
+         e_offset += *nedofs[i] * sizeof(double) * vdim;
+         f_offset += *nfdofs[i] * sizeof(double) * vdim;
+         r_offset += *nrdofs[i] * sizeof(double) * vdim;
+      }
+
+      int v_data_offset = 0;
+      int e_data_offset = v_data_offset + *nvdofs[file_rank] * vdim;
+      int f_data_offset = e_data_offset + *nedofs[file_rank] * vdim;
+      int r_data_offset = f_data_offset + *nfdofs[file_rank] * vdim;
+      
+      MPI_File_write_at_all(fh, v_offset, &data_[v_data_offset],
+                            *nvdofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+      MPI_File_write_at_all(fh, e_offset, &data_[e_data_offset],
+                            *nedofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+      MPI_File_write_at_all(fh, f_offset, &data_[f_data_offset],
+                            *nfdofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+      MPI_File_write_at_all(fh, r_offset, &data_[r_data_offset],
+                            *nrdofs[file_rank] * vdim, MPI_DOUBLE,
+                            MPI_STATUS_IGNORE);
+   }
+
+   MPI_File_close(&fh);
+   MPI_Comm_free(&file_comm);
+   
+   for (int i = 0; i < size; i++)
+   {
+      if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
+   }
+
+   delete[] dof_counts;
+   delete[] nv;
+   delete[] nvdofs;
+   delete[] nedofs;
+   delete[] nfdofs;
+   delete[] nrdofs;
+}
+
 void ParGridFunction::SaveAsOne(std::ostream &out)
 {
    int i, p;

fem/pgridfunc.hpp

@@ -83,6 +83,13 @@ public:
        constructed. The new ParGridFunction assumes ownership of both. */
    ParGridFunction(ParMesh *pmesh, std::istream &input);
 
+   /// Construct a ParGridFunction by loading a ParGridFunction saved using
+   /// ParGridFunction::Save(char *filename, int nfiles).
+   /** The parallel space @a *pf and the space used by the GridFunction saved
+       in @a *filename should match. The number of ranks used when loading the
+       ParGridFunction must be the same as when it was saved. */
+   ParGridFunction(ParFiniteElementSpace *pf, const char *filename);
+
    /// Copy assignment. Only the data of the base class Vector is copied.
    /** It is assumed that this object and @a rhs use ParFiniteElementSpace%s
        that have the same size.
@@ -324,6 +331,20 @@ public:
       const adios2stream::data_type type = adios2stream::data_type::point_data) const;
 #endif
 
+   /** Save the local grid functions to n number of files, where each file will
+       contain the grid functions from potentially multiple ranks. This is
+       similar to the syncIO approach from "Fu, Jing, et al. 'Scalable parallel
+       I/O alternatives for massively parallel partitioned solver systems.'
+       2010 IEEE International Symposium on Parallel & Distributed Processing,
+       Workshops and Phd Forum (IPDPSW). IEEE, 2010." 
+       @param[in] filename - filename for output files with extension
+       @param[in] nfiles - number of files to write using MPI-IO
+       @note - takes into account the signs of the local dofs.
+       @note - writes a binary file without the FESpace header; the saved file
+       should only be loaded by the accompanying constructor:
+       ParGridFunction(ParFiniteElementSpace *pf, const char *filename) */
+   void Save(const char *filename, const int nfiles = 1);
+
    /// Merge the local grid functions
    void SaveAsOne(std::ostream &out = mfem::out);