From e91aeb35f835d7e3d9cbef47838621502c2874b0 Mon Sep 17 00:00:00 2001
From: Matteo Cicuttin <datafl4sh@toxicnet.eu>
Date: Tue, 31 Aug 2021 17:06:22 +0200
Subject: [PATCH] Interprocess connectivity.
---
include/connectivity.h | 96 ++++++++++++++++++++++++++++++++++++++++++
include/gmsh_io.h | 9 ++++
src/connectivity.cpp | 21 +++++++++
src/gmsh_io.cpp | 83 +++++++++++++++++++++++++++++++++---
src/maxwell_common.cpp | 86 +++++++++++++++++++++++++------------
5 files changed, 263 insertions(+), 32 deletions(-)
diff --git a/include/connectivity.h b/include/connectivity.h
index 275ca63..23a4fa2 100644
--- a/include/connectivity.h
+++ b/include/connectivity.h
@@ -6,8 +6,13 @@
#include <map>
#include <set>
+#include "point.hpp"
#include "sgr.hpp"
+#ifdef USE_MPI
+#include "gmsh_mpi.h"
+#endif /* USE_MPI */
+
/***************************************************************************/
struct neighbour_descriptor
{
@@ -22,6 +27,25 @@ struct neighbour_descriptor
std::ostream& operator<<(std::ostream&, const neighbour_descriptor&);
+struct neighbour_descriptor_interprocess
+{
+ size_t iE; /* Entity offset in data structures */
+ size_t iT; /* Element offset in entity in data structures */
+ size_t via_iF; /* Element face index (i.e. 0,1,2 for triangles) */
+
+ int parent_tag; /* Tag of the parent entity */
+ int tag; /* Tag of the actual entity */
+ int partition; /* Partition on which the element resides */
+
+ point_3d bar; /* Barycenter, for materials. */
+
+ neighbour_descriptor_interprocess();
+ bool operator<(const neighbour_descriptor_interprocess&) const;
+
+};
+
+std::ostream& operator<<(std::ostream&, const neighbour_descriptor_interprocess&);
+
/***************************************************************************/
struct cell_descriptor
{
@@ -235,3 +259,75 @@ public:
}
};
+
+
+
+/***************************************************************************/
+#ifdef USE_MPI
+template<typename FaceKey>
+class interprocess_connectivity
+{
+public:
+ using ConnData = std::pair<FaceKey, neighbour_descriptor_interprocess>;
+ using NeighPair = neighbour_pair<neighbour_descriptor_interprocess>;
+
+private:
+ std::map<FaceKey, NeighPair> face_cell_adj;
+
+public:
+ interprocess_connectivity()
+ {}
+
+ void
+ connect(const FaceKey& fk, const neighbour_descriptor_interprocess& nd)
+ {
+ /* Insert the information that one the two neighbours around
+ * face `e' is `nd'. */
+ auto& nr = face_cell_adj[fk];
+ nr.add_neighbour(nd);
+ }
+
+ /* Find the neighbour of the element (iE,iT) via face fk */
+ std::pair<neighbour_descriptor_interprocess, bool>
+ neighbour_via(size_t iE, size_t iT, const FaceKey& fk) const
+ {
+ auto itor = face_cell_adj.find(fk);
+ if (itor == face_cell_adj.end())
+ return std::make_pair(neighbour_descriptor_interprocess(), false);
+
+ auto nr = (*itor).second;
+ assert( nr.num_neighbours() == 2 );
+
+ auto fst = nr.first();
+ auto snd = nr.second();
+
+ assert( (fst.iE == iE and fst.iT == iT) or
+ (snd.iE == iE and snd.iT == iT)
+ );
+
+ if ( fst.iE == iE and fst.iT == iT)
+ return std::make_pair(snd, true);
+
+ return std::make_pair(fst, true);
+ }
+
+ void clear(void)
+ {
+ face_cell_adj.clear();
+ }
+
+ void dump(void)
+ {
+ for (auto& [key, pair] : face_cell_adj)
+ {
+ std::cout << key << ": " << pair.first() << " - " << pair.second() << std::endl;
+ }
+ }
+
+ void mpi_bcast(int root, MPI_Comm comm)
+ {
+ priv_MPI_Bcast(face_cell_adj, root, comm);
+ }
+};
+
+#endif /* USE_MPI */
\ No newline at end of file
diff --git a/include/gmsh_io.h b/include/gmsh_io.h
index fd7f5c8..10a780f 100644
--- a/include/gmsh_io.h
+++ b/include/gmsh_io.h
@@ -111,6 +111,8 @@ class model
element_connectivity<element_key> conn;
#ifdef USE_MPI
+ interprocess_connectivity<element_key> ipconn;
+
/* Map from boundary tag to all the partitions it belongs to */
std::map<int, std::vector<int>> comm_map;
@@ -145,6 +147,7 @@ class model
void import_gmsh_entities_rankN(void);
int my_partition(void) const;
void make_partition_to_entities_map(void);
+ void map_interprocess_boundaries();
#endif /* USE_MPI */
void make_boundary_to_faces_map(void);
@@ -206,6 +209,12 @@ public:
entity& entity_world_lookup(int);
const entity& entity_world_lookup(int) const;
+
+ const interprocess_connectivity<element_key>& ip_connectivity() const
+ {
+ return ipconn;
+ }
+
#endif /* USE_MPI */
const element_connectivity<element_key>& connectivity() const
diff --git a/src/connectivity.cpp b/src/connectivity.cpp
index b42f125..a79245d 100644
--- a/src/connectivity.cpp
+++ b/src/connectivity.cpp
@@ -24,6 +24,27 @@ std::ostream& operator<<(std::ostream& os, const neighbour_descriptor& nd)
return os;
}
+
+
+
+neighbour_descriptor_interprocess::neighbour_descriptor_interprocess()
+{}
+
+bool
+neighbour_descriptor_interprocess::operator<(const neighbour_descriptor_interprocess& other) const
+{
+ std::array<size_t, 3> arr_me{iE, iT, via_iF};
+ std::array<size_t, 3> arr_other{other.iE, other.iT, other.via_iF};
+ return arr_me < arr_other;
+}
+
+std::ostream& operator<<(std::ostream& os, const neighbour_descriptor_interprocess& nd)
+{
+ os << "(" << nd.iE << ", " << nd.iT << ", " << nd.via_iF << ", PTAG: " << nd.parent_tag;
+ os << ", TAG: " << nd.tag << ", PART: " << nd.partition << ")";
+ return os;
+}
+
/***************************************************************************/
cell_descriptor::cell_descriptor()
{}
diff --git a/src/gmsh_io.cpp b/src/gmsh_io.cpp
index e403877..acc18cc 100644
--- a/src/gmsh_io.cpp
+++ b/src/gmsh_io.cpp
@@ -93,7 +93,9 @@ model::make_boundary_to_faces_map(void)
* used to determine who communicates with who. */
std::vector<int> parts;
gm::getPartitions(dim, tag, parts);
- comm_map[tag] = parts;
+ if (parts.size() == 2)
+ comm_map[tag] = parts;
+
if (pdim == dim)
surface_to_parent_map[tag] = ptag;
}
@@ -157,10 +159,6 @@ model::make_partition_to_entities_map()
void
model::update_connectivity(const entity& e, size_t entity_index)
{
-#ifdef USE_MPI
- //ASSERT_MPI_RANK_0;
-#endif /* USE_MPI */
-
for (size_t iT = 0; iT < e.num_cells(); iT++)
{
for (size_t iF = 0; iF < e.num_faces_per_elem(); iF++)
@@ -368,7 +366,7 @@ bool
model::is_interprocess_boundary(int tag)
{
if (num_partitions > 1)
- return comm_map.at(tag).size() == 2;
+ return comm_map.find(tag) != comm_map.end();
return false;
}
@@ -468,6 +466,73 @@ model::map_boundaries(void)
#endif
}
+#ifdef USE_MPI
+void
+model::map_interprocess_boundaries(void)
+{
+ /* Make a vector mapping element_key to entity tag */
+ using bfk_t = std::pair<element_key, int>;
+ std::vector<bfk_t> bfk;
+ for (auto& [tag, keys] : boundary_map)
+ {
+ if ( not is_interprocess_boundary(tag) )
+ continue;
+
+ for (auto& k : keys)
+ bfk.push_back( std::make_pair(k, tag) );
+ }
+ /* Sort it for lookups */
+ auto comp = [](const bfk_t& a, const bfk_t& b) -> bool {
+ return a.first < b.first;
+ };
+ std::sort(bfk.begin(), bfk.end(), comp);
+
+ auto lbcomp = [](const bfk_t& a, const element_key& b) -> bool {
+ return a.first < b;
+ };
+
+ auto connect = [&](const entity& e, int part) {
+ for (size_t iT = 0; iT < e.num_cells(); iT++)
+ {
+ auto nfpe = e.num_faces_per_elem();
+ for (size_t iF = 0; iF < nfpe; iF++)
+ {
+ auto fnum = iT*nfpe + iF;
+ element_key fk(e.face(fnum));
+
+ auto itor = std::lower_bound(bfk.begin(), bfk.end(), fk, lbcomp);
+ if ( (itor == bfk.end()) or not ((*itor).first == fk) )
+ continue;
+
+ neighbour_descriptor_interprocess ndi;
+ ndi.iE = e.number();
+ ndi.iT = iT;
+ ndi.via_iF = iF;
+ ndi.parent_tag = e.material_tag();
+ ndi.tag = e.gmsh_tag();
+ ndi.partition = part;
+ ndi.bar = e.cell(iT).barycenter();
+ ipconn.connect(fk, ndi);
+ }
+ }
+ };
+
+ for (auto& e : entities)
+ {
+ auto part = my_partition();
+ assert(part == 1);
+ connect(e, part);
+ }
+
+ for (auto& [part, res] : remote_entities)
+ {
+ assert(part != 1);
+ for (auto& e : res)
+ connect(e, part);
+ }
+}
+#endif
+
void
model::populate_from_gmsh_rank0()
{
@@ -512,6 +577,11 @@ model::populate_from_gmsh_rank0()
entities.clear();
import_gmsh_entities_rank0();
map_boundaries(); /* This must happen after import_gmsh_entities(). */
+
+#ifdef USE_MPI
+ map_interprocess_boundaries();
+ ipconn.mpi_bcast(0, MPI_COMM_WORLD);
+#endif /* USE_MPI */
}
#ifdef USE_MPI
@@ -539,6 +609,7 @@ model::populate_from_gmsh_rankN()
entities.clear();
import_gmsh_entities_rankN();
map_boundaries();
+ ipconn.mpi_bcast(0, MPI_COMM_WORLD);
}
#endif /* USE_MPI */
diff --git a/src/maxwell_common.cpp b/src/maxwell_common.cpp
index ba1e8db..e6a0feb 100644
--- a/src/maxwell_common.cpp
+++ b/src/maxwell_common.cpp
@@ -2,9 +2,10 @@
#include "maxwell_interface.h"
#include "maxwell_common.h"
-#define DIRICHLET 2.0
-#define ROBIN 1.0
-#define NEUMANN 0.0
+#define DIRICHLET 2.0
+#define ROBIN 1.0
+#define INTERPROCESS 1.0
+#define NEUMANN 0.0
namespace maxwell {
@@ -29,25 +30,32 @@ init_boundary_conditions(const model& mod, const parameter_loader& mpl,
continue;
double bc_coeff = DIRICHLET;
- switch ( mpl.boundary_type(bd.material_tag()) )
+ if (bd.type == face_type::INTERPROCESS_BOUNDARY)
{
- case boundary_condition::UNSPECIFIED:
- case boundary_condition::PEC:
- case boundary_condition::E_FIELD:
- bc_coeff = DIRICHLET;
- break;
-
- case boundary_condition::IMPEDANCE:
- case boundary_condition::PLANE_WAVE_E:
- case boundary_condition::PLANE_WAVE_H:
- bc_coeff = ROBIN;
- break;
-
- case boundary_condition::PMC:
- case boundary_condition::SURFACE_CURRENT:
- case boundary_condition::H_FIELD:
- bc_coeff = NEUMANN;
- break;
+ bc_coeff = INTERPROCESS;
+ }
+ else
+ {
+ switch ( mpl.boundary_type(bd.material_tag()) )
+ {
+ case boundary_condition::UNSPECIFIED:
+ case boundary_condition::PEC:
+ case boundary_condition::E_FIELD:
+ bc_coeff = DIRICHLET;
+ break;
+
+ case boundary_condition::IMPEDANCE:
+ case boundary_condition::PLANE_WAVE_E:
+ case boundary_condition::PLANE_WAVE_H:
+ bc_coeff = ROBIN;
+ break;
+
+ case boundary_condition::PMC:
+ case boundary_condition::SURFACE_CURRENT:
+ case boundary_condition::H_FIELD:
+ bc_coeff = NEUMANN;
+ break;
+ }
}
auto& pf = e.face(iF);
@@ -73,6 +81,9 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl,
pbH = vecxd::Zero( mod.num_fluxes() );
auto& conn = mod.connectivity();
+#ifdef USE_MPI
+ auto& ipconn = mod.ip_connectivity();
+#endif /* USE_MPI */
for (auto& e : mod)
{
for (size_t iT = 0; iT < e.num_cells(); iT++)
@@ -120,16 +131,39 @@ init_lax_milgram(const model& mod, const parameter_loader& mpl,
paH(glob_ofs+iD) = aH;
pbH(glob_ofs+iD) = bH;
}
+ continue;
}
- else
+#ifdef USE_MPI
+ auto [ip_neigh, has_ip_neigh] = ipconn.neighbour_via(e.number(), iT, fk);
+ if (has_ip_neigh)
{
+ auto ne_eps = mpl.epsilon(ip_neigh.parent_tag, ip_neigh.bar);
+ auto ne_mu = mpl.mu(ip_neigh.parent_tag, ip_neigh.bar);
+
+ auto ne_Z = std::sqrt(ne_mu/ne_eps);
+ auto ne_Y = 1./ne_Z;
+
+ auto aE = ne_Z/(my_Z + ne_Z);
+ auto bE = 1./(my_Z + ne_Z);
+ auto aH = ne_Y/(my_Y + ne_Y);
+ auto bH = 1./(my_Y + ne_Y);
+
for (size_t iD = 0; iD < rf.num_basis_functions(); iD++)
{
- paE(glob_ofs+iD) = 0.5;
- pbE(glob_ofs+iD) = 0.5/my_Z;
- paH(glob_ofs+iD) = 0.5;
- pbH(glob_ofs+iD) = 0.5/my_Y;
+ paE(glob_ofs+iD) = aE;
+ pbE(glob_ofs+iD) = bE;
+ paH(glob_ofs+iD) = aH;
+ pbH(glob_ofs+iD) = bH;
}
+ continue;
+ }
+#endif /* USE_MPI */
+ for (size_t iD = 0; iD < rf.num_basis_functions(); iD++)
+ {
+ paE(glob_ofs+iD) = 0.5;
+ pbE(glob_ofs+iD) = 0.5/my_Z;
+ paH(glob_ofs+iD) = 0.5;
+ pbH(glob_ofs+iD) = 0.5/my_Y;
}
}
}
--
GitLab