WIP: Re-introducing three point coupling

src/core/communication.cpp

@@ -159,7 +159,8 @@ int n_nodes = -1;
   CB(mpi_set_lb_fluid_counter)                                                 \
   CB(mpi_update_particle_slave)                                                \
   CB(mpi_bcast_lb_particle_coupling_slave)                                     \
-  CB(mpi_recv_lb_interpolated_velocity_slave)
+  CB(mpi_recv_lb_interpolated_velocity_slave)                                  \
+  CB(mpi_set_interpolation_order_slave)
 
 // create the forward declarations
 #define CB(name) void name(int node, int param);

src/core/grid_based_algorithms/lb_interpolation.cpp

@@ -1,17 +1,39 @@
+#include <boost/mpi/collectives.hpp>
+
 #include "communication.hpp"
 #include "config.hpp"
 #include "grid.hpp"
 #include "grid_based_algorithms/lattice.hpp"
+#include "grid_based_algorithms/lb_interpolation.hpp"
 #include "utils/Vector.hpp"
 
 #include "lb.hpp"
 #include "lb_interface.hpp"
 #include "lbgpu.hpp"
 
+namespace {
+
+InterpolationOrder interpolation_order = InterpolationOrder::linear;
+}
+
+void mpi_set_interpolation_order_slave(int, int) {
+  boost::mpi::broadcast(comm_cart, interpolation_order, 0);
+}
+
 #if defined(LB) || defined(LB_GPU)
 
-namespace {
+void lb_lbinterpolation_set_interpolation_order(
+    InterpolationOrder const &order) {
+  interpolation_order = order;
+  mpi_call(mpi_set_interpolation_order_slave, 0, 0);
+  boost::mpi::broadcast(comm_cart, interpolation_order, 0);
+}
 
+InterpolationOrder lb_lbinterpolation_get_interpolation_order() {
+  return interpolation_order;
+}
+
+namespace {
 template <typename Op>
 void lattice_interpolation(Lattice const &lattice, Vector3d const &pos,
                            Op &&op) {
@@ -72,17 +94,31 @@ lb_lbinterpolation_get_interpolated_velocity_global(const Vector3d &pos) {
   if (lattice_switch & LATTICE_LB_GPU) {
 #ifdef LB_GPU
     Vector3d interpolated_u{};
-    lb_get_interpolated_velocity_gpu(folded_pos.data(), interpolated_u.data(),
-                                     1);
+    switch (interpolation_order) {
+    case (InterpolationOrder::linear):
+      lb_get_interpolated_velocity_gpu<8>(folded_pos.data(),
+                                          interpolated_u.data(), 1);
+      break;
+    case (InterpolationOrder::quadratic):
+      lb_get_interpolated_velocity_gpu<27>(folded_pos.data(),
+                                           interpolated_u.data(), 1);
+      break;
+    }
     return interpolated_u;
 #endif
   } else if (lattice_switch & LATTICE_LB) {
 #ifdef LB
-    auto const node = map_position_node_array(folded_pos);
-    if (node == 0) {
-      return lb_lbinterpolation_get_interpolated_velocity(folded_pos);
-    } else {
-      return mpi_recv_lb_interpolated_velocity(node, folded_pos);
+    switch (interpolation_order) {
+    case (InterpolationOrder::quadratic):
+      throw std::runtime_error("The non-linear interpolation scheme is not "
+                               "implemented for the CPU LB.");
+    case (InterpolationOrder::linear):
+      auto const node = map_position_node_array(folded_pos);
+      if (node == 0) {
+        return lb_lbinterpolation_get_interpolated_velocity(folded_pos);
+      } else {
+        return mpi_recv_lb_interpolated_velocity(node, folded_pos);
+      }
     }
   }
 #endif
@@ -92,11 +128,18 @@ lb_lbinterpolation_get_interpolated_velocity_global(const Vector3d &pos) {
 #ifdef LB
 void lb_lbinterpolation_add_force_density(const Vector3d &pos,
                                           const Vector3d &force_density) {
-  lattice_interpolation(lblattice, pos,
-                        [&force_density](Lattice::index_t index, double w) {
-                          auto &field = lbfields[index];
-                          field.force_density += w * force_density;
-                        });
+  switch (interpolation_order) {
+  case (InterpolationOrder::quadratic):
+    throw std::runtime_error("The non-linear interpolation scheme is not "
+                             "implemented for the CPU LB.");
+  case (InterpolationOrder::linear):
+    lattice_interpolation(lblattice, pos,
+                          [&force_density](Lattice::index_t index, double w) {
+                            auto &field = lbfields[index];
+                            field.force_density += w * force_density;
+                          });
+    break;
+  }
 }
 #endif
 

src/core/grid_based_algorithms/lb_interpolation.hpp

@@ -1,6 +1,21 @@
 #ifndef LATTICE_INTERPOLATION_HPP
 #define LATTICE_INTERPOLATION_HPP
 
+#include "utils/Vector.hpp"
+/**
+ * @brief Interpolation order for the LB fluid interpolation.
+ * @note For the CPU LB only linear interpolation is available.
+ */
+enum class InterpolationOrder { linear, quadratic };
+
+/**
+ * @brief Set the interpolation order for the LB.
+ */
+void lb_lbinterpolation_set_interpolation_order(
+    InterpolationOrder const &interpolation_order);
+void mpi_set_interpolation_order_slave(int order, int);
+
+InterpolationOrder lb_lbinterpolation_get_interpolation_order();
 /**
  * @brief Calculates the fluid velocity at a given position of the
  * lattice.

src/core/grid_based_algorithms/lb_particle_coupling.cpp

@@ -168,71 +168,88 @@ void lb_lbcoupling_calc_particle_lattice_ia(bool couple_virtual) {
   ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
   if (lattice_switch & LATTICE_LB_GPU) {
 #ifdef LB_GPU
-    if (lb_particle_coupling.couple_to_md && this_node == 0)
-      lb_calc_particle_lattice_ia_gpu(couple_virtual,
-                                      lb_lbcoupling_get_gamma());
+    if (lb_particle_coupling.couple_to_md && this_node == 0) {
+      switch (lb_lbinterpolation_get_interpolation_order()) {
+      case (InterpolationOrder::linear):
+        lb_calc_particle_lattice_ia_gpu<8>(couple_virtual,
+                                           lb_lbcoupling_get_gamma());
+        break;
+      case (InterpolationOrder::quadratic):
+        lb_calc_particle_lattice_ia_gpu<27>(couple_virtual,
+                                            lb_lbcoupling_get_gamma());
+        break;
+      }
+    }
 #endif
   } else if (lattice_switch & LATTICE_LB) {
 #ifdef LB
     if (lb_particle_coupling.couple_to_md) {
-      using rng_type = r123::Philox4x64;
-      using ctr_type = rng_type::ctr_type;
-      using key_type = rng_type::key_type;
-
-      ctr_type c{{lb_particle_coupling.rng_counter_coupling.value(),
-                  static_cast<uint64_t>(RNGSalt::PARTICLES)}};
-
-      /* Eq. (16) Ahlrichs and Duenweg, JCP 111(17):8225 (1999).
-       * The factor 12 comes from the fact that we use random numbers
-       * from -0.5 to 0.5 (equally distributed) which have variance 1/12.
-       * time_step comes from the discretization.
-       */
-      auto const noise_amplitude = sqrt(12. * 2. * lb_lbcoupling_get_gamma() *
-                                        lb_lbfluid_get_kT() / time_step);
-      auto f_random = [&c](int id) -> Vector3d {
-        if (lb_lbfluid_get_kT() > 0.0) {
-          key_type k{{static_cast<uint32_t>(id)}};
-
-          auto const noise = rng_type{}(c, k);
-
-          using Utils::uniform;
-          return Vector3d{uniform(noise[0]), uniform(noise[1]),
-                          uniform(noise[2])} -
-                 Vector3d::broadcast(0.5);
-        } else {
-          return Vector3d{};
-        }
-      };
+      switch (lb_lbinterpolation_get_interpolation_order()) {
+      case (InterpolationOrder::quadratic):
+        throw std::runtime_error("The non-linear interpolation scheme is not "
+                                 "implemented for the CPU LB.");
+      case (InterpolationOrder::linear): {
+        using rng_type = r123::Philox4x64;
+        using ctr_type = rng_type::ctr_type;
+        using key_type = rng_type::key_type;
 
-      /* local cells */
-      for (auto &p : local_cells.particles()) {
-        if (!p.p.is_virtual or thermo_virtual or
-            (p.p.is_virtual && couple_virtual)) {
-          auto const force =
-              lb_viscous_coupling(&p, noise_amplitude * f_random(p.identity()));
-          /* add force to the particle */
-          p.f.f += force;
+        ctr_type c{{lb_particle_coupling.rng_counter_coupling.value(),
+                    static_cast<uint64_t>(RNGSalt::PARTICLES)}};
+
+        /* Eq. (16) Ahlrichs and Duenweg, JCP 111(17):8225 (1999).
+         * The factor 12 comes from the fact that we use random numbers
+         * from -0.5 to 0.5 (equally distributed) which have variance 1/12.
+         * time_step comes from the discretization.
+         */
+        auto const noise_amplitude = sqrt(12. * 2. * lb_lbcoupling_get_gamma() *
+                                          lb_lbfluid_get_kT() / time_step);
+        auto f_random = [&c](int id) -> Vector3d {
+          if (lb_lbfluid_get_kT() > 0.0) {
+            key_type k{{static_cast<uint32_t>(id)}};
+
+            auto const noise = rng_type{}(c, k);
+
+            using Utils::uniform;
+            return Vector3d{uniform(noise[0]), uniform(noise[1]),
+                            uniform(noise[2])} -
+                   Vector3d::broadcast(0.5);
+          } else {
+            return Vector3d{};
+          }
+        };
+
+        /* local cells */
+        for (auto &p : local_cells.particles()) {
+          if (!p.p.is_virtual or thermo_virtual or
+              (p.p.is_virtual && couple_virtual)) {
+            auto const force = lb_viscous_coupling(
+                &p, noise_amplitude * f_random(p.identity()));
+            /* add force to the particle */
+            p.f.f += force;
 #ifdef ENGINE
-          add_swimmer_force(p);
+            add_swimmer_force(p);
 #endif
+          }
         }
-      }
 
-      /* ghost cells */
-      for (auto &p : ghost_cells.particles()) {
-        /* for ghost particles we have to check if they lie
-         * in the range of the local lattice nodes */
-        if (in_local_domain(p.r.p)) {
-          if (!p.p.is_virtual || thermo_virtual) {
-            lb_viscous_coupling(&p, noise_amplitude * f_random(p.identity()));
+        /* ghost cells */
+        for (auto &p : ghost_cells.particles()) {
+          /* for ghost particles we have to check if they lie
+           * in the range of the local lattice nodes */
+          if (in_local_domain(p.r.p)) {
+            if (!p.p.is_virtual || thermo_virtual) {
+              lb_viscous_coupling(&p, noise_amplitude * f_random(p.identity()));
 #ifdef ENGINE
-            add_swimmer_force(p);
+              add_swimmer_force(p);
 #endif
+            }
           }
         }
+        break;
+      }
       }
-    }
 #endif
+    }
   }
 }
 

src/core/grid_based_algorithms/lbgpu.hpp

@@ -240,6 +240,8 @@ void lb_init_extern_nodeforcedensities_GPU(
     LB_parameters_gpu *lbpar_gpu);
 
 void lb_set_agrid_gpu(double agrid);
+
+template <std::size_t no_of_neighbours>
 void lb_calc_particle_lattice_ia_gpu(bool couple_virtual, double friction);
 
 void lb_calc_fluid_mass_GPU(double *mass);
@@ -273,8 +275,14 @@ void lb_lbfluid_particles_add_momentum(float velocity[3]);
 void lb_lbfluid_set_population(const Vector3i &, float[LBQ]);
 void lb_lbfluid_get_population(const Vector3i &, float[LBQ]);
 
+template <std::size_t no_of_neighbours>
 void lb_get_interpolated_velocity_gpu(double const *positions,
                                       double *velocities, int length);
+void linear_velocity_interpolation(double const *positions, double *velocities,
+                                   int length);
+void quadratic_velocity_interpolation(double const *positions,
+                                      double *velocities, int length);
+
 uint64_t lb_fluid_get_rng_state_gpu();
 void lb_fluid_set_rng_state_gpu(uint64_t counter);
 uint64_t lb_coupling_get_rng_state_gpu();