Interface Walberla LB

maintainer/CI/build_cmake.sh

@@ -101,6 +101,11 @@ fi
 if $with_walberla; then
   cmake_params="$cmake_params -DWITH_WALBERLA=ON"
 fi
+if [[ "$hide_gpu" == "true" ]]
+then
+  echo Hiding gpu from Cuda via CUDA_VISIBLE_DEVICES
+  export CUDA_VISIBLE_DEVICES=
+fi
 
 if $insource; then
     builddir=$srcdir

src/core/grid_based_algorithms/LbWalberla.cpp

@@ -44,21 +44,21 @@
 
 using namespace walberla;
 
-inline Vector3d to_vector3d(const Vector3<real_t> v) {
-  return Vector3d{v[0], v[1], v[2]};
+inline Utils::Vector3d to_vector3d(const Vector3<real_t> v) {
+  return Utils::Vector3d{v[0], v[1], v[2]};
 }
-Vector3<real_t> to_vector3(const Vector3d v) {
+Vector3<real_t> to_vector3(const Utils::Vector3d v) {
   return Vector3<real_t>{v[0], v[1], v[2]};
 }
 
 LbWalberla::LbWalberla(double viscosity, double agrid,
-                       const Vector3d &box_dimensions,
-                       const Vector3i &node_grid, double skin) {
+                       const Utils::Vector3d &box_dimensions,
+                       const Utils::Vector3i &node_grid, double skin) {
 
   m_skin = skin;
   m_agrid = agrid;
 
-  Vector3i grid_dimensions;
+  Utils::Vector3i grid_dimensions;
   for (int i = 0; i < 3; i++) {
     if (fmod(box_dimensions[i], agrid) >
         std::numeric_limits<double>::epsilon()) {
@@ -149,7 +149,7 @@ void LbWalberla::print_vtk_density(char *filename) {
 void LbWalberla::integrate() { m_time_loop->run(); }
 
 boost::optional<LbWalberla::BlockAndCell>
-LbWalberla::get_block_and_cell(const Vector3i &node) const {
+LbWalberla::get_block_and_cell(const Utils::Vector3i &node) const {
   // Get block and local cell
   Cell global_cell{uint_c(node[0]), uint_c(node[1]), uint_c(node[2])};
   auto block = m_blocks->getBlock(global_cell, 0);
@@ -163,8 +163,8 @@ LbWalberla::get_block_and_cell(const Vector3i &node) const {
   return {{block, local_cell}};
 }
 
-bool LbWalberla::set_node_velocity_at_boundary(const Vector3i node,
-                                               const Vector3d v) {
+bool LbWalberla::set_node_velocity_at_boundary(const Utils::Vector3i node,
+                                               const Utils::Vector3d v) {
   auto bc = get_block_and_cell(node);
   // Return if we don't have the cell.
   if (!bc)
@@ -181,9 +181,9 @@ bool LbWalberla::set_node_velocity_at_boundary(const Vector3i node,
   return true;
 }
 
-boost::optional<Vector3d>
-LbWalberla::get_node_velocity_at_boundary(const Vector3i &node) const {
-  boost::optional<Vector3d> res;
+boost::optional<Utils::Vector3d>
+LbWalberla::get_node_velocity_at_boundary(const Utils::Vector3i &node) const {
+  boost::optional<Utils::Vector3d> res;
   auto bc = get_block_and_cell(node);
   // return if we don't have the cell
   if (!bc)
@@ -196,14 +196,14 @@ LbWalberla::get_node_velocity_at_boundary(const Vector3i &node) const {
   if (!boundary_handling->isBoundary((*bc).cell))
     return res;
 
-  Vector3d v =
+  Utils::Vector3d v =
       to_vector3d(boundary_handling->getBoundaryCondition<UBB_t>(uid).getValue(
           (*bc).cell[0], (*bc).cell[1], (*bc).cell[2]));
   res = {v};
   return res;
 }
 
-bool LbWalberla::remove_node_from_boundary(const Vector3i &node) {
+bool LbWalberla::remove_node_from_boundary(const Utils::Vector3i &node) {
   auto bc = get_block_and_cell(node);
   if (!bc)
     return false;
@@ -215,7 +215,7 @@ bool LbWalberla::remove_node_from_boundary(const Vector3i &node) {
 }
 
 boost::optional<bool>
-LbWalberla::get_node_is_boundary(const Vector3i &node) const {
+LbWalberla::get_node_is_boundary(const Utils::Vector3i &node) const {
   auto bc = get_block_and_cell(node);
   if (!bc)
     return {boost::none};
@@ -256,8 +256,8 @@ LbWalberla::create_fluid_field_vtk_writer(
   return pdf_field_vtk_writer;
 }
 
-boost::optional<Vector3d>
-LbWalberla::get_node_velocity(const Vector3i node) const {
+boost::optional<Utils::Vector3d>
+LbWalberla::get_node_velocity(const Utils::Vector3i node) const {
   auto bc = get_block_and_cell(node);
   if (!bc)
     return {boost::none};
@@ -268,8 +268,8 @@ LbWalberla::get_node_velocity(const Vector3i node) const {
   return {to_vector3d(vel_adaptor->get((*bc).cell))};
 }
 
-boost::optional<Vector3d>
-LbWalberla::get_velocity_at_pos(const Vector3d &pos) const {
+boost::optional<Utils::Vector3d>
+LbWalberla::get_velocity_at_pos(const Utils::Vector3d &pos) const {
   auto block = m_blocks->getBlock(to_vector3(pos));
   if (!block)
     return {boost::none};
@@ -281,7 +281,8 @@ LbWalberla::get_velocity_at_pos(const Vector3d &pos) const {
   return {to_vector3d(v)};
 }
 
-bool LbWalberla::set_node_velocity(const Vector3i &node, const Vector3d v) {
+bool LbWalberla::set_node_velocity(const Utils::Vector3i &node,
+                                   const Utils::Vector3d v) {
   auto bc = get_block_and_cell(node);
   if (!bc)
     return false;
@@ -294,7 +295,8 @@ bool LbWalberla::set_node_velocity(const Vector3i &node, const Vector3d v) {
   return true;
 }
 
-bool LbWalberla::set_node_density(const Vector3i node, const double density) {
+bool LbWalberla::set_node_density(const Utils::Vector3i node,
+                                  const double density) {
   auto bc = get_block_and_cell(node);
   if (!bc)
     return false;
@@ -311,7 +313,7 @@ bool LbWalberla::set_node_density(const Vector3i node, const double density) {
 }
 
 boost::optional<double>
-LbWalberla::get_node_density(const Vector3i node) const {
+LbWalberla::get_node_density(const Utils::Vector3i node) const {
   auto bc = get_block_and_cell(node);
   if (!bc)
     return {boost::none};
@@ -329,13 +331,13 @@ void LbWalberla::set_viscosity(double viscosity) {
 double LbWalberla::get_viscosity() {
   return m_lattice_model->collisionModel().viscosity();
 }
-bool LbWalberla::node_in_local_domain(const Vector3i &node) const {
+bool LbWalberla::node_in_local_domain(const Utils::Vector3i &node) const {
   auto block = m_blocks->getBlock(
       Cell{uint_c(node[0]), uint_c(node[1]), uint_c(node[2])}, 0);
   return (block != nullptr);
 }
 
-bool LbWalberla::pos_in_local_domain(const Vector3d &pos) const {
+bool LbWalberla::pos_in_local_domain(const Utils::Vector3d &pos) const {
   auto block =
       m_blocks->getBlock(real_c(pos[0]), real_c(pos[1]), real_c(pos[2]));
   return (block != nullptr);

src/core/grid_based_algorithms/LbWalberla.hpp

@@ -9,6 +9,7 @@
 #undef PI
 #include "blockforest/StructuredBlockForest.h"
 #include "boost/optional.hpp"
+#include "boost/tuple/tuple.hpp"
 #include "boundary/BoundaryHandling.h"
 #include "core/mpi/Environment.h"
 #include "field/FlagField.h"
@@ -30,7 +31,7 @@ const walberla::FlagUID UBB_flag("velocity bounce back");
 class LbWalberla {
   double m_skin;
   double m_agrid;
-  Vector3d m_ext_force;
+  Utils::Vector3d m_ext_force;
 
   // Type definitions
   using vector_field_t =
@@ -43,10 +44,8 @@ class LbWalberla {
   using Pdf_field_t = walberla::lbm::PdfField<Lattice_model_t>;
 
   using UBB_t = walberla::lbm::UBB<Lattice_model_t, walberla::uint8_t>;
-  using Boundary_conditions_t = boost::tuples::tuple<UBB_t>;
   using Boundary_handling_t =
-      walberla::BoundaryHandling<Flag_field_t, Lattice_model_t::Stencil,
-                                 Boundary_conditions_t>;
+      walberla::BoundaryHandling<Flag_field_t, Lattice_model_t::Stencil, UBB_t>;
   class LB_boundary_handling {
   public:
     LB_boundary_handling(const walberla::BlockDataID &flag_field_id,
@@ -62,8 +61,7 @@ class LbWalberla {
 
       return new Boundary_handling_t(
           "boundary handling", flag_field, fluid,
-          boost::tuples::make_tuple(
-              UBB_t("velocity bounce back", UBB_flag, pdf_field, nullptr)));
+          UBB_t("velocity bounce back", UBB_flag, pdf_field, nullptr));
     }
 
   private:
@@ -77,50 +75,57 @@ class LbWalberla {
   };
 
 public:
-  LbWalberla(double viscosity, double agrid, const Vector3d &box_dimensions,
-             const Vector3i &node_grid, double skin);
+  LbWalberla(double viscosity, double agrid,
+             const Utils::Vector3d &box_dimensions,
+             const Utils::Vector3i &node_grid, double skin);
 
   void integrate();
-  boost::optional<Vector3d> get_node_velocity(const Vector3i node) const;
-  bool set_node_velocity(const Vector3i &node, const Vector3d v);
-
-  bool add_force_at_pos(const Vector3d &position, const Vector3d &force);
-  boost::optional<Vector3d> get_velocity_at_pos(const Vector3d &position) const;
-  boost::optional<double> get_density_at_pos(const Vector3d &position);
-  //  Vector3d get_stress_at_pos(const Vector3d& position);
-
-  boost::optional<double> get_node_density(const Vector3i node) const;
-  bool set_node_density(const Vector3i node, const double density);
-  boost::optional<Vector3d>
-  get_node_velocity_at_boundary(const Vector3i &node) const;
-  bool set_node_velocity_at_boundary(const Vector3i node, const Vector3d v);
-  bool remove_node_from_boundary(const Vector3i &node);
-  boost::optional<bool> get_node_is_boundary(const Vector3i &node) const;
+  boost::optional<Utils::Vector3d>
+  get_node_velocity(const Utils::Vector3i node) const;
+  bool set_node_velocity(const Utils::Vector3i &node, const Utils::Vector3d v);
+
+  bool add_force_at_pos(const Utils::Vector3d &position,
+                        const Utils::Vector3d &force);
+  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);
+
+  boost::optional<double> get_node_density(const Utils::Vector3i node) const;
+  bool set_node_density(const Utils::Vector3i node, const double density);
+  boost::optional<Utils::Vector3d>
+  get_node_velocity_at_boundary(const Utils::Vector3i &node) const;
+  bool set_node_velocity_at_boundary(const Utils::Vector3i node,
+                                     const Utils::Vector3d v);
+  bool remove_node_from_boundary(const Utils::Vector3i &node);
+  boost::optional<bool> get_node_is_boundary(const Utils::Vector3i &node) const;
 
   void print_vtk_velocity(char *filename);
   void print_vtk_density(char *filename);
   void print_vtk_boundary(char *filename);
 
-  void set_external_force(const Vector3d &ext_force) {
+  void set_external_force(const Utils::Vector3d &ext_force) {
     m_ext_force = ext_force;
   }
-  Vector3d get_external_force() { return m_ext_force; }
+  Utils::Vector3d get_external_force() { return m_ext_force; }
 
   void set_viscosity(double viscosity);
   double get_viscosity();
 
-  Vector3i get_grid_dimensions() {
-    return Vector3i{{int(m_blocks->getXSize()), int(m_blocks->getYSize()),
-                     int(m_blocks->getZSize())}};
+  Utils::Vector3i get_grid_dimensions() {
+    return Utils::Vector3i{{int(m_blocks->getXSize()),
+                            int(m_blocks->getYSize()),
+                            int(m_blocks->getZSize())}};
   };
 
   double get_grid_spacing() { return m_agrid; }
 
-  bool node_in_local_domain(const Vector3i &node) const;
-  bool pos_in_local_domain(const Vector3d &pos) const;
+  bool node_in_local_domain(const Utils::Vector3i &node) const;
+  bool pos_in_local_domain(const Utils::Vector3d &pos) const;
 
 private:
-  boost::optional<BlockAndCell> get_block_and_cell(const Vector3i &node) const;
+  boost::optional<BlockAndCell>
+  get_block_and_cell(const Utils::Vector3i &node) const;
   walberla::BlockDataID m_pdf_field_id;
   walberla::BlockDataID m_flag_field_id;
   walberla::BlockDataID m_force_field_id;

src/core/grid_based_algorithms/lb_interface.cpp

@@ -1,12 +1,10 @@
 #include "lb_interface.hpp"
-#include "LbWalberla.hpp"
 #include "communication.hpp"
 #include "config.hpp"
 #include "electrokinetics.hpp"
 #include "global.hpp"
 #include "grid.hpp"
 #include "lb.hpp"
-#include "lb_walberla_instance.hpp"
 #include "lbgpu.hpp"
 
 #include "utils/index.hpp"
@@ -174,7 +172,7 @@ void mpi_recv_fluid_populations(int node, int index, double *pop) {
 } // namespace
 #endif
 
-#if defined(LB) || defined(LB_GPU) || defined(LB_WALBERLA)
+#if defined(LB) || defined(LB_GPU)
 void lb_lbfluid_update() {
   if (lattice_switch == ActiveLB::CPU) {
 #ifdef LB
@@ -210,55 +208,19 @@ void lb_lbfluid_propagate() {
   }
 }
 
-/**
- * @brief Check the boundary velocities.
- * Sanity check if the velocity defined at LB boundaries is within the Mach
- * number limits of the scheme i.e. u < 0.3.
- */
-void lb_boundary_mach_check() {
-  // Boundary velocities are stored in MD units, therefore we need to scale them
-  // in order to get lattice units.
-  auto const conv_fac = lb_lbfluid_get_tau() / lb_lbfluid_get_agrid();
-  double constexpr mach_limit = 0.3;
-  using LBBoundaries::lbboundaries;
-  if (std::any_of(lbboundaries.begin(), lbboundaries.end(),
-                  [conv_fac, mach_limit](auto const &b) {
-                    return (b->velocity() * conv_fac).norm() >= mach_limit;
-                  })) {
-    runtimeErrorMsg() << "Lattice velocity exceeds the Mach number limit";
-  }
-}
-
-/**
- * @brief Perform LB parameter and boundary velocity checks.
- */
-void lb_lbfluid_sanity_checks() {
-  if (lattice_switch == ActiveLB::GPU) {
-#ifdef LB_GPU
-    if (this_node == 0) {
-      lb_GPU_sanity_checks();
-      lb_boundary_mach_check();
-    }
-#endif
-  } else if (lattice_switch == ActiveLB::CPU) {
-#ifdef LB
-    lb_sanity_checks();
-    lb_boundary_mach_check();
-#endif
-  }
-}
-
 void lb_lbfluid_on_integration_start() {
-  lb_lbfluid_sanity_checks();
   if (lattice_switch == ActiveLB::GPU) {
 #ifdef LB_GPU
-    if (this_node == 0 and lb_reinit_particles_gpu()) {
+    lb_GPU_sanity_checks();
+    if (this_node == 0 && lb_reinit_particles_gpu()) {
       lb_realloc_particles_gpu();
       lb_reinit_particles_gpu.validate();
     }
 #endif
   } else if (lattice_switch == ActiveLB::CPU) {
 #ifdef LB
+    lb_sanity_checks();
+
     halo_communication(&update_halo_comm,
                        reinterpret_cast<char *>(lbfluid[0].data()));
 #endif
@@ -1283,9 +1245,8 @@ const Utils::Vector3d lb_lbnode_get_velocity(const Utils::Vector3i &ind) {
   return {};
 }
 
-#if defined(LB) || defined(LB_GPU) // lbmodel not defined for LB_WALBERLA
-const Vector6d lb_lbnode_get_pi(const Vector3i &ind) {
-  Vector6d p_pi = lb_lbnode_get_pi_neq(ind);
+const Utils::Vector6d lb_lbnode_get_pi(const Utils::Vector3i &ind) {
+  Utils::Vector6d p_pi = lb_lbnode_get_pi_neq(ind);
 
   // Add equilibrium stress to the diagonal (in LB units)
   double const p0 = lb_lbfluid_get_density() * lbmodel.c_sound_sq;
@@ -1296,7 +1257,6 @@ const Vector6d lb_lbnode_get_pi(const Vector3i &ind) {
 
   return p_pi;
 }
-#endif
 
 const Utils::Vector6d lb_lbnode_get_pi_neq(const Utils::Vector3i &ind) {
   Utils::Vector6d p_pi{};
@@ -1587,5 +1547,4 @@ Utils::Vector3d lb_lbfluid_calc_fluid_momentum() {
   }
   return fluid_momentum;
 }
-
-#endif // end of LB guard
+#endif