Reimplement open boundary conditions

scripts/run_test.sh

@@ -1,6 +1,6 @@
 #!/bin/bash
 
-# Usage: scripts/run_test.sh name_of_test
+# Usage: scripts/run_test.sh name_of_test [gtest args]
 # Example: scripts/run_test.sh test_injector_boundary_inflow
 # Assumes you are in the root psc directory and have a build directory.
 
@@ -10,4 +10,4 @@ make $1
 mkdir -p runs
 cp ../bits/adios2cfg.xml runs/
 cd runs
-../src/libpsc/tests/$1 "${@:2}"
+mpirun -np 1 ../src/libpsc/tests/$1 "${@:2}"

src/include/bnd_particles_impl.hxx

@@ -7,7 +7,7 @@
 
 #include "balance.hxx"
 #include "ddc_particles.hxx"
-#include "bnd_particles.hxx"
+#include "particles.hxx"
 #include "particles_simple.hxx"
 
 extern int pr_time_step_no_comm;
@@ -17,10 +17,11 @@ extern double* psc_balance_comp_time_by_patch;
 // BndParticlesCommon
 
 template <typename MP>
-struct BndParticlesCommon : BndParticlesBase
+struct BndParticlesCommon
 {
   using Mparticles = MP;
   using real_t = typename Mparticles::real_t;
+  using Real3 = Vec3<real_t>;
   using ddcp_t = ddc_particles<Mparticles>;
   using BndBuffers = typename Mparticles::BndBuffers;
 
@@ -96,12 +97,9 @@ void BndParticlesCommon<MP>::process_patch(const Grid_t& grid,
   // New-style boundary requirements.
   // These will need revisiting when it comes to non-periodic domains.
 
+  const Int3& ldims = grid.ldims;
   const auto& gpatch = grid.patches[p];
-  const Int3& ldims = pi.ldims();
-  real_t xm[3];
-  for (int d = 0; d < 3; d++) {
-    xm[d] = gpatch.xe[d] - gpatch.xb[d];
-  }
+  Real3 patch_size = Real3(gpatch.xe - gpatch.xb);
 
   auto* dpatch = &ddcp->patches_[p];
   for (int dir1 = 0; dir1 < N_DIR; dir1++) {
@@ -121,12 +119,9 @@ void BndParticlesCommon<MP>::process_patch(const Grid_t& grid,
 
     Int3 pos = pi.cellPosition(xi);
 
-    if (pos[0] >= 0 &&
-        pos[0] < ldims[0] && // OPT, could be optimized with casts to unsigned
-        pos[1] >= 0 && pos[1] < ldims[1] && pos[2] >= 0 && pos[2] < ldims[2]) {
+    if (pi.isValidCellPosition(pos)) {
       // fast path
       // particle is still inside patch: move into right position
-      pi.validCellIndex(xi);
       buf[head++] = *prt;
       continue;
     }
@@ -139,7 +134,7 @@ void BndParticlesCommon<MP>::process_patch(const Grid_t& grid,
     for (int d = 0; d < 3; d++) {
       if (pos[d] < 0) {
         if (!grid.atBoundaryLo(p, d) || grid.bc.prt_lo[d] == BND_PRT_PERIODIC) {
-          xi[d] += xm[d];
+          xi[d] += patch_size[d];
           dir[d] = -1;
           int ci = pi.cellPosition(xi[d], d);
           if (ci >= ldims[d]) {
@@ -148,18 +143,26 @@ void BndParticlesCommon<MP>::process_patch(const Grid_t& grid,
           }
         } else {
           switch (grid.bc.prt_lo[d]) {
-            case BND_PRT_REFLECTING:
+            case BND_PRT_REFLECTING: {
               xi[d] = -xi[d];
               pxi[d] = -pxi[d];
               dir[d] = 0;
               break;
-            case BND_PRT_ABSORBING: drop = true; break;
+            }
+            case BND_PRT_OPEN: {
+              drop = true;
+              break;
+            }
+            case BND_PRT_ABSORBING: {
+              drop = true;
+              break;
+            }
             default: assert(0);
           }
         }
       } else if (pos[d] >= ldims[d]) {
         if (!grid.atBoundaryHi(p, d) || grid.bc.prt_hi[d] == BND_PRT_PERIODIC) {
-          xi[d] -= xm[d];
+          xi[d] -= patch_size[d];
           dir[d] = +1;
           int ci = pi.cellPosition(xi[d], d);
           if (ci < 0) {
@@ -168,7 +171,7 @@ void BndParticlesCommon<MP>::process_patch(const Grid_t& grid,
         } else {
           switch (grid.bc.prt_hi[d]) {
             case BND_PRT_REFLECTING: {
-              xi[d] = 2.f * xm[d] - xi[d];
+              xi[d] = 2.f * patch_size[d] - xi[d];
               pxi[d] = -pxi[d];
               dir[d] = 0;
               int ci = pi.cellPosition(xi[d], d);
@@ -177,7 +180,14 @@ void BndParticlesCommon<MP>::process_patch(const Grid_t& grid,
               }
               break;
             }
-            case BND_PRT_ABSORBING: drop = true; break;
+            case BND_PRT_OPEN: {
+              drop = true;
+              break;
+            }
+            case BND_PRT_ABSORBING: {
+              drop = true;
+              break;
+            }
             default: assert(0);
           }
         }
@@ -191,7 +201,7 @@ void BndParticlesCommon<MP>::process_patch(const Grid_t& grid,
           xi[d] = 0.f;
         }
         assert(xi[d] >= 0.f);
-        assert(xi[d] <= xm[d]);
+        assert(xi[d] <= patch_size[d]);
       }
     }
     if (!drop) {

src/include/boundary_injector.hxx

@@ -9,6 +9,7 @@
 #include "pushp.hxx"
 #include "dim.hxx"
 #include "setup_particles.hxx"
+#include "kg/VecRange.hxx"
 #include "../libpsc/psc_push_particles/inc_push.cxx"
 
 /// @brief A particle generator for use with @ref BoundaryInjector. Samples
@@ -72,15 +73,13 @@ public:
   using MfieldsState = typename PushParticles::MfieldsState;
   using AdvanceParticle_t = typename PushParticles::AdvanceParticle_t;
   using Current = typename PushParticles::Current;
-  using Dim = typename PushParticles::Dim;
   using real_t = typename PushParticles::real_t;
-  using Real3 = typename PushParticles::Real3;
-  using checks_order = typename PushParticles::checks_order;
+  using Real3 = Vec3<real_t>;
 
   BoundaryInjector(ParticleGenerator particle_generator, Grid_t& grid)
-    : partice_generator{particle_generator},
-      advance{grid.dt},
-      prts_per_unit_density{grid.norm.prts_per_unit_density}
+    : particle_generator_{particle_generator},
+      advance_{grid.dt},
+      prts_per_unit_density_{grid.norm.prts_per_unit_density}
   {}
 
   /// Injects particles at the lower y-bound as if there were a population of
@@ -97,84 +96,71 @@ public:
     const Grid_t& grid = mprts.grid();
     auto injectors_by_patch = mprts.injector();
 
-    Real3 dxi = Real3{1., 1., 1.} / Real3(grid.domain.dx);
+    Real3 dxi = grid.domain.dx_inv;
     Current current(grid);
 
     for (int patch_idx = 0; patch_idx < grid.n_patches(); patch_idx++) {
-      Int3 ilo = grid.patches[patch_idx].off;
-
-      if (ilo[INJECT_DIM_IDX_] != 0) {
+      if (!grid.atBoundaryLo(patch_idx, INJECT_DIM_IDX_)) {
         continue;
       }
 
-      Int3 ihi = ilo + grid.ldims;
+      Int3 ilo = {0, 0, 0};
+      Int3 ihi = grid.ldims;
+
+      ilo[INJECT_DIM_IDX_] = -1;
+      ihi[INJECT_DIM_IDX_] = 0;
 
       auto&& injector = injectors_by_patch[patch_idx];
       auto flds = mflds[patch_idx];
       typename Current::fields_t J(flds);
 
-      for (Int3 cell_idx = ilo; cell_idx[0] < ihi[0]; cell_idx[0]++) {
-        for (cell_idx[2] = ilo[2]; cell_idx[2] < ihi[2]; cell_idx[2]++) {
-          auto boundary_current_before =
-            flds.storage()(cell_idx[0] + grid.ibn[0], -1 + grid.ibn[1],
-                           cell_idx[2] + grid.ibn[2], JXI + 1);
-
-          cell_idx[INJECT_DIM_IDX_] = -1;
-
-          Real3 cell_corner =
-            grid.domain.corner + Double3(cell_idx) * grid.domain.dx;
-          int n_prts_to_try_inject =
-            get_n_in_cell(1.0, prts_per_unit_density, true);
-          real_t charge_injected = 0;
-
-          for (int prt_count = 0; prt_count < n_prts_to_try_inject;
-               prt_count++) {
-            psc::particle::Inject prt =
-              partice_generator.get(cell_corner, grid.domain.dx);
-
-            Real3 v = advance.calc_v(prt.u);
-            Real3 initial_x = prt.x;
-            advance.push_x(prt.x, v, 1.0);
-
-            if (prt.x[INJECT_DIM_IDX_] < grid.domain.corner[INJECT_DIM_IDX_]) {
-              // don't inject a particle that fails to enter the domain
-              continue;
-            }
-
-            injector(prt);
-
-            // Update currents
-            // Taken from push_particles_1vb.hxx PushParticlesVb::push_mprts()
+      for (Int3 initial_idx : VecRange(ilo, ihi)) {
+        Real3 cell_corner = Double3(initial_idx) * grid.domain.dx;
+        int n_prts_to_try_inject =
+          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);
+
+          Real3 v = advance_.calc_v(prt.u);
+          Real3 initial_x = prt.x;
+          advance_.push_x(prt.x, v);
+
+          if (prt.x[INJECT_DIM_IDX_] < 0.0) {
+            // don't inject a particle that fails to enter the patch
+            continue;
+          }
 
-            Real3 initial_normalized_pos = initial_x * dxi;
+          injector(prt);
 
-            Int3 final_idx;
-            Real3 final_normalized_pos;
-            find_idx_pos_1st_rel(prt.x, dxi, final_idx, final_normalized_pos,
-                                 real_t(0.));
+          // Update currents
+          // Taken from push_particles_1vb.hxx PushParticlesVb::push_mprts()
 
-            // CURRENT DENSITY BETWEEN (n+.5)*dt and (n+1.5)*dt
-            real_t qni_wni = grid.kinds[prt.kind].q * prt.w;
-            current.calc_j(J, initial_normalized_pos, final_normalized_pos,
-                           final_idx, cell_idx, qni_wni, v);
+          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;
 
-            charge_injected += qni_wni;
-          }
+          Real3 final_normalized_pos = prt.x * dxi;
+          Int3 final_idx = final_normalized_pos.fint();
 
-          // override whatever current was deposited in the boundary to be
-          // consistent with the charge having started completely out of the
-          // domain
-          flds.storage()(cell_idx[0] + grid.ibn[0], -1 + grid.ibn[1],
-                         cell_idx[2] + grid.ibn[2], JXI + 1) =
-            boundary_current_before + charge_injected * grid.domain.dx[1] /
-                                        grid.dt / prts_per_unit_density;
+          // CURRENT DENSITY BETWEEN (n+.5)*dt and (n+1.5)*dt
+          real_t qni_wni = grid.kinds[prt.kind].q * prt.w;
+          current.calc_j(J, initial_normalized_pos, final_normalized_pos,
+                         final_idx, initial_idx, qni_wni, v);
         }
       }
     }
   }
 
 private:
-  ParticleGenerator partice_generator;
-  AdvanceParticle<real_t, dim_y> advance;
-  real_t prts_per_unit_density;
+  ParticleGenerator particle_generator_;
+  AdvanceParticle_t advance_;
+  real_t prts_per_unit_density_;
 };

src/include/grid/domain.hxx

@@ -43,6 +43,7 @@ struct Domain
 
     ldims = gdims / np;
     dx = length / Real3(gdims);
+    dx_inv = Real3(gdims) / length;
   }
 
   void view() const
@@ -65,6 +66,8 @@ struct Domain
   Int3 ldims;
   /// @brief Side lengths of each grid cell, in physical units.
   Real3 dx;
+  /// @brief Inverses of side lengths of each grid cell, in physical units.
+  Real3 dx_inv;
 };
 
 template <typename R>

src/include/particle_indexer.hxx

@@ -66,7 +66,7 @@ struct ParticleIndexer
   using Real3 = Vec3<real_t>;
 
   ParticleIndexer(const Grid_t& grid)
-    : dxi_(Real3{1., 1., 1.} / Real3(grid.domain.dx)), ldims_(grid.ldims)
+    : dxi_(grid.domain.dx_inv), ldims_(grid.ldims)
   {
     n_cells_ = ldims_[0] * ldims_[1] * ldims_[2];
   }
@@ -98,10 +98,25 @@ struct ParticleIndexer
     return cellIndex(cpos);
   }
 
+  bool isValidCellPosition(const Int3& cpos) const
+  {
+    for (int d = 0; d < 3; d++) {
+      // optimization trick: if cpos[d]<0, it becomes larger than any positive
+      // int after casting both to uint
+      // TODO (cursed?) store ldims as uints, so this happens implicitly
+      if (uint(cpos[d]) >= ldims_[d]) {
+        return false;
+      }
+    }
+    return true;
+  }
+
   int validCellIndex(const real_t* pos) const
   {
     Int3 cpos = cellPosition(pos);
     for (int d = 0; d < 3; d++) {
+      // optimization trick: if cpos[d]<0, it becomes larger than any positive
+      // int after casting both to uint
       if (uint(cpos[d]) >= ldims_[d]) {
         printf("validCellIndex: cpos[%d] = %d ldims_[%d] = %d // pos[%d] = %g "
                "pos[%d]*dxi_[%d] = %g\n",

src/libpsc/integrate.cxx

@@ -5,7 +5,6 @@
 #include "push_particles.hxx"
 #include "sort.hxx"
 #include "collision.hxx"
-#include "bnd_particles.hxx"
 #include "checks_params.hxx"
 
 #include <mrc_common.h>

src/libpsc/psc_balance/psc_balance_impl.hxx

@@ -4,7 +4,6 @@
 #include "fields.hxx"
 #include "fields_traits.hxx"
 #include "balance.hxx"
-#include "bnd_particles.hxx"
 #include "bnd.hxx"
 #include "mpi_dtype_traits.hxx"