merge main

CHANGELOG.md

@@ -3,11 +3,11 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
-- [[PR 1171]](https://github.com/parthenon-hpc-lab/parthenon/pull/1171) Add PARTHENON_USE_SYSTEM_PACKAGES build option
-
+- [[PR 1161]](https://github.com/parthenon-hpc-lab/parthenon/pull/1161) Make flux field Metadata accessible, add Metadata::CellMemAligned flag, small perfomance upgrades
 
 ### Changed (changing behavior/API/variables/...)
-- [[PR1172]](https://github.com/parthenon-hpc-lab/parthenon/pull/1172) Make parthenon manager robust against external MPI init and finalize calls
+- [[PR 1171]](https://github.com/parthenon-hpc-lab/parthenon/pull/1171) Add PARTHENON_USE_SYSTEM_PACKAGES build option
+- [[PR 1172]](https://github.com/parthenon-hpc-lab/parthenon/pull/1172) Make parthenon manager robust against external MPI init and finalize calls
 
 ### Fixed (not changing behavior/API/variables/...)
 

doc/sphinx/src/interface/metadata.rst

@@ -158,9 +158,22 @@ classes may be allocated. The behaviours are the following:
   sense (e.g. if the ``WithFluxes`` variable has ``Metadata::Cell`` 
   set the new variable will have ``Metadata::Face``) will be created in
   the package with the name ``bnd_flux::name_of_original_variable`` and 
-  ``Metadata::Flux`` and ``Metadata::OneCopy``. When creating packs that 
-  include fluxes, the new flux field will be included in the flux portion 
-  of the pack if the parent field is in the pack. 
+  ``Metadata::Flux`` and ``Metadata::OneCopy``. Additionally, the flags 
+  ``Metadata::Sparse``, ``Metadata::Vector``, and ``Metadata::Tensor`` 
+  will propagate to the flux ``Metadata`` if they are set in the base field
+  ``Metadata``. By default, a flux field for a cell-centered field is 
+  built with ``Metadata::CellMemAligned`` flag set for backwards 
+  compatability. A shared pointer to the ``Metadata`` object for the flux
+  field can be accessed from the base ``Metadata`` with the method 
+  ``Metadata::GetSPtrFluxMetadata()``. This can be used to set flags other
+  than the defaults or set custom prolongation/restriction operations for
+  the fluxes. Note that calling `Metadata::RegisterRefinementOps<...>()`
+  on the base field propagates the registered refinement operations through
+  to the flux `Metadata` for backward compatibility. If separate operations 
+  are desired for the fluxes, the ordering of calls to `RegisterRefinementOps` 
+  on the base field and the flux field matters. When creating packs that 
+  include fluxes, the new flux field will be included in the flux portion of
+  the pack if the parent field is in the pack. 
 
 - If ``Metadata::Flux`` is set, this field is exchanged on shared elements 
   across fine-coarse boundaries when the flux correction tasks are called. 

src/bvals/comms/boundary_communication.cpp

@@ -213,7 +213,7 @@ TaskStatus ReceiveBoundBufs(std::shared_ptr<MeshData<Real>> &md) {
       [&all_received](auto pbuf) { all_received = pbuf->TryReceive() && all_received; });
 
   int ibound = 0;
-  if (Globals::sparse_config.enabled) {
+  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) {
           const std::size_t ibuf = cache.idx_vec[ibound];

src/interface/meshblock_data.hpp

@@ -77,6 +77,9 @@ class MeshBlockData {
   void SetAllowedDt(const Real dt) const { GetBlockPointer()->SetAllowedDt(dt); }
   Mesh *GetMeshPointer() const { return GetBlockPointer()->pmy_mesh; }
 
+  // This mirrors a MeshBlockData routine
+  int NumBlocks() const { return 1; }
+
   template <class... Ts>
   IndexRange GetBoundsI(Ts &&...args) const {
     return GetBlockPointer()->cellbounds.GetBoundsI(std::forward<Ts>(args)...);

src/interface/metadata.cpp

@@ -101,10 +101,13 @@ MetadataFlag Metadata::GetUserFlag(const std::string &flagname) {
 }
 
 namespace parthenon {
-Metadata::Metadata(const std::vector<MetadataFlag> &bits, const std::vector<int> &shape,
+Metadata::Metadata(const std::vector<MetadataFlag> &bits,
+                   const std::vector<MetadataFlag> &flux_bits,
+                   const std::vector<int> &shape,
                    const std::vector<std::string> &component_labels,
                    const std::string &associated,
-                   const refinement::RefinementFunctions_t ref_funcs_)
+                   const refinement::RefinementFunctions_t ref_funcs_,
+                   const refinement::RefinementFunctions_t flux_ref_funcs_)
     : shape_(shape), component_labels_(component_labels), associated_(associated) {
   // set flags
   for (const auto f : bits) {
@@ -164,6 +167,39 @@ Metadata::Metadata(const std::vector<MetadataFlag> &bits, const std::vector<int>
     deallocation_threshold_ = 0.0;
     default_value_ = 0.0;
   }
+
+  // Now create the flux metadata if required
+  if (IsSet(WithFluxes)) {
+    std::set<MetadataFlag> flux_flags;
+    for (const auto f : flux_bits)
+      flux_flags.insert(f);
+
+    // Set some standard defaults for the flux metadata if no
+    // flags were provided
+    if (flux_flags.size() == 0) {
+      flux_flags.insert(OneCopy);
+      if (IsSet(Fine)) flux_flags.insert(Fine);
+      if (IsSet(Cell)) flux_flags.insert(CellMemAligned);
+      if (IsSet(Sparse)) flux_flags.insert(Sparse);
+    }
+
+    // These flags are automatically propagated for fluxes
+    flux_flags.insert(Flux);
+    if (IsSet(Cell)) {
+      flux_flags.insert(Face);
+    } else if (IsSet(Face)) {
+      flux_flags.insert(Edge);
+    } else if (IsSet(Edge)) {
+      flux_flags.insert(Node);
+    }
+
+    if (IsSet(Tensor)) flux_flags.insert(Tensor);
+    if (IsSet(Vector)) flux_flags.insert(Vector);
+
+    std::vector<MetadataFlag> flux_flags_vec(flux_flags.begin(), flux_flags.end());
+    flux_metadata = std::make_shared<Metadata>(flux_flags_vec, shape, component_labels,
+                                               std::string(), flux_ref_funcs_);
+  }
 }
 
 std::ostream &operator<<(std::ostream &os, const parthenon::Metadata &m) {
@@ -252,6 +288,14 @@ bool Metadata::IsValid(bool throw_on_fail) const {
     }
   }
 
+  if (IsSet(FillGhost) && IsSet(CellMemAligned) && (!IsSet(Cell))) {
+    valid = false;
+    if (throw_on_fail) {
+      PARTHENON_THROW(
+          "Cannot communicate ghosts of non-cell fields that have cell aligned memory.");
+    }
+  }
+
   // Prolongation/restriction
   if (HasRefinementOps()) {
     if (refinement_funcs_.label().size() == 0) {
@@ -328,16 +372,12 @@ Metadata::GetArrayDims(std::weak_ptr<MeshBlock> wpmb, bool coarse) const {
       arrDims[i + 3] = 1;
     if (IsSet(Cell)) {
       arrDims[MAX_VARIABLE_DIMENSION - 1] = 1; // Only one cell center per cell
-    } else if (IsSet(Face) && IsSet(Flux)) {
-      // 3 directions but keep the same ijk shape as cell var for performance
-      arrDims[MAX_VARIABLE_DIMENSION - 1] = 3;
     } else if (IsSet(Face) || IsSet(Edge)) {
       arrDims[MAX_VARIABLE_DIMENSION - 1] = 3; // Three faces and edges per cell
-      arrDims[0]++;
-      if (arrDims[1] > 1) arrDims[1]++;
-      if (arrDims[2] > 1) arrDims[2]++;
     } else if (IsSet(Node)) {
       arrDims[MAX_VARIABLE_DIMENSION - 1] = 1; // Only one lower left node per cell
+    }
+    if (!IsSet(CellMemAligned) && !IsSet(Cell)) {
       arrDims[0]++;
       if (arrDims[1] > 1) arrDims[1]++;
       if (arrDims[2] > 1) arrDims[2]++;

src/interface/metadata.hpp

@@ -120,6 +120,9 @@
   PARTHENON_INTERNAL_FOR_FLAG(Fine)                                                      \
   /** this variable is the flux for another variable **/                                 \
   PARTHENON_INTERNAL_FOR_FLAG(Flux)                                                      \
+  /** Align memory of fields to cell centered memory                                     \
+      (Field will be missing one layer of ghosts if it is not cell centered) **/         \
+  PARTHENON_INTERNAL_FOR_FLAG(CellMemAligned)                                            \
   /************************************************/                                     \
   /** Vars specifying coordinates for visualization purposes **/                         \
   /** You can specify a single 3D var **/                                                \
@@ -325,28 +328,48 @@ class Metadata {
   // 4 constructors, this is the general constructor called by all other constructors, so
   // we do some sanity checks here
   Metadata(
-      const std::vector<MetadataFlag> &bits, const std::vector<int> &shape = {},
+      const std::vector<MetadataFlag> &bits, const std::vector<MetadataFlag> &flux_bits,
+      const std::vector<int> &shape = {},
       const std::vector<std::string> &component_labels = {},
       const std::string &associated = "",
       const refinement::RefinementFunctions_t ref_funcs_ =
+          refinement::RefinementFunctions_t::RegisterOps<
+              refinement_ops::ProlongateSharedMinMod, refinement_ops::RestrictAverage>(),
+      const refinement::RefinementFunctions_t flux_ref_funcs_ =
           refinement::RefinementFunctions_t::RegisterOps<
               refinement_ops::ProlongateSharedMinMod, refinement_ops::RestrictAverage>());
 
-  // 1 constructor
+  Metadata(
+      const std::vector<MetadataFlag> &bits, const std::vector<int> &shape = {},
+      const std::vector<std::string> &component_labels = {},
+      const std::string &associated = "",
+      const refinement::RefinementFunctions_t ref_funcs_ =
+          refinement::RefinementFunctions_t::RegisterOps<
+              refinement_ops::ProlongateSharedMinMod, refinement_ops::RestrictAverage>())
+      : Metadata(bits, {}, shape, component_labels, associated, ref_funcs_, ref_funcs_) {}
+
   Metadata(const std::vector<MetadataFlag> &bits, const std::vector<int> &shape,
            const std::string &associated)
       : Metadata(bits, shape, {}, associated) {}
 
-  // 2 constructors
   Metadata(const std::vector<MetadataFlag> &bits,
            const std::vector<std::string> component_labels,
            const std::string &associated = "")
       : Metadata(bits, {1}, component_labels, associated) {}
 
-  // 1 constructor
   Metadata(const std::vector<MetadataFlag> &bits, const std::string &associated)
       : Metadata(bits, {1}, {}, associated) {}
 
+  std::shared_ptr<Metadata> GetSPtrFluxMetadata() {
+    PARTHENON_REQUIRE(IsSet(WithFluxes),
+                      "Asking for flux metadata from metadata that doesn't have it.");
+    return flux_metadata;
+  }
+
+ private:
+  std::shared_ptr<Metadata> flux_metadata;
+
+ public:
   // Static routines
   static MetadataFlag AddUserFlag(const std::string &name);
   static bool FlagNameExists(const std::string &flagname);
@@ -536,6 +559,12 @@ class Metadata {
     refinement_funcs_ =
         refinement::RefinementFunctions_t::RegisterOps<ProlongationOp, RestrictionOp,
                                                        InternalProlongationOp>();
+    // Propagate refinement operations to flux metadata for backward compatibility
+    if (IsSet(WithFluxes)) {
+      flux_metadata->refinement_funcs_ =
+          refinement::RefinementFunctions_t::RegisterOps<ProlongationOp, RestrictionOp,
+                                                         InternalProlongationOp>();
+    }
   }
 
   // Operators
@@ -668,7 +697,6 @@ Set_t GetByFlag(const Metadata::FlagCollection &flags, NameMap_t &nameMap,
   return out;
 }
 } // namespace MetadataUtils
-
 } // namespace parthenon
 
 #endif // INTERFACE_METADATA_HPP_

src/interface/sparse_pack_base.cpp

@@ -67,19 +67,16 @@ SparsePackBase::GetAllocStatus(T *pmd, const PackDescriptor &desc,
                                const std::vector<bool> &include_block) {
   using mbd_t = MeshBlockData<Real>;
 
-  int nvar = desc.nvar_groups;
-
-  std::vector<int> astat;
+  const int nvar = desc.nvar_groups;
+  const int nblock = pmd->NumBlocks();
+  std::vector<int> astat(nblock * desc.nvar_tot);
+  int idx = 0;
   ForEachBlock(pmd, include_block, [&](int b, mbd_t *pmbd) {
     const auto &uid_map = pmbd->GetUidMap();
     for (int i = 0; i < nvar; ++i) {
       for (const auto &[var_name, uid] : desc.var_groups[i]) {
-        if (uid_map.count(uid) > 0) {
-          const auto pv = uid_map.at(uid);
-          astat.push_back(pv->GetAllocationStatus());
-        } else {
-          astat.push_back(-1);
-        }
+        astat[idx++] =
+            uid_map.count(uid) > 0 ? (uid_map.at(uid))->GetAllocationStatus() : -1;
       }
     }
   });

src/interface/sparse_pack_base.hpp

@@ -138,7 +138,7 @@ 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), identifier("") {}
+        coarse(false), flat(false), identifier(""), nvar_tot(0) {}
 
   template <class GROUP_t, class SELECTOR_t>
   PackDescriptor(StateDescriptor *psd, const std::vector<GROUP_t> &var_groups_in,
@@ -147,7 +147,7 @@ struct PackDescriptor {
         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)),
-        identifier(GetIdentifier()) {
+        identifier(GetIdentifier()), nvar_tot(GetNVarsTotal(var_groups)) {
     PARTHENON_REQUIRE(!(with_fluxes && coarse),
                       "Probably shouldn't be making a coarse pack with fine fluxes.");
   }
@@ -159,8 +159,18 @@ struct PackDescriptor {
   const bool coarse;
   const bool flat;
   const std::string identifier;
+  const std::size_t nvar_tot;
 
  private:
+  static int GetNVarsTotal(const std::vector<VariableGroup_t> &var_groups) {
+    int nvar_tot = 0;
+    for (const auto &group : var_groups) {
+      for (const auto &[a, b] : group) {
+        nvar_tot++;
+      }
+    }
+    return nvar_tot;
+  }
   std::string GetIdentifier() {
     std::string ident("");
     for (const auto &vgroup : var_groups) {

src/interface/sparse_pool.cpp

@@ -11,6 +11,8 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+#include <memory>
+
 #include "interface/sparse_pool.hpp"
 
 #include "interface/metadata.hpp"
@@ -48,44 +50,58 @@ SparsePool::SparsePool(const std::string &base_name, const Metadata &metadata,
   }
 }
 
-const Metadata &SparsePool::AddImpl(int sparse_id, const std::vector<int> &shape,
-                                    const MetadataFlag *vector_tensor,
-                                    const std::vector<std::string> &component_labels) {
-  PARTHENON_REQUIRE_THROWS(sparse_id != InvalidSparseID,
-                           "Tried to add InvalidSparseID to sparse pool " + base_name_);
-
+std::shared_ptr<Metadata>
+MakeSparseVarMetadataImpl(Metadata *in, const std::vector<int> &shape,
+                          const MetadataFlag *vector_tensor,
+                          const std::vector<std::string> &component_labels) {
   // copy shared metadata
-  Metadata this_metadata(
-      shared_metadata_.Flags(), shape.size() > 0 ? shape : shared_metadata_.Shape(),
-      component_labels.size() > 0 ? component_labels
-                                  : shared_metadata_.getComponentLabels(),
-      shared_metadata_.getAssociated(), shared_metadata_.GetRefinementFunctions());
+  auto this_metadata = std::make_shared<Metadata>(
+      in->Flags(), shape.size() > 0 ? shape : in->Shape(),
+      component_labels.size() > 0 ? component_labels : in->getComponentLabels(),
+      in->getAssociated(), in->GetRefinementFunctions());
 
-  this_metadata.SetSparseThresholds(shared_metadata_.GetAllocationThreshold(),
-                                    shared_metadata_.GetDeallocationThreshold(),
-                                    shared_metadata_.GetDefaultValue());
+  this_metadata->SetSparseThresholds(in->GetAllocationThreshold(),
+                                     in->GetDeallocationThreshold(),
+                                     in->GetDefaultValue());
 
   // if vector_tensor is set, apply it
   if (vector_tensor != nullptr) {
     if (*vector_tensor == Metadata::Vector) {
-      this_metadata.Unset(Metadata::Tensor);
-      this_metadata.Set(Metadata::Vector);
+      this_metadata->Unset(Metadata::Tensor);
+      this_metadata->Set(Metadata::Vector);
     } else if (*vector_tensor == Metadata::Tensor) {
-      this_metadata.Unset(Metadata::Vector);
-      this_metadata.Set(Metadata::Tensor);
+      this_metadata->Unset(Metadata::Vector);
+      this_metadata->Set(Metadata::Tensor);
     } else if (*vector_tensor == Metadata::None) {
-      this_metadata.Unset(Metadata::Vector);
-      this_metadata.Unset(Metadata::Tensor);
+      this_metadata->Unset(Metadata::Vector);
+      this_metadata->Unset(Metadata::Tensor);
     } else {
       PARTHENON_THROW("Expected MetadataFlag Vector, Tensor, or None, but got " +
                       vector_tensor->Name());
     }
   }
 
   // just in case
-  this_metadata.IsValid(true);
+  this_metadata->IsValid(true);
+
+  return this_metadata;
+}
+
+const Metadata &SparsePool::AddImpl(int sparse_id, const std::vector<int> &shape,
+                                    const MetadataFlag *vector_tensor,
+                                    const std::vector<std::string> &component_labels) {
+  PARTHENON_REQUIRE_THROWS(sparse_id != InvalidSparseID,
+                           "Tried to add InvalidSparseID to sparse pool " + base_name_);
+
+  auto this_metadata = MakeSparseVarMetadataImpl(&shared_metadata_, shape, vector_tensor,
+                                                 component_labels);
+  if (this_metadata->IsSet(Metadata::WithFluxes)) {
+    this_metadata->GetSPtrFluxMetadata() =
+        MakeSparseVarMetadataImpl(shared_metadata_.GetSPtrFluxMetadata().get(), shape,
+                                  vector_tensor, component_labels);
+  }
 
-  const auto ins = pool_.insert({sparse_id, this_metadata});
+  const auto ins = pool_.insert({sparse_id, *this_metadata});
   PARTHENON_REQUIRE_THROWS(ins.second, "Tried to add sparse ID " +
                                            std::to_string(sparse_id) +
                                            " to sparse pool '" + base_name_ +