Forest Of Octrees

CHANGELOG.md

@@ -25,6 +25,7 @@
 - [[PR978]](https://github.com/parthenon-hpc-lab/parthenon/pull/978) remove erroneous sparse check
 
 ### Infrastructure (changes irrelevant to downstream codes)
+- [[PR 1009]](https://github.com/parthenon-hpc-lab/parthenon/pull/1009) Move from a single octree to a forest of octrees
 - [[PR 1017]](https://github.com/parthenon-hpc-lab/parthenon/pull/1017) Make regression tests more verbose on failure
 - [[PR 1007]](https://github.com/parthenon-hpc-lab/parthenon/pull/1007) Split template instantiations for HDF5 Read/Write attributes to speed up compile times
 - [[PR 990]](https://github.com/parthenon-hpc-lab/parthenon/pull/990) Partial refactor of HDF5 I/O code for readability/extendability

doc/sphinx/src/boundary_communication.rst

@@ -155,17 +155,7 @@ In practice, we denote each channel by a unique key
 so that the ``Mesh`` can contain a map from these keys to communication
 channels. Then, at each remesh, sending blocks and blocks that are
 receiving from blocks on a different rank can create new communication
-channels and register them in this map. *Implementation detail:* To
-build these keys, we currently rely on the
-``MeshBlock::std::unique_ptr<BoundaryValues> pbval`` object to get
-information about the neighboring blocks and build the channels and
-keys. ``BoundaryValues`` has its own communication methods defined, but
-none of these are used for the sparse communication. We really only rely
-on the information stored in ``BoundaryBase`` (which contains general
-information about all of the neighboring blocks on the mesh), which
-``BoundaryValues`` inherits from. Eventually, I think ``pbval`` should
-be turned into a ``BoundaryBase`` object and ``BoundaryValues`` should
-be removed from the code base.
+channels and register them in this map. 
 
 MPI Communication IDs
 ~~~~~~~~~~~~~~~~~~~~~

doc/sphinx/src/interface/boundary.rst

@@ -12,53 +12,3 @@ BoundaryCommunication
 
 Pure abstract base class, defines interfaces for managing
 ``BoundaryStatus`` flags and MPI requests
-
-BoundaryBuffer
---------------
-
-Pure abstract base class, defines interfaces for managing MPI
-send/receive and loading/storing data from communication buffers.
-
-BoundaryVariable
-----------------
-
-**Derived from**: ``BoundaryCommunication`` and ``BoundaryBuffer``
-
-**Contains**: ``BoundaryData`` for variable and flux correction.
-
-**Knows about**: ``MeshBlock`` and ``Mesh``
-
-Still abstract base class, but implements some methods for sending and
-receiving buffers.
-
-BoundaryBase
-------------
-
-**Contains**: ``NeighborIndexes`` and ``NeighborBlock`` PER NEIGHBOR,
-number of neighbors, neighbor levels
-
-**Knows about**: ``MeshBlock``
-
-Implements ``SearchAndSetNeighbors``
-
-BoundaryValues
---------------
-
-**Derived from**: ``BoundaryBase`` and ``BoundaryCommunication``
-
-**Knows about**: ``MeshBlock``, all the ``BoundaryVariable`` connected
-to variables of this block
-
-Central class to interact with individual variable boundary data. Owned
-by ``MeshBlock``.
-
-CellCenteredBoundaryVariable
-----------------------------
-
-**Derived from**: ``BoundaryVariable``
-
-**Contains**: Shallow copies of variable data, coarse buffer, and fluxes
-(owned by ``Variable``)
-
-Owned by ``Variable``, implements loading and setting boundary data,
-sending and receiving flux corrections, and more.

example/particles/particles.cpp

@@ -529,8 +529,8 @@ TaskStatus StopCommunicationMesh(const BlockList_t &blocks) {
 #ifdef MPI_PARALLEL
   for (auto &block : blocks) {
     auto swarm = block->swarm_data.Get()->Get("my_particles");
-    for (int n = 0; n < block->pbval->nneighbor; n++) {
-      NeighborBlock &nb = block->pbval->neighbor[n];
+    for (int n = 0; n < block->neighbors.size(); n++) {
+      NeighborBlock &nb = block->neighbors[n];
       // TODO(BRR) May want logic like this if we have non-blocking TaskRegions
       // if (nb.snb.rank != Globals::my_rank) {
       //  if (swarm->vbswarm->bd_var_.flag[nb.bufid] != BoundaryStatus::completed) {
@@ -563,8 +563,8 @@ TaskStatus StopCommunicationMesh(const BlockList_t &blocks) {
     auto &pmb = block;
     auto sc = pmb->swarm_data.Get();
     auto swarm = sc->Get("my_particles");
-    for (int n = 0; n < swarm->vbswarm->bd_var_.nbmax; n++) {
-      auto &nb = pmb->pbval->neighbor[n];
+    for (int n = 0; n < pmb->neighbors.size(); n++) {
+      auto &nb = block->neighbors[n];
       swarm->vbswarm->bd_var_.flag[nb.bufid] = BoundaryStatus::waiting;
     }
   }

scripts/python/packages/parthenon_performance_app/parthenon_performance_app/githubapp.py

@@ -76,7 +76,6 @@ def printTree(self):
 
 
 class GitHubApp:
-
     """
     GitHubApp Class
 

src/CMakeLists.txt

@@ -108,14 +108,11 @@ add_library(parthenon
   bvals/boundary_conditions_generic.hpp
   bvals/boundary_conditions.cpp
   bvals/boundary_conditions.hpp
-
   bvals/bvals.cpp
   bvals/bvals.hpp
-  bvals/bvals_base.cpp
-  bvals/bvals_interfaces.hpp
+  bvals/neighbor_block.cpp
+  bvals/neighbor_block.hpp
   bvals/boundary_flag.cpp
-  bvals/bvals_var.cpp
-  bvals/bvals_swarm.cpp
 
   coordinates/coordinates.hpp
   coordinates/uniform_cartesian.hpp
@@ -161,6 +158,8 @@ add_library(parthenon
 
   mesh/amr_loadbalance.cpp
   mesh/domain.hpp
+  mesh/forest.cpp
+  mesh/forest.hpp
   mesh/logical_location.cpp
   mesh/logical_location.hpp
   mesh/mesh_refinement.cpp
@@ -170,8 +169,6 @@ add_library(parthenon
   mesh/mesh.hpp
   mesh/meshblock.hpp
   mesh/meshblock_pack.hpp
-  mesh/meshblock_tree.cpp
-  mesh/meshblock_tree.hpp
   mesh/meshblock.cpp
 
   outputs/ascent.cpp

src/bvals/boundary_conditions.cpp

@@ -16,7 +16,7 @@
 
 #include "bvals/boundary_conditions.hpp"
 #include "bvals/boundary_conditions_generic.hpp"
-#include "bvals/bvals_interfaces.hpp"
+#include "bvals/neighbor_block.hpp"
 #include "defs.hpp"
 #include "interface/meshblock_data.hpp"
 #include "mesh/domain.hpp"

src/bvals/boundary_flag.cpp

@@ -90,9 +90,8 @@ std::string GetBoundaryString(BoundaryFlag input_flag) {
 //! \fn CheckBoundaryFlag(BoundaryFlag block_flag, CoordinateDirection dir)
 //  \brief Called in each MeshBlock's BoundaryValues() constructor. Mesh() ctor only
 //  checks the validity of user's input mesh/ixn_bc, oxn_bc string values corresponding to
-//  a BoundaryFlag enumerator before passing it to a MeshBlock and then BoundaryBase
-//  object. However, not all BoundaryFlag enumerators can be used in all directions as a
-//  valid MeshBlock boundary.
+//  a BoundaryFlag enumerator before passing it to a MeshBlock. However, not all
+//  BoundaryFlag enumerators can be used in all directions as a valid MeshBlock boundary.
 
 void CheckBoundaryFlag(BoundaryFlag block_flag, CoordinateDirection dir) {
   std::stringstream msg;

src/bvals/bvals.cpp

@@ -14,8 +14,6 @@
 // license in this material to reproduce, prepare derivative works, distribute copies to
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
-//! \file bvals.cpp
-//  \brief constructor/destructor and utility functions for BoundaryValues class
 
 #include "bvals/bvals.hpp"
 
@@ -45,93 +43,131 @@
 
 namespace parthenon {
 
-// BoundaryValues constructor (the first object constructed inside the MeshBlock()
-// constructor): sets functions for the appropriate boundary conditions at each of the 6
-// dirs of a MeshBlock
-BoundaryValues::BoundaryValues(std::weak_ptr<MeshBlock> wpmb, BoundaryFlag *input_bcs,
-                               ParameterInput *pin)
-    : BoundaryBase(wpmb.lock()->pmy_mesh, wpmb.lock()->loc, wpmb.lock()->block_size,
-                   input_bcs),
-      pmy_block_(wpmb) {
-  // Check BC functions for each of the 6 boundaries in turn ---------------------
-  for (int i = 0; i < 6; i++) {
-    switch (block_bcs[i]) {
-    case BoundaryFlag::reflect:
-    case BoundaryFlag::outflow:
-      apply_bndry_fn_[i] = true;
-      break;
-    default: // already initialized to false in class
-      break;
-    }
-  }
-  // Inner x1
-  nface_ = 2;
-  nedge_ = 0;
-  CheckBoundaryFlag(block_bcs[BoundaryFace::inner_x1], CoordinateDirection::X1DIR);
-  CheckBoundaryFlag(block_bcs[BoundaryFace::outer_x1], CoordinateDirection::X1DIR);
-
-  std::shared_ptr<MeshBlock> pmb = GetBlockPointer();
-  if (!pmb->block_size.symmetry(X2DIR)) {
-    nface_ = 4;
-    nedge_ = 4;
-    CheckBoundaryFlag(block_bcs[BoundaryFace::inner_x2], CoordinateDirection::X2DIR);
-    CheckBoundaryFlag(block_bcs[BoundaryFace::outer_x2], CoordinateDirection::X2DIR);
-  }
+BoundarySwarm::BoundarySwarm(std::weak_ptr<MeshBlock> pmb, const std::string &label)
+    : bswarm_index(), pmy_block(pmb), pmy_mesh_(pmb.lock()->pmy_mesh) {
+#ifdef MPI_PARALLEL
+  swarm_comm = pmy_mesh_->GetMPIComm(label);
+#endif
+  InitBoundaryData(bd_var_);
+}
 
-  if (!pmb->block_size.symmetry(X3DIR)) {
-    nface_ = 6;
-    nedge_ = 12;
-    CheckBoundaryFlag(block_bcs[BoundaryFace::inner_x3], CoordinateDirection::X3DIR);
-    CheckBoundaryFlag(block_bcs[BoundaryFace::outer_x3], CoordinateDirection::X3DIR);
+void BoundarySwarm::InitBoundaryData(BoundaryData<> &bd) {
+  auto pmb = GetBlockPointer();
+  BufferID buffer_id(pmb->pmy_mesh->ndim, pmb->pmy_mesh->multilevel);
+  bd.nbmax = buffer_id.size();
+
+  for (int n = 0; n < bd.nbmax; n++) {
+    bd.flag[n] = BoundaryStatus::waiting;
+#ifdef MPI_PARALLEL
+    bd.req_send[n] = MPI_REQUEST_NULL;
+    bd.req_recv[n] = MPI_REQUEST_NULL;
+#endif
   }
-
-  // prevent reallocation of contiguous memory space for each of 4x possible calls to
-  // std::vector<BoundaryVariable *>.push_back() in Field, PassiveScalars
-  bvars.reserve(3);
 }
 
-// destructor
-
-//----------------------------------------------------------------------------------------
-//! \fn void BoundaryValues::SetupPersistentMPI()
-//  \brief Setup persistent MPI requests to be reused throughout the entire simulation
+void BoundarySwarm::SetupPersistentMPI() {
+#ifdef MPI_PARALLEL
+  std::shared_ptr<MeshBlock> pmb = GetBlockPointer();
 
-void BoundaryValues::SetupPersistentMPI() {
-  for (auto bvars_it = bvars.begin(); bvars_it != bvars.end(); ++bvars_it) {
-    (*bvars_it).second->SetupPersistentMPI();
+  // Initialize neighbor communications to other ranks
+  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) {
+      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) {
+        MPI_Request_free(&bd_var_.req_send[nb.bufid]);
+      }
+      if (bd_var_.req_recv[nb.bufid] != MPI_REQUEST_NULL) {
+        MPI_Request_free(&bd_var_.req_recv[nb.bufid]);
+      }
+    }
   }
+#endif
 }
 
-//----------------------------------------------------------------------------------------
-//! \fn void BoundaryValues::StartReceiving(BoundaryCommSubset phase)
-//  \brief initiate MPI_Irecv()
-
-void BoundaryValues::StartReceiving(BoundaryCommSubset phase) {
-  for (auto bvars_it = bvars.begin(); bvars_it != bvars.end(); ++bvars_it) {
-    (*bvars_it).second->StartReceiving(phase);
+// Send particle buffers across meshblocks. If different MPI ranks, use MPI, if same rank,
+// do a deep copy on device.
+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];
+    if (nb.snb.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],
+                                    swarm_comm, &(bd_var_.req_send[nb.bufid])));
+#endif // MPI_PARALLEL
+    } else {
+      MeshBlock &target_block = *pmy_mesh_->FindMeshBlock(nb.snb.gid);
+      std::shared_ptr<BoundarySwarm> ptarget_bswarm =
+          target_block.pbswarm->bswarms[bswarm_index];
+      if (send_size[nb.bufid] > 0) {
+        // Ensure target buffer is large enough
+        if (bd_var_.send[nb.bufid].extent(0) >
+            ptarget_bswarm->bd_var_.recv[nb.targetid].extent(0)) {
+          ptarget_bswarm->bd_var_.recv[nb.targetid] =
+              BufArray1D<Real>("Buffer", (bd_var_.send[nb.bufid].extent(0)));
+        }
+
+        target_block.deep_copy(ptarget_bswarm->bd_var_.recv[nb.targetid],
+                               bd_var_.send[nb.bufid]);
+        ptarget_bswarm->recv_size[nb.targetid] = send_size[nb.bufid];
+        ptarget_bswarm->bd_var_.flag[nb.targetid] = BoundaryStatus::arrived;
+      } else {
+        ptarget_bswarm->recv_size[nb.targetid] = 0;
+        ptarget_bswarm->bd_var_.flag[nb.targetid] = BoundaryStatus::completed;
+      }
+    }
   }
 }
 
-//----------------------------------------------------------------------------------------
-//! \fn void BoundaryValues::ClearBoundary(BoundaryCommSubset phase)
-//  \brief clean up the boundary flags after each loop
-
-void BoundaryValues::ClearBoundary(BoundaryCommSubset phase) {
-  // Note BoundaryCommSubset::mesh_init corresponds to initial exchange of conserved fluid
-  // variables and magentic fields
-  for (auto bvars_it = bvars.begin(); bvars_it != bvars.end(); ++bvars_it) {
-    (*bvars_it).second->ClearBoundary(phase);
+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];
+    if (nb.snb.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));
+        if (!static_cast<bool>(test)) {
+          bd_var_.flag[nb.bufid] = BoundaryStatus::waiting;
+        } else {
+          bd_var_.flag[nb.bufid] = BoundaryStatus::arrived;
+
+          // If message is available, receive it
+          PARTHENON_MPI_CHECK(
+              MPI_Get_count(&status, MPI_PARTHENON_REAL, &(recv_size[nb.bufid])));
+          if (recv_size[nb.bufid] > bd_var_.recv[nb.bufid].extent(0)) {
+            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));
+        }
+      }
+    }
   }
+#endif
 }
 
 // BoundarySwarms constructor (the first object constructed inside the MeshBlock()
 // constructor): sets functions for the appropriate boundary conditions at each of the 6
 // dirs of a MeshBlock
 BoundarySwarms::BoundarySwarms(std::weak_ptr<MeshBlock> wpmb, BoundaryFlag *input_bcs,
                                ParameterInput *pin)
-    : BoundaryBase(wpmb.lock()->pmy_mesh, wpmb.lock()->loc, wpmb.lock()->block_size,
-                   input_bcs),
-      pmy_block_(wpmb) {
+    : pmy_block_(wpmb) {
   // Check BC functions for each of the 6 boundaries in turn ---------------------
   // TODO(BRR) Add physical particle boundary conditions, maybe using the below code
   /*for (int i = 0; i < 6; i++) {
@@ -144,6 +180,9 @@ BoundarySwarms::BoundarySwarms(std::weak_ptr<MeshBlock> wpmb, BoundaryFlag *inpu
       break;
     }
   }*/
+  for (int i = 0; i < 6; ++i)
+    block_bcs[i] = input_bcs[i];
+
   // Inner x1
   nface_ = 2;
   nedge_ = 0;
@@ -167,7 +206,7 @@ BoundarySwarms::BoundarySwarms(std::weak_ptr<MeshBlock> wpmb, BoundaryFlag *inpu
 }
 
 //----------------------------------------------------------------------------------------
-//! \fn void BoundaryValues::SetupPersistentMPI()
+//! \fn void BoundarySwarms::SetupPersistentMPI()
 //  \brief Setup persistent MPI requests to be reused throughout the entire simulation
 
 void BoundarySwarms::SetupPersistentMPI() {