Refactor Lees-Edwards code

Capture Lees-Edwards parameters by reference, remove mention of Python
in the core, improve code coverage.

src/core/lb/particle_coupling.cpp

@@ -123,12 +123,12 @@ static void positions_in_halo_impl(Utils::Vector3d const &pos_folded,
             pos_folded + Utils::hadamard_product(box_geo.length(), shift);
 
         if (box_geo.type() == BoxType::LEES_EDWARDS) {
-          // note: Python-style modulo: 0<=offset <length
-          auto le = box_geo.lees_edwards_bc();
-          auto folded_offset =
-              std::fmod(le.pos_offset, box_geo.length()[le.shear_direction]);
+          // note: modulo convention: 0 <= offset < length
+          auto const &le = box_geo.lees_edwards_bc();
+          auto const length = box_geo.length()[le.shear_direction];
+          auto folded_offset = std::fmod(le.pos_offset, length);
           if (folded_offset < 0.) {
-            folded_offset += box_geo.length()[le.shear_direction];
+            folded_offset += length;
           }
           pos_shifted[le.shear_direction] +=
               shift[le.shear_plane_normal] * folded_offset;
@@ -166,10 +166,9 @@ std::vector<Utils::Vector3d> positions_in_halo(Utils::Vector3d const &pos,
   return res;
 }
 
-Utils::Vector3d
-lees_edwards_vel_shift(const Utils::Vector3d &pos_shifted_by_box_l,
-                       const Utils::Vector3d &orig_pos,
-                       const BoxGeometry &box_geo) {
+static auto lees_edwards_vel_shift(Utils::Vector3d const &pos_shifted_by_box_l,
+                                   Utils::Vector3d const &orig_pos,
+                                   BoxGeometry const &box_geo) {
   Utils::Vector3d vel_shift{{0., 0., 0.}};
   if (box_geo.type() == BoxType::LEES_EDWARDS) {
     auto le = box_geo.lees_edwards_bc();
@@ -291,7 +290,7 @@ void ParticleCoupling::kernel(std::vector<Particle *> const &particles) {
       auto const random_force = get_noise_term(p);
       force_on_particle = drag_force + random_force;
       ++it_interpolated_velocities;
-      ++it_positions_velocity_coupling++;
+      ++it_positions_velocity_coupling;
     }
 
     auto force_on_fluid = -force_on_particle;

src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp

@@ -31,16 +31,6 @@
 #include <stdexcept>
 #include <utility>
 
-namespace {
-template <typename T> T python_modulo(T a, T b) {
-  T res = std::fmod(a, b);
-  if (res < 0)
-    return res + b;
-  else
-    return res;
-}
-} // namespace
-
 namespace walberla {
 
 /**
@@ -118,26 +108,22 @@ class InterpolateAndShiftAtBoundary {
     // the target
     auto const prefactor =
         ((slab_dir == m_slab_max) ? FloatType{-1} : FloatType{1});
-    auto const offset = get_pos_offset() * prefactor;
-    auto const folded_offset = python_modulo(offset, length);
+    auto const offset = static_cast<FloatType>(get_pos_offset()) * prefactor;
+    auto const folded_offset = modulo(offset, length);
     // 0<=folded_offset<length
-    auto const weight1 = 1 - std::fmod(folded_offset, FloatType{1});
+    auto const weight1 = FloatType{1} - std::fmod(folded_offset, FloatType{1});
     auto const weight2 = std::fmod(folded_offset, FloatType{1});
     for (auto const &&cell : ci) {
       Cell source1 = cell;
       Cell source2 = cell;
-      int target_idx = cell[dir];
-      FloatType const source_pos = FloatType(target_idx) + folded_offset;
-      auto folded_source_pos = python_modulo(source_pos, length);
-      // 0<=folded_source_pos <length
+      auto const source_pos = static_cast<FloatType>(cell[dir]) + folded_offset;
+      auto const folded_source_pos = modulo(source_pos, length);
+      // 0 <= folded_source_pos < length
       source1[dir] = cell_idx_c(std::floor(folded_source_pos));
-      //  0<=source1[dir]<length, i.e. integer value sbetw. 0  and length-1
-      //  inclusive
-      source2[dir] =
-          cell_idx_c(python_modulo(FloatType(source1[dir] + 1), length));
-      // ineger values between 0 and length -1 inclusive
-
-      for (uint_t q = 0; q < FieldType::F_SIZE; ++q) {
+      // 0 <= source1[dir] < length, i.e. integer value up to length-1 inclusive
+      source2[dir] = cell_idx_c(modulo(FloatType(source1[dir] + 1), length));
+      // integer value between 0 and length -1 inclusive
+      for (uint_t q = 0u; q < FieldType::F_SIZE; ++q) {
         tmp_field->get(cell, q) =
             field->get(source1, q) * weight1 + field->get(source2, q) * weight2;
       }
@@ -152,6 +138,11 @@ class InterpolateAndShiftAtBoundary {
     }
   }
 
+  FloatType modulo(FloatType a, FloatType b) const {
+    auto const res = std::fmod(a, b);
+    return (res < FloatType{0}) ? res + b : res;
+  }
+
 private:
   std::shared_ptr<StructuredBlockForest> m_blocks;
   BlockDataID m_field_id;

testsuite/python/lb_lees_edwards_particle_coupling.py

@@ -54,28 +54,14 @@ def coupling_weight(pos_lb_units, node_idx, lb_shape):
     #    minimum image convention, node and coupling position can be
     #    at different sides of a periodic boundary
     dx = np.abs(min_image_dist(pos_lb_units, node_idx, lb_shape))
-    # If the coupling point is >=1 lattice constant away from the node, no coupling.
-    if np.any(dx >= 1): 
-        weight = 0
-    else:
-        # distance pos to node via module with lattice constant 1
-        weight = np.product(1 - dx)
+    # If the coupling point is >=1 lattice constant away from the node,
+    # no coupling. Otherwise, distance pos to node via module with lattice
+    # constant 1
+    weight = 0. if np.any(dx >= 1.) else np.product(1. - dx)
 
     return weight
 
 
-class MockLBF:
-    shape = np.array((10, 10, 10))
-    agrid = 1
-
-    def __getitem__(self, idx):
-        assert within_grid(idx, self.shape)
-        return f"node {idx}"
-
-
-mock_lbf = MockLBF()
-
-
 def le_aware_lb_nodes_around_pos(
         folded_pos, lbf, le_pos_offset, shear_direction, shear_plane_normal):
     """Returns LB node(s) relevant for interpolation around the given position"""
@@ -156,13 +142,16 @@ def lb_node(idx): return lbf[fold_index(idx, lbf.shape)]
 
 
 @utx.skipIfMissingFeatures("WALBERLA")
+@ut.skipIf(np.prod(system.cell_system.node_grid) != 1, "Requires 1 MPI rank")
 class LBLeesEdwardsParticleCoupling(ut.TestCase):
+    """Test LB Lees-Edwards corner cases with a random RNG seed (smoke test)"""
+
     def test_viscous_coupling_with_offset(self):
         system.lb = None
         system.time_step = 1
         system.cell_system.skin = 0.1
         system.cell_system.set_n_square()
-        offset = 2 * (np.random.random() - 1) * 3 * system.box_l[1]
+        offset = (np.random.random() - 1.) * 6. * system.box_l[1]
         protocol = lees_edwards.LinearShear(
             shear_velocity=0, initial_pos_offset=offset, time_0=0.)
         system.lees_edwards.set_boundary_conditions(
@@ -244,10 +233,10 @@ def v_x(x): return np.interp(
             if abs(y <= 0.5):
                 pref = -1.
                 dist_to_unshifted_lb_nodes = 0.5 - y
-            elif y >= system.box_l[1] - 0.5:
+            else:
+                assert y >= system.box_l[1] - 0.5
                 pref = 1.
                 dist_to_unshifted_lb_nodes = y - (system.box_l[2] - 0.5)
-            else: raise Exception()
             vel_shift = pref * shear_vel
             xs = 0.5 + np.arange(lbf.shape[0])
             ys = [v_x(x - pref * pos_offset) for x in xs]