FEAT: get grid to grid message passing to work

Author: aaronjridley

src/grid_match.cpp

@@ -62,11 +62,27 @@ bool exchange_information(int64_t *nPointsToPass,
   return true;
 }
 
-bool grid_match(Grid gGrid, 
-                Grid mGrid, 
+// -----------------------------------------------------------------------------
+// This function:
+//   on the requesting information side:
+//     - figures out which processor each point of the other grid is on
+//     - counts the points for each processor
+//     - exchanges how many points to pass for each processor
+//     - makes lists of coordinates to send to each processor
+//     - sends those lists
+//   on the interpolator side:
+//     - builds interpolators for the requested information
+// -----------------------------------------------------------------------------
+
+bool grid_match(Grid &gGrid, 
+                Grid &mGrid, 
                 Quadtree gQuadtree,
                 Quadtree mQuadtree) {
 
+  std::string function = "grid_match";
+  static int iFunction = -1;
+  report.enter(function, iFunction);
+
   // Let's do magnetic to geographic first:
 
   int64_t iX, mnX = mGrid.get_nX();
@@ -86,6 +102,7 @@ bool grid_match(Grid gGrid,
 
   // This is not the most efficient way to do this, but the first pass, let's
   // just count how many points we need to send to the other processors:
+  mGrid.gridToGridMap.set_size(mnX, mnY, mnZ);
   for (iX = mGCs; iX < mnX - mGCs; iX++) {
     for (iY = mGCs; iY < mnY - mGCs; iY++) {
       for (iZ = mGCs; iZ < mnZ - mGCs; iZ++) {
@@ -103,6 +120,7 @@ bool grid_match(Grid gGrid,
         if (jNode < 0 || jNode >= nGrids) {
             std::cout << "out of bounds!!! " << jNode << "\n";
         }
+        mGrid.gridToGridMap(iX, iY, iZ) = jNode;
         nPointsToPass[jNode] = nPointsToPass[jNode]+1;
         /* std::cout << "lon, lat, node: " << lon*cRtoD << " "
             << lat*cRtoD << " "
@@ -115,13 +133,13 @@ bool grid_match(Grid gGrid,
       }
     }
   }
-  std::cout << "made it here: " << iProc << "\n";
   MPI_Barrier(aether_comm);
 
-  for (jNode = 0; jNode < nGrids ; jNode++)
-    std::cout << "nPtsToPass : " << iProc << " " << nPointsToPass[jNode] << "\n";
-
-  std::cout << "sending number of points :\n";
+  if (report.test_verbose(3)) {
+    for (jNode = 0; jNode < nGrids ; jNode++)
+      std::cout << "nPtsToPass : " << iProc << " " << nPointsToPass[jNode] << "\n";
+    std::cout << "sending number of points :\n";
+  }
 
   // This section sends the number of points that need to be transfered to each processor.
   // Then the processor saves the number of points, so it can be remembered, and both the
@@ -135,27 +153,23 @@ bool grid_match(Grid gGrid,
     nPointsToReceive[jNode] = nPointsDummy[iGrid];
   }
 
-  MPI_Barrier(aether_comm);
-
-  for (jNode = 0; jNode < nGrids ; jNode++) {
-    std::cout << "nPtsToReceive : " << iProc << " " << jNode << " " << nPointsToReceive[jNode] << "\n";
-    MPI_Barrier(aether_comm);
+  if (report.test_verbose(3)) {
+    for (jNode = 0; jNode < nGrids ; jNode++) {
+      std::cout << "nPtsToReceive : " << iProc << " " << jNode << " " << nPointsToReceive[jNode] << "\n";
+    }
   }
 
   //  Now we need to create an array of send points and an array of receive points.
   std::vector<precision_t *> latsToPass(nGrids);
   std::vector<precision_t *> lonsToPass(nGrids);
   std::vector<precision_t *> altsToPass(nGrids);
-  for (jNode = 0; jNode < nGrids ; jNode++) {
-    latsToPass[jNode] = static_cast<precision_t*>(malloc(nPointsToPass[jNode] * sizeof(precision_t)));
-    lonsToPass[jNode] = static_cast<precision_t*>(malloc(nPointsToPass[jNode] * sizeof(precision_t)));
-    altsToPass[jNode] = static_cast<precision_t*>(malloc(nPointsToPass[jNode] * sizeof(precision_t)));
-  }
-
   std::vector<precision_t *> latsToInterTo(nGrids);
   std::vector<precision_t *> lonsToInterTo(nGrids);
   std::vector<precision_t *> altsToInterTo(nGrids);
   for (jNode = 0; jNode < nGrids ; jNode++) {
+    latsToPass[jNode] = static_cast<precision_t*>(malloc(nPointsToPass[jNode] * sizeof(precision_t)));
+    lonsToPass[jNode] = static_cast<precision_t*>(malloc(nPointsToPass[jNode] * sizeof(precision_t)));
+    altsToPass[jNode] = static_cast<precision_t*>(malloc(nPointsToPass[jNode] * sizeof(precision_t)));
     latsToInterTo[jNode] = static_cast<precision_t*>(malloc(nPointsToReceive[jNode] * sizeof(precision_t)));
     lonsToInterTo[jNode] = static_cast<precision_t*>(malloc(nPointsToReceive[jNode] * sizeof(precision_t)));
     altsToInterTo[jNode] = static_cast<precision_t*>(malloc(nPointsToReceive[jNode] * sizeof(precision_t)));
@@ -186,11 +200,123 @@ bool grid_match(Grid gGrid,
     }
   }
   bool didWork;
+  // Pass first coordinate (lons)
+  didWork = exchange_information(nPointsToPass,
+                                 lonsToPass,
+                                 nPointsToReceive,
+                                 lonsToInterTo);
+  // Pass second coordinate (lats)
   didWork = exchange_information(nPointsToPass,
                                  latsToPass,
                                  nPointsToReceive,
                                  latsToInterTo);
+  // Pass third coordinate (alts):
+  didWork = exchange_information(nPointsToPass,
+                                 altsToPass,
+                                 nPointsToReceive,
+                                 altsToInterTo);
 
+  if (report.test_verbose(2)) {
+    for (jNode = 0; jNode < nGrids ; jNode++) {
+      std::cout << "Received the following points from iGrid = " << jNode << "\n";
+      std::cout << " -> points received : " << nPointsToReceive[jNode] << "\n";
+      for (int64_t iPt = 0; iPt < nPointsToReceive[jNode]; iPt++) 
+        std::cout << "  -> " << iPt << " "
+                    << lonsToInterTo[jNode][iPt] << " "
+                    << latsToInterTo[jNode][iPt] << " "
+                    << altsToInterTo[jNode][iPt] << "\n";
+    }
+  }
 
-  return true;
+  struct grid_to_grid_t oneGrid;
+
+  int64_t nPts;
+  for (jNode = 0; jNode < nGrids ; jNode++) {
+    // These are backwards now, since we will switch sender and reciever:
+    oneGrid.nPts = nPointsToReceive[jNode];
+    oneGrid.nPtsReceive = nPointsToPass[jNode];
+    oneGrid.iProcTo = iMember * nGrids + jNode;
+    if (report.test_verbose(2))
+      std::cout << "Making interpolation coefficients for : " << jNode 
+        << "; points : " << oneGrid.nPts << "\n";
+    if (oneGrid.nPts > 0) {
+      // Interpolation function takes vectors, 
+      // so transfer these arrays to vectors:
+      std::vector<precision_t> Lons(oneGrid.nPts);
+      std::vector<precision_t> Lats(oneGrid.nPts);
+      std::vector<precision_t> Alts(oneGrid.nPts);
+      for (int64_t iPt = 0; iPt < oneGrid.nPts; iPt++) {
+        Lons[iPt] = lonsToInterTo[jNode][iPt];
+        Lats[iPt] = latsToInterTo[jNode][iPt];
+        Alts[iPt] = altsToInterTo[jNode][iPt];
+      }
+      oneGrid.interpCoefs = gGrid.get_interpolation_coefs(Lons, Lats, Alts);
+    }
+    gGrid.gridToGridCoefs.push_back(oneGrid);
+  }
+
+  report.exit(function);
+  return didWork;
 }
+
+bool get_data_from_other_grid(Grid &gGrid, 
+                              Grid &mGrid,
+                              arma_cube &gData,
+                              arma_cube &mData) {
+
+  std::string function = "get_data_from_other_grid";
+  static int iFunction = -1;
+  report.enter(function, iFunction);
+
+  int64_t jNode, iPt;
+  std::vector<precision_t *> dataToSend(nGrids);
+  std::vector<precision_t *> dataToReceive(nGrids);
+  int64_t *nPointsToSend = static_cast<int64_t*>(malloc(nGrids * sizeof(int64_t)));
+  int64_t *nPointsToReceive = static_cast<int64_t*>(malloc(nGrids * sizeof(int64_t)));
+
+  for (jNode = 0; jNode < nGrids ; jNode++) {
+    if (report.test_verbose(2))
+      std::cout << "nPts : " << jNode << " " << gGrid.gridToGridCoefs[jNode].nPts << "\n";
+    nPointsToSend[jNode] = gGrid.gridToGridCoefs[jNode].nPts;
+    nPointsToReceive[jNode] = gGrid.gridToGridCoefs[jNode].nPtsReceive;
+    dataToSend[jNode] = static_cast<precision_t*>(malloc(gGrid.gridToGridCoefs[jNode].nPts * sizeof(precision_t)));
+    dataToReceive[jNode] = static_cast<precision_t*>(malloc(gGrid.gridToGridCoefs[jNode].nPtsReceive * sizeof(precision_t)));
+    std::vector<precision_t> values = gGrid.get_interpolation_values(gData, gGrid.gridToGridCoefs[jNode].interpCoefs);
+    
+    for (iPt = 0; iPt < gGrid.gridToGridCoefs[jNode].nPts; iPt++) {
+      dataToSend[jNode][iPt] = values[iPt];
+      if (report.test_verbose(2))
+        std::cout << "datatosend : " << iPt << " " << dataToSend[jNode][iPt] << "\n";
+    }
+  }
+  bool didWork = exchange_information(nPointsToSend,
+                                      dataToSend,
+                                      nPointsToReceive,
+                                      dataToReceive);
+  int64_t iX, mnX = mGrid.get_nX();
+  int64_t iY, mnY = mGrid.get_nY();
+  int64_t iZ, mnZ = mGrid.get_nZ();
+  int64_t mGCs = mGrid.get_nGCs();
+  std::vector<int64_t> iCounter(nGrids);
+  for (jNode = 0; jNode < nGrids ; jNode++)
+    iCounter[jNode] = 0;
+
+  for (iX = mGCs; iX < mnX - mGCs; iX++) {
+    for (iY = mGCs; iY < mnY - mGCs; iY++) {
+      for (iZ = mGCs; iZ < mnZ - mGCs; iZ++) {
+        jNode = mGrid.gridToGridMap(iX, iY, iZ);
+        if (report.test_verbose(2)) {
+          std::cout << "unpacking point : " << iX << " " << iY << " " << iZ << " " << jNode << " " 
+          << iCounter[jNode] << " " << dataToReceive[jNode][iCounter[jNode]] << "\n";
+        }
+
+        mData(iX, iY, iZ) = dataToReceive[jNode][iCounter[jNode]];
+        iCounter[jNode] = iCounter[jNode]+1;
+      }
+    }
+  }
+
+  report.exit(function);
+  return true;
+
+}