Merge branch 'develop' into brryan/fix_vec_hst

CHANGELOG.md

@@ -28,6 +28,7 @@
 
 
 ### Incompatibilities (i.e. breaking changes)
+- [[PR 1042]](https://github.com/parthenon-hpc-lab/parthenon/pull/1042) Use Offset class and clean up of NeighborBlock
 - [[PR1019](https://github.com/parthenon-hpc-lab/parthenon/pull/1019) Remove support for file formats < 3
 
 

src/bvals/bvals.cpp

@@ -73,7 +73,7 @@ void BoundarySwarm::SetupPersistentMPI() {
   for (int n = 0; n < pmb->neighbors.size(); n++) {
     NeighborBlock &nb = pmb->neighbors[n];
     // Neighbor on different MPI process
-    if (nb.snb.rank != Globals::my_rank) {
+    if (nb.rank != Globals::my_rank) {
       send_tag[nb.bufid] = pmb->pmy_mesh->tag_map.GetTag(pmb.get(), nb);
       recv_tag[nb.bufid] = pmb->pmy_mesh->tag_map.GetTag(pmb.get(), nb);
       if (bd_var_.req_send[nb.bufid] != MPI_REQUEST_NULL) {
@@ -95,16 +95,16 @@ void BoundarySwarm::Send(BoundaryCommSubset phase) {
   pmb->exec_space.fence();
   for (int n = 0; n < pmb->neighbors.size(); n++) {
     NeighborBlock &nb = pmb->neighbors[n];
-    if (nb.snb.rank != Globals::my_rank) {
+    if (nb.rank != Globals::my_rank) {
 #ifdef MPI_PARALLEL
       PARTHENON_REQUIRE(bd_var_.req_send[nb.bufid] == MPI_REQUEST_NULL,
                         "Trying to create a new send before previous send completes!");
       PARTHENON_MPI_CHECK(MPI_Isend(bd_var_.send[nb.bufid].data(), send_size[nb.bufid],
-                                    MPI_PARTHENON_REAL, nb.snb.rank, send_tag[nb.bufid],
+                                    MPI_PARTHENON_REAL, nb.rank, send_tag[nb.bufid],
                                     swarm_comm, &(bd_var_.req_send[nb.bufid])));
 #endif // MPI_PARALLEL
     } else {
-      MeshBlock &target_block = *pmy_mesh_->FindMeshBlock(nb.snb.gid);
+      MeshBlock &target_block = *pmy_mesh_->FindMeshBlock(nb.gid);
       std::shared_ptr<BoundarySwarm> ptarget_bswarm =
           target_block.pbswarm->bswarms[bswarm_index];
       if (send_size[nb.bufid] > 0) {
@@ -133,14 +133,14 @@ void BoundarySwarm::Receive(BoundaryCommSubset phase) {
   const int &mylevel = pmb->loc.level();
   for (int n = 0; n < pmb->neighbors.size(); n++) {
     NeighborBlock &nb = pmb->neighbors[n];
-    if (nb.snb.rank != Globals::my_rank) {
+    if (nb.rank != Globals::my_rank) {
       // Check to see if we got a message
       int test;
       MPI_Status status;
 
       if (bd_var_.flag[nb.bufid] != BoundaryStatus::completed) {
         PARTHENON_MPI_CHECK(
-            MPI_Iprobe(nb.snb.rank, recv_tag[nb.bufid], swarm_comm, &test, &status));
+            MPI_Iprobe(nb.rank, recv_tag[nb.bufid], swarm_comm, &test, &status));
         if (!static_cast<bool>(test)) {
           bd_var_.flag[nb.bufid] = BoundaryStatus::waiting;
         } else {
@@ -153,8 +153,8 @@ void BoundarySwarm::Receive(BoundaryCommSubset phase) {
             bd_var_.recv[nb.bufid] = BufArray1D<Real>("Buffer", recv_size[nb.bufid]);
           }
           PARTHENON_MPI_CHECK(MPI_Recv(bd_var_.recv[nb.bufid].data(), recv_size[nb.bufid],
-                                       MPI_PARTHENON_REAL, nb.snb.rank,
-                                       recv_tag[nb.bufid], swarm_comm, &status));
+                                       MPI_PARTHENON_REAL, nb.rank, recv_tag[nb.bufid],
+                                       swarm_comm, &status));
         }
       }
     }

src/bvals/comms/bnd_info.cpp

@@ -80,7 +80,6 @@ void ProResCache_t::RegisterRegionHost(int region, ProResInfo pri, Variable<Real
 SpatiallyMaskedIndexer6D CalcIndices(const NeighborBlock &nb, MeshBlock *pmb,
                                      TopologicalElement el, IndexRangeType ir_type,
                                      bool prores, std::array<int, 3> tensor_shape) {
-  const auto &ni = nb.ni;
   const auto &loc = pmb->loc;
   auto shape = pmb->cellbounds;
   // Both prolongation and restriction always operate in the coarse
@@ -114,7 +113,7 @@ SpatiallyMaskedIndexer6D CalcIndices(const NeighborBlock &nb, MeshBlock *pmb,
   // the current block in some cases
   std::array<int, 3> top_offset{TopologicalOffsetI(el), TopologicalOffsetJ(el),
                                 TopologicalOffsetK(el)};
-  std::array<int, 3> block_offset = {ni.ox1, ni.ox2, ni.ox3};
+  std::array<int, 3> block_offset = nb.offsets;
 
   int interior_offset =
       ir_type == IndexRangeType::BoundaryInteriorSend ? Globals::nghost : 0;
@@ -184,9 +183,9 @@ SpatiallyMaskedIndexer6D CalcIndices(const NeighborBlock &nb, MeshBlock *pmb,
   // but it is possible to convince oneself that dealing with ownership in
   // only exterior index ranges works correctly
   if (ir_type == IndexRangeType::BoundaryExteriorRecv) {
-    int sox1 = -ni.ox1;
-    int sox2 = -ni.ox2;
-    int sox3 = -ni.ox3;
+    int sox1 = -block_offset[0];
+    int sox2 = -block_offset[1];
+    int sox3 = -block_offset[2];
     if (nb.loc.level() < loc.level()) {
       // For coarse to fine interfaces, we are passing zones from only an
       // interior corner of the cell, never an entire face or edge
@@ -214,10 +213,10 @@ int GetBufferSize(MeshBlock *pmb, const NeighborBlock &nb,
   const int isize = cb.ie(in) - cb.is(in) + 2;
   const int jsize = cb.je(in) - cb.js(in) + 2;
   const int ksize = cb.ke(in) - cb.ks(in) + 2;
-  return (nb.ni.ox1 == 0 ? isize : Globals::nghost + 1) *
-         (nb.ni.ox2 == 0 ? jsize : Globals::nghost + 1) *
-         (nb.ni.ox3 == 0 ? ksize : Globals::nghost + 1) * v->GetDim(6) * v->GetDim(5) *
-         v->GetDim(4) * topo_comp;
+  return (nb.offsets(X1DIR) == 0 ? isize : Globals::nghost + 1) *
+         (nb.offsets(X2DIR) == 0 ? jsize : Globals::nghost + 1) *
+         (nb.offsets(X3DIR) == 0 ? ksize : Globals::nghost + 1) * v->GetDim(6) *
+         v->GetDim(5) * v->GetDim(4) * topo_comp;
 }
 
 BndInfo BndInfo::GetSendBndInfo(MeshBlock *pmb, const NeighborBlock &nb,
@@ -239,13 +238,12 @@ BndInfo BndInfo::GetSendBndInfo(MeshBlock *pmb, const NeighborBlock &nb,
   auto elements = v->GetTopologicalElements();
   out.ntopological_elements = elements.size();
   auto idx_range_type = IndexRangeType::BoundaryInteriorSend;
-  if (std::abs(nb.ni.ox1) + std::abs(nb.ni.ox2) + std::abs(nb.ni.ox3) == 0)
-    idx_range_type = IndexRangeType::InteriorSend;
+  if (nb.offsets.IsCell()) idx_range_type = IndexRangeType::InteriorSend;
   for (auto el : elements) {
     int idx = static_cast<int>(el) % 3;
     out.idxer[idx] = CalcIndices(nb, pmb, el, idx_range_type, false, {Nt, Nu, Nv});
   }
-  if (nb.snb.level < mylevel) {
+  if (nb.loc.level() < mylevel) {
     out.var = v->coarse_s.Get();
   } else {
     out.var = v->data.Get();
@@ -278,13 +276,12 @@ BndInfo BndInfo::GetSetBndInfo(MeshBlock *pmb, const NeighborBlock &nb,
   auto elements = v->GetTopologicalElements();
   out.ntopological_elements = elements.size();
   auto idx_range_type = IndexRangeType::BoundaryExteriorRecv;
-  if (std::abs(nb.ni.ox1) + std::abs(nb.ni.ox2) + std::abs(nb.ni.ox3) == 0)
-    idx_range_type = IndexRangeType::InteriorRecv;
+  if (nb.offsets.IsCell()) idx_range_type = IndexRangeType::InteriorRecv;
   for (auto el : elements) {
     int idx = static_cast<int>(el) % 3;
     out.idxer[idx] = CalcIndices(nb, pmb, el, idx_range_type, false, {Nt, Nu, Nv});
   }
-  if (nb.snb.level < mylevel) {
+  if (nb.loc.level() < mylevel) {
     out.var = v->coarse_s.Get();
   } else {
     out.var = v->data.Get();
@@ -312,8 +309,7 @@ ProResInfo ProResInfo::GetInteriorRestrict(MeshBlock *pmb, const NeighborBlock &
 
   out.fine = v->data.Get();
   out.coarse = v->coarse_s.Get();
-  NeighborBlock nb(pmb->pmy_mesh, pmb->loc, Globals::my_rank, 0, 0, {0, 0, 0},
-                   NeighborConnect::none, 0, 0, 0, 0);
+  NeighborBlock nb(pmb->pmy_mesh, pmb->loc, Globals::my_rank, 0, {0, 0, 0}, 0, 0, 0, 0);
 
   auto elements = v->GetTopologicalElements();
   out.ntopological_elements = elements.size();
@@ -344,8 +340,7 @@ ProResInfo ProResInfo::GetInteriorProlongate(MeshBlock *pmb, const NeighborBlock
 
   out.fine = v->data.Get();
   out.coarse = v->coarse_s.Get();
-  NeighborBlock nb(pmb->pmy_mesh, pmb->loc, Globals::my_rank, 0, 0, {0, 0, 0},
-                   NeighborConnect::none, 0, 0, 0, 0);
+  NeighborBlock nb(pmb->pmy_mesh, pmb->loc, Globals::my_rank, 0, {0, 0, 0}, 0, 0, 0, 0);
 
   auto elements = v->GetTopologicalElements();
   out.ntopological_elements = elements.size();
@@ -378,7 +373,7 @@ ProResInfo ProResInfo::GetSend(MeshBlock *pmb, const NeighborBlock &nb,
 
   auto elements = v->GetTopologicalElements();
   out.ntopological_elements = elements.size();
-  if (nb.snb.level < mylevel) {
+  if (nb.loc.level() < mylevel) {
     for (auto el : elements) {
       out.idxer[static_cast<int>(el)] = CalcIndices(
           nb, pmb, el, IndexRangeType::BoundaryInteriorSend, true, {Nt, Nu, Nv});
@@ -416,7 +411,7 @@ ProResInfo ProResInfo::GetSet(MeshBlock *pmb, const NeighborBlock &nb,
   auto elements = v->GetTopologicalElements();
   out.ntopological_elements = elements.size();
   for (auto el : elements) {
-    if (nb.snb.level < mylevel) {
+    if (nb.loc.level() < mylevel) {
       out.refinement_op = RefinementOp_t::Prolongation;
     } else {
       if (restricted) {
@@ -436,7 +431,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.snb.level < mylevel) {
+  if (nb.loc.level() < mylevel) {
     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, el, IndexRangeType::BoundaryExteriorRecv, true, {Nt, Nu, Nv});
@@ -468,29 +463,20 @@ BndInfo BndInfo::GetSendCCFluxCor(MeshBlock *pmb, const NeighborBlock &nb,
   int si = ib.s;
   int ei = si + std::max((ib.e - ib.s + 1) / 2, 1) - 1;
 
-  if (nb.fid == BoundaryFace::inner_x1 || nb.fid == BoundaryFace::outer_x1) {
+  PARTHENON_REQUIRE(nb.offsets.IsFace(),
+                    "Flux corrections only occur on faces for CC variables.");
+  if (nb.offsets(X1DIR) != 0) {
     out.dir = X1DIR;
-    if (nb.fid == BoundaryFace::inner_x1)
-      si = ib.s;
-    else
-      si = ib.e + 1;
+    si = nb.offsets(X1DIR) < 0 ? ib.s : ib.e + 1;
     ei = si;
-  } else if (nb.fid == BoundaryFace::inner_x2 || nb.fid == BoundaryFace::outer_x2) {
+  } else if (nb.offsets(X2DIR) != 0) {
     out.dir = X2DIR;
-    if (nb.fid == BoundaryFace::inner_x2)
-      sj = jb.s;
-    else
-      sj = jb.e + 1;
+    sj = nb.offsets(X2DIR) < 0 ? jb.s : jb.e + 1;
     ej = sj;
-  } else if (nb.fid == BoundaryFace::inner_x3 || nb.fid == BoundaryFace::outer_x3) {
+  } else if (nb.offsets(X3DIR) != 0) {
     out.dir = X3DIR;
-    if (nb.fid == BoundaryFace::inner_x3)
-      sk = kb.s;
-    else
-      sk = kb.e + 1;
+    sk = nb.offsets(X3DIR) < 0 ? kb.s : kb.e + 1;
     ek = sk;
-  } else {
-    PARTHENON_FAIL("Flux corrections only occur on faces for CC variables.");
   }
 
   out.var = v->flux[out.dir];
@@ -526,45 +512,47 @@ BndInfo BndInfo::GetSetCCFluxCor(MeshBlock *pmb, const NeighborBlock &nb,
   int ek = kb.e;
   int ej = jb.e;
   int ei = ib.e;
-  if (nb.fid == BoundaryFace::inner_x1 || nb.fid == BoundaryFace::outer_x1) {
+  PARTHENON_REQUIRE(nb.offsets.IsFace(),
+                    "Flux corrections only occur on faces for CC variables.");
+  if (nb.offsets(X1DIR) != 0) {
     out.dir = X1DIR;
-    if (nb.fid == BoundaryFace::inner_x1)
+    if (nb.offsets(X1DIR) == -1)
       ei = si;
     else
       si = ++ei;
-    if (nb.ni.fi1 == 0)
+    if (nb.fi1 == 0)
       ej -= pmb->block_size.nx(X2DIR) / 2;
     else
       sj += pmb->block_size.nx(X2DIR) / 2;
-    if (nb.ni.fi2 == 0)
+    if (nb.fi2 == 0)
       ek -= pmb->block_size.nx(X3DIR) / 2;
     else
       sk += pmb->block_size.nx(X3DIR) / 2;
-  } else if (nb.fid == BoundaryFace::inner_x2 || nb.fid == BoundaryFace::outer_x2) {
+  } else if (nb.offsets(X2DIR) != 0) {
     out.dir = X2DIR;
-    if (nb.fid == BoundaryFace::inner_x2)
+    if (nb.offsets(X2DIR) == -1)
       ej = sj;
     else
       sj = ++ej;
-    if (nb.ni.fi1 == 0)
+    if (nb.fi1 == 0)
       ei -= pmb->block_size.nx(X1DIR) / 2;
     else
       si += pmb->block_size.nx(X1DIR) / 2;
-    if (nb.ni.fi2 == 0)
+    if (nb.fi2 == 0)
       ek -= pmb->block_size.nx(X3DIR) / 2;
     else
       sk += pmb->block_size.nx(X3DIR) / 2;
-  } else if (nb.fid == BoundaryFace::inner_x3 || nb.fid == BoundaryFace::outer_x3) {
+  } else if (nb.offsets(X3DIR) != 0) {
     out.dir = X3DIR;
-    if (nb.fid == BoundaryFace::inner_x3)
+    if (nb.offsets(X3DIR) == -1)
       ek = sk;
     else
       sk = ++ek;
-    if (nb.ni.fi1 == 0)
+    if (nb.fi1 == 0)
       ei -= pmb->block_size.nx(X1DIR) / 2;
     else
       si += pmb->block_size.nx(X1DIR) / 2;
-    if (nb.ni.fi2 == 0)
+    if (nb.fi2 == 0)
       ej -= pmb->block_size.nx(X2DIR) / 2;
     else
       sj += pmb->block_size.nx(X2DIR) / 2;

src/bvals/comms/build_boundary_buffers.cpp

@@ -64,7 +64,7 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
           })));
     }
 
-    const int receiver_rank = nb.snb.rank;
+    const int receiver_rank = nb.rank;
     const int sender_rank = Globals::my_rank;
 
     int tag = 0;

src/bvals/comms/bvals_utils.hpp

@@ -39,17 +39,17 @@ inline std::tuple<int, int, std::string, int>
 SendKey(const MeshBlock *pmb, const NeighborBlock &nb,
         const std::shared_ptr<Variable<Real>> &pcv) {
   const int sender_id = pmb->gid;
-  const int receiver_id = nb.snb.gid;
-  const int location_idx = (1 + nb.ni.ox1) + 3 * (1 + nb.ni.ox2 + 3 * (1 + nb.ni.ox3));
+  const int receiver_id = nb.gid;
+  const int location_idx = nb.offsets.GetIdx();
   return {sender_id, receiver_id, pcv->label(), location_idx};
 }
 
 inline std::tuple<int, int, std::string, int>
 ReceiveKey(const MeshBlock *pmb, const NeighborBlock &nb,
            const std::shared_ptr<Variable<Real>> &pcv) {
   const int receiver_id = pmb->gid;
-  const int sender_id = nb.snb.gid;
-  const int location_idx = (1 - nb.ni.ox1) + 3 * (1 - nb.ni.ox2 + 3 * (1 - nb.ni.ox3));
+  const int sender_id = nb.gid;
+  const int location_idx = nb.offsets.GetReverseIdx();
   return {sender_id, receiver_id, pcv->label(), location_idx};
 }
 

src/bvals/comms/tag_map.cpp

@@ -27,10 +27,10 @@ using namespace loops::shorthands;
 
 TagMap::rank_pair_t TagMap::MakeChannelPair(const MeshBlock *pmb,
                                             const NeighborBlock &nb) {
-  const int location_idx_me = (1 + nb.ni.ox1) + 3 * (1 + nb.ni.ox2 + 3 * (1 + nb.ni.ox3));
-  const int location_idx_nb = (1 - nb.ni.ox1) + 3 * (1 - nb.ni.ox2 + 3 * (1 - nb.ni.ox3));
+  const int location_idx_me = nb.offsets.GetIdx();
+  const int location_idx_nb = nb.offsets.GetReverseIdx();
   BlockGeometricElementId bgei_me{pmb->gid, location_idx_me};
-  BlockGeometricElementId bgei_nb{nb.snb.gid, location_idx_nb};
+  BlockGeometricElementId bgei_nb{nb.gid, location_idx_nb};
   return UnorderedPair<BlockGeometricElementId>(bgei_me, bgei_nb);
 }
 template <BoundaryType BOUND>
@@ -58,7 +58,7 @@ void TagMap::AddMeshDataToMap(std::shared_ptr<MeshData<Real>> &md) {
       return &(pmb->neighbors);
     }();
     for (auto &nb : *neighbors) {
-      const int other_rank = nb.snb.rank;
+      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];
       // Add channel key with an invalid tag
@@ -101,7 +101,7 @@ void TagMap::ResolveMap() {
 }
 
 int TagMap::GetTag(const MeshBlock *pmb, const NeighborBlock &nb) {
-  const int other_rank = nb.snb.rank;
+  const int other_rank = nb.rank;
   auto &pair_map = map_[other_rank];
   auto cpair = MakeChannelPair(pmb, nb);
   PARTHENON_REQUIRE(pair_map.count(cpair) == 1,