Merge pull request #390 from JamesMcClung/pr/marder-at-edges

Fix Marder correction at conductive and open edges

Author: JamesMcClung

src/include/boundary_injector.hxx

@@ -127,12 +127,6 @@ public:
           get_n_in_cell(1.0, prts_per_unit_density_, true);
 
         for (int prt_count = 0; prt_count < n_prts_to_try_inject; prt_count++) {
-          // FIXME #948112531345 (also see the other FIXMEs with this id)
-          // Depositing current in ghost corners leads to false-positive gauss
-          // errors. This could be avoided here by artificially constraining the
-          // trajectories of injected particles. However, assuming the particle
-          // boundary condition is "open", outflowing particles cause the same
-          // problem, and there's nothing to be done about that.
           psc::particle::Inject prt =
             particle_generator_.get(cell_corner, grid.domain.dx);
 
@@ -151,9 +145,6 @@ public:
           // Taken from push_particles_1vb.hxx PushParticlesVb::push_mprts()
 
           Real3 initial_normalized_pos = initial_x * dxi;
-          // pretend it came from the edge to inject proper current
-          initial_normalized_pos[INJECT_DIM_IDX_] = -1.0;
-
           Real3 final_normalized_pos = prt.x * dxi;
           Int3 final_idx = final_normalized_pos.fint();
 

src/libpsc/psc_checks/checks_impl.hxx

@@ -71,8 +71,23 @@ public:
     auto item_divj = Item_divj<MfieldsState>{};
 
     auto rho_p = psc::mflds::interior(grid, item_rho(mprts));
-    auto divj = psc::mflds::interior(grid, item_divj(mflds));
     auto d_rho = rho_p - rho_m_;
+
+    auto divj = psc::mflds::interior(grid, item_divj(mflds));
+
+    // account for insertion/deletion of incoming/outgoing particles
+    for (int p = 0; p < grid.n_patches(); p++) {
+      for (int d = 0; d < 3; d++) {
+        if (grid.atBoundaryLo(p, d) && grid.bc.fld_lo[d] == BND_FLD_OPEN) {
+          Int3 r = grid.ldims;
+          r[d] = 1;
+
+          divj.view(_s(0, r[0]), _s(0, r[1]), _s(0, r[2]), 0, p) =
+            -d_rho.view(_s(0, r[0]), _s(0, r[1]), _s(0, r[2]), 0, p) / grid.dt;
+        }
+      }
+    }
+
     auto dt_divj = grid.dt * divj;
 
     double local_err = gt::norm_linf(d_rho + dt_divj);
@@ -144,14 +159,12 @@ public:
     auto rho = psc::mflds::interior(grid, item_rho(mprts));
     auto dive = psc::mflds::interior(grid, item_dive(mflds));
 
+    // account for virtual charges implied by BCs
     for (int p = 0; p < grid.n_patches(); p++) {
       for (int d = 0; d < 3; d++) {
         if (grid.atBoundaryLo(p, d) &&
             (grid.bc.fld_lo[d] == BND_FLD_CONDUCTING_WALL ||
              grid.bc.fld_lo[d] == BND_FLD_OPEN)) {
-
-          // account for implicit surface charges
-
           Int3 r = grid.ldims;
           r[d] = 1;
 

src/libpsc/psc_push_fields/marder_impl.hxx

@@ -17,39 +17,6 @@ namespace psc
 namespace marder
 {
 
-namespace detail
-{
-
-inline void find_limits(const Grid_t& grid, int p, Int3& lx, Int3& rx, Int3& ly,
-                        Int3& ry, Int3& lz, Int3& rz)
-{
-  Int3 l_cc = {0, 0, 0}, r_cc = {0, 0, 0};
-  Int3 l_nc = {0, 0, 0}, r_nc = {0, 0, 0};
-  for (int d = 0; d < 3; d++) {
-    if (grid.bc.fld_lo[d] == BND_FLD_CONDUCTING_WALL &&
-        grid.atBoundaryLo(p, d)) {
-      l_cc[d] = -1;
-      l_nc[d] = -1;
-    }
-    if (grid.bc.fld_hi[d] == BND_FLD_CONDUCTING_WALL &&
-        grid.atBoundaryHi(p, d)) {
-      r_cc[d] = -1;
-      r_nc[d] = 0;
-    }
-  }
-  // FIXME, for conducting wall the signs here need checking...
-  lx = -Int3{l_cc[0], l_nc[1], l_nc[2]} + grid.ibn;
-  rx = Int3{r_cc[0], r_nc[1], r_nc[2]} + grid.ldims + grid.ibn;
-
-  ly = -Int3{l_nc[0], l_cc[1], l_nc[2]} + grid.ibn;
-  ry = Int3{r_nc[0], r_cc[1], r_nc[2]} + grid.ldims + grid.ibn;
-
-  lz = -Int3{l_nc[0], l_nc[1], l_cc[2]} + grid.ibn;
-  rz = Int3{r_nc[0], r_nc[1], r_cc[2]} + grid.ldims + grid.ibn;
-}
-
-} // namespace detail
-
 // ----------------------------------------------------------------------
 // correct
 //
@@ -70,53 +37,46 @@ inline void correct(const Grid_t& grid, E1& efield, const Int3& efield_ib,
   Real3 fac = .5f * real_t(grid.dt) * diffusion * Real3(grid.domain.dx_inv);
 
   for (int p = 0; p < grid.n_patches(); p++) {
-    Int3 lx, rx, ly, ry, lz, rz;
-    detail::find_limits(grid, p, lx, rx, ly, ry, lz, rz);
-
-    Int3 ls[3] = {lx, ly, lz};
-    Int3 rs[3] = {rx, ry, rz};
-
     auto res = mf.view(_all, _all, _all, 0, p);
     for (int d = 0; d < 3; d++) {
       if (grid.isInvar(d)) {
         continue;
       }
 
-      Int3 l = ls[d];
-      Int3 r = rs[d];
       auto e_comp = efield.view(_all, _all, _all, d, p);
 
-      gt::gslice s1x = _s(l[0], r[0]);
-      gt::gslice s1y = _s(l[1], r[1]);
-      gt::gslice s1z = _s(l[2], r[2]);
-
+      Int3 l = grid.ibn;
+      Int3 r = grid.ldims + grid.ibn;
+      gt::gslice s1[3] = {_s(l[0], r[0]), _s(l[1], r[1]), _s(l[2], r[2])};
       gt::gslice s2[3] = {_s(l[0], r[0]), _s(l[1], r[1]), _s(l[2], r[2])};
-      s2[d] = _s(l[d] + 1, r[d] + 1);
+      s2[d].start += 1;
+      s2[d].stop += 1;
 
-      e_comp.view(s1x, s1y, s1z) =
-        e_comp.view(s1x, s1y, s1z) +
-        (res.view(s2[0], s2[1], s2[2]) - res.view(s1x, s1y, s1z)) * fac[d];
+      e_comp.view(s1[0], s1[1], s1[2]) =
+        e_comp.view(s1[0], s1[1], s1[2]) +
+        (res.view(s2[0], s2[1], s2[2]) - res.view(s1[0], s1[1], s1[2])) *
+          fac[d];
     }
   }
 }
 
 #ifdef USE_CUDA
 
 template <typename E1, typename E2>
-inline void cuda_marder_correct_yz(E1& efield, E2& res, Float3 fac, Int3 ly,
-                                   Int3 ry, Int3 lz, Int3 rz)
+inline void cuda_marder_correct_yz(E1& efield, E2& res, Float3 fac, Int3 l,
+                                   Int3 r, Int3 l, Int3 r)
 {
   auto k_efield = efield.to_kernel();
   auto k_res = res.to_kernel();
   gt::launch<2>(
     {k_efield.shape(1), k_efield.shape(2)}, GT_LAMBDA(int iy, int iz) {
-      if ((iy >= ly[1] && iy < ry[1]) && (iz >= ly[2] && iz < ry[2])) {
+      if ((iy >= l[1] && iy < r[1]) && (iz >= l[2] && iz < r[2])) {
         k_efield(0, iy, iz, 1) =
           k_efield(0, iy, iz, 1) +
           fac[1] * (k_res(0, iy + 1, iz) - k_res(0, iy, iz));
       }
 
-      if ((iy >= lz[1] && iy < rz[1]) && (iz >= lz[2] && iz < rz[2])) {
+      if ((iy >= l[1] && iy < r[1]) && (iz >= l[2] && iz < r[2])) {
         k_efield(0, iy, iz, 2) =
           k_efield(0, iy, iz, 2) +
           fac[2] * (k_res(0, iy, iz + 1) - k_res(0, iy, iz));
@@ -126,30 +86,30 @@ inline void cuda_marder_correct_yz(E1& efield, E2& res, Float3 fac, Int3 ly,
 }
 
 template <typename E1, typename E2>
-inline void cuda_marder_correct_xyz(E1& efield, E2& res, Float3 fac, Int3 lx,
-                                    Int3 rx, Int3 ly, Int3 ry, Int3 lz, Int3 rz)
+inline void cuda_marder_correct_xyz(E1& efield, E2& res, Float3 fac, Int3 l,
+                                    Int3 r)
 {
   auto k_efield = efield.to_kernel();
   auto k_res = res.to_kernel();
   gt::launch<3>(
     {k_efield.shape(0), k_efield.shape(1), k_efield.shape(2)},
     GT_LAMBDA(int ix, int iy, int iz) {
-      if ((ix >= lx[0] && ix < rx[0]) && (iy >= lx[1] && iy < rx[1]) &&
-          (iz >= lx[2] && iz < rx[2])) {
+      if ((ix >= l[0] && ix < r[0]) && (iy >= l[1] && iy < r[1]) &&
+          (iz >= l[2] && iz < r[2])) {
         k_efield(ix, iy, iz, 0) =
           k_efield(ix, iy, iz, 0) +
           fac[0] * (k_res(ix, iy + 1, iz) - k_res(ix, iy, iz));
       }
 
-      if ((ix >= ly[0] && ix < ry[0]) && (iy >= ly[1] && iy < ry[1]) &&
-          (iz >= ly[2] && iz < ry[2])) {
+      if ((ix >= l[0] && ix < r[0]) && (iy >= l[1] && iy < r[1]) &&
+          (iz >= l[2] && iz < r[2])) {
         k_efield(ix, iy, iz, 1) =
           k_efield(ix, iy, iz, 1) +
           fac[1] * (k_res(ix, iy + 1, iz) - k_res(ix, iy, iz));
       }
 
-      if ((ix >= lz[0] && ix < rz[0]) && (iy >= lz[1] && iy < rz[1]) &&
-          (iz >= lz[2] && iz < rz[2])) {
+      if ((ix >= l[0] && ix < r[0]) && (iy >= l[1] && iy < r[1]) &&
+          (iz >= l[2] && iz < r[2])) {
         k_efield(ix, iy, iz, 2) =
           k_efield(ix, iy, iz, 2) +
           fac[2] * (k_res(ix, iy, iz + 1) - k_res(ix, iy, iz));
@@ -172,15 +132,15 @@ inline void correct(const Grid_t& grid, E1& efield, const Int3& efield_ib,
   assert(mf_ib == -grid.ibn);
   // OPT, do all patches in one kernel
   for (int p = 0; p < grid.n_patches(); p++) {
-    Int3 lx, rx, ly, ry, lz, rz;
-    detail::find_limits(grid, p, lx, rx, ly, ry, lz, rz);
+    Int3 l = grid.ibn;
+    Int3 r = grid.ibn + grid.ldims;
 
     auto p_efield = efield.view(_all, _all, _all, _all, p);
     auto p_res = mf.view(_all, _all, _all, 0, p);
     if (grid.isInvar(0)) {
-      cuda_marder_correct_yz(p_efield, p_res, fac, ly, ry, lz, rz);
+      cuda_marder_correct_yz(p_efield, p_res, fac, l, r);
     } else {
-      cuda_marder_correct_xyz(p_efield, p_res, fac, lx, rx, ly, ry, lz, rz);
+      cuda_marder_correct_xyz(p_efield, p_res, fac, l, r);
     }
   }
 }
@@ -201,7 +161,7 @@ public:
   using Bnd = BND;
   using real_t = typename storage_type::value_type;
 
-  // FIXME: checkpointing won't properly restore state
+  // FIXME: checkpointing won't properl restore state
 
   MarderCommon(const Grid_t& grid, real_t diffusion, int loop, bool dump)
     : diffusion_{diffusion}, loop_{loop}, dump_{dump}
@@ -262,6 +222,38 @@ public:
       Int3 res_ib = -grid.ibn;
       auto res = storage_type{psc::mflds::make_shape(grid, 1, res_ib)};
       psc::mflds::interior(grid, res) = dive - rho;
+
+      // Gauss' law is ostensibly violated at some boundaries, where virtual
+      // charges (i.e., charges that aren't associated with actual particles)
+      // implicitly shape the electric field. To account for virtual charges,
+      // simply set the error at those boundaries to 0.
+      for (int p = 0; p < grid.n_patches(); p++) {
+        for (int d = 0; d < 3; d++) {
+          if ((grid.bc.fld_lo[d] == BND_FLD_CONDUCTING_WALL ||
+               grid.bc.fld_lo[d] == BND_FLD_OPEN) &&
+              grid.atBoundaryLo(p, d)) {
+
+            gt::gslice slices[3] = {_s(grid.ibn[0], -grid.ibn[0]),
+                                    _s(grid.ibn[1], -grid.ibn[1]),
+                                    _s(grid.ibn[2], -grid.ibn[2])};
+            slices[d].stop = slices[d].start + 1;
+            res.view(slices[0], slices[1], slices[2], 0, p) = 0.0;
+          }
+
+          if ((grid.bc.fld_hi[d] == BND_FLD_CONDUCTING_WALL ||
+               grid.bc.fld_hi[d] == BND_FLD_OPEN) &&
+              grid.atBoundaryHi(p, d)) {
+            gt::gslice slices[3] = {_s(grid.ibn[0], -grid.ibn[0]),
+                                    _s(grid.ibn[1], -grid.ibn[1]),
+                                    _s(grid.ibn[2], -grid.ibn[2])};
+            // Note that upper edges are in the ghost region.
+            slices[d].start = slices[d].stop;
+            slices[d].stop += 1;
+            res.view(slices[0], slices[1], slices[2], 0, p) = 0.0;
+          }
+        }
+      }
+
       bnd_.fill_ghosts(grid, res, res_ib, 0, 1);
 
       print_progress(grid, rho, dive, res);

src/libpsc/psc_push_particles/push_particles_1vb.hxx

@@ -2,50 +2,6 @@
 // J. Villasenor and O. Buneman, "Rigorous charge conservation for local
 // electromagnetic field solvers", Computer Physics Communications 69 (1992) 306
 
-namespace
-{
-/**
- * @brief "Pushes" an exiting particle to the outer edge of the first layer of
- * ghost cells. This should happen before current deposition, so the full
- * exiting current is deposited. Note that higher order particles would need to
- * be pushed deeper into the ghost region.
- *
- * Actually, this function works with the grid-normalized particle position, not
- * the particle itself (or even its actual position). The grid-normalized
- * position is what's actually used to deposit current, and the particle is
- * dropped later, so its position doesn't need to be updated.
- * @tparam Real `float` or `double`
- * @param final_x_normed the particle's grid-normalized position, which is
- * possibly mutated
- * @param final_i3 the particle's 3d cell index
- * @param grid the grid
- * @param p the patch index
- */
-template <typename Real>
-void exit_to_edge(Vec3<Real>& final_x_normed, const Int3& final_i3,
-                  const Grid_t& grid, int p)
-{
-  // FIXME #948112531345 (also see the other FIXMEs with this id)
-  // Current deposited in ghost corners isn't sent to other patches, and
-  // thus leads to false-positive gauss errors.
-  // This is a general problem with non-periodic boundaries that have currents
-  // in ghost cells (i.e., just open BCs as of now), not just this function.
-  for (int d = 0; d < 3; d++) {
-    if (grid.bc.prt_lo[d] == BND_PRT_OPEN && grid.atBoundaryLo(p, d)) {
-      if (final_i3[d] < 0) {
-        final_x_normed[d] = -1.0;
-      }
-    }
-
-    if (grid.bc.prt_hi[d] == BND_PRT_OPEN && grid.atBoundaryHi(p, d)) {
-      if (final_i3[d] >= grid.ldims[d]) {
-        final_x_normed[d] = grid.ldims[d] + 1.0;
-      }
-    }
-  }
-}
-} // namespace
-
 // ======================================================================
 // PushParticlesVb
 
@@ -110,8 +66,6 @@ struct PushParticlesVb
         Real3 final_pos_normalized = prt.x() * dxi;
         Int3 final_index = final_pos_normalized.fint();
 
-        exit_to_edge(final_pos_normalized, final_index, grid, p);
-
         // CURRENT DENSITY BETWEEN (n+.5)*dt and (n+1.5)*dt
         Int3 initial_index;
         if (!Dim::InvarX::value) {

src/libpsc/tests/test_boundary_injector.cxx

@@ -89,12 +89,6 @@ struct ParticleGenerator
       uy = 0.0;
     }
 
-    // FIXME #948112531345 (also see the other FIXMEs with this id)
-    // If x[2] != 0, the injected particle deposits some
-    // current in the cells above its path. If one of those cells is in a ghost
-    // corner, the current deposited in that cell IS NOT sent across patches.
-    // That current is necessary to compute div E correctly. The run itself
-    // should be fine; it only throws off the gauss check.
     Real3 x = min_pos + pos_range * Real3{0, .999, 0.};
     Real3 u{0.0, uy, 0.0};
     Real w = 1.0;

src/libpsc/tests/test_push_fields.cxx

@@ -158,6 +158,11 @@ TYPED_TEST(PushFieldsTest, MarderCorrect)
 
   auto efield = mflds.storage().view(_all, _all, _all, _s(EX, EX + 3), _all);
   auto efield_ib = mflds.ib();
+
+  ASSERT_GT(gt::norm_linf(psc::mflds::interior(grid, mflds.gt()) -
+                          psc::mflds::interior(grid, mflds_ref.gt())),
+            1e-3);
+
   psc::marder::correct(grid, efield, efield_ib, mphi.storage(), mphi.ib(),
                        diffusion);