Fix and test LB Lees Edwards particle coupling (espressomd#4948)

* When finding ghosts of a position for applying forces to LB, take Lees-Edwards offset into account
* Correct interpolated fluid velocities for Lees-Edwards velocity shift when using an interpolation position across an LE boundary
* Simplify and correct the shifting and interpolation of LB PDF ghost layers, when a Lees-Edwards offset is applied
* Move application of forces from MD to the correct place for LB with pull scheme
* More extensive testing of 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) {
-          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;
+          // 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 += length;
+          }
+          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,25 @@ std::vector<Utils::Vector3d> positions_in_halo(Utils::Vector3d const &pos,
   return res;
 }
 
+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();
+    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 +264,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 +277,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

@@ -94,8 +94,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);
@@ -110,25 +108,24 @@ 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 offset = static_cast<FloatType>(get_pos_offset()) * prefactor;
+    auto const folded_offset = modulo(offset, length);
+    // 0<=folded_offset<length
+    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;
-      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);
-
-      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);
-
-      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;
+      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 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;
       }
       tmp_field->get(cell, m_shear_direction) -= prefactor * shift;
     }
@@ -141,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;

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();
   }