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,22 @@ 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,7 @@ 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>;
 
@@ -172,6 +172,7 @@ 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);

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

@@ -70,8 +70,8 @@ class LBBoundary : public AutoParameters<LBBoundary> {
         const auto rho = lb_lbfluid_get_density();
         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;
+        const double unit_conversion = agrid/tau/tau;
+            //agrid * agrid * agrid * agrid / rho / tau / tau;
         return m_lbboundary->get_force() * unit_conversion;
       }
       return none;

testsuite/python/lb_shear.py

@@ -139,6 +139,18 @@ 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)
+        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,12 +162,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)

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