WIP: Walberla boundary force

src/core/CMakeLists.txt

@@ -91,6 +91,7 @@ set(EspressoCore_SRC
   grid_based_algorithms/halo.cpp
   grid_based_algorithms/lattice.cpp
   grid_based_algorithms/lb_boundaries.cpp
+  grid_based_algorithms/lbboundaries/LBBoundary.cpp
   grid_based_algorithms/lb.cpp
   grid_based_algorithms/LbWalberla.cpp
   grid_based_algorithms/lb_walberla_instance.cpp

src/core/grid_based_algorithms/LbWalberla.cpp

@@ -55,8 +55,8 @@ LbWalberla::LbWalberla(double viscosity, double density, double agrid,
 
   Utils::Vector3i grid_dimensions;
   for (int i = 0; i < 3; i++) {
-    if (fabs(floor(box_dimensions[i] / agrid) * agrid - box_dimensions[i]) >
-        std::numeric_limits<double>::epsilon()) {
+    if (fabs((box_dimensions[i] / agrid) * agrid - box_dimensions[i]) >
+        2 * std::numeric_limits<double>::epsilon()) {
       throw std::runtime_error(
           "Box length not commensurate with agrid in direction " +
           std::to_string(i));
@@ -293,6 +293,24 @@ bool LbWalberla::add_force_at_pos(const Utils::Vector3d &pos,
   return true;
 }
 
+boost::optional<Utils::Vector3d>
+LbWalberla::get_node_boundary_force(const Utils::Vector3i node) const {
+  auto bc = get_block_and_cell(node, true); // including ghosts
+  if (!bc)
+    return {boost::none};
+  // Get boundary hanlindg
+  auto const &bh =
+      (*bc).block->getData<Boundary_handling_t>(m_boundary_handling_id);
+  auto const &ff = (*bc).block->getData<Flag_field_t>(m_flag_field_id);
+  if (!ff->isFlagSet((*bc).cell, ff->getFlag(UBB_flag)))
+    return {boost::none};
+
+  auto const uid = bh->getBoundaryUID(UBB_flag);
+  auto const &ubb = bh->getBoundaryCondition<UBB_t>(uid);
+  return {to_vector3d(
+      ubb.getForce((*bc).cell.x(), (*bc).cell.y(), (*bc).cell.z()))};
+}
+
 boost::optional<Utils::Vector3d>
 LbWalberla::get_node_velocity(const Utils::Vector3i node) const {
   auto bc = get_block_and_cell(node);

src/core/grid_based_algorithms/LbWalberla.hpp

@@ -112,7 +112,8 @@ class LbWalberla {
   using Flag_field_t = walberla::FlagField<walberla::uint8_t>;
   using Pdf_field_t = walberla::lbm::PdfField<Lattice_model_t>;
 
-  using UBB_t = walberla::lbm::UBB<Lattice_model_t, walberla::uint8_t>;
+  using UBB_t =
+      walberla::lbm::UBB<Lattice_model_t, walberla::uint8_t, true, true>;
   using Boundary_handling_t =
       walberla::BoundaryHandling<Flag_field_t, Lattice_model_t::Stencil, UBB_t>;
 
@@ -173,6 +174,8 @@ class LbWalberla {
   bool add_force_at_pos(const Utils::Vector3d &position,
                         const Utils::Vector3d &force);
   boost::optional<Utils::Vector3d>
+  get_node_boundary_force(const Utils::Vector3i node) const;
+  boost::optional<Utils::Vector3d>
   get_velocity_at_pos(const Utils::Vector3d &position) const;
   boost::optional<double> get_density_at_pos(const Utils::Vector3d &position);
   //  Utils::Vector3d get_stress_at_pos(const Utils::Vector3d& position);

src/core/grid_based_algorithms/lbboundaries/LBBoundary.cpp

@@ -0,0 +1,41 @@
+#include "LBBoundary.hpp"
+#include "config.hpp"
+#ifdef LB_BOUNDARIES
+#include "grid_based_algorithms/lb_interface.hpp"
+#include "grid_based_algorithms/lb_walberla_instance.hpp"
+namespace LBBoundaries {
+Utils::Vector3d LBBoundary::get_force() const {
+#if defined(LB_BOUNDARIES) || defined(LB_BOUNDARIES_GPU)
+  if (lattice_switch == ActiveLB::WALBERLA) {
+    auto const grid = lb_walberla()->get_grid_dimensions();
+    auto const agrid = lb_lbfluid_get_agrid();
+    Utils::Vector3d force{0, 0, 0};
+    for (auto index_and_pos : lb_walberla()->global_node_indices_positions()) {
+      // Convert to MD units
+      auto const index = index_and_pos.first;
+      auto const pos = index_and_pos.second * agrid;
+      if (calc_dist(pos) <= 0.)
+        for (int dx : {-1, 0, 1})
+          for (int dy : {-1, 0, 1})
+            for (int dz : {-1, 0, 1}) {
+              Utils::Vector3i shifted_index =
+                  index +
+                  Utils::Vector3i{dx * grid[0], dy * grid[1], dz * grid[2]};
+              auto node_force =
+                  lb_walberla()->get_node_boundary_force(shifted_index);
+              if (node_force) {
+                force += (*node_force);
+              }
+            } // Loop over cells
+    }         // loop over lb cells
+    boost::mpi::all_reduce(comm_cart, force, std::plus<Utils::Vector3d>());
+    return force;
+  } else {
+    return lbboundary_get_force(this);
+  }
+#else
+  throw std::runtime_error("Needs LB_BOUNDARIES or LB_BOUNDARIES_GPU.");
+#endif
+}
+} // namespace LBBoundaries
+#endif

src/core/grid_based_algorithms/lbboundaries/LBBoundary.hpp

@@ -49,6 +49,12 @@ class LBBoundary {
     m_shape->calculate_dist(pos, dist, vec);
   }
 
+  double calc_dist(const Utils::Vector3d &pos) const {
+    double tmp[3], dist;
+    calc_dist(pos, &dist, tmp);
+    return dist;
+  }
+
   void set_shape(std::shared_ptr<Shapes::Shape> const &shape) {
     m_shape = shape;
   }
@@ -64,13 +70,7 @@ class LBBoundary {
   Shapes::Shape const &shape() const { return *m_shape; }
   Utils::Vector3d &velocity() { return m_velocity; }
   Utils::Vector3d &force() { return m_force; }
-  Utils::Vector3d get_force() const {
-#if defined(LB_BOUNDARIES) || defined(LB_BOUNDARIES_GPU)
-    return lbboundary_get_force(this);
-#else
-    throw std::runtime_error("Needs LB_BOUNDARIES or LB_BOUNDARIES_GPU.");
-#endif
-  }
+  Utils::Vector3d get_force() const;
 
 #ifdef EK_BOUNDARIES // TODO: ugly. Better would be a class EKBoundaries,
                      // deriving from LBBoundaries, but that requires completely

src/script_interface/lbboundaries/LBBoundary.hpp

@@ -71,7 +71,8 @@ class LBBoundary : public AutoParameters<LBBoundary> {
         const auto agrid = lb_lbfluid_get_agrid();
         const auto tau = lb_lbfluid_get_tau();
         const double unit_conversion =
-            agrid * agrid * agrid * agrid / rho / tau / tau;
+            agrid / tau / tau * agrid * agrid * agrid;
+        // agrid * agrid * agrid * agrid / rho / tau / tau;
         return m_lbboundary->get_force() * unit_conversion;
       }
       return none;

testsuite/python/lb_shear.py

@@ -29,10 +29,10 @@
 """
 
 
-AGRID = 1 
+AGRID = .5
 VISC = 5.2
-DENS = 2.3
-TIME_STEP = 0.02
+DENS = 2.3 
+TIME_STEP = 0.03
 # Box size will be H +2 AGRID to make room for walls.
 # The number of grid cells should be devisible by four and 3 in all directions
 # for testing on multiple mpi nodes.
@@ -139,6 +139,24 @@ def check_profile(self, shear_plane_normal, shear_direction):
                     np.copy(v_measured),
                     np.copy(v_expected[j]) * shear_direction, atol=3E-3)
 
+        print(wall1.get_force())
+        print(wall2.get_force())
+        p_eq = DENS * AGRID**2 / TIME_STEP**2 / 3
+        print(
+            shear_direction *
+            SHEAR_VELOCITY /
+            H *
+            W**2 *
+            VISC +
+            p_eq *
+            W**2 *
+            shear_plane_normal)
+        np.testing.assert_allclose(
+            np.copy(wall1.get_force()),
+            -np.copy(wall2.get_force()),
+            atol=1E-4)
+        np.testing.assert_allclose(np.copy(wall1.get_force()), 
+                                   shear_direction * SHEAR_VELOCITY / H * W**2 * VISC, atol=2E-4)
         # Test steady state stress tensor on a node
         p_eq = DENS * AGRID**2 / TIME_STEP**2 / 3
         p_expected = np.diag((p_eq, p_eq, p_eq))
@@ -150,13 +168,6 @@ def check_profile(self, shear_plane_normal, shear_direction):
             np.testing.assert_allclose(node_stress,
                                        p_expected, atol=1E-5, rtol=5E-3)
 
-        np.testing.assert_allclose(
-            np.copy(wall1.get_force()),
-            -np.copy(wall2.get_force()),
-            atol=1E-4)
-        np.testing.assert_allclose(np.copy(wall1.get_force()), 
-                                   shear_direction * SHEAR_VELOCITY / H * W**2 * VISC, atol=2E-4)
-
     def test(self):
         x = np.array((1, 0, 0), dtype=float)
         y = np.array((0, 1, 0), dtype=float)

testsuite/python/lb_stokes_sphere.py

@@ -123,5 +123,12 @@ def setUp(self):
         self.lbf = espressomd.lb.LBFluid(**LB_PARAMS)
 
 
+@utx.skipIfMissingFeatures(['LB_WALBERLA', 'EXTERNAL_FORCES'])
+class LBWalberlaStokes(ut.TestCase, Stokes):
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluidWalberla(**LB_PARAMS)
+
+
 if __name__ == "__main__":
     ut.main()