WIP: LB coupling

maintainer/benchmarks/lb.py

@@ -62,7 +62,7 @@
 args = parser.parse_args()
 
 # process and check arguments
-n_iterations = 30
+n_iterations = 200
 assert args.volume_fraction > 0, "--volume_fraction must be a positive number"
 assert args.volume_fraction < np.pi / (3 * np.sqrt(2)), \
     "--volume_fraction exceeds the physical limit of sphere packing (~0.74)"

maintainer/benchmarks/lj.py

@@ -44,7 +44,7 @@
 
 # process and check arguments
 measurement_steps = int(np.round(5e6 / args.particles_per_core, -2))
-n_iterations = 30
+n_iterations = 2000 
 assert args.volume_fraction > 0, "volume_fraction must be a positive number"
 assert args.volume_fraction < np.pi / (3 * np.sqrt(2)), \
     "volume_fraction exceeds the physical limit of sphere packing (~0.74)"
@@ -119,10 +119,10 @@
 system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=42)
 
 # tuning and equilibration
-min_skin = 0.2
+min_skin = 0.1
 max_skin = 1.0
 print("Tune skin: {:.3f}".format(system.cell_system.tune_skin(
-    min_skin=min_skin, max_skin=max_skin, tol=0.05, int_steps=100)))
+    min_skin=min_skin, max_skin=max_skin, tol=0.005, int_steps=300)))
 print("Equilibration")
 system.integrator.run(min(5 * measurement_steps, 60000))
 print("Tune skin: {:.3f}".format(system.cell_system.tune_skin(

src/core/forces.cpp

@@ -51,8 +51,10 @@
 #include "short_range_loop.hpp"
 #include "system/System.hpp"
 #include "thermostat.hpp"
+
 #include "thermostats/langevin_inline.hpp"
 #include "virtual_sites/relative.hpp"
+#include <optional>
 
 #include <utils/Vector.hpp>
 #include <utils/math/sqr.hpp>
@@ -184,6 +186,9 @@ void System::System::calculate_forces() {
 #endif
   // Use combined function instead of two separate calls
 
+  // Prepare LB coupling state (initiates velocity interpolation)
+  std::optional<LB::ParticleCouplingState> lb_coupling_state;
+
   auto const elc_kernel = coulomb.pair_force_elc_kernel();
   auto const coulomb_kernel = coulomb.pair_force_kernel();
   auto const dipoles_kernel = dipoles.pair_force_kernel();
@@ -347,11 +352,6 @@ void System::System::calculate_forces() {
   // Must be done here. Forces need to be ghost-communicated
   immersed_boundaries->volume_conservation(*cell_structure);
 
-  if (thermostat->lb and (propagation->used_propagations &
-                          PropagationMode::TRANS_LB_MOMENTUM_EXCHANGE)) {
-    lb_couple_particles();
-  }
-
 #ifdef CUDA
 #ifdef CALIPER
   CALI_MARK_BEGIN("copy_forces_from_GPU");
@@ -362,6 +362,13 @@ void System::System::calculate_forces() {
 #endif
 #endif // CUDA
 
+  // Apply LB forces after GPU copy (hides interpolation latency)
+  if (thermostat->lb and (propagation->used_propagations &
+                          PropagationMode::TRANS_LB_MOMENTUM_EXCHANGE)) {
+    lb_coupling_state = lb_prepare_particle_coupling();
+    lb_apply_particle_forces(*lb_coupling_state);
+  }
+
 #ifdef VIRTUAL_SITES_RELATIVE
   if (propagation->used_propagations &
       (PropagationMode::TRANS_VS_RELATIVE | PropagationMode::ROT_VS_RELATIVE)) {

src/core/lb/Solver.cpp

@@ -39,6 +39,7 @@
 #include <cassert>
 #include <cmath>
 #include <functional>
+#include <future>
 #include <limits>
 #include <memory>
 #include <optional>
@@ -277,6 +278,14 @@ std::vector<Utils::Vector3d> Solver::get_coupling_interpolated_velocities(
       *impl->solver);
 }
 
+std::future<std::vector<Utils::Vector3d>>
+Solver::get_coupling_interpolated_velocities_async(
+    std::vector<Utils::Vector3d> const &pos) const {
+  return std::async(std::launch::async, [this, pos]() {
+    return get_coupling_interpolated_velocities(pos);
+  });
+}
+
 void Solver::add_forces_at_pos(std::vector<Utils::Vector3d> const &pos,
                                std::vector<Utils::Vector3d> const &forces) {
   std::visit(

src/core/lb/Solver.hpp

@@ -28,6 +28,7 @@
 #include <cassert>
 #include <cmath>
 #include <functional>
+#include <future>
 #include <memory>
 #include <optional>
 #include <vector>
@@ -190,6 +191,10 @@ struct Solver : public System::Leaf<Solver> {
   std::vector<Utils::Vector3d> get_coupling_interpolated_velocities(
       std::vector<Utils::Vector3d> const &pos) const;
 
+  std::future<std::vector<Utils::Vector3d>>
+  get_coupling_interpolated_velocities_async(
+      std::vector<Utils::Vector3d> const &pos) const;
+
   void add_forces_at_pos(std::vector<Utils::Vector3d> const &pos,
                          std::vector<Utils::Vector3d> const &forces);