Add Morton numbers and expand LogicalLocation functionality

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 900]](https://github.com/parthenon-hpc-lab/parthenon/pull/900) Add Morton numbers and expand functionality of LogicalLocation
 - [[PR 902]](https://github.com/parthenon-hpc-lab/parthenon/pull/902) Add ability to output NaNs for de-allocated sparse fields
 - [[PR 887]](https://github.com/parthenon-hpc-lab/parthenon/pull/887) Add ability to dump more types of params and read them from restarts
 - [[PR 884]](https://github.com/parthenon-hpc-lab/parthenon/pull/884) Add constant derivative BC and expose GenericBC

src/CMakeLists.txt

@@ -158,6 +158,7 @@ add_library(parthenon
 
   mesh/amr_loadbalance.cpp
   mesh/domain.hpp
+  mesh/logical_location.hpp
   mesh/mesh_refinement.cpp
   mesh/mesh_refinement.hpp
   mesh/mesh.cpp
@@ -229,6 +230,7 @@ add_library(parthenon
   utils/hash.hpp
   utils/indexer.hpp
   utils/loop_utils.hpp
+  utils/morton_number.hpp
   utils/mpi_types.hpp
   utils/multi_pointer.hpp
   utils/nan_payload_tag.hpp

src/amr_criteria/refinement_package.cpp

@@ -73,7 +73,7 @@ AmrTag CheckAllRefinement(MeshBlockData<Real> *rc) {
     for (auto &amr : desc->amr_criteria) {
       // get the recommended change in refinement level from this criteria
       AmrTag temp_delta = (*amr)(rc);
-      if ((temp_delta == AmrTag::refine) && pmb->loc.level >= amr->max_level) {
+      if ((temp_delta == AmrTag::refine) && pmb->loc.level() >= amr->max_level) {
         // don't refine if we're at the max level
         temp_delta = AmrTag::same;
       }

src/bvals/boundary_conditions.cpp

@@ -178,9 +178,9 @@ void ComputeProlongationBounds_(const std::shared_ptr<MeshBlock> &pmb,
     }
   };
 
-  getbounds(nb.ni.ox1, pmb->loc.lx1, pmb->c_cellbounds.GetBoundsI(interior), bi);
-  getbounds(nb.ni.ox2, pmb->loc.lx2, pmb->c_cellbounds.GetBoundsJ(interior), bj);
-  getbounds(nb.ni.ox3, pmb->loc.lx3, pmb->c_cellbounds.GetBoundsK(interior), bk);
+  getbounds(nb.ni.ox1, pmb->loc.lx1(), pmb->c_cellbounds.GetBoundsI(interior), bi);
+  getbounds(nb.ni.ox2, pmb->loc.lx2(), pmb->c_cellbounds.GetBoundsJ(interior), bj);
+  getbounds(nb.ni.ox3, pmb->loc.lx3(), pmb->c_cellbounds.GetBoundsK(interior), bk);
 }
 
 } // namespace boundary_cond_impl

src/bvals/bvals_base.cpp

@@ -51,9 +51,9 @@ NeighborIndexes BoundaryBase::ni[NMAX_NEIGHBORS];
 //                          int ibid, int itargetid, int ifi1=0, int ifi2=0)
 // \brief Set neighbor information
 
-void NeighborBlock::SetNeighbor(int irank, int ilevel, int igid, int ilid, int iox1,
-                                int iox2, int iox3, NeighborConnect itype, int ibid,
-                                int itargetid,
+void NeighborBlock::SetNeighbor(LogicalLocation inloc, int irank, int ilevel, int igid,
+                                int ilid, int iox1, int iox2, int iox3,
+                                NeighborConnect itype, int ibid, int itargetid,
                                 int ifi1, // =0
                                 int ifi2  // =0
 ) {
@@ -69,6 +69,7 @@ void NeighborBlock::SetNeighbor(int irank, int ilevel, int igid, int ilid, int i
   ni.fi2 = ifi2;
   bufid = ibid;
   targetid = itargetid;
+  loc = inloc;
   if (ni.type == NeighborConnect::face) {
     if (ni.ox1 == -1)
       fid = BoundaryFace::inner_x1;
@@ -304,12 +305,12 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
   Kokkos::Profiling::pushRegion("SearchAndSetNeighbors");
   MeshBlockTree *neibt;
   int myox1, myox2 = 0, myox3 = 0, myfx1, myfx2, myfx3;
-  myfx1 = ((loc.lx1 & 1LL) == 1LL);
-  myfx2 = ((loc.lx2 & 1LL) == 1LL);
-  myfx3 = ((loc.lx3 & 1LL) == 1LL);
-  myox1 = ((loc.lx1 & 1LL) == 1LL) * 2 - 1;
-  if (block_size_.nx2 > 1) myox2 = ((loc.lx2 & 1LL) == 1LL) * 2 - 1;
-  if (block_size_.nx3 > 1) myox3 = ((loc.lx3 & 1LL) == 1LL) * 2 - 1;
+  myfx1 = ((loc.lx1() & 1LL) == 1LL);
+  myfx2 = ((loc.lx2() & 1LL) == 1LL);
+  myfx3 = ((loc.lx3() & 1LL) == 1LL);
+  myox1 = ((loc.lx1() & 1LL) == 1LL) * 2 - 1;
+  if (block_size_.nx2 > 1) myox2 = ((loc.lx2() & 1LL) == 1LL) * 2 - 1;
+  if (block_size_.nx3 > 1) myox3 = ((loc.lx3() & 1LL) == 1LL) * 2 - 1;
 
   int nf1 = 1, nf2 = 1;
   if (pmy_mesh_->multilevel) {
@@ -324,7 +325,7 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
         nblevel[k][j][i] = -1;
     }
   }
-  nblevel[1][1][1] = loc.level;
+  nblevel[1][1][1] = loc.level();
 
   // x1 face
   for (int n = -1; n <= 1; n += 2) {
@@ -335,31 +336,31 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
     }
     if (neibt->pleaf_ != nullptr) { // neighbor at finer level
       int fface = 1 - (n + 1) / 2;  // 0 for BoundaryFace::outer_x1, 1 for inner_x1
-      nblevel[1][1][n + 1] = neibt->loc_.level + 1;
+      nblevel[1][1][n + 1] = neibt->loc_.level() + 1;
       for (int f2 = 0; f2 < nf2; f2++) {
         for (int f1 = 0; f1 < nf1; f1++) {
           MeshBlockTree *nf = neibt->GetLeaf(fface, f1, f2);
           int fid = nf->gid_;
-          int nlevel = nf->loc_.level;
+          int nlevel = nf->loc_.level();
           int tbid = FindBufferID(-n, 0, 0, 0, 0);
-          neighbor[nneighbor].SetNeighbor(ranklist[fid], nlevel, fid,
+          neighbor[nneighbor].SetNeighbor(nf->loc_, ranklist[fid], nlevel, fid,
                                           fid - nslist[ranklist[fid]], n, 0, 0,
                                           NeighborConnect::face, bufid, tbid, f1, f2);
           bufid++;
           nneighbor++;
         }
       }
     } else { // neighbor at same or coarser level
-      int nlevel = neibt->loc_.level;
+      int nlevel = neibt->loc_.level();
       int nid = neibt->gid_;
       nblevel[1][1][n + 1] = nlevel;
       int tbid;
-      if (nlevel == loc.level) { // neighbor at same level
+      if (nlevel == loc.level()) { // neighbor at same level
         tbid = FindBufferID(-n, 0, 0, 0, 0);
       } else { // neighbor at coarser level
         tbid = FindBufferID(-n, 0, 0, myfx2, myfx3);
       }
-      neighbor[nneighbor].SetNeighbor(ranklist[nid], nlevel, nid,
+      neighbor[nneighbor].SetNeighbor(neibt->loc_, ranklist[nid], nlevel, nid,
                                       nid - nslist[ranklist[nid]], n, 0, 0,
                                       NeighborConnect::face, bufid, tbid);
       bufid += nf1 * nf2;
@@ -380,31 +381,31 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
     }
     if (neibt->pleaf_ != nullptr) { // neighbor at finer level
       int fface = 1 - (n + 1) / 2;  // 0 for BoundaryFace::outer_x2, 1 for inner_x2
-      nblevel[1][n + 1][1] = neibt->loc_.level + 1;
+      nblevel[1][n + 1][1] = neibt->loc_.level() + 1;
       for (int f2 = 0; f2 < nf2; f2++) {
         for (int f1 = 0; f1 < nf1; f1++) {
           MeshBlockTree *nf = neibt->GetLeaf(f1, fface, f2);
           int fid = nf->gid_;
-          int nlevel = nf->loc_.level;
+          int nlevel = nf->loc_.level();
           int tbid = FindBufferID(0, -n, 0, 0, 0);
-          neighbor[nneighbor].SetNeighbor(ranklist[fid], nlevel, fid,
+          neighbor[nneighbor].SetNeighbor(nf->loc_, ranklist[fid], nlevel, fid,
                                           fid - nslist[ranklist[fid]], 0, n, 0,
                                           NeighborConnect::face, bufid, tbid, f1, f2);
           bufid++;
           nneighbor++;
         }
       }
     } else { // neighbor at same or coarser level
-      int nlevel = neibt->loc_.level;
+      int nlevel = neibt->loc_.level();
       int nid = neibt->gid_;
       nblevel[1][n + 1][1] = nlevel;
       int tbid;
-      if (nlevel == loc.level) { // neighbor at same level
+      if (nlevel == loc.level()) { // neighbor at same level
         tbid = FindBufferID(0, -n, 0, 0, 0);
       } else { // neighbor at coarser level
         tbid = FindBufferID(0, -n, 0, myfx1, myfx3);
       }
-      neighbor[nneighbor].SetNeighbor(ranklist[nid], nlevel, nid,
+      neighbor[nneighbor].SetNeighbor(neibt->loc_, ranklist[nid], nlevel, nid,
                                       nid - nslist[ranklist[nid]], 0, n, 0,
                                       NeighborConnect::face, bufid, tbid);
       bufid += nf1 * nf2;
@@ -422,31 +423,31 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
       }
       if (neibt->pleaf_ != nullptr) { // neighbor at finer level
         int fface = 1 - (n + 1) / 2;  // 0 for BoundaryFace::outer_x3, 1 for inner_x3
-        nblevel[n + 1][1][1] = neibt->loc_.level + 1;
+        nblevel[n + 1][1][1] = neibt->loc_.level() + 1;
         for (int f2 = 0; f2 < nf2; f2++) {
           for (int f1 = 0; f1 < nf1; f1++) {
             MeshBlockTree *nf = neibt->GetLeaf(f1, f2, fface);
             int fid = nf->gid_;
-            int nlevel = nf->loc_.level;
+            int nlevel = nf->loc_.level();
             int tbid = FindBufferID(0, 0, -n, 0, 0);
-            neighbor[nneighbor].SetNeighbor(ranklist[fid], nlevel, fid,
+            neighbor[nneighbor].SetNeighbor(nf->loc_, ranklist[fid], nlevel, fid,
                                             fid - nslist[ranklist[fid]], 0, 0, n,
                                             NeighborConnect::face, bufid, tbid, f1, f2);
             bufid++;
             nneighbor++;
           }
         }
       } else { // neighbor at same or coarser level
-        int nlevel = neibt->loc_.level;
+        int nlevel = neibt->loc_.level();
         int nid = neibt->gid_;
         nblevel[n + 1][1][1] = nlevel;
         int tbid;
-        if (nlevel == loc.level) { // neighbor at same level
+        if (nlevel == loc.level()) { // neighbor at same level
           tbid = FindBufferID(0, 0, -n, 0, 0);
         } else { // neighbor at coarser level
           tbid = FindBufferID(0, 0, -n, myfx1, myfx2);
         }
-        neighbor[nneighbor].SetNeighbor(ranklist[nid], nlevel, nid,
+        neighbor[nneighbor].SetNeighbor(neibt->loc_, ranklist[nid], nlevel, nid,
                                         nid - nslist[ranklist[nid]], 0, 0, n,
                                         NeighborConnect::face, bufid, tbid);
         bufid += nf1 * nf2;
@@ -467,30 +468,30 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
         int ff1 = 1 - (n + 1) / 2;    // 0 for BoundaryFace::outer_x1, 1 for inner_x1
         int ff2 = 1 - (m + 1) / 2;    // 0 for BoundaryFace::outer_x2, 1 for inner_x2
 
-        nblevel[1][m + 1][n + 1] = neibt->loc_.level + 1;
+        nblevel[1][m + 1][n + 1] = neibt->loc_.level() + 1;
         for (int f1 = 0; f1 < nf2; f1++) {
           MeshBlockTree *nf = neibt->GetLeaf(ff1, ff2, f1);
           int fid = nf->gid_;
-          int nlevel = nf->loc_.level;
+          int nlevel = nf->loc_.level();
           int tbid = FindBufferID(-n, -m, 0, 0, 0);
-          neighbor[nneighbor].SetNeighbor(ranklist[fid], nlevel, fid,
+          neighbor[nneighbor].SetNeighbor(nf->loc_, ranklist[fid], nlevel, fid,
                                           fid - nslist[ranklist[fid]], n, m, 0,
                                           NeighborConnect::edge, bufid, tbid, f1, 0);
           bufid++;
           nneighbor++;
         }
       } else { // neighbor at same or coarser level
-        int nlevel = neibt->loc_.level;
+        int nlevel = neibt->loc_.level();
         int nid = neibt->gid_;
         nblevel[1][m + 1][n + 1] = nlevel;
         int tbid;
-        if (nlevel == loc.level) { // neighbor at same level
+        if (nlevel == loc.level()) { // neighbor at same level
           tbid = FindBufferID(-n, -m, 0, 0, 0);
         } else { // neighbor at coarser level
           tbid = FindBufferID(-n, -m, 0, myfx3, 0);
         }
-        if (nlevel >= loc.level || (myox1 == n && myox2 == m)) {
-          neighbor[nneighbor].SetNeighbor(ranklist[nid], nlevel, nid,
+        if (nlevel >= loc.level() || (myox1 == n && myox2 == m)) {
+          neighbor[nneighbor].SetNeighbor(neibt->loc_, ranklist[nid], nlevel, nid,
                                           nid - nslist[ranklist[nid]], n, m, 0,
                                           NeighborConnect::edge, bufid, tbid);
           nneighbor++;
@@ -516,30 +517,30 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
       if (neibt->pleaf_ != nullptr) { // neighbor at finer level
         int ff1 = 1 - (n + 1) / 2;    // 0 for BoundaryFace::outer_x1, 1 for inner_x1
         int ff2 = 1 - (m + 1) / 2;    // 0 for BoundaryFace::outer_x3, 1 for inner_x3
-        nblevel[m + 1][1][n + 1] = neibt->loc_.level + 1;
+        nblevel[m + 1][1][n + 1] = neibt->loc_.level() + 1;
         for (int f1 = 0; f1 < nf1; f1++) {
           MeshBlockTree *nf = neibt->GetLeaf(ff1, f1, ff2);
           int fid = nf->gid_;
-          int nlevel = nf->loc_.level;
+          int nlevel = nf->loc_.level();
           int tbid = FindBufferID(-n, 0, -m, 0, 0);
-          neighbor[nneighbor].SetNeighbor(ranklist[fid], nlevel, fid,
+          neighbor[nneighbor].SetNeighbor(nf->loc_, ranklist[fid], nlevel, fid,
                                           fid - nslist[ranklist[fid]], n, 0, m,
                                           NeighborConnect::edge, bufid, tbid, f1, 0);
           bufid++;
           nneighbor++;
         }
       } else { // neighbor at same or coarser level
-        int nlevel = neibt->loc_.level;
+        int nlevel = neibt->loc_.level();
         int nid = neibt->gid_;
         nblevel[m + 1][1][n + 1] = nlevel;
         int tbid;
-        if (nlevel == loc.level) { // neighbor at same level
+        if (nlevel == loc.level()) { // neighbor at same level
           tbid = FindBufferID(-n, 0, -m, 0, 0);
         } else { // neighbor at coarser level
           tbid = FindBufferID(-n, 0, -m, myfx2, 0);
         }
-        if (nlevel >= loc.level || (myox1 == n && myox3 == m)) {
-          neighbor[nneighbor].SetNeighbor(ranklist[nid], nlevel, nid,
+        if (nlevel >= loc.level() || (myox1 == n && myox3 == m)) {
+          neighbor[nneighbor].SetNeighbor(neibt->loc_, ranklist[nid], nlevel, nid,
                                           nid - nslist[ranklist[nid]], n, 0, m,
                                           NeighborConnect::edge, bufid, tbid);
           nneighbor++;
@@ -560,30 +561,30 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
       if (neibt->pleaf_ != nullptr) { // neighbor at finer level
         int ff1 = 1 - (n + 1) / 2;    // 0 for BoundaryFace::outer_x2, 1 for inner_x2
         int ff2 = 1 - (m + 1) / 2;    // 0 for BoundaryFace::outer_x3, 1 for inner_x3
-        nblevel[m + 1][n + 1][1] = neibt->loc_.level + 1;
+        nblevel[m + 1][n + 1][1] = neibt->loc_.level() + 1;
         for (int f1 = 0; f1 < nf1; f1++) {
           MeshBlockTree *nf = neibt->GetLeaf(f1, ff1, ff2);
           int fid = nf->gid_;
-          int nlevel = nf->loc_.level;
+          int nlevel = nf->loc_.level();
           int tbid = FindBufferID(0, -n, -m, 0, 0);
-          neighbor[nneighbor].SetNeighbor(ranklist[fid], nlevel, fid,
+          neighbor[nneighbor].SetNeighbor(nf->loc_, ranklist[fid], nlevel, fid,
                                           fid - nslist[ranklist[fid]], 0, n, m,
                                           NeighborConnect::edge, bufid, tbid, f1, 0);
           bufid++;
           nneighbor++;
         }
       } else { // neighbor at same or coarser level
-        int nlevel = neibt->loc_.level;
+        int nlevel = neibt->loc_.level();
         int nid = neibt->gid_;
         nblevel[m + 1][n + 1][1] = nlevel;
         int tbid;
-        if (nlevel == loc.level) { // neighbor at same level
+        if (nlevel == loc.level()) { // neighbor at same level
           tbid = FindBufferID(0, -n, -m, 0, 0);
         } else { // neighbor at coarser level
           tbid = FindBufferID(0, -n, -m, myfx1, 0);
         }
-        if (nlevel >= loc.level || (myox2 == n && myox3 == m)) {
-          neighbor[nneighbor].SetNeighbor(ranklist[nid], nlevel, nid,
+        if (nlevel >= loc.level() || (myox2 == n && myox3 == m)) {
+          neighbor[nneighbor].SetNeighbor(neibt->loc_, ranklist[nid], nlevel, nid,
                                           nid - nslist[ranklist[nid]], 0, n, m,
                                           NeighborConnect::edge, bufid, tbid);
           nneighbor++;
@@ -608,12 +609,12 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
           int ff3 = 1 - (l + 1) / 2;    // 0 for BoundaryFace::outer_x3, 1 for inner_x3
           neibt = neibt->GetLeaf(ff1, ff2, ff3);
         }
-        int nlevel = neibt->loc_.level;
+        int nlevel = neibt->loc_.level();
         nblevel[l + 1][m + 1][n + 1] = nlevel;
-        if (nlevel >= loc.level || (myox1 == n && myox2 == m && myox3 == l)) {
+        if (nlevel >= loc.level() || (myox1 == n && myox2 == m && myox3 == l)) {
           int nid = neibt->gid_;
           int tbid = FindBufferID(-n, -m, -l, 0, 0);
-          neighbor[nneighbor].SetNeighbor(ranklist[nid], nlevel, nid,
+          neighbor[nneighbor].SetNeighbor(neibt->loc_, ranklist[nid], nlevel, nid,
                                           nid - nslist[ranklist[nid]], n, m, l,
                                           NeighborConnect::corner, bufid, tbid);
           nneighbor++;

src/bvals/bvals_interfaces.hpp

@@ -142,10 +142,10 @@ struct NeighborBlock { // aggregate and POD type. Inheritance breaks standard-la
 
   int bufid, eid, targetid;
   BoundaryFace fid;
-
-  void SetNeighbor(int irank, int ilevel, int igid, int ilid, int iox1, int iox2,
-                   int iox3, NeighborConnect itype, int ibid, int itargetid, int ifi1 = 0,
-                   int ifi2 = 0);
+  LogicalLocation loc;
+  void SetNeighbor(LogicalLocation inloc, int irank, int ilevel, int igid, int ilid,
+                   int iox1, int iox2, int iox3, NeighborConnect itype, int ibid,
+                   int itargetid, int ifi1 = 0, int ifi2 = 0);
 };
 
 //----------------------------------------------------------------------------------------

src/bvals/bvals_swarm.cpp

@@ -127,7 +127,7 @@ void BoundarySwarm::Send(BoundaryCommSubset phase) {
 void BoundarySwarm::Receive(BoundaryCommSubset phase) {
 #ifdef MPI_PARALLEL
   std::shared_ptr<MeshBlock> pmb = GetBlockPointer();
-  int &mylevel = pmb->loc.level;
+  const int &mylevel = pmb->loc.level();
   for (int n = 0; n < pmb->pbval->nneighbor; n++) {
     NeighborBlock &nb = pmb->pbval->neighbor[n];
     if (nb.snb.rank != Globals::my_rank) {

src/bvals/comms/bnd_info.cpp

@@ -106,7 +106,7 @@ Indexer6D CalcSetIndices(const NeighborIndexes &ni, LogicalLocation loc,
   // information for going from coarse to fine and the neighbor block
   // structure does not contain the logical location of the neighbor
   // block
-  std::array<std::int64_t, 3> logic_loc{loc.lx1, loc.lx2, loc.lx3};
+  std::array<std::int64_t, 3> logic_loc{loc.lx1(), loc.lx2(), loc.lx3()};
   std::array<int, 2> face_offset{ni.fi1, ni.fi2};
   std::array<int, 3> s, e;
 
@@ -172,7 +172,7 @@ BndInfo BndInfo::GetSendBndInfo(std::shared_ptr<MeshBlock> pmb, const NeighborBl
   int Nu = v->GetDim(5);
   int Nt = v->GetDim(6);
 
-  int mylevel = pmb->loc.level;
+  int mylevel = pmb->loc.level();
   out.coords = pmb->coords;
 
   if (pmb->pmr) out.coarse_coords = pmb->pmr->GetCoarseCoords();
@@ -236,7 +236,7 @@ BndInfo BndInfo::GetSetBndInfo(std::shared_ptr<MeshBlock> pmb, const NeighborBlo
   if (v->IsSet(Metadata::Edge)) elements = {TE::E1, TE::E2, TE::E3};
   if (v->IsSet(Metadata::Node)) elements = {TE::NN};
 
-  int mylevel = pmb->loc.level;
+  int mylevel = pmb->loc.level();
   out.coords = pmb->coords;
   if (pmb->pmr) out.coarse_coords = pmb->pmr->GetCoarseCoords();
   out.fine = v->data.Get();