core: Remove time_step from LB code

src/core/event.cpp

@@ -110,12 +110,14 @@ void on_integration_start(double time_step) {
   integrator_npt_sanity_checks();
 #endif
   interactions_sanity_checks();
-  lb_lbfluid_on_integration_start();
+  lb_lbfluid_sanity_checks(time_step);
 
   /********************************************/
   /* end sanity checks                        */
   /********************************************/
 
+  lb_lbfluid_on_integration_start();
+
 #ifdef CUDA
   MPI_Bcast(gpu_get_global_particle_vars_pointer_host(),
             sizeof(CUDA_global_part_vars), MPI_BYTE, 0, comm_cart);

src/core/grid_based_algorithms/lb.cpp

@@ -35,7 +35,6 @@
 #include "grid.hpp"
 #include "grid_based_algorithms/lb_boundaries.hpp"
 #include "halo.hpp"
-#include "integrate.hpp"
 #include "lb-d3q19.hpp"
 #include "random.hpp"
 
@@ -983,13 +982,10 @@ void lb_collide_stream() {
  *  for the lattice dynamics can be coarser than the MD time step, we
  *  monitor the time since the last lattice update.
  */
-void lattice_boltzmann_update(double time_step) {
-  auto const factor = static_cast<int>(round(lbpar.tau / time_step));
-
+void lattice_boltzmann_update(int lb_steps_per_md_step) {
   fluidstep += 1;
-  if (fluidstep >= factor) {
+  if (fluidstep >= lb_steps_per_md_step) {
     fluidstep = 0;
-
     lb_collide_stream();
   }
 }

src/core/grid_based_algorithms/lb.hpp

@@ -183,7 +183,7 @@ extern std::vector<LB_FluidNode> lbfields;
  *  collisions. If boundaries are present, it also applies the boundary
  *  conditions.
  */
-void lattice_boltzmann_update(double time_step);
+void lattice_boltzmann_update(int lb_steps_per_md_step);
 
 void lb_sanity_checks(const LB_Parameters &lb_parameters);
 

src/core/grid_based_algorithms/lb_interface.cpp

@@ -52,27 +52,27 @@ struct NoLBActive : public std::exception {
   const char *what() const noexcept override { return "LB not activated"; }
 };
 
-void lb_lbfluid_update(double time_step) {
+void lb_lbfluid_update(int lb_steps_per_md_step) {
   if (lattice_switch == ActiveLB::CPU) {
-    lattice_boltzmann_update(time_step);
+    lattice_boltzmann_update(lb_steps_per_md_step);
   } else if (lattice_switch == ActiveLB::GPU and this_node == 0) {
 #ifdef CUDA
 #ifdef ELECTROKINETICS
     if (ek_initialized) {
       ek_integrate();
     } else {
 #endif
-      lattice_boltzmann_update_gpu(time_step);
+      lattice_boltzmann_update_gpu(lb_steps_per_md_step);
 #ifdef ELECTROKINETICS
     }
 #endif
 #endif
   }
 }
 
-void lb_lbfluid_propagate(double time_step) {
+void lb_lbfluid_propagate(int lb_steps_per_md_step) {
   if (lattice_switch != ActiveLB::NONE) {
-    lb_lbfluid_update(time_step);
+    lb_lbfluid_update(lb_steps_per_md_step);
     if (lb_lbfluid_get_kT() > 0.0) {
       if (lattice_switch == ActiveLB::GPU) {
 #ifdef CUDA
@@ -104,8 +104,7 @@ void lb_boundary_mach_check() {
   }
 }
 
-void lb_lbfluid_sanity_checks() {
-  extern double time_step;
+void lb_lbfluid_sanity_checks(double time_step) {
   if (lattice_switch == ActiveLB::GPU && this_node == 0) {
 #ifdef CUDA
     lb_GPU_sanity_checks();
@@ -123,7 +122,6 @@ void lb_lbfluid_sanity_checks() {
 }
 
 void lb_lbfluid_on_integration_start() {
-  lb_lbfluid_sanity_checks();
   if (lattice_switch == ActiveLB::CPU) {
     halo_communication(update_halo_comm,
                        reinterpret_cast<char *>(lbfluid[0].data()));
@@ -418,18 +416,18 @@ void lb_lbfluid_set_tau(double tau) {
   }
 }
 
-void check_tau_time_step_consistency(double tau, double time_s) {
+void check_tau_time_step_consistency(double tau, double time_step) {
   auto const eps = std::numeric_limits<float>::epsilon();
-  if ((tau - time_s) / (tau + time_s) < -eps)
+  if ((tau - time_step) / (tau + time_step) < -eps)
     throw std::invalid_argument("LB tau (" + std::to_string(tau) +
                                 ") must be >= MD time_step (" +
-                                std::to_string(time_s) + ")");
-  auto const factor = tau / time_s;
+                                std::to_string(time_step) + ")");
+  auto const factor = tau / time_step;
   if (fabs(round(factor) - factor) / factor > eps)
     throw std::invalid_argument("LB tau (" + std::to_string(tau) +
                                 ") must be integer multiple of "
                                 "MD time_step (" +
-                                std::to_string(time_s) + "). Factor is " +
+                                std::to_string(time_step) + "). Factor is " +
                                 std::to_string(factor));
 }
 

src/core/grid_based_algorithms/lb_interface.hpp

@@ -37,7 +37,7 @@ extern ActiveLB lattice_switch;
 /**
  * @brief Propagate the LB fluid.
  */
-void lb_lbfluid_propagate(double time_step);
+void lb_lbfluid_propagate(int lb_steps_per_md_step);
 
 /**
  * @brief Event handler for integration start.
@@ -133,7 +133,7 @@ void lb_lbfluid_set_kT(double kT);
 /**
  * @brief Perform LB parameter and boundary velocity checks.
  */
-void lb_lbfluid_sanity_checks();
+void lb_lbfluid_sanity_checks(double time_step);
 
 /**
  * @brief Set the LB density for a single node.

src/core/grid_based_algorithms/lb_particle_coupling.cpp

@@ -122,6 +122,7 @@ namespace {
  * @brief Add a force to the lattice force density.
  * @param pos Position of the force
  * @param force Force in MD units.
+ * @param time_step MD time step.
  */
 void add_md_force(Utils::Vector3d const &pos, Utils::Vector3d const &force,
                   double time_step) {
@@ -137,9 +138,10 @@ void add_md_force(Utils::Vector3d const &pos, Utils::Vector3d const &force,
  *  Section II.C. @cite ahlrichs99a
  *
  *  @param[in] p             The coupled particle.
- *  @param[in]     f_random  Additional force to be included.
+ *  @param[in] f_random      Additional force to be included.
+ *  @param[in] time_step     MD time step.
  *
- *  @return The viscous coupling force plus f_random.
+ *  @return The viscous coupling force plus @p f_random.
  */
 Utils::Vector3d lb_viscous_coupling(Particle const &p,
                                     Utils::Vector3d const &f_random,

src/core/grid_based_algorithms/lbgpu.cpp

@@ -100,14 +100,9 @@ bool ek_initialized = false;
 /*-----------------------------------------------------------*/
 
 /** %Lattice Boltzmann update gpu called from integrate.cpp */
-void lattice_boltzmann_update_gpu(double time_step) {
-
-  auto factor = static_cast<int>(round(lbpar_gpu.tau / time_step));
-
+void lattice_boltzmann_update_gpu(int lb_steps_per_md_step) {
   fluidstep += 1;
-
-  if (fluidstep >= factor) {
-
+  if (fluidstep >= lb_steps_per_md_step) {
     fluidstep = 0;
     lb_integrate_GPU();
   }

src/core/grid_based_algorithms/lbgpu.hpp

@@ -150,7 +150,7 @@ void lb_GPU_sanity_checks();
 void lb_get_device_values_pointer(LB_rho_v_gpu **pointer_address);
 void lb_get_boundary_force_pointer(float **pointer_address);
 void lb_get_para_pointer(LB_parameters_gpu **pointer_address);
-void lattice_boltzmann_update_gpu(double time_step);
+void lattice_boltzmann_update_gpu(int lb_steps_per_md_step);
 
 /** Perform a full initialization of the lattice Boltzmann system.
  *  All derived parameters and the fluid are reset to their default values.

src/core/integrate.cpp

@@ -57,6 +57,10 @@
 
 #include <boost/range/algorithm/min_element.hpp>
 
+#include <algorithm>
+#include <cmath>
+#include <csignal>
+#include <functional>
 #include <stdexcept>
 
 #ifdef VALGRIND_INSTRUMENTATION
@@ -284,8 +288,12 @@ int integrate(int n_steps, int reuse_forces) {
 
     // propagate one-step functionalities
     if (integ_switch != INTEG_METHOD_STEEPEST_DESCENT) {
-      lb_lbfluid_propagate(time_step);
-      lb_lbcoupling_propagate();
+      if (lb_lbfluid_get_lattice_switch() != ActiveLB::NONE) {
+        auto const lb_steps_per_md_step =
+            static_cast<int>(std::round(lb_lbfluid_get_tau() / time_step));
+        lb_lbfluid_propagate(lb_steps_per_md_step);
+        lb_lbcoupling_propagate();
+      }
 
 #ifdef VIRTUAL_SITES
       virtual_sites()->after_lb_propagation(time_step);

src/python/espressomd/lb.pxd

@@ -66,7 +66,6 @@ cdef extern from "grid_based_algorithms/lb_interface.hpp":
     const Vector3d lb_lbfluid_get_ext_force_density() except +
     void lb_lbfluid_set_bulk_viscosity(double c_bulk_visc) except +
     double lb_lbfluid_get_bulk_viscosity() except +
-    void lb_lbfluid_sanity_checks() except +
     void lb_lbfluid_print_vtk_velocity(string filename) except +
     void lb_lbfluid_print_vtk_velocity(string filename, vector[int] bb1, vector[int] bb2) except +
     void lb_lbfluid_print_vtk_boundary(string filename) except +

src/python/espressomd/lb.pyx

@@ -54,7 +54,7 @@ cdef class HydrodynamicInteraction(Actor):
         Lattice constant. The box size in every direction must be an integer
         multiple of ``agrid``.
     tau : :obj:`float`
-        LB time step. The MD time step must be an integer multiple of ``tau``.
+        LB time step, must be an integer multiple of the MD time step.
     dens : :obj:`float`
         Fluid density.
     visc : :obj:`float`
@@ -157,7 +157,6 @@ cdef class HydrodynamicInteraction(Actor):
         if "gamma_even" in self._params:
             python_lbfluid_set_gamma_even(self._params["gamma_even"])
 
-        lb_lbfluid_sanity_checks()
         utils.handle_errors("LB fluid activation")
 
     def _get_params_from_es_core(self):