Cache PackDescriptors

CHANGELOG.md

@@ -8,6 +8,7 @@
 
 
 ### Changed (changing behavior/API/variables/...)
+- [[PR 1239]](https://github.com/parthenon-hpc-lab/parthenon/pull/1239) Automatically cache PackDescriptors in Mesh
 - [[PR 1242]](https://github.com/parthenon-hpc-lab/parthenon/pull/1242) Move to Kokkos 4.6.01 (for AMD APU support)
 - [[PR 1253]](https://github.com/parthenon-hpc-lab/parthenon/pull/1253) Add support for uint64 swarm variables and add default id
 

doc/sphinx/src/particles.rst

@@ -155,7 +155,7 @@ Similarly to grid variables, particle swarms support
 This also supports ``FlatIdx`` for indexing; see the
 ``particle_leapfrog`` example for usage.
 
-. _swarm_packs:
+.. _swarm_packs:
 
 ``SwarmPack``\ s
 ----------------

doc/sphinx/src/sparse_packs.rst

@@ -16,7 +16,7 @@ etc. within a kernel.
 
 ``SparsePack``\ s *work for all types of variables (both dense and sparse). They were originally implemented to support sparse variables and to supersede the older* ``VariablePack``\ s *and* ``VariableFluxPack``\ s *and picked up the* ``Sparse`` *modifier to differentiate them. The latter have not been removed from ``Parthenon`` because some downstream codes still rely on them, but they are deprecated and will be removed eventually.*
 
-*If you want to deal with particle fields, you will need to use* ``SwarmPack``\ s, *which are described at :ref:`swarm_packs`.*
+*If you want to deal with particle fields, you will need to use* ``SwarmPack``\ *s, which are described at* :ref:`swarm_packs`.
 
 Type-based Packing
 ------------------
@@ -64,7 +64,7 @@ Building and Using a ``SparsePack``
      return TaskStatus::complete;
    }
 
-``PackDescriptor``\ *s are somewhat expensive to build and are not currently cached. As a result, it is often beneficial to mark them as ``static`` where possible.*
+``PackDescriptor``\ s can be somewhat expensive to build because they require searching through all fields in simulation. Therefore, they are automatically cached in the ``StateDescriptor`` where possible. Additionally, it is often possible to declare ``PackDescriptors`` that are created in task functions to be ``static``.  
 
 ``PackDescriptor`` takes a ``std::set`` of ``PDOpt`` options to determine what to include in the pack:
 

doc/sphinx/src/utilities.rst

@@ -1,10 +1,13 @@
+Utilities
+=========
+
 ``TypeList``\ s
-===============
+---------------
 
 Provides a wrapper class around a variadic pack of types to simplify
 performing compile time operations on the pack. There are templated 
 types defined giving the type at a particular index in the pack, types 
-for extracting sub-``TypeList``s of the original type list, and ``constexpr``
+for extracting sub-``TypeList``\ s of the original type list, and ``constexpr``
 functions for getting the index of the first instance of a type in the 
 pack. Additionally it provides a capability for iterating an ``auto`` lambda
 over the type list, which can be useful for calling a ``static`` function 
@@ -19,20 +22,20 @@ lists of types associated with field variables that are used in type
 based ``SparsePack``\ s.  
 
 Robust
-======
+------
 
 Provides a number of functions for doing operations on floating point
 numbers that are bounded, avoid division by zero, etc. 
 
 C++11 Style Concepts Implementation
-===================================
+-----------------------------------
 
 *This documentation needs to be written (see issue #695), but there are
 extensive comments in src/utlils/concepts_lite.hpp and examples of
 useage in tst/unit/test_concepts_lite.hpp*
 
 ``Indexer``
-===========
+-----------
 
 Provides functionality for iterating over an arbitrary dimensional
 hyper-rectangular index space using a flattened loop. Specific 

example/fine_advection/advection_driver.cpp

@@ -103,19 +103,16 @@ TaskCollection AdvectionDriver::MakeTaskCollection(BlockList_t &blocks, const in
     // to the Stokes update routines, so that they can internally create variable
     // packs that operate on only the desired set of variables.
     using namespace advection_package::Conserved;
-    static auto desc = parthenon::MakePackDescriptor<phi>(
-        pmesh->resolved_packages.get(), {parthenon::Metadata::WithFluxes},
-        {parthenon::PDOpt::WithFluxes});
+    auto desc = parthenon::MakePackDescriptor<phi>(
+        pmesh, {parthenon::Metadata::WithFluxes}, {parthenon::PDOpt::WithFluxes});
     using pack_desc_t = decltype(desc);
 
-    static auto desc_fine = parthenon::MakePackDescriptor<phi_fine>(
-        pmesh->resolved_packages.get(), {parthenon::Metadata::WithFluxes},
-        {parthenon::PDOpt::WithFluxes});
+    auto desc_fine = parthenon::MakePackDescriptor<phi_fine>(
+        pmesh, {parthenon::Metadata::WithFluxes}, {parthenon::PDOpt::WithFluxes});
     using pack_desc_fine_t = decltype(desc_fine);
 
-    static auto desc_vec = parthenon::MakePackDescriptor<C, D>(
-        pmesh->resolved_packages.get(), {parthenon::Metadata::WithFluxes},
-        {parthenon::PDOpt::WithFluxes});
+    auto desc_vec = parthenon::MakePackDescriptor<C, D>(
+        pmesh, {parthenon::Metadata::WithFluxes}, {parthenon::PDOpt::WithFluxes});
 
     using TT = parthenon::TopologicalType;
     using TE = parthenon::TopologicalElement;

example/fine_advection/advection_package.cpp

@@ -130,7 +130,7 @@ void CheckRefinementMesh(MeshData<Real> *md, parthenon::ParArray1D<AmrTag> &amr_
   auto do_regular_advection = pkg->Param<bool>("do_regular_advection");
   if (do_regular_advection) {
     // refine on advected, for example.  could also be a derived quantity
-    static auto desc = parthenon::MakePackDescriptor<Conserved::phi>(md);
+    auto desc = parthenon::MakePackDescriptor<Conserved::phi>(md);
     auto pack = desc.GetPack(md);
 
     const auto &refine_tol = pkg->Param<Real>("refine_tol");
@@ -179,7 +179,7 @@ Real EstimateTimestep(MeshData<Real> *md) {
   const auto &vy = pkg->Param<Real>("vy");
   const auto &vz = pkg->Param<Real>("vz");
 
-  static auto desc = parthenon::MakePackDescriptor<Conserved::phi>(md);
+  auto desc = parthenon::MakePackDescriptor<Conserved::phi>(md);
   auto pack = desc.GetPack(md);
 
   IndexRange ib = md->GetBoundsI(IndexDomain::interior);
@@ -206,7 +206,7 @@ Real EstimateTimestep(MeshData<Real> *md) {
 }
 
 TaskStatus FillDerived(MeshData<Real> *md) {
-  static auto desc =
+  auto desc =
       parthenon::MakePackDescriptor<Conserved::phi_fine, Conserved::phi_fine_restricted,
                                     Conserved::C, Conserved::C_cc, Conserved::D,
                                     Conserved::D_cc, Conserved::divC, Conserved::divD>(

src/CMakeLists.txt

@@ -211,6 +211,7 @@ add_library(parthenon
 
   pack/block_selector.hpp
   pack/make_pack_descriptor.hpp
+  pack/make_pack_descriptor.cpp
   pack/make_swarm_pack_descriptor.hpp
   pack/pack_descriptor.cpp
   pack/pack_descriptor.hpp

src/amr_criteria/refinement_package.cpp

@@ -101,8 +101,7 @@ void FirstDerivative(const AMRBounds &bnds, MeshData<Real> *md, const std::strin
                      const int &idx, ParArray1D<AmrTag> &amr_tags,
                      const Real refine_criteria_, const Real derefine_criteria_,
                      const int max_level_) {
-  const auto desc =
-      MakePackDescriptor(md->GetMeshPointer()->resolved_packages.get(), {field});
+  const auto desc = MakePackDescriptor(md, {field});
   auto pack = desc.GetPack(md);
   const int ndim = md->GetMeshPointer()->ndim;
   const int nvars = pack.GetMaxNumberOfVars();
@@ -156,8 +155,7 @@ void SecondDerivative(const AMRBounds &bnds, MeshData<Real> *md, const std::stri
                       const int &idx, ParArray1D<AmrTag> &amr_tags,
                       const Real refine_criteria_, const Real derefine_criteria_,
                       const int max_level_) {
-  const auto desc =
-      MakePackDescriptor(md->GetMeshPointer()->resolved_packages.get(), {field});
+  const auto desc = MakePackDescriptor(md, {field});
   auto pack = desc.GetPack(md);
   const int ndim = md->GetMeshPointer()->ndim;
   const int nvars = pack.GetMaxNumberOfVars();

src/bvals/boundary_conditions_generic.hpp

@@ -137,40 +137,6 @@ void GenericSwarmBC(std::shared_ptr<Swarm> &swarm) {
       });
 }
 
-namespace impl {
-using desc_key_t = std::tuple<bool, bool, TopologicalType>;
-template <class... var_ts>
-using map_bc_pack_descriptor_t =
-    std::unordered_map<desc_key_t, typename SparsePack<var_ts...>::Descriptor,
-                       tuple_hash<desc_key_t>>;
-
-template <class... var_ts>
-map_bc_pack_descriptor_t<var_ts...>
-GetPackDescriptorMap(std::shared_ptr<MeshBlockData<Real>> &rc) {
-  std::vector<std::pair<TopologicalType, MetadataFlag>> elements{
-      {TopologicalType::Cell, Metadata::Cell},
-      {TopologicalType::Face, Metadata::Face},
-      {TopologicalType::Edge, Metadata::Edge},
-      {TopologicalType::Node, Metadata::Node}};
-  map_bc_pack_descriptor_t<var_ts...> my_map;
-  for (auto [tt, md] : elements) {
-    std::vector<MetadataFlag> flags{Metadata::FillGhost};
-    flags.push_back(md);
-    std::set<PDOpt> opts{PDOpt::Coarse};
-    my_map.emplace(std::make_pair(desc_key_t{true, false, tt},
-                                  MakePackDescriptor<var_ts...>(rc.get(), flags, opts)));
-    my_map.emplace(std::make_pair(desc_key_t{false, false, tt},
-                                  MakePackDescriptor<var_ts...>(rc.get(), flags)));
-    flags.push_back(Metadata::Fine);
-    my_map.emplace(std::make_pair(desc_key_t{true, true, tt},
-                                  MakePackDescriptor<var_ts...>(rc.get(), flags, opts)));
-    my_map.emplace(std::make_pair(desc_key_t{false, true, tt},
-                                  MakePackDescriptor<var_ts...>(rc.get(), flags)));
-  }
-  return my_map;
-}
-} // namespace impl
-
 template <CoordinateDirection DIR, BCSide SIDE, BCType TYPE, class... var_ts>
 void GenericBC(std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse,
                TopologicalElement el, Real val) {
@@ -183,10 +149,29 @@ void GenericBC(std::shared_ptr<MeshBlockData<Real>> &rc, bool coarse,
   constexpr bool X3 = (DIR == X3DIR);
   constexpr bool INNER = (SIDE == BCSide::Inner);
 
-  static auto descriptors = impl::GetPackDescriptorMap<var_ts...>(rc);
+  const auto ttFlag = [el] {
+    const auto tt = GetTopologicalType(el);
+    switch (tt) {
+    case (TopologicalType::Cell):
+      return Metadata::Cell;
+    case (TopologicalType::Face):
+      return Metadata::Face;
+    case (TopologicalType::Edge):
+      return Metadata::Edge;
+    case (TopologicalType::Node):
+      return Metadata::Node;
+    default:
+      PARTHENON_FAIL("Unknown topological type")
+    }
+  }();
+
+  // static auto descriptors = impl::GetPackDescriptorMap<var_ts...>(rc);
   for (auto fine : {false, true}) {
-    auto q = descriptors[impl::desc_key_t{coarse, fine, GetTopologicalType(el)}].GetPack(
-        rc.get());
+    std::vector<MetadataFlag> flags{Metadata::FillGhost, ttFlag};
+    if (fine) flags.push_back(Metadata::Fine);
+    std::set<PDOpt> opts = coarse ? std::set<PDOpt>{PDOpt::Coarse} : std::set<PDOpt>{};
+    const auto desc = MakePackDescriptor<var_ts...>(rc.get(), flags, opts);
+    auto q = desc.GetPack(rc.get());
     const int b = 0;
     const int lstart = q.GetLowerBoundHost(b);
     const int lend = q.GetUpperBoundHost(b);

src/interface/metadata.hpp

@@ -177,6 +177,7 @@ class MetadataFlag {
 
   std::string const &Name() const;
 
+  int Flag() const { return flag_; }
 #ifdef CATCH_VERSION_MAJOR
   // Should never be used for application code - only exposed for testing.
   constexpr int InternalFlagValue() const { return flag_; }
@@ -705,4 +706,11 @@ Set_t GetByFlag(const Metadata::FlagCollection &flags, NameMap_t &nameMap,
 } // namespace MetadataUtils
 } // namespace parthenon
 
+template <>
+struct std::hash<parthenon::MetadataFlag> {
+  std::size_t operator()(const parthenon::MetadataFlag &flag) const {
+    return flag.Flag();
+  }
+};
+
 #endif // INTERFACE_METADATA_HPP_

src/interface/state_descriptor.hpp

@@ -43,6 +43,7 @@ template <typename T>
 class MeshData;
 class AMRCriteria;
 class Packages_t;
+class PackDescriptorCacheBase;
 
 /// A little container class owning refinement function properties
 /// needed for the state descriptor.
@@ -403,6 +404,10 @@ class StateDescriptor {
   std::array<std::vector<BValFunc>, BOUNDARY_NFACES> UserBoundaryFunctions;
   std::array<std::vector<SBValFunc>, BOUNDARY_NFACES> UserSwarmBoundaryFunctions;
 
+  // Caches for PackDescriptors associated with this StateDescriptor
+  std::unordered_map<std::string, std::shared_ptr<PackDescriptorCacheBase>>
+      pack_desc_cache_map;
+
  protected:
   // internal function to add dense/sparse fields. Private because outside classes must
   // use the public interface below

src/interface/update.cpp

@@ -151,8 +151,7 @@ TaskStatus SparseDealloc(MeshData<Real> *md) {
   const IndexRange kb = md->GetBoundsK(IndexDomain::entire);
 
   auto control_vars = md->GetMeshPointer()->resolved_packages->GetControlVariables();
-  auto desc = MakePackDescriptor(md->GetMeshPointer()->resolved_packages.get(),
-                                 control_vars, {Metadata::Sparse});
+  auto desc = MakePackDescriptor(md, control_vars, {Metadata::Sparse});
   auto pack = desc.GetPack(md);
   auto packIdx = desc.GetMap();
 

src/interface/update.hpp

@@ -329,8 +329,8 @@ TaskStatus InitNewlyAllocatedVars(T *rc) {
     // This pack will always be freshly built, since we only get here if sparse data
     // was allocated and hasn't been initialized, which in turn implies the cached
     // pack must be stale.
-    auto desc = parthenon::MakePackDescriptor<variable_names::any>(
-        rc->GetMeshPointer()->resolved_packages.get(), {Metadata::Sparse});
+    auto desc =
+        parthenon::MakePackDescriptor<variable_names::any>(rc, {Metadata::Sparse});
     auto v = desc.GetPack(rc);
 
     Kokkos::parallel_for(

src/mesh/mesh.hpp

@@ -48,6 +48,7 @@
 #include "mesh/forest/forest_topology.hpp"
 #include "mesh/meshblock_pack.hpp"
 #include "outputs/io_wrapper.hpp"
+#include "pack/pack_descriptor.hpp"
 #include "parameter_input.hpp"
 #include "parthenon_arrays.hpp"
 #include "utils/communication_buffer.hpp"
@@ -68,6 +69,8 @@ class RestartReader;
 // Map from LogicalLocation to (gid, rank) pair of location
 using LogicalLocMap_t = std::map<LogicalLocation, std::pair<int, int>>;
 
+// Base class to allow cacheing of different types of PackDescriptors
+
 //----------------------------------------------------------------------------------------
 //! \class Mesh
 //  \brief data/functions associated with the overall mesh
@@ -240,8 +243,7 @@ class Mesh {
   using channel_key_t = std::tuple<int, int, std::string, int, int>;
   using comm_buf_t = CommBuffer<buf_pool_t<Real>::owner_t>;
   std::unordered_map<int, buf_pool_t<Real>> pool_map;
-  using comm_buf_map_t =
-      std::unordered_map<channel_key_t, comm_buf_t, tuple_hash<channel_key_t>>;
+  using comm_buf_map_t = std::unordered_map<channel_key_t, comm_buf_t>;
   comm_buf_map_t boundary_comm_map;
   TagMap tag_map;