Add and calculate coordinate container for off-rank neighbor nodes. (#…

…1081)

* Add and calculate coordinate container for off-rank neighbor nodes.

+ Create node-coordinate map that is 1-to-1 with map to node indexes in dual ghost layout.
+ Add specializations on double of allgatherv and determinate_allgatherv.
+ Use double determinate_allgatherv to get ghost neighbor node coordinates.
+ Use consistent loop variable names in loops serializing data for allgatherv.
+ Add tests of ghost coordinate values to tstDraco_Mesh_DD unit test.
+ Add mesh decomposition ASCII art schematics to tstDraco_Mesh_DD tests.

Note: this performs an additional allgatherv of four doubles between ranks
(two doubles for each coordinate adjacent to the node and its neighbor cell).

* Use soft_equiv for checking equality of doubles.

* Incorporate recommended improvements from KT.

+ Replace push_back with emplace_back in dual ghost map generation.
+ Add unit tests of [in]determinate_allgatherv specializations.
+ Fix else branch (missing) in logic for other gatherv unit tests.

* Convert type to auto in unit test, to address clang-tidy result.

src/c4/C4_MPI_gather_scatter_pt.cc

@@ -46,6 +46,9 @@ template int gatherv<char>(char *send_buffer, int send_size, char *receive_buffe
 template int allgatherv<unsigned>(unsigned *send_buffer, int send_size, unsigned *receive_buffer,
                                   int *receive_sizes, int *receive_displs);
 
+template int allgatherv<double>(double *send_buffer, int send_size, double *receive_buffer,
+                                int *receive_sizes, int *receive_displs);
+
 //------------------------------------------------------------------------------------------------//
 template int scatter<unsigned>(unsigned *send_buffer, unsigned *receive_buffer, int size);
 

src/c4/gatherv_pt.cc

@@ -51,6 +51,8 @@ template void indeterminate_allgatherv<unsigned>(std::vector<unsigned> &outgoing
 
 template void determinate_allgatherv<unsigned>(std::vector<unsigned> &outgoing_data,
                                                std::vector<std::vector<unsigned>> &incoming_data);
+template void determinate_allgatherv<double>(std::vector<double> &outgoing_data,
+                                             std::vector<std::vector<double>> &incoming_data);
 
 } // end namespace rtt_c4
 

src/c4/test/tstGatherScatter.cc

@@ -99,6 +99,7 @@ void tstIndeterminateGatherScatterv(UnitTest &ut) {
       for (unsigned p = 0; p < number_of_processors; ++p) {
         if (receive[p].size() != p) {
           FAILMSG("NOT correct number of elements in gatherv");
+        } else {
           for (unsigned i = 0; i < p; ++i) {
             if (receive[p][i] != p)
               FAILMSG("NOT correct values in gatherv");
@@ -137,6 +138,7 @@ void tstIndeterminateGatherScatterv(UnitTest &ut) {
       for (unsigned p = 0; p < number_of_processors; ++p) {
         if (receive[p].size() != p) {
           FAILMSG("NOT correct number of elements in gatherv");
+        } else {
           for (unsigned i = 0; i < p; ++i) {
             if (!rtt_dsxx::soft_equiv(receive[p][i], static_cast<double>(p)))
               FAILMSG("NOT correct values in gatherv");
@@ -175,6 +177,7 @@ void tstIndeterminateGatherScatterv(UnitTest &ut) {
       for (unsigned p = 0; p < number_of_processors; ++p) {
         if (receive[p].size() != p) {
           FAILMSG("NOT correct number of elements in gatherv");
+        } else {
           for (unsigned i = 0; i < p; ++i) {
             if (receive[p][i] != static_cast<int>(p))
               FAILMSG("NOT correct values in gatherv");
@@ -206,21 +209,15 @@ void tstIndeterminateGatherScatterv(UnitTest &ut) {
     indeterminate_gatherv(emptysend, emptyreceive);
     PASSMSG("No exception thrown for indeterminate_gatherv with empty containers.");
 
-    if (emptysend.size() != 0)
-      ITFAILS;
-    if (emptyreceive.size() != number_of_processors)
-      ITFAILS;
-    if (emptyreceive[pid].size() != 0)
-      ITFAILS;
+    FAIL_IF(emptysend.size() != 0);
+    FAIL_IF(emptyreceive.size() != number_of_processors);
+    FAIL_IF(emptyreceive[pid].size() != 0);
 
     indeterminate_scatterv(emptyreceive, emptysend);
 
-    if (emptysend.size() != 0)
-      ITFAILS;
-    if (emptyreceive.size() != number_of_processors)
-      ITFAILS;
-    if (emptyreceive[pid].size() != 0)
-      ITFAILS;
+    FAIL_IF(emptysend.size() != 0);
+    FAIL_IF(emptyreceive.size() != number_of_processors);
+    FAIL_IF(emptyreceive[pid].size() != 0);
   }
 
   return;
@@ -249,6 +246,7 @@ void tstDeterminateGatherScatterv(UnitTest &ut) {
       for (unsigned p = 0; p < number_of_processors; ++p) {
         if (receive[p].size() != p) {
           FAILMSG("NOT correct number of elements in gatherv");
+        } else {
           for (unsigned i = 0; i < p; ++i) {
             if (receive[p][i] != p)
               FAILMSG("NOT correct values in gatherv");
@@ -292,6 +290,7 @@ void tstDeterminateGatherScatterv(UnitTest &ut) {
       for (unsigned p = 0; p < number_of_processors; ++p) {
         if (receive[p].size() != p) {
           FAILMSG("NOT correct number of elements in gatherv");
+        } else {
           for (unsigned i = 0; i < p; ++i) {
             if (!rtt_dsxx::soft_equiv(receive[p][i], static_cast<double>(p)))
               FAILMSG("NOT correct values in gatherv");
@@ -335,6 +334,7 @@ void tstDeterminateGatherScatterv(UnitTest &ut) {
       for (unsigned p = 0; p < number_of_processors; ++p) {
         if (receive[p].size() != p) {
           FAILMSG("NOT correct number of elements in gatherv");
+        } else {
           for (unsigned i = 0; i < p; ++i) {
             if (receive[p][i] != static_cast<int>(p))
               FAILMSG("NOT correct values in gatherv");
@@ -378,6 +378,7 @@ void tstDeterminateGatherScatterv(UnitTest &ut) {
       for (unsigned p = 0; p < number_of_processors; ++p) {
         if (receive[p].size() != p) {
           FAILMSG("NOT correct number of elements in gatherv");
+        } else {
           for (unsigned i = 0; i < p; ++i) {
             if (receive[p][i] != 'A')
               FAILMSG("NOT correct values in gatherv");
@@ -414,8 +415,7 @@ void topology_report(UnitTest &ut) {
   // Look at the data found on the IO proc.
   if (my_mpi_rank == 0) {
 
-    if (procnames[my_mpi_rank].size() != namelen)
-      ITFAILS;
+    FAIL_IF(procnames[my_mpi_rank].size() != namelen);
 
     // Count unique processors
     vector<string> unique_processor_names;
@@ -434,8 +434,7 @@ void topology_report(UnitTest &ut) {
 
     for (size_t i = 0; i < mpi_ranks; ++i) {
       std::cout << "\n  - MPI rank " << i << " is on " << procnames[i];
-      if (procnames[i].size() < 1)
-        ITFAILS;
+      FAIL_IF(procnames[i].size() < 1);
     }
     std::cout << std::endl;
 
@@ -464,6 +463,97 @@ void topology_report(UnitTest &ut) {
   return;
 }
 
+//------------------------------------------------------------------------------------------------//
+void tstDeterminateAllGatherv(UnitTest &ut) {
+  unsigned const pid = node();
+  unsigned const number_of_processors = nodes();
+
+  { // T=unsigned
+    vector<unsigned> send(pid, pid);
+
+    // determinate_allgatherv already has the data sizes from the other MPI/C4 ranks/nodes
+    vector<vector<unsigned>> receive(number_of_processors);
+    for (unsigned p = 0; p < number_of_processors; ++p)
+      receive[p].resize(p);
+
+    // gather data from all nodes to all nodes at rank index in receive vector
+    determinate_allgatherv(send, receive);
+    PASSMSG("No exception thrown");
+
+    // check values gathered from each rank (and check on each rank of course)
+    for (unsigned p = 0; p < number_of_processors; ++p) {
+      for (unsigned i = 0; i < p; ++i) {
+        if (receive[p][i] != p)
+          FAILMSG("NOT correct values in allgatherv");
+      }
+    }
+  }
+
+  { // T=double
+    vector<double> send(pid, static_cast<double>(pid));
+
+    // determinate_allgatherv already has the data sizes from the other MPI/C4 ranks/nodes
+    vector<vector<double>> receive(number_of_processors);
+    for (unsigned p = 0; p < number_of_processors; ++p)
+      receive[p].resize(p);
+
+    // gather data from all nodes to all nodes at rank index in receive vector
+    determinate_allgatherv(send, receive);
+    PASSMSG("No exception thrown");
+
+    // check values gathered from each rank (and check on each rank of course)
+    for (unsigned p = 0; p < number_of_processors; ++p) {
+      const auto p_dbl = static_cast<double>(p);
+      for (unsigned i = 0; i < p; ++i) {
+        if (!rtt_dsxx::soft_equiv(receive[p][i], p_dbl))
+          FAILMSG("NOT correct values in allgatherv");
+      }
+    }
+  }
+
+  // successful test output
+  if (ut.numFails == 0)
+    PASSMSG("tstDeterminateAllGatherv tests ok.");
+  return;
+}
+
+//------------------------------------------------------------------------------------------------//
+void tstIndeterminateAllGatherv(UnitTest &ut) {
+  unsigned const pid = node();
+  unsigned const number_of_processors = nodes();
+
+  // T=unsigned
+  vector<unsigned> send(pid, pid);
+  vector<vector<unsigned>> receive;
+
+  // gather data from all nodes to all nodes at rank index in receive vector
+  indeterminate_allgatherv(send, receive);
+  PASSMSG("No exception thrown");
+
+  // check the size of the receiving vector is now the number of MPI/C4 ranks/nodes
+  if (receive.size() == number_of_processors)
+    PASSMSG("correct number of processors in allgatherv");
+  else
+    FAILMSG("NOT correct number of processors in allgatherv");
+
+  // check values gathered from each rank (and check on each rank of course)
+  for (unsigned p = 0; p < number_of_processors; ++p) {
+    if (receive[p].size() != p) {
+      FAILMSG("NOT correct number of elements in allgatherv");
+    } else {
+      for (unsigned i = 0; i < p; ++i) {
+        if (receive[p][i] != p)
+          FAILMSG("NOT correct values in allgatherv");
+      }
+    }
+  }
+
+  // successful test output
+  if (ut.numFails == 0)
+    PASSMSG("tstIndeterminateAllGatherv tests ok.");
+  return;
+}
+
 //------------------------------------------------------------------------------------------------//
 int main(int argc, char *argv[]) {
   rtt_c4::ParallelUnitTest ut(argc, argv, release);
@@ -473,6 +563,8 @@ int main(int argc, char *argv[]) {
     tstIndeterminateGatherScatterv(ut);
     tstDeterminateGatherScatterv(ut);
     topology_report(ut);
+    tstDeterminateAllGatherv(ut);
+    tstIndeterminateAllGatherv(ut);
   }
   UT_EPILOG(ut);
 }

src/mesh/Draco_Mesh.cc

@@ -478,6 +478,8 @@ void Draco_Mesh::compute_node_to_cell_linkage(
     return;
 
   //----------------------------------------------------------------------------------------------//
+  // \todo: Can this ghost data all gather procedure be split into sensible member functions?
+  //
   // When domain-decomposed, the following creates a map of local nodes to all ghost cells.
   // The procedure makes use of existing ghost data across cell faces, as follows:
   //
@@ -543,22 +545,35 @@ void Draco_Mesh::compute_node_to_cell_linkage(
   const size_t num_serial = cellnodes_per_serial.size() / 3;
 
   //----------------------------------------------------------------------------------------------//
-  // initialize a serial array of global node indices
+  // initialize a serial array of global node indices and neighbor node coordinates
   std::vector<unsigned> global_node_per_serial(num_serial);
+  std::vector<double> coord_nbrs_per_serial(4 * num_serial);
 
   // map global ghost node indices to serial index
   size_t serial_count = 0;
   for (const auto &global_node_cellnode_pair : global_node_to_local_cellnodes) {
 
     // get the number of local cells for this node
-    const size_t num_node_cells = global_node_cellnode_pair.second.size();
+    const size_t num_cell_nbrs = global_node_cellnode_pair.second.size();
+
+    // set the global and neighbor node coordinates per serial index
+    for (size_t j = 0; j < num_cell_nbrs; ++j) {
 
-    // set to the global node at the serial index for this local cell
-    for (size_t node_cell = 0; node_cell < num_node_cells; ++node_cell)
-      global_node_per_serial[serial_count + node_cell] = global_node_cellnode_pair.first;
+      // set to the global node at the serial index for this local cell
+      global_node_per_serial[serial_count + j] = global_node_cellnode_pair.first;
+
+      // set the neighbor node coordinates
+      const unsigned node1 = global_node_cellnode_pair.second[j].second[0];
+      const unsigned node2 = global_node_cellnode_pair.second[j].second[1];
+      const size_t cnbrs_offset = 4 * (serial_count + j);
+      coord_nbrs_per_serial[cnbrs_offset] = node_coord_vec[node1][0];
+      coord_nbrs_per_serial[cnbrs_offset + 1] = node_coord_vec[node1][1];
+      coord_nbrs_per_serial[cnbrs_offset + 2] = node_coord_vec[node2][0];
+      coord_nbrs_per_serial[cnbrs_offset + 3] = node_coord_vec[node2][1];
+    }
 
     // increment count over serial index
-    serial_count += num_node_cells;
+    serial_count += num_cell_nbrs;
   }
 
   // check that serial vector has been filled
@@ -586,10 +601,22 @@ void Draco_Mesh::compute_node_to_cell_linkage(
   // gather the global ghost node indices per serial index per rank
   rtt_c4::determinate_allgatherv(cellnodes_per_serial, cellnodes_per_serial_per_rank);
 
+  //----------------------------------------------------------------------------------------------//
+  // gather the local coordinate values for neighbor nodes per serial per rank
+
+  // resize gather target for neighbor node coordinates
+  std::vector<std::vector<double>> coord_nbrs_per_serial_per_rank(num_ranks);
+  for (unsigned rank = 0; rank < num_ranks; ++rank)
+    coord_nbrs_per_serial_per_rank[rank].resize(4 * global_node_per_serial_per_rank[rank].size());
+
+  // gather the global ghost node indices per serial index per rank
+  rtt_c4::determinate_allgatherv(coord_nbrs_per_serial, coord_nbrs_per_serial_per_rank);
+
   //----------------------------------------------------------------------------------------------//
   // merge global_node_per_serial_per_rank and cells_per_serial_per_rank into map per rank
 
   std::vector<std::map<unsigned, std::vector<CellNodes_Pair>>> ghost_dualmap_per_rank(num_ranks);
+  std::vector<std::map<unsigned, std::vector<Coord_NBRS>>> ghost_coord_nbrs_per_rank(num_ranks);
   for (unsigned rank = 0; rank < num_ranks; ++rank) {
 
     // short-cut to serialized vector size from rank
@@ -605,8 +632,15 @@ void Draco_Mesh::compute_node_to_cell_linkage(
       const std::array<unsigned, 2> node_nbrs = {cellnodes_per_serial_per_rank[rank][3 * i + 1],
                                                  cellnodes_per_serial_per_rank[rank][3 * i + 2]};
 
-      // accumulate local cells (for rank) adjacent to this global node
-      ghost_dualmap_per_rank[rank][global_node].push_back(std::make_pair(local_cell, node_nbrs));
+      // accumulate local cells and neighbor nodes (for rank) adjacent to this global node
+      ghost_dualmap_per_rank[rank][global_node].emplace_back(std::make_pair(local_cell, node_nbrs));
+
+      // accumulate neighbor node coordinates (in same order as node indices about global_node)
+      const std::array<double, 2> crd_nbr1 = {coord_nbrs_per_serial_per_rank[rank][4 * i],
+                                              coord_nbrs_per_serial_per_rank[rank][4 * i + 1]};
+      const std::array<double, 2> crd_nbr2 = {coord_nbrs_per_serial_per_rank[rank][4 * i + 2],
+                                              coord_nbrs_per_serial_per_rank[rank][4 * i + 3]};
+      ghost_coord_nbrs_per_rank[rank][global_node].emplace_back(std::make_pair(crd_nbr1, crd_nbr2));
     }
   }
 
@@ -616,10 +650,13 @@ void Draco_Mesh::compute_node_to_cell_linkage(
   for (unsigned node = 0; node < num_nodes; ++node)
     global_to_local_node[global_node_number[node]] = node;
 
+  //----------------------------------------------------------------------------------------------//
+  // generate dual ghost layout and coordinates by setting each ranks nodes back to local values
+
   // get this (my) rank
   const unsigned my_rank = rtt_c4::node();
 
-  // generate dual gost layout
+  // generate dual ghost layout
   for (unsigned rank = 0; rank < num_ranks; ++rank) {
 
     // exclude this rank
@@ -653,12 +690,17 @@ void Draco_Mesh::compute_node_to_cell_linkage(
 
       // append each local-cell-rank pair to dual ghost layout
       for (auto local_cellnodes : ghost_dualmap_per_rank[rank].at(gl_node))
-        node_to_ghost_cell_linkage[node].push_back(std::make_pair(local_cellnodes, rank));
+        node_to_ghost_cell_linkage[node].emplace_back(std::make_pair(local_cellnodes, rank));
+
+      // append each ghost coordinate pair bounding a ghost cell neighboring this node
+      for (auto coord_nbrs : ghost_coord_nbrs_per_rank[rank].at(gl_node))
+        node_to_ghost_coord_linkage[node].emplace_back(coord_nbrs);
     }
   }
 
   // since this mesh was constructed with ghost data, the resulting map must have non-zero size
   Ensure(node_to_ghost_cell_linkage.size() > 0);
+  Ensure(node_to_ghost_coord_linkage.size() > 0);
 }
 
 } // end namespace rtt_mesh

src/mesh/Draco_Mesh.hh

@@ -63,8 +63,10 @@ public:
 
   // e.g.: (key: node, value: vector of pairs of rank and local cell index on the rank)
   // vectors of pairs are ordered in increasing rank, by construction (see Draco_Mesh.cc)
+  using Coord_NBRS = std::pair<std::array<double, 2>, std::array<double, 2>>;
   using Dual_Ghost_Layout =
       std::map<unsigned int, std::vector<std::pair<CellNodes_Pair, unsigned int>>>;
+  using Dual_Ghost_Layout_Coords = std::map<unsigned int, std::vector<Coord_NBRS>>;
 
 protected:
   // >>> DATA
@@ -118,6 +120,9 @@ protected:
   // Node map to vector of ghost cells
   Dual_Ghost_Layout node_to_ghost_cell_linkage;
 
+  // Node map to vector of adjacent coordinates bounding adjacent ghost cells
+  Dual_Ghost_Layout_Coords node_to_ghost_coord_linkage;
+
 public:
   //! Constructor.
   Draco_Mesh(unsigned dimension_, Geometry geometry_,
@@ -158,6 +163,7 @@ public:
   Layout get_cg_linkage() const { return cell_to_ghost_cell_linkage; }
   Dual_Layout get_nc_linkage() const { return node_to_cellnode_linkage; }
   Dual_Ghost_Layout get_ngc_linkage() const { return node_to_ghost_cell_linkage; }
+  Dual_Ghost_Layout_Coords get_ngcoord_linkage() const { return node_to_ghost_coord_linkage; }
 
   // >>> SERVICES