Modifying model to allow multiple partitions.

CMakeLists.txt

@@ -71,6 +71,13 @@ set(SOL2_BUILD_LUA OFF CACHE BOOL "Override sol2 build lua option")
 add_subdirectory(contrib/sol2)
 include_directories(contrib/sol2/include)
 
+option(OPT_USE_MPI "Use MPI if available" OFF)
+find_package(MPI)
+if (OPT_USE_MPI AND MPI_FOUND)
+    add_definitions(-DUSE_MPI)
+    include_directories("${MPI_CXX_INCLUDE_DIRS}")
+    set(LINK_LIBS ${LINK_LIBS} ${MPI_CXX_LIBRARIES})
+endif()
 ######################################################################
 ## Use or not use GPU solvers
 option(OPT_ENABLE_GPU_SOLVER "Enable GPU solvers" OFF)
@@ -285,6 +292,10 @@ set(LIBGMSHDG_SOURCES   src/gmsh_io.cpp
                         src/param_loader.cpp
                     )
 
+if (OPT_USE_MPI)
+    set(LIBGMSHDG_SOURCES ${LIBGMSHDG_SOURCES} src/gmsh_mpi.cpp)
+endif()
+
 if (OPT_ENABLE_GPU_SOLVER)
     set(LIBGMSHDG_SOURCES ${LIBGMSHDG_SOURCES} src/kernels_gpu_mi.cpp)
 
@@ -303,6 +314,16 @@ else()
     target_link_libraries(gmshdg ${LINK_LIBS})
 endif()
 
+########## mpi testing stuff
+
+add_executable(testmpi src/testmpi.cpp)
+target_link_libraries(testmpi gmshdg)
+
+add_executable(testpart src/testpart.cpp)
+target_link_libraries(testpart gmshdg)
+
+add_executable(testmod src/testmod.cpp)
+target_link_libraries(testmod gmshdg)
 
 ########## maxwell solver
 set(MAXWELL_SOURCES     src/maxwell_interface.cpp

doc/changelog.md

@@ -23,5 +23,10 @@ Unavailable features:
 
 ## Version v0.3 (NOT RELEASED YET)
 - Volumetric sources on CPU & GPU (2ed2e9e7).
+- Fixed issues that prevented compilation on CECI clusters.
+- Minor performance optimizations
+- Line integration
+- Preliminary MPI support
+
 
 

include/gmsh_io.h

@@ -60,9 +60,11 @@ class model
 {
     int                                         geometric_order;
     int                                         approximation_order;
+    int                                         num_partitions;
 
     std::vector<entity>                         entities;
     std::map<size_t, std::vector<element_key>>  boundary_map;
+    std::map<size_t, std::vector<int>>          comm_map;
     std::vector<boundary_descriptor>            bnd_descriptors;
     element_connectivity<element_key>           conn;
 
@@ -73,7 +75,13 @@ class model
     void    import_gmsh_entities(void);
     void    update_connectivity(const entity&, size_t);
 
-    void    build_priv(void);
+    void    populate_from_gmsh_rank0(void);
+
+#ifdef HAVE_MPI
+    void    populate_from_gmsh_rankN(void);
+#endif
+
+    void    make_boundary_to_faces_map(void);
 
 public:
     model();
@@ -85,6 +93,11 @@ public:
     void build(void);
     void build(double);
 
+    void generate_mesh(void);
+    void generate_mesh(double);
+    void partition_mesh(int);
+    void populate_from_gmsh(void);
+
     const element_connectivity<element_key>& connectivity() const
     {
         return conn;

include/gmsh_mpi.h

+#pragma once
+
+#include <mpi/mpi.h>
+#include "gmsh.h"
+
+namespace g     = gmsh;
+namespace go    = gmsh::option;
+namespace gm    = gmsh::model;
+namespace gmm   = gmsh::model::mesh;
+namespace gmo   = gmsh::model::occ;
+namespace gmg   = gmsh::model::geo;
+namespace gv    = gmsh::view;
+
+using gvp_t = gmsh::vectorpair;
+
+/* Macros used to annotate GMSH integer inputs */
+#define DIMENSION(x) x
+
+void MPIGMSH_getEntities(gvp_t&, int, MPI_Comm);
+

src/gmsh_io.cpp

@@ -5,6 +5,10 @@
 #include "gmsh_io.h"
 #include "entity.h"
 
+#ifdef HAVE_MPI
+#include <mpi/mpi.h>
+#endif /* HAVE_MPI */
+
 std::string
 quadrature_name(int order)
 {
@@ -34,40 +38,26 @@ basis_grad_name(int order)
 }
 
 model::model()
-    : geometric_order(1), approximation_order(1)
+    : geometric_order(1), approximation_order(1), num_partitions(1)
 {}
 
 model::model(int pgo, int pao)
-    : geometric_order(pgo), approximation_order(pao)
+    : geometric_order(pgo), approximation_order(pao), num_partitions(1)
 {
     assert(geometric_order > 0);
     assert(approximation_order > 0);
 }
 
-model::model(const char *filename)
-    : geometric_order(1), approximation_order(1)
-{
-    gmsh::open( std::string(filename) );
-}
-
-model::model(const char *filename, int pgo, int pao)
-    : geometric_order(pgo), approximation_order(pao)
-{
-    assert(geometric_order > 0);
-    assert(approximation_order > 0);
-
-    gmsh::open( std::string(filename) );
-}
-
 model::~model()
 {
-    gmsh::clear();
+    //gmsh::clear();
 }
 
+
 void
-model::build_priv()
+model::make_boundary_to_faces_map(void)
 {
-    gmm::setOrder( geometric_order );
+    assert(num_partitions > 1);
 
     gvp_t gmsh_entities;
 
@@ -76,124 +66,38 @@ model::build_priv()
      * because otherwise you get spurious 2D entities (and this in
      * turn because the constructor of entity calls addDiscreteEntity()
      * to get all the element faces) .*/
-    std::map<size_t, std::vector<face_key>>   boundaries_elemTags;
     gm::getEntities(gmsh_entities, DIMENSION(2));
     for (auto [dim, tag] : gmsh_entities)
     {
-        std::vector<int> elemTypes;
-        gmm::getElementTypes(elemTypes, dim, tag);
-        assert(elemTypes.size() == 1);
-        for (auto& elemType : elemTypes)
-        {
-            element_key_factory ekf(DIMENSION(2), tag, elemType);
-            boundary_map[tag] = ekf.get_keys();
-        }
-    }
-
-    entities.clear();
-    import_gmsh_entities();
-    map_boundaries(); /* This must happen after import_gmsh_entities(). */
-}
-
-void
-model::build()
-{
-    gmm::generate( DIMENSION(3) );
-    build_priv();
-}
-
-void
-model::build(double sf)
-{
-    gmm::generate( DIMENSION(3) );
-
-    std::vector<double> tr = {sf,  0,  0,  0,
-                               0, sf,  0,  0,
-                               0,  0, sf,  0 };
-    gmm::affineTransform(tr);
-    build_priv();
-}
-
-std::vector<entity>::const_iterator
-model::begin() const
+        if (num_partitions > 1)
         {
-    return entities.begin();
-}
-
-std::vector<entity>::const_iterator
-model::end() const
-{
-    return entities.end();
-}
-
-std::vector<entity>::iterator
-model::begin()
-{
-    return entities.begin();
-}
-
-std::vector<entity>::iterator
-model::end()
-{
-    return entities.end();
-}
+            /* If the model is partitioned and entity has no parents,
+             * it means that there are no elements in the entity on which
+             * we will compute. Skip it. */
+            int pdim, ptag;
+            gmsh::model::getParent(dim, tag, pdim, ptag);
 
-entity&
-model::entity_at(size_t pos)
-{
-    return entities.at(pos);
-}
-
-const entity&
-model::entity_at(size_t pos) const
-{
-    return entities.at(pos);
-}
-
-size_t
-model::num_dofs(void) const
-{
-    size_t ret = 0;
-    for (const auto& e : entities)
-        ret += e.num_dofs();
+            if (ptag == -1)
+                continue;
 
-    return ret;
+            /* This stores which partitions share a boundary. It will be
+             * used to determine who communicates with who. */
+            std::vector<int> parts;
+            gm::getPartitions(dim, tag, parts);
+            comm_map[tag] = parts;
         }
 
-size_t
-model::num_fluxes(void) const
-{
-    size_t ret = 0;
-    for (const auto& e : entities)
-        ret += e.num_fluxes();
-
-    return ret;
-}
+        std::vector<int> elemTypes;
+        gmm::getElementTypes(elemTypes, dim, tag);
+        assert(elemTypes.size() == 1);
 
-size_t
-model::num_cells(void) const
+        /* Get all the face keys for a given boundary. */
+        for (auto& elemType : elemTypes)
         {
-    size_t ret = 0;
-    for (const auto& e : entities)
-        ret += e.num_cells();
-
-    return ret;
+            element_key_factory ekf(DIMENSION(2), tag, elemType);
+            boundary_map[tag] = ekf.get_keys();
         }
-
-size_t
-model::num_faces(void) const
-{
-    size_t ret = 0;
-    for (const auto& e : entities)
-        ret += e.num_faces();
-
-    return ret;
     }
-
-size_t
-model::num_entities(void) const
-{
-    return entities.size();
 }
 
 void
@@ -223,6 +127,18 @@ model::import_gmsh_entities(void)
     size_t index_base = 0;
     for (auto [dim, tag] : gmsh_entities)
     {
+        if (num_partitions > 1)
+        {
+            /* If the model is partitioned and entity has no parents,
+             * it means that there are no elements in the entity on which
+             * we will compute. Skip it. */
+            int pdim, ptag;
+            gmsh::model::getParent(dim, tag, pdim, ptag);
+
+            if (ptag == -1)
+                continue;
+        }
+
         std::vector<int> elemTypes;
         gmm::getElementTypes(elemTypes, dim, tag);
         
@@ -330,6 +246,171 @@ model::map_boundaries(void)
     */
 }
 
+void
+model::populate_from_gmsh_rank0()
+{
+    gmm::setOrder( geometric_order );
+
+    make_boundary_to_faces_map();
+
+    for (auto& [tag, parts] : comm_map)
+    {
+        std::cout << tag << ": ";
+        for (auto& part : parts)
+            std::cout << part << " ";
+        std::cout << std::endl;
+    }
+
+    entities.clear();
+    import_gmsh_entities();
+    map_boundaries(); /* This must happen after import_gmsh_entities(). */
+}
+
+#ifdef HAVE_MPI
+void
+model::populate_from_gmsh_rankN()
+{}
+#endif /* HAVE_MPI */
+
+void
+model::populate_from_gmsh()
+{
+#ifdef HAVE_MPI
+    int rank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    if (rank == 0)
+        populate_from_gmsh_rank0();
+    else
+        populate_from_gmsh_rankN();
+#else /* HAVE_MPI */
+    populate_from_gmsh_rank0();
+#endif /* HAVE_MPI */
+}
+
+void
+model::build()
+{
+    generate_mesh();
+    make_boundary_to_faces_map();
+    populate_from_gmsh_rank0();
+}
+
+void
+model::build(double sf)
+{
+    generate_mesh(sf);
+    make_boundary_to_faces_map();
+    populate_from_gmsh_rank0();
+}
+
+void
+model::generate_mesh()
+{
+    gmm::generate( DIMENSION(3) );
+}
+
+void
+model::generate_mesh(double sf)
+{
+    gmm::generate( DIMENSION(3) );
+
+    std::vector<double> tr = {sf,  0,  0,  0,
+                               0, sf,  0,  0,
+                               0,  0, sf,  0 };
+    gmm::affineTransform(tr);
+}
+
+void
+model::partition_mesh(int parts)
+{
+    gmm::partition(parts);
+    num_partitions = parts;
+}
+
+std::vector<entity>::const_iterator
+model::begin() const
+{
+    return entities.begin();
+}
+
+std::vector<entity>::const_iterator
+model::end() const
+{
+    return entities.end();
+}
+
+std::vector<entity>::iterator
+model::begin()
+{
+    return entities.begin();
+}
+
+std::vector<entity>::iterator
+model::end()
+{
+    return entities.end();
+}
+
+entity&
+model::entity_at(size_t pos)
+{
+    return entities.at(pos);
+}
+
+const entity&
+model::entity_at(size_t pos) const
+{
+    return entities.at(pos);
+}
+
+size_t
+model::num_dofs(void) const
+{
+    size_t ret = 0;
+    for (const auto& e : entities)
+        ret += e.num_dofs();
+
+    return ret;
+}
+
+size_t
+model::num_fluxes(void) const
+{
+    size_t ret = 0;
+    for (const auto& e : entities)
+        ret += e.num_fluxes();
+
+    return ret;
+}
+
+size_t
+model::num_cells(void) const
+{
+    size_t ret = 0;
+    for (const auto& e : entities)
+        ret += e.num_cells();
+
+    return ret;
+}
+
+size_t
+model::num_faces(void) const
+{
+    size_t ret = 0;
+    for (const auto& e : entities)
+        ret += e.num_faces();
+
+    return ret;
+}
+
+size_t
+model::num_entities(void) const
+{
+    return entities.size();
+}
+
+
+
 model::entofs_pair
 model::lookup_tag(size_t tag) const
 {

src/gmsh_mpi.cpp

+#include "gmsh_mpi.h"
+
+template<typename T1, typename T2>
+void zip(const std::vector<T1>& firsts, const std::vector<T2>& seconds,
+    std::vector<std::pair<T1, T2>>& pairs)
+{
+    pairs.clear();
+    auto len = std::min(firsts.size(), seconds.size());
+    for(size_t i = 0; i < len; i++)
+        pairs.push_back( std::make_pair(firsts[i], seconds[i]) );
+}
+
+template<typename T1, typename T2>
+void unzip(const std::vector<std::pair<T1, T2>>& pairs,
+    std::vector<T1>& firsts, std::vector<T2>& seconds)
+{
+    firsts.clear();
+    seconds.clear();
+    for (const auto& [fst, snd] : pairs)
+    {
+        firsts.push_back(fst);
+        seconds.push_back(snd);
+    }
+}
+
+void
+MPIGMSH_getEntities(gvp_t& gmsh_entities, int dim, MPI_Comm comm)
+{
+    int rank, size;
+    MPI_Comm_rank(comm, &rank);
+    MPI_Comm_size(comm, &size);
+
+    if (rank == 0)
+    {
+        gm::getEntities(gmsh_entities, dim);
+        size_t num_entities = gmsh_entities.size();
+
+        if (size > 1)
+        {
+            MPI_Bcast(&num_entities, 1, MPI_UNSIGNED_LONG_LONG, 0, comm);
+            std::vector<int> dims;
+            std::vector<int> tags;
+            unzip(gmsh_entities, dims, tags);
+            MPI_Bcast(dims.data(), num_entities, MPI_INT, 0, comm);
+            MPI_Bcast(tags.data(), num_entities, MPI_INT, 0, comm);
+        }
+    }
+    else
+    {
+        size_t num_entities;
+        MPI_Bcast(&num_entities, 1, MPI_UNSIGNED_LONG_LONG, 0, comm);
+        
+        std::vector<int> dims;
+        dims.resize(num_entities);
+        MPI_Bcast(dims.data(), num_entities, MPI_INT, 0, comm);
+
+        std::vector<int> tags;
+        tags.resize(num_entities);
+        MPI_Bcast(tags.data(), num_entities, MPI_INT, 0, comm);
+
+        zip(dims, tags, gmsh_entities);
+    }
+
+    return;
+}
+
+void
+MPIGMSH_getElementTypes(std::vector<int>& elemTypes, int dim, int tag, MPI_Comm comm)
+{
+    int rank, size;
+    MPI_Comm_rank(comm, &rank);
+    MPI_Comm_size(comm, &size);
+
+    if (rank == 0)
+    {
+        //gm::getEntities(elemTypes, dim, tag);
+        size_t num_eT = elemTypes.size();
+
+        if (size > 1)
+        {
+            MPI_Bcast(&num_eT, 1, MPI_UNSIGNED_LONG_LONG, 0, comm);
+            MPI_Bcast(elemTypes.data(), num_eT, MPI_INT, 0, comm);
+        }
+    }
+    else
+    {
+        size_t num_eT;
+        MPI_Bcast(&num_eT, 1, MPI_UNSIGNED_LONG_LONG, 0, comm);
+        elemTypes.resize(num_eT);
+        MPI_Bcast(elemTypes.data(), num_eT, MPI_INT, 0, comm);
+    }
+}
\ No newline at end of file

src/maxwell_solver.cpp

@@ -9,6 +9,10 @@
 #include <omp.h>
 #endif
 
+#ifdef USE_MPI
+#include <mpi/mpi.h>
+#endif
+
 #include <fstream>
 
 #include "gmsh.h"
@@ -275,10 +279,14 @@ register_lua_usertypes_bystate(maxwell::parameter_loader& mpl,
 
 int main(int argc, const char *argv[])
 {
-    gmsh::initialize();
-    //make_geometry(1, 0.1);
+#ifdef USE_MPI
+    MPI_Init(&argc, &argv);
 
-    gmsh::option::setNumber("General.Terminal", 0);
+    int rank, size;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    MPI_Comm_size(MPI_COMM_WORLD, &size);
+    std::cout << "Running with MPI: " << rank+1 << "/" << size << std::endl;
+#endif
 
     if (argc != 2)
     {
@@ -301,6 +309,8 @@ int main(int argc, const char *argv[])
     std::filesystem::create_directory( mpl.sim_name() );;
 #endif
 
+    gmsh::initialize();
+    gmsh::option::setNumber("General.Terminal", 0);
     gmsh::open( mpl.sim_gmshmodel() );
 
     model mod( mpl.sim_geomorder(), mpl.sim_approxorder() );
@@ -336,5 +346,9 @@ int main(int argc, const char *argv[])
 
     //gmsh::finalize();
 
+#ifdef USE_MPI
+    MPI_Finalize();
+#endif
+
     return 0;
 }

src/testmod.cpp

+#include <iostream>
+
+#ifdef USE_MPI
+#include <mpi/mpi.h>
+#endif
+
+#include "gmsh_io.h"
+
+std::pair<int, int>
+priv_MPI_Procinfo()
+{
+    int rank, size;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    MPI_Comm_size(MPI_COMM_WORLD, &size);
+    return std::make_pair(rank, size);
+}
+
+void
+init_model_rank0(model& mod)
+{
+    std::cout << "Rank 0 init" << std::endl;
+
+    auto [comm_rank, comm_size] = priv_MPI_Procinfo();
+
+    mod.generate_mesh();
+    if (comm_size > 1)
+        mod.partition_mesh(comm_size);
+
+    mod.populate_from_gmsh();
+
+}
+
+void
+init_model_rankN(model& mod)
+{
+    std::cout << "Rank N init" << std::endl;
+    mod.populate_from_gmsh();
+}
+
+int main(int argc, char **argv)
+{
+    if (argc < 2)
+    {
+        std::cout << "Please specify .geo" << std::endl;
+        return 1;
+    }
+
+    MPI_Init(&argc, &argv);
+
+    auto [rank, size] = priv_MPI_Procinfo();
+
+    model mod;
+
+    if (rank == 0)
+    {
+        gmsh::initialize();
+        gmsh::open(argv[1]);
+        init_model_rank0(mod);
+    }
+    else
+    {
+        init_model_rankN(mod);
+    }
+
+
+
+    if (rank == 0)
+    {
+        gmsh::finalize();
+    }
+
+    MPI_Finalize();
+    return 0;
+}
\ No newline at end of file

src/testmpi.cpp

+#include <iostream>
+#include <mpi/mpi.h>
+
+#include "gmsh_mpi.h"
+
+int main(int argc, char **argv)
+{
+    MPI_Init(&argc, &argv);
+
+    int rank, size;
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    MPI_Comm_size(MPI_COMM_WORLD, &size);
+    std::cout << "Running with MPI: " << rank+1 << "/" << size << std::endl;
+
+    gvp_t gmsh_entities;
+
+    if (rank == 0)
+    {
+        gmsh::initialize();
+        gmsh::open("yagi_prof.geo");
+        gmm::generate(3);
+        gmm::setOrder(1);
+    }
+
+    MPIGMSH_getEntities(gmsh_entities, 3, MPI_COMM_WORLD);
+
+    std::cout << "Rank: " << rank << " -> ";
+    for (auto& [dim, tag] : gmsh_entities)
+        std::cout << "(" << dim << ", " << tag << ") ";
+    std::cout << std::endl;
+
+    if (rank == 0)
+    {
+        //gmsh::finalize();
+    }
+
+    MPI_Finalize();
+    return 0;
+}

src/testpart.cpp

+#include <iostream>
+
+#include "gmsh.h"
+
+int main(int argc, char **argv)
+{
+    if (argc < 2)
+    {
+        std::cout << "Specify .geo" << std::endl;
+        return 1;
+    }
+
+    gmsh::initialize();
+    gmsh::open(argv[1]);
+
+    gmsh::vectorpair gmsh_entities;
+
+    gmsh::model::getEntities(gmsh_entities, 2);
+
+    for (auto [dim, tag] : gmsh_entities)
+    {
+        //std::vector<int> elemTypes;
+        //gmm::getElementTypes(elemTypes, dim, tag);
+
+        std::string etype;
+        gmsh::model::getType(dim, tag, etype);
+        std::cout << etype << ", tag: " << tag << std::endl;
+
+        int pdim, ptag;
+        gmsh::model::getParent(dim, tag, pdim, ptag);
+        std::cout << "  Parent: " << ptag << std::endl;
+
+        std::vector<int> parts;
+        gmsh::model::getPartitions(dim, tag, parts);
+        std::cout << "  Partitions: ";
+        for (const auto& part : parts)
+            std::cout << part << " ";
+        std::cout << std::endl;
+    }
+
+    gmsh::finalize();
+
+    return 0;
+}
\ No newline at end of file