Merge pull request #1083 from acreyes/acreyes/edge-flux-pack

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 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
 - [[PR 1073]](https://github.com/parthenon-hpc-lab/parthenon/pull/1073) Fix bug in AMR and sparse restarts

src/interface/variable_pack.hpp

@@ -452,9 +452,14 @@ class VariableFluxPack : public VariablePack<T> {
   // host only
   inline auto flux_alloc_status() const { return flux_alloc_status_; }
 
-  KOKKOS_FORCEINLINE_FUNCTION
-  const ViewOfParArrays<T> &flux(const int dir) const {
-    assert(dir > 0 && dir <= this->ndim_);
+  template <TopologicalType TT = TopologicalType::Face>
+  KOKKOS_FORCEINLINE_FUNCTION const ViewOfParArrays<T> &flux(const int dir) const {
+    assert(dir > 0);
+    if constexpr (TT == TopologicalType::Edge) {
+      assert(dir <= this->ndim_ || (this->ndim_ == 2 && dir == 3));
+    } else {
+      assert(dir <= this->ndim_);
+    }
     return f_[dir - 1];
   }
 
@@ -471,10 +476,11 @@ class VariableFluxPack : public VariablePack<T> {
   constexpr bool IsFluxAllocated(const int /*n*/) const { return true; }
 #endif
 
-  KOKKOS_FORCEINLINE_FUNCTION
-  T &flux(const int dir, const int n, const int k, const int j, const int i) const {
+  template <TopologicalType TT = TopologicalType::Face>
+  KOKKOS_FORCEINLINE_FUNCTION T &flux(const int dir, const int n, const int k,
+                                      const int j, const int i) const {
     assert(IsFluxAllocated(n));
-    return flux(dir)(n)(k, j, i);
+    return flux<TT>(dir)(n)(k, j, i);
   }
 
  private:
@@ -681,9 +687,15 @@ void FillFluxViews(const VariableVector<T> &vars, const int ndim,
         for (int i = 0; i < v->GetDim(4); i++) {
           host_al(vindex) = v->IsAllocated();
           if (v->IsAllocated()) {
-            host_f1(vindex) = v->data.Get(0, k, j, i);
-            if (ndim >= 2) host_f2(vindex) = v->data.Get(1, k, j, i);
-            if (ndim >= 3) host_f3(vindex) = v->data.Get(2, k, j, i);
+            if (v->IsSet(Metadata::Edge)) {
+              if (ndim >= 2) host_f3(vindex) = v->data.Get(2, k, j, i);
+              if (ndim >= 3) host_f2(vindex) = v->data.Get(1, k, j, i);
+              if (ndim >= 3) host_f1(vindex) = v->data.Get(0, k, j, i);
+            } else {
+              host_f1(vindex) = v->data.Get(0, k, j, i);
+              if (ndim >= 2) host_f2(vindex) = v->data.Get(1, k, j, i);
+              if (ndim >= 3) host_f3(vindex) = v->data.Get(2, k, j, i);
+            }
           }
 
           vindex++;