Fix errors in LB+Lees-Edwards

maintainer/benchmarks/lb.py

@@ -115,7 +115,7 @@
     system.non_bonded_inter[0, 0].lennard_jones.set_params(
         epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")
     system.part.add(pos=np.random.random((n_part, 3)) * system.box_l)
-    benchmarks.minimize(system, n_part / 2.)
+    benchmarks.minimize(system, n_part / 20.)
     system.integrator.set_vv()
     system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=42)
 

src/core/lb/particle_coupling.cpp

@@ -123,11 +123,20 @@ std::vector<Utils::Vector3d> positions_in_halo(Utils::Vector3d const &pos,
                                                double agrid) {
   auto const halo = 0.5 * agrid;
   auto const halo_vec = Utils::Vector3d::broadcast(halo);
-  auto const fully_inside_lower = local_box.my_left() + 2. * halo_vec;
-  auto const fully_inside_upper = local_box.my_right() - 2. * halo_vec;
-  if (in_box(pos, fully_inside_lower, fully_inside_upper)) {
-    return {pos};
+
+  if (box_geo.type() == BoxType::CUBOID) {
+    // If the particle is at least one agrid away from the node boundary
+    // any ghosts shifted by +- box_length cannot be in the lb volume
+    // accessible by this node (-agrid/2 to box_length +agrid/2)
+    auto const fully_inside_lower = local_box.my_left() + 2. * halo_vec;
+    auto const fully_inside_upper = local_box.my_right() - 2. * halo_vec;
+    if (in_box(pos, fully_inside_lower, fully_inside_upper)) {
+      return {pos};
+    }
   }
+
+  // For remaingin particles, positions shifted by +- box_length
+  // can potentially be in the locally accessible LB volume
   auto const halo_lower_corner = local_box.my_left() - halo_vec;
   auto const halo_upper_corner = local_box.my_right() + halo_vec;
 
@@ -138,25 +147,49 @@ std::vector<Utils::Vector3d> positions_in_halo(Utils::Vector3d const &pos,
         Utils::Vector3d shift{{double(i), double(j), double(k)}};
         Utils::Vector3d pos_shifted =
             pos + Utils::hadamard_product(box_geo.length(), shift);
-
+        // Apply additional shift from Lees-Edwards boundary conditions, when
+        // shifting across a periodic boundary
         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 const &le = box_geo.lees_edwards_bc();
+          auto const folded_offset =
+              std::fmod(le.pos_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)) {
-          res.push_back(pos_shifted);
+          res.emplace_back(pos_shifted);
         }
       }
     }
   }
   return res;
 }
 
+/* Determine Lees-Edwards velocity offset for positions shifted
+/  by +- box_length in one or more coordinates,
+/  i.e., those obtained form positoins_in_halo()
+*/
+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 {
@@ -184,7 +217,8 @@ void ParticleCoupling::kernel(Particle &p) {
 #endif
     for (auto const &pos : halo_pos) {
       if (in_local_halo(m_local_box, pos, agrid)) {
-        auto const vel_offset = lb_drift_velocity_offset(p);
+        auto const vel_offset = lb_drift_velocity_offset(p) +
+                                lees_edwards_vel_shift(pos, p.pos(), m_box_geo);
         auto const drag_force =
             lb_drag_force(m_lb, m_thermostat.gamma, p, pos, vel_offset);
         auto const random_force = get_noise_term(p);

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,25 @@ 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);
-
-      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;
+      int target_idx = cell[dir];
+      auto const source_pos = 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] = python_modulo(source1[dir]+1, length);
+      // ineger values between 0 and length -1 inclusive
+
+      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

@@ -457,14 +457,14 @@ 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);
     // Refresh ghost layers
+    integrate_reset_force(blocks);
     ghost_communication_pdfs();
   }
 
@@ -565,7 +565,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
         std::make_shared<InterpolateAndShiftAtBoundary<VectorField, FloatType>>(
             blocks, m_last_applied_force_field_id, m_vec_tmp_field_id,
             n_ghost_layers, shear_direction, shear_plane_normal,
-            m_lees_edwards_callbacks->get_pos_offset);
+            [this]() { return -1.0*m_lees_edwards_callbacks->get_pos_offset();});
   }
 
   void check_lebc(unsigned int shear_direction,

src/walberla_bridge/tests/LBWalberlaImpl_lees_edwards.cpp

@@ -92,15 +92,15 @@ BOOST_AUTO_TEST_CASE(test_transient_shear) {
   }
 }
 
-static auto setup_lb_with_offset(double offset) {
+static auto setup_lb_with_offset(double offset, double vel_shift = 0.) {
   using LBImplementation = walberla::LBWalberlaImpl<double, lbmpy::Arch::CPU>;
   auto density = 1.;
   auto viscosity = 1. / 7.;
   auto lattice =
       std::make_shared<LatticeWalberla>(Vector3i{10, 10, 10}, mpi_shape, 1);
   auto lb = std::make_shared<LBImplementation>(lattice, viscosity, density);
   auto le_pack = std::make_unique<LeesEdwardsPack>(
-      0u, 1u, [=]() { return offset; }, []() { return 0.0; });
+      0u, 1u, [=]() { return offset; }, [=]() { return vel_shift; });
   lb->set_collision_model(std::move(le_pack));
   lb->ghost_communication();
   return lb;
@@ -125,22 +125,148 @@ BOOST_AUTO_TEST_CASE(test_interpolation_force) {
   BOOST_CHECK_SMALL((laf - f1).norm(), 1E-10);
 }
 
-BOOST_AUTO_TEST_CASE(test_interpolation_velocity) {
-  auto const offset = 2;
-  auto lb = setup_lb_with_offset(offset);
+int fold(int i, int size) {
+  while (i < 0)
+    i += size;
+  while (i >= size)
+    i -= size;
+  return i;
+}
+
+std::vector<Utils::Vector3i>
+get_mirroring_xz_ghosts(const Utils::Vector3i &n, const Utils::Vector3i shape) {
+  std::vector<Utils::Vector3i> res;
+
+  auto in_ghost_layer = [](int i, int shape) {
+    return (i == -1 or i == shape);
+  };
+  auto in_box_or_ghost_layer = [](int i, int shape) {
+    return (i >= -1 and i <= shape);
+  };
+
+  for (int dx : {-shape[0], 0, shape[0]}) {
+    for (int dz : {-shape[2], 0, shape[2]}) {
+      if (dx == 0 and dz == 0) // no shift
+        continue;
+
+      auto shifted = Utils::Vector3i{n[0] + dx, n[1], n[2] + dz};
+      if ((in_ghost_layer(shifted[0], shape[0]) and
+           in_box_or_ghost_layer(shifted[2], shape[2])) or
+          (in_ghost_layer(shifted[2], shape[2]) and
+           in_box_or_ghost_layer(shifted[0], shape[0]))) {
+        res.push_back(shifted);
+      }
+    }
+  }
+  return res;
+}
+template <typename LB>
+void check_mirroring_ghost_vel(const LB &lb, Vector3i orig_ghost_node) {
+  auto const ref_vel = *(lb->get_node_velocity(orig_ghost_node, true));
   auto const shape = lb->get_lattice().get_grid_dimensions();
-  auto const xz = shape[0] / 2;
-  auto const y_max = shape[1] - 1;
+  //    std::cerr << "mirrors of "<<orig_ghost_node<<" v="<<ref_vel<<" lb
+  //    shape="<<shape<<std::endl;
+  for (auto const mirror_node :
+       get_mirroring_xz_ghosts(orig_ghost_node, shape)) {
+    auto const mirror_vel = *(lb->get_node_velocity(mirror_node, true));
+    //    std::cerr << "node "<<mirror_node<<" v="<<mirror_vel<<std::endl;
+    BOOST_CHECK_SMALL((mirror_vel - ref_vel).norm(), 1E-10);
+  }
+}
 
-  auto const source_node = Vector3i{xz, y_max, xz};
-  auto const v = Vector3d{0.3, -0.2, 0.3};
-  lb->set_node_velocity(source_node, v);
+BOOST_AUTO_TEST_CASE(test_interpolation_velocity_int_offset) {
+  auto const offset = 12;
+  auto vel_shift = 1.5;
+  Vector3d vel_shift_vec{vel_shift, 0., 0.};
+  auto lb = setup_lb_with_offset(offset, vel_shift);
+  auto const shape = lb->get_lattice().get_grid_dimensions();
+  for (int x = 0; x < shape[0]; x++) {
+    for (int z : {0, 3, shape[2] - 1}) {
+      auto const y_max = shape[1] - 1;
+      auto const v_upper = Vector3d{0.3 + x, -0.2 + y_max, 0.3 + z};
+      auto const v_lower = Vector3d{0.5 + x, -0.7, 0.9 + z};
+
+      auto const upper_source_node = Vector3i{x, y_max, z};
+      auto const lower_source_node = Vector3i{x, 0, z};
+      auto const ghost_of_upper_node =
+          Vector3i{fold(upper_source_node[0] - offset, shape[0]), -1,
+                   upper_source_node[2]};
+      auto const ghost_of_lower_node =
+          Vector3i{fold(lower_source_node[0] + offset, shape[0]), y_max + 1,
+                   lower_source_node[2]};
+
+      lb->set_node_velocity(upper_source_node, v_upper);
+      lb->set_node_velocity(lower_source_node, v_lower);
 
+      lb->ghost_communication();
+
+      auto const ghost_vel_for_upper =
+          *(lb->get_node_velocity(ghost_of_upper_node, true));
+      auto const ghost_vel_for_lower =
+          *(lb->get_node_velocity(ghost_of_lower_node, true));
+
+      BOOST_CHECK_SMALL((ghost_vel_for_upper + vel_shift_vec - v_upper).norm(),
+                        1E-10);
+      BOOST_CHECK_SMALL((ghost_vel_for_lower - vel_shift_vec - v_lower).norm(),
+                        1E-10);
+      check_mirroring_ghost_vel(lb, ghost_of_upper_node);
+      check_mirroring_ghost_vel(lb, ghost_of_lower_node);
+    }
+  }
+}
+
+
+Vector3d pos_from_node(Vector3i n) {
+  return {0.5+n[0],0.5+n[1],0.5+n[2]};
+}
+
+BOOST_AUTO_TEST_CASE(test_interpolation_velocity_non_int_offset) {
+  auto const offset = .5;
+  auto vel_shift = 0;
+  Vector3d vel_shift_vec{vel_shift, 0., 0.};
+  auto lb = setup_lb_with_offset(offset, vel_shift);
+  auto const shape = lb->get_lattice().get_grid_dimensions();
+  for (int x = 0; x < shape[0]; x++) {
+    for (int z : {0, 3, shape[2] - 1}) {
+      auto const y_max = shape[1] - 1;
+      auto const v_upper = Vector3d{x, -0.2 + y_max, 0.3 + z};
+      auto const v_lower = Vector3d{x, -0.7, 0.9 + z};
+
+      auto const upper_source_node = Vector3i{x, y_max, z};
+      auto const lower_source_node = Vector3i{x, 0, z};
+      lb->set_node_velocity(upper_source_node, v_upper);
+      lb->set_node_velocity(lower_source_node, v_lower);
+      auto upper_source_pos = pos_from_node(upper_source_node);
+      auto lower_source_pos = pos_from_node(lower_source_node);
+      BOOST_CHECK_SMALL((*(lb->get_velocity_at_pos(lower_source_pos)) - v_lower).norm(), 1E-8);
+      BOOST_CHECK_SMALL((*(lb->get_velocity_at_pos(upper_source_pos)) - v_upper).norm(), 1E-8);
+
+    }
+  }
+
   lb->ghost_communication();
+
+  for (int y: {0,10,9,-1}) {
+    std::cout << "y: "<<y<<"| ";
+    for (int x = -1; x <= shape[0]; x++) {
+       Vector3i node{x,y,0};
+       std::cout << (*(lb->get_node_velocity(node,true)))[0] << " ";
+    }
+    std::cout << std::endl;
+  }
+  return;
+  for (int x = 0; x < shape[0]; x++) {
+    int z=0;
+    Vector3i node = {x,-1,z};
+    auto vel = *(lb->get_node_velocity(node,true));
+    BOOST_CHECK_SMALL((*(lb->get_velocity_at_pos(pos_from_node(node),true)) - vel).norm(), 1E-8);
+    check_mirroring_ghost_vel(lb,node);
 
-  auto const ghost_node = Vector3i{source_node[0] - offset, -1, source_node[2]};
-  auto const ghost_vel = *(lb->get_node_velocity(ghost_node, true));
-  BOOST_CHECK_SMALL((ghost_vel - v).norm(), 1E-10);
+    node = {x,10,z};
+    vel = *(lb->get_node_velocity(node,true));
+    BOOST_CHECK_SMALL((*(lb->get_velocity_at_pos(pos_from_node(node), true)) - vel).norm(), 1E-8);
+    check_mirroring_ghost_vel(lb,node);
+  }
 }
 
 BOOST_AUTO_TEST_CASE(test_interpolation_pdf) {