Refactor waLBerla bridge

maintainer/walberla_kernels/templates/Boundary.tmpl.h

@@ -40,6 +40,7 @@
 #include <field/FlagField.h>
 #include <core/debug/Debug.h>
 
+#include <cassert>
 #include <functional>
 #include <memory>
 #include <vector>
@@ -197,8 +198,8 @@ class {{class_name}}
         auto * flagField = block->getData< FlagField_T > ( flagFieldID );
         {{additional_data_handler.additional_field_data|indent(4)}}
 
-        if( !(flagField->flagExists(boundaryFlagUID) && flagField->flagExists(domainFlagUID) ))
-            return;
+        assert(flagField->flagExists(boundaryFlagUID and
+               flagField->flagExists(domainFlagUID));
 
         auto boundaryFlag = flagField->getFlag(boundaryFlagUID);
         auto domainFlag = flagField->getFlag(domainFlagUID);

src/script_interface/walberla/LBFluidNode.cpp

@@ -47,13 +47,12 @@ Variant LBFluidNode::do_call_method(std::string const &name,
   if (name == "set_velocity_at_boundary") {
     if (is_none(params.at("value"))) {
       m_lb_fluid->remove_node_from_boundary(m_index);
-      m_lb_fluid->ghost_communication();
     } else {
       auto const u =
           get_value<Utils::Vector3d>(params, "value") * m_conv_velocity;
       m_lb_fluid->set_node_velocity_at_boundary(m_index, u);
-      m_lb_fluid->ghost_communication();
     }
+    m_lb_fluid->ghost_communication();
     return {};
   }
   if (name == "get_velocity_at_boundary") {

src/walberla_bridge/include/walberla_bridge/utils/boundary_utils.hpp

@@ -21,6 +21,8 @@
 
 #include "walberla_utils.hpp"
 
+#include <walberla_bridge/LatticeWalberla.hpp>
+
 #include <utils/Vector.hpp>
 
 #include <cassert>
@@ -81,30 +83,29 @@ void set_boundary_from_grid(BoundaryModel &boundary,
                             std::vector<int> const &raster_flat,
                             std::vector<DataType> const &data_flat) {
 
+  auto const &conv = es2walberla<DataType, typename BoundaryModel::value_type>;
   auto const grid_size = lattice.get_grid_dimensions();
   auto const offset = lattice.get_local_grid_range().first;
   auto const gl = static_cast<int>(lattice.get_ghost_layers());
   assert(raster_flat.size() == Utils::product(grid_size));
   auto const n_y = grid_size[1];
   auto const n_z = grid_size[2];
-  auto const off_i = offset[0];
-  auto const off_j = offset[1];
-  auto const off_k = offset[2];
 
-  auto const &blocks = lattice.get_blocks();
-  for (auto block = blocks->begin(); block != blocks->end(); ++block) {
-    auto const [size_i, size_j, size_k] = boundary.block_dims(*block);
+  for (auto const &block : *lattice.get_blocks()) {
+    auto const [size_i, size_j, size_k] = boundary.block_dims(block);
     // Get field data which knows about the indices
-    // In the loop, x,y,z are in block-local coordinates
-    for (int i = off_i - gl; i < size_i + off_i + gl; ++i) {
-      for (int j = off_j - gl; j < size_j + off_j + gl; ++j) {
-        for (int k = off_k - gl; k < size_k + off_k + gl; ++k) {
-          auto const node = Utils::Vector3i{{i, j, k}};
+    // In the loop, i,j,k are in block-local coordinates
+    for (int i = -gl; i < size_i + gl; ++i) {
+      for (int j = -gl; j < size_j + gl; ++j) {
+        for (int k = -gl; k < size_k + gl; ++k) {
+          auto const node = offset + Utils::Vector3i{{i, j, k}};
           auto const idx = (node + grid_size) % grid_size;
           auto const index = idx[0] * n_y * n_z + idx[1] * n_z + idx[2];
           if (raster_flat[index]) {
+            auto const &value = data_flat[index];
             auto const bc = get_block_and_cell(lattice, node, true);
-            boundary.set_node_value_at_boundary(node, data_flat[index], *bc);
+            assert(bc.has_value());
+            boundary.set_node_value_at_boundary(node, conv(value), *bc);
           }
         }
       }

src/walberla_bridge/src/BoundaryHandling.hpp

@@ -23,7 +23,10 @@
 #include <walberla_bridge/utils/walberla_utils.hpp>
 
 #include <blockforest/StructuredBlockForest.h>
+#include <domain_decomposition/BlockDataID.h>
+#include <domain_decomposition/IBlock.h>
 #include <field/FlagField.h>
+#include <field/FlagUID.h>
 
 #include <utils/Vector.hpp>
 
@@ -57,10 +60,9 @@ template <typename T, typename BoundaryClass> class BoundaryHandling {
       return get_value(global);
     }
 
-    template <typename U>
-    void set_node_boundary_value(Utils::Vector3i const &node, U const &val) {
+    void set_node_boundary_value(Utils::Vector3i const &node, T const &val) {
       auto const global = Cell(node[0], node[1], node[2]);
-      (*m_value_boundary)[global] = es2walberla<U, T>(val);
+      (*m_value_boundary)[global] = val;
     }
 
     void unset_node_boundary_value(Utils::Vector3i const &node) {
@@ -72,7 +74,7 @@ template <typename T, typename BoundaryClass> class BoundaryHandling {
     [[nodiscard]] auto
     get_node_boundary_value(Utils::Vector3i const &node) const {
       auto const global = Cell(node[0], node[1], node[2]);
-      return walberla2es(get_value(global));
+      return get_value(global);
     }
 
     bool node_is_boundary(Utils::Vector3i const &node) const {
@@ -105,17 +107,16 @@ template <typename T, typename BoundaryClass> class BoundaryHandling {
 
   BoundaryHandling(std::shared_ptr<StructuredBlockForest> blocks,
                    BlockDataID value_field_id, BlockDataID flag_field_id)
-      : m_blocks(std::move(blocks)), m_value_field_id(value_field_id),
-        m_flag_field_id(flag_field_id), m_callback(DynamicValueCallback()),
-        m_pending_changes(false) {
+      : m_blocks(std::move(blocks)), m_flag_field_id(flag_field_id),
+        m_callback(DynamicValueCallback()), m_pending_changes(false) {
     // reinitialize the flag field
     for (auto b = m_blocks->begin(); b != m_blocks->end(); ++b) {
       flag_reset_kernel(&*b);
     }
     // instantiate the boundary sweep
     std::function callback = m_callback;
     m_boundary =
-        std::make_shared<BoundaryClass>(m_blocks, m_value_field_id, callback);
+        std::make_shared<BoundaryClass>(m_blocks, value_field_id, callback);
   }
 
   void operator()(IBlock *block) { (*m_boundary)(block); }
@@ -129,15 +130,18 @@ template <typename T, typename BoundaryClass> class BoundaryHandling {
     return m_callback.get_node_boundary_value(node);
   }
 
-  template <typename U>
-  void set_node_value_at_boundary(Utils::Vector3i const &node, U const &v,
+  void set_node_value_at_boundary(Utils::Vector3i const &node, T const &v,
                                   BlockAndCell const &bc) {
     auto [flag_field, boundary_flag] = get_flag_field_and_flag(bc.block);
     m_callback.set_node_boundary_value(node, v);
     flag_field->addFlag(bc.cell, boundary_flag);
     m_pending_changes = true;
   }
 
+  void unpack_node(Utils::Vector3i const &node, T const &v) {
+    m_callback.set_node_boundary_value(node, v);
+  }
+
   void remove_node_from_boundary(Utils::Vector3i const &node,
                                  BlockAndCell const &bc) {
     auto [flag_field, boundary_flag] = get_flag_field_and_flag(bc.block);
@@ -163,7 +167,6 @@ template <typename T, typename BoundaryClass> class BoundaryHandling {
 
 private:
   std::shared_ptr<StructuredBlockForest> m_blocks;
-  BlockDataID m_value_field_id;
   BlockDataID m_flag_field_id;
   DynamicValueCallback m_callback;
   std::shared_ptr<BoundaryClass> m_boundary;

src/walberla_bridge/src/electrokinetics/EKinWalberlaImpl.hpp

@@ -195,7 +195,7 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
             m_density_field_id));
 
     // Synchronize ghost layers
-    (*m_full_communication)();
+    ghost_communication();
   }
 
   // Global parameters
@@ -453,7 +453,8 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
     if (!bc)
       return false;
 
-    m_boundary_flux->set_node_value_at_boundary(node, flux, *bc);
+    m_boundary_flux->set_node_value_at_boundary(
+        node, to_vector3<FloatType>(flux), *bc);
 
     return true;
   }
@@ -465,7 +466,7 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
     if (!bc or !m_boundary_flux->node_is_boundary(node))
       return std::nullopt;
 
-    return {m_boundary_flux->get_node_value_at_boundary(node)};
+    return {to_vector3d(m_boundary_flux->get_node_value_at_boundary(node))};
   }
 
   bool remove_node_from_flux_boundary(Utils::Vector3i const &node) override {
@@ -517,7 +518,8 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
             auto const bc = get_block_and_cell(lattice, node, false);
             auto const &opt = *it;
             if (opt) {
-              m_boundary_density->set_node_value_at_boundary(node, *opt, *bc);
+              m_boundary_density->set_node_value_at_boundary(
+                  node, FloatType_c(*opt), *bc);
             } else {
               m_boundary_density->remove_node_from_boundary(node, *bc);
             }
@@ -545,8 +547,8 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
           for (auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
             auto const node = local_offset + Utils::Vector3i{{x, y, z}};
             if (m_boundary_density->node_is_boundary(node)) {
-              out.emplace_back(
-                  m_boundary_density->get_node_value_at_boundary(node));
+              out.emplace_back(double_c(
+                  m_boundary_density->get_node_value_at_boundary(node)));
             } else {
               out.emplace_back(std::nullopt);
             }
@@ -575,7 +577,8 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
             auto const bc = get_block_and_cell(lattice, node, false);
             auto const &opt = *it;
             if (opt) {
-              m_boundary_flux->set_node_value_at_boundary(node, *opt, *bc);
+              m_boundary_flux->set_node_value_at_boundary(
+                  node, to_vector3<FloatType>(*opt), *bc);
             } else {
               m_boundary_flux->remove_node_from_boundary(node, *bc);
             }
@@ -603,8 +606,8 @@ class EKinWalberlaImpl : public EKinWalberlaBase {
           for (auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
             auto const node = local_offset + Utils::Vector3i{{x, y, z}};
             if (m_boundary_flux->node_is_boundary(node)) {
-              out.emplace_back(
-                  m_boundary_flux->get_node_value_at_boundary(node));
+              out.emplace_back(to_vector3d(
+                  m_boundary_flux->get_node_value_at_boundary(node)));
             } else {
               out.emplace_back(std::nullopt);
             }

src/walberla_bridge/src/electrokinetics/reactions/EKReactionImplIndexed.cpp

@@ -30,6 +30,7 @@
 #include <domain_decomposition/IBlock.h>
 #include <field/AddToStorage.h>
 
+#include <cassert>
 #include <cstddef>
 #include <memory>
 #include <optional>
@@ -68,9 +69,8 @@ void fillFromFlagField(IBlock *block, BlockDataID indexVectorID,
 
   auto *flagField = block->getData<FlagField>(flagFieldID);
 
-  if (!(flagField->flagExists(boundaryFlagUID) &&
-        flagField->flagExists(domainFlagUID)))
-    return;
+  assert(flagField->flagExists(boundaryFlagUID) and
+         flagField->flagExists(domainFlagUID));
 
   auto boundaryFlag = flagField->getFlag(boundaryFlagUID);
   auto domainFlag = flagField->getFlag(domainFlagUID);

src/walberla_bridge/src/lattice_boltzmann/LBWalberlaImpl.hpp

@@ -973,18 +973,17 @@ class LBWalberlaImpl : public LBWalberlaBase {
     if (!bc or !m_boundary->node_is_boundary(node))
       return std::nullopt;
 
-    return {m_boundary->get_node_value_at_boundary(node)};
+    return {to_vector3d(m_boundary->get_node_value_at_boundary(node))};
   }
 
   bool set_node_velocity_at_boundary(Utils::Vector3i const &node,
                                      Utils::Vector3d const &velocity) override {
     auto bc = get_block_and_cell(get_lattice(), node, true);
-    if (!bc)
-      return false;
-
-    m_boundary->set_node_value_at_boundary(node, velocity, *bc);
-
-    return true;
+    if (bc) {
+      m_boundary->set_node_value_at_boundary(
+          node, to_vector3<FloatType>(velocity), *bc);
+    }
+    return bc.has_value();
   }
 
   std::vector<std::optional<Utils::Vector3d>> get_slice_velocity_at_boundary(
@@ -1003,7 +1002,8 @@ class LBWalberlaImpl : public LBWalberlaBase {
           for (auto z = lower_cell.z(); z <= upper_cell.z(); ++z) {
             auto const node = local_offset + Utils::Vector3i{{x, y, z}};
             if (m_boundary->node_is_boundary(node)) {
-              out.emplace_back(m_boundary->get_node_value_at_boundary(node));
+              out.emplace_back(
+                  to_vector3d(m_boundary->get_node_value_at_boundary(node)));
             } else {
               out.emplace_back(std::nullopt);
             }
@@ -1032,7 +1032,8 @@ class LBWalberlaImpl : public LBWalberlaBase {
             auto const bc = get_block_and_cell(lattice, node, false);
             auto const &opt = *it;
             if (opt) {
-              m_boundary->set_node_value_at_boundary(node, *opt, *bc);
+              m_boundary->set_node_value_at_boundary(
+                  node, to_vector3<FloatType>(*opt), *bc);
             } else {
               m_boundary->remove_node_from_boundary(node, *bc);
             }
@@ -1054,12 +1055,10 @@ class LBWalberlaImpl : public LBWalberlaBase {
 
   bool remove_node_from_boundary(Utils::Vector3i const &node) override {
     auto bc = get_block_and_cell(get_lattice(), node, true);
-    if (!bc)
-      return false;
-
-    m_boundary->remove_node_from_boundary(node, *bc);
-
-    return true;
+    if (bc) {
+      m_boundary->remove_node_from_boundary(node, *bc);
+    }
+    return bc.has_value();
   }
 
   std::optional<bool>

src/walberla_bridge/src/lattice_boltzmann/generated_kernels/Dynamic_UBB_double_precision.h

@@ -31,6 +31,7 @@
 #include "field/FlagField.h"
 #include "field/GhostLayerField.h"
 
+#include <cassert>
 #include <functional>
 #include <memory>
 #include <vector>
@@ -138,9 +139,8 @@ class Dynamic_UBB_double_precision {
 
     auto *flagField = block->getData<FlagField_T>(flagFieldID);
 
-    if (!(flagField->flagExists(boundaryFlagUID) &&
-          flagField->flagExists(domainFlagUID)))
-      return;
+    assert(flagField->flagExists(boundaryFlagUID) and
+           flagField->flagExists(domainFlagUID));
 
     auto boundaryFlag = flagField->getFlag(boundaryFlagUID);
     auto domainFlag = flagField->getFlag(domainFlagUID);

src/walberla_bridge/src/lattice_boltzmann/generated_kernels/Dynamic_UBB_single_precision.h

@@ -31,6 +31,7 @@
 #include "field/FlagField.h"
 #include "field/GhostLayerField.h"
 
+#include <cassert>
 #include <functional>
 #include <memory>
 #include <vector>
@@ -138,9 +139,8 @@ class Dynamic_UBB_single_precision {
 
     auto *flagField = block->getData<FlagField_T>(flagFieldID);
 
-    if (!(flagField->flagExists(boundaryFlagUID) &&
-          flagField->flagExists(domainFlagUID)))
-      return;
+    assert(flagField->flagExists(boundaryFlagUID) and
+           flagField->flagExists(domainFlagUID));
 
     auto boundaryFlag = flagField->getFlag(boundaryFlagUID);
     auto domainFlag = flagField->getFlag(domainFlagUID);

src/walberla_bridge/tests/LBWalberlaImpl_unit_tests.cpp

@@ -435,7 +435,7 @@ BOOST_DATA_TEST_CASE(forces_interpolation, bdata::make(all_lbs()),
   // todo: check a less symmetrical situation, where the force is applied not
   // in the middle between the nodes
 
-  for (Vector3i n : all_nodes_incl_ghosts(lb->get_lattice())) {
+  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}};