initial commit with parRSB working

config/config.hpp.in

@@ -144,6 +144,9 @@
 // Enable MFEM functionality based on the GSLIB library
 // #define MFEM_USE_GSLIB
 
+// Enable MFEM functionality based on the PARRSB library
+// #define MFEM_USE_PARRSB
+
 // Build the NVIDIA GPU/CUDA-enabled version of the MFEM library.
 // Requires a CUDA compiler (nvcc).
 // #define MFEM_USE_CUDA

config/config.mk.in

@@ -48,6 +48,7 @@ MFEM_USE_CONDUIT       = @MFEM_USE_CONDUIT@
 MFEM_USE_PUMI          = @MFEM_USE_PUMI@
 MFEM_USE_HIOP          = @MFEM_USE_HIOP@
 MFEM_USE_GSLIB         = @MFEM_USE_GSLIB@
+MFEM_USE_PARRSB        = @MFEM_USE_PARRSB@
 MFEM_USE_CUDA          = @MFEM_USE_CUDA@
 MFEM_USE_HIP           = @MFEM_USE_HIP@
 MFEM_USE_RAJA          = @MFEM_USE_RAJA@

config/defaults.mk

@@ -141,6 +141,7 @@ MFEM_USE_CONDUIT       = NO
 MFEM_USE_PUMI          = NO
 MFEM_USE_HIOP          = NO
 MFEM_USE_GSLIB         = NO
+MFEM_USE_PARRSB         = NO
 MFEM_USE_CUDA          = NO
 MFEM_USE_HIP           = NO
 MFEM_USE_RAJA          = NO
@@ -411,6 +412,11 @@ GSLIB_DIR = @MFEM_DIR@/../gslib/build
 GSLIB_OPT = -I$(GSLIB_DIR)/include
 GSLIB_LIB = -L$(GSLIB_DIR)/lib -lgs
 
+# PARRSB library
+PARRSB_DIR = @MFEM_DIR@/../parRSB/build
+PARRSB_OPT = -I$(PARRSB_DIR)/include
+PARRSB_LIB = -L$(PARRSB_DIR)/lib -lparRSB
+
 # CUDA library configuration
 CUDA_OPT =
 CUDA_LIB = -lcusparse

data/inline-hex.mesh

@@ -1,9 +1,9 @@
 MFEM INLINE mesh v1.0
 
 type = hex
-nx = 4
-ny = 4
-nz = 4
+nx = 6
+ny = 6
+nz = 6
 sx = 1.0
 sy = 1.0
 sz = 1.0

data/inline-quad.mesh

@@ -1,7 +1,7 @@
 MFEM INLINE mesh v1.0
 
 type = quad
-nx = 4
-ny = 4
+nx = 114
+ny = 114
 sx = 1.0
 sy = 1.0

data/inline-tri.mesh

@@ -1,7 +1,7 @@
 MFEM INLINE mesh v1.0
 
 type = tri
-nx = 4
-ny = 4
+nx = 5
+ny = 6
 sx = 1.0
 sy = 1.0

fem/gslib.cpp

@@ -23,6 +23,7 @@
 namespace gslib
 {
 #include "gslib.h"
+#include "parRSBCPP.h"
 }
 
 #ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
@@ -1000,4 +1001,125 @@ void OversetFindPointsGSLIB::Interpolate(const Vector &point_pos,
 
 } // namespace mfem
 
+
+#ifdef MFEM_USE_PARRSB
+
+namespace mfem
+{
+
+MFEMPARRSB::~MFEMPARRSB()
+{
+   //delete gsl_comm;
+   //delete cr;
+}
+
+void MFEMPARRSB::GetPartitioningSerial(MPI_Comm comm, Array<int> &partition)
+{
+    MFEM_VERIFY(serial, " Initialize with serial mesh.");
+    Array<long long> vertglob;
+    int geomtype = m->GetElementBaseGeometry(0);
+    int nvert;
+    if (geomtype == Geometry::SQUARE) {
+        nvert = 4;
+    }
+    else if (geomtype == Geometry::TRIANGLE) {
+        nvert = 3;
+    }
+    else if (geomtype == Geometry::CUBE) {
+        nvert = 8;
+    }
+    else if (geomtype == Geometry::TETRAHEDRON) {
+        nvert = 4;
+    }
+    else {
+        MFEM_ABORT(" Geometry type not currently supported.");
+    }
+    Vector coord(m->GetNE()*nvert*m->Dimension());
+
+    int n = 0;
+    for (int i = 0; i < m->GetNE(); i++) {
+        Array<int> vert;
+        m->GetElementVertices(i, vert);
+        for (int j = 0; j < vert.Size(); j++) {
+            vertglob.Append(vert[j]);
+            double *pos = m->GetVertex(vert[j]);
+            for (int d = 0; d < m->Dimension(); d++) {
+                coord[n++] = pos[d];
+            }
+        }
+    }
+
+
+    int myid;
+    MPI_Comm_rank(comm, &myid);
+    int nel = m->GetNE();
+    partition.SetSize(m->GetNE());
+    if (myid != 0) {
+        //vertglob.SetSize(0);
+        //coord.SetSize(0);
+        nel = 0;
+        //partition.SetSize(0);
+    }
+
+    gslib::parrsb_options options = gslib::parrsb_default_options;
+    parRSB_partMesh(partition, NULL, vertglob.GetData(),
+                    coord.GetData(), nel, nvert, options, comm);
+    partition.SetSize(m->GetNE());
+    MPI_Bcast(partition.GetData(), partition.Size(), MPI_INT, 0, comm);
+}
+
+void MFEMPARRSB::GetPartitioningParallel(Array<int> &partition)
+{
+    MFEM_VERIFY(!serial, " Initialize with parallel mesh.");
+    Array<long long> vertglob;
+    int geomtype = pm->GetElementBaseGeometry(0);
+    int nvert;
+    if (geomtype == Geometry::SQUARE) {
+        nvert = 4;
+    }
+    else if (geomtype == Geometry::TRIANGLE) {
+        nvert = 3;
+    }
+    else if (geomtype == Geometry::CUBE) {
+        nvert = 8;
+    }
+    else if (geomtype == Geometry::TETRAHEDRON) {
+        nvert = 4;
+    }
+    else {
+        MFEM_ABORT(" Geometry type not currently supported.");
+    }
+    Vector coord(pm->GetNE()*nvert*pm->Dimension());
+
+    Array<HYPRE_Int> gi;
+    pm->GetGlobalVertexIndices(gi);
+    int n = 0;
+    for (int i = 0; i < pm->GetNE(); i++) {
+        Array<int> vert;
+        pm->GetElementVertices(i, vert);
+        for (int j = 0; j < vert.Size(); j++) {
+            vertglob.Append(gi[vert[j]]);
+            double *pos = pm->GetVertex(vert[j]);
+            for (int d = 0; d < pm->Dimension(); d++) {
+                coord[n++] = pos[d];
+            }
+        }
+    }
+
+    partition.SetSize(pm->GetNE());
+    partition = -1;
+    gslib::parrsb_options options = gslib::parrsb_default_options;
+    parRSB_partMesh(partition, NULL, vertglob.GetData(),
+                    coord.GetData(), pm->GetNE(), nvert, options ,pm->GetComm());
+}
+
+void MFEMPARRSB::FreeData()
+{
+
+}
+
+}
+
+#endif // MFEM_USE_PARRSB
+
 #endif // MFEM_USE_GSLIB

fem/gslib.hpp

@@ -14,6 +14,7 @@
 
 #include "../config/config.hpp"
 #include "gridfunc.hpp"
+#include "../mesh/pmesh.hpp"
 
 #ifdef MFEM_USE_GSLIB
 
@@ -248,6 +249,40 @@ public:
 
 } // namespace mfem
 
+#ifdef MFEM_USE_PARRSB
+
+namespace mfem
+{
+
+class MFEMPARRSB
+{
+public:
+   Mesh *m;
+#ifdef MFEM_USE_MPI
+   ParMesh *pm;
+#endif
+   bool serial;
+
+public:
+   MFEMPARRSB(Mesh &mesh) : m(&mesh), serial(true) { };
+
+#ifdef MFEM_USE_MPI
+   MFEMPARRSB(ParMesh &pmesh) : pm(&pmesh), serial(false) { };
+#endif
+
+   ~MFEMPARRSB();
+
+   void GetPartitioningSerial(MPI_Comm comm, Array<int> &partition);
+   void GetPartitioningParallel(Array<int> &partition);
+
+   void FreeData();
+};
+
+
+}
+
+#endif // MFEM_USE_PARRSB
+
 #endif // MFEM_USE_GSLIB
 
 #endif // MFEM_GSLIB

makefile

@@ -273,7 +273,7 @@ endif
 
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
- SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP GSLIB\
+ SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP GSLIB PARRSB\
  OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER
 
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
@@ -337,7 +337,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU\
  MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC\
  MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT\
- MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
+ MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_PARRSB MFEM_USE_CUDA MFEM_USE_HIP\
  MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE MFEM_USE_SIMD\
  MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_AMGX MFEM_USE_MUMPS\
  MFEM_USE_CALIPER MFEM_SOURCE_DIR MFEM_INSTALL_DIR
@@ -668,6 +668,7 @@ status info:
 	$(info MFEM_USE_PUMI          = $(MFEM_USE_PUMI))
 	$(info MFEM_USE_HIOP          = $(MFEM_USE_HIOP))
 	$(info MFEM_USE_GSLIB         = $(MFEM_USE_GSLIB))
+	$(info MFEM_USE_PARRSB        = $(MFEM_USE_PARRSB))
 	$(info MFEM_USE_CUDA          = $(MFEM_USE_CUDA))
 	$(info MFEM_USE_HIP           = $(MFEM_USE_HIP))
 	$(info MFEM_USE_RAJA          = $(MFEM_USE_RAJA))

miniapps/gslib/parrsb.cpp

@@ -0,0 +1,225 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// Compile with: make parrsb
+//
+// Sample runs:
+//    mpirun -np 4 parrsb -m ../../data/inline-quad.mesh
+//    mpirun -np 4 parrsb -m ../../data/inline-hex.mesh
+//    mpirun -np 4 parrsb -m ../../data/armadillo.vtk
+
+#include "../../mfem.hpp"
+
+using namespace mfem;
+using namespace std;
+
+int main (int argc, char *argv[])
+{
+   // Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // Set the method's default parameters.
+   const char *mesh_file = "../../data/rt-2d-q3.mesh";
+   int order             = 3;
+   int mesh_poly_deg     = 3;
+   int rs_levels         = 0;
+   int rp_levels         = 0;
+   bool visualization    = false;
+   int fieldtype         = 0;
+   int ncomp             = 1;
+
+   // Parse command-line options.
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&mesh_poly_deg, "-mo", "--mesh-order",
+                  "Polynomial degree of mesh finite element space.");
+   args.AddOption(&rs_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&rp_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&fieldtype, "-ft", "--field-type",
+                  "Field type: 0 - H1, 1 - L2, 2 - H(div), 3 - H(curl).");
+   args.AddOption(&ncomp, "-nc", "--ncomp",
+                  "Number of components for H1 or L2 GridFunctions");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (myid == 0) { args.PrintOptions(cout); }
+
+   // Initialize and refine the starting mesh.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
+   for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
+   const int dim = mesh->Dimension();
+   if (myid == 0)
+   {
+      cout << "Mesh curvature of the original mesh: ";
+      if (mesh->GetNodes()) { cout << mesh->GetNodes()->OwnFEC()->Name(); }
+      else { cout << "(NONE)"; }
+      cout << endl;
+   }
+
+   // Mesh bounding box (for the full serial mesh).
+   Vector pos_min, pos_max;
+   MFEM_VERIFY(mesh_poly_deg > 0, "The order of the mesh must be positive.");
+   mesh->GetBoundingBox(pos_min, pos_max, mesh_poly_deg);
+   if (myid == 0)
+   {
+      cout << "--- Generating equidistant point for:\n"
+           << "x in [" << pos_min(0) << ", " << pos_max(0) << "]\n"
+           << "y in [" << pos_min(1) << ", " << pos_max(1) << "]\n";
+      if (dim == 3)
+      {
+         cout << "z in [" << pos_min(2) << ", " << pos_max(2) << "]\n";
+      }
+   }
+
+   // Distribute the mesh.
+
+   //mesh->EnsureNCMesh();
+   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
+   int nel = pmesh.GetNE();
+   int nel_min, nel_max;
+   MPI_Allreduce(&nel, &nel_min, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
+   MPI_Allreduce(&nel, &nel_max, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
+   if (myid == 0) {
+       std::cout << nel_max << " " << nel_min << " METIS-NELMAX-NELMIN\n";
+   }
+
+   // Curve the mesh based on the chosen polynomial degree.
+   L2_FECollection fecm(0, dim);
+   ParFiniteElementSpace pfespace(&pmesh, &fecm);
+   ParGridFunction x(&pfespace);
+   for (int i = 0; i < pmesh.GetNE(); i++) {
+       x(i) = myid;
+   }
+
+   if (visualization) {
+       char vishost[] = "localhost";
+       int  visport   = 19916;
+       socketstream sout;
+       sout.open(vishost, visport);
+       sout << "parallel " << num_procs << " " << myid << "\n";
+       sout.precision(8);
+       sout << "solution\n" << pmesh << x;
+       if (myid == 0)
+       {
+          sout << "window_title 'Default partitioning'\n"
+               << "window_geometry "
+               << 00 << " " << 0 << " " << 400 << " " << 400 << "\n"
+               << "keys jRmclA" << endl;
+       }
+   }
+
+   MFEMPARRSB parrsb = MFEMPARRSB(pmesh);
+   Array<int> partition;
+   parrsb.GetPartitioningParallel(partition);
+   for (int i = 0; i < pmesh.GetNE(); i++) {
+       x(i) = partition[i];
+   }
+
+   if (visualization) {
+       char vishost[] = "localhost";
+       int  visport   = 19916;
+       socketstream sout;
+       sout.open(vishost, visport);
+       sout << "parallel " << num_procs << " " << myid << "\n";
+       sout.precision(8);
+       sout << "solution\n" << pmesh << x;
+       if (myid == 0)
+       {
+          sout << "window_title 'ParRSB'\n"
+               << "window_geometry "
+               << 400 << " " << 0 << " " << 400 << " " << 400 << "\n"
+               << "keys jRmclA" << endl;
+       }
+   }
+
+
+   MFEMPARRSB parrsb2 = MFEMPARRSB(*mesh);
+   parrsb2.GetPartitioningSerial(MPI_COMM_WORLD, partition);
+   ParMesh pmesh2(MPI_COMM_WORLD, *mesh, partition);
+   ParFiniteElementSpace pfespace2(&pmesh2, &fecm);
+   ParGridFunction x2(&pfespace2);
+   for (int i = 0; i < pmesh2.GetNE(); i++) {
+       x2(i) = myid;
+   }
+
+   if (visualization) {
+       char vishost[] = "localhost";
+       int  visport   = 19916;
+       socketstream sout;
+       sout.open(vishost, visport);
+       sout << "parallel " << num_procs << " " << myid << "\n";
+       sout.precision(8);
+       sout << "solution\n" << pmesh2 << x2;
+       if (myid == 0)
+       {
+          sout << "window_title 'ParRSB Serial'\n"
+               << "window_geometry "
+               << 800 << " " << 0 << " " << 400 << " " << 400 << "\n"
+               << "keys jRmclA" << endl;
+       }
+   }
+
+
+   //mesh->EnsureNCMesh();
+   partition.SetSize(mesh->GetNE());
+   for (int i = 0; i < mesh->GetNE(); i++) {
+       partition[i] = i % num_procs;
+   }
+   ParMesh pmesh3(MPI_COMM_WORLD, *mesh, partition);
+   MFEMPARRSB parrsb3 = MFEMPARRSB(pmesh3);
+   parrsb3.GetPartitioningParallel(partition);
+   ParFiniteElementSpace pfespace3(&pmesh3, &fecm);
+   ParGridFunction x3(&pfespace3);
+   for (int i = 0; i < pmesh3.GetNE(); i++) {
+       x3(i) = partition[i];
+   }
+
+   if (visualization) {
+       char vishost[] = "localhost";
+       int  visport   = 19916;
+       socketstream sout;
+       sout.open(vishost, visport);
+       sout << "parallel " << num_procs << " " << myid << "\n";
+       sout.precision(8);
+       sout << "solution\n" << pmesh3 << x3;
+       if (myid == 0)
+       {
+          sout << "window_title 'Manual partition then parRSB'\n"
+               << "window_geometry "
+               << 400 << " " << 0 << " " << 400 << " " << 400 << "\n"
+               << "keys jRmclA" << endl;
+       }
+   }
+
+
+
+
+
+
+
+
+
+
+   MPI_Finalize();
+
+
+
+
+
+   return 0;
+}