Fix regressions and improve unit tests

.gitlab-ci.yml

@@ -21,6 +21,7 @@ stages:
   image: ghcr.io/espressomd/docker/fedora:f7f8ef2c0ca93c67aa16b9f91785492fb04ecc1b
   variables:
     GIT_SUBMODULE_STRATEGY: none
+    GET_SOURCES_ATTEMPTS: 3
   before_script:
     - git config --global --add safe.directory ${CI_PROJECT_DIR}
   dependencies: []
@@ -32,6 +33,7 @@ stages:
 
 variables:
   GIT_SUBMODULE_STRATEGY: recursive
+  GET_SOURCES_ATTEMPTS: 3
   CCACHE_DIR: /cache
   CCACHE_MAXSIZE: 100G
   with_ccache: 'true'

src/core/electrostatics/elc.cpp

@@ -376,6 +376,7 @@ ElectrostaticLayerCorrection::z_energy(ParticleRange const &particles) const {
 
   if (elc.dielectric_contrast_on) {
     if (elc.const_pot) {
+      // metallic boundaries
       clear_vec(gblcblk, size);
       for (auto const &p : particles) {
         auto const z = p.pos()[2];
@@ -392,7 +393,7 @@ ElectrostaticLayerCorrection::z_energy(ParticleRange const &particles) const {
         }
       }
     } else {
-      // metallic boundaries
+      // dielectric boundaries
       auto const delta = elc.delta_mid_top * elc.delta_mid_bot;
       auto const fac_delta_mid_bot = elc.delta_mid_bot / (1. - delta);
       auto const fac_delta_mid_top = elc.delta_mid_top / (1. - delta);
@@ -447,6 +448,7 @@ void ElectrostaticLayerCorrection::add_z_force(
 
   if (elc.dielectric_contrast_on) {
     if (elc.const_pot) {
+      // metallic boundaries
       clear_vec(gblcblk, size);
       /* just counter the 2 pi |z| contribution stemming from P3M */
       for (auto const &p : particles) {
@@ -458,7 +460,7 @@ void ElectrostaticLayerCorrection::add_z_force(
           gblcblk[0] += elc.delta_mid_top * q;
       }
     } else {
-      // metallic boundaries
+      // dielectric boundaries
       auto const delta = elc.delta_mid_top * elc.delta_mid_bot;
       auto const fac_delta_mid_bot = elc.delta_mid_bot / (1. - delta);
       auto const fac_delta_mid_top = elc.delta_mid_top / (1. - delta);

src/walberla_bridge/include/walberla_bridge/lattice_boltzmann/LBWalberlaBase.hpp

@@ -59,6 +59,9 @@ class LBWalberlaBase : public LatticeModel {
   /** @brief Perform a ghost communication of the velocity field. */
   virtual void ghost_communication_vel() = 0;
 
+  /** @brief Perform a ghost communication of the last applied forces field. */
+  virtual void ghost_communication_laf() = 0;
+
   /** @brief Number of discretized velocities in the PDF. */
   virtual std::size_t stencil_size() const noexcept = 0;
 

src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp

@@ -640,7 +640,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
     m_pending_ghost_comm.set(GhostComm::VEL);
     m_pending_ghost_comm.set(GhostComm::LAF);
     // Refresh ghost layers
-    ghost_communication_pdfs();
+    ghost_communication_full();
   }
 
 protected:
@@ -668,11 +668,11 @@ class LBWalberlaImpl : public LBWalberlaBase {
   void ghost_communication() override {
     if (m_pending_ghost_comm.any()) {
       ghost_communication_boundary();
-      ghost_communication_pdfs();
+      ghost_communication_full();
     }
   }
 
-  void ghost_communication_pdf() override { // TODO: is this never traversed???
+  void ghost_communication_pdf() override {
     if (m_pending_ghost_comm.test(GhostComm::PDF)) {
       m_pdf_communicator->communicate();
       if (has_lees_edwards_bc()) {
@@ -694,7 +694,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
     }
   }
 
-  void ghost_communication_laf() {
+  void ghost_communication_laf() override {
     if (m_pending_ghost_comm.test(GhostComm::LAF)) {
       m_laf_communicator->communicate();
       if (has_lees_edwards_bc()) {
@@ -712,7 +712,7 @@ class LBWalberlaImpl : public LBWalberlaBase {
     }
   }
 
-  void ghost_communication_pdfs() {
+  void ghost_communication_full() {
     m_full_communicator->communicate();
     if (has_lees_edwards_bc()) {
       auto const &blocks = get_lattice().get_blocks();
@@ -975,15 +975,14 @@ class LBWalberlaImpl : public LBWalberlaBase {
     if (pos.empty()) {
       return {};
     }
+    std::vector<Utils::Vector3d> vel{};
     if constexpr (Architecture == lbmpy::Arch::CPU) {
-      std::vector<Utils::Vector3d> vel{};
       vel.reserve(pos.size());
       for (auto const &vec : pos) {
         auto res = get_velocity_at_pos(vec, true);
         assert(res.has_value());
         vel.emplace_back(*res);
       }
-      return vel;
     }
 #if defined(__CUDACC__)
     if constexpr (Architecture == lbmpy::Arch::GPU) {
@@ -1003,17 +1002,15 @@ class LBWalberlaImpl : public LBWalberlaBase {
       auto field =
           block.template uncheckedFastGetData<VectorField>(m_velocity_field_id);
       auto const res = lbm::accessor::Interpolation::get(field, host_pos, gl);
-      std::vector<Utils::Vector3d> vel{};
       vel.reserve(res.size() / 3ul);
       for (auto it = res.begin(); it != res.end(); it += 3) {
         vel.emplace_back(Utils::Vector3d{static_cast<double>(*(it + 0)),
                                          static_cast<double>(*(it + 1)),
                                          static_cast<double>(*(it + 2))});
       }
-      return vel;
     }
 #endif
-    return {};
+    return vel;
   }
 
   std::optional<Utils::Vector3d>

src/walberla_bridge/tests/LBWalberlaImpl_bspline_tests.cpp

@@ -144,8 +144,49 @@ BOOST_DATA_TEST_CASE(velocity_interpolation_bspline, bdata::make(all_lbs()),
   }
 }
 
-// TODO: check last applied force on a ghost node, i.e. when two forces
-// are applied at (agrid/2, 0, 0) and (box_l - agrid/2, 0, 0)
+BOOST_DATA_TEST_CASE(force_interpolation_ghosts, bdata::make(all_lbs()),
+                     lb_generator) {
+  auto lb = lb_generator(params);
+
+  auto const agrid = params.box_dimensions[0] / params.grid_dimensions[0];
+  auto const force = Vector3d{{1., 2., -3.}};
+  lb->add_force_at_pos(Vector3d{{-agrid / 2., 0., 0.}}, force);
+  lb->add_force_at_pos(Vector3d{{+agrid / 2., 0., 0.}}, force);
+  for (int x : {-1, 0, 1}) {
+    for (int y : {-1, 0, 1}) {
+      for (int z : {-1, 0, 1}) {
+        auto const check_node = Vector3i{{x, y, z}};
+        if (lb->get_lattice().node_in_local_halo(check_node)) {
+          auto const res = lb->get_node_force_to_be_applied(check_node);
+          if (x <= 0 and y <= 0 and z <= 0) {
+            BOOST_CHECK_SMALL(((*res) - force / 4.).norm(), 1E-10);
+          } else {
+            BOOST_CHECK_SMALL((*res).norm(), 1E-10);
+          }
+        }
+      }
+    }
+  }
+  lb->integrate();
+  lb->ghost_communication();
+  // last applied forces should now contain forces to be applied,
+  // and its ghost layer should have been zeroed out
+  for (int x : {-1, 0, 1}) {
+    for (int y : {-1, 0, 1}) {
+      for (int z : {-1, 0, 1}) {
+        auto const check_node = Vector3i{{x, y, z}};
+        if (lb->get_lattice().node_in_local_halo(check_node)) {
+          auto const res = lb->get_node_last_applied_force(check_node, true);
+          if (x == 0 and y == 0 and z == 0) {
+            BOOST_CHECK_SMALL(((*res) - force / 4.).norm(), 1E-10);
+          } else {
+            BOOST_CHECK_SMALL((*res).norm(), 1E-10);
+          }
+        }
+      }
+    }
+  }
+}
 
 int main(int argc, char **argv) {
   int n_nodes;

src/walberla_bridge/tests/LBWalberlaImpl_lees_edwards_tests.cpp

@@ -122,6 +122,16 @@ BOOST_AUTO_TEST_CASE(test_interpolation_force) {
   auto const ghost_node = Vector3i{force_node[0] - offset, -1, force_node[2]};
   auto const laf = *(lb->get_node_last_applied_force(ghost_node, true));
   BOOST_CHECK_SMALL((laf - f1).norm(), 1E-10);
+
+  // check setter
+  auto const f = Vector3d{{0.1, 0.2, -0.3}};
+  lb->set_node_last_applied_force(force_node, f);
+
+  lb->ghost_communication_laf();
+  lb->ghost_communication_vel();
+
+  auto const ghost_laf = *(lb->get_node_last_applied_force(ghost_node, true));
+  BOOST_CHECK_SMALL((ghost_laf - f).norm(), 1E-10);
 }
 
 BOOST_AUTO_TEST_CASE(test_interpolation_velocity) {
@@ -135,7 +145,8 @@ BOOST_AUTO_TEST_CASE(test_interpolation_velocity) {
   auto const v = Vector3d{0.3, -0.2, 0.3};
   lb->set_node_velocity(source_node, v);
 
-  lb->ghost_communication();
+  lb->ghost_communication_pdf();
+  lb->ghost_communication_vel();
 
   auto const ghost_node = Vector3i{source_node[0] - offset, -1, source_node[2]};
   auto const ghost_vel = *(lb->get_node_velocity(ghost_node, true));
@@ -158,7 +169,7 @@ BOOST_AUTO_TEST_CASE(test_interpolation_pdf) {
     x += .1;
   });
   lb->set_node_population(source_node, source_pop);
-  lb->ghost_communication();
+  lb->ghost_communication_pdf();
 
   auto const ghost_node = Vector3i{source_node[0] - offset, -1, source_node[2]};
   auto const ghost_pop = *(lb->get_node_population(ghost_node, true));

src/walberla_bridge/tests/LBWalberlaImpl_unit_tests.cpp

@@ -315,7 +315,8 @@ BOOST_DATA_TEST_CASE(velocity_at_node_and_pos, bdata::make(all_lbs()),
     }
   }
 
-  lb->ghost_communication();
+  lb->ghost_communication_pdf();
+  lb->ghost_communication_vel();
 
   // check velocities
   for (auto const &node : all_nodes_incl_ghosts(lb->get_lattice())) {
@@ -371,7 +372,7 @@ BOOST_DATA_TEST_CASE(density_at_pos, bdata::make(all_lbs()), lb_generator) {
     }
   }
 
-  lb->ghost_communication();
+  lb->ghost_communication_pdf();
 
   // check densities
   for (auto const &node : all_nodes_incl_ghosts(lb->get_lattice())) {
@@ -436,29 +437,45 @@ BOOST_DATA_TEST_CASE(forces_interpolation, bdata::make(all_lbs()),
                      lb_generator) {
   auto lb = lb_generator(params);
 
-  // todo: check a less symmetrical situation, where the force is applied not
-  // in the middle between the nodes
-
   for (auto const &n : all_nodes_incl_ghosts(lb->get_lattice())) {
     if (lb->get_lattice().node_in_local_halo(n)) {
-      auto const pos = 1. * n; // Mid point between nodes
-      auto const f = Vector3d{{1., 2., -3.5}};
-      lb->add_force_at_pos(pos, f);
-      lb->ghost_communication();
-      // Check neighboring nodes for force to be applied
-      for (int x : {0, 1})
-        for (int y : {0, 1})
-          for (int z : {0, 1}) {
-            auto const check_node = Vector3i{{n[0] - x, n[1] - y, n[2] - z}};
-            if (lb->get_lattice().node_in_local_halo(check_node)) {
-              auto const res = lb->get_node_force_to_be_applied(check_node);
-              BOOST_CHECK_SMALL(((*res) - f / 8.).norm(), 1E-10);
+      for (auto dir : {0u, 1u, 2u}) {
+        auto pos = 1. * n; // Mid point between nodes
+        pos[dir] += 0.25;
+        auto const f = Vector3d{{1., 2., -3.5}};
+        lb->add_force_at_pos(pos, f);
+        // Check neighboring nodes for force to be applied
+        for (int x : {0, 1}) {
+          for (int y : {0, 1}) {
+            for (int z : {0, 1}) {
+              auto const check_node = Vector3i{{n[0] - x, n[1] - y, n[2] - z}};
+              auto const weight = (check_node[dir] == n[dir]) ? 16. / 3. : 16.;
+              if (lb->get_lattice().node_in_local_halo(check_node)) {
+                auto const res = lb->get_node_force_to_be_applied(check_node);
+                BOOST_CHECK_SMALL(((*res) - f / weight).norm(), 1E-10);
+              }
             }
           }
-      // Apply counter force to clear force field
-      lb->add_force_at_pos(pos, -f);
-    } else {
-      lb->ghost_communication();
+        }
+        // Apply counter force to clear force field
+        lb->add_force_at_pos(pos, -f);
+      }
+    }
+  }
+}
+
+BOOST_DATA_TEST_CASE(last_applied_forces_setters, bdata::make(all_lbs()),
+                     lb_generator) {
+  auto lb = lb_generator(params);
+
+  for (auto const &n : all_nodes_incl_ghosts(lb->get_lattice(), false)) {
+    auto const f = Vector3d{{static_cast<double>(n[0] + 2), 2., -3.5}};
+    lb->set_node_last_applied_force(n, f);
+    lb->ghost_communication_laf();
+    lb->ghost_communication_vel();
+    if (lb->get_lattice().node_in_local_halo(n)) {
+      auto const res = lb->get_node_last_applied_force(n, true);
+      BOOST_CHECK_SMALL(((*res) - f).norm(), 1E-10);
     }
   }
 }