Draft: Test and fix reflective BCs

src/include/add_ghosts_reflecting.hxx

@@ -0,0 +1,71 @@
+#pragma once
+
+#include "../kg/include/kg/Vec3.h"
+#include "grid.hxx"
+
+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);
+  }
+}
+
+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)
+{
+  // FIXME only need to scan 1 cell into the ghost region for 1st-order moments
+
+  Int3 idx_end = ldims - ib;
+  Int3 idx_begin = 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;
+        }
+      }
+    }
+  }
+}

src/include/fields_item.hxx

@@ -7,6 +7,7 @@
 #include "fields3d.hxx"
 #include "particles.hxx"
 #include "psc_fields_c.h"
+#include "add_ghosts_reflecting.hxx"
 
 #include <mrc_profile.h>
 
@@ -53,74 +54,6 @@ 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);
-    }
-  }
-
-  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)
-  {
-    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);
-            }
-          }
-        }
-      }
-    } 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);
-            }
-          }
-        }
-      }
-    } else {
-      assert(0);
-    }
-  }
-
   template <typename FE>
   void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, const Int3& ib,
                            int p, int mb, int me)

src/libpsc/psc_checks/checks_impl.hxx

@@ -209,11 +209,10 @@ struct checks_order_2nd
   using Moment_rho_nc = Moment_rho_2nd_nc<S, D>;
 };
 
-template <typename MP, typename S, typename ITEM_RHO>
+template <typename S, typename ITEM_RHO>
 class ChecksCommon
 {
 public:
-  using Mparticles = MP;
   using storage_type = S;
   using item_rho_type = ITEM_RHO;
 
@@ -228,5 +227,5 @@ public:
 
 template <typename MP, typename MF, typename ORDER, typename D>
 using Checks_ =
-  ChecksCommon<MP, typename MF::Storage,
+  ChecksCommon<typename MF::Storage,
                typename ORDER::template Moment_rho_nc<typename MF::Storage, D>>;

src/libpsc/tests/CMakeLists.txt

@@ -21,6 +21,7 @@ add_psc_test(test_rng)
 add_psc_test(test_mparticles_cuda)
 add_psc_test(test_mparticles)
 add_psc_test(test_output_particles)
+add_psc_test(test_reflective_bcs)
 add_psc_test(test_mfields)
 add_psc_test(test_mfields_cuda)
 add_psc_test(test_bnd)