diff --git a/CMakeLists.txt b/CMakeLists.txt
index 80681c5e74e6e790604bd271e4285099dcac5bab..5415ec4976b207f91f3562dd889cc7c21509a14a 100644
--- a/CMakeLists.txt
+++ b/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
diff --git a/doc/changelog.md b/doc/changelog.md
index d3defac113566debd48c73ef83f93087fa36730b..a35cb5cc01ff7f2cc38ea75383b7f622f0ada7ac 100644
--- a/doc/changelog.md
+++ b/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
+
diff --git a/include/gmsh_io.h b/include/gmsh_io.h
index 0824294ad5f41df66b74ac9dd2fd91aacb61348b..8d32796a4be30b0230a54b18f343e721233a7ab2 100644
--- a/include/gmsh_io.h
+++ b/include/gmsh_io.h
@@ -58,11 +58,13 @@ struct boundary_descriptor
class model
{
- int geometric_order;
- int approximation_order;
+ 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;
diff --git a/include/gmsh_mpi.h b/include/gmsh_mpi.h
new file mode 100644
index 0000000000000000000000000000000000000000..c428de0fa77eab2de162603701b2e44202708e97
--- /dev/null
+++ b/include/gmsh_mpi.h
@@ -0,0 +1,20 @@
+#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);
+
diff --git a/src/gmsh_io.cpp b/src/gmsh_io.cpp
index bdddb6e911c5067c8cb38389327dacfd209962a6..7d8eebee3f68cd0597ef6410246627ced6bd6ba9 100644
--- a/src/gmsh_io.cpp
+++ b/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)
{
+ 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;
+
+ /* 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;
+ }
+
std::vector<int> elemTypes;
gmm::getElementTypes(elemTypes, dim, tag);
assert(elemTypes.size() == 1);
+
+ /* Get all the face keys for a given boundary. */
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
-{
- 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();
}
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
{
diff --git a/src/gmsh_mpi.cpp b/src/gmsh_mpi.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f644e3af0450807d511dc70507f5ee3482a9f8fc
--- /dev/null
+++ b/src/gmsh_mpi.cpp
@@ -0,0 +1,92 @@
+#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
diff --git a/src/maxwell_solver.cpp b/src/maxwell_solver.cpp
index 72ac3c5b636055881a8f57b437afcc12a013c5d0..ce3cd72d56c1847014b5f553b6a7aedf712ad9f0 100644
--- a/src/maxwell_solver.cpp
+++ b/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;
}
diff --git a/src/testmod.cpp b/src/testmod.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..74f9ce12ca4520abdb7ac33977c9b4c55a26f745
--- /dev/null
+++ b/src/testmod.cpp
@@ -0,0 +1,74 @@
+#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
diff --git a/src/testmpi.cpp b/src/testmpi.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bcfe50692465a49c88c220b98214dfa47032000d
--- /dev/null
+++ b/src/testmpi.cpp
@@ -0,0 +1,39 @@
+#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;
+}
diff --git a/src/testpart.cpp b/src/testpart.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ff467d3b0c0ba75cc38e36d040b9b7cfcc97a74d
--- /dev/null
+++ b/src/testpart.cpp
@@ -0,0 +1,44 @@
+#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