LBCPU: Fix rounding issue for grid size calculation (#3678)

Description of changes:
Replaces ceil by round in lbcpu n*agrid ~= box_l calculation + a test in which the old behavior resulted in a false positive

src/core/grid_based_algorithms/lattice.cpp

@@ -42,21 +42,23 @@ Lattice::Lattice(double agrid, double offset, int halo_size,
   /* determine the number of local lattice nodes */
   auto const epsilon = std::numeric_limits<double>::epsilon();
   for (int d = 0; d < 3; d++) {
-    grid[d] = static_cast<int>(ceil(local_box[d] / agrid));
+    grid[d] = static_cast<int>(round(local_box[d] / agrid));
     global_grid[d] = node_grid[d] * grid[d];
     local_index_offset[d] = node_pos[d] * grid[d];
   }
 
   // sanity checks
   for (int dir = 0; dir < 3; dir++) {
     // check if local_box_l is compatible with lattice spacing
-    if (fabs(local_box[dir] - grid[dir] * agrid) > epsilon * box_length[dir]) {
+    auto diff = fabs(local_box[dir] - grid[dir] * agrid);
+    if (diff > epsilon * box_length[dir]) {
       throw std::runtime_error(
           "Lattice spacing agrid[" + std::to_string(dir) +
           "]=" + std::to_string(agrid) + " is incompatible with local_box_l[" +
           std::to_string(dir) + "]=" + std::to_string(local_box[dir]) +
           " ( box_l[" + std::to_string(dir) +
-          "]=" + std::to_string(box_length[dir]) + " )");
+          "]=" + std::to_string(box_length[dir]) +
+          " ). Mismatch: " + std::to_string(diff));
     }
   }
 

testsuite/python/lb.py

@@ -315,6 +315,22 @@ def test_incompatible_agrid(self):
         print("End of LB error messages", file=sys.stderr)
         sys.stderr.flush()
 
+    def test_agrid_rounding(self):
+        """Tests agird*n ~= box_l for a case where rounding down is needed"""
+        system = self.system
+        old_l = system.box_l
+
+        n_part = 1000
+        phi = 0.05
+        lj_sig = 1.0
+        l = (n_part * 4. / 3. * np.pi * (lj_sig / 2.)**3 / phi)**(1. / 3.)
+        system.box_l = [l] * 3 * system.cell_system.node_grid
+        system.actors.add(self.lb_class(agrid=l / 31, dens=1,
+                                        visc=1, kT=0, tau=system.time_step))
+        system.integrator.run(steps=1)
+        system.actors.clear()
+        system.box_l = old_l
+
     @utx.skipIfMissingFeatures("EXTERNAL_FORCES")
     def test_viscous_coupling(self):
         v_part = np.array([1, 2, 3])