Small: Fix bug in sparse packs that include fluxes

CHANGELOG.md

@@ -21,6 +21,7 @@
 - [[PR 1004]](https://github.com/parthenon-hpc-lab/parthenon/pull/1004) Allow parameter modification from an input file for restarts
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 1088]](https://github.com/parthenon-hpc-lab/parthenon/pull/1088) Correctly fill fluxes for non-cell variables in SparsePacks
 - [[PR 1083]](https://github.com/parthenon-hpc-lab/parthenon/pull/1083) Correctly fill VariableFluxPack for edge fluxes in 2D
 - [[PR 1087]](https://github.com/parthenon-hpc-lab/parthenon/pull/1087) Make sure InnerLoopPatternTVR is resolved on device properly when it is the default loop pattern
 - [[PR 1071]](https://github.com/parthenon-hpc-lab/parthenon/pull/1070) Fix bug in static mesh refinement related to redefinition of Mesh::root_level

src/interface/sparse_pack.hpp

@@ -337,22 +337,22 @@ class SparsePack : public SparsePackBase {
   auto &flux(const int b, const int dir, const int idx) const {
     PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
     PARTHENON_DEBUG_REQUIRE(dir > 0 && dir < 4 && with_fluxes_, "Bad input to flux call");
-    return pack_(dir, b, idx);
+    return pack_(dir - 1 + flx_idx_, b, idx);
   }
 
   KOKKOS_INLINE_FUNCTION
   Real &flux(const int b, const int dir, const int idx, const int k, const int j,
              const int i) const {
     PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
     PARTHENON_DEBUG_REQUIRE(dir > 0 && dir < 4 && with_fluxes_, "Bad input to flux call");
-    return pack_(dir, b, idx)(k, j, i);
+    return pack_(dir - 1 + flx_idx_, b, idx)(k, j, i);
   }
 
   KOKKOS_INLINE_FUNCTION
   Real &flux(const int dir, const int idx, const int k, const int j, const int i) const {
     PARTHENON_DEBUG_REQUIRE(flat_, "Accessor must only be used for flat packs");
     PARTHENON_DEBUG_REQUIRE(dir > 0 && dir < 4 && with_fluxes_, "Bad input to flux call");
-    return pack_(dir, 0, idx)(k, j, i);
+    return pack_(dir - 1 + flx_idx_, 0, idx)(k, j, i);
   }
 
   KOKKOS_INLINE_FUNCTION
@@ -362,7 +362,7 @@ class SparsePack : public SparsePackBase {
     PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
     PARTHENON_DEBUG_REQUIRE(dir > 0 && dir < 4 && with_fluxes_, "Bad input to flux call");
     const int n = bounds_(0, b, idx.VariableIdx()) + idx.Offset();
-    return pack_(dir, b, n)(k, j, i);
+    return pack_(dir - 1 + flx_idx_, b, n)(k, j, i);
   }
 
   template <class TIn, REQUIRES(IncludesType<TIn, Ts...>::value)>
@@ -371,7 +371,7 @@ class SparsePack : public SparsePackBase {
     PARTHENON_DEBUG_REQUIRE(!flat_, "Accessor cannot be used for flat packs");
     PARTHENON_DEBUG_REQUIRE(dir > 0 && dir < 4 && with_fluxes_, "Bad input to flux call");
     const int vidx = GetLowerBound(b, t) + t.idx;
-    return pack_(dir, b, vidx)(k, j, i);
+    return pack_(dir - 1 + flx_idx_, b, vidx)(k, j, i);
   }
 
   template <class... VTs>

src/interface/sparse_pack_base.cpp

@@ -110,6 +110,7 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc,
   int max_size = 0;
   int nblocks = 0;
   bool contains_face_or_edge = false;
+  bool contains_face_with_fluxes = false;
   int size = 0; // local var used to compute size/block
   ForEachBlock(pmd, include_block, [&](int b, mbd_t *pmbd) {
     if (!desc.flat) {
@@ -122,8 +123,13 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc,
         if (uid_map.count(uid) > 0) {
           const auto pv = uid_map.at(uid);
           if (pv->IsAllocated()) {
-            if (pv->IsSet(Metadata::Face) || pv->IsSet(Metadata::Edge))
+            if (pv->IsSet(Metadata::Edge)) contains_face_or_edge = true;
+            if (pv->IsSet(Metadata::Face)) {
+              if (pv->IsSet(Metadata::WithFluxes) && desc.with_fluxes) {
+                contains_face_with_fluxes = true;
+              }
               contains_face_or_edge = true;
+            }
             int prod = pv->GetDim(6) * pv->GetDim(5) * pv->GetDim(4);
             size += prod;       // max size/block (or total size for flat)
             pack.size_ += prod; // total ragged size
@@ -138,7 +144,11 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc,
 
   // Allocate the views
   int leading_dim = 1;
-  if (desc.with_fluxes) {
+  pack.flx_idx_ = 1;
+  if (contains_face_with_fluxes) {
+    leading_dim += 5;
+    pack.flx_idx_ = 3;
+  } else if (desc.with_fluxes) {
     leading_dim += 3;
   } else if (contains_face_or_edge) {
     leading_dim += 2;
@@ -205,11 +215,11 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc,
 
                     if (pv->IsSet(Metadata::Face)) {
                       pack.pack_h_(0, b, idx).topological_element =
-                          TopologicalElement::E1;
+                          TopologicalElement::F1;
                       pack.pack_h_(1, b, idx).topological_element =
-                          TopologicalElement::E2;
+                          TopologicalElement::F2;
                       pack.pack_h_(2, b, idx).topological_element =
-                          TopologicalElement::E3;
+                          TopologicalElement::F3;
                     }
 
                   } else { // This is a cell, node, or a variable that doesn't have
@@ -221,13 +231,16 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc,
                     }
                     if (pv->IsSet(Metadata::Vector))
                       pack.pack_h_(0, b, idx).vector_component = v + 1;
+                  }
 
-                    if (desc.with_fluxes && pv->IsSet(Metadata::WithFluxes)) {
-                      pack.pack_h_(1, b, idx) = pvf->data.Get(0, t, u, v);
-                      pack.pack_h_(2, b, idx) = pvf->data.Get(1, t, u, v);
-                      pack.pack_h_(3, b, idx) = pvf->data.Get(2, t, u, v);
+                  if (desc.with_fluxes && pv->IsSet(Metadata::WithFluxes)) {
+                    pack.pack_h_(0 + pack.flx_idx_, b, idx) = pvf->data.Get(0, t, u, v);
+                    if (!pv->IsSet(Metadata::Edge)) {
+                      pack.pack_h_(1 + pack.flx_idx_, b, idx) = pvf->data.Get(1, t, u, v);
+                      pack.pack_h_(2 + pack.flx_idx_, b, idx) = pvf->data.Get(2, t, u, v);
                     }
                   }
+
                   for (auto el :
                        GetTopologicalElements(pack.pack_h_(0, b, idx).topological_type)) {
                     pack.pack_h_(static_cast<int>(el) % 3, b, idx).topological_element =

src/interface/sparse_pack_base.hpp

@@ -92,6 +92,7 @@ class SparsePackBase {
   bounds_h_t bounds_h_;
   coords_t coords_;
 
+  int flx_idx_;
   bool with_fluxes_;
   bool coarse_;
   bool flat_;

src/utils/interpolation.hpp

@@ -6,18 +6,16 @@
 // Interpolation copied/refactored from
 // https://github.com/lanl/phoebus and https://github.com/lanl/spiner
 //========================================================================================
-// © 2022. 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 for the U.S.
-// Department of Energy/National Nuclear Security Administration. All
-// rights in the program are reserved by Triad National Security, LLC,
-// and the U.S. Department of Energy/National Nuclear Security
-// Administration. The Government is granted for itself and others
-// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
-// 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.
+// (C) (or copyright) 2024. 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
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// 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.
 
 #ifndef UTILS_INTERPOLATION_HPP_
 #define UTILS_INTERPOLATION_HPP_

src/utils/robust.hpp

@@ -5,18 +5,16 @@
 //========================================================================================
 // Copied from https://github.com/lanl/phoebus
 //========================================================================================
-// © 2022. 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 for the U.S.
-// Department of Energy/National Nuclear Security Administration. All
-// rights in the program are reserved by Triad National Security, LLC,
-// and the U.S. Department of Energy/National Nuclear Security
-// Administration. The Government is granted for itself and others
-// acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
-// 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.
+// (C) (or copyright) 2024. 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
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// 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.
 
 #ifndef UTILS_ROBUST_HPP_
 #define UTILS_ROBUST_HPP_