WIP: Boundary support for Walberla

src/core/grid_based_algorithms/LbWalberla.cpp

@@ -4,6 +4,7 @@
 #include "communication.hpp"
 #include "lb_walberla_instance.hpp"
 #include "utils/Vector.hpp"
+#include "utils/math/make_lin_space.hpp"
 #include "utils/mpi/gatherv.hpp"
 
 #include "blockforest/Initialization.h"
@@ -69,7 +70,9 @@ LbWalberla::LbWalberla(double viscosity, double density, double agrid,
     }
   }
   m_grid_dimensions = grid_dimensions;
-  printf("grid: %d %d %d, node: %d %d %d\n",grid_dimensions[0],grid_dimensions[1],grid_dimensions[2],node_grid[0],node_grid[1],node_grid[2]);
+  printf("grid: %d %d %d, node: %d %d %d\n", grid_dimensions[0],
+         grid_dimensions[1], grid_dimensions[2], node_grid[0], node_grid[1],
+         node_grid[2]);
 
   m_blocks = blockforest::createUniformBlockGrid(
       uint_c(node_grid[0]), // blocks in x direction
@@ -108,7 +111,7 @@ LbWalberla::LbWalberla(double viscosity, double density, double agrid,
       LB_boundary_handling(m_flag_field_id, m_pdf_field_id),
       "boundary handling");
 
-  empty_flags_for_boundary(m_blocks, m_boundary_handling_id);
+  clear_boundaries();
 
   // 1 timestep each time the integrate function is called
   m_time_loop =
@@ -125,20 +128,19 @@ LbWalberla::LbWalberla(double viscosity, double density, double agrid,
       ResetForce<Pdf_field_t, vector_field_t, Boundary_handling_t>>(
       m_pdf_field_id, m_force_field_id, m_force_field_from_md_id,
       m_boundary_handling_id);
-  m_time_loop->add()
-      << // timeloop::BeforeFunction(communication, "communication")
-      timeloop::Sweep(
-          Boundary_handling_t::getBlockSweep(m_boundary_handling_id),
-          "boundary handling");
+  m_time_loop->add() << timeloop::BeforeFunction(communication, "communication")
+                     << timeloop::Sweep(Boundary_handling_t::getBlockSweep(
+                                            m_boundary_handling_id),
+                                        "boundary handling");
   m_time_loop->add() << timeloop::Sweep(makeSharedSweep(m_reset_force),
                                         "Reset force fields");
-  m_time_loop->add()
-      << timeloop::AfterFunction(communication, "communication")
-      << timeloop::Sweep(
-             domain_decomposition::makeSharedSweep(
-                 lbm::makeCellwiseSweep<Lattice_model_t, Flag_field_t>(
-                     m_pdf_field_id, m_flag_field_id, Fluid_flag)),
-             "LB stream & collide");
+  //  m_time_loop->add()
+  //      << timeloop::AfterFunction(communication, "communication")
+  m_time_loop->add() << timeloop::Sweep(
+      domain_decomposition::makeSharedSweep(
+          lbm::makeCellwiseSweep<Lattice_model_t, Flag_field_t>(
+              m_pdf_field_id, m_flag_field_id, Fluid_flag)),
+      "LB stream & collide");
 
   m_force_distributor_id =
       field::addDistributor<Vector_field_distributor_t, Flag_field_t>(
@@ -181,7 +183,7 @@ LbWalberla::get_block_and_cell(const Utils::Vector3i &node,
     // Try to find a block which has the cell as ghost layer
     for (auto b = m_blocks->begin(); b != m_blocks->end(); ++b) {
       if (b->getAABB()
-              .getExtended(m_skin / m_agrid)
+              .getExtended(1)
               .contains(real_c(node[0]), real_c(node[1]), real_c(node[2]))) {
         block = &(*b);
         break;
@@ -222,7 +224,7 @@ IBlock *LbWalberla::get_block(const Utils::Vector3d &pos,
 
 bool LbWalberla::set_node_velocity_at_boundary(const Utils::Vector3i node,
                                                const Utils::Vector3d v) {
-  auto bc = get_block_and_cell(node);
+  auto bc = get_block_and_cell(node,true);
   // Return if we don't have the cell.
   if (!bc)
     return false;
@@ -423,4 +425,56 @@ bool LbWalberla::pos_in_local_domain(const Utils::Vector3d &pos) const {
   return (block != nullptr);
 }
 
+std::vector<std::pair<Utils::Vector3i, Utils::Vector3d>>
+LbWalberla::node_indices_positions() {
+
+  std::vector<std::pair<Utils::Vector3i, Utils::Vector3d>> res;
+  for (auto block = m_blocks->begin(); block != m_blocks->end(); ++block) {
+    auto left = block->getAABB().min();
+    // Lattice constant is 1, node centers are offset by .5
+    Utils::Vector3d pos_offset =
+        to_vector3d(left) + Utils::Vector3d::broadcast(.5);
+
+    // Lattice constant is 1, so cast left corner position to ints
+    Utils::Vector3i index_offset =
+        Utils::Vector3i{int(left[0]), int(left[1]), int(left[2])};
+
+    // Get field data which knows about the indices
+    // In the loop, x,y,z are in block-local coordinates
+    auto pdf_field = block->getData<Pdf_field_t>(m_pdf_field_id);
+    for (int x : Utils::make_lin_space(int(0), int(pdf_field->xSize()),
+                                       int(pdf_field->xSize()), false)) {
+      for (int y : Utils::make_lin_space(int(0), int(pdf_field->ySize()),
+                                         int(pdf_field->ySize()), false)) {
+        for (int z : Utils::make_lin_space(int(0), int(pdf_field->zSize()),
+                                           int(pdf_field->zSize()), false)) {
+          res.push_back(
+              {index_offset + Utils::Vector3i{x, y, z},
+               pos_offset + Utils::Vector3d{double(x), double(y), double(z)}});
+        }
+      }
+    }
+  }
+  return res;
+}
+
+std::vector<std::pair<Utils::Vector3i, Utils::Vector3d>>
+LbWalberla::global_node_indices_positions() {
+
+  std::vector<std::pair<Utils::Vector3i, Utils::Vector3d>> res;
+  for (int x : Utils::make_lin_space(int(0), int(m_grid_dimensions[0]),
+                                     int(m_grid_dimensions[0]), false)) {
+  for (int y : Utils::make_lin_space(int(0), int(m_grid_dimensions[1]),
+                                     int(m_grid_dimensions[1]), false)) {
+  for (int z : Utils::make_lin_space(int(0), int(m_grid_dimensions[2]),
+                                     int(m_grid_dimensions[2]), false)) {
+          res.push_back(
+              {Utils::Vector3i{x, y, z},
+               Utils::Vector3d::broadcast(.5) +Utils::Vector3d{double(x), double(y), double(z)}});
+        }
+      }
+    }
+  return res;
+}
+
 #endif

src/core/grid_based_algorithms/LbWalberla.hpp

@@ -130,7 +130,10 @@ class LbWalberla {
       Flag_field_t *flag_field = block->getData<Flag_field_t>(m_flag_field_id);
       Pdf_field_t *pdf_field = block->getData<Pdf_field_t>(m_pdf_field_id);
 
-      const uint8_t fluid = flag_field->registerFlag(Fluid_flag);
+      // const uint8_t fluid = flag_field->registerFlag(Fluid_flag);
+      const auto fluid = flag_field->flagExists(Fluid_flag)
+                             ? flag_field->getFlag(Fluid_flag)
+                             : flag_field->registerFlag(Fluid_flag);
 
       return new Boundary_handling_t(
           "boundary handling", flag_field, fluid,
@@ -155,7 +158,7 @@ class LbWalberla {
   void integrate();
   std::pair<Utils::Vector3d, Utils::Vector3d> get_local_domain() {
     // We only have one block per mpi rank
-    assert(m_blocks->begin()++ == m_blocks->end());
+    assert(++(m_blocks->begin()) == m_blocks->end());
 
     auto const ab = m_blocks->begin()->getAABB();
     return {to_vector3d(ab.min()), to_vector3d(ab.max())};
@@ -261,19 +264,22 @@ class LbWalberla {
       walberla::field::TrilinearFieldInterpolator<VelocityAdaptor,
                                                   Flag_field_t>;
 
-  void empty_flags_for_boundary(
-      std::shared_ptr<walberla::StructuredBlockForest> &blocks,
-      const walberla::BlockDataID &boundary_handling_id) {
+public:
+  std::vector<std::pair<Utils::Vector3i, Utils::Vector3d>>
+  node_indices_positions();
+  std::vector<std::pair<Utils::Vector3i, Utils::Vector3d>>
+  global_node_indices_positions();
+  void clear_boundaries() {
     using namespace walberla;
     const CellInterval &domain_bb_in_global_cell_coordinates =
-        blocks->getDomainCellBB();
-    for (auto block = blocks->begin(); block != blocks->end(); ++block) {
+        m_blocks->getDomainCellBB();
+    for (auto block = m_blocks->begin(); block != m_blocks->end(); ++block) {
 
       Boundary_handling_t *boundary_handling =
-          block->getData<Boundary_handling_t>(boundary_handling_id);
+          block->getData<Boundary_handling_t>(m_boundary_handling_id);
 
       CellInterval domain_bb(domain_bb_in_global_cell_coordinates);
-      blocks->transformGlobalToBlockLocalCellInterval(domain_bb, *block);
+      m_blocks->transformGlobalToBlockLocalCellInterval(domain_bb, *block);
 
       boundary_handling->fillWithDomain(domain_bb);
     }

src/core/grid_based_algorithms/lb_boundaries.cpp

@@ -34,6 +34,7 @@
 #include "grid_based_algorithms/lattice.hpp"
 #include "grid_based_algorithms/lb.hpp"
 #include "grid_based_algorithms/lb_interface.hpp"
+#include "grid_based_algorithms/lb_walberla_instance.hpp"
 #include "grid_based_algorithms/lbgpu.hpp"
 #include "lbboundaries/LBBoundary.hpp"
 
@@ -271,9 +272,52 @@ void lb_init_boundaries() {
       }
     }
 #endif
-  }
+  } else if (lattice_switch == ActiveLB::WALBERLA) {
+#ifdef LB_WALBERLA
+#if defined(LB_BOUNDARIES)
+    Utils::Vector3i offset;
+    int the_boundary = -1;
+
+    lb_walberla()->clear_boundaries();
+
+    auto const agrid = lb_lbfluid_get_agrid();
+
+    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;
+
+      int n = 0;
+      for (auto it = lbboundaries.begin(); it != lbboundaries.end();
+           ++it, ++n) {
+        double dist;
+        Utils::Vector3d tmp;
+        (**it).calc_dist(pos, &dist, tmp.data());
+
+        if (dist <= 0) {
+
+          // Set boundaries on the ghost layers
+          auto const grid = lb_walberla()->get_grid_dimensions();
+          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]};
+                if (lb_walberla()->set_node_velocity_at_boundary(
+                  shifted_index, (**it).velocity() / lb_lbfluid_get_lattice_speed())) {
+          printf("Boundary: %d, %d %d %d, %g %g %g\n", n, shifted_index[0], shifted_index[1],
+                 shifted_index[2], pos[0], pos[1], pos[2]);
+            }
+            }
+            break;
+          } // if dist <=0
+        } // loop over boundaries
+      } // Loop over cells
+    } // lattice switch is WALBERLA
+#endif
+#endif
 }
 
+
 Utils::Vector3d lbboundary_get_force(LBBoundary const *lbb) {
   Utils::Vector3d force{};
 #if defined(LB_BOUNDARIES) || defined(LB_BOUNDARIES_GPU)

src/core/unit_tests/LbWalberla.cpp

@@ -75,6 +75,58 @@ BOOST_AUTO_TEST_CASE(boundary) {
   }
 }
 
+// BOOST_AUTO_TEST_CASE(boundary_flow_single_node) {
+//  Vector3d vel = {0.2, 3.8, 0};
+//  LbWalberla lb = LbWalberla(viscosity, density, agrid, tau, box_dimensions,
+//                             node_grid, skin);
+//    Vector3i node{5,5,5};
+//    if (lb.node_in_local_domain(node)) {
+//      BOOST_CHECK(lb.set_node_velocity_at_boundary(node, vel));
+//    }
+//    for( int i=0;i<10;i++) {
+//      lb.integrate();
+//    }
+//    for (int j=0;j<9;j++) {
+//        auto v = lb.get_node_velocity(Vector3i{j,node[1],node[2]});
+//        if (v) {
+//          if (j!=node[0]) {
+//            BOOST_CHECK((*v)[0]>1E-4);
+//            BOOST_CHECK((*v)[1]>1E-4);
+//            BOOST_CHECK(fabs((*v)[2])<1E-12);
+//            printf("%d,%g %g %g\n",j,(*v)[0],(*v)[1],(*v)[2]);
+//          }
+//        }
+//  }
+//}
+
+BOOST_AUTO_TEST_CASE(boundary_flow_shear) {
+  Vector3d vel = {0.2, -0.3, 0};
+  LbWalberla lb = LbWalberla(viscosity, density, agrid, tau, box_dimensions,
+                             node_grid, skin);
+  for (int x = 1; x < grid_dimensions[0] - 1; x++) {
+    for (int y = 1; y < grid_dimensions[1] - 1; y++) {
+      Vector3i node{x, y, 1};
+      if (lb.node_in_local_domain(node)) {
+        BOOST_CHECK(lb.set_node_velocity_at_boundary(node, vel));
+      }
+      node[2] = grid_dimensions[2] - 4;
+      if (lb.node_in_local_domain(node)) {
+        BOOST_CHECK(lb.set_node_velocity_at_boundary(node, -vel));
+      }
+    }
+  }
+
+  for (int i = 0; i < 150; i++) {
+    lb.integrate();
+  }
+  for (int j = 0; j < grid_dimensions[2]; j++) {
+    auto v = lb.get_node_velocity(Vector3i{0, 0, j});
+    if (v) {
+      printf("%d,%g %g %g\n", j, (*v)[0], (*v)[1], (*v)[2]);
+    }
+  }
+}
+
 BOOST_AUTO_TEST_CASE(velocity) {
   LbWalberla lb = LbWalberla(viscosity, density, agrid, tau, box_dimensions,
                              node_grid, skin);

testsuite/python/lb_shear.py

@@ -29,8 +29,8 @@
 """
 
 
-AGRID = 0.6
-VISC = 3.2
+AGRID = 1 
+VISC = 5.2
 DENS = 2.3
 TIME_STEP = 0.02
 # Box size will be H +2 AGRID to make room for walls.
@@ -79,7 +79,7 @@ def shear_flow(x, t, nu, v, h, k_max):
     return v * u
 
 
-class LBShearCommon(object):
+class LBShearCommon:
 
     """Base class of the test that holds the test logic."""
     lbf = None
@@ -115,22 +115,26 @@ def check_profile(self, shear_plane_normal, shear_direction):
         self.system.lbboundaries.add(wall2)
 
         t0 = self.system.time
-        sample_points = int(H / AGRID - 1)
+        sample_points = np.arange(H/AGRID+2)
 
         for i in range(9):
             self.system.integrator.run(50)
 
             v_expected = shear_flow(
-                x=(np.arange(0, sample_points) + .5) * AGRID,
+                x=(sample_points -.5) * AGRID,
                                 t=self.system.time - t0,
                                 nu=VISC,
                                 v=SHEAR_VELOCITY,
                                 h=H,
                                 k_max=100)
-            for j in range(1, sample_points):
-                ind = np.max(((1, 1, 1), shear_plane_normal * j + 1), 0)
+
+            # We omit the boundary nodes themselves, which have undefined velocities
+            for j in np.array(sample_points,dtype=int)[1:-1]:
+                ind = (1, 1, 1)
                 ind = np.array(ind, dtype=int)
+                ind[np.argmax(shear_plane_normal)]=j
                 v_measured = self.lbf[ind[0], ind[1], ind[2]].velocity
+                print(j,v_measured,v_expected[j])
                 np.testing.assert_allclose(
                     np.copy(v_measured),
                     np.copy(v_expected[j]) * shear_direction, atol=3E-3)
@@ -183,6 +187,15 @@ class LBGPUShear(ut.TestCase, LBShearCommon):
     def setUp(self):
         self.lbf = espressomd.lb.LBFluidGPU(**LB_PARAMS)
 
+@utx.skipIfMissingFeatures(['LB_WALBERLA', 'LB_BOUNDARIES'])
+class LBWalberlaShear(ut.TestCase, LBShearCommon):
+
+    """Test for the Walberla implementation of the LB."""
+
+    def setUp(self):
+        self.lbf = espressomd.lb.LBFluidWalberla(**LB_PARAMS)
+
+
 
 if __name__ == '__main__':
     ut.main()