enable RCBPtn to return multiple candidate ranks

CMakeLists.txt

@@ -158,6 +158,17 @@ if(BUILD_TESTING)
     NAME2 app PROCS2 1 EXE2 ./test_setup_classPtn ARGS2 0)
   add_exe(test_query test_query.cpp)
   mpi_test(test_query_1p 1 ./test_query)
+  # Unit tests for RCBPtn::GetCandidateRanks.
+  add_exe(test_candidateRanks test_candidateRanks.cpp)
+  mpi_test(test_candidateRanks_1p 1 ./test_candidateRanks)
+  # A standalone test of the boundary routing logic.  Routing uses the
+  # multi-candidate API (RCBPtn::GetCandidateRanks), which resolves boundary
+  # ambiguity within the rendezvous partition; this test still FAILS on rows
+  # whose true mesh owner is in a different rendezvous region (misalignment
+  # beyond adjacency).  WILL_FAIL keeps CI green while the known gap exists.
+  add_exe(test_boundaryRouting test_boundaryRouting.cpp)
+  mpi_test(test_boundaryRouting_1p 1 ./test_boundaryRouting)
+  set_tests_properties(test_boundaryRouting_1p PROPERTIES WILL_FAIL TRUE)
   add_exe(test_send test_send.cpp)
   dual_mpi_test(TESTNAME test_send_3p TIMEOUT ${test_timeout}
     NAME1 rdv PROCS1 4 EXE1 ./test_send ARGS1 1

redev.cpp

@@ -9,6 +9,8 @@
 #include <chrono>         // std::chrono::seconds
 #include <string>         // std::stoi
 #include <algorithm>      // std::find_if
+#include <cmath>          // std::abs, std::log2, std::pow
+#include <limits>         // std::numeric_limits
 
 namespace {
   //Wait for the file to be created by the writer.
@@ -236,6 +238,98 @@ namespace redev {
     return cuts;
   }
 
+  std::array<std::vector<redev::Real>, 3> RCBPtn::BuildDimCuts() const {
+    const auto len = static_cast<redev::LO>(cuts.size());
+    std::array<std::vector<redev::Real>, 3> dc;
+    for(redev::LO i=1; i<len; i++) {
+      const auto level = static_cast<redev::LO>(std::log2(i));
+      dc[level % dim].push_back(cuts[i]);
+    }
+    for(auto& list : dc) {
+      std::sort(list.begin(), list.end());
+      list.erase(std::unique(list.begin(), list.end()), list.end());
+    }
+    return dc;
+  }
+
+  redev::Real RCBPtn::LocalSpacing(
+      const std::array<std::vector<redev::Real>, 3>& dimCuts,
+      redev::LO d, redev::Real v, redev::Real fallback) {
+    const auto& list = dimCuts[d];
+    const auto it = std::lower_bound(list.begin(), list.end(), v);
+    if(it == list.end()) return fallback;
+    if(it != list.begin() && it + 1 != list.end())
+      return std::min(v - *(it-1), *(it+1) - v);
+    if(it != list.begin()) return v - *(it-1);
+    if(it + 1 != list.end()) return *(it+1) - v;
+    return fallback;
+  }
+
+  redev::Real RCBPtn::NodeTolerance(
+      const std::array<std::vector<redev::Real>, 3>& dimCuts,
+      redev::LO d, redev::Real v, redev::Real ptCoord,
+      redev::Real tolFactor) {
+    // hLocal = 0 when dim d has fewer than 2 distinct cuts: there are no
+    // adjacent sub-domains to broaden to on this axis, so the tolerance
+    // collapses to the machine-epsilon floor.
+    const auto hLocal = LocalSpacing(dimCuts, d, v, /*fallback=*/0);
+    const auto epsGuard =
+        defaultEpsFactor * std::numeric_limits<redev::Real>::epsilon() *
+        std::max({std::abs(ptCoord), std::abs(v), redev::Real(1)});
+    return std::max(tolFactor * hLocal, epsGuard);
+  }
+
+  redev::LOs RCBPtn::TraverseCandidates(
+      std::array<redev::Real, 3>& pt, redev::Real tolFactor,
+      const std::array<std::vector<redev::Real>, 3>& dimCuts) const {
+    REDEV_FUNCTION_TIMER;
+    assert(ranks.size() && cuts.size());
+    assert(dim>0 && dim<=3);
+    const auto levels = static_cast<redev::LO>(std::log2(cuts.size()));
+    const auto leafStart = static_cast<redev::LO>(std::pow(2, levels));
+    // Standard multi-path traversal of the RCB cut tree
+    struct Frame { redev::LO node; redev::LO d; };
+    std::vector<Frame> stack;
+    stack.reserve(static_cast<size_t>(levels) * 2);
+    redev::LOs cands;
+    stack.push_back({1, 0});
+    while(!stack.empty()) {
+      const auto frame = stack.back();
+      stack.pop_back();
+      if(frame.node >= leafStart) {
+        const auto rankIdx = frame.node - leafStart;
+        assert(rankIdx >= 0 && rankIdx < static_cast<redev::LO>(ranks.size()));
+        cands.push_back(ranks.at(rankIdx));
+        continue;
+      }
+      const auto v = cuts.at(frame.node);
+      const auto dNext = static_cast<redev::LO>((frame.d + 1) % dim);
+      const auto goLeft = pt[frame.d] < v;
+      const auto tol = NodeTolerance(dimCuts, frame.d, v, pt[frame.d],
+                                     tolFactor);
+      const auto nearPlane = std::abs(pt[frame.d] - v) <= tol;
+      stack.push_back({goLeft ? frame.node*2 : frame.node*2+1, dNext});
+      if(nearPlane)
+        stack.push_back({goLeft ? frame.node*2+1 : frame.node*2, dNext});
+    }
+    std::sort(cands.begin(), cands.end());
+    cands.erase(std::unique(cands.begin(), cands.end()), cands.end());
+    return cands;
+  }
+
+  redev::LOs RCBPtn::GetCandidateRanks(std::array<redev::Real,3>& pt) const {
+    REDEV_FUNCTION_TIMER;
+    const auto dimCuts = BuildDimCuts();
+    return TraverseCandidates(pt, defaultTolFactor, dimCuts);
+  }
+
+  redev::LOs RCBPtn::GetCandidateRanks(std::array<redev::Real,3>& pt,
+                                       redev::Real tolFactor) const {
+    REDEV_FUNCTION_TIMER;
+    const auto dimCuts = BuildDimCuts();
+    return TraverseCandidates(pt, tolFactor, dimCuts);
+  }
+
   void RCBPtn::Write(adios2::Engine& eng, adios2::IO& io) {
     REDEV_FUNCTION_TIMER;
     const auto len = ranks.size();

redev_partition.h

@@ -175,8 +175,62 @@ class RCBPtn {
    * are stored in the vector as (x0,y0(,z0),x1,y1(,z1),...xN-1,yN-1(,zN-1))
    */
   [[nodiscard]] std::vector<redev::Real> GetCuts() const;
+  /**
+   * Return the default tolerance factor applied to the LOCAL sub-domain scale
+   * to form the default search tolerance of GetCandidateRanks(pt) (no
+   * argument).  At each cut plane the local scale is the distance to the
+   * nearest other cut plane in the same dimension (zero when the dimension has
+   * only one cut); the tolerance is
+   * `tol = max(GetDefaultTolFactor() * localScale, machinePrecisionGuard)`,
+   * where machinePrecisionGuard keeps the tolerance above a few machine epsilons
+   * times the coordinate scale so floating-point representation error can
+   * never drive it to zero.
+   */
+  [[nodiscard]] static constexpr redev::Real GetDefaultTolFactor() noexcept {
+    return defaultTolFactor;
+  }
+  /**
+   * Return the ranks of the sub-domains that could own the given point under a
+   * (possibly different) boundary convention. The returned set is complete: every
+   * sub-domain whose closure is within tol of the point on every axis is
+   * included.  For points strictly interior to the owning sub-domain the
+   * returned vector is {GetRank(pt)}; on-boundary (and, with tol, near-
+   * boundary) points return 2..2^dim candidates (the returned vector is sorted
+   * and unique).
+   * The overload without arguments uses the default tolerance factor
+   * (GetDefaultTolFactor); pass tolFactor explicitly to override,
+   * e.g. `tolFactor = 0.0` for exact cut-plane detection only.
+   * @param[in] pt the cartesian point in space. The 3rd value is ignored for 2d
+   * domains.
+   * @param[in] tolFactor factor multiplied by the local sub-domain scale to
+   * form the search tolerance, clamped above a machine-epsilon guard
+   */
+  [[nodiscard]] redev::LOs GetCandidateRanks(std::array<redev::Real, 3>& pt) const;
+  [[nodiscard]] redev::LOs GetCandidateRanks(std::array<redev::Real, 3>& pt,
+                                             redev::Real tolFactor) const;
 
 private:
+  // Fraction of the local sub-domain scale used by the default tolerance.
+  static constexpr redev::Real defaultTolFactor = 1e-3;
+  // Machine-precision guard multiplier: floor tolerance at
+  // defaultEpsFactor * epsilon * max(|coord|, |v|, 1).
+  static constexpr redev::Real defaultEpsFactor = 64;
+
+  // Build per-dimension sorted-unique cut value arrays.
+  std::array<std::vector<redev::Real>, 3> BuildDimCuts() const;
+  // Distance from v to the nearest other cut in the same dimension.
+  static redev::Real LocalSpacing(
+      const std::array<std::vector<redev::Real>, 3>& dimCuts,
+      redev::LO d, redev::Real v, redev::Real fallback);
+  // Tolerance at node (d,v): max(tolFactor * localSpacing(d,v), epsGuard).
+  static redev::Real NodeTolerance(
+      const std::array<std::vector<redev::Real>, 3>& dimCuts,
+      redev::LO d, redev::Real v, redev::Real ptCoord,
+      redev::Real tolFactor);
+  // Fixed-radius multi-path traversal of the cut tree.
+  redev::LOs TraverseCandidates(std::array<redev::Real, 3>& pt,
+                                redev::Real tolFactor,
+                                const std::array<std::vector<redev::Real>, 3>& dimCuts) const;
   const std::string ranksVarName = "rcb partition ranks";
   const std::string cutsVarName = "rcb partition cuts";
   const std::string dimVarName = "rcb partition dim";

test_boundaryRouting.cpp

@@ -0,0 +1,162 @@
+#include <algorithm>  // std::count
+#include <array>    // std::array
+#include <cmath>    // std::nextafter
+#include <cstdlib>  // EXIT_FAILURE / EXIT_SUCCESS
+#include <iostream>
+#include "redev.h"
+
+/**
+ * Distributed point-search routing contract test.
+ *
+ * The scenario
+ * ============
+ * The rendezvous app owns the REDEZVOUS partition used to route query points to
+ * ranks (handcoded here as a 2D RCBPtn over [0,1]^2).  The receiving
+ * application owns the TRUE MESH partition (handcoded 2x2 blocks) whose cut
+ * positions (x=0.49, y=0.76) are deliberately OFFSET from the rendezvous cuts
+ * (x=0.5 and y=0.75/0.25) to model a realistic mesh whose blocks do not
+ * coincide with the rendezvous grid.
+ *
+ * Expected result (as of writing)
+ * ===============================
+ *  - Part 1 (rendezvous determinism): PASSES.  redev's answers are
+ *    deterministic and follow the documented half-open convention
+ *    (pt[d] < cut -> left, pt[d] >= cut -> right).
+ *  - Part 2 (candidate routing contract)      : FAILS on 4 of 14 rows.  The
+ *    destination set is now formed with RCBPtn::GetCandidateRanks (send-to-all-
+ *    candidates routing).  Rows that still FAIL are exactly those where the
+ *    true mesh owner's rendezvous sub-domain does NOT share a boundary with
+ *    the rendezvous cell containing the point - i.e. rendezvous/mesh
+ *    misalignment beyond the rendezvous adjacency structure.  This is the
+ *    intended signal: multi-candidate routing resolves boundary AMBIGUITY
+ *    within the rendezvous partition but cannot heal MISALIGNMENT between the
+ *    rendezvous and mesh partitions.
+ *
+ * Registered in CMakeLists.txt with WILL_FAIL TRUE so CI stays green while the
+ * known gap exists.  Remove WILL_FAIL (and expect this to pass) once the
+ * misaligned case is addressed, e.g. via cell->rank-overlap routing or
+ * rendezvous partitions aligned with the mesh.
+ */
+namespace {
+  /**
+   * Handcoded TRUE mesh partition: 2x2 axis-aligned blocks over [0,1]^2.
+   * Cut positions (x=0.49, y=0.76) are deliberately OFFSET from the rendezvous
+   * cuts.  Mesh convention: half-open blocks, x < 0.49 -> left column,
+   * y < 0.76 -> bottom row.  rank = 2*col + row.
+   */
+  struct MeshOwnership {
+    static redev::LO GetRank(const std::array<redev::Real, 3>& pt) {
+      redev::LO col = (pt[0] < 0.49) ? 0 : 1;
+      redev::LO row = (pt[1] < 0.76) ? 0 : 1;
+      return static_cast<redev::LO>(2 * col + row);
+    }
+  };
+
+  /** One handcoded query point and its hand-verified expected answers. */
+  struct Query {
+    const std::array<redev::Real, 3> pt;
+    const char* desc;
+    /** Hand-verified singular rendezvous answer (rdzv convention). */
+    redev::LO rdzvExpected;
+    /** Hand-verified true mesh owner (mesh convention). */
+    redev::LO meshExpected;
+  };
+
+  // 2D rendezvous partition (RCBPtn): regions are
+  //   rank0=[0,.5)x[0,.75)  rank1=[0,.5)x[.75,1]
+  //   rank2=[.5,1]x[0,.25)  rank3=[.5,1]x[.25,1]
+  // 2D mesh partition: regions are (xcut=0.49, ycut=0.76)
+  //   mesh0=[0,.49)x[0,.76)  mesh1=[0,.49)x[.76,1]
+  //   mesh2=[.49,1]x[0,.76)  mesh3=[.49,1]x[.76,1]
+  const std::vector<Query> queries = {
+      // --- Interior points (Part 2 must PASS today) ---
+      {{0.1, 0.1, 0.0}, "interior (both partitions)", 0, 0},
+      {{0.3, 0.3, 0.0}, "interior (both partitions)", 0, 0},
+      {{0.4, 0.8, 0.0}, "interior (both partitions)", 1, 1},
+      {{0.6, 0.1, 0.0}, "interior (both partitions)", 2, 2},
+      {{0.8, 0.9, 0.0}, "interior (both partitions)", 3, 3},
+      {{0.48, 0.3, 0.0}, "just inside both (left of mesh-cut 0.49)", 0, 0},
+      // --- Exactly on a rendezvous cut/corner (EXPECTED to FAIL) ---
+      {{0.5, 0.5, 0.0}, "EXACTLY on rdzv x-cut (0.5)", 3, 2},
+      {{0.25, 0.75, 0.0}, "EXACTLY on rdzv y-cut (0.75)", 1, 0},
+      {{0.75, 0.25, 0.0}, "EXACTLY on rdzv y-cut (0.25)", 3, 2},
+      {{0.5, 0.25, 0.0}, "EXACTLY on rdzv 4-region corner", 3, 2},
+      {{0.5, 0.75, 0.0}, "on rdzv x-cut AND y-cut", 3, 2},
+      // --- FP-slop / misaligned-strip probes (EXPECTED to FAIL) ---
+      {{0.495, 0.5, 0.0}, "rdzv-interior BUT in stagger strip (0.49,0.5)", 0, 2},
+      {{std::nextafter(0.5, 0.0), 0.5, 0.0},
+       "one ulp below rdzv x-cut (FP slop)", 0, 2},
+      // --- Educational: on a rdzv cut yet coincidentally consistent ---
+      {{0.5, 0.9, 0.0}, "on rdzv x-cut but mesh row agrees", 3, 3},
+  };
+} // namespace
+int main(int argc, char** argv) {
+  MPI_Init(&argc, &argv);
+  {
+    // Handcoded rendezvous partition.  cuts: slot 0 unused, then BFS order
+    // (root x=0.5, then y=0.75 at node 2, y=0.25 at node 3).
+    const redev::LO dim = 2;
+    redev::LOs ranks = {0, 1, 2, 3};
+    redev::Reals cuts = {0, 0.5, 0.75, 0.25};
+    redev::RCBPtn rdzvPtn(dim, ranks, cuts);
+    redev::Redev rdv(MPI_COMM_WORLD, redev::Partition{rdzvPtn},
+                     redev::ProcessType::Server, /*noClients=*/true);
+
+    std::cout << "test_boundaryRouting: handcoded rendezvous (RCBPtn) vs "
+                 "handcoded mesh (offset 2x2 blocks)\n"
+                 "row format: PASS/FAIL  candidates={dest ranks}  meshOwner=#  "
+                 "description\n\n";
+
+    // ---------------- Part 1: rendezvous determinism ----------------
+    // redev's GetRank must be deterministic and identical whether queried
+    // directly or through the partition stored inside Redev.
+    for (const auto& q : queries) {
+      std::array<redev::Real, 3> pt = q.pt;
+      const auto destLocal = rdzvPtn.GetRank(pt);
+      const auto destRedev =
+          std::get<redev::RCBPtn>(rdv.GetPartition()).GetRank(pt);
+      REDEV_ALWAYS_ASSERT(destLocal == q.rdzvExpected);
+      REDEV_ALWAYS_ASSERT(destRedev == destLocal);
+    }
+    std::cout << "Part 1 (rendezvous determinism + half-open convention): "
+                 "PASSED\n\n";
+
+    // ---------------- Part 2: candidate routing contract ----------------
+    // Contract: the destination-rank candidate set returned by
+    // GetCandidateRanks must contain the TRUE mesh owner.
+    std::size_t numUnrouted = 0;
+    for (const auto& q : queries) {
+      std::array<redev::Real, 3> pt = q.pt;
+      auto dest = std::get<redev::RCBPtn>(rdv.GetPartition())
+                      .GetCandidateRanks(pt, /*tolFactor=*/0.0);
+      const auto ownerMesh = MeshOwnership::GetRank(pt);
+      REDEV_ALWAYS_ASSERT(std::count(dest.begin(), dest.end(), q.rdzvExpected) != 0);
+      REDEV_ALWAYS_ASSERT(ownerMesh == q.meshExpected);
+      const bool routed = std::count(dest.begin(), dest.end(), ownerMesh) != 0;
+      if (!routed) ++numUnrouted;
+      std::cout << (routed ? "PASS" : "FAIL") << "  candidates={";
+      for (std::size_t i = 0; i < dest.size(); i++)
+        std::cout << (i ? "," : "") << dest[i];
+      std::cout << "}  meshOwner=" << ownerMesh << "  " << q.desc << "\n";
+    }
+
+    if (numUnrouted) {
+      std::cerr
+          << "\nPart 2 (candidate routing contract): " << numUnrouted << " of "
+          << queries.size()
+          << " query points have a TRUE mesh owner that is NOT in the candidate\n"
+          << "set returned by GetCandidateRanks.  That set contains only the\n"
+          << "rendezvous sub-domains sharing a boundary with the point, so these\n"
+          << "failures are rows where the mesh owner lies in a different\n"
+          << "rendezvous region (misalignment beyond adjacency).  This is the\n"
+          << "intended signal: multi-candidate routing resolves boundary\n"
+          << "AMBIGUITY within the rendezvous partition but cannot heal\n"
+          << "rendezvous/mesh MISALIGNMENT.\n";
+      MPI_Finalize();
+      return EXIT_FAILURE;
+    }
+    std::cout << "Part 2 (candidate routing contract) PASSED\n";
+  }
+  MPI_Finalize();
+  return EXIT_SUCCESS;
+}