Fix refinement next to refined blocks

src/bvals/bvals.hpp

@@ -20,6 +20,7 @@
 //  \brief defines BoundaryBase, BoundaryValues classes used for setting BCs on all data
 
 #include <memory>
+#include <set>
 #include <string>
 #include <vector>
 
@@ -78,7 +79,8 @@ class BoundaryBase {
   static int BufferID(int dim, bool multilevel);
   static int FindBufferID(int ox1, int ox2, int ox3, int fi1, int fi2);
 
-  void SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist, int *nslist);
+  void SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist, int *nslist,
+                             const std::set<LogicalLocation> &newly_refined = {});
 
  protected:
   // 1D refined or unrefined=2
@@ -90,7 +92,7 @@ class BoundaryBase {
   RegionSize block_size_;
   ParArrayND<Real> sarea_[2];
 
-  void SetNeighborOwnership();
+  void SetNeighborOwnership(const std::set<LogicalLocation> &newly_refined = {});
 
  private:
   // calculate 3x shared static data members when constructing only the 1st class instance

src/bvals/bvals_base.cpp

@@ -300,8 +300,8 @@ int BoundaryBase::CreateBvalsMPITag(int lid, int bufid) {
 
 // TODO(felker): break-up this long function
 
-void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
-                                         int *nslist) {
+void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist, int *nslist,
+                                         const std::set<LogicalLocation> &newly_refined) {
   Kokkos::Profiling::pushRegion("SearchAndSetNeighbors");
   MeshBlockTree *neibt;
   int myox1, myox2 = 0, myox3 = 0, myfx1, myfx2, myfx3;
@@ -368,7 +368,7 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
     }
   }
   if (block_size_.nx(X2DIR) == 1) {
-    SetNeighborOwnership();
+    SetNeighborOwnership(newly_refined);
     Kokkos::Profiling::popRegion(); // SearchAndSetNeighbors
     return;
   }
@@ -503,7 +503,7 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
   }
 
   if (block_size_.nx(X3DIR) == 1) {
-    SetNeighborOwnership();
+    SetNeighborOwnership(newly_refined);
     Kokkos::Profiling::popRegion(); // SearchAndSetNeighbors
     return;
   }
@@ -626,11 +626,11 @@ void BoundaryBase::SearchAndSetNeighbors(MeshBlockTree &tree, int *ranklist,
     }
   }
 
-  SetNeighborOwnership();
+  SetNeighborOwnership(newly_refined);
   Kokkos::Profiling::popRegion(); // SearchAndSetNeighbors
 }
 
-void BoundaryBase::SetNeighborOwnership() {
+void BoundaryBase::SetNeighborOwnership(const std::set<LogicalLocation> &newly_refined) {
   // Set neighbor block ownership
   std::set<LogicalLocation> allowed_neighbors;
   allowed_neighbors.insert(loc); // Insert the location of this block
@@ -642,7 +642,7 @@ void BoundaryBase::SetNeighborOwnership() {
   RootGridInfo rg_info = pmy_mesh_->GetRootGridInfo();
   for (int n = 0; n < nneighbor; ++n) {
     neighbor[n].ownership =
-        DetermineOwnership(neighbor[n].loc, allowed_neighbors, rg_info);
+        DetermineOwnership(neighbor[n].loc, allowed_neighbors, rg_info, newly_refined);
     neighbor[n].ownership.initialized = true;
   }
 }

src/mesh/amr_loadbalance.cpp

@@ -789,6 +789,7 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput
   RegionSize block_size = GetBlockSize();
 
   BlockList_t new_block_list(nbe - nbs + 1);
+  std::set<LogicalLocation> newly_refined;
   for (int n = nbs; n <= nbe; n++) {
     int on = newtoold[n];
     if ((ranklist[on] == Globals::my_rank) &&
@@ -810,6 +811,7 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput
       new_block_list[n - nbs] =
           MeshBlock::Make(n, n - nbs, newloc[n], block_size, block_bcs, this, pin, app_in,
                           packages, resolved_packages, gflag);
+      if (newloc[n].level() > loclist[on].level()) newly_refined.insert(newloc[n]);
     }
   }
 
@@ -913,29 +915,47 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput
     }
   }
   prolongation_cache.CopyToDevice();
+
   refinement::ProlongateShared(resolved_packages.get(), prolongation_cache,
                                block_list[0]->cellbounds, block_list[0]->c_cellbounds);
-  refinement::ProlongateInternal(resolved_packages.get(), prolongation_cache,
-                                 block_list[0]->cellbounds, block_list[0]->c_cellbounds);
-
-#ifdef MPI_PARALLEL
-  if (send_reqs.size() != 0)
-    PARTHENON_MPI_CHECK(
-        MPI_Waitall(send_reqs.size(), send_reqs.data(), MPI_STATUSES_IGNORE));
-#endif
-  Kokkos::Profiling::popRegion(); // AMR: Recv data and unpack
 
   // update the lists
   loclist = std::move(newloc);
   ranklist = std::move(newrank);
   costlist = std::move(newcost);
 
-  // re-initialize the MeshBlocks
+  // A block newly refined and prolongated may have neighbors which were
+  // already refined to the new level.
+  // If so, the prolongated versions of shared elements will not reflect
+  // the true, finer versions present in the neighbor block.
+  // We must create any new fine buffers and fill them from these neighbors
+  // in order to maintain a consistent global state.
+  // Thus we rebuild and synchronize the mesh now, but using a unique
+  // neighbor precedence favoring the "old" fine blocks over "new" ones
   for (auto &pmb : block_list) {
-    pmb->pbval->SearchAndSetNeighbors(tree, ranklist.data(), nslist.data());
+    pmb->pbval->SearchAndSetNeighbors(tree, ranklist.data(), nslist.data(),
+                                      newly_refined);
   }
   Initialize(false, pin, app_in);
 
+  // Internal refinement relies on the fine shared values, which are only consistent after
+  // being updated with any previously fine versions
+  refinement::ProlongateInternal(resolved_packages.get(), prolongation_cache,
+                                 block_list[0]->cellbounds, block_list[0]->c_cellbounds);
+
+  // Rebuild the ownership model, this time weighting the "new" fine blocks just like
+  // any other blocks at their level.
+  for (auto &pmb : block_list) {
+    pmb->pbval->SearchAndSetNeighbors(tree, ranklist.data(), nslist.data());
+  }
+
+#ifdef MPI_PARALLEL
+  if (send_reqs.size() != 0)
+    PARTHENON_MPI_CHECK(
+        MPI_Waitall(send_reqs.size(), send_reqs.data(), MPI_STATUSES_IGNORE));
+#endif
+  Kokkos::Profiling::popRegion(); // AMR: Recv data and unpack
+
   ResetLoadBalanceVariables();
 
   Kokkos::Profiling::popRegion(); // RedistributeAndRefineMeshBlocks

src/mesh/logical_location.hpp

@@ -251,15 +251,25 @@ struct block_ownership_t {
 inline block_ownership_t
 DetermineOwnership(const LogicalLocation &main_block,
                    const std::set<LogicalLocation> &allowed_neighbors,
-                   const RootGridInfo &rg_info = RootGridInfo()) {
+                   const RootGridInfo &rg_info = RootGridInfo(),
+                   const std::set<LogicalLocation> &newly_refined = {}) {
   block_ownership_t main_owns;
 
-  auto ownership_less_than = [](const LogicalLocation &a, const LogicalLocation &b) {
+  auto ownership_level = [&](const LogicalLocation &a) {
+    // Newly-refined blocks are treated as higher-level than blocks at their
+    // parent level, but lower-level than previously-refined blocks at their
+    // current level.
+    if (newly_refined.count(a)) return 2 * a.level() - 1;
+    return 2 * a.level();
+  };
+
+  auto ownership_less_than = [ownership_level](const LogicalLocation &a,
+                                               const LogicalLocation &b) {
     // Ownership is first determined by block with the highest level, then by maximum
     // Morton number this is reversed in precedence from the normal comparators where
     // Morton number takes precedence
-    if (a.level() == b.level()) return a.morton() < b.morton();
-    return a.level() < b.level();
+    if (ownership_level(a) == ownership_level(b)) return a.morton() < b.morton();
+    return ownership_level(a) < ownership_level(b);
   };
 
   for (int ox1 : {-1, 0, 1}) {

src/mesh/mesh.cpp

@@ -475,11 +475,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages,
   }
 
   ResetLoadBalanceVariables();
-
-  // Output variables in use in this run
-  if (Globals::my_rank == 0) {
-    std::cout << "#Variables in use:\n" << *(resolved_packages) << std::endl;
-  }
 }
 
 //----------------------------------------------------------------------------------------
@@ -738,11 +733,6 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, RestartReader &rr,
   }
 
   ResetLoadBalanceVariables();
-
-  // Output variables in use in this run
-  if (Globals::my_rank == 0) {
-    std::cout << "#Variables in use:\n" << *(resolved_packages) << std::endl;
-  }
 }
 
 //----------------------------------------------------------------------------------------