core: Fix LB/EK propagation counter

src/core/ek/Solver.hpp

@@ -73,10 +73,6 @@ struct Solver {
    */
   double get_tau() const;
 
-  auto get_steps_per_md_step(double time_step) const {
-    return static_cast<int>(std::round(get_tau() / time_step));
-  }
-
   /**
    * @brief Perform @c EK parameter checks.
    */

src/core/integrate.cpp

@@ -103,8 +103,8 @@ bool recalc_forces = true;
 /** Average number of integration steps the Verlet list has been re-using. */
 static double verlet_reuse = 0.0;
 
-static int fluid_step = 0;
-static int ek_step = 0;
+static int lb_skipped_md_steps = 0;
+static int ek_skipped_md_steps = 0;
 
 namespace {
 volatile std::sig_atomic_t ctrl_C = 0;
@@ -238,6 +238,10 @@ void walberla_agrid_sanity_checks(std::string method,
 }
 #endif
 
+static auto calc_md_steps_per_tau(double tau) {
+  return static_cast<int>(std::round(tau / ::time_step));
+}
+
 static void resort_particles_if_needed(ParticleRange const &particles) {
   auto const offset = LeesEdwards::verlet_list_offset(
       box_geo, cell_structure.get_le_pos_offset_at_last_resort());
@@ -446,39 +450,40 @@ int integrate(int n_steps, int reuse_forces) {
         auto &lb = system.lb;
         auto &ek = system.ek;
         // assume that they are coupled, which is not necessarily true
-        auto const lb_steps_per_md_step = lb.get_steps_per_md_step(time_step);
-        auto const ek_steps_per_md_step = ek.get_steps_per_md_step(time_step);
+        auto const md_steps_per_lb_step = calc_md_steps_per_tau(lb.get_tau());
+        auto const md_steps_per_ek_step = calc_md_steps_per_tau(ek.get_tau());
 
-        if (lb_steps_per_md_step != ek_steps_per_md_step) {
+        if (md_steps_per_lb_step != md_steps_per_ek_step) {
           runtimeErrorMsg()
               << "LB and EK are active but with different time steps.";
         }
 
-        // only use fluid_step in this case
-        assert(fluid_step == ek_step);
+        assert(lb_skipped_md_steps == ek_skipped_md_steps);
 
-        fluid_step += 1;
-        if (fluid_step >= lb_steps_per_md_step) {
-          fluid_step = 0;
+        lb_skipped_md_steps += 1;
+        ek_skipped_md_steps += 1;
+        if (lb_skipped_md_steps >= md_steps_per_lb_step) {
+          lb_skipped_md_steps = 0;
+          ek_skipped_md_steps = 0;
           lb.propagate();
           ek.propagate();
         }
         lb_lbcoupling_propagate();
       } else if (lb_active) {
         auto &lb = system.lb;
-        auto const lb_steps_per_md_step = lb.get_steps_per_md_step(time_step);
-        fluid_step += 1;
-        if (fluid_step >= lb_steps_per_md_step) {
-          fluid_step = 0;
+        auto const md_steps_per_lb_step = calc_md_steps_per_tau(lb.get_tau());
+        lb_skipped_md_steps += 1;
+        if (lb_skipped_md_steps >= md_steps_per_lb_step) {
+          lb_skipped_md_steps = 0;
           lb.propagate();
         }
         lb_lbcoupling_propagate();
       } else if (ek_active) {
         auto &ek = system.ek;
-        auto const ek_steps_per_md_step = ek.get_steps_per_md_step(time_step);
-        ek_step += 1;
-        if (ek_step >= ek_steps_per_md_step) {
-          ek_step = 0;
+        auto const md_steps_per_ek_step = calc_md_steps_per_tau(ek.get_tau());
+        ek_skipped_md_steps += 1;
+        if (ek_skipped_md_steps >= md_steps_per_ek_step) {
+          ek_skipped_md_steps = 0;
           ek.propagate();
         }
       }

src/core/lb/Solver.hpp

@@ -142,10 +142,6 @@ struct Solver {
   void add_force_density(Utils::Vector3d const &pos,
                          Utils::Vector3d const &force_density);
 
-  auto get_steps_per_md_step(double time_step) const {
-    return static_cast<int>(std::round(get_tau() / time_step));
-  }
-
   void on_boxl_change();
   void on_node_grid_change();
   void on_cell_structure_change();

src/core/unit_tests/ek_interface_test.cpp

@@ -138,9 +138,6 @@ BOOST_AUTO_TEST_CASE(ek_interface_walberla) {
     // tau setters and getters
     espresso::ek_container->set_tau(2.);
     BOOST_CHECK_EQUAL(ek.get_tau(), 2.);
-    BOOST_CHECK_EQUAL(ek.get_steps_per_md_step(1.), 2);
-    BOOST_CHECK_EQUAL(ek.get_steps_per_md_step(2.), 1);
-    BOOST_CHECK_EQUAL(ek.get_steps_per_md_step(5.), 0);
   }
 
   {
@@ -216,7 +213,6 @@ BOOST_AUTO_TEST_CASE(ek_interface_none) {
     BOOST_CHECK_THROW(ek.on_node_grid_change(), NoEKActive);
     BOOST_CHECK_THROW(ek.on_timestep_change(), NoEKActive);
     BOOST_CHECK_THROW(ek.on_temperature_change(), NoEKActive);
-    BOOST_CHECK_THROW(ek.get_steps_per_md_step(1.), std::exception);
     auto const err_msg = std::string(NoEKActive().what());
     auto const ref_msg = std::string("EK not activated");
     BOOST_CHECK_EQUAL(err_msg, ref_msg);

testsuite/python/save_checkpoint.py

@@ -81,7 +81,8 @@
 if lbf_class:
     lbf_cpt_mode = 0 if 'LB.ASCII' in modes else 1
     lbf = lbf_class(
-        lattice=lb_lattice, kinematic_viscosity=1.3, density=1.5, tau=0.01)
+        lattice=lb_lattice, kinematic_viscosity=1.3, density=1.5,
+        tau=system.time_step)
     wall1 = espressomd.shapes.Wall(normal=(1, 0, 0), dist=1.0)
     wall2 = espressomd.shapes.Wall(normal=(-1, 0, 0),
                                    dist=-(system.box_l[0] - 1.0))