Fix and test LB Lees-Edwards particle coupling

src/core/lb/particle_coupling.cpp

@@ -123,12 +123,15 @@ 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 normal_shift = (pos_shifted - pos_folded)[le.shear_plane_normal];
-          if (normal_shift > std::numeric_limits<double>::epsilon())
-            pos_shifted[le.shear_direction] += le.pos_offset;
-          if (normal_shift < -std::numeric_limits<double>::epsilon())
-            pos_shifted[le.shear_direction] -= le.pos_offset;
+          auto folded_offset =
+              std::fmod(le.pos_offset, box_geo.length()[le.shear_direction]);
+          if (folded_offset < 0.) {
+            folded_offset += box_geo.length()[le.shear_direction];
+          }
+          pos_shifted[le.shear_direction] +=
+              shift[le.shear_plane_normal] * folded_offset;
         }
 
         if (in_box(pos_shifted, halo_lower_corner, halo_upper_corner)) {
@@ -163,6 +166,26 @@ 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) {
+  Utils::Vector3d vel_shift{{0., 0., 0.}};
+  if (box_geo.type() == BoxType::LEES_EDWARDS) {
+    auto le = box_geo.lees_edwards_bc();
+    auto normal_shift =
+        (pos_shifted_by_box_l - orig_pos)[le.shear_plane_normal];
+    // normal_shift is +,- box_l or 0 up to floating point errors
+    if (normal_shift > std::numeric_limits<double>::epsilon()) {
+      vel_shift[le.shear_direction] -= le.shear_velocity;
+    }
+    if (normal_shift < -std::numeric_limits<double>::epsilon()) {
+      vel_shift[le.shear_direction] += le.shear_velocity;
+    }
+  }
+  return vel_shift;
+}
+
 namespace LB {
 
 Utils::Vector3d ParticleCoupling::get_noise_term(Particle const &p) const {
@@ -242,6 +265,7 @@ void ParticleCoupling::kernel(std::vector<Particle *> const &particles) {
   auto const &domain_upper_corner = m_local_box.my_right();
   auto it_interpolated_velocities = interpolated_velocities.begin();
   auto it_positions_force_coupling = positions_force_coupling.begin();
+  auto it_positions_velocity_coupling = positions_velocity_coupling.begin();
   auto it_positions_force_coupling_counter =
       positions_force_coupling_counter.begin();
   for (auto ptr : coupled_particles) {
@@ -254,11 +278,20 @@ void ParticleCoupling::kernel(std::vector<Particle *> const &particles) {
 #endif
     Utils::Vector3d force_on_particle = {};
     if (coupling_mode == particle_force) {
-      auto const &v_fluid = *it_interpolated_velocities;
+      auto &v_fluid = *it_interpolated_velocities;
+      if (m_box_geo.type() == BoxType::LEES_EDWARDS) {
+        // Account for the case where the interpolated velocity has been read
+        // from a ghost of the particle across the LE boundary (or vice verssa)
+        // Then the particle velocity is shifted by +,- the LE shear velocity
+        auto const vel_correction = lees_edwards_vel_shift(
+            *it_positions_velocity_coupling, p.pos(), m_box_geo);
+        v_fluid += vel_correction;
+      }
       auto const drag_force = lb_drag_force(p, m_thermostat.gamma, v_fluid);
       auto const random_force = get_noise_term(p);
       force_on_particle = drag_force + random_force;
       ++it_interpolated_velocities;
+      ++it_positions_velocity_coupling++;
     }
 
     auto force_on_fluid = -force_on_particle;

src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp

@@ -31,6 +31,16 @@
 #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 {
 
 /**
@@ -94,8 +104,6 @@ class InterpolateAndShiftAtBoundary {
     auto const dir = m_shear_direction;
     auto const dim = cell_idx_c(m_blocks->getNumberOfCells(*block, dir));
     auto const length = numeric_cast<FloatType>(dim);
-    auto const weight =
-        std::abs(std::fmod(get_pos_offset() + length, FloatType{1}));
 
     // setup slab
     auto field = block->template getData<FieldType>(m_field_id);
@@ -111,24 +119,27 @@ class InterpolateAndShiftAtBoundary {
     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);
+    // 0<=folded_offset<length
+    auto const weight1 = 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;
-      source1[dir] = cell_idx_c(std::floor(
-                         static_cast<FloatType>(source1[dir]) + offset)) %
-                     dim;
-      source1[dir] = cell_idx_c(static_cast<FloatType>(source1[dir]) + length);
-      source1[dir] = cell_idx_c(source1[dir] % dim);
-
+      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
+      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(std::ceil(static_cast<FloatType>(source2[dir]) + offset)) %
-          dim;
-      source2[dir] = cell_idx_c(static_cast<FloatType>(source2[dir]) + length);
-      source2[dir] = cell_idx_c(source2[dir] % dim);
+          cell_idx_c(python_modulo(FloatType(source1[dir] + 1), length));
+      // ineger values between 0 and length -1 inclusive
 
-      for (uint_t f = 0; f < FieldType::F_SIZE; ++f) {
-        tmp_field->get(cell, f) = field->get(source1, f) * (1 - weight) +
-                                  field->get(source2, f) * weight;
+      for (uint_t q = 0; q < FieldType::F_SIZE; ++q) {
+        tmp_field->get(cell, q) =
+            field->get(source1, q) * weight1 + field->get(source2, q) * weight2;
       }
       tmp_field->get(cell, m_shear_direction) -= prefactor * shift;
     }

src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp

@@ -535,13 +535,13 @@ class LBWalberlaImpl : public LBWalberlaBase {
 
   void integrate_pull_scheme() {
     auto const &blocks = get_lattice().get_blocks();
-    integrate_reset_force(blocks);
     // Handle boundaries
     integrate_boundaries(blocks);
     // LB stream
     integrate_stream(blocks);
     // LB collide
     integrate_collide(blocks);
+    integrate_reset_force(blocks);
     // Refresh ghost layers
     ghost_communication_pdfs();
   }