Keeping track of parent tag in boundary_descriptor.

CMakeLists.txt

@@ -292,10 +292,6 @@ 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)
 
@@ -316,9 +312,6 @@ 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)
 

include/gmsh_io.h

@@ -44,7 +44,8 @@ enum class face_type : int
 struct boundary_descriptor
 {
     face_type   type;
-    int         gmsh_entity;
+    int         gmsh_entity;        /* Actual GMSH entity */
+    int         parent_entity;      /* Parent GMSH entity in partitioned model */
 
     //auto operator <=> (const boundary_descriptor&) const = default;
 
@@ -60,8 +61,16 @@ struct boundary_descriptor
     }
 
     boundary_descriptor()
-        : type(face_type::NONE), gmsh_entity(0)
+        : type(face_type::NONE), gmsh_entity(0), parent_entity(-1)
     {}
+
+    int material_tag(void) const
+    {
+        if (parent_entity == -1)
+            return gmsh_entity;
+
+        return parent_entity;
+    }
 };
 
 
@@ -85,6 +94,9 @@ class model
 
     /* Map from the partition number to the entities belonging to that partition */
     std::map<int, std::vector<int>>             partition_map;
+    
+    /* In a partitioned model, map from 2D entity tag to parent tag */
+    std::map<int, int>                          surface_to_parent_map;
 
     /* Map from 3D entity tag to partition */
     std::map<int, int>                          partition_inverse_map;
@@ -107,7 +119,7 @@ class model
     void    populate_from_gmsh_rankN(void);
     void    import_gmsh_entities_rankN(void);
     int     my_partition(void) const;
-#endif
+#endif /* USE_MPI */
 
     void    make_boundary_to_faces_map(void);
     void    make_partition_to_entities_map(void);

include/gmsh_mpi.h

@@ -100,6 +100,46 @@ priv_MPI_Recv(T& elem, int src, MPI_Comm comm)
     MPI_Recv(&elem, sizeof(T), MPI_PACKED, src, 0, comm, &status);
 }
 
+template<typename T1, typename T2>
+void
+priv_MPI_Send(std::map<T1, T2>& map, int dst, MPI_Comm comm)
+{
+    static_assert(std::is_trivially_copyable<T1>::value, "Type must be trivially copyable");
+    static_assert(std::is_trivially_copyable<T2>::value, "Type must be trivially copyable");
+
+    std::vector<T1> left;
+    left.reserve( map.size() );
+    std::vector<T2> right;
+    right.reserve( map.size() );
+
+    for (auto& m : map)
+    {
+        left.push_back(m.first);
+        right.push_back(m.second);
+    }
+
+    priv_MPI_Send(left, dst, comm);
+    priv_MPI_Send(right, dst, comm);
+}
+
+template<typename T1, typename T2>
+void
+priv_MPI_Recv(std::map<T1, T2>& map, int src, MPI_Comm comm)
+{
+    static_assert(std::is_trivially_copyable<T1>::value, "Type must be trivially copyable");
+    static_assert(std::is_trivially_copyable<T2>::value, "Type must be trivially copyable");
+
+    std::vector<T1> left;
+    priv_MPI_Recv(left, src, comm);
+
+    std::vector<T2> right;
+    priv_MPI_Recv(right, src, comm);
+
+    assert(left.size() == right.size());
+
+    for (size_t i = 0; i < left.size(); i++)
+        map[ left[i] ] = right[i]; 
+}
 
 template<typename T>
 void
@@ -166,3 +206,43 @@ priv_MPI_Bcast(std::map<T1, std::vector<T2>>& map, int root, MPI_Comm comm)
     }
 }
 
+template<typename T1, typename T2>
+void
+priv_MPI_Bcast(std::map<T1, T2>& map, int root, MPI_Comm comm)
+{
+    static_assert(std::is_trivially_copyable<T1>::value, "Type must be trivially copyable");
+    static_assert(std::is_trivially_copyable<T2>::value, "Type must be trivially copyable");
+
+    int rank;
+    MPI_Comm_rank(comm, &rank);
+
+    if (rank == root)
+    {
+        std::vector<T1> left;
+        left.reserve( map.size() );
+        std::vector<T2> right;
+        right.reserve( map.size() );
+
+        for (auto& m : map)
+        {
+            left.push_back(m.first);
+            right.push_back(m.second);
+        }
+
+        priv_MPI_Bcast(left, root, comm);
+        priv_MPI_Bcast(right, root, comm);
+    }
+    else
+    {
+        std::vector<T1> left;
+        priv_MPI_Bcast(left, root, comm);
+
+        std::vector<T2> right;
+        priv_MPI_Bcast(right, root, comm);
+
+        assert(left.size() == right.size());
+
+        for (size_t i = 0; i < left.size(); i++)
+            map[ left[i] ] = right[i]; 
+    }
+}
\ No newline at end of file

src/gmsh_io.cpp

@@ -87,6 +87,7 @@ model::make_boundary_to_faces_map(void)
             std::vector<int> parts;
             gm::getPartitions(dim, tag, parts);
             comm_map[tag] = parts;
+            surface_to_parent_map[tag] = ptag;
         }
 
         std::vector<int> elemTypes;
@@ -300,7 +301,10 @@ model::import_gmsh_entities_rankN(void)
 bool
 model::is_interprocess_boundary(int tag)
 {
-    return comm_map.at(tag).size() == 2;
+    if (num_partitions > 1)
+        return comm_map.at(tag).size() == 2;
+    
+    return false;
 }
 
 #endif /* USE_MPI */
@@ -361,6 +365,8 @@ model::map_boundaries(void)
             }
             /* and store the gmsh tag in the descriptor. */
             bd.gmsh_entity = (*itor).second;
+            if (num_partitions > 1)
+                bd.parent_entity = surface_to_parent_map.at( (*itor).second );
 
             /* Finally, put the descriptor in the global array of faces. */
             bnd_descriptors.at(fbase + iF) = bd;
@@ -416,6 +422,7 @@ model::populate_from_gmsh_rank0()
     priv_MPI_Bcast(boundary_map, 0, MPI_COMM_WORLD);
     priv_MPI_Bcast(comm_map, 0, MPI_COMM_WORLD);
     priv_MPI_Bcast(partition_map, 0, MPI_COMM_WORLD);
+    priv_MPI_Bcast(surface_to_parent_map, 0, MPI_COMM_WORLD);
 
     for (auto& [tag, vec] : partition_map)
     {
@@ -443,6 +450,7 @@ model::populate_from_gmsh_rankN()
     priv_MPI_Bcast(boundary_map, 0, MPI_COMM_WORLD);
     priv_MPI_Bcast(comm_map, 0, MPI_COMM_WORLD);
     priv_MPI_Bcast(partition_map, 0, MPI_COMM_WORLD);
+    priv_MPI_Bcast(surface_to_parent_map, 0, MPI_COMM_WORLD);
 
     for (auto& [tag, vec] : partition_map)
     {

src/gmsh_mpi.cpp

-#include "gmsh_mpi.h"
-
-#if 0
-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);
-    }
-}
-
-#endif
\ No newline at end of file

src/maxwell_solver.cpp

@@ -5,10 +5,6 @@
 #include <filesystem>
 #endif
 
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
 #ifdef USE_MPI
 #include <mpi/mpi.h>
 #endif
@@ -203,10 +199,6 @@ void test_it(const model& mod, State& state, maxwell::parameter_loader& mpl)
     auto cycle_print_rate = mpl.postpro_cyclePrintRate();
     auto ti = mpl.sim_timeIntegrator();
 
-#ifdef _OPENMP
-    omp_set_num_threads(4);
-#endif
-
     std::cout << "    BEGINNING SIMULATION" << std::endl;
     std::cout << "I will do " << num_timesteps << " of " << state.delta_t;
     std::cout << "s each." << std::endl;

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;
-}