WIP: LB+Lees Edwards fixes

src/core/grid_based_algorithms/lb_particle_coupling.cpp

@@ -18,6 +18,7 @@
  */
 #include "LocalBox.hpp"
 #include "Particle.hpp"
+#include "algorithm/periodic_fold.hpp"
 #include "cells.hpp"
 #include "communication.hpp"
 #include "config/config.hpp"
@@ -162,27 +163,34 @@ bool in_local_halo(Utils::Vector3d const &pos) {
  * @brief Return a vector of positions shifted by +,- box length in each
  * coordinate
  */
-std::vector<Utils::Vector3d> positions_in_halo(Utils::Vector3d pos,
-                                               const BoxGeometry &box) {
-  std::vector<Utils::Vector3d> res;
+std::vector<std::pair<Utils::Vector3d, Utils::Vector3d>>
+positions_in_halo(Utils::Vector3d pos, const BoxGeometry &box) {
+  std::vector<std::pair<Utils::Vector3d, Utils::Vector3d>> res;
   for (int i : {-1, 0, 1}) {
     for (int j : {-1, 0, 1}) {
       for (int k : {-1, 0, 1}) {
         Utils::Vector3d shift{{double(i), double(j), double(k)}};
         Utils::Vector3d pos_shifted =
             pos + Utils::hadamard_product(box.length(), shift);
-
+        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 - pos)[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 = Algorithm::periodic_fold(
+              le.pos_offset, box_geo.length()[le.shear_direction]);
+          if (normal_shift > std::numeric_limits<double>::epsilon()) {
+            pos_shifted[le.shear_direction] += folded_offset;
+            vel_shift[le.shear_direction] -= le.shear_velocity;
+          }
+          if (normal_shift < -std::numeric_limits<double>::epsilon()) {
+            pos_shifted[le.shear_direction] -= folded_offset;
+            vel_shift[le.shear_direction] += le.shear_velocity;
+          }
         }
 
         if (in_local_halo(pos_shifted)) {
-          res.push_back(pos_shifted);
+          res.push_back({pos_shifted, vel_shift});
+          //          std::cout << "coupling at "<<pos_shifted<<std::endl;
         }
       }
     }
@@ -232,7 +240,9 @@ void add_swimmer_force(Particle const &p, double time_step) {
 
     // couple positions including shifts by one box length to add forces
     // to ghost layers
-    for (auto pos : positions_in_halo(source_position, box_geo)) {
+    for (auto shift : positions_in_halo(source_position, box_geo)) {
+      auto const pos = shift.first;
+
       add_md_force(pos, force, time_step);
     }
   }
@@ -260,10 +270,12 @@ void couple_particle(Particle &p, bool couple_virtual, double noise_amplitude,
 
   // Calculate coupling force
   Utils::Vector3d coupling_force = {};
-  for (auto pos : positions_in_halo(p.pos(), box_geo)) {
+  for (auto shift : positions_in_halo(p.pos(), box_geo)) {
+    auto const pos = shift.first;
+    auto const lees_edwards_vel = shift.second;
     if (in_local_halo(pos)) {
-      auto const drag_force =
-          lb_drag_force(p, pos, lb_particle_coupling_drift_vel_offset(p));
+      auto const drag_force = lb_drag_force(
+          p, pos, lb_particle_coupling_drift_vel_offset(p) + lees_edwards_vel);
       auto const random_force =
           noise_amplitude * lb_particle_coupling_noise(noise_amplitude > 0.0,
                                                        p.id(), rng_counter);
@@ -274,7 +286,8 @@ void couple_particle(Particle &p, bool couple_virtual, double noise_amplitude,
 
   // couple positions including shifts by one box length to add
   // forces to ghost layers
-  for (auto pos : positions_in_halo(p.pos(), box_geo)) {
+  for (auto shift : positions_in_halo(p.pos(), box_geo)) {
+    auto const pos = shift.first;
     if (in_local_domain(pos)) {
       /* Particle is in our LB volume, so this node
        * is responsible to adding its force */

src/core/grid_based_algorithms/lb_particle_coupling.hpp

@@ -96,8 +96,8 @@ Utils::Vector3d lb_particle_coupling_noise(bool enabled, int part_id,
                                            const OptionalCounter &rng_counter);
 
 // internal function exposed for unit testing
-std::vector<Utils::Vector3d> positions_in_halo(Utils::Vector3d pos,
-                                               const BoxGeometry &box);
+std::vector<std::pair<Utils::Vector3d, Utils::Vector3d>>
+positions_in_halo(Utils::Vector3d pos, const BoxGeometry &box);
 
 // internal function exposed for unit testing
 void couple_particle(Particle &p, bool couple_virtual, double noise_amplitude,

src/core/virtual_sites/VirtualSitesInertialessTracers.cpp

@@ -59,7 +59,8 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) {
       return;
     }
     if (should_be_coupled(p, coupled_ghost_particles)) {
-      for (auto pos : positions_in_halo(p.pos(), box_geo)) {
+      for (auto shift : positions_in_halo(p.pos(), box_geo)) {
+        auto const pos = shift.first;
         add_md_force(pos * to_lb_units, -p.force(), time_step);
       }
     }
@@ -71,7 +72,8 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) {
       return;
     }
     if (should_be_coupled(p, coupled_ghost_particles)) {
-      for (auto pos : positions_in_halo(p.pos(), box_geo)) {
+      for (auto shift : positions_in_halo(p.pos(), box_geo)) {
+        auto const pos = shift.first;
         add_md_force(pos * to_lb_units, -p.force(), time_step);
       }
     }

src/walberla_bridge/src/lattice_boltzmann/InterpolateAndShiftAtBoundary.hpp

@@ -31,6 +31,16 @@
 #include <stdexcept>
 #include <utility>
 
+namespace {
+double python_modulo(double a, double b) {
+  double 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,35 @@ 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);
+    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];
+      auto ghost_fold = [dim](FloatType x) {
+        //         std::cout << "ghost fold " << x<< " ";
+        while (x > dim) {
+          x -= dim;
+        };
+        while (x < -1) {
+          x += dim;
+        };
+        //         std::cout <<x<<std::endl;
+        return x;
+      };
+      source1[dir] =
+          cell_idx_c(std::floor(ghost_fold(target_idx + folded_offset)));
       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;
+          cell_idx_c(std::floor(ghost_fold(target_idx + folded_offset + 1)));
+      //      std::cout << offset<<"->"<<folded_offset<<" - " << target_idx<<",
+      //      " <<source1[dir]<<", " <<source2[dir]<<" | " <<weight1<<", "
+      //      <<weight2<<std::endl;
+
+      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

@@ -33,6 +33,7 @@
 #include <field/GhostLayerField.h>
 #include <field/vtk/FlagFieldCellFilter.h>
 #include <field/vtk/VTKWriter.h>
+#include <iostream>
 
 #include <field/AddToStorage.h>
 #include <field/FlagField.h>

testsuite/python/lb_lees_edwards_particle_coupling.py

@@ -23,18 +23,18 @@
 import numpy as np
 import unittest_decorators as utx
 
+system = espressomd.System(box_l=[10, 10, 10])
+
 
 @utx.skipIfMissingFeatures("WALBERLA")
 class LBLeesEdwardsParticleCoupling(ut.TestCase):
-    def test(self):
-        system = espressomd.System(box_l=[10, 10, 10])
-
+    def test_viscous_coupling_with_offset(self):
+        system.actors.clear()
         system.time_step = 1
         system.cell_system.skin = 0.1
         system.cell_system.set_n_square()
 
-        offset = 1
-        idx = int(offset)
+        offset = 0.9  # (np.random.random()-1/2) * 5*system.box_l[0]
         protocol = lees_edwards.LinearShear(
             shear_velocity=0., initial_pos_offset=offset, time_0=0.)
         system.lees_edwards.set_boundary_conditions(
@@ -45,49 +45,119 @@ def test(self):
         system.actors.add(lbf)
         system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1)
 
+        idx = int(offset % system.box_l[0])  # lb x index incl offset 
+        # the particle is placed in the center of 8 lb points, i.e., 
+        # 0.5 away from the planes with the grid points
         pos = [system.box_l[0] / 2., 0., system.box_l[0] / 2.]
         p = system.part.add(pos=pos)
-        v0 = np.array([1, 2, 3])
-        mid_x = lbf.shape[0] // 2
+        mid_x = lbf.shape[0] // 2  # lb index for the origin particle
         mid_z = lbf.shape[2] // 2
 
-        upper_y = lbf.shape[1] - 1
-        nodes = [lbf[mid_x - 1, 0, mid_z],
-                 lbf[mid_x, 0, mid_z - 1],
-                 lbf[mid_x - 1, 0, mid_z],
-                 lbf[mid_x, 0, mid_z],
-                 lbf[mid_x - 1 + idx, upper_y, mid_z],
-                 lbf[mid_x + idx, upper_y, mid_z - 1],
-                 lbf[mid_x - 1 + idx, upper_y, mid_z],
-                 lbf[mid_x + idx, upper_y, mid_z]]
+        upper_y = lbf.shape[1] - 1  # y is shear plane noremal
+        # LB Nodes surrounding the particle.
+        # On the particle's side 
+        unshifted = [lbf[mid_x - 1, 0, mid_z - 1], 
+                     lbf[mid_x, 0, mid_z - 1],
+                     lbf[mid_x - 1, 0, mid_z],
+                     lbf[mid_x, 0, mid_z]]
+        # across the Lees Edwads boundary conditoin
+        # if offset ^1<.5, it couples to the left neighbor
+        # otherwise to the right
+        if (offset % 1) >= .5: extra = 1
+        else: extra = 0
+        shifted_left = [lbf[(mid_x - 1 + idx + extra) % lbf.shape[0], upper_y, mid_z - 1],
+                        lbf[(mid_x - 1 + idx + extra) % lbf.shape[0], upper_y, mid_z]]
+        shifted_right = [lbf[(mid_x + idx + extra) % lbf.shape[0], upper_y, mid_z - 1],
+                         lbf[(mid_x + idx + extra) % lbf.shape[0], upper_y, mid_z]]
+        nodes = shifted_left + shifted_right + unshifted
+
+        v0 = np.array([1, 2, 3])
         for n in nodes:
             n.velocity = v0
 
         system.integrator.run(1)
+
+        # Gather forces applied to the LB by the particle coupling
         lb_forces = np.array([n.last_applied_force for n in nodes])
         lb_force = np.sum(lb_forces, axis=0)
+
+        # total force on lb = - force on particle?
         np.testing.assert_allclose(lb_force, -np.copy(p.f))
-        for f in lb_forces:
-            np.testing.assert_allclose(f, lb_forces[0])
+        # Particle couples to 8 nodes. On the sie of the particle
+        # each lb node should have 1/8 of the force
+        for n in unshifted:
+            np.testing.assert_allclose(
+                np.copy(n.last_applied_force), -np.copy(p.f) / 8)
+
+        # Across the lees edwards boundary, forces ahead and behind
+        # the particle in shear directiion are differently distributed
+        # For offset 0 (in the middle between two nodes)
+        # left and right weighs are equal (0.5)
+        weight_r = (offset + .5 - idx) % 1
+        weight_l = 1 - weight_r
+        np.testing.assert_allclose(weight_l + weight_r, 1)
+        for w in (weight_l, weight_r):
+            assert w >= 0 and w <= 1
+        # The total weight is the product of the weights
+        # in all 3 Cartesian directions. These are
+        # weight_l/r in the shear direction and 1/2 in the other two
+        for n in shifted_left:
+            np.testing.assert_allclose(
+                np.copy(n.last_applied_force), -np.copy(p.f) * 1 / 2 * 1 / 2 * weight_l)
+        for n in shifted_right:
+            np.testing.assert_allclose(
+                np.copy(n.last_applied_force), -np.copy(p.f) * 1 / 2 * 1 / 2 * weight_r)
 
+        # Check the LB velocity interpolation
+        # at a positoin, where the ghost layer of the
+        # lees edwards shear plane contributes
+        lbf[:, :, :].velocity = np.zeros(3)
         lbf[:, :, :].velocity = [0, 0, 0]
 
-        lower_nodes = nodes[:4]
-        upper_nodes = nodes[4:]
+        lower_nodes = unshifted
+        upper_nodes = shifted_left + shifted_right
         for n in lower_nodes:
             n.velocity = v0
         for n in upper_nodes:
             n.velocity = - v0
+        # When the offset modulo 1 is not zero
+        # in the ghost layer, neighboring cell swith zero velocity contribute
+        expected_v = 1 / 2 * v0 +\
+            1 / 4 * -v0 + 1 / 4 * (1 - abs(1 - (offset % 1))) * -v0
         p.update(dict(pos=pos, v=np.zeros(3)))
         np.testing.assert_allclose(
             np.copy(lbf.get_interpolated_velocity(pos=pos)),
-            np.zeros(3))
-        system.integrator.run(1)
-        np.testing.assert_allclose(np.copy(p.pos), pos)
-        np.testing.assert_allclose(np.copy(p.f), np.zeros(3))
-        for n in nodes:
-            np.testing.assert_allclose(
-                np.copy(n.last_applied_force), np.zeros(3))
+            expected_v)
+
+    def test_viscous_coupling_with_shear_vel(self):
+        # Places a co-moving particle close to the LE boundary
+        # in shear flow. chesk that it remains force free
+        # this is only the case, if the periodic imagesin the 
+        # halo regoin calculate the drag force including the LE
+        # shear velocity.
+        system.actors.clear()
+        system.time_step = 0.1
+        system.cell_system.skin = 0.1
+        system.cell_system.set_n_square()
+        v_shear = 0.5
+        protocol = lees_edwards.LinearShear(
+            shear_velocity=v_shear, initial_pos_offset=0, time_0=0.)
+        system.lees_edwards.set_boundary_conditions(
+            shear_direction="x", shear_plane_normal="y", protocol=protocol)
+
+        lbf = espressomd.lb.LBFluidWalberla(
+            agrid=1., density=1., kinematic_viscosity=1., tau=system.time_step)
+        system.actors.add(lbf)
+        system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=1)
+        system.integrator.run(200) 
+        p = system.part.add(
+            pos=(
+                0, 0, 0), v=lbf.get_interpolated_velocity(
+                pos=(
+                    0, 0, 0)))
+        for _ in range(100): 
+            system.integrator.run(1)
+            np.testing.assert_allclose(p.f, np.zeros(3), atol=1E-7)
 
 
 if __name__ == '__main__':