Merge branch 'lroberts36/add-multi-grid' into lroberts36/multigrid-ex…

…ample-update

CHANGELOG.md

@@ -26,14 +26,17 @@
 - [[PR 885]](https://github.com/parthenon-hpc-lab/parthenon/pull/885) Expose PackDescriptor and use uids in SparsePacks
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 947]](https://github.com/parthenon-hpc-lab/parthenon/pull/947) Add missing ForceRemeshComm dependencies
+- [[PR 928]](https://github.com/parthenon-hpc-lab/parthenon/pull/928) Fix boundary comms during refinement next to refined blocks
 - [[PR 937]](https://github.com/parthenon-hpc-lab/parthenon/pull/937) Fix multiple line continuations
 - [[PR 933]](https://github.com/parthenon-hpc-lab/parthenon/pull/933) Remove extraneous debug check
 - [[PR 917]](https://github.com/parthenon-hpc-lab/parthenon/pull/917) Update Iterative Tasking Infrastructure
 - [[PR 890]](https://github.com/parthenon-hpc-lab/parthenon/pull/890) Fix bugs in sparse communication and prolongation
 
 ### Infrastructure (changes irrelevant to downstream codes)
+- [[PR 944]](https://github.com/parthenon-hpc-lab/parthenon/pull/944) Move sparse pack identifier creation to descriptor
 - [[PR 904]](https://github.com/parthenon-hpc-lab/parthenon/pull/904) Move to prolongation/restriction in one for AMR and communicate non-cell centered fields
-- [[PR 918]](https://github.com/parthenon-hpc-lab/parthenon/pull/918) Refactor RegionSize 
+- [[PR 918]](https://github.com/parthenon-hpc-lab/parthenon/pull/918) Refactor RegionSize
 - [[PR 901]](https://github.com/parthenon-hpc-lab/parthenon/pull/901) Implement shared element ownership model
 
 ### Removed (removing behavior/API/varaibles/...)

doc/sphinx/src/mesh/mesh.rst

@@ -66,7 +66,7 @@ level that are neighbors to finer blocks, which implies that below the root grid
 level the blocks may not cover the entire mesh. For levels above the root grid, 
 blocks may change shape so that they only cover the domain of the root grid. Note 
 that leaf blocks may be contained in multiple blocklists, and the lists all point
-to the same block (not a separate copy).
+to the same block (not a separate copy). To be explicit, when ``parthenon/mesh/multigrid`` is set to ``true`` blocks corresponding to *all* internal nodes of the refinement tree are created, in addition to the leaf node blocks that are normally created.
 
 *GMG Implementation Note:*
 The reason for including two levels in the GMG block lists is for dealing with 
@@ -85,7 +85,7 @@ To work with these GMG levels, ``MeshData`` objects containing these blocks can
 be recovered from a ``Mesh`` pointer using 
 
 .. code:: c++
-  auto &md = pmesh->gmg_mesh_data[level].GetOrAdd(level, "base", i);
+  auto &md = pmesh->gmg_mesh_data[level].GetOrAdd(level, "base", partition_idx);
 
 This ``MeshData`` will include blocks at the current level and possibly some 
 blocks at the next coarser level. Often, one will only want to operate on blocks

src/bvals/bvals.hpp

@@ -21,6 +21,7 @@
 
 #include <memory>
 #include <string>
+#include <unordered_set>
 #include <vector>
 
 #include "basic_types.hpp"
@@ -78,7 +79,9 @@ class BoundaryBase {
   static int BufferID(int dim, bool multilevel);
   static int FindBufferID(int ox1, int ox2, int ox3, int fi1, int fi2);
 
-  void SearchAndSetNeighbors(Mesh *mesh, MeshBlockTree &tree, int *ranklist, int *nslist);
+  void
+  SearchAndSetNeighbors(Mesh *mesh, MeshBlockTree &tree, int *ranklist, int *nslist,
+                        const std::unordered_set<LogicalLocation> &newly_refined = {});
 
  protected:
   // 1D refined or unrefined=2
@@ -90,7 +93,8 @@ class BoundaryBase {
   RegionSize block_size_;
   ParArrayND<Real> sarea_[2];
 
-  void SetNeighborOwnership();
+  void
+  SetNeighborOwnership(const std::unordered_set<LogicalLocation> &newly_refined = {});
 
  private:
   // calculate 3x shared static data members when constructing only the 1st class instance

src/bvals/bvals_base.cpp

@@ -30,6 +30,7 @@
 #include <unordered_set>
 
 #include "globals.hpp"
+#include "mesh/logical_location.hpp"
 #include "mesh/mesh.hpp"
 #include "utils/buffer_utils.hpp"
 #include "utils/error_checking.hpp"
@@ -333,8 +334,9 @@ int BoundaryBase::CreateBvalsMPITag(int lid, int bufid) {
 
 // TODO(felker): break-up this long function
 
-void BoundaryBase::SearchAndSetNeighbors(Mesh *mesh, MeshBlockTree &tree, int *ranklist,
-                                         int *nslist) {
+void BoundaryBase::SearchAndSetNeighbors(
+    Mesh *mesh, MeshBlockTree &tree, int *ranklist, int *nslist,
+    const std::unordered_set<LogicalLocation> &newly_refined) {
   Kokkos::Profiling::pushRegion("SearchAndSetNeighbors");
   MeshBlockTree *neibt;
   int myox1, myox2 = 0, myox3 = 0, myfx1, myfx2, myfx3;
@@ -401,7 +403,7 @@ void BoundaryBase::SearchAndSetNeighbors(Mesh *mesh, MeshBlockTree &tree, int *r
     }
   }
   if (block_size_.nx(X2DIR) == 1) {
-    SetNeighborOwnership();
+    SetNeighborOwnership(newly_refined);
     Kokkos::Profiling::popRegion(); // SearchAndSetNeighbors
     return;
   }
@@ -536,7 +538,7 @@ void BoundaryBase::SearchAndSetNeighbors(Mesh *mesh, MeshBlockTree &tree, int *r
   }
 
   if (block_size_.nx(X3DIR) == 1) {
-    SetNeighborOwnership();
+    SetNeighborOwnership(newly_refined);
     Kokkos::Profiling::popRegion(); // SearchAndSetNeighbors
     return;
   }
@@ -659,11 +661,12 @@ void BoundaryBase::SearchAndSetNeighbors(Mesh *mesh, MeshBlockTree &tree, int *r
     }
   }
 
-  SetNeighborOwnership();
+  SetNeighborOwnership(newly_refined);
   Kokkos::Profiling::popRegion(); // SearchAndSetNeighbors
 }
 
-void BoundaryBase::SetNeighborOwnership() {
+void BoundaryBase::SetNeighborOwnership(
+    const std::unordered_set<LogicalLocation> &newly_refined) {
   // Set neighbor block ownership
   std::unordered_set<LogicalLocation> allowed_neighbors;
   allowed_neighbors.insert(loc); // Insert the location of this block
@@ -675,7 +678,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/bvals/comms/boundary_communication.cpp

@@ -155,7 +155,7 @@ template <BoundaryType bound_type>
 TaskStatus StartReceiveBoundBufs(std::shared_ptr<MeshData<Real>> &md) {
   Kokkos::Profiling::pushRegion("Task_StartReceiveBoundBufs");
   Mesh *pmesh = md->GetMeshPointer();
-  auto &cache = md->GetBvarsCache().GetSubCache(BoundaryType::flxcor_send, false);
+  auto &cache = md->GetBvarsCache().GetSubCache(bound_type, false);
   if (cache.buf_vec.size() == 0)
     InitializeBufferCache<bound_type>(md, &(pmesh->boundary_comm_map), &cache, ReceiveKey,
                                       false);

src/interface/sparse_pack_base.cpp

@@ -287,21 +287,28 @@ template SparsePackBase SparsePackBase::Build<MeshData<Real>>(MeshData<Real> *,
 template <class T>
 SparsePackBase &SparsePackCache::Get(T *pmd, const PackDescriptor &desc,
                                      const std::vector<bool> &include_block) {
-  std::string ident = GetIdentifier(desc, include_block);
-  if (pack_map.count(ident) > 0) {
-    auto &pack = pack_map[ident].first;
+  if (pack_map.count(desc.identifier) > 0) {
+    auto &cache_tuple = pack_map[desc.identifier];
+    auto &pack = std::get<0>(cache_tuple);
     auto alloc_status_in = SparsePackBase::GetAllocStatus(pmd, desc, include_block);
-    auto &alloc_status = pack_map[ident].second;
+    auto &alloc_status = std::get<1>(cache_tuple);
     if (alloc_status.size() != alloc_status_in.size())
-      return BuildAndAdd(pmd, desc, ident, include_block);
+      return BuildAndAdd(pmd, desc, include_block);
     for (int i = 0; i < alloc_status_in.size(); ++i) {
       if (alloc_status[i] != alloc_status_in[i])
-        return BuildAndAdd(pmd, desc, ident, include_block);
+        return BuildAndAdd(pmd, desc, include_block);
+    }
+    auto &include_status = std::get<2>(cache_tuple);
+    if (include_status.size() != include_block.size())
+      return BuildAndAdd(pmd, desc, include_block);
+    for (int i = 0; i < include_block.size(); ++i) {
+      if (include_status[i] != include_block[i])
+        return BuildAndAdd(pmd, desc, include_block);
     }
     // Cached version is not stale, so just return a reference to it
-    return pack_map[ident].first;
+    return std::get<0>(cache_tuple);
   }
-  return BuildAndAdd(pmd, desc, ident, include_block);
+  return BuildAndAdd(pmd, desc, include_block);
 }
 template SparsePackBase &SparsePackCache::Get<MeshData<Real>>(MeshData<Real> *,
                                                               const PackDescriptor &,
@@ -312,37 +319,17 @@ SparsePackCache::Get<MeshBlockData<Real>>(MeshBlockData<Real> *, const PackDescr
 
 template <class T>
 SparsePackBase &SparsePackCache::BuildAndAdd(T *pmd, const PackDescriptor &desc,
-                                             const std::string &ident,
                                              const std::vector<bool> &include_block) {
-  if (pack_map.count(ident) > 0) pack_map.erase(ident);
-  pack_map[ident] = {SparsePackBase::Build(pmd, desc, include_block),
-                     SparsePackBase::GetAllocStatus(pmd, desc, include_block)};
-  return pack_map[ident].first;
+  if (pack_map.count(desc.identifier) > 0) pack_map.erase(desc.identifier);
+  pack_map[desc.identifier] = {SparsePackBase::Build(pmd, desc, include_block),
+                               SparsePackBase::GetAllocStatus(pmd, desc, include_block),
+                               include_block};
+  return std::get<0>(pack_map[desc.identifier]);
 }
 template SparsePackBase &
 SparsePackCache::BuildAndAdd<MeshData<Real>>(MeshData<Real> *, const PackDescriptor &,
-                                             const std::string &,
                                              const std::vector<bool> &);
 template SparsePackBase &SparsePackCache::BuildAndAdd<MeshBlockData<Real>>(
-    MeshBlockData<Real> *, const PackDescriptor &, const std::string &,
-    const std::vector<bool> &);
-
-std::string SparsePackCache::GetIdentifier(const PackDescriptor &desc,
-                                           const std::vector<bool> &include_block) const {
-  std::string identifier("");
-  for (const auto &vgroup : desc.var_groups) {
-    for (const auto &[vid, uid] : vgroup) {
-      identifier += std::to_string(uid) + "_";
-    }
-    identifier += "|";
-  }
-  identifier += std::to_string(desc.with_fluxes);
-  identifier += std::to_string(desc.coarse);
-  identifier += std::to_string(desc.flat);
-  for (const auto b : include_block) {
-    identifier += std::to_string(b);
-  }
-  return identifier;
-}
+    MeshBlockData<Real> *, const PackDescriptor &, const std::vector<bool> &);
 
 } // namespace parthenon

src/interface/sparse_pack_base.hpp

@@ -54,6 +54,7 @@ class SparsePackBase {
   friend class SparsePackCache;
 
   using alloc_t = std::vector<int>;
+  using include_t = std::vector<bool>;
   using pack_t = ParArray3D<ParArray3D<Real, VariableState>>;
   using bounds_t = ParArray3D<int>;
   using bounds_h_t = typename ParArray3D<int>::HostMirror;
@@ -114,12 +115,10 @@ class SparsePackCache {
 
   template <class T>
   SparsePackBase &BuildAndAdd(T *pmd, const impl::PackDescriptor &desc,
-                              const std::string &ident,
                               const std::vector<bool> &include_block);
 
-  std::string GetIdentifier(const impl::PackDescriptor &desc,
-                            const std::vector<bool> &include_block) const;
-  std::unordered_map<std::string, std::pair<SparsePackBase, SparsePackBase::alloc_t>>
+  std::unordered_map<std::string, std::tuple<SparsePackBase, SparsePackBase::alloc_t,
+                                             SparsePackBase::include_t>>
       pack_map;
 
   friend class SparsePackBase;
@@ -136,15 +135,16 @@ struct PackDescriptor {
   // default constructor needed for certain use cases
   PackDescriptor()
       : nvar_groups(0), var_group_names({}), var_groups({}), with_fluxes(false),
-        coarse(false), flat(false) {}
+        coarse(false), flat(false), identifier("") {}
 
   template <class GROUP_t, class SELECTOR_t>
   PackDescriptor(StateDescriptor *psd, const std::vector<GROUP_t> &var_groups_in,
                  const SELECTOR_t &selector, const std::set<PDOpt> &options)
       : nvar_groups(var_groups_in.size()), var_group_names(MakeGroupNames(var_groups_in)),
         var_groups(BuildUids(var_groups_in.size(), psd, selector)),
         with_fluxes(options.count(PDOpt::WithFluxes)),
-        coarse(options.count(PDOpt::Coarse)), flat(options.count(PDOpt::Flatten)) {
+        coarse(options.count(PDOpt::Coarse)), flat(options.count(PDOpt::Flatten)),
+        identifier(GetIdentifier()) {
     PARTHENON_REQUIRE(!(with_fluxes && coarse),
                       "Probably shouldn't be making a coarse pack with fine fluxes.");
   }
@@ -155,8 +155,22 @@ struct PackDescriptor {
   const bool with_fluxes;
   const bool coarse;
   const bool flat;
+  const std::string identifier;
 
  private:
+  std::string GetIdentifier() {
+    std::string ident("");
+    for (const auto &vgroup : var_groups) {
+      for (const auto &[vid, uid] : vgroup) {
+        ident += std::to_string(uid) + "_";
+      }
+      ident += "|";
+    }
+    ident += std::to_string(with_fluxes);
+    ident += std::to_string(coarse);
+    ident += std::to_string(flat);
+    return ident;
+  }
   template <class FUNC_t>
   std::vector<PackDescriptor::VariableGroup_t>
   BuildUids(int nvgs, const StateDescriptor *const psd, const FUNC_t &selector) {

src/interface/variable.cpp

@@ -1,5 +1,5 @@
 //========================================================================================
-// (C) (or copyright) 2020-2022. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved.
 //
 // This program was produced under U.S. Government contract 89233218CNA000001 for Los
 // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -81,7 +81,8 @@ void Variable<T>::CopyFluxesAndBdryVar(const Variable<T> *src) {
     }
   }
 
-  if (IsSet(Metadata::FillGhost) || IsSet(Metadata::Independent)) {
+  if (IsSet(Metadata::FillGhost) || IsSet(Metadata::Independent) ||
+      IsSet(Metadata::ForceRemeshComm)) {
     // no need to check mesh->multilevel, if false, we're just making a shallow copy of
     // an empty ParArrayND
     coarse_s = src->coarse_s;
@@ -172,7 +173,8 @@ void Variable<T>::AllocateFluxesAndCoarse(std::weak_ptr<MeshBlock> wpmb) {
   }
 
   // Create the boundary object
-  if (IsSet(Metadata::FillGhost) || IsSet(Metadata::Independent)) {
+  if (IsSet(Metadata::FillGhost) || IsSet(Metadata::Independent) ||
+      IsSet(Metadata::ForceRemeshComm)) {
     if (wpmb.expired()) return;
     std::shared_ptr<MeshBlock> pmb = wpmb.lock();
 
@@ -205,7 +207,8 @@ std::int64_t Variable<T>::Deallocate() {
     }
   }
 
-  if (IsSet(Metadata::FillGhost) || IsSet(Metadata::Independent)) {
+  if (IsSet(Metadata::FillGhost) || IsSet(Metadata::Independent) ||
+      IsSet(Metadata::ForceRemeshComm)) {
     mem_size += coarse_s.size() * sizeof(T);
     coarse_s.Reset();
   }