Merge pull request #1367 from parthenon-hpc-lab/lroberts36/update-dempsey/bvals

Updates to #1365

Author: lroberts36

src/CMakeLists.txt

@@ -187,6 +187,7 @@ add_library(parthenon
   mesh/mesh-gmg.cpp
   mesh/mesh.cpp
   mesh/mesh.hpp
+  mesh/mesh_neighbors.hpp
   mesh/meshblock.hpp
   mesh/meshblock_pack.hpp
   mesh/meshblock.cpp

src/bvals/boundary_conditions.cpp

@@ -43,17 +43,7 @@ TaskStatus ApplyBoundaryConditionsOnCoarseOrFine(std::shared_ptr<MeshBlockData<R
 
   // We only need to call the BC on the coarse buffer if one of the neighbors
   // is at a coarser level than us.
-  if (coarse) {
-    bool has_coarser_neighbor = false;
-    const int mylevel = pmb->loc.level();
-    for (const auto &nb : pmb->neighbors) {
-      if (nb.origin_loc.level() < mylevel) {
-        has_coarser_neighbor = true;
-        break;
-      }
-    }
-    if (!has_coarser_neighbor) return TaskStatus::complete;
-  }
+  if (coarse && !pmb->HasCoarserNeighbors()) return TaskStatus::complete;
 
   auto &tree_bnd_func = pmesh->forest.GetTreePtr(pmb->loc.tree())->MeshBndryFnctn;
   auto &tree_bnd_func_user =

src/bvals/bvals.cpp

@@ -70,8 +70,9 @@ void BoundarySwarm::SetupPersistentMPI() {
   std::shared_ptr<MeshBlock> pmb = GetBlockPointer();
 
   // Initialize neighbor communications to other ranks
-  for (int n = 0; n < pmb->neighbors.size(); n++) {
-    NeighborBlock &nb = pmb->neighbors[n];
+  const auto &neighbors = pmb->GetNeighbors();
+  for (int n = 0; n < neighbors.size(); n++) {
+    const NeighborBlock &nb = neighbors[n];
     // Neighbor on different MPI process
     if (nb.rank != Globals::my_rank) {
       send_tag[nb.bufid] = pmb->pmy_mesh->tag_map.GetTag(pmb.get(), nb);
@@ -93,8 +94,9 @@ void BoundarySwarm::Send(BoundaryCommSubset phase) {
   std::shared_ptr<MeshBlock> pmb = GetBlockPointer();
   // Fence to make sure buffers are loaded before sending
   pmb->exec_space.fence();
-  for (int n = 0; n < pmb->neighbors.size(); n++) {
-    NeighborBlock &nb = pmb->neighbors[n];
+  const auto &neighbors = pmb->GetNeighbors();
+  for (int n = 0; n < neighbors.size(); n++) {
+    const NeighborBlock &nb = neighbors[n];
     if (nb.rank != Globals::my_rank) {
 #ifdef MPI_PARALLEL
       PARTHENON_REQUIRE(bd_var_.req_send[nb.bufid] == MPI_REQUEST_NULL,
@@ -131,8 +133,9 @@ void BoundarySwarm::Receive(BoundaryCommSubset phase) {
 #ifdef MPI_PARALLEL
   std::shared_ptr<MeshBlock> pmb = GetBlockPointer();
   const int &mylevel = pmb->loc.level();
-  for (int n = 0; n < pmb->neighbors.size(); n++) {
-    NeighborBlock &nb = pmb->neighbors[n];
+  const auto &neighbors = pmb->GetNeighbors();
+  for (int n = 0; n < neighbors.size(); n++) {
+    const NeighborBlock &nb = neighbors[n];
     if (nb.rank != Globals::my_rank) {
       // Check to see if we got a message
       int test;

src/bvals/comms/bnd_info.cpp

@@ -425,17 +425,11 @@ ProResInfo ProResInfo::GetSet(MeshBlock *pmb, const NeighborBlock &nb,
   // This will select a superset of the boundaries that actually need to be restricted,
   // more logic could be added to only restrict boundary regions that abut boundary
   // regions that were filled by coarser neighbors
-  bool restricted = false;
-  int mylevel = pmb->loc.level();
-  if (mylevel > 0) {
-    for (const auto &nb : pmb->neighbors) {
-      restricted = restricted || (nb.origin_loc.level() == (mylevel - 1));
-    }
-  }
+  bool restricted = pmb->HasCoarserNeighbors();
 
   for (auto el : v->GetTopologicalElements()) {
     out.IncludeTopoEl(el) = true;
-    if (nb.origin_loc.level() < mylevel) {
+    if (nb.origin_loc.level() < pmb->loc.level()) {
       out.refinement_op = RefinementOp_t::Prolongation;
     } else {
       if (restricted) {
@@ -455,7 +449,7 @@ ProResInfo ProResInfo::GetSet(MeshBlock *pmb, const NeighborBlock &nb,
   //      I doubt that the extra calculations matter, but the storage overhead could
   //      matter since each 6D indexer contains 18 ints and we are always carrying around
   //      10 indexers per bound info even if the field isn't allocated
-  if (nb.origin_loc.level() < mylevel) {
+  if (nb.origin_loc.level() < pmb->loc.level()) {
     for (auto el : {TE::CC, TE::F1, TE::F2, TE::F3, TE::E1, TE::E2, TE::E3, TE::NN})
       out.idxer[static_cast<int>(el)] =
           CalcIndices(nb, pmb, v, el, IndexRangeType::BoundaryExteriorRecv, true);

src/bvals/comms/boundary_communication.cpp

@@ -233,7 +233,7 @@ TaskStatus ReceiveBoundBufs(std::shared_ptr<MeshData<Real>> &md) {
   int ibound = 0;
   if (Globals::sparse_config.enabled && all_received) {
     ForEachBoundary<bound_type>(
-        md, [&](auto pmb, sp_mbd_t rc, nb_t &nb, const sp_cv_t v) {
+        md, [&](auto pmb, sp_mbd_t rc, const nb_t &nb, const sp_cv_t v) {
           const std::size_t ibuf = cache.idx_vec[ibound];
           auto &buf = *cache.buf_vec[ibuf];
 

src/bvals/comms/build_boundary_buffers.cpp

@@ -54,20 +54,23 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
   std::unordered_map<std::size_t, std::size_t>
       nbufs_allocated; // total that are actually allocated
 
-  ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, nb_t &nb, const sp_cv_t v) {
-    // Calculate the required size of the buffer for this boundary
-    int buf_size = GetBufferSize(pmb, nb, v);
-    //  LR: Multigrid logic requires blocks sending messages to themselves (since the same
-    //  block can show up on two multigrid levels). This doesn't require any data
-    //  transfer, so the message size can be zero. It is essentially just a flag to show
-    //  that the block is done being used on one level and can be used on the next level.
-    if (pmb->gid == nb.gid && nb.offsets.IsCell()) buf_size = 0;
-
-    nbufs[buf_size] += 1; // relying on value init of int to 0 for initial entry
-    nbufs_allocated[buf_size] += v->IsAllocated();
-  });
-
-  ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, nb_t &nb, const sp_cv_t v) {
+  ForEachBoundary<BTYPE>(
+      md, [&](auto pmb, sp_mbd_t /*rc*/, const nb_t &nb, const sp_cv_t v) {
+        // Calculate the required size of the buffer for this boundary
+        int buf_size = GetBufferSize(pmb, nb, v);
+        //  LR: Multigrid logic requires blocks sending messages to themselves (since the
+        //  same block can show up on two multigrid levels). This doesn't require any data
+        //  transfer, so the message size can be zero. It is essentially just a flag to
+        //  show that the block is done being used on one level and can be used on the
+        //  next level.
+        if (pmb->gid == nb.gid && nb.offsets.IsCell()) buf_size = 0;
+
+        nbufs[buf_size] += 1; // relying on value init of int to 0 for initial entry
+        nbufs_allocated[buf_size] += v->IsAllocated();
+      });
+
+  ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, const nb_t &nb,
+                                 const sp_cv_t v) {
     // Calculate the required size of the buffer for this boundary
     int buf_size = GetBufferSize(pmb, nb, v);
     // See comment above on the same logic.
@@ -139,18 +142,19 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
 template <BoundaryType BTYPE>
 void RegisterCoalescedCommsSubset(std::shared_ptr<MeshData<Real>> &md) {
   Mesh *pmesh = md->GetMeshPointer();
-  ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, nb_t &nb, const sp_cv_t v) {
-    const int receiver_rank = nb.rank;
-    const int sender_rank = Globals::my_rank;
-    if (receiver_rank != sender_rank) {
-      if constexpr (IsSender(BTYPE)) {
-        pmesh->pcoalesced_comms->AddSendBuffer(md->partition, pmb, nb, v, BTYPE);
-      }
-      if constexpr (IsReceiver(BTYPE)) {
-        pmesh->pcoalesced_comms->AddRecvBuffer(pmb, nb, v, BTYPE);
-      }
-    }
-  });
+  ForEachBoundary<BTYPE>(
+      md, [&](auto pmb, sp_mbd_t /*rc*/, const nb_t &nb, const sp_cv_t v) {
+        const int receiver_rank = nb.rank;
+        const int sender_rank = Globals::my_rank;
+        if (receiver_rank != sender_rank) {
+          if constexpr (IsSender(BTYPE)) {
+            pmesh->pcoalesced_comms->AddSendBuffer(md->partition, pmb, nb, v, BTYPE);
+          }
+          if constexpr (IsReceiver(BTYPE)) {
+            pmesh->pcoalesced_comms->AddRecvBuffer(pmb, nb, v, BTYPE);
+          }
+        }
+      });
 }
 
 } // namespace

src/bvals/comms/bvals_utils.hpp

@@ -99,16 +99,17 @@ void InitializeBufferCache(std::shared_ptr<MeshData<Real>> &md, COMM_MAP *comm_m
   std::vector<std::tuple<int, int, Mesh::channel_key_t>> key_order;
 
   int boundary_idx = 0;
-  ForEachBoundary<bound_type>(md, [&](auto pmb, sp_mbd_t rc, nb_t &nb, const sp_cv_t v) {
-    auto key = KeyFunc(pmb, nb, v, bound_type);
-    PARTHENON_DEBUG_REQUIRE(comm_map->count(key) > 0,
-                            "Boundary communicator does not exist");
-    // Create a unique index by combining receiver gid (second element of the key
-    // tuple) and geometric element index (fourth element of the key tuple)
-    int recvr_idx = 27 * GetReceiverGid(key) + GetLocIdx(key);
-    key_order.push_back({recvr_idx, boundary_idx, key});
-    ++boundary_idx;
-  });
+  ForEachBoundary<bound_type>(
+      md, [&](auto pmb, sp_mbd_t rc, const nb_t &nb, const sp_cv_t v) {
+        auto key = KeyFunc(pmb, nb, v, bound_type);
+        PARTHENON_DEBUG_REQUIRE(comm_map->count(key) > 0,
+                                "Boundary communicator does not exist");
+        // Create a unique index by combining receiver gid (second element of the key
+        // tuple) and geometric element index (fourth element of the key tuple)
+        int recvr_idx = 27 * GetReceiverGid(key) + GetLocIdx(key);
+        key_order.push_back({recvr_idx, boundary_idx, key});
+        ++boundary_idx;
+      });
 
   // If desired, sort the keys and boundary indices by receiver_idx
   // std::sort(key_order.begin(), key_order.end(),
@@ -155,7 +156,8 @@ inline auto CheckSendBufferCacheForRebuild(std::shared_ptr<MeshData<Real>> md) {
   bool rebuild = false;
   bool other_communication_unfinished = false;
   int nbound = 0;
-  ForEachBoundary<BOUND_TYPE>(md, [&](auto pmb, sp_mbd_t rc, nb_t &nb, const sp_cv_t v) {
+  ForEachBoundary<BOUND_TYPE>(md, [&](auto pmb, sp_mbd_t rc, const nb_t &nb,
+                                      const sp_cv_t v) {
     const std::size_t ibuf = cache.idx_vec[nbound];
     auto &buf = *(cache.buf_vec[ibuf]);
 
@@ -187,7 +189,8 @@ inline auto CheckReceiveBufferCacheForRebuild(std::shared_ptr<MeshData<Real>> md
   bool rebuild = false;
   int nbound = 0;
 
-  ForEachBoundary<BOUND_TYPE>(md, [&](auto pmb, sp_mbd_t rc, nb_t &nb, const sp_cv_t v) {
+  ForEachBoundary<BOUND_TYPE>(md, [&](auto pmb, sp_mbd_t rc, const nb_t &nb,
+                                      const sp_cv_t v) {
     const std::size_t ibuf = cache.idx_vec[nbound];
     auto &buf = *cache.buf_vec[ibuf];
     if (ibuf < cache.bnd_info_h.size()) {
@@ -237,19 +240,20 @@ inline void RebuildBufferCache(std::shared_ptr<MeshData<Real>> md, int nbound,
   cache.prores_cache.Initialize(nbound, pkg);
 
   int ibound = 0;
-  ForEachBoundary<BOUND_TYPE>(md, [&](auto pmb, sp_mbd_t rc, nb_t &nb, const sp_cv_t v) {
-    // bnd_info
-    const std::size_t ibuf = cache.idx_vec[ibound];
-    cache.bnd_info_h(ibuf) = BndInfoCreator(pmb, nb, v, cache.buf_vec[ibuf]);
-
-    // subsets ordering is same as in cache.bnd_info
-    // RefinementFunctions_t owns all relevant functionality, so
-    // only one ParArray2D needed.
-    cache.prores_cache.RegisterRegionHost(ibuf, ProResInfoCreator(pmb, nb, v), v.get(),
-                                          pkg);
-
-    ++ibound;
-  });
+  ForEachBoundary<BOUND_TYPE>(
+      md, [&](auto pmb, sp_mbd_t rc, const nb_t &nb, const sp_cv_t v) {
+        // bnd_info
+        const std::size_t ibuf = cache.idx_vec[ibound];
+        cache.bnd_info_h(ibuf) = BndInfoCreator(pmb, nb, v, cache.buf_vec[ibuf]);
+
+        // subsets ordering is same as in cache.bnd_info
+        // RefinementFunctions_t owns all relevant functionality, so
+        // only one ParArray2D needed.
+        cache.prores_cache.RegisterRegionHost(ibuf, ProResInfoCreator(pmb, nb, v),
+                                              v.get(), pkg);
+
+        ++ibound;
+      });
   Kokkos::deep_copy(cache.bnd_info, cache.bnd_info_h);
   cache.prores_cache.CopyToDevice();
 }

src/bvals/comms/tag_map.cpp

@@ -43,26 +43,26 @@ void TagMap::AddMeshDataToMap(std::shared_ptr<MeshData<Real>> &md) {
     auto pmb = rc->GetBlockPointer();
     // type_t var = []{...}() pattern defines and uses a lambda that
     // returns  to reduce initializations of var
-    auto *neighbors = [&pmb, &md] {
+    const auto &neighbors = [&pmb, &md] {
       if constexpr (BOUND == BoundaryType::gmg_restrict_send)
-        return pmb->loc.level() == md->grid.logical_level ? &(pmb->gmg_coarser_neighbors)
-                                                          : &(pmb->gmg_leaf_neighbors);
+        return pmb->loc.level() == md->grid.logical_level ? pmb->GetGMGCoarserNeighbors()
+                                                          : pmb->GetGMGLeafNeighbors();
       if constexpr (BOUND == BoundaryType::gmg_restrict_recv)
-        return pmb->gmg_finer_neighbors.size() > 0 ? &(pmb->gmg_finer_neighbors)
-                                                   : &(pmb->gmg_leaf_neighbors);
+        return pmb->GetGMGFinerNeighbors().size() > 0 ? pmb->GetGMGFinerNeighbors()
+                                                      : pmb->GetGMGLeafNeighbors();
       if constexpr (BOUND == BoundaryType::gmg_prolongate_send)
-        return &(pmb->gmg_finer_neighbors);
+        return pmb->GetGMGFinerNeighbors();
       if constexpr (BOUND == BoundaryType::gmg_prolongate_recv)
-        return &(pmb->gmg_coarser_neighbors);
+        return pmb->GetGMGCoarserNeighbors();
       if constexpr (BOUND == BoundaryType::gmg_prolongate_recv)
-        return &(pmb->gmg_coarser_neighbors);
+        return pmb->GetGMGCoarserNeighbors();
       if constexpr (BOUND == BoundaryType::gmg_same)
         return pmb->loc.level() == md->grid.logical_level
-                   ? &(pmb->gmg_same_neighbors)
-                   : &(pmb->gmg_composite_finer_neighbors);
-      return &(pmb->neighbors);
+                   ? pmb->GetGMGSameNeighbors()
+                   : pmb->GetGMGCompositeFinerNeighbors();
+      return pmb->GetNeighbors();
     }();
-    for (auto &nb : *neighbors) {
+    for (const auto &nb : neighbors) {
       const int other_rank = nb.rank;
       if (map_.count(other_rank) < 1) map_[other_rank] = rank_pair_map_t();
       auto &pair_map = map_[other_rank];