Merge branch 'develop' into lroberts36/multigrid-reduce-block-size-first

Author: lroberts36

.github/workflows/ci-macos.yml

@@ -30,7 +30,16 @@ jobs:
 
     - name: Install dependencies
       run: |
-        brew install openmpi hdf5-mpi adios2 || true
+        brew install openmpi adios2 || true
+
+    # this is a workaround because homebrew's stable is hdf5 2.0.0 and cmake's
+    # find_package(HDF5) doesn't properly determine HDF5_IS_PARALLEL
+    # this *should* be fixed in cmake 4.3
+    - name: Install hdf5 from old formula
+      run: |
+        mkdir -p $(brew --repository)/Library/Taps/user/homebrew-repo
+        curl -o $(brew --repository)/Library/Taps/user/homebrew-repo/hdf5-mpi.rb https://raw.githubusercontent.com/Homebrew/homebrew-core/d03791212a4670e68f21c66939c5d348bc1d4bef/Formula/h/hdf5-mpi.rb
+        brew install user/repo/hdf5-mpi
 
     - name: Install parthenon_tools
       run: |

CHANGELOG.md

@@ -4,13 +4,15 @@
 
 ### Added (new features/APIs/variables/...)
 - [[PR 1346]](https://github.com/parthenon-hpc-lab/parthenon/pull/1346) Allow for user defined inter-level restrictions in multigrid 
+- [[PR 1308]](https://github.com/parthenon-hpc-lab/parthenon/pull/1308) Add additional indexing options for (b, type, ...indices) meshdata fields. Helpful for lower dimensional Metadata::None fields.
 - [[PR 1344]](https://github.com/parthenon-hpc-lab/parthenon/pull/1344) Add option to communicate single layer of ghosts, only communicate required two-level composite boundaries
 
 
 ### Changed (changing behavior/API/variables/...)
 - [[PR 1355]](https://github.com/parthenon-hpc-lab/parthenon/pull/1355) Allow disabling format and lint targets
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 1365]](https://github.com/parthenon-hpc-lab/parthenon/pull/1365) Fix boundary condition being called with coarse=true but no coarse neighbors.
 - [[PR 1345]](https://github.com/parthenon-hpc-lab/parthenon/pull/1345) Coalesce dot product reductions and speed up kernel
 - [[PR 1360]](https://github.com/parthenon-hpc-lab/parthenon/pull/1360) Fix boundary cache clearing in different MeshData partitions
 

external/Kokkos

@@ -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/CMakeLists.txt

@@ -41,6 +41,10 @@ TaskStatus ApplyBoundaryConditionsOnCoarseOrFine(std::shared_ptr<MeshBlockData<R
   Mesh *pmesh = pmb->pmy_mesh;
   const int ndim = pmesh->ndim;
 
+  // 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 && !pmb->HasCoarserNeighbors()) return TaskStatus::complete;
+
   auto &tree_bnd_func = pmesh->forest.GetTreePtr(pmb->loc.tree())->MeshBndryFnctn;
   auto &tree_bnd_func_user =
       pmesh->forest.GetTreePtr(pmb->loc.tree())->UserBoundaryFunctions;

src/bvals/boundary_conditions.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/bvals.cpp

@@ -441,13 +441,7 @@ 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;

src/bvals/comms/bnd_info.cpp

@@ -234,7 +234,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/boundary_communication.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.
@@ -143,18 +146,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/build_boundary_buffers.cpp

@@ -105,7 +105,7 @@ 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) {
+  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, md->GetBoundBufferId(bound_type));
     PARTHENON_DEBUG_REQUIRE(comm_map->count(key) > 0,
                             "Boundary communicator does not exist");
@@ -159,7 +159,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]);
 
@@ -191,7 +192,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()) {
@@ -241,19 +243,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();
 }