Fix reflective BC moment calculations

src/include/add_ghosts_reflecting.hxx

@@ -0,0 +1,150 @@
+#pragma once
+
+#include "../kg/include/kg/Vec3.h"
+#include "grid.hxx"
+
+template <typename FE>
+void add_ghosts_reflecting_lo_cc(const Int3& ldims, FE& mres_gt, const Int3& ib,
+                                 int p, int d, int mb, int me)
+{
+  // FIXME only need to scan 1 cell into the ghost region for 1st-order moments
+
+  Int3 idx_begin = ib;
+  Int3 idx_end = ldims - ib;
+
+  idx_begin[d] = 0;
+  idx_end[d] = -ib[d];
+
+  Int3 idx;
+
+  for (int m = mb; m < me; m++) {
+    for (idx[2] = idx_begin[2]; idx[2] < idx_end[2]; idx[2]++) {
+      for (idx[1] = idx_begin[1]; idx[1] < idx_end[1]; idx[1]++) {
+        for (idx[0] = idx_begin[0]; idx[0] < idx_end[0]; idx[0]++) {
+          Int3 idx_reflected = idx;
+          idx_reflected[d] = -idx[d] - 1;
+
+          auto val_reflected =
+            mres_gt(idx_reflected[0] - ib[0], idx_reflected[1] - ib[1],
+                    idx_reflected[2] - ib[2], m, p);
+
+          mres_gt(idx[0] - ib[0], idx[1] - ib[1], idx[2] - ib[2], m, p) +=
+            val_reflected;
+        }
+      }
+    }
+  }
+}
+
+template <typename FE>
+void add_ghosts_reflecting_hi_cc(const Int3& ldims, FE& mres_gt, const Int3& ib,
+                                 int p, int d, int mb, int me)
+{
+  // FIXME only need to scan 1 cell into the ghost region for 1st-order moments
+
+  Int3 idx_begin = ib;
+  Int3 idx_end = ldims - ib;
+
+  idx_begin[d] = ldims[d] + ib[d];
+  idx_end[d] = ldims[d];
+
+  Int3 idx;
+
+  for (int m = mb; m < me; m++) {
+    for (idx[2] = idx_begin[2]; idx[2] < idx_end[2]; idx[2]++) {
+      for (idx[1] = idx_begin[1]; idx[1] < idx_end[1]; idx[1]++) {
+        for (idx[0] = idx_begin[0]; idx[0] < idx_end[0]; idx[0]++) {
+          Int3 idx_reflected = idx;
+          idx_reflected[d] = 2 * ldims[d] - idx[d] - 1;
+
+          auto val_reflected =
+            mres_gt(idx_reflected[0] - ib[0], idx_reflected[1] - ib[1],
+                    idx_reflected[2] - ib[2], m, p);
+
+          mres_gt(idx[0] - ib[0], idx[1] - ib[1], idx[2] - ib[2], m, p) +=
+            val_reflected;
+        }
+      }
+    }
+  }
+}
+
+template <typename FE>
+void add_ghosts_reflecting_lo_nc(const Int3& ldims, FE& mres_gt, const Int3& ib,
+                                 int p, int d, int mb, int me)
+{
+  // When lower boundary has 2 ghosts, upper boundary of nc grid only has 1
+  // because nc has +1 point compared to cc, and takes it from upper ghost
+  // region. Thus, ignore a ghost from lower ghost region to balance it out.
+
+  Int3 unused_ghost = {!!ib[0], !!ib[1], !!ib[2]};
+
+  Int3 idx_begin = ib + unused_ghost;
+  Int3 idx_end = ldims - ib;
+
+  // Also note that the boundary nc point is the point about which reflections
+  // occur, and isn't itself reflected.
+
+  idx_begin[d] = 1;
+  idx_end[d] = 1 - ib[d] - unused_ghost[d];
+
+  Int3 idx;
+
+  for (int m = mb; m < me; m++) {
+    for (idx[2] = idx_begin[2]; idx[2] < idx_end[2]; idx[2]++) {
+      for (idx[1] = idx_begin[1]; idx[1] < idx_end[1]; idx[1]++) {
+        for (idx[0] = idx_begin[0]; idx[0] < idx_end[0]; idx[0]++) {
+          Int3 idx_reflected = idx;
+          idx_reflected[d] = -idx[d];
+
+          auto val_reflected =
+            mres_gt(idx_reflected[0] - ib[0], idx_reflected[1] - ib[1],
+                    idx_reflected[2] - ib[2], m, p);
+
+          mres_gt(idx[0] - ib[0], idx[1] - ib[1], idx[2] - ib[2], m, p) +=
+            val_reflected;
+        }
+      }
+    }
+  }
+}
+
+template <typename FE>
+void add_ghosts_reflecting_hi_nc(const Int3& ldims, FE& mres_gt, const Int3& ib,
+                                 int p, int d, int mb, int me)
+{
+  // When lower boundary has 2 ghosts, upper boundary of nc grid only has 1
+  // because nc has +1 point compared to cc, and takes it from upper ghost
+  // region. Thus, ignore a ghost from lower ghost region to balance it out.
+
+  Int3 unused_ghost = {!!ib[0], !!ib[1], !!ib[2]};
+
+  Int3 idx_begin = ib + unused_ghost;
+  Int3 idx_end = ldims - ib;
+
+  // Also note that the extra nc point is the point about which reflections
+  // occur, and isn't itself reflected.
+
+  idx_begin[d] = ldims[d] + ib[d] + unused_ghost[d];
+  idx_end[d] = ldims[d];
+
+  Int3 idx;
+
+  for (int m = mb; m < me; m++) {
+    for (idx[2] = idx_begin[2]; idx[2] < idx_end[2]; idx[2]++) {
+      for (idx[1] = idx_begin[1]; idx[1] < idx_end[1]; idx[1]++) {
+        for (idx[0] = idx_begin[0]; idx[0] < idx_end[0]; idx[0]++) {
+          Int3 idx_reflected = idx;
+          idx_reflected[d] = 2 * ldims[d] - idx[d];
+
+          auto val_reflected =
+            mres_gt(idx_reflected[0] - ib[0], idx_reflected[1] - ib[1],
+                    idx_reflected[2] - ib[2], m, p);
+
+          mres_gt(idx[0] - ib[0], idx[1] - ib[1], idx[2] - ib[2], m, p) +=
+            val_reflected;
+        }
+      }
+    }
+  }
+}

src/include/fields_item.hxx

@@ -3,10 +3,12 @@
 
 #include "../libpsc/psc_bnd/psc_bnd_impl.hxx"
 #include "bnd.hxx"
+#include "centering.hxx"
 #include "fields.hxx"
 #include "fields3d.hxx"
 #include "particles.hxx"
 #include "psc_fields_c.h"
+#include "add_ghosts_reflecting.hxx"
 
 #include <mrc_profile.h>
 
@@ -32,118 +34,89 @@ inline std::vector<std::string> addKindSuffix(
 // ======================================================================
 // ItemMomentBnd
 
-template <typename S, typename Bnd>
-class ItemMomentBnd
+template <centering::Centering C>
+class DomainBoundary
 {
 public:
-  using storage_type = S;
-
-  ItemMomentBnd(const Grid_t& grid) {}
-
-  void add_ghosts(const Grid_t& grid, storage_type& mres_gt, const Int3& ib)
-  {
-    for (int p = 0; p < mres_gt.shape(4); p++) {
-      add_ghosts_boundary(grid, mres_gt, ib, p, 0, mres_gt.shape(3));
-    }
-
-    bnd_.add_ghosts(grid, mres_gt, ib, 0, mres_gt.shape(3));
-  }
-
-private:
-  // ----------------------------------------------------------------------
-  // boundary stuff FIXME, should go elsewhere...
-
   template <typename FE>
-  void add_ghosts_reflecting_lo(const Int3& ldims, FE& mres_gt, const Int3& ib,
-                                int p, int d, int mb, int me)
-  {
-    int bx = ldims[0] == 1 ? 0 : 1;
-    if (d == 1) {
-      for (int iz = -1; iz < ldims[2] + 1; iz++) {
-        for (int ix = -bx; ix < ldims[0] + bx; ix++) {
-          int iy = 0;
-          {
-            for (int m = mb; m < me; m++) {
-              mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) +=
-                mres_gt(ix - ib[0], iy - 1 - ib[1], iz - ib[2], m, p);
-            }
-          }
-        }
-      }
-    } else if (d == 2) {
-      for (int iy = 0 * -1; iy < ldims[1] + 0 * 1; iy++) {
-        for (int ix = -bx; ix < ldims[0] + bx; ix++) {
-          int iz = 0;
-          {
-            for (int m = mb; m < me; m++) {
-              mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) +=
-                mres_gt(ix - ib[0], iy - ib[1], iz - 1 - ib[2], m, p);
-            }
-          }
-        }
-      }
-    } else {
-      assert(0);
-    }
-  }
+  static void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt,
+                                  const Int3& ib, int p, int mb, int me);
+};
 
+template <>
+class DomainBoundary<centering::CC>
+{
+public:
   template <typename FE>
-  void add_ghosts_reflecting_hi(const Int3& ldims, FE& mres_gt, const Int3& ib,
-                                int p, int d, int mb, int me)
+  static void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt,
+                                  const Int3& ib, int p, int mb, int me)
   {
-    int bx = ldims[0] == 1 ? 0 : 1;
-    if (d == 1) {
-      for (int iz = -1; iz < ldims[2] + 1; iz++) {
-        for (int ix = -bx; ix < ldims[0] + bx; ix++) {
-          int iy = ldims[1] - 1;
-          {
-            for (int m = mb; m < me; m++) {
-              mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) +=
-                mres_gt(ix - ib[0], iy - 1 + ib[1], iz - ib[2], m, p);
-            }
-          }
-        }
+    // lo
+    for (int d = 0; d < 3; d++) {
+      // FIXME why reflect for open BCs?
+      if (grid.atBoundaryLo(p, d) && (grid.bc.prt_lo[d] == BND_PRT_REFLECTING ||
+                                      grid.bc.prt_lo[d] == BND_PRT_OPEN)) {
+        add_ghosts_reflecting_lo_cc(grid.ldims, mres_gt, ib, p, d, mb, me);
       }
-    } else if (d == 2) {
-      for (int iy = 0 * -1; iy < ldims[1] + 0 * 1; iy++) {
-        for (int ix = -bx; ix < ldims[0] + bx; ix++) {
-          int iz = ldims[2] - 1;
-          {
-            for (int m = mb; m < me; m++) {
-              mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) +=
-                mres_gt(ix - ib[0], iy - ib[1], iz + 1 - ib[2], m, p);
-            }
-          }
-        }
+    }
+    // hi
+    for (int d = 0; d < 3; d++) {
+      // FIXME why reflect for open BCs?
+      if (grid.atBoundaryHi(p, d) && (grid.bc.prt_hi[d] == BND_PRT_REFLECTING ||
+                                      grid.bc.prt_hi[d] == BND_PRT_OPEN)) {
+        add_ghosts_reflecting_hi_cc(grid.ldims, mres_gt, ib, p, d, mb, me);
       }
-    } else {
-      assert(0);
     }
   }
+};
 
+template <>
+class DomainBoundary<centering::NC>
+{
+public:
   template <typename FE>
-  void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, const Int3& ib,
-                           int p, int mb, int me)
+  static void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt,
+                                  const Int3& ib, int p, int mb, int me)
   {
     // lo
     for (int d = 0; d < 3; d++) {
-      if (grid.atBoundaryLo(p, d)) {
-        if (grid.bc.prt_lo[d] == BND_PRT_REFLECTING ||
-            grid.bc.prt_lo[d] == BND_PRT_OPEN) {
-          add_ghosts_reflecting_lo(grid.ldims, mres_gt, ib, p, d, mb, me);
-        }
+      // FIXME why reflect for open BCs?
+      if (grid.atBoundaryLo(p, d) && (grid.bc.prt_lo[d] == BND_PRT_REFLECTING ||
+                                      grid.bc.prt_lo[d] == BND_PRT_OPEN)) {
+        add_ghosts_reflecting_lo_nc(grid.ldims, mres_gt, ib, p, d, mb, me);
       }
     }
     // hi
     for (int d = 0; d < 3; d++) {
-      if (grid.atBoundaryHi(p, d)) {
-        if (grid.bc.prt_hi[d] == BND_PRT_REFLECTING ||
-            grid.bc.prt_hi[d] == BND_PRT_OPEN) {
-          add_ghosts_reflecting_hi(grid.ldims, mres_gt, ib, p, d, mb, me);
-        }
+      // FIXME why reflect for open BCs?
+      if (grid.atBoundaryHi(p, d) && (grid.bc.prt_hi[d] == BND_PRT_REFLECTING ||
+                                      grid.bc.prt_hi[d] == BND_PRT_OPEN)) {
+        add_ghosts_reflecting_hi_nc(grid.ldims, mres_gt, ib, p, d, mb, me);
       }
     }
   }
+};
+
+// ======================================================================
+// ItemMomentBnd
+
+template <typename S, typename Bnd, centering::Centering C>
+class ItemMomentBnd
+{
+public:
+  using storage_type = S;
+
+  ItemMomentBnd(const Grid_t& grid) {}
+
+  void add_ghosts(const Grid_t& grid, storage_type& mres_gt, const Int3& ib)
+  {
+    for (int p = 0; p < mres_gt.shape(4); p++) {
+      DomainBoundary<C>::add_ghosts_boundary(grid, mres_gt, ib, p, 0,
+                                             mres_gt.shape(3));
+    }
+
+    bnd_.add_ghosts(grid, mres_gt, ib, 0, mres_gt.shape(3));
+  }
 
 private:
   Bnd bnd_;
@@ -189,6 +162,6 @@ public:
   }
 
 private:
-  ItemMomentBnd<storage_type, Bnd> bnd_;
+  ItemMomentBnd<storage_type, Bnd, moment_type::CENTERING> bnd_;
   std::vector<std::string> comp_names_;
 };