Fix regression in CPU LB thermalization.

src/core/grid_based_algorithms/lb.cpp

@@ -853,7 +853,7 @@ lb_thermalize_modes(Lattice::index_t index, const std::array<T, 19> &modes,
         rng_type{}(c, {{static_cast<uint64_t>(index), 2ul}}),
         rng_type{}(c, {{static_cast<uint64_t>(index), 3ul}})};
 
-    auto rng = [&](int i) { return uniform(noise[i / 4][i % 4]); };
+    auto rng = [&](int i) { return uniform(noise[i / 4][i % 4]) - 0.5; };
 
     return {/* conserved modes */
             {modes[0], modes[1], modes[2], modes[3],

testsuite/python/CMakeLists.txt

@@ -123,7 +123,10 @@ python_test(FILE widom_insertion.py MAX_NUM_PROC 1)
 python_test(FILE constant_pH.py MAX_NUM_PROC 4)
 python_test(FILE writevtf.py MAX_NUM_PROC 4)
 python_test(FILE lb_stokes_sphere.py MAX_NUM_PROC 4 LABELS gpu long)
-python_test(FILE lb_pressure_tensor_acf.py MAX_NUM_PROC 4 LABELS gpu long)
+if(NOT WITH_COVERAGE)
+  python_test(FILE lb_pressure_tensor.py MAX_NUM_PROC 1 LABELS gpu long)
+  python_test(FILE wang_landau_reaction_ensemble.py MAX_NUM_PROC 1 LABELS long)
+endif()
 python_test(FILE ek_fluctuations.py MAX_NUM_PROC 1 LABELS gpu)
 python_test(FILE ek_charged_plate.py MAX_NUM_PROC 1 LABELS gpu)
 python_test(FILE ek_eof_one_species_x.py MAX_NUM_PROC 1 LABELS gpu)
@@ -164,9 +167,6 @@ python_test(FILE analyze_distance.py MAX_NUM_PROC 1)
 python_test(FILE comfixed.py MAX_NUM_PROC 2)
 python_test(FILE rescale.py MAX_NUM_PROC 2)
 python_test(FILE accumulator.py MAX_NUM_PROC 4)
-if(NOT WITH_COVERAGE)
-  python_test(FILE wang_landau_reaction_ensemble.py MAX_NUM_PROC 1 LABELS long)
-endif()
 python_test(FILE array_properties.py MAX_NUM_PROC 4)
 python_test(FILE analyze_distribution.py MAX_NUM_PROC 1)
 python_test(FILE observable_profile.py MAX_NUM_PROC 4)

testsuite/python/lb_pressure_tensor_acf.py → testsuite/python/lb_pressure_tensor.py

@@ -32,7 +32,7 @@
 KT = 0.8
 
 
-class TestLBPressureACF:
+class TestLBPressureTensor:
     """Tests that the thermalized LB pressure auto correlation function
     is consistent with the chosen viscosity
     """
@@ -46,52 +46,108 @@ def tearDown(self):
         self.system.actors.clear()
         self.system.thermostat.turn_off()
 
-    def test(self):
-        # setup
+    def sample_pressure_tensor(self):
+        # Setup
         system = self.system
         lb = self.lb_class(agrid=AGRID, dens=DENS, visc=VISC,
                            tau=TAU, kT=KT, seed=SEED)
         system.actors.add(lb)
-        system.thermostat.set_lb(LB_fluid=lb, seed=2)
+        system.thermostat.set_lb(LB_fluid=lb, seed=SEED + 1)
 
         # Warmup
         system.integrator.run(500)
 
-        # sampling
-        steps = 50000
-        p_global = np.zeros((steps, 3, 3))
-        p_node = np.zeros((steps, 3, 3))
+        # Sampling
+        self.p_global = np.zeros((self.steps, 3, 3))
+        self.p_node0 = np.zeros((self.steps, 3, 3))
+        self.p_node1 = np.zeros((self.steps, 3, 3))
 
-        node = lb[0, 0, 0]
+        node0 = lb[0, 0, 0]
+        node1 = lb[3 * [N_CELLS // 2]]
 
-        for i in range(steps):
-            p_node[i] = node.pressure_tensor
-            p_global[i] = lb.pressure_tensor
+        for i in range(self.steps):
+            self.p_node0[i] = node0.pressure_tensor
+            self.p_node1[i] = node1.pressure_tensor
+            self.p_global[i] = lb.pressure_tensor
 
             system.integrator.run(2)
 
+    def assert_allclose_matrix(self, x, y, atol_diag, atol_offdiag):
+        assert x.shape == y.shape
+        n = min(x.shape)
+        mask_offdiag = ~np.identity(n, dtype=bool)
+
+        np.testing.assert_allclose(np.diag(x), np.diag(y), atol=atol_diag)
+        np.testing.assert_allclose(
+            x[mask_offdiag],
+            y[mask_offdiag],
+            atol=atol_offdiag)
+
+    def test_averages(self):
+        # Sound speed for D3Q19 in LB lattice units
+        c_s_lb = np.sqrt(1 / 3)
+        # And in MD units
+        c_s = c_s_lb * AGRID / TAU
+
+        # Test time average of pressure tensor against expectation ...
+        p_avg_expected = np.diag(3 * [DENS * c_s**2 + KT / AGRID**3])
+
+        # ... globally,
+        self.assert_allclose_matrix(
+            np.mean(self.p_global, axis=0),
+            p_avg_expected, atol_diag=c_s_lb**2 / 6, atol_offdiag=c_s_lb**2 / 9)
+
+        # ... for two nodes.
+        for time_series in [self.p_node0, self.p_node1]:
+            self.assert_allclose_matrix(
+                np.mean(time_series, axis=0),
+                p_avg_expected, atol_diag=c_s_lb**2 * 10, atol_offdiag=c_s_lb**2 * 6)
+
         # Test that <sigma_[i!=j]> ~=0 and sigma_[ij]=sigma_[ji]
-        tol_global = 4 / np.sqrt(steps)
+        tol_global = 4 / np.sqrt(self.steps)
         tol_node = tol_global * np.sqrt(N_CELLS**3)
 
-        # check single node
+        # check for the two sampled nodes
         for i in range(3):
             for j in range(i + 1, 3):
-                avg_ij = np.average(p_node[:, i, j])
-                avg_ji = np.average(p_node[:, i, j])
-                self.assertEqual(avg_ij, avg_ji)
+                avg_node0_ij = np.average(self.p_node0[:, i, j])
+                avg_node0_ji = np.average(self.p_node0[:, j, i])
+                avg_node1_ij = np.average(self.p_node1[:, i, j])
+                avg_node1_ji = np.average(self.p_node1[:, j, i])
+
+                self.assertEqual(avg_node0_ij, avg_node0_ji)
+                self.assertEqual(avg_node1_ij, avg_node1_ji)
 
-                self.assertLess(avg_ij, tol_node)
+                self.assertLess(avg_node0_ij, tol_node)
+                self.assertLess(avg_node1_ij, tol_node)
 
         # check system-wide pressure
         for i in range(3):
             for j in range(i + 1, 3):
-                avg_ij = np.average(p_global[:, i, j])
-                avg_ji = np.average(p_global[:, i, j])
+                avg_ij = np.average(self.p_global[:, i, j])
+                avg_ji = np.average(self.p_global[:, j, i])
                 self.assertEqual(avg_ij, avg_ji)
 
                 self.assertLess(avg_ij, tol_global)
 
+
+class TestLBPressureTensorCPU(TestLBPressureTensor, ut.TestCase):
+
+    def setUp(self):
+        self.lb_class = espressomd.lb.LBFluid
+        self.steps = 5000
+        self.sample_pressure_tensor()
+
+
+@utx.skipIfMissingGPU()
+class TestLBPressureTensorGPU(TestLBPressureTensor, ut.TestCase):
+
+    def setUp(self):
+        self.lb_class = espressomd.lb.LBFluidGPU
+        self.steps = 50000
+        self.sample_pressure_tensor()
+
+    def test_gk_viscosity(self):
         # Check that stress auto correlatin matches dynamic viscosity
         # eta = V/kT integral(stress acf)
         all_viscs = []
@@ -100,8 +156,9 @@ def test(self):
 
                 # Calculate acf
                 tmp = np.correlate(
-                    p_global[:, i, j], p_global[:, i, j], mode="full")
-                acf = tmp[len(tmp) // 2:] / steps
+                    self.p_global[:, i, j],
+                    self.p_global[:, i, j], mode="full")
+                acf = tmp[len(tmp) // 2:] / self.steps
 
                 # integrate first part numerically, fit exponential to tail
                 t_max_fit = 50 * TAU
@@ -116,7 +173,7 @@ def f(x, a, b): return a * np.exp(-b * x)
 
                 integral = numeric_integral + tail
 
-                measured_visc = integral * system.volume() / KT
+                measured_visc = integral * self.system.volume() / KT
 
                 self.assertAlmostEqual(
                     measured_visc, VISC * DENS, delta=VISC * DENS * .15)
@@ -126,21 +183,6 @@ def f(x, a, b): return a * np.exp(-b * x)
         self.assertAlmostEqual(np.average(all_viscs),
                                VISC * DENS, delta=VISC * DENS * .07)
 
-# DISABLE CPU TEST UNTIL #3804 IS SOLVED
-# class TestLBPressureACFCPU(TestLBPressureACF, ut.TestCase):
-#
-#    def setUp(self):
-#        self.lb_class = espressomd.lb.LBFluid
-#
-#
-
-
-@utx.skipIfMissingGPU()
-class TestLBPressureACFGPU(TestLBPressureACF, ut.TestCase):
-
-    def setUp(self):
-        self.lb_class = espressomd.lb.LBFluidGPU
-
 
 if __name__ == "__main__":
     ut.main()